Introduction In electrical engineering and computer science, image processing refers to any form of signal processing in which the input is an image and the output can be either an image or a set of parameters related to the image. Generally, image processing includes image enhancement, restoration and reconstruction, edge and boundary detection, classification and segmentation, object recognition and identification, compression and communication, etc. Among them, image restoration is a classical problem and is generally a preprocessing stage of higher level processing. In many applications, the measured images are degraded by blurs; e.g. the optical system in a camera lens may be out of focus, so that the incoming light is smeared out, and in astronomical imaging the incoming light in the telescope has been slightly bent by turbulence in the atmosphere. In addition, images that occur in practical applications inevitably suffer from noise, which arise from numerous sources such as radiation scatter from the surface before the image is sensed, electrical noise in the sensor or camera, transmission errors, and bit errors as the image is digitized, etc. In such situations, the image formation process is usually modeled by the following equation f ( x ) = ( k * u ¯ ) ( x ) + ω ( x ) , x ∈ Ω , where u¯(x) is an unknown clean image over a region Ω⊂R2, “*" denotes the convolution operation, k(x),n(x) and f(x) are real-valued functions from R2 to R representing, respectively, convolution kernel, additive noise, and the blurry and noisy observation. Usually, the convolution process neither absorbs nor generates optical energy, i.e., ∫Ωk(x)dx=1, and the additive noise has zero mean. Deblurring or decovolution aims to recover the unknown image u¯(x) from f(x) and k(x) based on (). When k(x) is unknown or only an estimate of it is available, recovering u¯(x) from f(x) is called blind deconvolution. Throughout this module, we assume that k(x) is known and ω(x) is either Gaussian or impulsive noise. When k(x) is equal to the Dirac delta, the recovery of u¯(x) becomes a pure denoising problem. In the rest of this section, we review the TV-based variational models for image restoration and introduce necessary notation for analysis.

Total Variation for Image Restoration The TV regularization was first proposed by Rudin, Osher and Fatemi in for image denoising, and then extended to image deblurring in . The TV of u is defined as TV ( u ) = ∫ Ω ∥ ∇ u ( x ) ∥ d x . When ∇u(x) does not exist, the TV is defined using a dual formulation , which is equivalent to () when u is differentiable. We point out that, in practical computation, discrete forms of regularization are always used where differential operators are replaced by ceratin finite difference operators. We refer TV regularization and its variants as TV-like regularization. In comparison to Tikhonov-like regularization, the homogeneous penalty on image smoothness in TV-like regularization can better preserve sharp edges and object boundaries that are usually the most important features to recover. Variational models with TV regularization and ℓ2 fidelity has been widely studied in image restoration; see e.g. , and references therein. For ℓ1 fidelity with TV regularization, its geometric properties are analyzed in , , . The superiority of TV over Tikhonov-like regularization was analyzed in , for recovering images containing piecewise smooth objects. Besides Tikhonov and TV-like regularization, there are other well studied regularizers in the literature, e.g. the Mumford-Shah regularization . In this module, we concentrate on TV-like regularization. We derive fast algorithms, study their convergence, and examine their performance.

Discretization and Notation As used before, we let ∥·∥ be the 2-norm. In practice, we always discretize an image defined on Ω, and vectorize the two-dimensional digitalized image into a long one-dimensional vector. We assume that Ω is a square region in R2. Specifically, we first discretize u(x) into a digital image represented by a matrix U∈Rn×n. Then we vectorize U column by column into a vector u∈Rn2, i.e. u i = U p q , i = 1 , ... , n 2 , where ui denotes the ith component of u, Upq is the component of U at pth row and qth column, and p and q are determined by i=(q-1)n+p and 1≤q≤n. Other quantities such as the convolution kernel k(x), additive noise ω(x), and the observation f(x) are all discretized correspondingly. Now we present the discrete forms of the previously presented equations. The discrete form of () is f = K u ¯ + ω , where in this case, u¯,ω,f∈Rn2 are all vectors representing, respectively, the discrete forms of the original image, additive noise and the blurry and noisy observation, and K∈Rn2×n2 is a convolution matrix representing the kernel k(x). The gradient ∇u(x) is replaced by certain first-order finite difference at pixel i. Let Di∈R2×n2 be a first-order local finite difference matrix at pixel i in horizontal and vertical directions. E.g. when the forward finite difference is used, we have D i u = u i + n - u i u i + 1 - u i ∈ R 2 , for i=1,...,n2 (with certain boundary conditions assumed for i>n2-n). Then the discrete form of TV defined in () is given by TV ( u ) = ∑ i = 1 n 2 ∥ D i u ∥ . We will refer to min u T V ( u ) + μ 2 ∥ K u - f ∥ 2 with discretized TV regularization () as TV/L2. For impulsive noise, we replace the ℓ2 fidelity by ℓ1 fidelity and refer to the resulted problem as TV/L1. Now we introduce several more notation. For simplicity, we let ∑i be the summation taken over all pixels. The two first-order global finite difference operators in horizontal and vertical directions are, respectively, denoted by D(1) and D(2) which are n2-by-n2 matrices (boundary conditions are the same as those assumed on Di). As such, it is worth noting that the two-row matrix Di is formed by stacking the ith row of D(1) on that of D(2). For vectors v1 and v2, we let v=(v1;v2)≜(v1⊤,v2⊤)⊤, i.e. v is the vector formed by stacking v1 on the top of v2. Similarly, we let D=(D(1);D(2))=(D(1))⊤,(D(2))⊤⊤. Given a matrix T, we let diag (T) be the vector containing the elements on the diagonal of T, and F(T)=FTF-1, where F∈n2×n2 is the 2D discrete Fourier transform matrix.

Existing Methods Since TV is nonsmooth, quite a few algorithms are based on smoothing the TV term and solving an approximation problem. The TV of u is usually replaced by TV ϵ ( u ) = ∑ i ∥ D i u ∥ 2 + ϵ , where ϵ>0 is a small constant. Then the resulted approximate TV/L2 problem is smooth and many optimization methods are available. Among others, the simplest method is the gradient descent method as was used in . However, this method suffers slow convergence especially when the iterate point is close to the solution. Another important method is the linearized gradient method proposed in for denoising and in for deblurring. Both the gradient descent and the linearized gradient methods are globally and at best linearly convergent. To obtain super linear convergence, a primal-dual based Newton method was proposed in . Both the linearized gradient method and this primal-dual method need to solve a large system of linear equations at each iteration. When ϵ is small and/or K becomes more ill-conditioned, the linear system becomes more and more difficult to solve. Another class of well-known methods for TV/L2 are the iterative shrinkage/thresholding (IST) based methods . For IST-based methods, a TV denoising problem needs to be solved at each iteration. Also, in the authors transformed the TV/L2 problem into a second order cone program and solved it by interior point method.

A New Alternating Minimization Algorithm In this section, we derive a new algorithm for the TV/L2 problem min u ∑ i ∥ D i u ∥ + μ 2 ∥ K u - f ∥ 2 . In (), the fidelity term is quadratic with respect to u. Moreover, K is a convolution matrix and thus can be easily diagonalized by fast transforms (with proper boundary conditions assumed on u). Therefore, the main difficulty in solving () is caused by the nondifferentiability and the universal coupling of variables of the TV term. Our algorithm is derived from the well-known variable-splitting and penalty techniques in optimization. First, we introduce an auxiliary variable wi∈R2 at pixel i to transfer Diu out of the nondifferentiable term ∥·∥. Then we penalize the difference between wi and Diu quadratically. As such, the auxiliary variables wi's are separable with respect to one another. For convenience, in the following we let w≜[w1,...,wn2]. The approximation model to () is given by min w , u ∑ i ∥ w i ∥ + β 2 ∑ i ∥ w i - D i u ∥ 2 + μ 2 ∥ K u - f ∥ 2 , where β≫0 is a penalty parameter. It is well known that the solution of () converges to that of () as β→∞. In the following, we concentrate on problem ().

Basic Algorithm The benefit of () is that while either one of the two variables u and w is fixed, minimizing the objective function with respect to the other has a closed-form formula that we will specify below. First, for a fixed u, the first two terms in () are separable with respect to wi, and thus the minimization for w is equivalent to solving min w i ∥ w i ∥ + β 2 ∥ w i - D i u ∥ 2 , i = 1 , 2 , ... , n 2 . It is easy to verify that the unique solutions of () are w i = max ∥ D i u ∥ - 1 β , 0 D i u ∥ D i u ∥ , i = 1 , ... , n 2 , where the convention 0·(0/0)=0 is followed. On the other hand, for a fixed w, () is quadratic in u and the minimizer u is given by the normal equations ∑ i D i ⊤ D i + μ β K ⊤ K u = ∑ i D i ⊤ w i + μ β K ⊤ f . By noting the relation between D and Di and a reordering of variables, () can be rewritten as D ⊤ D + μ β K ⊤ K u = D ⊤ w + μ β K ⊤ f , where w ≜ w 1 w 2 ∈ R 2 n 2 and w j ≜ ( w 1 ) j ⋮ ( w n 2 ) j , j = 1 , 2 . The normal equation () can also be solved easily provided that proper boundary conditions are assumed on u. Since both the finite difference operations and the convolution are not well-defined on the boundary of u, certain boundary assumptions are needed when solving (). Under the periodic boundary conditions for u, i.e. the 2D image u is treated as a periodic function in both horizontal and vertical directions, D(1), D(2) and K are all block circulant matrices with circulant blocks; see e.g. , . Therefore, the Hessian matrix on the left-hand side of () has a block circulant structure and thus can be diagonalized by the 2D discrete Fourier transform F, see e.g. . Using the convolution theorem of Fourier transforms, the solution of () is given by u = F - 1 F D ⊤ w + ( μ / β ) K ⊤ f diag F ( D ⊤ D + ( μ / β ) K ⊤ K ) , where the division is implemented by componentwise. Since all quantities but w are constant for given β, computing u from () involves merely the finite difference operation on w and two FFTs (including one inverse FFT), once the constant quantities are computed. Since minimizing the objective function in () with respect to either w or u is computationally inexpensive, we solve () for a fixed β by an alternating minimization scheme given below. Algorithm : Input f, K and μ>0. Given β>0 and initialize u=f. While “not converged”, Do Compute w according to () for fixed u. Compute u according to () for fixed w (or equivalently w). End Do

Optimality Conditions and Convergence Results To present the convergence results of Algorithm "Basic Algorithm" for a fixed β, we make the following weak assumption. Assumption 1 N(K)∩N(D)={0}, where N(·) represents the null space of a matrix. Define M = D ⊤ D + μ β K ⊤ K and T = D M - 1 D ⊤ . Furthermore, we will make use of the following two index sets: L = i : ∥ D i u * ∥ < 1 β and E = { 1 , ... , n 2 } ∖ L . Under Assumption 1, the proposed algorithm has the following convergence properties. Theorem 1 For any fixed β>0, the sequence {(wk,uk)} generated by Algorithm "Basic Algorithm" from any starting point (w0,u0) converges to a solution (w*,u*) of (). Furthermore, the sequence satisfies ∥ w E k + 1 - w E * ∥ ≤ ρ ( ( T 2 ) E E ) ∥ w E k - w E * ∥ ; ∥ u k + 1 - u * ∥ M ≤ ρ ( T E E ) ∥ u k - u * ∥ M ; for all k sufficiently large, where TEE=[Ti,j]i,j∈E∪(n2+E) is a minor of T, ∥v∥M2=v⊤Mv and ρ(·) is the spectral radius of its argument.

Extensions to Multichannel Images and TV/L The alternating minimization algorithm given in "A New Alternating Minimization Algorithm" can be extended to solve multichannel extension of () when the underlying image has more than one channels and TV/L1 when the additive noise is impulsive.

Multichannel image deconvolution Let u¯=[u¯(1);...;u¯(m)]∈Rmn2 be an m-channel image, where, for each j, u¯(j)∈Rn2 represents the jth channel. An observation of u¯ is modeled by (), in which case f=[f(1);...;f(m)] and ω=[ω(1);...;ω(m)] have the same size and the number of channels as u¯, and K is a multichannel blurring operator of the form K = K 11 K 12 ⋯ K 1 m K 21 K 22 ⋯ K 2 m ⋮ ⋮ ⋱ ⋮ K m 1 K m 2 ⋯ K m m ∈ R m n 2 × m n 2 , where Kij∈Rn2×n2, each diagonal submatrix Kii defines the blurring operator within the ith channel, and each off-diagonal matrix Kij, i≠j, defines how the jth channel affects the ith channel. The multichannel extension of () is min u ∑ i ∥ ( I m ⊗ D i ) u ∥ + μ 2 ∥ K u - f ∥ 2 , where Im is the identity matrix of order m, and “⊗" is the Kronecker product. By introducing auxiliary variables wi∈R2m, i=1,...,n2, () is approximated by min w , u ∑ i ∥ w i ∥ + β 2 ∑ i ∥ w i - ( I m ⊗ D i ) u ∥ 2 + μ 2 ∥ K u - f ∥ 2 . For fixed u, the minimizer function for w is given by () in which Diu should be replaced by (Im⊗Di)u. On the other hand, for fixed w, the minimization for u is a least squares problem which is equivalent to the normal equations I 3 ⊗ ( D ⊤ D ) + μ β K ⊤ K u = ( I 3 ⊗ D ) ⊤ w + μ β K ⊤ f , where w is a reordering of variables in a similar way as given in (). Under the periodic boundary condition, () can be block diagonalized by FFTs and then solved by a low complexity Gaussian elimination method.

Deconvolution with Impulsive Noise When the blurred image is corrupted by impulsive noise rather than Gaussian, we recover u¯ as the minimizer of a TV/L1 problem. For simplicity, we again assume u¯∈Rn2 is a single channel image and the extension to multichannel case can be similarly done as in "Multichannel image deconvolution". The TV/L1 problem is min u ∑ i ∥ D i u ∥ + μ ∥ K u - f ∥ 1 . Since the data-fidelity term is also not differentiable, in addition to w, we introduce z∈Rn2 and add a quadratic penalty term. The approximation problem to () is min w , z , u ∑ i ∥ w i ∥ + β 2 ∥ w i - D i u ∥ 2 + μ ∥ z ∥ 1 + γ 2 ∥ z - ( K u - f ) ∥ 2 , where β,γ≫0 are penalty parameters. For fixed u, the minimization for w is the same as before, while the minimizer function for z is given by the famous one-dimensional shrinkage: z = max | K u - f | - 1 γ , 0 · sgn ( K u - f ) . On the other hand, for fixed w and z, the minimization for u is a least squares problem which is equivalent to the normal equations D ⊤ D + μ γ β K ⊤ K u = D ⊤ w + μ γ β K ⊤ ( f + z ) . Similar to previous arguments, () can be easily solved by FFTs.

Experiments In this section, we present the practical implementation and numerical results of the proposed algorithms. We used two images, Man (grayscale) and Lena (color) in our experiments, see . The two images are widely used in the field of image processing because they contain nice mixture of detail, flat regions, shading area and texture.

Test images: Man (left, 1024×1024) and Lena (right, 512×512). We tested several kinds of blurring kernels including Gaussian, average and motion. The additive noise is Gaussian for TV/L2 problems and impulsive for TV/L1 problem. The quality of image is measured by the signal-to-noise ratio (SNR) defined by SNR ≜ 10 * log 10 ∥ u ¯ - E ( u ¯ ) ∥ 2 ∥ u ¯ - u ∥ 2 , where u¯ is the original image and E(u¯) is the mean intensity value of u¯. All blurring effects were generated using the MATLAB function “imfilter " with periodic boundary conditions, and noise was added using “imnoise ". All the experiments were finished under Windows Vista Premium and MATLAB v7.6 (R2008a) running on a Lenovo laptop with an Intel Core 2 Duo CPU at 2 GHz and 2 GB of memory.

Practical Implementation Generally, the quality of the restored image is expected to increase as β increases because the approximation problems become closer to the original ones. However, the alternating algorithms converge slowly when β is large, which is well-known for the class of penalty methods. An effective remedy is to gradually increase β from a small value to a pre-specified one. compares the different convergence behaviors of the proposed algorithm when with and without continuation, where we used Gaussian blur of size 11 and standard deviation 5 and added white Gaussian noise with mean zero and standard deviation 10-3.

Continuation vs. no continuation: u* is an “exact” solution corresponding to β=214. The horizontal axis represents the number of iterations, and the vertical axis is the relative error ek=∥uk-u*∥/∥u*∥. In this continuation framework, we compute a solution of an approximation problem which used a smaller beta, and use the solution to warm-start the next approximation problem corresponding to a bigger β. As can be seen from , with continuation on β the convergence is greatly sped up. In our experiments, we implemented the alternating minimization algorithms with continuation on β, which we call the resulting algorithm “Fast Total Variation de-convolution” or FTVd, which, for TV/L2, the framework is given below. [FTVd]: Input f, K and μ>0. Given βmax>β0>0. Initialize u=f, up=0, β=β0 and ϵ>0. While β≤βmax, Do Run Algorithm "Basic Algorithm" until an optimality condition is met. β←2*β. End Do

SNRs of images recovered from () for different β.

Results recovered from TV/L2. Image Man is blurred by a Gaussian kernel, while image Lena is blurred by a cross-channel kernel. Gaussian noise with zero mean and standard deviation 10-3 is added to both blurred images. The left images are the blurry and noisy observations, and the right ones are recovered by FTVd. Generally, it is difficult to determine how large β is sufficient to generate a solution that is close to be a solution of the original problems. In practice, we observed that the SNR values of recovered images from the approximation problems are stabilized once β reached a reasonably large value. To see this, we plot the SNR values of restored images corresponding to β=20,21,⋯,218 in . In this experiment, we used the same blur and noise as we used in the testing of continuation. As can be seen from , the SNR values on both images essentially remain constant for β≥27. This suggests that β need not to be excessively large from a practical point of view. In our experiments, we set β0=1 and βmax=27 in Algorithm "Practical Implementation". For each β, the inner iteration was stopped once an optimality condition is satisfied. For TV/L1 problems, we also implement continuation on γ, and used similar settings as used in TV/L2.

Recovered Results In this subsection, we present results recovered from TV/L2 and TV/L1 problems including (), () and their multichannel extensions. We tested various of blurs with different levels of Gaussian noise and impulsive noise. Here we merely present serval test results. gives two examples of blurry and noisy images and the recovered ones, where the blurred images are corrupted by Gaussian noise, while gives the recovered results where the blurred images are corrupted by random-valued noise. For TV/L1 problems, we set γ=215 and β=210 in the approximation model and implemented continuation on both β and γ.

Results recovered from TV/L1. Image Lena is blurred by a cross-channel kernel and corrupted by 40% (left) and 50% (right) random-valued noise. The top row contains the blurry and noisy observations and the bottom row shows the results recovered by FTVd.

Concluding Remarks We proposed, analyzed and tested an alternating algorithm FTVd which for solving the TV/L2 problem. This algorithm was extended to solve the TV/L1 model and their multichannel extensions by incorporating an extension of TV. Cross-channel blurs are permitted when the underlying image has more than one channels. We established strong convergence results for the algorithms and validated a continuation scheme. Numerical results are given to demonstrate the feasibility and efficiency of the proposed algorithms.

Acknowledgements This Connexions module describes work conducted as part of Rice University's VIGRE program, supported by National Science Foundation grant DMS-0739420.