A Fast Algorithm for Deblurring Models with Neumann Boundary Conditions.
Blur removal is an important problem in signal and image processing. The blurring matrices obtained by using the zero boundary condition (corresponding to assuming dark background outside the scene) are Toeplitz matrices for one-dimensional problems and block-Toeplitz--Toeplitz-block matrices for two-dimensional cases. They are computationally intensive to invert especially in the block case. If the periodic boundary condition is used, the matrices become (block) circulant and can be diagonalized by discrete Fourier transform matrices. In this paper, we consider the use of the Neumann boundary condition (corresponding to a reflection of the original scene at the boundary). The resulting matrices are (block) Toeplitz-plus-Hankel matrices. We show that for symmetric blurring functions, these blurring matrices can always be diagonalized by discrete cosine transform matrices. Thus the cost of inversion is significantly lower than that of using the zero or periodic boundary conditions. We also show that the use of the Neumann boundary condition provides an easy way of estimating the regularization parameter when the generalized cross-validation is used. When the blurring function is nonsymmetric, we show that the optimal cosine transform preconditioner of the blurring matrix is equal to the blurring matrix generated by the symmetric part of the blurring function. Numerical results are given to illustrate the efficiency of using the Neumann boundary condition.
Introduction
A fundamental issue in signal and image processing is blur removal. The signal or image obtained
from a point source under the blurring process is called the impulse response function or the
point spread function. The observed signal or image g is just the convolution of this blurring
function h with the "true" signal or image f . The deblurring problem is to recover f from the
blurred function g given the blurring function h. This basic problem appears in many forms in
signal and image processing [2, 5, 12, 14].
In practice, the observed signal or image g is of finite length (and width) and we use it to
recover a finite section of f . Because of the convolution, g is not completely determined by
f in the same domain where g is defined. More precisely, if a blurred signal g is defined on
the interval [a; b] say, then it is not completely determined by the values of the true signal f
on [a; b] only. It is also affected by the values of f close to the boundary of [a; b] because of
the convolution. The size of the interval that affects g depends on the support of the blurring
function h. Thus in solving f from a finite length g, we need some assumptions on the values of
f outside the domain where g is defined. These assumptions are called the boundary conditions.
The natural and classical approach is to use the zero (Dirichlet) boundary condition [2,
pp.211-220]. It assumes that the values of f outside the domain of consideration are zero. This
results in a blurring matrix which is a Toeplitz matrix in the 1-dimensional case and a block-
Toeplitz-Toeplitz-block matrix in the 2-dimensional case, see [2, p.71]. However, these matrices
are known to be computationally intensive to invert, especially in the 2-dimensional case, see [2,
p.126]. Also ringing effects will appear at the boundary if the data are indeed not close to zero
outside the domain.
One way to alleviate the computational cost is to assume the periodic boundary condition,
i.e., data outside the domain of consideration are exact copies of data inside [12, p.258]. The
resulting blurring matrix is a circulant matrix in the 1-dimensional case and a block-circulant-
circulant-block matrix in the 2-dimensional case. These matrices can be diagonalized by discrete
Fourier matrices and hence their inverses can easily be found by using the Fast Fourier Transforms
see [12, p.258]. However, ringing effects will also appear at the boundary unless
f is close to periodic, and that is not common in practice.
In the image processing literature, other methods have also been proposed to assign boundary
values, see Lagendijk and Biemond [18, p.22] and the references therein. For instance, the
boundary values may be fixed at a local image mean, or they can be obtained by a model-based
extrapolation. In this paper, we consider the use of the Neumann (reflective) boundary condition
for image restoration. It sets the data outside the domain of consideration as reflection of the
data inside. The Neumann boundary condition has been studied in image restoration [21, 3, 18]
and in image compression [25, 20]. In image restoration, the boundary condition restores a
balance that is lost by ignoring the energy that spreads outside of the area of interest [21],
and also minimizes the distortion at the borders caused by deconvolution algorithms [3]. This
approach can also eliminate the artificial boundary discontinuities contributed to the energy
compaction property that is exploited in transform image coding [20].
The use of the Neumann boundary condition results in a blurring matrix that is a Toeplitz-
plus-Hankel matrix in the 1-dimensional case and a block Toeplitz-plus-Hankel matrix with
Toeplitz-plus-Hankel blocks in the 2-dimensional case. Although these matrices have more
complicated structures, we show that they can always be diagonalized by the discrete cosine
transform matrix provided that the blurring function h is symmetric. Thus their inverses can be
obtained by using fast cosine transforms (FCTs). Because FCT requires only real multiplications
and can be done at half of the cost of FFT, see [23, pp.59-60], inversion of these matrices is
faster than that of those matrices obtained from either the zero or periodic boundary conditions.
We also show that the use of the Neumann boundary condition provides an easy way of
estimating the regularization parameter when using the generalized cross-validation. We remark
that blurring functions are usually symmetric, see [14, p.269]. However, in the case where the
blurring function is nonsymmetric, we show that the optimal cosine transform preconditioner
[6] of the blurring matrix is generated by the symmetric part of the blurring function. Thus if
the blurring function is close to symmetric, the optimal cosine transform preconditioner should
be a good preconditioner.
The outline of the paper is as follows. In x2, we introduce the three different boundary con-
ditions. In x3, we show that symmetric blurring matrices obtained from the Neumann boundary
condition can always be diagonalized by the discrete cosine transform matrix. In x4, we show
that using the Neumann boundary condition, the generalized cross-validation estimate of the
regularization parameter can be done in a straightforward way. In x5, we give the construction
of the optimal cosine transform preconditioners for the matrices generated by nonsymmetric
blurring functions. In x6, we illustrate by numerical examples from image restorations that our
algorithm is efficient. Concluding remarks are given in x7.
2 The Deblurring Problem
For simplicity, we begin with the 1-dimensional deblurring problem. Consider the original signal
and the blurring function given by
The blurred signal is the convolution of h and ~ f , i.e., the i-th entry g i of the blurred signal is
given by
The deblurring problem is to recover the vector given the blurring function
h and a blurred signal of finite length. From (2), we
f \Gammam+1
f \Gammam+2
Thus the blurred signal g is determined not by f only, but by (f
. The linear system (3) is underdetermined. To overcome this, we make certain assumptions
(called boundary conditions) on the unknown data f \Gammam+1
so as to reduce the number of unknowns.
Before we discuss the boundary conditions, let us first rewrite (3) as
where
f \Gammam+1
f \Gammam+2
hm
. h \Gammam
2.1 The Zero (Dirichlet) Boundary Condition
The zero (or Dirichlet) boundary condition assumes that the signal outside the domain of the
observed vector g is zero [2, pp.211-220], i.e.,
the zero vector. The matrix system in (4) becomes
We see from (6) that the coefficient matrix T is a Toeplitz matrix.
There are many iterative or direct Toeplitz solvers that can solve the Toeplitz system (8)
with costs ranging from O(n log n) to O(n 2 ) operations, see for instance [19, 16, 1, 7]. In the 2-
dimensional case, the resulting matrices will be block-Toeplitz-Toeplitz-block matrices. Inversion
of these matrices is known to be very expensive, e.g. the fastest direct Toeplitz solver is of O(n 4 )
operations for an n 2 -by-n 2 block-Toeplitz-Toeplitz-block matrix, see [17].
2.2 The Periodic Boundary Condition
For practical applications, especially in the 2-dimensional case, where we need to solve the system
efficiently, one usually resort to the periodic boundary condition. This amounts to setting
in (3), see [12, p.258]. The matrix system in (4) becomes
are n-by-n Toeplitz matrices obtained by augmenting (n \Gamma m) zero
columns to T l and T r respectively.
The most important advantage of using the periodic boundary condition is that B so obtained
is a circulant matrix. Hence B can be diagonalized by the discrete Fourier matrix and (9) can
be solved by using three fast Fourier transforms (FFTs) (one for finding the eigenvalues of the
matrix B and two for solving the system, cf (15) below). Thus the total cost is of O(n log n)
operations.
In the 2-dimensional case, the blurring matrix is a block-circulant-circulant-block matrix
and can be diagonalized by the 2-dimensional FFTs (which are tensor-products of 1-dimensional
FFTs) in O(n 2 log n) operations.
2.3 The Neumann Boundary Condition
For the Neumann boundary condition, we assume that the data outside f are reflection of the
data inside f . More precisely, we set
and
in (3). Thus (4) becomes
where J is the n-by-n reversal matrix.
We remark that the coefficient matrix A in (10) is neither Toeplitz nor circulant. It is a
Toeplitz-plus-Hankel matrix. Although these matrices have more complicated structures, we
will show in x3 that the matrix A can always be diagonalized by the discrete cosine transform
matrix provided that the blurring function h is symmetric, i.e., h
It follows that (10) can be solved by using three fast cosine transforms (FCTs) in O(n log n)
operations, see (15) below. This approach is computationally attractive as FCT requires only
real operations and is about twice as fast as the FFT, see [23, pp.59-60]. Thus solving a problem
with the Neumann boundary condition is twice as fast as solving a problem with the periodic
boundary condition.
We will establish similar results in the 2-dimensional case, where the blurring matrices will
be block Toeplitz-plus-Hankel matrices with Toeplitz-plus-Hankel blocks.
3 Diagonalization of the Neumann Blurring Matrices
3.1 One-Dimensional Problems
We first review some definitions and properties of the discrete cosine transform matrix. Let C
be the n-by-n discrete cosine transform matrix with entries
r
is the Kronecker delta, see [14, p.150]. We note that C is orthogonal, i.e., C t
Also, for any n-vector v, the matrix-vector multiplications Cv and C t v can be
computed in O(n log n) real operations by FCTs; see [23, pp.59-60].
Let C be the space containing all matrices that can be diagonalized by C, i.e.
is an n-by-n real diagonal matrixg: (12)
to both sides of we see that
the eigenvalues [ ] i;i of Q are given by
Hence, the eigenvalues of Q can be obtained by taking an FCT of the first column of Q. In
particular, any matrix in C is uniquely determined by its first column.
Next we give a characterization of the class of matrices C. Let us define the shift of any vector
(v) to be the n-by-n symmetric
Toeplitz matrix with v as the first column and H(x; y) to be the n-by-n Hankel matrix with x
as the first column and y as the last column.
Lemma 1 (Chan, Chan, and Wong [6], Kailath and Olshevsky [15], Martucci [22]
and Sanchez et al. [24]) Let C be the class of matrices that can be diagonalized by the discrete
cosine transform matrix C. Then
It follows from Lemma 1 that matrices that can be diagonalized by C are some special
Toeplitz-plus-Hankel matrices.
Theorem 1 Let the blurring function h be symmetric, i.e., h Then the matrix
A given in (10) can be written as
where In particular, A can be diagonalized by C.
Proof: By (10), . By (6), it is clear that T is equal to T (u). From
the definitions of T l and T r in (5) and (7), it is also obvious that
Hence
By Theorem 1, the solution f of (10) is given by
where   is the diagonal matrix holding the eigenvalues of A. By (13),   can be obtained in one
FCT. Hence f can be obtained in three FCTs.
We remark that from (14), it is straightforward to construct the Neumann blurring matrix
A from the Dirichlet blurring matrix (6). All we need is to reflect the first column
of T to get the Hankel matrix H(oe(u); J oe(u)) and add it to T . Clearly the storage requirements
of both matrices A and T are the same - we only need to store the first column.
3.2 Two-Dimensional Problems
The results of x3.1 can be extended in a natural way to 2-dimensional image restoration problems.
In this case, one is concerned with solving a least squares problem similar to that in (3), except
that the matrix is now a block matrix. For the zero boundary condition, the resulting blurring
matrix is a block-Toeplitz-Toeplitz-block matrix of the form
where each block T (j) is a Toeplitz matrix of the form given in (7). The first column and row
of T in (16) are completely determined by the blurring function of the blurring process.
With the Neumann boundary condition, the resulting matrix A is a block Toeplitz-plus-
Hankel matrix with Toeplitz-plus-Hankel blocks. More precisely,
A (m) . A (\Gammam)
A (m\Gamma1) A (\Gammam+1)
with each block A (j) being an n-by-n matrix of the form given in (10). We note that the A (j) in
(17) and the T (j) in (16) are related by (14). Thus again it is straightforward to construct the
blurring matrix A from the matrix T or from the blurring function directly. Obviously, storage
requirements of A and T are the same.
We next show that for a symmetric blurring function, the blurring matrix A in (17) can be
diagonalized by the 2-dimensional discrete cosine transform matrix. Hence inversion of A can
be done by using only three 2-dimensional FCTs.
Theorem 2 If the blurring function h i;j is symmetric, i.e.,
for all i and j, then A can be diagonalized by the 2-dimensional discrete cosine transform matrix
C\Omega C,
where\Omega is the tensor product.
Proof: We note that
(C\Omega I)(I\Omega C). Since each block A (j) in (17) is of the form given
by (14), by Theorem 1, A (j) can be diagonalized by C, i.e.,
It follows that
(C\Omega C)A(C
(C\Omega I)(I\Omega C)A(I\Omega C t )(C
t\Omega I)
(C\Omega I ) (C
t\Omega I)
where
(m) .   (m)
Let P be the permutation matrix that satisfies
i.e. the (i; j)th entry of the (k; ')th block in   is permuted to the (k; ')th entry of the (i; j)th
block. Then we have P t
(I\Omega C) and
~
A (1) 0
~
A (2)
From (18), each matrix ~
A (j) has the same form as A in (14). In particular, for all j, C ~
~
(j) , a diagonal matrix. Thus
(C\Omega C)A(C
(C\Omega I ) (C
t\Omega I)
~
~
(2)
which is a permutation of a diagonal matrix and hence is still a diagonal matrix.
4 Estimation of Regularization Parameters
Besides the issue of boundary conditions, it is well-known that blurring matrices are in general
ill-conditioned and deblurring algorithms will be extremely sensitive to noise [12, p.282]. The
ill-conditioning of the blurring matrices stems from the wide range of magnitudes of their eigen-values
[10, p.31]. Therefore, excess amplification of the noise at small eigenvalues can occur. The
method of regularization is used to achieve stability for deblurring problems. In the classical
regularization [10, p.117], stability is attained by introducing a regularization operator
which restricts the set of admissible solutions. More specifically, the regularized solution f (-)
is computed as the solution to
min
The term kDf(-)k 2
2 is added in order to regularize the solution. The regularization parameter
- controls the degree of regularity (i.e., the degree of bias) of the solution.
One can find the solution f (-) in (19) by solving the normal equations
Usually, kDfk 2 is chosen to be the L 2 norm kfk 2 or the H 1 norm kLfk 2 where L is the first
order difference operator matrix, see [14, 8, 12]. Correspondingly, the matrix D t D in (20) is
the identity matrix or the discrete Laplacian matrix with some boundary conditions. In the
latter case, if the zero boundary condition is imposed, D t D is just the discrete Laplacian with
the Dirichlet boundary condition. For the periodic boundary condition, D t D is circulant and
can be diagonalized by the FFTs, see for instance [12, p.283]. For the Neumann boundary
condition, D t D is the discrete Laplacian with the Neumann boundary condition, which can be
diagonalized by the discrete cosine transform matrix, see for instance [4]. Thus if we use the
Neumann boundary condition for both the blurring matrix A and the regularization operator
D, then the matrix in (20) can be diagonalized by the discrete cosine transform matrix and
hence its inversion can still be done in three FCTs for any fixed -, cf. (15).
Another difficulty in regularization is the choice of -. Generalized cross-validation [11] is
a technique that estimates - directly without requiring an estimate of the noise variance. It
is based on the concept of prediction errors. For each
- be the vector that
minimizes the error measure:
is the ith element of Af and [g] i is the ith element of g. If - is such that f k
- is a good
estimate of f , then [Af k
should be a good approximation of [g] k on average. For a given -, the
average squared error between the predicted value [Af k
and the actual value [g] k is given byn
The generalized cross-validation (GCV) is a weighted version of the above error:
where m jj (-) is the (j; j)th entry of the so-called influence matrix
In [11], Golub et al. have shown that v(-) can be written as
The optimal regularization parameter is chosen to be the - that minimizes v(-). Since v(-)
is a nonlinear function, the minimizer usually cannot be determined analytically. However, if
the Neumann boundary condition is used for both A and D t D, we can rewrite v(-) as
where ff i and fi i represent the eigenvalues of A and D t D respectively. We recall that ff i and
can be obtained by taking the FCT of the first column of A and D t D respectively since the
matrices can be diagonalized by the discrete cosine transform matrix C. Thus the GCV estimate
of - can be computed in a straightforward manner, see [13].
For the periodic boundary condition, the GCV estimate can also be computed by a similar
procedure. However, if we use the zero boundary condition, determining the GCV estimate of -
will require the inversion of a large matrix which is clearly an overwhelming task for any images
of reasonable size.
5 Optimal Cosine Transform Preconditioners
Because all matrices in C are symmetric (see (12)), discrete cosine transform matrices can only
diagonalized blurring matrices from symmetric blurring functions. For nonsymmetric blurring
functions, matrices in C may be used as preconditioners to speed up the convergence of iterative
methods such as the conjugate gradient method. Given a matrix A, we define the optimal cosine
transform preconditioners c(A) to be the minimizer of kQ \Gamma
is the Frobenius norm.
In [6, 16], c(A) are obtained by solving linear systems. Here we give a simple approach for
finding c(A).
Theorem 3 Let h be an arbitrary blurring function and A be the blurring matrix of h with the
Neumann boundary condition imposed. Then the optimal cosine transform preconditioner c(A)
of A can be found as follows:
1. In the one-dimensional case, c(A) is the blurring matrix corresponding to the symmetric
blurring function s with the Neumann boundary condition imposed.
2. In the 2-dimensional case, c(A) is the blurring matrix corresponding to the symmetric
blurring function given by s i;j j (h i;j +h i;\Gammaj +h \Gammai;j +h \Gammai;\Gammaj )=4 with the Neumann boundary
condition imposed.
Proof: We only give the proof for the one-dimensional case. The proof for the two-dimensional
case is similar. We first note that if U and V are symmetric and skew-symmetric matrices
respectively, then kU
F . Hence, for any
Since the second term in the right hand side above does not affect the diagonal matrix  , the
minimizer   is given by the diagonal of
It is easy to check that2
where
and
We claim that the diagonal entries [CH(v; \GammaJ v)]C t this is true, then the
minimizer   is given by
and Theorem 3 then follows directly from Theorem 1.
To prove our claim, we first note that by (11),
r
cos
Also it is clear that T j JH(v; \GammaJ v) is a skew-symmetric Toeplitz matrix. Therefore
In view of Theorem 3 and the results we have in x3, it is easy to find c(A) for blurring
matrices generated by nonsymmetric blurring functions. We just take the symmetric part of the
blurring functions and form the (block) Toeplitz-plus-Hankel matrices as given in (10) or (17).
From Theorem 3, we also see that if the blurring function is close to symmetric, then c(A) will
be a good approximation (hence a good preconditioner) to A, see the numerical results in x6.
6 Numerical Experiments
In this section, we illustrate the efficiency of employing the Neumann boundary condition over
the other two boundary conditions for image restoration problems. All our tests were done by
Matlab. The data source is a photo from the 1964 Gatlinburg Conference on Numerical Algebra
taken from Matlab. From (4), we see that to construct the right hand side vector g correctly, we
need the vectors f l and f r , i.e., we need to know the image outside the given domain. Thus we

Figure

1: The "Gatlinburg Conference" image.1020150.0050.0150.0250.0351020515250.010.020.03
(a) (b) (c) (d)

Figure

2: (a) Gaussian blur; (b) out-of-focus blur; (c) noisy and blurred image by Gaussian blur
and (d) by out-of-focus blur.
start with the 480-by-640 image of the photo and cut out a 256-by-256 portion from the image.

Figure

1 gives the 256-by-256 image of this picture.
We consider restoring the "Gatlinburg Conference" image blurred by the following two blurring
functions, see [14, p.269]:
(i) a truncated Gaussian blur:
ae
ce
(ii) an out-of-focus blur:
ae c; if
where h i;j is the jth entry of the first column of T (i) in (16) and c is the normalization constant
such that
We remark that the Gaussian blur is symmetric and separable whereas
the out-of-focus blur is symmetric but not separable, see Figures 2(a) and 2(b). A Gaussian
white noise n with signal-to-noise ratio of 50dB is then added to the blurred images. The noisy
blurred images are shown in Figures 2(c) and 2(d). We note that after the blurring, the cigarette
held by Prof. Householder (the rightmost person) is not clearly shown, cf. Figure 1.
We remark that the regularization parameter based on the generalized cross-validation
method (see x4) is not suitable when the zero boundary condition is imposed. The method
will require the inversion of many large matrices. As a comparison of the cost among different
boundary conditions, we chose the optimal regularization parameter -   such that it minimizes
the relative error of the reconstructed image f (-) which is defined as kf \Gamma f (-)k 2 =kfk 2 where f
is the original image. The optimal -   , accurate up to one significant digit, is obtained by trial
and error.
In

Figures

3 and 4, we present the restored images for the three different boundary conditions,
the optimal -   and the relative errors. We used the L 2 norm as the regularization functional
here. We see from the figures that by imposing the Neumann boundary condition, the relative
error and the ringing effect are the smallest. Also the cigarette is better reconstructed by using
the Neumann boundary condition than by those using the other two boundary conditions.
rel.

Figure

3: Restoring Gaussian blur with zero (left), periodic (middle) and Neumann (right)
boundary conditions.
rel.

Figure

4: Restoring out-of-focus blur with zero (left), periodic (middle) and Neumann (right)
boundary conditions.
Next let us consider the cost. Recall that for each -, we only need three 2-dimensional FFTs
and FCTs to compute the restored images for the periodic and the Neumann boundary conditions
respectively. Thus the costs for both approaches are about O(n 2 log n) operations though the
Neumann one is twice as fast because FCT requires only real multiplications [23, pp.59-60]. For
the zero boundary condition, we have to solve a block-Toeplitz-Toeplitz-block system for each -.
The fastest direct Toeplitz solver requires O(n 4 ) operations, see [17]. In our tests, the systems
Blurring function 3 \Theta
Gaussian 67
Out-of-focus 81 73 48 43 26
Blurring function 3 \Theta
Gaussian 17 17 13 12 7 7
Out-of-focus 22 21 15 14 9 9

Table

1: The numbers of iterations required for using the zero boundary condition.
are solved by the preconditioned conjugate gradient method with circulant preconditioners [7].

Table

1 shows the numbers of iterations required for the two blurring functions for different -.
The stopping tolerance is 10 \Gamma6 . We note that the cost per iteration is about four 2-dimensional
FFTs. Thus the cost is extremely expensive especially when - is small. In conclusion, we see
that the cost of using the Neumann boundary condition is lower than that of using the other
two boundary conditions.
Finally we illustrate the effectiveness of the optimal cosine transform preconditioners for
blurring functions that are close to symmetric. More general tests on the preconditioners are
given in [6]. We consider a 2-dimensional deconvolution problem arising in the ground-based
atmospheric imaging. Figure 5(a) gives the 256-by-256 blurred and noisy image of an ocean
reconnaissance satellite observed by a ground-based imaging system and Figure 5(b) is a 256-
by-256 image of a guide star observed under similar circumstances, see [8]. The discrete blurring
function h is given by the pixel values of the guide star image. The blurring matrix A is obtained
as in (17) by imposing the Neumann boundary condition.
50 100 150 200 2500.40.8

Figure

5: (a) Observed image, (b) the guide star image and (c) a cross-section of the blurring
function.
We note that the blurring function is not exactly symmetric in this case, see Figure 5(b).
However, from the cross-sections of the blurring function (see for instance Figure 6(c)), we know
that it is close to symmetric. Therefore, we used the preconditioned conjugate gradient algorithm
with the optimal cosine transform preconditioner to remove the blurring, see x5. Again we use
the L 2 norm as the regularization functional here. In Figure 6(a), we present the restored image
with the optimal -   . The original image is given in Figure 6(b) for comparison.

Figure

(a) Restored image (-  (b) the true image.

Table

2: The number of iterations for convergence.
Using the optimal cosine transform preconditioner, the image is restored in 4 iterations for
a stopping tolerance of 10 \Gamma6 . If no preconditioner is used, acceptable restoration is achieved
after 134 iterations, see Table 2. We remark that the cost per iteration for using the optimal
cosine transform preconditioner is almost the same as that with no preconditioner: they are
about 1:178 \Theta 10 8 and 1:150 \Theta 10 8 floating point operations per iteration respectively. Thus
we see that the preconditioned conjugate gradient algorithm with the optimal cosine transform
preconditioner is an efficient and effective method for this problem.
7 Concluding Remarks
In this paper, we have shown that discrete cosine transform matrices can diagonalize dense
Toeplitz-plus-Hankel blurring matrices arising from using the Neumann (reflective)
boundary condition. Numerical results suggest that the Neumann boundary condition provides
an effective model for image restoration problems, both in terms of the computational cost
and of minimizing the ringing effects near the boundary. It is interesting to note that discrete
sine transform matrices can diagonalize Toeplitz matrices with at most 3 bands (such as the
discrete Laplacian with zero boundary conditions) but not dense Toeplitz matrices in general,
see [9] for instance.



--R

Superfast Solution of Real Positive Definite Toeplitz Systems

IEEE Signal Processing Mag- azine
Fast Computation of a Discretized Thin-plate Smoothing Spline for Image Data

Based Preconditioners for Total Variation Deblurring
Conjugate Gradient Methods for Toeplitz Systems
Generalization of Strang's Preconditioner with Applications to Toeplitz Least Squares Problems
Sine Transform Based Preconditioners for Symmetric Toeplitz Systems
Regularization of Inverse Problems
Generalized Cross-validation as a Method for Choosing a Good Ridge Parameter
New York
a Matlab Package for Analysis and Solution of Discrete Ill-Posed Problems
Fundamentals of

Theory and Applications
Fast Algorithms for Block Toeplitz Matrices with Toeplitz Entries
Iterative Identification and Restoration of Images
The Wiener RMS (Root Mean Square) Error Criterion in Filter Design and Prediction

Reducing Boundary Distortion in Image Restoration
Symmetric Convolution and the Discrete Sine and Cosine Transforms

Diagonalization Properties of the Discrete Cosine Transforms
The
--TR

--CTR
Ben Appleton , Hugues Talbot, Recursive filtering of images with symmetric extension, Signal Processing, v.85 n.8, p.1546-1556, August 2005
Marco Donatelli , Claudio Estatico , Stefano Serra-Capizzano, Boundary conditions and multiple-image re-blurring: the LBT case, Journal of Computational and Applied Mathematics, v.198 n.2, p.426-442, 15 January 2007
Michael K. Ng , Andy M. Yip, A Fast MAP Algorithm for High-Resolution Image Reconstruction with Multisensors, Multidimensional Systems and Signal Processing, v.12 n.2, p.143-164, April 2001
Daniela Calvetti, Preconditioned iterative methods for linear discrete ill-posed problems from a Bayesian inversion perspective, Journal of Computational and Applied Mathematics, v.198 n.2, p.378-395, 15 January 2007
Jian-Feng Cai , Raymond H. Chan , Carmine Fiore, Minimization of a Detail-Preserving Regularization Functional for Impulse Noise Removal, Journal of Mathematical Imaging and Vision, v.29 n.1, p.79-91, September 2007
Michael K. Ng , Andy C. Yau, Super-Resolution Image Restoration from Blurred Low-Resolution Images, Journal of Mathematical Imaging and Vision, v.23 n.3, p.367-378, November  2005
toeplitz matrix;circulant matrix;boundary conditions;cosine transform;deblurring;hankel matrix