Modeling nonstationary lens blur using eigen blur kernels for restoration

: Images acquired through a lens show nonstationary blur due to defocus and optical aberrations. This paper presents a method for accurately modeling nonstationary lens blur using eigen blur kernels obtained from samples of blur kernels through principal component analysis. Pixelwise variant nonstationary lens blur is expressed as a linear combination of stationary blur by eigen blur kernels. Operations that represent nonstationary blur can be implemented efficiently using the discrete Fourier transform. The proposed method provides a more accurate and efficient approach to modeling nonstationary blur compared with a widely used method called the efficient filter flow, which assumes stationarity within image regions. The proposed eigen blur kernel-based modeling is applied to total variation restoration of nonstationary lens blur. Accurate and efficient modeling of blur leads to improved restoration performance. The proposed method can be applied to model various nonstationary degradations of image acquisition processes, where degradation information is available only at some sparse pixel locations.


Introduction
Image acquisition processes can be modeled by the following system of linear equations: where b, x, and v are vector representations of the observed, original, and noise images, respectively. Matrix A represents the degradation introduced by an image acquisition process.
In this paper, we address blur introduced by a lens due to defocus and optical aberrations. The severity and characteristics of the blur can change based on the focusing accuracy and sensor pixel locations. Defocus and optical aberrations (e.g., spherical aberrations, coma, astigmatism, and field curvature) are responsible for the nonstationarity of the blur [1]. Modern lenses are designed to reduce optical aberrations through careful arrangements of optical elements and use of aspherical elements [2]. However, images acquired using modern lenses still suffer from lens blur-even when images are properly focused and acquired with no camera or subject motion [3][4][5][6][7][8][9][10]. Thus, unlike blur due to camera shake or subject motion, for which nonstationary blur is determined at exposure by the amount of motion or shake, [11,12], lens blur can be precalibrated once and subsequently used to restore photographed images [5,7,8,10]. A nonstationary blur operation can be written as a filtering operation using spatially variant blur kernels. However, constructing a matrix A that represents the filtering with pixelwise variant blur kernels is difficult. Moreover, matrix A that represents nonstationary blur is neither a Toepliz nor a block circulant matrix. Matrix A is not diagonalized by the discrete Fourier transform (DFT) [13,14]; consequently, any efficient restoration algorithms that rely on the diagonalization of matrix A by the DFT [15][16][17][18][19] cannot be applied to restore nonstationary blur. In many practical situations, nonstationary blur information can be acquired through the available blur kernels, estimation, or measurements at sparse pixel locations [4,5,7,8]. The so-called efficient filter flow (EFF) [20,21] method is often adopted to model nonstationary blur by assuming stationarity inside the regions of an image where blur kernels are available. Because the stationary blur operation inside each region can be implemented efficiently with the DFT, the EFF has been widely used to model nonstationary blur for efficient restoration [4,5,7,8,10,22,23]. However, the following problems have not been addressed: (1) the accuracy with which the EFF can model nonstationary blur and (2) how the modeling accuracy of the EFF affects restoration performance.
This paper presents an efficient modeling method to model nonstationary blur. Blur at each pixel location is expressed as a linear combination of eigen blur kernels, which are obtained through principal component analysis (PCA) of blur kernels available at sparse pixel locations. Blur operations using the eigen blur kernels can be implemented efficiently with the DFT. Moreover, the proposed method requires no region-wise processing, which can affect modeling accuracy and processing efficiency. PCA has been used to model point spread functions (psf's) of optical systems. In [24,25] a three dimensional (3D) depth-variant psf of a microscope is modeled using principal components of 3D blur kernels measured at various object distances. In astronomical imaging, psf at a star location is modeled as a linear combination of point source images, basis functions, or principal components of point source images [26,27]. We use the PCA to model spatially variant nonstationary lens blur caused by defocus and optical aberrations. Through controlled experiments, in which pixelwise nonstationary blur models were available, we showed that nonstationary lens blur could be modeled more accurately and more efficiently when using the eigen blur kernels than when using the EFF. We applied the proposed modeling method to total variation (TV) restoration [28][29][30] and showed that more accurate and efficient modeling leads to improved restoration performance. We also applied the proposed method to model nonstationary blur for real-world lenses to restore defocus and optical aberrations.
The proposed eigen blur kernel-based modeling method can be applied to model various nonstationary degradations that occur during image acquisition. When the degradation changes characteristics in a pixelwise fashion, and information concerning such degradation is estimated or measured at sparse pixel locations for practical consideration, the proposed method provides a method to accurately model and efficiently implement the degradation processes for accurate and efficient restoration.
The remainder of this paper is organized as follows. In Section 2.1, nonstationary lens blur due to defocus and optical aberrations is introduced. In Section 2.2, we present the modeling of nonstationary lens blur by the eigen blur kernels. We introduce TV restoration in Section 3.1. Section 3.2 describes how to implement nonstationary blur operations used in restoration efficiently with the proposed eigen blur kernels. In Sections 4.1 and 4.2, the modeling accuracy and TV restoration performance of the proposed method are analyzed through controlled experiments using Gaussian blur and thin lens blur models. Section 4.3 presents the modeling and restoration of real-world lenses. Finally, Section 5 concludes the paper.

Nonstationary lens blur
The image for a point light source not at the focal distance forms either in front of or behind a sensor. Consequently, the point source will appear as a circular image called the "circle of confusion." Any object that does not provide an acceptable circle of confusion will appear blurry. The point spread function (psf) that represents the image shape of a point light source changes depending on the pixel locations, aperture size, and focusing accuracy. For example, the psf of a thin lens with a circular exit pupil of radius a at a sensor located at distance z is given by [1] where (r, θ) are the locations at the sensor in polar coordinates, R represents the distance from the exit pupil to the focused image, and ρ = r/a. Function Φ(ρ, ϕ) is called a phase aberration, which is the sum of the primary optical aberration and defocus aberration. The primary optical aberration, which includes spherical aberration, coma, astigmatism, field curvature, and distortion, changes based on the pixel location (r, θ). The defocus aberration changes with the location of sensor z relative to the focused image location R. The overall shape of the psf changes with the defocus, the position of z relative to R, and the spatial location of a pixel on the sensor, (r, θ). Hence, the blur introduced by psfs is nonstationary with respect to the defocus and pixel locations. Nonstationary blur can be represented by a filtering operation with spatially variant blur kernels [31]. Let h(i, j, ξ, ζ; d) be the blur kernel at the pixel location (i, j) in Cartesian coordinates, when the camera focus is behind the scene at distance d. Variables (ξ, ζ) are local spatial variables for the individual blur kernel. An image acquired through a lens can be modeled as follows: where b, x, and v are the observed, original, and noise images, respectively. The acquisition model is a nonstationary process that involves pixelwise variant blur kernels. When the blur is stationary, the acquisition model in (3) is reduced to the convolution by a single spatially invariant blur kernel with additive noise.

Modeling using eigen blur kernels
Accurate modeling of a nonstationary image acquisition process requires knowledge of the pixelwise variant blur kernels. In many practical situations, the blur of an image acquisition system is measured or estimated at some sparse pixel locations. In [4,5,7,8,32], the blur kernels are estimated on grid points using a test pattern. In this section, we present a modeling method that can construct an accurate nonstationary blur model from blur kernels measured or estimated at any focal distance and pixel locations. Figure 1 shows a scheme of the image acquisition system used in our modeling. A test image, a part of which is shown in Fig. 1(a), consists of repeated circular patterns. The test image is photographed using a lens, for which we build a blur model. A camera is positioned perpendicular to the test image and is focused at the focal distance s. The aperture is fixed at a given f-number. We take multiple photographs of the test image after placing the test image at a defocusing distance d without changing the focus. For example, to model the commercially available lens in Section 4.3, we take three photographs with the test image: at the focal distance, in front of the focal distance, and behind the focal distance. This procedure allows us to collect blur kernels at different defocus and pixel locations and construct a nonstationary lens blur model that accounts for both the defocus and primary optical aberrations. Patches of photographed test images containing a group of circular patterns are extracted, and blur kernels are estimated for each patch. The photographed test images are divided into patches of circular patterns, as specified by the four red dots shown in the image in Fig. 1(a). Blur is estimated in a patchwise manner, and we assume that the blur inside the patches is stationary. Under this assumption, the acquired image in a patch is the convolution of a blur kernel and the original test pattern. For the rth patch, the observed image b r is obtained by applying a homography transform to align the four dots of the acquired image to the four dots of the original test pattern x p . The original test pattern x p is available from the definition of the pattern. Blur of each patch image can be estimated from the acquired and original image using a blur estimation method [4,5,7,8,10,33]. Once the blur kernels at various locations are obtained, the proposed method can be applied to model pixelwise variation of nonstationary blur. Subsequently, the modeling of nonstationary blur is independent of a patchwise blur estimation method used. In this work, the nonparametric blur kernels are found by solving the following optimization problem [4,10]: minimize where b r and h r are the vector representations of the observed image b r and the blur kernel h r , respectively. The elements of matrix X p are the pixel values of x p rearranged to represent the convolution. Matrices C d are the horizontal, vertical, and diagonal difference operations, and weights µ d are chosen to be small values. A solution is then found using the MATLAB optimization toolbox. The estimated blur kernels are postprocessed with thresholding and renormalization such that all the values are positive, and small values due to noise are removed. A model that characterizes the shapes and spread of blur kernels of a lens can be constructed from the estimated blur kernels h r . The principal components of the blur kernels estimated through (4) are computed using PCA. Let e k (ξ, ζ) be the kth principal component, which we call an eigen blur kernel. We assume that eigenvalues λ k corresponding to the eigen blur kernels e k (ξ, ζ) are sorted in decreasing order. The blur kernel at the (i, j)th pixel location defocused by d can be approximated as a linear combination of the K eigen blur kernels: where c k (i, j; d) are the weights for the kth eigen blur kernel e k . The approximation of blur kernels by the eigen blur kernels in (5) is available only for the pixel locations (i, j), where we estimated the blur kernel h r . However, assuming that the blur of a lens does not change abruptly, a blur kernel at any pixel location can be represented as a weighted average of the available blur kernels in the neighborhood. Thus, the blur kernel at any (i, j)th pixel location can be written as where ϕ(s, t) is the interpolation weight, and N (i, j) is a set of indices for the neighbors of (i, j) for which estimated blur kernels are available. Substituting (5) into (6), we obtain where Hence, the blur kernel at any pixel location can be found by using the interpolated weights c k (i, j; d) and the eigen blur kernels e k (ξ, ζ).

Total variation restoration
Image restoration is an inverse problem: to obtain the original image x from a degraded image b, which is often ill-posed [34,35]. An optimization problem with constraints that reflect our prior knowledge about images can be used to find a meaningful solution for the ill-posed problem. The solutions can be found by various methods [16,[28][29][30].
In this paper, we consider the restoration problem given by the following optimization problem: where TV(x) is the TV. We consider three restoration methods: the primal-dual algorithm (PA) [29], NESTA [28], and KADMM [30]. The most computationally expensive algorithm operations are Ax and A H x operations. Hence, the way these operations are implemented is vital for efficient restoration algorithms. In this study, we consider the TV restoration approach. However, the following development that shows efficient implementations of Ax and A H x operations can be applied to any restoration methods that require these operations. For example, restoration algorithms that utilize linear operators such as ∥x∥ or ∥Dx∥ instead of the TV also require operations Ax and A H x. https://www.overleaf.com/project/5f20edabff82f10001467968 When blur is stationary, matrix A represents the convolution by a stationary blur kernel and becomes Toeplitz. Moreover, given a periodic boundary condition, matrix A becomes block circulant. Hence, operations Ax and A H x can be implemented using the DFT matrix [13,14]. The operations are where A is a diagonal matrix, A * is the complex conjugate of A, and F is the DFT matrix. When blur is nonstationary, matrix A represents the filtering by the spatially variant blur kernels, and matrix A is neither Toeplitz nor block circulant. Note that the efficient implementations of Ax and A H x using the diagonalization by the DFT matrix are not available for nonstationary blur. The EFF method is often used to implement the operations efficiently when blur is nonstationary [20,21]. Here, the blur in regions of an image is assumed to be stationary, and operations Ax and A H x are approximated by where the sum is for the regions of image x. Matrix W l is a diagonal matrix that picks up pixels in the lth region with an overlapping window, P l crops the pixels in the lth region, F l A l F H l applies the stationary blur to the cropped lth region, and P T l inserts the results. Within the EFF, A eff x and A H eff x operations are implemented with the DFT (applied in a regionwise fashion). Matrix W l is used as a window function to prevent ringing due to the regionwise applications of the DFT and is also responsible for smooth transitions of the blur kernels across the region boundaries. The elements in W l 's are prepared so that the sum of the elements operating on a pixel equals one. Because matrix W l serves these two purposes, designing W l to accurately approximate matrix A is not straightforward.

Nonstationary blur operations using eigen blur kernels
The blur kernel at any pixel location can be written as a linear combination of K eigen blur kernels, as shown in (7), and it can be written in matrix-vector notation as follows: where C k is a diagonal matrix whose elements are c k (i, j; d), and E k represents the convolution by the stationary eigen blur kernel e k . Given a periodic boundary condition, matrix E k is Toeplitz and block circulant. Matrix E k is diagonalized by the DFT matrix as Hence, the multiplication operations can be written as follows: Because C k and E k are diagonal matrices, operations A K x and A H K x can be implemented efficiently using the DFT.
Matrix A can be approximated more accurately by the eigen blur kernel-based implementation when more blur kernels are available. The interpolation of the coefficients in (8) does not require regularly sampled blur kernels. Hence, as we estimate or measure blur kernels at more pixel points, we can acquire more accurate approximations of the nonstationarity. The EFF works with blur kernels at regularly divided regions of an image. Adding more blur kernels results in smaller size regions. Overlapping the smaller regions with the window function W l , which is responsible for both overlapping and the transitions of blur kernels across the region boundaries, may hinder the approximation accuracy.
The dominant arithmetic operations of A K x are the computations of the DFT, which are relative to the size of the image. The computational complexity of using K eigen blur kernels for M × N size images is on the order of K MN log(MN).
For the EFF using S i × S j blur kernels for M × N size images, the dominant arithmetic operations of A eff x are the DFT computations of size where o is the number of overlapping pixels. The computational complexity of the EFF is on the order of The computational complexity of both methods is on the order of MN log q with different q values. However, in terms of implementation, the eigen blur kernel-based method offers simpler data access than the EFF. All the operations involved in the eigen kernel-based approximation are imagewise operations; there are no regionwise accesses. In contrast, the EFF requires regionwise accesses to and from various regions of an image, which complicates the processing flow considerably.
The eigen blur kernel-based operations A K x and A H K x in (13) can be applied to the three TV restoration algorithms that we consider (PA, NESTA, and KADMM) to perform efficient nonstationary blur restoration.

Experimental results and discussions
We investigated the modeling and restoration of nonstationary blur using the proposed method through two controlled experiments and an experiment using real-world lenses. For the controlled experiments, we introduced pixelwise nonstationary blur to images using matrix A, which was constructed following the analytic blur models. Subsequently, the blurred images were restored using the proposed eigen blur kernel-based modeling method using matrix A K . The controlled experiments allow us to evaluate the modeling accuracy quantitatively because the analytic expressions for the pixelwise nonstationary blur are available. Moreover, it allows us to analyze how modeling accuracy affects the restoration performance because both the pixelwise blur models and the original images are available for quantitative evaluation. For the experiments using real-world lenses, we acquired images degraded optically by lenses using a camera. Subsequently, we restored the acquired blurry images using the eigen blur kernel-based modeling, which uses matrix A K . The experiments with real-world lenses allow us to validate the application of the proposed method to calibrate nonstationary lens blur for restoration purposes. For the first controlled experiment, we consider pixelwise variant rotated Gaussian blur. A two-dimensional Gaussian distribution and its variants were used to model nonstationary lens blur [4,9,10]. The Gaussian blur kernels were prepared with the standard deviations in the principal axis set as follows: where r is the radius from the center of the image normalized by half of the diagonal sensor length, and γ(d) determines the severity of nonstationarity depending on the defocus d. Here, α i , α j , β i , and β j are constants. The principal axes of the Gaussian kernels are rotated to align with the radial and tangential directions. The blur kernel h(i, j, ξ, ζ; d) is found for every pixel index (i, j) ∈ [1, M] × [1, N] and for the defocusing distance d.
We simulate a practical situation, in which blur kernels are available through measurements or estimations only at grid points, and a nonstationary blur model is constructed from the available blur kernels. The pixelwise variant Gaussian blur kernels are sampled at S i × S j grid points for every defocusing distance d. Each sampled blur kernel represents the blur inside a bs × bs region of an image. Figure 2 shows the blur kernels sampled at 8 × 8 grid points. The blur kernels simulate the lens blur at three different defocusing distances with γ(d) = 0.5, 1.0, and 2.0. It can be seen that blur kernels change their shapes depending on the defocus and the pixel locations. The blur kernels and their contours at the two locations, indicated with the red squares, are shown for a close inspection.
Eigen blur kernels are computed from the 3 × 8 × 8 sampled Gaussian blur kernels shown in Fig. 2. The PCA is applied, and the eigen blur kernels that correspond to the K largest eigenvalues are selected. The results are evaluated using the mean square error (MSE). Using K eigen blur kernels, the nonstationary blur is modeled using (12). For the interpolated weights in C k , bicubic interpolation is used throughout the experiments.
The pixelwise spatially variant blur matrix A is prepared from the analytic expression for the blur kernels h(i, j, ξ, ζ; d) at a chosen distance d. The blur matrix A serves as a benchmark by which the modeling accuracy of the EFF (using A eff ) and the eigen blur kernels (using A K ) are evaluated. Table 1 shows the modeling accuracy of the EFF method, ∥A − A eff ∥ F and that of the proposed eigen blur kernel-based method using the K eigen blur kernels ∥A − A 16 ∥ F . We consider the cases with image sizes of 512 × 512, 1024 × 1024, and 2560 × 2048 and where blur kernels are available for each block of size 16 × 16, 32 × 32, 64 × 64, 128 × 128, and 256 × 256. We assume that at least 8 samples exist in the shorter dimension. In all the cases, A 16 approximated A more accurately than A eff .  Figure 3(a) shows MSE from approximating the sampled blur kernels using the K eigen blur kernels. As shown in Fig. 3(b), 16 eigen blur kernels accurately approximate the sampled blur kernels.  The Gaussian blur at any pixel location can be approximated using linear combinations of the K eigen blur kernels by (7). Figure 4 shows the interpolated weights c k (i, j; d) computed for the case in which γ(d) = 1.0, and the image size is 512 × 512 for k = 1, 2, . . . , 16. Table 2 compares the average processing time for A eff x and A 16 x operations in terms of clock time. The average of 100 multiplications is shown. The code for the EFF implementation was provided by [12], which we run with our nonstationary lens blur kernels. The eigen blur kernel-based method executes in less time than that required by the EFF in all the cases. Moreover, it executes considerably faster when a large number of blur kernel samples are used. In contrast, the regionwise data accesses complicate the processing flow of the EFF considerably. In particular, the execution of FFT requires the number of data to be the power of two. In some cases, the use of larger number of blur kernel samples increases the number of block accesses without decreasing the size of the blocks. In those cases, the computational complexity increases with the number of blur kernel samples. The computational complexity of the eigen blur kernel-based method does not change based on the number of blur kernel samples, S i × S j ; instead, it changes based on the number of eigen blur kernels, K.

Restoration
We tested TV restoration of nonstationary blur using nonstationary Gaussian blur. The Lena (512 × 512), man (1, 024 × 1024), and cafe (2, 560 × 2, 048) images are used as test images. Test images are shown in Fig. 5(a), (b), and (c). The test images were blurred using matrix A, which represented the pixelwise variant blur. Noise was added such that the blurred signal-to-noise ratio (BSNR) given by where var(Ax) is the variance of the blurred image, and σ 2 v is the noise variance, is equal to 60 dB. Subsequently, TV was restored by PA, NESTA and KADMM. All three restoration algorithms terminated when the relative change became smaller than a threshold, which was set to 10 −5 . The operations in the restoration algorithms were implemented using A, A eff , and A 16 . The restoration using A represented a case where we had exact knowledge and implementation of nonstationary blur for restoration. The restoration results using A served as a benchmark for evaluating how modeling accuracy affected the restoration performance. The restoration using A eff and A 16 represented a case where only blur models built from sparse measurements or estimations for restoration were available. Note that we apply pixelwise variant blur using A to obtain observed images and the models using A eff and A 16 for restoration. Hence, any modeling inaccuracy would degrade the restoration performance. Table 3 shows the performance of the three restoration algorithms. For all cases of restoration algorithms, blur nonstationarity γ, image sizes M × N, and the number of blur kernel samples S i × S j , the MSE values are lower when the restoration algorithms are executed using A 16 than when they are executed using A eff . Moreover, using more blur kernel samples increases the restoration performance when the algorithms are run with A 16 ; however, that is not always true when they are run with A eff . For small images, adding more blur kernels increases overlapping of the EFF window to multiple patches, which may degrade the restoration performance. For example, for 512 × 512 images, modeling by EFF with 32 × 32 blur kernels provided lower     performance than modeling by EFF with 16 × 16 blur kernels. These results indicate that the more accurate modeling by the eigen blur kernel-based method improves the restoration performance. Figure 6 shows the restoration results using the cafe image with γ = 2.0 using 8 × 8 blur kernel samples. The figure shows portions of the image from the lower right side. The improvement in MSEs in this case corresponds to 1.45 dB, 1.33 dB, and 2.50 dB in the peak signal-to-noise ratio (PSNR) for PA, NESTA, and KADMM, respectively. Errors between the original images and the restored images using KADMM are also shown. According to the results, the eigen blur kernel based method restored the edges and details in the original images more closely. Figure 7 shows the evolution of cost functions vs. processing time for the tested restoration algorithms. Note that PA and NESTA run considerably faster when using A 16 than when using A eff . PA requires one Ax operation per iteration, and NESTA requires one Ax and one A H x operation per iteration. The use of the more efficient implementations of A 16 x and A H 16 x operations results in faster restoration. KADMM reduces the solution space to the nth order Krylov subspace; therefore, it requires only n multiplications by A 16 or A eff , where n is typically a small number. The processing times for the two implementations of KADMM are similar. In all cases, the restoration accuracy is higher when the restoration algorithms are run with A 16 than when they are run with A eff .

Modeling
For the second controlled experiment, we consider a blur by a thin lens. We prepared pixelwise variant blur kernels for 256 × 256-pixel images. A blur kernel at every pixel location is obtained through the numerical integration given in (2). The parameters from the numerical example in [1] are used, and a matrix A representing the pixelwise variant blur of a thin lens was prepared.
We simulated a practical situation, in which the blur kernels are available only at grid points. The blur kernels of a thin lens were sampled at 8 × 8 grid points. Each sampled blur kernel represents the blur inside a 32 × 32-pixel region of an image. Figure 8 shows the sampled blur kernels for three different defocusing distances: z = 8.375, 8.750, and 9.125. As shown, the blur kernels change their shapes depending on the defocus distance and the pixel locations. Moreover, the blur kernel shapes are more complicated than the simple Gaussian model we considered in Section 4.1. The blur kernels and their contours at the two locations, indicated with the red squares, are shown for a closer inspection.
Eigen blur kernels were computed from the 3 × 8 × 8 sampled blur kernels, and PCA was applied to find the eigen blur kernels. Figure 9(a) shows the MSE from approximating the sampled blur kernels using the K eigen blur kernels. According to the results, 16 eigen blur kernels accurately approximate the sampled blur kernels. The 16 eigen blur kernels are shown in Fig. 9(b). Table 4 shows the modeling accuracy of the EFF method, ∥A − A eff ∥ F and of the proposed eigen blur kernels using the K eigen blur kernels ∥A − A 16 ∥ F . In all cases, A 16 approximated A more accurately than A eff .      larger support for z = 8.375 and 9.125. The EFF implementation requires more padding to accommodate the larger blur kernels. The processing time of the A eff x operation is affected by the larger size of the blur kernels. The A 16 x operation, which is not affected by the size of blur kernels, is considerably faster than A eff x.

Restoration
We applied TV restoration to restore the blur of a thin lens using the baboon image (256 × 256) shown in Fig. 5(d). The test image is blurred using matrix A prepared using (2) for each pixel, and noise is added for BSNR at 60 dB. Subsequently, TV was restored by NESTA and KADMM. The restoration algorithms terminated when the relative change was smaller than a threshold, which was set to 10 −5 . Table 6 shows the performances of the two restoration methods. For all the tested restoration algorithms and defocusing distances, the MSE values are lower when the restoration algorithms are executed with A 16 than A eff .  Figure 10 shows restoration results using the baboon image with z = 9.135 using 8 × 8 blur kernel samples. The improvement in MSEs in this case corresponds to 1.85 dB and 1.73 dB in the PSNR for NESTA, and KADMM, respectively. Errors between the original images and the restored images using KADMM are also shown. It can be seen that the eigen blur kernel based method restored the edges and details in the original images more closely.

Modeling
We also characterized the blur of a commercially available lens-a Nikon AF-G 50 mm f/1.8 lens. Blur kernels at f-number 1.8 are modeled for a 4016 × 6016-pixel sensor. The test image is photographed at three different distances from the camera. When a camera with a lens with focal length f and aperture A (in f-number) is focused at s, any objects between s − d near and s + d far provides an acceptable circle of confusion c, where  We took three photographs of the test image at three distances, placing the test image at s − 2d near , s, and s + 2d far . This setup allows us to model nonstationary lens blur not only when it is properly focused and but also when it is slightly defocused.
We estimated the blur kernels at 6 × 10 locations and 3 defocus distances using (4) (the estimated blur kernels are shown in Fig. 11). The blur kernels at different defocusing distances and different pixel locations show a variety of shapes and spreads and include some shapes quite different from those of the Gaussian model. The AF 50 mm f/1.8 lens consists of seven lens elements in six groups designed to provide excellent optical properties. However, the lens still suffers from nonstationary lens blur, even when properly focused, revealing blur kernels with complicated shapes.
Eigen blur kernels of the AF-G 50 mm f/1.8 lens are found through PCA applied to the 3 × 6 × 10 blur kernels shown in Fig. 11. Figure 12(a) shows the MSE between the estimated blur kernels and their approximations by the K eigen blur kernels. The average MSE of the blur kernels for a given K is shown. The figure shows that the blur kernels can be approximated accurately using only a small number of eigen blur kernels. Figure 12(b) shows the 16 eigen blur kernels found from the estimated blur kernels.

Restoration
For commercial lens blur, matrix A that represents the pixelwise variant lens blur is not available. However, any image photographed using a lens is an observed image b that is degraded optically by pixelwise variant lens blur and observation noise. We use KADMM for the restoration with A eff and A 16 . Because no original image free of lens blur and noise is available, it is difficult to evaluate the restoration performance quantitively. Therefore, for quantitative evaluation, we use the same photographed test images used to estimate the blur kernels for restoration. We extracted the patches inside of the four red dots of the restored images and compared them with the definition of the circular patterns. We measured the MSE after compensating for the dynamic ranges of the patches. Table 7 shows the average MSE of the patches. For the photographed test images at all three distances, the average MSE of the patches restored using A 16 is higher than that using A eff . Figure 13 shows examples of the restored images with the DCU3 test chart photographed at the three distances. It can be seen blur due to defocus and optical aberrations can be restored using precalibrated nonstationary blur models.

Conclusion
In this study, nonstationary lens blur was modeled by eigen blur kernels, which are the principal components of sampled blur kernels. The proposed eigen blur kernel-based method provides a more accurate and efficient method for modeling nonstationary lens blur compared with a widely used modeling method that assumes regionwise stationary blur. The improved restoration performance achieved by the proposed modeling method was validated with TV restoration. The proposed method can be applied to model various nonstationary degradations of image acquisition processes in practical situations, where degradation information is available only at some sparse pixel locations. We are currently working on modeling nonstationary blur due to camera shake using eigen blur kernels.

Funding
Ulsan National Institute of Science and Technology (1.200033.01); National Research Foundation of Korea (2016R1D1A1B01016041).