An approach to computational problems associated with generation and estimation of large Gaussian fields is studied. Fast algorithms for matrix operations on circulant matrices are presented, and a connection between such matrices and covariance matrices of Gaussian fields is established. Based on this approach, a model for discrete spatial data is introduced, extending the work of Nott and Wilson (1997). We assume that discrete random fields are obtained by clipping a stationary zero mean Gaussian random field at several fixed levels. The model is defined by this set of levels, a choice of a family of covariance functions for the Gaussian field, and a parameter vector specifying a particular covariance function within the family. For this model, the Stochastic Expectation-Maximization algorithm for estimating the covariance parameter vector is presented. The algorithm includes conditional generation of Gaussian fields given that components fall within specified intervals; this is achieved by the Gibbs sampler--a Markov Chain Monte Carlo technique. The precision of the algorithm--understood in terms of the variance of the resulting estimator of the correlation function--is compared to that of estimating the parameter directly from the Gaussian data by the Maximum Likelihood method. For this purpose, the Fisher information matrix in the Gaussian model is computed, the asymptotic distribution of the MLE estimator of the correlation parameter is established, and simulations are performed to compare the empirical variances of the MLE and several SEM estimators (the latter based on various quantizations) to the variance predicted by the theory, and to each other.