BIT Numer Math DOI 10.1007/s10543-013-0464-y

Rotational image deblurring with sparse matrices Per Christian Hansen · James G. Nagy · Konstantinos Tigkos

Received: 22 July 2013 / Accepted: 4 December 2013 © Springer Science+Business Media Dordrecht 2013

Abstract We describe iterative deblurring algorithms that can handle blur caused by a rotation along an arbitrary axis (including the common case of pure rotation). Our algorithms use a sparse-matrix representation of the blurring operation, which allows us to easily handle several different boundary conditions. We also include robust stopping rules for the iterations. The performance of our algorithms is illustrated with examples. Keywords Image deblurring · Iterative algorithms · Stopping rules · Sparse matrices · Boundary conditions Mathematics Subject Classification (2010)

65F22 · 65F10

Communicated by Rosemary Renaut. This work was supported by Grant No. 274-07-0065 from the Danish Research Council for Technology and Production Sciences, Grant No. DMS-1115627 from the US National Science Foundation, and Grant No. AF9550-12-1-0084 from the US Air Force Office of Scientific Research. P. C. Hansen (B) Department of Applied Mathematics and Computer Science, Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark e-mail: [email protected] J. G. Nagy Department of Mathematics and Computer Science, Emory University, Atlanta, GA 30322, USA e-mail: [email protected] K. Tigkos Fraunhofer Institute for Integrated Circuits, Am Wolfsmantel 33, 91058 Erlangen, Germany e-mail: [email protected]

123

P. C. Hansen et al.

1 Introduction Image deblurring is an important image restoration problem in its own right, as well as a component in navigation systems, robotics, etc, and a large number of specialized algorithms have been developed for this task. A common cause of blurring is motion of either the recording device or the object. Examples include movement of the heart during cardiac imaging processes, movement of a low flying aircraft in aerial imaging, and people or vehicle movement in surveillance imaging. Linear motion with constant speed and direction is well understood, but relatively little work has been done in the development of mathematical models and efficient algorithms for restoration of images degraded by nonlinear and nonuniform motion blur where the blurring is spatially variant [4,10,11,19,28,34,40,41]. This work focuses on a particular type of motion blur where the object or scene is rotated during the recording of the image, leading to spatially variant blur. An important special case is when the rotation axis of the scene points towards the camera, and this case has already been studied in several papers, e.g., [10,32,35]. In this case, if the rotation also has constant velocity, then it is possible to transform the blurring to polar coordinates, which results in a shift-invariant operator [10]. Although this is mathematically convenient, it is not possible to obtain a purely shift-invariant operator for more complicated, nonuniform motion blurs. In addition, the transformation between rectangular and polar coordinates can introduce severe interpolation errors in the reconstruction. We are concerned with algorithms that avoid such transformations. One approach would be to use general algorithms for spatially variant blur [22,40] that use a sectioning approach that assumes the motion blur is approximately spatially invariant in various sections of the image. However, this only works well if the motion blur changes smoothly (and slightly) between adjacent sections. The same is true for the approach in [38]. A completely general treatment of motion blur, including rotational blur around an arbitrary axis, is given in [39] which uses an iterative reconstruction algorithm. Here the “forward computation” of blurring is achieved by summing a sequence of rotated and/or translated images. This matrix-free approach is storage efficient, but can be slow because an accurate representation of the blur requires many small steps to represent the motion accurately. Our algorithms are also based on iterative algorithms, but instead of the matrix-free approach we form a sparse matrix that represents the spatially variant rotational blur. The matrix needs only be constructed once, and we can easily incorporate suitable boundary conditions into the matrix. The advantages are fast execution, at the expense of storage requirements for the sparse matrix, and flexible treatment of boundary conditions. To make our algorithms attractive we consider robust stopping rules, we include nonnegativity and box constraints on the pixel values, and we incorporate boundary conditions in order to reduce artifacts from the edges. Few of the algorithms presented in the literature cover all these important aspects. Sections 2 and 3 describe the underlying geometry and how we model the forward problem and create the sparse matrix, including a discussion of boundary conditions. In Sect. 4 we review the iterative reconstruction methods, and in Sect. 5 we demonstrate

123

Rotational image deblurring

their performance for various rotational deblurring problems, including a discussion of the influence of boundary conditions. Section 6 discusses the stopping rules and our implementation of these rules. Finally, in Sect. 7 we illustrate the performance of our algorithms on a realistic test problem. 2 Setting the stage Our geometry setup is as follows: the origin of the coordinate system is the camera’s focal point. The object or scene being recorded is assumed to be a plane, and this plane is projected onto the image plane during the rotation. The rotation is always around a single fixed rotation axis that goes through the focal point. We consider the following three cases: – Pure rotation. The rotation axis points toward the camera (i.e., the scene rotates around its center) and is orthogonal to the image plane; hence we have simple rotation around a single point in the image. – Pure tilt. The rotation axis is orthogonal to the direction of the scene (i.e., the scene moves, or spins, on a circle round the focal point at a fixed distance). – Tilted rotation. Similar to the pure rotation, except that the rotation axis is not orthogonal to the image plane. The first case of pure rotation has been dealt with in several papers, e.g., [32]. To our knowledge the latter two cases have not been covered in the literature, except in the general treatment of Tai et al. [39]. We assume throughout our paper that the directionof-rotation as well as the angle of rotation are known – it is outside the scope of this work to determine these parameters, and instead we refer to [5,7,38]. Previous work for pure-rotation deblurring has mainly been based on a transformation to polar coordinates, where the motion blur is linear and for which compensation is easy. This approach, however, has two main drawbacks: (1) interpolation errors are introduced in the transformations to and from polar coordinates, and (2) the transformation to polar coordinates cannot handle a general direction of rotation. Hence we need a more general formulation of the deblurring problem. Our approach is to set up a constrained linear least squares problem associated with the deblurring problem:   (2.1) min  A f − g 2 , f∈C

in which the vectors f and g represent the image to be reconstructed and the noisy blurred image, respectively; both are assumed to be N × N images, so f, g ∈ Rn with n = N 2 . The n × n matrix A represents the spatially variant blurring; this matrix is sparse, but without any further structure. The set C is either Rn (no constraint), the first orthant Rn+ (nonnegativity constraints), or the n-dimensional box [0, L]n , where L is the dynamical range (box constraints on all pixels). We consider two iterative methods for solving the constrained least squares problem: the projected Landweber (PL) method with box constraints [6,27], and the modified residual norm steepest descent (MRNSD) algorithm [3,24] with nonnegativity constraints. These methods are chosen for their robustness and their robust convergence

123

P. C. Hansen et al.

behavior [23], and the nonnegativity/box constraints are known to have a stabilizing effect [2]. We remark that the well-known Richardson–Lucy algorithm is not used in our work because MRNSD can be used for applications where Poisson statistical models are necessary, and MRNSD generally converges much faster than Richardson– Lucy; see [3]. As stopping rules we consider generalized cross validation (GCV) and the normalized cumulative periodogram (NCP) criterion, both based on statistical considerations about the residual vector [14]. We did not implement the predictive risk estimator from [1] because its performance is similar to that of GCV. 3 Modeling of the forward problem We model the forward problem as g = Af + ε ,

(3.1)

where f represents the true, unknown image, g regresents the observed blurred and noisy image, and A is a large sparse matrix that models the blurring operation. It is important to note that the blur is highly spatially variant, so it is not possible to use a convolution model, and the matrix A does not have the associated convenient Toeplitz structure. As we describe below, it is possible to construct a sparse matrix A that models the blurring operation. 3.1 Towards the discretized problem The forward model is a mathematical description of the blurring process, here the blurring due to the rotation. Assume that we perform a rotation α with a fixed rotation axis; this is a projective transformation in which the original image coordinates ξ = (x, y, z) and the transformed image coordinates ξ  (α) are related by ξ  (α) = H (α) ξ ,

(3.2)

where the 3 × 3 matrix H (α) represents a homography [17] that can describe such planar projective transformations. This model is well suited for planar scenes without out-of-plane components (which would be distorted by the projective transformation). Thus we can compute the rotated image by using the above coordinate transformation. The blurring through an arbitrary angle is described by an integral of the rotated images. For generality, let the rotation angle α(t) depend on a parameter t (which could be time), going from α(0) = 0 to α(1) = θ . If F is a function that represents the sharp image then we can write the blurred image as a function G given by 1 G(ξ ) =

  F H (α(t)) ξ dt,

(3.3)

0

where H (α(t)) is the above-mentioned homography matrix that describes the image coordinate transformation for rotation angle α(t).

123

Rotational image deblurring

Fig. 1 The principle behind creating the blurred image as a sum of images that undergo a sequence of small rotations

We model (3.3) as a finite sum over rotated images through M − 1 small rotations α0 = 0, α1 , α2 , . . ., α M−1 = θ , see Fig. 1, and we obtain: G(ξ ) ≈

M−1  1   F H (αi ) ξ . M

(3.4)

i=0

Note that in the special case αi = i·θ/M, i = 0, 1, . . . , M − 1 we have H (αi ) = H (θ/M)i . In practice, we work with digital images F and G instead of the continuous images F and G. Utilizing (3.4) we can create a black-box operator which, given a sharp image F and the geometric information about the blurring, computes an approximation to the blurred image G via M − 1 interpolations in the sharp image. With the notation f = vec(F) and g = vec(G) (where “vec” denotes stacking the matrix columns into a long vector), this operation represents the forward computation g = A f in the notation of (2.1). We emphasize that in this approach, the matrix A does not need to be explicitly formed, so it has the advantage that the memory requirement is very low—but for each forward computation in the iterative reconstruction algorithm M − 1 interpolations must be performed for each pixel, which can greatly slow down the reconstruction algorithm. 3.2 Setting up the sparse matrix As discussed in the previous section, the basic idea is to assume the motion-blurred image is the scaled sum of images at incremental times, or positions, during acquisition: G=

M−1 1  F( ξ i ). M

(3.5)

i=0

– Each F( ξ i ) represents the ith snapshot of the object, – ξ i , i = 0, 1, 2, . . . , M − 1 are position vectors of the ith transformed images, where ξ i = H (αi )ξ 0 and θi is the angle at which the original image is rotated to obtain the transformed image, and – ξ 0 is the position vector corresponding to the original image. To represent the digital image, we must determine approximations of F at the pixel centers; we do this using bilinear interpolation. Specifically, each pixel center is approx-

123

P. C. Hansen et al. Fig. 2 Illustration of bilinear interpolation to approximate FPC using the known values FNW , FNE , FSW , and FSE . The approximation is given by (3.6), using x and y in this figure

imated by a weighted average of the four transformed values that surround it, and the weights are proportional to the distance from the pixel center. For example, consider the case of a pure rotation. Without loss of generality we can assume that the distance between pixel centers is one, and we consider the illustration in Fig. 2 where FPC is a particular pixel center with surrounding transformed values FNW , FNE , FSW and FSE . Then bilinear interpolation approximates FPC as FPC ≈ (1 − x)(1 − y) FSW + xy FNE +(1 − x)y FNW + x(1 − y) FSE .

(3.6)

While the weights in (3.6) are determined under the assumption that each pixel is a square with sides of length one, it is not difficult to modify the weights for pixels of arbitrary size. Since f = vec(F) we can write vec(F( ξ i )) = Ai f , where Ai is a sparse matrix whose rows contain the interpolation weights, (1 − x)(1 − y), (1 − x)y, x(1 − y), xy, for the ith transformed image. We emphasize that by using a sparse data format (e.g., compressed row storage [8]) to represent Ai , we need only keep track of the nonzero entries and their locations in the matrix Ai . As previously discussed, we assume the observed motion blurred image is the scaled sum of images at incremental times during the acquisition. From (3.5) we thus obtain g=

M−1 M−1 1  1  vec(F( ξ )) = Ai f = A f, M M i=0

i=0

where the matrix that models the motion blur is A=

M−1 1  Ai . M i=0

123

(3.7)

Rotational image deblurring

Incorporating additive noise in the model, we arrive at the linear inverse problem given in (3.1). While it may not be possible to precisely know the relative position of the object and camera, there are applications in which this information can be either approximated [40], or measured to high accuracy [11,28]. Sparseness of A decreases (A becomes more dense) as more motion information is used, or if a higher-order interpolation method (e.g., cubic) is used to construct each Ai . For large rotational blurs, there is a tradeoff between accurately modeling the motion blur and computational cost. For example, in [28], the authors apply motion correction to a brain imaging application, with an approach similar to what is considered in this paper. Although cubic interpolation was attempted, they eventually used constant (nearest neighbor) interpolation, stating that the “insignificant” improvement in the restored images, when using cubic interpolation, did not justify the substantial increase in computational cost. We feel that bilinear interpolation is a good tradeoff, giving more accuracy than constant interpolation, and yet the construction of A is generally not too time consuming, and the sparsity can be efficiently exploited when performing matrix-vector products in the iterative image deblurring methods.

3.3 Incorporating boundary conditions In most image deblurring problems, it is essential to somehow account for information that is close to the image borders but outside of the captured image. Since this information is usually not available, one uses instead information from inside the image, a technique typically referred to as incorporating boundary conditions (BC) in the reconstruction; see, for example, [16,30,36]. Some of these techniques can be implemented very efficiently with fast transforms (e.g., FFT or DCT) if the blur is spatially invariant, but may require additional assumptions on the point spread function (PSF), such as strong symmetry or very small support of the PSF. However, for iterative methods all of these techniques, as well as many more sophisticated approaches (e.g., [12]), can be used by simply padding the images with appropriate pixels values. In spatially invariant deblurring, boundary condition mismatch often creates ringing artifacts near the image borders. As we shall illustrate below, the influence of the boundary conditions in connection with rotational blur, which is highly spatially variant, has a different and even more severe effect on the reconstruction quality. During the acquisition of an image subject to rotational motion, information that was not part of the initial frame (before the motion took place) is recorded—and at the same time information that was part of the initial frame is lost. Thus, as the rotation angle increases, we can say that we have more “data loss” and also more information from outside the original scene being recorded in the blurred image. The aspect of “data loss” is reflected in the forward computation (3.4): for each step of the spatial transformation, a part of the transformed image moves out of the frame, and another part enters the image. A good reconstruction model must reflect this, in order to reduce artifacts at the corners of the reconstructed image. This is achieved by augmenting the forward model with specific assumptions about information outside the borders of the frame prior to the motion, i.e., by specifying the boundary conditions to be included in the forward model. These boundary conditions are defined and

123

P. C. Hansen et al.

integrated in an intuitive way and, as seen in the experimental results, this approach gives good results. Specifically, the boundary conditions come into play by “filling in” the missing information in each step of the transformation, thus trying to account for the “unknown” information that was not part of the initial image (prior to the rotation). We consider four possible boundary conditions: – Zero: the scene outside the initial image is black. – Replicate: each pixel at the border of the initial image is infinitely replicated outside the image. – Periodic: the scene outside the initial image is a periodic replication of the image. – Reflective: the scene outside the initial image is a reflection of the image along the image border. The pixel values corresponding to coordinates ξ i outside the initial image domain are thus determined by the imposed boundary conditions. Concerning the sparse matrix approximation, not every row of the Ai matrices will necessarily contain four interpolation weights. Instead, the values and column indices of the interpolation weights for rows that correspond to pixels mapped outside the initial image domain should be determined by the boundary conditions. 4 Iterative reconstruction methods The linear system in (3.1) is an example of an ill-posed inverse problem [14]. Regularization is needed to suppress noise amplification in the reconstructed image, and because A is a very large, sparse matrix, it is essential to use iterative methods. Here we consider two algorithms: the projected Landweber method [27] and the MRNSD method [24]. Both methods have iterates of the form f (k+1) = f (k) + τk d (k) ,

(4.1)

where k denotes the iteration index, d (k) is the step direction vector, and the scalar τk is the step length parameter. 4.1 Projected Landweber (PL) method The Landweber method uses the step direction vector   d (k) = A T g − A f (k)

(4.2)

which is the steepest descent direction for the associated least squares minimization problem. The step length parameter τk = τ remains constant for each iteration, and to ensure convergence it must satisfy 0 < τ < 2/A22 , where A2 is the 2-norm of A. The regularization properties of the Landweber method are well understood [14]. While its convergence may be slow, it has the advantage that it is very easy to introduce

123

Rotational image deblurring

box constraints by simply projecting the approximate solution at each iteration:    f (k+1) = P f (k) + τ A T g − A f (k) . Here P represents the orthogonal projection corresponding to the box constraints: (k+1)

f min ≤ f i

≤ f max , i = 1, 2, . . . n,

where f min and f max are, respectively, lower and upper bounds on the pixel values of the image. For properties of the projected Landweber method, including preconditioning to accelerate convergence, see [6,27]. Properties of the more general projected steepest (gradient) descent method can be found in one of many references on numerical optimization, such as [21,25], and a discussion in the context of image restoration can be found in [42]. 4.2 Modified residual norm steepest descent (MRNSD) The MRNSD algorithm was developed in order to provide a method with better convergence behavior. It uses the step direction vector   d (k) = Dk A T g − A f (k) , where Dk is a diagonal matrix with diagonal entries given by f (k) . The step length τk is chosen to ensure that f (k+1) has nonnegative entries. The derivation of MRNSD can be found in [3], and further details of the method can be found in [20,24]. We remark that if Dk = I (the identity matrix) and τk = τ (a constant), then MRNSD is simply the Landweber method. The presence of Dk acts like a preconditioner, and is related to scaled projected gradient descent methods [21], where Dk is used as a crude approximation to the reduced Hessian. MRNSD is an example of a covariance preconditioned iterative method [3], which incorporates statistical information in terms of an approximation of the noise covariance matrix. In fact, as is shown in [3], the well known Richardson–Lucy (or EM) algorithm is another example in this class of iterative methods, but the Richardson– Lucy method uses a different approximation of the noise covariance matrix, as well as a constant step size at each iteration, which results in much slower convergence than MRNSD. Although MRNSD enforces a nonnegative constraint, the general method does not allow for enforcing an upper bound constraint. In our numerical experiments, if we want to enforce box constraints (that is, both lower and upper bounds) on the pixel values, we apply the upper bound constraint after the final solution is computed. This simple approach does not solve the actual box-constrained optimization problem (see Fig. 2.1 in [9] for an illustation); but the solutions are improved, and it is outside the scope of this work to incorporate box constraints in MRNSD.

123

P. C. Hansen et al.

4.3 Semi-convergence Both methods considered here exhibit semi-convergence behavior, where the early iterations produce iteration vectors that converge towards the desired noise-free solution, while later iterations produce noisy vectors that diverge from from this solution. It is easiest to see this in the basic Landweber iteration without projection. n Suppose, u i σi viT without loss of generality, that the initial guess is f (0) = 0, and A = i=1 is the singular value decomposition (SVD) of A. Then it can be shown [14] that each iteration of the Landweber iteration can be written as f (k) =

n  i=1

(k)

where φi

T (k) u i g

φi

σi

vi ,

(4.3)

are filter factors for the kth iteration, given by (k)

φi

= 1 − (1 − τ σi2 )k , i = 1, 2, . . . , n .

(4.4)

For the largest singular values we have φi(k) ≈ 1, while for the smallest singular values φi(k) ≈ kτ σi2 , i.e., they are proportional to σi2 and therefore damp the corresponding SVD components in the sum (4.3). As the number of iterations k increases, more filter factors become closer to 1, and thus more SVD components are included in f (k) (and it approaches the noise-free solution); but when too many SVD components are included the reconstructions start to become noisy. This semi-convergence behavior acts as a regularization process, where k plays the role of the regularization parameter. We can see from (4.4) that it may take many iterations to reconstruct SVD components corresponding to intermediate singular values. For example, if σ1 = τ = 1 then for σi = 0.1 we obtain the corresponding filter (k) factors φi = 1 − (1 − τ σi )k = 1 − (0.99)k , and thus k needs to be very large to (k) obtain φi ≈ 1. With a faster converging method such as MRNSD, more SVD components are reconstructed at the early iterations. This also means that the noise begins to corrupt solutions sooner, and it is therefore essential to have stopping rules that determine when the noise begins to corrupt the solution. We return to stopping rules in Sect. 6. 5 Comparison of performance In this section we present and discuss some experimental results that illustrate the overall performance of our algorithms. All test problems in this section are synthetically generated by using the same forward model that is used in the reconstruction (this practice is known as “inverse crime” and is only used for studying and understanding the properties of the reconstruction algorithms). The exact test image can be seen in Fig. 5. Given the noise-free blurred image A f exact , we add noise to obtain the noisy blurred image g. Two types of noise commonly present in images, Gaussian and Poisson, were

123

Rotational image deblurring

Rel. error ρ

0.2 0.15 0.1 0.05 0

200

400

600

800 1000

PL MRNSD

0.29 0.28 0.27 0.26 0.25 0.24 0.23 0.22 0.21 0.2 0

50 100 150 200 250 300 350 400

2% Poisson noise Rel. error ρ

PL MRNSD

0.25

Rel. error ρ

2% Gaussian noise

No noise

0.3

PL MRNSD

0.29 0.28 0.27 0.26 0.25 0.24 0.23 0.22 0.21 0.2 0

Iteration k

Iteration k

50 100 150 200 250 300 350 400

Iteration k

Fig. 3 Error histories for “inverse crime” experiments with no noise, as well as with 2 % Gaussian and 2 % Poisson noise. The semi-convergence of both methods is evident in the case of noise. Note the different scalings of the axes

used in our experiments. The noise level in all experiments of this section is 2 %, i.e., g − A f exact 2 /A f exact 2 = 0.02. To evaluate the reconstruction algorithms we use two quality measures to compare the reconstructed image Frec with the exact one Fex . The relative RMS error is defined as (5.1) ρ = Fex − Frec F /Fex F , where FF =  f 2 is the RMS of the pixel values; this is a simple measure that gives the average relative difference of two images. A more sophisticated measure that tries to reflect how the human visual system perceives structure is the Structural SIMilarity (SSIM) index [43], which separates luminance and contrast discrepancies from structural ones, enabling one to penalize each of them. The standard expression is: S=

(2μex μrec + C1 )(2σex,rec + C2 ) , 2 + σ2 + C ) (μ2ex + μ2rec + C1 )(σex 2 rec

(5.2)

where μex and μrec are the mean intensities of the two images, σex and σrec their 2 is their covariance. We use the default values C1 = standard deviations, and σex,rec 2 (0.01L) and C2 = (0.03L)2 where L is the dynamical range of the pixels. In practice, it is advised to compute the SSIM index as a sum of local averages across the entire image; the implementation used in this work is the one accompanying [43]. 5.1 Reconstruction results Figure 3 (left) shows the error histories, ρ versus k, of both reconstruction algorithms, for a 20◦ pure rotation with no noise and with zero BC. It is evident that both algorithms converge to the exact solution. Figure 3 (middle and right) also shows the error histories for the same test problem but with Gaussian and Poisson noise, respectively. Now the iterates approach the exact solution in the early iterations until a minimum error is reached, after which the iterates start to diverge as the reconstruction gets dominated by noise. Figure 4 shows the optimal reconstructions for Gaussian noise, obtained at the point of semi-convergence, i.e., at the minimum of the error history. The first row shows a

123

P. C. Hansen et al.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Fig. 4 Reconstructions for “inverse crime” experiments with 2 % Gaussian noise. Each picture shows the reconstruction at the point of semi-convergence, together with the corresponding values of the relative error ρ and the SSIM index S. a Pure rotation, b PL, ρ = 0.21, S = 0.72; c MRNSD, ρ = 0.20, S = 0.75; d pure tilt, e PL, ρ = 0.23, S = 0.65; f MRNSD, ρ = 0.22, S = 0.68; g tilted rotation, h PL, ρ = 0.22, S = 0.71; i MRNSD, S = 0.22, S = 0.73

case of 20◦ pure rotation, the second row a case of 20◦ pure tilt, and the third row a case of tilted rotation, with a 20◦ rotation where the planar scene has been tilted by 20◦ before recording the image. Similar results were Table 1.

123

Rotational image deblurring Table 1 Reconstruction results for “inverse crime” experiments with 2 % Poisson noise The table shows the minimum relative error ρ and SSIM index S at the point of semi-convergence

Rotation

PL ρ

MRNSD ρ

S

S

Pure rotation

0.2149

0.7328

0.2117

0.7448

Pure tilt

0.2320

0.6635

0.2273

0.6761

Tilted rotation

0.2218

0.7234

0.2240

0.7300

(a)

(c)

(b)

(d)

(e)

Fig. 5 Illustration of the role of boundary conditions. The reflective BC in d match the sharp padded image, while this is not true for the zero BC in e. a Sharp padded image, b blurred padded image, c cropped, d rec., reflexive BC, e rec., zero BC

5.2 Influence of boundary conditions An “inverse crime” experiment was conducted in order to illustrate the importance of boundary conditions and the resulting effect of boundary condition mismatch for the particular problem. We use 20◦ pure rotation blur and no noise. The test image, marked by the white box in Fig. 5a, was padded to simulate a case of reflective BC, i.e., assuming that the scene around the captured image is a reflection of the image. This

123

P. C. Hansen et al.

image is then blurred (cf. Fig. 5b) and cropped around the white box to the dimensions of the original image shown (Fig. 5c). Figure 5d shows a reconstruction where reflective BC were incorporated in the matrix A; the image has been very well reconstructed. Figure 5e shows a similar reconstruction, but with zero BC; in this case the image is very well reconstructed inside a circle while the corners of the image are ruined. This experiment illustrates that our treatment of boundary conditions by incorporating different kinds of “fill” can be expected to have a dramatic influence on the reconstruction. For other types of rotation, there is a similar effect when there is a BC mismatch, but the shape and size of the reconstructible area depends on the particular rotation. 5.3 The “ghosts” Image deblurring often introduces artifacts in the reconstruction image. For example, in the case of out-of-focus blur or Gaussian blur, ringing artifacts tend to arise around edges in the image—this is a Gibbs phenomenon due to the damping of high spatial frequencies in the reconstructed image. In case of motion blur, the PSF can be considered as a piecewise constant function, and it tends to introduces replicas of sharp edges in the reconstruction. We refer to these artifacts as “ghosts”, and they are plainly visible in reconstructions b and c in Fig. 4; they are also visible in several works on motion deblurring, e.g., [37,39,44]. To illustrate how these “ghosts” arise in an iterative reconstruction method, consider the basic Landweber algorithm without projection. It follows from (4.1) and (4.2) that each iterate is given by a matrix polynomial in A T A times the vector A T g, and hence the iterate f (k) lies in the Krylov subspace spanned by the vectors A T g, A T A A T g, (A T A)2 A T g, . . ., (A T A)k−1 g. If the vectors w 1k, w2 , . . . , wk form an orthonormal T (k) are basis for this Krylov subspace then f (k) = i=1 ci wi , where ci = wi f coefficients. For illustration, we consider a simple 1D example. Figure 6 shows a rectangular PSF, an exact piecewise function, the blurred version of this function, together with the first nine components ci wi for i = 1, 2, . . . , 9. In this example, the first four components are smooth, while the fifth component c5 u 5 has large and spiky features that are located precisely at the edges of the exact function—and since the reconstruction has this vector as one of its components, the oscillations will also be present in the reconstruction. Thus, the “ghosts” are in fact due to the repeated multiplications with A T A that underlie all the iterative methods: they correspond to convolutions of the image by a piecewise constant point-spread function, and these convolutions give rise to certain basis vectors for the solution that have distinct localized features at the positions of sharp edges in the blurred image g. 6 Stopping rules In order to provide a more complete algorithm, we integrated stopping rules in the reconstruction algorithms. Ideally, this allows the reconstruction algorithm to auto-

123

Rotational image deblurring

PSF

Exact function

0.2 0.1 0

1

0.5

0.5

0 0

10

20

30

0 0

10

i=1 0.6 0.4 0.2 0 −0.2

20

0.6 0.4 0.2 0 −0.2

0.6 0.4 0.2 0 −0.2

30

0

i=2

i=4

20

30

0.6 0.4 0.2 0 −0.2

0.6 0.4 0.2 0 −0.2

i=6 0.6 0.4 0.2 0 −0.2

i=8 0.6 0.4 0.2 0 −0.2

10

i=3

i=5

i=7 0.6 0.4 0.2 0 −0.2

Blurred function

1

i=9 0.6 0.4 0.2 0 −0.2

Fig. 6 1D illustration of the appearance of “ghosts”. The three bottom rows show the components ci wi of the reconstruction in the orthonormal basis vectors of the Krylov subspace. Spikes are clearly visible for i =5

matically stop at a point where the solution is closest to the exact one, i.e., at the point of semi-convergence. A stopping rule that often works well, in theory, is the “discrepancy principle” which amounts to stopping when the RMS of the residual vector approximately equals the RMS of the errors in the data. In practise, however, this method only works well given a good estimate of the variance of the noise, which is not available in many applications. Hence, we focus on stopping rules that are independent of this information, to make the algorithms more generally useful. See [31] for a recent survey of stopping rules for inverse problems. 6.1 Normalized cumulative periodogram Assuming that the blurred image is corrupted by white noise, at some iteration the residual vector r (k) = g − A f (k) should have the spectral characteristics of this noise type, i.e., its power spectrum should be approximately uniform. At earlier iterations, not all available information has been extracted from the data, and the residual still resembles a signal; at later iterations the residual becomes dominated by the higher frequency noise components. The point of semi-convergence is precisely when the residual resembles white noise the most, for then we have extracted all available signal from the data.

123

P. C. Hansen et al.

In order to practically verify whether the spectral characteristics of the residual are similar to those of white noise one can make use of the NCP. Here we summarize the derivation of the NCP criterion from [15,33]. Let R (k) denote the N × N residual image obtained by reshaping the residual vector r (k) into the original image dimensions, and let P (k) denote the q × q matrix (q = N /2 + 1) with the corresponding 2D power spectrum of R (k) . The elements of P (k) are:

  2 (k) Pi j = dft2 R (k) i j , i, j = 1, . . . , q = N /2 + 1, where “dft2” denotes the 2D discrete Fourier transform. First the elements of the q 2 vector vec(P (k) ) must then be reordered in order of increasing spatial frequency. For this purpose a new q × q matrix is introduced with elements i j = i 2 + j 2 , and then we define the q 2 × q 2 permutation matrix Π such that the elements of Π vec( ) appear in nondecreasing order. Then the permuted vector pˆ (k) = Π vec(P (k) ) holds all the power spectrum elements in order of increasing spatial frequency. We can now define the NCP as the vector with elements  (k)   pˆ (2: j+1) (k)  1 , j = 1, . . . , q 2 − 1. (6.1) c(R ) j =  (k)  pˆ (2:q 2 ) 1

If the elements of R (k) are white noise then c(R (k) ) will resemble a straight line when plotted. Specifically, if the NCP curve for the residual lies within the so-called Kolmogorov–Smirnoff limits ±1.35q of the straight line that represents the ideal white-noise NCP, then we accept the given NCP c(R (k) ) as representing white noise with a 5 % significance level and stop the iterations. 6.2 Monte-Carlo GCV This stopping rule is based on the classical method of cross validation. The underlying idea is to stop the iteration at a point where the reconstructed data A f (k) can best predict the exact data A f ex . The method of GCV (see, e.g., [16,29]) is a robust way to implement this approach, where we must the k that minimizes the GCV function:   g − A f (k) 2 2 (k) =   2 , trace I − A A#k

(6.2)

in which I the identity matrix and A#k the influence matrix such that f (k) = A#k g. Explicitly working with the matrix A#k would require a formidable computational effort, and thus the trace term has to be estimated. Girard [13] and Hutchinson [18] proposed two similar methods for estimating this term based on a sequence of MonteCarlo experiments. In this work, the Hutchinson method was used as it seemed to be slightly more robust. According to this method, the trace of a matrix M can be estimated as y T M y, where y is a vector whose elements are samples of a discrete random variable that assumes

123

Rotational image deblurring −4

x 10

−4

Landweber GCV

6

5.8

6

5.7

5.8 5.6

GCV envelope

5.6 5.5 5.4

5.4 5.2 10

MRNSD GCV

5.9

6.2

GCV(k)

GCV(k)

6.4

x 10

5.3 5.2 20

30

40

50

60

70

80

Iterations (k)

10

20

30

40

50

60

70

80

Iterations (k)

Fig. 7 Monte-Carlo GCV functions for 5 % Gaussian noise

values of ±1 with equal probability. In order to estimate the GCV denominator, we need an additional vector z (k) = A#k y, which can be interpreted as a “reconstructed” version of y at the kth iteration. An update for the z (k) vector can be calculated in a similar manner to the iteration vector f (k) ; see [26] for Landweber and [1] for MRNSD. Using (6.2) the GCV function can then be estimated by the Monte-Carlo GCV function:   g − A f (k) 2 2 MC (k) =  (6.3) 2 . y T y − y T A z (k) Typically, the Monte-Carlo GCV stopping rule increases the computational requirements of the reconstruction algorithm by a factor of two. In some cases this might pose a serious disadvantage; however, for our sparse blurring matrix model the increase in computational complexity is acceptable. 6.3 Performance studies Here we present experimental results regarding the stopping rules. In order to better study their behavior the “inverse crime” approach was used to generate the blurred image which, as in the previous experiments, is corrupted by a 20◦ pure rotational blur. In order to obtain a better understanding of the robustness of the stopping rules, two noise levels were used, 2 and 5 % of Gaussian and Poisson noise. A first observation about the Monte-Carlo GCV can be made by studying its behavior shown in Fig. 7. In both cases the GCV function MC (k) has a minimum, but for MRNSD there were always some small oscillations in MC (k) close to the global minimum. This behavior was also observed in [26]. Since the minimum in this case is fairly close to the optimal stopping point, we disregard the oscillations by considering the “envelope” of the GCV function. Specifically, we interpolate all the peaks of the oscillations to produce a smooth envelope of the GCV function. This can be done for every iteration with a relatively negligible computational cost, and it enables us to approximate the global GCV minimum quite well in order to stop the algorithm at this point.

123

P. C. Hansen et al.

0.36

Pure rotation 5% Gaussian noise

0.35 0.34

error history min. error GCV NCP

0.32

0.32 0.31 0.3

0.31 0.3 0.29

0.29

0.28

0.28

0.27 0.26

0.27 20

40

60

80

100

10

20

30

40

50

60

70

80

Iterations (k)

Iterations (k)

Pure rotation 2% Gaussian noise

Pure rotation 2% Gaussian noise

error history min. error GCV

0.224 0.222

0.225

0.22 0.218

0.22 0.215 0.21

0.216

0.205

0.214 100

error history min. error GCV NCP

0.23

RE(k)

RE(k)

Pure rotation 5% Gaussian noise

0.33

RE(k)

0.33

RE(k)

0.34

error history min. error GCV NCP

150

200

250

300

0.2

350

60

80

100

120

140

160

180

Iterations (k)

Iterations (k)

Fig. 8 Stopping rule performance for 5 and 2 % Gaussian noise for PL (left) and MRNSD (left). The markers show the minima and the points chosen by the stopping rules Table 2 Stopping rule performance for Poisson noise

Algorithm

Min error ρ

k

NCP

GCV ρ

k

ρ

k

5 % Poisson noise PL

46

0.2716

46

0.2716

57

0.2729

MRNSD

33

0.2681

35

0.2682

64

0.2960

212

0.2149

492

0.2359

200

0.2150

95

0.2117

221

0.2385

130

0.2193

2 % Poisson noise PL MRNSD

Figure 8 shows the error histories for 5 and 2 % Gaussian noise. The optimal stopping points and the ones chosen by the stopping rules are marked on the curves. In the case of 5 % noise, we can see fairly good performance of both stopping rules. For a noise level of 2 % the NCP criterion stops MRNSD too late, and it fails to stop the PL iterations. This might be due to the fact that for lower noise levels more noise components are missing from the residual, thus making its spectrum differ substantially from that of white noise. On the other hand, the Monte-Carlo GCV rule seems to be quite robust, stopping the iterations fairly close to the optimal point.

123

Rotational image deblurring

The corresponding experiments with Poisson noise are summarized in Table 2, and the results are very similar in most cases. For higher levels of Poisson noise the NCP criterion seems to perform better. We conclude that the overall performance of the Monte-Carlo GCV criterion is superior and quite satisfactory in most cases for these series of “inverse crime” experiments. 7 Simulation studies In this final section we present results using more realistic simulations, i.e., without committing “inverse crime”. 7.1 Creating realistic test problems In order to test the robustness of our algorithms and study their behavior under realistic circumstances we developed a simulator that produces realistic test data for rotational motion blur. The goal is to simulate the capturing of a continuous scene under a rotational motion blur scheme by creating a “high resolution” simulation of the forward problem. This means that we want to take into account (1) the discretization effects when capturing a continuous scene, (2) the temporal integration of the captured image described by (3.3), and (3) the scene outside the recorded image. In order to achieve these requirements, the realistic blur simulation consists of the following steps: – We assume a continuous planar scene in the form of a high-resolution test image. A resolution of 2,000 × 2,000 is sufficient to approximate the discretization effect for a 400 × 400 target image. – The high-resolution test image is blurred using a matrix-free forward model. We applied a pure rotation of 20◦ ; a step of 0.5◦ was sufficient to approximate the temporal integration of the captured image. – In order to simulate the effect of the scene surrounding the image, the high√ resolution blurred image is cropped by a factor of 2 which, for a pure rotation, eliminates all edge effects in the cropped blurred image. – To approximate the discretization of the continuous scene being recorded, the blurred and cropped high-resolution image is down-sampled by averaging neighborhoods of pixel values to produce the down-sampled image. – Gaussian or Poisson noise is added to the down-sampled image to simulate the noise usually present in captured images. 7.2 Simulation results As test image we used a high-resolution image STS105_707_019HR.JPG of the International Space Station from http://cadmos.cnes.fr/modules/resources/download/ phototheque/ISS. The image is fairly complex, containing both high and low frequency detail, and it was blurred using the simulator described in the previous section.

123

P. C. Hansen et al.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Fig. 9 Simulation experiments. All results are for 2 % Gaussian noise except (d) which is for 2 % Poisson noise. Images e–l show the best reconstructions within the circular region for different BCs, while c and d show the best reconstructions when the complete region is considered. a Exact image, b blurred image, c Gaussian noise, ρ = 0.2602, s = 0.7194; d Poisson noise, ρ = 0.2623, s = 0.7082; e PL, zero BC, ρ = 0.2362, S = 0.7571; f PL, replicate BC, ρ = 0.2415, s = 0.7550; g PL, reflective BC, ρ = 0.2428, s = 0.7430; h PL, periodic BC, ρ = 0.2594, s = 0.7242; i MRNSD, zero BC, ρ = 0.2250, S = 0.7957; j MRNSD, replic. BC, ρ = 0.22384, S = 0.7963; k MRNSD, reflect. BC, ρ = 0.2231, S = 0.7784; l MRNSD, periodic BC, ρ = 0.2492, S = 0.7610

In Sect. 5.2 we illustrated the effects of a boundary condition mismatch. For a pure rotation, when the boundary conditions do not reflect the actual scene outside the frame, the reconstructible area is limited to a circle—while the corners of the reconstructed image are practically ruined. Since the boundary conditions only try to approximate the scene outside the frame, it is obvious that a good reconstruction cannot be guaranteed in the entire image domain. This justifies the use of a pure rotational motion blur for the simulation experiments as we know that the reconstructible area is a circle. Therefore, in this section we focus on evaluating the quality of the reconstruction within this circle, by applying a circular binary mask to both the reconstructed image and the exact image, before computing the quality measures ρ and S. Figure 9 shows the reconstructed images within the circular mask using the four different BCs. We observe that different boundary conditions can have an impact even

123

Rotational image deblurring Table 3 Simulation results for 2 % Poisson noise

BC

PL ρ

MRNSD S

ρ

S

Zero

0.2368

0.7073

0.2288

0.7370

Replicate

0.2420

0.7011

0.2274

0.7382

Reflective

0.2436

0.6871

0.2366

0.7140

Periodic

0.2602

0.6622

0.2518

0.6926

within the circle defined by the mask. In particular, the replicate BC seems to give best results in most cases. Similar results, obtained with Poisson noise, are summarized in Table 3. We note that the best overall reconstruction for this image is obtained with the MRNSD algorithm using the replicate BC. This holds true for both Gaussian and Poisson noise and the full reconstructions are shown in Fig. 9c, d; note that the corners of the images (outside the circle) are surprisingly well reconstructed.

8 Conclusion We developed iterative deblurring algorithms with robust stopping rules that can handle rotational blur along an arbitrary axis (very few other algorithms can do that). Efficiency of our algorithms is achieved through the use of a sparse-matrix representation of the blurring operators, which allows us to easily incorporate different types of boundary conditions. Our experiments demonstrate that we are able to achieve good reconstructions, and that the use of boundary conditions is essential for this. Acknowledgments We thank Henrik Aanæs for his invaluable help with the geometric aspects, and John Bardsley for his help with the implementation of the randomized GCV method.

References 1. Bardsley, J.M.: Stopping rules for a nonnegatively constrained iterative method for ill-posed Poisson imaging problems. BIT 48, 651–664 (2008) 2. Bardsley, J.M., Merikoski, J.K., Vio, R.: The stabilizing properties of nonnegativity constraints in least-squares image reconstruction. Int. J. Pure Appl. Math. 43, 95–109 (2008) 3. Bardsley, J.M., Nagy, J.G.: Covariance-preconditioned iterative methods for nonnegatively constrained astronomical imaging. SIAM J. Matrix Anal. Appl. 27, 1184–1197 (2006) 4. Ben-Ezra, M., Nayar, S.K.: Motion based motion deblurring. IEEE Trans. Pattern Anal. Mach. Int. 26(6), 689–698 (2004) 5. Boracchi, G., Caglioti, V., Danese, A.: Estimating camera rotation parameters from a blurred image. In: Proceedings of the 3rd International Conference on Computer Vision Theory and Applications (VISAPP 2008), Funchal (2008) 6. Brianzi, P., Di Bendetto, F., Estatico, C.: Improvement of space-invariant image deblurring by preconditioned Landweber iterations. SIAM J. Sci. Comput. 30, 1430–1458 (2008) 7. Dai, S., Wu, Y.: Motion from blur. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 1–8 (2008) 8. Davis, T.A.: Direct Methods for Sparse Linear Systems. SIAM, Philadelphia (2006)

123

P. C. Hansen et al. 9. Elfving, T., Hansen, P.C., Nikazad, T.: Semi-convergence and relaxation parameters for projected SIRT algorithms. SIAM J. Sci. Comp. 34, A2000–A2017 (2012) 10. Estatico, C., Di Bendetto, F.: Shift-invariant approximations of structured shift-variant blurring matrices. Numer. Algorithms 62, 615–635 (2013) 11. Faber, T.L., Raghunath, N., Tudorascu, D., Votaw, J.R.: Motion correction of PET brain images through deconvolution: I. Theoretical development and analysis in software simulations. Phys. Med. Biol. 54(3), 797–811 (2009) 12. Fan, Y.W., Nagy, J.G.: Synthetic boundary conditions for image deblurring. Linear Algebra Appl. 434, 2244–2268 (2010) 13. Girard, A.: A fast Monte-Carlo cross-validation procedure for large least squares problems with noisy data. Numerische Mathematik 56, 1–23 (1989) 14. Hansen, P.C.: Discrete Inverse Problems: Insight and Algorithms. SIAM, Philadelphia (2010) 15. Hansen, P.C., Kilmer, M.E., Kjeldsen, R.H.: Exploiting residual information in the parameter choice for discrete ill-posed problems. BIT 46, 41–59 (2006) 16. Hansen, P.C., Nagy, J.G., O’Leary, D.P.: Deblurring Images: Matrices, Spectra, and Filtering. SIAM, Philadelphia (2006) 17. Hartley, R., Zisserman, A.: Multiple View Geometry in Computer Vision. Cambridge University Press, New York (2003) 18. Hutchinson, M.: A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines. Commun. Stat. Simul. Comput. 18, 1059–1076 (1989) 19. Kang, S.K., Min, J.H., Paik, J.K.: Segmentation-based spatially adaptive motion blur removal and its application to surveillance systems. In: Proceedings of the International Conference Image Processing, vol. 1, pp. 245–248 (2001) 20. Kaufman, L.: Maximum likelihood, least squares, and penalized least squares for PET. IEEE Trans. Med. Imaging 12, 200–214 (1993) 21. Kelley, C.T.: Iterative Methods for Optimization. SIAM, Philadelphia (1999) 22. Nagy, J.G., O’Leary, D.P.: Restoring images degraded by spatially variant blur. SIAM J. Sci. Comput. 19(4), 1063–1082 (1998) 23. Nagy, J.G., Palmer, K.: Steepest descent, cg, and iterative regularization of ill-posed problems. BIT Numer. Math. 43(5), 1003–1017 (2003) 24. Nagy, J.G., Strakoš, Z.: Enforcing nonnegativity in image reconstruction algorithms. In: Wilson, D.C., et al. (eds.) Proceedings of SPIE. Mathematical modeling, estimation, and imaging, vol. 4121, pp. 182–190. (2000). doi:10.1117/12.402439 25. Nocedal, J., Wright, S.: Numerical Optimization. Springer, New York (1999) 26. Perry, K., Reeves, S.: A practical stopping rule for iterative signal restoration. IEEE Trans. Signal Proces. 42(7), 1829–1833 (1994) 27. Piana, M., Bertero, M.: Projected Landweber method and preconditioning. Inverse Probl. 13, 441–463 (1997) 28. Raghunath, N., Faber, T.L., Suryanarayanan, S., Votaw, J.R.: Motion correction of PET brain images through deconvolution: II. Practical implementation and algorithm optimization. Phys. Med. Biol. 54(3), 813–829 (2009) 29. Reeves, S.J.: Generalized cross-validation as a stopping rule for the Richardson–Lucy algorithm. Int. J. Imaging Syst. Technol. 6(4), 387–391 (1995) 30. Reeves, S.J.: Fast image restoration without boundary artifacts. IEEE Trans. Image Process. 14, 1448– 1453 (2005) 31. Reichel, L., Rodriguez, G.: Old and new parameter choice rules for discrete ill-posed problems. Numer. Algorithms 63, 65–87 (2013) 32. Ribaric, S., Milani, M., Kalafatic, Z.: Restoration of images blurred by circular motion. In: Proceedings of the First International Workshop on Image and Signal Processing and Analysis (IWISPA 2000), pp. 53–60. IEEE (2000) 33. Rust, B.W., O’Leary, D.P.: Residual periodograms for choosing regularization parameters for ill-posed problems. Inverse Probl. 24, 034005 (2008). doi:10.1088/0266-5611/24/3/034005 34. Sawchuk, A.A.: Space-variant image motion degradation and restoration. Proc. IEEE 60, 854–861 (1972) 35. Sawchuk, A.A.: Space-variant image restoration by coordinate transformations. J. Opt. Soc. Am. 64, 138–144 (1974)

123

Rotational image deblurring 36. Serra-Capizzano, S.: A note on antireflective boundary conditions and fast deblurring models. SIAM J. Sci. Comput. 25(4), 1307–1325 (2003) 37. Shan, Q., Jia, J., Agarwala, A.: High-quality motion deblurring from a single image. ACM Trans. Graph. (SIGGRAPH 2008) 27(3) (2008). Article 73. doi:10.1145/1399504.1360672 38. Shan, Q., Xiong, W., Jia, J.: Rotational motion deblurring of a rigid object from a single image. In: Proceedings of the 11th International Conference on Computer Vision (ICCV 2007), pp. 1–8. IEEE (2007) 39. Tai, Y.W., Tan, P., Brown, M.: Richardson–Lucy deblurring for scenes under a projective motion path. IEEE Trans. Pattern Anal. Mach. Intell. 33(8), 1603–1618 (2011) 40. Trussell, H.J., Fogel, S.: Identification and restoration of spatially variant motion blurs in sequential images. IEEE Trans. Image Process. 1, 123–126 (1992) 41. Tull, D.L., Katsaggelos, A.K.: Iterative restoration of fast-moving objects in dynamic image sequences. Opt. Eng. 35(12), 3460–3469 (1996) 42. Vogel, C.R.: Computational Methods for Inverse Problems. SIAM, Philadelphia (2002) 43. Wang, Z., Bovik, A., Sheikh, H., Simoncelli, E.: Image quality assessment: from error visibility to structural similarity. IEEE Trans. Image Process. 13(4), 600–612 (2004) 44. Whyte, O., Sivic, J., Zisserman, A., Ponce, J.: Non-uniform deblurring for shaken images. Int. J. Comput. Vis. 98, 168–186 (2012)

123