Jean-Michel Morel

Pablo Musé

Andrés Almansa

Resumen de la Tesis

Mauricio Delbracio Bentancor Dos problemas en la formaci ón de imágenes digitales: la estimaci ón de la funci ón de dispersi ón de punto de una cámara fotográfica y la aceleraci ón de renderers estocásticos por filtrado auto-

Introduction

Motivation

Digital images are generated by using physical acquisition devices, such as digital cameras, but also by simulating light propagation through environmental models. This thesis deals with two fundamental problems of digital image generation: the modeling and estimation of the camera's intrinsic blur presented in all optical digital devices and the generation of fast and realistic synthetic images.

Mostly due to advances in technology, pushed by a massive market penetration, digital image quality has significantly improved in the past two decades. This has caused that certain characteristic problems, caused by the nature of the acquisition devices, that were not critical until recently, are now becoming the bottleneck to get further improvements. Blur produced by light diffraction, geometrical distortions caused by the use of low-cost lenses and thermal noise due to electronic circuits are examples of this kind of outcomes.

Image blur can be observed when the camera's focus is not correctly adjusted by the user, when the objects in the scene appear at different depths, or when the relative motion between the camera and the scene is faster than the shutter speed (motion blur). In addition to these sources of blur, even in ideal acquisition conditions, there is a permanent intrinsic physical camera blur due to light diffraction, sensor resolution, lens aberration, and anti-aliasing filters. The first part of this dissertation, addresses the problem of accurately estimating the Point Spread Function (psf), that models the intrinsic camera blur. This function can be locally interpreted as the response of the camera to a point light source.

At the other end, synthesizing high quality realistic images in a reasonable amount of time remains a major challenge in computer graphics. The aim of realistic image synthesis is to generate new images from a complete three-dimensional description of a virtual scene. The scene description should contain at least the geometry, location and properties of objects, the camera viewpoint and a characterization of light sources. The generated picture should be as photorealistic as possible: if the three-dimensional scene is constructed and a photograph is taken from the same camera's point of view, the difference should be negligible. Of course this requires a perfect knowledge of how the light interacts with the environment and extremely accurate material models; oversimplifications must be avoided.

The seminal paper by Kajiya [] presented the rendering equation, an integral equation modeling the steady-state light distribution in a scene. Except for very simple scenes, analytical solutions are impossible to obtain, so most typical approaches are based on Monte Carlo numerical integration techniques. Image pixels are formed by averaging the contribution of stochastic rays cast from a virtual camera through the scene. The principal problem of Monte Carlo ren- ders is that the variance of the estimator decreases linearly with the number of stochastic samples. Thus the root mean square error to an ideal image decreases as the square root of the number of samples. Several hours or even days may be necessary to produce noiseless realistic images. Indeed, at present, the final image quality is indirectly topped by the available production time and computational resources.

This makes this problem interesting not only from an academic point of view. Indeed, the motivation for this problem came through a partnership with the French company e-on software  , specialist in modeling and rendering natural environments. The company's ceo, Nicholas Phelps proposed a concrete well defined problem:

Is there any way of generating realistic synthetic images, more quickly, with the same amount of computational resources and without loosing quality?

The second part of this dissertation analyzes some ideas and proposes a new method to accelerate Monte Carlo renderers, which are the most popular realistic renderers that are currently used.

Part I -Recovering the Camera Point Spread Function

The point spread function (psf) describes the distribution of light in the camera focal plane for a point light source. This function is strictly related to the resolution and blur of an optical device. Most medium to high quality digital cameras (dslrs) acquire images at a spatial rate which is below the Nyquist rate. For this reason only aliased versions of the camera point-spread function can be directly observed. In addition, since the acquisition system is only locally stationary, the psf estimation must be local.

psf estimation methods can be classified as blind or non-blind, depending on whether they use or not snapshots of a specially designed calibration pattern. Blind approaches try to estimate the psf from photographs of an unknown scene. They do assume, however, that the scene involved in the estimation follows some statistical model of sharp images, or includes a significant amount of geometric cues such as sharp edges. Most of these psf estimation approaches attempt to detect edges, which are modeled as pure step-edge functions convolved with the psf kernel [Chalmond ; Luxen and Förstner ; Capel ; Smith ]. In this setting, the estimation is very illposed; to solve the inverse problem, the solution space has to be constrained by considering kernels with a parametric model or with strong regularity assumptions. Therefore, such blind estimation techniques do not normally lead to accurate psf estimates and are only used in image restoration problems, where precision is not the main objective. For this reason, most accurate psf estimation procedures rely on the use of specially designed calibration patterns. A local kernel estimation is performed by comparing the ideal calibration pattern to its photographs.

In this thesis, two fully automatic algorithms for the psf estimation are introduced. The first one performs an estimation from a photograph of a specially designed calibration pattern. As we will show, this non-blind algorithm achieves high accuracy. The second one uses two photographs of a planar scene, taken at two different distances with the same camera configuration. This algorithm lays in an intermediate category: semi-blind. Although, it does not make use of a psf estimation pattern, the photographed scene should be planar and textured to get accurate results.

Pursuit of the Optimal Calibration Pattern

Up to the present, even non-blind subpixel psf estimation methods reported in the literature led to ill-posed inverse problems. The inversion required the imposition of simple psf parametric models, or other regularity or symmetry priors. In Chapter  we show that such a priori assumptions on the psf are actually unnecessary and jeopardize the estimation accuracy.

A mathematical digital image formation model that takes into account blur due to intrinsic and extrinsic phenomena, geometrical lens distortions, sampling and noise is presented in Section .. One of the main problems for the psf estimation, is that digital cameras may capture images that are undersampled according to Shannon theory. The mathematical model, built on Lemma , considers this issue specifically. In this discrete model, images are sampled in a virtual superresolved grid. This high-resolution lattice will be for example, 4 × 4 times the resolution of the physical camera allowing a 4× psf estimation.

Assume we can unveil exactly the latent sharp pattern image that produced the blurry aliased observation. Then, solving for the psf amounts to solve an inverse problem governed by the image formation model. To achieve this several problems need to be addressed:

• how to choose a good psf calibration pattern;

• how to estimate the geometric deformation between the pattern and the acquired image;

• how to estimate the non-uniform illumination;

• how to numerically solve the inverse problem.

If the pattern and the captured image are perfectly registered and its non-uniform illumination compensated, the accuracy of the psf estimation depends on how well we can invert an operator. This operator mainly depends on the calibration pattern. In Section . we present a quality measure (that we call the γ value) of a given pattern view in terms of the well-posedness of the psf estimation problem. Formalizing this well-posedness measure allows us to define an optimal digital pattern for the subpixel psf estimation (Definition ). A theoretical bound on the optimality of a pattern is presented in Proposition . This mathematical bound gives a limit of the performance that we can theoretically achieve. Several numerical experiments conducted in Section . allow us to conclude that in realistic conditions, near-optimal quality measure values are reached with a pattern created from a realization of a random field of independent Bernoulli black or white pixels.. The mathematical modeling of the psf estimation problem together with the near-optimality of the Bernoulli pattern is one of the main contributions of this dissertation.

Is it feasible to construct such an ideal calibration pattern?

Figure . shows the proposed Bernoulli noise pattern, compared to the pattern designed by Joshi et al. [] consisting of 120 • arc step edges. The proposed noise pattern consists of 256 × 256 small black/white squares drawn independently from a Bernoulli equiprobable distribution. This central region is surrounded by checkerboard-like marks and black/white squares for alignment and illumination estimation purposes.The fact that the pattern is black and white avoids to calibrate the printer used for generating the pattern.



Slant-edge pattern by Joshi et al. [] Bernoulli pattern (proposed) The well-posedness measure for the corresponding patterns and the theoretical optimal values are shown in Table .. The Bernoulli noise pattern significantly outperforms the one from Joshi et al. [].

The pattern is designed to cover a small region of 100 × 100 pixels in the acquired image. This permits an estimation which is both local and well-posed. An example of an acquired image is shown in Figure .. This image shows one of the green channels of the Bayer raw camera output.

The proposed non-blind psf estimation algorithm is fully automatic. The captured image is precisely aligned to the analytic pattern by means of the surrounding checkerboard markers. Nonuniform illumination is corrected from the acquired image by using the auxiliary black and white flat regions. Non-linear sensor response is also roughly estimated based on the fact that the central noise part of the pattern should have a perfect average of black and white. Of course, since we work with raw camera output, this response is almost linear.

Once these intermediate steps are performed, the local psf is directly computed by inverting a linear system. Since the psf must be non-negative, as we justify in the image formation model, a numerical constraint is imposed to enforce this non-negativity. Notwithstanding, the experimental section shows that not enforcing the kernel to be non-negative essentially yields the same results. In fact, this serves as a sanity check on the proposed method.

In order to validate the calibration procedure and the quality of the Bernoulli pattern we conducted several synthetic and real camera experiments. A comparison to the commercial software  Imatest [LLC ] and the state-of-the-art method of Joshi et al. is presented in Figure .. This figure shows the horizontal profile of the modulus of the psf Fourier spectrum. The Imatest estimation is performed from a slanted-edge image and only gives an estimate of the point spread function at the direction orthogonal to the slanted-edge. The Joshi et al. method is forced to use a penalty term on the norm of the psf gradient, since the inverse problem using their slant-edge pattern is ill-posed. In the low frequencies Joshi et al. and the proposed method yield very similar results. However, for higher frequencies the result of Joshi et al. is strongly dependent on the regularization level. The Imatest estimate is quite noisy and does not resolve frequencies above twice the sampling rate. The proposed algorithm based on the Bernoulli pattern generates much more information than the typical slanted-edge psf calibration. Joshi et al. [2008] , Imatest [LLC 2010] commercial software and the proposed Bernoulli pattern algorithm. 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. The Imatest software produces a quite noisy estimation.

Figure . shows the result of estimating a 4× psf for one of the green channels at different image locations using the proposed Bernoulli pattern. Kernels closer to image borders are larger and more asymmetrical than the kernel at the image center. This agrees with the expected result, since lenses are designed to minimize aberrations near the optical center. Note that although no regularization is imposed, the resulting psfs are smooth.

Avoiding the Use of a Calibration Pattern

Although very precise, the use of a calibration pattern can be sometimes tedious and impractical: these approaches rely on a careful setup, and the calibration grid has to be properly assembled, whereby a good quality print is essential. Therefore, we explored the feasibility of obtaining accurate psf estimates, while avoiding the explicit use of a calibration pattern.

Chapter  proves that, instead of using a photograph of a known calibration pattern, two photographs of the same scene acquired at different distances with fixed camera configuration are For each pixel at the sensor resolution we estimate 4 × 4 samples of the PSF. The 1× resolution (camera grid) is shown in green dotted lines. Kernels closer to image borders are larger and more asymmetrical than the kernel at the image center. This is mainly due to lens aberration. Although no regularization is imposed, the resulting kernels are smooth.

enough to recover a regularization-free subpixel psf. The mathematical relation between these two fronto-parallel images allows us to introduce the concept of inter-image kernel between them. This kernel can be interpreted as the blur that should be applied to the closest image (followed by the necessary zoom-out) to produce the farthest image (Definition ). An interesting observation, which is proved in Lemma , is that the inter-image kernel k, and the camera psf h are closely related through the following equation

H λ h * k = h,
where H λ is a zoom of the necessary factor to put the two views in the same scale. The derivation of the camera psf from the previous equation is not straightforward. However, as proved in Proposition ,

h = lim n→∞ H λ n-1 k * H λ n-2 k * • • • * H λ k * k.
Thus, it is possible to recover the camera psf from the inter-image kernel k. One of the main contributions of this dissertation is the mathematical development that proves that such an estimation is possible.

In practice, we deal with discrete images, that may suffer from noise and aliasing, so the interimage kernel estimation may be biased. Based on a mathematical analysis, we give some ideas on how to mitigate the impact of these problems and to increase the accuracy of the inter-image kernel estimation. This is done by properly choosing the scene and the distance between the acquired images. Indeed, the inter-image kernel is estimated by solving a least squares problem similar to the one for the pattern based psf estimation. As shown in Chapter  the inverse problem is wellposed as long as the photographed scene presents textured characteristics (similar to white noise). In this setting, the closest image plays a similar role as the one played by the calibration pattern in a traditional non-blind estimation.

We introduce an algorithm that is completely automatic. It gives a subpixel estimate of the camera psf from two photographs of a textured planar object taken at different distances. An Then, a least-squares problem is solved to get the inter-image kernel estimation. Finally, the camera psf is computed from the inter image-kernel through an iterative procedure. Experimental evidence shows the well-posedness of the problem and the convergence of the proposed algorithm to the camera in-focus psf. The mathematical development is made under a technical assumption (e.g., re-focusing does not change the psf) that may not strictly hold in practice. Nevertheless, as shown in Figure . the psf estimated with the Bernoulli pattern and the estimation with the two scaled photographs are convincingly close.

The proposed acquisition procedure is simple and handy in comparison to a non-blind approach. The choice of the photographed scene is important but not critical. For a wide range of everyday textured scenes, the acquired image pairs lead to well posed inversions and accurate results. Figure . gives an idea of the sharpness of the method. In this experiment, the psf at four times the camera resolution of the four color Bayer channels (typical color filter array in raw camera output consists of two green channels, one blue and one red) are estimated. Notice that the red channel psf is wider than the green and the blue one, as expected from the physics of diffraction-limited optical systems, since the wavelengths associated to red light are larger than the rest. The differences between the dominant orientations of the red/blue and green psf spectra can be explained by the sensor shape and layout. In fact, each sensor active zone is usually L-shaped  (see for example Yadid-Pecht []), and the red and blue sensors are rotated 90 • w.r.t. the green ones . These rotations are consistently observed in the psfs and Fourier spectra estimated with the proposed two photographs method. This illustrates the precision of the proposed approach.

 -2 -1 0 1 2 -2 -1 0 1 2 green  -2 -1 0 1 2 -2 -1 0 1 2 blue -2 -1 0 1 2 -2 -1 0 1 2 red -2 -1 0 1 2 -2 -1 0 1 2 Figure 1
.7: Top row: the 4× PSF estimated for the four Bayer channels (two green channels, one red and one blue). Bottom row: their corresponding Fourier spectrum modulus. The red PSF is larger than the blue and the green ones. This is consistent with the diffraction phenomenon: the red wavelengths are larger than the rest, thus its diffraction kernel is wider. Also note the differences between the shape of the red/blue and green PSF spectra (bottom row). Red and blue MTFs are rotated 90 • with respect to the green ones. This symmetric behavior is consistent with the layout of L-shaped sensors Yadid-Pecht [2000].

Part II -Accelerating Realistic Image Synthesis

In order to synthesize an image with global illumination, a radiance value must be assigned to each pixel in the image. Path tracing (and more generally ray-tracing) is a popular technique for  resolving the rendering equation ruling the steady state equilibrium of light in a scene. In a raytracing scenario, this value is computed as a weighted average of radiance values incident on the image plane, along light rays coming from the light sources, bouncing in the scene, passing trough the pixel, and pointing to the virtual camera.

Unfortunately, only a finite number of rays can be cast, so the radiance value is computed only approximately. To avoid artifacts, rays are cast randomly. Mathematically, this is equivalent to solving the rendering equation through a Monte Carlo numerical integration procedure. The main problem of Monte Carlo rendering is that the variance of the estimator converges only linearly with the number of random samples. An example of an scene rendered with a varying number of rays per pixel is shown in There are mainly two approaches to accelerate the convergence of Monte Carlo rendering to obtain good quality images. One of these approaches is adaptive sampling. This class of algorithms locally adapt the number of samples cast per pixel. The idea is to increase the number of rays in complex parts of the scene while maintaining a reduced number in simple parts, such as flat regions. Complex textures or defocused zones are typical elements that require large amounts of rays to be properly rendered. Hachisuka et al. [] proposed to adaptively distribute a set of samples in the full, multidimensional sampling domain where the rendering equation is computed. However, as more Monte Carlo effects are considered (e.g. depth-of-field, motion blur, area lighting) the dimension of this space will be larger and thus will suffer from the curse of dimensionality. One of the most significant adaptive sampling algorithms is certainly the Adaptive Wavelet Rendering by Overbeck et al. []. This method adaptively distributes Monte Carlo samples in the screen space to reduce the variance of a wavelet basis scale coefficients. Then, the image is reconstructed from these non-uniformly distributed samples by using a suitable wavelet approximation.

The other approach is adaptive filtering. In this family of algorithms, the existing set of samples are combined to produce a better estimator of the pixel color using ray information in a pixel and in its neighbors. Adaptive filtering may take place at sample level (i.e., primarily filtering the ray colors) or at pixel level (i.e., primarily filtering pixel color values). The simplest adaptive filters act at pixel level, like any filter used in classical image processing [Jensen and Christensen ; The majority of these methods can actually be written as generalized versions of the bilateral



filter (or the sigma-filter [Lee ]) applying a weighted average of the samples (resp. of the pixels) in a neighborhood. This general bilateral filter obeys the law of joint destiny (Gesetz des gemeinsamen Schicksals) introduced in Gestalt psychology by Wertheimer []. This law states that similar pixels are grouped by our perception. Its obvious generalization in image processing is to say that similar pixels must be denoised jointly, being different samples of the same model. The law of common destiny is implicitly used by the sigma-filter and by the nl-means algorithm [Buades et al. ]. The main disadvantage of traditional bilateral filters is that by comparing noisy pixel values, they cannot easily distinguish noise from intrinsic pixel variability. Thus, the clustering of similar pixels is potentially subject to errors and the filtering will result in a significantly biased image.

As we will present in Chapter , in computer graphics, the statistics of ray samples permits to identify much more rigorously than in classic image processing the pixels sharing the same model. Indeed, all ray samples hitting a given pixel and its neighbors can be used for that purpose.

Similar pixels can be detected by comparing their empirical ray color distributions using an adequate histogram distance. Since the order in which the samples are calculated is irrelevant, the sample color empirical distribution appears as a natural and complete descriptor of the compared sets. Figure . shows a small region of a Monte Carlo rendered image where two pixels are singled out. Although both pixels have different colors, their color distributions are strikingly similar and can be fused. The difference in the pixel colors may be the consequence of the presence of a single very bright ray sample in one of the distributions. By comparing the ray color distributions, it is nevertheless possible to conclude that both pixels are from the same "nature", while this conclusion could not be reached by comparing the pixel values.

Figure 1.9: Monte Carlo rendered pixels can be grouped very efficiently by comparing their ray color distribution. Left: a crop of a Monte Carlo rendered image where two pixels with different colors are singled out. The difference in color is due to a poor estimate from a low number of rays cast at each pixel. Right: the color sample distributions of each pixel. The color sample distribution is represented in the RGB color box, where the color of each of the rays cast at a pixel is one point. The color distributions are strikingly similar and can be fused, which is the principle of the proposed algorithm.

We propose a simple but powerful filtering algorithm that uses exclusively the colors and positions of the cast ray samples. It can be thus coupled with any Monte Carlo sampler keeping a record of rays. The algorithm does not assume any noise model. It generalizes the nl-means denoiser, and shares with it an artifact-free record. The cornerstone of the proposed algorithm is to find and average the most similar patches by comparing the ray color sample distributions of each of its pixels. Inspired by this concept we name the algorithm ray histogram fusion. Figure . shows a running example.

In a pure Monte Carlo rendering the estimation error presents white noise characteristics.



This means that all frequencies are equally contaminated by noise. To this purpose, we introduce a multi-scale implementation that sequentially decomposes the input noisy image at each scale, filters each scale and reconstructs the multi-scale filtered image.

As shown in the experimental section of Chapter , the proposed filter is consistent. As the number of samples increases, more evidence is required to average two pixels. In the limit two pixels will be averaged only if their color distributions are the same. Therefore, in practice, as the number of samples grows the method converges to the expected solution. The acceleration factor depends on the degree of self similarity of the scene, which fortunately is usually high [Lebrun et al. ]. The algorithm provides a psnr gain of  to  decibels, or equivalently accelerates the rendering process by using  to  times fewer samples without observable bias. It is immediately extendable to synthetic movies. Being based on the ray color values only, it can be combined with all rendering effects.

Boosting Monte Carlo renderers

We propose a new multi-scale filter to be used as an acceleration of Monte Carlo renderers. Each image pixel is represented by the colors of the rays that reach its surface. A robust histogram distance compares the empirical color distribution associated with each pixel at each scale and decides whether two pixels can share their rays. This simple and easily reproducible algorithm provides a significant gain in psnr, or equivalently accelerates the rendering process by using fewer samples without introducing observable bias. The algorithm is universal in the sense that it can be combined with any rendering effect. It is consistent and does not assume a particular noise model and is immediately extendable to synthetic movies.

List of Publications

The work in this thesis has led to the publication of the following articles: 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 subpixel resolution, by a numerical method. Since the acquisition system is only locally stationary, this psf estimation must be local. This chapter presents a theoretical study proving that the subpixel 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(.) random noise. The physical realization of this psf estimation method is demonstrated in many comparative experiments. We use an algorithm to accurately estimate the pattern position and its illumination conditions. Once this accurate registration is obtained, the local psf can be directly computed by inverting a well conditioned linear system. The psf estimates reach stringent accuracy levels with a relative error in the order of % to %. To the best of our knowledge, such a regularization-free and model-free subpixel psf estimation scheme is the first of its kind.

• M.

. Introduction

Image blur can be observed when the camera focus 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 function -psf, 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 superresolution, 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 subpixel estimation of the psf is therefore crucial to evaluate the image quality in terms of a trade-off between sharpness and aliasing effects.

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

Existing psf estimation methods can be classified as blind or non-blind, parametric or nonparametric. 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 [Chalmond ; Luxen and 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 slanted-edge patterns to random noise images. The subpixel psf estimation problem is generally treated as ill-posed. Most non-blind methods 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 subpixel-shifted versions of the spot image would be needed. Bar or sine patterns can also help sample the mtf, but only up to the Nyquist frequency.

The iso  standard [ISO ] 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 [Reichenbach et al. ], 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. Zandhuis et al. [] propose a slanted-edge non-parametric subpixel psf estimation method that admits geometrical distortions. A parametric and non-parametric edge spread function estimation procedure is proposed by Claxton and Staunton []. 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.

Before this thesis, the recent method by Joshi et al. [] arguably represented the state-ofthe-art of slanted-edge methods. It proposes a flexible blind and non-blind non-parametric local psf estimation algorithm. Its approach is based on the ability to detect edges with subpixel accuracy. In order to get a precise local psf a specially designed pattern formed by -degrees-arcstep-edges is used. The method directly solves the de-convolution and superresolution problem  for a bi-dimensional subpixel psf. To reach a subpixel 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 ., 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 subpixel 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 by Daniels et al. []; Levy et al. []; Backman et al. [, ]. 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 subpixel accuracy, and aliasing effects are not taken into account. Brauers et al. [] also used a random noise pattern, but in a completely different approach. The acquired image is registered to match the pattern. Then, by doing de-convolution with the almost flat spectrum noise pattern, this method succeeds in characterizing locally the psf. However, the method assumes that the camera over-samples 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 . 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), subpixel (sp). The main systematic perturbations in imaging (optical distortion (D), 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 was the Joshi et al. method, but it includes a regularization which will be shown fatal to the high frequency kernel content. The method proposed here has all "good" features. Indeed, it will be shown 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 pattern, 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 nonuniform 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 subpixel 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 chapter is organized as follows. Section . describes the general mathematical digital camera model used for psf estimation method. Section . proposes a mathematical theory of

Brauers et al. [] NB R K 1-2D-L G Proposed NB R - 1-2D-L-sp DIG



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 . describes all the steps of the proposed psf estimation protocol. In Section . experimental results generated with both simulated and real camera data are presented, cross-validated, and compared with the results of state of the art previous methods. Section . is a final discussion.

. 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 three-dimensional 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, optical anti-aliasing filters, 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

v = S 1 g u • H) * h ex • F * h + n, (M ′ )
where F (•) is the geometric distortion operator, h is the convolution kernel due to all intrinsic psflike effects, h ex is the convolution due to extrinsic blurring effects that occur outside the camera (like motion blur and atmospheric turbulence), and g(•) is a monotone non-decreasing function that describes the non-linear sensor response (camera response function -crf). The operator S 1 is the bi-dimensional ideal sampling operator due to the sensor array, and n(u) models the sensor noise. Physical models of digital camera sensors, both for ccd and cmos sensors, suggest that the readout noise n(u) is a mixture of luminance independent (Gaussian, thermal) noise, and luminance dependent (Poisson or photon counting) noise [Healey and Kondepudy ; Tian et al. ; Marion ]. In fact, the noise can be modeled as white Gaussian noise with luminancedependent variance. For the purposes of this study, however, precise statistics of noise are not critical, and only the global snr is significant, so we shall stick to the more traditional and simpler model of image independent white Gaussian noise.

Furthermore, as stated in the introduction, we will only deal here with intrinsic psf-like effects at the camera focal plane. Therefore we assume and that the experimental setup is capable of avoiding motion blur, atmospheric turbulence and out-of-focus blur as much as possible, thus permitting to neglect the effect of h ex . Strictly speaking out-of-focus blur should be included in the intrinsic camera blur kernel h, but since it can become negligible under controlled experimental setups, we chose to exclude this kind of psf-like effect from our study.



As a result of the previous discussion we shall simplify model (M ′ ) and use

v = S 1 g u • F • H * h + n, (M ′′ )
where h models the camera-intrinsic psf-like blurring effects, and n is an image-independent Gaussian white noise.

The blur kernel h is space variant, but it varies 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 implicit in the combined deformation operator F •H can be considered as a whole: there is no need to estimate separately the projective and non-projective parts. We shall therefore denote by D the whole geometric transformation, and the image formation model becomes

v = S 1 g (u D * h) + n, (M ′′′ )
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

v = S 1 (u D * h) + n, (M)
and explain later on how g 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 [Goodman ]. As a trivial consequence the psf diffraction kernel is always non-negative.

In the reasonable, though inexact hypothesis that the aperture is circular, the diffraction kernel writes

h diff (x) = A λf 2 • 2J 1 (r) r 2 ,
with r = πD|x| λf and A = π(D/2) 2 . The function J 1 (r) = 1 π π 0 cos(θr sin θ)dθ is the Bessel function of the first kind and order one, f is the lens focal length, D the aperture diameter, and λ the wavelength. In the case of a circular aperture, the cutoff frequency of a diffraction limited system is ρ c = 2πD λf . This frequency depends only on the so called F-number = f D and the wavelength λ.

Optical aberrations degrade this ideal system where only diffraction is considered, producing larger kernels [Williams and Becklund ]. In addition, optical anti-aliasing filters -olpf 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 [Zhao et al. ].



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. Yadid-Pecht [] 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. Most digital cameras have only one sensor array. In order to acquire color information, each photo-sensor is filtered to capture only wavelengths of a particular band for the red, green, or blue channels which is done by a color filter array (cfa). The most popular cfa is the Bayer filter mosaic, which covers the sensor plane with 50% of green filters and 25% of blue and red filters respectively (see Figure .). The image formed by the data as it comes directly off the sensor array is called raw image. In a typical configuration of f /D = 5.6 the diffraction cutoff frequency for the green light (λ = 530 nm) is ρ c = 2π 530×10 -9 ×5.6 . Hence, to avoid aliasing the inter pixel distance should be at least δ c = π ρc = 530×10 -9 ×5.6

2

. Assume we have a digital camera with a Bayer sensor of size 22.2 mm × 14.8 mm with a resolution of 3888 × 2592 pixels. This leads to an inter pixel distance δ s = 22.2×10 -3 3888/2 . Thus, in order to avoid aliasing in this diffraction limited system we would need to get samples at s = δ s /δ c ≈ 7.7× the camera resolution. In practice the psf cutoff frequency will be much smaller due to the anti-aliasing filter and the light integration in the sensor.

. 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 re-sampling 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. ), 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 as a white noise image is chosen as the calibration pattern.

.. Inverse problem statement in terms of the re-sampling operator and the calibration pattern

In the following, F denotes the Fourier transform and f the Fourier transform of a function f . The s-Shannon-Whittaker interpolator defined as

I s u(x) = k u(k) sinc(s -1 x -k) is denoted by I s , S s is the s-over-sampling operator S s u(k) = u(s -1 k
) and lpf w is the frequency cut-off low pass filter that cuts the spectrum of a signal to [-wπ, wπ] 2 . Suppose that h is band-limited within supp( ĥ) = [-δπ, δπ] 2 . 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.

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

u * h = I s (ũ * h) ,
where h = S s h and ũ = S s lpf s u. -sπ,sπ] 2 . This implies that u * h = ũ * h. Indeed,

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

ũ * h = I s S s (ũ) * I s S s (h) = I s (ũ) * I s (h) = I s (ũ * h).
Remark . Note that u does not need to be band-limited, only h. Notwithstanding, if we can find the spectral cut-off ũ 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--sub-sampling operator

S s = S 1 I s .
It follows from Lemma  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 over-sampling 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 (.)
Here, ũD is 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 pattern;

• to estimate the function g(•), the non-linear sensor 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 (.) can be rewritten in matrix form,

arg min h S s C[ũ D ]h -v 2 2 , (P)
where C[ũ D ] 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 C[ũ D ] is of size M × N , with N = r × r.
The solution of Problem (P) is easily obtained using a least squares estimation procedure, and is given by

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

(S s C[ũ D ])h + n = v,



the estimation error is given by n e = (S s C[ũ D ]) + 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 C[ũ D ]) + n 2 2 = E    M j=1 N i=1 (S s C[ũ D ]) + ij n i 2    = M j=1 N i=1 N k=1 (S s C[ũ D ]) + ij (S s C[ũ D ]) + 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 C[ũ D ]) + ij 2 σ 2 n = (S s C[ũ D ]) + 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 C[ũ D ] are non zero, the singular values of (S s C[ũ D ]) + are the inverses of the singular values of (S s C[ũ D ]). 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 C[ũ D ]. Then (S s C[ũ D ]) + 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 C[ũ D ]) := N i=1 σ -2 i
It should be pointed out that γ depends on the rate s and on the samples ũD . The superresolution 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× subpixel 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  (Optimal digital pattern). Given a kernel support N = r × r and a window observation size M = m × n, the optimal pattern for the s× subpixel psf estimation is the digital calibration pattern u * such that u * = arg min

a≤u ij ≤b γ(S s C[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 C[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 . 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. (.)
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 . (.)



Thus, by (.) and (.), 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  (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 , for any matrix Φ with entries in [a, b], and in particular for the ones that attain

γ * = 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 ∈[a,b] γ(Φ) ≥ 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 . 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 . 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 . Notice also that in Lemma  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 handy than the bound that was obtained in Lemma , and is worthless since both bounds are extremely close, as shown in Proposition  (Non-asymptotic bound for optimal patterns).

Let u = {u ij } be a ms × ns digital image with all its values in [a, b]. Let S s C[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≤u ij ≤b γ(S s C[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 , since the operator

S s C[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(.) stationary noise. It will be shown that this calibration pattern is so close to the γ(S s C[u]) lower bound given by Lemma  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(.) 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 C[u]) t (S s C[u]) = M m 2 u  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  t + σ 2 u I has only two different eigenvalues:

σ 1 = M (N 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 C[u] when using a Bernoulli(.) 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 . Indeed, for

M ≥ N ≫ 1 we have γ * -γ ≈ 4 M (b-a) 2 . This small difference is illustrated in Figure . for the particular case a = 0, b = 1. Notice also that since E (S s C[u]) t (S s C[u]) = lim M →∞ (S s C[u]) t (S s C[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(.) 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(.), and evaluated γ(S s C[u]) for fixed down-sampling rate s = 4. Figure . shows that the obtained γ is very close to the non-asymptotic lower bound (Lemma ), 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 C[ũ D ] 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 0 ≤ u ij ≤ 1.

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 is almost perfectly satisfied if the views of the pattern are 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 setup 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 random variables u ∈ {0, 1} with Bernoulli(.) distribution. Yet, while this perfect geometric situation is unattainable, the experiments show that γ stays close to its optimal value when s/t is between . and , as it is shown in in

.. 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 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 zoomin. 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. 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 pattern, and to correct non-uniform illumination. 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.

. The Complete PSF Estimation Procedure

In this section we describe the steps that lead to a local subpixel psf estimates. The complete chain is summarized in the block diagram of

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. [Harris and Stephens ], [Cheng et al. ], [Lucchese and Mitra ]), 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 [], that allows us to iteratively refine the detected corner positions to reach subpixel accuracy.

Geometric transform estimation

The estimation of the psf does not require a decomposition of the distortion into its homography and non-homography parts, as it is done in classical geometric camera calibration [Zhang ], 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-plate splines were originally conceived as an exact interpolation method [Bookstein ] they can be easily extended to the approximation problem [Sprengel et al. ]. 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. Accurate geometric transform estimation is essential for good performance. Although there are no "control points" inside the random pattern (only on borders as depicted in Fig. .(a)), as the pattern is designed to cover about 100 × 100 pixels in an image, the local geometrical distortion inside such a small region will be practically affine. Hence, by using thin-plate splines we can achieve the necessary registration accuracy.

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. Since we are working with the raw data and out of the saturation region of the sensors, the sensor response should be almost linear. Hence, for simplicity we model the camera response function as a polynomial of order no larger than . 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 ( and  respectively, after normalization). Consequently, since the psf has unit area, the mean gray value within the observed image should be .. The solution is defined as a parabolic function u → αu 2 + (1α)u where α is chosen so that 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 band-limited in [sπ, sπ]. (Remember that the camera resolution is 1×. Thus the digital pattern has an s× over-sampling). The procedure is:

. 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;

. Frequencies higher than sπ are cut off from the digital pattern u to get ũ;

. 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 non-negative kernels. Thus the estimated psf must be non-negative. We can therefore constrain the solution to be nonnegative, thus reducing the space of solutions. Section ., Fig.

. shows that not imposing this non-negativity 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 to simply solve a least squares 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× subpixel resolution. Also suppose that the estimated support of the psf is inside a r × r image. The matrix S s C[ũ D ] 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

arg min h S s C[ũ D ]h -v 2 subject to h i ≥ 0, i = 1, . . . , r 2 . (P ′ )
Problem (P ′ ) can be solved using standard convex optimization solvers such as cvx [Grant and Boyd ]. A simpler Newton interior point algorithm proposed in [Portugal et al. ] was used and always converged rapidly.

. Experimental Results

This section is dedicated to the evaluation of the proposed non-blind subpixel psf estimation method, and to the comparison of its performance with two state of the art proposed approaches. A complete algorithmic description, an online demo facility and a source code can be found in our IPOL publication Delbracio et al. [b].

We selected a method recently reported in the literature by Joshi et al. [], and a mtf commercial software, Imatest [LLC ]. 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 down-sampled 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 . 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 outperforms Joshi et al.'s algorithm, achieving a much less noisy 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.

.. 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 D camera provided with a Tamron AF -mm F/. XR Di-II lens was used. The focal length was fixed at . 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. . shows the diffraction-limited mtf for a circular f/. aperture and green monochromatic light (See the end of Sec. ..). The estimated response for our camera-lens system at aperture f/. 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.

Different apertures

Estimation of camera response function

This experiment evaluates the impact of the nonlinearity 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.

Figure .(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. 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. 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 axis-symmetric. 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. .12: Diffraction-limited System. Theoretical diffraction MTF for monochromatic green light with circular f/5.6 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.



-2 -1 0 1 2 -2 -1 0 1 2 f/. -2 -1 0 1 2 -2 -1 0 1 2 f/. -2 -1 0 1 2 -2 -1 0 1 2 f/. -2 -1 0 1 2 -2 -1 0 1 2 f/. -2 -1 0 1 2 -2 -1 0 1 2 f/. (





-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 module real-part -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -4 -2 0 2 4 Figure 2
.13: OTF phase. Estimation done for the green channel g1, f/5.6 at the center of the sensor array. The Figure on the top shows the modulus of a horizontal profile of the optical transfer function -OTF 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 π.



• estimation without any crf correction from the raw observed values (psf-nocrf).

• 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 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.

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/.. The results for the 4× psf estimation located in the image center are shown in Fig. .. 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 accept 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 . displays the 4× psf estimates for one of the green channels, at different image locations, for f/.. 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/.), from similar viewpoints. The estimation was performed using the random pattern captured at aperture f/5.6. The red PSF is larger than the green and 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 note 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.



-2 -1 0 1 2 -2 -1 0 1 2 red -2 -1 0 1 2 -2 -1 0 1 2 blue -2 -1 0 1 2 -2 -1 0 1 2 green  -2 -1 0 1 2 -2 -1 0 1 2 green  (a) -2 -1.5 -1 -0.5 0 0.5 1

 left center right

Figure 2.16: Different locations. Taken at f/5.6 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.



For each acquisition, the 4× psf of one of the green channels at the image center was estimated. Figure . 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 Figure 2.17: Stability of the estimation. Average and standard deviation statistics were generated from 13 estimations computed at f/5.6 (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 method stability.

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 . . 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. Figure . 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 slantededge mtf calibration.

Comparison of several methods

. Discussion

The work presented in this chapter is an attempt to define an optimal non-blind subpixel psf estimation method from a single aliased image. The method is successful, but its setup is tight. The pattern must be large enough (some cm 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 well-conditioning 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 Joshi et al. [2008] 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.6 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.



estimate of the psf, as amply shown in the various comparative experiments, with quite stringent accuracy levels (relative error in the order of % to %). 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 regularizationfree subpixel 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  or even . 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 subpixel 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  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. A possible solution to avoid these many photographs would be to print a very large random pattern covering a whole wall, that would cover the whole visual field of the camera. While this is not easy to implement, it is indeed doable in lab conditions.



PSF Estimation from Two Photographs at Different Distances

In most digital cameras, and even in high-end digital single lens reflex cameras, the acquired images are sampled at rates below the Nyquist critical rate, causing aliasing effects. This chapter introduces an algorithm for the subpixel estimation of the point spread function of a digital camera from aliased photographs. The numerical procedure simply uses two fronto-parallel photographs of any planar textured scene at different distances. The mathematical theory developed herein proves that the camera psf can be derived from these two images, under reasonable conditions. Mathematical proofs supplemented by experimental evidence shows the well-posedness of the problem and the convergence of the proposed algorithm to the camera in-focus psf.

An experimental comparison of the resulting psf estimates shows that the proposed algorithm reaches the accuracy levels of the best non-blind state-of-the-art methods.

. Introduction

Light diffraction, lens aberrations, sensor averaging and anti-aliasing filters are some of the inherent camera factors that unavoidably introduce blur in photographs. The blur that results from the combination of all these factors can be modeled locally as a convolution kernel known as point spread function (psf), which corresponds to the space variant impulse response of the whole camera, including the sensor, before the final sampling.

The area enclosed by the first zero crossing of the psf, usually called Airy pattern, is arguably the most reasonable characterization of the optical system resolution. Top camera/lens manufacturers use charts based on the psf Fourier spectrum modulus (the modulated transfer function, mtf) to describe their products. But accurate knowledge of the psf is not limited to quality assessment of optical devices, and it proves to be extremely useful or even necessary for several image processing tasks such as deblurring [Ng et In most typical digital cameras, both compact and high-end dslrs, images are sampled at frequencies below the Nyquist critical rate. Consequently, only aliased versions of the camera psf can be directly observed. Yet, to fully characterize the psf, it is necessary to recover it at a subpixel resolution.

psf estimation methods can be classified as blind or non-blind, depending on whether they use or not snapshots of a specially designed calibration pattern. Blind approaches try to estimate the psf from photographs of an unknown scene. They do assume, however, that the scene involved  in the estimation follows some statistical model of sharp images or include a significant amount of geometric cues such as sharp edges. Most of these psf approaches attempt to detect edges, which are modeled as pure step-edge functions convolved with the psf kernel [Chalmond ; Luxen and Förstner ; Capel ; Smith ]. In this setting, the estimation is very ill-posed; to solve the inverse problem, the solution space has to be constrained by considering kernels with a parametric model or with strong regularity assumptions. Therefore, such blind estimation techniques do not lead to accurate psf estimates and are consequently constrained to image restoration problems, where precision is not the main objective. For this reason, accurate psf estimation procedures rely on the use of specially designed calibration patterns. A local kernel estimation is performed by comparing the ideal calibration pattern to its photographs.

As Even non-blind subpixel psf estimation methods reported in the literature led to ill-posed inverse problems. The inversion required the imposition of simple psf parametric models or other regularity or symmetry priors. In Chapter  we have shown that such a priori assumptions on the psf are actually unnecessary and jeopardize the estimation accuracy. More precisely, by carefully modeling the image acquisition system, a calibration pattern made of a white Bernoulli noise realization is nearly optimal in terms of well-conditioning of the problem. This procedure leads to very accurate regularization-free subpixel psf estimation.

The purpose of the present work is to explore the feasibility of obtaining accurate psf estimates, while avoiding the explicit use of a calibration pattern. The motivation comes from the fact that, although very precise, the use of a calibration pattern can be sometimes tedious and impractical: these approaches rely on a careful setup, and the calibration grid has to be properly assembled, whereby a good quality print is essential.

We show that, instead of using a photograph of a known calibration pattern, two photographs of the same scene acquired at different distances with fixed camera configuration are enough to recover a regularization-free subpixel psf. The proposed acquisition procedure is simple and handy in comparison to a non-blind approach. Experimental evidence will show that the resulting estimates do not exhibit any significant accuracy loss compared to their best non-blind competitors. The choice of the photographed scene is important but not critical. For a wide range of everyday textured scenes, the acquired image pairs lead to well posed inversions and highly accurate results. The proposed method can be used with the previously introduced Bernoulli pattern as well, the difference being that in that case the quality requirement for the printed pattern are no longer stringent.

This chapter is written with a dual public in mind: mathematicians and/or image processing specialists. We have tried to define accurately all mathematical objects necessary to deal rigorously with image formation. An accurate formalism is needed to justify the somewhat intricate interlacement of sampling and convolution operations. This forces one to check on the compatibility of all function or distribution spaces to which the objects belong and to verify that the formulas are mathematically consistent. Nevertheless, the application-oriented reader can skip the proofs and the functional space details at a first reading, and simply focus on the standard image processing formalism and algorithms. Most proofs are placed at the end of the chapter. A glossary is appended to display all notation in a single place.

The chapter is organized as follows: Section . presents a mathematical model of the digital image acquisition system. This model is used in Section ., where it is shown that the camera psf  can be recovered from a pair of unknown scaled images. We define the notion of blur between such a pair of images, and we propose a method to perform its estimation. Then we prove that the camera psf can be recovered from this inter-image blur. Section . presents an algorithm that implements the complete psf estimation procedure described in Section .. In Section . we discuss a series of experiments on real and simulated images. Finally, Section . closes with a brief recapitulation and conclusions. The details of the adopted notation and the complete mathematical proofs are presented in Appendices .A and .B respectively.

. Image Formation Model

.. Generalized digital pin-hole camera

An accurate estimation of the psf requires a proper modeling of the digital image formation process. The geometric component of this process is most often modeled in computer vision by a pin-hole camera. An ideal pin-hole camera with focal length f , shooting at a planar scene u from a distance d and at fronto-parallel pose, will produce an image w(x) = u(λx) which is just an homothecy of scale factor λ = d f of the original planar scene u. If the pose is not perfectly fronto-parallel or the pin-hole camera presents non-canonical internal calibration parameters, w and u are related by a planar homography D, i.e., w = u • D. In a more accurate camera model the distortion D takes the form of a more general (but regular) diffeomorphism. This is required when the scene is a regular close-to-planar surface (as it is assumed here) or when the geometric distortion due to the optical system is taken into account as suggested in Chapter .

For the purpose of psf estimation this simple model needs to be augmented with an accurate radiometric component, comprising at least the following elements.

Blurring

The psf kernel h models blur due to intrinsic camera characteristics, such as diffraction when light goes through a finite aperture, light averaging within the sensor and lens aberration. Other blur sources such as motion, atmospheric turbulence or defocus blur, that may change from one snapshot to another, will be minimized by the experimental procedure, and it is not the goal of the present work to estimate them. Another implicit assumption that is usually made is that as long as the camera is in focus, the psf is independent of the focus position, i.e., the relative distance between the sensor array and the lens system. Therefore in focus images captured with the same camera configuration are affected by the same psf.

The diffraction kernel is determined by the shape and size of the aperture, the focal length, and the wavelength of the considered monochromatic light. Under the Fraunhofer far-field approximation, for incoherent light this kernel is the squared Fourier transform modulus of the camera's aperture indicator function [Goodman ]. It follows that the psf diffraction kernel is always non-negative and band-limited.

Besides the kernel due to diffraction, other sources of blur inherent to the optical system are present in real cameras. These are mainly optical aberrations, and anti-aliasing filters (which reduce aliasing but do not completely cancel it) introduced in the system prior to sampling [Williams and Becklund ; Zhao et al. ]. The sampling process also introduces blur. Indeed, each photo-sensor in the rectangular sampling grid integrates the light arriving at a particular exposure time. This corresponds to a convolution with the indicator function of the photo-sensor active  area. To sum up, the unknown psf results basically from the convolution of three non-negative kernels (diffraction, aberrations and anti-aliasing filters, and sensor averaging), one of them being band-limited. No parametrical model on the psf will be adopted here. Nonetheless the physical modeling justifies our assumption that the psf is band-limited and non-negative.

Sampling

We model the continuous to digital conversion at the image plane by the introduction of an ideal sampling operator S 1 and additive noise n due to measurement uncertainties. Physical models of digital camera sensors, both for ccd and cmos sensors, suggest that the readout noise n is a mixture of luminance independent (Gaussian, thermal) noise, and luminance dependent (Poisson or photon counting) noise [Healey and Kondepudy ; Tian et al. ; Marion ]. A usual simplification of this model, which we follow here, assumes the noise is image independent, white and Gaussian, with constant variance.

The whole image formation process can then be summarized in a single equation:

ṽ = g (S 1 ((u • D) * h)) + n,
where g(•) is a monotone non-decreasing function that describes the non-linear sensor response. If the camera is working outside the saturation zone, in raw images this response can be reasonably assumed to be linear (see Figure . in the Experimental section of Chapter ). This boils down to a rescaling of the dynamics of u and therefore disappears with out loss of generality from the model. Hence, in what follows, the image formation model will be

ṽ = S 1 ((u • D) * h) + n.
(M)

.. Inverse problem statement in terms of digital sequences

Since in practice our data consist exclusively of discrete sequences (or digital images), the image formation model will be rewritten in terms of discrete sequences. This requires the introduction of additional notation, summarized in Table . (a more precise definition of each term is presented in Appendix .A). It would be cumbersome to verify systematically all regularity requirements on all functions and distributions needed in the proofs. Thus, all necessary results are given in a precise form in the appendices. They will be invoked in the proofs, and the reader is invited to check that their use was licit. Suppose that the psf h is s-band-limited, that is, supp( ĥ) = [-sπ, sπ] 2 . Then, if sampled at a rate s, the Nyquist sampling theorem guarantees perfect reconstruction of h from its samples h = S 1 H 1 s h. We are actually interested in the case s > 1, usual for digital cameras. This means that the images obtained from (M) may be subject to aliasing.

Proposition  (discrete camera model). Let u ∈ BL 2 0 and h ∈ L 1 ∩ BL 2 0 , band-limited in [-sπ, sπ] 2 . Then S 1 (u * h) = S s (ū * h), (.)
where we have called

ū = S 1 W 1 H 1 s u and h = S 1 H 1 s h.

 u, v

Images defined on continuous domain x ∈ R 2 u, v

Digital images are sampled on a discrete grid k ∈ Z 2 F Fourier transform f Fourier transform of a function f I 1

Shannon-Whittaker interpolator:

I 1 u(x) = k u(k)sinc(x -k) S 1 1-sampling operator: u(k) = (S 1 u)(k) = u(k) W w
Ideal low-pass filter cuts the spectrum of continuous signals to [-wπ, wπ] 2

S s

The s-to--resampling operator

S s = S 1 H s I 1 H λ Continuous homothecy: H λ u(x, y) = λ 2 u(λx, λy). (λ < 1 dilation) H α
Digital Nyquist homothecy operator of parameter α:

H α u := S 1 W 1 H α I 1 u C[u] Linear map associated to the convolution with a digital image u L * Adjoint of a linear operator L L + Pseudo-inverse L + := (L * L) -1 L * of a linear operator L L 1 Integrable functions on R 2 (L 1 (R 2 )) L 2 Square integrable functions (L 2 (R 2 )) BL 2 L 2 functions, band-limited in [-π, π] 2 BL 2 0
L 2 functions with compactly supported Fourier transform Proof. We first derive the expression and then justify the application of each result. The set of used properties are detailed in the Appendices .A and .B.

S 1 (u * h) = S 1 H s H 1 s (u * h) (.) = S 1 H s H 1 s (W s u * h) (.) = S 1 H s (H 1 s W s u * H 1 s h) (.) = S 1 H s (W 1 H 1 s u * H 1 s h) (.) = S 1 H s I 1 S 1 (W 1 H 1 s u * H 1 s h) (.) = S 1 H s I 1 (S 1 W 1 H 1 s u * S 1 H 1 s h) def = S 1 H s I 1 (ū * h) def = S s (ū * h).
First note that as u ∈ BL 2 0 and h ∈ L 1 are band-limited in [-sπ, sπ] 2 , we can apply (.) and (.) directly. As W 1 u is in BL 2 we can apply (.). The Nyquist theorem (.) is valid since

u ∈ L 2 and h ∈ L 1 , then W 1 H 1 s u * H 1 s h belongs to BL 2 . Both W 1 H 1 s
u and H 1 s h are band-limited finite energy functions so we are free to apply (.).

Since the sequence (ū * h) is the sampling of the band-limited

L 2 function W 1 H 1 s u * H 1 s h, it belongs to ℓ 2 (Lemma ). Finally, the interpolation I 1 (ū * h) is well defined.
The previous proposition shows that the image formation model (M) can be written in terms of discrete sequences.

ṽ = S s (ū D * h) + n = S s C[ū D ]h + n. (.) The digital image ūD = S 1 W 1 H 1 s u D is a well-sampled version of the distorted image u D = u•D.
The value s is the resampling rate from the high resolution lattice s×, where the psf estimation will take place, to the 1× sensor grid.



The numerical method will recover only a finite number of samples h of h. Strictly speaking h, being band-limited, cannot be compactly supported. Nonetheless, the error introduced by assuming that the support of h is bounded will prove negligible in comparison to the other sources of error: image noise, quantization, slight estimation errors of D, etc. Indeed, the retrieved solution h will prove to be experimentally independent from variations of its assumed support as long as it is large enough for errors to be negligible, and small enough for the operator to still be well conditioned.

When n is a zero-mean white discrete Gaussian noise, it follows from the previous formula that h e = (S s C[ū D ]) + ṽ is an unbiased estimator of h, as long as the linear operator S s C[ū D ] is injective. It can be shown that the estimator variance is proportional to the Hilbert-Schmidt norm of (S s C[ū D ]) (for matrices, the Frobenius norm), and that it is nearly minimal when ūD is a white noise realization (see Section .).

. PSF Estimation from an Unknown Pair of Scaled Images

Assume that we have perfect knowledge of the latent sharp image u that produced the blurry aliased observation ṽ. Under this non-blind assumption, solving for the psf amounts to solving an inverse problem governed by the image formation model (M). Of course, this would require the use of a specially designed calibration pattern. We are now interested in investigating to what extent the use of such pattern could be circumvented. We will propose a method that allows us to accurately estimate the psf from a pair of snapshots of the same scene, captured from different distances. In this method, the closest image will play a role similar to that of a calibration pattern in a classical non-blind approach.

In Chapter  we have shown that the highest accuracy in the psf estimation is obtained by using a realization of a Bernoulli white noise as calibration pattern. However, many highly textured scenes do exist in nature which, while not being optimal, may still lead to a well-posed inverse problem. In what follows, we prove that from two far apart snapshots of this kind of scene, complete recovery of the camera psf is theoretically possible based on the estimation of the blur between this pair.

.. Relative blur between two images: the inter-image kernel

Consider two digital images ṽ1 , ṽ2 of the same planar scene u, captured from different distances in a fronto-parallel position with negligible rotation around the optical axis. Let λ 1 and λ 2 denote the corresponding scale factors between the scene and each of the images. Then,

ṽi = S 1 H λ i u * h + n i for i = 1, 2 (.) = S 1 v i + n i = v i + n i ,
where v i := H λ i u * h and v i := S 1 v i . We will realistically assume that h ∈ L 1 ∩ BL 2 0 is nonnegative with h L 1 = 1, and u ∈ BL 2 0 (details on the appropriateness of these assumptions are given in Appendix .A). Also, it will be assumed that the acquisition distances are such that sλ 1 < λ 2 ; the importance of this assumption will soon become clear.



Definition . Let v 1 , v 2 ∈ BL 2 0 be two fronto-parallel continuous views of the same scene, acquired from different distances λ 1 and λ 2 respectively. We define an inter-image kernel between v 1 and v 2 as any kernel k ∈ BL 2 0 satisfying v 2 = Hλ 2/λ 1 v 1 * k.

The following lemma provides a characterization of the inter-image kernel.

Lemma . Let h ∈ L 1 ∩ BL 2 0 be non-negative, band-limited with supp( ĥ) ⊂ [-sπ, sπ] 2 and ĥ(0) = 1. Let ρ be the largest positive number such that | ĥ(ζ)| > 0 for every ζ ∞ < ρπ and assume that λ 2 ρ > sλ 1 . Then there is an inter-image kernel k ∈ BL 2 0 with support in [-sπ, sπ] 2 between (fronto-parallel views) v 1 and v 2 that satisfies

H λ h * k = h, where λ = λ 2 λ 1 . (.) If û does not vanish inside [-s π λ 2 , s π λ 2 ] 2
, then the inter image-kernel is unique and depends only on h and λ.

Proof. If k is an inter-image kernel between v 1 and v 2 , according to Definition  it must satisfy

F(H λ v 1 )(ζ) k(ζ) = v2 (ζ).
Since v i := H λ i u * h, the right-hand side of the previous equation is given by

v2 (ζ) = û ( ζ /λ 2 ) ĥ (ζ) .
In the same way, for the left-hand side,

H λ v 1 = H λ (H λ 1 u * h) (.) = H λ 2 u * H λ h, i.e., F(H λ v 1 )(ζ) k(ζ) = û ( ζ /λ 2 ) ĥ ( ζ /λ). Hence, û ( ζ /λ 2 ) ĥ ( ζ /λ) k(ζ) = û ( ζ /λ 2 ) ĥ (ζ) . (.)
It follows that a sufficient condition for k to be an inter-image kernel is ĥ ( ζ /λ) k(ζ) = ĥ (ζ) . Since h ∈ L 1 , ĥ is continuous. It follows that ρ is necessarily positive, since ĥ(0) = 1 > 0. In addition, as λ > s ρ by hypothesis, F(H λ h)(ζ) = ĥ( ζ /λ) does not vanish inside [-sπ, sπ] 2 and

k(ζ) = ĥ(ζ) ĥ( ζ /λ) (.)
is well defined all over its support, supp( k) ⊂ [-sπ, sπ] 2 . Finally, if û ( ζ /λ 2 ) does not vanish within the support of ĥ, from Eq. (.) k is unique.

Remark . In Lemma  it is assumed that the psf h is the same for the two images. This has at least two practical implications. First, we assume that both images are taken in perfect focus through proper refocusing. The only camera parameter allowed to change is the focus (aperture and focal distance remain unchanged). Second, the common area between v 1 and v 2 covers an important part of v 1 , and consequently its psf may exhibit some space variance that may degrade the estimation. Indeed, if v 1 is acquired through another psf h ′ such that ĥ′ (ζ) = ĥ(ζ) in [-s π λ , s π λ ] 2 , then Equation (.) will still be valid. Thus, the real requirement is that the low frequencies of the psf (i.e., frequencies in [-s π λ , s π λ ] 2 ) do not change.



.. Estimation of the inter-image kernel

The next goal is to estimate the inter-image kernel k. Since k is an s-band-limited function, we will work with its s× samples k = S 1 H 1 s k. We will show that under reasonable conditions, k can be recovered from the noisy aliased observations ṽ1 and ṽ2 . Let us first build up some intuition on how to derive the proposed estimator. In what follows, v1 = S 1 W 1 H λ s v 1 denotes a well sampled homothecy of parameter λ /s of v 1 .

Proposition . Under the assumptions of Lemma ,

v 2 = (S s C[v 1 ])k. (.)
Proof. Being k an inter-image kernel between v 1 and v 2 , it satisfies Eq. (.). Then,

v 2 = S 1 (v 2 ) = S 1 (H λ v 1 * k). Since k is s-band-limited, it follows that v 2 (.) = S 1 (W s H λ v 1 * W s k). (.)
Using the Nyquist-Shannon theorem for a band-limited signal and a set of properties detailed in Appendices .A and .B, yields

v 2 s>0 = S 1 H s H 1 s (W s H λ v 1 * W s k) (.) = S 1 H s (H 1 s W s H λ v 1 * H 1 s W s k) (.) = S 1 H s (W 1 H λ s v 1 * W 1 H 1 s k) (.) = S 1 H s I 1 S 1 (W 1 H λ s v 1 * W 1 H 1 s k) (.) = S 1 H s I 1 (S 1 W 1 H λ s v 1 * S 1 W 1 H 1 s k) def = S 1 H s I 1 (v 1 * k) def = S s (v 1 * k) def = S s C[v 1 ]k.
Of course, in practice we do not have access to v1 or to v 2 , but only to their noisy, aliased versions ṽ1 and ṽ2 . Thus k cannot be directly estimated from Eq. (.). However, a relationship between v1 and ṽ1 can be established as follows:

H λ s ṽ1 = H λ s (v 1 + n 1 ) + v1 -v1 = v1 + S 1 W 1 H λ s (I 1 v 1 -v 1 ) r +H λ s n 1 , (.)
where the last equality results from the definition of the discrete homothecy operator. The term r is a consequence of aliasing when sampling v 1 and introduces an unknown bias in the estimation of k. While this bias cannot be fully controlled, its impact can be mitigated. Indeed, since

r = S 1 W 1 H λ s (I 1 v 1 -v 1 ) (.) = S 1 H λ s W s λ (I 1 v 1 -v 1 ),



the aliasing term r will be non-zero only if there are aliasing components in the frequency interval s λ π, s λ π 2 . This allows us to choose

v 1 = H λ 1 u such that supp(v 1 ) ⊂ -2π + s λ π, 2π -s λ π 2 (see Figure .)
. Thus, to minimize the impact of the aliasing term the images should be acquired from a pair of fronto-parallel locations as far as possible one from the other, since that amounts to increasing the value of λ.

aliasing terms

Figure 3.1: Neglecting the aliasing. The estimation will be affected by aliasing only if there are aliasing components in the interval -s λ π, s λ π 2 . Hence, to avoid aliasing one can choose

v 1 = H λ1 u such that supp( v1 ) ⊂ -2π + s λ π, 2π -s λ π 2 .
From now on, we assume that the snapshots are acquired following the previous considerations. Therefore, we can ignore the aliasing term in (.), which leads to

v 2 = (S s C[v 1 ])k = (S s C[H λ s v 1 -H λ s n 1 ])k, that is (S s C[H λ s ṽ1 -H λ s n 1 ])k = ṽ2 -n 2 .
One could be tempted to solve for k in the previous equation using a total least squares based approach:

arg min k,δ,ǫ δ + κ ǫ subject to S s C[H λ s ṽ1 + δ]k = ṽ2 + ǫ. (TLS)
However, the particular structure of the operator S s C[H λ s ṽ1 + δ] makes this problem a difficult one. Instead we prefer to follow a simpler approach, which results from neglecting the noise term H λ s n 1 . This yields to the least squares estimation problem

arg min k,ǫ ǫ subject to S s C[H λ s ṽ1 ]k = ṽ2 + ǫ, (LS)
whose solution is given by

k e = S s C[H λ s ṽ1 ] + ṽ2 . (.)



If the noise n 1 is small compared to v 1 , this solution would be very close to the one that would be obtained from Problem (TLS). If, in addition, n 2 is small compared to v 2 , both solutions would be close to the actual inter-image kernel k = (S s C[v 1 ]) + v 2 . This follows directly from the continuity and injectivity assumptions on S s C[v 1 ], as a consequence of Lemma . This being said, we will consider the estimator of the inter-image kernel in Eq. (.).

Remark . If λ < s, the convolution between k and H λ s ṽ1 is not invertible so the operator S s C[H λ s ṽ1 ] will not be injective. This constraint on λ is necessary but not sufficient to make S s C[H λ s ṽ1 ] invertible. In addition, it is required that the spectrum of the image H λ s ṽ1 exhibits slow decay. Indeed, as shown in Section ., the flatter the spectrum of the image scene is, the better conditioned is the inverse problem. For that reason, in order to obtain accurate estimates of k, it is desirable that the chosen scene u exhibits white noise characteristics.

.. From relative to absolute blur

Now that we have a method for estimating the inter-image kernel k, we will concentrate on how to recover the camera psf. Notice that h is related to k by H λ h * k = h, and therefore its derivation is not straightforward. However, it holds that

h = lim n→∞ H λ n-1 k * H λ n-2 k * • • • * H λ k * k, (.)
as proved in the following proposition.

Proposition . Let h ∈ L 1 ∩ BL 2 0 and k ∈ BL 2 0 such that k(ζ) = h(ζ) h( ζ λ )
. Assume λ large enough to ensure that ĥ( ζ /λ) does not vanish in the support of k. Then if λ > 1, we have

lim n→∞ H λ n-1 k * H λ n-2 k * • • • * H λ k * k = h,
where the limit is in

L 2 ∩ C 0 . Proof. Let us call u n = H λ n-1 k * • • • * H λ k * k.
Then in the Fourier domain we have

lim n→∞ ûn (ζ) = lim n→∞ n-1 i=0 k ζ λ i = lim n→∞ ĥ(ζ) ĥ( ζ /λ n ) .
Since h ∈ L 1 , then ĥ ∈ C 0 and we have

lim n→∞ ĥ( ζ /λ n ) = ĥ(0) = 1
The convergence is uniform on a fixed compact set because ĥ is continuous and compactly supported. This implies that the convergence holds in L 1 and L 2 . Therefore

H λ n-1 k * H λ n-2 k * • • • * H λ k * k L 2 ∩C 0 -→ h.



This limit shows that it is possible to recover the camera psf h from the inter-image kernel k. Recall that in practice we have access only to discrete sequences; therefore it is convenient to derive a discrete equivalent of the previous limit. Since k is s-band-limited,

S 1 H 1 s (H λ k * k) (.) = S 1 H 1 s (W s H λ k * k) (.) = S 1 (H 1 s W s H λ k * H 1 s k) (.) = S 1 (W 1 H λ s k * H 1 s k) (.) = S 1 W 1 H λ s k * S 1 H 1 s k def = H λ k * k.
Iteratively applying this result to Eq (.) yields

h = lim n→∞ H λ n-1 k * H λ n-2 k * • • • * H λ k * k.
(.)

. The Complete PSF Estimation Procedure

This section describes the algorithmic steps that lead to local subpixel psf estimates.

Image alignment

In order to estimate the geometric transformation between both images, they need to be precisely aligned. This alignment can be obtained by matching sift descriptors [Lowe ], which have the advantage of being scale invariant.

Geometric transform estimation

The complete geometric transformation from one image to the other was approximated with a homography from the matched sift pairs. This permits the correction for deviations from the fronto-parallel assumption in the acquisition. Of course, if the  distortion is significant the assumed inter-image kernel Eq. (.) will not be accurate. The relative scale λ = (λ x , λ y ) is taken directly from the estimated homography. The IPOL implementation by Moisan et al. [] was chosen because of the efficiency of the optimized random sampling algorithm (orsa) rejection of false matches.

Gray level adjustment

Both snapshots should be acquired with exactly the same camera configuration and constant scene illumination. This ensures that there is no contrast change between them.

Resampling and distortion correction of ṽ1

The generation of the rescaled samples H λ s ṽ1 requires the interpolation of ṽ1 at the desired scale λ /s. This is done by using the estimated geometric transformation with bicubic interpolation. Notice that since ṽ1 is not very aliased, one can correctly interpolate it without introducing artifacts.

Numerical methods for inter-image kernel estimation

Suppose that the image ṽ2 has size m×n. The goal is to estimate k at s× the resolution of ṽ2 (camera sensor resolution). Also suppose that the estimated support of the inter-image kernel k is contained in an r × r patch. Then the matrix S s C[H λ s ṽ1 ] is of size mn × r 2 . A simple least squares procedure yields the inter-image kernel estimator:

k e = arg min k S s C[H λ s ṽ1 ]k -ṽ2 2 .
Transforming the kernel: from k to h Recovering the samples of the camera psf h amounts to evaluate the limit in Eq (.). Directly working with the digital sequences requires some care in how the successive convolutions are computed. Since λ > 1, the application of H λ would require a low-pass filter to avoid aliasing artifacts. To bypass this inconvenience one can restate the limit convolution as follows:

h = lim n→∞ H λ n (k * H 1 λ k * • • • * H 1 λ n-1 k * H 1 λ n k).
If implemented in this way, the successive discrete convolutions can be computed without any special care. To apply the discrete homothecy operator to k, we need to resample k using the Shannon-Whittaker interpolator. Because of its slow decay, in order to reduce ringing and other windowing effects, we opted to use bicubic interpolation. We get h by an iterative procedure that converges after a few iterations since λ n grows very fast. See Appendix A. for the algorithmic details.

In theory, as we already stated, the estimated psf should be non-negative. In practice, small negative values may be observed, due to deviations from model assumptions and numerical artifacts. To correct for these deviations, we simply set all negative values to zero.

. Experimental Results

Since there is no psf ground truth available, the validation of the proposed method was carried out by simulations and by comparing the results with state-of-the-art methods Joshi et al. []; LLC [] and with the proposed non-blind psf estimation algorithm presented in Chapter . Comparison was made only to non-blind, pattern based methods, as the accuracy of blind methods  is significantly lower. A complete algorithmic description, an online demo facility, and a reference source code can be found at the IPOL workshop by Delbracio et al. [a].

.. Simulations as a sanity check

A synthetic random image u was generated and re-interpolated 4× in order to get the "continuous" sharp homothecy of the image u. Next both images were convolved with a psf-like kernel (in this case a Gaussian isotropic kernel), and down-sampled to get the respective observed digital images at the camera resolution (i.e., 1×). The kernel was chosen so that the low resolution image presents aliasing artifacts. By generating the views of u in this way, there are no aliasing artifacts in the closest image. This experiment was done as a sanity check of the proposed method. A 4× kernel was estimated from the observed image pair. The results are shown in Figures . and ..

The procedure was tested for both automatic sift-based registration and the ideal (known) alignment. Although both estimates are significantly accurate, the automatic registration introduces a small misalignment, as shown in the difference images. See the caption of Figures . and . for details.

.. Real camera examples

The behavior of the proposed approach was tested for several different image pairs and for superresolution estimations ranging from 1× to 4×. The experiments were performed using a a Canon EOS D camera equipped with a Tamron AF -mm F/. XR Di-II lens. The focal length was fixed to . mm.

Two-scale versus non-blind pattern based method

In Chapter  we presented a non-blind method that uses a realization of white noise as calibration pattern. It was proved that, until to now, this method is the one tat estimates the psf with highest accuracy. Therefore, the psf resulting from this method will be used here as ground truth. Figure . shows the 4× psf estimated by the proposed two-scale method from a pair of views of a wall shown in Figure .. The estimation was conducted for one of the green channels (half of the green pixels of the Bayer matrix), with the camera aperture set to f/.. The estimated psf is quite close to the one obtained by using our Bernoulli pattern based approach. In particular their sizes are similar, and their corresponding mtfs present zeros at the same locations. Notice that the red channel psf is wider than the green and the blue one, as expected from the physics of diffraction-limited optical systems, since the wavelengths associated to red light are larger than the rest. The differences between the dominant orientations of the red, green and blue psfs spectra can be explained by the sensor shape and layout. In fact, each sensor active zone is usually L-shaped, and the red and blue sensors are rotated 90 • w.r.t. the green ones (see, for example, Yadid-Pecht []). These rotations are consistently observed in the psfs and mtfs estimated with our two-scale method. This clearly illustrates the accuracy of the proposed approach.

Color filter array estimations



Different kinds of scenes

The wall images in the previous experiments are well adapted for our two-scale psf estimation method, since their spectra show slow decay. A priori one would think that images from pure white noise would yield better estimates, since this is what happens in our previous pattern based approach. But for our two-scale approach, this would be true if both snapshots could be precisely aligned, which is not the case in practice. Indeed, sift descriptors are not stable in the presence of aliasing. Hence, there is a trade-off between having accurate sift matches and textures with high frequency information. The texture shown in Figure . is an example of an appropriate trade-off.

Figure . shows two snapshots of a photograph in a magazine, with the corresponding 1× to 4× psf estimations for the first green channel. The estimation was performed at the image center for the camera working at f/. aperture. All the subpixel estimations are consistent: their mtfs exhibit good overlap in common regions. While these newspaper images produce accurate sift points, their spectra decay faster than those of the wall images. Consequently, the high frequencies in the psf estimate are noisier. This can be readily seen by comparing both estimates at 4× resolution.

What kind of textures should be used? It follows from the previous analysis that, in order to simultaneously produce good sift points and a sufficiently slow frequency decay, textures composed of elements with different sizes are to be preferred. Three-dimensional (3d) textures like those shown in Figure . can be problematic for this approach. Even though they respect the two previous conditions, their 3d nature produces disparities, and occlusions which change the image beyond a simple zoom. Likewise, non-Lambertian surfaces and dynamical scenes are not appropriated either.

Comparison to other methods

In this experiment we compare the performance of the twoscale method proposed here with three state-of-the-art non-blind methods: that of Joshi et al.

[], Imatest commercial software LLC [], and our previous Bernoulli pattern method. All the estimates were computed at the image center, with aperture f/.. For the two-scale approach, we used the wall image pair shown in green g

-2 -1 0 1 2 -2 -1 0 1 2 green g -2 -1 0 1 2 -2 -1 0 1 2 blue b -2 -1 0 1 2 -2 -1 0 1 2 red r -2 -1 0 1 2 -2 -1 0 1 2 Figure 3
.7: Different color channels. PSF estimation at 4× resolution for the four Bayer channels (two greens, red and blue). Top row: two distant, parallel views of a concrete wall (close ups, the full images are shown in Figure 1.5). Middle row: the 4× PSF estimated for the four channels. Bottom row: their corresponding Fourier spectrum modulus. The estimation was performed with images captured at aperture f/5.6. The red PSF is larger than the green and blue ones. This is consistent with the diffraction phenomenon: the red wavelengths are larger than the rest, thus the diffraction kernel is wider. Also note the differences between the shape of the red, blue and green PSF spectra (bottom row). Red and blue MTFs are rotated 90 • with respect to the green ones. This symmetric behavior is consistent with the layout of L-shaped sensors Yadid-Pecht [2000].



. Conclusion

In this chapter we presented an algorithm for the subpixel estimation of the point spread function (psf) of a digital camera from aliased photographs. The procedure is based on taking two frontoparallel photographs of the same flat textured scene, from different distances leading to different geometric scales, and then estimating the kernel blur between them.

The estimation method is regularization-free. In that sense, the technique is closely related to the non-blind estimation method presented in Chapter , which uses a random noise pattern. The main difference is that non-blind methods can directly estimate the psf using the perfect knowledge of the pattern. In the proposed two-scale method the question is far more intricate because only the blur between the acquisitions can be estimated. Thus a mathematical analysis and new algorithms have been introduced proving how the psf can be recovered from the interimage kernel.

To reach high accuracy, images of textured scenes with sufficient flat spectra are preferred. It was experimentally verified that many textures found in nature are well adapted to these requirements. A comparison of the resulting psf estimates with other subpixel psf estimation methods shows that the proposed algorithm reaches accuracy levels similar to those of state-of-the-art methods, with the advantage of not requiring any special acquisition setup or calibration pattern and thus being much more practical.

.A Mathematical Framework and Physical Modeling Functional spaces and other notation

• R 2 is the set of pairs of real numbers x = (x 1 , x 2 ), and Z 2 the set of pairs of integers Joshi et al. [2008], Imatest commercial software LLC [2010], our previous Bernoulli pattern method, and the two-scale method proposed in this work (applied to the images of the wall shown in Fig. 3.7). On the low frequencies all algorithms produced very similar estimates, while on the higher frequencies the Joshi et al. estimation depends strongly on the regularization level. Although much effort was made 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. The proposed two-scale algorithm is the one presenting an estimation closest to the non-blind estimation presented in Chapter 2, considered as ground truth by virtue of its high accuracy.

k = (k 1 , k 2 ). L 1 (R 2 ) is the set of integrable functions on R 2 , L 2 (R 2 )
 integrable functions, C 0 b (R 2 ) is the set of continuous bounded functions, C ∞ (R 2
) is the set of infinitely differentiable functions, S(R 2 ) is the Schwartz class of C ∞ functions whose derivatives of all orders have fast decay, S ′ (R 2 ) is its dual, the space of tempered distributions, E ′ the subset of S ′ (R 2 ) of compactly supported distributions. We shall use the properties of the convolution

L 1 * L 2 ⊂ L 2 , L 1 * L 1 ⊂ L 1 , L 2 * L 2 ⊂ C 0 , E ′ * S ′ ⊂ S ′ .
• We denote by BL 2 (R 2 ) (or BL 2 for short) the set of L 2 functions that are band-limited inside [-π, π] 2 . More generally, BL 2 0 denotes the space of L 2 functions with compactly supported Fourier transform.

The following conventions and notations will be used in what follows:

• F is the Fourier transform operator defined on S ′ ; F(f

)(ζ) = f (ζ) = e -ix•ζ f (x)dx defines it for a function f ∈ L 1 (R 2 ) in a point ζ = (ζ 1 , ζ 2
). This formula is still valid for functions belonging to L p (R 2 ) with 1 < p ≤ 2 (see, e.g., the reference books by Stein and Weiss []; Bony []).

• Continuous images are defined for x ∈ R 2 , whereas digital images are sampled on a discrete grid k ∈ Z 2 . realistic assumption is to consider them non-negative.

• S 1 :

C 0 b → ℓ ∞ (Z 2 ) is the 1-sampling operator such that u(k) = (S 1 u)(k).
From the distribution viewpoint S 1 is the product by a Dirac comb Π s := k δ sk with s = 1, namely S 1 u = Π 1 .u where u must be a continuous function. Both representations of the sampling operator will be identified, and it will be clear from the context which representation is intended.

• A digital image u will be represented either as a sequence (u(k)) k in ℓ ∞ (Z 2 ) or as the corresponding Dirac comb u := k∈Z 2 u(k)δ k .

• The operator I 1 : ℓ 2 (Z 2 ) → BL 2 (R 2 ) denotes the Shannon-Whittaker interpolator, defined by I 1 u(x) = k∈Z 2 u(k)sinc(xk), where sinc(x) = sin(πx) πx sin(πy)

πy . We therefore have

I 1 u = F -1 ( k u(k)e -ik•ξ 1 [-π,π] 2 ).
When u ∈ ℓ 2 , F(I 1 u) belongs to L 2 and is compactly supported. Thus I 1 u ∈ BL 2 , and we have S 1 I 1 = Id.

• The filter

W w u = F -1 (û • 1 [-wπ,wπ] 2 ) is an ideal low-pass filter that cuts the spectrum of u to [-wπ, wπ] 2 . It is defined if û is a function. Note that if u ∈ L 1 ∪ L 2 , then W 1 u is in BL 2 .
• H λ u(x) = λ 2 u(λx) is the continuous homothecy (i.e., λ > 1 is a contraction); the rationale for its normalization is to preserve the image mean (its zero-frequency coefficient).

In the Fourier domain

F(H λ u )(ζ) = û( ζ λ ), so if u is α-band-limited, then H 1 α u is band- limited.
• S s : ℓ 2 (Z 2 ) → ℓ 2 (Z 2 ) denotes the s-to--resampling operator S s = S 1 H s I 1 (i.e., s > 1 is a subsampling by s).



• C[u] : ℓ 2 (Z 2 ) → ℓ 2 (Z 2 ) denotes the linear map associated to the convolution with a digital image u. The convolved sequence belongs to ℓ 2 (Z 2 ) which in general is satisfied if u ∈ ℓ 1 (Z 2 ).

• The digital Nyquist homothecy operator

H α : ℓ 2 (Z 2 ) → ℓ 2 (Z 2 ) is defined by H α u := S 1 W 1 H α I 1 u. It is a digital contraction if α > 1.
• Let L be a bounded linear operator over a Hilbert space. L * is its adjoint and L + (if it exists) is its pseudo-inverse, i.e., the minimum-norm solution of (L * L)L + := L * .

Physical and mathematical modeling of continuous images

Continuous images will be assumed to be functions in BL 2 0 (R 2 ). This choice is consistent, since these functions are continuous (actually C ∞ ) and the sampling is well defined. Moreover, as suggested by Morel and Ladjal [] and later by Almansa et al. [, Appendix A] this choice is sufficiently general to model the continuous landscape observed by a camera just before sampling takes place at the sensors.

In fact, even if the raw physical image before blur and sampling is, realistically, a positive Radon measure O (due to the photon-counting nature of sensitive digital systems) with compact support (imposed by the finite number of photons), it will still be blurred by the camera psf h which will be regular enough for h * O to be in BL 2 0 . How regular can it realistically be assumed to be? The kernel h originates in several physical phenomena from diffraction, anti-aliasing filtering and sensor integration. Each one of these phenomena, and their combination as well, lead to model h as a nonnegative function with finite mass h = 1 (normalized to ). In addition the diffraction part ensures that ĥ is compactly supported. From this one deduces that h ∈ BL 2 0 ∩ L 1 . We now turn to the problem of simplifying O to a more manageable function u, which is indistinguishable from O after convolution with the psf h. Let B = supp( ĥ) be the (compact) spectral support of the psf h. Hence h can be idempotently written as h = h * h 0 , where h 0 ∈ S ′ has a compactly supported spectrum satisfying ĥ0 (η) = 1 for η ∈ B. The function ĥ0 can easily constructed by an explicit formula as a C ∞ and compactly supported function satisfying ĥ0 (η) = 1 on B. Then its inverse Fourier transform has all required properties.

So we have

v = h * O = h * u, where u = h 0 * O.
Consequently, the observed landscape can be assumed without loss of generality to be u = h 0 * O instead of O. Being the convolution of a compactly supported positive Radon measure O ∈ E ′ with h 0 ∈ BL 2 0 ∩ L 1 , u also belongs to BL 2 0 , and its convolution with

h ∈ BL 2 0 ∩ L 1 is the observed image v ∈ BL 2 0 .

Standard results from Fourier analysis

The following two main results from standard Fourier analysis and distribution theory are stated without proof. The reader is referred to [Stein and Weiss ; Hörmander ] for the proofs in the particular setting chosen here.

Proposition  (convolution through Fourier transform). The relation [Stein and Weiss , Theorem .]).

F(f * g) = F(f ) • F(g) (.)  is valid in either of these cases . g ∈ L 1 (R 2 ) and f ∈ L p (R 2 ) for 1 ≤ p ≤ 2. Then f * g belongs to L p (R 2 ) (see
. g ∈ E ′ and f ∈ S ′ . Then f * g belongs to S ′ (see [Hörmander , Theorem ..]).

Applying the Fourier transform on both sides of (.) and recalling that the squared Fourier transform operator F 2 (u) = (2π) 2 [x → u(-x)] is almost the identity (except for flipping and a constant factor), we obtain the following:

Corollary  (product through Fourier transform). The relation

F(f • g) = 1 (2π) 2 F(f ) * F(g) F -1 (f • g) = F -1 (f ) * F -1 (g) (.) holds when ĝ ∈ E ′ and f ∈ S ′ . Then f • g belongs to S ′ . Proposition  (Poisson formula in R 2 for tempered distributions [Hörmander ]). Π1 = (2π) 2 Π 2π . (.) Lemma . If û ∈ E ′ , then F(Π 1 • u) = Π 2π * û. (.)
Proof. We can apply the first form of Corollary  where f = Π 1 ∈ S ′ and ĝ = û ∈ E ′ to obtain

F(Π 1 • u) = (2π) -2 Π1 * û = Π 2π * û
where the last equality is deduced from the Poisson formula (.).

The Shannon-Whittaker sampling theorem is then a direct consequence of the two previous results.

Proposition  (Nyquist-Shannon theorem). If u ∈ BL 2 (R 2 ), then

u = I 1 S 1 u. (.)
Proof. We can apply Lemma . Multiplying both sides of Equation (.) by

F(sinc) = ✶ [B] we obtain F(sinc) • F[S 1 u] = F(sinc) • [Π 2π * û] = k∈Z 2 û(• + 2πk)✶ [B] = û,
where in the right-hand side the only non-null term is for k = 0 because u is band-limited in

B = [-π, π] 2 and F(sinc) = ✶ [B]
. Finally, using the second form of Corollary , we obtain

sinc * (S 1 u) = u,
and the left term is by definition I 1 S 1 u.

Corollary . If u ∈ L 2 is s-band-limited, then u = H s I 1 S 1 H 1 s u. (.)



.B Proof of Auxiliary Results

Common hypotheses According to the discussion in Appendix .A, and in order to justify all the lemmas and propositions we will require that h ∈ BL 2 0 ∩ L 1 (R 2 ), non-negative ĥ(0) = 1;

-u ∈ BL 2 0 .
This ensures that the convolution u * h = v is well defined with u ∈ BL 2 0 . For the uniqueness of the inter-image kernel we shall additionally assume that û does not vanish inside [-s λ 2 π, s λ 2 π] 2 .

We now prove several properties that are used throughout this chapter.

Lemma . If u ∈ BL 2 0 , then S 1 u ∈ ℓ 2 (Z 2 ).
Proof. As u in BL 2 0 there exists s > 0 such that û ⊂ [-sπ, sπ] 2 . Furthermore, since û ∈ E ′ applying (.) we have F(S 1 u) = (2π) 2 Π 2π * û. Since u belongs to L 2 , then û is again in L 2 . Thus, Π 2π * û is the 2π-periodic version of a L 2 function in [-sπ, sπ] 2 . Consequently the inverse Fourier transform of Π 2π * û is a Dirac comb whose coefficients are the Fourier series coefficients of û. Thus the coefficients of S 1 u form an ℓ 2 sequence.

Proposition . Let h ∈ L 1 (R 2 ) and u, v ∈ L 1 ∪ L 2 (R 2 ).
The following equalities hold:

W 1 (h * v) = W 1 h * v = h * W 1 v, (.) W 1 H λ v = H λ W 1 λ v, (.) H α (u * v) = H α u * H α v. (.)
Proof. This is the proof of (.). In the Fourier domain,

F(W 1 (h * v)) def = F(h * v) • 1 [-π,π] 2 (.) = F(h) • F(v) • 1 [-π,π] 2 .
Thus,

F(h) • F(v) • 1 [-π,π] 2 = F(h) • 1 [-π,π] 2 • F(v) • 1 [-π,π] 2 ,
and all results are deduced from this last statement.

Proof. This is the proof of (.). Since

F(H λ v) = λ 2 F(v(λ•)) = λ 2 1 λ 2 v( . λ ) = λ 2 H 1 λ v,
we have

F(W 1 H λ v) (.) = F(H λ v) • 1 [-π,π] 2 = λ 2 H 1 λ v • 1 [-π,π] 2 .
On the other hand,

F(H λ W 1 λ u) = λ 2 H 1 λ F(W 1 λ v) (.) = λ 2 H 1 λ (û • 1 [-π λ , π λ ] ) = λ 2 (H 1 λ v) • 1 [-π,π] 2 .



Proof. This is the proof of (.). The proof is a mere change of variables:

H α (u * v)(x) = α 2 u(s)v(αx -s)ds = α 4 u(αs)v(αx -αs)ds = (H α u * H α v)(x). Lemma . Let u, v ∈ BL 2 0 (R 2 ). If either u or v is band-limited, then S 1 (u * v) = S 1 ū * S 1 v, (.)
where we have called ū = W 1 u and v = W 1 v .

Proof. We will prove this statement in the tempered distribution sense. We will consider

S 1 u = Π 1 • u = k δ k • u as a Dirac comb.
The application of S 1 to ū, v and u * v is well defined as all functions are in BL 2 (R 2 ) and by consequence they are in

C ∞ . Recall that if u ∈ D ′ and f is C ∞ then f • u ∈ D ′
thus in this framework we need a function to be in C ∞ to be sampled.

From Lemma  we know that the sequence of coefficients from S 1 ū and S 1 v are in ℓ 2 (Z 2 ). Thus (S 1 ū) * (S 1 v) is a bounded sequence and therefore every term is well defined.

Finally

F(S 1 (u * v)) = Π 2π * (û.v) = (Π 2π * û) • (Π 2π * v)
is true because all considered functions happen to be 2π-periodizations of compactly supported functions in (-π, π) 2 , namely û, v, and their product.

Lemma  (stability of the inter image kernel estimation). Let A be an injective bounded linear operator (iblo) defined on a Banach space X and ∆A a perturbation of A such that A + ∆A is also iblo and A ∆A < 1. Let b ∈ X and ∆b be a perturbation of b. Then, the solution of x = A + b and x * = (A + ∆A) + (b + δb) satisfy:

x * -x x ≤ cond(A) 1 -A + ∆A δb b + ∆A A , (.)
where cond(A) = A A + .

Proof. First note that as A is full rank the pseudo-inverse is the left inverse of A, namely A + A = I. Since A ∆A < 1 we have that (A + ∆A) + = (I + A + ∆A) -1 A + and we also have

(I + A + ∆A) -1 = (A + ∆A) k ≤ (A + ∆A) k = 1 1 -A + ∆A .
Hence,

x * -x = (A + ∆A) + (b + δb) -A + b = (I + A + ∆A) -1 A + (b + δb) -A + b; therefore (I + A + ∆A)(x * -x) = A + (b + δb) -A + b -A + ∆AA + b = A + (δb -∆Ax);
 and then

x * -x x ≤ A + 1 -A + ∆A δb + ∆Ax x = cond(A) 1 -A + ∆A δb + ∆Ax A x ≤ cond(A) 1 -A + ∆A δb Ax + ∆A x A x ≤ cond(A) 1 -A + ∆A δb b + ∆A A .





Part II

Accelerating realistic image synthesis



This chapter describes a new multi-scale filter for accelerating Monte Carlo renderers. Each pixel in the image is characterized by the colors of the rays that reach its surface.

The proposed filter uses a statistical distance to compare the ray color distribution associated with each pixel at each scale. Based on this distance, it decides whether two pixels can share their rays. This simple and easily reproducible algorithm provides a psnr gain of  to  decibels, or equivalently accelerates the rendering process by using  to  times fewer samples without observable bias. The algorithm is consistent, does not assume a particular noise model, and is immediately extendable to synthetic movies. Being based on the ray color values only, it can be combined with all rendering effects.

. Introduction

In computer graphics, producing high quality realistic images in a reasonable amount of time is still a major challenge. The goal of a global illumination algorithm is to estimate the light distribution in a scene. The color of each pixel in the image results from the superposition of light rays transported by an infinite number of paths that lead to it, either directly from light sources, or indirectly after bouncing in the scene. The light distribution in a scene can be obtained as a solution of the rendering equation, an integral equation that models the radiance equilibrium as a light transport in a scene [Kajiya ]. Solving this equation for real scenes is an intractable problem. Approximate solutions are usually obtained by Monte Carlo numerical integration techniques where image pixels are formed by averaging the contribution of stochastic rays cast from the camera through the scene. The main limitation of Monte Carlo rendering is that the variance of the estimator decreases linearly with the number of stochastic samples. Thus the root mean squared error of the estimated image decreases as the square root of the number of primary rays cast from the camera (which we call samples from now on). While variance reduction techniques such as importance sampling, Russian roulette, or Markov Chain Monte Carlo methods can be used to accelerate convergence, still several hours or even days may be necessary to produce noiseless photorealistic images.

To reduce the time required by Monte Carlo rendering to produce good quality images, two main strategies have been proposed, that may be called adaptive rendering and rendering postprocessing. In the first strategy, the idea is to act during the rendering process by locally adapting the number of rays cast per pixel, depending on the complexity of certain zones. Post-processing is applied once rendering has been completed, and mainly consists of filtering or interpolating either samples or pixels. Both strategies can be combined.

It is worth noting that the target quality of the images may vary depending on the application.



The quality required for pre-visualization, where important time constraints have to be met, is clearly not as high as for applications where photo-realistically rendered scenes are an objective by themselves. In pre-visualization scenarios, computational time reduction is obtained by using the renderer to produce only a very small number of samples (say,  to  samples per pixel). In order to produce images of high enough quality from such a sparse and noisy data, it is necessary to filter or re-synthesize samples using as much information as possible. Indeed, for each sample the rendering system keeps track of relevant information associated to the ray path: geometric, color and texture features, object and material properties, Monte Carlo random parameters, etc. Using these fat samples, state of the art methods such as the ones proposed by Sen and Darabi [] or by Lehtinen et al. [] produce spectacular results. However, the quality of the results obtained by this approach remains scene dependent, being potentially affected by the strong under-sampling of high dimensional data. The larger the number of effects that are simultaneously present, the higher the risks of this under-sampling. Proper up-sampling or interpolation of the sample space is therefore only possible under strong regularity conditions on the fat samples distribution. This explains why the best performances are observed for highly diffusive scenes (where impressive results are obtained from only one sample per pixel). As pointed out by Lehtinen et al. [], poor performance is instead expected when the scene contains high frequency illumination effects, incompatible with a low sampling rate.

In short, the generation of high quality images, specially when simulating complex effects such as anti-aliasing, indirect illumination, depth-of-field, motion, requires a large number of rendered paths to correctly sample the path space. The required number of fat samples is certainly too large, not only because of the computational time that would be required to process them but, most fundamentally, because its memory storage would exceed any reasonable capacity limit (more than 100 bytes per fat sample [Sen and Darabi ; Lehtinen et al. ]). The natural alternative is to give up using fat samples, and to store only part of their information. In the limit, the information can be reduced to color samples, that is the final color transported by each ray when hitting the screen. In this case, we say that the method works on the screen space, as opposed to the previous methods which work in the space of paths. Working on the screen space allows one to avoid memory saturation, while keeping a number of samples which may be large enough to capture the sample space variability. The works by Rousselle et al. [], Dammertz et al. [] and Overbeck et al. [] are among the most representative ones of this kind of approach.

In the present work we propose and study a new, intermediate, filtering approach that works on the screen space but keeps and uses the color samples at each pixel. Thus, it can be coupled with any Monte Carlo renderer keeping a record of the samples color. The cornerstone of the proposed method is to measure the similarity between any two pixels as the statistical distance between the histograms of rays color that hit them. If the comparison is positive, the ray color histograms of the similar pixels can be fused. The final color of a pixel is then obtained as the average of a all rays color of all similar pixels. This fusion is made still more reliable by comparing patches instead of pixels, and by allowing long range interaction by a multiscale procedure. The ray color histogram characterizes better the physical and geometric properties of a pixel than just its color or the color of its neighbors. The proposed approach is related to bilateral filters, which were first applied to denoise Monte Carlo rendered scenes by Xu and Pattanaik []. The idea of comparing patches instead of individual pixels goes back to Buades et al.'s Non-Local Means . Our approach still presents a fundamental difference with Non-Local Means or any classic variation of bilateral filters proposed in image processing or computer graphics: instead of defining similarity by computing distances between pixels color, we compute distances between color distributions and fuse them when the distance is small enough. For this reason, we will call this method ray histogram fusion  rhf. Distances based on distributions are of course much more informative than comparing just their averages, as bilateral filters do.

rhf is simple, easy to implement, and therefore fully reproducible. It is independent of the sample generation process. It can be easily coupled with any renderer and even with any other acceleration method. Most importantly, the method does not make any particular assumption on the scene. As will be demonstrated by our experiments, it therefore copes with a wide range of scenes and multiple simultaneous effects. Finally, its time and memory complexities grow linearly with the image size and are independent of the number of input samples.

The limitations of rhf are also clear. Its performance depends on the degree of self similarity of the scene, which fortunately is usually high [Lebrun et al. ], and the price to pay for its generality is the requirement of a relatively large number of input samples.

The plan of the chapter is as follows. In Section . we review relevant previous work. Section . defines a pixel similarity measure based on the corresponding cast rays color, and discusses its statistical interpretation. Section . describes the rhf algorithm, and shows how it successfully makes use of the whole ray color histogram information. Section . reports quantitative and qualitative results on the algorithm performance. We close with Section ., discussing limitations of our approach and outlining future work and conclude in Section ..

. Previous Work

A thorough analysis of Monte Carlo rendering is far beyond the scope of the present work. The interested reader may consult the introductory book by Dutré et al. [] and the one by Pharr and Humphreys []. However, for what follows it is enough to note that there are mainly two approaches to reduce the time required by Monte Carlo rendering to obtain good quality images.

One of these approaches is adaptive sampling. This class of algorithms locally adapt the number of rays cast per pixel. The idea is to increase the number of samples in complex parts of the scene while maintaining a reduced number in simple parts, such as flat regions. Complex textures or defocused zones are typical elements that require large amounts of rays to be properly rendered. Hachisuka et al. [] (mdas) proposed to adaptively distribute a set of samples in the full, multidimensional sampling domain where the rendering equation is computed. However, as more Monte Carlo effects are considered (e.g. depth of field, motion blur, area lighting, etc.) the dimension of this space will be larger and thus will suffer from the curse of dimensionality. One of the most significant adaptive sampling algorithms is certainly the Adaptive Wavelet Rendering (awr) by Overbeck et al. []. This method adaptively distributes Monte Carlo samples in the screen space to reduce the variance of a wavelet basis scale coefficients. Then, the image is reconstructed from these non-uniformly distributed samples by using a suitable wavelet approximation. Soler et al. [] proposed to analyze the depth of field effect in the Fourier domain. By properly predicting the local bandwidth their algorithm adaptively samples the multidimensional domain. In a similar fashion, Egan et al. ; b; a addressed motion blur, soft shadows and directional occlusions respectively, by adaptively sampling the multidimensional domain followed by a sheared reconstruction. This allowed reusing samples between pixels in specific effects.

The reconstruction scheme proposed by Rousselle et al. [] attempts to minimize the mean squared error. The idea is that given the current distribution of samples in the screen space, the algorithm chooses the best reconstruction filter (among a set of predefined filters -e.g. Gaussian filters) at each pixel to minimize some error criterion. Next, given the current filter selection, new samples are distributed to minimize the error. Thus, this algorithm is both an adaptive sampling  and reconstruction filtering. This state of the art algorithm will be used in the experimental section for comparison.

The other approach is denoising or adaptive filtering. In this family of algorithms, the existing set of samples are combined to produce a better estimator of the pixel color using sample information in a pixel. Adaptive filtering may take place at sample level (i.e., primarily filtering samples) or at pixel level (i.e., primarily filtering pixel values). The majority of these methods can actually be written as generalizations of the bilateral filter [Paris et al. ] applying a weighted average of the samples (resp. of the pixels) in a neighborhood. The complexity of the method depends on whether it is applied at a pixel or sample level and how deep the method digs into the rendering information (e.g., information about each sample history: color, normal, object of the last impact; information about the random parameters used to calculate the sample) in order to compute the weights of the samples. In order to show how the conception of the proposed filter appears as the natural evolution of the previous work, we briefly present significant contributions, from a historical perspective. We will see that the general trend in this evolution is to rely more strongly in the auxiliary information available from the rendering system.

The simplest adaptive filters act at pixel level, like any filter used in classical image processing. Lee and Redner [] presented a seminal work defining an alpha trimmed filter (a generalization of the median/mean). Jensen and Christensen [] proposed to apply Gaussian or median filters with 3 × 3 pixels support to light having been reflected diffusely at least twice. The trilateral filter of Choudhury and Tumblin [] involves an adaptive neighborhood function and the image gradient. Again a classic image bilateral filter was proposed by Xu and Pattanaik []. Notice that unlike the work by Lee and Redner [], classical bilateral filters cannot remove outliers. To overcome this limitation, the weights of the bilateral filter by Xu and Pattanaik are computed based on a denoised version of the original image.

More complex filters make use of sample information available from the renderer in order to filter still at a pixel level. Rushmeier and Ward [] proposed to spread out noisy pixels (e.g. pixels whose variance is larger than a threshold after a fixed number of iterations) into a region of influence. A noisy pixel will contribute to several denoised output pixels, and since the filter is normalized no energy will be leaked. McCool [] proposed another classical filter that uses pixels geometric information. It is an anisotropic diffusion (of the Perona-Malik type) removing noise from Monte Carlo rendering. The conductance function that models the strength of the diffusion scheme in a pixel is estimated from a coherence map using depth and normal information gathered during rendering (contained in the G-buffer) along with a color coherence map. More recently Dammertz et al. [] presented a fast wavelet filtering scheme designed for ray traced Monte Carlo global illumination images. For that purpose the filter uses rt-buffer information about direct or indirect illumination, and the buffer information on normals and position. The bilateral filter is also invoked by Xu et al. [] to denoise images created with complex light paths in smoke or fog. In this work, additional bilateral weights based on the path gradient direction are used to better guide the denoising scheme.

The last class of filters uses the additional sample information to adaptively filter the sample values. Shirley et al. [] addressed the question of noise in defocused or motion blurred regions. The image filter is adapted to the a priori knowledge of the kind of blur in a given image region. This is a very natural and successful ad hoc strategy for these regions. Probably the most impressive results are those recently reported by Sen and Darabi []. This method uses the whole information of the rendering process and the whole information on each numerical photon to denoise by bilateral filtering. The bilateral filter takes simultaneously into account in its weights the sample position and spatial neighborhood in the image, the random synthesis parameters, the scene features (nor- mal, world space position, texture values) and finally the sample color. It computes as a mutual information the statistical dependence on the random generation parameters of the pixels sharing the same features and colors. Although the results are outstanding at very low samples per pixel, the complexity of the method makes it not scalable to generate high quality images from a large number of input samples.

Lehtinen et al.

[] described a reconstruction technique that allows rendering a combination of motion blur, depth of field and soft shadows by exploiting the anisotropy in the temporal light field. The effective sampling rate is increased by a large factor by efficiently reusing samples between pixels. Recently Lehtinen et al. [] generalized these ideas to deal with indirect illumination. By contemplating the properties of diffuse surfaces, their algorithm permits to interpolate the light field to produce results similar to those that would have been obtained by rendering a much larger number of samples. For instance, from an input image of  samples per pixel, they synthesize images of  samples per pixel, whose quality is similar to  samples per pixel generated by standard path tracing. While the quality increase is impressive, the noise level in these images is still too strong for applications requiring high quality images. Our algorithm is somehow complementary to this approach. Indeed, it can be used to boost the performance of a pure Monte Carlo renderer or any other set of samples like the ones generated by Lehtinen et al. [], and can perfectly deal with a number of samples in this order of magnitude.

The above bibliographical analysis has shown that most Monte Carlo denoising methods are generalizations of the bilateral filter (or sigma-filter [Lee ]). The general principle behind the bilateral filter is that similar pixels must be denoised jointly, being different samples of the same model. This is also implicitly used by the sigma-filter and by the nl-means algorithm [Buades et al. ]. In computer graphics, ray information permits to identify still more rigorously than in classic image processing the pixels sharing the same model. Indeed, all ray samples hitting a given pixel and its neighbors can be used for that purpose.

. Proposed Approach

.. Rationale

In contrast to classical photography where only the energy arriving at the sensor plane can be measured, in a rendering scenario much more information about pixel formation is available. In particular, the light contribution and the screen position of each path can be stored, as well as the associated geometrical and scene information about the objects encountered along the ray path.

As pointed out by Veach [], the light transport problem can be stated in the space of paths, and global illumination can be estimated by computing a transport measure over each individual path. Under this path integral formulation, each pixel color u(x) = (u R (x), u G (x), u B (x)) is given by the integral over all possible light paths

u(x) = Ωx f (p)dµ(p),
where Ω x is the space of paths originated at pixel x, p is a path of any length, and dµ(p) is a measure in the path-space. The function f (p) describes the energy contribution through a path p and is the product of several scene factors due to the interaction of light within the path plus initial self-emitted radiance and importance distributions. Thanks to this formulation, the image color at pixel x can be estimated from n x random paths p 1

x , . . . , p nx x , generated by an appropriate Monte



Carlo sampling procedure. If c j x denotes the color transported by random path p j x (for instance, in path tracing c j x = f (p j x )), the Monte Carlo approximation of u(x) is computed as

ũ(x) = 1 n x nx j=1 c j x . (.)
Consider now the Monte Carlo approximation error n(x), given by

n(x) = ũ(x) -u(x). (.)
The Monte Carlo approximation is asymptotically unbiased, but the mean squared error E[n 2 (x)] decays linearly with the number of samples n x . Consequently, unless the rendering system spends several hours or even days producing samples, the resulting images will be contaminated by white noise. This is a consequence of the fact that the Monte Carlo samples are independent and therefore the random process {n(x), x image pixels} is white.

One possibility to reduce the approximation error while keeping the rendering time reasonable is to render fewer samples, and to filter the pixel values afterwards. Filtering will always result in a significant variance reduction, however, it may also severely increase the approximation bias. The only filtering processes that do not introduce bias are those that combine pixels of the same "nature", that is pixels x having the same ideal value u(x). While identifying two similar pixels x and y based on the unknown pixel values u(x) and u(y) is of course impossible, it is reasonable to expect that their samples color {c 1

x , . . . , c nx x } and {c 1 y , . . . , c ny y } will follow similar distributions. Moreover, if N pixels share the same sample color distribution, the union of the samples can be seen as an N times larger super-set following the underlying distribution. By simply averaging them the variance reduction is increased by a factor of N .

The cornerstone of the proposed approach is to find the most similar pixels to each given pixel by comparing their underlying sample color distributions. This is the object of the next section.

.. Distribution-driven pixel similarity

Consider the empirical distribution of the samples color at a given pixel. Figure . depicts this distribution for five different pixels on two different scenes, for samples generated by a Monte Carlo path-tracing algorithm. In the first example (top row) the three pixels were selected because their colors are extremely close. A quick visual inspection shows immediately that the samples of the two edge pixels follow roughly the same color distribution, and that this distribution is considerably different from the one of the third pixel. This example illustrates to what extent the information provided by the sample color distribution can help discriminate pixels of different nature, even when their pixels color are similar.

In the following, we denote by C x = {c 1 x , . . . , c nx x } the set of the color of samples cast from pixel x, and by h(x) the corresponding empirical color distribution. To measure pixel similarity we propose to use the binned empirical distributions as pixel descriptors. Since in general we deal with tri-stimulus color images, we can choose to build this descriptor either as a single histogram in the three-dimensional color space, or as three one-dimensional histograms (one per color channel).

Given the samples color C x and C y at pixels x and y, and their corresponding n b -binned distributions h(x) = (h 1 (x), . . . , h n b (x)) and h(y) = (h 1 (y), . . . , h n b (y)), the Chi-Square distance is given by where n x = k h k (x) and n y = k h k (y) are the total number of samples in each set.

d χ 2 (C x , C y ) = n b k=1 ny nx h k (x) -nx ny h k (y) 2 h k (x) + h k (y) , (.)



In order to take into account spatial coherence, the previous pixel-wise distance can be extended to patches of half-size w centered at x and y as follows,

d χ 2 (P x , P y ) = |t|≤w d χ 2 (C x+t , C y+t ). (.)
Since the order in which the samples are calculated is irrelevant, the sample color distribution appears as a natural and complete descriptor of the compared sets. There are different ways of measuring the similarity between two distributions depending on the data type. In the case of continuous data, the Cramer-von Mises [Anderson ; Anderson and Darling ], the Kolmogorov-Smirnov [Stephens ; Press et al. ] or the Kantorovich-Mallows-Monge-Wasserstein distances (also known as the Earth Mover's Distance [Rubner et al. ]) are all accepted ways to compare distributions. These three similarity measures are computed as L p distances between the two cumulative distributions (L ∞ , L 2 and L 1 respectively). For categorical data, the most popular measure to compare distributions is the χ 2 distance previously defined in (.).

By discretizing the data in a fixed number of histogram bins, the computational complexity of measuring the similarity between two data sets can be kept bounded and independent of the  number of samples. This is important since this similarity measure is evaluated a large number of times. Thus, the color space will be divided into fixed bins, and the χ 2 distance fits well to this form for the data. However, if an image is rendered with very few samples, one of the other two metrics would be preferable.

Since we are only interested in similar distributions, we set a threshold κ on the Chi-Square distance divided by the number of non-empty bins k,

κ = d min/k.
The normalization by the number of non-empty bins is necessary since only the bins carrying information should be considered in the comparison.

Remark: comparing pixel values versus comparing distributions State of the art image

denoising algorithms measure pixel similarity by comparing pixel colors. Indeed, the bilateral filter and NL-Means replace each noisy pixel by a weighted average of the most similar ones. In the case of nl-means, the pixel comparison is performed with patches centered around each pixel. For a very recent review on patch based denoising methods, we refer to [Lebrun et al. ] and for a fast implementation to [Adams et al. ].

Nevertheless, image denoising algorithms must know or measure the noise variance to evaluate properly the similarity of noisy samples. Fortunately, Monte Carlo rendering is an almost ideal situation where mean and variance values of the rays cast from each pixel can be estimated.

The main disadvantage of this formulation is that it cannot distinguish noise from intrinsic pixel variability. As a first example, suppose that a pixel is situated on an edge. In that case the sample color distribution will be at least bi-modal. Thus, it will probably have a large variance. This variance will result in a large tolerance to differences in the means, and consequently different pixel types may be wrongly mixed up. A case of this type is shown in Figure . (top row).

On the other hand, by directly comparing distributions, pixels lit from several sources can be better clustered. In the case of the histogram comparison, we will need no implicit nor explicit noise model assumption.

The bottom row of Figure . shows two pixels with very different pixel values. This is the consequence of the presence of a single very bright ray sample in one of the distributions. By comparing the ray color distributions, it is nevertheless possible to conclude that both pixels are from the same "nature", while this conclusion could not be reached by comparing the averages.



.. Distribution-driven average

For each pixel x, we define N ǫ (x) as the set of pixels y whose centered patches P y are such that d χ 2 (P x , P y ) ≤ ǫ. Then, if ǫ is such that these pixels are of the same nature as x, the maximum likelihood estimator of the noiseless pixel color is simply their arithmetic mean

ū(x) = 1 |N ǫ(x)| y∈Nǫ(x) ũ(y).
Unlike the previous estimator, where only the center of the patch is averaged, we can proceed to denoise the whole patch, and to denoise the image patchwise. Let us denote by V x the color values of a denoised patch centered at pixel x. Similarly, this denoised patch can be computed by averaging the patches which are at a Chi-square distance smaller than ǫ:

V x = 1 |N ǫ(x)| y∈Nǫ(x) ũ(P y ),
where we use the convention that ũ(P y ) is the evaluation of u on each pixel in patch P y .

In this way, by applying this aggregation procedure for all patches in the image, we shall dispose of (2w + 1) 2 possible estimates for each pixel. These estimates can be finally averaged at each pixel location in order to build the final denoised image.

ũ(x) = 1 (2w + 1) 2 |y-x|≤w V y (y -x).
This patchwise implementation is the one considered in this work.

.. Removing low-frequency noise

As already mentioned, in a pure Monte Carlo scenario the approximation error is characterized by a white random noise. This means that all frequencies are equally contaminated by noise. The proposed filtering procedure described so far filters noise at patch scale. Long wavelength noise cannot be eliminated by this procedure, because long wavelength structures cannot be captured by small patches. Removing noise at lower frequencies requires a (straightforward) multi-scale extension of the method. Let us define two useful operators, the s× Gaussian downsampling D s u(x) := (G 2 s σ * u)(2 s x) and U s the s× bicubic interpolator. Now, for each scale s, the corresponding histograms h s (x) have to be computed. Since each pixel at scale s results from the fusion of a set of neighboring pixels in the original finer scale, the new histograms are obtained by fusing the color histograms of all pixels in the same neighborhood. To obtain h s (x), the same down-sampling operator D s can be applied to the original color distribution h(x). Then, at each scale, the resulting histograms are re-normalized so that the sum of their areas is preserved across scales (thus preserving the original total number of samples in the finer scale).

Given a noisy image input ũ and its respective pixel color distribution h(x) the multi-scale histogram fusion proceeds as follows:

. Generate the Gaussian multi-scale sequence: ũ0 = ũ, ũs = D s ũ, s = 1, . . . , N , and their respective sample color distributions.

. Apply the denoising algorithm separately to each scale to recover ū0 , ū1 , . . . , ūN . The multi-scale approach eliminates low-frequency noise, as can be seen in the second and third row, and in the profile for a particular line shown in the bottom row.



. Implementation Details

As previously stated, our approach builds on two basic blocks: the estimation of the sample color distribution at each pixel, and a non-local multi-scale filtering based on averaging pixels sharing similar sample color distributions. This requires two kinds of data from the rendering system: the noisy Monte Carlo image ũ(x) and the associated sample color histograms h(x).

A fundamental aspect of our method is that sample color histograms can be computed on the fly, in parallel with the Monte Carlo rendering process. This is extremely important, since it makes the memory requirements independent from the number of rendered samples. The memory complexity bounds are fixed by the number of pixels, and therefore, as pointed out in the introduction, the input Monte Carlo images may consist of a large enough number of samples to produce high quality images.

This section gives the implementation details to estimate the sample color distribution and to perform the non-local multi-scale filtering.

.. Computing the color distribution of samples

To approximate a distribution using a histogram, one has to divide the range of possible values into discrete bins and count the number of samples within each bin. Smoother estimates can be obtained using kernel density estimation, by interpolating the contribution of each sample using a kernel. In this work, we used a triangular kernel to linearly interpolate the contribution of each sample color value to its adjacent bins.

We recall that pixel values are obtained by averaging sample color values. Hence, despite the saturation value for pixels (perfect white) is one, the range spanned by the sample values is much larger. In order to take into account the fact that high-energy (bright) samples are less frequent than low-energy ones, the bins are designed so that their sizes increase with the sample value, following an exponential law of exponent 2.2. The range covered by the histograms is set to [0, 7.5], and all samples exceeding this range are counted assigned to the last bin. It is worth mentioning that although histogram comparison is not particularly sensitive to these parameters, they must be chosen to cover the dynamic range adequately.

In general, sample values have a tri-stimulus color representation. Therefore we can either compute one three-dimensional distribution (3d) where bins are boxes in the full 3d color space, or compute three one dimensional distributions, one for each color. Although distributions in the 3d color space can capture inter-color correlations, a much larger number of bins are required to keep the same quantization step, and consequently a larger number of samples. In Section . we present a comparison of both strategies. This comparison shows that there is no advantage in using the full 3d color space.

The estimation of the color distribution can be done on the fly while samples are being computed. This keeps the memory bounded and independent of the number of samples per pixel produced by the rendering system.

.. Filtering

The implementation of the rhf filter is straightforward. In addition to the parameters needed to compute the histogram, four parameters are involved in the algorithm: the number of scales n s , half the patch size w, half the search block size b, and the χ 2 distance threshold.



The search of similar patches is restricted to a block of size (2b+1)×(2b+1). This is reasonable since the probability that two patches are similar will be smaller if one is distant from the other. The threshold is directly set on the Chi-square distance as a product of κ (the user parameter) by the number of bins k where both histograms are non-empty. A pseudo-code of both the filtering at each scale and the multi-scale generalization are presented in Algorithms  and , respectively. In algorithm , the denoised version of patch P i is obtained by averaging all patches Q j such that

d 2 χ (P i , Q j ) < k • κ.
Note that the only user parameter is κ. This parameter controls the amount of noise that is removed, or in other words the trade-off between image smoothness and noise. Its optimal choice depends mostly on the samples generation process (the considered renderer). The dependence on the rendered scene is actually very weak, as will be demonstrated by experiments in Section .. A simple intuitive explanation for the dependence of the optimal κ on the rendering method comes from the observation that the value of κ is related to the confidence associated to the color samples.

If the samples values are computed with low confidence, the distance threshold should be less restrictive. For instance, in pure Monte Carlo path tracing, each sample carries the energy of a single light path, while in volumetric ray tracing each sample value is computed as the average of several light paths. Therefore, the samples generated with pure path tracing have lower confidence, and this explains why the threshold should be less restrictive.

The practical implications of this fact is that, once a rendering method has been chosen, the value of κ can be safely fixed once for all. Moreover this tuning is not time consuming: indeed, since the distance between patches (the heaviest computational task) is independent of the parameter, its computation can be first performed and then several values of the parameter can be tested with practically no additional cost.

The multi-scale implementation in Algorithm , as detailed in Section .., sequentially decomposes the input noisy image at each scale, filters each scale using Algorithm  and reconstructs the multi-scaled filtered image.

.. Time complexity

The complexity of the filtering at each scale is O(N wbn b ) where N is the number of pixels. Note that the computational cost is independent of the number of samples.

In the case that two scales are used the computational cost increases by about %, the lowfrequency noise filtering being done on a four times smaller image. If n s scales are used the computational cost is bounded from above by % of the filtering time at the finest resolution.

. Experimental Set-up and Results

Different types of scenes containing complex geometries, indirect illumination, depth-of-field and other effects were rendered using the software pbrt-v2 [Pharr and Humphreys ]. The color distribution estimation stage was implemented on top of pbrt, so the color histograms were produced online as the samples were computed. The filtering-reconstruction stage was implemented in a stand-alone application which makes use of a multichannel image that to each pixel associates its corresponding sample color histogram, and the noisy Monte Carlo image generated with a box filter.

We compared the proposed algorithm to three different methods, both regarding image quality and execution time. The first one is a pure Monte Carlo rendering (mc): this is the basic approach  Algorithm : Single-Scale Ray Histogram Fusion Input: mc image ũ, corresponding histograms h, patch size w, search block size b , distance threshold κ Output: Filtered image ū : ū ← 0 : n ← 0 //auxiliary counter at each pixel in the image : for every pixel i do :

P i ← patch centered at pixel i : W i ← search block for pixel i : c ← 0 and V ← 0 : for every j ∈ W i do : Q j ← patch centered at pixel j : d ← ChiSquareDistance(h(P i ), h(Q j )) : k ← Number of non-empty bins in h(P i ) + h(Q j ) : if d < κ • k then : V ← V + ũ(Q j ) : c ← c + 1 : end if : end for : V ← V /c :
n(P i ) ← n(P i ) + 1 // +1 for each pixel in P i :

ū(P i ) ← ū(P i ) + Vū(P i ) ./n(P i ) : end for Notation convention: ũ(P i ) is the evaluation of ũ on each pixel in patch P i (the same applies for ū, n, h). The operator ./ (line ) represents element-wise division.

Algorithm : Ray Histogram Fusion

Input: mc image ũ, corresponding histograms h, patch size w, search block size b , distance threshold κ, number of scales n s . Output:

Filtered image ū = ū0 : s ← n s -1 : n T ← x,k h k (x)
// total number of samples : while s ≥ 0 do :

u s ← D s (ũ) : h s ← D s (h), n s T ← x,k h s k (x), h s ← n T n s T h s : ūs ← rhf (u s , h s , w, b, κ) : if s < n s -1 then : ūs ← ūs -U 1 D 1 ūs + U 1 ūold : end if :
ūold ← ūs :

s ← s -1 : end while  to generate photorealistic images. This technique is asymptotically unbiased but the variance shows slow (linear) decay with the number of samples.

The second algorithm chosen for comparison is an adaptation of the classic nl-means [Buades et al. ]. In image processing, nl-means performs denoising by averaging similar patches. In the rendering scenario, this method is obviously valid and can be improved by considering the noise level at each pixel, estimated from the variance of the samples that are cast from it. The main difference is in the way similar patches are identified. The performance comparison with nl-means will show that the knowledge of the sample color distribution adds a very significant amount of information, not yet contained in the patch colors.

Finally, we also consider comparison with the Adaptive Sampling and Reconstruction technique asr by Rousselle et al. []. As already discussed in Section ., this Monte Carlo based method can estimate the reconstruction error and control the number of samples cast from each pixel to reduce it. This methods is similar to ours in the sense that it does not rely on fat samples, and uses only the final color of each rendered sample. As such, the method scales well with the number of samples and can be used to produce high quality renderings of complex scenes. asr is a state of the art algorithm in this class of methods. Comparison is made using the code provided by the authors, and manually setting the parameters to produce the highest possible image quality, while matching the execution time of our algorithm (including both the samples rendering time and the filtering stage).

The success criterion is to get an image that is very close to the ground truth in a much shorter time. Image quality is assessed by comparing results to reference images, generated by pure Monte Carlo rendering with a very large number of samples per pixel. The performance measure is the standard peak-signal-to-noise ratio (psnr) calculated as psnr = 10 log 1 MSE where mse is the mean square error to the reference image. The psnr is a reliable criterion to characterize the quality of the reconstruction. It will nonetheless be complemented by some close-ups of difficult image details. All experiments were performed on a × Intel Xeon cpu X @ .GHz ( cores) with GB of RAM.

All the algorithms were run on several scenes from the pbrt software, simulating various effects with varying complexity levels.

In all cases three independent histograms were calculated, one for each channel (R, G, B) with nbins=. The search for similar patches was limited to a 13 × 13 window centered on the filtered pixel. The κ threshold (the user parameter) was manually set to produce a good balance between smoothness and remaining noise. As previously explained, the optimal value for this trade-off depends mostly on the rendering method. The values of κ that were chosen in the experiments are shown in Table .. Note that for all the renderings performed with a pure path tracing, we set κ = 1, while for the scene rendered using volumetric ray tracing (plants-dusk), κ = 0.37. This is consistent with the fact that the color samples generated with volumetric ray tracing result from an average of several light paths, and therefore are more precise than in pure path tracing.

A summary indicating all the considered effects, rendering method and image size is shown in Table ..

D versus 3× D Histograms. Table . illustrates the performance of the method as a function of the number of bins of the rays color histogram. The experiments do not support the use of three-dimensional color space bins. Thus, independent histograms were generated for the R, G, B channels. The number of bins must be large enough to capture the histogram structure, but not too large to grant a robust histogram comparison. Average of the k = 15 closest neighbors, the performance metric is the psnr with respect to a ground truth image. Two different ways of estimating the ray color histogram: bins in the 3d color space or three one-dimensional histograms one for each color. In most cases, taking more bins increases the performance, but also the computational cost of the algorithm. Estimating three independent histograms one for each color gives better results than estimating the histogram in the original 3d color space.



Robustness: Comparing Means versus Comparing Distributions.

Suppose that an external Oracle tells us the exact number of closest patches that should be averaged in each pixel to minimize the mse. By comparing the resulting psnr, we can compare the maximum theoretical performance of the method. A comparison for different patch sizes is given in Table .. The results show that using histogram information to compare patches permits to significantly reduce the patch size. This is an advantage since it permits the algorithm to be more local, hence better preserving small details. As for the psnr, the experiment shows that using the color distribution information leads to much better results than using just the pixels colors, as nl-means does.

Comparisons for several scenes. rhf systematically outperforms asr, as shown in Figures

.-.. Even if in some scenes both algorithms reach very similar psnrs, the proposed algorithm does not introduce artifacts while asr often fails to capture the geometry and causes spots. The psnr gain by rhf filtering is spectacularly larger than the one that would be obtained by generating more Monte Carlo samples using the same time span. Indeed, a db psnr gain by a pure mc  algorithm requires to double the number of samples. Instead, the filtering increases the psnr by  to  decibels. This amounts to decreasing the overall sampling time by a factor ranging from  to .

In the case of the yeahright image in Figures ., . and . the asr algorithm produces a slightly better psnr. This scene is best suited for this algorithm, because several regions are flat and asr can distribute most of the samples in the problematic parts. Nevertheless, in the shadows rhf produces a more natural smooth result. As previously commented, the nl-means based approach cannot distinguish between a large histogram variance due to pixel complexity, from a variance due to mc noise. This fact is well illustrated in the metal edge of Figure ., which is completely removed by nl-means, while it is well preserved by rhf.

The plants-dusk scene with participating media, in Figures ., . and . is a very challenging one. Here, the principal problem is the complex geometry of the vegetation. The proposed algorithm tends to blur and to slightly flatten some texture details. Nevertheless, contrarily to asr, no artifacts are introduced.

To illustrate the fact that the proposed denoising method is independent from the rendering system, in Figures . and . we present a filtering experiment that runs on a image generated by photon-mapping and final gathering. This scene comes from [Rousselle et al. ]. The noise has been properly removed, and no artifacts are observed.

Extension to animated sequences. The ideas behind this approach can be immediately extended to video sequences where pixels on every frame will be candidates in the search. The supplementary video shows the result of denoising an animated sequence, by implementing this simple generalization. Similar patches are searched within a temporal window of size  (namely in the previous, actual and next frame). Although no explicit temporal correlation is enforced, the filtered sequence does not show significant flicker. This is a consequence of the stability of the proposed multi-scale filter. .17: Results in plants-dusk scene (close-ups). In all cases, the PSNR values are given for the whole image. The comparison is done in such a way that the ASR computational time matches the Monte Carlo samples generation + RHF filtering time. This image presents a very complex fine geometry which is very difficult to capture with few samples. Nonetheless, the proposed algorithm produce acceptable quality and PSNR. For NL-means and RHF the indicated time follows the format total time (filtering time).











. Discussion, Limitations and Future Work

The maximum distance authorized between two patches plays an important role in the biasvariance tradeoff of the method. If the threshold is set in a conservative way, then very few pixels will be averaged. Thus, the filtering stage will not introduce bias, but the variance reduction will be low. On the other hand, if set too large many pixels of different nature would be considered similar, and averaged by error. Then the resulting image would be smooth but also biased (see Figure .).

If we accept that the selection of the most similar pixels for each noisy input pixel is independent of the number of samples, then the gain in psnr when casting more samples is only due to the averaging of less noisy pixels. By the randomness of the MC rendering, the noise of the input pixels is reduced by +db/octave, thus the ideal (best) slope should be +db/octave. This is the ideal, because it assumes that there is no error in the selection of similar pixels. Therefore, we can consider the difference in slope to the ideal +db/octave as a measure of experimental bias (introduced error). While the proposed algorithm rhf has an experimental bias of .db, the nl-means bias is three times larger. More important, the relative bias to the mc gain is 0.2/3 ≈ 0.07, which demonstrates that the proposed algorithm makes very few wrong ray color attributions.

Moreover, the proposed algorithm is consistent up to the discretization of the color distribution. As the number of samples increases, more evidence is required to average two pixels. In the limit two pixels will be averaged only if their color histograms are the same. Therefore, in practice, as the number of samples grows the method converges to the expected solution, as illustrated by the experiment in Figure ..

The acceleration factor depends on the degree of self-similarity of the scene, which fortunately is usually high [Lebrun et al. ]. Besides, in order to capture details, pixels need a large enough number of color samples to be well characterized. This is actually a design decision: we wanted our method to produce unbiased high quality images for any kind of scenes and complex effects, and this naturally requires a proper sampling of the light field. If this requirement is not met, the algorithm may not properly cluster similar pixels and details may be removed due to over-blur, as it happens with some details in the plants image. In the case of very low sample numbers, if the path-space is regular enough to be well described by sparse sampling, methods based on strong scene hypothesis that use fat samples [Sen and Darabi ; Lehtinen et al. ] are certainly much more adapted.

. Conclusion

In this work we have introduced rhf, an adaptive filtering scheme that accelerates Monte Carlo renderers. In the proposed approach, each pixel in the image is characterized by the collection of rays that reach its surface. The proposed filter uses a distance based on the sample color distribution of each pixel, to decide whether two pixels can share their samples. This permits to boost the performance of a Monte-Carlo render by reusing samples without introducing significant bias.

We have presented several experiments showing that rhf achieves artifact-free high quality noise reduction on a variety of scenes, and is able to cope with multiple simultaneous effects. The method is not only capable of removing high frequency noise: thanks to its natural multi-scale design, it can also successfully remove low-frequency noise. The proposed method can be easily extended to process animated sequences.

The method is independent of the rendering system and can be applied to samples generated by different methods, such as pure Monte Carlo path-tracing or photon-mapping with final gathering. It could also be potentially applied to post-process other methods that re-synthesize samples If κ is too small the test on the similarity is excessively conservative, and the noise is not reduced (high variance). If κ is large, too many pixels are averaged and the image is blurred (high bias). The results were calculated on the toasters scene generated with 256 spp.



using information from the scene, like the one recently proposed by Lehtinen et al. []. An advantage of the proposed filter is that its time and memory complexities do not depend on the number of input samples, and scale linearly with the image size. Finally, since a direct output of our method is the number of similar pixels for each given pixel, a decision on where to distribute new samples can be adopted. This may lead to an adaptive rendering version of the proposed filtering approach.

 

Conclusions and Perspectives

This dissertation concerns the digital image formation process and the quest for (i) mathematical models to better describe the kind of images that are obtained by this process, and

(ii) algorithms to improve image quality based on such models and to estimate the parameters of those models.

More precisely, we focus on the problems of denoising and blur estimation and address them for images that have been generated via two fundamentally different formation processes. In part I we deal with optically generated images, that have been produced by real photons after conversion to digital information by the sensors in a real camera. In contrast, part II is concerned with images that have been synthetically generated by computer simulations of virtual light hitting a digital 3d model and captured by a virtual camera.

We have introduced a mathematical model of an optical digital camera that provides a framework to evaluate the quality of psf estimation patterns. The key idea is to find the pattern that makes the psf estimation inverse problem as well-posed as possible. Or in other words, which is the pattern that carries the most information for the psf estimation problem. Indeed, the theoretical analysis concludes that a random pattern made from Bernoulli noise is nearly optimal.

An algorithm for estimating the psf from a photograph of this pattern is proposed. The method is very strict in the precautions to compensate for variations in illumination and possible geometrical distortions between the ideal pattern and the observed image. The procedure is successful, but its setup is somehow tedious since a calibration pattern must be printed and properly assembled. Very accurate estimates of the psf are achieved, as it is amply shown in the various comparative experiments (the relative error is of the order of 2% to 5%). As usual, a localityaccuracy trade-off had to be resolved. Locality of an order of a few hundred pixels can be achieved under common noise conditions.

Within the same mathematical framework, we confirmed that traditional state-of-the-art slantedge based methods lead to inverse problems that are ill-posed and require some kind of regularization. This is a major inconvenient, since introducing regularization considerably biases the estimation.

To avoid the use of a pattern, we introduced a semi-blind algorithm for the subpixel estimation of the camera psf that uses aliased photographs. The procedure is based on taking two frontoparallel photographs of the same flat textured scene, from different distances leading to different geometric scales, and then estimating the blur between them.

The psf estimation method is regularization-free, being closely related to the Bernoulli pattern algorithm. Nevertheless, the question is far more intricate since only the relative blur between the  acquisitions can be estimated. The algorithm is based on a mathematical analysis that proves that the psf can be recovered from the relative blur.

Images of planar textured scenes with a flat enough spectrum are necessary to reach high accuracy. Fortunately, as was experimentally shown, many textures found in nature are well adapted to these requirements. This semi-blind algorithm reaches similar accuracy levels to the one with the near-optimal Bernoulli pattern, with the advantage of not requiring any special acquisition setup or calibration pattern, thus being much more practical.

The experiments here suggest that for typical dslr cameras, each color channel is undersampled 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 subpixel estimation procedure.

Accurate psf estimates are important for several classical image processing applications, such as image superresolution, image de-blurring, and camera quality evaluation. In what follows we present some insights of possible future research directions regarding these and other applications.

In the current conditions, the camera psf acts as a strong low-pass filter attenuating significantly all frequencies beyond two times the camera's channel sampling frequency. From the point of view of the Shannon-Nyquist theory, this sets an upper limit in the gain of multi-image superresolution that is not too attractive.

In a modern dslr camera with no significant lens aberration, working at a wide aperture (e.g., f-number ≤ 5.6), the psf is mainly due to the optical anti-aliasing filter and the light integration in the sensor array. Hence, it becomes fundamental to remove this optical filter, that has been introduced to mitigate the aliasing artifacts, to unveil the latent superresolution potential that a digital camera may have. Some companies in the market, such as max-max  , offer to remove the camera optical anti-aliasing filter and to replace it by an ordinary glass with the same refractive index. Recently, Nikon has introduced a new model (d800e) of its d800, which is essentially the same but does not apply the optical anti-aliasing filtering. Although not a majority, some camera manufacturers are already considering the possibility of doing in-camera superresolution. Hasselblad developed a camera (h4d-200ms) capable of producing images of 4× the sensor resolution by fusing several snapshots with the sensor subpixically shifted.

The problem of image superresolution and deconvolution with a spatially variable psf are intertwined. Traonmilin et al. [] analyzed the problem of superresolution, and showed as an example how to make independent color channel demosaicking under the simplistic hypothesis that the deconvolution and superresolution are decoupled. However, as the authors claim, this is only accurate if the psf is space invariant and the the snapshots are purely translational motions. Having a tool to estimate locally and accurately the psf can unlock this problem and allow a multiimage camera to be used as a true measuring instrument without aliasing.

Plenoptic cameras, which sample the 4d light field, have proved to be useful for single image digital refocusing or depth estimation (e.g., the Lytro camera  ). However, this richer sampling comes at a significant associated cost in spatial resolution. Some researchers [Bishop et al. ; Bouman et al. ; Bishop and Favaro ]) have recently explored to extend the resolution of these cameras combining all the information captured from the light field. A mathematical modeling of the acquisition process, along with a local estimate of the psf is essential to achieve this goal.

Another related problem is the design of fast autofocus algorithms. It then becomes essential to make a very fast blur measure, especially in the case of autofocus in video, and to decide in which  www.maxmax.com  www.lytro.com  direction and how far the focal point should be moved to capture an in-focus image. In addition, a precise calibration of the out-of-focus psf could be very useful for the problem of estimating the depth from defocus. Currently, the in-focus kernels are oversimplistically modeled as Gaussian or circular kernels. However, it is well known by experimented photographers that these kernels, which give rise to the popularly known bokeh effect, have much more complex shapes. Codedaperture and coded-exposure imaging (see e.g., [Raskar et al. ; Levin et al. ]) allow the point spread function to be engineered for better conditionate specific applications such as motion deblurring or depth-from-defocus. Maximizing the quality of the restored image requires an accurate estimation of the resulting psf. We would like to explore and quantify the impact of the psf estimation in these problems.

Another useful application and a possible future research direction of having accurate estimate of the psf is to know the types of aberrations that an optical camera introduces. The proposed accurate non-parametric psf estimate can be matched to parametrical simulations of different types of aberrations to conclude which type of aberrations the camera presents (such as, comma, astigmatism, chromatic blur). This information can be helpful to quantitatively assess the quality of a device, but also to better design the optical system.

Detecting the parts of an image that were poorly sampled (leading to aliasing) and restore them is an important problem in image processing. The work of Coulange and Moisan [] proposed to detect, by analyzing the Fourier spectrum, those harmonics that are likely to be aliased. One possible avenue of research is trying to use the this technique along with the work of Almansa et al. [] for image restoration with samples in non-rectangular lattices to make single-image antialiasing for traditional and non-traditional light sensors.

Last, but not least, online psf estimation is of significant importance in the calibration of modern astronomical telescopes. Most of these telescopes have a large mirror (of several meters of diameter) which is exposed to variable stress due to unpredicted weather conditions. The only way to build and maintain in service such a large system is to work in an active optics principle: the mirror has actuators that modify its shape (and thus the wave front and the psf) proactively to prevent deformations. Image analysis and continuous psf re-estimation becomes necessary to adapt the optical shape and obtain a psf as sharp as possible.

The second part of this thesis addresses a denoising problem in image synthesis. Although there have been several breakthroughs very recently, synthesizing high quality realistic images in a reasonable amount of time remains a major challenge in computer graphics. In this dissertation, we introduced an adaptive filtering algorithm for accelerating Monte Carlo renderers. In the proposed approach, each image pixel is characterized by the set of rays that reach its surface. The algorithm, that we named ray histogram fusion (rhf), uses a similarity measure on the empirical ray color distribution of each pixel, to decide whether two pixels can be fused.

This simple procedure permits to boost the performance of any stochastic renderer by reusing samples without introducing significant bias. The proposed method achieves artifact-free high quality noise reduction on a variety of scenes, and is able to cope with multiple simultaneous rendering effects. The method is renderer independent and can be applied to samples generated by different methods, such as pure Monte Carlo path-tracing or photon-mapping with final gathering.

Thanks to its natural multi-scale design, it can successfully remove noise at all scales. Moreover, the algorithm is easily extended to process videos from animated sequences. One major advantage of the ray histogram fusion filter is that its time and memory complexities are independent of the number of input samples, and it scales linearly with the image size.



As a future work, we would like to investigate how rhf can be applied to post-process other methods that re-synthesize samples using information from the scene, like the one recently proposed by Lehtinen et al. []. Also, since a direct output of the method is the number of similar pixels that each pixel has, a decision on where to distribute new samples can be adopted. This may be the basis of an adaptive rendering version of the proposed filtering.

The proposed rhf filter assumes that the pixels grouped as similar have exactly the same expected value. In practice, this does not strictly hold and can therefore lead to the introduction of a small bias. We would like to explore whether the use of more general models (e.g., affine or more complex statistical models) can improve the performance and keep bias controlled.

Besides, it would also be interesting to explore different ways of reducing the computational cost of the rhf filter. Recently, Gastal and Oliveira [] have introduced a technique for accelerating filters based on the auto-similarity principle. Their algorithm, which reaches outstanding results, learns a set of manifolds that well capture the image structure, and then filters each of them separately. Hence, this is a natural research direction to reduce the computational cost of the rhf filter.



A Detailed Description of the PSF Estimation Algorithms

This section presents a complete detailed algorithmic description of the proposed non-blind and two-scale photographs psf estimation procedures. A demo facility and a reference source code can be found in the respective IPOL publications [Delbracio et al. a,b].

A. PSF Estimation from a Calibration Pattern Image

The captured image is precisely aligned to the analytic pattern by means of the surrounding checkerboard markers. Non-uniform illumination and non-linear camera response function impact -crf are corrected from the captured image to allow an artifact-free superresolved psf estimation.

In the next paragraphs we present a brief summary for each block.

Pattern detection

In order to detect the pattern in the image we use the line segment detector (lsd) algorithm by Grompone von Gioi et al. []. The idea is to detect the segments that are present in the pattern structure as can be seen in Figure A.. By using the detected line segments (in Figure A. are shown in red, green and blue for illustration purposes) we can approximate the checkerboard corners. These initial corners locations are re-adjusted by a subsequent subpixel stage. The pattern detection procedure can be summarized into the following steps:

. Detect all lines segments in the image. Each segment is represented by its two extreme points.

. For each detected segment s i do:

(a) Take another detected segment s j , j = i.

(b) Calculate s ij , the coordinates of segment s j in a new coordinate system relative to segment s i . This new coordinate system maps the first extreme point of segment s i to (0, 0) and the second one to (0, 1). The idea is to detect all the red and blue segments in relative position to the green one. i. Compute the distance between s ij and the ideal relative position of every searched segments (red and blue segments).

ii. If it is up to a tolerance and the distance is less than in the previous detection of the same segment, update the segment and the distance to the optimal position . If there is one and only one segment that has all associated segments we consider that the pattern has been detected.

Once the pattern with the aforementioned segments has been detected, an initial estimate of the positions of O, P, Q, R is computed as the middle points of segments ,, and .

Pattern subpixel alignment

In order to deal with geometric distortions the ideal pattern and its observation have to be precisely aligned locally. For that reason we have introduced checkerboard corners. Several methods to detect checkerboard corners have been reported in the Computer Vision literature ranging from differential operators such as Harris detector to more specific correlation methods. In this work we used a Harris-Stephens based corner detector in which we iteratively refine the detected corner positions to reach subpixel accuracy. The present algorithm is included in the OpenCV library [Bradski ]. The subpixel corner detector is based on the fact that the image gradient at a point close to the center is orthogonal to every vector from the center to that point.

 C

This can be mathematically expressed as

∇I t P • CP = 0
for every P in a neighborhood of C. As in practice the acquired image will be contaminated with noise, this expression will not be strictly zero. This last expression can be written as (∇I P • ∇I For all our experiments we set tol = 10 -5 , max iter = 200, R = 3.

Geometric transformation estimation

For our purpose of psf estimation we do not need to decompose the distortion into its homography and non-homography parts, as it is done in classical geometric camera calibration. Instead we use thin-plate splines [Sprengel et al. ] to model the whole deformation. Since we have previously detected the {p i } checkerboard corners centers from the observed pattern image and we know exactly their ideal corresponding locations {P i }, we can use these correspondences to find a smooth mapping from the non-distorted to the distorted space.



The thin-plates are found by minimizing the functional:

E = i f (P i ) -pi 2 + λ f 2 xx + 2f 2 xy + f 2 yy dxdy,
where λ controls the amount of regularization. As we will show if λ is zero the thin-plates interpolate the example points and if it goes to infinity the mapping becomes a pure affine transformation. The solution of this functional is of the form

f i (x, y) = d i0 + d i1 x + d i2 y + n j=1 c ij Φ( P j -(x, y) ) i = 1, 2,
where Φ(x) = x 2 log x , n is the number of points (in our case n is the number of checkerboard corners i.e., ). A nice property of the thin-plates is that it can always be decomposed into an affine and a non-affine component. The  coefficients (d ij ) form the affine part and the n×2 coefficients (c ij ) form the non-affine part. We now show how both set of coefficients can be obtained. In the following we represent points in homogeneous coordinates (x, y, 1). Let us call Y and X the n × 3 concatenated versions of the point coordinates P i and pi respectively. We denote by Φ the n × n matrix formed from Φ( P ipj ), i, j = 1..n, representing the thin-plate kernel.

It can be shown (see the work by Sprengel et al. []) that the solution is given by

ĉ = Q 2 (Q T 2 ΦQ 2 + λI (n-3) ) -1 Q T 2 Y d = R -1 Q T 1 (Y -Φĉ),
where Q 1 and Q 2 are obtained from the qr decomposition of matrix X,

X = Q 1 Q 2 R 0 .
In this application as we always have a fixed number of points -12 -and the support of the pattern in the image is very small (about 100 × 100 pixels), the geometrical distortion is expected to be minimal and thus very well approximated by an affinity. In practice, setting λ = 10 proved to be an appropriate trade-off.

Illumination estimation and normalization

In order to recover the psf, the gray levels in the sharp image pattern and those in the observed image have to be matched. For that purpose, the proposed pattern has black and white square regions to estimate the mapping between black and white colors and the respective observed gray values. As there are several of these constant regions located at different places, we can estimate a black (white) image that models the black (white) intensity level at each pixel. We model this black (white) image by a second order polynomial:

I black (x, y) = ax 2 + by 2 + cx + dx + e.

The estimation of the coefficients is done by least squares from the known pairs (value, position). Then, assume v(x, y) is the observed image, it can be corrected to get v c according to: v c (x, y) = v(x, y) -I black (x, y) / I white (x, y) -I black (x, y)



CRF estimation g(•)

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

The solution crf is defined as a parabolic function u → αu 2 + (1α)u where α is chosen so that the mean of the pattern after the correction is .. Hence, α can be calculated as:

α = 0.5 -u (u 2 -u) ,
where α is well defined as long as u is not a binary image which will never happen in our case.

Pattern rasterization

Rasterization is the procedure of converting a vector image into a raster pixel image. As we have a vector description of the pattern image we can rasterize it at any desired resolution. For that, it must be interpolated at the desired resolution by taking into account the estimated geometrical transformation. Then, we also need to cut the spectrum of the rasterized image to be band-limited at twice the desired frequency at which the psf is estimated. We do this by the following procedure:

. The continuous pattern u is sampled at a very high resolution, e.g. each of the flat black or white squares is sampled with a 4 × 4 block of pixels. We are going to work directly with these samples u (digital image), instead of the continuous pattern.

. Frequencies higher than sπ are cut off from the digital pattern u to get ũ. The dct of u is computed and the dct coefficients larger than (m × s, n × s) are set to zero (m and n are the number of rows and columns of the noise part in the captured image). Then we recover the filtered version ũ by applying the Inverse dct. Note that π represents the camera sensor sampling frequency that is why the factors m and n are included.

. Using the previously computed geometric distortion, the filtered pattern ũ is scaled to the desired resolution s by bicubic interpolation. The resultant image ũD is of size (s × (m -1) + 1, s × (n -1) + 1).

Building the linear system

In order to find the psf, we generate the following linear system: . The s-down-sampling matrix S s takes one sample per each block of s × s pixels.



. Finally a mask is applied that puts to zero all values that are outside the region of interest.

The region of interest consists of a trapezoid mask that restricts the convolution to the noise part of the pattern. Also this mask is eroded by a factor (max(p, q) -1)/2s to avoid boundary problems due to the convolution of finite support sequences.

We rewrite the observed image v in vector form to be consistent with the matrix formulation of the system. The mask restricting to the region of interest is also applied to the observation v.

Numerical methods for PSF estimation

Finally, we need to solve for a non-negative psf, arg min

h i S s C[ũ D ]h -v 2
2 subject to h i ≥ 0, i = 1, . . . , r 2 . (P)

To directly solve the non-negative least squares problem (P) we used the Newton interior point algorithm proposed by Portugal et al. [] and if we release the non-negative hypothesis the solution is simply found by a least squares algorithm. Another option is to threshold the least squares solution to be non-negative. The reference source code can operate in any of these three options.

A. PSF estimation from Two Photographs at Different Distances

The input of the algorithm are the two digital images: ṽ1 , ṽ2 , the superresolution factor s and the kernels (inter-image kernel and psf) support size: p × p at the s× superresolved grid. The output of the algorithm are a s× sampling of the inter-image kernel k and a s× sampling of the camera psf h. Both images are of size p × p.

Image subpixel alignment and geometric transformation estimation

In order to align both images and to estimate the geometric transformation from one to the other we use sift points and the orsa-Homography subroutine by Moisan et al. []. These subroutines may be replaced by any other accurate subpixel registration method. The important output of this stage is that given the two images we have a function D that maps one to the other. Suppose that the common parts of ṽ1 and ṽ2 are respectively of size m ′ × n ′ and m × n. Then

D : [0, m ′ -1] × [0, n ′ -1] → [0, m -1] × [0, n -1].
In the case of a homography, D can be expressed in homogeneous coordinates as linear transform represented by the 3 × 3 matrix:

D =   h 0,0 h 0,1 h 0,2 h 1,0 h 1,1 h 1,2 h 2,0 h 2,1 1   .



Image interpolation: H λ/s ṽ1

In order to generate the rescaled samples H λ/s ṽ1 from the digital image ṽ1 we need to interpolate it at the desired scale λ/s . This is done by using the estimated geometric transformation D.

From now on we will consider that 1× is the camera sampling frequency and the equivalent frequency band [-π, π] 2 .

The spectrum of the resampled image should be cut to be band-limited at [-sπ, sπ] 2 before resampling it. This is necessary to avoid aliasing artifacts. We do this by the following procedure: . Frequencies higher than sπ are cut-off from the zoom-in image ṽ1 . The dct of ṽ1 is computed and the dct coefficients larger than (m × s, n × s) are set to zero (m and n are the number of rows and columns of the common region part in the captured zoom-out image).

Then we recover the filtered version by applying the inverse dct. Note that 2π represents the camera sensor sampling frequency that is why the factors m and n are included.

. Using the previously computed geometric transformation D, the filtered zoom-in image is bicubically interpolated at the desired resolution λ/s. We do this by estimating the image values at a regular s× grid: [0, m ′ -1] × [0, n ′ -1] (i.e., the step size is /s). The resultant image H λ/s ṽ1 is of size (s × (m -1) + 1, s × (n -1) + 1) .

Solving for the inter-image kernel

First, the following linear system is built:

arg min k MS s C[H λ/s ṽ1 ]k -Mṽ 2 2 2 , (Pi)
where the matrix MS s C[H λ/s ṽ1 ] is composed by: . The matrix C[H λ/s ṽ1 ] associated to the D convolution operator with the interpolated image H λ/s ṽ1 . The convolution is done with a kernel of size p × p.

. The s-down-sampling matrix S s takes one sample per each block of s × s pixels.

. Finally a mask M is applied that puts to zero all values that are outside the region of interest.

The region of interest consists of a trapezoidal mask that restricts the convolution to the common part in the two images. Also this mask is eroded by a factor p-1 2s to avoid boundary problems due to the convolution of finite support sequences.

We rewrite the zoom-out observed image as a vector to be consistent with the matrix formulation of the system. The mask M restricting the image to the region of interest is also applied to the observation ṽ2 . Next, we need to solve Problem (Pi) which is simply a least squares algorithm.

From the inter-image kernel to the the PSF

In order to recover the camera psf h we need to compute (see Section .)

h = lim n→∞ H λ n (k * H 1 λ k * • • • * H 1 λ n k).
The value λ = (λ x , λ y ) is estimated from the geometric transformation D. In the case D is estimated as a homography, the scale values are taken to λ = (h 0,0 , h 1,1 ). This corresponds to the  situation where D is a pure zoom, and is a good approximation to the fronto-parallel acquisition with negligible rotation. We proceed as follows:

. Initialize u 0 = k, n = 1.

. Compute H 1/λ n k by using λ = (λ x , λ y ) (bicubic interpolation).

. Calculate u n = H 1/λ n k * u n-1 .

. If min{λ n x , λ n y } > λ max go to . Else update n := n + 1 and repeat from .

. Calculate h = H λ n u n (bicubic interpolation).

The algorithm converges after a few iterations since λ n grows very fast. We set λ max = 50 , since the convolution with the inter-image kernel zoomed-out 50× or greater produces a negligible change in the final result.

Since negative light does not exist the estimated psf should be positive. We can therefore constrain the solution to be non-negative by projecting the result of step  to the non-negative half-space.

General tips for the set-up

. The scene should be as planar and as textured as possible.

. The photographs should be taken between  and  relative distance. The possible superresolution factor is always less than the relative distance, so for 4× estimation relative distance should be higher than .

. To produce accurate estimations it is highly recommended to use a tripod to avoid handheld shake.

. Both photographs should be taken with the same camera parameters. The only exception is camera focus, that should be re-set to have both images in focus.

. The images should be taken with the same illumination conditions.

. Both images should be recorded in raw format (no compression, no post-processing: no demosaicking, no denoising, no enhancing, etc). raw conversion is camera-dependent and not provided by our demo. In our examples a suitable conversion of the raw format to a pgm image containing the Bayer pattern could be achieved by the command "dcraw - -D input.raw output.pgm"  . A single color channel should then be extracted and use as the input to our algorithm.

 dcraw is the Dave Coffin's utility for decoding raw digital photos in Linux.



  Several patterns have been used for psf estimation, ranging from pin-hole, slanted-edge [ISO ; Reichenbach et al. ; Zandhuis et al. ; Claxton and Staunton ], or arc-step-edge patterns [Joshi et al. ] to random noise images [Daniels et al. ; Levy et al. ; Backman et al. , ; Brauers et al. ].

Figure 1 . 1 :

 11 Figure 1.1: Different calibration patterns for local PSF estimation. On left the Joshi et al. pattern consisting of 120 • arc step edges. On the right the proposed noise pattern consists of 256 × 256 small black/white squares drawn independently from a Bernoulli equiprobable distribution. The lateral checkerboard-like marks and black/white squares are introduced for alignment and illumination estimation purposes.

Figure 1 . 2 :

 12 Figure 1.2: Example of an acquired image of the Bernoulli pattern. On the right, a crop of the central part containing the pattern. As is shown on the left, the pattern should cover only a small region of the image (roughly 100 × 100 pixels). This allows for a local PSF estimation.

  Figure 1.3: Horizontal profile of the PSF Fourier spectrum modulus for different state-of-the-art methods: Joshi et al. [2008] , Imatest [LLC 2010] commercial software and the proposed Bernoulli pattern algorithm. 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. The Imatest software produces a quite noisy estimation.

Figure 1

 1 Figure 1.4: The result of estimating a 4× PSF for one of the green channels at different image locations.For each pixel at the sensor resolution we estimate 4 × 4 samples of the PSF. The 1× resolution (camera grid) is shown in green dotted lines. Kernels closer to image borders are larger and more asymmetrical than the kernel at the image center. This is mainly due to lens aberration. Although no regularization is imposed, the resulting kernels are smooth.

  Figure 1.5: An example of a pair of digital images that allow to estimate the PSF. Top: two distant, parallel views of a wall. Bottom: the inter-image kernel between these two views. The inter-image kernel models the necessary blur that should be applied to the closest image to produce the farthest image (with the necessary zooming). The camera PSF is recovered from the inter-image kernel through an iterative procedure. The estimated inter-image kernel and camera PSF are obtained at 4× the camera resolution for the blue channel. Although no regularization is imposed, the kernels are smooth.

  Figure 1.6: Comparison of the PSF estimation from the two-scaled photographs and the non-blind estimation using the Bernoulli pattern. The estimations are for the same camera and configuration at 4× the resolution. Both estimations are significantly close and the difference is mostly due to noise present in the two-scale estimation.

  green

Figure

  Figure 1.8: Example of an image rendered with Monte Carlo path-tracing. In a pure MC scenario the square error decreases linearly with the number of samples per pixel (spp), thus the convergence is quite slow.

  Choudhury and Tumblin ; Xu and Pattanaik ]. More complex filters make use of ray information available from the renderer in order to filter also at pixel level [Rushmeier and Ward ; McCool ; Dammertz et al. ; Xu et al. ]. The most sophisticated filters, use the additional ray information to adaptively filter the sample rays [Shirley et al. ; Sen and Darabi ; Rousselle et al. ; Lehtinen et al. ].

Figure 1 .

 1 Figure 1.10: A running example of the ray histogram fusion algorithm. The filter increased the PSNR of this Monte Carlo rendered image by +11.6 decibels. To get an equal PSNR with pure Monte Carlo, 15× more samples would have been needed.

Figure 2

 2 Figure 2.1: Typical Bayer pattern. Image taken from Wikipedia [2012].

Figure 2 . 2 :

 22 Figure 2.2: 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 the plot on the left, as a function of M . The plot at right shows their difference, also as a function of M .

Figure 2 . 3 :

 23 Figure 2.3: Comparison of the lower bound given by Lemma 3 and the gamma obtained from the expected S s C[u] operator when using a Bernoulli(0.5) random noise pattern. Both bounds are shown in the plot on the left, as a function of M whereas the plot at right shows their difference, also as a function of M .

  Figure 2.4: Reaching theoretical bounds. A random Bernoulli binary image is used to generate the S s C[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.

  Figure .. The resulting ũD ij for distances in a range from s/t = 1 produce γ values close to the γ bound for entries in [0, 1].

  Figure2.5: 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 zoomin. 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 . shows the proposed random pattern, compared to a pattern designed by Joshi et al. consisting of 120 • arc step edges.

Fig

  Fig. . shows the eigenvalues of the S s C[ũ D ] 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 Bernoulli pattern. The random pattern produces secondary eigenvalues very similar in contrast to the fast decay shown by the eigenvalues of the slanted-edge Joshi pattern. The γ values for the corresponding patterns are shown in Table .. In all cases, the random pattern significantly outperforms the Joshi et. al. pattern. The γ bound value was computed by taking into account the effective observed window size, that is, leaving out the auxiliary region with the checkerboard and flat regions.

Figure

  Figure 2.6: Calibration patterns for local PSF estimation

  Fig. .. The next paragraphs present brief summaries for each block. A detailed description is given in Appendix A.

Figure 2

 2 Figure 2.8: 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 impact -CRF are corrected from the captured image to allow an artifact-free s× PSF estimation.

  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 .. The estimations were performed at the image center from a window of size 90 × 90 pixels. Notice that kernels at apertures f/ and f/ 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 .. Fig. 

Figure 2

 2 Figure 2.9: 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.

Figure 2

 2 Figure 2.10: 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 axis-symmetric. 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.

Figure 2 .

 2 Figure 2.11: Real camera example. Taken at aperture f/5.6. 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).

  Figure2.12: Diffraction-limited System. Theoretical diffraction MTF for monochromatic green light with circular f/5.6 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.

  Figure 2.14: Dependence on the CRF correction, for a 4×PSF estimation of the green channel, at aperture f/5.6. 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.

  Figure2.15: 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 moduli. Bottom row: MTF horizontal and vertical profiles. The estimation was performed using the random pattern captured at aperture f/5.6. The red PSF is larger than the green and 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 note 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. . caption.

  Figure2.18: Non-negative Constraint. This experiment analyzes how the PSF estimation changes by not assuming the non-negative hypothesis. On the left we show both estimations: the no-constrained and the non-negative 4× PSF for the green channel, f/5.6 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.

Figure 2 Figure 2

 22 Figure2.19: 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 .

  al. ], superresolution [Park et al. ; Robinson and Milanfar ] or shape from defocus [Chaudhuri and Rajagopalan ].

  we have presented in the previous chapter, several patterns have been used for psf estimation, ranging from pin-hole, slanted-edge [ISO ; Reichenbach et al. ; Zandhuis et al. ; Claxton and Staunton ], or arc-step-edge patterns Joshi et al. [] to random noise images [Daniels et al. ; Levy et al. ; Backman et al. , ; Brauers et al. ].

  Two pictures of another textured wall shown in Figure . were used to estimate the psf of the four color Bayer channels (raw camera output). This wall texture presents characteristics similar to those of white noise. The results for the 4× psf estimated at the image center are shown in Figure ..

Figure 3 . 3 :Figure 3 . 4 :

 3334 Figure 3.3: Synthetic example: 4× PSF estimation for simulated data. Top row: the closest and farthest images.Middle row: the simulated PSF (ground truth) and the respective PSF estimations using the automatic SIFT points / homography alignment and the ideal alignment. Both estimations are accurate. However, as shown in the difference images the automatic registration introduces a small misalignment. This can also be seen in the phase and modulus of the PSF Fourier transform vertical profile, shown in the bottom row. Bottom row (right): comparison of the inter-image and PSF kernels. Since both input images are simulated at distances in a ratio of λ = 4×, h is very close to k.

  Figure3.5: Wall image: two-scale versus Bernoulli white noise pattern based estimation. Estimation at 4× PSF resolution for one of the green channels from the camera RAW output. Top row: two distant, parallel views of a textured wall. Middle row: the PSF estimated with the proposed algorithm and the one estimated using the Bernoulli pattern method. Bottom row: vertical profile of the MTF. Both estimations are close. In particular the associated airy disks have similar sizes, and the MTFs vanish in the same locations.

Figure 3 . 8 :

 38 Figure 3.8: Magazine image: 1×, 2×, 3× and 4× estimations for the first green channel from a pair of photographs of a newspaper image. The estimation was done at the image center for the camera working at an f/5.6 aperture. All the estimations are consistent: their MTFs show good overlap. The 4× PSF estimation is noisier than the one produced from the wall images. The main reason is that the spectrum of the magazine image decays faster.

Figure

  Figure . shows the mtf profiles of the obtained psf estimates. The proposed two-scale method performs at least as well as the non-blind methods under comparison. The method of Joshi and colleagues shows similar performance for a carefully, manually chosen regularization parameter. See caption for details.

Figure 3

 3 Figure3.10: Comparison of PSF/MTF estimation methods. Our implementation of the PSF estimation algorithm byJoshi et al. [2008], Imatest commercial software LLC[2010], our previous Bernoulli pattern method, and the two-scale method proposed in this work (applied to the images of the wall shown in Fig.3.7). On the low frequencies all algorithms produced very similar estimates, while on the higher frequencies the Joshi et al. estimation depends strongly on the regularization level. Although much effort was made 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. The proposed two-scale algorithm is the one presenting an estimation closest to the non-blind estimation presented in Chapter 2, considered as ground truth by virtue of its high accuracy.

Figure 4 . 1 :

 41 Figure 4.1: The top row singles out three pixels in Cornell Box scene and their sample color distributions. (The samples with color values falling out of the [0, 1] 3 -box are by convention colored in red.) The first pixel, situated on the brown wall, has a unimodal sample color distribution. The other two pixels belong to an occlusion boundary showing a bimodal green-brown distribution. This feature is shared by many pixels on the same boundary, which can therefore share their samples. The bottom row shows two pixels of the toasters scene with different colors. Their sample color distributions are nevertheless very similar and will therefore be merged as well.

.

  Compute the final image ū = û0 by the recursion ûi = ūi -U 1 D 1 ūi + U 1 ûi+1 initialized with ûN = ūN for i = N .

Figure

  Figure . shows the importance of dealing with noise at multiple scales. When filtering only at a single fine scale, conspicuous low frequency noise remains. This noise is almost completely eliminated by the multi-scale procedure with three scales.

Figure 4 . 3 :

 43 Figure 4.3: Results in sibenik scene. The result of RHF is shown on the left while the noisy input is shown on the right. For close-ups on difficult parts see Figures 4.4 and 4.5.

Figure 4 . 4 :

 44 Figure 4.4: Results in sibenik scene (close-ups). In all cases, the PSNR values are given for the whole image. In general ASR tends to create artifacts near edges. RHF produces the best PSNR, with no visible artifacts. For NL-means and RHF the indicated time follows the format total time (filtering time).

Figure 4 . 5 :

 45 Figure 4.5: Results in sibenik scene (close-ups). In all cases, the PSNR values are given for the whole image. In general ASR tends to create artifacts near edges. RHF produces the best PSNR, with no visible artifacts. For NL-means and RHF the indicated time follows the format total time (filtering time).



Figure 4

 4 Figure 4.6: Results in cornell-box scene. The result of RHF is shown on the left while the noisy input is shown on the right. For close-ups on difficult parts see Figures 4.7 and 4.8.

Figure 4 . 9 :

 49 Figure 4.9: Results in toasters scene. The result of RHF is shown on the left while the noisy input is shown on the right. For close-ups on difficult parts see Figures 4.10and 4.11.



Figure 4 .

 4 Figure 4.12: Results in yeahright scene. Fine geometry details, glossy surfaces and indirect illumination presented in the yeahright scene are rendered with the PBRT-V2 path-tracing algorithm. The result of RHF on the left hand part of the image is shown on the left while the noisy input of the right part is shown on the right. For close-ups on difficult parts see Figures 4.13and 4.14.

Figure 4 .

 4 Figure 4.15: Results in plants-dusk scene. The result of RHF is shown on the left while the noisy input is shown on the right. For close-ups on difficult parts see Figures 4.16 and 4.17.

Figure 4 .

 4 Figure 4.18: Light interaction with participating media rendered through a photon mapping algorithm shows the generality of the proposed filtering. The comparison is done in such a way that the ASR computational time matches the Monte Carlo samples generation + RHF filtering time. For close-ups on difficult part see Figure 4.19.

Figure 4 .

 4 Figure 4.19: Results in dragon-fog scene (close-ups). The comparison is done in such a way that the ASR computational time matches the Monte Carlo samples generation + RHF filtering time. Since Photon Mapping is a biased rendering algorithm no reference image (and therefore no PSNR) was computed.

Figure 4

 4 Figure4.20: Changing the distance threshold κ: it fixes the maximum distance that two color distributions can differ. (a): a small detail in the toasters image filtered with the RHF algorithm with growing κ values. The MSE presents a minimum for κ = 0.7 -1.0 (b): If κ is too small the test on the similarity is excessively conservative, and the noise is not reduced (high variance). If κ is large, too many pixels are averaged and the image is blurred (high bias). The results were calculated on the toasters scene generated with 256 spp.

Figure 4

 4 Figure 4.21: Top: A comparison of the PSNR for RHF, NL-means, and the pure MC pathtracing algorithm on the toasters scene.As the number of samples per pixel increases the PSNR increases. In a pure MC scenario the square error decreases linearly with the number of samples (which is a trivial consequence of averaging independent samples). Thus, duplicating the number of samples produces a 3db gain. Although for the proposed algorithm the slope gain is a little smaller (2.8db/octave), RHF reduces the error significantly in comparison to MC. +15.5dB is a huge difference; it permits to reach the same image quality with 35× fewer samples. On the bottom we show a close up, generated with a varying number of samples through pure MC (top) and filtered with RHF (bottom).

  Figure A.1: The Bernoulli pattern and the line segments used for the detection of the pattern in an input image (shown in red, green and blue).

  t P ) • (P -C) = 0 and since this expression holds for P in a neighborhood of C,P ∈N (C) (∇I P • ∇I t P ) • P = P ∈N (C) (∇I P • ∇I t P ) • CIn order to give more importance to points closer to the center C, a weighting functionW P = f ( P -C ) is included.Then, C is found by least squares idea, the checkerboard detector algorithm iteratively runs as follows:. Given C (with possible subpixel precision), compute N (C) as a square of size 2R+1×2R+ 1 centered in C. The non-integer pixel values are calculated using bilinear interpolation.. The gradient ∇I P is calculated by finite differences in every point P of N (C).. ∇I P • ∇I t P W P is computed whereW P = f ( P -C ) is an isotropic Gaussian function centered in C. The variance of the Gaussian function is σ 2 W P = (2R+1) 2 2 .. Set C old = C and compute the new C by solving (A.).. If C -C old < tol or ++iter > max iter exit, else go to .

  matrix S s C[ũ D ] is composed by: . The matrix C[ũ D ] associated to the D convolution operator with the image ũD . The convolution is done with a kernel of size p × q.

  Delbracio, P. Musé, A. Almansa and J.M. Morel. The Non-parametric Sub-pixel Local Point Spread Function Estimation Is a Well Posed Problem. International Journal of Computer Vision, vol. , no. , pp. -, January .

• M. Delbracio, A. Almansa, J.M. Morel and P. Musé. Subpixel Point Spread Function Estimation from Two Photographs at Different Distances. SIAM Journal on Imaging Sciences, vol. , no. , pp. -, November .  Part I Recovering the Camera Point Spread Function  2 PSF Estimation from a Calibration Pattern Image

  Förstner ; Capel ; Smith ]. However, in real images, a step-edge convolved with the psf kernel is generally not a good model of the observed edges as noted by Ladjal [, chapter ]. Other blind approaches try to estimate the psf based on statistical models of sharp images [Chalmond ; Rooms et al. ; Zhang and Cham ; Šroubek et al. ]. 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.

Table 2

 2

	.1: PSF estimation algorithm summary.

  While this may not be new (noise  patterns have been used in non-blind psf estimation, see e.g. [Daniels et al. ; Levy et al. ; Backman et al. , ; Brauers et al. ]

Table 3 .

 3 1: Summary of the notation used in this chapter. A more precise definition of each term is presented in Appendix 3.A

  The complete chain is summarized in the block diagram of Figure .. The next paragraphs present brief summaries for each block. A complete algorithmic description is given in Appendix A.. Then the interpolated image H λ s ṽ1 and image ṽ2 are compared to obtain the inter-image kernel k, which is later iteratively updated to obtain the absolute camera PSF h.

Figure 3.2: Algorithm description. Both captured images are aligned via SIFT feature matching followed by the estimation of a homography. The relative geometric transformation and gray-level corrections are applied to a low-pass (unaliased) version of the finest scale image ṽ1 .

Table 4 .

 4 1: Summary of the tested scenes. The scenes cornell-box, plants-dusk and sibenik are from [Pharr and Humphreys ] while the scenes dragons and toasters were taken from [Rousselle et al. ].

	Scene	Effects	Size		Generation	κ-rhf
	cornell-box toasters plants-dusk ALPD AI AILD sibenik AILD 1024 × 1024 256 × 256 512 × 512 800 × 400 yeahright AI 800 × 800 dragons AILDP 512 × 512	path tracing path tracing ray tracing path tracing path tracing photon mapping + final gathering	. . . . . .
	Considered effects: anti-aliasing (A); indirect illumination (I); area lights (L); depth-
	of-field (D); participating media (P). Table 4.2: Performance comparison: estimating the histogram.
			3d			3 × 1D
		3 3	4 3	5 3	3×5 3×10 3×15 3×20 3×25
	cornell-box . . .	. . . . .
	toasters	. . .	. . . . .
	plants-dusk . . .	. . . . .
	yeahright . . .	. . . . .
	sibenik	. . .	. . . . .

Table 4 .

 4 3: Oracle performance comparison.For each pixel the ideal number of closest patches to minimize the error with respect to the ground truth image has been computed and fixed. The table compares nl-means and rhf for different patch sizes. It shows that rhf permits to reduce significantly the patch size.

		nl-means	rhf
		1×1 3×3 5×5	1×1 3×3 5×5
	cornell-box . . .	. . .
	toasters	. . .	. . .
	plants-dusk . . .	. . .
	yeahright	. . .	. . .
	sibenik	. . .	. . .

Acknowledgments

First and foremost I would like to express my heartfelt gratitude to my advisors Andrés Almansa, Jean-Michel Morel and Pablo Musé, with whom I have had the tremendous privilege to work for the past four years. The enormous willingness of Andrés to discuss ideas, his creativity and scientific rigor have profoundly marked this thesis. His unfailing good mood and constant optimism have made this work very enjoyable. It is hard to imagine a more generous advisor than Jean-Michel, from whom I not only have learned a lot of technical insights but also a research methodology based on mathematical thinking. His wise advice during our frequent meetings allowed me to work with total freedom in a very comfortable environment. My dear friend Pablo, who guided me and continuously transmitted his personal experience and humanism over these years. His continuous and total support along with his scientific advice and knowledge allowed me to always have confidence in our work.