Open Access
Add to BibSonomy
Abstract
We quantitatively investigate how spatially varying deblurring algorithms can improve the imaging performance of hybrid optical/digital systems affected by field aberrations. To this end, we validate a theoretical model of the maximal gain that linear and spatially varying deblurring can bring to any given lens, and derive a practical algorithm to implement this type of deblurring with low computational complexity. The results demonstrate the usefulness to properly coordinate and balance the roles of the imaging optical system and raw image post-processing: optimal final imaging quality can be obtained by a lens that has been optically designed to reduce field aberrations at the price of lower average raw optical quality, associated with a fast and “slightly” spatially varying piecewise Wiener deconvolution algorithm.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Lenses with large field of view (FoV) are affected by field aberrations which make image quality spatially inhomogeneous. This type of field aberrations could be reduced by adding a certain number of optical elements, but this may increase the cost significantly. However, nowadays, it is unusual to operate an optical system without digital post-processing [1–7]. Digital processing can thus be leveraged to take on a portion of the burden of making image quality homogeneous over the FoV, i.e. to relax the quality level required on lens optical performance which may reduce at the same time the lens complexity and cost [8–14]. The question is to properly coordinate and balance the roles of the optical part and the processing part towards this goal.
In this article, we address this issue by considering the example of two relatively simple Cooke triplet lenses with large field of view (FoV) that have been optimized with the commercial lens design software CodeV in two different ways. The first triplet is optimized with a classical lens design criterion, and presents significant field aberrations, but a good quality on-axis. The second one is co-designed by minimizing the mean square error (MSE) of the final image after post-processing with a single deconvolution filter. The quality of the raw image is more uniform, but lower on average over the FoV, and in particular on-axis. Our goal is to improve the imaging performance of these two lenses by compensating for raw image field aberrations with deblurring algorithms adapted to spatially varying point spread function (PSF). In particular, we propose a theoretical model of the maximal gain that linear and spatially varying deblurring can bring to any given lens, and derive a practical algorithm to implement this type of deblurring with low computational complexity.
The article is organized as follows. In Sec. 2, we describe the two triplet lenses studied in this paper and theoretically evaluate the imaging performance provided by an ideal reconstruction algorithm taking fully into account the continuous spatial variation of the PSF. We compare this performance with that of a spatially uniform average Wiener filter for the two considered triplets. In Sec. 3, we introduce and theoretically evaluate the idea of piecewise Wiener deconvolution that makes it possible to get a good approximation to the ideal spatially varying deblurring algorithm while keeping a short execution time. Finally, in Sec. 4, we validate our theoretical models by comparing them with practical implementations of piecewise Wiener deconvolution and fully spatially varying deblurring applied to realistically simulated images.
2. Theoretical comparison of uniform and spatially varying deblurring
A hybrid imaging system consists of a lens with its sensor and a digital processing algorithm. In this paper, we consider two Cooke triplets optimized in two different ways. They are described in Sec. 2.1. Then, in Sec.2.2, we evaluate and compare the imaging performance obtained with hybrid systems composed of these two lenses combined with various deblurring algorithms. The imaging performance is evaluated using a theoretical model of the MSE between an ideal image of the scene and the final reconstructed image.
2.1 Presentation of co-designed triplets
The first considered optical system is a Cooke triplet optimized in a conventional way by minimizing the quadratic sum of the spot diagram diameters at selected positions in the FoV. Its maximum half-FoV is $20\deg$, with a focal length of $50$ mm and a $F/4$ aperture. It will be referred to as “conventional triplet” in the following. Its modulation transfer functions (MTFs) at different positions in the FoV are represented in Fig. 1(a), up to the Nyquist frequency of the sensor (taken here at $100$ lp.mm$^{-1}$, i.e. a $5$ $\mathrm {\mu }$m sensor pitch). They are calculated from the diffractive Point Spread Function (PSF) extracted from CodeV, which performs a physical calculation involving both diffraction and real aberrations, based on real ray-tracing and thorough Fourier Optics [15,16]. As most optical systems optimized in a standard way, it presents a good performance on axis but the performance degrades at peripheral positions in the FoV (Fig. 1(a)).
Fig. 1. (a) MTFs of the conventional triplet (b) MTFs of the MMSE triplet. The legend of (a) is the same as the legend of (b).
The second Cooke triplet we will consider has been optimized following the procedure described in [17]. In short, it has been co-optimized by minimizing the global averaged MSE between the deconvolved image and the ideal image of the scene denoted $MSE_G$), assuming that the deconvolution algorithm was a single linear filter common to all positions in the FoV. It will thus be referred to as “MMSE triplet” in the following. The MSE at the position $\psi$ in the FoV can be expressed, in the Fourier domain, in the following way [18]:
2.2 Influence of the deblurring algorithm on image quality
Let us now evaluate the imaging performance obtained by hybrid systems composed of these two triplets and 3 different post-processings, namely, no deblurring, deconvolution with the averaged Wiener filter (for which the MMSE triplet has been optimized) and deblurring with an ideal locally adaptive algorithm that takes into account the spatial variations of the MTF over the FoV.
When there is no processing, $\tilde {w}(\nu )$ can be considered as a unitary filter (filled with “$1$” until the Nyquist frequency, and “$0$” beyond) in Eq. (1). The expression of the local MSE then becomes:
Fig. 2. $MSE(\psi )$ (Eq. (1)) as a function of the field position parameterized by its angle when there is no deconvolution, there is a deconvolution of the whole FoV with a single average Wiener filter, and when there is an adapted deconvolution. (a) For the conventional triplet (b) For the MMSE triplet.
Let us now assume that the images are processed with the average Wiener filter (Eq. (3)), which is the filter for which the MMSE triplet has been co-optimized. In this case, the value of the local MSE is obtained by substituting the expression of the average Wiener filter (see Eq. (3)) for $\tilde {w}(\nu )$ in Eq. (1). We then obtain, for the two triplets, the purple curves on Fig. 2. We observe that in the case of the conventional triplet (Fig. 2(a)), deconvolution brings almost no improvement at the center of the FoV, and only a slight one at peripheral positions. Conversely, the improvement brought to the MMSE triplet is much larger since the MSE is divided by about three at all positions in the FoV (Fig. 2(b)). This strong improvement is due to the fact that, contrary to the conventional triplet, the MMSE triplet was co-optimized with the average Wiener filter. This filter is not suitable to work with an optical system whose OTFs varies greatly within the FoV - which is the case of the conventional triplet.
In order to improve the results of the conventional triplet after deconvolution, one must use a spatially varying deblurring algorithm that adapts to the local characteristics of the OTFs. For that purpose, let us assume an “ideal” algorithm that would implement, at each FoV position, the Wiener filter corresponding to the exact OTF at this position. The MSE that would be obtained with this purely theoretical method, that we will refer to as “adapted deconvolution”, can be deduced (from Eq. (1)) by substituting for $\tilde {w}(\nu )$ the expression of the standard Wiener filter adapted to the local MTF $\tilde {h}_\psi (\nu )$:
The values of the MSE defined in Eq. (5) are plotted as blue curves on Fig. 2. In the case of the conventional triplet, the performance is significantly improved with respect to the average Wiener filter (Fig. 2(a)). On the other hand, in the case of the MMSE triplet, the gain exists but is limited (Fig. 2(b)). This was expected since the benefit of using such an adaptive deconvolution method is greater when the optical properties have significant spatial variations. We may also notice that for both triplets, the adapted deconvolution (blue curve) does not lead to a constant MSE over the FoV: it presents globally the same variations as the MSE curve obtained without deconvolution (red curve), but attenuated. In fact, a signal too attenuated at some spatial frequencies (e.g. Figure 1(a) with poor raw MTFs which exhibit nullings), may not be recovered at all by deconvolution, and especially by Wiener filtering that avoids, by construction, to amplify the signal and also the noise at spatial frequencies where the raw $SNR(\nu )$ ($S_{oo}(\nu )/S_{nn}(\nu )$) is low. This means that even if the adapted deconvolution improves the overall image quality, the variability of the optical performance of the lens with the position in the FoV observed before deconvolution remains present after deconvolution. Thus, if a specific characteristic is desired for the hybrid optical/digital imaging system (for example, homogenization of performance in the field, or on the contrary, an excellent performance in the center specifically), it is necessary that the optical system itself already carries this characteristic, and thus be optimized in this sense.
3. Theoretical impact of piecewise Wiener deconvolution on image quality
We have considered in the previous section two “extreme cases” for digital post-processing: the average Wiener filter, which is the simplest and fastest one, and a hypothetical “perfectly adaptive” one, which takes into account the spatial variation of the MTF of the optical system but which implementation will inevitably be much more complex and computationally intensive. In this section, we propose a digital post-processing with a philosophy intermediate between these two extreme cases, which has the computational complexity of a linear shift-invariant filter but takes into account the spatial variations of the MTF. We investigate its theoretical behavior and compare its performance with the two “extreme case” algorithms. Implementation and performance evaluation on images will be discussed in Sec. 4.
3.1 Theoretical principle of “piecewise Wiener deconvolution”
The idea of piecewise Wiener deconvolution is to partition the entire FoV into square zones of equal areas, and to deconvolve each area with a specific average Wiener filter. We consider $N_Z$ zones $Z_1$, $Z_2$,…, $Z_{N_Z}$, described respectively by the sets of field positions $\Psi _{Z_1}$, $\Psi _{Z_2}$,…, $\Psi _{Z_{N_Z}}$. Each set contains an equal number K of field positions. We also define the set $\Psi _{\rm {FoV}} = \bigcup \Psi _{Z_n}$ corresponding to all FoV positions. An example of such a partition is given on Fig. 3 for a division into $N_Z=3\times 3=9$ square zones. Each area is deconvolved with an average Wiener filter (Eq. (3)) built from the $K$ OTFs characterizing each area. The obtained expression of the $MSE_G(\Psi )$ for a given region $\Psi$, defined in Eq. (2), can be written in the following way [23]:
withFig. 3. Schematic drawing illustrating the FoV positions describing the whole FoV and an example of division into $N_Z=9$ square zones.
3.2 Evolution of piecewise Wiener deconvolution performance with the number of zones
Figure 4 represents the value of the square roots of $MSE_I(\Psi _{Z_n})$, $MSE_D(\Psi _{Z_n})$, and $MSE_G(\Psi _{Z_n})$ for a division of the FoV into a variable number $N_Z$ of zones in the case of the conventional triplet. This number $N_Z$ of zones ranges from $1$ (which amounts to using a single deconvolution filter for the entire field) in the first row, to $25$ in the last row (the FoV is divided into $5 \times 5$ zones, and thus $25$ different average Wiener filters are used). When $N_Z=1$, the term $MSE_D$, that corresponds to disparity between MTFs (Fig. 4(1b)), is significantly larger that the term $MSE_I$, that corresponds to intrinsic deconvolution difficulty (Fig. 4(1a)). This was expected since we have seen in Fig. 1 that the MTFs of this triplet are quite disparate.
Fig. 4. For the conventional triplet used with piecewise Wiener deconvolution, square root of the average MSE of each area. (First column) Intrinsic difficulty of deconvolution $\sqrt {MSE_I}$. (Second column) Disparity between MTFs $\sqrt {MSE_D}$. (Third column) Global averaged MSE $\sqrt {MSE_G}$. The number of areas used for piecewise Wiener deconvolution depends on the row : (1) Deconvolution with a single average Wiener filter on the whole FoV; (2) four areas; (3) nine areas; (4) sixteen areas; (5) twenty-five areas. The colorbar is the same for all graphs.
As the number of zones increases, the two components take on values that are closer to each other and decrease on average. For this reason, the global averaged MSE (Eq. (6)) displayed in the third column of Fig. 4 decreases with the number of zones, which demonstrates the benefit of using spatially adaptive deconvolution. It should be noted, however, that this decrease is heterogeneous in the FoV. For example, the variation of the term $MSE_I$ is similar to the variation of the local MSE observed in Fig. 2(a)): it is larger at the edges of the FoV because the MTFs are lower there (Fig. 1(a)). Similarly, the term of disparity between MTFs, $MSE_D$, is also larger at the edges, especially even when the partition is done with a large number of zones, since the MTFs are also more disparate in the peripheral positions of the FoV (Fig. 4(5b)).
Figure 5 displays the same data as Fig. 4, but for the MMSE triplet. By comparing Fig. 4(1c) and Fig. 5(1c), it is seen that when deconvolving with a single filter, the MMSE triplet shows a better performance (i.e., a lower $MSE_G$) than the conventional triplet. This observation was already made in Fig. 2, but Fig. 5 allows us to detail the evolution of the two terms of $MSE_G$ with the number of zones. This evolution is not the same as for the conventional triplet: the disparity term $MSE_D$ is much lower because, by construction, the MTFs are close to each other (Fig. 1(b)), but the intrinsic deconvolution difficulty $MSE_I$ is a bit higher because the MTFs are lower. Moreover, by increasing the number of zones, $MSE_I$ (Fig. 5(a)) becomes slightly heterogeneous but less than for the conventional triplet since the MTFs are much less disparate over the FoV. On the other hand, $MSE_D$ (Fig. 5(b)) decreases until almost zero when the number of zones increases, since decreasing the size of the zones makes the MTFs within each zone even more similar to each other. Its value is also quite constant across the FoV. It follows that for a large number of zones (Fig. 5(5c)), the global averaged MSE for each zone is lower than with the conventional triplet and more homogeneous across the FoV.
Fig. 5. For the MMSE triplet used with piecewise Wiener deconvolution, square root of the average MSE of each area. (First column) Intrinsic difficulty of deconvolution $\sqrt {MSE_I}$. (Second column) Disparity between MTFs $\sqrt {MSE_D}$. (Third column) Global averaged MSE $\sqrt {MSE}$. The number of areas used for piecewise Wiener deconvolution depends on the row : (1) Deconvolution with a single average Wiener filter on the whole FoV; (2) four areas; (3) nine areas; (4) sixteen areas; (5) twenty-five areas. The colorbar is the same for all graphs.
3.3 Theoretical image quality with piecewise Wiener deconvolution
Let us now evaluate the global performance of the hybrid systems composed of the two triplets and piecewise Wiener deconvolution. The average MSE over the whole FoV is therefore:
where $MSE_G(\Psi _{Z_n})$ is calculated with Eq. (2), by substituting Eq. (3) for $\tilde {w}(\nu )$ and $\Psi _{Z_n}$ set for the $\Psi$ set into Eq. (1). Note that this $MSE_{FoV}(N_Z)$ expression is a theoretical limit as it is a somewhat optimistic simple theoretical approximation of the whole image MSE (in particular, it does not take into account edge effects that may occur, and a possible geometrical distortion is not considered). Instead of using directly $MSE_{FoV}$, we will quantify it, in an equivalent way, with the image quality (IQ): In this definition, we consider that the statistical variance of the image is normalized to one ($\int S_{oo}(\nu ) d\nu = 1$).This quantity is represented in Fig. 6 as a function of $N_Z$, for the two triplets. First, it can be seen, as in Fig. 2, that the performance before deconvolution (corresponding to $N_Z = 0$ on Fig. 6) is much better (about $+2.5$ dB) for the conventional triplet than for the MMSE triplet. After deconvolution with a single average Wiener filter (i.e., $N_Z=1$), the gain for the conventional triplet is quite low because the large variations of the MTFs over the FoV prevent deconvolution with a single filter from working properly. Then, when the number of deconvolution zones $N_Z$ increases, the $IQ$ of the conventional triplet increases slowly. The MMSE triplet has a very different behavior: most of the improvement is obtained with a single average Wiener filter ($N_Z=1$). Then, its performance slowly increases with the number of zones, but quickly reaches a horizontal asymptotic behavior. Thus, piecewise Wiener deconvolution allows to increase the average performance of the two systems, but in a different way according to the criterion used for their optimization.
Fig. 6. Theorical $IQ$ as a function of the number of areas for piecewise Wiener deconvolution, for the conventional triplet and the MMSE triplet. $N_Z = 0$ means that there is no deconvolution. The performance in case of adapted deconvolution is indicated by a dotted line.
For both triplets, piecewise Wiener deconvolution tends toward a maximal value, which is the performance of the adapted deconvolution. This performance can be evaluated by summing the local MSE defined in Eq. (5) over all the FoV positions. The corresponding $IQ$ values are represented by dotted lines on Fig. 6. We can see that they have similar values, even if the performance of the conventional triplet is slightly lower: $17$ dB instead of $17.3$ dB for the MMSE triplet. This similarity of adapted deconvolution performance for the two triplets could be interpreted as follows. These systems have globally the same “amount of aberrations” but these aberrations are simply distributed in a different way over the FoV: for the conventional triplet, they are low in the center of the FoV but high at the edges, and for the MMSE triplet, they are more homogeneously distributed. It follows from this difference in distribution that the “path” to achieve the limit performance by increasing the number of zones (thus reducing their size) is quite different for the two systems.
Therefore, to conclude, it is the MMSE triplet that benefits the most from piecewise Wiener deconvolution, since this algorithm is already effective with a small number of zones. This is mainly due to the much faster decrease of the $MSE_D$ disparity term with optical systems co-optimized for a single deconvolution filter (as can be seen in Fig. 5 and [14]). Of course, the performance of the conventional triplet is also gradually improved with such a method, but it tends much more slowly to the performance of the adapted deconvolution as the number of zones increases.
4. Practical implementation of algorithms taking into account the spatial variation of the PSF
We have so far evaluated the theoretical performance of two hybrid systems using piecewise Wiener deconvolution and adapted deconvolution. This study was based on simple image models, which has made it possible to derive theoretical, closed-form expressions of the MSE obtained with these algorithms. These expressions are interesting since they allow us to study the performance of any hybrid system and outline its characteristics with simple and fast computation. However, in practice, these deconvolution methods have to be implemented. The goal of this section is to propose algorithms that allow to reach the previously outlined theoretical performance and trends. These algorithms will be validated on realistically simulated images from the two considered triplets. First, we briefly present the direct model used for simulating images with spatially varying PSFs (Sec. 4.1). Then, we detail the implementation of the piecewise Wiener deconvolution and adapted reconstruction algorithms (Sec. 4.2). Finally, we evaluate their performance on simulated images (Sec. 4.3) and compare it with the theoretical performance analyzed in Sections 2. and 3.
4.1 Direct model of spatially varying PSF
To simulate images with PSFs continuously varying across the FoV, we use bilinearly interpolated PSFs with the fast formulation based on image weighting and convolutions proposed by Denis et al. and precisely described in [11]. The spatially varying blur operator is computed from an array representing a finite number of PSFs describing the optical system. This operator consists in the (well chosen) weighting of the ideal image around the positions corresponding to each PSF, and then in the convolution by each PSF of the different weighted image portions, and their appropriate juxtaposition. We adapt the code so that it can use the PSFs of the two studied triplets generated by the CodeV software (taking into account of the sensor pixel size). For this study, the image of Barbara (Fig. 7) - presenting many areas with periodic elements - will be used. Note that we have tested other images (having different PSDs) and obtained similar results, as can be seen in the Supplemental document.
This image is first degraded with a spatially varying blur applied from $31\times 31$ simulated PSFs (computed with CodeV). Once degraded by the simulated optical system, it is corrupted with noise according to a realistic model of modern CMOS sensors, combining Poisson noise and additive Gaussian noise, and leading to a SNR close to $36$ dB (the SNR depending slightly on the used triplet due to the presence of Poisson noise).
4.2 Practical implementation of the spatially varying deblurring algorithms
In this section, we describe the implementation of the two considered spatially varying reconstruction algorithms, namely, piecewise Wiener deconvolution and adapted deconvolution.
4.2.1 Piecewise Wiener deconvolution
Piecewise deconvolution is a method that comes quite naturally when one wishes to reconstruct a spatially varying PSF with several filters (both in depth of field [24,25] or in field of view [26]). To implement piecewise Wiener deconvolution in practice, we use the following method. It is illustrated in Fig. 8, for the case of a division of the FoV in $3 \times 3$ zones.
- • Step 2: $N_z$ “enlarged subparts” of the raw image (denoted $I_1$, $I_2$,…, $I_{N_Z}$) slightly larger than each zone (at least the size of the area to which is added a few pixels corresponding to the global size of the PSF, except on the edges of the whole image where the area is not enlarged) are selected and separated from each other. In the example, the enlarged subpart around zone $5$ is shown in dotted lines : it is $I_5$. In our case, we implement a border of size $10 \%$ the total size of the zone. Each enlarged subpart is deconvolved with the deconvolution filter whose size is adapted to the subpart : $\tilde {\hat {O}}_n(\nu ) = \tilde {I}_n(\nu ) \times \tilde {w}_{\Psi _{Z_n}}(\nu )$.
- • Step 3: A weighting $P_n(r)$ (with $r$ the spatial coordinates) is associated to each enlarged subpart. In order to build $P_n(r)$ we convolve a matrix $P_{I_n}$ with a Gaussian kernel g (of standard deviation $3 \%$ of the size of the zone), so that $\sum _{n=1}^{N_Z} P_n(r) = 1$ for any $r$ (for the boundaries, input array values outside the bounds of the array are assumed to be equal to the nearest array border value). Here, we use a matrix $P_{I_n}$ which is equal to $1$ all inside the zone, and $0$ everywhere else (edges of the enlarged zone). The example in Fig. 8 shows the weighting associated with the zone 5 obtained in this way.
- • Step 4: the total image is reconstructed by weighting each deconvolved enlarged subpart by the weights created in the previous step : $\hat {O}(r) = \sum _{n=1}^{N_Z} \hat {O}_n(r) \times P_n(r)$. This weighting varies gradually from the inside to the outside of each “guard zone” surrounding each zone (as described by hatching for zone n$^{\circ }$5 represented in step 4 in Fig. 8). A blending of the edges of the different independently deconvoluted areas is thus performed to avoid having too abrupt transitions in the overall image. This amounts to interpolating the edges of each area so that the “inter-area” edges are not visible. In the example, the reconstruction of the hatched part of zone 5 depends on the reconstruction of all the adjacent zones.
4.2.2 Spatially varying deblurring
In order to implement the “adapted deconvolution” as described in Sec. 2.2 (asymptotic behavior of the curves in Fig. 6), we use a spatially varying deblurring algorithm that takes into account continuous variations of the PSF across the FoV and is described in [11]. This article describes several interesting state-of-the-art algorithms and Refs. [27,28]. The algorithm we employ makes use of the direct imaging model (via a data fidelity term) and of a regularization term. We considered quadratic regularization (not total variation regularization as in [11]) in order to mimic the use of an average Wiener filter. The estimate of the image is then:
where $||.||$ denotes the 2-norm, y the observed image $I$ in vector form, x the image of interest $O$ in vector form, H the observation operator (implemented by the direct model), D the derivative operator and $\mu$ the regularization parameter. We choose by hand the regularization parameter $\mu$ that maximizes the final IQ. This convex optimization problem is solved using an iterative quasi-Newton method [29], and requires at least one hundred estimations of the direct model. Its computation time is thus much higher than piecewise Wiener deconvolution. The direct model is computed as described in the previous section, with a dense a priori of $31\times 31$ known PSFs in order to be as close as possible to the best achievable performance with this algorithm.4.3 Performance comparison of piecewise Wiener deconvolution and spatially varying deblurring
Let us now evaluate the performance of piecewise Wiener deconvolution and spatially varying deblurring on simulated images, and compare them with the theoretical results obtained in Sec. 2. and 3. The ideal scene is displayed in Fig. 7. Figure 9 shows the center subpart of images of this scene provided by hybrid systems based on conventional (first column) and MMSE (second column) triplets, when there is no deconvolution (first row), with piecewise Wiener deconvolution by $1\times 1$, $2\times 2$ and $5\times 5$ zones (second to fourth rows), and with spatially varying deblurring (fifth row). The third column is the subpart of the ideal image, for comparison. Figure 10 shows the same data for a subpart located at one edge of the image.
Fig. 9. Center subpart of the image of Barbara. (First column) Conventional triplet. (Second column) MMSE triplet. (a;f) Without deconvolution (b;g) Deconvolution with a single average Wiener filter (c;h) $2\times 2$ piecewise Wiener deconvolution (d;i) $5\times 5$ piecewise Wiener deconvolution (e;j) spatially varying deblurring.(k) Ideal image.
Fig. 10. Edge subpart of the image of Barbara. (First column) Conventional triplet. (Second column) MMSE triplet. (a;f) Without deconvolution (b;g) Deconvolution with a single average Wiener filter (c;h) $2\times 2$ piecewise Wiener deconvolution (d;i) $5\times 5$ piecewise Wiener deconvolution (e;j) spatially varying deblurring. (k) Ideal image.
Let us first consider piecewise Wiener deconvolution. We can notice that for the conventional triplet (Fig. 9, first column), the visual quality slowly improves as the number of zones increases. This is especially evident in the scarf pattern just below the woman’s chin. This observation confirms the theoretical predictions of Fig. 4. The improvement is less visible at the edges of the image (Fig. 10(a-d)). Indeed, at the periphery of the FoV, the spatial variation of the PSF shape is faster, resulting in a larger term of disparity between MTFs (Fig. 4(1b-3b)). To sufficiently reduce this disparity, it is thus necessary to reduce the size of the zones. As a consequence, we can observe that a minimum of 25 deconvolution zones is required for the reconstruction to appear visually acceptable (Fig. 10(d)). The MMSE triplet has a quite different behavior. First, the image without deconvolution is very blurry in the center (Fig. 9(f)), compared to the conventional triplet’s case (Fig. 9(a)). The whole image being blurred homogeneously, the edge subpart is also very blurry (Fig. 10), but one can observe that the details are not completely lost. As predicted theoretically (Fig. 5), the reconstruction with a single average Wiener filter is already very efficient, for the whole FoV (Fig. 9(g) and Fig. 10(g)). Then, the deconvolution by a higher number of zones improves the visual quality rather moderately, even if some details appear sharper and sharper (for example, the stripes on the bottom left of the image on the axis, Fig. 10(g-i)).
Let us now analyze the results of spatially varying deblurring. They are displayed in the fifth rows of Fig. 9 and Fig. 10. Visually, in the case of the MMSE triplet, the image quality is very close to that obtained with piecewise Wiener deconvolution with $5 \times 5$ zones, both on the axis (Fig. 9(i-j)) and on the edges (Fig. 10(i-j)). In sharp contrast, there is a wide performance gap between both algorithms with the conventional triplet, especially at the edge of the image (Fig. 10(d-e)). Nevertheless, at this position in the FoV, the image obtained with the conventional triplet and processed by spatially varying deblurring (Fig. 10(e)) has a worse quality than the one obtained with the MMSE triplet and piecewise Wiener deconvolution (Fig. 10(i)): a better processing cannot completely compensate for a significantly lower MTF at this position in the FoV. This effect was also clearly observed from the theoretical study of Sec. 2, by comparing the blue curves in Fig. 2(a) and Fig. 2(b). We can also observe another important point: whatever the post-processing algorithm used, the differences that initially exist between the conventional triplet and the MMSE triplet persist after processing. For example, the conventional triplet remains excellent in the center but less efficient on the edges. This can be seen both theoretically (compare red curve and blue curves in Fig. 2(a)) and in the image simulations (Fig. 9 and Fig. 10).
Let us now consider the global IQ values obtained with the different deblurring methods applied to simulated images from both triplets. They are represented in Fig. 11 as a function of the number of zones. It is observed that the graph is quite similar to Fig. 6, which was obtained from the theoretical model. For the conventional triplet, the IQ obtained with piecewise Wiener deconvolution (red solid line) increases progressively with the number of zones, which is consistent with our analyses of Fig. 9 and Fig. 10. For the MMSE triplet, the IQ curve (blue solid line) follows the expected theoretical behavior as well, with a sharp increase followed by an asymptotic behavior. Of course, the absolute IQ values are not the same as in Fig. 6 since the realistic image and the noise model used in these simulations are different from the theoretical generic PSD models used in Sec. 2.
Fig. 11. $IQ$ on the simulated image (Barbara) as a function of the number of areas for piecewise Wiener deconvolution, for the conventional triplet and the MMSE triplet. $N_Z = 0$ means that there is no deconvolution. The performance in case of spatially varying deblurring is indicated by a dotted line.
The IQ values obtained with spatially varying deblurring are plotted with dotted lines. They are close to each other, and the value obtained with the conventional triplet is slightly lower than that obtained with the MMSE triplet ($-0.5$ dB). This fits well the theoretical prediction (Fig. 6), except for the fact that the difference between the IQ obtained with spatially varying deblurring (dotted lines) and piecewise Wiener deconvolution (solid lines) is larger than theoretically expected. This can be explained by the fact that the average Wiener filters used in each region by piecewise Wiener deconvolution do not fully correspond to the PSD of the object and noise in simulated image, and that this algorithm presents some remaining edge effects at the borders of the zones. To reduce them, one could look at other edge conditions for processing: for example, using symmetric edge conditions rather than periodic edge conditions [30]. These defects prevent piecewise Wiener deconvolution from getting closer to the performance of spatially varying deblurring.
In order to verify the generality of our conclusions, we have conducted the same simulations and analyses on 3 totally different images (see Supplemental document). It is observed that the IQ curves are very similar, and that the analyses of the center and edge subpart of the simulated images lead to the same conclusions. This confirms that the theoretical MSE expressions developed in Sections 2. and 3. can be safely used to characterize the properties of a given lens under various post-processing algorithms.
In summary, we have described two algorithms that implement respectively piecewise Wiener deconvolution and spatially varying deblurring, and evaluated their performance for post-processing images acquired by the conventional and the MMSE triplets. We have shown that in the case of the MMSE triplet, piecewise Wiener deconvolution yields a performance slightly inferior but comparable to spatially varying deblurring, but for a much lower computational cost. Moreover, the simulation results confirm an important conclusion drawn from the theoretical studies in Sections 2. and 3.: to reach a given imaging goal with a co-designed hybrid system, it is preferable that the lens already have optical properties adapted to this goal. Digital processing cannot compensate for all the optical defects. In our case, the goal was to reach a good imaging performance over the whole FoV, and we have shown that for that purpose, the MMSE triplet, which was optimized with a true (MSE) co-design criterion and has nearly uniform MTFs without nulling, yields better average performance than the conventional triplet after deconvolution.
5. Conclusion
We have studied the contribution of spatially varying deblurring algorithms to improve the imaging performance of lenses affected by field aberrations. We first developed a simple model based on closed-form expressions in the Fourier domain of the MSE that enables us to evaluate, in a fast and efficient manner, the performance gain brought by such algorithms to any given lens defined by a set of OTFs over the FoV. This model has been validated on realistically simulated images post-processed by three different deblurring algorithms. In particular, we have proposed a piecewise Wiener deconvolution algorithm that takes efficiently into account the spatial variations of the PSF and has a computational complexity similar to a single deconvolution filter. Moreover, this algorithm does not require precise calibration of the spatially varying PSF thanks to the “averaging” nature of the filters over each zone.
The performance of an imaging system relies on the quality of both the lens and the digital post-processing [31]. One of the main conclusions of this article is that the contribution of these two components must be balanced. The processing alone cannot completely compensate for the variations of the optical quality over the FoV. In our case, the best solution is provided by a lens which has been optically designed to reduce field aberrations at the price of lower average optical quality, and a fast and “slightly” spatially varying piecewise Wiener deconvolution algorithm that efficiently enhances the global image quality and compensates for the residual non-homogeneity over the FoV.
The two lenses considered in this article were optimized with the CodeV lens design software by assuming no processing (conventional triplet) or deconvolution with an average Wiener filter (MMSE triplet). The spatially varying deblurring algorithms were used only to improve the performance of these already optimized lenses. The most interesting perspective to this work would be to perform lens optimization by taking explicitly into account, right from the optimization process in the lens design software, the fact that the deblurring algorithm is spatially varying. The closed-form MSE models developed in Sec. 2. and 3. provide a practical way to do so. Analyzing the influence of the SNR value on lens optimization and deblurring performance (in particular, higher noise levels corresponding to more challenging imaging conditions) is another important perspective for future work.
Acknowledgments
The work reported in this study was supported in part by the Agence de l’Innovation de Défense (AID) that provides half of a PhD fellowship to Alice Fontbonne. We also deeply thank Loïc Denis for providing us the codes of [11], and for his useful remarks on this article.
Disclosures
The authors declare no conflicts of interest.
Data availability
The exhaustive description of the two lenses analyzed in this article (Conventional and MMSE triplet) is available in the Supplemental document of [17]. The code underlying the results presented in this article is not publicly available at this time but may be obtained from the authors upon reasonable request.
Supplemental document
See Supplement 1 for supporting content.
References
1. E. R. Dowski and W. T. Cathey, “Extended depth of field through wave-front coding,” Appl. Opt. 34(11), 1859–1866 (1995). [CrossRef]
2. M.-A. Burcklen, H. Sauer, F. Diaz, and F. Goudail, “Joint digital-optical design of complex lenses using a surrogate image quality criterion adapted to commercial optical design software,” Appl. Opt. 57(30), 9005 (2018). [CrossRef]
3. D. G. Stork and M. D. Robinson, “Theoretical foundations for joint digital-optical analysis of electro-optical imaging systems,” Appl. Opt. 47(10), B64–B75 (2008). [CrossRef]
4. T. Vettenburg and A. R. Harvey, “Holistic optical-digital hybrid-imaging design:wide-field reflective imaging,” Appl. Opt. 52(17), 3931–3936 (2013). [CrossRef]
5. J. Wang, L. Wang, Y. Yang, R. Gong, X. Shao, C. Liang, and J. Xu, “An integral design strategy combining optical system and image processing to obtain high resolution images,” in Remotely Sensed Data Compression, Communications, and Processing XII, B. Huang, C.-I. Chang, and C. Lee, eds. (SPIE, 2016).
6. J. Portilla and S. Barbero, “Hybrid digital-optical imaging design for reducing surface asphericity cost while keeping high performance,” in Computational Optics II, D. G. Smith, F. Wyrowski, and A. Erdmann, eds. (SPIE, 2018).
7. O. Lévêque, C. Kulcsár, and F. Goudail, “Comparison of linear and nonlinear deconvolution algorithms for co-optimization of depth-of-field enhancing binary phase masks,” OSA Continuum 4(2), 589 (2021). [CrossRef]
8. A. W. Lohmann and D. P. Paris, “Space-variant image formation*,” J. Opt. Soc. Am. 55(8), 1007 (1965). [CrossRef]
9. J. G. Nagy and D. P. O’Leary, “Restoring images degraded by spatially variant blur,” SIAM J. Sci. Comput. 19(4), 1063–1082 (1998). [CrossRef]
10. M. Hirsch, S. Sra, B. Scholkopf, and S. Harmeling, “Efficient filter flow for space-variant multiframe blind deconvolution,” in 2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, (IEEE, 2010).
11. L. Denis, E. Thiébaut, F. Soulez, J.-M. Becker, and R. Mourya, “Fast approximations of shift-variant blur,” Int. J. Comput. Vis. 115(3), 253–278 (2015). [CrossRef]
12. F. Sroubek, J. Kamenicky, and Y. M. Lu, “Decomposition of space-variant blur in image deconvolution,” IEEE Signal Process. Lett. 23(3), 346–350 (2016). [CrossRef]
13. J. Zhang, J. Pan, J. Ren, Y. Song, L. Bao, R. W. Lau, and M.-H. Yang, “Dynamic scene deblurring using spatially variant recurrent neural networks,” in 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, (IEEE, 2018).
14. A. Fontbonne, H. Sauer, and F. Goudail, “How to integrate digital post-processing algorithms in professional optical design software for co-designing complex optical systems?” in Optical Design and Engineering VIII, J. Babington, U. Fuchs, and L. Mazuray, eds. (SPIE, 2021).
15. J. Goodman, Introduction to Fourier Optics, 4th Ed. (W. H. Freeman, Macmillan Learning, 2017).
16. H. Hopkins and M. J. Yzuel, “The computation of diffraction patterns in the presence of aberrations,” Optica Acta: Int. J. Opt. 17(3), 157–182 (1970). [CrossRef]
17. A. Fontbonne, H. Sauer, and F. Goudail, “Comparison of methods for end-to-end co-optimization of optical systems and image processing with commercial lens design software,” Opt. Express 30(8), 13556–13571 (2022). [CrossRef]
18. D. Robinson and D. G. Stork, “Joint design of lens systems and digital image processing,” in International Optical Design, (Optical Society of America, 2006), p. WB4.
19. D. L. Ruderman, “Origins of scaling in natural images,” Vision Res. 37(23), 3385–3398 (1997). [CrossRef]
20. A. van der Schaaf and J. van Hateren, “Modelling the power spectra of natural images: Statistics and information,” Vision Res. 36(17), 2759–2770 (1996). [CrossRef]
21. R. N. Clark, “Digital camera reviews and sensor performance summary,” http://www.clarkvision.com/articles/digital.sensor.performance.summary (2016).
22. F. Diaz, F. Goudail, B. Loiseaux, and J.-P. Huignard, “Increase in depth of field taking into account deconvolution by optimization of pupil mask,” Opt. Lett. 34(19), 2970–2972 (2009). [CrossRef]
23. D. Robinson and D. G. Stork, “Extending depth-of-field: Spherical coding versus asymmetric wavefront coding,” in Frontiers in Optics 2009/Laser Science XXV/Fall 2009 OSA Optics & Photonics Technical Digest, (Optical Society of America, 2009), p. CThB3.
24. E. Maalouf, B. Colicchio, and A. Dieterlen, “Fast deconvolution with non-invariant PSF for 3-D fluorescence microscopy,” in Optical and Digital Image Processing, vol. 7000P. Schelkens, T. Ebrahimi, G. Cristóbal, and F. Truchetet, eds., International Society for Optics and Photonics (SPIE, 2008), pp. 478–485.
25. E. Maalouf, B. Colicchio, and A. Dieterlen, “Fluorescence microscopy three-dimensional depth variant point spread function interpolation using zernike moments,” J. Opt. Soc. Am. A 28(9), 1864–1870 (2011). [CrossRef]
26. T. P. Costello and W. B. Mikhael, “Efficient restoration of space-variant blurs from physical optics by sectioning with modified wiener filtering,” Digital Signal Process. 13(1), 1–22 (2003). [CrossRef]
27. P. Perona, “Deformable kernels for early vision,” IEEE Trans. Pattern Anal. Machine Intell. 17(5), 488–499 (1995). [CrossRef]
28. D. Miraut and J. Portilla, “Efficient shift-variant image restoration using deformable filtering (part i),” EURASIP J. Adv. Signal Process. 2012(1), 100 (2012). [CrossRef]
29. J. Nocedal, “Updating quasi-newton matrices with limited storage,” Math. Comp. 35(151), 773–782 (1980). [CrossRef]
30. M. K. Ng, R. H. Chan, and W.-C. Tang, “A fast algorithm for deblurring models with neumann boundary conditions,” SIAM J. Sci. Comput. 21(3), 851–866 (1999). [CrossRef]
31. E. Tseng, A. Mosleh, F. Mannan, K. St-Arnaud, A. Sharma, Y. Peng, A. Braun, D. Nowrouzezahrai, J.-F. Lalonde, and F. Heide, “Differentiable compound optics and processing pipeline optimization for end-to-end camera design,” ACM Trans. Graph. 40(2), 1–19 (2021). [CrossRef]
