Antireflective Boundary Conditions for Deblurring Problems

Abstract

This survey paper deals with the use of antireflective boundary conditions for deblurring problems where the issues that we consider are the precision of the reconstruction when the noise is not present, the linear algebra related to these boundary conditions, the iterative and noniterative regularization solvers when the noise is considered, both from the viewpoint of the computational cost and from the viewpoint of the quality of the reconstruction. In the latter case, we consider a reblurring approach that replaces the transposition operation with correlation. For many of the considered items, the anti-reflective algebra coming from the given boundary conditions is the optimal choice. Numerical experiments corroborating the previous statement and a conclusion section end the paper.

1. Introduction

Formation of a blurred signal/image is typically modeled as a linear system:𝐴𝐟=𝐠, (1) where 𝐟 is the true object, 𝐠 is the blurred object, and 𝐴 models the blurring process. Obtaining an accurate deblurring model (i.e., the matrix 𝐴) requires essentially two main pieces of information: (1)identification of the blur operator, called a point spread function (PSF),(2)choosing an appropriate boundary condition (BC), assuming that the observed image is always finite.

The identification of the blur operator is related to the infinite dimensional problem and it decides the essential structure of the matrix 𝐴. In the spatially invariant case, due to the shift invariance of the blurring, for image deblurring we deduce a two-level Toeplitz structure (i.e., a two-level Toeplitz matrix or in a different language a block Toeplitz matrix with Toeplitz blocks). The choice of the BCs influences small sections of 𝐴 by a low-rank plus low-norm term. However, these correction matrices have a substantial impact in two important directions: (a)precision of the reconstruction especially close to the boundaries (presence of ringing effects); (b)cost of the computation for recovering the “true” image from the blurred one with or without noise.

On the other hand, the involved correction matrices do not modify significantly the spaces of ill-conditioning (where the coefficient matrix 𝐴 shows small eigenvalues). This happen because the global matrix is of convolution type. Therefore, the eigenvectors look globally like the Fourier vectors 𝐦_𝑘 (they are exactly the Fourier vectors when periodic BCs are imposed). Conversely, the small rank matrices induced by the chosen BCs are not generic since they are necessarily localized in space. More specifically, in one dimension there exists 𝑘≪𝑛, 𝑛 being the size of the matrix, such that the low-rank correction belongs to the span of canonical matrices 𝐸_𝑖,𝑗=𝐞_𝑖𝐞^𝑇_𝑗 with 𝑖,𝑗≤𝑘 and/or 𝑖,𝑗≥𝑛+1−𝑘, 𝐞_𝑘 being the 𝑘th vector of the canonical basis. In that case, we observe ‖𝐸_𝑖,𝑗𝐦_𝑘‖₂≪‖𝐦_𝑘‖₂ with ‖⋅‖₂ being the classical Euclidean norm and hence the term 𝐸_𝑖,𝑗 is unable to modify significantly the characterization in frequencies of the eigenspaces. In actuality, due to the relation ‖𝐸_𝑖,𝑗𝐦_𝑘‖₂≪‖𝐦_𝑘‖₂, the effect induced by any 𝐸_𝑖,𝑗 is more or less uniformly distributed in all the Fourier directions.

For image deblurring, a classical choice is to use zero (Dirichlet) BCs (see Andrews and Hunt [1], pp. 211–220) where we assume that the values of the image 𝐟 outside the domain of consideration are zero. The corresponding blurring matrix is two-level Toeplitz, which is known to be computationally expensive to invert when direct solvers are used (see for instance the review [2]): the fastest available direct Toeplitz solver is of 𝑂(𝑛⁴) operations for an 𝑛²-by-𝑛² two-level Toeplitz matrix [3]. In addition, we have no good news when applying iterative methods due to the negative results in [4–6], where it is proved that we should not expect the classical matrix algebra preconditioners to lead to optimal iterative solvers in the ill-conditioned case. More specifically, we lose the strong cluster at unity of the preconditioned matrix sequence [4, 5] and the spectral boundedness of the preconditioned matrix sequence and its inverse [6]. Therefore, by the convergence analysis of Axelsson and coworkers (see e.g., the pioneering paper [7]), we cannot expect optimal convergence (i.e., a convergence rate independent of the matrix size) both for Krylov type methods and for classical stationary solvers (with the partial exception of multigrid type techniques, see [8, 9] for a recent adaptation in the context of image deblurring). In image restoration problems, we consider regularized problems so that the previous statement is less relevant; however, for linear systems associated to the Tikhonov regularization, the optimality could guarantee a convergence speed independent of the regularizing parameter (for a discussion on this kind of robustness see [9, 10]). Moreover, an image has boundaries, and if the true image is not close to zero at the boundaries, it means that Dirichlet BCs may introduce an artificial discontinuity, which in turn results in the well-known Gibbs phenomenon (observed as a ringing effect in the reconstructed image). Therefore, we face substantial problems both with respect to (a) and (b). A possibility for drastically reducing the computational cost is to impose periodic BCs since this implies that 𝐴 is a two-level circulant matrix which can be efficiently inverted by means of few fast Fourier transforms (FFTs). In this case, we have a good solution to problem (b) (see Gonzalez and Woods, [11] p. 258) but often the image at the right (upper) boundary has nothing to do with the image at the left (lower) boundary. Again this implies that we are introducing an artificial discontinuity which in turn leads to ringing effects. An important proposal in the direction of reducing the ringing effect and of designing fast algorithms is the one by Ng et al. [12] (see also Strang, [13] page 145, and Chan et al. [14]). In particular, they consider Neumann (reflective) BCs, which means that the data outside the domain of consideration is taken as a reflection of the data inside. The approach is fast for circularly symmetric PSFs since the corresponding blurring matrix belongs to a special two-level algebra simultaneously diagonalized by the fast cosine transform (DCT III). Because a two-level DCT III requires only real multiplications and can be done at half of the cost of a single FFT (see Rao and Yip, [15] pp. 59-60), inversion of these matrices is substantially faster than those obtained from zero BCs and is comparable to (or little better than) the case of periodic BCs. Thus, for circularly symmetric PSFs, requirement (b) is satisfied when using reflective BCs. Moreover, the artificial boundary discontinuities are eliminated, and therefore the ringing effects are strongly reduced. However, although reflection guarantees continuity, it generally fails to guarantee continuity of the normal derivative except for the nongeneric case where the normal derivative is zero at the boundary. Therefore if the image is smooth, then the reflective BCs only save the continuity but introduce an artificial discontinuity of the normal derivative (though requirement (b) is substantially improved with respect to periodic and Dirichlet BCs). In the image processing literature, other methods have also been proposed to assign boundary values as a local image mean, or are obtained by extrapolation methods in order to insure continuity (see Lagendijk and Biemond [16], p. 22, and the references therein). However, these techniques that meet the approximation requirement (a), generally do not satisfy the computational requirement (b) so we still have to face a difficult linear algebra problem.

Here we follow the approach in [17] in which new BCs are proposed that can further reduce ringing effects and that, under the assumption of symmetry of the PSF, leads to a matrix 𝐴 belonging to the sine algebra of type I for which fast (real) algorithms are available. More precisely, we consider the use of antireflective BCs (AR-BCs) where the data outside the domain of consideration are taken as an antireflection of the data inside. Consider the original signal 𝐟̃=(…,𝑓_−𝑚+1,…,𝑓₀,𝑓₁,…,𝑓_𝑛,𝑓_𝑛+1,…,𝑓_𝑛+𝑚,…)^𝑇 and the normalized blurring function given byℎ=(…,0,0,ℎ_−𝑚,ℎ_−𝑚+1,…,ℎ₀,…,ℎ_𝑚−1,ℎ_𝑚,0,0,…)^𝑇,_{𝑚𝑗=−𝑚}ℎ_𝑗=1. (2) The blurred signal 𝐠 is the convolution of ℎ and 𝐟̃ and consequently 𝑔_𝑖=∑^∞_𝑗=−∞ℎ_𝑗𝑓_𝑖−𝑗. The deblurring problem is to recover the vector 𝐟=(𝑓₁,…,𝑓_𝑛)^𝑇 given the blurring function ℎ and a blurred signal 𝐠=(𝑔₁,…,𝑔_𝑛)^𝑇 of finite length, that is,𝐴̃ℎ𝐠=𝐟,𝐴=^𝑚⋯ℎ₀⋯ℎ_−𝑚ℎ_𝑚ℎ₀ℎ_−𝑚0⋱⋱⋱⋱⋱⋱⋱⋱⋱⋱0ℎ_𝑚ℎ₀ℎ_−𝑚ℎ_𝑚⋯ℎ₀⋯ℎ_−𝑚. (3) Thus the blurred signal 𝐠 is determined not by 𝐟 only, but also by (𝑓_−𝑚+1,…,𝑓₀)^𝑇 and (𝑓_𝑛+1,…,𝑓_𝑛+𝑚)^𝑇 and the linear system (3) is underdetermined. To overcome this, we make certain assumptions (called BCs) on the unknown boundary data 𝑓_−𝑚+1,…,𝑓₀ and 𝑓_𝑛+1,…,𝑓_𝑛+𝑚 in such a way that the number of unknowns equals the number of equations.

In terms of the objects defined so far, we recall that zero (Dirichlet) BCs means 𝑓_𝑗=𝑓_𝑛+𝑗=0 for all 𝑗 in (3) so that a Toeplitz structure is encountered. If we consider periodic BCs we set 𝑓_𝑗=𝑓_𝑛+𝑗, for all 𝑗 in (3) and the matrix system in (3) becomes 𝑛-by-𝑛 circulant so that it can be diagonalized by the discrete Fourier matrix and the above system can be solved by using three FFTs (one for finding the eigenvalues of the blurring matrix and two for solving the system). For the Neumann BCs, we assume that the data outside 𝑓 are a reflection of the data inside 𝑓 (refer to [12]) so that 𝑓_1−𝑗=𝑓_𝑗 for 𝑗=1,…,𝑚 and 𝑓_𝑛+𝑗=𝑓_{𝑛+1−𝑗} for all 𝑗=1,…,𝑚 in (3). Thus (3) becomes 𝐶𝑓=𝑔, where 𝐶 is neither Toeplitz nor circulant but it is a special 𝑛-by-𝑛 Toeplitz plus Hankel matrix which is diagonalized by the discrete cosine transform provided that the blurring function ℎ is symmetric, that is, ℎ_𝑗=ℎ_−𝑗 for all 𝑗 in (2). It follows that the above system can be solved by using three transforms DCT III in 𝑂(𝑛log𝑛) operations (refer to [12]). The latter approach is computationally interesting since a DCT III requires only real operations and is about twice as fast as the FFT (see [15], pp. 59-60), and this is true in two dimensions as well. With the help of a different Toeplitz plus Hankel structure, we establish similar results for the antireflective BCs both in one dimension and two dimensions. It is worth finally to remark that La Spina has analyzed [18] the BCs associated with all the known trigonometric algebras: the interesting facts are two, first some matrix algebras do not lead to any boundary condition, second the highest precision is reached by the classical DCT III algebra and by the classical sine transform algebra of type I which ensures only the continuity of the signal. Therefore, in this sense we can claim that among the known algebras the antireflective one is that related to the most precise reconstructions.

The paper is organized as follows. In the Section 2 we examine the antireflective BCs, the related algebra, and the related transforms also from a computational viewpoint. At the end of Section 2 we briefly consider the multilevel extension while in Section 3 we study some regularization techniques specifically adapted to the features of the antireflective algebra. Finally, Section 4 contains numerical experiments both with Tikhonov like solvers (with reblurring) and with Krylov techniques with early termination as stopping criterion. A final section of conclusions and future problems end the paper.

2. The Algebra of AR Matrices

In this section we describe the AR-BCs and the algebra 𝒜ℛ_𝑛≡𝒜ℛ, 𝑛≥3, of matrices arising from the imposition of AR-BCs.

2.1. The Antireflective BCs

When defining the antireflective boundary conditions, we assume that the data outside 𝐟 are an antireflection of the data inside 𝐟. More precisely, if 𝑥 is a point outside the domain and 𝑥^∗ is the closest boundary point, then we have 𝑥=𝑥^∗−𝛿𝑥 and the quantity 𝑓(𝑥) is approximated by 𝑓(𝑥^∗)−(𝑓(𝑥^∗+𝛿𝑥)−𝑓(𝑥^∗)) so that we impose𝑓_1−𝑗=𝑓₁−(𝑓_𝑗+1−𝑓₁)=2𝑓₁−𝑓_𝑗+1,𝑓forall𝑗=1,…,𝑚,_𝑛+𝑗=𝑓_𝑛−(𝑓_𝑛−𝑗−𝑓_𝑛)=2𝑓_𝑛−𝑓_𝑛−𝑗,forall𝑗=1,…,𝑚 (4) in (3). If we define the vector 𝐳 whose components are 𝑧_𝑗=2∑^𝑚_𝑘=𝑗ℎ_𝑘 for 𝑗≤𝑚 and zero otherwise and if we define the vector 𝐰 whose components are 𝑤_{𝑛+1−𝑗}=2∑^𝑚_𝑘=𝑗ℎ_−𝑘 for 𝑗≤𝑚 and zero otherwise, then (3) becomes𝐴𝐟̃=[𝐳𝐞^𝑇₁−(0∣𝑇_𝑙)𝐽̃+𝑇−(𝑇_𝑟∣0)𝐽̂+𝐰𝐞^𝑇_𝑛]𝐟̃=𝐠, (5) where 𝐞_𝑘 is the 𝑘th vector of the canonical basis, 𝐽̃=(⁰⁰_0𝐽), 𝐽̂=(^𝐽0₀₀) with 𝐽 denoting the 𝑛−1 dimensional flip matrix having entries [𝐽]_𝑠,𝑡=1 if 𝑠+𝑡=𝑛+1 and zero otherwise, and where the matrices 𝑇_𝑙,𝑇, and 𝑇_𝑟 are given byℎ_𝑚⋯ℎ₁ℎ⋱⋱_𝑚0,ℎ₀⋯ℎ_−𝑚0ℎ⋮⋱⋱⋱_𝑚⋱⋱⋱ℎ_−𝑚⋱⋱⋱⋮0ℎ_𝑚⋯ℎ₀,0ℎ_−𝑚ℎ⋮⋱₋₁⋯ℎ_−𝑚. (6) Therefore, the matrices (0∣𝑇_𝑙)𝐽̃ and (𝑇_𝑟∣0)𝐽̂ involved in (5) take the form0ℎ₁⋯ℎ_𝑚ℎ0⋯0⋮⋰_𝑚,ℎ⋮𝑂⋮0⋯00⋯0⋮𝑂⋮_−𝑚⋰⋮0⋯0ℎ_−𝑚⋯ℎ₋₁0. (7)

We remark that the coefficient matrix 𝐴 in (5) is neither Toeplitz nor circulant, but it is a Toeplitz plus Hankel plus 2 rank correction matrix, where the correction is placed at the first and the last column. We will show that the linear system (5) can be reduced to an (𝑛−2) by (𝑛−2) new system whose coefficient matrix can always be diagonalized by the discrete sine transform DST I (associated with the 𝜏 class) matrix provided that the blurring function ℎ is symmetric, that is, ℎ_𝑗=ℎ_−𝑗 for all 𝑗 in (2). It follows that (5) can be solved by using three FSTs in 𝑂(𝑛log𝑛) operations. This approach is computationally attractive as FST requires only real operations and is about twice as fast as the FFT and hence solving a problem with the AR BCs is twice as fast as solving a problem with the periodic BCs, has the same cost as solving a problem with the reflective BCs, and one gains one order of precision due to the 𝐶¹ continuity. Moreover, all these remarks stand in two dimensions as well and indeed an abstract treatment of the two-level and multilevel settings is reported in Section 2.5.

Finally it is worth mentioning that the use of AR BCs has been considered by several authors (see e.g., [18–26]). In particular in [24] the mean BCs have been introduced as a slight variation of the AR BCs. The order of approximation is the same but the hidden constants are smaller so that the mean BCs are generally more precise than the AR BCs. However, the resulting matrices do not form an algebra so that we cannot define a new transform and this is a confirmation of the negative result found in [18].

2.2. The 𝜏 Algebra

Let 𝑄 be the type I sine transform matrix of order 𝑛 (see [27]) with entries[𝑄]_𝑖,𝑗=2𝑛+1sin(𝑗𝑖𝜋𝑛+1),𝑖,𝑗=1,…,𝑛. (8) It is known that the real matrix 𝑄 is orthogonal and symmetric (𝑄⁻¹=𝑄^𝑇=𝑄). For any 𝑛-dimensional real vector 𝐯, the matrix-vector multiplication 𝑄𝐯 (DST-I transform) can be computed in 𝑂(𝑛log(𝑛)) real operations by using the algorithm FST-I. Let 𝜏 be the space of all the matrices that can be diagonalized by 𝑄:𝜏={𝑄𝐷𝑄∶𝐷isarealdiagonalmatrixofsize𝑛}. (9) Let 𝑋=𝑄𝐷𝑄∈𝜏, then 𝑄𝑋=𝐷𝑄. Consequently, if we let 𝐞₁ denote the first column of the identity matrix, then the relationship 𝑄𝑋𝐞₁=𝐷𝑄𝐞₁ implies that the eigenvalues [𝐷]_𝑖,𝑖 of 𝑋 are given by [𝐷]_𝑖,𝑖=[𝑄(𝑋𝐞₁)]_𝑖/[𝑄𝐞₁]_𝑖,𝑖=1,…,𝑛. Therefore, the eigenvalues of 𝑋 can be obtained by applying a DST-I transform to the first column of 𝑋 and, in addition, any matrix in 𝜏 is uniquely determined by its first column.

Now we report a characterization of the 𝜏 class, which is important for analyzing the structure of AR-matrices. Let us define the shift of any vector 𝐡=[ℎ₀,…,ℎ_𝑛−1]^𝑇 as 𝜎(𝐡)=[ℎ₁,ℎ₂,…,ℎ_𝑛−1,0]^𝑇. According to a Matlab-like notation, we define 𝑇(𝐱) to be the 𝑛-by-𝑛 symmetric Toeplitz matrix whose first column is 𝐱 and 𝐻(𝐱,𝐲) to be the 𝑛-by-𝑛 Hankel matrix whose first and last column are 𝐱 and 𝐲, respectively. Every matrix of the class (9) can be written as (see [27])𝑇(𝐡)−𝐻(𝜎²(𝐡),𝐽𝜎²(𝐡)), (10) where 𝐡=[ℎ₀,…,ℎ_𝑛−1]^𝑇∈ℝ^𝑛 and 𝐽 is the flip matrix. The structure defined by (10) means that matrices in the 𝜏 class are special instances of Toeplitz-plus-Hankel matrices.

Moreover, the eigenvalues of the 𝜏 matrix in (10) are given by the cosine function ℎ(𝑦) evaluated at the grid points vector 𝐺_𝑛=[𝑘𝜋/(𝑛+1)]^𝑛_𝑘=1, whereℎ(𝑦)=_{|𝑗|≤𝑛−1}ℎ_𝑗exp(𝐢𝑗𝑦), (11)𝐢²=−1, and ℎ_𝑗=ℎ_|𝑗| for |𝑗|≤𝑛−1. The 𝜏 matrix in (10) is usually denoted by 𝜏(ℎ) and is called the 𝜏 matrix generated by the function or symbol ℎ=ℎ(⋅) (see the seminal paper [27] where this notation was originally proposed).

2.3. The AR-Algebras 𝒜ℛ

Let ℎ=ℎ(⋅) be a real-valued cosine polynomial of degree 𝑙 and let 𝜏_𝑘(ℎ)≡𝑄diag(ℎ(𝐺_𝑘))𝑄 (note that 𝜏_𝑘(ℎ) coincides with the matrix in (10)-(11), when 𝑙≤𝑘−1). Hence the Fourier coefficients of ℎ are such that ℎ_𝑖=ℎ_−𝑖∈ℝ with ℎ_𝑖=0 if |𝑖|>𝑙, and for 𝑘=𝑛−2 we can define the one-level AR_𝑛(⋅) operatorAR_𝑛(ℎ)=𝐯ℎ(0)_𝑛−2(ℎ)𝜏_𝑛−2(ℎ)𝐽𝐯_𝑛−2(ℎ)ℎ(0), (12) where 𝐯_𝑛−2(ℎ)=𝜏_𝑛−2((𝜙(ℎ))(⋅))𝐞₁ and(𝜙(ℎ))(𝑦)=ℎ(𝑦)−ℎ(0).2(cos(𝑦)−1) (13) It is interesting to observe that ℎ(𝑦)−ℎ(0) has a zero of order at least 2 at zero (since ℎ is a cosine polynomial) and therefore 𝜙(ℎ)=(𝜙(ℎ))(⋅) is still a cosine polynomial of degree 𝑙−1, whose value at zero is −ℎ^′′(0)/2 (in other words the function is well defined at zero).

As proved in [28], with the above definition of the operator AR_𝑛(⋅), we have

(1) 𝛼AR_𝑛(ℎ₁)+𝛽AR_𝑛(ℎ₂)=AR_𝑛(𝛼ℎ₁+𝛽ℎ₂),

(2) AR_𝑛(ℎ₁)AR_𝑛(ℎ₂)=AR_𝑛(ℎ₁ℎ₂),

for real 𝛼 and 𝛽 and for cosine functions ℎ₁=ℎ₁(⋅) and ℎ₂=ℎ₂(⋅).

These properties allow us to define 𝒜ℛ as the algebra (closed under linear combinations, product, and inversion) of matrices AR_𝑛(ℎ), with ℎ being a cosine polynomial. By standard interpolation arguments it is easy to see that 𝒜ℛ can be defined as the set of matrices AR_𝑛(ℎ), where ℎ is a cosine polynomial of degree at most 𝑛−3. Therefore, it is clear that dim(𝒜ℛ)=𝑛−2. Moreover, the algebra 𝒜ℛ is commutative thanks to item 2, since ℎ₁(𝑦)ℎ₂(𝑦)=ℎ₂(𝑦)ℎ₁(𝑦) at every 𝑦. Consequently, if matrices of the form AR_n(ℎ) are diagonalizable, then they must have the same set of eigenvectors [29]. This means there must exist an “antireflective transform" that diagonalizes the matrices in 𝒜ℛ. Unfortunately this transform fails to be unitary, since the matrices in 𝒜ℛ are generically not normal. However the AR-transform and its inverse are close in rank to orthogonal linear mappings and only moderately ill conditioned.

Following the analysis given in [17], if the blurring function (the PSF) 𝐡 is symmetric (i.e., ℎ_𝑖=ℎ_−𝑖,forall𝑖∈ℤ), if ℎ_𝑖=0 for |𝑖|≥𝑛−2 (degree condition), and if 𝐡 is normalized so that ∑^𝑚_𝑖=−𝑚ℎ_𝑖=1, then the structure of the 𝑛×𝑛 antireflective blurring matrix 𝐴 is𝑧𝐴=⁰0^𝑇0𝑧₁𝑧_𝑚⋮𝑧𝐴⋮_𝑚𝑧₁00^𝑇𝑧₀, (14) where 𝐴_1,1=𝐴_𝑛,𝑛=1, 𝑧_𝑖=ℎ_𝑖+2∑^𝑚_𝑘=𝑖+1ℎ_𝑘, 𝐴 has order 𝑛−2 and𝐴=𝑇(𝐡)−𝐻(𝜎²(𝐡),𝐽𝜎²(𝐡)). (15) According to the brief discussion of Section 2.2, relation (15) implies that 𝐴=𝜏_𝑛−2(ℎ) with ℎ(𝑦)=ℎ₀+2∑^𝑚_𝑘=1ℎ_𝑘cos(𝑘𝑦) (see (10) and (11)). Moreover in [28] it is proved that 𝐴=AR_𝑛(ℎ).

2.4. Eigenvalues and Eigenvectors of AR-BC Matrices

We first describe the spectrum of AR-BC matrices, under the usual mild degree condition (i.e., the PSF 𝐡 has finite support), with symmetric, normalized PSFs. Then we describe the eigenvector structure and we introduce the AR-transform.

The spectral structure of any AR-BC matrix, with symmetric PSF 𝐡, is concisely described in the following result.

Theorem 1 (see [28]). Let the blurring function (PSF) 𝐡 be symmetric (i.e., ℎ_𝑠=ℎ_−𝑠), normalized, and satisfying the usual degree condition. As a consequence the eigenvalues of the 𝑛×𝑛 AR-BCs blurring matrix 𝐴 in (14), 𝑛≥3, are given by ℎ(0)=1 with multiplicity two and ℎ(𝐺_𝑛−2).

The proof can be easily derived by (12) which shows that the eigenvalues of AR_𝑛(ℎ) are ℎ(0) with multiplicity 2 and those of 𝜏_𝑛−2(ℎ), that is, ℎ(𝐺_𝑛−2), with multiplicity 1 each.

Here we will determine the eigenvectors of every matrix AR_𝑛(ℎ). In particular, we show that every AR-BCs matrix is diagonalizable, and we demonstrate independence of the eigenvectors from the symbol ℎ. With reference to the notation in (8)–(11), calling 𝐪_𝑗^(𝑛−2) the 𝑗 th column of 𝑄_𝑛−2, and 𝑦_𝑗^(𝑛−2) the 𝑗th point of 𝐺_𝑛−2, 𝑗=1,…,𝑛−2, we haveAR_𝑛0(ℎ)𝐪_𝑗^(𝑛−2)0=𝐯ℎ(0)_𝑛−2(ℎ)𝜏_𝑛−2(ℎ)𝐽𝐯_𝑛−20(ℎ)𝐪ℎ(0)0^{𝑗(𝑛−2)}=0𝜏_𝑛−2(ℎ)𝐪_𝑗^(𝑛−2)0=ℎ(𝑦_𝑗^(𝑛−2))0𝐪_𝑗^(𝑛−2)0, (16) since 𝐪_𝑗^(𝑛−2) is an eigenvector of 𝜏_𝑛−2(ℎ) and ℎ(𝑦_𝑗^(𝑛−2)) is the related eigenvalue. Due to the centrosymmetry of the involved matrix, if [1,𝐩^𝑇,0]^𝑇 is an eigenvector of AR_𝑛(ℎ) related to the eigenvalue ℎ(0), then the other is its flip, that is, [0,(𝐽𝐩)^𝑇,1]^𝑇. Let us look for this eigenvector by imposing the equalityAR_𝑛1𝐩01𝐩0(ℎ)=ℎ(0), (17) which is equivalent to seeking a vector 𝐩 that satisfies𝐯_𝑛−2(ℎ)+𝜏_𝑛−2(ℎ)𝐩=ℎ(0)𝐩. (18) Since 𝐯_𝑛−2(ℎ)=𝜏_𝑛−2(𝜙(ℎ))𝐞₁ by definition of the operator 𝐯_𝑛−2(⋅) (see (12) and the lines below), and because of the algebra structure of 𝜏_𝑛−2 and thanks to the above relation, we deduce that the vector 𝐩 satisfies the relation𝜏_𝑛−2(ℎ−ℎ(0))[−𝐿⁻¹_𝑛−2𝐞₁+𝐩]=0, (19) where 𝐿_𝑛−2 is the discrete one-level Laplacian, that is, 𝐿_𝑛−2=𝜏_𝑛−2(2−2cos(⋅)). Therefore, by (19), the solution is given by 𝐩=𝐿⁻¹_𝑛−2𝐞₁. If 𝜏_𝑛−2(ℎ−ℎ(0)) is invertible (as it happens for every nontrivial PSF, since the unique maximum of the function is reached at 𝑦=0, which is not a grid point of 𝐺_𝑛−2), then the solution is unique. Hence, independently of ℎ, we haveAR_𝑛1(ℎ)𝐩𝑄_𝑛−21=1𝐽𝐩𝐩𝑄_𝑛−21𝐽𝐩ℎ(0)diag(ℎ(𝐺_𝑛−2)).ℎ(0) (20) Now we observe that the 𝑗th eigenvector is unitary, 𝑗=2,…,𝑛−1, because 𝑄_𝑛−2 is unitary: we wish to impose the same condition on the first and the last eigenvector. The interesting fact is that 𝐩 has an explicit expression. By using standard finite difference techniques, it follows that 𝑝_𝑗=1−𝑗/(𝑛−1) so that the first eigenvector is exactly the sampling of the function 1−𝑥 on the grid 𝑗/(𝑛−1) for 𝑗=0,…,𝑛−1. Its Euclidean norm is 𝛼_𝑛=∑^𝑛−1_𝑗=0𝑗²/(𝑛−1)∼𝑛/3, where, for nonnegative sequences 𝛽_𝑛, 𝛾_𝑛, the relation 𝛾_𝑛∼𝛽_𝑛 means 𝛾_𝑛=𝛽_𝑛(1+𝑜(1)). In this way, the (normalized) AR-transform can be defined as 𝑇_𝑛=𝛼_𝑛⁻¹𝛼_𝑛⁻¹𝐏𝑄_𝑛−2𝛼_𝑛⁻¹𝛼𝐽𝐏_𝑛−1. (21)

Remark 2. With the normalization condition in (21), all the columns of 𝑇_𝑛 are unitary. However orthogonality is only partially fulfilled since it holds for the central columns, while the first and last columns are not orthogonal to each other, and neither one is orthogonal to the central columns. We can solve the first problem: the sum of the first and of the last column (suitably normalized) and the difference of the first and the last column (suitably normalized) become orthonormal, and are still eigenvectors related to the eigenvalue ℎ(0). However, since 𝐪₁^(𝑛−2) has only positive components and the vector space generated by the first and the last column of 𝑇_𝑛 contains positive vectors, it follows that 𝑇_𝑛 cannot be made orthonormal just by operating on the first and the last column. Indeed, we do not want to change the central block of 𝑇_𝑛 since it is related to a fast 𝑂(𝑛log(𝑛)) real transform and hence, necessarily, we cannot get rid of this quite mild lack of orthogonality.

Remark 3. There is a suggestive functional interpretation of the transform 𝑇_𝑛. When considering periodic BCs, the transform of the related matrices is the Fourier transform: its 𝑗th column vector, up to a normalizing scalar factor, can be viewed as a sampling, over a suitable uniform gridding of [0,2𝜋], of the frequency function exp(−𝐢𝑗𝑦). Analogously, when imposing reflective BCs with a symmetric PSF, the transform of the related matrices is the cosine transform: its 𝑗 th column vector, up to a normalizing scalar factor, can be viewed as a sampling, over a suitable uniform gridding of [0,𝜋], of the frequency function cos(𝑗𝑦). Here the imposition of the antireflective BCs can be functionally interpreted as a linear combination of sine functions and of linear polynomials (whose use is exactly required for imposing 𝐶¹ continuity at the borders).

The previous observation becomes evident in the expression of 𝑇_𝑛 in (21). Indeed, by defining the one-dimensional grid 𝐺_𝑛=[0,𝐺^𝑇_𝑛−2,𝜋]^𝑇=[𝑗𝜋/(𝑛−1)]^𝑛−1_𝑗=0, which is a subset of [0,𝜋], we infer that the first column of 𝑇_𝑛 is given by 𝛼_𝑛⁻¹𝐺(1−𝑦/𝜋)|_𝑛, the 𝑗th column of 𝑇_𝑛, 𝑗=2,…,𝑛−1, is given by 𝐺2/(𝑛−1)(sin(𝑗𝑦))|_𝑛, and finally the last column of 𝑇_𝑛 is given by 𝛼_𝑛⁻¹𝐺(𝑦/𝜋)|_𝑛, that is,𝑇_𝑛𝑦=[1−𝜋𝑦,sin(𝑦),…,sin((𝑛−2)𝑦),𝜋𝐺]|_𝑛⋅Δ_𝑛,Δ_𝑛=diag(𝛼_𝑛⁻¹,2𝐼𝑛−1^𝑛−2,𝛼_𝑛⁻¹). (22)

Finally, it is worth mentioning that the inverse transform is also described in terms of the same block structure since 𝑇_𝑛⁻¹=𝛼_𝑛−𝑄_𝑛−2𝐏𝑄_𝑛−2−𝑄_𝑛−2𝛼𝐽𝐏_𝑛 (23)

Theorem 4 (𝐴𝑅_𝑛(⋅) Jordan Canonical Form). With the notation and assumptions of Theorem 1, the 𝑛×𝑛 AR-BCs blurring matrix 𝐴 in (14), 𝑛≥3, coincides with𝐴𝑅_𝑛(ℎ)=𝑇_𝑛diag(ℎ(𝐺_𝑛))𝑇_𝑛⁻¹, (24) where 𝑇_𝑛 and 𝑇_𝑛⁻¹ are defined in (22) and (23), while 𝐺_𝑛=[0,𝐺^𝑇_𝑛−2,0]^𝑇.

2.5. Multilevel Extension

Here we provide some comments on the extension of our findings to 𝑑-dimensional objects with 𝑑>1. When 𝑑=1, 𝐡 is a vector, when 𝑑=2, 𝐡 is a 2D array, when 𝑑=3, 𝐡 is a 3D tensor and so forth.

With reference to Section 2.3 we propose a (canonical) multidimensional extension of the algebras 𝒜ℛ and of the operators AR_𝑛(⋅), 𝑛=(𝑛₁,…,𝑛_𝑑): the idea is to use tensor products. If ℎ=ℎ(⋅) is 𝑑-variate real-valued cosine polynomial, then its Fourier coefficients form a real 𝑑-dimensional tensor which is strongly symmetric, that is symmetric with respect to every direction, that is, ℎ_𝑗=ℎ_|𝑗| for all 𝑗∈ℤ^𝑑. In addition, ℎ(𝐲), 𝐲=(𝑦₁,…,𝑦_𝑑), can be written as a linear combination of independent terms of the form 𝑚(𝐲)=∏^𝑑_𝑗=1cos(𝛼_𝑗𝑦_𝑗) where any 𝛼_𝑗 is a nonnegative integer. Therefore, we defineAR_𝑛(𝑚(𝐲))=AR_𝑛1(cos(𝛼₁𝑦₁))⊗⋯⊗AR_𝑛𝑑(cos(𝛼_𝑑𝑦_𝑑)), (25) where ⊗ denotes Kronecker product, and we forceAR_𝑛(𝛼ℎ₁+𝛽ℎ₂)=𝛼AR_𝑛(ℎ₁)+𝛽AR_𝑛(ℎ₂) (26) for every real 𝛼 and 𝛽 and for every 𝑑-variate real-valued cosine polynomials ℎ₁=ℎ₁(⋅) and ℎ₂=ℎ₂(⋅). It is clear that the request that AR_𝑛(⋅) is a linear operator (for 𝑑>1, we impose this property in (26) by definition) is sufficient for defining completely the operator in the 𝑑-dimensional setting.

With the above definition of the operator AR_𝑛(⋅), we have (1)𝛼AR_𝑛(ℎ₁)+𝛽AR_𝑛(ℎ₂)=AR_𝑛(𝛼ℎ₁+𝛽ℎ₂), (2)AR_𝑛(ℎ₁)AR_𝑛(ℎ₂)=AR_𝑛(ℎ₁ℎ₂),

for real 𝛼 and 𝛽 and for cosine functions ℎ₁=ℎ₁(⋅) and ℎ₂=ℎ₂(⋅).

The latter properties of algebra homomorphism allows to define a commutative algebra 𝒜ℛ of the matrices AR_𝑛(ℎ), with ℎ(⋅) being a 𝑑-variate cosine polynomial. By standard interpolation arguments it is easy to see that 𝒜ℛ can be defined as the set of matrices AR_𝑛(ℎ), where ℎ is a 𝑑-variate cosine polynomial of degree at most 𝑛_𝑗−3 in the 𝑗th variable for every 𝑗 ranging in {1,…,𝑑}: we denote the latter polynomial set by 𝒫^(𝑑,even)_𝑛−2𝑒, with 𝑒 being the vector of all ones. Here we have to be a bit careful in specifying the meaning of algebra when talking of polynomials. More precisely, for ℎ₁,ℎ₂∈𝒫^(𝑑,even)_𝑛−2𝑒 the product ℎ₁⋅ℎ₂ is the unique polynomial ℎ∈𝒫^(𝑑,even)_𝑛−2𝑒 satisfying the following interpolation condition:ℎ(𝑦)=𝑧_𝑦,𝑧_𝑦≡ℎ₁(𝑦)ℎ₂(𝑦),∀𝑦∈𝐺^(𝑑)_𝑛−2. (27) If the degree of ℎ₁ plus the degree of ℎ₂ in the 𝑗 th variable does not exceed 𝑛_𝑗−2, 𝑗=1,…,𝑑, then the uniqueness of the interpolant implies that ℎ coincides with the product between polynomials in the usual sense. The uniqueness holds also for 𝑑≥2 thanks to the tensor form of the grid 𝐺^(𝑑)_𝑛−2 (see [28] for more details). The very same idea applies when considering inversion. In conclusion, with this careful definition of the product/inversion and with the standard definition of addition, 𝒫^(𝑑,even)_𝑛−2𝑒 has become an algebra, showing the vector-space dimension equal to (𝑛₁−2)⋅(𝑛₂−2)⋯(𝑛_𝑑−2) which coincides with that of 𝒜ℛ_𝑛.

Without loss of generality and for the sake of notational clarity, in the following we assume 𝑛_𝑗=𝑛 for 𝑗=1,…,𝑑. Thanks to the tensor structure emphasized in (25)–(26), and by using Theorem 4 for every term AR_𝑛(cos(𝛼_𝑗𝑦_𝑗)), 𝑗=1,…,𝑑, of AR_𝑛(𝑚) the 𝑑-level extension of such a theorem easily follows. More precisely, if ℎ is a 𝑑-variate real-valued cosine symbol related to a 𝑑-dimensional strongly symmetric and normalized mask 𝐡, thenAR_𝑛(ℎ)=𝑇_𝑛^(𝑑)𝐷_𝑛(𝑇_𝑛^(𝑑))⁻¹,𝑇_𝑛^(𝑑)=𝑇_𝑛⊗⋯⊗𝑇_𝑛, (28) (𝑑 times) where 𝐷_𝑛 is the diagonal matrix containing the eigenvalues of AR_𝑛(ℎ). The description of 𝐷_𝑛 in 𝑑 dimensions is quite involved when compared with the case 𝑑=1, implicitly reported in Theorem 1.

For a complete analysis of the spectrum of AR_𝑛(ℎ) we refer the reader to [28]. Here we give details on a specific aspect. More precisely we attribute a correspondence in a precise and simple way among eigenvalues and eigenvectors, by making recourse only to the main 𝑑-variate symbol ℎ(⋅). Let 𝐱_𝑛=𝐱_𝑛⁽¹⁾⊗𝐱_𝑛⁽²⁾⊗⋯⊗𝐱_𝑛^(𝑑) be a column of 𝑇_𝑛^(𝑑), with 𝐱_𝑛^(𝑗)∈{𝛼_𝑛⁻¹[1,𝐩^𝑇,0]^𝑇,𝛼_𝑛⁻¹[0,(𝐽𝐩)^𝑇,1]^𝑇} or 𝐱_𝑛^(𝑗)=[0,𝐪^𝑇_𝑠𝑗,0]^𝑇, 1≤𝑠_𝑗≤𝑛−2 and 𝐪_𝑠𝑗 is the (𝑠_𝑗) th column of 𝑄_𝑛−2, for 𝑗=1,…,𝑑. Let ℱ_𝐱𝑛={𝑗∣𝐱_𝑛^(𝑗)=𝛼_𝑛⁻¹[1,𝐩^𝑇,0]^𝑇or𝐱_𝑛^(𝑗)=𝛼_𝑛⁻¹[0,(𝐽𝐩)^𝑇,1]^𝑇}⊂{1,…,𝑑}, (29) with 𝐱_𝑛 being the generic eigenvector, that is, the generic column of 𝑇_𝑛^(𝑑). The eigenvalue related to 𝐱_𝑛 is𝜆=ℎ(𝑦₁^(𝑛),…,𝑦_𝑑^(𝑛)), (30) where 𝑦_𝑗^(𝑛)=0 for 𝑗∈ℱ_𝐱𝑛 and 𝑦_𝑗^(𝑛)=𝜋𝑣_𝑗/𝑛−1 for 𝑗∉ℱ_𝐱𝑛. Defining the 𝑑-dimensional gridG_𝑛^(𝑑)=𝐺_𝑛⊗⋯⊗𝐺_𝑛,𝑑times, (31) we can evaluate the 𝑑-variate function ℎ on 𝐺_𝑛^(𝑑) by ℎ(reshape(𝐺_𝑛^(𝑑),𝑛)), wherereshape(𝑋,𝑛) (32) arranges the entries of 𝑋 in a 𝑑-dimensional array of length 𝑛 along each direction according to Matlab notation. Using this notation the following compact and elegant result can be stated (its proof is omitted since it is simply the combination of the eigenvalue analysis in [28], of Theorem 4, and of the previous tensor arguments).

Theorem 5 (𝐴𝑅_𝑛(⋅) Jordan Canonical Form). The 𝑛^𝑑×𝑛^𝑑 AR-BCs blurring matrix 𝐴, obtained when using a strongly symmetric 𝑑-dimensional mask ℎ such that ℎ_𝑖=0 if |𝑖_𝑗|≥𝑛−2 for some 𝑗∈{1,…,𝑑} (the 𝑑-dimensional degree condition), 𝑛≥3, coincides with AR_𝑛(ℎ)=𝑇_𝑛^(𝑑)diag(reshape(ℎ(reshape(𝐺_𝑛^(𝑑),𝑛)),𝑛^𝑑))(𝑇_𝑛^(𝑑))⁻¹, (33) where 𝑇_𝑛^(𝑑) and 𝐺_𝑛^(𝑑) are defined in (28) and (31).

It is worth observing that the matrix AR_𝑛(ℎ) spectrally analyzed in the previous theorem is exactly the same matrix arising from the imposition of AR-BCs applied separately in every direction, when ℎ is the multivariate cosine symbol coming from the 𝑑-D tensor mask 𝐡 defining the shift-invariant 𝑑-dimensional blurring operator.

3. Regularization by Reblurring

When the observed signal (or image) is noise-free, then there is a substantial gain of the reflective boundary conditions (R-BCs) with respect to both the periodic and zero Dirichlet BCs and, at the same time, there is a significant improvement of the AR-BCs with regard to the R-BCs (see [17, 30]). Since the deconvolution problem is ill posed (components of the solution related to high frequency data errors are greatly amplified) regardless of the chosen BCs, it is evident that we have to regularize the problem. Two classical methods, that is, Tikhonov regularization, with direct or iterative solution of the Tikhonov linear system, and regularization iterative solvers, with early termination, for normal equations (CG [31] or Landweber method [32]) can be used. We observe that in both the cases, the coefficient matrix involves a shift of 𝐴^𝑇𝐴 and that the righthand-side is given by 𝐴^𝑇𝐠. Quite surprisingly, the AR-BCs may be spoiled by this approach at least for 𝑑=1 and if we compute explicitly the matrix 𝐴^𝑇𝐴 and the vector 𝐴^𝑇𝐠, see [33]: more in detail, even in presence of a moderate noise and a strongly symmetric PSF, the accuracy of AR-BCs restorations becomes worse in some examples than the accuracy of R-BCs restorations (see [33]). The reason of this fact relies upon the properties of the matrix 𝐴^𝑇, and we give some insights in the following.

3.1. The Reblurring Operator

A key point is that, for zero Dirichlet, periodic and reflective BCs, when the kernel ℎ is symmetric, the matrix 𝐴^𝑇 is still a blurring operator since 𝐴^𝑇=𝐴, while, in the case of the AR-BCs matrix, 𝐴^𝑇 cannot be interpreted as a blurring operator. A (normalized) blurring operator is characterized by nonnegative coefficients such that every row sum is equal to 1 (and it is still acceptable if it is not higher than 1): in the case of 𝐴^𝑇 with AR-BCs the row sum of the first and of the last row can be substantially larger than 1. This means that modified signal 𝐴^𝑇𝐠 may have artifacts at the borders and this seems to be a potential motivation for which a reduction of the reconstruction quality occurs.

Furthermore, the structure of the matrix 𝐴^𝑇𝐴 is also spoiled and, in the case of images (𝑑=2) we lose the 𝑂(𝑛²log(𝑛)) computational cost for solving a generic system 𝐴^𝑇𝐴𝐱=𝐛. The cost of solving such a type of linear systems is proportional to 𝑛³ by using for example, Shermann-Morrison formulae (which by the way can be numerically unstable [34]). Dealing with higher dimensions, the scenario is even worse [35], because in the 𝑑-dimensional setting the solution of the normal equation linear system is asymptotic to 𝑛^3(𝑑−1), where 𝑛^𝑑 is the size of the matrix 𝐴. In order to overcome these problems (which arise only with the most precise AR-BCs for strongly symmetric PSFs), we replace 𝐴^𝑇 by 𝐴^′, where 𝐴^′ is the matrix obtained by imposing the current BCs to the center-flipped PSF (i.e., in the 2D case, to the PSF rotated by 180 degrees).

For the sake of simplicity we first consider a strongly symmetric PSF so that the associated normal equations can be read as 𝐴²𝐟=𝐴𝐠, whenever zero Dirichlet, periodic or reflective BCs are considered. Therefore, the observed image 𝐠 is reblurred. The reblurring is the key of the success of classical regularization techniques (Tikhonov or CG, Landweber for the normal equations) since also the noise is blurred and this makes the contribution of the noise less evident. We notice that the two systems 𝐴²𝐟=𝐴𝐠 and 𝐴𝐟=𝐠 are algebraically equivalent if 𝐴 is invertible: in any case, if 𝐴 is also positive definite, the first represents the minimization of the functional ‖𝐴𝐟−𝐠‖²₂ while the second represents the minimization of the functional ‖𝐴^1/2𝐟−𝐴^−1/2𝐠‖²₂ so that the first can be considered the blurred version of the second and in fact the first approach is uniformly better than the second. On these grounds, in the case of AR-BCs, since 𝐴^𝑇≠𝐴, we can replace 𝐴^𝑇 by 𝐴 which is again a low-pass filter (see [33]). In this way, we overcome also the computational problems.

In order to provide a general reblurring approach also in the case of nonsymmetric PSFs, we consider the correlation operation instead of the transposition (see [36]). In the discrete 2D case, the correlation performs the same operation of the convolution, but rotating the mask (the PSF in our case) of 180 degrees. Moreover, we note that in the continuous case over an infinite domain, the correlation and the adjoint are exactly the same operation, provided that the convolution kernel is real. Indeed, let 𝐾 be the convolution operator with shift-invariant kernel 𝑘(𝑠), then [𝐾𝑓](𝑥)=∫𝑘(𝑥−𝑦)𝑓(𝑦)d𝑦. Since the PSF (and then 𝑘) is real (and then real valued), the adjoint operator 𝐾^∗ is [𝐾^∗𝑓](𝑥)=∫𝑘(𝑦−𝑥)𝑓(𝑦)d𝑦 which is a correlation operator. We remark that here the convolution and the correlation use the same kernel except for the sign of the variable (i.e., 𝑘(𝐬) vs 𝑘(−𝐬)), and, in the 2D case, the change of sign in the variable 𝐬 of the kernel can be viewed as a 180 degrees rotation of the PSF mask.

By virtue of these arguments, in order to overcome the problems arising with the normal equations approach for AR-BCs 2D restorations, we propose to replace 𝐴^𝑇 by 𝐴^′ (the reblurring matrix), where 𝐴^′ is the matrix obtained by imposing the BCs to the PSF rotated by 180 degrees. Using Matlab notation, if 𝐻 is a 𝑞×𝑞 PSF, its 180 degrees rotated version is 𝐻^′=fliplr(flipud(𝐻))=𝐽_𝑞𝐻𝐽_𝑞, where 𝐽_𝑞 is the flip matrix defined as [𝐽_𝑞]_𝑖,𝑗=1 if 𝑖+𝑗=𝑞+1 for 𝑖,𝑗=1,…,𝑞, and zero otherwise. For a 𝑑-dimensional problem, 𝐴^′ is obtained by imposing the BCs to the PSF flipped with respect to the origin, or, in other words, to the new PSF where all the coefficients are flipped with respect to every variable.

In this way 𝐴^′ has the same computational properties of 𝐴 (it belongs to 𝒜ℛ in the case of AR BCs) and it is certainly a low-pass filter. In the reblurring approach the normal equations are replaced by𝐴^′𝐴𝐟=𝐴^′𝐠. (34) Furthermore, using the Tikhonov regularization in the reblurred version, we propose to use[𝐴^′𝐴+𝜇𝐿^′𝐿]𝐟=𝐴^′𝐠, (35) with 𝐿 being any discrete. differential operator with AR-BCs. In general 𝐴^′𝐴 is nonsymmetric, but the asymptotic eigenvalue analysis in [28, Section 3.3] has shown that the spectrum is clustered around the positive real interval given by the range of |ℎ|² if ℎ is the symbol associated to the PSF. Such a cluster is strong if the decay of the PSF coefficients is fast enough, as it occurs in real world PSFs. The latter statements give a reasonable support to the applicability of the CGLS when 𝐴^𝑇 is replaced by 𝐴^′ and, indeed, in our numerical tests such a regularizing method has always worked perfectly. From the viewpoint of the modeler, the previous considerations can be summarized in the following motivation. The image restoration problem is the restriction in the FOV of an infinite dimensional problem. We can follow two ways to design the linear system to solve: (1)to impose BCs and then to look at a least-squares solution, (2)to formulate a least-squares solution on the infinite dimensional problem, and then to impose the BCs to the two infinite dimensional operators 𝐾 and 𝐾^∗, separately.

A third possibility is to formulate a least-squares solution on the infinite dimensional problem, and then to impose the BCs to this minimum problem: a difficulty in this case is that, even without noise, the resulting system is not equivalent in an algebraic sense to the original equations 𝐴𝐟=𝐠. In the first case we resort to the normal equations in the finite dimensional space. On the contrary, in the second case we apply the BCs to 𝐾 and 𝐾^∗ in the infinite dimensional normal equations (where the adjoint operator 𝐾^∗ performs a correlation operation) and then we obtain (34). More precisely, the discretization of 𝐾 and 𝐾^∗ in the continuous equation 𝐾^∗𝐾𝑓=𝐾^∗𝑔 with any fixed BCs gives (34).

3.2. Linear Algebra and Computational Issues

We note that in the 1D case 𝐴^′_𝑛=𝐽_𝑛𝐴_𝑛𝐽_𝑛. In the 𝑑-dimensional case, let 𝑛=(𝑛₁,𝑛₂,…,𝑛_𝑑) be the partial dimensions of the problem, whose total size is ∏^𝑑_𝑖=1𝑛_𝑖. By flipping each variable, we obtain𝐴^′_𝑛=𝐽_𝑛𝐴_𝑛𝐽_𝑛,𝐽_𝑛=^𝑑_𝑖=1𝐽_𝑛𝑖. (36) For the analysis of properties of the reblurring scheme (34) with respect to all the different BCs, we now study the discretization of the continuous operator 𝐾. Let us consider the Toeplitz 𝑑-level matrix 𝑇_𝑛(𝜙) of partial dimensions 𝑛=(𝑛₁,…,𝑛_𝑑)∈ℕ^𝑑₊ and generating function 𝜙 [2], which is defined as 𝑇_𝑛=(𝜙)_{|𝑗|≤𝑛−𝑒}𝑎_𝑗𝑍_𝑛^[𝑗]=_{|𝑗1|<𝑛1}⋯_{|𝑗𝑑|<𝑛𝑑}𝑎_{(𝑗1,…,𝑗𝑑)}𝑍^[𝑗1]_𝑛1⊗⋯⊗𝑍^[𝑗𝑑]_𝑛𝑑 (37) (𝑒=(1,…,1)∈ℕ^𝑑₊) by means of the Fourier coefficients of 𝜙𝑎_𝑘=1(2𝜋)^{𝑑[−𝜋,𝜋]𝑑}𝜙(𝑥)exp^{−𝐢⟨𝑘∣𝑥⟩}dx,𝐢²=−1,𝑘∈ℤ^𝑑. (38) Here ⟨𝑘∣𝑥⟩=∑^𝑑_𝑖=1𝑘_𝑖𝑥_𝑖 and for 𝑗∈ℤ, 𝑚∈𝑁₊, 𝑍_𝑚^[𝑗]∈ℝ^𝑚×𝑚 is the matrix whose (𝑠,𝑟) th entry is 1 if 𝑠−𝑟=𝑗, and 0 elsewhere. As it is well known for multilevel Toeplitz matrices 𝑇^𝐻_𝑛(𝜙)=𝑇_𝑛(𝜙), where 𝜙 is the conjugate of the function 𝜙, and the Fourier coefficients of 𝜙 are the same of 𝜙, but conjugated and flipped. Moreover, since 𝑇_𝑛(𝜙)=𝐽_𝑛𝑇_𝑛(𝜙)𝐽_𝑛, if 𝐴_𝑛=𝑇_𝑛(𝜙) is real then 𝐴^𝑇_𝑛=𝐽_𝑛𝐴_𝑛𝐽_𝑛=𝐴^′_𝑛. This means that for Dirichlet BCs (D-BCs) and periodic BCs (P-BCs) the reblurring approach is exactly equal to the classical normal equations approach, since in these two cases the corresponding blurring matrix 𝐴_𝑛 is multilevel Toeplitz: indeed, concerning P-BCs, we notice that the resulting multilevel circulant structure is a special instance of the multilevel Toeplitz case. Unfortunately, in the case of Hankel matrices (or multilevel mixed Toeplitz-Hankel matrices with at least one Hankel level) this is no longer true in general. However, a sufficient and necessary condition to have 𝐴^′_𝑛=𝐴^𝑇_𝑛 is 𝐽_𝑛𝐴_𝑛=(𝐽_𝑛𝐴_𝑛)^𝑇 (or equivalently 𝐴_𝑛𝐽_𝑛=(𝐴_𝑛𝐽_𝑛)^𝑇), which is a multilevel antidiagonal symmetry called persymmetry. Therefore in the case of R-BCs, where the matrix 𝐴_𝑛 involves nested Hankel parts, in general 𝐴^′_𝑛≠𝐴^𝑇_𝑛, while 𝐴^′_𝑛=𝐴^𝑇_𝑛 only when the PSF is strongly symmetric since in this case 𝐽_𝑛𝐴_𝑛=(𝐽_𝑛𝐴_𝑛)^𝑇. Dealing with the AR-BCs, the situation is even more involved, since 𝐴^′_𝑛≠𝐴^𝑇_𝑛 also for strongly symmetric PSFs, owing to the low-rank correction term. Hence, we can state that the reblurring is a new proposal not only for the AR-BCs, but also for all those BCs for which 𝐽_𝑛𝐴_𝑛≠(𝐽_𝑛𝐴_𝑛)^𝑇. As a nontrivial and unexpected example, it is important to stress that the imposition of R-BCs with nonstrong symmetric PSFs implies 𝐽_𝑛𝐴_𝑛≠(𝐽_𝑛𝐴_𝑛)^𝑇, that is, 𝐴^′_𝑛≠𝐴^𝑇_𝑛.

We provide now a computational motivation for the choice of using 𝐴^′ as an alternative to 𝐴^𝑇∶𝐴^′ is the usual operation which has to be implemented to perform the adjoint operation in the Fourier domain. Indeed, the convolution with prescribed BCs can be implemented by first enlarging the image according to the considered BCs and then by computing the matrix vector product by a simple circular convolution operation, see [37]. More precisely, let 𝑋 and 𝐻 be two matrices such that 𝑋 represents an 𝑛×𝑛 image and 𝐻 is the discrete 𝑞×𝑞 2D PSF, 𝑞=2𝑚+1, which leads to the matrix blurring 𝐴. By using the Matlab notation 𝐱=𝑋(∶) (i.e., the vector 𝐱 is the column-stacked version of 𝑋), the product 𝐴𝐱 in the 2D case can be implemented (1)by using an enlarged image 𝑋, which is the (𝑛+𝑞−1)×(𝑛+𝑞−1) image 𝑋 extended at the boundaries according to the imposed BCs, (2)computing 𝐻∗𝑋, where the symbol “∗” denotes the circular convolution operator (𝐻 should be zero padded to have the same size of 𝑋), (3)and then taking the inner part of the result having the same size of 𝑋.

The circular convolution can be computed using the 2D discrete Fourier transform (DFT2) and its inverse (IDFT2), since we have𝐻∗𝑋=IDFT2(DFT2(𝐻)⊙DFT2(𝑋)), (39) where “⊙” is the componentwise product. If 𝒞_𝑘(𝑓) denotes the block circulant matrix with circulant blocks, of block size 𝑘 with blocks of size 𝑘 and generating function 𝑓, then (39) represents the same operation as 𝒞_{𝑛+𝑞−1}(𝜙)𝐱, where 𝐱=𝑋(∶) (according to a 2D ordering). Conversely, it is well known that the operation corresponding to the product with the adjoint operator, in the Fourier domain gives rise to𝑌=𝐻^′∗𝑋=IDFT2(𝐻DFT2(𝐻)⊙DFT2(𝑋)),_′𝑖,𝑗=𝐻_{−𝑖,−𝑗},𝑖,𝑗∈ℤ, (40) where the overline symbol denotes the complex conjugation. As a result, since the transform 𝐻↦𝐻^′ is equivalent to the transform 𝜙↦𝜙 (because 𝜙(𝑥,𝑦)=∑_{𝑖,𝑗∈ℤ}𝐻_𝑖,𝑗exp^{𝐢⟨(𝑖,𝑗)∣(𝑥,𝑦)⟩}), and since 𝒞_{𝑛+𝑞−1}(𝜙)=𝒞^𝑇_{𝑛+𝑞−1}(𝜙), if 𝐲=𝑌(∶), then it follows that 𝐲=𝒞^𝑇_{𝑛+𝑞−1}(𝜙)𝐱. Therefore, for any of the considered BCs, the inner part of 𝐲 is exactly 𝐴^′𝐱. Here it is worthwhile to specify exactly what we mean for inner part: if the vector 𝐲 is partitioned in 𝑛+𝑞−1 blocks of size 𝑛+𝑞−1, 𝑞=2𝑚+1, then for inner part we mean that we are excluding the first and the last 𝑚 blocks and, in any of the remaining blocks, we are deleting the first and the last 𝑚 entries. More generally, if the PSF is arbitrary (e.g. nonsymmetric) that is, the nonzero coefficients of the PSF have first index belonging to [−𝑚⁻₁,𝑚⁺₁] and second index in the range [−𝑚⁻₂,𝑚⁺₂], then we have to delete the first 𝑚⁻₁ and the last 𝑚⁺₁ blocks and, in any of the other blocks, we have to exclude the first 𝑚⁻₂ and the last 𝑚⁺₂ entries.

Since the DFT and its inverse can be computed in 𝑂(𝑛²log(𝑛)) arithmetic operations using FFTs, the previous approach is implemented in the Matlab toolbox RestoreTools [37]. We have added the AR-BCs in such a toolbox for the matrix vector product, suitable for iterative regularizing methods. This code has been used for the numerical tests of Section 3 and it is downloadable from the homepage “http://scienze-como.uninsubria.it/mdonatelli/software.html”.

3.3. Filtering Methods for AR-BCs Matrices

As mentioned in Section 1, regardless of the imposed BCs, matrices 𝐴 that arise in signal and image restoration are typically severely ill conditioned, and regularization is needed in order to compute a stable approximation of the solution of (1). A class of regularization methods is obtained through spectral filtering [38, 39]. Specifically, if the spectral decomposition of 𝐴 is 𝐴=𝑇_𝑛diag(𝐝)𝑇_𝑛⁻¹,𝑇_𝑛=𝐭₁𝐭₂⋯𝐭_𝑛,𝑇_𝑛⁻¹=𝐭^𝑇₁̃𝐭^𝑇₂̃⋮𝐭^𝑇_𝑛̃, (41) with 𝐝=ℎ(𝐺_𝑛), then a spectral filter solution is given by𝐟_reg=^𝑛_𝑖=1𝜙_𝑖𝐭^𝑇_𝑖̃𝐠𝑑_𝑖𝐭_𝑖, (42) where 𝜙_𝑖 are filter factors that satisfy 𝜙_𝑖≈1,if𝑑_𝑖islarge,0,if𝑑_𝑖issmall. (43) The small eigenvalues correspond to eigenvectors with high frequency components, and are typically associated with the noise space, while the large eigenvalues correspond to eigenvectors with low-frequency components, and are associated with the signal space. Thus filtering methods attempt to reconstruct signal space components of the solution, while avoiding reconstruction of noise space components.

For example, the filter factors for two well-known filtering methods, truncated spectral value decomposition (TSVD) and Tikhonov regularization, are𝜙^tsvd_𝑖=1,if𝑑_𝑖≥𝛿,0,if𝑑_𝑖𝜙<𝛿,^tik𝑖=𝑑²_𝑖𝑑²_𝑖+𝜆,𝜆>0, (44) where the problem dependent regularization parameters 𝛿 and 𝜆 must be chosen [39]. Several techniques can be used to estimate appropriate choices for the regularization parameters when the SVD is used for filtering (i.e., 𝑑_𝑖 are the singular values), including generalized cross validation (GCV), L-curve, and the discrepancy principle [38, 40, 41].

In our case, the notation in (44) defines a slight abuse of notation, because the eigenvalues 𝑑_𝑖 are not the singular values: in fact the Jordan canonical form (CF) in (24) is different from the singular value decomposition (SVD), since the transform 𝑇_𝑛 is not orthogonal (indeed it is a rank-2 correction of a symmetric orthogonal matrix). Therefore, note that the use of 𝜙^tsvd_𝑖 in (42) defines the filtering of the eigenvalues in the Jordan CF instead of the more classical filtering of the singular values in the SVD. However, we note that in general computing the SVD can be computationally very expensive, especially in the multidimensional case and also in the strongly symmetric case. Moreover, quite surprisingly, a recent and quite exhaustive set of numerical tests, both in the case of signals and images (see [42]), has shown that the truncated Jordan CF is more or less equivalent to the truncated SVD in terms of quality of the restored object: indeed this is a delicate issue that deserves more attention in the future.

Our final aim is to compute (42) in a fast and stable way. This is the classic approach implemented for instance with periodic BCs by using three FFTs. In our case we employ the AR-transform (21), its inverse (23), and a fast algorithm for computing the eigenvalues described in [28].

Algorithm 6. (1) 𝐠=𝑇_𝑛⁻¹𝐠.
(2) 𝐝=[ℎ(0),𝐝^𝑇,ℎ(0)]^𝑇, where 𝐝=[𝑑₂,…,𝑑_𝑛−1]^𝑇 are the eigenvalues of 𝜏_𝑛−2(ℎ) that can be computed by a fast sine transform (FST).
(3) 𝐟̃=(𝜙./𝐝)⨀𝐠, where the dot operations are component-wise.
(4) 𝐟_reg=𝑇_𝑛𝐟.̃

The product 𝑇_𝑛𝐟̃ can be clearly computed in a fast and stable way by one FST. Indeed for all 𝐱∈ℝ^𝑛𝑇_𝑛𝐱=𝛼_𝑛⁻¹𝑥₁1𝐩0+0𝑄_𝑛−20𝐱(2∶𝑛−1)+𝛼_𝑛⁻¹𝑥_𝑛01𝐽𝐩, (45) where 𝐱(2∶𝑛−1) in Matlab notation is the vector 𝐱 with components indexed from 2 to 𝑛−1. A similar strategy can be followed for computing the matrix-vector product 𝑇_𝑛⁻¹𝐠. Instead of 𝛼_𝑛⁻¹𝐩 there is 𝐮=−𝑄_𝑛−2𝐩 and instead of 𝛼_𝑛⁻¹𝐽𝐩 there is 𝐰=−𝑄_𝑛−2𝐽𝐩. Recalling that 𝐩=𝐿⁻¹_𝑛−2𝐞₁ the two vectors 𝐮 and 𝐰 can be explicitly computed obtaining 𝑢_𝑖=(2𝑛−2)^−1/2cot(𝑖𝜋/(2𝑛−2)), for 𝑖=1,…,𝑛−2 and 𝐰=diag_{𝑖=1,…,𝑛−2}(−1)^𝑖+1𝐮.

Remark 7. As discussed in Remark 3, there is a natural interpretation in terms of frequencies when considering one-dimensional periodic and reflective BCs. The eigenvalue obtained as a sampling of the symbol ℎ at a grid-point close to zero, that is, close to the maximum point of ℎ, has an associated eigenvector that corresponds to low-frequency (signal space) information, while the eigenvalue obtained as a sampling of the symbol ℎ at a grid-point far away from zero (and, in particular, close to 𝜋), has an associated eigenvector that corresponds to high-frequency (noise space) information. Concerning antireflective BCs, the same situation occurs when dealing with the frequency eigenvectors 𝐺2/(𝑛−1)(sin(𝑗𝑦))|_𝑛, 𝑗=2,…,𝑛−1. The other two exceptional eigenvectors generate the space of linear polynomials and therefore they correspond to low-frequency information: this intuition is well supported by the fact that the related eigenvalue is ℎ(0), that is, the maximum and the infinity norm of ℎ, and by the fact that AR-BCs are more precise than other classical BCs.

Remark 8. For 𝑑=1 and with reference to the previous sections, we have proved that, thanks to the definition of a (fast) AR-transform, it is possible to define a truncated spectral decomposition which is computationally very effective and, surprisingly enough, quite competitive when compared with the celebrated but costly truncated SVD in terms of restoration quality. However, we are well aware that the real challenge is represented by a general extension to the multidimensional setting. Indeed, except for more involved multi-index notations, all the techniques are plainly generalized in the multilevel setting, maintaining a cost proportional to three 𝑑-level FSTs of size (𝑛−2)^𝑑, and the key tool is the simplified eigenvalue-eigenvector correspondence concisely indicated in Theorem 5. In reality, regarding the previous Algorithm 6 the only difficult task is the computation in step (2), where we have to compute the eigenvalues in the right order. For this task we refer to [28], where an algorithm is proposed and studied: more specifically the related procedure in [28] is based on a single 𝑑-level FST of size (𝑛−2)^𝑑 plus lower order computations.

3.4. Convergence of the Reblurring Approach

We remark that the antireflective transform can be defined also by the eigenvector matrix𝑉_𝑛=10𝑄0𝓘, (46) where 11=√𝑛11⋮1√,𝓘=3𝑛(𝑛²−1)⋮1−𝑛3−𝑛𝑛−1. (47)

Note that 1 and 𝓘 differ from the corresponding eigenvectors [1,𝐩^𝑇,0]^𝑇 and [0,𝐽𝐩^𝑇,1]^𝑇 used in Section 2.4, but they are a linear combination of them. They have been chosen here to form an orthonormal basis of the grid samples of all linear polynomials and this property will be useful in the following.

According to Theorem 4 the spectral decomposition of AR_𝑛(ℎ) can now be written asAR_𝑛(ℎ)=𝑉_𝑛Λ𝑉_𝑛⁻¹, (48) where the diagonal entries 𝜆_𝑗𝑗 of Λ are given by𝜆_𝑗𝑗=ℎ(𝑗−1𝑛−1𝜋),1≤𝑗<𝑛,ℎ(0),𝑗=𝑛. (49) Here we prove that the reblurring approach for Tikhonov regularization in (35) where 𝐿=𝐼 is a regularization method. For a complete analysis we need to compute the SVD of 𝑉_𝑛.

We use the notation 𝐲≐𝐳 if the two vectors 𝐲 and 𝐳 depend on 𝑛 and for each entry 𝑦_𝑖 of 𝐲 and the corresponding entry 𝑧_𝑖 of 𝐳 there holds 𝑦_𝑖/𝑧_𝑖→1 as 𝑛→∞.

Theorem 9 (see [43]). The two dominant singular values of 𝑉_𝑛 are given by 𝜎₁≐𝜎₂≐√2, (50) where √2 is, in fact, a strict upper bound, and the two minimal singular values are given by 𝜎_𝑛−1≐√3√𝑛and𝜎_𝑛≐1√𝑛, (51) respectively. Fix 1^′=1(2∶𝑛−1), 𝓘^′=𝓘(2∶𝑛−1), 𝐚=𝑄1^′, and 𝐛=𝑄𝓘^′. The corresponding right singular vectors are 𝐯₁≐1√21𝐚0,𝐯₂≐1√20𝐛1,𝐯_𝑛−1≐1√20𝐛−1,𝐯_𝑛≐1√2𝐚0,−1 (52) and the left singular vectors are 𝐮₁≐1(2√1𝑛)_′1(2√𝑛),𝐮₂≐𝑙₁2𝓘^′𝑙_𝑛2,𝐮_𝑛−1≐1√213𝑛𝓘^′−1,𝐮_𝑛≐1√21−1_′√𝑛,−1 (53) respectively. The remaining singular values are equal to one, and the corresponding left and right singular vectors have homogeneous boundary values.

Now we are in the position to determine the condition number of the antireflective transform to first order.

Corollary 10. The condition number of the antireflective transform satisfies 𝜇(𝑉_𝑛)≐√2𝑛,𝑛→∞. (54)

Remark 11. It is important to note that the ill-conditioned subspace of 𝑉_𝑛⁻¹ has dimension two, independent of 𝑛, since 𝑉_𝑛 has two singular values that decay like 1/√𝑛 while all others are between one and two. Also, 𝑉_𝑛⁻¹ only amplifies vectors that fail to be orthogonal to 𝒰=span{𝐮_𝑛−1,𝐮_𝑛}. According to Theorem 9 the vectors from 𝒰 are essentially zero, except for their two boundary values.

Convergent bounds for the Tikhonov reblurring approach, can be obtained with the usual remedy from the theory of ill-posed problems, which consists in so-called smoothness assumptions, the most simple one being as follows.

Assumption. Let 𝑓 be itself a blurred version of a continuous signal 𝑤, that is, 𝑓(𝑥)=_ℝ𝑘(𝑥−𝑦)𝑤(𝑦)𝑑𝑦,𝑥∈ℝ, (55) where 𝑤 satisfies the same BCs as 𝑓 (i.e., periodic, reflective, or antireflective ones).

On the grounds of Assumption 12 we may therefore assume that𝐟=𝐴𝐰,forsome𝐰∈ℝⁿ (56) with a moderate bound|‖𝐰‖|≤‖𝐰‖_∞≤𝜚, (57) where |‖𝐱‖|∶=‖𝐱‖₂/√𝑛,𝐱∈ℝ^𝑛. Since the observed object is usually affected by noise, instead of the blurred object 𝐠 we have to work with the blurred and noisy object 𝐠=𝐠+𝐞, where ‖𝐞‖_∞<𝜀. In this way, for 𝐿=𝐼 (35) becomes (𝐴^′𝐴+𝜇𝐼)𝐟=𝐴^′𝐠.

Theorem 13. Let the exact solution 𝐟 of (1) satisfy (56) with (57). Furthermore the total error of the reblurring strategy with AR BCs satisfies |‖𝐟^𝜀_𝛼−𝐟‖|=𝑂(√𝜀𝜚), (58) for 𝛼=𝛼(𝜀)=𝜀/𝜚, where the constant in the 𝑂(⋅)-notation is independent of the dimension 𝑛.

Note that the upper bound from Theorem 13 is the same as for Tikhonov regularization with reflective or periodic BCs; only the constant hidden in the 𝑂(⋅)-notation may be somewhat larger for the reblurring strategy.

4. Numerical Results

The section is devoted to numerical experiments, with reference to the Tikhonov regularization in the reblurred version and to the classical conjugate gradient (CG) regularization with early termination. In both cases the use of antireflective BCs improves the quality of the restored image, without penalizing the computational cost. Another promising approach not discussed here both from the viewpoint of the quality and of the complexity is that based on a regularized version of multigrid-type techniques (see [44, 45]): also this idea can be successfully implemented in combination with the choice of AR BCs.

In our numerical experiments we use Matlab 6.5 and the toolbox RestoreTools [37] suitably extended for dealing with AR-BCs. The relative restoration error (RRE) is ‖𝐟̂−𝐟‖₂/‖𝐟‖₂, where 𝐟̂ is the computed approximation of the true image 𝐟. The signal-to-noise ratio (SNR) is computed as 20log₁₀‖𝐠_𝑏‖₂/‖𝝂‖₂, where 𝐠_𝑏 is the blurred image without noise and 𝝂 is the noise vector [32].

4.1. Reblurring and Tikhonov Regularization

Let us begin with an example illustrating the approach discussed in Section 3.3 for a 2-dimensional imaging problem. We do not take an extensive comparison of the AR-BCs with other classic BCs, like periodic or reflective, since the topic and related issues have been already widely discussed in several works (see e.g., [19, 33, 46]), where the advantage on some classes of images, in terms of the restored image quality, of the application of AR-BCs has been emphasized. Here we present only a 2D image deblurring example with Gaussian blur and various levels of white Gaussian noise.

The true and the observed images are in Figure 1, where the observed image is affected by a Gaussian blur and 1% noise. We compare the AR-BCs only with the reflective BCs since for this test other BCs like periodic or Dirichlet do not produce satisfactory restorations. In Figure 2 we observe a better restoration and reduced ringing effects at the edges for AR-BCs with respect to reflective BCs. Restored images in Figure 2 are obtained with the minimum relative restoration error varying several values of the regularization parameter 𝜆.

(a)

(b)

Figure 1 

Test problem with Gaussian blur and 1% white Gaussian noise.

(a)

(b)

Figure 2 

Restored images for the test problem in Figure 1.

From Table 1, we note that for the 10% noise case, all of the approaches give comparable restorations. On the other hand, decreasing the noise, that is, passing to 1% and then to 0.1% noise, the AR-BCs improve the restoration while the reflective BCs are not able to do that, due to the barrier of the ringing effects.

Table 1 

RRE for the test problem in Figure 1.

4.2. Reblurring and CG Regularization

In the following tests, the reblurring strategy will be applied with R-BCs and AR-BCs when the PSF is not necessarily strongly symmetric. Indeed, in the case of D-BCs and P-BCs, the reblurring approach is equal to the classical normal equations, while, in the case of R-BCs and AR-BCs, this is no longer true in general.

We provide two problems using iterative regularization by CG.

Test I: Cameraman
The first test is reported in Figure 3. The true image is a cameraman and the 61×61 PSF is associated with a Gaussian distribution in the square domain [−30,30]×[−30,30], with variance equal to four. We add a Gaussian noise.

(a)

(b)

(c)

Figure 3 

Test I: true image, Gaussian PSF, and the blurred image without noise.

Test II: Astronomy
We are dealing with a nonsymmetric experimental 256×256 PSF developed by US Air Force Phillips Laboratory, Lasers and Imaging Directorate, Kirtland Air Force Base, New Mexico, widely used in literature (see e.g. [12, 47]). The true object is the image of Saturn in Figure 4; a Poissonian noise is added, as it is customary when dealing with astronomical images.

(a)

(b)

(c)

Figure 4 

Test II: true image, PSF, and the blurred image without noise.

We show the results corresponding only to P-BCs, R-BCs, and AR-BCs. For shortness, we do not report the reconstructions coming from D-BCs, since the related restorations are usually not better than those with P-BCs.

Tables 2 and 3 show the best RREs for various levels of noise. In Table 2 we also report the 𝐿² norm of the residuals, that is, ‖𝐠−𝐴𝐟‖̂₂, where 𝐠 is the observed image, 𝐟̂ is the computed approximation, and 𝐴 is the coefficient matrix constructed according to the considered BCs: the latter measure is the sum of square errors and it represents, up to the scaling of the variance, the 𝜒² statistical measure of the error. As already pointed out, the choice of the BCs is important mainly for low levels of noise, that is, for high values of SNR. Indeed, in the last row of these tables (SNR =10), the errors due to noise start to dominate the restoration process and therefore the choice of particular BCs is not relevant for the restoration accuracy. In the other rows, where the noise is lower, the choice of the BCs becomes crucial. In particular, the AR-BCs improve substantially the quality of the restorations with respect to the other BCs. This is especially evident in Test II (see Table 3). The reason of the observed high improvement is due to the shape of the PSF, since, basically, the more the support of the PSF is large, the more the ringing effects (and hence the BCs) become dominating.

Table 2 

Test I: best RREs and 𝐿² norm of the residuals for CG with reblurring varying the SNR (SNR =∞ means 0% of noise).

Table 3 

Test II: best RREs for CG with reblurring varying the SNR (SNR =∞ means 0% of noise) within 200 iterations.

To emphasize the quality of the restored images, we consider the reconstruction in the case of Test I for a fixed SNR equal to 40. In Figure 5, we report the restored images and in Figure 6 the residuals of the computed solutions for each pixel divided by the variance of the noise. The last one should have a normal distribution in the case of a good restoration since we add a Gaussian noise. In Figure 5 is evident the reduction of the ringing effects passing from P-BCs to R-BCs (ringings such as the horizontal white line on the top and the horizontal black line on the bottom disappear) and from R-BCs to AR-BCs (ringings such as the two vertical white lines in the bottom left disappear). Indeed, in the same figure we note a higher level of detail in the case of AR-BCs, especially concerning the face of the cameramen. The image restored with R-BCs is smoother when compared with the one restored with AR-BCs, where we can see the “freckles” effect, typical of the 𝐿² norm restoration. Indeed the CGLS method computes the least-squares solution that is well known to be affected by such a phenomenon [39]. When passing to the R-BCs, the considered effect is less evident: indeed it seems that the slightly greater ringing effects smooth the image reducing the “freckles”, but also reducing some details like, for example, the eye of the cameraman. The previous comments suggest that using AR-BCs in connection with regularization methods related to other norms, like the Total Variation [48], could lead to a reconstruction with sharper edges. In Figure 6, we observe a normal distribution of the scaled residuals only with the AR-BCs, while, also with the R-BCs, some further errors corresponding to the ringing effects emerge at the boundary and at the edges of the image: this means that the imposition of the R-BCs is not good enough as model at least for this example. On the other hand, with the AR-BCs, this kind of error seems to disappear and the scaled residual seems to follow a normal distribution. This confirms the goodness of the restoration obtained with the AR-BCs.

(a)

(b)

(c)

Figure 5 

Restored images with SNR=40.

(a)

(b)

(c)

Figure 6 

Residuals divided by the variance of the noise for the restored images with SNR=40.

Two convergence histories, that is, the RREs at any iteration, are plotted in Figure 7 for two different values of the SNR. It should be stressed that the AR-BCs give the best results and the lowest level of RRE. Such behavior is again more evident considering Test II. Indeed, in such case the RREs with AR-BCs continue to decrease even after the first 200 iterations. On the other hand, the restorations with R-BCs start to deteriorate after the very first iterations. In addition, we notice that AR-BCs curves are in general quite flat. This is a very useful feature since the estimation of the optimal stopping value, which is well known to be a crucial and difficult task, can be done with low precision. Indeed, in order to stress the applicability of the AR-BCs to real problems, we consider for the Test I the discrepancy principle widely used with iterative methods [31]. Since we know the 𝐿² norm of the error, we stop the CG when the 𝐿² norm of the residual becomes lower than the 𝐿² norm of the error. Such a criterion seems to work quite well for the AR-BCs as we can see in Figure 8. The restored images are good enough with respect to the optimal solution and also the stopping iteration, at least in this example, is close to the optimal one. On the other hand, such criterion is not always effective for the other BCs in this example. For instance, the stopping iteration for R-BCs is greater than 1000 in the case of SNR = 50 and it is 13 for SNR = 30. However, an analysis of the stopping criterion, in connection with AR-BCs, should be further investigated in the future.

(a)

(b)

Figure 7 

Test I: RREs at each CG iteration with different BCs: dotted line: Periodic, dashed line: Reflective, solid line: Antireflective.

(a)

(b)

(c)

(d)

Figure 8 

Discrepancy principle: images restored at the optimum value iteration and at which computed with the discrepancy principle.

Finally, we remark that also for Test II the CG applied to the linear system, (34) works without breakdown, both with R-BCs and AR-BCs. Therefore, it is possible that the applicability of the CG to (34) is a general property, which does not depend on the particular choice of the BCs.

5. Conclusions

In this contribution we have dealt with the use of antireflective BCs for deblurring problems where the considered issues have been: the precision of the reconstruction when the noise is not present, the linear algebra related to these BCs, the computational costs, and the reconstruction quality associated with iterative and noniterative regularizing solvers when the noise is considered. For many of the considered items, the antireflective algebra coming from the given BCs is the optimal choice: for instance in the work of La Spina it is proven that no one of the trigonometric algebras can be associated with BCs of the same precision as the antireflective. Numerical experiments corroborating the previous statements have been reported and discussed: in this direction it remains an open problem to understand why the CG works without any numerical problem even if the antireflective structure is non symmetric (an event not normal as emphasized by the Jordan decomposition).

Acknowledgment

The paper was partially supported by MIUR 2008, Grant number 20083KLJEZ.