# Antireflective Boundary Conditions for Deblurring Problems

Authors: Marco Donatelli; Stefano Serra-Capizzano
Journal: Journal of Electrical and Computer Engineering
 (2010)

---

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


---

## Body

## 1. Introduction

Formation of a blurred signal/image is typically modeled as a linear system:(1)Af=g,
where f is the true object, g is the blurred object,  and A models the blurring process. Obtaining an accurate deblurring model (i.e., the matrix A) 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 matrixA. 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 A 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 matrixA shows small eigenvalues). This happen because the global matrix is of convolution type. Therefore, the eigenvectors look globally like the Fourier vectors mk (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 k≪n, n being the size of the matrix, such that the low-rank correction belongs to the span of canonical matrices Ei,j=eiejT with i,j≤k and/or i,j≥n+1-k, ek being the kth vector of the canonical basis. In that case, we observe ∥Ei,jmk∥2≪∥mk∥2 with ∥·∥2 being the classical Euclidean norm and hence the term Ei,j is unable to modify significantly the characterization in frequencies of the eigenspaces. In actuality, due to the relation ∥Ei,jmk∥2≪∥mk∥2, the effect induced by any Ei,j 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 f 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 O(n4) operations for an n2-by-n2 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 A 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 A 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 f̃=(…,f-m+1,…,f0,f1,…,fn,fn+1,…,fn+m,…)T and the normalized blurring function given by(2)h=(…,0,0,h-m,h-m+1,…,h0,…,hm-1,hm,0,0,…)T,∑j=-mmhj=1.
The blurred signal g is the convolution of h and f̃ and consequently gi=∑j=-∞∞hjfi-j. The deblurring problem is to recover the vector f=(f1,…,fn)T given the blurring function h and a blurred signal g=(g1,…,gn)T of finite length, that is,(3)g=Ãf̃,Ã=(hm⋯h0⋯h-mhmh0h-m0⋱⋱⋱⋱⋱⋱⋱⋱⋱⋱0hmh0h-mhm⋯h0⋯h-m).
Thus the blurred signal g is determined not by f only, but also by (f-m+1,…,f0)T and (fn+1,…,fn+m)T and the linear system (3) is underdetermined. To overcome this, we make certain assumptions (called BCs) on the unknown boundary data f-m+1,…,f0 and fn+1,…,fn+m 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 meansfj=fn+j=0 for all j in (3) so that a Toeplitz structure is encountered. If we consider periodic BCs we set fj=fn+j, for all j in (3) and the matrix system in (3) becomes n-by-n 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 f are a reflection of the data inside f (refer to [12]) so that f1-j=fj for j=1,…,m and fn+j=fn+1-j for all j=1,…,m in (3). Thus (3) becomes Cf=g, where C is neither Toeplitz nor circulant but it is a special n-by-n Toeplitz plus Hankel matrix which is diagonalized by the discrete cosine transform provided that the blurring function h is symmetric, that is, hj=h-j for all j in (2). It follows that the above system can be solved by using three transforms DCT III in O(nlogn) 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 Section2 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𝒜ℛn≡𝒜ℛ, n≥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 outsidef are an antireflection of the data inside f. More precisely, if x is a point outside the domain and x* is the closest boundary point, then we have x=x*-δx and the quantity f(x) is approximated by f(x*)-(f(x*+δx)-f(x*)) so that we impose(4)f1-j=f1-(fj+1-f1)=2f1-fj+1,forallj=1,…,m,fn+j=fn-(fn-j-fn)=2fn-fn-j,forallj=1,…,m
in (3). If we define the vector z whose components are zj=2∑k=jmhk for j≤m and zero otherwise and if we define the vector w whose components are wn+1-j=2∑k=jmh-k for j≤m and zero otherwise, then (3) becomes(5)Ãf̃=[ze1T-(0∣Tl)J̃+T-(Tr∣0)Ĵ+wenT]f̃=g,
where ek is the kth vector of the canonical basis, J̃=(000J), Ĵ=(J000) with J denoting the n-1 dimensional flip matrix having entries [J]s,t=1 if s+t=n+1 and zero otherwise, and where the matrices Tl,T, and Tr are given by(6)(hm⋯h1⋱⋱hm0),(h0⋯h-m0⋮⋱⋱⋱hm⋱⋱⋱h-m⋱⋱⋱⋮0hm⋯h0),(0h-m⋮⋱h-1⋯h-m).
Therefore, the matrices (0∣Tl)J̃ and (Tr∣0)Ĵ involved in (5) take the form(7)(0h1⋯hm0⋯0⋮⋰hm⋮O⋮0⋯0),(0⋯0⋮O⋮h-m⋰⋮0⋯0h-m⋯h-10).We remark that the coefficient matrixA 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 (n-2) by (n-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 h is symmetric, that is, hj=h-j for all j in (2). It follows that (5) can be solved by using three FSTs in O(nlogn) 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 C1 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

LetQ be the type I sine transform matrix of order n (see [27]) with entries(8)[Q]i,j=2n+1sin(jiπn+1),i,j=1,…,n.
It is known that the real matrix Q is orthogonal and symmetric (Q-1=QT=Q). For any n-dimensional real vector v, the matrix-vector multiplication Qv (DST-I transform) can be computed in O(nlog(n)) real operations by using the algorithm FST-I. Let τ be the space of all the matrices that can be diagonalized by Q:(9)τ={QDQ:Disarealdiagonalmatrixofsizen}.
Let X=QDQ∈τ, then QX=DQ. Consequently, if we let e1 denote the first column of the identity matrix, then the relationship QXe1=DQe1 implies that the eigenvalues [D]i,i of X are given by [D]i,i=[Q(Xe1)]i/[Qe1]i,i=1,…,n. Therefore, the eigenvalues of X can be obtained by applying a DST-I transform to the first column of X 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 h=[h0,…,hn-1]T as σ(h)=[h1,h2,…,hn-1,0]T. According to a Matlab-like notation, we define T(x) to be the n-by-n symmetric Toeplitz matrix whose first column is x and H(x,y) to be the n-by-n Hankel matrix whose first and last column are x and y, respectively. Every matrix of the class (9) can be written as (see  [27])(10)T(h)-H(σ2(h),Jσ2(h)),
where h=[h0,…,hn-1]T∈ℝn and J 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 h(y) evaluated at the grid points vector Gn=[kπ/(n+1)]k=1n, where(11)h(y)=∑|j|≤n-1hjexp(ijy),i2=-1, and hj=h|j| for |j|≤n-1. The τ matrix in (10) is usually denoted by τ(h) and is called the τ matrix generated by the function or symbolh=h(·) (see the seminal paper [27] where this notation was originally proposed).


### 2.3. The AR-Algebras𝒜ℛ

Leth=h(·) be a real-valued cosine polynomial of degree l and let τk(h)≡Qdiag(h(Gk))Q (note that τk(h) coincides with the matrix in (10)-(11), when l≤k-1). Hence the Fourier coefficients of h are such that hi=h-i∈ℝ with hi=0 if |i|>l, and for k=n-2 we can define the one-level ARn(·) operator(12)ARn(h)=[h(0)vn-2(h)τn-2(h)Jvn-2(h)h(0)],
where vn-2(h)=τn-2((ϕ(h))(·))e1 and(13)(ϕ(h))(y)=h(y)-h(0)2(cos(y)-1).
It is interesting to observe that h(y)-h(0) has a zero of order at least 2 at zero (since h is a cosine polynomial) and therefore ϕ(h)=(ϕ(h))(·) is still a cosine polynomial of degree l-1, whose value at zero is -h′′(0)/2 (in other words the function is well defined at zero).As proved in [28], with the above definition of the operator ARn(·), we have(1)αARn(h1)+βARn(h2)=ARn(αh1+βh2),(2)ARn(h1)ARn(h2)=ARn(h1h2),for realα and β and for cosine functions h1=h1(·) and h2=h2(·).These properties allow us to define𝒜ℛ as the algebra (closed under linear combinations, product, and inversion) of matrices ARn(h), with h being a cosine polynomial. By standard interpolation arguments it is easy to see that 𝒜ℛ can be defined as the set of matrices ARn(h), where h is a cosine polynomial of degree at most n-3. Therefore, it is clear that dim(𝒜ℛ)=n-2. Moreover, the algebra 𝒜ℛ is commutative thanks to item 2, since h1(y)h2(y)=h2(y)h1(y) at every y. Consequently, if matrices of the form ARn(h) 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) h is symmetric (i.e., hi=h-i,foralli∈ℤ), if hi=0 for |i|≥n-2 (degree condition), and if h is normalized so that ∑i=-mmhi=1, then the structure of the n×n antireflective blurring matrix A is(14)A=[z00T0z1zm⋮Â⋮zmz100Tz0],
where A1,1=An,n=1, zi=hi+2∑k=i+1mhk, Â has order n-2 and(15)Â=T(h)-H(σ2(h),Jσ2(h)).
According to the brief discussion of Section 2.2, relation (15) implies that Â=τn-2(h) with h(y)=h0+2∑k=1mhkcos(ky) (see (10) and (11)). Moreover in [28] it is proved that A=ARn(h).


### 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 PSFh 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 PSFh, is concisely described in the following result.Theorem 1 (see [28]).
Let the blurring function (PSF)h be symmetric (i.e., hs=h-s), normalized, and satisfying the usual degree condition. As a consequence the eigenvalues of the n×n AR-BCs blurring matrix A in (14), n≥3, are given by h(0)=1 with multiplicity two and h(Gn-2).The proof can be easily derived by (12) which shows that the eigenvalues of ARn(h) are h(0) with multiplicity 2 and those of τn-2(h), that is, h(Gn-2), with multiplicity 1 each.Here we will determine the eigenvectors of every matrixARn(h). In particular, we show that every AR-BCs matrix is diagonalizable, and we demonstrate independence of the eigenvectors from the symbol h. With reference to the notation in (8)–(11), calling qj(n-2) the j th column of Qn-2, and yj(n-2) the jth point of Gn-2, j=1,…,n-2, we have(16)ARn(h)[0qj(n-2)0]=[h(0)vn-2(h)τn-2(h)Jvn-2(h)h(0)][0qj(n-2)0]=[0τn-2(h)qj(n-2)0]=h(yj(n-2))[0qj(n-2)0],
since qj(n-2) is an eigenvector of τn-2(h) and h(yj(n-2)) is the related eigenvalue. Due to the centrosymmetry of the involved matrix, if [1,pT,0]T is an eigenvector of ARn(h) related to the eigenvalue h(0), then the other is its flip, that is, [0,(Jp)T,1]T. Let us look for this eigenvector by imposing the equality(17)ARn(h)[1p0]=h(0)[1p0],
which is equivalent to seeking a vector p that satisfies(18)vn-2(h)+τn-2(h)p=h(0)p.
Since vn-2(h)=τn-2(ϕ(h))e1 by definition of the operator vn-2(·) (see (12) and the lines below), and because of the algebra structure of τn-2 and thanks to the above relation, we deduce that the vector p satisfies the relation(19)τn-2(h-h(0))[-Ln-2-1e1+p]=0,
where Ln-2 is the discrete one-level Laplacian, that is, Ln-2=τn-2(2-2cos(·)). Therefore, by (19), the solution is given by p=Ln-2-1e1. If τn-2(h-h(0)) is invertible (as it happens for every nontrivial PSF, since the unique maximum of the function is reached at y=0, which is not a grid point of Gn-2), then the solution is unique. Hence, independently of h, we have(20)ARn(h)[1pQn-2Jp1]=[1pQn-2Jp1][h(0)diag(h(Gn-2))h(0)].
Now we observe that the jth eigenvector is unitary, j=2,…,n-1, because Qn-2 is unitary: we wish to impose the same condition on the first and the last eigenvector. The interesting fact is that p has an explicit expression. By using standard finite difference techniques, it follows that pj=1-j/(n-1) so that the first eigenvector is exactly the sampling of the function 1-x on the grid j/(n-1) for j=0,…,n-1. Its Euclidean norm is αn=∑j=0n-1j2/(n-1)~n/3, where, for nonnegative sequences βn, γn, the relation γn~βn means γn=βn(1+o(1)). In this way, the (normalized) AR-transform can be defined as
(21)Tn=[αn-1αn-1PQn-2αn-1JPαn-1].Remark 2.
With the normalization condition in (21), all the columns of Tn 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 h(0). However, since q1(n-2) has only positive components and the vector space generated by the first and the last column of Tn contains positive vectors, it follows that Tn 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 Tn since it is related to a fast O(nlog(n)) 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 transformTn. When considering periodic BCs, the transform of the related matrices is the Fourier transform: its jth 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(-ijy). Analogously, when imposing reflective BCs with a symmetric PSF, the transform of the related matrices is the cosine transform: its j 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(jy). 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 C1 continuity at the borders).The previous observation becomes evident in the expression ofTn in (21). Indeed, by defining the one-dimensional grid G̃n=[0,Gn-2T,π]T=[jπ/(n-1)]j=0n-1, which is a subset of [0,π], we infer that the first column of Tn is given by αn-1(1-y/π)|G̃n, the jth column of Tn, j=2,…,n-1, is given by 2/(n-1)(sin(jy))|G̃n, and finally the last column of Tn is given by αn-1(y/π)|G̃n, that is,(22)Tn=[1-yπ,sin(y),…,sin((n-2)y),yπ]|G̃n·Δn,Δn=diag(αn-1,2n-1In-2,αn-1).Finally, it is worth mentioning that the inverse transform is also described in terms of the same block structure since(23)Tn-1=[αn-Qn-2PQn-2-Qn-2JPαn]Theorem 4 (ARn(·) Jordan Canonical Form).
With the notation and assumptions of Theorem1, the n×n AR-BCs blurring matrix A in (14), n≥3, coincides with(24)ARn(h)=Tndiag(h(Ĝn))Tn-1,
where Tn and Tn-1 are defined in (22) and (23), while Ĝn=[0,Gn-2T,0]T.


### 2.5. Multilevel Extension

Here we provide some comments on the extension of our findings tod-dimensional objects with d>1. When d=1, h is a vector, when d=2, h is a 2D array, when d=3, h is a 3D tensor and so forth.With reference to Section2.3 we propose a (canonical) multidimensional extension of the algebras 𝒜ℛ and of the operators ARn(·), n=(n1,…,nd): the idea is to use tensor products. If h=h(·) is d-variate real-valued cosine polynomial, then its Fourier coefficients form a real d-dimensional tensor which is strongly symmetric, that is symmetric with respect to every direction, that is, hj=h|j| for all j∈ℤd. In addition, h(y), y=(y1,…,yd), can be written as a linear combination of independent terms of the form m(y)=∏j=1dcos(αjyj) where any αj is a nonnegative integer. Therefore, we define(25)ARn(m(y))=ARn1(cos(α1y1))⊗⋯⊗ARnd(cos(αdyd)),
where ⊗ denotes Kronecker product, and we force(26)ARn(αh1+βh2)=αARn(h1)+βARn(h2)
for every real α and β and for every d-variate real-valued cosine polynomials h1=h1(·) and h2=h2(·). It is clear that the request that ARn(·) is a linear operator (for d>1, we impose this property in (26) by definition) is sufficient for defining completely the operator in the d-dimensional setting.With the above definition of the operatorARn(·), we have  (1)
αARn(h1)+βARn(h2)=ARn(αh1+βh2),(2)
ARn(h1)ARn(h2)=ARn(h1h2),for realα and β and for cosine functions h1=h1(·) and h2=h2(·).The latter properties of algebra homomorphism allows to define a commutative algebra𝒜ℛ of the matrices ARn(h), with h(·) being a d-variate cosine polynomial. By standard interpolation arguments it is easy to see that 𝒜ℛ can be defined as the set of matrices ARn(h), where h is a d-variate cosine polynomial of degree at most nj-3 in the jth variable for every j ranging in {1,…,d}: we denote the latter polynomial set by 𝒫n-2e(d,even), with e 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 h1,h2∈𝒫n-2e(d,even) the product h1·h2 is the unique polynomial h∈𝒫n-2e(d,even) satisfying the following interpolation condition:(27)h(y)=zy,zy≡h1(y)h2(y),∀y∈Gn-2(d).
If the degree of h1 plus the degree of h2 in the j th variable does not exceed nj-2, j=1,…,d, then the uniqueness of the interpolant implies that h coincides with the product between polynomials in the usual sense. The uniqueness holds also for d≥2 thanks to the tensor form of the grid Gn-2(d) (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, 𝒫n-2e(d,even) has become an algebra, showing the vector-space dimension equal to (n1-2)·(n2-2)⋯(nd-2) which coincides with that of 𝒜ℛn.Without loss of generality and for the sake of notational clarity, in the following we assumenj=n for j=1,…,d. Thanks to the tensor structure emphasized in (25)–(26), and by using Theorem 4 for every term ARn(cos(αjyj)), j=1,…,d, of ARn(m) the d-level extension of such a theorem easily follows. More precisely, if h is a d-variate real-valued cosine symbol related to a d-dimensional strongly symmetric and normalized mask h, then(28)ARn(h)=Tn(d)Dn(Tn(d))-1,Tn(d)=Tn⊗⋯⊗Tn,
(d times) where Dn is the diagonal matrix containing the eigenvalues of ARn(h). The description of Dn in d dimensions is quite involved when compared with the case d=1, implicitly reported in Theorem 1.For a complete analysis of the spectrum ofARn(h) 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 d-variate symbol h(·). Let xn=xn(1)⊗xn(2)⊗⋯⊗xn(d) be a column of Tn(d), with xn(j)∈{αn-1[1,pT,0]T,αn-1[0,(Jp)T,1]T} or xn(j)=[0,qsjT,0]T, 1≤sj≤n-2 and qsj is the (sj) th column of Qn-2, for j=1,…,d. Let (29)ℱxn={j∣xn(j)=αn-1[1,pT,0]Torxn(j)=αn-1[0,(Jp)T,1]T}⊂{1,…,d},
with xn being the generic eigenvector, that is, the generic column of Tn(d). The eigenvalue related to xn is(30)λ=h(y1(n),…,yd(n)),
where yj(n)=0 for j∈ℱxn and yj(n)=πvj/n-1 for j∉ℱxn. Defining the d-dimensional grid(31)Ĝn(d)=Ĝn⊗⋯⊗Ĝn,dtimes,
we can evaluate the d-variate function h on Ĝn(d) by h(reshape(Ĝn(d),n)), where(32)reshape(X,n)
arranges the entries of X in a d-dimensional array of length n 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 (ARn(·) Jordan Canonical Form).
Thend×nd AR-BCs blurring matrix A, obtained when using a strongly symmetric d-dimensional mask h such that hi=0 if |ij|≥n-2 for some j∈{1,…,d} (the d-dimensional degree condition), n≥3, coincides with
(33)ARn(h)=Tn(d)diag(reshape(h(reshape(Ĝn(d),n)),nd))(Tn(d))-1,
where Tn(d) and Ĝn(d) are defined in (28) and (31).It is worth observing that the matrixARn(h) spectrally analyzed in the previous theorem is exactly the same matrix arising from the imposition of AR-BCs applied separately in every direction, when h is the multivariate cosine symbol coming from the d-D tensor mask h defining the shift-invariant d-dimensional blurring operator.


## 2.1. The Antireflective BCs

When defining the antireflective boundary conditions, we assume that the data outsidef are an antireflection of the data inside f. More precisely, if x is a point outside the domain and x* is the closest boundary point, then we have x=x*-δx and the quantity f(x) is approximated by f(x*)-(f(x*+δx)-f(x*)) so that we impose(4)f1-j=f1-(fj+1-f1)=2f1-fj+1,forallj=1,…,m,fn+j=fn-(fn-j-fn)=2fn-fn-j,forallj=1,…,m
in (3). If we define the vector z whose components are zj=2∑k=jmhk for j≤m and zero otherwise and if we define the vector w whose components are wn+1-j=2∑k=jmh-k for j≤m and zero otherwise, then (3) becomes(5)Ãf̃=[ze1T-(0∣Tl)J̃+T-(Tr∣0)Ĵ+wenT]f̃=g,
where ek is the kth vector of the canonical basis, J̃=(000J), Ĵ=(J000) with J denoting the n-1 dimensional flip matrix having entries [J]s,t=1 if s+t=n+1 and zero otherwise, and where the matrices Tl,T, and Tr are given by(6)(hm⋯h1⋱⋱hm0),(h0⋯h-m0⋮⋱⋱⋱hm⋱⋱⋱h-m⋱⋱⋱⋮0hm⋯h0),(0h-m⋮⋱h-1⋯h-m).
Therefore, the matrices (0∣Tl)J̃ and (Tr∣0)Ĵ involved in (5) take the form(7)(0h1⋯hm0⋯0⋮⋰hm⋮O⋮0⋯0),(0⋯0⋮O⋮h-m⋰⋮0⋯0h-m⋯h-10).We remark that the coefficient matrixA 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 (n-2) by (n-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 h is symmetric, that is, hj=h-j for all j in (2). It follows that (5) can be solved by using three FSTs in O(nlogn) 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 C1 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

LetQ be the type I sine transform matrix of order n (see [27]) with entries(8)[Q]i,j=2n+1sin(jiπn+1),i,j=1,…,n.
It is known that the real matrix Q is orthogonal and symmetric (Q-1=QT=Q). For any n-dimensional real vector v, the matrix-vector multiplication Qv (DST-I transform) can be computed in O(nlog(n)) real operations by using the algorithm FST-I. Let τ be the space of all the matrices that can be diagonalized by Q:(9)τ={QDQ:Disarealdiagonalmatrixofsizen}.
Let X=QDQ∈τ, then QX=DQ. Consequently, if we let e1 denote the first column of the identity matrix, then the relationship QXe1=DQe1 implies that the eigenvalues [D]i,i of X are given by [D]i,i=[Q(Xe1)]i/[Qe1]i,i=1,…,n. Therefore, the eigenvalues of X can be obtained by applying a DST-I transform to the first column of X 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 h=[h0,…,hn-1]T as σ(h)=[h1,h2,…,hn-1,0]T. According to a Matlab-like notation, we define T(x) to be the n-by-n symmetric Toeplitz matrix whose first column is x and H(x,y) to be the n-by-n Hankel matrix whose first and last column are x and y, respectively. Every matrix of the class (9) can be written as (see  [27])(10)T(h)-H(σ2(h),Jσ2(h)),
where h=[h0,…,hn-1]T∈ℝn and J 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 h(y) evaluated at the grid points vector Gn=[kπ/(n+1)]k=1n, where(11)h(y)=∑|j|≤n-1hjexp(ijy),i2=-1, and hj=h|j| for |j|≤n-1. The τ matrix in (10) is usually denoted by τ(h) and is called the τ matrix generated by the function or symbolh=h(·) (see the seminal paper [27] where this notation was originally proposed).


## 2.3. The AR-Algebras𝒜ℛ

Leth=h(·) be a real-valued cosine polynomial of degree l and let τk(h)≡Qdiag(h(Gk))Q (note that τk(h) coincides with the matrix in (10)-(11), when l≤k-1). Hence the Fourier coefficients of h are such that hi=h-i∈ℝ with hi=0 if |i|>l, and for k=n-2 we can define the one-level ARn(·) operator(12)ARn(h)=[h(0)vn-2(h)τn-2(h)Jvn-2(h)h(0)],
where vn-2(h)=τn-2((ϕ(h))(·))e1 and(13)(ϕ(h))(y)=h(y)-h(0)2(cos(y)-1).
It is interesting to observe that h(y)-h(0) has a zero of order at least 2 at zero (since h is a cosine polynomial) and therefore ϕ(h)=(ϕ(h))(·) is still a cosine polynomial of degree l-1, whose value at zero is -h′′(0)/2 (in other words the function is well defined at zero).As proved in [28], with the above definition of the operator ARn(·), we have(1)αARn(h1)+βARn(h2)=ARn(αh1+βh2),(2)ARn(h1)ARn(h2)=ARn(h1h2),for realα and β and for cosine functions h1=h1(·) and h2=h2(·).These properties allow us to define𝒜ℛ as the algebra (closed under linear combinations, product, and inversion) of matrices ARn(h), with h being a cosine polynomial. By standard interpolation arguments it is easy to see that 𝒜ℛ can be defined as the set of matrices ARn(h), where h is a cosine polynomial of degree at most n-3. Therefore, it is clear that dim(𝒜ℛ)=n-2. Moreover, the algebra 𝒜ℛ is commutative thanks to item 2, since h1(y)h2(y)=h2(y)h1(y) at every y. Consequently, if matrices of the form ARn(h) 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) h is symmetric (i.e., hi=h-i,foralli∈ℤ), if hi=0 for |i|≥n-2 (degree condition), and if h is normalized so that ∑i=-mmhi=1, then the structure of the n×n antireflective blurring matrix A is(14)A=[z00T0z1zm⋮Â⋮zmz100Tz0],
where A1,1=An,n=1, zi=hi+2∑k=i+1mhk, Â has order n-2 and(15)Â=T(h)-H(σ2(h),Jσ2(h)).
According to the brief discussion of Section 2.2, relation (15) implies that Â=τn-2(h) with h(y)=h0+2∑k=1mhkcos(ky) (see (10) and (11)). Moreover in [28] it is proved that A=ARn(h).


## 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 PSFh 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 PSFh, is concisely described in the following result.Theorem 1 (see [28]).
Let the blurring function (PSF)h be symmetric (i.e., hs=h-s), normalized, and satisfying the usual degree condition. As a consequence the eigenvalues of the n×n AR-BCs blurring matrix A in (14), n≥3, are given by h(0)=1 with multiplicity two and h(Gn-2).The proof can be easily derived by (12) which shows that the eigenvalues of ARn(h) are h(0) with multiplicity 2 and those of τn-2(h), that is, h(Gn-2), with multiplicity 1 each.Here we will determine the eigenvectors of every matrixARn(h). In particular, we show that every AR-BCs matrix is diagonalizable, and we demonstrate independence of the eigenvectors from the symbol h. With reference to the notation in (8)–(11), calling qj(n-2) the j th column of Qn-2, and yj(n-2) the jth point of Gn-2, j=1,…,n-2, we have(16)ARn(h)[0qj(n-2)0]=[h(0)vn-2(h)τn-2(h)Jvn-2(h)h(0)][0qj(n-2)0]=[0τn-2(h)qj(n-2)0]=h(yj(n-2))[0qj(n-2)0],
since qj(n-2) is an eigenvector of τn-2(h) and h(yj(n-2)) is the related eigenvalue. Due to the centrosymmetry of the involved matrix, if [1,pT,0]T is an eigenvector of ARn(h) related to the eigenvalue h(0), then the other is its flip, that is, [0,(Jp)T,1]T. Let us look for this eigenvector by imposing the equality(17)ARn(h)[1p0]=h(0)[1p0],
which is equivalent to seeking a vector p that satisfies(18)vn-2(h)+τn-2(h)p=h(0)p.
Since vn-2(h)=τn-2(ϕ(h))e1 by definition of the operator vn-2(·) (see (12) and the lines below), and because of the algebra structure of τn-2 and thanks to the above relation, we deduce that the vector p satisfies the relation(19)τn-2(h-h(0))[-Ln-2-1e1+p]=0,
where Ln-2 is the discrete one-level Laplacian, that is, Ln-2=τn-2(2-2cos(·)). Therefore, by (19), the solution is given by p=Ln-2-1e1. If τn-2(h-h(0)) is invertible (as it happens for every nontrivial PSF, since the unique maximum of the function is reached at y=0, which is not a grid point of Gn-2), then the solution is unique. Hence, independently of h, we have(20)ARn(h)[1pQn-2Jp1]=[1pQn-2Jp1][h(0)diag(h(Gn-2))h(0)].
Now we observe that the jth eigenvector is unitary, j=2,…,n-1, because Qn-2 is unitary: we wish to impose the same condition on the first and the last eigenvector. The interesting fact is that p has an explicit expression. By using standard finite difference techniques, it follows that pj=1-j/(n-1) so that the first eigenvector is exactly the sampling of the function 1-x on the grid j/(n-1) for j=0,…,n-1. Its Euclidean norm is αn=∑j=0n-1j2/(n-1)~n/3, where, for nonnegative sequences βn, γn, the relation γn~βn means γn=βn(1+o(1)). In this way, the (normalized) AR-transform can be defined as
(21)Tn=[αn-1αn-1PQn-2αn-1JPαn-1].Remark 2.
With the normalization condition in (21), all the columns of Tn 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 h(0). However, since q1(n-2) has only positive components and the vector space generated by the first and the last column of Tn contains positive vectors, it follows that Tn 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 Tn since it is related to a fast O(nlog(n)) 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 transformTn. When considering periodic BCs, the transform of the related matrices is the Fourier transform: its jth 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(-ijy). Analogously, when imposing reflective BCs with a symmetric PSF, the transform of the related matrices is the cosine transform: its j 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(jy). 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 C1 continuity at the borders).The previous observation becomes evident in the expression ofTn in (21). Indeed, by defining the one-dimensional grid G̃n=[0,Gn-2T,π]T=[jπ/(n-1)]j=0n-1, which is a subset of [0,π], we infer that the first column of Tn is given by αn-1(1-y/π)|G̃n, the jth column of Tn, j=2,…,n-1, is given by 2/(n-1)(sin(jy))|G̃n, and finally the last column of Tn is given by αn-1(y/π)|G̃n, that is,(22)Tn=[1-yπ,sin(y),…,sin((n-2)y),yπ]|G̃n·Δn,Δn=diag(αn-1,2n-1In-2,αn-1).Finally, it is worth mentioning that the inverse transform is also described in terms of the same block structure since(23)Tn-1=[αn-Qn-2PQn-2-Qn-2JPαn]Theorem 4 (ARn(·) Jordan Canonical Form).
With the notation and assumptions of Theorem1, the n×n AR-BCs blurring matrix A in (14), n≥3, coincides with(24)ARn(h)=Tndiag(h(Ĝn))Tn-1,
where Tn and Tn-1 are defined in (22) and (23), while Ĝn=[0,Gn-2T,0]T.


## 2.5. Multilevel Extension

Here we provide some comments on the extension of our findings tod-dimensional objects with d>1. When d=1, h is a vector, when d=2, h is a 2D array, when d=3, h is a 3D tensor and so forth.With reference to Section2.3 we propose a (canonical) multidimensional extension of the algebras 𝒜ℛ and of the operators ARn(·), n=(n1,…,nd): the idea is to use tensor products. If h=h(·) is d-variate real-valued cosine polynomial, then its Fourier coefficients form a real d-dimensional tensor which is strongly symmetric, that is symmetric with respect to every direction, that is, hj=h|j| for all j∈ℤd. In addition, h(y), y=(y1,…,yd), can be written as a linear combination of independent terms of the form m(y)=∏j=1dcos(αjyj) where any αj is a nonnegative integer. Therefore, we define(25)ARn(m(y))=ARn1(cos(α1y1))⊗⋯⊗ARnd(cos(αdyd)),
where ⊗ denotes Kronecker product, and we force(26)ARn(αh1+βh2)=αARn(h1)+βARn(h2)
for every real α and β and for every d-variate real-valued cosine polynomials h1=h1(·) and h2=h2(·). It is clear that the request that ARn(·) is a linear operator (for d>1, we impose this property in (26) by definition) is sufficient for defining completely the operator in the d-dimensional setting.With the above definition of the operatorARn(·), we have  (1)
αARn(h1)+βARn(h2)=ARn(αh1+βh2),(2)
ARn(h1)ARn(h2)=ARn(h1h2),for realα and β and for cosine functions h1=h1(·) and h2=h2(·).The latter properties of algebra homomorphism allows to define a commutative algebra𝒜ℛ of the matrices ARn(h), with h(·) being a d-variate cosine polynomial. By standard interpolation arguments it is easy to see that 𝒜ℛ can be defined as the set of matrices ARn(h), where h is a d-variate cosine polynomial of degree at most nj-3 in the jth variable for every j ranging in {1,…,d}: we denote the latter polynomial set by 𝒫n-2e(d,even), with e 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 h1,h2∈𝒫n-2e(d,even) the product h1·h2 is the unique polynomial h∈𝒫n-2e(d,even) satisfying the following interpolation condition:(27)h(y)=zy,zy≡h1(y)h2(y),∀y∈Gn-2(d).
If the degree of h1 plus the degree of h2 in the j th variable does not exceed nj-2, j=1,…,d, then the uniqueness of the interpolant implies that h coincides with the product between polynomials in the usual sense. The uniqueness holds also for d≥2 thanks to the tensor form of the grid Gn-2(d) (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, 𝒫n-2e(d,even) has become an algebra, showing the vector-space dimension equal to (n1-2)·(n2-2)⋯(nd-2) which coincides with that of 𝒜ℛn.Without loss of generality and for the sake of notational clarity, in the following we assumenj=n for j=1,…,d. Thanks to the tensor structure emphasized in (25)–(26), and by using Theorem 4 for every term ARn(cos(αjyj)), j=1,…,d, of ARn(m) the d-level extension of such a theorem easily follows. More precisely, if h is a d-variate real-valued cosine symbol related to a d-dimensional strongly symmetric and normalized mask h, then(28)ARn(h)=Tn(d)Dn(Tn(d))-1,Tn(d)=Tn⊗⋯⊗Tn,
(d times) where Dn is the diagonal matrix containing the eigenvalues of ARn(h). The description of Dn in d dimensions is quite involved when compared with the case d=1, implicitly reported in Theorem 1.For a complete analysis of the spectrum ofARn(h) 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 d-variate symbol h(·). Let xn=xn(1)⊗xn(2)⊗⋯⊗xn(d) be a column of Tn(d), with xn(j)∈{αn-1[1,pT,0]T,αn-1[0,(Jp)T,1]T} or xn(j)=[0,qsjT,0]T, 1≤sj≤n-2 and qsj is the (sj) th column of Qn-2, for j=1,…,d. Let (29)ℱxn={j∣xn(j)=αn-1[1,pT,0]Torxn(j)=αn-1[0,(Jp)T,1]T}⊂{1,…,d},
with xn being the generic eigenvector, that is, the generic column of Tn(d). The eigenvalue related to xn is(30)λ=h(y1(n),…,yd(n)),
where yj(n)=0 for j∈ℱxn and yj(n)=πvj/n-1 for j∉ℱxn. Defining the d-dimensional grid(31)Ĝn(d)=Ĝn⊗⋯⊗Ĝn,dtimes,
we can evaluate the d-variate function h on Ĝn(d) by h(reshape(Ĝn(d),n)), where(32)reshape(X,n)
arranges the entries of X in a d-dimensional array of length n 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 (ARn(·) Jordan Canonical Form).
Thend×nd AR-BCs blurring matrix A, obtained when using a strongly symmetric d-dimensional mask h such that hi=0 if |ij|≥n-2 for some j∈{1,…,d} (the d-dimensional degree condition), n≥3, coincides with
(33)ARn(h)=Tn(d)diag(reshape(h(reshape(Ĝn(d),n)),nd))(Tn(d))-1,
where Tn(d) and Ĝn(d) are defined in (28) and (31).It is worth observing that the matrixARn(h) spectrally analyzed in the previous theorem is exactly the same matrix arising from the imposition of AR-BCs applied separately in every direction, when h is the multivariate cosine symbol coming from the d-D tensor mask h defining the shift-invariant d-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 ATA and that the righthand-side is given by ATg. Quite surprisingly, the AR-BCs may be spoiled by this approach at least for d=1 and if we compute explicitly the matrix ATA and the vector ATg, 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 AT, 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 kernelh is symmetric, the matrix AT is still a blurring operator since AT=A, while, in the case of the AR-BCs matrix, AT 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 AT 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 ATg 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 matrixATA is also spoiled and, in the case of images (d=2) we lose the O(n2log(n)) computational cost for solving a generic system ATAx=b. The cost of solving such a type of linear systems is proportional to n3 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 d-dimensional setting the solution of the normal equation linear system is asymptotic to n3(d-1), where nd is the size of the matrix A. In order to overcome these problems (which arise only with the most precise AR-BCs for strongly symmetric PSFs), we replace AT by A′, where A′ 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 asA2f=Ag, whenever zero Dirichlet, periodic or reflective BCs are considered. Therefore, the observed image g 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 A2f=Ag and Af=g are algebraically equivalent if A is invertible: in any case, if A is also positive definite, the first represents the minimization of the functional ∥Af-g∥22 while the second represents the minimization of the functional ∥A1/2f-A-1/2g∥22 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 AT≠A, we can replace AT by A 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 K be the convolution operator with shift-invariant kernel k(s), then [Kf](x)=∫k(x-y)f(y)dy. Since the PSF (and then k) is real (and then real valued), the adjoint operator K* is [K*f](x)=∫k(y-x)f(y)dy 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., k(s) vs k(-s)), and, in the 2D case, the change of sign in the variable s 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 replaceAT by A′ (the reblurring matrix), where A′ is the matrix obtained by imposing the BCs to the PSF rotated by 180 degrees. Using Matlab notation, if H is a q×q PSF, its 180 degrees rotated version is H′=fliplr(flipud(H))=JqHJq, where Jq is the flip matrix defined as [Jq]i,j=1 if i+j=q+1 for i,j=1,…,q, and zero otherwise. For a d-dimensional problem, A′ 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 wayA′ has the same computational properties of A (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)A′Af=A′g.
Furthermore, using the Tikhonov regularization in the reblurred version, we propose to use(35)[A′A+μL′L]f=A′g,
with L being any discrete.  differential operator with AR-BCs. In general A′A 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 |h|2 if h 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 AT is replaced by A′ 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 operatorsK and K*, 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 equationsAf=g. 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 K and K* in the infinite dimensional normal equations (where the adjoint operator K* performs a correlation operation) and then we obtain (34). More precisely, the discretization of K and K* in the continuous equation K*Kf=K*g with any fixed BCs gives (34).


### 3.2. Linear Algebra and Computational Issues

We note that in the 1D caseAn′=JnAnJn. In the d-dimensional case, let n=(n1,n2,…,nd) be the partial dimensions of the problem, whose total size is ∏i=1dni. By flipping each variable, we obtain(36)An′=JnAnJn,Jn=⨂i=1dJni.
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 K. Let us consider the Toeplitz d-level matrix Tn(ϕ) of partial dimensions n=(n1,…,nd)∈ℕ+d and generating function ϕ [2], which is defined as (37)Tn(ϕ)=∑|j|≤n-eajZn[j]=∑|j1|<n1⋯∑|jd|<nda(j1,…,jd)Zn1[j1]⊗⋯⊗Znd[jd]
(e=(1,…,1)∈ℕ+d) by means of the Fourier coefficients of ϕ(38)ak=1(2π)d∫[-π,π]dϕ(x)exp-i〈k∣x〉dx,i2=-1,k∈ℤd.
Here 〈k∣x〉=∑i=1dkixi and for j∈ℤ, m∈N+, Zm[j]∈ℝm×m is the matrix whose (s,r) th entry is 1 if s-r=j, and 0 elsewhere. As it is well known for multilevel Toeplitz matrices TnH(ϕ)=Tn(ϕ¯), where ϕ¯ is the conjugate of the function ϕ, and the Fourier coefficients of ϕ¯ are the same of ϕ, but conjugated and flipped. Moreover, since Tn(ϕ¯)=JnTn(ϕ)¯Jn, if An=Tn(ϕ) is real then AnT=JnAnJn=An′. 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 An 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 An′=AnT is JnAn=(JnAn)T (or equivalently AnJn=(AnJn)T), which is a multilevel antidiagonal symmetry called persymmetry. Therefore in the case of R-BCs, where the matrix An involves nested Hankel parts, in general An′≠AnT, while An′=AnT only when the PSF is strongly symmetric since in this case JnAn=(JnAn)T. Dealing with the AR-BCs, the situation is even more involved, since A′n≠AnT 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 JnAn≠(JnAn)T. As a nontrivial and unexpected example, it is important to stress that the imposition of R-BCs with nonstrong symmetric PSFs implies JnAn≠(JnAn)T, that is, An′≠AnT.We provide now a computational motivation for the choice of usingA′ as an alternative to AT:A′ 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 X and H be two matrices such that X represents an n×n image and H is the discrete q×q 2D PSF, q=2m+1, which leads to the matrix blurring A. By using the Matlab notation x=X(:) (i.e., the vector x is the column-stacked version of X), the product Ax in the 2D case can be implemented  (1)
by using an enlarged imageX̃, which is the (n+q-1)×(n+q-1) image X extended at the boundaries according to the imposed BCs,(2)
computingH*X̃, where the symbol “*” denotes the circular convolution operator (H should be zero padded to have the same size of X̃),(3)
and then taking the inner part of the result having the same size ofX.The circular convolution can be computed using the 2D discrete Fourier transform (DFT2) and its inverse (IDFT2), since we have(39)H*X̃=IDFT2(DFT2(H)⊙DFT2(X̃)),
where “⊙” is the componentwise product. If 𝒞k(f) denotes the block circulant matrix with circulant blocks, of block size k with blocks of size k and generating function f, then (39) represents the same operation as 𝒞n+q-1(ϕ)x̃, where x̃=X̃(:) (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(40)Y=H′*X̃=IDFT2(DFT2(H)¯⊙DFT2(X̃)),H'i,j=H-i,-j,i,j∈ℤ,
where the overline symbol denotes the complex conjugation. As a result, since the transform H↦H′ is equivalent to the transform ϕ↦ϕ¯ (because ϕ(x,y)=∑i,j∈ℤHi,jexpi〈(i,j)∣(x,y)〉), and since 𝒞n+q-1(ϕ¯)=𝒞n+q-1T(ϕ), if y=Y(:), then it follows that y=𝒞n+q-1T(ϕ)x̃. Therefore, for any of the considered BCs, the inner part of y is exactly A′x. Here it is worthwhile to specify exactly what we mean for inner part: if the vector y is partitioned in n+q-1 blocks of size n+q-1, q=2m+1, then for inner part we mean that we are excluding the first and the last m blocks and, in any of the remaining blocks, we are deleting the first and the last m entries. More generally, if the PSF is arbitrary (e.g. nonsymmetric) that is, the nonzero coefficients of the PSF have first index belonging to [-m1-,m1+] and second index in the range [-m2-,m2+], then we have to delete the first m1- and the last m1+ blocks and, in any of the other blocks, we have to exclude the first m2- and the last m2+ entries.Since the DFT and its inverse can be computed inO(n2log(n)) 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 Section1, regardless of the imposed BCs, matrices A 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 A is (41)A=Tndiag(d)Tn-1,Tn=[t1t2⋯tn],Tn-1=[t̃1Tt̃2T⋮t̃nT],
with d=h(Ĝn), then a spectral filter solution is given by(42)freg=∑i=1nϕit̃iTgditi,
where ϕi are filter factors that satisfy (43)ϕi≈{1,ifdiislarge,0,ifdiissmall.
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(44)ϕitsvd={1,ifdi≥δ,0,ifdi<δ,ϕitik=di2di2+λ,λ>0,
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., di 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 di 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 Tn is not orthogonal (indeed it is a rank-2 correction of a symmetric orthogonal matrix). Therefore, note that the use of ϕitsvd 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)g̃=Tn-1g.
(2)d=[h(0),d̂T,h(0)]T, where d̂=[d2,…,dn-1]T are the eigenvalues of τn-2(h) that can be computed by a fast sine transform (FST).
(3)f̃=(ϕ./d)⨀g̃, where the dot operations are component-wise.
(4)freg=Tnf̃.The productTnf̃ can be clearly computed in a fast and stable way by one FST. Indeed for all x∈ℝn(45)Tnx=αn-1x1[1p0]+[0Qn-2x(2:n-1)0]+αn-1xn[0Jp1],
where x(2:n-1) in Matlab notation is the vector x with components indexed from 2 to n-1. A similar strategy can be followed for computing the matrix-vector product Tn-1g. Instead of αn-1p there is u=-Qn-2p and instead of αn-1Jp there is w=-Qn-2Jp. Recalling that p=Ln-2-1e1 the two vectors u and w can be explicitly computed obtaining ui=(2n-2)-1/2cot(iπ/(2n-2)), for i=1,…,n-2 and w=diagi=1,…,n-2(-1)i+1u.Remark 7.
As discussed in Remark3, 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 h at a grid-point close to zero, that is, close to the maximum point of h, has an associated eigenvector that corresponds to low-frequency (signal space) information, while the eigenvalue obtained as a sampling of the symbol h 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/(n-1)(sin(jy))|G̃n, j=2,…,n-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 h(0), that is, the maximum and the infinity norm of h, and by the fact that AR-BCs are more precise than other classical BCs.Remark 8.
Ford=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 d-level FSTs of size (n-2)d, 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 d-level FST of size (n-2)d 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(46)Vn=[1|0Q0|l],
where (47)1=1n[11⋮1],l=3n(n2-1)[1-n3-n⋮n-1].Note that1 and l differ from the corresponding eigenvectors [1,pT,0]T and [0,JpT,1]T 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 Theorem4 the spectral decomposition of ARn(h) can now be written as(48)ARn(h)=VnΛVn-1,
where the diagonal entries λjj of Λ are given by(49)λjj={h(j-1n-1π),1≤j<n,h(0),j=n.
Here we prove that the reblurring approach for Tikhonov regularization in (35) where L=I is a regularization method. For a complete analysis we need to compute the SVD of Vn.We use the notationy≐z if the two vectors y and z depend on n and for each entry yi of y and the corresponding entry zi of z there holds yi/zi→1 as n→∞.Theorem 9 (see [43]).
The two dominant singular values ofVn are given by
(50)σ1≐σ2≐2,
where 2 is, in fact, a strict upper bound, and the two minimal singular values are given by
(51)σn-1≐3nandσn≐1n,
respectively. Fix 1′=1(2:n-1), l′=l(2:n-1), a=Q1′, and b=Ql′. The corresponding right singular vectors are
(52)v1≐12[1a0],v2≐12[0b1],vn-1≐12[0b-1],vn≐12[-1a0],
and the left singular vectors are
(53)u1≐[1(2n)1′1(2n)],u2≐[l12l′ln2],un-1≐12[13nl′-1],un≐12[-11′n-1],
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(54)μ(Vn)≐2n,n→∞.Remark 11.
It is important to note that the ill-conditioned subspace ofVn-1 has dimension two, independent of n, since Vn has two singular values that decay like 1/n while all others are between one and two. Also, Vn-1 only amplifies vectors that fail to be orthogonal to 𝒰=span{un-1,un}. 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.
Letf be itself a blurred version of a continuous signal w, that is,
(55)f(x)=∫ℝk(x-y)w(y)dy,x∈ℝ,
where w satisfies the same BCs as f (i.e., periodic, reflective, or antireflective ones).On the grounds of Assumption12 we may therefore assume that(56)f=Aw,forsomew∈ℝn
with a moderate bound(57)|∥w∥|≤∥w∥∞≤ϱ,
where |∥x∥|:=∥x∥2/n,x∈ℝn. Since the observed object is usually affected by noise, instead of the blurred object g we have to work with the blurred and noisy object g̃=g+e, where ∥e∥∞<ε. In this way, for L=I (35) becomes (A'A+μI)f=A'g.Theorem 13.
Let the exact solutionf of (1) satisfy (56) with (57). Furthermore the total error of the reblurring strategy with AR BCs satisfies
(58)|∥fαε-f∥|=O(εϱ),
for α=α(ε)=ε/ϱ, where the constant in the O(·)-notation is independent of the dimension n.Note that the upper bound from Theorem13 is the same as for Tikhonov regularization with reflective or periodic BCs; only the constant hidden in the O(·)-notation may be somewhat larger for the reblurring strategy.


## 3.1. The Reblurring Operator

A key point is that, for zero Dirichlet, periodic and reflective BCs, when the kernelh is symmetric, the matrix AT is still a blurring operator since AT=A, while, in the case of the AR-BCs matrix, AT 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 AT 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 ATg 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 matrixATA is also spoiled and, in the case of images (d=2) we lose the O(n2log(n)) computational cost for solving a generic system ATAx=b. The cost of solving such a type of linear systems is proportional to n3 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 d-dimensional setting the solution of the normal equation linear system is asymptotic to n3(d-1), where nd is the size of the matrix A. In order to overcome these problems (which arise only with the most precise AR-BCs for strongly symmetric PSFs), we replace AT by A′, where A′ 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 asA2f=Ag, whenever zero Dirichlet, periodic or reflective BCs are considered. Therefore, the observed image g 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 A2f=Ag and Af=g are algebraically equivalent if A is invertible: in any case, if A is also positive definite, the first represents the minimization of the functional ∥Af-g∥22 while the second represents the minimization of the functional ∥A1/2f-A-1/2g∥22 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 AT≠A, we can replace AT by A 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 K be the convolution operator with shift-invariant kernel k(s), then [Kf](x)=∫k(x-y)f(y)dy. Since the PSF (and then k) is real (and then real valued), the adjoint operator K* is [K*f](x)=∫k(y-x)f(y)dy 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., k(s) vs k(-s)), and, in the 2D case, the change of sign in the variable s 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 replaceAT by A′ (the reblurring matrix), where A′ is the matrix obtained by imposing the BCs to the PSF rotated by 180 degrees. Using Matlab notation, if H is a q×q PSF, its 180 degrees rotated version is H′=fliplr(flipud(H))=JqHJq, where Jq is the flip matrix defined as [Jq]i,j=1 if i+j=q+1 for i,j=1,…,q, and zero otherwise. For a d-dimensional problem, A′ 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 wayA′ has the same computational properties of A (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)A′Af=A′g.
Furthermore, using the Tikhonov regularization in the reblurred version, we propose to use(35)[A′A+μL′L]f=A′g,
with L being any discrete.  differential operator with AR-BCs. In general A′A 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 |h|2 if h 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 AT is replaced by A′ 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 operatorsK and K*, 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 equationsAf=g. 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 K and K* in the infinite dimensional normal equations (where the adjoint operator K* performs a correlation operation) and then we obtain (34). More precisely, the discretization of K and K* in the continuous equation K*Kf=K*g with any fixed BCs gives (34).


## 3.2. Linear Algebra and Computational Issues

We note that in the 1D caseAn′=JnAnJn. In the d-dimensional case, let n=(n1,n2,…,nd) be the partial dimensions of the problem, whose total size is ∏i=1dni. By flipping each variable, we obtain(36)An′=JnAnJn,Jn=⨂i=1dJni.
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 K. Let us consider the Toeplitz d-level matrix Tn(ϕ) of partial dimensions n=(n1,…,nd)∈ℕ+d and generating function ϕ [2], which is defined as (37)Tn(ϕ)=∑|j|≤n-eajZn[j]=∑|j1|<n1⋯∑|jd|<nda(j1,…,jd)Zn1[j1]⊗⋯⊗Znd[jd]
(e=(1,…,1)∈ℕ+d) by means of the Fourier coefficients of ϕ(38)ak=1(2π)d∫[-π,π]dϕ(x)exp-i〈k∣x〉dx,i2=-1,k∈ℤd.
Here 〈k∣x〉=∑i=1dkixi and for j∈ℤ, m∈N+, Zm[j]∈ℝm×m is the matrix whose (s,r) th entry is 1 if s-r=j, and 0 elsewhere. As it is well known for multilevel Toeplitz matrices TnH(ϕ)=Tn(ϕ¯), where ϕ¯ is the conjugate of the function ϕ, and the Fourier coefficients of ϕ¯ are the same of ϕ, but conjugated and flipped. Moreover, since Tn(ϕ¯)=JnTn(ϕ)¯Jn, if An=Tn(ϕ) is real then AnT=JnAnJn=An′. 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 An 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 An′=AnT is JnAn=(JnAn)T (or equivalently AnJn=(AnJn)T), which is a multilevel antidiagonal symmetry called persymmetry. Therefore in the case of R-BCs, where the matrix An involves nested Hankel parts, in general An′≠AnT, while An′=AnT only when the PSF is strongly symmetric since in this case JnAn=(JnAn)T. Dealing with the AR-BCs, the situation is even more involved, since A′n≠AnT 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 JnAn≠(JnAn)T. As a nontrivial and unexpected example, it is important to stress that the imposition of R-BCs with nonstrong symmetric PSFs implies JnAn≠(JnAn)T, that is, An′≠AnT.We provide now a computational motivation for the choice of usingA′ as an alternative to AT:A′ 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 X and H be two matrices such that X represents an n×n image and H is the discrete q×q 2D PSF, q=2m+1, which leads to the matrix blurring A. By using the Matlab notation x=X(:) (i.e., the vector x is the column-stacked version of X), the product Ax in the 2D case can be implemented  (1)
by using an enlarged imageX̃, which is the (n+q-1)×(n+q-1) image X extended at the boundaries according to the imposed BCs,(2)
computingH*X̃, where the symbol “*” denotes the circular convolution operator (H should be zero padded to have the same size of X̃),(3)
and then taking the inner part of the result having the same size ofX.The circular convolution can be computed using the 2D discrete Fourier transform (DFT2) and its inverse (IDFT2), since we have(39)H*X̃=IDFT2(DFT2(H)⊙DFT2(X̃)),
where “⊙” is the componentwise product. If 𝒞k(f) denotes the block circulant matrix with circulant blocks, of block size k with blocks of size k and generating function f, then (39) represents the same operation as 𝒞n+q-1(ϕ)x̃, where x̃=X̃(:) (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(40)Y=H′*X̃=IDFT2(DFT2(H)¯⊙DFT2(X̃)),H'i,j=H-i,-j,i,j∈ℤ,
where the overline symbol denotes the complex conjugation. As a result, since the transform H↦H′ is equivalent to the transform ϕ↦ϕ¯ (because ϕ(x,y)=∑i,j∈ℤHi,jexpi〈(i,j)∣(x,y)〉), and since 𝒞n+q-1(ϕ¯)=𝒞n+q-1T(ϕ), if y=Y(:), then it follows that y=𝒞n+q-1T(ϕ)x̃. Therefore, for any of the considered BCs, the inner part of y is exactly A′x. Here it is worthwhile to specify exactly what we mean for inner part: if the vector y is partitioned in n+q-1 blocks of size n+q-1, q=2m+1, then for inner part we mean that we are excluding the first and the last m blocks and, in any of the remaining blocks, we are deleting the first and the last m entries. More generally, if the PSF is arbitrary (e.g. nonsymmetric) that is, the nonzero coefficients of the PSF have first index belonging to [-m1-,m1+] and second index in the range [-m2-,m2+], then we have to delete the first m1- and the last m1+ blocks and, in any of the other blocks, we have to exclude the first m2- and the last m2+ entries.Since the DFT and its inverse can be computed inO(n2log(n)) 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 Section1, regardless of the imposed BCs, matrices A 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 A is (41)A=Tndiag(d)Tn-1,Tn=[t1t2⋯tn],Tn-1=[t̃1Tt̃2T⋮t̃nT],
with d=h(Ĝn), then a spectral filter solution is given by(42)freg=∑i=1nϕit̃iTgditi,
where ϕi are filter factors that satisfy (43)ϕi≈{1,ifdiislarge,0,ifdiissmall.
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(44)ϕitsvd={1,ifdi≥δ,0,ifdi<δ,ϕitik=di2di2+λ,λ>0,
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., di 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 di 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 Tn is not orthogonal (indeed it is a rank-2 correction of a symmetric orthogonal matrix). Therefore, note that the use of ϕitsvd 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)g̃=Tn-1g.
(2)d=[h(0),d̂T,h(0)]T, where d̂=[d2,…,dn-1]T are the eigenvalues of τn-2(h) that can be computed by a fast sine transform (FST).
(3)f̃=(ϕ./d)⨀g̃, where the dot operations are component-wise.
(4)freg=Tnf̃.The productTnf̃ can be clearly computed in a fast and stable way by one FST. Indeed for all x∈ℝn(45)Tnx=αn-1x1[1p0]+[0Qn-2x(2:n-1)0]+αn-1xn[0Jp1],
where x(2:n-1) in Matlab notation is the vector x with components indexed from 2 to n-1. A similar strategy can be followed for computing the matrix-vector product Tn-1g. Instead of αn-1p there is u=-Qn-2p and instead of αn-1Jp there is w=-Qn-2Jp. Recalling that p=Ln-2-1e1 the two vectors u and w can be explicitly computed obtaining ui=(2n-2)-1/2cot(iπ/(2n-2)), for i=1,…,n-2 and w=diagi=1,…,n-2(-1)i+1u.Remark 7.
As discussed in Remark3, 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 h at a grid-point close to zero, that is, close to the maximum point of h, has an associated eigenvector that corresponds to low-frequency (signal space) information, while the eigenvalue obtained as a sampling of the symbol h 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/(n-1)(sin(jy))|G̃n, j=2,…,n-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 h(0), that is, the maximum and the infinity norm of h, and by the fact that AR-BCs are more precise than other classical BCs.Remark 8.
Ford=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 d-level FSTs of size (n-2)d, 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 d-level FST of size (n-2)d 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(46)Vn=[1|0Q0|l],
where (47)1=1n[11⋮1],l=3n(n2-1)[1-n3-n⋮n-1].Note that1 and l differ from the corresponding eigenvectors [1,pT,0]T and [0,JpT,1]T 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 Theorem4 the spectral decomposition of ARn(h) can now be written as(48)ARn(h)=VnΛVn-1,
where the diagonal entries λjj of Λ are given by(49)λjj={h(j-1n-1π),1≤j<n,h(0),j=n.
Here we prove that the reblurring approach for Tikhonov regularization in (35) where L=I is a regularization method. For a complete analysis we need to compute the SVD of Vn.We use the notationy≐z if the two vectors y and z depend on n and for each entry yi of y and the corresponding entry zi of z there holds yi/zi→1 as n→∞.Theorem 9 (see [43]).
The two dominant singular values ofVn are given by
(50)σ1≐σ2≐2,
where 2 is, in fact, a strict upper bound, and the two minimal singular values are given by
(51)σn-1≐3nandσn≐1n,
respectively. Fix 1′=1(2:n-1), l′=l(2:n-1), a=Q1′, and b=Ql′. The corresponding right singular vectors are
(52)v1≐12[1a0],v2≐12[0b1],vn-1≐12[0b-1],vn≐12[-1a0],
and the left singular vectors are
(53)u1≐[1(2n)1′1(2n)],u2≐[l12l′ln2],un-1≐12[13nl′-1],un≐12[-11′n-1],
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(54)μ(Vn)≐2n,n→∞.Remark 11.
It is important to note that the ill-conditioned subspace ofVn-1 has dimension two, independent of n, since Vn has two singular values that decay like 1/n while all others are between one and two. Also, Vn-1 only amplifies vectors that fail to be orthogonal to 𝒰=span{un-1,un}. 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.
Letf be itself a blurred version of a continuous signal w, that is,
(55)f(x)=∫ℝk(x-y)w(y)dy,x∈ℝ,
where w satisfies the same BCs as f (i.e., periodic, reflective, or antireflective ones).On the grounds of Assumption12 we may therefore assume that(56)f=Aw,forsomew∈ℝn
with a moderate bound(57)|∥w∥|≤∥w∥∞≤ϱ,
where |∥x∥|:=∥x∥2/n,x∈ℝn. Since the observed object is usually affected by noise, instead of the blurred object g we have to work with the blurred and noisy object g̃=g+e, where ∥e∥∞<ε. In this way, for L=I (35) becomes (A'A+μI)f=A'g.Theorem 13.
Let the exact solutionf of (1) satisfy (56) with (57). Furthermore the total error of the reblurring strategy with AR BCs satisfies
(58)|∥fαε-f∥|=O(εϱ),
for α=α(ε)=ε/ϱ, where the constant in the O(·)-notation is independent of the dimension n.Note that the upper bound from Theorem13 is the same as for Tikhonov regularization with reflective or periodic BCs; only the constant hidden in the O(·)-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 ∥f̂-f∥2/∥f∥2, where f̂ is the computed approximation of the true image f. The signal-to-noise ratio (SNR) is computed as 20log10∥gb∥2/∥ν∥2, where gb 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 Section3.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 Figure1, 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 λ.Test problem with Gaussian blur and1% white Gaussian noise.
(a)(b)Restored images for the test problem in Figure1.
(a)(b)From Table1, 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 Figure1.
NoiseReflectiveAR10%0.12840.12611%0.11880.10340.1%0.11860.0989


### 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 Figure3. 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.Test I: true image, Gaussian PSF, and the blurred image without noise.
(a)(b)(c)Test II: Astronomy
We are dealing with a nonsymmetric experimental256×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.Test II: true image, PSF, and the blurred image without noise.
(a)(b)(c)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.Tables2 and 3 show the best RREs for various levels of noise. In Table 2 we also report the L2 norm of the residuals, that is, ∥g-Af̂∥2, where g is the observed image, f̂ is the computed approximation, and A 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 χ2 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 andL2 norm of the residuals for CG with reblurring varying the SNR (SNR =∞ means 0% of noise).
SNRRelative restoration errorsL2 norm of the residuals (∥g-Af̂∥2)PeriodicReflectiveAntireflectivePeriodicReflectiveAntireflective∞0.22750.19930.18311.23630.01130.0004500.22760.19960.18501.27280.04800.0374400.22780.20070.19211.60780.37660.3654300.23000.20880.20514.93953.73463.7133200.24870.23820.237838.235737.210037.2149100.38360.38140.3823379.2445376.0303376.0419Table 3
Test II: best RREs for CG with reblurring varying the SNR (SNR=∞ means 0% of noise) within 200 iterations.
SNRPeriodicReflectiveAntireflective∞0.24150.15290.0604500.24160.15300.0616400.24180.15320.0737300.24300.15750.1052200.25740.18380.1605100.36890.33290.3289To 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 Figure5, 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 L2 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.Restored images with SNR=40.
(a)(b)(c)Residuals divided by the variance of the noise for the restored images with SNR=40.
(a)(b)(c)Two convergence histories, that is, the RREs at any iteration, are plotted in Figure7 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 L2 norm of the error, we stop the CG when the L2 norm of the residual becomes lower than the L2 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.Test I: RREs at each CG iteration with different BCs: dotted line: Periodic, dashed line: Reflective, solid line: Antireflective.
(a)(b)Discrepancy principle: images restored at the optimum value iteration and at which computed with the discrepancy principle.
(a)(b)(c)(d)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.


## 4.1. Reblurring and Tikhonov Regularization

Let us begin with an example illustrating the approach discussed in Section3.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 Figure1, 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 λ.Test problem with Gaussian blur and1% white Gaussian noise.
(a)(b)Restored images for the test problem in Figure1.
(a)(b)From Table1, 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 Figure1.
NoiseReflectiveAR10%0.12840.12611%0.11880.10340.1%0.11860.0989


## 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 Figure3. 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.Test I: true image, Gaussian PSF, and the blurred image without noise.
(a)(b)(c)Test II: Astronomy
We are dealing with a nonsymmetric experimental256×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.Test II: true image, PSF, and the blurred image without noise.
(a)(b)(c)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.Tables2 and 3 show the best RREs for various levels of noise. In Table 2 we also report the L2 norm of the residuals, that is, ∥g-Af̂∥2, where g is the observed image, f̂ is the computed approximation, and A 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 χ2 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 andL2 norm of the residuals for CG with reblurring varying the SNR (SNR =∞ means 0% of noise).
SNRRelative restoration errorsL2 norm of the residuals (∥g-Af̂∥2)PeriodicReflectiveAntireflectivePeriodicReflectiveAntireflective∞0.22750.19930.18311.23630.01130.0004500.22760.19960.18501.27280.04800.0374400.22780.20070.19211.60780.37660.3654300.23000.20880.20514.93953.73463.7133200.24870.23820.237838.235737.210037.2149100.38360.38140.3823379.2445376.0303376.0419Table 3
Test II: best RREs for CG with reblurring varying the SNR (SNR=∞ means 0% of noise) within 200 iterations.
SNRPeriodicReflectiveAntireflective∞0.24150.15290.0604500.24160.15300.0616400.24180.15320.0737300.24300.15750.1052200.25740.18380.1605100.36890.33290.3289To 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 Figure5, 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 L2 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.Restored images with SNR=40.
(a)(b)(c)Residuals divided by the variance of the noise for the restored images with SNR=40.
(a)(b)(c)Two convergence histories, that is, the RREs at any iteration, are plotted in Figure7 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 L2 norm of the error, we stop the CG when the L2 norm of the residual becomes lower than the L2 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.Test I: RREs at each CG iteration with different BCs: dotted line: Periodic, dashed line: Reflective, solid line: Antireflective.
(a)(b)Discrepancy principle: images restored at the optimum value iteration and at which computed with the discrepancy principle.
(a)(b)(c)(d)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).



---

*Source: 241467-2010-09-21.xml*