130 CHAPTER 8. DECONVOLUTION METHODS
All matrices in this set are of rank at most ` and have a prescribed eigenspace. The minimization problem (8.5) withA:“A`simplifies to the`-dimensional minimization problem with diagonal matrices,
min
y`PR`
Λ`PR`ˆ`
diagonal
!
}Λ`y`´˜bδ`}2F `α1}y`´y0,`}2F `α2}Λ`´Λ0,`}2F )
. (8.10)
Letty˚`,Λ˚`u minimize (8.10). The restored image then is given by
x˚“V`U0,`y˚`. (8.11)
In many applications, the matrix A0 has a structure that allows fast evaluation of matrix-vector products. For instance, when bδ represents an image in two space-dimensions that has been contaminated by space-invariant Gaussian blur, thenA0 can be chosen to be a symmetric block Toeplitz matrix with Toeplitz blocks. The number of steps`of the Lanczos process should be large enough so that the largest eigenvalues and eigenvectors of A0 can be approximated fairly accurately by Ritz values and Ritz vectors defined by the Lanczos decomposition (8.6).
Minimizing (8.10) leads to ` decoupled one-dimensional minimization problems (see [35]). From these problems, we can deduce λj in dependence of y. Under certain conditions, substituting the expression for λ in the corresponding equation leads to an equation ppyq “0 for a polynomial of degree five. Having computed its zeros ypkqj , k “ 1,2, . . . ,5, we determine the associated λ-values λpkqj . We are only interested in the real zeros, since our image is real valued. The pointspyjpkq, λpkqj qassociated with the real zeros yjpkq of p are possible minima of our one-dimensional function. Evaluating the function with all these points determines all solutions of (8.10) with α1 “ α and α2 “αt.
CHAPTER 8. DECONVOLUTION METHODS 131
Figure 8.1: Log plot of both the true star cluster image from [113] and the same image blurred using the true PSF associated to an observation. The log plot is used to exaggerate the star brightness.
package [113].
Let us recall that the true PSF is unknown and only an estimate can be obtained through the reconstruction methods presented in Chapter 6 and 7. Usually, the PSF is not given as a matrix, but the action of convolution with the PSF induces a matrix with dimensions corresponding to the image dimensions. Therefore, we need to com-pute the corresponding matrix.
We note that in the following experiments, the initial approximation for the image is always the zero vector, i.e., an empty image. Thus for a reconstructed image x, wer denote as the relative error norm}xtrue´xr}{}xtrue}. For the operator reconstruction, the initial approximation is the given operator A0. Thus for a reconstructed operator Arwe denote as the relative operator reconstruction error}Ar´Atrue}{}A0´Atrue}.
8.4.1 Experiment setting
Before being able to use our algorithm for image improvement, we have to artificially blur an image with the true PSF corresponding to the observation. As in the numerical tests of Section 6.5 and 7.4, we use as a benchmark a high quality approximation of the PSF influenced by atmospheric turbulence from OCTOPUS, and we use this high quality PSF as the true one. Before we are able to use the PSFs in a (de-)convolution step, we need to make sure that the size of one pixel in the science image is the same as the size of one pixel in both PSFs, the true and the reconstructed one. This can be achieved by using theMATLAB-functioninterp2. For the planned MICADO instrument on the ELT, the resolution of the science images will be 4 mas{pixel. We scale the used PSFs to this resolution and assume that one pixel in the science image is 4mas
132 CHAPTER 8. DECONVOLUTION METHODS
wide.
As the Lanczos-process is a matrix-based method, we need to transform the PSF into the induced matrix. This matrix can be obtained by convolving the PSF with all canonical basis vectors from Rpq (i.e., the space of pˆq images). With the help of this step, we noticed immediately that the PSFs from this application do not actually induce symmetric matrices. This is due to the presence of small perturbations from noise in the PSF.
However, one could think of using a radial average to get a symmetrized version of the PSF, e.g., using the MATLAB-package [108]. The underlying formula is given as
PSFsymprq “ 1 2π
ż2π 0
PSFpr, θqdθ, @r ą0,
where PSFpr, θq is a representation of PSF in polar coordinates. This process is a natural extension of our reconstruction method and it is also common to look only at the radial average of a PSF in real applications (cf., e.g., [107]). The idea why this might be natural is that in the used AO system all wavefront distortions are small and have an expected value of zero, clearly indicating that the PSF associated to a run of the AO system should not be too far of from the optimal PSF of the telescope which is symmetric. Acquiring an image takes a long time compared to measuring the wave-front distortions, thus speckles occurring due to the distortions in a short time frame and measurement errors will average out. The reconstructed PSF is in center direction of the telescope, thus no built-in asymmetric behavior will show up in the center of our image. Additionally, as the field of view is small, the PSF will be (almost) stable across the field. Last, but not least, the reconstructed PSF is an estimate for the true PSF anyway, so radially averaging it will introduce some (small) additional error. As we aim for a blind deconvolution, this additional error just changes our initial guess slightly and hopefully does not change the result of the deconvolution process. Note that this process is applied before the matrix representation of the PSF is computed.
Computing the matrix first and then symmetrizing it would result in different errors.
In Figure 8.3, we show the symmetrized versions of the PSFs shown in Figure 8.2.
The distortion created by the action of each radially symmetrized PSF can be seen in Figure 8.4.
8.4.2 Comparing the induced matrices
Having at hand our true and false (approximated) radially symmetric PSFs, we can induce two symmetric matrices from R1024ˆ1024 and study their properties. One sees immediately from a plot of their entries in Figure 8.5 that the induced matrices are highly structured. Indeed, a look at the entries indicates that they are Toeplitz matri-ces. This makes sense, as the matrices are induced by discrete convolution operations.
CHAPTER 8. DECONVOLUTION METHODS 133
Figure 8.2: The true PSF obtained by OCTOPUS (left) and the PSF reconstructed from WFS data computed in Section 6.5 (right).
Figure 8.3: Radially symmetrized versions of the true PSF obtained by OCTOPUS (left) and the PSF reconstructed from WFS data computed in Section 6.5 (right), using the the PSFs of Figure 8.2.
Figure 8.6 displays the log-plot of the positive eigenvalues in descending order of the matrices displayed in Figure 8.5, and Figure 8.7 shows four normed eigenvectors associated with the largest eigenvalues of both matrices. One can see that the dominant eigenvectors look fairly similar and that the eigenvalues decay rapidly. We deduce that already a few iterations with the Lanczos process applied to this blurring matrix has the potential of giving accurate approximations of the dominant eigenvectors of the matrix induced by the true PSF. Thus the method proposed in this section can be used to determine suitable eigenvalues.
134 CHAPTER 8. DECONVOLUTION METHODS
Figure 8.4: Log plot of the star cluster image, from [113], blurred using the true and the reconstructed PSF. Note that the left image is the same as the right one in Figure 8.1 for demonstration purpose.
200 400 600 800 1000
100 200 300 400 500 600 700 800 900
1000 -7
-6 -5 -4 -3 -2 -1
200 400 600 800 1000
100 200 300 400 500 600 700 800 900 1000
-6 -5 -4 -3 -2 -1
Figure 8.5: Log plot of the absolute values of the matrices determined by radially symmetrized true and reconstructed PSFs.
8.4.3 Deconvolving images from AO assisted observations
We want to illustrate the effect of our algorithm on the image. For this illustration, we take a symmetrized version of the true PSF from an SCAO simulation run in OC-TOPUS and convolve a star image, taken from [113]. From the same SCAO run in OCTOPUS, we get an estimate of the PSF through the reconstruction methods of Chapter 6, which we symmetrize for the reasons discussed above. The matrix induced by this symmetrized reconstructed PSF is then used as matrixA0 in our algorithm.
For the SCAO run, we use the setting from Section 6.5.3. We consider a high guide star flux, i.e.,nph “1000, and the PSF reconstructed as in Section 6.5.4, withρ“0.125m.
CHAPTER 8. DECONVOLUTION METHODS 135
0 200 400 600 800 1000
10-4 10-3 10-2 10-1 100 101 102
Postive eigenvalues of induced matrices
True PSF False PSF
Figure 8.6: Log plot of the positive eigenvalues of the matrices induced by applications of the true and false PSFs..
0 200 400 600 800 1000 1200
-0.05 -0.045 -0.04 -0.035 -0.03 -0.025 -0.02 -0.015
-0.01 First eigenvector
True PSF False PSF
0 200 400 600 800 1000 1200
-0.06 -0.04 -0.02 0 0.02 0.04
0.06 Second eigenvector
True PSF False PSF
0 200 400 600 800 1000 1200
-0.06 -0.04 -0.02 0 0.02 0.04
0.06 Third eigenvector
True PSF False PSF
0 200 400 600 800 1000 1200
-0.06 -0.04 -0.02 0 0.02 0.04
0.06 Fourth eigenvector
True PSF False PSF
Figure 8.7: Comparison of first four unit eigenvectors of matrices induced by the two PSFs.
As the deviation between our reconstructed PSF and the true PSF is only 2% in the peak, we expect only a small deviation in the matrices induced by symmetrized ver-sions of these PSFs. We start from a 2D image of the PSF calculated by the means of Chapter 6. Using theMATLAB function radialavg.m from [108], we obtain a radially symmetrized version of our PSF. This radially symmetric PSF can then be transformed into the matrixA0 by applying the MATLAB-in-built function imfilter.
For these PSFs and the blurred image in Figure 8.1 (right), we carry out tests with 10 and 20 Lanczos iterations, and investigate for various values oft and α the quality of the reconstructed image and low-rank operator reconstruction. For the image, we simply compute the relative error for all tests and these are shown in Figure 8.8.
Judging the quality of the low-rank operator reconstruction is a bit more difficult.
There is no explicitly available true operator with which to compare, and for an im-age of this size, the induced matrix representation of the PSF applications would be prohibitively expensive to compute.
In Figure 8.9, we present a reconstruction usingj “5 Lanczos iterations and
param-136 CHAPTER 8. DECONVOLUTION METHODS
j= 10
t
0 1000 2000 3000 4000 5000
α
0 100 200 300 400 500 600 700 800 900
1000 0.978
0.98 0.982 0.984 0.986 0.988 0.99 0.992 0.994 0.996 0.998
j= 20
t
0 1000 2000 3000 4000 5000
α
0 100 200 300 400 500 600 700 800 900 1000
0.98 0.985 0.99 0.995
Figure 8.8: Relative image reconstruction norm errors for different combinations oft and α for 10 and 20 Lanczos iterations.
eters α1 “ 0.1, α2 “ tα1 and t “ 10. For comparison, we compute the relative error norm
}Itrue´Irec}2
}Itrue´Io}2
,
where Itrue is the original image, Irec is the reconstruction obtained by applying our method to the observed, blurred image Io. For j “5 Lanczos iterations this ratio is 0.98, i.e., the norm error is improved by only 2 % even though the visual impression of the reconstructed image is much better. Increasing toj “20 gives a ratio of 0.65, i.e., an improvement by 35 %, but the reconstructed image shows higher values in the parts of the image which are supposed to be black. This will require further investigation in order to find the best possible setting for real life applications.