Mauricio Delbracio

Pablo Musé

Andrés Almansa
email: andres.almansa@telecom-paristech.fr

Jean-Michel Morel
email: morel@cmla.ens-cachan.fr

The non-parametric sub-pixel local point spread function estimation is a well posed problem

Keywords: subpixel convolution kernel estimation, aliasing, inverse problems, camera quality assessment

Most medium to high quality digital cameras (DSLRs) acquire images at a spatial rate which is several times below the ideal Nyquist rate. For this reason only aliased versions of the cameral point-spread function (psf) can be directly observed. Yet, it can be recovered, at a sub-pixel resolution, by a numerical method. Since the acquisition system is only locally stationary, this psf estimation must be local. This paper presents a theoretical study proving that the sub-pixel psf estimation problem is well-posed even with a single well chosen observation. Indeed, theoretical bounds show that a near-optimal accuracy can be achieved with a calibration pattern mimicking a Bernoulli(0.5) random noise. The physical realization of this psf estimation method is demonstrated in many comparative experiments. They use an algorithm estimating accurately the pattern position and its illumination conditions. Once this accurate registration is obtained, the local psf can be directly computed by inverting a well condi-

Introduction

Extrinsic image blur can be observed when the camera focal distance was wrong, when there are different objects at different depths, or when there is a motion blur. But there is a permanent intrinsic physical camera blur due to light diffraction, sensor resolution, lens aberration, and anti-aliasing filters. Our goal here is to accurately estimate the point spread functionpsf, that models the intrinsic camera blur. This function can be locally interpreted as the response of the camera to a point light.

There are several key applications of psf estimation, among them image super-resolution, image de-blurring and camera quality evaluation. Traditionally sharp psfs are considered to lead to better images, but too sharp psfs (containing significant frequency components beyond the Nyquist frequency) cause aliasing effects that may also affect the quality of digital images. An accurate sub-pixel estimation of the psf is therefore crucial to evaluate the image quality in terms of a trade-off between sharpness and aliasing effects.

Image super-resolution is the longstanding problem of increasing the resolution of an aliased imaging system by interpolating a single-frame, or by fusing together several low-resolution images. For this difficult super-resolution process, an accurate psf is fundamental. Surprisingly, there are many more works on blind de-convolution associated to image restoration or on super-resolution, than on the accurate psf estimation.

Existing psf estimation methods can be classified as blind or non-blind, parametric or non-parametric. Blind methods estimate the psf from a single image or from a set of acquired images, without any knowledge of the scene. On the contrary, non-blind methods use a specially designed calibration pattern. Blind methods endeavor to model features of the latent sharp image and to find by optimization the most suitable kernel that predicts them from the blurry observation. Most of them attempt to detect edges in the blurred image, modeling them as the result of blurring pure step-edge functions [START_REF] Chalmond | Psf estimation for image deblurring[END_REF][START_REF] Luxen | Characterizing image quality: Blind estimation of the point spread function from a single image[END_REF][START_REF] Capel | Image Mosaicing and Super-Resolution (Cphc/Bcs Distinguished Dissertations[END_REF][START_REF] Smith | Psf estimation by gradient descent fit to the esf[END_REF]. However, in real images, a step-edge convolved with the psf kernel is generally not a good model of the observed edges ( [START_REF] Ladjal | Flou et quantification dans les images numriques[END_REF], Chapter 4). Other blind approaches try to estimate the psf based on statistical models of sharp images [START_REF] Chalmond | Psf estimation for image deblurring[END_REF][START_REF] Rooms | Parametric psf estimation via sparseness maximization in the wavelet domain[END_REF][START_REF] Zhang | A single image based blind super-resolution approach[END_REF][START_REF] Sroubek | A unified approach to superresolution and multichannel blind deconvolution[END_REF]. Since the blind estimation is an ill-posed problem (blind source separation), strong kernel smoothness assumptions or, equivalently, very simple parametric models are necessary. These inaccurate approaches are necessary to characterize and to blindly restore images affected by contingent motion or out of focus blur.

Non-blind methods instead address the problem of estimating accurately the inherent camera blur. They rely on photographs of calibration patterns to estimate the psf. These patterns range from pin-hole or slantededge patterns to random noise images. The sub-pixel psf estimation problem is generally treated as ill-posed. Most of its algorithms therefore introduce a psf model constraining the space of possible solutions. Parametric models, priors on the regularity of the psf or on its symmetry are the most current assumptions. However, these a priori assumptions can jeopardize the estimation accuracy.

The ideal calibration pattern that comes to mind would be a perfect pin-hole image simulating a Dirac delta impulse, permitting to directly observe samples of the psf. However, in such an observation the signal to noise ratio would be very low, the spot support being ideally infinitesimal. Furthermore, for producing subpixel psf estimates several sub-pixel-shifted versions of the spot image would be needed. Bar or sine patterns can also help to sample the mtf, but only up to the Nyquist frequency.

The iso 12233 standard [START_REF]Photography -electronic stillpicture cameras -resolution measurements[END_REF] gives a normalized pattern and a procedure for measuring the one-dimensional mtf, i.e. the modulus of the Fourier transform of the system's impulse response (psf) in a particular orientation. This standard is based on the slanted-edge method [START_REF] Reichenbach | Characterizing digital image acquisition devices[END_REF], which is an extension of the step-edge technique to achieve sub pixel resolution on the estimation. By aligning the step-edge slightly off the orthogonal scan direction the effective sampling rate is increased. Also, scan-line averaging successfully suppresses noise and increases signal-to-noise ratio making the estimation more stable. In [START_REF] Zandhuis | Sub-pixel non-parametric psf estimation for image enhancement[END_REF] the authors propose a slantededge non-parametric sub-pixel psf estimation method that admits geometrical distortions. A parametric and non-parametric edge spread function estimation procedure is proposed in [START_REF] Claxton | Measurement of the point-spread function of a noisy imaging system[END_REF]. Non-uniform illumination is also taken into account. However, the differentiation step that gives back the psf requires regularization and therefore loses accuracy. Since the previous methods are based on estimating several one-dimensional responses, several images or symmetry assumptions are needed to reconstruct a full bi-dimensional psf.

The recent method by Joshi et al. [START_REF] Joshi | Psf estimation using sharp edge prediction[END_REF][START_REF] Joshi | Enhancing photographs using content-specific image priors[END_REF] arguably represents the current state of the art of slanted-edge methods. It proposes a flexible blind and non-blind nonparametric local psf estimation algorithm. Its approach is based on the ability to detect edges with sub-pixel accuracy. In order to get a precise local psf a specially designed pattern formed by 120-degrees-arc-step-edges is used. The method directly solves the de-convolution and superresolution problem for a bi-dimensional subpixel psf. To reach a sub-pixel accuracy a penalty term on the norm of the psf gradient is introduced, the inverse problem being ill-posed. As we shall see in Section 5, this penalty causes inaccurate estimates in the high frequency components of the psf. If the observed image is under-sampled, which is highly probable and the reason why a sub-pixel psf estimation will be proposed here, interpolating it tramples high frequency information.

As we shall try to prove, there are two main possible improvements to the Joshi et al. method, and they are linked: one is the use of a random noise pattern and the other is the removal of any regularity term, thus transforming the psf estimation problem into a well-posed problem. The use of random noise patterns with known power spectral density has been explored for mtf estimation in [START_REF] Daniels | Random transparency targets for modulation transfer function measurement in the visible and infrared regions[END_REF]18,[START_REF] Backman | Fast lens testing using random targets[END_REF][START_REF] Backman | Random target method for fast mtf inspection[END_REF]. In an ideal situation, the power spectral density psd(f ) of the observed digital image at frequency f is equal to the input power spectral density psd i (f ) times the squared mtf(f ). The advantage of this procedure is that the mtf can be directly calculated. It does not require knowledge of the particular noise realization, relying only on statistical assumptions. A strong limitation of this approach is that the estimation is done up to half the sampling frequency. Consequently it does not reach a sub-pixel accuracy, and aliasing effects are not taken into account.

In [START_REF] Brauers | Direct PSF estimation using a random noise target[END_REF] a random noise pattern is also used, but in a completely different approach. The acquired image is registered to match the target. Then, by doing deconvolution with the almost flat spectrum noise target, this method succeeds in characterizing locally the psf. However, the method assumes that the camera oversamples the signal, which is a correct hypothesis for the particular multi-spectral-camera-lens system, but unrealistic for a classical optical camera. This method contemplates the possibility of a non-linear light sensor response, but does not correct the non-projective distortion. Again, the question is treated as an ill-posed problem and noise-free kernels are produced by regularization.

Table 1 summarizes some of the existing algorithms for psf estimation. It first gives the abbreviations for the five criteria characterizing calibration methods. The above analysis suggests that an ideal method must be non-blind (NB), with no regularization. The kernel estimation must be 2D, local (L), sub-pixel (sp). The main systematic perturbations in imaging (optical distortion (I), non uniform illumination (I), non linear sensor response (G)) must be corrected when comparing the ideal pattern to the photographed one. In short, an ideal method must be (NB,R,2D-L-sp, DIG) with no (C,P,K). The closest to this ideal in the state of the art is the Joshi et al. method, but it includes a regularization which be shown fatal to the high frequency kernel content. The proposed method here has all "good" features. It shows mathematically and practically that an adequate noise pattern permits to avoid any regularization. The camera kernel is directly recovered from the comparison of the ideal noise pattern to the observed one by the inversion of a well-conditioned matrix. We will also verify that this is not possible with an edge based pattern.

By correctly choosing the calibration target, a subpixel psf estimation is therefore feasible without a priori kernel model, without regularization, and with a single aliased input image capture. Nevertheless, this requires the careful correction of the geometrical distortion, of the non-uniform illumination, and of the non-linearity of the sensor response. In short, with a noise pattern, and thanks to this careful elimination of all bias, the psf sub-pixel estimation becomes wellposed. Theoretical bounds will also demonstrate the quasi-optimality of white noise calibration patterns to that purpose. Given that the psf is space variant, due to lens aberrations or sensor non uniformity, the estimation must be done as local as possible.

This article is organized as follows. Section 2 describes the general mathematical digital camera model used for psf estimation method. Section 3 proposes a mathematical theory of optimal patterns. It studies the optimality of the calibration pattern in terms of the well-posedness of the psf estimation problem, and concludes with the proposition of a near optimal and physically feasible random noise pattern. Section 4 describes all the steps of the proposed psf estimation protocol. In Section 5 experimental results generated with both simulated and real camera data are presented, crossvalidated, and compared with the results of state of the art previous methods. Section 6 is a final discussion.

B N C-P 1-2D - [7] B E-N K 1-2D - [29] NB E - k-2D-L-sp D [9] NB E P-K k-2D-sp I [23] NB E - 1-1D-sp - [16] NB E K 1-2D-L-sp D [10, 18, 2, 1] NB R C 1-1D - [5] NB R K 1-2D-L G Proposed NB R - 1-2D-L-sp DIG

Image Formation Model

An accurate estimation of the psf requires a proper modeling of the digital image formation process. The basic pin-hole camera model consists of a perspective projection of the 3D world scene into the focal plane.

In real cameras, a system of lenses is needed to concentrate the light rays toward the focal point, passing through a finite but non pin-hole aperture. Hence, the perspective projection is followed by geometric distortions, which are always present in any camera/lens system. This process can be faithfully modeled as a diffeomorphism from the focal plane into itself. The blur of the resulting image in the focal plane is modeled by a kernel that captures all psf like effects (diffraction due to finite aperture, lens aberration, olpf, sensor light integration, etc). Finally the resulting analog image is sampled into a discrete image by the sensor array. If we consider that the observed scene is a planar scene u, the perspective projection is reduced to a planar homography that will be denoted by H. The whole image formation process can therefore be summarized in a single equation

(M') v = S 1 (g (F (H(u)) * h)) + n,
where F (•) is the geometric distortion field, h is the convolution kernel due to all psf like effects and g(•) is a monotone non-decreasing function that describes the non-linear sensor response (camera response functioncrf). The operator S 1 is the bi-dimensional ideal sampling operator due to the sensor array, and n models the CCD noise.

The blur kernel h is space variant, but evolves smoothly. Thus, the symbol * is understood as a local convolution product, the kernel h varying smoothly with the position in the image domain.

The model can be further simplified by noticing that, in order to estimate h, the geometric transformation F • H can be considered as a whole: there is not need to estimate separately the projective and nonprojective parts. We shall therefore denote by D the whole geometric transformation, and the image formation model finally is

(M") v = S 1 g (u D * h) + n,
where u D is the geometrically transformed image, namely u D (x) = u(D(x)). This model can be further simplified. Indeed the sampling and the contrast change g commute, so that S 1 g (u D * h) = g (S 1 u D * h) . As we shall see, the contrast change g can be recovered from the image samples. Thus we shall first focus on the simplified formation model

(M) v = S 1 (u D * h) + n,
and explain later on how the g term can be eliminated. The next section discusses the structure of the optical kernel h.

Diffraction-Limited Optical Systems

Ideal optical systems present psfs only caused by the optical light diffraction. In the case where there are no aberrations the diffraction kernel is determined by the shape and size of the aperture, the focal length, and the wavelength of the considered monochromatic light. If the shape and size of the aperture is known, the far field approximation (Fraunhofer diffraction) can be explicitly computed as the square of the Fourier transform modulus of the aperture function [START_REF] Goodman | Introduction To Fourier Optics[END_REF]. As a trivial consequence the psf diffraction kernel is always nonnegative.

Optical aberrations degrade this ideal system where only diffraction is considered, producing larger kernels [START_REF] Williams | Introduction to the Optical Transfer Function[END_REF]. In addition, optical anti-aliasing filtersolpf may be introduced in the camera before sampling. They are typically made of several birefringent crystals that separate a light spot into several divergent light spots, leading to an effect similar to having a larger pixel pitch. An analysis of the filters commonly used in digital cameras can be found in [START_REF] Zhao | Characteristicanalysis of optical low pass filter used in digital camera[END_REF].

In most cameras, the digitization process is performed by a rectangular grid of photo-sensors (ccd or cmos) located on the focal plane. Each photo-sensor integrates the light arriving at a particular exposure time. This sensor light integration can be modeled by a convolution with a kernel h sensor = 1 C , the indicator function of the photo-sensor region C. In [START_REF] Yadid-Pecht | Geometrical modulation transfer function for different pixel active area shapes[END_REF] the author performs a theoretical analysis of the mtf for the active area shape and deduces explicit formulas for the transfer function for cmos pixel arrays with square, rectangular and L shaped active areas, which are regularly used. In conclusion, the unknown kernel h results from the convolution of some three different kernels, all nonnegative.

Optimality criterion and quality measure for calibration patterns

Assume we can unveil exactly the latent sharp image that produced the blurry aliased observation. Then, solving for the psf amounts to solve an inverse problem governed by the image formation model (M). The first step toward solving this problem is to carefully model the re-sampling operator that produced an aliased observation. The inverse problem to be solved can be stated in terms of the re-sampling rate and of the observed pattern image u D , which is a function of the calibration pattern. It follows, as will soon become clear, that the accuracy of the estimation of h depends on how well we can invert an operator that depends on the resampling operator and on the calibration pattern. In this section we show that a nearly optimal conditioning is obtained when the calibration pattern is a realization of a white noise. While this may not be new (noise patterns have been used in non-blind psf estimation, see e.g. [START_REF] Daniels | Random transparency targets for modulation transfer function measurement in the visible and infrared regions[END_REF]18,[START_REF] Backman | Random target method for fast mtf inspection[END_REF][START_REF] Backman | Fast lens testing using random targets[END_REF][START_REF] Brauers | Direct PSF estimation using a random noise target[END_REF]), the novelty presented in this section is that the use of white noise patterns allows one to solve for super-resolved psfs without the need for any regularization, and without any prior model for h. In other words, the system is well posed as long a white noise image is chosen as the calibration pattern.

Inverse problem statement in terms of the re-sampling operator and the calibration pattern

Suppose that h is band-limited within supp( ĥ) = [-δπ, δπ]. If the psf is sampled at a rate s, where s > δ, the Nyquist sampling theorem guarantees a perfect signal reconstruction. We will consider the case where δ > 1, which corresponds to aliased images, as in practice most digital cameras introduce aliasing.

In the following, F denotes the Fourier Transform. We denote by f the Fourier Transform of a function f . The s-Shannon-Whittaker interpolator defined as

I s u(x) = n u(n)sinc( x s -n) is denoted by I s , S s is the s-over- sampling operator S s u(n) = u(s -1 n
) and lpf w is the frequency cut-off low pass filter that cuts the spectrum of a signal to [-wπ, wπ].

Lemma 1 (Discrete Convolution) Let u and h be images in L 2 (R 2 ) such that h is band-limited, i.e. supp( ĥ) = [-sπ, sπ]. Then

u * h = I s (ũ * h) ,
where h = S s h and

ũ = S s lpf s u. Proof Set ũ := lpf s u = F -1 û • 1 [-sπ,sπ] , so that ũ = S s ũ and û = û • 1 [-sπ,sπ] . This implies that u * h = ũ * h. Indeed, F -1 ûĥ = F -1 û • ĥ • 1 [-sπ,sπ] = F -1 ûĥ .
Now, since both û and ĥ are supported in [-sπ, sπ], it follows that

ũ * h = I s S s (ũ) * I s S s (h) = I s (ũ) * I s (h) = I s (ũ * h).

⊓ ⊔

Remark 1 Note that u does not need to be band-limited, only h. Notwithstanding, if we can find the spectral cutoff ũ of u, then this lemma implies that the continuous convolution u * h can be simulated exactly with a discrete set of samples.

Let us denote by S s the s-to-1-sub-sampling operator S s = S 1 I s .

It follows from Lemma 1 that the image formation model (M) can be rewritten in terms of discrete sequences as

v = S s ũD * h + n,
where h and ũD are sampled at rate s such that s > δ for h to be well sampled. The value s is the oversampling rate to the high resolution lattice, where the psf estimation is going to take place, from the 1× sensor grid.

Assuming that n is a zero-mean stationary white Gaussian noise, the kernel samples h can be obtained by solving

arg min h S s ũD * h -v 2 2 (1)
Here, ũD is the the result of the Shannon-sampling on the s× grid of the distorted continuous pattern signal ũD = S s lpf s u(D(x)), and v the blurred degraded digital observation on the camera 1× sensor grid.

As inferred by the above discussion, to estimate the psf by a non-blind method raises the following issues:

to choose a good psf characterization target; -to estimate the function g(•), the non-linear ccd response; -to estimate the geometric deformation D(•); -to generate ũD from the sharp latent pattern image u; -to find numerical algorithms calculating the psf.

So far h is only assumed to be band-limited. The numerical method will recover only a finite number of samples of h, which is well localized, and therefore in practice compactly supported. Strictly speaking h being band limited cannot be compactly supported. However, the error introduced by a restriction on the support will prove negligible in comparison to the other sources of error: image noise, quantization, slight estimation errors of g, D,... The found solution h is experimentally independent from variations of its assumed support.

The problem in (1) can be rewritten in matrix form,

(P) arg min h S s Uh -v 2 2 ,
where U is the convolution matrix by ũD . (This matrix is applied to the sample vector h). Assuming that the observed image v is of size m × n, the sizes of ũD and h are ms × ns and r × r, respectively. The matrix S s is the downsampling matrix of size M × M s 2 , where M = m × n. As mentioned above, we need s > δ to recover h from its samples. Thus, s is an integer greater than δ, which facilitates the construction of the subsampling matrix

(S s u)(m, n) = u(ms, ns). Then S s U is of size M × N , with N = r × r.
The solution of (P) is easily obtained using a least squares estimation procedure, and is given by

h e = (S s U) + v, where (S s U) + = ((S s U) t (S s U))
-1 (S s U) t is the Moore-Penrose pseudo-inverse of (S s U). Depending on the condition number of this matrix, the inversion would be well-posed and the solution would be unique. Since

(S s U)h + n = v,
the estimation error is given by n e = (S s U) + n. The noise has zero-mean, thus the estimator h e is unbiased and its variance is

E n e 2 2 = E (S s U) + n 2 2 = E    M j=1 N i=1 (S s U) + ij n i 2    = M j=1 N i=1 N k=1 (S s U) + ij (S s U) + kj E {n i n k } .
Since n is white and stationary, with zero mean, it follows that

E n e 2 2 = M j=1 N i=1 (S s U) + ij 2 σ 2 n = (S s U) + 2 F σ 2 n ,
where σ 2 n denotes the noise variance, and • F is the Frobenius norm of a matrix.

If all singular values of S s U are non zero, the singular values of (S s U) + are the inverses of the singular values of (S s U). If some singular value is zero, the system is ill posed and the estimation problem cannot be solved, unless some kind of regularization on h is imposed.

Let {σ 1 , σ 2 , . . . , σ N } be the singular values of S s U. Then

(S s U) + 2 F = N i=1 σ -2 i .
In order to minimize the variance of the estimator h e (i.e. to minimize the noise amplification), one has to minimize the function

γ(S s U) := N i=1 σ -2 i
It should be pointed out that γ depends on the rate s and on the samples ũD . The super-resolution rate s is determined by the spectral support of the psf. The sequence ũD depends on the adopted continuous pattern u, on the geometric transformation D (that includes the perspective projection associated to the particular pattern's view) and also to other possible distortions presented in the camera-lens system. Hence, for the s× sub-pixel psf estimation problem, γ measures the quality of any given view of a calibration pattern.

In order to find the best ideal pattern independently of the view and distortion, we will consider first the discrete problem of finding the best sequence ũD , minimizing the γ value. To simplify the notation we write u ij = (ũ D ) ij . This motivates the following definition.

Definition 1 (Optimal digital pattern) Given a kernel support N = r × r and a window observation size M = m × n, the optimal pattern for the s× sub-pixel psf estimation is the digital calibration pattern u * such that u * = arg min

a≤uij ≤b γ(S s U).
where the constraints on u ij are linked to the physical realization of the pattern and to the sensibility of the sensors. (The conclusions of the analysis will prove independent of the particular value of these bounds.)

Characterization of optimal digital calibration patterns

In this section, we derive a lower bound for γ(S s U) that will be used to design calibration patterns. Indeed, it will then be shown that for a realization of white stationary Bernoulli noise, the γ value is so close to this bound, that in practice these patterns can be considered to be optimal.

Lemma 2 Let Φ be a M × N matrix, M > N , with all its entries in [a, b]. Let σ 1 ≥ σ 2 ≥ • • • ≥ σ N denote its singular values. Then N i=2 σ 2 i ≤ (a + b) √ M N σ 1 -σ 2 1 -abM N.
Proof Let ϕ ij be the (i, j) entry of Φ, and ϕ j its j-th column. Let also φj = 1 M M i=1 ϕ ij denote each column's mean, and φj = ϕ jφj 1.

The Frobenius norm of Φ can be expressed as

Φ 2 F = trace(Φ t Φ) = N j=1 ϕ t j ϕ j . Since φij := ϕij -a b-a ∈ [0, 1],
we have φ2 ij ≤ φij , and then

ϕ t j ϕ j = M i=1 ϕ 2 ij ≤ M i=1 (a+b)ϕ ij -abM = M (a+b) φj -abM Thus, Φ 2 F ≤ M (a + b) N j=1 φj -abM N. (2)
On the other hand, for all

x such that x = 1, Φ 2 ≥ Φx . Let us take x = 1 √ N 1. Then Φ 2 2 ≥ Φx 2 2 = 1 N N j=1 ϕ j 2 = 1 N N j=1 φj + N j=1 φj 1 2 = 1 N   N j=1 φj 2 + N j=1 φj 1 2   + 1 N   N j=1 φj   t   N j=1 φj 1   = 1 N   N j=1 φj 2 + N j=1 φj 1 2   ≥ 1 N N j=1 φj 1 2 = M N   N j=1 φj   2 . (3)
Thus, by Eqs. ( 2) and (3), we have:

Φ 2 ≥ Φ 2 F + abM N √ M N (a + b)
.

Then, since Φ 2 = σ 1 , σ 1 ≥ Φ 2 F + abM N √ M N (a + b) = N i=1 σ 2 i + abM N √ M N (a + b) Finally, N i=2 σ 2 i ≤ (a + b) √ M N σ 1 -σ 2 1 -abM N. ⊓ ⊔ Lemma 3 (A bound on γ) Let Φ be a M ×N matrix, M > N , with all its entries ϕ ij in [a, b]. Then min ϕij ∈[a,b] γ(Φ) ≥ 1 M N 1 b 2 + 4(N -1) 2 (b -a) 2 .
Proof According to Lemma 2, for any matrix Φ with entries in [a, b], and in particular for the ones that attains

γ * = min σ1,...,σN N i=1 σ -2 i , the inequality N i=2 σ 2 i ≤ (a + b) √ M N σ 1 -σ 2 1 -abM N holds. Thus min ϕij ∈[0,1] γ(Φ) ≥ min σ∈D f (σ),
where σ = (σ 1 , . . . , σ N ),

D := {σ | σ i ≥ 0, N i=1 σ 2 i -(a+b) √ M N σ 1 +abM N ≤ 0},
and

f (σ) := N i=1 σ -2 i .
The function f being strictly convex on D, which is itself a convex and compact domain, it follows that the minimum of f on D is attained at a unique point. D and f being invariant by any permutation of σ 2 , . . . , σ N , the minimum point being unique satisfies

σ 2 = • • • = σ N . Since this minimum belongs to D, N i=2 σ 2 i = (N -1)σ 2 2 ≤ (a + b) √ M N σ 1 -σ 2 1 -abM N.
By noting that the maximum value of

σ 1 → (a + b) √ M N σ 1 -σ 2 1 -abM N is ( b-a 2 ) 2 M N , it follows that σ 2 2 ≤ b -a 2 2 M N N -1 .
On the other hand for any point of D we have

σ 2 1 ≤ (a + b) √ M N σ 1 -abM N. Then, it follows that σ 2 1 ≤ b 2 M N . Consequently, min σ∈D f (σ) = N i=1 σ -2 i ≥ 1 M N 1 b 2 + 4(N -1) 2 (b -a) 2 .

⊓ ⊔

Remark 2 It should be noted that in the proof of the previous lemma, the condition that the entries of Φ belong to [a, b] was replaced by the weaker condition given by the inequality proved in Lemma 2. This amounts to enlarge the space of matrices that was originally considered, thus the real optimum that can be attained by matrices with entries in [a, b] will necessarily lead to higher values of γ.

Remark 3 Notice also that in Lemma 3 we did not solve the complete constrained optimization problem min σ1,...,σN

N i=1 σ -2 i subject to N i=2 σ 2 i ≤ (a+b) √ M N σ 1 - σ 2
1 -abM N . While this problem can be solved via the Karush-Kuhn-Tucker conditions, according to the previous remark it would still lead to a lower bound on γ. The solution of this constrained minimization problem leads to a closed form which is significantly less handier than the bound that was obtained in Lemma 3, and is worthless since both bounds are extremely close, as shown in Figure 1 for a = 0, b = 1. Let S s U be the operator associated to the convolution of the r × r kernel with the image u, followed by the downsampling operator of rate s. Then

min a≤uij ≤b γ(S s U) ≥ 1 M N 1 b 2 + 4(N -1) 2 (b -a) 2 .
where M = m × n is the observation window size and N = r × r is the kernel size .

Proof The result follows directly from Lemma 3, since the operator S s U associated to u is a M × N matrix with all its entries in [a, b].

⊓ ⊔

We will propose as calibration pattern a realization of a white Bernoulli(0.5) stationary noise. It will be shown that this calibration pattern is so close to the γ(S s U) lower bound given by Lemma 3 that for practical calibration purposes, it can considered to be optimal.

The motivation for choosing stationary white noise patterns is not new: white noise has been widely used for system identification applications, since it optimizes the minimum variance of unbiased estimators. Now, the choice of Bernoulli(0.5) distribution can be explained as follows. Suppose u = {u ij }, where u ij ∈ [a, b] are mutually independent random variables, identically distributed with mean m u and variance σ 2 u . In this case, it can easily be shown that

E (S s U) t (S s U) = M m 2 u 11 t + σ 2 u I .
This is a direct consequence of the non-correlated nature of u and that subsamples of white noise remain white noise. Observe that M m 2 u 11 t + σ 2 u I has only two different eigenvalues:

σ 1 = M (m 2 u + σ 2 u ) and σ 2 = • • • = σ N = M σ 2
u . Thus, its γ value is

γ = 1 M 1 N m 2 u + σ 2 u + N -1 σ 2 u .
On the one hand, in order to minimize γ, m u and σ 2 u values should be as large as possible. On the other hand there is a trade-off between both values and they cannot be simultaneously maximized. Indeed, any random variable with support [a, b] satisfies

σ 2 u ≤ (m u -a)(b -m u ).
Nonetheless, the equality holds for the Bernoulli distribution. Hence, from now on we restrict the analysis to the Bernoulli case which, from the previous reason, is optimal. Therefore we can express γ as

γ = 1 M 1 N m 2 u + (m u -a)(b -m u ) + N -1 (m u -a)(b -m u )
.

It can be shown that the m u ∈ [a, b] value where γ attains its minimum is always very close to m u = a+b 2 . This happens independently of M and N . However, the exact value depends on N . It is therefore convenient, to avoid dependence on N , to fix m u = a+b 2 by using an equiprobable Bernoulli distribution. Finally, the γ value for the expected operator S s U when using a Bernoulli(0.5) pattern is

γ = 4 M 1 N (a + b) 2 + (b -a) 2 + N -1 (b -a) 2 .
This value is very close to the bound provided by Lemma 3. Indeed, for

M ≥ N ≫ 1 we have γ * -γ ≈ 4 M(b-a) 2
. This small difference is illustrated in Figure 2 for the particular case a = 0, b = 1. Notice also that since E {(S s U) t (S s U)} = lim M→∞ (S s U) t (S s U), large M values may be required in order to reach the optimal γ. However, this is clearly not our case of interest, our goal being to perform a local kernel estimation. Nevertheless, we may still be interested in exploring the use of a realization of white stationary Bernoulli(0.5) noise as calibration pattern, for finite and realistic values of M and N (the non-asymptotic case). In order to show that the choice of such a calibration pattern can be considered to be optimal for practical psf estimation, we generated a white random binary image u ij ∈ {0, 1}, Bernoulli(0.5), and evaluated γ(S s U) for fixed down-sampling rate s = 4. Figure 3 shows that the obtained γ is very close to the nonasymptotic lower bound (Lemma 3), indicating that this pattern is optimal in a practical sense. Concluding Remark The mathematical argument and experiments above show that near-optimal γ values are reached with a Bernoulli random noise pattern for reasonable observation, kernel and pattern sizes. Slightly better γ values could be achieved if we allowed the pattern to adapt to the kernel size. This is nevertheless not practical. The payoff would be a negligible improvement of the well-posedness, and the exact psf support size being anyway a priori unknown.

From continuous patterns to digital patterns

Based on the previous section it comes into view that good psf estimation patterns are those that produce very contrasted random ũD sequences. However, we cannot choose directly the values inside the S s U operator. Indeed, the γ value depends on ũD , obtained by sampling on the s× grid the distorted continuous pattern image.

Consider the set of analogical patterns formed of constant u ij gray value squares regions, u(x) = i,j u ij 1 x-(i,j) ≤ 1 2 . Since signals in optical systems are non-negative in nature and bounded, we can assume w.l.o.g. that 1 ≥ u ij ≥ 0.

For the mathematical exploration of optimal patterns, we will restrict ourselves to the case where the geometrical transformation D is a zoom-out with factor t -1 , Z t -1 . This assumption will almost be verified, the views of the pattern taken frontally. Notice that the s-sampling operator can be written as

S s = S 1 Z s . Thus, ũD = S s lpf s u D = S 1 lpf 1 Z s u D = S 1 lpf 1 Z s Z t -1 u = c * u
where c is the digital filter

c i,j = 1 2 -1 2 1 2 -1 2 sinc( s t ζ -i)sinc( s t η -j)dζdη.
As mentioned earlier the goal is to produce values (ũ D ) ij as independent and contrasted as possible. This motivates the following simplification. Suppose that the set up realizes t = s. An ideal unattainable situation would be that the re-sampling operator and the low-pass filter do not produce inter-symbol interference (i.e. the discrete filter c does not change the input signal u). Then each of the square gray values would be equal to the sample after low-pass filtering ũD ≈ u ij . In this particular case we would have a perfect one-to-one correspondence between the gray values of the pattern and the ũD digital signal which would be a Bernoulli pattern. Due to the constraints on u ij the best we can do is to choose iid Bernoulli(0.5) u ∈ {0, 1}. Yet, while this perfect geometric situation is unattainable, the experiments show that γ stays close to its optimal value when s/t is between 0.7 and 2, as it is shown in in Figure 4. The resulting ũDij for distances in a range from s/t = 1 produce γ values close to the γ bound for entries in [0, 1]. . The zoom factor s/t is closely related to the focal distance and to the distance from the camera to the pattern. For example if the distance to the pattern is too small (small s/t value) the pattern will look like a step-edge pattern because of the zoom-in. The corresponding γ value will be higher than the optimal. On the other hand, if the distance to the pattern is large then the γ value will also be larger than the optimal one, because of the contrast loss due to the zoom-out. In agreement with the theoretical study, the views with zoom factors close to one (i.e. t ≈ s) produce the best γ values.

Comparison of calibration patterns

The γ factor introduced above permits to compare the suitability of different patterns for the psf estimation problem. Since the noise amplification is governed by the sum of the inverses of the singular values, it is desirable to use patterns that produce singular values that are all as large as possible. For this purpose, and justified by the previous theory, we shall use a binary random pattern. The proposed noise pattern consists of a matrix of 256 × 256 black and white random squares generated from an equiprobable Bernoulli distribution. The pattern was printed at a high enough resolution so that artifacts introduced by the printer could be neglected. Several cross marks and white/black flat regions were added, to easily align the acquired image with the target, and to correct non-uniform illumination.

The complete psf Estimation Procedure

In this section we describe the steps that lead to a local sub-pixel psf estimates. The complete chain is summarized in the block diagram of Fig. 7. The next paragraphs present brief summaries for each block.

Feature Detection

In order to deal with geometric distortions the ideal pattern and its observation have to be precisely locally aligned. To that purpose checkerboard corners were introduced along the boundary of the noise calibration pattern. Assuming that the psf is (approximately) symmetric, these x-corners will not suffer from shrinkage. Several methods to detect checkerboard corners have been reported in the Computer Vision literature (e.g. [START_REF] Harris | A combined corner and edge detector[END_REF], [START_REF] Cheng | Evaluating harris method in camera calibration[END_REF], [START_REF] Lucchese | Using saddle points for subpixel feature detection in camera calibration targets[END_REF]), ranging from differential operators such as the Harris detector to more specific correlation methods. We used a Harris-Stephens based corner detector implemented by Bouguet [START_REF] Bouguet | Camera calibration toolbox for matlab[END_REF], that allows to iteratively refine the detected corner positions to reach sub-pixel accuracy.

Geometric Transform Estimation

The estimation of the psf does not require a decomposition of the distortion into its homography and nonhomography parts, as it is done in classical geometric camera calibration [START_REF] Zhang | A flexible new technique for camera calibration[END_REF], where a global radial lens distortion model is usually adopted. In order to avoid that computation and to utilize a more flexible model that may capture local lens distortion, the complete geometric distortion was approximated with thin-plate splines. While thin-plates splines were originally conceived as an exact interpolation method [START_REF] Bookstein | Principal warps: Thin-plate splines and the decomposition of deformations[END_REF] they can be easily extended to the approximation problem [START_REF] Sprengel | Thin-plate spline approximation for image registration[END_REF]. The mapping from the non-distorted to the distorted space is estimated from the detected corners {p i }, and their correspondences in the ideal pattern {P i }, whose coordinates are perfectly known.

Illumination Estimation and Normalization

In order to match the gray levels in the sharp pattern to those in the observed image, black and white square flat regions were included along the boundary of the noise pattern. These regions permit to estimate the mapping between black and white colors and the corresponding observed gray level values. The presence of these constant regions all around the pattern permit to estimate a black (white) image that models the black (white) intensity level at each pixel. These light images have been modeled by second order polynomials whose coefficients are estimated by least squares from the known pairs (value, position). In continuation each pixel value in the observed image is linearly rescaled within the range [0, 1], by considering the respective estimated black and white values.

CRF Estimation g(•)

Once the nonuniform illumination has been compensated, the camera response function can finally be estimated and the non-linear response of the sensors corrected. The estimation and correction procedure is based on a strong property of our pattern: the white noise pattern was generated assigning equal probabilities to black and white values (0 and 1 respectively, after normalization). Consequently, since the psf has unit area, the mean gray value within the observed image should be 0.5.

The solution is defined as a parabolic function u → αu 2 + (1α)u where α is chosen so that the the mean of the pattern after the correction is 1/2.

Pattern Rasterization

In order to generate the samples ũD from the ideal continuous pattern image u, we need to sample this image at the desired s resolution after deforming it by the estimated geometric transformation. For that purpose the distorted continuous pattern u D must be low pass filtered to be bandlimited in [sπ, sπ]. (Remember that the camera resolution is 1×. Thus the digital pattern has an s× over-sampling). The procedure is:

1. The continuous pattern u is sampled at a very high resolution. From the vectorial description of the pattern a digital image is generated (this procedure is called rasterization) by replacing each one of the flat squares by a 4 × 4 block of pixels with the same gray value. The re-sampling starts from these samples u instead of the continuous pattern; 2. Frequencies higher than sπ are cut off from the digital pattern u to get ũ; 3. By help of the previously computed geometric distortion the filtered pattern ũ is bi-cubically interpolated at the desired resolution s× ũD .

Numerical Methods for psf Estimation

We have seen that light diffraction, optical low pas filtering, and sensor light integration all produce nonnegative kernels. Thus the estimated psf must be nonnegative. We can therefore constrain the solution to be non negative, thus reducing the space of solutions. Section 5, Fig. 17 shows that not imposing this nonnegativity assumption yields essentially the same results, which in fact verifies the correctness of the proposed image formation model. Hence, we can opt to solve a non-negative least squares, or simply solve a least problem and then threshold the solution to eliminate very little components. Suppose that the local grid pattern observation v has size m × n and that we want to estimate a psf at s× sub-pixel resolution. Also suppose that the estimated support of the psf is inside a r × r image. The matrix S s U corresponding to the s-down-sampling of the convolution with the distorted calibration pattern, has size mn × r 2 . Thus, the problem to be solved can be formally written as

(P ) arg min h S s Uh -v 2 subject to h i ≥ 0, i = 1, . . . , r 2 .
Problem (P ) can be solved using standard convex optimization solvers such as cvx [START_REF] Grant | Cvx: Matlab software for disciplined convex programming (web page and software[END_REF]. A simpler Newton interior point algorithm proposed in [START_REF] Portugal | A comparison of block pivoting and interior-point algorithms for linear least squares problems with nonnegative variables[END_REF] was used and always converged rapidly.

Experimental Results

This section is dedicated to the evaluation of the proposed non-blind sub-pixel psf estimation method, and to the comparison of its performance with two state of the art proposed approaches. We selected a method recently reported in the literature by Joshi et al. [START_REF] Joshi | Psf estimation using sharp edge prediction[END_REF][START_REF] Joshi | Enhancing photographs using content-specific image priors[END_REF], and a mtf commercial software, Imatest [19]. Since we do not have real camera ground truth for the psf, the performance evaluation was first carried out on simulated data. A real psf estimation on real cameras was in continuation tried under varying acquisition conditions. Particular attention was paid to the aliasing effect caused by sampling under the Nyquist frequency.

Simulations for Objective Evaluations

The simulation of the camera acquisition process was as follows. The grid pattern was rasterized at a very high resolution (i.e. 8×), convolved with a psf like kernel (in this case a Gaussian isotropic kernel), and downsampled to get the observed digital image at the camera resolution (i.e. 1×). The kernel was chosen so that the low resolution image presented aliasing artifacts. We also added white Gaussian noise of standard deviation σ = 0.02. We compared the performance of the proposed approach to that of Joshi et al. using their calibration pattern and our implementation of their approach, with different regularization levels. A 4× kernel was estimated for both algorithms from the observed window of size 110 × 110 pixels.

Figure 8 shows the results for 4× psf estimation from the simulated observation. Solutions with the Joshi et al. method with three levels of regularization are presented, along with the proposed approach (which is regularization-free). In this experiment the proposed method significantly outperform's Joshi et al.'s algorithm, achieving a much less noisier estimation. Joshi's algorithm needs a strong regularization to stabilize the estimation and to avoid an amplification of high frequency noise. Consequently, its estimation tends to penalize high frequency components and to produce a biased kernel with amplified lower frequency components. See caption for details. The observed window has 110×110 pixels. The top row shows the kernel estimation and the middle row the difference image between the estimation and ground truth for one of the realizations. The proposed method significantly outperforms the Joshi et al. algorithm, achieving a much less noisier estimation as shown by the difference images and by the peak signal to noise ratios. The bottom row shows central horizontal profiles for all the estimated mtfs (0.5 is the Nyquist frequency). Notice that in the Joshi et al. method the estimation is unstable. The estimates show extremely noisy components for frequencies higher than the sampling frequency, when the amount of kernel regularization is too small. On the other hand, if a strong regularization is imposed, the penalization of the kernel gradient adopted by Joshi et al. tends to produce kernels with under-estimated high frequency components. The method proposed here does not rely on a regularization and produces nonetheless noiseless and unbiased results.

Experiments with Real Camera Images

In this section we present several local 4× psf estimation examples from real camera acquisitions. In all cases a Canon EOS 400D camera provided with a Tamron AF 17-50mm F/2.8 XR Di-II lens was used. The focal length was fixed at 50.0 mm. Based on these experiments the behavior of the proposed method was analyzed with varying camera aperture. The impact of the crf estimation/correction was evaluated, and the psf estimates obtained for the four color channels in the Bayer pattern compared. Variations of the kernel estimates depending on their location in the image were also explored. This was followed by an evaluation of the stability of the estimation procedure, and of the influence of the kernel support size. Finally the results were again compared with the Joshi et al. algorithm and with Imatest, applied to real cameras.

Different Apertures

The estimation was conducted using the proposed random pattern captured at five different apertures. For each acquisition, a 4× psf estimation for one of the green channels ((half of the green pixels of the Bayer matrix) was performed. Results are shown in Figure 9. The estimations were performed at the image center from a window of size 90 × 90 pixels. Notice that kernels at apertures f/32 and f/16 are significantly larger than the rest, as predicted by diffraction theory (see caption for details). An example of the acquired blurry image is shown in Figure 12.

Fig. 11 shows the diffraction-limited mtf for a circular f/5.7 aperture and green monochromatic light (See the end of Sec. 2.1). The estimated response for our camera-lens system at aperture f/5.7 and for the green channel is under the ideal diffraction-limited response. This can be a consequence of the light integration in the sensor array but also of the optical low pass filter specifically included to avoid aliasing.

Estimation of Camera Response Function

This experiment evaluates the impact of the non-linearity of the camera sensors response. To conduct this experiment, the camera response curve was computed using a specially designed pattern for crf estimation. In order to assess the impact of the crf on the psf estimation, the observed image was corrected using the special purpose crf estimate, to compare the results that yield the psf estimation algorithm.

Fig. 13(a) compares the crf estimated using the special purpose pattern with the crf estimate embedded in the proposed psf estimation algorithm. Notice that both estimates are hardly non-linear and extremely close to each other, so the psf estimation algorithm seems to be capable of giving a reasonable crf estimation.

Fig. 13(b) shows the mtfs obtained under four different situations:

psf estimation with embedded crf correction from the raw observed values (psf-crf). -psf estimation without any crf correction from the raw observed values (psf-nocrf). -psf estimation without embedded crf correction from the adjusted values after correction via the special purpose crf estimate (psf-nocrf-eq).

In all cases, the estimation yielded very similar results. The top and bottom rows show respectively the estimated psfs and a few level lines of the corresponding mtfs that prove that the kernels are not exactly axissymmetric. The kernels at apertures f/32 and f/16 are considerably larger than the rest in agreement with diffraction theory. This phenomenon also stands out in the modulus of the estimated psf spectra, which also shows that the psfs/mtfs are not axis symmetric. Figure (b): Vertical cuts of the spectrum modulus. The camera seems to have the sharpest response from apertures f/3.2 to f/12.9. At apertures f/32 and f/16 the camera cuts high frequencies significantly more than the rest, as predicted by diffraction theory. Notice that in all cases, except at aperture f/32, the mtf at the Nyquist frequency (f = 0.5) is significantly greater than zero. Hence, the camera introduces aliasing.

Color Estimation The goal of this experiment is to compare the psf estimates for all four channels from the Bayer raw camera output (two greens, red and blue).

The estimation was performed using the random pattern captured at apertures f/5.7. The results for the 4× psf estimation located in the image center are shown in Fig. 14. It is easily seen that the red psf is larger than the green and the blue one (i.e. produces blurry images). This is reasonable, since the wavelengths associated to red are smaller than the rest. Hence the red diffraction kernel will be larger than the green and blue kernels for the same camera configuration. The differences between the shapes of the red/blue and green psf spectra can be explained by the sensor shape. If we ac- The estimated response for our camera-lens system is under the ideal diffraction-limited response. This can be consequence of the light integration in the sensor array, but also of the optical low pass filter specifically included to avoid aliasing.

cept that the sensor active zone is L-shaped, then by the red/blue sensors in the Bayer pattern will have the same sensor term mtf and will be rotated 45 • with respect to the green channels.

Location Figure 15 displays the 4× psf estimates for one of the green channels, at different image locations, for f/5.7. Kernels closer to image borders are larger and more asymmetrical than the kernel at the image center. This seems to be a consequence of lens aberrations that deteriorate the system performance.

Stability of the estimation procedure

A set of thirteen images of the noise calibration pattern were acquired with exactly the same camera configuration (f/5.7), from similar viewpoints. For each acquisition, the 4× psf of one of the green channels at Both curves coincide, implying that the otf is real and thus the psf is symmetric. This is also seen in the bottom figure that shows that the otf phase is 0 or π. The estimates from the raw gray values with and without crf estimation/compensation (psf-crf and psf-nocrf, resp.) gave very similar results. After compensation of the gray values using an external special purpose estimation of the crf, the psf estimation procedure (psf-nocrf-eq) also led to very similar results.

the image center were estimated. Figure 16 shows the average mtf vertical profile, and its standard deviation band. It is clear from the small value of the standard deviation that the estimation method is highly stable, in agreement with the fact that the corresponding linear system to be inverted is very well-posed. More details are given in Fig. 16 caption. The red psf is larger than the green and the blue ones. Since the wavelengths associated to red are smaller than the rest, the diffraction components for the red channel will be larger than those for green and blue for the same camera configuration. Also notice the differences between the shape of the red/blue and green psf spectra (bottom row). Red and blue mtf seem to be 45 • rotated with respect to the green ones. This symmetrical behavior is plausible for an L-shaped active zone sensor array.

Support We can consider that the proposed approach has only one main parameter: the kernel support size. The choice of this size implies a trade-off between the model validity and the feasibility of the estimation. On the one hand, if the support is too large the kernel estimation will be very noisy, since the γ factor increases with the support size. On the other hand, if the kernel support is too small the considered image formation model will not be accurate.

Fig. 18 shows the 4× psf estimation for various kernel support sizes. All estimations for the supports 17 × 17, 25 × 25 and 33 × 33 turn out to be very close to each other. Nevertheless, the 9 × 9 kernel support does not seem to be large enough to correctly model the psf. Hence, as soon as the support size exceeds such a lower bound, the proposed algorithm does not appear to be sensitive to this parameter. We define the relative mtf sensitivity in a region Ω as: s(Ω) = mean(std(mtf))/mean(mtf) where the mean values are computed inside the region Ω. In this case we construct Ω(threshold) = {x : mtf(x) ≥ threshold}. The relative sensitivity in the whole spectrum does not exceed 0.08 and what is more if the mtf values smaller than 5% are not considered, then the relative sensitivity is less than 3%. The small standard deviation and sensitivity demonstrate the algorithm stability.

Comparison of several methods This section ends up with a comparison between the Joshi et al. method, Imatest, and the proposed approach to non-blind subpixel psf estimation [19]. Imatest is a commercial mtf estimation software. The Imatest estimation is performed from a slanted-edge image and only gives an estimate of the mtf at the direction orthogonal to the slantededge. The estimation was conducted with images taken at aperture f/5.7 with patterns located at the center of the image. Figure 19 shows the horizontal mtf profiles obtained with the Joshi et al. method using various regularization levels, with Imatest, and with the proposed approach for one of the green channels. In the low frequencies Joshi and the proposed approach yield very similar results. However, for higher frequencies the Joshi et al. results vary strongly with the regularization level. The Imatest estimate is quite noisy and does not resolve frequencies above twice the sampling rate. The proposed random pattern algorithm generates much more information than the typical slanted-edge mtf calibration. for one of the green channels. On the low frequencies all algorithms gave very similar estimations, while on the higher frequencies the Joshi et al. estimation depends on the regularization level. Although we did our best to get a noise free mtf estimation from the Imatest software, the final estimation is quite noisy. The Imatest estimation is done from a slanted-edge image and only gives an estimation for the mtf at the slanted-edge orthogonal direction.

Discussion

This work is an attempt to define an optimal nonblind sub-pixel psf estimation method from a single aliased image. The method is successful, but its set up is tight. The pattern must be large enough (some 70cm in our experiments), printed with good quality ink. The random squares must be large enough to avoid any ink soaking bias, and a good quality print is recommended. The mathematical analysis demonstrated that a Bernoulli pattern is nearly optimal in terms of wellconditioning of the matrix to inverse. The pattern was therefore placed in an approximately frontal position. The photographs were taken at the right distance to ensure that the camera sampling grid and the pattern grid had similar meshes. These position requirements are not strict, though, the experiments showing only a slow degradation of the results when the distance varies around the optimal position. The method is also very strict in the precautions to compensate for the variations in illumination and to estimate the exact deformation between the ideal pattern and the observed one.

Nevertheless, the pay off of this careful procedure is high. The method delivers a very accurate estimate of the psf, as amply shown in the manifold comparative experiments, with quite stringent accuracy levels (relative error in the order of 2-5%). It remains to wonder why the former methods added regularizing terms or a priori models if these were not needed. Yet, the numerical experiments have confirmed that the inverse estimation problem is indeed ill-posed with slanted edge patterns, which accounts for the necessity of regularization terms for such patterns. Although random noise patterns have been widely used in the past, up to our knowledge no regularization-free sub-pixel psf estimation scheme had been previously proposed. For these previous methods with noise patterns, the lack of a careful correction for all perturbations may explain the need for a regularization or an a priori model. The experiments here have confirmed that for typical DSLR cameras, each color channel is under-sampled with respect to the ideal Nyquist rate given by the psf, by a factor of 2 or even 4. This fact was confirmed, even with DSLR models including an optical anti-aliasing filter on the sensor. This more than justifies a posteriori the need of a sub-pixel estimation procedure. As usual, a locality-accuracy trade-off had to be resolved. The locality of the order of a few hundred pixels can be achieved under common noise conditions.

Of course a wholesome local camera calibration remains a heavy procedure. According to the above setting, some 100 snapshots of the pattern are needed to cover the whole image domain to get an accurate enough psf estimate everywhere. Indeed, the experiments show that this kernel varies significantly, particularly near the image boundaries. The only solution to avoid these many photographs would be to print a very large random pattern covering a whole wall, which is actually quite doable in lab conditions.

Fig. 1 :

 1 Fig. 1: Comparison of the lower bound given by Lemma 3 and the one obtained by solving the KKT conditions for the case a = 0, b = 1. Both bounds are shown in (a), as a function of M . The plot in (b) shows their difference, also as a function of M .

Fig. 2 :

 2 Fig. 2: Comparison of the lower bound given by Lemma 3 and the gamma obtained from the expected S s U operator when using a Bernoulli(0.5) random noise target. Both bounds are shown in (a), as a function of M whereas the Plot (b) shows their difference, also as a function of M .

Fig. 3 :

 3 Fig. 3: Reaching theoretical bounds. A random Bernoulli binary image is used to generate the S s U. We set s = 4 and estimate γ for different observed image sizes (m, n values) (a) and different kernel sizes (r value) (b). The proximity between the obtained γ and the theoretical bound shows the tightness of the derived γ lower bound.

Fig. 4 :

 4 Fig.4: Random Pattern Analysis. Sensitivity of the γ value to the kernel support size (a) and to the t/s zoom factor (b) s = 4. The larger the support of the kernel, the noisier the estimation when the gamma value increases with the kernel support size (a). The zoom factor s/t is closely related to the focal distance and to the distance from the camera to the pattern. For example if the distance to the pattern is too small (small s/t value) the pattern will look like a step-edge pattern because of the zoom-in. The corresponding γ value will be higher than the optimal. On the other hand, if the distance to the pattern is large then the γ value will also be larger than the optimal one, because of the contrast loss due to the zoom-out. In agreement with the theoretical study, the views with zoom factors close to one (i.e. t ≈ s) produce the best γ values.

  Fig 5 shows the proposed random pattern, compared to a pattern designed by Joshi et al. consisting of 120 • arc step edges.Suppose we want to do a s = 4× psf estimation. As shown in the previous section, the pattern should be photographed at such a distance that the pattern covers more or less 256/4 × 256/4 pixels. In practice, this permits a very local psf estimation.Fig.6shows the eigenvalues of the S s U matrix for s = 4, an observed window with size 80 × 80, and varying kernel sizes, for Joshi et al.'s pattern and for the proposed random target. The random pattern produces

Fig. 5 :

 5 Fig. 5: Calibration patterns for local psf estimation

Fig. 7 :Fig. 6 :

 76 Fig.7: Algorithm Description. The captured image is precisely aligned to the analytic pattern through intentionally inserted checkerboard markers. Non-uniform illumination and non-linear camera response function impactcrf are corrected from the captured image allow an artifact-free s× psf estimation.

Fig. 8 :

 8 Fig.8: Synthetic example I. Performance comparison for simulated data. A 4× kernel is estimated using the Joshi et al. algorithm, with varying regularization level, and the proposed approach. The observed window has 110×110 pixels. The top row shows the kernel estimation and the middle row the difference image between the estimation and ground truth for one of the realizations. The proposed method significantly outperforms the Joshi et al. algorithm, achieving a much less noisier estimation as shown by the difference images and by the peak signal to noise ratios. The bottom row shows central horizontal profiles for all the estimated mtfs (0.5 is the Nyquist frequency). Notice that in the Joshi et al. method the estimation is unstable. The estimates show extremely noisy components for frequencies higher than the sampling frequency, when the amount of kernel regularization is too small. On the other hand, if a strong regularization is imposed, the penalization of the kernel gradient adopted by Joshi et al. tends to produce kernels with under-estimated high frequency components. The method proposed here does not rely on a regularization and produces nonetheless noiseless and unbiased results.

Fig. 9 :

 9 Fig. 9: Different apertures. Taken at different apertures, green channel g1, 100 iso, 50mm. All estimated 4× kernels are quite smooth. Fig. (a):The top and bottom rows show respectively the estimated psfs and a few level lines of the corresponding mtfs that prove that the kernels are not exactly axissymmetric. The kernels at apertures f/32 and f/16 are considerably larger than the rest in agreement with diffraction theory. This phenomenon also stands out in the modulus of the estimated psf spectra, which also shows that the psfs/mtfs are not axis symmetric. Figure (b): Vertical cuts of the spectrum modulus. The camera seems to have the sharpest response from apertures f/3.2 to f/12.9. At apertures f/32 and f/16 the camera cuts high frequencies significantly more than the rest, as predicted by diffraction theory. Notice that in all cases, except at aperture f/32, the mtf at the Nyquist frequency (f = 0.5) is significantly greater than zero. Hence, the camera introduces aliasing.

Fig. 10 :

 10 Fig. 10: Real camera example. Taken at f/5.7. An example image, to show how local the psf estimation is (left), and a zoom of the observed window of size 110 × 110 pixels (right).

Fig. 11 :

 11 Fig.11: Diffraction-limited System. Theoretical diffraction mtf for monochromatic green light with circular f/5.7 aperture and the estimation for the green channel at the same aperture. The estimated response for our camera-lens system is under the ideal diffraction-limited response. This can be consequence of the light integration in the sensor array, but also of the optical low pass filter specifically included to avoid aliasing.

Fig. 12 :

 12 Fig.12: otf phase. Estimation done for the green channel g1, 100 iso, 50mm, f/5.7 at the center of the sensor array. The Figure on the top shows the modulus of a horizontal profile of the optical transfer functionotf and its real component. Both curves coincide, implying that the otf is real and thus the psf is symmetric. This is also seen in the bottom figure that shows that the otf phase is 0 or π.

Fig. 13 :

 13 Fig.13: Dependence on the crf correction, for a 4×psf estimation of the green channel, at aperture f/5.7. Figure (a): crf estimates obtained with the estimation embedded in the proposed psf estimation algorithm (crf-psf estimation), and with the one generated independently from a special purpose crf calibration pattern (crf-pattern). Both estimates are very similar and hardly non-linear. Figure (b): vertical profile of mtfs. The estimates from the raw gray values with and without crf estimation/compensation (psf-crf and psf-nocrf, resp.) gave very similar results. After compensation of the gray values using an external special purpose estimation of the crf, the psf estimation procedure (psf-nocrf-eq) also led to very similar results.

Fig. 14 :

 14 Fig.14: Different color channels. 4× psf estimation for the four Bayer pattern channels (two greens, red and blue) from a camera raw output. Top row: psf estimation. Middle row: the corresponding Fourier spectrum modulus. Bottom row: mtf horizontal and vertical profiles. The estimation was performed using the random pattern captured at aperture f/5.7. The red psf is larger than the green and the blue ones. Since the wavelengths associated to red are smaller than the rest, the diffraction components for the red channel will be larger than those for green and blue for the same camera configuration. Also notice the differences between the shape of the red/blue and green psf spectra (bottom row). Red and blue mtf seem to be 45 • rotated with respect to the green ones. This symmetrical behavior is plausible for an L-shaped active zone sensor array.

Fig. 15 :

 15 Fig.15: Different locations. Taken at f/5.7 for one of the green channels. The psfs estimated far from image center are larger and more asymmetrical than the one estimated at the center. This is certainly due to lens aberrations, which are more significant near the image borders.

Fig. 16 :

 16 Fig.16: Stability of the estimation. Average and standard deviation statistics were generated from 13 estimations computed at f/5.7 (for one of the green channels). The small standard deviation in the vertical profile of the Fourier spectrum modulus is shown in (a). The relative mtf sensitivity vs region threshold is shown in (b). We define the relative mtf sensitivity in a region Ω as: s(Ω) = mean(std(mtf))/mean(mtf) where the mean values are computed inside the region Ω. In this case we construct Ω(threshold) = {x : mtf(x) ≥ threshold}. The relative sensitivity in the whole spectrum does not exceed 0.08 and what is more if the mtf values smaller than 5% are not considered, then the relative sensitivity is less than 3%. The small standard deviation and sensitivity demonstrate the algorithm stability.

Fig. 17 :

 17 Fig.17: Non-negative Constraint. This experiment analyzes how the psf estimation changes by not assuming the nonnegative hypothesis. On the left we show both estimations: the no-constrained and the non-negative 4× psf for the green channel, f/5.7 at the center of the image. Since there is no structure in the image produced by subtracting both estimations and since the relation between the energy of the image difference and the energy of the non-negative estimation is 0.001, we can conclude that both estimations are extremely close. This is confirmed by observing in the left figure a horizontal profile of the mtf for both estimations.

Fig. 18 :

 18 Fig. 18: Changing support size. This experiment analyzes how the psf estimation changes with the desired psf support size. Several psf estimations for various kernels support sizes (left). Only the central 9×9 regions are shown. All the estimates are very close, specially 17 × 17, 25 × 25 and 33 × 33.However, the 9 × 9 kernel support seems to be hardly sufficient for correctly modeling the psf, as indicated by the mtf vertical profiles on the right. The proposed algorithm does not appear to be sensitive to this parameter as soon as the kernel support exceeds this minimal size .

Fig. 19 :

 19 Fig.19: Comparison of psf/ mtf estimation methods applied to a real camera. Our implementation of Joshi et al. psf estimation algorithm, the Imatest commercial software and the proposed random patter algorithm. All estimations are done at the center of the image with a camera at aperture f/5.7 for one of the green channels. On the low frequencies all algorithms gave very similar estimations, while on the higher frequencies the Joshi et al. estimation depends on the regularization level. Although we did our best to get a noise free mtf estimation from the Imatest software, the final estimation is quite noisy. The Imatest estimation is done from a slanted-edge image and only gives an estimation for the mtf at the slanted-edge orthogonal direction.

Table 1 :

 1 psf estimation algorithm summary.

Table 2

 2

			. In all cases, the random
	pattern significantly outperforms the Joshi et. al. pat-
	tern. The γ bound value was computed by taking into
	account the effective observed window size, that is, leav-
	ing out the auxiliary region with the checkerboard and
	flat regions.
		9 × 9 17 × 17 25 × 25	33 × 33
	Joshi	99.44 1133.05 6445.87 58419.08
	Random	0.19	0.69	1.54	2.98
	Bound	0.10	0.35	0.70	1.15

Table 2 :

 2 A comparison of pattern realizations through the γ value. The random pattern produces significantly smaller γ values than the slanted-edge Joshi pattern.

Acknowledgements The authors would like to thank Rafael Grompone von Gioi and Saïd Ladjal for fruitful comments and discussions. Research partially funded by the Uruguayan Agency for Research and Innovation (ANII) under grant PR-POS-2008-003, the Centre National d'Etudes Spatiales (R&T), the European Research Council, ECOS-Sud Project number U06E01, and the Office of Naval Research (grant N00014-97-1-0839).