Review Article  Open Access
Marco Donatelli, Stefano SerraCapizzano, "Antireflective Boundary Conditions for Deblurring Problems", Journal of Electrical and Computer Engineering, vol. 2010, Article ID 241467, 18 pages, 2010. https://doi.org/10.1155/2010/241467
Antireflective Boundary Conditions for Deblurring Problems
Abstract
This survey paper deals with the use of antireflective boundary conditions for deblurring problems where the issues that we consider are the precision of the reconstruction when the noise is not present, the linear algebra related to these boundary conditions, the iterative and noniterative regularization solvers when the noise is considered, both from the viewpoint of the computational cost and from the viewpoint of the quality of the reconstruction. In the latter case, we consider a reblurring approach that replaces the transposition operation with correlation. For many of the considered items, the antireflective algebra coming from the given boundary conditions is the optimal choice. Numerical experiments corroborating the previous statement and a conclusion section end the paper.
1. Introduction
Formation of a blurred signal/image is typically modeled as a linear system: where is the true object, is the blurred object, and models the blurring process. Obtaining an accurate deblurring model (i.e., the matrix ) requires essentially two main pieces of information: (1)identification of the blur operator, called a point spread function (PSF),(2)choosing an appropriate boundary condition (BC), assuming that the observed image is always finite.
The identification of the blur operator is related to the infinite dimensional problem and it decides the essential structure of the matrix . In the spatially invariant case, due to the shift invariance of the blurring, for image deblurring we deduce a twolevel Toeplitz structure (i.e., a twolevel Toeplitz matrix or in a different language a block Toeplitz matrix with Toeplitz blocks). The choice of the BCs influences small sections of by a lowrank plus lownorm 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 illconditioning (where the coefficient matrix shows small eigenvalues). This happen because the global matrix is of convolution type. Therefore, the eigenvectors look globally like the Fourier vectors (they are exactly the Fourier vectors when periodic BCs are imposed). Conversely, the small rank matrices induced by the chosen BCs are not generic since they are necessarily localized in space. More specifically, in one dimension there exists , being the size of the matrix, such that the lowrank correction belongs to the span of canonical matrices with and/or , being the th vector of the canonical basis. In that case, we observe with being the classical Euclidean norm and hence the term is unable to modify significantly the characterization in frequencies of the eigenspaces. In actuality, due to the relation , the effect induced by any is more or less uniformly distributed in all the Fourier directions.
For image deblurring, a classical choice is to use zero (Dirichlet) BCs (see Andrews and Hunt [1], pp. 211–220) where we assume that the values of the image outside the domain of consideration are zero. The corresponding blurring matrix is twolevel Toeplitz, which is known to be computationally expensive to invert when direct solvers are used (see for instance the review [2]): the fastest available direct Toeplitz solver is of operations for an by twolevel 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 illconditioned 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 wellknown Gibbs phenomenon (observed as a ringing effect in the reconstructed image). Therefore, we face substantial problems both with respect to (a) and (b). A possibility for drastically reducing the computational cost is to impose periodic BCs since this implies that is a twolevel 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 twolevel algebra simultaneously diagonalized by the fast cosine transform (DCT III). Because a twolevel 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. 5960), inversion of these matrices is substantially faster than those obtained from zero BCs and is comparable to (or little better than) the case of periodic BCs. Thus, for circularly symmetric PSFs, requirement (b) is satisfied when using reflective BCs. Moreover, the artificial boundary discontinuities are eliminated, and therefore the ringing effects are strongly reduced. However, although reflection guarantees continuity, it generally fails to guarantee continuity of the normal derivative except for the nongeneric case where the normal derivative is zero at the boundary. Therefore if the image is smooth, then the reflective BCs only save the continuity but introduce an artificial discontinuity of the normal derivative (though requirement (b) is substantially improved with respect to periodic and Dirichlet BCs). In the image processing literature, other methods have also been proposed to assign boundary values as a local image mean, or are obtained by extrapolation methods in order to insure continuity (see Lagendijk and Biemond [16], p. 22, and the references therein). However, these techniques that meet the approximation requirement (a), generally do not satisfy the computational requirement (b) so we still have to face a difficult linear algebra problem.
Here we follow the approach in [17] in which new BCs are proposed that can further reduce ringing effects and that, under the assumption of symmetry of the PSF, leads to a matrix belonging to the sine algebra of type I for which fast (real) algorithms are available. More precisely, we consider the use of antireflective BCs (ARBCs) where the data outside the domain of consideration are taken as an antireflection of the data inside. Consider the original signal and the normalized blurring function given by The blurred signal is the convolution of and and consequently . The deblurring problem is to recover the vector given the blurring function and a blurred signal of finite length, that is, Thus the blurred signal is determined not by only, but also by and and the linear system (3) is underdetermined. To overcome this, we make certain assumptions (called BCs) on the unknown boundary data and in such a way that the number of unknowns equals the number of equations.
In terms of the objects defined so far, we recall that zero (Dirichlet) BCs means for all in (3) so that a Toeplitz structure is encountered. If we consider periodic BCs we set , for all in (3) and the matrix system in (3) becomes by circulant so that it can be diagonalized by the discrete Fourier matrix and the above system can be solved by using three FFTs (one for finding the eigenvalues of the blurring matrix and two for solving the system). For the Neumann BCs, we assume that the data outside are a reflection of the data inside (refer to [12]) so that for and for all in (3). Thus (3) becomes , where is neither Toeplitz nor circulant but it is a special by Toeplitz plus Hankel matrix which is diagonalized by the discrete cosine transform provided that the blurring function is symmetric, that is, for all in (2). It follows that the above system can be solved by using three transforms DCT III in 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. 5960), and this is true in two dimensions as well. With the help of a different Toeplitz plus Hankel structure, we establish similar results for the antireflective BCs both in one dimension and two dimensions. It is worth finally to remark that La Spina has analyzed [18] the BCs associated with all the known trigonometric algebras: the interesting facts are two, first some matrix algebras do not lead to any boundary condition, second the highest precision is reached by the classical DCT III algebra and by the classical sine transform algebra of type I which ensures only the continuity of the signal. Therefore, in this sense we can claim that among the known algebras the antireflective one is that related to the most precise reconstructions.
The paper is organized as follows. In the Section 2 we examine the antireflective BCs, the related algebra, and the related transforms also from a computational viewpoint. At the end of Section 2 we briefly consider the multilevel extension while in Section 3 we study some regularization techniques specifically adapted to the features of the antireflective algebra. Finally, Section 4 contains numerical experiments both with Tikhonov like solvers (with reblurring) and with Krylov techniques with early termination as stopping criterion. A final section of conclusions and future problems end the paper.
2. The Algebra of AR Matrices
In this section we describe the ARBCs and the algebra , , of matrices arising from the imposition of ARBCs.
2.1. The Antireflective BCs
When defining the antireflective boundary conditions, we assume that the data outside are an antireflection of the data inside . More precisely, if is a point outside the domain and is the closest boundary point, then we have and the quantity is approximated by so that we impose in (3). If we define the vector whose components are for and zero otherwise and if we define the vector whose components are for and zero otherwise, then (3) becomes where is the th vector of the canonical basis, , with denoting the dimensional flip matrix having entries if and zero otherwise, and where the matrices , and are given by Therefore, the matrices and involved in (5) take the form
We remark that the coefficient matrix in (5) is neither Toeplitz nor circulant, but it is a Toeplitz plus Hankel plus 2 rank correction matrix, where the correction is placed at the first and the last column. We will show that the linear system (5) can be reduced to an by new system whose coefficient matrix can always be diagonalized by the discrete sine transform DST I (associated with the class) matrix provided that the blurring function is symmetric, that is, for all in (2). It follows that (5) can be solved by using three FSTs in operations. This approach is computationally attractive as FST requires only real operations and is about twice as fast as the FFT and hence solving a problem with the AR BCs is twice as fast as solving a problem with the periodic BCs, has the same cost as solving a problem with the reflective BCs, and one gains one order of precision due to the continuity. Moreover, all these remarks stand in two dimensions as well and indeed an abstract treatment of the twolevel and multilevel settings is reported in Section 2.5.
Finally it is worth mentioning that the use of AR BCs has been considered by several authors (see e.g., [18–26]). In particular in [24] the mean BCs have been introduced as a slight variation of the AR BCs. The order of approximation is the same but the hidden constants are smaller so that the mean BCs are generally more precise than the AR BCs. However, the resulting matrices do not form an algebra so that we cannot define a new transform and this is a confirmation of the negative result found in [18].
2.2. The Algebra
Let be the type I sine transform matrix of order (see [27]) with entries It is known that the real matrix is orthogonal and symmetric (). For any dimensional real vector , the matrixvector multiplication (DSTI transform) can be computed in real operations by using the algorithm FSTI. Let be the space of all the matrices that can be diagonalized by : Let , then . Consequently, if we let denote the first column of the identity matrix, then the relationship implies that the eigenvalues of are given by . Therefore, the eigenvalues of can be obtained by applying a DSTI transform to the first column of and, in addition, any matrix in is uniquely determined by its first column.
Now we report a characterization of the class, which is important for analyzing the structure of ARmatrices. Let us define the shift of any vector as . According to a Matlablike notation, we define to be the by symmetric Toeplitz matrix whose first column is and to be the by Hankel matrix whose first and last column are and , respectively. Every matrix of the class (9) can be written as (see [27]) where and is the flip matrix. The structure defined by (10) means that matrices in the class are special instances of ToeplitzplusHankel matrices.
Moreover, the eigenvalues of the matrix in (10) are given by the cosine function evaluated at the grid points vector , where, and for . The matrix in (10) is usually denoted by and is called the matrix generated by the function or symbol (see the seminal paper [27] where this notation was originally proposed).
2.3. The ARAlgebras
Let be a realvalued cosine polynomial of degree and let (note that coincides with the matrix in (10)(11), when ). Hence the Fourier coefficients of are such that with if , and for we can define the onelevel operator where and It is interesting to observe that has a zero of order at least at zero (since is a cosine polynomial) and therefore is still a cosine polynomial of degree , whose value at zero is (in other words the function is well defined at zero).
As proved in [28], with the above definition of the operator , we have
(1) ,
(2) ,
for real and and for cosine functions and .
These properties allow us to define as the algebra (closed under linear combinations, product, and inversion) of matrices , with being a cosine polynomial. By standard interpolation arguments it is easy to see that can be defined as the set of matrices , where is a cosine polynomial of degree at most . Therefore, it is clear that . Moreover, the algebra is commutative thanks to item 2, since at every . Consequently, if matrices of the form 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 ARtransform and its inverse are close in rank to orthogonal linear mappings and only moderately ill conditioned.
Following the analysis given in [17], if the blurring function (the PSF) is symmetric (i.e., ), if for (degree condition), and if is normalized so that , then the structure of the antireflective blurring matrix is where , , has order and According to the brief discussion of Section 2.2, relation (15) implies that with (see (10) and (11)). Moreover in [28] it is proved that .
2.4. Eigenvalues and Eigenvectors of ARBC Matrices
We first describe the spectrum of ARBC matrices, under the usual mild degree condition (i.e., the PSF has finite support), with symmetric, normalized PSFs. Then we describe the eigenvector structure and we introduce the ARtransform.
The spectral structure of any ARBC matrix, with symmetric PSF , is concisely described in the following result.
Theorem 1 (see [28]). Let the blurring function (PSF) be symmetric (i.e., ), normalized, and satisfying the usual degree condition. As a consequence the eigenvalues of the ARBCs blurring matrix in (14), , are given by with multiplicity two and .
The proof can be easily derived by (12) which shows that the eigenvalues of are with multiplicity and those of , that is, , with multiplicity each.
Here we will determine the eigenvectors of every matrix . In particular, we show that every ARBCs matrix is diagonalizable, and we demonstrate independence of the eigenvectors from the symbol . With reference to the notation in (8)–(11), calling the th column of , and the th point of , , we have since is an eigenvector of and is the related eigenvalue. Due to the centrosymmetry of the involved matrix, if is an eigenvector of related to the eigenvalue , then the other is its flip, that is, . Let us look for this eigenvector by imposing the equality which is equivalent to seeking a vector that satisfies Since by definition of the operator (see (12) and the lines below), and because of the algebra structure of and thanks to the above relation, we deduce that the vector satisfies the relation where is the discrete onelevel Laplacian, that is, . Therefore, by (19), the solution is given by . If is invertible (as it happens for every nontrivial PSF, since the unique maximum of the function is reached at , which is not a grid point of ), then the solution is unique. Hence, independently of , we have Now we observe that the th eigenvector is unitary, , because is unitary: we wish to impose the same condition on the first and the last eigenvector. The interesting fact is that has an explicit expression. By using standard finite difference techniques, it follows that so that the first eigenvector is exactly the sampling of the function on the grid for . Its Euclidean norm is , where, for nonnegative sequences , , the relation means . In this way, the (normalized) ARtransform can be defined as
Remark 2. With the normalization condition in (21), all the columns of are unitary. However orthogonality is only partially fulfilled since it holds for the central columns, while the first and last columns are not orthogonal to each other, and neither one is orthogonal to the central columns. We can solve the first problem: the sum of the first and of the last column (suitably normalized) and the difference of the first and the last column (suitably normalized) become orthonormal, and are still eigenvectors related to the eigenvalue . However, since has only positive components and the vector space generated by the first and the last column of contains positive vectors, it follows that cannot be made orthonormal just by operating on the first and the last column. Indeed, we do not want to change the central block of since it is related to a fast real transform and hence, necessarily, we cannot get rid of this quite mild lack of orthogonality.
Remark 3. There is a suggestive functional interpretation of the transform . When considering periodic BCs, the transform of the related matrices is the Fourier transform: its th column vector, up to a normalizing scalar factor, can be viewed as a sampling, over a suitable uniform gridding of , of the frequency function . Analogously, when imposing reflective BCs with a symmetric PSF, the transform of the related matrices is the cosine transform: its th column vector, up to a normalizing scalar factor, can be viewed as a sampling, over a suitable uniform gridding of , of the frequency function . Here the imposition of the antireflective BCs can be functionally interpreted as a linear combination of sine functions and of linear polynomials (whose use is exactly required for imposing continuity at the borders).
The previous observation becomes evident in the expression of in (21). Indeed, by defining the onedimensional grid , which is a subset of , we infer that the first column of is given by , the th column of , , is given by , and finally the last column of is given by that is,
Finally, it is worth mentioning that the inverse transform is also described in terms of the same block structure since
Theorem 4 ( Jordan Canonical Form). With the notation and assumptions of Theorem 1, the ARBCs blurring matrix in (14), , coincides with where and are defined in (22) and (23), while .
2.5. Multilevel Extension
Here we provide some comments on the extension of our findings to dimensional objects with . When , is a vector, when , is a 2D array, when , is a 3D tensor and so forth.
With reference to Section 2.3 we propose a (canonical) multidimensional extension of the algebras and of the operators , : the idea is to use tensor products. If is variate realvalued cosine polynomial, then its Fourier coefficients form a real dimensional tensor which is strongly symmetric, that is symmetric with respect to every direction, that is, for all . In addition, , , can be written as a linear combination of independent terms of the form where any is a nonnegative integer. Therefore, we define where denotes Kronecker product, and we force for every real and and for every variate realvalued cosine polynomials and . It is clear that the request that is a linear operator (for , we impose this property in (26) by definition) is sufficient for defining completely the operator in the dimensional setting.
With the above definition of the operator , we have (1), (2),
for real and and for cosine functions and .
The latter properties of algebra homomorphism allows to define a commutative algebra of the matrices , with being a variate cosine polynomial. By standard interpolation arguments it is easy to see that can be defined as the set of matrices , where is a variate cosine polynomial of degree at most in the th variable for every ranging in : we denote the latter polynomial set by , with being the vector of all ones. Here we have to be a bit careful in specifying the meaning of algebra when talking of polynomials. More precisely, for the product is the unique polynomial satisfying the following interpolation condition: If the degree of plus the degree of in the th variable does not exceed , , then the uniqueness of the interpolant implies that coincides with the product between polynomials in the usual sense. The uniqueness holds also for thanks to the tensor form of the grid (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, has become an algebra, showing the vectorspace dimension equal to which coincides with that of .
Without loss of generality and for the sake of notational clarity, in the following we assume for . Thanks to the tensor structure emphasized in (25)–(26), and by using Theorem 4 for every term , , of the level extension of such a theorem easily follows. More precisely, if is a variate realvalued cosine symbol related to a dimensional strongly symmetric and normalized mask , then ( times) where is the diagonal matrix containing the eigenvalues of . The description of in dimensions is quite involved when compared with the case , implicitly reported in Theorem 1.
For a complete analysis of the spectrum of we refer the reader to [28]. Here we give details on a specific aspect. More precisely we attribute a correspondence in a precise and simple way among eigenvalues and eigenvectors, by making recourse only to the main variate symbol . Let be a column of , with or , and is the th column of , for . Let with being the generic eigenvector, that is, the generic column of . The eigenvalue related to is where for and for . Defining the dimensional grid we can evaluate the variate function on by , where arranges the entries of in a dimensional array of length along each direction according to Matlab notation. Using this notation the following compact and elegant result can be stated (its proof is omitted since it is simply the combination of the eigenvalue analysis in [28], of Theorem 4, and of the previous tensor arguments).
Theorem 5 ( Jordan Canonical Form). The ARBCs blurring matrix , obtained when using a strongly symmetric dimensional mask such that if for some (the dimensional degree condition), , coincides with where and are defined in (28) and (31).
It is worth observing that the matrix spectrally analyzed in the previous theorem is exactly the same matrix arising from the imposition of ARBCs applied separately in every direction, when is the multivariate cosine symbol coming from the D tensor mask defining the shiftinvariant dimensional blurring operator.
3. Regularization by Reblurring
When the observed signal (or image) is noisefree, then there is a substantial gain of the reflective boundary conditions (RBCs) with respect to both the periodic and zero Dirichlet BCs and, at the same time, there is a significant improvement of the ARBCs with regard to the RBCs (see [17, 30]). Since the deconvolution problem is ill posed (components of the solution related to high frequency data errors are greatly amplified) regardless of the chosen BCs, it is evident that we have to regularize the problem. Two classical methods, that is, Tikhonov regularization, with direct or iterative solution of the Tikhonov linear system, and regularization iterative solvers, with early termination, for normal equations (CG [31] or Landweber method [32]) can be used. We observe that in both the cases, the coefficient matrix involves a shift of and that the righthandside is given by . Quite surprisingly, the ARBCs may be spoiled by this approach at least for and if we compute explicitly the matrix and the vector , see [33]: more in detail, even in presence of a moderate noise and a strongly symmetric PSF, the accuracy of ARBCs restorations becomes worse in some examples than the accuracy of RBCs restorations (see [33]). The reason of this fact relies upon the properties of the matrix , and we give some insights in the following.
3.1. The Reblurring Operator
A key point is that, for zero Dirichlet, periodic and reflective BCs, when the kernel is symmetric, the matrix is still a blurring operator since , while, in the case of the ARBCs matrix, cannot be interpreted as a blurring operator. A (normalized) blurring operator is characterized by nonnegative coefficients such that every row sum is equal to (and it is still acceptable if it is not higher than ): in the case of with ARBCs the row sum of the first and of the last row can be substantially larger than . This means that modified signal may have artifacts at the borders and this seems to be a potential motivation for which a reduction of the reconstruction quality occurs.
Furthermore, the structure of the matrix is also spoiled and, in the case of images () we lose the computational cost for solving a generic system . The cost of solving such a type of linear systems is proportional to by using for example, ShermannMorrison formulae (which by the way can be numerically unstable [34]). Dealing with higher dimensions, the scenario is even worse [35], because in the dimensional setting the solution of the normal equation linear system is asymptotic to , where is the size of the matrix . In order to overcome these problems (which arise only with the most precise ARBCs for strongly symmetric PSFs), we replace by , where is the matrix obtained by imposing the current BCs to the centerflipped PSF (i.e., in the 2D case, to the PSF rotated by 180 degrees).
For the sake of simplicity we first consider a strongly symmetric PSF so that the associated normal equations can be read as , whenever zero Dirichlet, periodic or reflective BCs are considered. Therefore, the observed image is reblurred. The reblurring is the key of the success of classical regularization techniques (Tikhonov or CG, Landweber for the normal equations) since also the noise is blurred and this makes the contribution of the noise less evident. We notice that the two systems and are algebraically equivalent if is invertible: in any case, if is also positive definite, the first represents the minimization of the functional while the second represents the minimization of the functional 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 ARBCs, since , we can replace by which is again a lowpass filter (see [33]). In this way, we overcome also the computational problems.
In order to provide a general reblurring approach also in the case of nonsymmetric PSFs, we consider the correlation operation instead of the transposition (see [36]). In the discrete 2D case, the correlation performs the same operation of the convolution, but rotating the mask (the PSF in our case) of 180 degrees. Moreover, we note that in the continuous case over an infinite domain, the correlation and the adjoint are exactly the same operation, provided that the convolution kernel is real. Indeed, let be the convolution operator with shiftinvariant kernel , then . Since the PSF (and then ) is real (and then real valued), the adjoint operator is which is a correlation operator. We remark that here the convolution and the correlation use the same kernel except for the sign of the variable (i.e., vs ), and, in the 2D case, the change of sign in the variable of the kernel can be viewed as a 180 degrees rotation of the PSF mask.
By virtue of these arguments, in order to overcome the problems arising with the normal equations approach for ARBCs 2D restorations, we propose to replace by (the reblurring matrix), where is the matrix obtained by imposing the BCs to the PSF rotated by 180 degrees. Using Matlab notation, if is a PSF, its 180 degrees rotated version is , where is the flip matrix defined as if for , and zero otherwise. For a dimensional problem, is obtained by imposing the BCs to the PSF flipped with respect to the origin, or, in other words, to the new PSF where all the coefficients are flipped with respect to every variable.
In this way has the same computational properties of (it belongs to in the case of AR BCs) and it is certainly a lowpass filter. In the reblurring approach the normal equations are replaced by Furthermore, using the Tikhonov regularization in the reblurred version, we propose to use with being any discrete. differential operator with ARBCs. In general is nonsymmetric, but the asymptotic eigenvalue analysis in [28, Section ] has shown that the spectrum is clustered around the positive real interval given by the range of if is the symbol associated to the PSF. Such a cluster is strong if the decay of the PSF coefficients is fast enough, as it occurs in real world PSFs. The latter statements give a reasonable support to the applicability of the CGLS when is replaced by and, indeed, in our numerical tests such a regularizing method has always worked perfectly. From the viewpoint of the modeler, the previous considerations can be summarized in the following motivation. The image restoration problem is the restriction in the FOV of an infinite dimensional problem. We can follow two ways to design the linear system to solve: (1)to impose BCs and then to look at a leastsquares solution, (2)to formulate a leastsquares solution on the infinite dimensional problem, and then to impose the BCs to the two infinite dimensional operators and , separately.
A third possibility is to formulate a leastsquares solution on the infinite dimensional problem, and then to impose the BCs to this minimum problem: a difficulty in this case is that, even without noise, the resulting system is not equivalent in an algebraic sense to the original equations . In the first case we resort to the normal equations in the finite dimensional space. On the contrary, in the second case we apply the BCs to and in the infinite dimensional normal equations (where the adjoint operator performs a correlation operation) and then we obtain (34). More precisely, the discretization of and in the continuous equation with any fixed BCs gives (34).
3.2. Linear Algebra and Computational Issues
We note that in the 1D case . In the dimensional case, let be the partial dimensions of the problem, whose total size is . By flipping each variable, we obtain For the analysis of properties of the reblurring scheme (34) with respect to all the different BCs, we now study the discretization of the continuous operator . Let us consider the Toeplitz level matrix of partial dimensions and generating function [2], which is defined as () by means of the Fourier coefficients of Here and for , , is the matrix whose () th entry is if , and elsewhere. As it is well known for multilevel Toeplitz matrices , where is the conjugate of the function , and the Fourier coefficients of are the same of , but conjugated and flipped. Moreover, since , if is real then . This means that for Dirichlet BCs (DBCs) and periodic BCs (PBCs) the reblurring approach is exactly equal to the classical normal equations approach, since in these two cases the corresponding blurring matrix is multilevel Toeplitz: indeed, concerning PBCs, 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 ToeplitzHankel matrices with at least one Hankel level) this is no longer true in general. However, a sufficient and necessary condition to have is (or equivalently ), which is a multilevel antidiagonal symmetry called persymmetry. Therefore in the case of RBCs, where the matrix involves nested Hankel parts, in general , while only when the PSF is strongly symmetric since in this case . Dealing with the ARBCs, the situation is even more involved, since also for strongly symmetric PSFs, owing to the lowrank correction term. Hence, we can state that the reblurring is a new proposal not only for the ARBCs, but also for all those BCs for which . As a nontrivial and unexpected example, it is important to stress that the imposition of RBCs with nonstrong symmetric PSFs implies that is, .
We provide now a computational motivation for the choice of using as an alternative to is the usual operation which has to be implemented to perform the adjoint operation in the Fourier domain. Indeed, the convolution with prescribed BCs can be implemented by first enlarging the image according to the considered BCs and then by computing the matrix vector product by a simple circular convolution operation, see [37]. More precisely, let and be two matrices such that represents an image and is the discrete 2D PSF, , which leads to the matrix blurring . By using the Matlab notation (i.e., the vector is the columnstacked version of ), the product in the 2D case can be implemented (1)by using an enlarged image , which is the image extended at the boundaries according to the imposed BCs, (2)computing , where the symbol “” denotes the circular convolution operator ( should be zero padded to have the same size of ), (3)and then taking the inner part of the result having the same size of .
The circular convolution can be computed using the 2D discrete Fourier transform (DFT2) and its inverse (IDFT2), since we have where “” is the componentwise product. If denotes the block circulant matrix with circulant blocks, of block size with blocks of size and generating function , then (39) represents the same operation as , where (according to a 2D ordering). Conversely, it is well known that the operation corresponding to the product with the adjoint operator, in the Fourier domain gives rise to where the overline symbol denotes the complex conjugation. As a result, since the transform is equivalent to the transform (because ), and since , if , then it follows that . Therefore, for any of the considered BCs, the inner part of is exactly . Here it is worthwhile to specify exactly what we mean for inner part: if the vector is partitioned in blocks of size , , then for inner part we mean that we are excluding the first and the last blocks and, in any of the remaining blocks, we are deleting the first and the last entries. More generally, if the PSF is arbitrary (e.g. nonsymmetric) that is, the nonzero coefficients of the PSF have first index belonging to and second index in the range , then we have to delete the first and the last blocks and, in any of the other blocks, we have to exclude the first and the last entries.
Since the DFT and its inverse can be computed in arithmetic operations using FFTs, the previous approach is implemented in the Matlab toolbox RestoreTools [37]. We have added the ARBCs 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://scienzecomo.uninsubria.it/mdonatelli/software.html”.
3.3. Filtering Methods for ARBCs Matrices
As mentioned in Section 1, regardless of the imposed BCs, matrices that arise in signal and image restoration are typically severely ill conditioned, and regularization is needed in order to compute a stable approximation of the solution of (1). A class of regularization methods is obtained through spectral filtering [38, 39]. Specifically, if the spectral decomposition of is with , then a spectral filter solution is given by where are filter factors that satisfy 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 lowfrequency 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 wellknown filtering methods, truncated spectral value decomposition (TSVD) and Tikhonov regularization, are where the problem dependent regularization parameters and must be chosen [39]. Several techniques can be used to estimate appropriate choices for the regularization parameters when the SVD is used for filtering (i.e., are the singular values), including generalized cross validation (GCV), Lcurve, and the discrepancy principle [38, 40, 41].
In our case, the notation in (44) defines a slight abuse of notation, because the eigenvalues are not the singular values: in fact the Jordan canonical form (CF) in (24) is different from the singular value decomposition (SVD), since the transform is not orthogonal (indeed it is a rank correction of a symmetric orthogonal matrix). Therefore, note that the use of 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 ARtransform (21), its inverse (23), and a fast algorithm for computing the eigenvalues described in [28].
Algorithm 6. (1)
(2) , where are the eigenvalues of that can be computed by a fast sine transform (FST).
(3) /, where the dot operations are componentwise.
(4)
The product can be clearly computed in a fast and stable way by one FST. Indeed for all where in Matlab notation is the vector with components indexed from 2 to . A similar strategy can be followed for computing the matrixvector product . Instead of there is and instead of there is . Recalling that the two vectors and can be explicitly computed obtaining , for and .
Remark 7. As discussed in Remark 3, there is a natural interpretation in terms of frequencies when considering onedimensional periodic and reflective BCs. The eigenvalue obtained as a sampling of the symbol at a gridpoint close to zero, that is, close to the maximum point of , has an associated eigenvector that corresponds to lowfrequency (signal space) information, while the eigenvalue obtained as a sampling of the symbol at a gridpoint far away from zero (and, in particular, close to ), has an associated eigenvector that corresponds to highfrequency (noise space) information. Concerning antireflective BCs, the same situation occurs when dealing with the frequency eigenvectors , . The other two exceptional eigenvectors generate the space of linear polynomials and therefore they correspond to lowfrequency information: this intuition is well supported by the fact that the related eigenvalue is , that is, the maximum and the infinity norm of , and by the fact that ARBCs are more precise than other classical BCs.
Remark 8. For and with reference to the previous sections, we have proved that, thanks to the definition of a (fast) ARtransform, 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 multiindex notations, all the techniques are plainly generalized in the multilevel setting, maintaining a cost proportional to three level FSTs of size , and the key tool is the simplified eigenvalueeigenvector correspondence concisely indicated in Theorem 5. In reality, regarding the previous Algorithm 6 the only difficult task is the computation in step , where we have to compute the eigenvalues in the right order. For this task we refer to [28], where an algorithm is proposed and studied: more specifically the related procedure in [28] is based on a single level FST of size 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 where
Note that and differ from the corresponding eigenvectors and used in Section 2.4, but they are a linear combination of them. They have been chosen here to form an orthonormal basis of the grid samples of all linear polynomials and this property will be useful in the following.
According to Theorem 4 the spectral decomposition of can now be written as where the diagonal entries of are given by Here we prove that the reblurring approach for Tikhonov regularization in (35) where is a regularization method. For a complete analysis we need to compute the SVD of .
We use the notation if the two vectors and depend on and for each entry of and the corresponding entry of there holds as .
Theorem 9 (see [43]). The two dominant singular values of are given by where is, in fact, a strict upper bound, and the two minimal singular values are given by respectively. Fix , , , and . The corresponding right singular vectors are and the left singular vectors are 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
Remark 11. It is important to note that the illconditioned subspace of has dimension two, independent of , since has two singular values that decay like while all others are between one and two. Also, only amplifies vectors that fail to be orthogonal to . 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 illposed problems, which consists in socalled smoothness assumptions, the most simple one being as follows.
Assumption. Let be itself a blurred version of a continuous signal , that is, where satisfies the same BCs as (i.e., periodic, reflective, or antireflective ones).
On the grounds of Assumption 12 we may therefore assume that with a moderate bound where . Since the observed object is usually affected by noise, instead of the blurred object we have to work with the blurred and noisy object , where . In this way, for (35) becomes
Theorem 13. Let the exact solution of (1) satisfy (56) with (57). Furthermore the total error of the reblurring strategy with AR BCs satisfies for , where the constant in the notation is independent of the dimension .
Note that the upper bound from Theorem 13 is the same as for Tikhonov regularization with reflective or periodic BCs; only the constant hidden in the notation may be somewhat larger for the reblurring strategy.
4. Numerical Results
The section is devoted to numerical experiments, with reference to the Tikhonov regularization in the reblurred version and to the classical conjugate gradient (CG) regularization with early termination. In both cases the use of antireflective BCs improves the quality of the restored image, without penalizing the computational cost. Another promising approach not discussed here both from the viewpoint of the quality and of the complexity is that based on a regularized version of multigridtype 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 ARBCs. The relative restoration error (RRE) is , where is the computed approximation of the true image . The signaltonoise ratio (SNR) is computed as , where is the blurred image without noise and is the noise vector [32].
4.1. Reblurring and Tikhonov Regularization
Let us begin with an example illustrating the approach discussed in Section 3.3 for a 2dimensional imaging problem. We do not take an extensive comparison of the ARBCs 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 ARBCs has been emphasized. Here we present only a 2D image deblurring example with Gaussian blur and various levels of white Gaussian noise.
The true and the observed images are in Figure 1, where the observed image is affected by a Gaussian blur and noise. We compare the ARBCs 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 ARBCs with respect to reflective BCs. Restored images in Figure 2 are obtained with the minimum relative restoration error varying several values of the regularization parameter .
(a)
(b)
(a)
(b)
From Table 1, we note that for the noise case, all of the approaches give comparable restorations. On the other hand, decreasing the noise, that is, passing to and then to noise, the ARBCs improve the restoration while the reflective BCs are not able to do that, due to the barrier of the ringing effects.

4.2. Reblurring and CG Regularization
In the following tests, the reblurring strategy will be applied with RBCs and ARBCs when the PSF is not necessarily strongly symmetric. Indeed, in the case of DBCs and PBCs, the reblurring approach is equal to the classical normal equations, while, in the case of RBCs and ARBCs, this is no longer true in general.
We provide two problems using iterative regularization by CG.
Test I: Cameraman
The first test is reported in Figure 3. The true image is a cameraman and the PSF is associated with a Gaussian distribution in the square domain , with variance equal to four. We add a Gaussian noise.
(a)
(b)
(c)
Test II: Astronomy
We are dealing with a nonsymmetric experimental PSF developed by US Air Force Phillips Laboratory, Lasers and Imaging Directorate, Kirtland Air Force Base, New Mexico, widely used in literature (see e.g. [12, 47]). The true object is the image of Saturn in Figure 4; a Poissonian noise is added, as it is customary when dealing with astronomical images.
(a)
(b)
(c)
We show the results corresponding only to PBCs, RBCs, and ARBCs. For shortness, we do not report the reconstructions coming from DBCs, since the related restorations are usually not better than those with PBCs.
Tables 2 and 3 show the best RREs for various levels of noise. In Table 2 we also report the norm of the residuals, that is, , where is the observed image, is the computed approximation, and is the coefficient matrix constructed according to the considered BCs: the latter measure is the sum of square errors and it represents, up to the scaling of the variance, the statistical measure of the error. As already pointed out, the choice of the BCs is important mainly for low levels of noise, that is, for high values of SNR. Indeed, in the last row of these tables (SNR ), 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 ARBCs 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.


To emphasize the quality of the restored images, we consider the reconstruction in the case of Test I for a fixed SNR equal to 40. In Figure 5, we report the restored images and in Figure 6 the residuals of the computed solutions for each pixel divided by the variance of the noise. The last one should have a normal distribution in the case of a good restoration since we add a Gaussian noise. In Figure 5 is evident the reduction of the ringing effects passing from PBCs to RBCs (ringings such as the horizontal white line on the top and the horizontal black line on the bottom disappear) and from RBCs to ARBCs (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 ARBCs, especially concerning the face of the cameramen. The image restored with RBCs is smoother when compared with the one restored with ARBCs, where we can see the “freckles” effect, typical of the norm restoration. Indeed the CGLS method computes the leastsquares solution that is well known to be affected by such a phenomenon [39]. When passing to the RBCs, 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 ARBCs 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 ARBCs, while, also with the RBCs, 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 RBCs is not good enough as model at least for this example. On the other hand, with the ARBCs, 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 ARBCs.
(a)
(b)
(c)
(a)
(b)
(c)
Two convergence histories, that is, the RREs at any iteration, are plotted in Figure 7 for two different values of the SNR. It should be stressed that the ARBCs 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 ARBCs continue to decrease even after the first 200 iterations. On the other hand, the restorations with RBCs start to deteriorate after the very first iterations. In addition, we notice that ARBCs 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 ARBCs to real problems, we consider for the Test I the discrepancy principle widely used with iterative methods [31]. Since we know the norm of the error, we stop the CG when the norm of the residual becomes lower than the norm of the error. Such a criterion seems to work quite well for the ARBCs 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 RBCs 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 ARBCs, should be further investigated in the future.
(a)
(b)
(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 RBCs and ARBCs. Therefore, it is possible that the applicability of the CG to (34) is a general property, which does not depend on the particular choice of the BCs.
5. Conclusions
In this contribution we have dealt with the use of antireflective BCs for deblurring problems where the considered issues have been: the precision of the reconstruction when the noise is not present, the linear algebra related to these BCs, the computational costs, and the reconstruction quality associated with iterative and noniterative regularizing solvers when the noise is considered. For many of the considered items, the antireflective algebra coming from the given BCs is the optimal choice: for instance in the work of La Spina it is proven that no one of the trigonometric algebras can be associated with BCs of the same precision as the antireflective. Numerical experiments corroborating the previous statements have been reported and discussed: in this direction it remains an open problem to understand why the CG works without any numerical problem even if the antireflective structure is non symmetric (an event not normal as emphasized by the Jordan decomposition).
Acknowledgment
The paper was partially supported by MIUR 2008, Grant number 20083KLJEZ.
References
 H. Andrews and B. Hunt, Digital Image Restoration, PrenticeHall, Englewood Cliffs, NJ, USA, 1977.
 R. H. Chan and M. Ng, “Conjugate gradient methods for Toeplitz systems,” SIAM Review, vol. 38, no. 3, pp. 427–482, 1996. View at: Google Scholar
 N. Kalouptsidis, G. Carayannis, and D. Manolakis, “Fast algorithms for block Toeplitz matrices with Toeplitz entries,” Signal Processing, vol. 6, no. 1, pp. 77–81, 1984. View at: Google Scholar
 S. SerraCapizzano and E. Tyrtyshnikov, “Any circulantlike preconditioner for multilevel matrices is not superlinear,” SIAM Journal on Matrix Analysis and Applications, vol. 21, no. 2, pp. 431–439, 2000. View at: Google Scholar
 S. SerraCapizzano, “Matrix algebra preconditioners for multilevel Toeplitz matrices are not superlinear,” Linear Algebra and Its Applications, vol. 343/344, pp. 303–319, 2002. View at: Google Scholar
 D. Noutsos, S. SerraCapizzano, and P. Vassalos, “Matrix algebra preconditioners for multilevel Toeplitz systems do not insure optimal convergence rate,” Theoretical Computer Science, vol. 315, no. 23, pp. 557–579, 2004. View at: Publisher Site  Google Scholar
 O. Axelsson and G. Lindskog, “On the rate of convergence of the preconditioned conjugate gradient method,” Numerische Mathematik, vol. 48, no. 5, pp. 499–523, 1986. View at: Publisher Site  Google Scholar
 G. Fiorentino and S. SerraCapizzano, “Multigrid methods for symmetric positive definite block toeplitz matrices with nonnegative generating functions,” SIAM Journal of Scientific Computing, vol. 17, no. 5, pp. 1068–1081, 1996. View at: Google Scholar
 M. Donatelli, “A multigrid for image deblurring with Tikhonov regularization,” Numerical Linear Algebra with Applications, vol. 12, no. 8, pp. 715–729, 2005. View at: Publisher Site  Google Scholar
 R. H. Chan, M. Donatelli, S. SerraCapizzano, and C. TablinoPossio, “Application of multigrid techniques to image restoration problems,” in Advanced Signal Processing Algorithms, Architec