Iterative Signal Processing in CommunicationsView this Special Issue
Review Article | Open Access
Marco Donatelli, Stefano Serra-Capizzano, "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
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.
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 two-level Toeplitz structure (i.e., a two-level Toeplitz matrix or in a different language a block Toeplitz matrix with Toeplitz blocks). The choice of the BCs influences small sections of by a low-rank plus low-norm term. However, these correction matrices have a substantial impact in two important directions: (a)precision of the reconstruction especially close to the boundaries (presence of ringing effects); (b)cost of the computation for recovering the “true” image from the blurred one with or without noise.
On the other hand, the involved correction matrices do not modify significantly the spaces of ill-conditioning (where the coefficient matrix shows small eigenvalues). This happen because the global matrix is of convolution type. Therefore, the eigenvectors look globally like the Fourier vectors (they are exactly the Fourier vectors when periodic BCs are imposed). Conversely, the small rank matrices induced by the chosen BCs are not generic since they are necessarily localized in space. More specifically, in one dimension there exists , being the size of the matrix, such that the low-rank correction belongs to the span of canonical matrices with and/or , 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 , pp. 211–220) where we assume that the values of the image outside the domain of consideration are zero. The corresponding blurring matrix is two-level Toeplitz, which is known to be computationally expensive to invert when direct solvers are used (see for instance the review ): the fastest available direct Toeplitz solver is of operations for an -by- two-level Toeplitz matrix . 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 . Therefore, by the convergence analysis of Axelsson and coworkers (see e.g., the pioneering paper ), we cannot expect optimal convergence (i.e., a convergence rate independent of the matrix size) both for Krylov type methods and for classical stationary solvers (with the partial exception of multigrid type techniques, see [8, 9] for a recent adaptation in the context of image deblurring). In image restoration problems, we consider regularized problems so that the previous statement is less relevant; however, for linear systems associated to the Tikhonov regularization, the optimality could guarantee a convergence speed independent of the regularizing parameter (for a discussion on this kind of robustness see [9, 10]). Moreover, an image has boundaries, and if the true image is not close to zero at the boundaries, it means that Dirichlet BCs may introduce an artificial discontinuity, which in turn results in the well-known Gibbs phenomenon (observed as a ringing effect in the reconstructed image). Therefore, we face substantial problems both with respect to (a) and (b). A possibility for drastically reducing the computational cost is to impose periodic BCs since this implies that is a two-level circulant matrix which can be efficiently inverted by means of few fast Fourier transforms (FFTs). In this case, we have a good solution to problem (b) (see Gonzalez and Woods,  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.  (see also Strang,  page 145, and Chan et al. ). 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,  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 , 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  in which new BCs are proposed that can further reduce ringing effects and that, under the assumption of symmetry of the PSF, leads to a matrix belonging to the sine algebra of type I for which fast (real) algorithms are available. More precisely, we consider the use of antireflective BCs (AR-BCs) where the data outside the domain of consideration are taken as an antireflection of the data inside. Consider the original signal 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 ) 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 ). The latter approach is computationally interesting since a DCT III requires only real operations and is about twice as fast as the FFT (see , 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  the BCs associated with all the known trigonometric algebras: the interesting facts are two, first some matrix algebras do not lead to any boundary condition, second the highest precision is reached by the classical DCT III algebra and by the classical sine transform algebra of type I which ensures only the continuity of the signal. Therefore, in this sense we can claim that among the known algebras the antireflective one is that related to the most precise reconstructions.
The paper is organized as follows. In the Section 2 we examine the antireflective BCs, the related algebra, and the related transforms also from a computational viewpoint. At the end of Section 2 we briefly consider the multilevel extension while in Section 3 we study some regularization techniques specifically adapted to the features of the antireflective algebra. Finally, Section 4 contains numerical experiments both with Tikhonov like solvers (with reblurring) and with Krylov techniques with early termination as stopping criterion. A final section of conclusions and future problems end the paper.
2. The Algebra of AR Matrices
In this section we describe the AR-BCs and the algebra , , of matrices arising from the imposition of AR-BCs.
2.1. The Antireflective BCs
When defining the antireflective boundary conditions, we assume that the data outside are an antireflection of the data inside . More precisely, if is a point outside the domain and is the closest boundary point, then we have and the quantity is approximated by so that we impose 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 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  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 .
2.2. The Algebra
Let be the type I sine transform matrix of order (see ) with entries It is known that the real matrix is orthogonal and symmetric (). For any -dimensional real vector , the matrix-vector multiplication (DST-I transform) can be computed in real operations by using the algorithm FST-I. 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 DST-I transform to the first column of and, in addition, any matrix in is uniquely determined by its first column.
Now we report a characterization of the class, which is important for analyzing the structure of AR-matrices. Let us define the shift of any vector as . According to a Matlab-like notation, we define to be the -by- symmetric Toeplitz matrix whose first column is and to be the -by- Hankel matrix whose first and last column are and , respectively. Every matrix of the class (9) can be written as (see ) where and is the flip matrix. The structure defined by (10) means that matrices in the class are special instances of Toeplitz-plus-Hankel matrices.
Moreover, the eigenvalues of the matrix in (10) are given by the cosine function evaluated at the grid points vector , 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  where this notation was originally proposed).
2.3. The AR-Algebras
Let be a real-valued 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 one-level 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 , with the above definition of the operator , we have
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 . 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 , 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  it is proved that .
2.4. Eigenvalues and Eigenvectors of AR-BC Matrices
We first describe the spectrum of AR-BC matrices, under the usual mild degree condition (i.e., the PSF has finite support), with symmetric, normalized PSFs. Then we describe the eigenvector structure and we introduce the AR-transform.
The spectral structure of any AR-BC matrix, with symmetric PSF , is concisely described in the following result.
Theorem 1 (see ). Let the blurring function (PSF) be symmetric (i.e., ), normalized, and satisfying the usual degree condition. As a consequence the eigenvalues of the AR-BCs blurring matrix in (14), , 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 AR-BCs 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 one-level 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) AR-transform 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 one-dimensional 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
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 real-valued cosine polynomial, then its Fourier coefficients form a real -dimensional tensor which is strongly symmetric, that is symmetric with respect to every direction, that is, for all . In addition, , , can be written as a linear combination of independent terms of the form where any is a nonnegative integer. Therefore, we define where denotes Kronecker product, and we force for every real and and for every -variate real-valued 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  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 vector-space 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 real-valued 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 . 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 , of Theorem 4, and of the previous tensor arguments).
Theorem 5 ( Jordan Canonical Form). The AR-BCs blurring matrix , obtained when using a strongly symmetric -dimensional mask such that 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 AR-BCs applied separately in every direction, when is the multivariate cosine symbol coming from the -D tensor mask defining the shift-invariant -dimensional blurring operator.
3. Regularization by Reblurring
When the observed signal (or image) is noise-free, then there is a substantial gain of the reflective boundary conditions (R-BCs) with respect to both the periodic and zero Dirichlet BCs and, at the same time, there is a significant improvement of the AR-BCs with regard to the R-BCs (see [17, 30]). Since the deconvolution problem is ill posed (components of the solution related to high frequency data errors are greatly amplified) regardless of the chosen BCs, it is evident that we have to regularize the problem. Two classical methods, that is, Tikhonov regularization, with direct or iterative solution of the Tikhonov linear system, and regularization iterative solvers, with early termination, for normal equations (CG  or Landweber method ) can be used. We observe that in both the cases, the coefficient matrix involves a shift of and that the righthand-side is given by . Quite surprisingly, the AR-BCs may be spoiled by this approach at least for and if we compute explicitly the matrix and the vector , see : 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 ). The reason of this fact relies upon the properties of the matrix , and we give some insights in the following.
3.1. The Reblurring Operator
A key point is that, for zero Dirichlet, periodic and reflective BCs, when the kernel is symmetric, the matrix is still a blurring operator since , while, in the case of the AR-BCs matrix, cannot be interpreted as a blurring operator. A (normalized) blurring operator is characterized by nonnegative coefficients such that every row sum is equal to (and it is still acceptable if it is not higher than ): in the case of with AR-BCs 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, Shermann-Morrison formulae (which by the way can be numerically unstable ). Dealing with higher dimensions, the scenario is even worse , 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 AR-BCs for strongly symmetric PSFs), we replace by , where is the matrix obtained by imposing the current BCs to the center-flipped PSF (i.e., in the 2D case, to the PSF rotated by 180 degrees).
For the sake of simplicity we first consider a strongly symmetric PSF so that the associated normal equations can be read as , whenever zero Dirichlet, periodic or reflective BCs are considered. Therefore, the observed image is reblurred. The reblurring is the key of the success of classical regularization techniques (Tikhonov or CG, Landweber for the normal equations) since also the noise is blurred and this makes the contribution of the noise less evident. We notice that the two systems and are algebraically equivalent if is invertible: in any case, if is also positive definite, the first represents the minimization of the functional while the second represents the minimization of the functional so that the first can be considered the blurred version of the second and in fact the first approach is uniformly better than the second. On these grounds, in the case of AR-BCs, since , we can replace by which is again a low-pass filter (see ). 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 ). In the discrete 2D case, the correlation performs the same operation of the convolution, but rotating the mask (the PSF in our case) of 180 degrees. Moreover, we note that in the continuous case over an infinite domain, the correlation and the adjoint are exactly the same operation, provided that the convolution kernel is real. Indeed, let be the convolution operator with shift-invariant kernel , then . 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 AR-BCs 2D restorations, we propose to replace by (the reblurring matrix), where is the matrix obtained by imposing the BCs to the PSF rotated by 180 degrees. Using Matlab notation, if is a PSF, its 180 degrees rotated version is , 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 low-pass 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 AR-BCs. 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 least-squares solution, (2)to formulate a least-squares solution on the infinite dimensional problem, and then to impose the BCs to the two infinite dimensional operators and , separately.
A third possibility is to formulate a least-squares solution on the infinite dimensional problem, and then to impose the BCs to this minimum problem: a difficulty in this case is that, even without noise, the resulting system is not equivalent in an algebraic sense to the original equations . In the first case we resort to the normal equations in the finite dimensional space. On the contrary, in the second case we apply the BCs to and in the infinite dimensional normal equations (where the adjoint operator performs a correlation operation) and then we obtain (34). More precisely, the discretization of and in the continuous equation with any fixed BCs gives (34).
3.2. Linear Algebra and Computational Issues
We note that in the 1D case . In the -dimensional case, let be the partial dimensions of the problem, whose total size is . 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 , 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 (D-BCs) and periodic BCs (P-BCs) the reblurring approach is exactly equal to the classical normal equations approach, since in these two cases the corresponding blurring matrix is multilevel Toeplitz: indeed, concerning P-BCs, we notice that the resulting multilevel circulant structure is a special instance of the multilevel Toeplitz case. Unfortunately, in the case of Hankel matrices (or multilevel mixed Toeplitz-Hankel matrices with at least one Hankel level) this is no longer true in general. However, a sufficient and necessary condition to have is (or equivalently ), which is a multilevel antidiagonal symmetry called persymmetry. Therefore in the case of R-BCs, where the matrix involves nested Hankel parts, in general , while only when the PSF is strongly symmetric since in this case . Dealing with the AR-BCs, the situation is even more involved, since also for strongly symmetric PSFs, owing to the low-rank correction term. Hence, we can state that the reblurring is a new proposal not only for the AR-BCs, but also for all those BCs for which . As a nontrivial and unexpected example, it is important to stress that the imposition of R-BCs with nonstrong symmetric PSFs implies that is, .
We provide now a computational motivation for the choice of using as an alternative to is the usual operation which has to be implemented to perform the adjoint operation in the Fourier domain. Indeed, the convolution with prescribed BCs can be implemented by first enlarging the image according to the considered BCs and then by computing the matrix vector product by a simple circular convolution operation, see . 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 column-stacked 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 . We have added the AR-BCs in such a toolbox for the matrix vector product, suitable for iterative regularizing methods. This code has been used for the numerical tests of Section 3 and it is downloadable from the homepage “http://scienze-como.uninsubria.it/mdonatelli/software.html”.
3.3. Filtering Methods for AR-BCs Matrices
As mentioned in Section 1, regardless of the imposed BCs, matrices that arise in signal and image restoration are typically severely ill conditioned, and regularization is needed in order to compute a stable approximation of the solution of (1). A class of regularization methods is obtained through spectral filtering [38, 39]. Specifically, if the spectral decomposition of is 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 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 where the problem dependent regularization parameters and must be chosen . Several techniques can be used to estimate appropriate choices for the regularization parameters when the SVD is used for filtering (i.e., are the singular values), including generalized cross validation (GCV), L-curve, and the discrepancy principle [38, 40, 41].
In our case, the notation in (44) defines a slight abuse of notation, because the eigenvalues are not the singular values: in fact the Jordan canonical form (CF) in (24) is different from the singular value decomposition (SVD), since the transform is not orthogonal (indeed it is a rank- 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 ), 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 .
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 component-wise.
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 matrix-vector 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 one-dimensional periodic and reflective BCs. The eigenvalue obtained as a sampling of the symbol at a grid-point close to zero, that is, close to the maximum point of , has an associated eigenvector that corresponds to low-frequency (signal space) information, while the eigenvalue obtained as a sampling of the symbol at a grid-point far away from zero (and, in particular, close to ), has an associated eigenvector that corresponds to high-frequency (noise space) information. Concerning antireflective BCs, the same situation occurs when dealing with the frequency eigenvectors , . 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 , that is, the maximum and the infinity norm of , and by the fact that AR-BCs are more precise than other classical BCs.
Remark 8. For and with reference to the previous sections, we have proved that, thanks to the definition of a (fast) AR-transform, it is possible to define a truncated spectral decomposition which is computationally very effective and, surprisingly enough, quite competitive when compared with the celebrated but costly truncated SVD in terms of restoration quality. However, we are well aware that the real challenge is represented by a general extension to the multidimensional setting. Indeed, except for more involved multi-index notations, all the techniques are plainly generalized in the multilevel setting, maintaining a cost proportional to three -level FSTs of size , 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 , where we have to compute the eigenvalues in the right order. For this task we refer to , where an algorithm is proposed and studied: more specifically the related procedure in  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 ). 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 ill-conditioned 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 ill-posed problems, which consists in so-called smoothness assumptions, the most simple one being as follows.
Assumption. Let be itself a blurred version of a continuous signal , that is, 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 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  suitably extended for dealing with AR-BCs. The relative restoration error (RRE) is , where is the computed approximation of the true image . The signal-to-noise ratio (SNR) is computed as , where is the blurred image without noise and is the noise vector .
4.1. Reblurring and Tikhonov Regularization
Let us begin with an example illustrating the approach discussed in Section 3.3 for a 2-dimensional imaging problem. We do not take an extensive comparison of the AR-BCs with other classic BCs, like periodic or reflective, since the topic and related issues have been already widely discussed in several works (see e.g., [19, 33, 46]), where the advantage on some classes of images, in terms of the restored image quality, of the application of AR-BCs has been emphasized. Here we present only a 2D image deblurring example with Gaussian blur and various levels of white Gaussian noise.
The true and the observed images are in Figure 1, where the observed image is affected by a Gaussian blur and 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 .
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 AR-BCs 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 R-BCs and AR-BCs when the PSF is not necessarily strongly symmetric. Indeed, in the case of D-BCs and P-BCs, the reblurring approach is equal to the classical normal equations, while, in the case of R-BCs and AR-BCs, this is no longer true in general.
We provide two problems using iterative regularization by CG.
Test I: Cameraman
The first test is reported in Figure 3. The true image is a cameraman and the PSF is associated with a Gaussian distribution in the square domain , with variance equal to four. We add a Gaussian noise.
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.
We show the results corresponding only to P-BCs, R-BCs, and AR-BCs. For shortness, we do not report the reconstructions coming from D-BCs, since the related restorations are usually not better than those with P-BCs.
Tables 2 and 3 show the best RREs for various levels of noise. In Table 2 we also report the norm of the residuals, that is, , where is the observed image, is the computed approximation, and is the coefficient matrix constructed according to the considered BCs: the latter measure is the sum of square errors and it represents, up to the scaling of the variance, the statistical measure of the error. As already pointed out, the choice of the BCs is important mainly for low levels of noise, that is, for high values of SNR. Indeed, in the last row of these tables (SNR ), 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.
To emphasize the quality of the restored images, we consider the reconstruction in the case of Test I for a fixed SNR equal to 40. In Figure 5, we report the restored images and in Figure 6 the residuals of the computed solutions for each pixel divided by the variance of the noise. The last one should have a normal distribution in the case of a good restoration since we add a Gaussian noise. In Figure 5 is evident the reduction of the ringing effects passing from P-BCs to R-BCs (ringings such as the horizontal white line on the top and the horizontal black line on the bottom disappear) and from R-BCs to AR-BCs (ringings such as the two vertical white lines in the bottom left disappear). Indeed, in the same figure we note a higher level of detail in the case of AR-BCs, especially concerning the face of the cameramen. The image restored with R-BCs is smoother when compared with the one restored with AR-BCs, where we can see the “freckles” effect, typical of the norm restoration. Indeed the CGLS method computes the least-squares solution that is well known to be affected by such a phenomenon . 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 , 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.
Two convergence histories, that is, the RREs at any iteration, are plotted in Figure 7 for two different values of the SNR. It should be stressed that the AR-BCs give the best results and the lowest level of RRE. Such behavior is again more evident considering Test II. Indeed, in such case the RREs with AR-BCs continue to decrease even after the first 200 iterations. On the other hand, the restorations with R-BCs start to deteriorate after the very first iterations. In addition, we notice that AR-BCs curves are in general quite flat. This is a very useful feature since the estimation of the optimal stopping value, which is well known to be a crucial and difficult task, can be done with low precision. Indeed, in order to stress the applicability of the AR-BCs to real problems, we consider for the Test I the discrepancy principle widely used with iterative methods . Since we know the norm of the error, we stop the CG when the norm of the residual becomes lower than the norm of the error. Such a criterion seems to work quite well for the AR-BCs as we can see in Figure 8. The restored images are good enough with respect to the optimal solution and also the stopping iteration, at least in this example, is close to the optimal one. On the other hand, such criterion is not always effective for the other BCs in this example. For instance, the stopping iteration for R-BCs is greater than 1000 in the case of SNR = 50 and it is 13 for SNR = 30. However, an analysis of the stopping criterion, in connection with AR-BCs, should be further investigated in the future.
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.
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).
The paper was partially supported by MIUR 2008, Grant number 20083KLJEZ.
- H. Andrews and B. Hunt, Digital Image Restoration, Prentice-Hall, 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.
- 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.
- S. Serra-Capizzano and E. Tyrtyshnikov, “Any circulant-like preconditioner for multilevel matrices is not superlinear,” SIAM Journal on Matrix Analysis and Applications, vol. 21, no. 2, pp. 431–439, 2000.
- S. Serra-Capizzano, “Matrix algebra preconditioners for multilevel Toeplitz matrices are not superlinear,” Linear Algebra and Its Applications, vol. 343/344, pp. 303–319, 2002.
- D. Noutsos, S. Serra-Capizzano, and P. Vassalos, “Matrix algebra preconditioners for multilevel Toeplitz systems do not insure optimal convergence rate,” Theoretical Computer Science, vol. 315, no. 2-3, pp. 557–579, 2004.
- 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.
- G. Fiorentino and S. Serra-Capizzano, “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.
- M. Donatelli, “A multigrid for image deblurring with Tikhonov regularization,” Numerical Linear Algebra with Applications, vol. 12, no. 8, pp. 715–729, 2005.
- R. H. Chan, M. Donatelli, S. Serra-Capizzano, and C. Tablino-Possio, “Application of multigrid techniques to image restoration problems,” in Advanced Signal Processing Algorithms, Architec