Restoration Based on Adaptively Weighted Regularization

Since linear shift invariant algorithms only lead to an unsatisfactory compromise between resolution and smoothness of the estimate (Andrews and Hunt 1977), we derive an adaptive method. This is done by incorporation of the inverse weighting matrices and into the quadratic error term . The Toeplitz matrix determines a linear shift invariant regularization function:

The optimization of and is carried out by minimization of the expectation of the squared distance between the original signal and the estimated signal:

Application of matrix differentiation rules and matrix algebra leads to the same matrix equation in both cases which means that the results are not independent of each other:

Since there are infinitely many solutions we are free to select one of the matrices such that the simplest solution is obtained:

The matrices are determined by the statistics of the original signal and the noise which are described by the autocorrelation matrices. If Eq. (7) is inserted into Eq. (3), the resulting estimate is similar to the estimate provided by a shift variant Wiener filter (Andrews and Hunt 1977). If the signal and the noise are Gaussian distributed it is also similar to a Maximum a Posteriori estimate (Andrews and Hunt 1977). Theoretically this is the optimal least squares solution of our restoration problem. But there are two reasons why it does not work in practical applications. First, the autocorrelation matrix of the original signal must be known quite accurately. Therefore, the algorithm is extremely sensitive against modeling errors and usually can not be estimated from the data at hand with sufficient accuracy. Second, its size is prohibitively large. In the case of a 512512 pixel image the size of the matrix is 262144262144, which exceeds even the memory capacity of today's computers.

A suboptimal but practical algorithm is obtained if we replace the general weighting matrices and by the diagonal matrices and . Now the number of parameters of the statistical model is drastically reduced and this makes the algorithm much more robust. Eq. (2) and Eq. (3) must be replaced by

The optimization of and is again carried out by minimization of the expectation of the squared distance between the original signal and the estimated signal:

Application of matrix differentiation rules and some matrix algebra show that the diagonal elements of the optimal diagonal matrices must be equal to the diagonal elements of the previously determined optimal general weighting matrices so that

Now we can derive a linear adaptively regularized restoration algorithm. Eq. (3) is rearranged, and using components we obtain a linear system of equations:

is determined from the local noise variance. Basically this means that in the case of HST data, the Poisson distributed noise is approximated by Gaussian distributed noise with signal dependent variance:

is estimated from the measured data. Since noise is Poisson distributed its variance is proportional to the local signal value. The HST point spread function leads to strongly smoothed images so that some additional smoothing using the operator sm causes only minor errors in the estimation of the underlying noise-free signal. Therefore we estimate:

The reciprocals of the diagonal elements of the inverse weighting matrix determine the optimal weighting function :

Obviously the weighting function is matched to the regularization function in a remarkably simple way. Unless a simulation is performed, is not known a priori! Therefore, an estimate must be determined from the available data. Before the iteration starts the measured data themselves provide the best available estimate of the original signal. If we assume that during the iteration the estimate converges to the original signal, estimates of the optimal weighting function can be determined successively from intermediate results of the iteration. In practical applications this tricky approach converges. The algorithm is self-updating at least up to a certain degree (Bundschuh 1992):