Josef Sivic

Andrew Zisserman

Jean Ponce

I would like to thank my supervisors

Chapter 1 Introduction

With the explosion of digital photography in recent years, many of us take large numbers of digital photos with cameras or camera-phones. When we review our photos later however, there is sometimes a divergence between what we remember seeing at the time, and what our cameras recorded. This disparity can perhaps be explained by the old adage that we see with our brains, not with our eyes. Whether our photos contain a luridly-dressed tourist that our brain had filtered out at the time, or our photos appear blurry due to camera shake in low light, it is not uncommon that the photos we find ourselves with do not capture what we wanted to record.

In this thesis, we develop models and methods aimed at "restoring" photographs to bring them closer to the images we hoped to record at the time of their taking.

For many people, existing photos capture important and fleeting moments that it may be impossible to recapture. As opposed to proposing new hardware or methods for taking future photographs, we are predominantly concerned with handling images that have already been captured. By using accurate models of the image formation process, and incorporating strong prior information on the images we would like to recover, we endeavour to make software post-processing the "brain" to the digital "eyes" of our cameras. Specifically, our objective in this thesis is to automatically restore photographs, when these photographs contains unwanted occluders as shown in Figure 1.2, or when Removing camera shake blur from photos. The blur caused by camera shake, such as in the blurry image shown here (a), is typically spatially-variant. Most previous work on removing camera shake has assumed a spatially-invariant blur model, leading to deblurred images such as that shown in (b). Using the model for spatially-variant camera shake blur proposed in Chapter 3 and the blind deblurring algorithms described in Chapter 4, we are able to model the spatially-variant blur correctly and obtain superior deblurred results (c). Despite the large blur (around 30 pixels), much of the text that is illegible in the input image, and which is not restored sufficiently with the spatially-invariant blur model, can be read clearly in our deblurred image.

they are blurry due to camera shake, as shown in Figure 1.1. Besides the emotional motivation for improving peoples' photos, more concrete motivations for these tasks might be the need to recover visual information (e.g. vehicle license plate numbers) from blurry photos, to reduce the cost of obtaining high-quality sharp images (by using cheap computational power instead of expensive camera hardware), or to remove trademarked or otherwise sensitive objects from photos before publishing them.

Problem Statement

The types of image degradation considered in this work lead to inherently ill-posed image restoration problems. Starting from only a degraded image, we wish to recover a good image of the same scene. If there is a loss of information, or we have more unknowns than observations, there may be a large family of valid solutions, which we must somehow choose amongst when producing the "restored" image. This is the case in both the problems discussed here. For deblurring images, where the blur is unknown, we must estimate both the parameters describing the blur and the sharp image. Since the sharp image has the same number of pixels as the blurry image, we evidently have more unknowns than observations. When we wish to remove occluders from photographs, the recorded image contains no information about what is behind the occluder; it could be concealing a building, a tree, a patch of grass, or a group of people. Making the right choice and producing visually-pleasing results requires good models of the image formation process, and equally importantly, good prior information about the unknowns.

Contributions

The main contributions of this thesis can be divided broadly into deblurring of camera shake, and removing large occluders from photographs.

Restoring Photographs Blurred Due to Camera Shake

The main contributions with respect to deblurring camera shake are first to demonstrate that camera shake is mainly caused by 3D rotation, as opposed to 3D translation of the camera, causing spatially-variant blur, and second to derive a geometric model for the blur process based on this. We propose a formulation of this model which is directly applicable in existing blind deblurring algorithms, and consequently demonstrate the ability to remove spatially-variant camera shake blur from photographs. Figure 1.1 demonstrates an example result on a real image blurred by camera shake, and shows that the spatially-variant blur model allows us to recover a significantly better result than only modelling spatially-invariant blur.

We propose an efficient approximation for this model which significantly reduces the computational burden associated with using it, and makes spatially-variant blind deblurring of camera shake practical for real images. Finally we address deblurring of images containing clipped, or saturated pixels. We propose a forward model that includes sensor saturation, and propose a non-blind deblurring algorithm that incorporates this model while preventing artefacts from appearing in the deblurred results. shows a result deblurring a real image, with saturated highlights in the background.

The image is taken from Canon's "Second Shot" advertising campaign1 , in which the owners of badly blurred or degraded images were given the chance to go back and retake their photographs with a new camera. Our result is of course not as good as Canon's retaken photo, but is a little more achievable for those without the means or desire to return to the scene of every blurry photo they ever took.

Removing Occluders from Photographs

Our contribution with respect to removing large occluders from photographs is to propose an automatic system for replacing occluded regions of photographs using photographs of the same scene, retrieved from the Internet. Our system is able to take an input image with the occluders marked by a user, and return a restored image, where the occluded region has been seamlessly replaced with realistic image content that corresponds to the true underlying scene. Figure 1.2 demonstrates an example result produced by our system. Note that the underlying scene is realistically rendered, without any knowledge of the 3D structure of the scene or of the environmental conditions.

Thesis Outline

We begin in Chapter 2 with some background on image restoration and a discussion of relevant work from the literature.

In Chapter 3 we examine blur caused by camera shake. We begin by considering the geometric relationship between motion of the camera and the apparent motion of the image. Following this, we propose and discuss a novel formulation for camera shake blur, which handles spatially-variant blur naturally. We present the discrete equivalent of our model, and provide a comparison to other models for camera shake blur. Finally we discuss some practical implementation considerations.

In Chapter 4 we demonstrate blind deblurring of shaken images, demonstrating that our model is easily applicable within existing deblurring paradigms. We first apply our model within two algorithms for single-image deblurring, using two different approaches to estimating the parameters of the blur. We present results of these algorithms and discuss some of the advantages and disadvantages of using our model instead of the classical spatially-invariant blur model. Secondly we demonstrate deblurring of images when there is a additional noisy image available of the scene. We conclude with some implementation issues relevant to all the algorithms presented.

In Chapter 5 we introduce an efficient approximation scheme for our spatially-variant blur model. We begin by examining the bottlenecks in the deblurring process, before presenting our approximation scheme and showing how our approximation greatly reduces the computational burden of these steps. We present results using the approximation, and compare to those produced using the exact model.

In Chapter 6 we approach deblurring of images containing clipped / saturated pixels.

We discuss possible ways of handling such pixels in the deblurring process, and propose a modification to the image formation model which incorporates the saturation process.

Further to this, we examine how saturated pixels cause artefacts in the deblurred results, and propose an algorithm which handles these clipped pixels explicitly to prevent visible artefacts being introduced. Finally we present results of our method and compare to 1.4 Publications other methods.

In Chapter 7, we approach the removal of large occluders from photographs of famous landmarks. We discuss the visual search engine which we use to retrieve a set of example images of the same scene, and how we register these images geometrically and photometrically to an input image. Following this we present our algorithm for replacing the occluders, by combining pixels from several of the example images. We present results of our system on a range of input images.

Finally in Chapter 8 we summarise the contributions of this thesis and discuss potential directions for future work.

Publications

Parts of the work in this thesis have appeared in the following publications:

O. Whyte, J. Sivic, and A. Zisserman. Get out of my picture! Internet-based inpainting.

In Chapter 2

Background and Related Work

Image restoration is a vast field of research, having been studied for many decades since the advent of digital imagery. No such thing as a perfect camera exists, and all images are, to some extent, affected by noise and blur. Classical work on digital image restoration tackled these problems, and are covered in depth by [START_REF] Gonzalez | Digital Image Processing[END_REF]. Recently, the field has expanded to include more diverse sources of degradation such as chromatic aberrations or missing pixels. Without doubt, the use of a priori information has been crucial in many recent advances, often allowing surprising amounts of of information to be recovered (or perhaps, hallucinated) from images degraded by noise, blur, or other sources of corruption. This information may concern the statistics of undegraded images, the properties of the underlying scene, or other factors such as idiosyncrasies of human perception. In this chapter we recap some of the most relevant work to the subjects covered in this thesis: deblurring and inpainting.

The problem of restoring images can be broken down into several components. First, a model is needed to relate the undegraded, "ideal" image to the observed image produced by the camera. Second, the parameters of this model must be estimated, and finally the sharp image can be estimated, given the model and the estimated parameters.

We begin by introducing a general forward model for image degradation in Sec-9 tion 2.1, and discuss the concepts which underpin successful image restoration algorithms. Following from a probabilistic model, which relates the unknown (latent) ideal image with the observations and our prior knowledge, we discuss the two general approaches to image restoration by maximising or marginalising this probability. In Section 2.2 we discuss existing algorithms for deblurring images when the blur is known.

Following this, in Section 2.3 we recap some of the existing work for deblurring images when the blur is unknown. In Section 2.4 we introduce the problem of parameterising spatially-variant blur, and discuss some recent work on this problem. Finally in Section 2.5 we address the relevant work and background to the problem of "inpainting".

Notation

In this thesis we will use some notation consistently, which we introduce here for reference. We use bold lower-case letters to denote vectors, e.g. f , and bold upper-case letters to denote matrices, e.g. A. We use subscripts on non-bold letters to index into vectors and matrices, e.g. f j indicates the j th element of the vector f , while A ij indicates element (i, j) of matrix A.

Calligraphic letters such as U denote sets or domains, depending on the context. We denote functions with lower-case letters, i.e. f : R → R + .

To denote 2D discrete images we use 1D vectors, e.g. f ∈ R N for an image of height H and width W , where N = H × W . In this context we denote 2D convolutions between two images as operating directly on the 1D vectors, i.e. u * v, and likewise we denote the 2D discrete Fourier transform simply as a function mapping a vector to a vector F : R N → C N . We denote the Hadamard (element-wise) product between two vectors by u • v. We abuse standard notation somewhat and write the element-wise division of two vectors simply using a fraction u v . When several objects of the same type are collected into a set, we use numeric superscripts in parentheses to index into this set, e.g. d (q) ∈ {d (q) }. When objects of the same type are identified by a symbolic property, we use superscripts without parentheses, e.g. d x , d y to indicate derivative filters in the x and y directions.

Image Degradation and Restoration

To frame the problems discussed in this thesis, we begin with a general forward model of image degradation. There are four important components in the image degradation process:

• The (observed) degraded image g : Ω → R + , which is the image output by the camera's sensor. The domain Ω ⊂ R 2 is the 2D plane of the camera's sensor.

• The (latent) sharp image f : Ω ′ → R + . This is the underlying, ideal, sharp image of the scene, which we would like to recover. In most cases, the domain Ω ′ ⊂ R 2 is chosen to be the same as Ω.

• The degradation operator H that acts on f , and which describes how light from the sharp image f is distributed in the observed image g. Depending on the situation, H may be known in advance (e.g. from optical properties of the camera) or may be unknown (e.g. arbitrary camera motion during camera shake).

• The random noise N that perturbs the recorded image after the sharp image has been degraded by H.

These components are combined in the generic image degradation model g(x) = N (Hf )(x) ,

(2.1) where x is a point (x, y) ∈ Ω. We will denote the "noiseless" degraded image Hf by g * , and also assume that the degradation operator H is linear, allowing us to write

g * (x) = (Hf )(x) = Ω ′ h(x, x ′ )f (x ′ ) dx ′ , (2.2)
where x ′ is a point (x ′ , y ′ ) ∈ Ω ′ , and h : Ω×Ω ′ → R + is referred to as the impulse response of H, or the point-spread function (PSF). For a point light source with magnitude 1 at a particular point x ′ 0 ∈ Ω ′ , the 2D function h( • , x ′ 0 ) is the response produced by H, and 11 describes how the light from this point is spread over the observed image [START_REF] Gonzalez | Digital Image Processing[END_REF]. In Section 2.1.2 we discuss possible forms for the PSF, and in Section 2.1.4

we discuss several noise models that are useful in practice.

It is important to note that in general, for a given degraded image g, both the sharp image f and the PSF h may be unknown. This makes the image restoration problem particularly difficult. The forward model plays an important role in solving this problem, but some additional information is typically needed to recover a good estimate of the latent image. This additional information can come in the form of statistical priors, which encourage the latent image to look realistic, and improve the conditioning of the problem. In Section 2.1.5 we discuss priors for the latent image and for the PSF.

Spatially-invariant Blur

A common assumption when modelling image blur is that the PSF is spatially-invariant, which is to say that there exists a function a : R 2 → R + such that h(x, x ′ ) = a(xx ′ ).

The function a is typically referred to as a convolution kernel, and in this case the dimensionality of the PSF is reduced from four to two, and Equation (2.2) reduces to a 2D convolution of f with a:

g * (x) = (a * f )(x) = Ω ′ a(x -x ′ )f (x ′ ) dx ′ .
(2.3)

In general, blur may be spatially-variant, and we will return to this issue in Section 2.4, and again in Chapter 3.

The Discrete Setting

Real cameras are equipped with a discrete set of pixels, and output a discrete set of samples of the degraded image, denoted by the vector g ∈ R N + , where N = H × W pixels for an image with H rows and W columns. We consider the sharp image also to be discrete: f ∈ R N + . We use i to index into the degraded image g, i.e. g i = g(x i ), where x i ∈ Ω is the coordinate of the i th pixel. Likewise, we use j to index into the sharp image f , such that f j = f (x ′ j ) for a coordinate x ′ j ∈ Ω ′ . Finally, we note that to evaluate an image at arbitrary (sub-pixel) locations, we interpolate from nearby pixels.

In this work, we use linear interpolation schemes, whereby sub-pixel values of an image, say g(x) are interpolated as a linear combination of nearby pixels: (2.4) where the coefficients b(x, x i ) are calculated using a standard method such as bilinear or bicubic interpolation.

g(x) = i b(x, x i )g i ,
In this discrete setting, we can write Equation (2.2) as (2.5) or in matrix-vector notation,

g * i = j A ij f j ,
g * = Af , (2.6)
where the N × N matrix A captures the discrete PSF. Each column of the matrix A contains the PSF for the corresponding pixel in the latent image f . In most cases of blur, the light from each pixel in f is spread over a relatively small number of nearby pixels in g. As a result, the PSF matrix A for an image is usually sparse (contain a relatively small number of non-zero values). When the PSF is spatially-invariant, we denote the discrete convolution kernel by a, and write g * = a * f .

(2.7)

Types of Degradation

In this thesis we are concerned with two types of PSF: those arising from camera shake, which causes image blur, and those arising when some pixels from the observed image are missing or deleted.

• Image Blur. When an image is degraded by blur, light from a single point in f is spread across a region in g. In this case, the PSF A will have many non-zeros in 13 each column. When the PSF is unknown, the problem is typically referred to as "blind" deblurring. If the PSF is known a priori , or has already been estimated, the problem of recovering the sharp image is, by contrast, referred to as "non-blind" deblurring.

• Missing or Deleted Pixels. When the image has been corrupted in certain regions, or contains occluders that we wish to remove, the PSF is simply the identity, i.e. A = I, since we assume that no blurring has occurred. We model the missing pixels as noise with very high uncertainty, such that their observed intensity is unrelated to their latent intensity.

Probabilistic Formulation

Given the general model of image degradation discussed above, a natural starting point for image restoration algorithms is to write down a probabilistic generative model for the observed image. If we know the type of noise affecting the observed image, we can write down the likelihood of the observed image p(g|f , A), which is the probability density of g, conditioned on f and A. If we then wish to find the latent image and PSF which best match the observed image g, an obvious choice would be to find the f and A which maximise this likelihood. Due to the loss of information that occurs in the degradation process, however, image restoration algorithms which simply maximise this likelihood are known to be ill-conditioned. Algorithms of this sort are prone to producing results containing artefacts, and in which noise is amplified.

If we have some prior knowledge about the latent variables being recovered, then using Bayes' rule, we can formulate the posterior distribution for the unknowns (the sharp image f and also, if unknown, the PSF A). The posterior incorporates the likelihood, which arises from the random noise in the observed image, and also prior knowledge about the unknowns, and in our case is given by

p(f , A|g) ∝ p(g|f , A)p(f )p(A)
for an unknown PSF (2.8)

p(f |g, A) ∝ p(g|f , A)p(f )
for a known PSF.

(2.9)

The priors p(f ) and p(A) can either be manually defined, or learnt from example data.

By incorporating this prior knowledge about the unknown variables, the posterior ameliorates the instability of estimating f or A from the likelihood alone. We will discuss several popular choices for these terms in the following sections.

Given the posterior distribution, it is the task of image restoration algorithms to find the "most probable" sharp image f . The idea of what is "most probable" however is not clear: is it the image f that maximises the posterior probability over all possible f and A, or is it the expected value of f under the distribution? In subsequent sections we will describe different approaches to the image restoration problem, based on different answers to this question. Finally we note that although not all image restoration algorithms have a probabilistic interpretation, this formulation can nevertheless provide some useful intuition into why they succeed.

Maximim A Posteriori

One popular method of finding the "most probable" values of the unknown variables in a system is to find those values which maximise the posterior probability. In this approach, for a known PSF, the estimated latent image f is found as

f = arg max f p(f |g, A), (2.10)
while for unknown PSF, the latent image and PSF are estimated simultaneously

{ f , Â} = arg max f ,A p(f , A|g). (2.11)
This maximisation is typically addressed by first transforming the the probability maximisation problem into an energy minimisation problem. The forms of the likelihood and priors in Equation (2.8) are generally chosen such that the posterior can be written as a Gibbs distribution, with the form

p(f , A|g) = 1 Z exp -1 T U (f , A) , (2.12)
where U is an energy function which depends on g, and Z is a normalising constant. For this distribution, the sharp image f and PSF A which maximise the posterior probability are those which minimise the energy function U . Given that U (f , A) ∝log p(f , A|g) + const., we can write the MAP problem for the posterior in Equation (2.8) as min

f ,A F g, Af + αρ f (f ) + βρ A (A), , (2.13)
where the function F is derived as the negative log-likelihood, and penalises latent images or PSFs which do not agree with the observed data, while the functions ρ f and ρ A are derived as the negative log-priors, and penalise latent images or PSFs which are unlikely under those priors. The function F is referred to as the data fidelity term, while ρ f and ρ A are referred to as the regularisation terms.

The difficulty in the MAP approach is that the search space is very large, with potentially millions of unknowns. Furthermore the energy function may not be convex, in which case it is not generally possible to reach the globally optimal solution.

Marginalisation

An alternative approach to maximising the posterior probability is to attempt to find the expected (mean) values of the unknowns under the posterior distribution. In this approach, for a known PSF the estimated latent image is calculated as

f = f p(f |g, A) df .
(2.14) while for unknown PSF, the latent image and PSF would be estimated as (2.16) This approach sometimes has advantages over maximising the posterior, since the MAP solution may find a peak which has high probability density but has very little probability mass below it. Marginalising will find a solution which is a combination of all posssible solutions, weighted by their probability density, and so a wide peak with a lower maximum density may have more influence than a narrow but higher one.

f = f p(f , A|g) df dA (2.15) Â = A p(f , A|g) df dA.
The difficulty in estimating the unknowns by marginalisation is that the expectations in Equations (2.14) to (2.16) are analytically intractable. To evaluate these integrals some approximation strategy must be used, of which there are in general two types: stochastic approximations and parametric approximations. Stochastic approximation methods, such as a Markov Chain Monte Carlo [START_REF] Neal | Probabilistic inference using markov chain monte carlo methods[END_REF], attempt to evaluate the integrals stochastically by drawing samples from the true posterior distribution. Such algorithms expend most of their effort on drawing samples in a way that ensures convergence to the true distribution. Parametric approximation methods, such as ensemble learning [START_REF] Lappalainen | Ensemble learning[END_REF], attempt to find a parametric approximation of the posterior for which the integrals become tractable. These algorithms spend most of their computational effort estimating the parameters of the approximating distribution, while the final marginalisation under the parametric distribution is often trivial or easy to compute.

Noise Models

The likelihood term p(g|f , A) in the posterior distribution is defined by the type of noise that is present in the observed image. Here we recap three noise models which are used in this work: Gaussian noise, Poisson noise, and uniform noise. In all cases we assume that all pixels g i in the observed image are independent, conditioned on the sharp image f and the PSF A, i.e. p(g|f , A) = i p(g i |f , A). In deblurring, Gaussian and Poisson noise are two widely applicable models [START_REF] Boncelet | Image noise models[END_REF]. Uniform noise on the other hand, is useful for modelling corrupted or deleted pixels, such as in inpainting.

Poisson Noise

Poisson noise is a realistic noise model to use for images, since the arrival of photons at the sensor is naturally modelled as a Poisson process. Indeed, in low light, "quantum noise" (noise caused by the random arrival times of individual photons) dominates over other types of noise in a digital camera. For a Poisson random variable z with mean λ, we write z ∼ Pois(λ), and z has the probability density function (PDF)

p(z) = λ z e -λ z! .
(2.17)

Applying this in our degradation model, we have g i ∼ Pois g * i :

p(g i |g * i ) ∝ g * i g i e -g * i g i ! , (2.18) and -log p(g i |g * i ) = -g i log g * i -g * i -log g i ! + const. (2.19)
Figure 2.1 plots the likelihood p(g i |g * i ) as a function of the unknown g * i , and the corresponding negative log-likelihood, which is used as the data fidelity term in the MAP problem in Equation (2.13). As can be seen, the peaks of p(g i |g * i ) become wider at higher g i , indicating greater uncertainty in g * i . In image restoration tasks, where we wish to estimate g * i , this corresponds to penalising errors less at bright pixels than at dark ones. Also, as g * i → 0,log p(g i |g * i ) → ∞, implicitly enforcing positivity on g * i . Finally, we note that the negative log-likelihood (Equation (2.19) and Figure 2.1 (b)) is convex, allowing a global minimum to be found if we use this noise model in a MAP approach.

Gaussian Noise

Gaussian noise is perhaps the most common noise model used in image processing, due to its tractability and wide applicability. For a Gaussian random variable z with mean µ and variance σ 2 , we write z ∼ N (µ, σ 2 ), and z has the PDF

p(z) = 1 √ 2πσ exp (z -µ) 2 2σ 2 . (2.20)
Applying this in our degradation model, we have

g i ∼ N g * i , σ 2 : p(g i |g * i ) = 1 √ 2πσ exp g i -g * i 2 2σ 2 (2.21) and -log p(g i |g * i ) = 1 2σ 2 (g i -g * i ) 2 + const. (2.22)
Figure 2.2 plots the likelihood p(g i |g * i ) for several values of g i and the corresponding negative log-likelihoods. In contrast to Poisson noise, the variance remains constant at all brightnesses, and the likelihood does not enforce positivity on g * i . Minimisation of the negative log-likelihood is a linear least-squares problem, which is typically easier to solve than that in Equation (2.19), and for which many good algorithms exist. Furthermore it is convex, making it a good candidate for MAP approaches.

Uniform Noise

As well as the random fluctuations in the light arriving at the sensor, we can also model corrupted or deleted pixels in the image as noise with a uniform distribution. For a uniform random variable z, we write z ∼ Unif(a, b), with a < b, and z has the PDF

p(z) =              0 z < a 1 b-a a ≤ z ≤ b 0 b < z (2.23)
Applying this in our degradation model, we have g i ∼ Unif(0, 256), and

-log p(g i |g * i ) =              ∞ g i < 0 const. 0 ≤ g i ≤ 256 ∞ 256 < g i . (2.24)
Under this noise model, an observed pixel g i is independent of the noiseless pixel g * i , and is thus decoupled from the latent image f . This noise model is thus applicable when the degraded image contains pixels that we would like to replace, since we would like the recovered image to be independent of the replaced pixels. Recently, [START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF] have used the uniform noise model in deblurring, to handle outlier pixels in the blurry image, e.g. due to "hot pixels", or salt and pepper noise.

Priors

As mentioned earlier, image restoration problems are frequently ill-posed, such that methods which use the data fidelity term alone may amplify noise in the solutions, or fail completely to find a good solution. As a result, some prior knowledge of the unknowns is often included in the restoration problem to improve the conditioning and stabilise the solution. In this section we discuss several such priors for the latent image and the PSF.

Priors for the Sharp Image

An important aspect of image restoration is the prior on the latent image p(f ). Many different image priors have been proposed for image reconstruction tasks such as deblurring, denoising, and super-resolution, often based on the statistics of natural images. Early work on the statistics of natural images observed that the amplitude spectrum, when radially averaged over orientation, approximately follows a power law [START_REF] Field | Relations between the statistics of natural images and the response properties of cortical cells[END_REF], or that the entropy of sharp images is higher than that of blurred images [START_REF] Gull | Maximum entropy method in image processing[END_REF].

Recent work has focused particularly on the distribution of the responses to dictionaries of filters. For a set of filters {d (q) }, these priors model the statistics of the image responses {f (q) }, where f (q) = d (q) * f . The simplest variant of these considers the response of the image to horizontal and vertical derivative filters {d x , d y }, essentially modelling the strength of the correlation between neighbouring pixel values. This use of spatial correlations has a long history [START_REF] Molina | Using spatial models as priors in astronomical image analysis[END_REF], however early work used Gaussian distributions to model these correlations, which leads to overly-smoothed results when applied to image restoration. It has been observed that the derivatives of natural images in fact generally follow heavy-tailed sparse distributions, rather than a Gaussian distribution. These distributions have a large probability mass near zero, reflecting the fact that most pixels are strongly correlated with their neighbours. However, the heavy tail is also necessary to model sharp edges, which have large gradient magnitudes, and which are not adequately captured by a Gaussian distribution. Beginning with the use of total variation and related norms [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], priors of this nature have been used extensively in image restoration tasks [START_REF] Fergus | Removing camera shake from a single photograph[END_REF][START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF][START_REF] Krishnan | Blind deconvolution using a normalized sparsity measure[END_REF][START_REF] Levin | Image and depth from a conventional camera with a coded aperture[END_REF][START_REF] Shan | High-quality motion deblurring from a single image[END_REF][START_REF] Tappen | Exploiting the sparse derivative prior for super-resolution and image demosaicing[END_REF]. Priors of this type have proven to be effective at suppressing noise in restored images while encouraging sharp edges to appear. The disadvantage of these priors is that the corresponding regularisers ρ f in the MAP problem (Equation (2.13)) are generally non-convex. This implies that we cannot find the globally optimal latent image. However, in practice this does not generally cause visible degradation in the restored image. [START_REF] Fergus | Removing camera shake from a single photograph[END_REF] proposed a prior on the sharp image which assumes that the derivatives for all pixels are independent and follow a mixture of C zero-mean Gaussians with mixture weights π c and variances v c (also known as a Gaussian scale mixture):

p(f (q) ) = j C c=1 π c exp   - f (q) j 2 2v c   .
(2.25) Other authors have used priors with the form While priors on the derivatives of images concern pairs of neighbouring pixels, other authors have focussed on the statistics of small image patches. Some authors have proposed to use simple sparse distributions with large dictionaries of filters to model this distribution, such as wavelets [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] or learnt dictionaries [START_REF] Aharon | K-SVD: an algorithm for designing overcomplete dictionaries for sparse representation[END_REF]Mairal et al., 2008;[START_REF] Olshausen | Emergence of simple-cell receptive field properties by learning a sparse code for natural images[END_REF]. Others have proposed to use more 23 complex distributions and learn their parameters [START_REF] Roth | Fields of experts: a framework for learning image priors[END_REF][START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF].

p(f (q) ) ∝ j exp -φ(f (q) j ) (2.

Priors for the PSF

When we wish to estimate the PSF for a blurry image, we also need priors for the parameters of the PSF. For example, when the blur is spatially-invariant and parameterised by a convolution kernel a, we use the prior p(a). Since in this work we are primarily concerned with blur caused by camera shake, we will discuss priors related to this type of PSF here.

Successful algorithms for estimating camera shake PSFs generally share the same two pieces of prior information about the blur kernel being estimated. First, all its elements are non-negative, since the image formation process is additive, with sensor elements accumulating photons, i.e. a k ≥ 0 for k = 1, . . . , K. 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. Figure 2.4 shows some typical camera shake kernels, which clearly exhibit these properties.

The nonnegative sparsity prior has been a prominent feature of previous camera shake removal algorithms, and has also been leveraged for the alignment of blurred image with non-blurred images (Yuan et al., 2007a). [START_REF] Fergus | Removing camera shake from a single photograph[END_REF] (2.28)

Figure 2.4 (b) plots these priors along with a histogram of values from some real camera shake blur kernels.

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", i.e. a = Φb where Φ is the matrix of curvelets, and place an exponential prior on the coefficients b of the curvelets, rather than on the kernel elements directly. By forcing the kernel to be constructed from a set short curves, this approach encourages the kernel to contain thin linear structures, which intuitively correspond to the camera's motion path. [START_REF] Cho | Fast motion deblurring[END_REF] and Yuan et al. (2007b) use simple Gaussian priors to regularise the kernel values, followed by thresholding such that most of the kernel elements are set to zero.

Algorithms for Non-Blind Image Deblurring

When the PSF for a blurry image is known, the problem of estimating the sharp image f is typically referred to as "non-blind" image restoration. This problem has also been studied for a long time, e.g. [START_REF] Wiener | Extrapolation, Interpolation, and Smoothing of Stationary Time Series[END_REF]. Many algorithms exist for non-blind deblurring, perhaps most famously the Richardson-Lucy algorithm [START_REF] Lucy | An iterative technique for the rectification of observed distributions[END_REF][START_REF] Richardson | Bayesian-based iterative method of image restoration[END_REF], which recovers the maximum-likelihood estimate of the sharp image under a Poisson noise model. For a Gaussian noise model, the maximum-likelihood solution can be found using standard least-squares methods such as conjugate gradient descent (Shewchuk, 1994), and for this reason, Gaussian noise has been favoured in the recent literature. Much work has revolved around the incorporation of strong priors, such as those described in Section 2.1.5, to suppress noise in the output while encouraging realistic image features to appear. In theory, the image priors are independent of the noise model used, so we can mix and match. In practice, particular authors have considered certain combinations. In general, maximum likelihood and maximum a posteriori algorithms are used to estimate the latent image, as opposed to marginalisation algorithms, and in this section we do not consider marginalisation approaches to non-blind deblurring, although this is feasible, as mentioned later in Section 4.3.3.

Non-Blind Deblurring With Poisson Noise

One very popular algorithm for non-blind deblurring under a Poisson noise model is the Richardson-Lucy algorithm [START_REF] Lucy | An iterative technique for the rectification of observed distributions[END_REF][START_REF] Richardson | Bayesian-based iterative method of image restoration[END_REF]. This algorithm has been shown to converge to the maximum-likelihood solution under this noise model [START_REF] Shepp | Maximum likelihood reconstruction for emission tomography[END_REF]. Variants of this algorithm have also appeared for general inversion and matrix-factorisation problems [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]. We can write the likelihood of the blurry image g as a function of f as follows:

p(g|f ) = i g * i g i e -g * i g i !
, where (2.29)

g * i = j A ij f j .
(2.30)

The negative log-likelihood is then:

L pois = i g * i -g i log g * i + log(g i !).
(2.31)

As a maximum-likelihood algorithm, the Richardson-Lucy algorithm minimises L pois over f . The Karush-Kuhn-Tucker conditions for the minimum are, for all f j :

f j ∂ ∂f j L pois = 0, (2.32) and ∂ ∂f j L pois ≥ 0 if f j = 0.
(2.33)

The Richardson-Lucy update rule can be derived as a fixed-point iteration from the condition in Equation (2.32):

f j i ∂ ∂f j (g * i -g i log g * i + log(g i !)) = 0 (2.34) f j i ∂g * i ∂f j - g i g * i ∂g * i ∂f j = 0 (2.35) f j i A ij = f j i A ij g i g * i ,
(2.36)

Assuming that the blur is conservative, i.e. i A ij = 1, we obtain the update rule by replacing the unknowns f j and g * i with our estimates fj and (A f ) i = j ′ A ij ′ fj ′ , and adding iteration indices:

f t+1 j = f t j i A ij g i (A f t ) i , (2.37)
where t indexes the iteration. In matrix-vector form, this update equation is

f t+1 = f t • A ⊤ g A f t , (2.38)
where • represents the Hadamard (element-wise) product and the fraction represents element-wise division of two vectors. Note that the convergence of the Richardson-Lucy algorithm can be improved significantly with acceleration schemes, such as that of [START_REF] Biggs | Acceleration of iterative image restoration algorithms[END_REF], however for simplicity we do not discuss these here.

Incorporating Image Priors

In the form in Equation (2.38), the Richardson-Lucy algorithm corresponds to a maximumlikelihood algorithm, and as such noise is amplified in the deblurred result. The most popular way of mitigating this problem is to stop the algorithm before it has converged, typically after 20-50 iterations. This "early-stopping" technique is a form of regularisation, although it does not correspond to a probabilistic prior that can be written down [START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF]. Probabilistic priors however can also be incorporated in the algorithm by starting from the posterior p(f |g, A) = p(g|f , A)p(f ), which leads to the following regularised update rule:

f t+1 = f t 1 -α∇ρ f ( f t ) • A ⊤ g A f t , (2.39)
where α is the regularisation weight, and ∇ρ f ( f t ) is the vector of partial derivatives of the regulariser ρ f with respect to each latent pixel, evaluated at the current estimate f t , i.e. ∇ρ f ( f t ) j = ∂ρ f ∂f j f = f t . Using this formulation, [START_REF] Tai | Richardson-Lucy deblurring for scenes under a projective motion path[END_REF] have investigated the use of various types of regularisation with the Richardson-Lucy algorithm in detail.

Figure 2.5 shows a synthetic example of an image blurred and corrupted by Poisson noise, along with the non-blind deblurring results using both the standard and regularised Richardson-Lucy algorithm. Note that this method limits the magnitude of the regularisation weight α, as the update vector may become negative if α is too large. [START_REF] Welk | Robust variational approaches to positivity-constrained image deconvolution[END_REF] avoids this limitation by splitting the regularisation gradient into positive and negative contributions: where [z] ± = 1 2 (z ± |z|). This update rule ensures that the non-negativity constraint on f is not violated, regardless of the value of α.

f t+1 = f t 1 -α[∇ρ f ( f t )] - • A ⊤ g A f t + α[∇ρ f ( f t )] + , ( 2
The Richarson-Lucy is an iterative algorithm and may require a large number of iterations to converge. As an alternative method for quickly finding the MAP estimate of the latent image with total-variation (TV) regularisation, Getreuer (2010) propose a fast algorithm using a "split Bregman" method. [START_REF] Yuan | Progressive inter-scale and intra-scale nonblind image deconvolution[END_REF] propose a method for suppressing ringing artefacts in non-blind deblurring.

Non-Blind Deblurring With Gaussian Noise

Many approaches to non-bind deblurring are based on an additive uniform Gaussian noise model, rather than a Poisson noise model, since it is very amenable to optimization, and is a realistic model in practice. For Gaussian noise, the likelihood is

p(g|f ) = 1 √ 2πσ i exp (g i -g * i ) 2 2σ 2 .
(2.41)

The negative log-likelihood is then

L gauss = log √ 2πσ + i (g i -g * i ) 2 2σ 2 .
(2.42)

Since ĝ is linear in the sharp image f , minimising L gauss is a linear least-squares minimisation problem, which is convex and can be solved with almost any standard least-squares algorithm, e.g. conjugate gradient descent. The solution to this problem, which is the maximum-likelihood sharp image, can be written using the pseudo-inverse of A as

f = (A ⊤ A) -1 A ⊤ g (2.43)
One particular case of interest under a Gaussian noise model is that of spatiallyinvariant blur, i.e. g * = a * f . In this case, using Parseval's theorem, L gauss can be transformed into N independent 1D quadratic minimisations in the frequency domain, allowing the maximum-likelihood solution to be obtained directly by pixel-wise division in the frequency domain [START_REF] Gamelin | The dimensionality of scene appearance[END_REF]:

f = F -1 F (a) * • F (g) F (a) * • F (a) .
(2.44)

Incorporating Image Priors

The maximum-likelihood solution for Gaussian noise in non-blind deblurring is typically very ill-conditioned. This is most clearly seen in Equation (2.44), where the denominator is the power spectrum of the blur kernel a. If this power spectrum has some small or zero values, those frequencies will be unstable and may be amplified in the solution. If the observed image were noise-free and we had infinite computational precision available, this would not be a problem, however in practice this solution is ill-conditioned. To avoid this, most non-blind deblurring uses some form of regularisation to smooth the solution. The resulting regularised least-squares problem is

min f Af -g 2 2 + αρ f (f ), (2.45)
where α is a weight that balances the data-fidelity term and the regularisation term.

The simplest form of regularisation to add to a Gaussian noise model is ℓ 2 (Tykhonov) regularisation directly on the elements of f , such that ρ f (f ) = f 2 2 . This produces a cost function which remains quadratic, and is thus straightforward to solve. Unfortunately, in the case of image restoration, such a regulariser does not correspond to a realistic prior for natural images.

Many authors have proposed algorithms to exploit the priors discussed in Section 2.1.5 to regularise least-squares solutions, which amounts to adding a smoothness penalty to the negative log likelihood. One of the most popular regularisers is derived from the sparse gradient prior of Equation (2.26), which penalises the derivatives of the deblurred image, which we write as the set of convolutions {d (q) * f }. In this notation, the set of filters {d (q) } = {d x , d y } are 1 st order horizontal and vertical finite-difference derivative filters. While it is possible to extend the regularisation to higher-order derivatives d xx , d xy , d yy , etc. this is not generally done in practice. This kind of regularisation can be written as ρ f (f ) = Φ {d (q) * f } , where the function Φ penalises the derivatives of the deblurred image. Using priors of the form in Equation (2.26), Φ has the form of a sum of 1D functions of each derivative at each pixel: Φ({v (q) }) = q j φ(v (q) j ). [START_REF] Levin | Image and depth from a conventional camera with a coded aperture[END_REF] use φ(x) = |x| p with p = 0.8, and use iteratively-reweighted least squares (IRLS) to minimise Equation (2.45) with this regulariser. This method involves solving a series of weighted least-squares problems, e.g. using conjugate gradient descent (CG).

Although the minimisation problem is non-convex for p < 1, the algorithm typically produces good results.

One problem with the IRLS method is that it typically takes a long time to run, sometimes tens of minutes for a megapixel image, as the conjugate gradient descent must be performed multiple times. To mitigate this, several authors have proposed to use so called "half-quadratic" optimisation schemes, whereby auxiliary variables are introduced to represent the derivatives of the deblurred image, and an alternating minimisation can be performed to obtain the deblurred result quickly [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF]. For the sparse gradient prior, a pair of auxiliary images {v (q) } = {v x , v y } is introduced, corresponding to the derivatives of f , and Equation (2.45) (where f is regularised directly) is replaced with the following problem, wherein f is indirectly regularised via the auxiliary variables:

min f ,{v (q) } Af -g 2 2 + λ q d (q) * f -v (q) 2 2 + αΦ {v (q) } . (2.46)
For sufficiently large λ, this problem is approximately equal to the standard regularised least squares problem in Equation (2.45) with ρ f (f ) = Φ {d (q) * f } . The optimisation alternates between solving for f while holding {v (q) } fixed, and vice versa, gradually increasing λ. For fixed {v (q) }, the problem is quadratic in f (hence the name "halfquadratic"), and for spatially-invariant blur can be solved directly in the frequency domain:

f = F -1 F (a) * • F (g) + λ q F d (q) * • F v (q) F (a) * • F (a) + λ q F d (q) * • F d (q) .
(2.47)

Holding f fixed, each optimal v (q) can be found by 2N independent 1D minimisations (1 for each pixel of each auxiliary image), of the form

min v (q) j λ (d (q) * f ) j -v (q) j 2 + αφ v (q) j .
(2.48)

For the hyper-Laplacian prior, where φ(x) = |x| p , [START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF] demonstrated analytic solutions to these minimisations for particular exponents p = 1 2 and p = 2 3 , and showed how to use a lookup table to solve this problem quickly for arbitrary exponents. [START_REF] Shan | High-quality motion deblurring from a single image[END_REF] also use the half-quadratic method for their piecewise prior φ(x) = min(k |x| , ax 2 + b). Using this choice for φ, Equation (2.48) can be minimised separately for each piece of φ, before taking the minimum amongst them. The deblurred results are produced using the algorithms of (c) [START_REF] Levin | Image and depth from a conventional camera with a coded aperture[END_REF], and (d) [START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF].

In a related approach, [START_REF] Wang | A new alternating minimization algorithm for total variation image reconstruction[END_REF] use the half-quadratic method to apply total variation regularisation to the deblurring problem, where Φ({v (q) }) does not penalise each derivative independently but couples them such that Φ({v

(q) }) = j v x j 2 + v y j 2 .
In this case, the optimal {v (q) } are found as the solutions of N independent 2D minimisations, and can still be found in closed-form. Recently, [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] have applied the half-quadratic method with a regulariser that models the statistics of small patches of the image rather than the derivatives. Several authors have investigated variants of these split-variable approaches, and accelerated algorithms for solving them [START_REF] Afonso | Fast image recovery using variable splitting and constrained optimization[END_REF][START_REF] Osher | Fast linearized bregman iteration for compressive sensing and sparse denoising[END_REF].

Non-Blind Deblurring With Other Noise Models

As seen in previous sections, the data fidelity term used for non-blind deblurring can be derived as the negative log-likelihood for a particular noise distribution. However, when the appropriate noise model is not known a priori , or there are outliers to the real noise distribution in the blurry image, artefacts can appear in the deblurred image.

In pursuit of deblurring algorithms that are robust against such problems other authors have proposed algorithms using robust data terms. These data fidelity terms reduce the risk of over-penalising a large data fidelity error. For examples, Xu and Jia (2010) use an ℓ 1 data term in non-blind deblurring to be robust to non-Gaussian noise, with TV regularisation, while [START_REF] Welk | Robust variational approaches to positivity-constrained image deconvolution[END_REF] proposes a robust version of the Richardson-Lucy algorithm.

Algorithms for Blind PSF Estimation

In general, the PSF for a blurry image may not be known. For some causes of blur, e.g. optical aberrations, the parameters may be estimated once per camera, however for others such as camera shake, the PSF is unique to each image and must be estimated on a per-image basis. "Blind" estimation of the PSF parameters directly from images has a long history, including methods based on entropy maximisation [START_REF] Gull | Maximum entropy method in image processing[END_REF], iterative projection on convex sets [START_REF] Ayers | Iterative blind deconvolution method and its applications[END_REF], and modification of the Richardson-Lucy algorithm [START_REF] Fish | Blind deconvolution by means of the Richardson-Lucy algorithm[END_REF]. The particular task of deblurring "shaken" images has received considerable attention recently. [START_REF] Fergus | Removing camera shake from a single photograph[END_REF] proposed the first robust algorithm for blind estimation of a camera shake blur kernel directly from a single blurry image, using an approximate marginalisation approach based on the method of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF] (also closely related to the more general blind deblurring algorithm of [START_REF] Molina | Blind deconvolution using a variational approach to parameter, image, and blur estimation[END_REF]). [START_REF] Shan | High-quality motion deblurring from a single image[END_REF] and [START_REF] Cho | Fast motion deblurring[END_REF] subsequently proposed maximum a posteriori (MAP) type algorithms capable of tackling this difficult problem in a much shorter amount of time. Several authors have also proposed similar MAP-type algorithms [START_REF] Cai | Blind motion deblurring from a single image using sparse approximation[END_REF][START_REF] Gupta | Single image deblurring using motion density functions[END_REF][START_REF] Krishnan | Blind deconvolution using a normalized sparsity measure[END_REF]Xu and Jia, 2010), while most recently [START_REF] Levin | Efficient marginal likelihood optimization in blind deconvolution[END_REF] propose a more efficient formulation of the marginalisation approach.

While non-blind deblurring may be posed as a convex problem, or at least a wellconditioned one, the problem of finding the sharp image and PSF that best reconstruct the observed image is in general ill-posed, since we have more unknowns than equations.

In fact, for a given blurry image g, there are an infinite number of sharp image / PSF pairs that can reconstruct g equally well, as demonstrated in Figure 2.7. To obtain a useful solution, it is necessary to add some regularisation and/or constraints on both the sharp image and the PSF.

For now, we will assume that the blur is spatially-invariant, such that the PSF A is parameterised by the single convolution kernel a.

Single-Image PSF Estimation

As discussed in Section 2.1.5, good priors for f and a are necessary for blind deblurring to be successful, so approaches to single-image PSF estimation generally take the posterior distribution for f and a as their starting point:

p(f , a|g) ∝ p(g|f , a)p(f )p(a), (2.49)
where the likelihood is derived from an assumption of uniform Gaussian noise, as discussed in Section 2.1.4:

p(g|f , a) ∝ i exp - ((a * f ) i -g i ) 2 2σ 2 , (2.50)
where σ is the standard deviation of the noise, and the priors p(f ) and p(a) may be taken from those discussed in to find the values of f and a that jointly maximise p(f , a|g). In this section we will review two popular approaches to blind PSF estimation, which follow these two broad directions. Figure 2.8 demonstrates an example result deblurring an image with these two algorithms.

The Marginalisation Approach

Fergus et al. (2006) proposed a marginalisation-based algorithm for blind deconvolution of a single image, blurred by camera shake. 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.

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 , a|g) by a simpler, factorised distribution. The factorised form of this distribution means that it is straightforward to marginalise out the sharp image f and find the expected value of a. [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.

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. They use the mixture of exponentials prior on the kernel described in Equation (2.27), and the mixture of Gaussians prior from Equation (2.25) on the sharp image. 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 }. 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 denote the set of latent variables f ∂ , a, and β σ as the "ensemble" Θ. The aim is to find the factorised distribution

q(Θ) = q(β σ ) j q(f ∂ j ) k q(a k ) (2.51)
that best approximates the true posterior p(Θ|g). Note that in this context, the latent image pixels f ∂ j , kernel elements a k , and noise precision β σ are random variables. Since the distribution is factorised, the posterior for a, p(a|g) is approximated simply by q(a). Thus, having found the optimal q(Θ), the expected value of a under q(a) is taken to be the optimal blur kernel, i.e., â = a q(a) . Since q(a) is itself factorised, i.e. q(a) = k q(a k ), this expectation is generally trivial to compute. Fergus et al. discard the latent image distribution q(f ) and estimate the sharp image using a nonblind deblurring algorithm (although as shown in Figure 4.2 (d), f q(f ) does in fact provide a useful estimate of the sharp image).

To find the optimal q(Θ), Miskin and MacKay minimise the following functional (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Θ. (2.52)
Minimising this functional is equivalent to minimising the Kullback-Leibler divergence between the posterior and the approximating distribution [START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF], and this is tackled by first using the calculus of variations to derive the optimal forms of q(f ∂ j ), q(a k ) and q(β σ ), then iteratively optimising their parameters. We omit a complete explanation of the details of the optimisation, however the main feature of the iterative process is as follows: Each distribution, for example q(a k ), is parameterised by two sufficient statistics, which we will denote a (1) k anda (2) k , essentially a mean and variance for a k . [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF] first derive equations for the optimal values of every unknown (i.e. a

(1)

k , a (2) k , f (1) k , f (2)
k , etc.) in terms of the other unknowns. By applying these optimality conditions as fixed-point update rules, [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF] aim to decrease the Kullback-Leibler divergence between the true and approximate posterior.

The Maximum a Posteriori Approach

In contrast to the marginalisation approach of the previous section, several authors have proposed algorithms which estimate a convolution kernel a by maximising the posterior p(f , a|g) [START_REF] Cai | Blind motion deblurring from a single image using sparse approximation[END_REF][START_REF] Cho | Fast motion deblurring[END_REF][START_REF] Shan | High-quality motion deblurring from a single image[END_REF]Xu and Jia, 2010). These maximum a posteriori (MAP) approaches vary in the priors they use, but generally follow the same structure, alternating between estimating a and f and proceeding from coarse to fine in a multi-scale framework.

Updating the blur kernel

Cho and [START_REF] Cho | Fast motion deblurring[END_REF] propose a method of kernel estimation which uses simple Gaussian priors on the kernel elements. Simply optimising a using this simple prior would almost certainly fail to produce a reasonable result, since it is not realistic. 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 2.1.5.

To update the kernel, [START_REF] Cho | Fast motion deblurring[END_REF] first filter the current estimate f of the sharp image to predict step edges. A bilateral filter [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF] is applied to denoise the image, followed by a shock filter [START_REF] Olshausen | Emergence of simple-cell receptive field properties by learning a sparse code for natural images[END_REF] to enhance sharp edges, producing an image f ′ . The derivatives d x * f ′ and d y * f ′ of this filtered image are then computed. Finally, these derivatives are thresholded to produce sparse derivatives maps p x and p y . The threshold is chosen so as to retain only a small number of non-zero entries, while ensuring that sufficient non-zeros are retained from all edge orientations.

Having predicted and retained salient sharp edges from the latent image, the set of sparse derivatives maps is augmented with higher order derivatives to form the set

{p (q) } = {p x , p y , d x * p x , 1 2 (d x * p y + d y * p x ), d y * p y },
ω q a * p (q) -d (q) * g 2 2 + β a 2 2 , (2.53)
where d (q) ∈ {d x , d y , d xx , d xy , d yy }, β is a weight for the Tykhonov regularisation, and the weights ω q weight each partial derivative. This is simply a linear least squares problem, which can be solved directly in the frequency domain, analogously to Equation (2.47):

â = F -1 q w q F p (q) * • F g (q)
q w q F p (q) * • F p (q) + β .

(2.54)

Having found the kernel â which minimises Equation (2.53), 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.

Xu and Jia (2010) take essentially the same approach to kernel estimation as [START_REF] Cho | Fast motion deblurring[END_REF], with the addition of an extra filtering step aimed at suppressing unhelpful non-zeros in the sparse derivative maps {p x , p y }. Xu and Jia show that non-zeros in p x and p y that arise from narrow peaks in the image (as opposed to step edges) actually reduce the quality of the estimated kernel, and as such should be removed before solving Equation (2.53).

In a contrasting approach, [START_REF] Shan | High-quality motion deblurring from a single image[END_REF] do not perform any image filtering, and replace the Tykhonov regularisation in Equation (2.53) with ℓ 1 regularisation, solving min

a q ω q a * d (q) * f -d (q) * g 2 2 + β a 1 . (2.55)
Here, instead of the Tykhonov regularisation used by Cho and Lee, ℓ 1 regularisation is used, corresponding to an exponential prior and encouraging sparsity directly. Using this regulariser however, precludes direct minimisation in the frequency domain, and Shan et al. use an interior point method [START_REF] Kim | An interior-point method for large-scale ℓ 1 -regularized least squares[END_REF] to solve Equation (2.55).

Updating the latent image

In order to update the latent image, the current estimate of the blur kernel â is used to deconvolve the blurry image and obtain an improved estimate of the sharp image. [START_REF] Cho | Fast motion deblurring[END_REF] do this by solving min

f q ω q â * d (q) * f -d (q) * g 2 2 + α d x * f 2 2 + α d y * f 2 2 , (2.56)
where α is a regularisation weight, and now, the partial derivatives include the zeroth order:

d (q) ∈ {1, d x , d y , d xx , d xy , d yy }.
Note that in this problem, the data term involves higher-order derivatives of the image, although the regularisation does not. This is solved in closed-form by pixel-wise division in the frequency domain. [START_REF] Shan | High-quality motion deblurring from a single image[END_REF] update the latent image using their non-blind deblurring algorithm mentioned in Section 2.2.2, which uses a sparsity prior on the image gradients. Xu and Jia (2010) choose to add ℓ 2 regularisation that encourages the latent image gradients to be close to the sparse derivative maps {p (q) } predicted from the previous iteration:

min f q ω q â * d (q) * f -d (q) * g 2 2 + α d x * f -p x 2 2 + α d y * f -p y 2 2 .
(2.57)

This reduces the problem of over-smoothing caused by the ℓ 2 regularisation, since sharp edges are encouraged to appear where they exist in {p (q) }, while still permitting closedform minimisation in the frequency domain.

As shown by [START_REF] Shan | High-quality motion deblurring from a single image[END_REF], the appropriate values for the relative weights of the different orders of partial derivatives are ω q = 1 2 κ(q) ω 0 , where κ(q) denotes the order of derivative d q .

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 5-7 iterations are typically sufficient.

Clearly the effectiveness of the MAP approach is not dependent on the exact priors used to regularise the least-squares estimation steps, since [START_REF] Cho | Fast motion deblurring[END_REF] use very weak Gaussian priors. However, what is important is that sharp edges are encouraged to appear in the latent image, and that the kernel is encouraged to be sparse. While [START_REF] Shan | High-quality motion deblurring from a single image[END_REF] do indeed achieve this using regularisation, [START_REF] Cho | Fast motion deblurring[END_REF] avoid the costly optimisation problems and achieve the same aim instead with simple non-linear filtering and thresholding steps.

Recently, [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF] have analysed the differences between the marginalisation and MAP approaches to blind deblurring, and conclude that in order to estimate the PSF, it is advantageous to marginalise over the sharp image f , while maximising over the kernel a. In a subsequent work, [START_REF] Levin | Efficient marginal likelihood optimization in blind deconvolution[END_REF] propose an algorithm which does this explicitly, maximising the posterior for the kernel alone p(a|g), which is obtained by marginalising over the unknown sharp image f , i.e. p(a|g) = p(a, f |g) df . To do this, [START_REF] Levin | Efficient marginal likelihood optimization in blind deconvolution[END_REF] adopt a strategy that is similar in spirit to the marginalisation algorithm of [START_REF] Fergus | Removing camera shake from a single photograph[END_REF], described in the previous section. Yuan et al. (2007b) proposed a method for blind deblurring where, in addition to a blurry image g, a sharp but noisy image f N of the same scene is also available. 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 by denoising, it can initially be used as a proxy for the sharp image, allowing us to estimate the PSF. Following this, with the PSF fixed, the noisy image can also be used to improve the final non-blind deconvolution step, using a modified version of the Richardson-Lucy algorithm, which uses f N to suppress ringing artefacts in the result.

Deblurring With Noisy / Blurry Image Pairs

In order to estimate the blur kernel, Yuan et al. (2007b) solve

min a a * f N -g 2 2 + β a 2 2 .
(2.58)

Having estimated the PSF, Yuan et al. propose several modifications to the Richardson-Lucy (RL) algorithm to perform the non-blind deblurring, 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" artefacts -low-frequency ripples spreading across the image, such as in Figure 2.9 (c)

-but using the denoised image it is possible to disambiguate the true low frequencies from these artefacts, 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., 2007b) for full details of the augmented RL algorithm. Figure 2.9

shows an example result deblurred using this approach.

Modelling Spatially-Variant Blur

Until this point, we have only considered spatially-invariant blur, where the image formation model in Equation (3.13) can be written as a 2D convolution. In most of the work mentioned above, and in general through the deblurring literature, convolution is used as a simple model for image blur. In general however, the blur in an image may be spatially-variant, and the matrix A will contain a different PSF for each pixel of the image. Figure 2.10 shows an example of spatially-variant blur due to camera shake.

For real images of more than a few hundred pixels across, it is not practical to compute and store every element of the large matrix A explicitly, and instead a more compact parameterisation of the matrix is required. For spatially-invariant blur, such a parameterisation comes in the form of a convolution kernel. For spatially-variant blur, recent work has investigated analogous parameterisations of the matrix A (a notable exception is the work of [START_REF] Seitz | Filter flow[END_REF], who attempt to estimate the elements of A itself, with no underlying parameterisation). The first distinction to draw is between spatially-variant blur caused by moving objects in the scene, and that caused by the image capture process. When the blur is caused by moving objects, the spatial-variance is largely due to the fact that different objects move differently, and the blur changes abruptly at object boundaries. Within each object, simple linear blurs are often assumed [START_REF] Chakrabarti | Analyzing spatially-varying blur[END_REF][START_REF] Cho | Removing non-uniform motion blur from images[END_REF][START_REF] Levin | Blind motion deblurring using image statistics[END_REF], such that the matrix A is parameterised by a line segment representing the blur, and a segmentation delineating the object.

When the blur is caused by the image capture process, rather than scene motion, the blur is more likely to vary smoothly across the image. Previous approaches to modelling this kind of blur can largely be divided into two categories: global models, which use a single set of parameters to describe the blur at all points in the image, and local models, which divide the image into a set of regions and use a separate set of blur parameters for each region. Although in this thesis we are concerned with blur caused by camera motion, we note that models for spatially-variant blur caused by optical aberrations of the camera have also recently been proposed by [START_REF] Kee | Modeling and removing spatially-varying optical blur[END_REF][START_REF] Schuler | Non-stationary correction of optical aberrations[END_REF].

Global Models

On of the contributions of this thesis (see Chapter 3) is to propose a model for spataillyvarying camera shake blur. Concurrently with this work, several models have been proposed to parameterise spatially-variant blur caused by unconstrained camera motion [START_REF] Gupta | Single image deblurring using motion density functions[END_REF][START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF][START_REF] Tai | Richardson-Lucy deblurring for scenes under a projective motion path[END_REF]. These approaches are derived from the relationship between motion of the camera and transformations of the sharp image. The general form of these models is

g * = K k=1 w k T (k) f , (2.59)
where each matrix T (k) applies some transformation to the sharp image, which then contributes to the blurry image with a weight w k . In this case, the weights (which we concatenate into the K element vector w) provide the compact parameterisation of the PSF.

The models proposed by various authors vary in their choice of the set of image transformations {T (k) }. [START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF] model the image motion using homographies which correspond to a planar scene. These homographies are parameterised by the pose of the camera, which is recorded using inertial measurement sensors attached to the camera. [START_REF] Tai | Richardson-Lucy deblurring for scenes under a projective motion path[END_REF] also model the image motion using homographies, although without constraining them to correspond to scene planes or camera rotations. In this case, the homographies are estimated by applying optical flow to the output of a second, high-speed camera attached to the main camera. Both [START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF] and [START_REF] Tai | Richardson-Lucy deblurring for scenes under a projective motion path[END_REF] apply their models only to non-blind deblurring, where the parameters w are estimated using additional hardware attached to the camera. [START_REF] Gupta | Single image deblurring using motion density functions[END_REF] propose to use the set of image-plane translations and rotations, and estimate the weights w directly from the blurry image, i.e. in a blind deblurring setting.

While the models described above are largely agnostic of the scene content, [START_REF] Šorel | Space-variant restoration of images degraded by camera motion blur[END_REF] tackled simultaneous blind deblurring and depth map estimation from multiple blurry images. Global models of spatially-variant blur have been also applied under various constrained motion models [START_REF] Klein | A single-frame visual gyroscope[END_REF][START_REF] Sawchuk | Space-variant image restoration by coordinate transformations[END_REF][START_REF] Shan | Rotational motion deblurring of a rigid object from a single image[END_REF]Tai et al., 2010a), and although global descriptors are used to parameterise these continuously-varying blurs, the constraints are often restrictive.

In Chapter 3 we propose a global model for camera shake blur with the form in Equation (2.59). Our method is distinguished from previous approaches by its foundation in the geometry of camera rotations, which provides a minimal parameterisation of the blur while remaining accurate at all possible camera settings.

Local Models

An alternative method for parameterising a spatially-variant blur matrix A is to assume that instead of each pixel having a unique PSF, nearby pixels can (approximately) share a single PSF. In this way, a spatially-variant blur can be parameterised using only a small number of convolution kernels, each assigned to a different region of the image. In contrast to "global" models such as those mentioned above, which use a complex motion model and a single set of parameters to describe the blur at every point in the image, these approaches are local, using a simple motion model but separate sets of parameters for each region in the image. Much previous work has relied on this idea, modelling the blur as being locally uniform [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF]Nagy and O'Leary, 1998;Tai et al., 2010b;[START_REF] Vio | Multiple image deblurring with spatially variant PSFs[END_REF].

The two general forms for these local models are

g * = P r=1 m (r) • (a (r) * f ) (2.60) g * = P r=1 a (r) * (m (r) • f ), (2.61)
where the blur within each region r is modelled using the convolution filter a (r) and a windowing function m (r) . The windowing functions m (r) are usually fixed manually in advance, and the set of kernels {a (r) } parameterises the PSF. In the first form, m (r) controls the region of influence of the filter a (r) (i.e. which parts of the blurry image are affected by the filter a (r) ), while in the second form, m (r) controls its region of dependence (i.e. which parts of the sharp image are blurred with the filter a (r) ). This locally-uniform method has proven to be powerful enough to model complex spatially-variant blur with reasonable accuracy, while enjoying computational advantages over global models, since the FFT can be used to compute the convolutions efficiently. In both cases, the windows r) are typically chosen to sum to unity everywhere, i.e. r m (r) = 1, to ensure that the blur is conservative.

m (
Nagy and O'Leary (1998) consider local blur models of the form in Equation (2.60), using either rectangular or triangular regions, for both the case where the regions are disjoint, and where they overlap. Tai et al. (2010b) also adopt a spatially-variant blur model of the form in Equation (2.60). However, instead of fixing the regions of influence r) and estimating the local kernels a (r) which best fit the observed blur, they fix the kernels a (r) using a pre-determined basis, and estimate the regions of influence m (r) which best reconstruct the observed blur. r) for each region r. This gives a smooth fall-off from one region to the next, allowing the blur to vary smoothly across the image, instead of changing abruptly at region boundaries. The fact that each rectangular region is small and has a spatially-invariant blur allows the forward mode to be computed using P small convolutions. To take advantage of the small support of each region, it is possible to define a matrix C (r) to crop out the r th rectangular region, and re-write the forward model in (2.61) as

m (
g * = P r=1 C (r) ⊤ a (r) * (m • C (r) f ) . (2.62)
Here, all regions use the same Bartlett-Hann windowing function m, and we note that the matrix C (r) ⊤ re-inserts the r th region at its correct location. This model can be computed efficiently in the frequency domain:

g * = P r=1 C (r) ⊤ F -1 F a (r) • F m • C (r) f , (2.63)
where F ( • ) takes the 2D discrete Fourier transform (computed using the fast Fourier transform), and F -1 ( • ) the inverse Fourier transform.

This region-based approach to convolution was originally proposed by [START_REF] Stockham | Image completion with structure propagation[END_REF], as a means of computing spatially-invariant convolutions and correlations more efficiently. When the filters are all identical, i.e. a (r) = a for all r and have small support, then the following equality holds exactly:

P r=1 C (r) ⊤ a * (m • C (r) f ) = a * f . (2.64)
The two main advantages of computing a spatially-invariant convolution in sections is that the entire input signal f does not all need to be stored in memory at once, and that it is cheaper to compute several small FFTs than to compute one large FFT. The FFT has computational complexity O N log N for a signal of length N , and thus computing two

FFTs of size N/2 has complexity O 2N/2 log N/2 = O N log N/2
. This advantage in speed also holds true for the spatially-variant case, allowing spatially-variant blur to be computed in a similar amount of time as spatially-invariant blur.

Say how the priors used on spatially-invariant kernels a can sometimes apply to the spatially-variant parameterisation.

Inpainting

The task of replacing missing or undesirable pixels in an image is known as "inpainting" [START_REF] Bertalmío | Image inpainting[END_REF], or sometimes "image completion" or "scene completion". When we are concerned with replacing missing or deleted pixels in an image, the restoration problem is quite different from deblurring, where all of the latent pixels are, in some sense, observed. Here, the missing latent pixels are completely decoupled from the observed image, and so cannot be recovered at all using maximum-likelihood techniques which estimate the latent image from the data alone. Instead, we require informative priors on the latent image.

In the problem of inpainting, we are generally not concerned with blur, so the PSF is simpyl the identity, i.e. A = I. The occlusion / corruption that we wish to replace is modelled as noise with a uniform distribution, and the other observed pixels are assumed to be noiseless, i.e. exactly equal to their counterparts in the latent image. We refer to the unknown / corrupted region to be filled as the target region, and denote this by the set Ψ, such that g i ∼ Unif(0, 256) for i ∈ Ψ. The known / good region is denoted by the set M, and with no noise, we have f M = g M , as shown in Figure 2.11. The pixels in Ψ which lie along the boundary with M are denoted by δΨ. In a probabilistic formulation,

the likelihood here provides no information about the unknown pixels f Ψ , as they are decoupled from the known pixels in M. Thus we require some prior information about the latent image before we can begin to replace the missing pixels.

Previous work on the problem of inpainting is diverse, varying significantly according to the size of the region to be filled, the content of the image surrounding it, and the amount of user interaction required.

Early work on the problem relied on the notion that the latent image should belong to certain convex sets. If projection operators P for these sets are available, then identities of the form f = Pf can be used to derive inpainting algorithms which find a solution which lies in the intersection of all the sets. Ferreira and Pinho (1994) used this approach to solve the inpainting problem as interpolation of a band-limited signal, and proposed an algorithm of this type capable of filling small holes, while [START_REF] Hirani | Combining frequency and spatial domain information for fast interactive image noise removal[END_REF] proposed a projection-based method capable of filling a textured region containing straight-line edge structures, by constraining the power spectrum of the region to match a manually-selected reference region.

The problem of filling small regions has generally been approached using local models of images, i.e. models which account for the behaviour of images in areas a few pixels across. [START_REF] Bertalmío | Image inpainting[END_REF] approached the problem using the equations of fluid dynamics, following the intuition that lines of constant colour (isophotes) should continue into the region to be filled and connect smoothly within it. Their algorithm iteratively propagates colour into the target region from around the boundary, leading to a smooth continuation of colours.Recent work has focussed on the use of learnt priors on the statistics of local patches (Mairal et al., 2008;[START_REF] Roth | Fields of experts: a framework for learning image priors[END_REF][START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]. The strength of these approaches is that the local structure in natural images is reasonably predictable, such that a pixels' value is closely related to its neighbours' values. The general form for these approaches is essentially a MAP problem with no data fidelity term (since the observed pixels contain no information about the missing pixels):

min f -log p(f ) (2.65)
In this formulation the unknown pixels are decoupled from the observed pixels in the PSF. The priors provide this coupling, and essentially provide a way of producing statistically probable (and thus hopefully convincing) hallucinations for the filled pixels. By looking at the pixels around the missing region, the priors are able to "suggest" likely content in the region.

Except in some cases, these local approaches fail when the region is larger than a few pixels across. When the pixels in the middle of the region are remote from those at the edge, it is no longer possible to infer their likely content from the surrounding pixels. One exception is when filling a region with homogeneous texture, for which a local texture synthesis algorithm can be applied [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF].

When filling a region with large-scale structures however, purely-local methods fail to capture both the structure and the fine-scale texture. To cope with this, subsequent work has considered combining texture synthesis methods with higher-level structural information about strong edges inside the region. This structural information is either inferred from the surrounding image [START_REF] Bertalmío | Simultaneous structure and texture image inpainting[END_REF][START_REF] Criminisi | Object removal by exemplar-based inpainting[END_REF][START_REF] Jia | Image repairing: robust image synthesis by adaptive N D tensor voting[END_REF], or manually provided by the user (Sun et al., 2005). In all these works, the texture synthesis algorithm replaces individual pixels or patches one at a time, selecting a pixel (or patch) from δΨ using some criteria, then replacing it. By combining texture synthesis with information about the large-scale content of the region, more complex image content can be produced. However, these algorithms are still constrained to draw their information from the image itself, or from simple manually-input guides.

Thus there is a limit to the kind of regions they can reconstruct.

While methods of texture-synthesis and edge-continuation aim to hallucinate realistic low-level details of an image and simple structures, a separate line of work has investigated synthesising images from multiple input images of the same scene. [START_REF] Agarwala | Interactive digital photomontage[END_REF][START_REF] Rasmussen | Spatiotemporal inpainting for recovering texture maps of partially occluded building facades[END_REF][START_REF] Wilczkowiak | Two-phase kernel estimation for robust motion deblurring[END_REF] propose methods for combining multiple images, captured in quick succession by the same camera. An integral component of these methods is "Poisson blending" [START_REF] Pérez | Poisson image editing[END_REF], a method for seamlessly stitching together regions from different images without introducing visible boundaries.

On the other hand, there has been significant interest recently in organising, visualising, and searching the growing number of images available in online photo collections [START_REF] Chum | Total recall: automatic query expansion with a generative feature model for object retrieval[END_REF][START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF][START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF][START_REF] Snavely | Photo tourism: exploring photo collections in 3D[END_REF].

For popular locations, several authors have demonstrated that there exist sufficient images to perform accurate 3D reconstructions [START_REF] Furukawa | Towards internet-scale multi-view stereo[END_REF][START_REF] Goesele | Multi-view stereo for community photo collections[END_REF] and 3D scene visualisations [START_REF] Shahrokni | From visual query to visual portrayal[END_REF]. [START_REF] Hays | Scene completion using millions of photographs[END_REF] first proposed to exploit this abundance of information to perform inpainting, using a large photo collection retrieved from the Internet to remove large objects from photographs. Their algorithm finds images with similar semantic content to that being restored, and directly copies the contents into the target image using Poisson blending. While the results are often convincing, they do not necessarily correspond to the true underlying scene, since the source images may depict an entirely different part of the world. [START_REF] Gamelin | The dimensionality of scene appearance[END_REF] also demonstrated the use of online photo collections for inpainting, by performing a 51 dense 3D reconstruction of the scene and learning an appearance model for its surface.

These approaches to inpainting can be distinguished from the approaches discussed above by their use of direct exemplars, as opposed to statistical priors on the local behaviour of natural images. This idea goes back at least to [START_REF] Hirani | Combining frequency and spatial domain information for fast interactive image noise removal[END_REF], who require the user to select a "reference" patch, which is then used as a constraint to ensure that the filled region of the latent image to has a similar spectrum to the reference region.

In Chapter 7 we propose a system which builds on this work, by retrieving reference images of the same scene as the image to be restored. We replace the occluded region using these exemplars, without an explicit 3D reconstruction of the scene or any knowledge of the camera positions. As opposed to filling the missing region with content synthesised from statistical information or semantically similar scenes, we use these reference images to produce a realistic output that corresponds to the true occluded scene.

Chapter 3

Modelling Spatially-Variant Camera Shake Blur

Introduction

Everybody is familiar with camera shake, since the resulting blur spoils many photos taken in low-light conditions. Often in such conditions, the shutter speed selected by the camera is slow, and unless the photographer is able to hold the camera perfectly still, the camera's motion causes the photograph to be blurry. While significant progress has been made recently towards removing such blur from images, most approaches to this problem assume that the blurred image can be modelled as the 2D convolution of a sharp image with a spatially-invariant 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., 2007b), as discussed in Section 2.3. However, real camera shake does not, in general, cause such spatially-invariant blur [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF]. Figure 3.1

shows an example of a typical spatially-variant blur caused by camera shake.

In this chapter we propose a geometrically motivated model of spatially-variant image blur due to camera shake. We develop a global descriptor for such blur, analogous to a convolution kernel, and demonstrate its ability to model a more general class of blurs than previous approaches, including spatially-invariant blur as a special case. In this chapter, we limit our scope to describing and discussing the model itself. We defer the application of the model to the following and later chapters.

We begin by considering the relative blurring effect of different camera motions, and deriving a geometric model for camera shake in Section 3.2. In Section 3.3 we develop this into a practical model for deblurring real images, and provide a comparison to other related models proposed recently. We conclude in Section 3.4.

A Geometric Model for Camera Shake

When we say that a photograph is blurry due to camera shake, it is generally understood that the camera moved during the exposure, and that this motion is to blame for the blur. More precisely, while the shutter is open the camera passes through a sequence of different poses, each of which gives a different view of the scene. The sensor accumulates all of these views, summing them up to form the recorded image. Intuitively, summing up many different views of a scene will result in a blurry image, and the more the views differ, the larger the blur will be. However, a camera may move in several different ways, and it is not necessarily obvious which kinds of motion cause large changes in the view (and hence a large blur), and which cause relatively small changes. Furthermore, even if the camera's motion is fully known for a given photograph, we need a model to translate this physical 3D motion into image-domain motion before we can begin to deblur the photograph.

In this section we address the challenges mentioned above. We examine the causes of camera shake blur, and derive a geometric model for the process, which fits into the linear degradation model presented in the previous chapter. In this work we limit our scope to photographs of static scenes, i.e. the blur is solely due to the motion of the camera.

Components of Camera Motion

To begin our discussion of camera motion, we note that the pose of a camera incorporates two components: position and orientation. Intuitively, the position tells us where the camera is, while the orientation tells us which way it is pointing. Both of these components may vary while the camera's shutter is open. In this section, we discuss the contribution of each component to the image blur, and conclude that in most cases of camera shake, the changes in orientation (rotation) of the camera during exposure have a significantly larger blurring effect than the changes in position (translation).

Consider the simplified case shown in Figure 3.2 of a scene point P , at a distance D from the camera, being imaged at the centre of the camera's retina / sensor. During the exposure the image of the point is blurred through a distance δ pixels. In (a) the camera translates through a distance X parallel to the image plane, while in (b) the camera rotates through an angle θ about its optical centre. By simple trigonometry, we can see that in (a) the camera must translate by

X = δ F D, (3.1)
where F is the camera's focal length, while in (b) the camera must rotate through an angle

θ = tan -1 δ F . (3.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 Equation (3.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 Equation (3.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 56 3.2 A Geometric Model for Camera Shake 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).

motion is at least an order of magnitude smaller than for a translation of the camera causing an equivalent amount of blur.

In reality, both the position and orientation of the camera vary simultaneously during the exposure. However, if 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. Figure 3.3 supports this assumption, showing some shaken images captured by [START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF],

for which the motion of the camera during the exposure was recorded. As can be seen, the camera translates by a few millimetres, while rotating by up to 0.5 • . For these magnitudes of motion, the blur for subjects more than 1 m away from the camera can almost entirely be attributed to the rotation.

From now on, we assume that the translational component of camera motion does not cause any blur. Furthermore, we assume that all rotations occur about the camera's optical centre. Note, however, that a camera rotation about a centre that is not the optical centre can be written as a rotation about the optical centre composed with a translation; these translations should generally be small for rotation centres that are not far from the optical centre.

Motion Blur and Homographies

It is well known that under a pinhole camera model, and assuming that the scene being photographed is static, rotations of a camera about its optical centre induce projective 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 centre, and the pixels are assumed to be square.

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.3)
where the 3 × 3 matrix R is a rotation matrix describing the orientation of the camera, and K is the camera's internal calibration matrix [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. In this work, we assume that the calibration matrix K is known (see Section 3.2.3).

The rotation matrix R has only 3 parameters. We adopt here the "angle-axis" parameterisation, in which a rotation is described by the 3-vector θ = (θ X , θ Y , θ Z ). R is then given by the matrix exponential

R θ = e [θ] × , where (3.4) [θ] × =      0 -θ Z θ Y θ Z 0 -θ X -θ Y θ X 0      . (3.5)
We fix our 3D coordinate frame to have its origin at the camera's optical centre. 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.4 (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.

g * (x) = T 0 f H t x dt, (3.6)
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 . Under this model, the apparent motion of scene points may vary significantly across the image. Figure 3.4 demonstrates this, showing the paths followed by points in an image under a Y -axis rotation, a Z-axis rotation, or an arbitrary sequence of rotations 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 nonuniformity 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.4. In addition, it is common for camera shake to include an in-plane rotational motion. From this, it is clear that modelling camera shake as a convolution with a spatially-invariant kernel is insufficient to fully describe its effects (see also Figure 3.1).

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

g * (x) = f H θ x w(θ) dθ, (3.7)
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, and the value of the function at a point θ along the trajectory corresponds to the duration the camera spent at the orientation θ.

Camera Calibration

In order to compute the homography in Equation (3.3) that is induced by a particular rotation of the camera, we need to know the camera's calibration matrix K. For the results shown in this thesis, we assume that K takes the standard form

K =      F 0 x 0 0 F y 0 0 0 1      . (3.8)
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. 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 centre of the image.

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 software1 , 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.

In practice, there are several potential pitfalls in relying on the EXIF information to recover the internal calibration of the camera. For example, if the image has been cropped or resized, the focal length and pixel size recovered from the EXIF tags will no longer be correct. In addition, many image editors delete or overwrite the EXIF tags, and while almost all cameras include the focal length in the tags, some e.g. cameraphones, do not include sufficient information to estimate the pixel size. In the case that the EXIF tags are incorrect or incomplete, we set the focal length to the width of the image, and the principal point to the center of the image. Although this is unlikely to be an accurate estimate, it is typically sufficient for deblurring purposes.

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 (3.4) can be approximated by discarding the 2 nd and higher order terms:

R θ ≈      1 0 θ Y 0 1 -θ X -θ Y θ X 1      . (3.9)
Combining this with Equations (3.3) and (3.8), it can be shown that as F → ∞,

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

A Computational Model for Camera Shake

So far, our model has been defined in terms of the continuous functions f and g, and the weight function w. We discretise the camera orientation space into a 3D volumetric grid of size N X × N Y × N Z , and assign each orientation θ (k) a weight w k , for k ∈ {1, . . . , K},

where

K = N X N Y N Z .
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. 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 modelled as a weighted sum of a set of projectively-transformed versions of f :

g * = k w k T (k) f , (3.11)
where T (k) is the N × N matrix which applies homography H k to the sharp image f . The matrix T (k) is very sparse. For example, if bilinear interpolation is used when transforming the image, each row has only 4 non-zero elements. Section 3.3.2 describes how to calculate the matrix T (k) . By writing out Equation (3.11) for a single pixel, we obtain the discrete analog of Equation (3.7):

g * i = k w k j T (k) ij f j , (3.12)
where i and j index the pixels of the blurry image and the sharp image, respectively. For a blurry pixel g * i with coordinate vector x i , the sum j T (k) ij f j interpolates the value of f (H k x i ). Figure 3.5 shows an example of this, where a blurry pixel with homogeneous coordinate vector x i is mapped under a homography H k to the point H k x i in the sharp image. The value of f at the point H k x i is then interpolated as a weighted sum of the nearby pixels of f . Section 3.3.2 gives further details of this interpolation.

Due to the bilinear form of Equation (3.12), note that when either the blur kernel or the sharp image is known, the blurry image is linear in the remaining unknowns, i.e.

given w, g * = Af , where

A ij = k T (k) ij w k , (3.13) given f , g * = Bw, where B ik = j T (k) ij f j . (3.14)
In the first form, A ∈ R N ×N is a large sparse matrix, whose rows each contain a local

g i x i H k → f 1 f 3 f 2 f 4 H k x i
Blurry image g Sharp image f

Figure 3.5. Interpolation of sub-pixel locations in the sharp image. In general, a homography H k will not map a pixel (e.g. x i ) 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 i is interpolated as a weighted sum of nearby pixels. Using bilinear interpolation, the value of f at H k x i is computed as a linear combination of f 1 , f 2 , f 3 , and f 4 .

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

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 spatially-variant blur. In common with our approach, 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 image motion using a temporally-ordered sequence of unconstrained 8-parameter homographies, which is known in advance. [START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF] also use a temporally-ordered sequence of 6-parameter homographies which correspond to a fronto-parallel planar scene, parameterised by the camera's pose. [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. [START_REF] Gupta | Single image deblurring using motion density functions[END_REF] provide some analysis of the accuracy of their model as the focal length of the camera varies, and show that as the focal length decreases, their model becomes less accurate. Starting from our model, this is not unexpected, as it is possible to show that their model is a special case of ours, obtained by letting F → ∞. By following essentially the same derivation as in Section 3.2.4, except with θ Z = 0, it is possible to show that for small θ X and θ Y , as F → ∞, Of course, the need to know the focal length for an image is an additional requirement, compared to the model of [START_REF] Gupta | Single image deblurring using motion density functions[END_REF]. However, given that we recover their model by setting F to a large number, and that it is often possible to extract the focal length of a photo from the EXIF tags, we do not believe this is a significant disadvantage of our model.

H θ →      cos θ Z -sin θ Z F θ Y sin θ Z cos θ Z -F θ X 0 0 1      , ( 3

Computation of Interpolation Coefficients

Here we give details of how to calculate the matrix T (k) that applies a homography H k to an image, using bilinear interpolation. Note that these are standard interpolation weights, and are not specific to this deblurring application.

Using bilinear interpolation, the value of an image f at a sub-pixel location (x ′ , y ′ ) is interpolated from the 4 pixels surrounding (x ′ , y ′ ) as follows:

δ x = x ′ -⌊x ′ ⌋ (3.16) δ y = y ′ -⌊y ′ ⌋ (3.17) f (x ′ , y ′ ) = (1 -δ x )(1 -δ y ) f ⌊x ′ ⌋, ⌊y ′ ⌋ + (1 -δ x )δ y f ⌊x ′ ⌋, ⌊y ′ ⌋ + 1 + δ x (1 -δ y ) f ⌊x ′ ⌋ + 1, ⌊y ′ ⌋ + δ x δ y f ⌊x ′ ⌋ + 1, ⌊y ′ ⌋ + 1 (3.18)
where ⌊ • ⌋ takes the integer part of a positive scalar.

In Equation (3.12), we write this interpolation operation as the sum

f (H k x i ) = j T (k) ij f j . (3.19)
Using bilinear interpolation as described above, we can see that the elements T (k) ij take the following values:

j T (k) ij sub2ind(⌊x ′ ⌋, ⌊y ′ ⌋) (1 -δ x )(1 -δ y ) sub2ind(⌊x ′ ⌋, ⌊y ′ ⌋ + 1) (1 -δ x )δ y sub2ind(⌊x ′ ⌋ + 1, ⌊y ′ ⌋) δ x (1 -δ y )
sub2ind(⌊x ′ ⌋ + 1, ⌊y ′ ⌋ + 1) δ x δ y all others 0

where the Matlab function sub2ind provides the correspondence between the (x, y)

coordinates of a pixel and its index j.

Sampling the Set of Rotations

One important detail to consider is how finely to discretise 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 neighbours. Setting the grid spacing in terms of pixels also has the advantage that our 3D blur kernels are defined 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 • ].

Conclusion

In this chapter we have proposed a geometrically-derived model for blur caused by camera shake. For a static scene and a camera with known focal length, we have shown that the blur caused by camera rotation can be modelled using a weighted set of homographies, and have proposed a practical formulation of this model in which the blurry image is bilinear in the sharp image and the weights. The model assumes that the motion of the camera during exposure is limited to rotations about its optical centre, and is temporally-agnostic to the distribution over camera orientations. Our model is not applicable for non-static scenes, or nearby scenes with large camera translations where parallax effects may become significant.

Chapter 4

Estimating and Removing Spatially-Variant Camera Shake Blur

Introduction

In this chapter, we demonstrate the effectiveness of the spatially-variant blur model presented in Chapter 3 by using it to replace the uniform blur model in three existing approaches to camera shake removal. We show quantitative and qualitative improvements in the deblurred results.

In Section 4.3 we consider single-image deblurring, where only a blurry image is available, applying our model within the algorithms proposed by [START_REF] Fergus | Removing camera shake from a single photograph[END_REF] (Section 4.3.1) and Cho and Lee (2009) (Section 4.3.2). We compare our results to those of the original algorithms (which handle only spatially-invariant blur), and discuss the advantages and limitations of using our model. In Section 4.4 we apply our model to the case where an additional sharp but noisy image of the same scene is available, as proposed by Yuan et al. (2007b). In Section 4.5 we discuss some implementation considerations of using our model in these algorithms, before concluding in Section 4.6.

Application to Existing Deblurring Algorithms

The fact that Equation (3.12) 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 (3.12) without significant modification to the algorithm. Here we discuss some of the general issues related to achieving this.

The constraint that all the PSF parameters be positive is intuitively applicable to our model, since each kernel element w k corresponds directly to the length of time spent at a camera orientation θ k . If the camera passed through orientation θ k during the exposure, w k will be positive, and if not, w k = 0.

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

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, e.g. the Wiener filter [START_REF] Wiener | Extrapolation, Interpolation, and Smoothing of Stationary Time Series[END_REF]. 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 (see Chapter 5 for an efficient approximation to our model which benefits from this frequency-domain property).

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

The Marginalisation 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, introduced in Section 2.3.2. 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 optimisation process, and we show in Section 4.3.3 that doing so improves the deblurred results.

Recall from Section 2.3.2 that the aim is to find the factorised distribution

q(Θ) = q(β σ ) j q(f ∂ j ) k q(w k ) (4.1)
that best approximates the true posterior p(Θ|g). In Section 2.3.2, the "ensemble" Θ incorporated the latent image gradients f ∂ , a convolution kernel a and the noise precision

β σ .
In our case we replace the convolution kernel with the blur kernel w for our model.

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 Appendix A). 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 T (k) ij f ∂ j 2 q(f ∂ ) (4.2) w (1) k w (2) k = β σ i g i j T (k) ij f ∂ j q(f ∂ j ) - k ′ =k j T (k) ij f ∂ j j T (k ′ ) ij f ∂ j q(f ∂ ) w k ′ q(w k ′ ) (4.3) f (2) j = β σ i k T (k) ij w k 2 q(w) (4.4) f (1) j f (2) j = β σ i g i k T (k) ij w k q(w k ) - j ′ =j f ∂ j ′ q(f ∂ j ′ ) k T (k) ij ′ w k k T (k) ij w k q(w) , (4.5)
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. For spatially-invariant blur (i.e. when the matrices T (k) apply 2D translations of the image), these equations reduce to those given by [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF]. Note that these equations cannot be implemented directly in this form, and must be re-written as computations involving only the parameters f

(1) j , f (2) j , w (1)
k , and w

(2) k , as follows. These are exactly the parameter updates computed in our implementation of the marginalisation algorithm described in Chapter 2. First, let

v w k = w 2 k -w k 2 (4.6) v f j = f ∂ j 2 -f ∂ j 2 (4.7) A ij = k T (k) ij w k (4.8) B ik = j T (k) ij f ∂ j (4.9) g * i = k j T (k) ij f ∂ j w k . (4.10) Then, w (2)
k = β σ i,j T (k) ij 2 v f j + β σ i B ik 2 (4.11) w (1) k w (2) k = β σ i B ik (g i -g * i ) -β σ i,j T (k) ij A ij v f j + w k w (2) k (4.12) f (2) j = β σ i,k T (k) ij 2 v w k + β σ i A ij 2 (4.13) f (1) j f (2) j = β σ i A ij (g i -g * i ) -β σ i,k T (k) ij B ik v w k + f ∂ j f (2) j . (4.14)
Finally, in order to compute the expectation β σ of the noise variance, it is necessary to evaluate the expected value of the squared reconstruction error (Miskin and MacKay, 2000, Eqn. (40)), which can be computed as follows:

(g i -g * i ) 2 = (g i -g * i ) 2 + j,k T (k) ij 2 v f j v w k + j A ij 2 v f j + k B ik 2 v w k . ( 4

The Maximum a Posteriori Approach

In this section we adapt the single image PSF estimation algorithm of [START_REF] Cho | Fast motion deblurring[END_REF], introduced in Section 2.3.3, to use our spatially-variant blur model. Like the marginalisation algorithm of the previous section, this algorithm can be readily adapted to handle non-uniform blur, substituting our model in place of convolution.

The main changes to the algorithm described in Chapter 2 are simply to replace the convolution in Equation (2.53) with our blur model, so that instead of estimating a convolution kernel â by solving Equation (2.53), we instead estimate the spatially-variant blur descriptor ŵ by solving min

w q ω q k w k T (k) p (q) -d (q) * g 2 2 + β w 2 2 , (4.16)
where ω q are weights for the data fidelity term, p (q) are thresholded derivative maps obtained by non-linear filtering of the current estimate of the sharp image, and d (q) are derivative filters. Similarly, when updating the latent image, we modify Equation (2.56) to be min

f q ω q k ŵk T (k) (d (q) * f ) -d (q) * g 2 2 + α d x * f 2 2 + α d y * f 2 2 .
(4.17)

Since the spatially-variant blur model is linear in w and f , these are still linear least squares problems. Although we are not able to take full advantage of the speed optimisations proposed by [START_REF] Cho | Fast motion deblurring[END_REF], 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 marginalisation algorithm of Section 4.3.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. This can be seen by considering a vertical blur at the left or right-hand 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 in order to best constrain the kernel, we must ensure that the pixels used to estimate the kernel (the non-zeroes 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 well-conditioned. 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 in-plane rotation of the camera. As can be seen in Figure 3.4, a rotation of the camera about its Z axis causes a very small displacement for pixels towards the centre 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 minimising Equation (4.16) with our model generally contain many non-zeros spread smoothly throughout, and do not produce good deblurred outputs (see Figure 4.5).

If instead of the ℓ 2 regularisation in Equation (4.16), we apply ℓ 1 regularisation combined with non-negativity constraints, the optimisation 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 is equivalent to an exponential prior on the kernel values. In this case, the kernel estimation problem becomes min

w q ω q k w k T (k) p (q) -d (q) * g 2 2 + β k w k s.t. ∀k = 1, . . . , K, w k ≥ 0. (4.18)
This is an instance of the lasso problem [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], for which efficient optimisation algorithms exist [START_REF] Efron | Least angle regression[END_REF][START_REF] Kim | An interior-point method for large-scale ℓ 1 -regularized least squares[END_REF][START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF]. This problem is convex, so that we can be sure of obtaining a global minimum. The different results obtained using ℓ 2 and ℓ 1 regularisation are discussed in Section 4.3.3. With the use of the ℓ 1 regularisation, 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 work using ℓ 1 regularisation, we set β = 0.01. In the remainder of the thesis, we refer to the original algorithm of Cho and Lee as MAP-ℓ 2 , and our ℓ 1 regularised version as MAP-ℓ 1 .

Single-Image Deblurring Results

We show in this section results of single-image deblurring using the algorithms described in Section 4.3 to estimate the spatially-variant PSF, 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. The estimated PSF is used with the non-blind deblurring algorithm of [START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF] to estimate the final deblurred image. This algorithm is easily adapted to non-uniform blur since it involves repeated minimisations of quadratic cost functions of the form (see Equation (2.46))

min f Af -g 2 2 + α d x * f -v x 2 2 + α d y * f -v y 2 2 , (4.19)
where v x and v y are intermediate variables of the optimisation scheme. For our nonuniform blur model, we use the conjugate gradient algorithm to minimise this cost function.

Figure 3.1 in Chapter 3 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 any least-squares based algorithms, such as that of [START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF]. We examine non-blind deblurring for images containing saturated pixels in Chapter 6. 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 4.2 (d), we also demonstrate the use of the variational marginalisation algorithm of Fergus et al. to produce the deblurred output, as opposed to the algorithm of [START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF], used in Figure 4.2 (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. The result of blind deblurring on a real camera shake image (a), captured with a shutter speed of 1 2 second, using both the marginalisation 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 marginalisation 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. 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 summarised 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. 4.2. Shown is the estimated kernel and deblurred result when using our non-uniform blur model in the algorithm of Cho and Lee with ℓ 2 regularisation on the kernel. As can be seen, the ℓ 2 regularisation is not sufficient to produce a good estimate of the kernel, and results in a deblurred output containing many artefacts. 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 4.4 shows results for blind deblurring of synthetically blurred images using the two methods (marginalisation and MAP), 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 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 marginalisation 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 marginalisation approach. In Figure 4.6, we compare our approach to that of [START_REF] Fergus | Removing camera shake from a single photograph[END_REF] 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 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 4.7 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 ignoring the camera's translation is reasonable in practice.

Note also that we estimate the kernel directly from the blurry image, whereas Joshi et al.

(2010) use additional hardware to record the camera's motion.

Discussion

Besides comparing the results of a given algorithm with either a uniform or non-uniform blur model, we can also compare the marginalisation 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 daytime image of a street scene, so the application of this prior to an image with a very different distribution of intensities (e.g. a photo taken at night) 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 3.1, the marginalisation approach (using the street scene prior) fails to find a useful kernel, while the MAP approach finds a good kernel, as demonstrated by the deblurred result in Figure 3.1 (b).

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 marginalisation approach has almost no parameters to tune. The problem of setting these parameters automatically merits future study, as it often requires some amount of manual tuning before a good result is achieved. A principled way of setting α and β would link their values to the noise variance in the image, however we defer this to future work. In this thesis, α and β were set manually, although we have found that α = 0.0005 and β = 0.01 to be generally good values. 2010), with the deblurred results from the original work (b) and using the MAP-ℓ 1 method with our blur model (c). The blur in this image is approximately 35 pixels long, and 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. To estimate the kernel for this image, we found that was necessary to reduce the threshold on the kernel values (Section 2.3.3) to 1 100 the maximum value, instead of 1 20 , in order to capture the long thin structure of the non-zeros. Without this, the vertical line of non-zeros was set to zero early in the MAP-ℓ 1 algorithm and was not subsequently recovered (see Section 4.5.1 for more details).

Convergence. Since both algorithms (marginalisation and MAP) attempt to solve nonconvex minimisation problems, neither can guarantee the ability to arrive at a globally optimal solution. However, in practice we have found them both to perform reliably. This is due in large part to the multi-scale approach; 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 minimisations 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 edge prediction step helps direct the optimisation process towards a desirable minimum, by encouraging the sharp image to contain step edges. In the marginalisation algorithm, we have found that the algorithm converges equally reliably for both the uniform model and our model, although a greater number of iterations are typically required when using our model, due to the greater number of parameters to estimate.

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 marginalisation algorithm, in Equations (4.11) to (4.15), are computationally expensive. Due to this expense, and the larger number of iterations required for our model compared to the uniform model, the marginalisation algorithm with our non-uniform model can take several hours to deblur an image of several hundred pixels across under Matlab on a Linux workstation with an Intel Xeon 2.93GHz CPU and 8GB of memory. Deblurring larger images with this method is not currently practical.

Our Matlab and C implementation of the MAP algorithm, on the other hand, can deblur the same images in under an hour, and megapixel images in 1 -3 hours, depending on the size of the blur (see Chapter 5 for an approximation which allows us to substantially decrease the blind deblurring time further).

Limitations. Both of these algorithms are capable of removing large blurs -the blurs removed in this thesis are up to 35 pixels wide (e.g. Figures 1.1 and 6.1), 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 1000 pixels.

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 for images which contain only fine-scale texture, or where the blur is too large to allow this.

Deblurring With Noisy / Blurry Image Pairs

As for the single-image deblurring algorithms discussed in the previous sections, we directly substitute our blur model in place of convolution in the original kernel estimation problem of Yuan et al. (2007b) in Equation (2.58). As discussed in Figure 2.3, some prior knowledge must be applied to recover a good kernel estimate. In their algorithm, Yuan et al. (2007b) constrain the kernel to have non-negative elements and unit ℓ 1 norm, however they simultaneously penalise the ℓ 2 norm of the kernel, reducing the sparsity-inducing effect of the constraint. To help find a sparse kernel, they propose a thresholding scheme with hysteresis 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 solve the following problem min

w k w k T (k) f N -g 2 2 s.t. ∀ k = 1, . . . , K, w k ≥ 0 and k w k = 1. (4.20)
Similar to Equation (4.18), this least-squares formulation with non-negative ℓ 1 constraints can be solved efficiently [START_REF] Kim | An interior-point 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.

such that we can be sure that we have reached the globally-optimal deblurred image. The solutions are also not particularly sensitive to the values of regularisation parameters, since the problems are well-conditioned, with many more observations than unknowns.

However, this approach is only applicable when a noisy image f N is available, which precludes most existing blurry images. As a method of capturing a good image with a standard camera in low-light however, the noisy / blurry image pair is a good technique for photographers to keep in mind.

Implementation

The implementation of the variational kernel estimation method presented in Section 4.3.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] 1 . 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 Equations (4.11) to (4.15). A package containing our code is available online2 . The implementation of the non-blind deblurring algorithm of [START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF] is also based on Matlab code made available online by the authors3 . 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. (2007b) are our own, and we use these implementations for both uniform and non-uniform 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.

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 re-fining 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 (4.20) directly at full resolution would involve transforming f N for every possible rotation under consideration and storing all the copies simultaneously in a matrix B N .

For a blur of size 10 pixels in an image 1000 pixels wide, this could be over 1000 copies of f 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 chapter, which estimate the kernel iteratively, we use our current estimate of the kernel ŵs at a scale s to constrain our estimate at the next iteration. To do this, we define an "active region" where ŵs is non-zero, and constrain the non-zeros at the next iteration to lie within this region. By fixing many kernel elements to zero, we eliminate a large amount of computation and memory requirements associated with estimating those elements' values. In the example of the previous paragraph, this corresponds to discarding many columns of B N , clearly reducing the size of the problem. 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 iteration, we find the non-zero elements of our current estimate of the kernel ŵs , and dilate this region using a 3 × 3 × 3 cube to define the active region for the next iteration. When moving from one scale s to the next scale s + 1, we upsample ŵs using bilinear interpolation, find the non-zero elements of this upsampled kernel, and as before, dilate this region using a 3 × 3 × 3 cube. This initialises the active region for our next estimate ŵs+1 . We repeat this process at each scale, until we have found the optimal kernel at the finest scale.

This approach is generally effective at reducing the computational burden of the kernel estimation without reducing accuracy, however, problems may occur if the blur kernel contains long faint structures, as it is possible for these to be clamped to zero at a coarse scale and never to be recovered, as mentioned in Figure 4.7.

Conclusion

We have applied our new model for spatially-variant camera shake blur within the frameworks of several 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. In Section 5.2 we examine some of the bottlenecks in the spatially-variant blind deblurring process, before introducing an efficient approximation in Section 5.3 to the model proposed in Chapter 3. We compare results and deblurring times using the approximation, to those using the exact model, and conclude in Section 5.4.

Bottlenecks in Spatially-Variant Blind Deblurring

Although the model for a spatially-variant PSF introduced in Chapter 3 is parameterised by the vector w, which may have a similar number of elements to a traditional convolution kernel, the computational cost of using this model in blind deblurring is substantially higher. This is due to several factors, for example that the transformations {T (k) } involve sub-pixel interpolations, that the set of matrices {T (k) } is too large to fit in memory and must be re-computed on-the-fly frequently, and that there is no equivalent to the frequency-domain convolution which allows spatially-invariant blur (2D convolution) to be computed in O N log N time instead of O N 2 , where N is the number of pixels in the image. These factors mean that blind PSF estimation using this model is significantly slower than for spatially-invariant blur. In this section we discuss the cost of the most computationally-intensive steps in the MAP blind PSF estimation algorithm.

Recall from Chapter 4 that the MAP algorithm has two major steps: updating the blur kernel and updating the sharp image; we discuss each of these in turn.

Updating the kernel

To update the kernel in the MAP-ℓ 1 blind deblurring algorithm, recall from Chapter 4 that we solve Equation (4.18) (p. 75). For ease of exposition, we note that the problem can be written in the general form

min w Pw -g ∂ 2 2 + β k w k s.t. ∀k = 1, . . . , K, w k ≥ 0, (5.1)
where K is the number of non-zero weights in the kernel, the vector g ∂ is the concatenation of all the derivatives of g, and the matrix P incorporates all the transformed versions of the predicted sparse derivative maps {p (q) } (as discussed in Section 2. 3.3 (p. 38)). The matrix P has the form

P =      T (1) p (1) T (2) p (1) • • • T (K-1) p (1) T (K) p (1) T (1) p (2) T (2) p (2) • • • T (K-1) p (2) T (K) p (2) . . . . . . . . . . . . . . .      , (5.2)
where in the algorithm of [START_REF] Cho | Fast motion deblurring[END_REF], p (1) = p x , p (2) = p y , and p (q) for 3 ≤ q ≤ 5 contain second-order derivatives.

In iterative minimisation algorithms used for solving such least-squares problems, e.g. LARS-lasso [START_REF] Efron | Least angle regression[END_REF] or conjugate-gradient descent (Shewchuk, 1994),

we must repeatedly compute matrix-vectors products such as Px, P ⊤ Px and P ⊤ y for arbitrary vectors x and y. The 5N × K matrix P is generally too large to fit into memory (e.g. for a megapixel image P could have size 5 • 10 6 × 10 3 elements, requiring 37 GB of memory at double precision), so we must generate its elements on-the-fly when computing these matrix-vector products. Computing each element of P is essentially equivalent to performing a sub-pixel interpolation, which is computationally expensive.

As a result, these matrix-vector products are by far the biggest bottleneck in the kernel estimation step, since they might involve computing a large number of image warps and multiplications at each iteration.

To solve Equation (5.1), we use the LARS-lasso algorithm of [START_REF] Efron | Least angle regression[END_REF], since it involves relatively few products of the type described above per iteration, and is guaranteed to reach the global minimum in a finite number of steps. To perform a single iteration of this algorithm, we have the choice of computing three Px-type products and one P ⊤ y-type product, or computing one P ⊤ Px-type product. Here, note that since K << N , it is generally feasible to store the K × K matrix P ⊤ P, referred to as the Gram matrix, in memory. This matrix can be computed without requiring any element of P more than once, using the fact that P ⊤ P = i p i p ⊤ i , where p ⊤ i is the i th row of P. Thus we can perform the entire LARS-lasso kernel estimation process while computing each element of P only once. With P ⊤ P stored in memory, the P ⊤ Px product at each iteration can be performed with optimised linear algebra libraries. We have found empirically that this approach (pre-computing P ⊤ P, as opposed to computing the Px and P ⊤ y products on-the-fly) gives the shortest kernel estimation times. However, the pre-computation of P ⊤ P may still by time-consuming. For megapixel images, or images with many non-zeros in {p (q) }, the pre-computation may take minutes, or even hours.

Updating the sharp image

To update the sharp image in the MAP-ℓ 1 blind deblurring algorithm, we use conjugate gradient descent (CG) to solve the minimisation problem in Equation (4.17) (p. 73), which essentially has the form (5.3) where A = k w k T k . To compute the gradient of the cost function, we need to evaluate products of the type A ⊤ Af at each iteration. As in the kernel estimation step, the image transformations T (k) f involved in computing these products are the most computationally-expensive operations. Unlike for the kernel estimation step, the N × N Gram matrix A ⊤ A is too large to store, so cannot be pre-computed. However, we typically perform a small, fixed number of iterations of CG (e.g. 20 iterations), and as such the time spent performing this optimisation is typically much less than the kernel update step. Nevertheless, this CG minimisation may still take several minutes for a megapixel image, as shown in Figure 5.3 (b).

min f Af -g 2 2 + α d x * f 2 2 + α d y * f 2 2 ,

Locally-Uniform Approximation

In this section we describe an approximation to the forward model described in Chapter 3 that provides a fast means of calculating the spatially-variant forward model. More importantly however, we describe how this approximation can be used to (a) compute the kernel Gram matrix P ⊤ P very quickly, allowing the kernel update step to be performed in seconds rather than hours, and (b) approximately minimise Equation (5.3) directly in two steps, allowing the sharp image to be updated in seconds rather than minutes.

This approximation is based on the "Efficient Filter Flow" model for spatially-variant blur proposed by [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF], discussed earlier in Section 2.4.2. [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF] observe that in some cases of spatially-variant image blur, the blur may vary slowly and smoothly across the image. As a result, the PSFs of nearby pixels may be very similar, and it is reasonable to approximate spatially-variant blur as being locally-uniform. Following this observation, they propose an approximation of the forward model, whereby the sharp image f is covered with a coarse grid of P overlapping patches, each of which is modelled as having a spatially-invariant blur. The overlap between patches ensures that the blur varies smoothly across the image, rather than changing abruptly at the boundary between two patches. The fact that each patch has a spatially-invariant blur allows the forward model to be computed using P small convolutions. [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF] assign each patch r a spatially-invariant blur filter a (r) , and the forward model is given by:

g * = P r=1 C (r) ⊤ a (r) * (m • C (r) f ) , (5.4)
where C (r) is a matrix that crops out the r th patch from the image f (and thus C (r) ⊤ reinserts it at its correct location). The vector m is a windowing function, e.g. the Bartlett-Hann window, which produces the smooth transition between neighbouring patches, and • • • represents the Hadamard (element-wise) product between two vectors.

This model can be computed efficiently in the frequency domain:

g * = P r=1 C (r) ⊤ F -1 F a (r) • F m • C (r) f , (5.5)
where F ( • ) takes the 2D discrete Fourier transform (computed using the fast Fourier transform), and F -1 ( • ) the inverse Fourier transform. The windowing function m is chosen such that the total weight at each pixel is 1, i.e. r C (r) ⊤ m = 1.

A Globally-Consistent Approximation for Camera Shake

In their original work, [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF] estimate a separate filter a (r) for each patch r. Likewise Harmeling et al. (2010b), who apply this approximation in the context of single-image camera shake removal, also estimate a separate filter per patch, using the MAP algorithm of [START_REF] Cho | Fast motion deblurring[END_REF]. One weakness of this approach is that in regions with little texture or strong edges, it may not be possible to estimate a good kernel.

Since the purely-local model has no notion of global consistency between the different regions' kernels, heuristics are needed to encourage neighbouring filters to be similar, and to detect failed regions and hypothesise kernels for those regions, based on their neighbours. However, given the forward blur model for camera shake in Equation (3.11), which is parameterised by a single set of weights w, we can in fact write each a (r) in terms of w. For each patch r, we choose a (r) to be the point spread function for the central pixel i r , which is given by the i th r row of A. Since A is linear in w, we can construct a matrix J (r) such that a (r) = C (r) J (r) w. The elements of each matrix J (r) are simply a re-arrangement of the elements of the matrices T (k) : J 5.1 shows how the quality of the approximation varies with the number of patches being used, compared to the exact model. In all our experiments, we use a grid of 6×8 patches.

(r) jk = T (k) irj . Figure
Having written each filter a (r) in terms of w, we can then substitute this into Equation (5.5) to obtain the following approximation of the forward model in Equation (3.11):

g * ≃ P r=1 C (r) ⊤ F -1 F C (r) J (r) w • F m • C (r) f .
(5.6)

These equations allow the forward model to be computed quickly using only a handful of frequency-domain convolutions. Furthermore, the derivatives of g * with respect to f and w can also be computed using a small number of frequency-domain convolutions and correlations. As discussed in the previous section, computing the forward model and its derivatives are bottlenecks in both the blind PSF estimation algorithm of [START_REF] Cho | Fast motion deblurring[END_REF].

Equation (5.6) provides us with a fast approximate way of computing the forward model for g * = Af = Bw. By re-arranging this approximation, we can compute B ⊤ y for arbitrary y and B ⊤ B using 3 FFTs per patch as follows:

P ⊤ y ≃ P r=1 J (r) ⊤ C (r) ⊤ F -1 F m • C (r) f * • F C (r) y
(5.7) r) .

P ⊤ P ≃ P r=1 J (r) ⊤ C (r) ⊤ XCorrMatrix F -1 F m • C (r) f * • F C (r) f C (r) J (
(5.8)

The function XCorrMatrix constructs the full cross-correlation matrix from a crosscorrelation vector between two signals. To clarify, consider a pair of 1D signals u and v.

We can compute the cross-correlation vector c for all relative offsets between them as c =

F -1 F (u) * • F (v) .
If we define the matrix S k which shifts a signal by k elements, then the element c k stores the cross-correlation between S k u and v, i.e. c k = (S k u) ⊤ v. With appropriate boundary conditions, we also have

c k = (S k-1 u) ⊤ (S 1 v) = (S k-2 u) ⊤ (S 2 v),
and so on. Thus to compute the cross-correlation matrix C with elements

C k 1 k 2 = (S k 1 u) ⊤ (S k 2 v)
, we need only re-arrange and duplicate the elements of the vector c. This argument extends naturally to 2D, and this action is denoted by C = XCorrMatrix(c).

Using the approximation described above, computation of

P ⊤ P costs O K+N log(N ) .
This is in comparison to the cost of O N K 2 to compute P ⊤ P using the exact model.

For camera shake, the blur kernel tends to consist of a line of non-zeros, so that the number of kernel elements K grows (at best) linearly with the length of the blur. Since the length of the blur will grow linearly with the image width W , we can assert that K grows as than √ N . Thus the cost of computing P ⊤ P with the exact model grows at least as O N 2 , while the approximate model grows as O N log N .

Fast Independent Non-Blind Deconvolution of Patches

Although the use of a locally-uniform approximation allows the forward model to be computed much faster using the FFT, we cannot immediately take full advantage of working in the frequency domain when it comes to non-blind deconvolution. For example, in the inner loop of the blind deblurring algorithm of [START_REF] Cho | Fast motion deblurring[END_REF], we must repeatedly solve problems of the form in Equation (5.3). When the blur is spatiallyinvariant, Af = a * f , and Equation (5.3) can be minimised directly in a single-step by pixel-wise division in the frequency domain, as discussed in Section 2.2.2. However, Blind deblurring time using the exact and approximate model. This log-log plot shows the time taken to perform blind deblurring of several images at various different sizes, using either the exact forward model (solid lines) or the approximate forward model (dashed lines). The blind deblurring algorithm is implemented in Matlab, and the function to compute the exact forward model is implemented in compiled C, while the approximate model is implemented in Matlab. The blind deblurring time using the exact model grows significantly faster with respect to the image size than using the approximate model, and for images 1024 pixels wide already takes almost 2 orders of magnitude longer. The images for these times are from the following figures: Istanbul -Figure 6 for spatially-variant blur, the locally-uniform approximation does not permit this. Even though each patch has a spatially-invariant blur, the fact that the patches overlap means that they cannot be deconvolved independently. This can be seen by expanding the nonblind deblurring problem in Equation (5.3), where the sum over patches lies inside the data fidelity term of the cost function:

181
min f P r=1 C (r) ⊤ a (r) * (m • C (r) f ) -g 2 2 + α d x * f 2 2 + α d y * f 2 2 .
(5.9)

If we note however that extracting and recombining the patches of an image has no 101 effect on an image, or on a filtering operation (Equation (2.64)), i.e.

g = P r=1 C (r) ⊤ (m • C (r) g)
(5.10)

d x * f = P r=1 C (r) ⊤ (d x * m • C (r) f ) , (5.11)
we can rewrite Equation (5.9) as

min f P r=1 C (r) ⊤ a (r) * (m • C (r) f ) -m • C (r) g 2 2 + α P r=1 C (r) ⊤ d x * (m • C (r) f ) 2 2 + α P r=1 C (r) ⊤ d y * (m • C (r) f ) 2 2 .
(5.12)

Although this cost function still does not permit a closed-form solution, we can obtain the following upper bound (up to a constant scale factor of P ) using Jensen's inequality [START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF]:

min f P r=1 C (r) ⊤ a (r) * (m • C (r) f ) -m • C (r) g 2 2 + α P r=1 C (r) ⊤ d x * (m • C (r) f ) 2 2 + α P r=1 C (r) ⊤ d y * (m • C (r) f ) 2 2 . (5.13)
This formulation penalises each patch independently, and can be minimised using a splitvariable approach, iterating between (a) estimating a deblurred patch f (r) for each patch r, and (b) estimating the complete deblurred image f that best matches these patches.

In practice a single iteration is sufficient to obtain a good result. As we show in the following, step (a) can be performed independently for each patch, and step (b) can be performed independently for each pixel of the deblurred image.

Specifically, we define the set of blurry patches g (r) such that (5.14) and note that the multiplication by C (r) ⊤ at the beginning of each term in Equation (5.13) only zero-pads an image, so has no effect on the cost. We first estimate a set of deblurred patches { f (r) } by solving min

g (r) = m • C (r) g,
f (r) a (r) * f (r) -g (r) 2 2 + α d x * f (r) 2 2 + α d y * f (r) 2 2 (5.15)
for each patch. This can be done using the direct method of Equation (2.44). Secondly, to estimate the full deblurred image that best matches the deblurred patches f (r) , we solve

min f P r=1 m • C (r) f -f (r) 2 2 .
(5.16)

To minimise this, we first note that we can multiply through by C (r) ⊤ without altering the cost (since C (r) ⊤ simply applies zero-padding), and that

C (r) ⊤ (m • C (r) f ) = (C (r) ⊤ m) • f .
Thus we obtain the following equivalent problem:

min f P r=1 (C (r) ⊤ m) • f -C (r) ⊤ f (r) 2 2 .
(5.17)

This problem can be solved independently for each pixel of f , yielding the following solution:

f = r C (r) ⊤ (m • f (r) ) r C (r) ⊤ (m • m)
.

(5.18)

Algorithm 1 summarizes the steps to performing this fast direct deconvolution, and in Figure 5.3 we compare this method of non-blind deblurring to other possibilities for solving Equation (5.3). The fast independent method produces results visually very similar to the exact model, in a significantly shorter amount of time. We use this direct deconvolution method to update the latent image in blind deblurring, as discussed in Section 5.2.2. This provides an additional speed improvement in blind deblurring, compared to the use of conjugate-gradient descent with the approximate forward model. CG with the approximate forward model, and (d) direct deconvolution using the approach described in Section 5.3.2. The results are visually similar using all three methods. Using CG with the approximate forward model is much faster than with the exact model, however the direct approach takes only a fraction of the time of either of these.

Algorithm 1: Fast non-blind deconvolution of spatially-varying blur.

Input: Blurry image g, blur descriptor w Output: Deblurred image f

1 Decompose g into patches g (r) using (5.14)

2 Deconvolve each patch to estimate f (r) by solving (5.15) in the frequency domain

3 Combine patches to estimate sharp image f using (5.18)

Conclusion

We have described how an efficient approximation for spatially-variant blur can be used to reduce the computational cost of computing the spatially-variant forward blur model from Chapter 3. In addition, we have proposed an approximate direct method for spatially-varying non-blind deblurring, which is able to perform least-squares non-blind deblurring in a fraction of the time required for iterative methods such as conjugate gradient descent. As a result of the approximations described, we are able to obtain a substantial speed-up in both the blind PSF estimation and the non-blind deblurring steps over the exact model, reducing the time required to completely deblur an image from several hours to a few minutes, with no considerable reduction in quality.

For a 1024 × 768 image, our C implementation of the exact model in Equation (3.11) takes approximately 5 seconds to compute, compared to 2 seconds for our Matlab implementation of the approximate forward model in Equation (5.6), on an Intel Xeon 2.93GHz CPU.

Chapter 6

Handling Saturation in Non-Blind Deblurring

Introduction

When the PSF for a blurry image is known, it is the goal of "non-blind" deblurring to recover an estimate of the sharp image. This problem has a long history, having been studied since the 1940s [START_REF] Wiener | Extrapolation, Interpolation, and Smoothing of Stationary Time Series[END_REF]. However, one feature of blurry images, and in particular "shaken" images, that is problematic for non-blind deblurring but that has received very little attention is the presence of saturated pixels. These are caused when the radiance of the scene exceeds the range of the camera's sensor, leaving bright highlights clipped at the maximum output value (e.g. 255 for an 8-bit image). To anyone who has attempted to take hand-held photographs at night, this effect should be familiar as the conspicuous bright streaks left by electric lights, such as in Figure 6.1 (a). These bright pixels, with their clipped values, violate the assumption made by many algorithms that the image formation process is linear, and as a result can cause obtrusive artefacts in the deblurred images. This can be seen in the deblurred images in Figures 6.1 (b) and 6.1 (c).

Saturation has not received wide attention in the literature, although it has been cited as the cause of artefacts in the deblurred outputs from deconvolution algorithms.

For example, [START_REF] Fergus | Removing camera shake from a single photograph[END_REF], [START_REF] Cho | Fast motion deblurring[END_REF] and [START_REF] Tai | Richardson-Lucy deblurring for scenes under a projective motion path[END_REF] mention the fact that saturated pixels cause problems, sometimes showing their effect on the deblurred output, but leave the problem to be addressed in future work. An exception is Harmeling et al. (2010a), who address the issue in the setting of multi-frame blind deblurring by thresholding the blurry image to detect saturated pixels, and ignoring these in the deblurring process. When many blurry images are available, these pixels can be safely discarded, since there will generally remain unsaturated pixels covering the same area in other images.

In this chapter we address the problem of non-blind deblurring for images that contain saturated pixels. By handling such pixels explicitly, we are able to produce significantly better results than existing methods. Figure 6.1 (d) shows the output of the proposed algorithm, which contains far fewer artefacts than the two existing algorithms shown for comparison. We begin in Section 6.2 by considering possible ways of explicitly handling saturated pixels, and propose modified versions of the Richardson-Lucy algorithm for each. In Section 6.3 we augment the Richardson-Lucy algorithm with an additional modification, which prevents the propagation of ringing artefacts. We discuss the implementation of the algorithm in Section 6.4, and show results in Section 6.5.

In Section 6.6 we look at the different causes of ringing in non-blind deblurring. We conclude in Section 6.7.

Explicitly Handling Saturated Pixels

We model sensor saturation as follows: the sensor outputs pixel values which are proportional to the scene radiance, up to some limit, beyond which the pixel value is clipped at the maximum output level. This model is supported by the data in Figure 6. of the exposure increases, more pixels saturate.

This suggests two possible ways of handling saturation when performing non-blind deblurring: (a) discard the clipped pixels, so that we only use data which follows the linear model, or (b) modify the forward model to take into account this non-linear relationship. We describe both of these approaches in the following, by modifying the popular Richardson-Lucy algorithm, and demonstrate that the second approach is better, as it does not discard any potentially-useful data.

Discarding Saturated Pixels

It is possible to estimate which blurry pixels are saturated by defining a threshold T , above which a blurry pixel is considered to be saturated, and therefore an outlier to the linear model. If we discard these pixels, the problem of deblurring with saturated pixels becomes deblurring with missing data.

Richardson-Lucy With Missing Data

We begin by defining the set of unsaturated (inlier) blurry pixels Z = {i|g i ≤ T }. We also define the corresponding binary mask z, where each element z i = 1 if blurry pixel g i ≤ T , and z i = 0 if g i > T . The image formation model for a single blurry pixel g i in Z is

g * i = j A ij f j (6.1) ∂g * i ∂f j = A ij . (6.2)
Recall from Section 2.1.4 that for a Poisson noise model, the likelihood for the inlier blurry pixels is

p(g Z |w, f ) = i∈Z g * i g i e -g * i g i ! (6.3)
The negative log-likelihood is then:

L Z = - i∈Z log g * i g i e -g * i g i ! (6.4) = i∈Z g * i -g i log g * i + log(g i !), (6.5)
The Richardson-Lucy algorithm minimises L Z over the sharp image f . The Karush-Kuhn-Tucker conditions for the optimum of this likelihood are, for all f j :

f j ∂ ∂f j L Z = 0, (6.6) and ∂ ∂f j L Z ≥ 0 if f j = 0. (6.7)
The Richardson-Lucy update rule is a fixed-point iteration derived from Equation (6.6):

f j i∈Z ∂ ∂f j (g i log g * i -g * i -log(g i !)) = 0 (6.8) f j i∈Z g i g * i ∂g * i ∂f j - ∂g * i ∂f j = 0. (6.9)
Substituting in Equation ( 6.2) we get:

f j i∈Z A ij = f j i∈Z A ij g i g * i . (6.10)
We can replace the sum over Z with a weighted sum over all blurry pixels, using the mask elements z i as the weights:

f j i A ij z i = f j i A ij g i g * i z i . (6.11)
Here, in order to obtain an expression for f j , we could divide both sides by i A ij z i , however this would be unstable in regions where z i ≃ 0. Thus we add f j (1

-i A ij z i )
to both sides:

f j = f j i A ij g i g * i z i + 1 - i A ij z i . (6.12)
Assuming that the blur is conservative, i.e. i A ij = 1, we can factorise this and obtain the update rule by replacing the unknowns f j and g * i with our estimates fj and (A f ) i = j ′ A ij ′ fj ′ , and adding iteration indices:

f t+1 j = f t j i A ij g i (A f t ) i z i + 1 -z i , (6.13)
In matrix-vector form, this update equation is

f t+1 = f t • A ⊤ g • z A f t + 1 -z , (6.14)
where 1 is a vector of ones. For an unsaturated pixel g i , the mask z i = 1, and the term in parentheses is the same as for the standard RL update. For a saturated (outlier/missing) pixel, z i = 0, so the term in parentheses is equal to unity. Since the update is multiplicative, this means that the saturated observation g i has no influence on the estimated latent image f .

The choice of threshold T can be problematic however; there may not be a single threshold that correctly separates the saturated pixels from the unsaturated ones. A low threshold may discard large numbers of inlying pixels from g, causing some parts of f to become decoupled from the data, and remain unchanged from their initialisation by the algorithm. A high threshold, on the other hand, may treat some saturated pixels as inliers, causing artefacts in the deblurred result. Figure 6.3 shows the result of deblurring using Equation (6.14) for different values of threshold T . As is visible in the figure, no particular threshold produces a result free of artefacts. At high values of T , the building is deblurred well, but artefacts appear around the lights. At the lowest value of T , the lights are deblurred reasonably well, but the face of the building is mistakenly discarded and thus remains blurry in the output. As the threshold decreases, the artefacts around the bright lights are reduced compared to the standard RL result in Figure 6.1 (b). At the lowest threshold (d) the fewest artefacts appear, but parts of the church are also discarded, hence remain blurred.

A Forward Model for Saturation

Instead of attempting to segment the blurry image into saturated and unsaturated regions, we may instead modify our forward model to include the saturation process. This avoids making a priori decisions about which data are inliers or outliers. To this end, we modify the original image degradation model in Equation (2.1) by introducing a response function R : R + → R + , which models the non-linear response of the sensor. This function is applied to the degraded image before it is output by the sensor, i.e. (6.15) where, recalling the notation of Section 2.1, H is the blurring operator and N represents the perturbation by random noise. This formulation however, is not particularly tractable, except for trivial response functions R. Since the noise lies inside R, the likelihood is distorted and can generally no longer be written in closed-form.

g(x) = R N (Hf )(x) ,
Instead, we use the following, more tractable, formulation: (6.16) which allows us to incorporate the response into the formation of the "noiseless" blurry image g * (x) = R (Hf )(x) , (6.17) (6.18) where the function R is applied element-wise to a vector.

g(x) = N R (Hf )(x) ,
g * = R (Af ) ,

Richardson-Lucy With a Non-Linear Response

We can incorporate this non-linear image formation model into the Richardson-Lucy algorithm. From Equation (6.18) we have (6.19) Substituting this into Equation (6.9) and following a similar derivation as for the Richardson-Lucy with missing data, we arrive at the update rule

∂g * i ∂f j = A ij R ′ (Af ) i .
f t+1 j = f t j i A ij g i R (A f t ) i R ′ (A f t ) i + 1 -R ′ (A f t ) i .
(6.20)

In matrix-vector form, this update equation is

f t+1 = f t • A ⊤ g • R ′ (A f t ) R(A f t ) + 1 -R ′ (A f t ) . (6.21)

Choosing the Response Function

One choice for R would be simply to truncate the linear model at 1 (the maximum pixel value), i.e. R(x) = min(x, 1). This choice is empirically justified, as can be seen in Figure 6.2. However, this function is non-differentiable at x = 1, i.e. R ′ (1) is not defined.

We thus use a smooth approximation [START_REF] Chen | A class of smoothing functions for nonlinear and mixed complementarity problems[END_REF], where

R(x) = x - 1 a log 1 + exp a(x -1) (6.22) R ′ (x) = 1 1 + exp a(x -1) , (6.23)
where R ′ represents the first derivative of R. The parameter a controls the smoothness of the approximation, and we have found a = 50 to be a suitable compromise between smoothness and accuracy. Figure 6.5 shows the shape of R and R ′ compared to the simple truncated linear model. In all our experiments a = 50.

Given the shape of R, Equation (6.21) can easily be interpreted: in the linear portion R(x) ≃ x and R ′ (x) ≃ 1, so that the term in parentheses is the same as for the standard RL algorithm, while in the saturated portion R(x) ≃ 1 and R ′ (x) ≃ 0, so that the term in parentheses is equal to unity and has no influence on f . Considering this shape of R ′ , and by comparing the update equation in Equation (6.21) to the "masked RL" update in Equation (6.14), we can see that R ′ (Af ) play a role very similar to the mask z, downweighting saturated observations. However, it is important to note that this behaviour arises naturally from the forward model, and is not detected from the blurry image. Given our current estimate of the latent image f , the forward model tells us which blurry pixels are saturated, and automatically downweights those pixels in the update equation. We refer to the algorithm using this update rule as "saturated RL".

Figure 6.6 demonstrates the advantage of this method over the standard RL algorithm on a synthetic 1D example.

Preventing the Propagation of Errors

It is important to note that even when we take account of saturation in our image formation model, we are not necessarily able to estimate every pixel in f accurately. In the blurring process, each latent pixel f j in the sharp image is blurred across multiple pixels in the blurry image g. If some (or all) of these are saturated, we are left with an incomplete set of data concerning f j , and our estimate of f j is likely to be less accurate than if we had a full set of unsaturated observations available. This mis-estimation is one source of "ringing" artefacts in the deblurred output. An over-estimate at one pixel must be balanced out by an under-estimate at a neighbouring pixel, which must in turn be balanced out by another over-estimate. In this way, an error at one pixel spreads outwards in waves across the image. In order to mitigate this effect, we propose a second modification to the Richardson-Lucy algorithm to stop the propagation of these errors.

0 0.5 1 1.5 2 0 0.5 1 x R(x) R ′ (x) (
Following the arguments above, we see that it is the brightest pixels in f , and their neighbours, that we are likely to estimate poorly (since they cause saturation), and that we would like to prevent the propagation of errors from these poorly estimated pixels to their neighbours. To this end, we segment f into two disjoint regions: S, which includes the bright pixels and their neighbours that we are unlikely to estimate accurately, and U , which covers the rest of the image and which we can estimate accurately. We decompose the latent image correspondingly:

f = f U + f S .
Our aim is then to prevent the propagation of errors from f S to f U . To achieve this, we propose to estimate f U using only data which is not influenced by any pixels from S. To this end, we first define the region (denoted by V) of the blurry image which is independent of f S , by eroding U using the non-zero elements of the PSF: (6.24) where U T (k) denotes the set U transformed by T (k) . By taking the intersection of all the transformed versions of U , we ensure that V contains only those blurry pixels that are completely independent of S. We can then estimate f U using only the data in V, by defining the binary mask v which corresponds to V and adapting the update equation from Equation (6.14) for Richardson-Lucy with missing data (Section 6.2.1):

V = k:w k >0 U T (k) ,
f t+1 U = f t U • A ⊤ g • R ′ (A f t ) • v R(A f t ) + 1 -R ′ (A f t ) • v . (6.25)
We estimate f S using the previously defined "saturated RL" algorithm:

f t+1 S = f t S • A ⊤ g • R ′ (A f t ) R(A f t ) + 1 -R ′ (A f t ) . (6.26)
Since the Richardson-Lucy algorithm is an iterative process, we do not know beforehand which parts of f belong in U and which in S. We thus perform the segmentation at each iteration t using a threshold on the latent image:

U = j f t j ≤ ϕ . (6.27)
We decompose f t according to (6.29) where u is the binary mask corresponding to U . We then compute V, update fU and fS using Equations (6.25) and (6.26), and recombine them to form our new estimate of the latent image f t+1 = f t+1 U + f t+1 S . We refer to this algorithm as "combined RL", and Figure 6.6 shows the results of applying it to a synthetic 1D example, demonstrating the advantage over the standard RL and "saturated RL" algorithms. Algorithm 2 summarises the method. Compute set V of blurry pixels independent of S using (6.24)

f t U = u • f t (6.28) f t S = f t -f t U ,
5
Compute f t+1 U using only data from V (6.25)

6
Compute f t+1 S using all data (6.26)

7 f t+1 = f t+1 U + f t+1 S 8 end
Although this combined RL algorithm involves the use of a threshold to segment the image (Equation (6.27)), its effect is less dramatic than in Section 6.2.1. In this case, the threshold only determines whether a given pixel f j should be updated using all the available data, or a subset of it. This is in contrast to Section 6.2.1, where parts of the data are discarded and never used again. Since our aim is to ensure that no large errors are introduced in f U , we set the threshold low enough that most potentially-bright pixels are assigned to S. Empirically, we choose ϕ = 0.9 for the results in this paper.

Implementation

In this section we describe some of the implementation details of the proposed algorithm. When segmenting the current estimate of the latent image in the combined RL algorithm, we take additional steps to ensure that we make a conservative estimate of which pixels can be estimated accurately. First, after thresholding the latent image in Equation (6.27), we perform a binary erosion on U , such that (6.30) where ⊖ denotes binary erosion, and the structural element M used for erosion is a disk of radius 3 pixels. This ensures that all poorly-estimated pixels are correctly assigned to S (perhaps at the expense of wrongly including some well-estimated pixels too).

U = j f j ≤ ϕ ⊖ M,
Performing this step improves the deblurred results, since it is not only the bright pixels whose value is likely to be inaccurate due to saturation, but their neighbours too, and fewer artefacts arise from wrongly assigning a well-estimated pixel into S than the other way around. Second, in order to avoid introducing visible boundaries between the two regions, we blur the mask u slightly using a Gaussian filter with standard deviation 3 pixels to produce a smoother set of weights when extracting f t U and f t S from the current latent image f t in Equations (6.28) and (6.29).

Results

Figures 6.1 and 6.7 to 6.9 show results of non-blind deblurring using the proposed "combined RL" algorithm described in Section 6.3 on real hand-held photographs. The (spatially-variant) PSFs for these images were estimated from the blurry images themselves using the MAP-ℓ 1 algorithm described in Chapter 4. The only modification required to handle saturated images using this algorithm is to discard potentially saturated regions of the blurry image using a threshold. Since in this case the aim is only to estimate the PSF (and not a complete deblurred image), we can safely discard all of these pixels, since the number of saturated pixels in an image is typically small compared to the total number of pixels. There will typically remain sufficient unsaturated pixels from which to estimate the PSF.

Note that the standard Richardson-Lucy algorithm and the approach of [START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF] produce large amounts of ringing around the saturated regions, while the proposed "combined RL" algorithm avoids this with no loss of quality elsewhere. In Comparison to the method of [START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF]. The blurry images, the spatially-invariant PSFs and the results of their method were provided by the authors. spatial frequencies, and numerical round-off and noise may dominate over the true signal. This immediately gives rise to ringing, however this kind of ringing can often be suppressed with regularisation.

Regularisation in non-blind deblurring is able to alleviate ringing to some extent, however since most regularisers only concerns image derivatives, they are only capable of reducing ringing at high spatial frequencies. [START_REF] Yuan | Progressive inter-scale and intra-scale nonblind image deconvolution[END_REF] observe this fact and propose a multi-scale non-blind deblurring algorithm capable of preventing ringing caused by noise and numerical inaccuracies. However their method does not handle ringing caused by outliers in the data, shown in Figure 6.12.

The other factor to consider is that the zero singular values represent spatial frequencies which are very poorly constrained by the data. Thus if the blurry image contains outliers which do not fit the linear model, any disruption to the deblurred image will happen in these spatial frequencies first, since they are the cheapest place to hide the errors.

Ringing Due to Outliers

We can see how outliers in the data cause ringing by decomposing the observed image g into an "ideal" noiseless component g * , a noise component n, and an outlier component e. For simplicity, we will consider spatially-invariant blur with a blur kernel a, such that the noiseless blurry image g * is related to the true sharp image f by g * = a * f . We decompose the observed blurry image as The component fg * , which is deconvolved from g * will be close to the true sharp image, apart from some smoothing due to the regularisation. The magnitude of the noise n is typically small relative to the image, and so the component fn will also have a fairly small amplitude in the deblurred image. The outliers e however, may have a large magnitude, and are also likely to have a high bandwidth -e.g. will contain impulses or step edges which have large frequency components at a wide range of frequencies. This is a significant source of ringing, since those spatial frequencies for which the power spectrum of a is small (in the denominator of Equation (6.33)) may be amplified. Since the outliers have such a large bandwidth, it is very likely that a large frequency component of the outliers will coincide with a small frequency of the kernel. Figure 6.13 demonstrates this visually with a synthetic example.

g = g * + n + e. ( 6
This also explains why ringing appears around sharp edges, since such features in an image will have a broad frequency bandwidth, and even small errors may cause large amplification of the zero frequencies of a.

Ringing Due to Kernel Errors

When the kernel contains errors, a second source of ringing is introduced. We denote the estimated kernel as â, and the true kernel as a, and again write Fourier transform of the deblurred image f as (6.36) In this case, we need only examine fg * , since the effect of deconvolving n and e with the incorrect kernel is essentially the same as deconvolving them with the correct kernel.

F( f ) = F(â) * F(â) * • F(â) + φ • F (g * ) + F (n) + F (e) (6.35) f = fg * + fn + fe .
F( fg * ) = F (f ) • F(â) * • F(a) F(â) * • F(â) + φ (6.37) = F (f ) • 1 + ∆a k F(â) * • F(â) F(â) * • F(â) + φ - F(â) * • F(∆a) F(â) * • F(â) + φ . (6.38)
We can write the deconvolved image as the sum of two terms, from the two terms in the parentheses, i.e.

F( f (1) g * ) = F (f ) • 1 + ∆a k F(â) * • F(â) F(â) * • F(â) + φ (6.39) F( f (2) g * ) = -F (f ) • F(â) * • F(∆a) F(â) * • F(â) + φ (6.40) fg * = f (1) g * + f (2)
g * (6.41) Examining the spectra of the deconvolved image fg * , and the two components f (1) g * and f (2) g * , we can see that f (1) g * is very similar to the original sharp image (slightly brighter due to the (1 + ∆a k ) scale factor), while the ringing is caused by f (2) g * , which amplifies small frequency components of â:

(d) |F ( fg * )| (e) |F ( f (1) g * )| (f) |F( f (2) g * )| (g) Deconvolved image fg * (h) f (1) g * (i) f (2)
g * Figure 6.14. Synthetic example showing ringing due to kernel errors. By decomposing the deconvolved image as described in Equations (6.39) to (6.41), we can see that the ringing is caused when the kernel error ∆a contains frequency components which are small in the estimated kernel â. .

The first component f (1)

g * is close to the true sharp image scaled by (1 + ∆a k ), except with some frequency components reduced due to the regularisation φ. The second component f (2) g * is the source of the ringing, which arises when the error ∆a has a broad spectrum (e.g. delta functions) which aligns with the small frequency components of the corrupted kernel â. Figure 6.14 demonstrates ringing caused by kernel errors on a synthetic example.

Implications for Blind and Non-Blind Deblurring

In the preceding sections we have examined ringing caused by outliers in the blurry image, and errors in the kernel. What do these insights mean for deblurring in practice?

In terms of errors in the kernel, the good news is that unless the errors are large (i.e. the kernel is completely wrong) the ringing will have a small amplitude. If we want to improve our estimate of the kernel, the deblurred image will typically contain the correct edges and texture for us to do so. This can be seen in Figure 6.14 (g), where the step edges and ridges that are used by the blind PSF estimation algorithm of [START_REF] Cho | Fast motion deblurring[END_REF] are clearly visible, and the ringing does not interfere too much with this.

In an iterative PSF estimation algorithm, the next estimate of the kernel would almost certainly be closer to the true kernel.

For ringing caused by outliers in the data, our conclusions must depend on the sources of these outliers. When the true forward model is non-linear, and we assume a linear forward model (as in this section and in most work on deblurring), then outliers will occur anywhere the true forward model is non-linear. Since these outliers are caused by an error in the model rather than a random noise process, the chances of the outliers having structure, and thus potentially a wide bandwidth are much greater. In order to reduce the ringing caused, we can either improve our models to take account of the non-linearities, or attempt to find and discard the non-linear data. These approaches are described in Section 6.2, and this is the approach to handling saturation proposed by [START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF]. As shown in Section 6.3 however, this is typically insufficient, since the loss of information will cause inaccuracies in the deblurred image, which then propagate as ringing. Thus in our "combined RL" algorithm we take the additional step of attempting to determine which pixels can be estimated accurately, and decoupling their estimates from the poorly-estimated pixels.

When outliers in the data are caused by noise such as "hot pixels", dirt or scratches, or deleted pixels, a non-linear model will not help to prevent ringing. In these cases, the outliers must be prevented from affecting our estimate of the sharp image, and this is the approach taken by [START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF]. In this chapter we have not addressed deblurring in the presence of such outliers, however we envision that our "combined RL" approach could be extended to handle such phenomena.

Conclusion

In this chapter we have developed an approach for the problem of blind deblurring for images blurred with camera shake and suffering from saturation. We have proposed a non-blind deblurring algorithm, based on the Richardson-Lucy algorithm, which is able to deblur images containing saturated regions without introducing artifacts, without sacrificing detail in the result.

Chapter 7

Removing Occluders from Photos of Famous Landmarks

Introduction

Often when we review our holiday photos, we notice things that we wish we could have avoided, such as vehicles, construction work, or simply other tourists. To go back and retake the photograph is impossible, so what can we do if we want to remove these things from our photos? We wish to replace the offending pixels in as convincing a way as possible, a problem often referred to as image completion or inpainting, which has been widely studied recently. However, most existing methods are only applicable to small image regions, or require significant user intervention. We would ideally like to have no limit on the complexity of the image in the regions we are replacing, and would like to be able to fill pixels with whatever would actually have been observed there had the occlusions not occurred, with minimal user interaction. With the growth of online photo sharing websites, we have millions of photos of the world available to us, probably including many of the very same place as our own photo, so how can we use them to fix up our snaps?

In the same way that strong priors are essential to successful blind deblurring algo- rithms, it is not possible to replace a region of an image convincingly without some prior knowledge of the replaced pixels' appearance. While previous authors have focussed on statistical approaches or intuition to model the local behaviour of images (as discussed in Section 2.5), in general we may wish to replace a region of an image that contains structures that are too complex to be reproduced using these approaches. In this chapter we leverage recent advances in viewpoint invariant image search [START_REF] Chum | Total recall: automatic query expansion with a generative feature model for object retrieval[END_REF][START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF] to find other images of the same scene. Beginning with a query image g containing a target region Ψ to be replaced, our aim is to replace the pixels in the target region Ψ in a way that is convincing and corresponds to the true underlying scene. We first use an online image search engine to retrieve images of the same scene, which we refer to as a set of oracles. We use this set of oracles to propose a set of solutions {p (q) } to the filling problem. Finally we solve a labelling problem to decide how best to combine the proposals into a single output result f . We begin by describing the viewpoint invariant image search used in Section 7.2, before describing our method of geometrically and photometrically registering each search result to the query image in Section 7.3. In Section 7.4 we discuss how to combine the proposals generated from all the retrieved images into a single output image. We present some results and failure cases for our system in Section 7.5, and conclude in Section 7.6

Retrieving Oracle Images

Recent works [START_REF] Chum | Total recall: automatic query expansion with a generative feature model for object retrieval[END_REF][START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF][START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF] have demonstrated the feasibility of large-scale image-based search engines, allowing the user to input an image region as a search query, and be presented with other images depicting the same object or scene. Such viewpoint invariant image search is now becoming commercially available, e.g. Google goggles1 . In this work we use the "Oxford Buildings

Search"2 online demonstration provided by [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF]. The system uses a set of 100,000 images collected from Flickr3 , for which an indexing structure has already been computed. Both the query images and oracle images from which we construct the solution come from this set. Figure 7.2 shows a typical query image and the top 30 results returned by the search engine. These search results are the oracles which are used to replace the target region in the query image g. For popular landmarks it is reasonable to expect the top retrieved results to contain the building in the query image, and for some to have been taken from a similar viewpoint. These search results will be used to replace the target region of our query image.

Geometric and Photometric Registration

Having retrieved a set of oracle images depicting the same place as the query image, we first need to deduce their geometrical and photometric relationship to the query image. We would also like to reason jointly about their registration, in order to be able to combine them into a single, coherent, occlusion-free image. In this section, we first explain the use of homographies to register the oracle images to a query image, and the extension to multiple homographies that is necessary to register scenes with multiple planes. We then find regions which have been registered well by each homography, and use these regions to estimate a global transformation on the colour channels of each oracle. Finally, we group together homographies which are likely to have registered the same scene plane, and for each group compute an unoccluded "median" image, which will be used to guide the region selection process in Section 7.4. The output of this part In each column, the first row shows a query photograph and another photograph of the same scene returned by the image search engine. Below, we show the second image registered to the first using a homography, and the two images overlaid. The last row shows a close-up of the overlaid images. Left column: A pair of images with approximately the same camera centre -the transformation applies almost equally well to all parts of the image. Right column: A pair of images of a piecewise planar scene -the homography registers the dominant wall, where the variation in depth is small, but is no longer valid for points not on that plane, such as the protruding wall on the right.

Multiple homographies

For a scene containing multiple planes, a single homography will in general be insufficient to register an oracle to the whole query image. Thus we allow multiple homographies to be detected for each oracle. After estimating a homography between the two images, we remove all the inlying interest points from consideration, and run RANSAC again on the remaining putative matches. We repeat this process until the last estimated homography is rejected by the criteria above.

Ground plane registration

While simply removing inliers to previous homographies is sufficient for RANSAC to register the different planes of buildings, it often fails to register the ground plane, where there are generally very few interest points, which can get lost amongst the large number of putative matches. However, given that many occlusions are caused by people or objects standing on the ground, it may be important for the ground plane to be well registered. In the situation where the ground plane is not registered and causes significant problems, we allow the user to interactively segment the ground region in the query and oracle images, and run RANSAC again using only the interest points in that region. This process is illustrated in Figure 7.4, where removing the putative matches from other parts of the image enables RANSAC to find the homography which registers the ground plane. Given the mask of well-registered pixels, we estimate the linear photometric correction by taking the median ratio of gradient magnitudes between the query image g and the registered oracle o, in the well-registered region, i.e. median i ( (∇g) i / (∇o) i ),

Photometric registration

where (∇g) i is the magnitude of the gradient of image g at pixel i. The correction is performed separately in each colour channel.

Grouping homographies

The goal here is to group homographies corresponding to a particular scene plane, and subsequently obtain an unoccluded "median" image for each group, which will be used to guide selection of unoccluded and well-registered oracles for inpainting the target region (Section 7.4).

In order to group together homographies which correspond to a particular scene plane, we consider each oracle-homography combination as a node of a graph, and for two homographies having as inliers the sets of interest points S and Q, place an edge between the two nodes if the overlap of these sets is above a threshold: |S∩Q| |S∪Q| > 0.1. We consider each connected component of the graph as a group of homographies likely to register the same scene plane.

To compute the median image for a homography group, we follow the approach proposed by Weiss ( 2001) for computing occlusion-free "intrinsic images". For a group of homographies G, we first compute the x and y derivatives of all the registered oracles in that group. At each pixel we take the median x derivative and the median y derivative over all the registered oracles, before using Poisson blending [START_REF] Pérez | Poisson image editing[END_REF] The query image with each homography group's inlying interest points marked, some of the registered oracles from each group, with the regions considered to be well-registered highlighted, and (far right) the "median" image for each group within the target region. Note that for each group, the median image provides a sharp, unoccluded estimate of the relevant plane, while it is blurry elsewhere. Thus, the difference between a registered oracle and this image will be small where the oracle is well-registered and unoccluded, but large elsewhere .

to reconstruct the "median" proposal p median(G) i for that group. combined into a single result which could not have been achieved using any single oracle.

Reasons for this might be occlusions in the oracles themselves, or the fact that a single homography may not be able to register the whole target region. Figure 7.6 shows the advantage of combining multiple proposals; a single oracle provides most of the result but requires other oracles to provide some small parts, partly due to the mis-registration of the ground plane.

In order to decide which proposal should be used at which pixel, we want to design and optimise an energy function which encourages each pixel to choose well, but to regularise this with the idea that pixels should agree with their neighbours about what their neighbourhood should look like. We can consider this a labelling problem, where the label l i for a pixel i corresponds to which proposal is used there. This can be formulated as a multi-label conditional random field (CRF), where we wish to find the optimal label configuration l by solving a problem of the form min

l i∈Ψ E 1 i, l i + (i,i ′ )∈E E 2 i, i ′ , l i , l i ′ , (7.1)
where Ψ indicates the set of pixels in the region being solved, (i, i ′ ) indicates a pair of neighbouring pixels (4-neighbours), with E being the set of all such pairs in the region being solved. E 1 i, l is the "cost" of using the proposal p (l) at pixel i, encoding our wishes for individual pixels, while E 2 i, i ′ , l i , l i ′ is the cost of using proposals p (l i ) and p (l i ′ ) at neighbouring pixels i and i ′ , encoding the way we wish neighbouring pixels to agree with each other.

Pixels outside the target region should look similar to the original query image, since they lie outside the region originally specified for replacement. Pixels on the inside of the target region however should be similar to some robust estimate of the unoccluded scene, to avoid inserting new occlusions into the image. To achieve these two goals we choose E 1 to have the form

E 1 i, l i = k query m i p (l i ) i -g i + k median m i p (l i ) i -p median(G(l i )) i , (7.2)
where m is a binary mask indicating the target region, and m is its logical negation.

Outside the target region, where m i = 0, the cost depends on the difference between p (l i ) i

(the colour of proposal l i at pixel i) and g i (the colour of the query image g at the pixel). This term penalises any differences between the result and the input query image, while allowing the optimisation to choose the best location for the boundary of the replaced region. Inside the target region, where m i = 1, the cost depends on the difference between p

(l i ) i and p median(G(l i )) i
, which is the "median" image for that oracle's homography group G(l i ) as described in Section 7.3.5. This term serves a dual purpose, helping both to avoid inserting new occlusions into the result, and avoid using any proposals outside the region where they are registered, since in both these cases, the deviation p

(l i ) i -p median(G(l i )) i
should be large. Further to this we set E 1 i, l i to a large number if proposal l i does not cover pixel i. This is effectively a hard constraint which prevents a registered oracle being used for regions outside its bounds. The parameters k query and k median weight the terms according to their relative importance. The norm p ig i is simply the Euclidean distance in RGB space.

The purpose of E 2 is to encourage a few large regions to be combined instead of many small ones, and to ensure that boundaries between regions from different proposals occur at places where they will be least obvious. For this we use the "gradient" cost function suggested by [START_REF] Agarwala | Interactive digital photomontage[END_REF], where E 2 = 0 if l i = l i ′ , and otherwise (7.3) where (∇p) i is the concatenation of the image gradients at pixel i in all colour channels, i.e. a 6D vector. The first term penalises the difference between the two proposals' gradients at pixel i, and the second term penalises the difference between the two proposals' gradients at pixel i ′ . This cost helps prevent boundaries between regions using different proposals, if the two proposals' image gradients differ at that location. By encouraging these boundaries to fall in places where the image gradients in the two proposals match well, the transitions between regions are hidden. For the results in this paper we used k query = 1, k median = 1, and k grad = 10, and for computational purposes we used at most 10 proposals from each homography group (obtained as in Section 7.3.5), ranked by the number of homography inliers.

E 2 i, i ′ , l i , l i ′ = k grad (∇p (l i ) ) i -(∇p (l i ′ ) ) i + (∇p (l i ) ) i ′ -(∇p (l i ′ ) ) i ′ ,
Finally, in order to achieve a good transition between the region that has been filled and the query image surrounding it, the query image itself is included as a proposal, and the CRF is solved over a larger region than the original target region. By generating the proposals such that they extend outside the target region, the optimisation may choose the best place in this boundary region at which to switch from the query image to the proposals. To optimise the CRF described above, we use tree-reweighted belief propagation [START_REF] Kolmogorov | Convergent tree-reweighted message passing for energy minimization[END_REF][START_REF] Szeliski | A comparative study of energy minimization methods for Markov random fields[END_REF][START_REF] Wainwright | MAP estimation via agreement on trees: message-passing and linear-programming approaches[END_REF], using the software made available online5 by [START_REF] Szeliski | A comparative study of energy minimization methods for Markov random fields[END_REF].

Results and Discussion

In this section, we demonstrate our method on several query images from the Oxford Buildings image set, with various difficulties which our method is able to overcome.

In Figure 7.1, the region to be replaced is spanned almost entirely by a single scene plane. In this case, the image search returned 49 photos of the same scene, of which 15 were accurately geometrically registered to the query image using single homographies.

The result comes mostly from a single oracle, which was automatically selected by our CRF formulation since it provided a consistent, occlusion free replacement for the target region.

Figure 7.6 shows the advantage of our method's ability to combine multiple proposals.

The search returned 48 photos of the same scene, of which 9 were registered to the query image, providing 16 proposals due to multiple homographies. Most of the result comes from a single proposal, but the other proposals are used to fill some regions where this fails demonstrating the benefit of using multiple oracles. The final composite obtained by our method is significantly better than results obtained by the approaches of Criminisi 2003) and (i) the result using the method of [START_REF] Hays | Scene completion using millions of photographs[END_REF]. In our result, we obtain a much better image by combining regions from several proposals, and choosing the boundaries to be as inconspicuous as possible. The result of Criminisi et al.

propagates strong edges into the target region, but cannot reason about where best to connect them or insert details that are hidden (e.g. the grate at the base of the wall). The method of Hays and Efros produces a result that fits well with the surrounding image, but does not correspond to the actual underlying scene. We are grateful to Hays and Efros for running their system on our input image to produce the result in (i). The result shown in Figure 7.9 demonstrates the effect of the unary term E 1 defined in Equation (7.2). When the weight k median is low, the pairwise smoothness term has a greater effect, and the optimisation chooses to combine fewer proposals, while following the median proposal less closely. When the weight k median is high, the optimisation finds a solution using more proposals, and avoiding introducing new occlusions into the result.

Figure 7.10 shows a typical failure case for our system, which tends to occur when we have only a limited number of oracles available. In this case, there may not be sufficient images to construct an occlusion-free result. Even if there are enough oracles, the median image used to guide the final solution may be misleading, since it will have been computed from relatively few samples, which may include new occlusions of misregistered planes.

Conclusion

We have demonstrated an inpainting method which is able to combine images taken from different viewpoints and under different lighting conditions to convincingly replace large regions of the query photograph containing complex image structures, thereby significantly extending the capabilities of local patch-based inpainting methods [START_REF] Criminisi | Object removal by exemplar-based inpainting[END_REF][START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF].

The approach is mainly applicable to tourist snapshots of frequently visited landmarks, which have many pictures taken by other people available on the Internet. Although results are visually pleasing, in some cases under a close inspection subtle artefacts remain in the final composites. These are mainly due to photometric and resolution issues such as differences in the length of shadows and different image resolution/focus. Chapter 8

Perspectives

In this chapter we review the contributions of this thesis, and discuss possible directions for future research.

Contributions

In Chapter 3 we presented a new model for spatially-variant blur caused by camera shake.

We showed that camera shake blur is predominantly caused by rotations of the camera during the exposure, as opposed to translations, and that such blur is generally spatiallyvariant. Starting from the geometry of camera rotations, we modelled the blurry image as a weighted combination of projectively-transformed versions of the sharp image. We demonstrated a practical, discrete version of this model, which parameterises the pointspread function (PSF) of a "shaken" image using a single set of weights.

In Chapter 4 we demonstrated that the model of spatially-variant blur proposed in Chapter 3 can be applied within existing deblurring algorithms, originally designed to estimate only spatially-invariant PSFs. We show that the resulting algorithms are able to estimate spatially-variant PSFs, and that by applying our model in non-blind deblurring algorithms, we can deblur a wider range of images than using a spatiallyinvariant blur model. We successfully applied our model to both a marginalisation and a MAP algorithm for blind PSF estimation, as well as in the case where an additional noisy (but sharp) image of the scene is available.

In Chapter 5 we discussed the increase in computational cost incurred by using a spatially-variant blur model instead of a spatially-invariant one. While spatiallyinvariant blur can be computed quickly in the frequency domain using the fast Fourier transform (FFT), there is no equivalent for our proposed blur model. We proposed an approximation scheme for our spatially-variant blur model that divides the blurry image into a small number of overlapping sub-images, each of which is treated as having spatially-invariant blur. Using this scheme we provide the equations for quickly approximating the most intensive computations in the blind deblurring process. We demonstrated an order of magnitude speed-up of blind deblurring compared to the exact model, while retaining the accuracy and consistency of our global model from Chapter 3.

In Chapter 6 we extended our forward model for camera shake blur to include sensor saturation, which introduces a non-linearity into the image formation process. We derived a modified version of the Richardson-Lucy algorithm which incorporates this non-linear forward model and is able to recover intensities that lie outside the camera's dynamic range. We proposed a second modification to the Richardson-Lucy algorithm which explicitly handles saturated pixels separately from unsaturated ones, and prevents ringing artefacts by decoupling the update equations for the two sets of pixels. Using the proposed algorithm, we demonstrated non-blind deblurring of images with saturated pixels without introducing large ringing artefacts.

Finally in Chapter 7 we tackled the problem of removing large occluders from photographs of popular landmarks. We proposed a system that for popular locations, is able to realistically hallucinate the occluded scene content, without any user input and without an explicit 3D representation of the scene. Our system processes "oracle" images retrieved from a visual search engine by geometrically and photometrically registering them to the target image, before collecting a set of proposals for how to fill the occluded region. In the final step, the proposals are combined using a CRF formulation, which encourages a result free from occlusions, and without any visible region boundaries.

Future Work

One important avenue for future research is to improve the robustness and reliability of blind deblurring algorithms. The automatic algorithms used in Chapter 4 for blind PSF estimation sometimes fail, for example if the blurry image contains too much noise, or if the sharp image does not match the prior that is used. One potential line of work is to investigate the use of more powerful image priors, for example that proposed by [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF], which models the joint distributions of pixels in image patches, rather than just pairs of pixels. Another interesting possibility is to consider the role of user interaction in blind deblurring. For many people, it is natural to interact with an image restoration process, and such a system would have the potential to be more robust to difficult images. Sometimes in these images, the shape of the PSF or the location of strong edges can be perceived visually, and it may be possible for a user to provide some assistance, allowing the PSF to be estimated. The best way of integrating user

interaction is an open question, however some potentially-natural modes of interaction might be to (a) draw the approximate shape of the PSF at one or more locations in the image, (b) draw the locations of strong step edges, or (c) draw a rectangle around locations where the shape of the PSF is visible (e.g. trails left by bright spots).

Another possibility for increasing the robustness of the blind deblurring process is to take advantage of the trails left by bright lights in the image. Many images taken at night include bright electric lights, and in Chapter 6 we addressed non-blind deblurring of images that contain such saturated regions. However, the trails left by bright lights in blurry images often give a clear outline of the PSF (although the intensity is lost due to saturation), and could provide a powerful constraint for blind PSF estimation.

In this thesis, we have considered camera shake blur using a model with 3 degrees of freedom, leading to 3-dimensional blur kernels. However, in some cases it may be useful to consider additional degrees of freedom, such as camera translation (for example when photographing from a fast-moving vehicle) or varying focal length (for example if the photographer accidentally nudged a manual zoom lens while taking a photograph). In these cases the dimensionality of the kernel would grow, and it may not be practical to estimate the high-dimensional kernel explicitly. Once the kernel dimension reaches a certain level, it may no longer provide the most compact parameterisation of the PSF.

Instead, it might be easier to parameterise the PSF using a set of local filters, such as in Chapter 5. The challenge to making such a parameterisation effective is the need to enforce global consistency between the filters, without knowing the global kernel itself.

This would be done with pairwise constraints between all pairs of filters, and could be approached in a purely numerical way, simply exploiting the fact that each filter is related to the kernel as a (r) = J r w. Alternatively, it may be possible to utilise the underlying geometry, for example to find a sort of epipolar geometry between pairs of filters.

variable, the parameters of its distribution depend on the distributions of all the other latent variables, e.g. w

(1) k and w

(2) k depend on q(β σ ), q(f j ) for all j and q(w k ′ ) for all k ′ = k. For our non-uniform blur model, we find the following optimal values for the parameters, as given earlier in Equations (4.2) to (4.5) (cf . Equations (46-49) of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF]):

w (2) k = β σ i j T (k) ij f j 2 q(f ) (A.4) w (1)
k w

(2)

k = β σ i g i j T (k) ij f j q(f j ) - k ′ =k j T (k) ij f j j T (k ′ ) ij f j q(f ) w k ′ q(w k ′ ) (A.5) f (2) j = β σ i k T (k) ij w k 2 q(w) (A.6) f (1) j f (2) j = β σ i g i k T (k) ij w k q(w k ) - j ′ =j f j ′ q(f j ′ ) k T (k) ij ′ w k k T (k) ij w k q(w)
.

(A.7)

The details of the derivation are given next.

A.1 Variational method

For convenience, we will collect the latent variables f , w, and β σ into the "ensemble"

Θ. The aim is to approximate the true posterior p(Θ|g) with a simpler factorized distribution q(Θ|g), denoted for simplicity as q(Θ) = q(β σ ) j q(f j ) k q(w k ). Our model, from Equation (3.12), provides the likelihood p(g|Θ):

g * i = j,k w k T (k) ij f j , (A.8) p(g|Θ) = i G g i ; g * i , β -1 σ ,
see Equation ( 7) of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF] (A.9)

Similarly, since q(f ) = j q(f j ), q(w) = k p(w k ), p(f ) = j p(f j ), and p(w) = k p(w k ),

C KL = j q(f j ) ln q(f j ) p(f j ) df j + k q(w k ) ln q(w k ) p(w k ) dw k + q(β σ ) ln q(β σ ) p(β σ ) dβ σ -q(Θ) ln p(g|Θ) dΘ. (A.16)

Finally, we expand the last term: .16) and ignoring terms independent of Θ, C KL = j q(f j ) ln q(f j ) p(f j ) df j + k q(w k ) ln q(w k ) p(w k ) dw k + q(β σ ) ln q(β σ ) p(β σ ) dβ σ -1 2 q(β σ ) i ln β σ -β σ q(f )q(w)(g i -g * i ) 2 df dw dβ σ (A.22)

p(g|Θ) = i G(g i ; g * i ,

A.3 Optimal Distributions

A.3.1 Optimal q(β σ )

To derive the optimal form of q(β σ ), we ignore terms in C KL independent of β σ , add a Lagrange multiplier for the constraint that q(β σ ) dβ σ = 1, and differentiate with respect to q(β σ ) :

C KL q(β σ ) = q(β σ ) ln q(β σ ) p(β σ ) -1 2 i ln β σ -β σ (g i -g * i ) 2 dβ σ + λ σ q(β σ ) dβ σ -1 , (A.23)

where • denotes the expectation under the approximating distribution q(Θ).

C KL q(β σ ) + δq(β σ ) = q(β σ ) ln q(β σ ) + δq(β σ ) p(β σ ) -1 2 i ln β σ -β σ (g i -g * i ) 2 dβ σ + δq(β σ ) ln q(β σ ) + δq(β σ ) p(β σ ) -1 2 i ln β σ -β σ (g i -g * i ) 2 dβ σ + λ σ q(β σ ) dβ σ + δq(β σ ) dβ σ -1 (A.24) ln q(β σ ) + δq(β σ ) ≃ ln q(β σ ) + δq(β σ ) q(β σ ) to first order, so (A.25)

C KL q(β σ ) + δq(β σ ) = C KL q(β σ ) + δq(β σ ) dβ σ + δq(β σ ) ln q(β σ ) p(β σ ) + δq(β σ ) q(β σ ) -1 2 i ln β σ -β σ (g i -g * i ) 2 dβ σ + λ σ δq(β σ ) dβ σ .

(A.26)

Discarding higher order terms in δq, δC KL = δq(β σ ) 1 + ln q(β σ ) p(β σ ) -1 2 i ln β σ -β σ (g i -g * i ) 2 + λ σ dβ σ (A.27)

∂C KL ∂q(β σ ) = 1 + ln q(β σ ) p(β σ ) - 1 2 i ln β σ -β σ (g i -g * i ) 2 + λ σ . (A.28)
For convenience, we partition f into f j , the pixel of interest, and f j * , the remaining pixels:

C KL q(f j ) = q(f j ) ln q(f j ) p(f j ) + 1 2 β σ i q(f j * )q(w)(g ig * i ) 2 df j * dw df j + λ j q(f j ) df j -1 (A.34)

C KL q(f j ) + δq(f j ) = q(f j ) ln q(f j ) p(f j ) + δq(f j ) q(f j ) + 1 2 β σ i q(f j * )q(w)(g ig * i ) 2 df j * dw df j + δq(f j ) ln q(f j ) p(f j ) + δq(f j ) q(f j ) + 1 2 β σ i q(f j * )q(w)(g ig * i ) 2 df j * dw df j + λ j q(f j ) df j + δq(f j ) df j -1 (A.35) = C KL q(f j ) + (1 + λ j ) δq(f j ) df j + δq(f j ) ln q(f j ) p(f j ) + 1 2 β σ i q(f j * )q(w)(g ig * i ) 2 df j * dw df j (A.36)

∂C KL ∂q(f j ) = ln q(f j ) p(f j ) + 1 2 β σ i q(f j * )q(w)(g ig * i ) 2 df j * dw + 1 + λ j .

(A.37)

Setting this equal to zero, we obtain the optimal form ln q(f j ) = ln p(f j ) -1 2 β σ i q(f j * )q(w)(g ig * i ) 2 df j * dw -1λ j .

(A.38)

Here we need to make some simplifications to obtain a function of f j . For convenience we re-write the forward model in Equation (3.11) (p. 64) as (A.39)

g * i = f ⊤ C i w,

  Figure1.1. Removing camera shake blur from photos. The blur caused by camera shake, such as in the blurry image shown here (a), is typically spatially-variant. Most previous work on removing camera shake has assumed a spatially-invariant blur model, leading to deblurred images such as that shown in (b). Using the model for spatially-variant camera shake blur proposed in Chapter 3 and the blind deblurring algorithms described in Chapter 4, we are able to model the spatially-variant blur correctly and obtain superior deblurred results (c). Despite the large blur (around 30 pixels), much of the text that is illegible in the input image, and which is not restored sufficiently with the spatially-invariant blur model, can be read clearly in our deblurred image.

  Figure 1.2. Removing occluders from photos. This figure shows an example result produced by the system described in Chapter 7. The tourists are removed and replaced with a faithful rendition of the underlying scene, as can be seen by comparison with another image of the scene, shown in (d).

Figure

  Figure 1.3. Deblurring a real saturated image. Taken from Canon's "Your Second Shot" advertising campaign. The original image (a) is very blurry, however using blind PSF estimation (Chapter 4) to estimate how the image was blurred, followed by non-blind deblurring, handling the saturated pixels (Chapter 6) we are able to obtain a much better image (b). For the advertising campaign, Canon flew the couple back to Barcelona to take the photo again (c) with a new camera.

  Proceedings of the 20th British Machine Vision Conference. London, 2009 O. Whyte, J. Sivic, A. Zisserman, and J. Ponce. Non-uniform deblurring for shaken images. In Proceedings of the 23rd IEEE Conference on Computer Vision and Pattern Recognition. San Francisco, CA, 2010 O. Whyte, J. Sivic, and A. Zisserman. Deblurring shaken and partially saturated images. In Proceedings of the IEEE Workshop on Color and Photometry in Computer Vision (CPCV 2011), with ICCV 2011. Barcelona, Spain, 2011a O. Whyte, J. Sivic, A. Zisserman, and J. Ponce. Non-uniform deblurring for shaken images. International Journal of Computer Vision. Accepted, 2011b

  Figure 2.1. The Poisson noise distribution. The likelihood (a) and negative log-likelihood (b) of different observed pixel values g i , as a function of the noiseless value g * i .

  Figure 2.2. The Gaussian noise distribution. The likelihood (a) and negative log-likelihood (b) of different observed pixel values g i , as a function of the noiseless value g * i .

  Figure 2.3. Statistics of, and priors for image gradients in photographs. (a) shows a photograph of a street scene, from Fergus et al. (2006). In (b), the logarithm of the histogram of gradients is plotted, along with several sparse-derivatives priors.

  26)[START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF];[START_REF] Levin | Image and depth from a conventional camera with a coded aperture[END_REF];[START_REF] Tappen | Exploiting the sparse derivative prior for super-resolution and image demosaicing[END_REF] explore the use of the hyper-Laplacian distribution where φ(x) = |x| p , for some exponent p, typically chosen to be 0.5 ≤ p ≤ 0.8.[START_REF] Shan | High-quality motion deblurring from a single image[END_REF] use a piecewise prior where φ(x) = min(k |x| , ax 2 + b), where the parameters k, a, and b are set manually. Figure 2.3 plots the derivatives of a natural image, and various distributions which have been used to fit this empirical distribution. The histogram is characterised by a sharp peak at zero, with a slow fall-off.

  Figure 2.4. Statistics of, and priors for camera shake blur kernels. Levin et al. (2009) estimated a set of "true" camera shake convolution kernels (a) using sharp / blurry pairs of images. In (b), the logarithm of the histogram of all the kernel values in (a) is plotted, along with two sparsity-inducing priors.

  encourage sparsity by placing a mixture of exponentials prior on the kernel values, such that p(a) = k D d=1 π d exp(-λ d a k ), (2.27) where π d is the mixture weight for the d th exponential, with rate parameter λ d . Fergus et al. (2006) learn the parameters {π d } and {λ d } from a set of real blur kernels. Shan et al. (2008) use a single exponential distribution for each of the kernel elements, such that p(a) = k λ exp(-λa k ).

  Figure 2.5. Non-blind deblurring with Poisson noise. In this synthetic example, the original sharp image (a) is blurred and corrupted with Poisson noise to produce the degraded image (b).The deblurred result (c) is produced using the Richardson-Lucy algorithm. Since no regularisation is used, noise is amplified and artefacts appear in the deblurred result (e.g. the vertical stripes to the right of the table). In the regularised Richardson-Lucy result (d), the noise and artefacts are reduced, although some image features are also suppressed.

  Figure 2.6. Non-blind deblurring with Gaussian noise. In this synthetic example, the original sharp image (a) is blurred and corrupted with Gaussian noise to produce the degraded image (b).The deblurred results are produced using the algorithms of (c)[START_REF] Levin | Image and depth from a conventional camera with a coded aperture[END_REF], and (d)[START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF].

  Figure 2.7. Infinite solutions to the blind deblurring problem. Given a blurry image (a), there are an infinite number of possible PSFs and sharp images (b), (c), (d) that could have produced the blurry image. It is particularly challenging to find the correct solution amongst all these possible solutions, and strong priors on both the sharp image and PSF are required.

  Figure 2.8. An example of blind deblurring by marginalisation or MAP algorithms. For the blurry image (a), the single-image deblurring algorithms of Fergus et al. (2006) and Cho and Lee (2009) estimate a convolution kernel, and subsequently deblur the image using a nonblind deconvolution algorithm, such as Richardson-Lucy. For this image the results produced by the marginalisation approach of Fergus et al. and the MAP approach of Cho and Lee produce similar results.

  Figure 2.9. An example result of the algorithm of Yuan et al. (2007b). Using a blurry image (a) and a noisy image (b) as inputs, Yuan et al. (2007b) estimate a spatially-invariant blur kernel (shown in (d)). Performing non-blind deconvolution using the standard Richardson-Lucy algorithm produces a deblurred image (c) containing visible ringing artefacts. Yuan et al. (2007b) propose several modifications to the Richardson-Lucy algorithm, using the noisy image as a guide to produce a deblurred image (d) free from ringing.

  Figure 2.10. An example of spatially-variant blur due to camera shake. This real handheld image exhibits a blur due to camera shake which is spatially-variant. This can be seen in the close-ups (b).

  [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF] propose a model of the form in Equation (2.61), using a coarse grid of overlapping rectangular regions to cover the image. They propose to use the Bartlett-Hann window as m (

Figure 2 .

 2 Figure 2.11. The inpainting problem. This image shows the notation for the inpainting problem. The known region of the images is denoted by M, while the unknown target region is denoted by Ψ. The pixels in Ψ that lie along the boundary with M are denote by δΨ.

  Figure 3.1. Modelling 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 work, we can describe this blur using a single global descriptor (c), which in this case has been estimated from the blurry image itself, simply by modifying existing algorithms for blind deblurring (see Chapter 4 for details). Close-ups of different parts of the image (d) 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 the deblurred image in (b) (see Chapter 6 for details on deblurring images like this with saturated pixels) and the close-ups in (e), different parts of the image, blurred in different ways, can be deblurred to recover a sharp image.

  Figure3.2. Blur 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 centre. 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).

  Figure 3.3. Real measurements of camera motion during a long exposure. These photographs, from[START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF], are taken with a 1/3 second exposure, and show how the pose of the camera varied during the exposure. The motion is recovered using an additional high-speed camera attached to the main camera -by performing structure-from-motion on the high-speed frames, the camera's pose can be recovered. From the motion plots, we see that the camera rotates between ±0.5 • , and translates between ±4 mm while the shutter is open.

  Figure3.4. Our coordinate frame with respect to initial camera orientation, and the paths followed by image points under different camera rotations. We define our coordinate frame (a) to have its origin at the camera's optical centre, 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 centre, and the pixels are assumed to be square.

  Having defined the type of image transformations we expect to occur while the shutter is open, we can write out the image degradation model. Let T denote the exposure time of the photograph. While the shutter is open, the camera passes through a sequence of orientations θ t , t ∈ [0, T ]. As discussed above, at each pose θ t , the sensor is exposed to a projectively transformed version of the sharp image f , where the projective transformation H t is given by Equations (3.3) to (3.5). The noiseless blurry image g * is then modelled as the integral over the exposure time T of all the transformed versions of f :

Figure 3

 3 Figure 3.1 (c) shows a visualisation 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 visualisation, with white corresponding to a large value, and black corresponding to

  .15) which covers the set of 2D translations and rotations considered in their work. By including the focal length in our model, we ensure accuracy even at short focal lengths.

  uniform blur, allowing them to handle non-uniform blur: those of Fergus et al. (2006) (Section 4.3.1) and Cho and Lee (2009) (Section 4.3.2).

  results of the marginalisation algorithm using these modified update equations and our spatially-variant blur model, compared to the original algorithm of[START_REF] Fergus | Removing camera shake from a single photograph[END_REF].

Figures 4 .

 4 Figures 4.1 and 4.2 show blind deblurring results on images blurred by real camera shake. Our model, used in both the marginalisation and MAP algorithms, is able to capture and remove the blur, while the original algorithms of Fergus et al. and Cho

Figure 4 .

 4 Figure 4.3 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

  Figure 4.2. Blind deblurring of real camera shake, example 2. A hand-held image with camera shake (a), captured with a shutter speed of 1 second, with the results of blind deblurring using the marginalisation 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 estimated kernels for the two models are shown inset in the deblurred results.The variational output (d) is estimated using the marginalisation 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 artefacts 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) and (f) cover ±0.7 • in θ X and θ Y and ±1.4 • in θ Z .

Figure 4 . 4 .

 44 Figure 4.4. Blind deblurring of synthetic single-axis blurs. A sharp image (top left) withexamples 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 summarised 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.

Figure 4 . 5 .

 45 Figure 4.5. Poor performance of MAP-ℓ 2 with non-uniform blur model. The corresponding blurry image can be seen in Figure4.2. Shown is the estimated kernel and deblurred result when using our non-uniform blur model in the algorithm of Cho and Lee with ℓ 2 regularisation on the kernel. As can be seen, the ℓ 2 regularisation is not sufficient to produce a good estimate of the kernel, and results in a deblurred output containing many artefacts.

Figure 4 .

 4 Figure 4.5 shows the failure of the MAP algorithm to produce a good result (using the blurry image from Figure 4.2 (a)) when using the original ℓ 2 regularisation proposed by Cho and Lee (2009) with our non-uniform blur model. As discussed in Section 4.3.2, the kernel produced is highly non-sparse despite the thresholding step, and the deconvolved output correspondingly exhibits many artefacts compared to the MAP-ℓ 1 result in Figure 4.2 (f).

  Figure 4.6. Blind deblurring of a real uniform blur. A real camera shake blur (a)-(b) from the dataset of[START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF], deblurred using kernels estimated with the marginalisation algorithm. We show deblurred results and kernels for four cases; (c) uniform blur using original algorithm ofFergus 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 (e), while not perfect, does have the same diagonal shape as the true blur, with the non-zeros centered around a single value of θ Z .

Figure 4 . 7 .

 47 Figure 4.7.Comparison of MAP-ℓ 1 method with that of[START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF]. A handheld 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). The blur in this image is approximately 35 pixels long, and 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. To estimate the kernel for this image, we found that was necessary to reduce the threshold on the kernel values (Section 2.3.3) to 1 100 the maximum value, instead of

  Figure 4.8. Deblurring 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 deblurred results for the uniform (c) and non-uniform (d) blur models. As can be seen in the close-ups (bottom two rows), our result contains more details and fewer artefacts than when using the uniform blur model, and reveals features not visible in either the noisy or the blurry image. The non-uniform kernel in (d) covers ±3 • along each dimension.

Figure 5 .

 5 Figure 5.2 plots the blind PSF estimation times for several images with different sizes of blur, using both the exact forward model and the fast approximate forward model.

Figure

  Figure 5.2. Blind deblurring time using the exact and approximate model. This log-log plot shows the time taken to perform blind deblurring of several images at various different sizes, using either the exact forward model (solid lines) or the approximate forward model (dashed lines). The blind deblurring algorithm is implemented in Matlab, and the function to compute the exact forward model is implemented in compiled C, while the approximate model is implemented in Matlab. The blind deblurring time using the exact model grows significantly faster with respect to the image size than using the approximate model, and for images 1024 pixels wide already takes almost 2 orders of magnitude longer. The images for these times are from the following figures: Istanbul -Figure6.8 (p. 125), LivingRoom -Figure 5.3 (p. 104), Pantheon -Figure 4.2 (p. 78), Book -Figure 4.3 (p. 79).

  Figure 5.2. Blind deblurring time using the exact and approximate model. This log-log plot shows the time taken to perform blind deblurring of several images at various different sizes, using either the exact forward model (solid lines) or the approximate forward model (dashed lines). The blind deblurring algorithm is implemented in Matlab, and the function to compute the exact forward model is implemented in compiled C, while the approximate model is implemented in Matlab. The blind deblurring time using the exact model grows significantly faster with respect to the image size than using the approximate model, and for images 1024 pixels wide already takes almost 2 orders of magnitude longer. The images for these times are from the following figures: Istanbul -Figure6.8 (p. 125), LivingRoom -Figure 5.3 (p. 104), Pantheon -Figure 4.2 (p. 78), Book -Figure 4.3 (p. 79).

  Figure 5.3. Least-squares non-blind deblurring using the exact and approximate forward models. Given a blurry image of size 1024 × 768 and blur kernel (a), this figure shows the results and computation times for least-squares deconvolution with ℓ 2 gradient regularisation (Equation (5.3)), using (b) conjugate gradient descent (CG) with the exact forward model, (c)CG with the approximate forward model, and (d) direct deconvolution using the approach described in Section 5.3.2. The results are visually similar using all three methods. Using CG with the approximate forward model is much faster than with the exact model, however the direct approach takes only a fraction of the time of either of these.

  Figure 6.1. Deblurring in the presence of saturation. Existing deblurring methods, such as those in (b) and (c), do not take account of saturated pixels. This leads to large and unsightly artefacts in the results, such as the "ringing" around the bright lights in the zoomed section. Using the proposed method (d), the ringing is greatly reduced and the quality of the deblurring improved. The PSF for this image causes a blur about 35 pixels in width, and was estimated using the MAP-ℓ 1 algorithm from Chapter 4.

Figure 6 . 3 .

 63 Figure 6.3. Ignoring saturated pixels using a threshold. A simple way to handle saturation is to threshold the blurry image at some level T , and discard the blurry pixels above this threshold. Shown are the results of running the Richardson-Lucy algorithm for different thresholds.As the threshold decreases, the artefacts around the bright lights are reduced compared to the standard RL result in Figure6.1 (b). At the lowest threshold (d) the fewest artefacts appear, but parts of the church are also discarded, hence remain blurred.

FigureFigure 6 . 4 .

 64 Figure 6.4. Diagram of image formation process. We model the response function as acting on the noiseless blurry image, before the noise is added.

  Figure 6.5. Modelling the saturated sensor response. (a) Ideal clipped linear response function (solid blue line) and its derivative (dashed red line). The derivative is not defined at x = 1. (b) Smooth and differentiable approximation to (a) defined in Equation (6.22). The derivative is also smooth and defined everywhere.

Figure 6 . 6 .

 66 Figure 6.6. Synthetic example of blur and saturation. Each column shows a sharp "tophat" signal, blurred using the box filter shown at the top left. Gaussian noise is added and the blurred signal is clipped, to model saturation. The kernel is also degraded with noise and one large error to produce a "perturbed" kernel which is used for deconvolution, to simulate errors in the kernel estimation step. The last three rows show the deblurred outputs for three algorithms discussed in Section 6.2. (a) With no saturation, all three algorithms produce similar results. (b) When some of the blurred signal is saturated (region B), the standard RL algorithm produces an output with large ringing artifacts. Although region A is not itself saturated, the ringing propagates outwards from B & C across the whole signal. The "saturated RL" algorithm reduces the ringing and correctly estimates the height of the top-hat at its edges (region C), where there are some unsaturated observations available. In region B all information about the height of the sharp signal is lost, and the output takes a sensible value close to 1. (c) When the blurred top-hat is completely saturated, it is no longer possible to estimate its true height anywhere. The saturated RL result accurately locates the top-hat, but contains ringing. The proposed method (combined RL) mitigates this by preventing the propagation of errors to the non-saturated region (compare D to E).

  Figure 6.7. Deblurring saturated images. Note that the ringing around saturated regions, visible in columns (b)-(c) is removed by our method (d), without causing any loss in visual quality elsewhere.

  .31) Using direct deconvolution in the frequency domain with an ℓ 2 regularisation term φ (e.g. φ = αF (d x ) * • F (d x ) + αF (d y ) * • F (d y )), we can write Fourier transform of the deblurred image f asF( f ) = F (a) * • F (g) F (a) * • F (a) + φ (6.32) = F (a) * F (a) * • F (a) + φ • F (g * ) + F (n) + F (e) (6.33) f = fg * + fn + fe .(6.34)

  Figure 6.12. Comparison to Cho et al. (2011) and Yuan et al. (2008). This figure compares non-blind deblurring results, on saturated images, for (b) the algorithm of Cho et al. (2011), (c) the Richardson-Lucy-based algorithm of Yuan et al. (2008), and (d) our proposed "combined RL" algorithm. While both our algorithm and the algorithm of Yuan et al. are based on the Richardson-Lucy algorithm, our results contain much less ringing, due to handling the saturated pixels explicitly. Compared to the method of Cho et al., our results contain similar or less ringing. Results in (b) and (c) provided by Cho et al.

  Figure 6.13. Synthetic example showing ringing due to outliers. By decomposing the deconvolved image into the contributions from the "noiseless" blurry image g * , the noise n, and outliers e, we can see that the ringing is caused entirely by the outliers. By examining the spectra of the outliers and the blur kernel, we can trace the causes of the ringing to places in the frequency domain where small components of the kernel coincide with large components of the outliers.

  Figure 7.1. An example result from our algorithm. The replaced region is consistent with the rest of the original image and the boundary between is effectively hidden, producing a convincing result. Note the complexity of image structures on the inpainted facade.

  Figure 7.1 shows an example result obtained from our system.Image completion using images of the same scene from the Internet was independently of our work approached by[START_REF] Amirshahi | Photo completion using images from internet photo sharing sites[END_REF];[START_REF] Amirshahi | An image completion algorithm using occlusion-free images from internet photo sharing sites[END_REF], who 7.2 Retrieving Oracle Images obtain relevant images using a text-based image search, then combine unoccluded blocks from them to replace the target region using a greedy algorithm. Our novel contributions are the following: (a) candidate images are obtained automatically using viewpoint invariant image search, (b) geometric registration is performed using multiple homographies, which allows registering significantly more complex scenes captured from varying camera viewpoints, and (c) occlusion reasoning and seamless combination of multiple registered oracle images is formulated as a single labelling problem efficiently solved using tree-reweighted belief propagation.

Figure 7 . 2 .

 72 Figure 7.2. Two example queries and the first 30 results returned by the viewpoint invariant image search engine.For popular landmarks it is reasonable to expect the top retrieved results to contain the building in the query image, and for some to have been taken from a similar viewpoint. These search results will be used to replace the target region of our query image.

Figure 7 . 3 .

 73 Figure 7.3. Pairs of images related by homographies.In each column, the first row shows a query photograph and another photograph of the same scene returned by the image search engine. Below, we show the second image registered to the first using a homography, and the two images overlaid. The last row shows a close-up of the overlaid images. Left column: A pair of images with approximately the same camera centre -the transformation applies almost equally well to all parts of the image. Right column: A pair of images of a piecewise planar scene -the homography registers the dominant wall, where the variation in depth is small, but is no longer valid for points not on that plane, such as the protruding wall on the right.

  Due to the diverse origins of the images retrieved from the Internet, they may have very different lighting conditions from the query image. To reduce the effect of lighting variations when combining oracle images with the query, we work in the gradient domain, and furthermore estimate a global linear correction on the gradients of each colour channel of each oracle image. The correction is robustly estimated by first estimating which regions are well-registered in each oracle, in a lighting-invariant way. To do this we compute the normalised cross-correlation (NCC) between the square (15 × 15) patch around each pixel in the query image with the corresponding square patch in the transformed oracle image. To convert the resulting NCC scores into a mask of well-registered regions, we use a binary Markov random field (MRF) with unary costs which are linear in the NCC scores, and a Potts potential on adjacent pixels, and minimise the cost function using graph-cuts. Figure7.5 shows typical binary masks indicating the well-registered regions for each oracle.

  Figure 7.4. Semi-automatic ground plane registration. This figure shows (a) a query image, (b) an oracle image, (c) the first homography extracted with inliers shown, and (d) the second homography extracted after the user manually indicates ground plane region (below horizontal line). The first homography extracts the dominant plane, and by manually indicating the ground plane region RANSAC is able to register the ground plane in the second homography.

  Figure 7.5 shows an example of two homography groups extracted for a scene with two dominant planes, and the median proposals for each.

  Figure 7.6. Combining multiple proposals. The inputs to the system (a) & (b) and the output labels (c) showing the combination of proposals in the final result (colours correspond to the borders of the images below). (d) The top 5 original (unregistered) automatically retrieved oracle images used in the result. (e) The geometrically registered oracles. (f) The proposals generated by each oracle. Note that none of the individual proposals covers the entire target region and provides a satisfactory result. Last row: (g) Our final result, (h) the result using the algorithm of Criminisi et al. (2003) and (i) the result using the method of[START_REF] Hays | Scene completion using millions of photographs[END_REF]. In our result, we obtain a much better image by combining regions from several proposals, and choosing the boundaries to be as inconspicuous as possible. The result ofCriminisi et al. propagates strong edges into the target region, but cannot reason about where best to connect them or insert details that are hidden (e.g. the grate at the base of the wall). The method of Hays and Efros produces a result that fits well with the surrounding image, but does not correspond to the actual underlying scene. We are grateful to Hays and Efros for running their system on our input image to produce the result in (i).

  Figure 7.7. Example Result 1.An example query and the corresponding result from our system. The CRF optimisation correctly chooses oracles registered to the ground plane where relevant (red, blue), and combines other oracles registered to the wall (green, yellow, magenta) in order to complete the result .

  Figure 7.9. Example result demonstrating the effect of unary term. In this example, it is possible to see the effect of the unary term, which encourages the solution to look similar to the median proposal. With the unary term at its default value (d)-(f) many proposals are combined and the alcove shown in the close-up is empty. When we lower the weight k median (g)-(i), the pairwise smoothness term has a greater effect, and the optimisation prefers to use fewer proposals, leading to a person appearing in the alcove, transferred from one of the proposals.

Figure 7

 7 Figure7.10. A failure case of the system. In this example, there are only 4 oracles (e)-(h) registered correctly in the target region. As a result, the solution cannot contain a good unoccluded view of the region, and instead contains a combination of the scaffolding from (e) and some mis-registered oracles (not shown).

http://yoursecondshot.usa.canon.com/

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/

http://cg.postech.ac.kr/research/deconv_outliers/ (accessed November 12,

2011).

http://www.google.com/mobile/goggles

http://www.robots.ox.ac.uk/ ~vgg/research/oxbuildings/

http://www.flickr.com

http://vision.middlebury.edu/MRF/code/

Acknowledgements

Deblurring With Noisy / Blurry Image Pairs

For comparison, we have also implemented this algorithm for uniform blurs, using convolution in Equation (4.20).

To perform the non-blind deblurring using the estimated kernel, we apply the augmented RL algorithm of Yuan et al. (2007b) directly, replacing convolution with our spatially-variant blur model.

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., 2007b).

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 .

Results and Discussion

In this section, we present results with noisy / blurry image pairs, and refer the reader to Section 4.5 for implementation details. Figures 4.8 and 4.9 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. (2007b). As can be seen from the deblurred images obtained with the two models, our results exhibit more detail and fewer artefacts than those using the uniform blur model.

Blind deblurring with the aid of a noisy image is significantly more robust than single image deblurring. This is because it involves solving two convex optimisation problems, all results in this paper we performed 50 iterations of all Richardson-Lucy variants. [START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF] Figure 6.11 shows the results of non-blind deblurring using our algorithm, alongside those of the recently-proposed algorithm of [START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF]. The blurry images and their spatially-invariant PSFs are provided by Cho et al., along with their deblurred results 1 .

Comparison to

The results of [START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF] are smooth from heavy sparse-gradient regularisation.

Our results are produced without any regularisation, however, despite this, they exhibit fewer visible artefacts in most cases than those of [START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF].

Perspective: The Causes of Ringing

In this section we offer some analysis of the causes of "ringing" artefacts in deblurred images -medium-frequency ripples spreading across the image. Ringing is common in deblurred images, and has been explained variously as due to the Gibbs phenomenon (Yuan et al., 2007b) (arising from the inability of a finite Fourier basis to reconstruct perfect step edges), or a combination of inaccurately-modelled image noise and errors in the estimated PSF [START_REF] Shan | High-quality motion deblurring from a single image[END_REF]. We have found the latter to be accurate, and provide some analysis of this here.

Ringing is caused by the interaction between zeros in the PSF and (a) outliers in the data (e.g. caused by noise, non-linearities in the camera response, non-static scene) and (b) errors in the estimated blur kernel. By zeros in the PSF, we are referring to singular values of A that are zero or very small. Each singular value of A corresponds approximately to particular spatial frequencies (they correspond exactly for spatiallyinvariant blur, where the singular value decomposition (SVD) of A is equivalent to the discrete Fourier transform of a). Thus small or zero singular values correspond to spatial frequencies that are annihilated by the blur.

When we attempt to deblur the image, there is very little information about these of the system is a set of geometrically and photometrically registered oracles {o (q) }, and a set of unoccluded "median" images which are used as guides in the final part of the system. Note that since we allow multiple homographies, each retrieved oracle may appear multiple times in the set of registered oracles {o (q) }.

Registering two images with a homography [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]) is valid when the optical centres of the two cameras coincide, or when the scene being imaged is planar, as shown in Figure 7.3. In outdoor urban scenes, the second condition is approximately satisfied very frequently. Additionally, when we have a large number of photographs of a scene taken by many photographers, we find that the first condition can also sometimes be approximately satisfied. Although we are unlikely to find a photo taken from exactly the same location, we are likely to find many for which the effects of parallax will be small.

Homography estimation

We use the standard method [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] of determining putative point matches between the two images (query and oracle) and estimating the inliers and homography simultaneously using RANSAC [START_REF] Fischler | Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography[END_REF]. Here we use Harris-Affine [START_REF] Mikolajczyk | Scale & affine invariant interest point detectors[END_REF]) and SIFT [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] interest points (typically 5000-15000 per image) and SIFT descriptors. 4 Having estimated the homography using the inlying interest points, we check that the number of inliers is sufficient for it to be reliable, discarding any homographies with fewer than 50 inliers. We also check whether the line at infinity in the oracle becomes visible under the homography, and discard it if so, since we deem the transformation to be too extreme in this case. Finally, we retain the oracle only if it covers at least 25% of the target region after warping, since an oracle covering less than this is unlikely to be much use for replacing the region's contents.

Generating and Combining Proposals

Once we have an oracle geometrically and photometrically registered, we would like to use each registered oracle o (q) to generate a proposal p (q) for how the target region should be filled. The most direct way of doing this would be to simply copy the pixels from the oracle into the region, i.e. p (q) = m • o (q) + m • g, where m is a binary mask corresponding to the target region Ψ. However in practice the variations in lighting mean that this approach will produce poor quality results, with clear boundaries at the edge of the region. The problem of how best to combine two images whose properties do not necessarily match has been approached in many ways, from methods which aim to conceal boundaries between regions, such as Burt and Adelson's multiresolution spline [START_REF] Burt | A multiresolution spline with application to image mosaics[END_REF] and Poisson blending [START_REF] Pérez | Poisson image editing[END_REF], to methods that try and find the best place to locate the boundary itself, such as the dynamic programming of [START_REF] Efros | Image quilting for texture synthesis and transfer[END_REF] and the graph-cut technique of [START_REF] Kwatra | Graphcut textures: image and video synthesis using graph cuts[END_REF].

In this work we use Poisson blending to combine the images, whereby instead of combining pixels from the two images, their gradient fields are combined to form a composite gradient field ∇p (q) = m • ∇o (q) + m • ∇g. The composite gradient field ∇p (q) can then be reconstructed into an image p (q) by solving Poisson's equation. The query image provides Dirichlet boundary conditions (which constrain the colour of the solution) for the equation around the target region. If the transformed oracle does not span the entire target region, pixels bordering the remaining unfilled region (where we have no colour information) take Neumann boundary conditions (which constrain the gradient of the solution), in order to reduce colour artefacts. Figure 7.6 (f) shows some of the individual proposals generated for a query image using Poisson blending.

Combining Multiple Proposals

Following the steps described in the previous sections, we have a set of proposals {p (q) }, where each proposal is generated using a single registered oracle. However, it may be that individual oracles cannot provide the best solution when taken alone, but may be

Appendix A Parameter Update Derivation for Marginalisation Algorithm

In this appendix we derive the optimal forms and parameters of the approximating distributions q(f ), q(w) and q(β σ ), used in the marginalisation algorithm for spatiallyvariant blind deblurring in Chapter 4. For our spatially-variant blur model, the optimal distributions for the latent variables are the same as for spatially-invariant blur (cf .

Equations (42,43,17) of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF]):

where w

(1)

k , f

(1) j , f

(2) j are parameters of the distributions, g * i is the "noiseless" value of blurry pixel i, related to the unknown latent image f and kernel w through the forward model in Equation (3.11) (p. 64). Note that f and w are random variables, so in this context g * i is also a random variable. N is the number of observed blurry pixels, and • q represents the expectation with respect to the the distribution q. For each latent where G( • ; µ, σ 2 ) is a Gaussian with mean µ and variance σ 2 . In order to get to the posterior, we also need a prior p(Θ) for our latent variables. The latent variables are assumed to be independent, so that the prior factorizes:

and furthermore the elements of both f and w are assumed to be independent and identically-distributed, i.e.

From Equation ( 10) of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF], we wish to minimize the following cost function, first using the calculus of variations to find the optimal form of the approximate distributions, then iteratively optimizing their parameters, which is equivalent to minimizing the Kullback-Leibler (KL) divergence between the posterior and the approximating distribution (see Bishop, 2006, Equation. (10.3), page 463):

ln p(g|Θ) dΘ. (A.13)

A.2 Inside the Cost Function

Since q(Θ) = q(f )q(w)q(β σ ) and p(Θ) = p(f )p(w)p(β σ ), .15) Setting this derivative to zero, we obtain an relation similar to Equation ( 14) in [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF]:

Thus the optimal distribution is

which, given that p(ln β σ ) = 1, which implies that p(β σ ) = Γ(β σ ; ǫ, ǫ) with ǫ → 0 (Equation ( 8) of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF]), gives Equation ( 17) of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF]:

where the Γ distribution is given by Equation ( 15) of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF]:

Starting from Equation (A.22), and isolating the relevant terms, .33) where the N ×K matrix C i is obtained by re-arranging the elements of the transformation matrices T (k) . We split g * i into the contribution from j and j * :

where c ⊤ ij is the j th row of C i , and C ij * is C i with this row removed. .48) which is just a quadratic in f j . Replacing the coefficients with a j and b j ,

cf . Equation ( 43) of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF] where

cf . Equation ( 48) of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF] =

cf . Equation ( 49) of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF] =

.

(A.59)

A.3.3 Optimal q(w k )

We proceed much the same as for q(f j ), starting from Equation (A.22), and isolating the relevant terms,

similarly, we partition w into w k , the element of interest, and w k * , the remaining elements:

and obtain the optimal form ln q(w k ) = ln p(w k ) -1 2

Here we need to make some simplifications to obtain a function of w k , as in Appendix A.3.2.

where c ik is the k th column of C i , and C ik * is C i with this column removed. .68) q(w k ) ∝ p(w k ) exp -1 2 w

(2) k w kw

(1) k 2 (A.69) cf . Equation ( 42) of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF] where w

(2)

cf . Equation ( 46) of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF] = β σ i j

k w

(2)

cf . Equation ( 47) of [START_REF] Miskin | Ensemble learning for blind image separation and deconvolution[END_REF]. .73)