Total Variation for Image Restoration

The TV regularization was first proposed by Rudin, Osher and Fatemi in [12] for image denoising, and then extended to image deblurring in [11]. The TV of uu is defined as

When ∇u(x)∇u(x) does not exist, the TV is defined using a dual formulation [18], which is equivalent to (Equation 2) when uu 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ℓ2 fidelity has been widely studied in image restoration; see e.g. [4], [3] and references therein. For ℓ1ℓ1 fidelity with TV regularization, its geometric properties are analyzed in [2], [16], [17]. The superiority of TV over Tikhonov-like regularization was analyzed in [1], [5] 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 [9]. 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 R2R2. Specifically, we first discretize u(x)u(x) into a digital image represented by a matrix U∈Rn×nU∈Rn×n. Then we vectorize UU column by column into a vector u∈Rn2u∈Rn2, i.e.

where uiui denotes the iith component of uu, UpqUpq is the component of UU at ppth row and qqth column, and pp and qq are determined by i=(q-1)n+pi=(q-1)n+p and 1≤q≤n1≤q≤n. Other quantities such as the convolution kernel k(x)k(x), additive noise ω(x)ω(x), and the observation f(x)f(x) are all discretized correspondingly. Now we present the discrete forms of the previously presented equations. The discrete form of (Equation 1) is

where in this case, u¯,ω,f∈Rn2u¯,ω,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×n2K∈Rn2×n2 is a convolution matrix representing the kernel k(x)k(x). The gradient ∇u(x)∇u(x) is replaced by certain first-order finite difference at pixel ii. Let Di∈R2×n2Di∈R2×n2 be a first-order local finite difference matrix at pixel ii in horizontal and vertical directions. E.g. when the forward finite difference is used, we have

for i=1,...,n2i=1,...,n2 (with certain boundary conditions assumed for i>n2-ni>n2-n). Then the discrete form of TV defined in (Equation 2) is given by

We will refer to

with discretized TV regularization (Equation 6) as TV/L22. For impulsive noise, we replace the ℓ2ℓ2 fidelity by ℓ1ℓ1 fidelity and refer to the resulted problem as TV/L11.

Now we introduce several more notation. For simplicity, we let ∑i∑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)D(1) and D(2)D(2) which are n2n2-by-n2n2 matrices (boundary conditions are the same as those assumed on DiDi). As such, it is worth noting that the two-row matrix DiDi is formed by stacking the iith row of D(1)D(1) on that of D(2)D(2). For vectors v1v1 and v2v2, we let v=(v1;v2)≜(v1⊤,v2⊤)⊤v=(v1;v2)≜(v1⊤,v2⊤)⊤, i.e. vv is the vector formed by stacking v1v1 on the top of v2v2. Similarly, we let D=(D(1);D(2))=(D(1))⊤,(D(2))⊤⊤D=(D(1);D(2))=(D(1))⊤,(D(2))⊤⊤. Given a matrix TT, we let diag (T) diag (T) be the vector containing the elements on the diagonal of TT, and F(T)=FTF-1F(T)=FTF-1, where F∈n2×n2F∈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 uu is usually replaced by

where ϵ>0ϵ>0 is a small constant. Then the resulted approximate TV/L22 problem is smooth and many optimization methods are available. Among others, the simplest method is the gradient descent method as was used in [12]. 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 [14] for denoising and in [15] 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 [13]. 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 KK becomes more ill-conditioned, the linear system becomes more and more difficult to solve. Another class of well-known methods for TV/L22 are the iterative shrinkage/thresholding (IST) based methods [8]. For IST-based methods, a TV denoising problem needs to be solved at each iteration. Also, in [7] the authors transformed the TV/L22 problem into a second order cone program and solved it by interior point method.

Basic Algorithm

The benefit of (Equation 10) is that while either one of the two variables uu and ww 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 uu, the first two terms in (Equation 10) are separable with respect to wiwi, and thus the minimization for ww is equivalent to solving

It is easy to verify that the unique solutions of (Equation 11) are

where the convention 0·(0/0)=00·(0/0)=0 is followed. On the other hand, for a fixed ww, (Equation 10) is quadratic in uu and the minimizer uu is given by the normal equations

By noting the relation between DD and DiDi and a reordering of variables, (Equation 13) can be rewritten as

where

The normal equation (Equation 14) can also be solved easily provided that proper boundary conditions are assumed on uu. Since both the finite difference operations and the convolution are not well-defined on the boundary of uu, certain boundary assumptions are needed when solving (Equation 14). Under the periodic boundary conditions for uu, i.e. the 2D image uu is treated as a periodic function in both horizontal and vertical directions, D(1)D(1), D(2)D(2) and KK are all block circulant matrices with circulant blocks; see e.g. [10], [6]. Therefore, the Hessian matrix on the left-hand side of (Equation 14) has a block circulant structure and thus can be diagonalized by the 2D discrete Fourier transform FF, see e.g. [6]. Using the convolution theorem of Fourier transforms, the solution of (Equation 14) is given by

where the division is implemented by componentwise. Since all quantities but ww are constant for given ββ, computing uu from (Equation 16) involves merely the finite difference operation on ww and two FFTs (including one inverse FFT), once the constant quantities are computed.

Since minimizing the objective function in (Equation 10) with respect to either ww or uu is computationally inexpensive, we solve (Equation 10) for a fixed ββ by an alternating minimization scheme given below.

Algorithm :

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}N(K)∩N(D)={0}, where N(·)N(·) represents the null space of a matrix.

Define

Furthermore, we will make use of the following two index sets:

Under Assumption 1, the proposed algorithm has the following convergence properties.

Theorem 1 For any fixed β>0β>0, the sequence {(wk,uk)}{(wk,uk)} generated by Algorithm "Basic Algorithm" from any starting point (w0,u0)(w0,u0) converges to a solution (w*,u*)(w*,u*) of (Equation 10). Furthermore, the sequence satisfies

for all kk sufficiently large, where TEE=[Ti,j]i,j∈E∪(n2+E)TEE=[Ti,j]i,j∈E∪(n2+E) is a minor of TT, ∥v∥M2=v⊤Mv∥v∥M2=v⊤Mv and ρ(·)ρ(·) is the spectral radius of its argument.

Multichannel image deconvolution

Let u¯=[u¯(1);...;u¯(m)]∈Rmn2u¯=[u¯(1);...;u¯(m)]∈Rmn2 be an mm-channel image, where, for each jj, u¯(j)∈Rn2u¯(j)∈Rn2 represents the jjth channel. An observation of u¯u¯ is modeled by (Equation 4), in which case f=[f(1);...;f(m)]f=[f(1);...;f(m)] and ω=[ω(1);...;ω(m)]ω=[ω(1);...;ω(m)] have the same size and the number of channels as u¯u¯, and KK is a multichannel blurring operator of the form

where Kij∈Rn2×n2Kij∈Rn2×n2, each diagonal submatrix KiiKii defines the blurring operator within the iith channel, and each off-diagonal matrix KijKij, i≠ji≠j, defines how the jjth channel affects the iith channel.

The multichannel extension of (Equation 9) is

where ImIm is the identity matrix of order mm, and “⊗⊗" is the Kronecker product. By introducing auxiliary variables wi∈R2mwi∈R2m, i=1,...,n2i=1,...,n2, (Equation 22) is approximated by

For fixed uu, the minimizer function for ww is given by (Equation 12) in which DiuDiu should be replaced by (Im⊗Di)u(Im⊗Di)u. On the other hand, for fixed ww, the minimization for uu is a least squares problem which is equivalent to the normal equations

where ww is a reordering of variables in a similar way as given in (Equation 15). Under the periodic boundary condition, (Equation 24) 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¯u¯ as the minimizer of a TV/L11 problem. For simplicity, we again assume u¯∈Rn2u¯∈Rn2 is a single channel image and the extension to multichannel case can be similarly done as in "Multichannel image deconvolution". The TV/L11 problem is

Since the data-fidelity term is also not differentiable, in addition to ww, we introduce z∈Rn2z∈Rn2 and add a quadratic penalty term. The approximation problem to (Equation 25) is

where β,γ≫0β,γ≫0 are penalty parameters. For fixed uu, the minimization for ww is the same as before, while the minimizer function for zz is given by the famous one-dimensional shrinkage:

On the other hand, for fixed ww and zz, the minimization for uu is a least squares problem which is equivalent to the normal equations

Similar to previous arguments, (Equation 28) can be easily solved by FFTs.

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. Figure 2 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-310-3.

Figure 2: Continuation vs. no continuation: u*u* is an “exact” solution corresponding to β=214β=214. The horizontal axis represents the number of iterations, and the vertical axis is the relative error ek=∥uk-u*∥/∥u*∥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 Figure 2, 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/L22, the framework is given below.

[FTVd]:

Figure 3: SNRs of images recovered from () for different ββ.

Figure 4: Results recovered from TV/L22. 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-310-3 is added to both blurred images. The left images are the blurry and noisy observations, and the right ones are recovered by FTVd.

(a) (b)

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β=20,21,⋯,218 in Figure 3. In this experiment, we used the same blur and noise as we used in the testing of continuation. As can be seen from Figure 3, the SNR values on both images essentially remain constant for β≥27β≥27. This suggests that ββ need not to be excessively large from a practical point of view. In our experiments, we set β0=1β0=1 and βmax=27βmax=27 in Algorithm "Practical Implementation". For each ββ, the inner iteration was stopped once an optimality condition is satisfied. For TV/L11 problems, we also implement continuation on γγ, and used similar settings as used in TV/L22.

Recovered Results

In this subsection, we present results recovered from TV/L22 and TV/L11 problems including (Equation 9), (Equation 25) 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. Figure 4 gives two examples of blurry and noisy images and the recovered ones, where the blurred images are corrupted by Gaussian noise, while Figure 5 gives the recovered results where the blurred images are corrupted by random-valued noise. For TV/L11 problems, we set γ=215γ=215 and β=210β=210 in the approximation model and implemented continuation on both ββ and γγ.

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

(a) (b)