Point Spread Function Filtering for Radially Variant Blur Restoration

In this paper, we present a point spread function (PSF) filtering technique for solving the radially variant blur restoration problem. Radially variant blur is generated by a spherical single-element lens imaging system (SSLIS) that is embedded in an experimental camera module. The restoration of this category of blur is carried out in a polar coordinate system using polar PSFs at different fields of view (FOVs). However, restoration using large PSFs tends to introduce severe ringing artifacts in the restored image owing to the nonsparse nature of these PSFs. We show in this paper that the PSF filtering technique can effectively minimize ringing artifacts by filtering out some PSF pixels with an intensity lower than the threshold intensity. As a result, a nonsparse PSF becomes a sparse PSF, which is for good restoration results. The effectiveness of the PSF filtering technique was validated by visual comparison using three test images captured by the SSLIS camera module. In addition, a systematic way to determine the optimal filtering coefficient for a PSF at any FOV within the FOV range is also introduced.


Introduction
The point spread function (PSF) plays a very important role in the formation and restoration of an image. In an image formation model, the final recorded scene on the image sensor is a function of the PSF. (1)(2)(3)(4)(5)(6) The final image quality greatly depends on the size, shape, and intensity of PSFs distributed across the image area. In the literature of image restoration, a PSF serves as a restoration kernel that can be used to restore the recorded blurred image to a deblurred one similar to the original scene. (1)(2)(3)(4)(5)(6) The image restoration process involves two-dimensional deconvolution between the PSF kernel and the recorded image, which is usually discretized and carried out in matrix form. If the PSF is spatially invariant, the matrix of the PSF takes the form of a block Toeplitz matrix with Toeplitz blocks (BTTB). (7)(8)(9)(10)(11) It has been proved in our previous paper that a sparse PSF BTTB results in good restored image quality. (12) The sparsity of a BTTB can be defined by the ratio of the bandwidth of the BTTB to the number of Toeplitz blocks in the BTTB, denoted β/n. The bandwidth β is determined by the PSF size, while n is determined by the recorded image size. If the PSF size is less than the image size, the ratio β/n will be small, which indicates a sparse PSF BTTB. However, in many practical restoration problems, we need to deal with a large nonsparse matrix. Restoration using a nonsparse BTTB directly introduces strong boundary ringing artifacts across the restored image. Researchers in the field of image restoration have proposed iterative algorithms for reducing the ringing artifacts and gradually improve a restored image until it is similar to the original scene without changing the BTTB sparsity. For example, Reeves proposed an iterative method to gradually approach the true solution of the padded elements of a lexicographically ordered column vector of the degraded image (13) (referred to as "the captured image" in our paper) such that the ringing artifacts are markedly minimized. Although the iterative method is very effective and efficient for image restoration using a sparse PSF, it is time-consuming when the PSF matrix is not sparse.
We proposed a spherical single-element lens imaging system (SSLIS) in our previous research, (12,(14)(15)(16)(17) aiming to reduce the optical complexity and manufacturing costs of a digital camera. This imaging system introduces a radially expanding blur from the center to the four vertices of the image plane because optical aberrations such as field curvature become severe as the field of view (FOV) increases. Since the imaging system has only a single-element lens, conventional optical means of reducing the blur using multiple lenses are not possible. Therefore, we proposed an image processing method of reducing the blur that uses a polar domain deconvolution algorithm. (12,14) This algorithm converts the recorded image and PSF to a polar coordinate system and carries out 2D deconvolution in the polar domain. The PSF matrix used in the algorithm takes the form of a block Toeplitz matrix with circulant blocks (BTCB) rather than a BTTB because pixels on the left and right boundaries of the converted panoramic polar image are neighbors (see ref. 12 for details). The aforementioned β/n ratio is also a key factor that affects the restored image quality. Similarly to most of the spatially invariant blur restoration problems, a nonsparse matrix with a large β/n ratio will introduce severe ringing artifacts across the radially restored image, which may reduce the image quality.
In this paper, we propose a PSF filtering technique for solving the radially variant blur restoration problem. This technique increases the sparsity of a PSF BTCB by filtering out PSF pixels whose intensities are lower than the threshold intensity so that some low intensity entries in the BTCB become zeros. The restored image is therefore visually superior owing to the decrease in the β/n ratio. This technique does not markedly modify the original PSF shape and keeps the PSF pixel values that are higher than the threshold intensity unchanged. Only the size of the original PSF is reduced. This method does not involve an iterative algorithm, meaning that it is not computationally demanding. We believe that the proposed technique can also be applied to other image restoration problems. The effectiveness of this technique will be demonstrated in § 3.

Spherical single-element lens imaging system (SSLIS) and radially expanding blur
The SSLIS camera module is shown in Fig. 1(a). It consists of a spherical singleelement lens fixed inside a lens holder and a CMOS image sensor IC chip. The lens has a double-convex shape and a fixed focal length of 10.0 mm. Its refractive index is 1.673 at the design wavelength of 587.6 nm. The lens holder is screwed into the CMOS image sensor IC chip and the distance between the lens and the image sensor can be adjusted. The image sensor has 2048 [H] × 1536 [V] pixels and the image area is 6.6 × 4.9 mm 2 . The maximum FOV can be calculated using the equation y diag = ƒ tan(θ max /2), where ƒ denotes the focal length and y diag is the semi-diagonal distance. The maximum FOV θ max is 44.6° for our SSLIS.
The PSF distribution of the SSLIS across an image is shown in Fig. 1(b). All the PSFs shown in the figure are drawn at the same scale, and it can be observed that the size of the PSF increases from low FOVs to high FOVs, indicating that the blur at high FOVs is severer than that at low FOVs. As the imaging system is rotationally symmetric, some of the FOV corresponds to positions with the same distance from the image center. The PSFs at these positions can be regarded as the same PSFs in the polar coordinate system. All the PSFs distributed across the image are radially expanding but rotationally invariant.

Image restoration model and structure of PSF BTCB
The polar domain deconvolution for radially variant blur restoration can be expressed as the following matrix-vector multiplication: where (H * H + αL * L) −1 H * is a constrained least-squares filter. The recorded scene on the image sensor and the restored image in the polar coordinate system are represented by the lexicographically ordered column vectors g and f , respectively. H and H * denote the PSF BTCB and its conjugate transpose, respectively. The term αL * L in eq. (1) represents Tikhonov regularization, where α is the regularization parameter and L is the regularization operator.
If the size of the polar image is m × n, then the column vectors g and f will have nm × 1 entries and the BTCB H as well as L will have nm × nm entries. The BTCB structure and its relation to the polar image and polar PSF are shown in Fig. 2. As shown in Figs. 2(a) and 2(b), the BTCB consists of n circulant blocks. Those with nonzero entries are denoted by H i, j (i = 1...n, j = 1...n) and those with all zero entries are represented by a space in Fig. 2. The row and column bandwidths are marked in Fig. 2(a), which are both β. Essentially, the bandwidths are determined by the polar PSF size along the vertical direction, which has 2β + 1 rows, as shown in Fig. 2

(d).
Here, The entries of each circulant block H i, j are shown in Fig. 2(b). Note that each row or column has m entries, in which the number of nonzero entries is 2α + 1. α is determined by the polar PSF size along the horizontal direction, which can be confirmed from Fig.  2(d). The panoramic polar image converted and stretched from the recorded image is given in Fig. 2(c), which has m columns and n rows. The line in the recorded image represents the location where the left and right boundaries of the panoramic polar image are connected, which is the reason why the PSF matrix has circulant rather than Toeplitz blocks. It is then easy to find that the BTCB matrix size and the size of each circulant block are determined by the polar image size. If the polar PSF is very small compared with the polar image (i.e., 2β + 1 rows in the polar PSF is very small compared with n rows in the polar image, and 2α + 1 columns in the polar PSF is also very small compared with m columns in the polar image), the circulant block and BTCB will become sparse matrices. The trick in our PSF filtering technique is to reduce the polar PSF size without destroying its original shape, while maintaining the polar image size so as to increase the sparsity of the BTCB. Details will be given in § 2.3.

PSF filtering technique
The fundamental principle of the PSF filtering technique can be summarized by the equation: where I t denotes the threshold intensity level and I m denotes the maximum intensity of the PSF. p is a filtering coefficient determining the threshold value I t . p is a value between 0 and 1 expressed as a percentage; thus, eq. (2) indicates the percentage of the maximum intensity of the polar PSF that is assigned to the threshold intensity. PSF pixels whose intensities are below this threshold are filtered out. The criteria for setting a reasonable coefficient p are different between different FOVs. For the radially expanding blur generated by our SSLIS, the PSFs at low FOVs require a small p such that I t is also small with respect to I m . This is because the majority of the nonzero pixels of a low-FOV PSF are distributed around a very small area, and only a small fraction of the nonzero pixels, whose intensities account for a very small percentage of I m , are scattered outside this area. Therefore, this small fraction of nonzero pixels can be filtered out by setting a small p. As a result, the number of nonzero entries in the BTCB is very small compared with the numbers of rows and columns of the polar image; thus, the BTCB is a sparse matrix. In contrast, the nonzero pixels of high-FOV PSFs whose intensities are not small with respect to to I m are scattered around a large area. In this case, p should be sufficiently large to filter out those pixels below the threshold intensity in order to reduce the area of nonzero pixels, thus increasing the sparsity of the corresponding BTCB matrix. This can be explained visually by observing filtering for PSFs at semi-FOVs of 3 and 14°, and the change in the sparsity of their BTCB matrices before and after filtering in Fig. 3.
The black areas shown in Figs. 3(a) and 3(b) indicate PSF pixels with zero intensity, while the other areas represent nonzero pixels. The brighter the nonzero pixels, the nonzero pixels whose higher their intensities with respect to the maximum intensity I m . It can be confirmed from Fig. 3(a) that, in the case of the PSF at 3°, the nonzero pixels are almost completely distributed around the polar PSF center and only a small fraction of pixels with a very low intensity are scattered elsewhere. We assigned a small coefficient such as p = 0.8% so that pixels whose intensities are below p = 0.8% of the maximum intensity were filtered out. The β/n ratios before and after PSF filtering are 8.37 and 1.08%, respectively, which are marked in Fig. 3(a). This suggests that the sparsity of the BTCB matrix of the polar PSF at a semi-FOV of 3° was increased effectively by setting a small coefficient. As the semi-FOV increases to 14°, we note from Fig. 3(b) that there are more bright nonzero pixels in the polar PSF at 14° that are scattered across a wider area than they are in the polar PSF at 3°. The intensities of these pixels are actually not low compared with I m . We have to set a large coefficient such as p = 10% to filter out some nonzero pixels this time, otherwise the sparsity of the corresponding BTCB matrix cannot be increased effectively. The values of sparsity before and after PSF filtering are marked in Fig. 3(b), and are 7.49 and 2.16%, respectively.
Note that the β/n ratio of the polar PSF at 14° before filtering is smaller than that of the polar PSF at 3° because there are a large number of nonzero pixels with a very low intensity in the latter, which contribute more to the bandwidth β than those in the former. Since the polar image size does not change for different PSFs during image restoration, the β/n ratio of the latter is larger than that of the former.
The polar PSFs before and after filtering are presented in a 3D view in Figs. 3(c) and 3(d), respectively. We give a 3D view because the changes in PSF shape and intensity due to PSF filtering cannot be observed directly in a 2D view. In Figs. 3(c) and 3(d), the PSFs before and after filtering have the same scale in the x-, y-, and z-directions. We can observe that the original shapes of the polar PSFs at 3 and 14° did not change significantly; only the pixels with intensity below the threshold were filtered out. The pixels with intensities above the threshold were unmodified.
p can be chosen systematically as follows: 1) Select PSFs at some discrete points of semi-FOV within the FOV range of the SSLIS. 2) For each PSF selected in 1), restore the blurred image using this PSF by setting different values of p. Obtain an optimal p by trial and error. The optimal p can be considered as a balance point where the error due to the nonsparsity in the form of ringing artifacts is not high, and the degree of smoothness of the restored image due to weakening of the PSF restoration kernel by modifying the original PSF shape is not high either. 3) Quantify the original PSF size by measuring the difference between the pixel values of this PSF and a black background of the same resolution whose pixel values are all zeros. This can be calculated using root mean square error (RMSE), which is defined as where o(i, j) and b(i, j) denote pixels of the original PSF and the black background with zero entries, respectively. N and M define the resolution of the PSF as well as the black background. A PSF at low FOV will result in a small RMSE because it is small and vice versa. 4) Calculate the size shrinkage after filtering the PSF using the optimal p. This time RMSE can be defined as where ƒ(i, j) represents pixels of the filtered PSF. This equation defines the similarity between before and after PSF filtering. A small RMSE indicates that the optimal p is small so that the size shrinkage is also small after filtering and vice versa. 5) Fit the data from the selected discrete PSFs to a polynomial function. The horizontal axis is the original size of the PSF quantified by RMSE in step 3), the vertical axis is the size shrinkage after filtering by the optimal p, which is calculated in step 4). 6) Use the polynomial function obtained in 5) to determine an optimal p for PSF at any FOV according to the original size of the PSF.
In the next section, we will investigate the effect of p on the restored image quality for different FOVs and give an example on how to determine an optimal p for PSF at any FOV by using the six steps described above.

Results
To evaluate the PSF filtering technique, we used three test images captured by the SSLIS, as shown in Fig. 4.
We selected PSFs at three semi-FOVs of 3, 7, and 14° to restore each captured image. For each PSF, we set three values for the filtering coefficient p to determine its effects on the sparsity of the PSF BTCB matrix and the PSF shape. Table 1 gives the three values of p for these PSFs and the β/n ratios before and after PSF filtering. Two-dimensional views of the PSFs are plotted in Fig. 5. The images resteored obtained using PSFs at the Fig. 4. Three test images captured by the SSLIS camera module.  Table 1 and Fig. 5 suggest that, for PSFs at 3 and 7°, setting a large filtering coefficient p such as 10 or 20% will destroy the Gaussian-shaped PSF, leaving only a "cylinder"-shaped PSF, even though the β/n ratio is significantly reduced from 8.37 to 0.14% and 0 for PSFs at 3 and 7°, respectively. This weakens the effect of the PSF as a deblurring (restoration) kernel. As a result, the images resteored show little improvement in terms of deblurring compared with the captured image, which can be confirmed from (a), (c), (d), (f), and (g) in Figs. 6-8.  It can be observed from Fig. 5 that the shape of the PSFs at 3 and 7° after filtering using p = 0.8% did not change significantly. Filtering reduced only the number of nonzero pixels that are below 0.8% of the maximum intensity I m of these PSFs. We can confirm from Table 1 that the β/n ratios are reduced from 8.37 to 1.08% and from 7.83 to 2.36% for the PSFs at 3 and 7°, respectively, indicating that the sparsity after PSF filtering greatly increased. The images restored using p = 0.8% are also superior to those using other values of p, which can be observed from (b) and (e) in Figs. 6-8. In contrast, setting a small p such as p = 0.8% for the PSF at 14° changed neither the sparsity nor shape of the PSF, which can be confirmed from Table 1 and Fig. 5.
As a result, the images shown in (h) in Figs. 6-8 were not well restored. The image quality was reduced by severe ringing artifacts due to the large β/n ratio. It can be observed from Fig. 5 that the intensities of the PSF pixels away from the PSF center are very low for the PSFs at 3 and 7°, but not for the PSF at 14° before filtering. We have measured and marked in Fig. 5 the intensities of those nonzero pixels away from the PSF center and found that these pixels of the PSF at 14° account for 2 to 8% of the maximum intensity I m , which is the reason why a small p failed to increase its sparsity. The sparsity was increased by setting a large p such as 10 or 20%, as shown in Table 1, which reduced the β/n ratio from 7.49 to 2.16 and 0.95%, respectively. However, we can see from  that the shape of the PSF at 14° is greatly changed from its original Gaussian shape when we set p = 20%. This also results in poor restoration results, as shown in (j) in Figs. 6 -8. The best results for the PSF at 14° are obtained by setting p = 10% and the images resteored are presented in (i) in Figs. 6-8. We can conclude that the PSF filtering technique is effective for reducing the β/n ratio and hence increasing the sparsity of a PSF BTCB matrix. The filtering coefficient p is a function of FOV. The restored image quality can be improved if p is suitably set for different FOVs. An acceptable p should not only increase the sparsity but also not markedly modify the original PSF shape.
In the following paragraphs, we show an example on how to obtain an optimal p for a PSF at any FOV using the six steps described at the end of § 2. Test image 1 will be used in this example.
In step 2), we aim to obtain the optimal p for each of the six PSFs by trial and error. We found that the optimal p values for PSFs at semi-FOVs of 3, 7, 10, and 14° are approximately 0.8, 1.5, 7, and 10%, respectively. The images resteored using different p values for the four PSFs are presented in Figs. 9-12. The images resteored using Fig. 9. Determine the optimal filtering coefficient p at semi-FOV of 3° by trial and error. the optimal p demonstrate fewer ringing artifacts than those restored using a smaller p because of an increased sparsity, and a better restoration result than in the case of using a larger p is obtained because the latter greatly modifies the original PSF shape (we show the 2D view of the PSF before and after filtering when using a larger p in Figs. [9][10][11][12], which weakens the effect of the PSF restoration kernel and smoothens the resulting image. As to PSF at 16 or 18°, we cannot find an optimal p by trial and error because the images resteored either present significant ringing artifacts across the image or demonstrates no improvement in terms of restoration (See Figs. 13 and 14).
In step 4), we calculated the RMSEs for the four PSFs after filtering using the optimal to be 0.1255, 0.3680, 1.5230, and 4.3728, respectively, using eq. (4).
In step 5), we drew a graph using the RMSE data of the four PSFs obtained in steps 3) and 4). It is found that the relationship can be represented by a linear function, which is marked in Fig. 15.
In step 6), we can obtain an optimal p for a PSF at any FOV within the range of 0 to 14° using the linear function obtained in step 5). For example, we want to determine the optimal p for the PSFs at semi-FOVs of 8 and 12°. The original PSF sizes by RMSE are 1.8751 and 5.4549, respectively. We can calculate from the linear function that the size shrinkages for the PSFs at 8 and 12° are 0.7811 and 2.6877, respectively. The corresponding optimal p values are approximately 5.2 and 8.4%, respectively. Figures  16 and 17 give the images resteored using the calculated optimal p.

Discussion
The blur restoration in this paper used a single PSF to restore the whole test image, which is not desirable when we want to obtain a good restored image. The radial blur of an SSLIS should be restored by multiple PSFs at different FOVs because of its radially expanding nature, which can be achieved by segmenting the panoramic polar image into different regions in which the PSF applied to each region can be regarded as locally invariant. Note that, in (i) in Figs. 6-8, the regions near the image center were poorly restored by using the PSF at 14° because these regions are formed by low-FOV PSFs in the image formation process. Therefore, these regions should be restored by low-FOV PSFs. Details of restoration by segmenting a polar image and using multiple PSFs are documented in our previous works. (12,14) The focus of this work is to demonstrate the effectiveness of the PSF filtering technique; thus, multiple PSF restoration by image   segmentation is not discussed here. The six steps discussed above can also be used to determine the optimal p for each PSF used in the multiple PSF restoration process.
On the other hand, although we are eager to compare the images resteored with the original scene using a quantitative means of evaluation such as the RMSE, it is difficult to find an original scene that is free of any distortion or blur. Therefore, we merely evaluated the restored image visually in this paper.   Furthermore, different imaging systems have different PSF distributions and hence the relationships between the original PSF size and the size shrinkage are also different. However, we believe that the six steps described in § 2 are not limit to the radially variant blur produced by the SSLIS, but can also be applied to other imaging systems with an different PSF distributions. We will validate this in the future.
Finally, PSF filtering may introduce an abrupt change in intensity at the threshold intensity level I t , such as for the PSF at 14° filtered by p = 10 and 20% in Fig. 5. This abrupt change in intensity can be considered to be at the edges of these filtered PSFs. When I t is not small with respect to I m , these edges also introduce ringing artifacts to the restored image. We could use a ramp or higher-order function to treat the PSF pixels with an intensity below I t instead of treating them as zeros such that the edges of the filtered PSF can be smoothed.

Conclusions
In this paper, we introduced a PSF filtering technique for restoring a radially expanding blur generated by a SSLIS. The proposed technique increases the sparsity of a PSF BTCB matrix by filtering out PSF pixels whose intensities are below the threshold intensity. The effectiveness of this technique was proved by the restoration of three images captured by the SSLIS camera module. The investigation of different filtering coefficients for PSFs at different FOVs suggested that good restoration can be obtained when the original PSF shape is not markedly modified. A systematic way to determine an optimal filtering coefficient for a PSF at any FOV was introduced as well. Our future work will focus on the evaluation of this technique for other imaging systems that have different PSF distributions and involve a nonsparse BTTB or BTCB matrix.