where L is a nonsingular univariate discrete Laplacian, L3 = L ⊗ I3 × 3, where ⊗ denotes the Kronecker product, W is a certain weight matrix defined in the weighted minimum norm solution, Λ is a diagonal matrix of regularizing parameters, and parameters τ and α are introduced. τ controls the amount of smoothness and α the relative importance of each grid point. Estimators are calculated iteratively, starting with a given initial estimate D0 (which may be taken to be D^LOR), Λ i is estimated from Di - 1, then D i from Λ i until one of them converges.