Oliver Whyte
email: oliver.whyte@ens.fr

Josef Sivic

Andrew Zisserman

Jean Ponce

Non-uniform Deblurring for Shaken Images

Keywords: Motion blur, Blind deconvolution, Camera shake, Non-uniform / spatially-varying blur

Photographs taken in low-light conditions are often blurry as a result of camera shake, i.e. a motion of the camera while its shutter is open. Most existing deblurring methods model the observed blurry image as the convolution of a sharp image with a uniform blur kernel. However, we show that blur from camera shake is in general mostly due to the 3D rotation of the camera, resulting in a blur that can be significantly non-uniform across the image. We propose a new parametrized geometric model of the blurring process in terms of the rotational motion of the camera during exposure. This model is able to capture non-uniform blur in an image due to camera shake using a single global descriptor, and can be substituted into existing deblurring algorithms with only small modifications. To demonstrate its effectiveness, we apply this model to two deblurring problems; first, the case where a single blurry image is available, for which we examine both an approximate marginalization approach and a maximum a posteriori approach, and second, the case where a sharp but noisy image of the scene is available in addition to the blurry image. We show that our approach makes it possible to model and remove a wider class of blurs than previous approaches, including uniform

Introduction

Everybody is familiar with camera shake, since the resulting blur spoils many photos taken in low-light conditions. While significant progress has been made recently towards removing this blur from images, almost all current approaches to this problem model the blurred image as the convolution of a sharp image with a spatially uniform filter [START_REF] Chan | Total variation blind deconvolution[END_REF][START_REF] Fergus | Removing camera shake from a single photograph[END_REF][START_REF] Shan | Rotational motion deblurring of a rigid object from a single image[END_REF]Yuan et al. 2007a). However, real camera shake, which we show can be (mostly) attributed to the rotation of the camera during exposure, does not in general cause uniform blur, as illustrated by Figure 1.

In this paper we propose a geometrically motivated model of non-uniform image blur due to camera shake. By showing that such blur can be mainly attributed to the rotation (as opposed to translation) of the camera during exposure, we develop a global descriptor for such parametrically nonuniform blur, derived from the geometry of camera rotations about a fixed center. Our descriptor can be seen as a generalization of a convolution kernel, and as such our model includes uniform blur as a special case. We demonstrate the effectiveness of our model by using it to replace the uniform blur model in three existing approaches to camera shake removal, and show quantitative and qualitative improvements in the results.

Specifically, we consider the problems of "blind" deblurring, where only a single blurry image is available, and the case where an additional sharp but noisy image of the same scene is available. To approach these two problems, we apply our model within the frameworks proposed by [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF] and [START_REF] Fergus | Removing camera shake from a single photograph[END_REF], and by [START_REF] Cho | Fast motion deblurring[END_REF] for the blind case, and Yuan et al. (2007a) for the case of a noisy / blurry image pair.

Related Work

The problem of modeling non-uniform blur is not new, and previous approaches to this problem are diverse, as are the many possible causes of such blurs. Much of the previous work has relied on the local uniformity of the blur, for example, modeling blur due to moving objects as piecewise uniform [START_REF] Levin | Blind motion deblurring using image statistics[END_REF][START_REF] Cho | Removing non-uniform motion blur from images[END_REF][START_REF] Chakrabarti | Analyzing spatiallyvarying blur[END_REF], or approximating a continuously varying blur by a spatially varying combination of localized uniform blurs [START_REF] Nagy | Restoring images degraded by spatially variant blur[END_REF][START_REF] Vio | Multiple image deblurring with spatially variant PSFs[END_REF][START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF]Tai et al. 2010a). Models of non-uniform blur that do not rely on assumptions of local uniformity have been applied under various constrained motion models [START_REF] Shan | Rotational motion deblurring of a rigid object from a single image[END_REF][START_REF] Sawchuk | Space-variant image restoration by coordinate transformations[END_REF][START_REF] Klein | A single-frame visual gyroscope[END_REF]Tai et al. 2010b), and while (as in this work) global models are used to describe these continuously varying blurs, the constraints are often restrictive.

Recent work has investigated the automatic estimation of global descriptors for the non-uniform blur caused by camera shake. [START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF] use inertial measurement sensors to estimate the motion of the camera over the course of the exposure. This information can then be used to effectively deblur the captured image. [START_REF] Gupta | Single image deblurring using motion density functions[END_REF] propose a model similar to ours, where the blur-causing motions are approximated using image-plane translations and rotations, as opposed to the 3D camera rotations used in this work. [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF] note that the assumption of uniform blur made by most algorithms is often violated, but do not address this fact.

If the point spread function (PSF) for each pixel in an image is known, the problem of recovering a sharp image is generally referred to as non-blind deconvolution, for which standard techniques such as the Wiener filter (for uniform blur) or the Richardson-Lucy algorithm (for general blur) exist [START_REF] Richardson | Bayesian-based iterative method of image restoration[END_REF][START_REF] Lucy | An iterative technique for the rectification of observed distributions[END_REF][START_REF] Banham | Digital image restoration[END_REF][START_REF] Puetter | Digital image reconstruction: Deblurring and denoising[END_REF]. Sophisticated algorithms for non-blind deconvolution have recently been proposed [START_REF] Dabov | Image restoration by sparse 3D transform-domain collaborative filtering[END_REF][START_REF] Shan | High-quality motion deblurring from a single image[END_REF][START_REF] Yuan | Progressive inter-scale and intra-scale non-blind image deconvolution[END_REF][START_REF] Couzinie-Devy | Dictionary learning for deblurring and digital zoom[END_REF]), but their application has generally been limited to the case of uniform blur. Note however that [START_REF] Tai | Richardson-Lucy deblurring for scenes under projective motion path[END_REF] propose a modified version of the Richardson-Lucy algorithm for deblurring scenes under general projective motion, where the temporal sequence of projective transformations which caused the blur is known.

The task of recovering a sharp image when the pixels' PSFs are unknown, so-called blind deblurring, is a difficult problem. Existing approaches for uniform blur, where a single PSF, or "blur kernel", describes the blur everywhere typically proceed by first estimating this kernel, then applying a non-blind deconvolution algorithm to estimate the sharp image. For uniform blur, [START_REF] Fergus | Removing camera shake from a single photograph[END_REF] estimate the kernel by first applying the variational algorithm of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF] to approximate the posterior for the kernel and sharp image with a simpler distribution. This distribution is chosen such that it is then trivial to estimate the kernel by marginalizing over all possible sharp images. Many different maximum a posteriori (MAP) formulations have also been proposed, such as those of [START_REF] Shan | High-quality motion deblurring from a single image[END_REF], [START_REF] Cho | Fast motion deblurring[END_REF], [START_REF] Cai | Blind motion deblurring from a single image using sparse approximation[END_REF], and [START_REF] Xu | Two-phase kernel estimation for robust motion deblurring[END_REF]. These algorithms typically use an alternation scheme, updating the estimate of the blur kernel at one step, and of the sharp image at the next. The algorithm proposed by [START_REF] Gupta | Single image deblurring using motion density functions[END_REF] for their non-uniform blur model also follows this paradigm. To simplify the deblurring problem, others have considered using additional information in the form of additional blurry images [START_REF] Rav-Acha | Two motion-blurred images are better than one[END_REF][START_REF] Chen | Robust dual motion deblurring[END_REF], or a sharp but noisy image of the same scene (Yuan et al. 2007a;[START_REF] Lim | Estimation and removal of motion blur by capturing two images with different exposures[END_REF]. The rest of this paper is organized as follows. Section 2 presents our geometric model. Section 3 presents a discrete version of this model. In Section 4, we demonstrate its application within two existing algorithms for deblurring a single blurry image, both an approximate marginalization approach and a maximum a posteriori (MAP) approach, and in Section 5 we compare the results obtained these two algorithms. In Section 6 we examine a second deblurring problem, where a sharp but noisy image of the same scene is available, in addition to the blurry image. In Section 7 we describe some of the implementation details for the presented algorithms, and in Section 8 we conclude with a discussion of some limitations of our work and potential future research directions.

Geometric Model

To motivate our approach, let us begin by noting that the blur in a "shaken" image is caused by the motion of the camera during the exposure, i.e. changes in the pose of the camera. The pose of a camera can be split into two components: position and orientation, and in this section, we argue that in most cases of camera shake, the changes in orientation (rotations) of the camera during exposure have a significantly larger effect than the changes in position (translations). Consider the simplified case shown in Figure 2 of a scene point P , at a distance D from the camera, being imaged at the center of the camera's retina. During the exposure the image of the point is blurred through a distance δ pixels, either by (a) translating the camera through a distance X parallel to the image plane, or (b) rotating the camera through an angle θ We consider the scene point P at a distance D from the camera, whose image is blurred by δ pixels as a result of either of the two motions. In most cases, for a given blur size δ the rotation θ constitutes a significantly smaller motion of the photographer's hands than the translation X (see text for details)

about its optical center. By simple trigonometry, we can see that in the case of translation the camera must move by

X = δ F D, (1)
where F is the camera's focal length, while for the rotation, the camera must move through an angle

θ = tan -1 δ F . ( 2
)
If we make the common assumption that the camera's focal length F is approximately equal to the width of the sensor, say 1000 pixels, then to cause a blur of δ = 10 pixels by translating the camera, we can see from (1) that X = 1 100 D. Thus the required translation grows with the subject's distance from the camera, and for a subject just 1 metre away, we must move the camera by X = 1 cm to cause the blur. When photographing a subject 30 metres away, such as a large landmark, we would have to move the camera by 30 cm! To cause the same amount of blur by rotating the camera, on the other hand, we can see from (2) that we would need to rotate the camera by θ = tan -1 1 100 ≈ 0.6 • , independent of the subject's distance from the camera. To put this in terms of the motion of the photographer's hands, then for example if the camera body is 10 cm wide, such a rotation could be caused by moving one hand just 1 mm forwards or backwards relative to the other. Provided the subject is more than 1 metre from the camera, this motion is at least an order of magnitude smaller than for a translation of the camera.

In reality, both the position and orientation of the camera vary simultaneously during the exposure. However, if we assume that the camera only undergoes small changes in position (translations), then following the discussion above, we can assert that the variations in the camera's orientation (rotations) are the only significant cause of blur. We do this from now on, and ignore the translational component of camera motion. We consider all rotations to occur about the camera's optical center, and although this may not be the case, we note that rotations about a different center can be written as rotations about the optical center, plus translations.

Motion Blur and Homographies

Assuming that the scene being photographed is static, it is well known that rotations of a camera about its optical center induce projective transformations of the image being observed, assuming a pinhole camera model. That is to say that, excluding boundary effects, the image at one camera orientation is related to the image at any other by a 2D projective transformation, or homography. For an uncalibrated camera, this is a general 8-parameter homography, but for a camera with known internal parameters, the homography H is given by

H = KRK -1 , (3)
where the 3 × 3 matrices R and K are, respectively, a rotation matrix describing the motion of the camera, and the camera's internal calibration matrix [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF].

The matrix R has only 3 parameters. We adopt here the "angle-axis" parameterization, in which a rotation is described by the angle θ moved about an axis a (a unit-norm 3-vector). This can be summarized in a single vector θ = θa = (θ X , θ Y , θ Z ) ⊤ . R is then given by the matrix exponential

R θ = e [θ]× , where (4)
[θ] × =   0 -θ Z θ Y θ Z 0 -θ X -θ Y θ X 0   .
(5)

We fix our 3D coordinate frame to have its origin at the camera's optical center. The axes are aligned with the camera's initial orientation, such that the XY -plane is aligned with the camera sensor's coordinate frame and the Z-axis is parallel to the camera's optical axis, as shown in Figure 3 (a). In this configuration, θ X describes the "pitch" of the camera, θ Y the "yaw", and θ Z the "roll", or in-plane rotation, of the camera.

In this work, we assume that the calibration matrix K is known and takes the standard form

K =   F 0 x 0 0 F y 0 0 0 1   . (6)
This corresponds to a camera whose sensor has square pixels, and whose optical axis intersects the sensor at (x 0 , y 0 ), referred to as the principal point. Section 2.3 describes how we estimate K in practice.

Having defined the type of image transformation we expect, we now assume that when the shutter of the camera opens, there is a sharp image f : R 2 → R of a static scene that we would like to capture. The camera's sensor accumulates photons while the shutter is open, and outputs an observed image g : R 2 → R. In the ideal pinhole case, each point on the sensor sees a single scene point throughout the exposure, giving us a sharp image. However if, while the shutter is open, the camera undergoes a sequence of rotations R t , the sensor is exposed to a sequence of projectively transformed versions of the sharp image f . For each point on the sensor (a point in the observed blurry image g), denoted by the homogeneous vector x, we can trace the sequence of points x ′ t in the ideal sharp image f which were visible there during the exposure:

x ′ t ∼ H t x, (7)
where H t is the homography induced by the rotation R t , and ∼ denotes equality up to scale. The observed image g is then the integral over the exposure time T of all the projectively-transformed versions of f , plus some observation noise ε:

g(x) = T 0 f H t x dt + ε, (8)
where, with a slight abuse of notation, we use g(x) to denote the value of g at the 2D image point represented by the homogeneous vector x, and similarly for f . According to this model, the apparent motion of scene points may vary significantly across the image. Figure 3 demonstrates this, showing the paths followed by points in an image under rotations about either the Y or Z axis of the camera. Under the (in-plane) Z-axis rotation, the paths vary significantly across the image. Under the (out-of-plane) rotation about the Y -axis, the paths, while varying considerably less, are still non-uniform. It should be noted that the degree of non-uniformity of this out-of-plane motion is dependent on the focal length of the camera, decreasing as the focal length increases. However, it is typical for consumer cameras to have focal lengths of the same order as their sensor width, as is the case in Figure 3. In addition, it is common for camera shake to include an in-plane rotational motion. Fig. 3 Our coordinate frame with respect to initial camera orientation, and the paths followed by image points under single-axis rotations.

We define our coordinate frame (a) to have its origin at the camera's optical center, with the X and Y axes aligned with those of the camera's sensor, and the Z axis parallel to the camera's optical axis. Under single-axis rotations of the camera, for example about its Y -axis (b), or its Z-axis (c), the paths traced by points in the image are visibly curved and non-uniform across the image. This non-uniformity remains true for general camera shakes (d), which do not follow such simple single-axis rotations, but rather take arbitrary paths through camera pose space. The focal length of the camera in this simulation is equal to the width of the image, the principal point is at the image's center, and the pixels are assumed to be square

From this, it is clear that modeling camera shake as a convolution with a spatially invariant kernel is insufficient to fully describe its effects (see also Figure 1).

In general, a blurry image has no temporal information associated with it, so it is convenient to replace the temporal integral in (8) by a weighted integral over a set of camera orientations:

g(x) = f H θ x w(θ) dθ + ε, (9)
where the weight function w(θ) encodes the camera's trajectory in a time-agnostic fashion. The weight will be zero everywhere except along the camera's trajectory, while the value of the function along that trajectory corresponds (inversely) to the camera's rotational speed, i.e. if the camera moves slowly through a certain orientation, the weight will be large there, and vice versa.

Uniform Blur as a Special Case

One consequence of our model for camera shake is that it includes uniform blur as a special case, and thus gives the conditions under which a uniform blur model is applicable.

From the definition of the matrix exponential, e A = I + A + 1 2! A 2 + . . ., we can see that if θ Z = 0 and θ X , θ Y are small, Equation (4) can be approximated by discarding the 2 nd and higher order terms:

R θ ≈   1 0 θ Y 0 1 -θ X -θ Y θ X 1   . ( 10
)
Combining this with Equations (3) and ( 6), it can be shown that as F → ∞,

H θ →   1 0 F θ Y 0 1 -F θ X 0 0 1   , (11)
which simply amounts to a translation in the image plane of (F θ Y , -F θ X ) ⊤ . Noting that for typical camera shakes, θ X and θ Y will indeed be small, we can see that if the focal length of the camera is large and there is no in-plane rotation, a uniform blur model may be sufficient to describe the blur.

Camera Calibration

In order to compute the homography in Equation (3) that is induced by a particular rotation of the camera, we need to know the camera's calibration matrix K, as given by Equation (6). To estimate K, we recover the pixel size and focal length of the camera from the image's EXIF tags, and assume that the principal point is at the center of the image.

Note that this assumes the image has not been cropped or resized beforehand, as these operations will generally invalidate the EXIF information.

The radial distortion present in many consumer-grade digital cameras can represent a significant deviation from the pinhole camera model. Rather than incorporating the distortion explicitly into our model, we pre-process images with the commercially available PTLENS tool1 , which uses a database of lens and camera parameters to correct for the distortion.

A second distortion present in many digital images comes from the fact that the pixel values stored in, for example, a jpeg file, do not correspond linearly to the scene radiance. Most cameras apply a compression curve before storing the values, sometimes referred to as "gamma correction". Where possible we avoid this problem by using raw camera output images, such that the pixel values correspond linearly to scene radiance. In other cases, where the compression curve is known (e.g. having been calibrated), we preprocess the blurry images with the inverse of this curve to recover the linear values, and where it is unknown, we apply a generic sRGB curve.

Restoration Model

So far, our model has been defined in terms of the continuous functions f and g, and the weight function w. Real cameras are equipped with a discrete set of pixels, and output an observed blurry image g ∈ R N , where N = H × W pixels for an image with H rows and W columns. We consider g to be generated by a sharp image f ∈ R N and a set of weights

w ∈ R K , whose size K = N X × N Y × N Z is controlled
by the number of rotation steps about each axis that we consider. The set of weights w forms a global descriptor for the camera shake blur in an image, and by analogy with convolutional blur, we refer to w as the blur kernel. Figure 1 (b) shows a visualization of w, where the cuboidal volume of size N X × N Y × N Z is shown, with the yellow points inside representing the non-zero elements of w in 3D. The kernel has also been projected onto the 3 back faces of the cuboid to aid visualization, with white corresponding to a large value, and black corresponding to zero.

Each element w k corresponds to a camera orientation θ k , and consequently to a homography H k , so that in the discrete setting, the blurry image g is modeled as a weighted sum of a set of projectively transformed versions of f :

g = k w k C k f + ε, ( 12
)
where C k is the N × N matrix which applies homography H k to the sharp image f . The matrix C k is very sparse. For example, if bilinear interpolation is used when transforming the image, each row has only 4 non-zero elements. Expanding Equation ( 12), we obtain the discrete analog of Equation (9):

g i = k w k j C ijk f j + ε, (13)
where i and j index the pixels of the observed image and the sharp image, respectively. Appendix B describes how to calculate the coefficients C ijk . For an observed pixel g i with coordinate vector x i , the sum j C ijk f j interpolates the point H k x i in the sharp image. Figure 4 shows an example of this, where a blurry pixel g 5 , with homogeneous coordinate vector x 5 , is mapped under a homography H k to the point H k x 5 in the sharp image. The value of f at the point H k x 5 is then interpolated as a weighted sum of the pixels of f . Due to the bilinear form of Equation ( 13), note that when either the blur kernel or the sharp image is known, the blurry

g 1 g 2 g 3 g 4 g 5 g 6 g 7 g 8 g 9 x 5 Blurry image g H k → f 4 f 5 f 6 f 7 f 8 f 9 f 10 f 11 f 12 H k x 5 Sharp image f
Fig. 4 Interpolation of sub-pixel locations in the sharp image. In general, a homography H k will not map a pixel (e.g. x 5 ) in the blurry image g to a single pixel in the sharp image x. Instead, the value of f at the point H k x 5 is interpolated as a weighted sum of nearby pixels. Using bilinear interpolation, the value of f at H k x 5 will be interpolated from the pixels f 7 , f 8 , f 10 , and f 11 image is linear in the remaining unknowns, i.e.

g = Af + ε, or (14)
g = Bw + ε, (15)
where

A ij = k C ijk w k , and B ik = j C ijk f j .
In the first form, A ∈ R N ×N is a large sparse matrix, whose rows each contain a local blur filter acting on f to generate a blurry pixel. In the second form, when the sharp image is known, each column of B ∈ R N ×K contains a projectively transformed copy of the sharp image. We will use each of these forms in the following.

Comparison to Other Non-uniform Blur Models

Concurrently with our proposal of this model for camera shake blur, several other authors have proposed global models of non-uniform blur. In common with our model, they generally model the blurry image as a sum of transformed versions of the sharp image. [START_REF] Tai | Richardson-Lucy deblurring for scenes under a projective motion path[END_REF] model the blur process using a temporally-ordered sequence of homographies, which is known in advance, for example by using additional hardware attached to the camera. [START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF] also recover a temporally ordered sequence of homographies by using inertial measurement sensors attached to the camera to recover the camera's path through 6D (rotation and translation) pose space, and assuming that the scene can be modelled as a single fronto-parallel plane (although the authors note that the blur is generally independent of depth for objects further than 1 metre from the camera). [START_REF] Gupta | Single image deblurring using motion density functions[END_REF] propose a model which is similar in spirit to our own, recovering a set of weights over a 3D parameter space which describes transformations of the sharp image. However, they consider image plane translations and rotations, rather than the camera pose-induced homographies used in this work.

Application to Existing Deblurring Algorithms

The fact that Equation ( 13) is bilinear in the sharp image and blur kernel is the key feature that allows our model to be applied within existing deblurring algorithms previously applied only to uniform blur. Since convolution is also a bilinear operation on the sharp image and the blur kernel, it can often be replaced with the general bilinear form in Equation (13) without significant modification to the algorithm. The remainder of the paper demonstrates this, first in Sections 4 and 5 for the problem of single-image deblurring, and second in Section 6 for the case where a sharp but noisy image of the scene is also available.

It should be noted that an important case where our model cannot easily be substituted in place of convolution is when an algorithm relies on the ability to work in the frequency domain. When taking the Fourier transform, convolution becomes an element-wise multiplication of the frequency components of the image and kernel, however this is not the case for the more general bilinear form in our model.

The aim of deblurring algorithms is to recover an estimate f of the true sharp image f . Generally, the approach taken is to also estimate the blur kernel ŵ such that together, f and ŵ are able to accurately reconstruct the observed blurry image g. We denote this reconstruction as ĝ( f , ŵ), where, under our model, a blurry pixel is reconstructed as:

ĝi = k ŵk j C ijk fj . ( 16
)
The problem of finding the sharp image and blur kernel that best reconstruct the observed image is in general ill-posed, since we have fewer equations than parameters. In fact, for a given g, there are an infinite number of ( f , ŵ) pairs that can reconstruct g equally well. To obtain a useful solution, it is thus necessary to add some regularization and/or constraints on both the sharp image and the kernel. Successful algorithms for deblurring camera shake generally share the same two pieces of prior information about the blur kernel being estimated, which we mention here. First, all its elements are non-negative, since the image formation process is additive, with sensor elements accumulating photons. This constraint is equally applicable to our model, since each kernel element w k corresponds to a camera orientation, so that if that if the camera passed through orientation θ k during the exposure, w k will be positive, and if not, w k = 0.

The second, and arguably more important fact to observe about a blur kernel for camera shake is that it should be sparse, i.e. contain relatively few non-zero elements. This sparsity prior has been a prominent feature of previous camera shake removal algorithms, and has also been leveraged for the alignment of blurred / non-blurred images (Yuan et al. 2007b). [START_REF] Fergus | Removing camera shake from a single photograph[END_REF] encourage sparsity by placing a mixture of exponentials prior on the kernel values, while [START_REF] Cho | Fast motion deblurring[END_REF] and Yuan et al. (2007a) proceed by thresholding the estimated kernel values such that most of the kernel is set to zero. In a contrasting approach, [START_REF] Cai | Blind motion deblurring from a single image using sparse approximation[END_REF] choose to construct the blur kernel as a linear combination of a predefined set of "curvelets", and place the sparsity prior on the coefficients of the curvelets, rather than on the kernel elements directly. This sparsity prior is intuitively applicable to blur kernels for our model too, since the camera follows a path θ(t) through the space of camera orientations, and thus will only pass through a small subset of all possible orientations while the shutter is open.

Many different image priors/regularizers have been proposed for image reconstruction tasks such as deblurring, denoising, and super-resolution, often based on the statistics of natural images. The most commonly used priors for deblurring are those which encourage the image's response to a set of derivative filters to follow heavy-tailed distributions, e.g. [START_REF] Fergus | Removing camera shake from a single photograph[END_REF][START_REF] Shan | High-quality motion deblurring from a single image[END_REF][START_REF] Krishnan | Fast image deconvolution using hyperlaplacian priors[END_REF], which have been shown to be effective at suppressing noise while preserving sharp edges in the reconstructed images. When substituting our blur model into the algorithms in Sections 4 and 6, we use the image regularizers suggested in the original works. In order to compare the final results of the two methods for single-image deblurring, shown in Section 5, we use the algorithm of [START_REF] Krishnan | Fast image deconvolution using hyperlaplacian priors[END_REF], adapted for our non-uniform blur model.

Single-Image Deblurring

In this section, we examine the case where we have only a single blurry input image g from which to estimate f . We substitute our model into two successful algorithms for uniform blur, allowing them to handle non-uniform blur: those of [START_REF] Fergus | Removing camera shake from a single photograph[END_REF] and [START_REF] Cho | Fast motion deblurring[END_REF]. As discussed in the previous section, good priors for f and w are necessary for blind deblurring to be successful, so both approaches take the posterior distribution for f and w as their starting point:

p(f , w|g) ∝ p(g|f , w)p(f )p(w), (17)
where the likelihood is derived from an assumption of isotropic Gaussian noise:

p(g|f , w) ∝ i exp - (ĝ i (f , w) -g i ) 2 2σ 2 , ( 18
)
where σ is the standard deviation of the noise, and the individual priors used by the two algorithms will be discussed in the following sections. Both of the algorithms are mainly concerned with estimating the blur kernel ŵ, and after the termination of this process, a final non-blind image reconstruction step is performed using the estimate of ŵ to produce the deblurred output f . To estimate the kernel, Fergus et al. ( 2006) use the variational inference approach of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF] to perform approximate marginalization of the posterior over f , while [START_REF] Cho | Fast motion deblurring[END_REF] use alternating optimizations to maximize the posterior over both f and w.

The Marginalization Approach

In this section we adapt the algorithm proposed by [START_REF] Fergus | Removing camera shake from a single photograph[END_REF] for blind deconvolution of a single image. The algorithm is based on the variational inference approach of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF], originally designed for simultaneous deblurring and source separation of cartoon images. We show that the convolutional blur model in the original algorithm can be replaced with our non-uniform blur model, leading to new update equations for the optimization process, and we show in Section 5 that doing so improves the deblurred results.

The algorithm proposed by [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF] attempts to approximate the posterior distribution for both the kernel and the sharp image p(f , w|g) by a simpler, factorized distribution using a variational method. The factorized form of this distribution means that it is straightforward to marginalize over the sharp image in order to produce an estimate ŵ of the kernel. [START_REF] Fergus | Removing camera shake from a single photograph[END_REF] successfully adapted this algorithm to the removal of camera shake blur from photographs by applying it in the gradient domain, within a multiscale framework. They use a prior on the kernel which assumes that each element w k is independent, and follows a mixture of exponential distributions with mixture weights π d and decay parameters λ d :

p(w) = k D d=1 π d exp(-λ d w k ). ( 19
)
By working in the gradient domain, the latent variable f j for the intensity of a pixel is replaced by the x and y derivatives f x j and f y j at that pixel, which are treated as separate variables. Fergus et al. use a prior on the sharp image which assumes that the derivatives for all pixels are independent and follow a mixture of zero-mean Gaussians with mixture weights π c and variances v c :

p(f x ) = j C c=1 π c exp - f x j 2 2v c , (20)
and likewise for the y derivatives. For simplicity, within the context of this algorithm, we use f to denote the concatenation of the derivative images f x and f y , and use j to index over this, i.e. j ∈ {1, . . . , 2N }. Fergus et al. learn the parameters π d , λ d , π c and v c from real data, and we use the values provided by them directly. Finally, to free the user from manually tuning the noise variance σ 2 , the inverse variance β σ = σ -2 is also considered as a latent variable. Following [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF], we collect the latent variables f , w, and β σ into an "ensemble" Θ. The aim is to find the factorized distribution

q(Θ) = q(β σ ) j q(f j ) k q(w k ) (21)
that best approximates the true posterior p(Θ|g), by minimizing the following cost function (Miskin and MacKay 2000, Eqn. ( 10)) over both the form and the parameters of q(Θ):

C KL = q(Θ) ln q(Θ) p(Θ) -ln p(g|Θ) dΘ. ( 22
)
Minimizing this cost function is equivalent to minimizing the Kullback-Leibler divergence between the posterior and the approximating distribution [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF], and this is tackled by first using the calculus of variations to derive the optimal forms of q(f j ), q(w k ) and q(β σ ), then iteratively optimizing their parameters. For our blur model, the optimal q(Θ) has the same form as in [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF]. However the equations for the optimal parameter values differ significantly and we have calculated these afresh (the derivation is provided in the supplementary material).

For our non-uniform blur model, we find the following optimal values for the parameters, cf . (Miskin and MacKay 2000, Eqns. 46-49):

w (2) k = β σ i j C ijk f j 2 q(f ) (23) w (1) k w
(2)

k = β σ i g i j C ijk f j q(fj ) (24)
-

k ′ =k j C ijk f j j C ijk ′ f j q(f ) w k ′ q(w k ′ ) f (2) j = β σ i k C ijk w k 2 q(w) (25)
f (1) j f (2) j = β σ i g i k C ijk w k q(w k ) (26) - j ′ =j f j ′ q(f j ′ ) k C ij ′ k w k k C ijk w k q(w)
, where w

(1)

k and w

(2)

k are the parameters of q(w k ), f

(1) j and f

(2) j are the parameters of q(f j ), q(f ) = j q(f j ), q(w) = k q(w k ), and • q represents the expectation with respect to the distribution q. Note that in this context, the latent image pixels f j , kernel elements w k , and noise precision β σ are random variables. Note also that these equations cannot be implemented directly in this form, as they involve expectations over combinations of random variables. However, they may be rewritten in terms of the mean and variance of individual variables, and we provide these expanded versions in Appendix A. Having found the optimal q(Θ), the expectation of q(w) is taken to be the optimal blur kernel, i.e., ŵ = w q(w) . Fergus et al. choose to discard the latent image distribution q(f ), although as shown in Figure 8 (d), this may in fact provide a useful estimate of the sharp image.

The Maximum A Posteriori Approach

Cho and [START_REF] Cho | Fast motion deblurring[END_REF] proposed an effective single image deblurring algorithm, optimized for speed on uniform blurs. Again, we show that this algorithm can be readily adapted to handle non-uniform blur, substituting our model in place of convolution. The algorithm can be considered to perform alternating maximum a posteriori estimation of the blur kernel, using Gaussian priors on the kernel elements and on the latent image gradients. Simply performing an alternating optimization of f and w using these priors would almost certainly not produce any reasonable result, however the introduction of non-linear filtering and thresholding steps into the process encourages the algorithm to find a latent image with sparse gradients and a blur kernel with sparse non-zero elements, such as discussed in Section 3.2. The algorithm proceeds by iterating over three main steps, of which we give a brief overview here.

The first step takes the current estimate f of the sharp image and aims to predict strong edges, which are useful for the kernel estimation step, by applying a bilateral filter [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF] followed by a shock filter (Osher and Rudin 1990). The derivatives of this filtered image are computed, then thresholded to produce sparse gradient maps {p x , p y } which contain only the most salient edges. The threshold is chosen so as to retain only a small number of non-zero gradients, while ensuring that all orientations are well-represented.

In the second step, the gradient maps are used to estimate the blur kernel by minimizing the energy function

E w (w) = (p * ,g * ) ω * ĝ(p * , w) -g * 2 2 + β w 2 2 , ( 27
)
where the weights ω * ∈ {ω 1 , ω 2 } weight each partial derivative, (p * , g * ) varies among (p x , D x g), (p y , D y g), (D x p x , D xx g), (D y p y , D yy g), ( 1 2 (D x p y + D y p x ), D xy g) , and β is the regularization weight. Since ĝ(p * , w) is linear in w, this is simply a linear least squares problem, which can be solved efficiently using a conjugate gradient method. Having found the kernel that minimizes ( 27), the values are thresholded, such that any element whose value is smaller than 1 20 the largest element's value is set to zero. This encourages sparsity in the kernel, and ensures that all the elements are positive, as discussed in Section 3.2.

In the third step, the current estimate of the blur kernel ŵ is used to deconvolve the blurry image and obtain an improved estimate of the sharp image. This is performed by minimizing the energy function

E f (f ) = D * ω * ĝ(D * f , ŵ) -D * g 2 2 + α ∇f 2 2 , ( 28
)
where α is the regularization weight, and now, the partial derivatives include the zeroth order: D * ∈ {I, D x , D y , D xx , D xy , D yy }, where I is the identity, and ω * ∈ {ω 0 , ω 1 , ω 2 }. The use of the partial derivatives D * in the data terms of ( 27) and ( 28), as suggested by [START_REF] Shan | High-quality motion deblurring from a single image[END_REF], has the effect of improving the conditioning and regularizing the solutions. These steps are applied iteratively, working from coarse to fine in a multi-scale framework. The iterative process generally converges quickly at each scale, and 7 iterations are typically sufficient. For the parameters ω 0 , ω 1 , ω 2 , α, and β, we use the values given by [START_REF] Cho | Fast motion deblurring[END_REF]. Although we are not able to take full advantage of the speed optimizations proposed by Cho and Lee, due to their use of Fourier transforms to compute convolutions, the algorithm is generally able to estimate a blur kernel in a much shorter time than the marginalization algorithm of Section 4.1.

Modification for Non-uniform Blur

When applying our model within this algorithm, we must take into account some important differences between our 3D kernels and 2D convolution kernels. First, we note that the point spread function (PSF) of a single pixel does not uniquely determine the full 3D kernel, i.e. for every PSF there are many different kernels that could explain it. This can be seen by considering a vertical blur at the left or righthand side of the image. Such a blur could be explained either by a rotation of the camera about its X axis, a rotation about its Z axis, or some combination of the two. Thus we must ensure that the pixels used to estimate the kernel (the nonzeroes in {p x , p y }) do not only come from a small region of the image, in order for the kernel estimation step to be wellconditioned. To achieve this, we simply subdivide the image into 3 × 3 regions, and apply the gradient thresholding step independently on each. This ensures that we retain a set of gradients that are well distributed over both orientation and location.

A second observation is that our 3D kernels contain a certain degree of redundancy, arising largely from the inplane rotation of the camera. As can be seen in Figure 3, a rotation of the camera about its Z axis causes a very small displacement for pixels towards the center of the image. Thus, in the kernel estimation step, the information provided by these pixels will be ambiguous with respect to this component of the camera's motion. Only pixels near the edge of the image will be able to provide detailed information concerning this motion. While the spatial binning mentioned above goes some way to ensuring that these pixels from the edge of the image are present in {p x , p y }, they may be greatly outnumbered by pixels from the interior. As a result, the kernels recovered by minimizing (27) with our model generally contain many non-zeros spread smoothly throughout, and do not produce good deblurred outputs (see Figure 9).

If instead of the ℓ 2 regularization in ( 27), we apply ℓ 1 regularization combined with non-negativity constraints, the optimization is encouraged to find a sparse kernel and is more likely to choose between ambiguous camera orientations, as opposed to spreading non-zero values across all of them. This type of ℓ 1 kernel regularization was previously applied by [START_REF] Shan | High-quality motion deblurring from a single image[END_REF] for uniform blur. In our case, the energy function becomes

E w (w) = (p * ,g * ) ω * ĝ(p * , w) -g * 2 2 + β k w k s.t. ∀k = 1, . . . , K, w k ≥ 0. ( 29
)
This is an instance of the lasso problem [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], for which efficient optimization algorithms exist [START_REF] Efron | Least angle regression[END_REF][START_REF] Kim | An interiorpoint method for large-scale ℓ 1 -regularized least squares[END_REF][START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF]). The different results obtained using ℓ 2 and ℓ 1 regularization are discussed in Section 5. With the use of the ℓ 1 regularization, we found that the best results were obtained with a lower value of β than that given by Cho and Lee, and for the results in this paper using ℓ 1 regularization, we set β = 0.1. In the remainder of the paper, we refer to the original algorithm of Cho and Lee as MAP-ℓ 2 , and our ℓ 1 -regularized version as MAP-ℓ 1 .

Image Reconstruction

Having estimated the blur kernel ŵ for the blurry image, we wish to invert Equation ( 14) in order to estimate the sharp image f . This process is often referred to as deconvolution, and while many algorithms exist for this process [START_REF] Banham | Digital image restoration[END_REF][START_REF] Puetter | Digital image reconstruction: Deblurring and denoising[END_REF][START_REF] Dabov | Image restoration by sparse 3D transform-domain collaborative filtering[END_REF], they are often applicable only to uniform blur, since they typically rely on convolutions or the ability to work in the Fourier domain.

For the results of single-image deblurring in Section 5, we have adapted the deconvolution algorithm of [START_REF] Krishnan | Fast image deconvolution using hyperlaplacian priors[END_REF], which performs MAP estimation of the sharp image using a hyper-Laplacian prior on the image gradients. Specifically, it attempts to maximize the following posterior over f :

p(f |g, w) ∝ p(g|f , w)p(f ), (30)
with the prior

p(f ) = j exp(-λ|f x j | p ) exp(-λ|f y j | p ), (31)
where the exponent p is chosen to be less than one, to encourage sparsity on the sharp image gradients. In this work we use p = 0.5. We refer the reader to the original work for full details, but note that the algorithm is easily adapted to non-uniform blur since it involves repeated minimizations of quadratic cost functions of the form

E(f ) = Af -g 2 2 + α D x f -v x 2 2 + α D y f -v y 2 2 , (32)
where v x and v y are intermediate variables of the optimization scheme, used to decouple the exact form of the prior from the main image reconstruction step. For our non-uniform blur model, we use the conjugate gradient algorithm to minimize this cost function.

Another method frequently used for deconvolution is the Richardson-Lucy algorithm [START_REF] Richardson | Bayesian-based iterative method of image restoration[END_REF][START_REF] Lucy | An iterative technique for the rectification of observed distributions[END_REF]. Although originally proposed for convolutional blur, this algorithm can equally be used to invert general linear systems [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]. Using the notation of Equation ( 14) for a known blur, the algorithm iteratively improves the estimate f using the following update equation:

f ← f ⊙ A ⊤ g ⊘ A f , ( 33
)
where g is the observed blurry image, and the matrix A depends on the estimated non-uniform blur. Here, ⊙ represents the element-wise product and ⊘ the element-wise division of two vectors. We have found that for images containing saturated regions (pixels where the signal is clipped and the linear model is no longer valid), such as in Figure 1, the Richardson-Lucy algorithm gives better results, with fewer artifacts around saturated regions such as the bright streetlights.

Single-Image Deblurring Results

We show in this section results of single-image deblurring using the algorithms described in Section 4, with comparisons to results obtained with the original algorithms of [START_REF] Fergus | Removing camera shake from a single photograph[END_REF] and [START_REF] Cho | Fast motion deblurring[END_REF] on both synthetic and real data. Implementation details are discussed in Section 7.

Figure 1 shows a result on a real camera shake blur, using the MAP-ℓ 1 algorithm to estimate the kernel, and the Richardson-Lucy algorithm to perform the final deblurring. The blurry image has many saturated regions (e.g. the bright street lights), and in such cases we found the Richardson-Lucy algorithm to produce significantly better results than True sharp image RMS errors between deblurred results and true sharp image (with ratios to the error obtained with ground-truth kernel in parentheses) We compare the results of blind deblurring for two sizes of blur and three noise levels, and the reconstruction errors are summarized in the table at the top. For each single-axis blur, the table contains the root-mean-square (RMS) errors between the deblurred results and the ground-truth sharp image for blurs with a maximum size of 10 or 20 pixels in the image, using our non-uniform model and the uniform model. In each cell we also show, in parentheses, the ratio between the RMS error and the corresponding error for that blurry image deblurred with the ground-truth kernel. Note that to facilitate comparison without the influence of image priors, the deblurred images were all produced using the Richardson-Lucy algorithm any least-squares based algorithms, such as that of [START_REF] Krishnan | Fast image deconvolution using hyperlaplacian priors[END_REF].

Figures 5 and8 show blind deblurring results on images blurred by real camera shake. Our model, used in both the marginalization and MAP algorithms, is able to capture and remove the blur, while the original algorithms of Fergus et al. and Cho and Lee, using a uniform blur model, fail to find meaningful kernels or good deblurred results. This is explained by both the wide field of view, and the fact that the kernels estimated using our algorithm exhibit significant in-plane rotation.

In Figure 8 (d), we also demonstrate the use of the variational marginalization algorithm of Fergus et al. to produce the deblurred output, as opposed to the Richardson-Lucy algorithm used in Figure 8 (c). Although, for computational simplicity, the kernel estimation step uses a grayscale image, at the convergence of this process the distributions q(w) and q(β σ ) for the kernel and noise variance can be fixed. The variational algorithm can then be run again to estimate q(f c ) for each color channel c separately. In the final step, each color channel can be reconstructed from q(f c ) using Poisson reconstruction [START_REF] Pérez | Poisson image editing[END_REF], before matching the color histogram to that of the blurry image. As can be seen, a good deblurred image is produced, underlining the fact that our blur model is valid throughout the image, and that the kernel produced provides a good description of the true non-uniform blur in the image.

Figure 6 shows a third result of single-image deblurring, using the MAP algorithm. While the uniform blur kernel provides a reasonable estimate of the true blur, and allows us to resolve some of the text on the book's cover, the use of our non-uniform blur model provides a clear improvement, and permits almost all of the text to be read.

Figure 7 shows results for blind deblurring of synthetically blurred images using the two methods, and demonstrates two important points: first, small out-of-plane (e.g. Y -axis) components of a blur are sufficiently uniform that the two models both perform well, although the rotational The marginalization approach, when using our model (d) recovers a useful kernel and a good deblurred image, but when using the uniform model (c) does not. Using the MAP approach, the uniform model (e) finds a reasonable approximation to the non-uniform blur, which is valid on the left side of the image. However, on the right side, the error in the kernel leaves diagonal streaks on the deblurred output. Using our non-uniform model (f), however, avoids this problem. The blur kernels for our model in (d) and (f) cover ±1.3 • along each dimension. We encourage the reader to examine the figures in more detail in the digital version of the paper model performs better. Second, our approach is the only one capable of removing in-plane (Z-axis) blurs, which cannot be represented as convolutions. In this case, and also for the largest out-of-plane blurs, we are able to recover a good sharp image, whereas the uniform approach breaks down due to the blur's non-uniformity. The MAP and marginalization algorithms exhibit similar performance across the different blur sizes and noise levels, although as demonstrated by the displayed kernels, the MAP-ℓ 1 approach tends to find sparser, less contiguous kernels than the marginalization approach.

Figure 9 shows the failure of the MAP algorithm to produce a good result (using the blurry image from Figure 8 (a)) when using the original ℓ 2 regularization proposed by [START_REF] Cho | Fast motion deblurring[END_REF] with our non-uniform blur model. As discussed in Section 4.2.1, the kernel produced is highly non- Fig. 6 Blind deblurring of real camera shake, example 2. The result of blind deblurring on a real camera shake image (a), captured with a shutter speed of 1 second, using the MAP approach of Cho and Lee with both the uniform and non-uniform blur models. Also shown in (b) are some of the local PSFs generated from the blur kernel in (d) at various points in the image. In the blurry image, most of the text on the book cover is too blurred to read. Deblurring the image with the uniform blur model (c) allows some of the text on the cover of the book to be read, however, after deblurring with our non-uniform model (d), all but the smallest text becomes legible. The blur kernel in (d) covers ±0.4 • in θ X and θ Y , and ±0.9 • in θ Z sparse despite the thresholding step, and the deconvolved output correspondingly exhibits many artifacts compared to the MAP-ℓ 1 result in Figure 8 (f).

In Figure 10, we compare our approach to that of Fergus et al. ( 2006) on a real, uniformly blurred image, taken from the dataset of [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF], where the true blur is known, and also known to be uniform. This demonstrates the fact that our model includes uniform blur as a special case; by setting the focal length to be large and applying the constraint that θ Z = 0, we obtain results indistinguishable from those of [START_REF] Fergus | Removing camera shake from a single photograph[END_REF]. When we do not apply 2009), deblurred using kernels estimated with the marginalization algorithm. We show deblurred results and kernels for four cases; (c) uniform blur using original algorithm of Fergus et al., (d) our model with a large focal length F and no in-plane rotation (θ Z = 0), (e) our approach with a large focal length F but with θ Z unconstrained, and (f) the ground-truth (uniform) kernel, provided with the dataset. Note that (d) is indistinguishable from (c), apart from a translation, and that the kernel in (d), while not perfect, does have the same diagonal shape as the true blur, with the non-zeros centered around a single value of θ Z the constraint on θ Z , our algorithm still produces a good result, but unsurprisingly does not perform as well, since the number of kernel elements to be estimated is much larger (K is increased by a factor of 8).

Figure 11 shows the result of the non-uniform MAP-ℓ 1 approach on an image of [START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF]. Although the scene is close to the camera, we are able to obtain a comparable result to that of Joshi et al. without considering the camera's translational motion. This suggests that our decision to ignore the camera's translation is reasonable in practice.

Besides comparing the results of a given algorithm with either a uniform or non-uniform blur model, we can also compare the marginalization and MAP approaches for a given model. In our experiments, we have observed that the MAP algorithm is generally more robust to the level of contrast in the input image. The parameters of the image prior provided by [START_REF] Fergus | Removing camera shake from a single photograph[END_REF] are learnt from a single image of a street scene, so the application of this prior to an image with a very different distribution of pixels is liable to produce poor results. The MAP algorithm however only relies on the ability to predict step edges from a blurry image, and adapts its threshold for predicting these edges depending on the contrast of the image. On an image containing only low-contrast edges then, such as in Figure 1, the marginalization approach (using the street scene prior) fails to find a useful kernel, while the MAP approach finds a good kernel, as demonstrated by the comparison in Figure 1 (c). On the other hand, as discussed by [START_REF] Cho | Fast motion deblurring[END_REF], the performance of the MAP approach is sensitive to the values of the parameters α and β, which must be manually specified, while the marginalization approach has almost no parameters to tune.

Convergence. Neither of the algorithms can guarantee the ability to arrive at a globally optimal solution. However, in practice we have found them both to perform reliably. By finding a sequence of solutions at increasingly fine resolutions, the large scale structures in the blur kernel and sharp image are resolved before the fine details. In the case of the MAP algorithm, each of the individual minimizations over the sharp image f and the blur kernel w is convex, ensuring convergence to a local minimum, even though the overall problem is not jointly convex in both f and w. The gradient prediction step helps direct the optimizations process towards a desirable minimum by encouraging the sharp image to contain step edges. In the marginalization algorithm, we have found that the algorithm converges equally reliably for both the uniform model and our model in a similar number of iterations.

Running time. An important difference between the two approaches is that the MAP algorithm typically takes a much shorter amount of time to run, since the parameter updates for the marginalization algorithm, given in Equations (23)(24)(25)(26), are computationally expensive. Due to this expense, and the larger number of iterations required for our model compared to the uniform model, the marginalization algorithm with our non-uniform model can take several hours to deblur an image of several hundred pixels across on a modern workstation. Deblurring larger images with this method is not currently practical, whereas the MAP algorithm can deblur the same images in under an hour.

Limitations. Both of these algorithms are capable of removing large blurs -up to around 50 pixels across, for both uniform and non-uniform blur. For our model, this corresponds to around 3 • -5 • of rotation around each axis for a photograph whose width and focal length are both 1000px. Since we have assumed that camera translation has a negligible blurring effect, our model (and in general the uniform model too) is unlikely to produce good results on images for which this is not true, due to the depth-dependent blur produced. Another typical failure case for the MAP algorithm comes from the fact that it relies on the ability to predict sharp step edges from blurry ones, which may not be the case on images which contain only fine-scale texture, or where the blur is too large to allow this.

Deblurring with Noisy / Blurry Image Pairs

In this section, we apply our model to the case where, in addition to g, we have a sharp but noisy image f N of the same scene, as proposed by Yuan et al. (2007a). The motivation for this is that in low light, blurry images occur at long shutter speeds, however it is often also possible to use a short exposure at a high ISO setting to obtain a sharp but noisy image of the same scene. While the noisy image may be degraded too badly to allow the direct recovery of a good sharp image, it can initially be used as a proxy for the sharp image, allowing us to estimate the blur kernel ŵ by solving Equation (15). Following this, the kernel is assumed to be known, and used to deblur g, solving (14). The noisy image can also be used to improve this deconvolution step, and Yuan et al. propose a modified version of the Richardson-Lucy algorithm, which uses f N to suppress artifacts in the result.

Kernel Estimation

As discussed in Section 3.2, some prior knowledge must be applied to recover a good kernel estimate. In their algorithm, Yuan et al. (2007a) constrain the kernel to have non-negative elements and unit ℓ 1 norm, however they simultaneously [START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF]. A hand-held image with camera shake (a), from [START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF], with the deblurred results from the original work (b) and using the MAP-ℓ 1 method with our blur model (c). We obtain a comparable result, without the use of additional hardware and without considering the camera's translation during the exposure, despite the scene being close to the camera.

penalize the ℓ 2 norm of the kernel, reducing the sparsityinducing effect of the constraint. To help find a sparse kernel, they propose a thresholding scheme which sets some kernel elements to zero at each iteration. In our approach, we opt to use the ℓ 1 and positivity constraints alone, since they lead naturally to a sparse kernel [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], a fact also exploited by [START_REF] Shan | Rotational motion deblurring of a rigid object from a single image[END_REF] for blur kernel estimation.

In order to estimate the blur kernel, we minimize the following energy function:

E N (w) = ĝ(f N , w) -g 2 2 s.t. ∀ k = 1, . . . , K, w k ≥ 0 and k w k = 1, (34)
where, by analogy with Equation (15), ĝ(f N , w) = B N w, and B N is the matrix whose columns each contain a projectively transformed copy of f N . Similar to Equation (29), this least-squares formulation with non-negative ℓ 1 constraints can be solved efficiently [START_REF] Kim | An interiorpoint method for large-scale ℓ 1 -regularized least squares[END_REF][START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF]. Since the energy function is convex with convex constraints, we can be sure of reaching the global minimum.

For comparison, we have also implemented this algorithm for uniform blurs, using a matrix B N in Equation (34) whose columns contained translated versions of f N , rather than projectively transformed versions.

Image Reconstruction

Having estimated the blur kernel, Yuan et al. (2007a) propose several modifications to the Richardson-Lucy (RL) algorithm, which take advantage of the fact that it is possible to recover much of the low-frequency content of f from a denoised version of f N . Images deblurred with the standard RL algorithm often exhibit "ringing" artifacts -low-frequency ripples spreading across the image, such as in Figure 9 -but using the denoised image it is possible to disambiguate the true low frequencies from these artifacts, and largely remove them from the result. Doing this significantly improves the deblurred results compared to the standard RL algorithm. We refer the reader to (Yuan et al. 2007a) for full details of the augmented RL algorithm, omitted here for brevity. We have adapted the algorithm for our non-uniform blur model, along the same lines as for the standard RL algorithm in Section 4.3.

Results

In this section, we present results with noisy / blurry image pairs, and refer the reader to Section 7 for implementation details. Figures 12 and13 show a comparison between the uniform model and ours, using the algorithm described above to estimate the blur kernels. Having estimated the kernel, we deblur the blurred images using the augmented RL algorithm of Yuan et al. (2007a). As can be seen from the deblurred images obtained with the two models, our results exhibit more detail and fewer artifacts than those using the uniform blur model.

Implementation

The implementation of the variational kernel estimation method presented in Section 4.1 is based on the code made available by [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF] and by [START_REF] Fergus | Removing camera shake from a single photograph[END_REF] 2 . We have modified the algorithm to use our blur model and replaced the parameter update equations with the corresponding versions derived for our bilinear blur model in . A package containing our code is available online3 . The implementation of the image reconstruction algorithm of [START_REF] Krishnan | Fast image deconvolution using hyperlaplacian priors[END_REF] is also based on MATLAB code made available online by the authors 4 The implementations of the Richardson-Lucy algorithm, the algorithm of [START_REF] Cho | Fast motion deblurring[END_REF], and the augmented RL algorithm of Yuan et al. (2007a) are our own, and we use these implementations for both uniform and nonuniform blur models when comparing results. A binary executable for Cho and Lee's algorithm is available, however we did not observe an improvement in the results obtained, and thus use our own implementation to permit a fairer comparison between the results from the uniform and non-uniform blur models.

Sampling the Set of Rotations

One important detail to consider is how finely to discretize the orientation parameter θ. Undersampling the set of orientations will affect our ability to accurately reconstruct the blurred image, but sampling it too finely will lead to unnecessary calculations. Since the kernel is defined over the 3 parameters θ X , θ Y and θ Z , doubling the sampling resolution increases the number of kernel elements by a factor of 8. In practice, we have found that a good choice of grid spacing is that which corresponds to a maximum displacement of 1 pixel in the image. Since we are fundamentally limited by the resolution of our images, reducing the spacing further leads to redundant orientations, which are indistinguishable from their neighbors. Setting the grid spacing in terms of pixels also has the advantage that our 3D blur kernels are de-fined on a grid which allows direct comparison to the pixel grid of the image. We set the size of our kernel along each dimension in terms of the size of the blur we need to model, typically a few degrees along each dimension of θ, e.g.

[-5 • , 5 • ].

Multiscale Implementation

All of the kernel estimation algorithms presented here are applied within a multiscale framework, starting with a coarse representation of image and kernel, and repeatedly refining the estimated kernel at higher resolutions. In the case of single-image deblurring, this is essential to avoid poor local minima, however it is also important for computational reasons in both the single-image and noisy / blurry image pair cases. The kernel at the original image resolution may have thousands or tens of thousands of elements, however very few of these should have non-zero values. For example to solve Equation (34) directly at full resolution would involve transforming f N for every possible rotation under consideration and storing all the copies simultaneously in B N . This represents a significant amount of redundant computation, since most of these copies will correspond to zeros in the kernel, and furthermore B N may have too many columns to fit in the computer's memory. The effect on the computation and memory requirements for single-image deblurring is comparable. Thus, in all of the applications presented in this paper, we use the solution ŵs at each scale s to constrain the solution at the next scale ŵs+1 , by defining an "active region" where ŵs is non-zero, and constraining the non-zeros at the next scale to lie within this region. In the example above, this corresponds to discarding many columns of B N , reducing both the computation and memory demands of the algorithm. We first build Gaussian pyramids for the blurred image (and noisy image, if applicable), and at the coarsest scale s = 0, define the active region to cover the full kernel. At each scale s, we find the optimal kernel ŵs for that scale. We then upsample ŵs to the next scale (s + 1) using bilinear interpolation, find the non-zero elements of this upsampled kernel, and dilate this region using a 3 × 3 × 3 cube. When finding the optimal kernel ŵs+1 , we fix all elements outside the active region to zero. We repeat this process at each scale, until we have found the optimal kernel at the finest scale.

Geometric and Photometric Registration

For the case of noisy / blurry image pairs, the two images are simply taken one after the other with a hand-held camera, so they may not be registered with each other. Thus, we estimate an approximate registration θ 0 between them at the coarsest scale, using an exhaustive search over a large set of rotations, for example ±10 • about all 3 axes using the same step size as for the blur kernel, and we remove this mis-registration from the noisy image. When applying the uniform blur model in this case, we manually estimate the in-plane rotation to best register the two images, as in (Yuan et al. 2007a).

To compensate for the difference in exposure between the noisy and blurry images, at each scale s, after computing ŵs for that scale, we estimate a linear rescaling a by computing the linear least-squares fit between the pixels of g s and those of ĝs ( ŵs , f N,s ), and apply this to the noisy image, i.e. f N ← af N .

Conclusion

We have proposed a new model for camera shake, derived from the geometric properties of cameras, and applied it to two deblurring problems within the frameworks of existing camera shake removal algorithms. We have validated the model with experiments on real and synthetic data, demonstrating superior results compared to the uniform blur model. The model assumes that the motion of the camera during exposure is limited to rotations about its optical center, and is temporally-agnostic to the distribution over camera orientations. Note, however, that camera rotations that are off the optical center can be modeled by camera rotations about the optical center together with translation; these translations should generally be small for rotation centers that are not far from the optical center. The model is not applicable for non-static scenes, or nearby scenes with large camera translations where parallax effects may become significant.

In the future, we plan to investigate the use of our general bilinear model to other non-uniform blurs. We also plan to investigate means of reducing the computational overhead of the model, for example with the use of a suitable approximation strategy, such as [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF]).

Acknowledgements

We are grateful for discussions with Bryan Russell, and comments from Fredo Durand and the reviewers. Thank you to James Miskin and Rob Fergus for making their code available. Financial support was provided by ONR MURI N00014-07-1-0182, ERC grant VisRec no. 228180, the MSR-INRIA laboratory, ANR grant HFIBMR (ANR-07-BLAN-0331-01), the EIT-ICT labs (activity 10863), and the ERC grant Vide-oWorld.

A Parameter update equations for marginalization algorithm

Equations (23-26) cannot be evaluated directly, since they involve expectations over combinations of variables. For implementation they must be expanded and written in terms of the mean and variance of individual variables. Here we give these expansions, which are exactly the parameter updates computed in our implementation. First, let

v w k = w 2 k -w k 2 (35) v f j = f 2 j -f j 2 (36) A ij = k C ijk w k (37) B ik = j C ijk f j (38) ĝi = k j C ijk f j w k . (39)
Then, w

k = βσ i,j C 2 ijk v f j + βσ i B ik 2 (2)
w

(1) k w

(2)

k = βσ i B ik g i -ĝi -βσ i,j C ijk A ij v f j + w k w (2) k (41) f (2) j = βσ i,k C 2 ijk v w k + βσ i A ij 2 (42) f (1) j f
(2) j

= βσ i A ij g i -ĝi -βσ i,k C ijk B ik v w k + f j f
(2) j .

(43)

Finally, in evaluating the parameters of the distribution q(βσ) for the variance of the noise, (Miskin and MacKay 2000, Eqn. 40), it is necessary to evaluate the following quantity:

(g i -ĝi ) 2 = (g i -ĝi ) 2 + j,k C 2 ijk v f j v w k + j A ij 2 v f j + k B ik 2 v w k .
(44)

B Computation of interpolation coefficients

Here we give details of how the values of C ijk can be calculated if bilinear interpolation is used. Note that these are standard interpolation weights, and are not specific to this deblurring application. We consider a pixel g i in the blurry image which is mapped under a homography H k to a point x ′ in the sharp image, i.e. H k x i ∼ (x ′ , y ′ , 1) ⊤ . The point x ′ will, in general, be a sub-pixel location, and the value of the sharp image f at x ′ is interpolated from the 4 pixels surrounding x ′ , using the following weights:

(x j , y j ) C ijk ⌊x ′ ⌋, ⌊y ′ ⌋ ⌊x ′ ⌋, ⌊y ′ ⌋ + 1 ⌊x ′ ⌋ + 1,⌊y ′ ⌋ ⌊x ′ ⌋ + 1,⌊y ′ ⌋ + 1 ⌊x ′ ⌋ + 1-x ′ ⌊y ′ ⌋ + 1-y ′ ⌊x ′ ⌋ + 1-x ′ y ′ -⌊y ′ ⌋ x ′ -⌊x ′ ⌋ ⌊y ′ ⌋ + 1-y ′ x ′ -⌊x ′ ⌋ y ′ -⌊y ′ ⌋
where ⌊ • ⌋ takes the integer part of a positive scalar. The correspondence between a pixel f j 's index j and coordinates (x j , y j ) can be obtained using, for example, the MATLAB functions ind2sub / sub2ind.

  Fig. 1 Modeling non-uniform blur in a shaken image. The blurry image (a) clearly exhibits blur which is non-uniform, as highlighted at different locations in the image. Using the model proposed in this paper, we can describe this blur using a single global descriptor (b), which in this case has been estimated from the blurry image itself, simply by modifying existing algorithms for blind deblurring (see Section 3 for a complete explanation). Close-ups of different parts of the image (c) show the variation in the shape of the blur, which can be accurately reproduced using our model, as shown by the local point spread functions generated from it. As can be seen in (d) and the close-ups in (e), different parts of the image, blurred in different ways, can be deblurred to recover a sharp image

Fig. 2

 2 Fig.2Blur due to translation or rotation of the camera. In this simplified example, we consider capturing a blurry image by either (a) translating the camera through a distance X parallel to the image plane, or (b) rotating the camera through an angle θ about its optical center. We consider the scene point P at a distance D from the camera, whose image is blurred by δ pixels as a result of either of the two motions. In most cases, for a given blur size δ the rotation θ constitutes a significantly smaller motion of the photographer's hands than the translation X (see text for details)

  Y -axis rotation of the camera (c) Z-axis rotation of the camera (d) Arbitrary sequence of rotations

Fig. 7

 7 Fig.7Blind deblurring of synthetic single-axis blurs. A sharp image (top left) with examples of synthetic blur by rotation of the camera about its Y and Z-axis, and the kernels and deblurred results for different cases. We compare the results of blind deblurring for two sizes of blur and three noise levels, and the reconstruction errors are summarized in the table at the top. For each single-axis blur, the table contains the root-mean-square (RMS) errors between the deblurred results and the ground-truth sharp image for blurs with a maximum size of 10 or 20 pixels in the image, using our non-uniform model and the uniform model. In each cell we also show, in parentheses, the ratio between the RMS error and the corresponding error for that blurry image deblurred with the ground-truth kernel. Note that to facilitate comparison without the influence of image priors, the deblurred images were all produced using the Richardson-Lucy algorithm

  Fig.5Blind deblurring of real camera shake, example 1. The result of blind deblurring on a real camera shake image (a), captured with a shutter speed of 1 2 second, using both the marginalization algorithm of Fergus et al. and the MAP approach of Cho and Lee with both the uniform and non-uniform blur models. Also shown in (b) are some of the local PSFs generated from the blur kernel in (d) at various points in the image. The marginalization approach, when using our model (d) recovers a useful kernel and a good deblurred image, but when using the uniform model (c) does not. Using the MAP approach, the uniform model (e) finds a reasonable approximation to the non-uniform blur, which is valid on the left side of the image. However, on the right side, the error in the kernel leaves diagonal streaks on the deblurred output. Using our non-uniform model (f), however, avoids this problem. The blur kernels for our model in (d) and (f) cover ±1.3 • along each dimension. We encourage the reader to examine the figures in more detail in the digital version of the paper

  MAP-ℓ2, uniform (detail) (b) Some local PSFs for (d), magnified (d) MAP-ℓ1, non-uniform (d) MAP-ℓ1, non-uniform (detail)

Fig. 8

 8 Fig.8Blind deblurring of real camera shake, example 3. A hand-held image with camera shake (a), captured with a shutter speed of 1 second, with the results of blind deblurring using the marginalization algorithm of Fergus et al. under both a uniform (b) and non-uniform (c-d) blur model, and the MAP algorithm of Cho and Lee with a uniform (e) and non-uniform (f) blur model. The variational output (d) is estimated using the marginalization algorithm for the non-uniform case (calculated as f q(f ) then converted from gradients to intensities using Poisson reconstruction[START_REF] Pérez | Poisson image editing[END_REF]). The results using our blur model show more detail and fewer artifacts than those using the uniform blur model, as can be seen in the zoomed-in portions shown in the last row. The rotational blur kernels in (c) (f) cover ±0.7 • in θ X and θ Y and ±1.4 • in θ Z

  Fig.11Blind deblurring of an image from[START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF]. A hand-held image with camera shake (a), from[START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF], with the deblurred results from the original work (b) and using the MAP-ℓ 1 method with our blur model (c). We obtain a comparable result, without the use of additional hardware and without considering the camera's translation during the exposure, despite the scene being close to the camera.

  Fig.12Deblurring real camera shake blur using a noisy / blurry image pair. A noisy image (a) and a blurry image (b) captured with a hand-held camera, with the estimated kernels (c-d) and deblurred results (f-g) for the uniform and non-uniform blur models. Also shown for illustration are a selection of the local PSFs generated by the rotational kernel (e). As can be seen in the close-ups (h-k), our result (k) contains more details and fewer artifacts than when using the uniform blur model (j), and reveals features not visible in either the noisy or the blurry image. The non-uniform kernel in (d) covers ±3 • along each dimension. We encourage the reader to examine the images in the digital version of the paper

http://epaperpress.com/ptlens/

http://cs.nyu.edu/~fergus/research/deblur.html

http://www.di.ens.fr/willow/research/deblurring/

http://cs.nyu.edu/~dilip/research/fast-deconvolution/