Oliver Whyte

Josef Sivic

Andrew Zisserman

Jean Ponce

Fergus

Hertzmann Roweis & Singh

Freeman

Jia & Agarwala

Cho Lee

Sivic Zisserman & Whyte

Ponce

Sivic Whyte

Efficient, Blind, Spatially-Variant Deblurring for Shaken Images

In this chapter we discuss modeling and removing spatially-variant blur from photographs. We describe a compact global parameterization of camera shake blur, based on the 3D rotation of the camera during the exposure. Our model uses three-parameter homographies to connect camera motion to image motion and, by assigning weights to a set of these homographies, can be seen as a generalization of the standard, spatially-invariant convolutional model of image blur. As such we show how existing algorithms, designed for spatially-invariant deblurring, can be "upgraded" in a straightforward manner to handle spatially-variant blur instead. We demonstrate this with algorithms working on real images, showing results for blind estimation of blur parameters from single images, followed by non-blind image restoration using these parameters. Finally, we introduce an efficient approximation to the global model, which significantly reduces the computational cost of modeling the spatially-variant blur. By approximating the blur as locally-uniform, we can take advantage of fast Fourier-domain convolution and deconvolution, reducing the time required for blind deblurring by an order of magnitude.

1.1

Introduction

Everybody is familiar with camera shake, since the resulting blur spoils many photos taken in low-light conditions. Camera shake blur is caused by motion of the camera during the exposure; 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, which is blurred as a result. We would like to be able to deblur such images to recover the underlying sharp image, which we would have captured if the camera had not moved.

In general, the problem of restoring an image after it has suffered some degradation can be broken down into three stages: first, a generative model is needed to relate the undegraded, "ideal" or latent image (that we would like to recover) to the observed image produced by the camera. Second, the parameters of this tions, and deriving a geometric model for camera shake. In Section 1.3 we develop this into a practical model for deblurring real images. 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

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, and both 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 In this simplified example, we consider capturing a blurred 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). Reproduced from [START_REF] Whyte | Non-uniform deblurring for shaken images[END_REF]) with permission, © Springer 2012.

(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 1.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 camera moves, and 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, (1.1)
where F is the camera's focal length, while in (b) the camera must rotate through an angle

θ = tan -1 δ F . (1.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 (1.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 (1.2) that we would need to rotate the camera by θ = tan -1 1 100 ≃ 0.6 • , independent of the subject's distance from the camera. To put this in terms of the motion of the photographer's hands, then for example if the camera body is 10 cm wide, such a rotation could be caused by moving one hand just 1 mm forwards or backwards relative to the other. Provided the subject is more than 1 metre from the camera, this motion is at least an order of magnitude smaller than for a translation of the camera 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.

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 will generally be small if the centre of rotation is not far from the optical centre. As we shall see in the following section, this model of camera motion leads to spatially-variant blur.

Motion blur and homographies

Under a pinhole camera model, and assuming that the scene being photographed is static, rotations of a camera about its optical centre induce projective transformations of the image being observed. In other words, the image observed at one camera orientation is related to the image at any other by a 2D projective transformation, or homography. For an uncalibrated camera, this is a general 8-parameter homography, but for a camera with known internal parameters, the homography H is specified by three parameters and is given by

H = KRK -1 , (1.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 1.2.3). The rotation matrix R has only three parameters. We adopt here the "angleaxis" parameterisation, in which a rotation is described by the angle θ moved about an axis a (a unit-norm 3-vector). This can be summarised in a single 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. Reproduced from [START_REF] Whyte | Non-uniform deblurring for shaken images[END_REF]) with permission, © Springer 2012.

3-vector θ = θa = (θ X , θ Y , θ Z ). R is then given by the matrix exponential

R θ = e [θ]× , where (1.4) [θ] × =   0 -θ Z θ Y θ Z 0 -θ X -θ Y θ X 0   . (1.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 XYplane 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 1. 3 (a). In this configuration, θ X describes the "pitch" of the camera, θ Y the "yaw", and θ Z the "roll", or in-plane rotation, of the camera.

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 (1.3) to (1.5). The noiseless blurred image g * is then modelled as the integral over the exposure time T of all the transformed versions of f :

g * (x) = T 0 f H t x dt, (1.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 1.3 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 non-uniformity of this out-of-plane motion is dependent on the focal length of the camera, decreasing as the focal length increases. However, it is typical for consumer cameras to have focal lengths of the same order as their sensor width, as is the case in Figure 1.3. 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 1.1).

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

g * (x) = f H θ x w(θ) dθ, (1.7)
where the weight function w(θ) encodes the camera's trajectory in a timeagnostic 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 (1.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 chapter, we assume that K takes the standard form

K =   F 0 x 0 0 F y 0 0 0 1   . (1.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. We recover the focal length of the camera from the image's EXIF tags, and assume that the principal point is at the centre of the image (which is typically sufficient for modelling homographies caused by camera rotations [START_REF] Szeliski | Image alignment and stitching: A tutorial[END_REF]).

The radial distortion present in many digital cameras can represent a significant deviation from the pinhole camera model. Rather than modelling the distortion explicitly, we pre-process images with the commercially available PTLens software1 to remove it. A second distortion present in many digital images comes from a non-linear intensity mapping applied by the camera storing the image, sometimes referred to as "gamma correction". In this work, we either use raw camera output images which do not have this non-linearity, or preprocess the blurred images with the inverse mapping.

1.2.4

Uniform blur as a special case

One consequence of our model for camera shake is that it includes spatiallyinvariant blur as a special case, and thus gives the conditions under which such a blur model is applicable. Given the definitions of R, K and H in Equations (1.3), (1.4) and (1.8), and assuming θ Z = 0, it can be shown [START_REF] Whyte | Removing Camera Shake Blur and Unwanted Occluders from Photographs[END_REF]) that for small θ X , θ Y , as F → ∞,

H θ →   1 0 F θ Y 0 1 -F θ X 0 0 1   , (1.9)
which simply amounts to a translation of the image. 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, the blur which results from camera motion will be approximately spatially-invariant.

The computational model

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 bilinear interpolation, whereby sub-pixel values of an image, say g(x) are interpolated as a linear combination of the four nearest pixels.

In this discrete setting, we write the blurred image as a linear function of the sharp image:

g * i = j A ij f j , (1.10) or in matrix-vector notation, g * = Af , (1.11)
where the N × N matrix A captures the discrete point spread function (PSF).

Each row of the matrix A corresponds to a single blurred pixel in g, and captures the fact that each blurred pixel is a weighted sum of multiple sharp pixels. In most cases of blur, the light received by each pixel in g comes from a relatively small number of nearby pixels in f . As a result, the PSF matrix A for an image is usually sparse (contains a relatively small number of non-zero values).

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.

Figure 1.1 (c) shows a visualisation of w, where the cuboidal volume of size N X ×N Y ×N Z is shown, with the points inside representing the non-zero elements of w in 3D. The kernel has also been projected onto the three back faces of the cuboid to aid visualisation, with white corresponding to a large value, and black corresponding to zero.

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

g * = k w k T (k) f , (1.12)
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 four non-zero elements corresponding to the interpolation weights for each pixel. By writing out Equation (1.12) for a single pixel, we obtain the discrete analog of Equation (1.7):

g * i = k w k j T (k) ij f j , (1.13)
where i and j index the pixels of the blurred image and the sharp image, respectively. For a blurred pixel g * i with coordinate vector x i , the sum j T (k) ij f j interpolates the value of the sub-pixel location f (H k x i ). We will return to how H k , and the corresponding matrices T (k) , are sampled in Section 1.7.

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

given w, g * = Af , where

A ij = k T (k) ij w k , (1.14) given f , g * = Bw, where B ik = j T (k) ij f j . (1.15)
In the first form, A ∈ R N ×N is a large sparse matrix, whose rows each contain a local blur filter acting on f to generate a blurred 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 sections.

In contrast to this model, when the PSF is spatially-invariant, we denote the discrete convolution kernel by a, and write

g * = a * f .
(1.16)

Blind estimation of blur from a single image

In this section, we demonstrate the effectiveness of the spatially-variant blur model presented in Section 1.3 by using it for blind deblurring of images acquired under camera shake. We consider single-image deblurring, where only a blurred image is available, and use our model to replace the spatially-invariant blur model in an existing algorithm for blind deblurring; specifically, applying our model within the algorithm proposed by [START_REF] Cho | Fast motion deblurring[END_REF]. "Blind" estimation of PSFs directly from images has a long history [START_REF] Gull | Maximum entropy method in image processing[END_REF][START_REF] Ayers | Iterative blind deconvolution method and its applications[END_REF][START_REF] Fish | Blind deconvolution by means of the Richardson-Lucy algorithm[END_REF], however the particular case of camera-shake blur has attracted significant attention recently [START_REF] Fergus | Removing camera shake from a single photograph[END_REF][START_REF] Shan | High-quality motion deblurring from a single image[END_REF][START_REF] Cai | Blind motion deblurring from a single image using sparse approximation[END_REF][START_REF] Cho | Fast motion deblurring[END_REF][START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF][START_REF] Xu | Two-phase kernel estimation for robust motion deblurring[END_REF][START_REF] Gupta | Single image deblurring using motion density functions[END_REF], Harmeling, Hirsch & Schölkopf 2010[START_REF] Krishnan | Blind deconvolution using a normalized sparsity measure[END_REF].

The method of [START_REF] Cho | Fast motion deblurring[END_REF] is similar in spirit to many other MAPtype algorithms recently proposed [START_REF] Shan | High-quality motion deblurring from a single image[END_REF][START_REF] Cai | Blind motion deblurring from a single image using sparse approximation[END_REF][START_REF] Xu | Two-phase kernel estimation for robust motion deblurring[END_REF][START_REF] Krishnan | Blind deconvolution using a normalized sparsity measure[END_REF], and proceeds by alternately updating the blur kernel and latent image, in a multi-scale framework. The algorithm iterates over three main steps. In the first step, a set of non-linear filters are applied to the current estimate f of the sharp image, in order to predict the locations and magnitudes of step edges in the sharp image. This is done using a bilateral filter [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF] to denoise the image, followed by a shock filter [START_REF] Osher | Feature oriented image enhancement using shock filters[END_REF] to enhance sharp edges, producing an image f ′ . Finally, the derivatives of f ′ are computed, and thresholded based on their orientation and magnitude, to produce a set of sparse edge maps {p (q) }.

In the second step, these predicted edges are used to estimate the blur kernel by solving a linear least-squares problem. Having found the kernel â that minimises this problem, any elements whose value are below a threshold are set to zero. This encourages sparsity in the kernel, and ensures that all the elements are positive. In the third step the sharp image is estimated, using the current estimate of the blur kernel â to deconvolve the blurred image and obtain an improved estimate of the sharp image f . This step is also performed by solving a linear least-squares problem.

These three 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. The kernel-update and imageupdate sub-problems involved in the algorithm of [START_REF] Cho | Fast motion deblurring[END_REF] are linear least-squares by virtue of the fact that convolution is a bilinear operation on the sharp image and the blur kernel. The fact that our blur model, given in Equation (1.13) is bilinear in the sharp image and blur kernel is the key feature that allows it to be applied within this, and other deblurring algorithms.

Although we have chosen the algorithm of [START_REF] Cho | Fast motion deblurring[END_REF] to apply our blur model, we note that many MAP-style blind deblurring algorithms consist of the same basic blocks, alternating between updating the kernel and the latent image by solving linear least-squares sub-problems, often in a multi-scale framework [START_REF] Shan | High-quality motion deblurring from a single image[END_REF][START_REF] Xu | Two-phase kernel estimation for robust motion deblurring[END_REF][START_REF] Krishnan | Blind deconvolution using a normalized sparsity measure[END_REF]). As such, the following discussion can apply equally well to a number of algorithms.

Updating the blur kernel

The general form of the spatially-invariant kernel update problem (the second step) in the algorithm of [START_REF] Cho | Fast motion deblurring[END_REF] is

min a a * f -g 2 2 + β a 2 2 , (1.17)
where the first data reconstruction term measures how well the blur kernel a applied to the current estimate of the latent image f reconstructs the observed blurred image g, and the second term, weighted by a scalar parameter β, is the ℓ 2 regularisation of the blur kernel a. Note that in [START_REF] Cho | Fast motion deblurring[END_REF] formulation, the data term replaces f with the sparse edge maps p (q) produced by the nonlinear filtering of the current estimate of the latent image f , i.e. has the form a * p (q)d (q) * g 2 2 , where q indexes over the first and second-order spatial derivatives, and d (q) are spatial derivative filters.

When applying our spatially-variant blur model, we modify the problem given by Equation (1.17) to update the spatially variant kernel w as

min w Bw -g 2 2 + β k w k , s.t. ∀k w k ≥ 0, (1.18)
where B is given by Equation (1.15). Here we have substituted our blur model into the data-reconstruction term, and replaced the ℓ 2 Tykhonov regularisation with ℓ 1 regularisation and non-negativity constraints on the kernel. Now, instead of simply a linear least-squares problem to update w, we have 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]. Similarly to [START_REF] Cho | Fast motion deblurring[END_REF] we use sparse edge maps p (q) which results, after expanding B using Equation (1.15), in a data term of the form

k w k T (k) p (q) -d (q) * g 2 2 .
The need for the ℓ 1 (instead of the simpler ℓ 2 ) regularisation arises from the differences between our kernel and convolution kernels. Fundamentally, our kernels cover a larger space of image transformations than convolution kernels (3parameter homographies instead of 2D image translations), and we must estimate more parameters from the same amount of data. As a result, the PSF estimation process is liable to become poorly conditioned, due to an increased amount of ambiguity in the data-reconstruction term. We have observed that when simply replacing the convolution in Equation (1.17) with our model, but without changing the regularisation, the resulting 3D kernels contain many non-zeros spread smoothly throughout, and do not produce good deblurred outputs (see Figure 1.5). On the other hand, with ℓ 1 regularisation, the optimisation is more likely to choose between ambiguous camera orientations than spreading non-zero values across all of them. In the remainder of the chapter, we refer to the original algorithm of Cho & Lee as MAP-ℓ 2 , and our ℓ 1 regularised version as MAP-ℓ 1 .

In addition to the use of ℓ 1 regularisation, we note that in order to constrain the 3D kernel adequately, we require data from all regions of the image. 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, or a rotation about its Z axis. In order to resolve this ambiguity, we would need to look at other regions of the image. To ensure that we use observations from all parts of the image when updating the kernel, we modify the construction of the edge-maps p (q) from the filtered sharp image f ′ . We simply subdivide the image into 3 × 3 regions, and apply the gradient thresholding step independently on each.

Updating the latent image

The latent image update (the third step) in the algorithm of [START_REF] Cho | Fast motion deblurring[END_REF] is performed in a similar way to many non-blind deblurring algorithms [START_REF] Xu | Two-phase kernel estimation for robust motion deblurring[END_REF][START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF]) by solving a linear least-squares problem of the form

min f a * f -g 2 2 + α d x * f 2 2 + α d y * f 2 2 , (1.19)
where d x and d y represent horizontal and vertical derivative filters, and α is a regularisation weight. Note that in [START_REF] Cho | Fast motion deblurring[END_REF] formulation, the data term also takes into account the derivatives of f and g as well as the intensities (i.e. a * (d (q) * f )d (q) * g 2 2 is also penalised, for various derivative filters d (q) ).

To apply our blur model, we simply replace the convolution in the data reconstruction term and update the latent image as

min f Af -g 2 2 + α d x * f 2 2 + α d y * f 2 2 , (1.20)
where A is given by Equation (1.14). We minimise Equation (1.20) using conjugategradient descent [START_REF] Shewchuk | An introduction to the conjugate gradient method without the agonizing pain[END_REF]. Note that at this point, 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 and to perform direct minimisation of Equation (1.19) in the frequency domain. However, in the following section we will describe an efficient approximation of the spatially-variant blur model which enables this.

Efficient computation of the spatially-variant model

Due to the additional computational expense incurred by using a spatiallyvariant blur model instead of a spatially-invariant one, both blind and non-blind deblurring under this model can be very time consuming. As seen in the previous section, MAP-type deblurring algorithms typically involve solving linear least-squares problems, minimising Afg 2 2 with respect to f , and Bwg with respect to to w. As defined in Equations (1.14) and (1.15), the matrices A and B involve sums over a set of homographies {T (k) }, where the size K of the set can be large, on the order of hundreds or thousands.

In iterative minimisation algorithms for solving such least-squares problems, we must repeatedly compute arbitrary matrix-vector multiplications involving A and B. For images with millions of pixels these matrices are generally too large to fit into memory, and the matrix-vector products must be computed on-the-fly, by explicitly warping the image for each homography T (k) . This is by far the biggest bottleneck in the use of our model, and means that both blind PSF estimation and non-blind deblurring are significantly slower than for spatially-invariant blur.

To reduce the running time of the whole deblurring process, in this section we propose an efficient approximation to the blur model from Equation (1.12), based on the locally-uniform "Efficient Filter Flow" proposed by [START_REF] Hirsch | Efficient filter flow for spacevariant multiframe blind deconvolution[END_REF]. We begin by describing the approximation in Section 1.5.1, before demonstrating in Section 1.5.2 how it can be used to compute the matrix-vector products necessary for the kernel update step (Equation (1.18)) very quickly. In Section 1.5.3 we describe how this approach allows us to update the sharp image (Equation (1.20)) directly in the frequency domain (instead of using iterative methods). The approximation allows blind deblurring to be performed an order of magnitude faster than when using the exact forward model. Note that concurrently with our work, Hirsch, Schuler, [START_REF] Hirsch | Fast removal of nonuniform camera shake[END_REF] proposed a similar model for efficiently computing spatiallyvariant camera-shake blur.

1.5.1

A locally-uniform approximation for camera shake [START_REF] Hirsch | Efficient filter flow for spacevariant 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. In these cases, it is reasonable to approximate spatially-variant blur as locally-uniform. Following this observation, they propose a model for spatially-variant blur, 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 overlapping patches ensure the blur varies smoothly, while allowing the forward model to be computed using P small convolutions. [START_REF] Hirsch | Efficient filter flow for spacevariant multiframe blind deconvolution[END_REF] assign each patch r a spatiallyinvariant blur filter a (r) , and their model is given by:

g * = P r=1 C (r) ⊤ a (r) * m • C (r) f , (1.21)
where C (r) is a matrix that crops the r th patch of the image f (thus C (r) ⊤ reinserts it via zero-padding). The vector m is the Bartlett-Hann window, and In their original work, [START_REF] Hirsch | Efficient filter flow for spacevariant multiframe blind deconvolution[END_REF] parameterise the blur using a separate filter a (r) for each patch r. Likewise [START_REF] Harmeling | Space-variant single-image blind deconvolution for removing camera shake[END_REF], who apply this model to 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 textureless regions, the algorithm of Cho & Lee may fail, and so heuristics are needed to encourage similarity between neighbouring filters, and also to detect and replace failed local kernel estimates.

•
Given the forward blur model for camera shake in Equation (1.12), 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 J (r) are simply a re-arrangement of the elements of the matrices T (k) ; J

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

g * = Af = Bw ≃ P r=1 C (r) ⊤ C (r) J (r) w * m • C (r) f . (1.22)
This equation allows the forward model to be computed quickly using only a handful of small convolutions, which can be performed efficiently in the frequency domain. Figure 1.4 shows how the quality of the locally uniform 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.

Updating the blur kernel

In iterative algorithms for estimating the blur kernel w from Equation (1.18), we typically need to compute the gradient of Bwg 2 2 with respect to w, which involves computing B ⊤ y for arbitrary y and B ⊤ B. In addition to being able to compute the forward model quickly, Equation (1.22) provides us with a fast approximate way of computing these products using one convolution per patch:

B ⊤ y ≃ P r=1 J (r) ⊤ C (r) ⊤ m • C (r) f ⊗ C (r) y (1.23) B ⊤ B ≃ P r=1 J (r) ⊤ C (r) ⊤ XCorrMatrix m • C (r) f ⊗ C (r) f C (r) J (r) , (

Updating the latent image: fast, non-iterative non-blind deconvolution

The locally-uniform approximation allows the forward model and its derivatives to be computed much faster using the FFT, such that we can quickly compute the A ⊤ y and A ⊤ A products needed to perform gradient-based optimisation of the sharp image in Equation (1.20). However, for spatially-invariant blur, the fastest way of updating the sharp image (in Equation (1.19)) is not using iterative methods such as conjugate-gradient descent, but rather using non-iterative frequency-domain deconvolution, which we outline below. In this section we extend this method to handle spatially-variant blur, via the locally-uniform approximation.

When blur is spatially-invariant, using Parseval's theorem, Equation (1.19) can be transformed into N independent 1D quadratic minimisations in the frequency domain, allowing the solution to be obtained directly by pixel-wise division in the frequency domain [START_REF] Gamelin | Complex Analysis[END_REF]):

f = F -1 F (a) * • F (g) F (a) * • F (a) + α F (d x ) * • F (d x ) + F (d y ) * • F (d y )
, (1.25) where F ( • ) takes the 2D discrete Fourier transform (computed using the fast Fourier transform), and F -1 ( • ) the inverse Fourier transform. However, for spatially-variant blur, the locally-uniform approximation does not immediately permit this. Even though each patch has a spatially-invariant blur, the fact that the patches overlap means that the reconstruction of any blurred pixel will involve several patches. This can be seen by inserting the locally-uniform model into the non-blind deblurring problem in Equation (1.20), and seeing that the sum over patches lies inside the data-reconstruction term:

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 . (1.26)
Thus, this equation cannot be minimised independently for each patch.

In order to be able to estimate each patch independently, one simple solution is to use Jensen's inequality to obtain an upper bound on Equation (1.26), taking the sum over patches outside the data-reconstruction term, and expanding the regularisation terms to also penalise individual patches. Having done this, we minimise the upper bound instead of the original cost. We define the set of blurred patches {g (r) } and sharp patches {f (r) } such that g (r) = m • C (r) g and f (r) = m • C (r) f , and estimate each deblurred patch f (r) by solving min

f (r) a (r) * f (r) -g (r) 2 2 + α d x * f (r) 2 2 + α d y * f (r) 2 2 (1.27)
for each patch. This can be done using the direct frequency-domain method of Equation (1.25). Second, we estimate the full deblurred image f that best matches the deblurred patches { f (r) } by minimising

P r=1 m • C (r) f -f (r) 2
2 . This problem can be solved independently for each pixel, yielding the following solution:

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

(1.28)

In Figure 1.7 we compare this method of non-blind deblurring to other possibilities for solving Equation (1.20) with spatially-variant camera shake blur. The fast independent method produces results which are visually very similar to those obtained using the exact model, in a significantly shorter amount of time. Using this direct deconvolution method to update the latent image via Equation (1.20) during blind deblurring provides an additional speed improvement, compared to the use of conjugate-gradient descent with the approximate forward model.

As a demonstration of the speed-up achievable with this approximation, the examples in Figures 1.5 to 1.7 all took more than 3 hours for blind deblurring using the exact model (implemented in Matlab and C), compared to under 6 minutes using the approximation presented in this section (implemented in Matlab). The results shown in this chapter were all produced using the approximation. The reduction in computational complexity can be quantified by comparing the exact model in Equation (1.13) with the approximate model in Equation (1.22). Evaluating Equation (1.13) requires O N K operations, where N is the number of pixels in the image and K is the number of non-zeros in w, whereas Equation (1.22) requires O N log N operations (due to the use of the FFT to compute the convolution). For a camera shake blur, O K ≥ O √ N , and thus the approximation provides a significant reduction in computational complexity. For a complete discussion of the computational complexity, as well as a full derivation of the equations presented here, please refer to [START_REF] Whyte | Removing Camera Shake Blur and Unwanted Occluders from Photographs[END_REF].

Application to other non-blind deblurring algorithms

The fast, non-iterative method for non-blind deblurring presented in this section can be used as a building-block in more sophisticated non-blind deblurring algorithms. A number of recent algorithms, which place sparse-gradient priors on the deblurred image [START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF][START_REF] Wang | A new alternating minimization algorithm for total variation image reconstruction[END_REF][START_REF] Shan | High-quality motion deblurring from a single image[END_REF], or which use more robust data-reconstruction terms [START_REF] Yan | An efficient TVL1 algorithm for deblurring multichannel images corrupted by impulsive noise[END_REF][START_REF] Xu | Two-phase kernel estimation for robust motion deblurring[END_REF], involve solving sub-problems of the form of Equation (1.19), which has a quadratic data-reconstruction term and quadratic regularisation. As such, the non-iterative method presented in this section can be used to extend these more sophisticated algorithms to handle spatially-variant blur, without a substantial increase in processing time. The results shown in Section 1.6 use the algorithm of [START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF], modified in this way, to perform the final non-blind deblurring.

Single-image deblurring results

In this section we present results of single-image deblurring using the MAPℓ 1 algorithm to estimate the spatially-variant PSF, with comparisons to results obtained with the original algorithm of [START_REF] Cho | Fast motion deblurring[END_REF] on real data. Having estimated the PSF, we apply 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 similar to Equation (1.20). Figure 1.5 shows a blind deblurring result on an image blurred by real camera shake. Our model is able to model and remove the blur, while the results with the original algorithm contain visible artefacts. 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. Also shown is the result of using ℓ 2 regularisation on the spatially-variant kernel. As discussed in Section 1.4.1, the kernel produced is not sparse, and as a result the deconvolved output exhibits many artefacts compared to the MAP-ℓ 1 result.

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

Additional

Limitations and failures

Since the MAP approach to blind deblurring attempts to solve a non-convex minimisation problem, it is not possible to guarantee a globally optimal solution. However, in practice we have found the MAP-ℓ 1 algorithm to be capable of deblurring a wide range of images with large blurs -the blurs removed in this chapter are up to 35 pixels wide (e.g. Figure 1.6) -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. 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. Nevertheless, failures do occur and Figure 1.8 shows an example case. We shows 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. Also shown is the result when using our blur model with ℓ2 regularisation on the kernel (c). 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. The rotational blur kernels in (c-d) cover ±0.7 • in θX and θY and ±1.4 • in θZ .

have observed failures caused by several factors, including a high level of noise in the input images, an excessively large blur, or an unknown non-linear camera response function e.g. Figure 1.8 (bottom row). Also, the algorithm may fail when the edge-based heuristics which guide the blind deblurring do not match

d) direct deconvolution using the approach described in Section 1.5.3. 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.

the particular image being estimated, e.g. in an image that contains only finescale texture with no large-scale step edges. Since we have assumed that camera translation has a negligible blurring effect, our model (and 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. As discussed in Section 1.2, this is unlikely to be a problem on most shaken images, except for close-up photos where the subject is less than about 1 m from the camera.

Implementation

The implementation of the algorithm of [START_REF] Cho | Fast motion deblurring[END_REF] is our own, and we use this implementation for both uniform and non-uniform blur models when comparing results. A binary executable for [START_REF] Cho | Fast motion deblurring[END_REF] 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. The implementation of the non-blind deblurring algorithm of [START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF] is based on Matlab code made available online by the authors2 .

Sampling the set of rotations One important detail to consider is how finely to discretise the orientation parameter θ. In the discrete case, each grid point θ (k) corresponds to a transformation matrix T (k) in the sum in Equation (1.12). Undersampling the space of orientations will affect our ability to accurately reconstruct the blurred image, but sampling it too finely will lead to unnecessary computation. 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. We set the size of the 3D kernel in terms of the size of the blur we are attempting to remove, typically a few degrees along each dimension of θ, e.g. [-5 • , 5 • ].

Multiscale implementation

Most successful blind kernel estimation algorithms here are applied within a multiscale framework, starting with a coarse representation of image and kernel, and repeatedly refining the estimated kernel at higher resolutions. In the case of single-image deblurring, this is essential to avoid poor local minima, however in our case, it is also important for computational reasons. At the original image resolution, the kernel may have thousands or tens of thousands of elements, however very few of these should have non-zero values.

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 clamping many kernel elements to zero, we eliminate a large amount of computation and memory requirements associated with estimating those elements' values. We first build Gaussian pyramids for the blurred image, 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.

Running time

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

Our implementation, in Matlab and C, of the MAP-ℓ 1 algorithm for spatiallyvariant blur takes over 3 hours to deblur a 1 megapixel (1024 × 768) image, depending on the size of the blur. This is a significant departure from the spatiallyinvariant algorithm of [START_REF] Cho | Fast motion deblurring[END_REF], who report deblurring times of under one minute using their C++ implementation. Using the efficient approximation described in Section 1.5, we are able to perform blind deblurring of the same images, using our spatially-variant blur model in under 6 minutes.

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 blurred image is bilinear in the sharp image and the weights.

We have applied our model for spatially-variant camera shake blur within an existing camera shake removal algorithms, and validated the model with experiments demonstrating superior results compared to the spatially-invariant blur model. We have also described how an efficient approximation for spatiallyvariant blur can be used to reduce the computational cost of computing the spatially-variant forward blur model.

Although we have demonstrated our model only in the algorithm of [START_REF] Cho | Fast motion deblurring[END_REF] for blind deblurring, it could equally be used with more recent algorithms, which have shown superior results, such as those of [START_REF] Xu | Two-phase kernel estimation for robust motion deblurring[END_REF] and [START_REF] Krishnan | Blind deconvolution using a normalized sparsity measure[END_REF]. Equally, faster or more sophisticated methods for non-blind deblurring, such as those of [START_REF] Afonso | Fast image recovery using variable splitting and constrained optimization[END_REF] and [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF], could be extended to use our spatially-variant blur model.

Interested readers can access an online demonstration of the blind deblurring algorithm at http://www.di.ens.fr/willow/research/deblurring/.

  (a) An image blurred by camera shake (b) Deblurred using the estimated blur descriptor

Figure 1 . 1

 11 Figure 1.1 Modelling non-uniform blur in a shaken image. The blurred 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 blurred image itself, simply by modifying existing algorithms for blind deblurring (see Section 1.4 for details). Having estimated the blur, the standard Richardson-Lucy algorithm is used to estimate the sharp image. 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) and the close-ups in (e), different parts of the image, blurred in different ways, can be deblurred to recover a sharp image. Reproduced from (Whyte et al. 2012) with permission, ©Springer 2012.

Figure 1 . 2

 12 Figure1.2 Blur due to translation or rotation of the camera. In this simplified example, we consider capturing a blurred 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). Reproduced from[START_REF] Whyte | Non-uniform deblurring for shaken images[END_REF]) with permission, © Springer 2012.

  Figure1.3 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. Reproduced from[START_REF] Whyte | Non-uniform deblurring for shaken images[END_REF] with permission, © Springer 2012.

  • • represents the Hadamard (element-wise) product. Note that this model can be computed very efficiently by computing the discrete Fourier transforms of each patch and filter using the fast Fourier transform (FFT), multiplying them element-wise in the frequency domain, and then taking the inverse discrete Fourier transform of the result. Under varying assumptions, different authors have also proposed locally-uniform models of spatially-variant blur, which take similar forms to Equation (1.21)[START_REF] Nagy | Restoring images degraded by spatially variant blur[END_REF][START_REF] Vio | Multiple image deblurring with spatially variant PSFs[END_REF][START_REF] Tai | Correction of spatially varying image and video motion blur using a hybrid camera[END_REF].

  1.24) where • ⊗ • represents 2D correlation, and the function XCorrMatrix constructs the full cross-correlation matrix from a cross-correlation vector between two signals by replicating and shifting elements. A cross-correlation matrix M for a pair of signals u and v stores their inner-product at all possible translations of both signals. Assuming appropriate boundary conditions, each row of M is simply a shifted version of the cross-correlation u ⊗ v.

Figure 1 . 4

 14 Figure 1.4 Approximating spatially-variant blur by combining uniformly-blurred, overlapping patches. Using the model described in Section 1.5.1, we can efficiently compute approximations to the spatially-variant blur model in Equation (1.12). With a small number of patches (a), the PSF at each pixel is visibly the sum of different blurs from overlapping patches. As more patches are used (b-c), the approximation becomes increasingly close to the exact model (d) -at 12 × 16 patches it is almost indistinguishable.

  blind deblurring results of real camera shakes are shown in figures Figure 1.1 and Figure 1.7.

Figure 1 . 5

 15 Figure1.5 Blind deblurring of real camera shake. A hand-held image with camera shake (a), captured with a shutter speed of 1 second, with the results of blind deblurring using the algorithm of Cho & Lee with (b) a uniform blur model and (c-d) our spatially-varying blur model. The estimated kernels are shown inset in the deblurred results. The result using our blur model in the MAP-ℓ1 algorithm (d) shows 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. Also shown is the result when using our blur model with ℓ2 regularisation on the kernel (c). 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. The rotational blur kernels in (c-d) cover ±0.7 • in θX and θY and ±1.4 • in θZ .

Figure 1 . 6 (Figure 1 . 7

 1617 Figure 1.6 Blind deblurring of real camera shake. The result of blind deblurring on a real camera shake image (a), captured with a shutter speed of 1 second, using the MAP approach of Cho & Lee with both the uniform and non-uniform blur models. Also shown in (b) are some of the local PSFs generated from the blur kernel in (d) at various points in the image. In the blurred image, most of the text on the book cover is too blurred to read. Deblurring the image with the uniform blur model (c) allows some of the text on the cover of the book to be read, however, after deblurring with our non-uniform model (d), all but the smallest text becomes legible. The estimated kernels for the two models are shown inset in the deblurred results. The blur kernel in (d) covers ±0.4 • in θX and θY , and ±0.9 • in θZ . Reproduced from (Whyte et al. 2012) with permission, © Springer 2012.

Figure 1 . 8

 18 Figure 1.8 Blind deblurring failures. A blurred image (a) with the deblurred result (b) and kernel (c) estimated by the MAP-ℓ1 algorithm.This image contains a large amount of noise, and was captured with a camera-phone whose response function is unknown. As a result, the algorithm fails to estimate a good kernel, and instead returns a kernel which is close to a delta function.

http://epaperpress.com/ptlens/

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

Acknowledgements

This work was supported in part by the MSR-INRIA laboratory, the EIT ICT labs, ERC grants VisRec and VideoWorld, and by the Institut Universitaire de France.


Nicolas Sutton-Charani

Sébastien Destercke

Thierry Denoeux

Application of E 2 M decision trees to rubber quality prediction

Keywords: classification, decision trees, rubber quality, hevea, belief functions, algorithm EM

In many applications, data are often imperfect, incomplete or more generally uncertain. This imperfection has to be integrated into the learning process as an information in itself. The E 2 M decision trees is a methodology that provides predictions from uncertain data modelled by belief functions. In this paper, the problem of rubber quality prediction is presented with a belief function modelling of some data uncertainties. Some resulting E 2 M decision trees are presented in order to improve the interpretation of the tree compared to standard decision trees.

Introduction

Learning a classifier from uncertain data necessitates an adequate modelling of this uncertainty, however learning with uncertain data is rarely straightforward. As data uncertainty is of epistemic nature, the standard probabilistic framework is not necessarily the best framework to deal with it. More general frameworks have therefore been proposed [START_REF] Zadeh | Fuzzy sets as a basis for theory of possibility[END_REF][START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF][START_REF] Walley | Statistical reasoning with imprecise probabilities[END_REF] that provide more adequate model for this type of uncertainty. Different classifier learning techniques [START_REF] Kanj | Evidential multi-label classification using the random k-label sets approach[END_REF][START_REF] Sutton-Charani | Classification trees based on belief functions[END_REF][START_REF] Janikow | Fuzzy decision trees: issues and methods[END_REF] using these models have then been developed.

In this paper, our goal is to learn a model from agronomic data. More precisely, we want to predict natural rubber quality from data concerning latex culture and natural rubber maturation. Generally speaking, uncertain measurements and expert assessments are common in agronomy and life science, mainly due to field and economy constraints. They are therefore domains where data uncertainty happens a lot. We retain the belief-function theory [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF][START_REF] Shafer | A mathematical theory of evidence[END_REF], as it is flexible enough to model a large variety of data uncertainties. The chosen classifier is the E 2 M decision tree [START_REF] Sutton-Charani | Learning decision trees from uncertain data with an evidential em approach[END_REF], for it is usually efficient and interpretable (an essential feature for agronomic experts).

After a short overview on the necessary background in Section 2, we detail in Section 3 the application context as well as the uncertainty models we used. We conclude Section 3 by comparing the results of the obtained E 2 M decision trees with classical ones.

Background

We briefly recall the elements needed in the application.

Formalism

As in any classification problem, the aim is to predict a class label Y from a set of attributes (or features) X. The classifier is learnt on a learning dataset LD containing samples of (X,Y ). The classifier is then evaluated by measuring its accuracy on a test dataset T D, comparing the predicted class labels with the real ones.

The attributes X = (X 1 , . . . , X J ) take their values on

Ω X = Ω X 1 ו • •×Ω X J , the class Y on Ω Y = {ω 1 , . . . , ω K }.
That is, K different classes are predicted using J different attributes (either categorical or real-valued).

A precise dataset containing N samples is a set of observations of (X,Y ) and is denoted by

D =    x 1 , y 1 . . . x N , y N    =    x 1 1 , . . . , x J 1 , y 1 . . . x 1 N , . . . , x J N , y N    .
Samples are here assumed to be i.i.d (independant and identically distributed).

Belief-function theory

The theory of belief functions (T BF), also called evidence theory or Dempster-Shafer theory was first presented by Dempster [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF] in a statistical approach. The very basis of the T BF is here presented, with a special focus on the evidential likelihood proposed by Denoeux [START_REF] Denoeux | Maximum likelihood estimation from uncertain data in the belief function framework[END_REF].

Generalities Assume we have an uncertain observation of a variable W defined on a finite space Ω W . We model this observation by a belief mass m W : 2

Ω W → [0, 1] verifying ∑ B∈2 Ω W m W (B) = 1.
We assume here that m W ( / 0) = 0. A focal element A ∈ 2 Ω W is a set such that m W (A) > 0. From this mass, the belief and plausibility functions are defined by:

Bel W (A) = ∑ B⊆A m W (B), Pl W (A) = ∑ B∩A = / 0 m W (B)
Bel W (A) measures the amount of information that implies W ∈ A, and is a measure of certainty, while Pl W (A) measures the amount of information that does not conflict with W ∈ A, and is a measure of plausibility. We naturally have Bel W (A) ≤ Pl W (A) with the two being equal in the specific case of probabilities.

The particular cases of precise data, imprecise data, missing data, probabilities and possibilities can all be modelled by belief functions:

precise data : m W ({w}) = 1 imprecise data : m W (A) = 1 missing data : m W (Ω W ) = 1 (complete ignorance) probabilities : m W (A) > 0 if |A| = 1 consonnant mass functions : m W (A) > 0 and m W (B) > 0 only if A ⊂ B or B ⊂ A
In our classification context, an evidential dataset ED will be of the form

ED = m X,Y =    m X,Y 1 . . . m X,Y N    =    m X 1 1 • • • m X J 1 m Y 1 . . . . . . . . . . . . m X 1 N • • • m X J N m Y N   
where m X,Y i describes the i th sample with its uncertainty.

Evidential Likelihood Assume now we want to fit a parametric model with parameter θ to the data. Likelihood maximisation often provides a good estimator θ of the unknown parameter θ . When data are uncertain, the likelihood can be re-written in the following way:

precise likelihood: L(θ ; w) = P θ (W = w) imprecise likelihood: L(θ ; A) = ∑ w∈A L(θ ; w) (1)
evidential likelihood: L(θ ; m W ) = ∑ A⊆Ω W m W (A i )L(θ ; A i )
As shown in [START_REF] Denoeux | Maximum likelihood estimation from uncertain data in the belief function framework[END_REF], this evidential likelihood can be maximised by the adaptation of the EM algorithm to belief functions: the E 2 M algorithm. This algorithm is quite similar to the EM and is guaranteed to converge towards a local maximum. The main difference is at the Expectation step (E), where the measure used to compute the expectation is the conjunctive combination of P θ and m W .

Decision trees

Decision trees are basic classifiers widely used in many areas such as machine learning, statistics, data mining etc. They are usually built from precise datasets by partitioning the attribute space in a set of leaves, each leave being thus attached to some conditions on the attribute values and to a predicted class.

As each leave is characterized by the proportion of the population "falling" into it as well as the frequencies of the class within this population, a decision tree can be viewed as a multinomial mixture whose parameters are the leaves probabilities (corresponding to the mixture coefficients) and the class probabilities inside the leaves (multinomial).

Example 1. Let us consider a data set with N = 12 data items, J = 2 attributes and K = 3 classes. Spaces are described as follows: Figures 1 and 2 respectively represent a possible decision tree and its corresponding partition of the attribute space with the learning dataset represented in it. To leaf t 2 are associated the values

Ω X 1 = [1, 10], Ω X 2 = [13, 150], Ω Y = {a, b, c}. X 1 y = a X 1 < 3 X 2 y = b X 2 < 45 y = c X 2 ≥ 45 X 1 ≥ 3 t 1 t 2 t 3
A 1 2 × A 2 2 = [3, 10] × [13, 45[
, its estimated probability is 3 /12 and the class estimated probabilities are ( 1 /3, 2 /3, 0), and the prediction for any test data item in t 2 (for example, x 1 = 6 and x 2 = 30) will be y = b.

The growing of a decision tree is recursive, and aims at separating classes inside leaves. This separation is usually measured by an impurity measure i computed from the leaves probabilities and the class distributions inside leaves. From a root node containing all the samples of the learning dataset, all possible splits (on all attributes and for all values) are explored, and the one with the highest purity gain (i.e., highest impurity reduction) is chosen. Then, for each newly created leave and the sample "falling" into it, the process is iterated until some stopping criteria is reached.

The E 2 M decision trees

To learn decision trees from uncertain data (where potentially both attributes and classes can be uncertain), we proposed [START_REF] Sutton-Charani | Learning decision trees from uncertain data with an evidential em approach[END_REF] to learn the multinomial mixture of the tree through the E 2 M algorithm. Périnel [START_REF] Périnel | Construire un arbre de discrimination binaire à partir de données imprécises[END_REF] proposed a similar idea, yet only dealt with uncertain attributes and probabilistic data. We refer to [START_REF] Sutton-Charani | Learning decision trees from uncertain data with an evidential em approach[END_REF] for technical details about the learning process of E 2 M decision trees, as well as for some experiments on benchmark data comparing E 2 M and CART decision trees and showing the potential interest of E 2 M decision trees in terms of accuracy in noisy environments.

Application: the rubber quality problem

We first present the application and two uncertainty models we used on the data, before showing some interesting results. Note that the uncertainty models could be re-used in similar situations.

Problem description

Natural rubber is the result of a milky fluid transformation: the latex that is extracted from the hevea tree. Compared to synthetic rubber, natural rubber presents some unique physical properties but suffers from a pretty high variability in terms of quality, that experts are not fully able to explain and to control.

In order to better control this quality, one of the biggest rubber company has put some effort to study this variability in one of its plantations by collecting cultural data. This plantation is located in a part of Brazil where the climate has the particularity to be more variable than in the rest of Brazil (the natural origin of the hevea tree is in the Brazilian forests). However, due to the size of the plantation (approximatively 70 hectars) and various factors (e.g., untracked delay of collection/process due to weather conditions, tanks mixing productions of many parcels), some variables are subject to high uncertainties.

Data are constituted of many variables summarised in Table 1 (no Unit means a dimensionless variable). Meteorological data may influence the latex during three different periods: the latex fabrication by the tree (one week before tapping), the tapping day during which latex is collected, and the latex maturation in tanks (five days). For the temperature and the relative humidity, the minimum, median and maximum values are computed for each day.

The data set contains 3053 examples described by 106 attributes. The quality is measured by the P 30 index which is an elasticity index. In order to use the E 2 M decision trees methodology, the P 30 was discretised into 5 equiprobable classes. This discretisation is presented in Table 2.

Data uncertainty modelling

Two types of uncertainty were modelled in this application: one relative to the rainfall, and one due to parcel mixture in tanks.

Rainfall uncertainty

The rain is a phenomenon that is geographically very variable, especially in tropical areas. In the plantation, all the rain data come from a single meteorological station located inside the plantation. Since the plantation area is very large, it is sensible to make the hypothesis that the farther is located a parcel in the plantation from the meteorological station, the more uncertain is its rainfall data. This uncertainty is non-probabilistic and progressive, so we chose to model it with a consonant mass function. Moreover, as more rainfall implies more uncertainty, it is logical to assume that imprecision of focal elements increases multiplicatively (i.e. more rainfall measured by the station implies wider focal elements).

To keep the complexity of the uncertainty reasonable, we limited the mass to five focal elements of the form: [w(1δ ), w(1 + δ )] where w is the original precise rainfall data and where δ ∈ ∆ = {0, 0.25, 0.5, 0.75, 1}. The proposed model is easy to expand to more than five focal elements and can therefore accommodate various levels of complexity (depending on the available computational power and on possible time constraints). We define a function g : [0, d max ] × ∆ → [0, 1] (d max being the maximal distance between a location of interest and the measurement station) such that the rainfall w of a parcel located at a distance d from the meteorological is characterized by g

(d, δ ) = m W ([w(1 -δ ), w(1 + δ )]) for all δ ∈ ∆ .
We distinguish two types of focal elements, the most precise ones (δ < 0.5), and the most imprecise ones (δ ≥ 0.5). In the first case (precise ones), we assume that the farther was the parcel from the meteorological station (d increasing), the smaller had to be the masses assigned to those focal elements. In the second, we want to assign bigger masses to the farther parcels. Such assumptions can be translated in the following constraints:

δ < 0.5 → ∂ g ∂ d < 0 δ ≥ 0.5 → ∂ g ∂ d > 0 (2)
that merely translate the assumption into derivative behaviors. As function g is used to define mass functions, we must add the following constraints

∀(d, δ ) ∈ [0, d max ] × ∆ , g(d, δ ) ≥ 0 (3) ∀d ∈ [0, d max ], ∑ δ ∈∆ g(d, δ ) = 1 (4)
that simply ensure one that the obtained mass functions will be well-defined, i.e. that it will be positive (constraint ( 3)) and will sum up to one (constraint ( 4)).

One simple solution of this problem is to use two linear functions, one increasing for the most precise focal elements (δ < 0.5) and one decreasing for the most imprecise ones (δ ≥ 0.5), and to use a convex sum of those two functions. We obtain:

g(d, δ ) = δ ( 2d 5d max ) + (1 -δ )( 2 5 - 2d 5d max ) (5)
Example 2. Consider three rainfall measurements w 1 = 0, w 2 = 10 and w 3 = 30 from the station, and for each of these measurements some corresponding parcels of interest respectively distant of 20km, 50km and 2km from the station. Assuming that d max = 80, we obtain m w 1 ({0}) = 1, given the multiplicative uncertainty, and

           m w 2 ({10}) = 0.15 m w 2 ([7.5, 12.5]) = 0.175 m w 2 ([5, 15]) = 0.200 m w 2 ([2.5, 17.5]) = 0.225 m w 2 ([0, 20]) = 0.250            m w 3 ({30}) = 0.39 m w 3 ([22.5, 37.5]) = 0.295 m w 3 ([15, 45]) = 0.2 m w 3 ([7.5, 52.5]) = 0.105 m w 3 ([0, 60]) = 0.01.
The absence of rain is thus considered certain (m w 1 ) whereas positive rainfall data masses are concentrated on imprecise focal elements when coming from distant parcels (m w 2 ) and on more precise ones when coming from parcels close from the meteorological station (m w 3 ).

Parcelar mixtures uncertainty During the harvest, the latex coming from many parcels is usually mixed in some tank. All the parcel variables (i.e., most agronomical variables of Table 1) we have are therefore subject to uncertainty as the amount of latex coming from each parcel in tanks is not tracked. During a pre-treatment of the data, we therefore split all those parcel variables into rough and uncertain proportions (due to latex production high variability) computed from the weight of latex produced annually by each parcel (shifting from 18 original attributes to 106, with all the split ones being in range [0, 1]).

For example, if 25% of a tank content comes from clone A parcels and 75% from clone B parcels, the actual amount of clones A and B latex in the tank may be quite different, as each clone has variable production capacities (that may depend differently on the weather, soil conditions, etc.). We model this uncertainty such that the more balanced are the parcel proportions in a tank, the more uncertain become those proportions: proportions of a tank with latex from only one pure parcel should remain certain, while proportions of a tank with equal amounts of latex coming from different parcels should be maximally uncertain.

To do that we used simple intervals around computed crisp proportions, with the interval width increasing as proportions become uniform. To measure this uniformity we simply used the Shannon entropy denoted ent computed on the set of parcel proportions of the tanks. The obtained model is the following: for each parcel variable X j having r modalities with the positive proportions {p 1 , . . . , p r },

       m j 1 ([max(p 1 -ent(p 1 ,...,p r ) r , 0), min(p 1 + ent(p 1 ,...,p r ) r , 1)]) = 1 . . .

Experiments

In order to see the consequences of integrating data uncertainty the way we described in Section 3.2, we perform some experiments comparing standard CART trees and E 2 M trees on the original precise dataset and its corresponding evidential dataset obtained with our uncertainty models, respectively.

For both methodologies, the stopping criteria is a maximum of 5 leaves (to preserve a high interpretability) and a relative purity gain of 0.05.The error rates were computed as the proportion of misclassified examples in the test dataset. Their means were computed from ten 3-fold cross validations. It is noticeable that we used standard (precise error rates) even for the E 2 M decision trees for comparison purposes. Given the small tree size, no pruning is done. As shown in Table 3, the accuracies of the two methodologies are quite similar, even if the E 2 M accuracy is slightly better. Let us now shift to the main interesting part for the experts (and hence for the application goal): the interpretation of the trees.

Interpretation

In a pure knowledge discovery concern, we learnt CART and E 2 M decision trees on the whole dataset in order to compare the informations they provide about the natural rubber quality explanation. Within such a goal, it should also be noted that integrating data uncertainty makes the result somehow more faithful to our actual information (and therefore more reliable for the experts). The learning parameters were exactly the same as in Section 3.3. In Figures 3 and4 that show the learning results, UR min matur is the minimum of relative humidity during the maturation of the latex in buckets, Temp med (resp. min) tapping is the medium (resp. minimum) temperature during the tapping day, # of tapped tree per Ha is the number of tapped trees per hectare.

As shown in those figures, integrating data uncertainty makes the number of tapped trees per hectare appear in the tree. Not only the E 2 M decision tree suggests a high relation between quality and productivity, but it also provides a density acceptable bound inside the plantation. According to the experts, this issue points out the potential need to investigate further this path. We can also notice that this data uncertainty lessen the role of clone PR255 proportion in the rubber quality explanation.

Discussion

In the rubber application, the prediction results from CART and E 2 M decision trees do not seem very different in term of accuracy (even if E 2 M has a slight advantage), but the interpretation of the trees can be quite different; this latter interpretation being as important as accuracy in many applications. Indeed, we observe that some variables may play a significant role in rubber quality explanation once their uncertainty is modelled (here by belief functions). Modelling uncertainty also provides more trustful results for the experts.

Perspectives concerning the application include the modelling of additional uncertainties (in particular with respect to the tapping day) as well as additional studies involving more variables (e.g., tree diseases, soil conditions) or other prediction goal (e.g., quantity of production, which is also important from an economical viewpoint).

Finally, as the uncertainty models we have introduced may be useful in other areas (particularly the distance-based model of rainfall), we plan to generalize them and study in more details their properties.

Fig. 1 : 3 Fig. 2 :

 132 Fig. 1: Decision tree illustration

, 1 )]) = 1 ( 6 )

 116 m j r ([max(p r -ent(p 1 ,...,p r ) r , 0), min(p r + ent(p 1 ,...,p r ) r with m j i modelling the uncertainty about the j i th proportion. Example 3. Let us consider a tank containing 75% of clone A and 25% of clone B. The entropy on those proportions is equal to 0.8113. The obtained masses are therefore m clone A ([34.43%, 100%]) = 1 m clone B ([0%, 65.57%]) = 1

Fig. 4 :

 4 Fig. 3: CART decision tree

Table 2 :

 2 P 30 discretisation

	variable		type	category	unit
	period	categorical agronomical
	season	categorical climatic
	weight	numerical agronomical	kg
	X (latitude)	numerical geographical	km
	Y (longitude)	numerical geographical	km
	clone	categorical agronomical
	panel	categorical agronomical
	tapping system	categorical agronomical
	surface	numerical agronomical	hectar
	planting year	numerical agronomical	year
	first tapping year	numerical agronomical	year
	tapping age	numerical agronomical	year
	annual number of tapped trees numerical agronomical
	tapped tree per hectar	numerical agronomical
	temperature	numerical climatic celcius degrees
	relative humidity	numerical climatic
	sun hours	numerical climatic	hours
	rainfall	numerical climatic	mm
	P30	numerical agronomical
	Table 1: variables characteristics
	class labels P 30 range
	very bad [1.87 ; 14.7[
	bad		[14.7 ; 21.6[
	medium [21.6 ; 27.4[
	good	[27.4 ; 32.9[
	very good [32.9 ; 49.5[

Table 3 :

 3 Results

	methodology mean error rates 95% confidence interval
	CART	0.6607	[0.6314 ; 0.6900]
	E 2 M	0.6560	[0.6266 ; 0.6854]

Acknowledgments. This work has been done thanks to CIRAD fundings and LASB laboratory support. We also wish to greatly thank Dr. Patrice Loisel for his great help on the rainfall uncertainty model.