The reader will find extensive studies in Lannes et al. (1987a-b) where the method summarized here is presented in a closed form, and some astrophysical applications in Fraix-Burnet et al. (1989), Roques et al. (1993) and Bouyoucef et al. (1994). To a first approximation, the experimental data are related to the original object , the intensity of the source at some high level of resolution, by an experimental transfer relation: where is the Point Spread Function. The error term includes random or systematic errors (e.g., errors on the determination of the PSF, linearity assumption of the imaging relation, image sampling) and signal-uncorrelated random noise (telescope, detectors, atmosphere, guiding, etc.).
The accuracy of the approximation considered when writing the convolution equation is controlled in Fourier space by a pointwise image-error bound such that . Moreover, can be regarded as a member of a certain family of objects. For each spatial frequency , it is possible to exhibit a suitable upper bound of the Fourier transform of : .
Due to the noise-type and systematic errors characterized by the pointwise Signal-to-Noise Ratio (SNR) in the frequency space: it is preferable to give up trying to determine at its highest level of resolution. One then defines the ``object to be reconstructed'' as a smoothed version of by a relation of the form , where is some synthetic transfer function small in the mean square sense outside some extent . This domain is the frequency coverage to be synthesized and regularizes the effective support of the transfer function . The choice of its diameter is of fundamental importance because it is closely related to the resolution limit of the reconstruction process. Intuitively, the greater is with respect to , the greater is the gain in resolution, but at the same time, the less stable is the deconvolution process. The problem is then to define the best compromise between the resolution to reach and the stability of the solution. The support of is contained in some finite region a priori known, which size and shape, determined in an interactive manner, will prove to play, together with , an essential part.
This filter fulfills three conditions:
A first approximation of of the object to be reconstructed is: where is some threshold value of the order of unity.
In the frequency domain where the SNR is considered as good (), is more or less reliable. One then assigns to the data a weight depending on the values of SNR. Let be the threshold value beyond which the whole information contained in is considered as entirely reliable. Let be the ``weight function'' characterizing the weight attached to this information. The function is defined as follows: One supposes that outside . So defined, characterizes the way of grounding the reconstruction on the information. Our deterministic procedure based on a least squares minimization defines the reconstructed object as the function which minimizes the functional:
The stability of the reconstruction problem is conditioned by the smallest eigenvalue of the imaging operator. It can be analytically estimated by examining some physical parameters: , characteristic functions of , and related to the choice of and to SNR. Indeed, this eigenvalue is a function of the ``interpolation parameter'' : characterizing the amount the interpolation to be performed both is real and in Fourier spaces. One has the following relation: where the ' depend on and and are of the kind ``moment of inertia'' relatively to and . This equation provides useful approximation of the minimum eigenvalue of the imaging operator occuring in the expression of an upper limit of the reconstruction error:
When this analysis is implemented before the reconstruction, one has only an estimation of and the estimation of . The upper bound of the quadratic reconstruction error is given by: In these two expressions, and can be regarded as errors resulting from the noise analysis and from systematic errors related to the choice of .
Note that these majorants correspond to the case where the overall error would be concentrated in the eigenspace associated with the minimum eigenvalue, what is not the case in practice (the error is generally distributed amongst all the eigenspaces).
Once the stability conditions are fulfilled, the solution can then be obtained by several well known minimization iterative methods. We apply the conjugate gradients method (Hestenes et al., 1952) because a particular implementation allows to compute the spectral decomposition of the imaging operator (then a more precise error can be calculated). The convergence is superlinear, and the least squares solution is reached in a number of iterations of the order of the number of degrees of freedom of the reconstruction process.