Restoration of Images with Wavefront Aberrations

TL;DR

This paper introduces a new wavefront correction algorithm for coherent image restoration, effectively handling severe distortions in optical systems like astronomy and microscopy, outperforming traditional methods.

Contribution

A novel wavefront correction algorithm based on evolving Gerchberg-Saxton and Hybrid-Input-Output methods for coherent image restoration.

Findings

High-quality restoration with severe wavefront distortions

Effective on simulated and real microscopic images

Improves coherence-based optical imaging quality

Abstract

This contribution deals with image restoration in optical systems with coherent illumination, which is an important topic in astronomy, coherent microscopy and radar imaging. Such optical systems suffer from wavefront distortions, which are caused by imperfect imaging components and conditions. Known image restoration algorithms work well for incoherent imaging, they fail in case of coherent images. In this paper a novel wavefront correction algorithm is presented, which allows image restoration under coherent conditions. In most coherent imaging systems, especially in astronomy, the wavefront deformation is known. Using this information, the proposed algorithm allows a high quality restoration even in case of severe wavefront distortions. We present two versions of this algorithm, which are an evolution of the Gerchberg-Saxton and the Hybrid-Input-Output algorithm. The algorithm is verified on simulated and real microscopic images.

Equations

$$M=I\otimes B+N.$$

$$A_m=U\otimes B_a+N_a,$$

$$B=(|B_a|{)}^2.$$

$$o(n)=\mathcal{F}(A)(n)=\int A(k)\exp (-j2\pi kn)\mathrm{d}k.$$

$$A_i(k)$$

$$A_i^{\prime }(k)$$

$$o_i^{\prime }(n)$$

$$o_{i+1}(n)$$

$$o_{i+1}(n)=\{\begin{array}{ll}{o}_i^{\prime }(n)& n\in V\\ o_i(n)-\beta o_i^{\prime }(n)& n\notin V,\end{array}$$

$$P_A(y)=\{x\in A,\Vert x-y\Vert =\Vert y-A\Vert \}.$$

$$o_{i+1}(n)=\{\begin{array}{ll}{P}_A(o_i)(n)& n\in V\\ o_i(n)-\beta P_A(o_i)(n)& n\notin V.\end{array}$$

$$p_d=\exp (j\frac{k}{2r}(x^2+y^2)),$$

$$A_i(k)$$

$$f_i(n)$$

$$f_i^{\prime }(n)$$

$$A_i^{\prime }(k)$$

$$o_i(n)$$

$$o_{i+1}(n)$$

$$E_{o,i}^2={\int }_n(o_{i+1}(n)-o_i^{\prime }(n){)}^2\mathrm{d}n.$$

$$E_{f,i}^2={\int }_n(f_i(n)-f_i^{\prime }(n){)}^2\mathrm{d}n.$$

$$E_{f,i}^2={\int }_k(A_i(k)p_s^{-1}-A_i^{\prime }(k)p_s^{-1}{)}^2\mathrm{d}k.$$

$$E_{f,i}^2={\int }_k(A_i(k)-A_i^{\prime }(k){)}^2\mathrm{d}k.$$

$$E_{f,i}^2={\int }_n(o_i(n)-o_i^{\prime }(n){)}^2\mathrm{d}n.$$

$$(o_i(n)-o_i^{\prime }(n){)}^2\le (o_{i+1}(n)-o_i^{\prime }(n){)}^2,$$

$$E_{o,i}^2\le E_{f,i}^2.$$

$$(f_{i+1}(n)-f_{i+1}^{\prime }(n){)}^2\le (f_i(n)-f_i^{\prime }(n){)}^2.$$

$$E_{f,i+1}^2\le E_{0,i}^2.$$

$$E_{f,i+1}^2\le E_{0,i}^2\le E_{f,i}^2$$

$$E_{f,i+1}^2\le E_{f,i}^2.$$

$$o_{i+1}(n)$$

$$=1_V{P}_f(o_i)(n)+(1-1_V(o_i-\beta P_f(o_i))(n)$$

$$=1_V((1+\beta )P_f(o_i)-o_i)+(I-P_f)(o_i)(n)$$

$$=(P_o((1+\beta )P_f-I)+I-P_f)(o_i)(n)$$

Full text

In this paper a novel wavefront correction algorithm is presented, which allows image restoration under coherent conditions. In most coherent imaging systems, especially in astronomy, the wavefront deformation is known. Using this information, the proposed algorithm allows a high quality restoration even in case of severe wavefront distortions. We present two versions of this algorithm, which are an evolution of the Gerchberg-Saxton and the Hybrid-Input-Output algorithm. The algorithm is verified on simulated and real microscopic images.

I Introduction

In this paper an algorithm based on the Gerchberg-Saxton algorithm is presented, which restores images disturbed by optical aberrations in case of coherent or partially coherent illumination. Prior deconvolution algorithm are not capable of this task. This algorithm, the wavefront correction algorithm (WFC), can be applied directly on imagery from coherent microscopy and astronomy and deal with the restoration of images blurred by wavefront errors, caused frequently by imperfect optical elements, alignment or by defocus. An example of blur caused by defocus (a spherical wavefront aberration) and the restored image is shown in Figure 1.

The removal of these aberrations in the context of coherent illumination is important, as coherent imaging conditions occur even in standard microscopic settings (see experimental Section III-B) and the light from the stars in astronomic imaging fulfills the conditions of spatial coherence.

Perfect aberration free optics are very costly and hard to design as the initially blurry images of the Hubble space telescope have shown [1]. The wavefront aberration in diffraction limited, sometimes called ”Fourier optics”, needs to be less than $\lambda/5$. Sometimes a total aberration of even less than $\lambda/10$ is required. This is hard to achieve, especially in systems composed of many elements.

The proposed algorithm can restore high quality images from imperfect coherent optical systems, such as Hubble space telescope before repairs. Restoration is possible even if the aberrations are constantly changing, like in earth-bound astronomy, where complex adaptive optics with deformable mirrors are necessary to compensate atmospheric turbulence [2]. A major advantage is that this algorithmic approach is not limited in modulation speed and magnitude in contrast to deformable mirrors.

Standard image processing deconvolution algorithms such as the Wiener restoration filter [3], the Richardson-Lucy algorithm [4] [5] and modern blind deconvolution algorithms [6] are based on the following image formation model, where an image $M$ is formed from the undisturbed image $I$, blur kernel $B$ and additional noise $N$

[TABLE]

This model describes the case of incoherent imaging well.

In case of coherent imaging, as in astronomy or coherent microscopy, the image formation is described by a different model [7]

[TABLE]

where $U$ is the undisturbed complex amplitude and $B_{a}$ and $N_{a}$ amplitude blur kernel and noise. Cameras acquire an image by measuring the light intensity distribution over an image sensor or historically a photosensitive film. The magnitude of the amplitude distribution is the square root of the light intensity distribution. $A_{m}$ is the complex amplitude distribution of the optical wave in the image plane. The recorded intensity image contains no phase information, therefore the complex part of the amplitude is unknown and standard algorithms cannot be applied, as they cannot deal with interference (see Figure 5). Clearly incoherent imaging is linear in intensity, while coherent imaging is linear in amplitude. The blur kernels are related [7]:

[TABLE]

In the following we shall deal with coherent imaging only. Phase retrieval algorithms solve the task of restoring the unknown phase of a wavefront from measured intensity distributions in Fourier and object plane with additional spatial constraints in the object plane. A more detailed overview of phase retrieval algorithms is given in [8] and a performance comparison of various algorithms in [9]. Notable recent developments are [10] and [11]. In the area of coherent optics phase retrieval algorithms are widely used in many applications from crystallography, holography to ptychography and astronomy. An overview of the applications of the Gerchberg-Saxton (GS) algorithm and its variations is given in [12]. The classic phase retrieval algorithm is the GS-algorithm [13]. It is discussed in the following section, together with the Hybrid-Input-Output (HIO) algorithm for phase retrieval by Fienup [14].

They form the base of our algorithm, however it is important to note that this paper deals with a different task, namely the correction of wavefront distortions solely from a measured intensity distribution in the image plane. Additionally some estimate of the wavefront distortion is needed. This cannot be accomplished with phase retrieval or the prior non-coherent image restoration algorithms (see Figure 5).

Furthermore we use different constraints than prior phase retrieval algorithms, as we introduce an additional virtual focus plane, whose amplitude distribution we require to be real and positive. Good theoretical convergence properties of the algorithm were confirmed by its application on simulated and real complex images. A comparison between simulated and real images is performed on a standard bright field microscope with a microscopic lens. The illumination lamp was replaced by a standard LED, which showed satisfactory wavelength coherence. The required spatial coherence could be achieved by adjusting the illumination diaphragm.

I-A The Phase Retrieval Problem

The amplitude distributions across the real-space object plane $o$ and a Fourier-space plane $A$ are linked by the Fourier transform operator $\mathcal{F}$:

[TABLE]

The goal is to recover the phase of $o$ from a set of measurements, while fulfilling the constraints on the object plane and on the Fourier plane. The choice of constraints by Gerchberg in [13] is restricting the object magnitude in real- and Fourier-space to the measured magnitudes $|o_{0}(n)|$ and $|A_{0}(k)|$.

The algorithm starts with an arbitrary phase distribution in $o$ and loops over real-space and Fourier-space, while enforcing the constraints in every iteration $i$, (see Figure 2):

[TABLE]

A disadvantage of this basic algorithm is its slow convergence. This was greatly improved by Fienup in [14]. His Hybrid-Input-Output (HIO) algorithm is widely used today. This algorithm treats the application of Fourier domain constraints as a feedback system with the input $o_{i}$ in Eq. 5 and output $o^{\prime}_{i}(n)$ in Eq. 7. Furthermore an additional support constraint is added, which restricts the size of the reconstruction to a prior defined set of points. The object-space update Eq. 8 is changed to:

[TABLE]

where $V$ is a set of points, where the object constraints are valid and $\beta\in\mathrm{R}$ is a parameter influencing the convergence speed. A value commonly used is $0.7$ [12], we follow this convention. For an evaluation of the influence of this parameter see [14].

From an input image, that is violating a constraint, the enforcement of a constraint outputs a result which is the first valid input closest to that violating input. This application of a constraint can also be expressed in the framework of projections [15], where the application of a projection means the enforcement of a constraint. Let $A$ be the set of values for which a constraint holds, then we define a projection $P$ of $y$ onto $A$ as the set of values $x\in A$, which are closest to $y$:

[TABLE]

For practical purposes if the set has more than one value, we choose one of these points. For conditions for single-valuedness of the Fienup constraints and a more detailed introduction see [16]. The application of the constraint in Eq. 8 is defined as $P_{o}$ and the subsequent application of equations 5, 6 and 7 as $P_{A}$. The HIO-algorithm now reads:

[TABLE]

II Wavefront Correction Algorithm

The goal is to develop an image restoration algorithm, based solely on the knowledge of the intensity distribution of the distorted image.

The light distribution in the focal plane of a thin lens is the Fourier transform of the light distribution of the aperture (Fourier plane). Thus assuming that the input image is captured in perfect focus we can switch between this focal image plane and the Fourier plane by applying the Fourier transform and the inverse Fourier transform, as in the classic phase retrieval algorithms see Section I-A. However, if the image is not captured in perfect focus, this condition no longer holds. We need to include the additional wavefront deformation introduced by defocus or other optical aberrations.

For coherent light incident on a convex thin lens with radius $r$ (assuming the Fresnel approximation [17]) the phase delay is

[TABLE]

where $k$ is the wave number and $x$ horizontal and $y$ vertical distance from the center of the lens perpendicular to the optical axis. Thus the effect of a thin lens is equivalent to a spherical wavefront deformation, which has the same form as the defocus error. We generally describe these wavefront aberrations using Zernike polynomials [18].

In some applications, as in a typical microscopic setting, we have no sensors to measure the Fourier space magnitude.However we know that in an ideal coherent optical system the image of a real object is a real positive function with planar phase. This is a consequence of the linearity of an aberration free imaging system.

Therefore we can use two constraints: The first constraint is the measured amplitude distribution, which is the square root of the intensity distribution, the second one is the requirement for a real and positive amplitude distribution in the virtual ideal imaging plane, where all aberrations are removed. The application of phase aberration function $p_{s}$ or its inverse $p_{s}^{-1}$ means a phase shift on the corresponding complex amplitudes. It is defined analogous to the specific defocus phase shift in Eq. 12.

II-A WFC-GS

We introduce a virtual focus plane $f$, which can be reached from Fourier-space $A$ via the inverse Fourier transform. The following equations define the Gerchberg-Saxton (GS) version of the Wavefront Correction Algorithm (WFC-GS) in iteration step $i$, illustrated by Figure 3:

[TABLE]

The algorithm starts similar to the GS-algorithm with the measured $o_{0}$ as $o_{i}$. In the first step (Equation 13) the Fourier transform is applied. Then the wavefront deformation is applied by multiplying the Fourier space distribution with phase delay $p_{s}$ of this wavefront. Transferred to the virtual focal plane (14), the real and positive constraint of the distribution is enforced (15). Then we transform back into object-space (16, 17), enforce the object magnitude constraint (18) and start iteration $i+1$.

The algorithm starts in object space with the measured magnitude distribution.

For the practical application of the algorithm note that the knowledge of the exact wave number of Eq. 12 is not necessary, as it is only a linear factor with the radius. More practical is a measure of maximum phase difference dependent on the image size.

In [14] the weak convergence for the GS-algorithm is shown using an error reduction argument, we expand on this to prove the same property for the GS-version of the wavefront correction algorithm (WFC-GS).

In iteration $i$ the squared error against the real-space plane constraint over the entire image is

[TABLE]

Similar for the focus plane constraint the squared error is

[TABLE]

The squared error may also be seen as the energy of the error. With Parseval’s theorem it follows that

[TABLE]

Because of the energy conservation of optical light transport, we obtain

[TABLE]

This allows a second application of Parseval’s theorem:

[TABLE]

Because for any $n$

[TABLE]

it follows with equations 19 and 23 that

[TABLE]

With the same arguments for $f$, clearly

[TABLE]

From this it follows with equations 20 and 23 that

[TABLE]

Combined with equations 25 we see

[TABLE]

and

[TABLE]

Therefore the error decreases or stagnates with every iteration. This result shows the weak convergence, but does not guarantee that the algorithm always converges to a global minimum.

II-B WFC-HIO

As a further step we adapt the HIO-algorithm to focal-space constraints and wavefront correction. We denote the application of Eq. 18 as $P_{o}$ and the application of Eq. 13 to 17 as $P_{f}$. Because we lack information about the object support, we cannot directly apply the standard HIO-definition Eq. 11. Let $\boldsymbol{1}_{V}$ be the characteristic function of the set of points $V$ where the constraints (such as object support) are valid and $\boldsymbol{1}_{\boldsymbol{C}V}$ of its complement. With [16] we write the HIO-version of the wavefront correction algorithm (WFC-HIO) as:

[TABLE]

and can apply it on images with wavefront aberration in the following section, without requiring object support information.

III Experimental Results

To verify the algorithm both in correctness of the results and applicability, we apply it both on synthetic images with simulated aberrations and to real microscopic images.

III-A Evaluation on Synthetic Images

For the visual verification of restoration quality, we choose a test image with sufficiently high contrast and fine structure (see Figure 4a). As a first test synthetic defocus aberration is applied on the image by convolving it with the PSF generated from the spherical wavefront deformation in Figure 7a. The task of the algorithm is now to restore the undisturbed image from this aberrated image. We initialize the unknown phase in real-space with zero and apply the WFC-GS algorithm using the known value of the optical aberration. The synthetic image together with the very convincing restoration is depicted in Figure 1. The results in Figure 1 show a near perfect restoration even for this very strong aberration and thereby verify that the algorithm works as theoretically predicted.

Above we have shown that the GS-version of wavefront correction image restoration algorithm (WFC-GS) converges in the weak sense. For the Gerchberg-Saxton algorithm for phase retrieval quite high iteration counts of over $2000$ have been reported to overcome stagnation in the results [14], which makes the GS-algorithm unpractical and was the primary reason for the improvements developed by Fienup. To test the convergence of the presented algorithm, we measure the peak signal to noise ratio (PSNR) between the focal image in every iteration and the ideal original.

To verify that the algorithm is not limited to defocus, we apply a combination of synthetic spherical aberration, astigmatism and coma original image, by convolving it with the PSF generated from the Zernike polynomials of the aberration [18]. We start with this strongly distorted image in Figure 4b, where the wavefront aberration is depicted in Figure 7b. Very high PSNR values of over $40\,\text{dB}$ prove that the result of the algorithm in Figure 6a is not only visually convincing, but also leads to correct results, even in case of complex wavefront aberrations. This is better than the PSNR values of $35\,\text{dB}$ reported by [20] for high quality compression algorithms.

Comparing both algorithms Figure 8, note the logarithmic scale, we see that the WFC-HIO shows the expected superior convergence. The WFC-HIO reaches a PSNR of $38.3\,\text{dB}$ in iteration 100, while the GS-version requires twice the number of iterations for this value. The difference in iteration $200$ is $1.5\,\text{dB}$. Overall both algorithms show good performance and results.

The curves show a sub-linear rate of convergence and that many iterations are necessary for near-perfect results, however even after only 10 iterations the WFC-HIO shows good results, with a PSNR of $33.2\,\text{dB}$, see Figure 6b. For a $512\times 512$ pixel image, on an $4\,\text{GHz}$ desktop processor the time spent per iteration is $14\,\text{ms}$, dominated by the fast Fourier transform.

The results of the Richardson-Lucy algorithm and a modern deconvolution algorithm by Krishnan and Fergus [6] in Figure 5b confirm that on images with such strong aberrations incoherent deconvolution algorithms are not applicable.

III-B Evaluation on Real Images

For the evaluation on real data, we choose a standard transmissive bright field microscopic setting. A standard USAF 1951 target is back-illuminated with a $650\,\text{nm}$ LED. The transmitted light is captured by a standard microscopic lens, the eyepiece was replaced by a digital camera. Before the application of the algorithm, the gamma curve of the images is linearized and edge tapering to reduce ringing artifacts is applied. We use the Matlab edge tapering function [19], the same edge tapering technique as in [21]. Also the pixel offset and scale need to be adjusted.

For this imaging setup we acquired a series of images with $50\,\mu\text{m}$ steps of the sharpest image on the optical axis. In Figure 9a the sharp image of the target is shown, in Figure 9c the image at a distance of $1950\,\mu\text{m}$ from the focal plane. We choose such a strong defocus, to demonstrate the effectiveness of the algorithm.

The restoration is depicted in Figure 9d. The algorithm achieves very good results even for this very strong wavefront aberration. None of the numbers are identifiable in the input image, but good contrast and high readability in the restoration is achieved. We consider acquisition noise and non-perfect coherence of the LED illumination as the source of the ringing artifacts.

Although an LED is far less coherent than a laser, our simple experiment shows that it can be used for a coherent illumination in microscopic applications. For a detailed discussion on spatial and frequency coherence see [7].

A synthetic defocus image with the same defocus parameters as used for the restoration is shown in Figure 9b. This image is at first glance almost identical to the real measurement in Figure 9c, which again verifies the correctness of the measurement.

IV Conclusion

The presented algorithm allows the restoration of coherent images even in case of strong wavefront aberrations. This should permit an essential cost reduction in coherent optical systems, for example in astronomy and coherent microscopy. Superb images can be achieved even with low quality optical components.

In most applications the wavefront error of the optical system is known, it can be measured using an interferometer or a point light source. In astronomy, it is common to use an artificial laser guide star to measure the wavefront aberrations induced by the turbulence of the atmosphere [2]. Adaptive optics, such as deformable mirrors are then used to compensate these wavefront deformations. With the new algorithm the deformable mirrors and their limitations can be avoided. Microscopic applications of adaptive optics with deformable mirrors [22] should also benefit from these improvements. Moreover, the WFC-algorithm can improve the performance of optical systems, in all cases where the implementation of a deformable mirror is not feasible.

The high iteration speed, together with the restoration result with $10$ iterations of the WFC-HIO (Figure 6b) demonstrate that if decent, yet non-perfect restorations are acceptable, real-time application of this algorithm is possible. This speed already allows parameter optimization and fine-tuning of the wavefront aberration by measuring the image sharpness. There are many possibilities for a faster implementation of this algorithm, e.g. on a GPU. By changing the value of the spherical aberration, we should be able to simulate autofocus. Furthermore, scanning through different values of the spherical aberration could be an interesting feature in the evaluation of coherent microscopic samples. An obvious choice as general optimization parameters would be the coefficients of the Zernike polynomials.

Our future work will be focused on further development, optimization and practical implementation of the algorithm.

Acknowledgment

This work has partly been supported by the German Research Foundation (DFG) Cluster of Excellence FUTURE OCEAN under proposals CP1331 and CP1525.