Optical aberrations Correction in Postprocessing using Imaging Simulation

For photographic objective systems with large FOVs and apertures, it is inaccurate to derive the PSF analytic solution by Gaussian paraxial approximation (Born and Wolf, 2013) or fast Fourier transform (FFT) (Manual, 2009) (we recommend readers to supplementary materials for relevant information). In other words, for various commercial optical designs, sequential raytracing and coherent superposition are the only universal methods to obtain precise point spread functions on the imaging plane.

In the image formation process, every point of the object can be regarded as a coherent light source $\mathcal{S}$. For a particular light source $\mathcal{S}$, a grid of rays is launched through the optical system, which forms the wavefront $\mathcal{\sum}$ at the exit pupil. According to the Huygens principle, every point on a wavefront is the source of spherical wavelets, and the secondary wavelets emanating from different points mutually interfere (Baker and Copson, 2003; Paul, 2014). In other words, the diffraction intensity at any point on the image plane is the complex sum of all these wavelets, squared. Therefore, in our proposed optical PSF model, the PSF on the image grid is computed in two stages, the first stage uses the sequential raytracing method to compute the wavefront distribution at the exit pupil plane, and the second stage uses coherent superposition to output the PSF calculation.

Each pixel in a conventional camera sensor is covered by a single red, green, or blue color filter (Brooks et al., 2019). These color filters have different responses to light waves of different wavelengths, which can be formulated as $\mathit{C_{wav}(\lambda)}$. Moreover, in the imaging formation process, the off-axis field will suffer a loss of illumination, which is caused by lens shading and can be formulated as $\mathit{C_{ill}(fov,\lambda)}$. Therefore, for a precise simulation, the PSFs of different wavelengths and FOVs are multiplied by the specific distribution of wavelength responses and lens shading parameters. The influence of the wave response and relative illumination can be easily modeled by a coefficient:

Here $fov$ is the normalized FOV, and $\lambda$ is the wavelength of object point $\mathcal{S}$.

Owing to the differences between the dynamic range of raw sensor data and the sensitivity of the humans eye, imaging signal processing systems are typically used to convert the raw sensor energy data to photographs (Hasinoff et al., 2016). In our imaging simulation pipeline, we transform the natural object image $L$ to the energy domain data $L_{e}$ with the energy domain transformation, which is shown in Figure 5. Specifically, for a general natural image, we first perform gamma decompression with the standard gamma curve, while fixing the input to the gamma curve with $\epsilon=10^{-8}$ to prevent numerical instability. We refer the readers to the work by (Brooks et al., 2019) for a more in-depth discussion of gamma decompression. After gamma decomposition, we convert the sRGB image to camera RGB data with the inverse of the fixed $3\times 3$ color correction matrix (CCM). The CCM is determined by the camera sensor. Then, the white balance (WB) of the image is inverted by dividing the white balance value at random color temperature. The WB values of different color temperatures are obtained by calibration. In this way, we obtain the energy domain image $L_{e}$.

Based on the optical PSFs $I_{x^{\prime\prime}y^{\prime\prime}}$ calculated in Section. 4.1, we degrade the energy domain image $L_{e}$ with partitioned convolution. Partitioned convolution can be divided into three parts: Firstly we crop the image into uniform small pieces, then these pieces are convoluted with the PSFs of the corresponding FOVs, and finally the degraded pieces are spliced together. It is worth mentioning that the optical PSF calculated in Algorithm 1 is the PSF of one point $(x_{0},y_{0},z_{0})$. To precisely simulate the imaging results, we uniformly sample the object plane and perform the PSF calculation in Algorithm 1 on every sampled point. The overall PSF can be represented by $I(h,w,\lambda)$. Here, $h$, $w$ are the coordinates on the sensor plane (corresponding to the sampled point on the object plane), and $\lambda$ is the wavelength of the light source. For a captured image with $n\times m$ resolution, let $J_{e}$, $L_{e}$, $N_{e}$ $\in\mathbb{R}^{n\times m}$ be the observed energy image, the underlying sharp image, and the measurement noise image, respectively. The formation of the blurred observation $J_{e}$ can then be formulated as:

Here, $fov$ can be calculated by $(h,w)$ and the full FOV of the lens. It is important to note that $J_{e},L_{e}$ are raw-like images that indicate the energy received by the sensor.

After generating the blurred observation $J_{e}$ with Eq.(11), we invert the energy domain transformation by sequentially applying the operations in the lower-right corner of Figure 5. To be more specific, because the red, green, and blue color filters are arranged in a Bayer pattern, such as R-G-G-B, we omit two of its three color values according to the Bayer filter pattern of the camera. After acquiring the mosaiced raw energy data, we approximate the shot and read noise together as a single heteroscedastic Gaussian distribution and add it to each channel of the raw image (Brooks et al., 2019). We refer the readers to (Brooks et al., 2019) for details on synthetic noise implementation. Moreover, we sequentially apply the variation in the AHD demosaic algorithm (Lian et al., 2007), white balance, color correction matrix and the inverse of the gamma operator to $J_{e}$. The pipeline of the imaging simulation system is shown in Figure 5.

In conclusion, there are two significant differences between our proposed method and the imaging simulation methods that have been used in $Zemax^{\circledR}$ or $CODEV^{\circledR}$. First, our method calculates the superposition of every spherical wavelet, where the phase term of the ray is considered. Therefore, the proposed method can generate more accurate PSFs when compared with the traditional raytracing or FFT methods. Second, we combine the sensor information and the ISP parameters into the imaging simulation system. In this way, our method not only simulates the imaging process from the physical process but also obtains very similar imaging results visually.

We propose a novel network for the retrieval of the latent image $L$ from corrupted data $J$. Because we generate the training data by simulation, the input and the ground truth of our network are pixel-to-pixel aligned. This characteristic greatly improves the stability of training. In addition, the aberrations have a high correlation with the FOV, and we compute the $x$ and $y$ coordinates of each pixel and concatenate the FOV image with the aberration inputs (as shown in Figure 6). This operation enables our network to directly perceive the FOV information, which can be regarded as another kind of natural prior. The proposed network is shown in Figure 6. Specifically, we adopt a variant of the UNet architecture (Ronneberger et al., 2015).

In regard to the loss function, the strictly paired data for supervised training allow us to use fidelity loss to learn robust texture details, thereby preventing fake information in the recovered image $L$. Specifically, the fidelity loss function is:

where $\theta$ denotes the learned parameters in the network. $I_{n}^{ab}$ are the aberration inputs, and $I_{n}^{gt}$ are the corresponding ground truths of $I_{n}^{ab}$. We note that many loss functions (such as perceptual loss (Johnson et al., 2016) or GAN loss (Ledig et al., 2017)) will generate sharper visual results on natural images. These loss functions also introduce fake texture to the restorations and are not good at suppressing noise. Because the noise level of mobilephone varies greatly while illumination and color temperature are changing. Therefore, we adopt the MSELoss, which is robust to different noise levels.

When the field of view increases, the higher-order aberration coefficient deviates from that computed by the Gaussian paraxial model. The different sizes and shapes of off-axis PSFs are caused by these deviations. However, the traditional convolution operation strictly adopts the feature of fixed locations around its center, which will introduce irrelevant features into the output and discard the information truly related to PSFs.

The size of optical PSFs ranges from 100 pixels to 1,600 pixels (squared), which means that our network must extract the features based on the FOV information. Therefore, we introduce the deformable convolution in the decorder path (Dai et al., 2017; Zhu et al., 2019). While different from the method in (Chang et al., 2020), we directly concatenate the input and the upsampled offset. After the concatenation, we refine the output offset with a convolution layer. In this way, the shape and the size of the convolution kernel can be changed adaptively according to the field-of-view information.

Multiscale information is important for handling large PSFs. However, multiple downsamplings will increase the overhead of the network and destroy image structure. Therefore, we introduce a context block (Chang et al., 2020) into the minimum scale between encoder and decoder, which increases the receptive field and reconstructs multiscale information without further downsampling. As shown in Figure 6, several dilated convolutions are used to extract features from different receptive fields, which are concatenated later to estimate the output.

As shown in Figure 7, the simulation result of our simulation system and the corresponding digital image are the input and ground truth of our network, respectively. In the partitioned convolution operation in Figure 5, we calculate 240,000 ($400\times 600$) PSFs for the customized DSLR camera lens and 120,000 ($300\times 400$) PSFs for HUAWEI HONOR 20. The number of calculated PSFs is decided by the resolution of the sensor, i.e., the resolution of the sensor in the DSLR camera is $4,000\times 6,000$, and we compute PSFs for every $10\times 10$ patch. We adopt DIV2K (Martin et al., 2001) which contains 800 images of 2K resolution as the digital images to generate training data pairs.

In regard to the PSF calculation, we set the fixed object distance $z_{0}$ for each camera. The shooting distance is $1750mm$ for the DSLR camera, and it is $600mm$ for HUAWEI HONOR 20. For the DSLR camera, we compute 31 PSFs with wavelengths varying from 400nm to 700nm with a 10nm interval. For HUAWEI HONOR 20, we calculate 34 PSFs varying from 400nm to 730nm. The wave distribution of PSFs is determined by the upper limit and the bottom limit of the sensor’s wave response. Finally, the PSFs $I(h,w,\lambda)$ are reunited by the wave response functions of the sensor and lens shading parameters $\mathcal{C}_{e}(fov,\lambda)$.

To evaluate the effect of our trained network, we photograph 100 images with a DSLR camera and HUAWEI HONOR 20, respectively. After acquiring the raw data taken by the camera, we sequentially apply the steps in Figure 3.b to obtain the sRGB image. These sRGB photographs make up the testing datasets. We note that the sharpen operation in classical ISP will affect optical aberration degradation, therefore, we use the self-defined postprocessing pipeline to transform raw data into sRGB images.

In this way, we construct the data pairs for our network. In our network architecture, the channel numbers of each layer are marked at the bottom of the corresponding blocks (as shown in Figure 6). In terms of the hyperparameter of training, our settings are basically the same as (Chang et al., 2020). In addition, the $H,W$ of our training input is $256\times 256$, which allows the network to receive a large amount of spatially variant degradation.

The optical PSF distributions of the customized DSLR camera lens and HUAWEI HONOR 20 are sketched on the top of Figure 8, which are centrally symmetric and have a high correlation with the FOV. We note that the optical PSFs of HUAWEI HONOR 20 reveal a wide range of spatial variation over the region, while the PSFs of the customized DSLR camera are relatively consistent. Our network architecture takes advantage of these characteristics. It engages with the FOV information to perform spatial adaptive reconstruction for degraded images. It is worth mentioning that our method directly calculates the resulting PSFs, which is different from the optimization method of PSF estimation (Mosleh et al., 2015).

We plot the Strehl ratio (Mahajan, 1982) curves and the measured patches in $0.1$, $0.3$, $0.7$, and $0.96$ FOVs of two optical designs, see the bottom of Figure 8. The Strehl ratio can be used to evaluate the energy diffusion in different FOVs, which has some correlation to the degradation of the image. As shown in Figure 8, the Strehl ratio turns lower from the center to the whole field, resulting in more blurred photography. Moreover, unlike the customized DSLR camera lens, the shape and size of PSFs in HUAWEI HONOR 20 abruptly change in regard to the $0.9$ FOV, which is common to camera lenses. Thus, the PSF characteristics of different devices are totally different, which means that it is significant to generate targeted datasets for a specific lens.

To prove the accuracy of the proposed imaging simulation system, we designed a two-pronged experiment. One prong is to verify the accuracy of the calculated PSFs: First, we calibrate the PSFs of the DSLR camera with the virtual point-like source (Jemec et al., 2017). Then, the PSFs of the same camera are estimated by the optimization method used in (Mosleh et al., 2015) (we use the chart of the noise pattern). Finally, we calculate the PSFs by $Zemax^{\circledR}$’s raytracing. The comparison is shown in Figure 9, and we postprocess all the PSFs by the same postprocessing operations for visualization. When the FOV increases, the PSFs estimated by (Mosleh et al., 2015) (Figure 9.(a)) are becoming more blurry and the $PSNRs$ decrease greatly. We note that from the center to the edge of the image, the relative illumination decreases a lot. Therefore, the edge’s results of (Mosleh et al., 2015) are not ideal because the method has difficulty converging when the degraded information is masked by noise. The marginal PSFs calibrated by (Jemec et al., 2017) (Figure 9.(b)) also suffer from this problem. Different form these methods that need to shoot targets, our method directly calculate the PSFs from the lens prescriptions and will not be affected by noise. By comparing the $PSNRs$ of the proposed method and $Zemax^{\circledR}$, we prove that our results are more accurate and more similiar to the real image. Moreover, these results also validate the necessity of the coherent superposition.

Another method for verifying the authenticity of our method is to compare the simulation results. We shoot the checkerboard chart with the corresponding camera and simulate the imaging results with $Zemax^{\circledR}$, where the parameters used are the same as those in our simulation. Moreover, we also perform an ablation study on the proposed imaging simulation system. The visual comparison of the results is shown in Figure 10. Obviously, the imaging results of our method is the closest to the real data. We note that the image generated by $Zemax^{\circledR}$ has no combination of the sensor’s spectral information and the ISP parameters, accurately. This might be the reason that $\textit{Zemax}^{\circledR}$ produces inaccurate results. In addition, we compare the modulation transfer function (MTF) (Schroeder, 1981) in different FOVs, which shows the accuracy of the proposed method through the whole region of the sensor (as shown in Table 1).

In summary, we prove the authenticity of the proposed imaging simulation system mainly from two aspects. First, we demonstrate that our PSF calculation method is more accurate than the optimization method. When compared with the calibration methods, the proposed method is more practical and will not be affected by noise. The comparison between $Zemax^{\circledR}$ and the proposed method proves the accuracy of our method and validates the necessity of coherent superposition. Second, we compare the simulation results with the image of $Zemax^{\circledR}$’s simulation and the photograph of the camera, which verifies the authenticity of the imaging simulation system and every component is necessary.

We perform the ablation study of our network architecture on the simulation dataset of the DSLR camera lens. The results are shown in Table 2, and the evaluation metric is PSNR.

Ablation on FOV Input We concatenate the FOV image with the aberrations input image. It is obvious that the FOV image helps the network determine the spatially variant blurring and displacement of optical aberrations.

Ablation on Deformable ResBlock Without the deformable ResBlock, the network loses the spatial adaptability to handle different PSFs of optical aberrations. If we replace the deformable ResBlock with ResBlock, the performance of our networks significantly decreases, which indicates its effect.

Ablation on Context Block The context block helps the network capture larger field information while not damaging the structure of the feature map. The comparison validates its effectiveness.

Because of the spatial variation over aberrant PSFs, we design a brand new architecture to perform spatially adaptive recovery on the degraded image. In this section, we compare our algorithm with state-of-the-art deblurring methods.

We compare the chromatism and the MTF of degraded image and recovered image (shown in Figure 14). Though only trained with synthetic data, the proposed method can restore perfect knife edge and greatly improve the MTF results. Moreover, the pixel values of the R, G, and B channels of the reconstructed image are closer to each other, which indicates that our method successfully handles the chromatism between different channels. The lower restored CA curve also proves that the shift between RGB channels is smaller. Moreover, the enhancement of line-pair per inch (LPI) is shown in Figure 12, where the line-pair resolution is improved in the center of the ISO12233 chart. We note that the restored MTF results approach the diffraction limit of the lens, which is the physical limit resolution of an imaging device.

The ultimate goal of the proposed work is to integrate the module of aberration correction into an ISP system and to produce nearly perfect images. In today’s handcrafted ISPs, sharpening is the core method for alleviating the blurring caused by optical aberrations. Unfortunately, the sharpening operation is globally consistent and cannot solve the spatially variant blur. Therefore, it introduces unsolved blur in the boundary of FOV and ”ringing” effect in the center of figure (as shown in Figure 2c). In Figure 13, we compare the results of our proposed pipeline (as shown in Figure 3) and the images generated by HUAWEI ISP (for more results, please refer to the supplementary materials). Obviously, our proposed method adapts to the spatially variant degradation and the results are visually superior to HUAWEI’s ISP. Together with accurate synthetic data and the customized restoration algorithm, we demonstrate that it is possible to handle optical aberrations in postprocessing. This also encourages us to use a lightweight lens for high-quality imaging.