Free Essay

In: Film and Music

Submitted By kxsurat

Words 12310

Pages 50

Words 12310

Pages 50

c 2008 Society for Industrial and Applied Mathematics

A New Alternating Minimization Algorithm for Total Variation Image Reconstruction∗

Yilun Wang†, Junfeng Yang‡, Wotao Yin†, and Yin Zhang†

Abstract. We propose, analyze, and test an alternating minimization algorithm for recovering images from blurry and noisy observations with total variation (TV) regularization. This algorithm arises from a new half-quadratic model applicable to not only the anisotropic but also the isotropic forms of TV discretizations. The per-iteration computational complexity of the algorithm is three fast Fourier transforms. We establish strong convergence properties for the algorithm including ﬁnite convergence for some variables and relatively fast exponential (or q-linear in optimization terminology) convergence for the others. Furthermore, we propose a continuation scheme to accelerate the practical convergence of the algorithm. Extensive numerical results show that our algorithm performs favorably in comparison to several state-of-the-art algorithms. In particular, it runs orders of magnitude faster than the lagged diﬀusivity algorithm for TV-based deblurring. Some extensions of our algorithm are also discussed. Key words. half-quadratic, image deblurring, isotropic total variation, fast Fourier transform AMS subject classiﬁcations. 68U10, 65J22, 65K10, 65T50, 90C25 DOI. 10.1137/080724265

1. Introduction. In this paper, we propose a fast algorithm for reconstructing images from blurry and noisy observations. For simplicity, we assume that the underlying images have square domains, but all discussions can be equally applied to rectangle domains. Let 2 2 2 u0 ∈ Rn be an original n×n grayscale image, K ∈ Rn ×n represent a blurring (or convolution) 2 2 operator, ω ∈ Rn be additive noise, and f ∈ Rn be an observation which satisﬁes the relationship (1.1) f = Ku0 + ω.

Given f and K, the image u0 is recovered from the model n2 (1.2)

∗

min u i=1

Di u +

μ Ku − f 2

2 2,

Received by the editors July 3, 2007; accepted for publication (in revised form) May 16, 2008; published electronically August 20, 2008. http://www.siam.org/journals/siims/1-3/72426.html † Department of Computational and Applied Mathematics, Rice University, 6100 Main, Houston, TX, 77005 (yilun.wang@rice.edu, wotao.yin@rice.edu, yzhang@rice.edu). The work of the ﬁrst and third authors was supported by NSF Career Grant DMS-0748839 and the latter author’s internal faculty research grant from the Dean of Engineering at Rice University. The fourth author’s work was supported in part by NSF grant DMS-0405831. ‡ Department of Mathematics, Nanjing University, 22 Hankou Road, Nanjing, Jiangsu Province, 210093, People’s Republic of China (jfyang2992@yahoo.com.cn). This author’s work was supported by the Chinese Scholarship Council during the author’s visit to Rice University. 248

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

249

where Di u ∈ R2 denotes the discrete gradient of u at pixel i and the sum Di u is the discrete total variation (TV) of u (the issue of boundary conditions will be discussed later). The origin of (1.2) and related results will be reviewed brieﬂy in section 1.2. Model (1.2) is also widely referred to as the TV/L2 model. In (1.2), the TV is isotropic if the norm · in the summation is the 2-norm, and anisotropic if it is the 1-norm. We emphasize here that, unlike most previous work, our method applies to both isotropic and anisotropic TV. However, we will treat only the isotropic case, · = · 2 , in detail because the treatment for the anisotropic case is completely analogous. Our algorithm is derived from the well-known variable-splitting and penalty techniques in optimization; speciﬁcally, at each pixel an auxiliary variable wi ∈ R2 is introduced to transfer Di u out of the nondiﬀerentiable term · 2 , and the diﬀerence between wi and Di u is penalized, yielding the following approximation model to (1.2): (1.3) min w,u i

wi

2

+

β 2

wi − Di u i 2 2

+

μ Ku − f 2

2 2,

with a suﬃciently large penalty parameter β. This type of quadratic penalty approach can be traced back as early as Courant [13] in 1943. It is well known that the solution of (1.3) converges to that of (1.2) as β → ∞. Clearly, the objective function in (1.3) is convex in (u, w). The motivation for this formulation is that, while either one of the two variables u and w is ﬁxed, minimizing the function with respect to the other has a closed-form formula with low computational complexity and high numerical stability. We will show that, for any ﬁxed β, the algorithm of minimizing u and w alternately has ﬁnite convergence for some variables and q-linear convergence rates for the rest. Moreover, the overall convergence of this algorithm can be signiﬁcantly accelerated by a theoretically well-justiﬁed continuation approach on the penalty parameter β. The objective function in (1.3) is quadratic in u and separable in w; therefore, problem (1.3) is half-quadratic according to Geman and Reynolds [17] and Geman and Yang [18]. Following the framework proposed in [17, 18], a number of half-quadratic models have been derived and analyzed. However, model (1.3) cannot be derived from (1.2) by the constructions given in [17, 18]; instead, a new approximate TV function must ﬁrst be introduced before applying the construction of [18] (see section 2.3 for more details). The approximate TV model (1.3) and the resulting alternating minimization algorithm were ﬁrst proposed in [42] without a convergence analysis. A similar split method has recently been proposed that uses Bregman iterations [20]. We have implemented the proposed algorithm in MATLAB and compared it with a C++ implementation of the lagged diﬀusivity algorithm [41], which is one of the fastest algorithms for solving (1.2). Numerical results presented in section 5 show that our algorithm is much faster than the lagged diﬀusivity algorithm especially when the blurring kernel is relatively large. Compared with a few other deblurring algorithms which solve diﬀerent models, including ForWaRD (a Fourier and wavelet shrinkage method [32]) and MATLAB Image Processing Toolbox functions “deconvwnr” and “deconvreg,” our algorithm consistently generates higher quality images in comparable running times. In the rest of this section, we will introduce the notation, give a brief review of regularization approaches, and then summarize the contributions and organization of this paper.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

250

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

2 2

1.1. Notation. Let D(j) ∈ Rn ×n , j = 1, 2, 3, . . . , be ﬁnite diﬀerence matrices and, in particular, let D(1) and D(2) be the ﬁrst-order ﬁnite diﬀerence matrices in the horizontal and vertical directions, respectively. For j ≥ 3, D(j) represents higher-order ﬁnite diﬀerence matrices. For scalars αi , vectors vi , and matrices Ai of appropriate sizes, i = 1, 2, we let T T a = (α1 ; α2 ) (α1 , α2 )T , v = (v1 ; v2 ) (v1 , v2 )T , and A = (A1 ; A2 ) (AT , AT )T . As is used 1 2 2×n2 is a two-row matrix formed by stacking the ith rows of D (1) and D (2) . in (1.2), Di ∈ R Let ρ(T ) be the spectral radius of matrix T and let P(·) be a projection operator. From here on, the norm · refers to the 2-norm unless otherwise indicated. Additional notation will be introduced as the paper progresses. 1.2. Existing regularization approaches. There are diﬀerent approaches based on statistics [14, 5], Fourier and/or wavelet transforms [25, 31], or variational analysis [37, 7, 11] for image deblurring. Among them the simplest is the maximum likelihood estimation, which solves the least squares problem minu Ku − f 2 and is also known as the inverse ﬁlter. However, the solution of this inverse ﬁlter, though it best matches the probabilistic behavior of the data, is often unacceptable because the image deblurring problems are usually severely ill-conditioned and the observed data usually contains noise. To stabilize the restoration process, some a priori knowledge about the unknown image is utilized through the addition of a regularization term, resulting in the model (1.4) min J(u) = Φreg (u) + u μ Ku − f 2

2

,

where Φreg (u) is the regularizer that enforces the a priori knowledge and the parameter μ is used to balance the two terms. Two classes of regularizers are well known. One is the Tikhonov-like class [40] including Φreg (u) = j D(j) u 2 , where D(j) ’s stand for some ﬁnite diﬀerence operators. In these cases, since the resulting objective functions J(u) are quadratic, it is relatively inexpensive to minimize them by solving linear systems of equations. However, since Tikhonov-like regularizers tend to make images overly smooth, they often fail to adequately preserve important image attributes such as sharp edges. Another well-known class of regularizers are based on TV, which was ﬁrst proposed for image denoising by Rudin, Osher and Fatemi in [38], and then extended to image deblurring in [37]. In comparison to Tikhonov-like regularizers, TV regularizers can better preserve sharp edges or object boundaries that are usually the most important features to recover. The TV/L2 model (1.2) has been widely studied (see, for example, [7, 8] and the references therein). The superiority of TV over Tikhonov-like regularization was analyzed in [1, 15] for recovering images containing piecewise smooth objects. However, due to the nondiﬀerentiability and nonlinearity of the TV functions, the TV/L2 model (1.2) is computationally more diﬃcult to solve. Despite numerous eﬀorts over the years, algorithms for solving the isotropic TV/L2 model are still much slower than those for solving Tikhonov-like regularization models. In this work, we aim to overcome this diﬃculty. The most important structure of the TV/L2 model is that both the blurring operator and the ﬁnite-diﬀerence operators can be treated as discrete convolutions under suitable boundary conditions. Determining how to best exploit this structure is crucial for achieving high computational eﬃciency.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

251

1.3. Contributions. This paper has three major contributions: the derivation of a new half-quadratic model and the accompanying algorithm, the analysis of the algorithm’s convergence properties, and the evaluation of its practical performance. 1. To the best of our knowledge, (1.3) is the ﬁrst half-quadratic method for the isotropic TV/L2 model whose leading computational requirement is three fast Fourier transforms (FFTs) per iteration. Previously, this level of complexity was achievable only for the anisotropic case. Our extension not only enables better image quality, but, more importantly, can be generalized to color (or multichannel) image reconstruction with TV regularization. 2. We established convergence results for our algorithm that are stronger than those proved for existing half-quadratic methods. We established ﬁnite convergence for some auxiliary variables and strong q-linear convergence rates for the other variables. The q-linear rates are stronger than usual because they depend on the spectral radii of submatrices rather than of the whole matrix as in usual cases. 3. Based on our convergence results, we propose a continuation scheme on the penalty parameter with a warm start. Our extensive numerical results have conﬁrmed that the proposed algorithm can be orders of magnitude faster than the lagged diﬀusivity method—one of the state-of-the-art methods for solving the TV/L2 model (1.2). In addition, we can extend our method to other problems involving TV regularization such as compressive sensing. These extensions will be discussed in section 6. 1.4. Organization. The rest of this paper is organized as follows. In section 2, we propose our alternating minimization algorithm, study optimality conditions, and describe its relations with previous works. Convergence properties of the proposed algorithm are analyzed in section 3. In section 4, we propose a continuation scheme, demonstrate its eﬃciency, and describe the implementation of our complete algorithm. Numerical results in comparison with the lagged diﬀusivity (LD) and some other typical image restoration methods are presented in section 5. In section 6, we discuss how our method can be extended to more general inverse problems, including some compressed sensing formulations. Finally, concluding remarks are given in section 7. 2. Basic algorithm, optimality, and related works. We ﬁrst introduce two auxiliary 2 vectors w1 , w2 ∈ Rn to approximate D(1) u and D(2) u, respectively, where we recall that 2 2 D(1) , D(2) ∈ Rn ×n represent the two ﬁrst-order forward ﬁnite diﬀerence operators with ap2 propriate boundary conditions. We denote w = (w1 ; w2 ) ∈ R2n and D = (D(1) ; D(2) ) ∈ 2 ×n2 , where the latter is the total ﬁnite diﬀerence operator. For i = 1, . . . , n2 , we let R2n wi = ((w1 )i ; (w2 )i ) ∈ R2 , which is an approximation of Di u = (D(1) u)i ; (D(2) u)i ∈ R2 . To keep wi close to Di u, we apply a quadratic penalty term, which is the second term in (1.3). 2.1. An alternating minimization algorithm. It is easy to minimize the objective function in (1.3) with respect to either u or w. For a ﬁxed u, the ﬁrst two terms in (1.3) are separable with respect to wi , so minimizing (1.3) for w is equivalent to solving, for i = 1, 2, . . . , n2 , (2.1) min wi wi +

β wi − Di u 2

2

,

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

252

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

for which the unique minimizer is given by the following two-dimensional (2D) shrinkage formula: Di u 1 , i = 1, . . . , n2 , (2.2) wi = max Di u − , 0 β Di u where the convention 0 · (0/0) = 0 is followed. For the anisotropic case · = · 1 , each component of wi is given by the simpler one-dimensional (1D) shrinkage formula: wi = max{|Di u| − 1/β, 0} sgn(Di u), where all operations are done componentwise. On the other hand, for a ﬁxed w, (1.3) is quadratic in u, and the minimizer u is given by the normal equations

T Di Di + i

μ T K K β

u= i T Di wi +

μ T K f β

or, equivalently, (2.3) (D(1) )T D(1) + (D(2) )T D(2) + μ T μ K K u = (D(1) )T w1 + (D(2) )T w2 + K T f. β β

Under the periodic boundary condition for u, (D(1) )T D(1) , (D(2) )T D(2) , and K T K are all block circulant (see [33, 21], for example). Therefore, the Hessian matrix on the left-hand side of (2.3) can be diagonalized by 2D discrete Fourier transform F. Using the convolution theorem of Fourier transforms, we can write (2.4) u = F −1 F(D(1) )∗ ◦ F(w1 ) + F(D(2) )∗ ◦ F(w2 ) + (μ/β)F(K)∗ ◦ F(f ) F(D(1) )∗ ◦ F(D(1) ) + F(D(2) )∗ ◦ F(D(2) ) + (μ/β)F(K)∗ ◦ F(K) ,

where “∗” denotes complex conjugacy, “◦” denotes componentwise multiplication, and the division is componentwise as well. With a slight abuse of notation, we have used F(K) for the Fourier transform of the function represented by K in the convolution Ku (and similarly for D(1) and D(2) ). Since all quantities but w1 and w2 are constant, computing u from (2.4) involves two FFTs and one inverse FFT, once the constant quantities are computed. Alternatively, under the Neumann boundary condition and assuming that the blurring kernel K is symmetric, the left-hand side of (2.3) is a block Toeplitz-plus-Hankel matrix (see [33]), so (2.3) can be diagonalized by 2D discrete cosine transforms (DCTs), and solving (2.3) requires three DCT calls. In our numerical experiments, we used the periodic boundary condition and FFTs. Since minimizing the objective function in (1.3) with respect to either w or u is computationally inexpensive, we solve (1.3) for a ﬁxed β by an alternating minimization scheme given below. Algorithm 1. Input f , K, μ > 0, and β > 0. Initialize u = f . While “not converged,” Do (1) Compute w according to (2.2) for ﬁxed u. (2) Compute u according to (2.4) for ﬁxed w. End Do More details about this algorithm will be presented later. Next, we derive optimality conditions for (1.3), based on which we will specify the stopping criterion (2.10) for Algorithm 1.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

253

2.2. Optimality conditions. To study the optimality conditions of (1.3) for a ﬁxed β, we need the following proposition. Ax is Proposition 2.1. For any A ∈ Rp×n , the subdiﬀerential of f (x) (2.5) ∂f (x) = {AT Ax/ Ax } AT h : h ≤ 1, h ∈ Rp if Ax = 0, otherwise.

Proof. By the deﬁnition of subdiﬀerential for a convex function, (2.6) ∂f (x) {g ∈ Rn : Ay − Ax ≥ g T (y − x) ∀ y ∈ Rn }.

If Ax = 0, then f is diﬀerentiable at x, and clearly (2.5) holds. Now assume Ax = 0. In this case, after setting y = x + d, we have ∂f (x) = {g ∈ Rn : Ad ≥ g T d ∀ d ∈ Rn }. We will show that ∂f (x) ≡ S AT h : h ≤ 1, h ∈ Rp . First, for any AT h ∈ S, (AT h)T d ≤ Ad by the Cauchy–Schwarz inequality, implying S ⊂ ∂f (x). Next, we show ∂f (x) ⊂ S by contradiction. Suppose that g ∈ ∂f (x), but g ∈ S. Since S is closed and convex, by the well/ known separation theorem, there exist d ∈ Rn and τ ∈ R such that the hyperplane dT x = τ separates g and S so that dT g > τ > dT p for all p ∈ S. It follows that dT g > τ ≥ sup{dT AT h : h ≤ 1, h ∈ Rp } = Ad , contradicting the assumption that g ∈ ∂f (x). Therefore, we have ∂f (x) = S, and the result is proved. Since the objective function in (1.3) is convex, a pair (w, u) solves (1.3) if and only if the subdiﬀerential of the objective function at (w, u) contains the origin, which gives the following optimality conditions in light of Proposition 2.1 with A = I: (2.7) (2.8) wi / wi + β(wi − Di u) = 0, β Di u ≤ 1, i ∈ I1 i ∈ I2 {i : wi = 0}, {i : wi = 0},

βDT (Du − w) + μK T (Ku − f ) = 0.

Our stopping criterion for Algorithm 1 is based on the optimality conditions (2.7) and (2.8). Let ⎧ i ∈ I1 , ⎨ r1 (i) (wi / wi )/β + wi − Di u, Di u − 1/β, i ∈ I2 , r2 (i) ⎩ T (Du − w) + μK T (Ku − f ), r3 βD

(2.9)

where I1 and I2 are deﬁned as in (2.7). Algorithm 1 is terminated once (2.10) Res max max{ r1 (i) }, max{r2 (i)}, r3 i∈I1 i∈I2 ∞

≤

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

254

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

is met, where Res measures the total residual and > 0 is a prescribed tolerance. In (2.9), condition (2.7) is scaled by a factor of 1/β, but (2.8) is not, because in practice the latter can be met quickly even without this scaling. Combining (2.7) and (2.8) to eliminate the w-variable, we can derive (2.11) i∈I1 T Di

Di u + Di u

T Di hi + μK T (Ku − f ) = 0, i∈I2

∗ where hi βDi u satisﬁes hi ≤ 1. Let u∗ be any solution of (1.2). Deﬁne I1 = {i, Di u∗ = 0} ∗ = {1, . . . , n2 } \ I ∗ . Then, from Proposition 2.1, there exist {h∗ ∈ R2 : h∗ ≤ 1, i ∈ I ∗ } and I2 1 2 i i such that

(2.12)

∗ i∈I1

T Di

Di u∗ + Di u∗

T Di h∗ + μK T (Ku∗ − f ) = 0. i

∗ i∈I2

Equation (2.11) diﬀers from (2.12) only in the index sets over which summations are taken. ∗ As β increases, I1 should approach I1 . In section 3, we will study the convergence properties of Algorithm 1. 2.3. Related works. While (1.3) is an instance of the classical quadratic penalty method, it is also a half-quadratic model according to Geman and Reyonlds [17] and Geman and Yang [18]. However, this model cannot be directly derived from the original frameworks in [17, 18] without introducing a new approximation to the isotropic TV function. Consider the following general framework of recovering an image u from its corrupted measurements f : (2.13)

2

min u i

T φ(gi u) +

μ Ku − f 2

2

,

T where gi ∈ Rn is a local ﬁnite diﬀerence operator, φ(gi ·) is convex and edge-preserving, and K is a convolution operator. Instead of solving (2.13) directly, the authors of [17] and [18] propose solving the equivalent problem

(2.14)

min u,b i

T Q(gi u, bi ) + ψ(bi ) +

μ Ku − f 2

2

,

where Q(t, s) and ψ(s) are chosen such that Q(t, s) is quadratic in t and (2.15) φ(t) = min(Q(t, s) + ψ(s)) ∀t ∈ R. s∈R The name “half-quadratic” comes from the fact that the objective function in (2.14) is quadratic in u but not in b. As such, under certain conditions such as convexity (2.14) can be solved by minimizing with respect to u and b alternately. Computationally, it is important that the constructed Q and ψ allow fast solutions for u and b, respectively. For this purpose, two forms of half-quadratic formulations have been widely studied: the multiplicative form [17] in which Q(t, s) = 1 t2 s and the additive form [18] in which Q(t, s) = 1 (t − s)2 . For 2 2

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

255

comparisons of these two forms of methods in both theory and practice, see [24, 2, 34, 12] and the references therein. Unfortunately, the constructions given in [17] and [18] cannot be directly applied to either the anisotropic or the isotropic TV because some technical conditions fail to hold. Over the years, numerous half-quadratic models have been proposed to approximate the TV functions. Some of these models are of the additive form [24, 34], while others are of the multiplicative form [7, 24]. However, in the half-quadratic models of the multiplicative form, the partial Hessians with respect to u depend on b and thus vary from one iteration to another. In addition, such partial Hessians do not have a block circulant structure. As a result, halfquadratics of the multiplicative form cannot be eﬃciently minimized by FFTs. Our half-quadratic model in (1.3) approximates the isotropic TV function and is of the additive form. The Hessian of the quadratic is constant and has a block circulant structure (under appropriate boundary conditions). The only other half-quadratic model for TV approximation that shares these attributes is the one by Nikolova and Ng [34], but it is for the anisotropic rather than the isotropic TV function. In [34], the authors use the Huber function

1 2 2 t

(2.16)

φ(t) =

|t| −

2

if |t| ≤ , otherwise,

where 0 < 1, to approximate φ(t) = |t| in (2.13) and apply the additive form halfquadratic construction. The resulting problem (2.14) is quadratic in u and separable in b, and the Hessian of the quadratic is constant and has a block circulant structure, enabling the use of fast transforms. However, [34] does not include a convergence analysis or extensive numerical experiments on the half-quadratic model for this anisotropic TV approximation. In [34], two approaches are proposed to approximate the isotropic TV function. One is based on the multiplicative construction; the other employs the additive construction but does not yield a half-quadratic model. Neither approach leads to a fast alternating minimization algorithm. Thus far all of the edge-preserving functions φ used in half-quadratic constructions that we are aware of take only scalar variables. In order to approximate the isotropic TV function, one should allow φ to take vector variables; then our model (1.3) can be derived from a similar additive construction, which we describe below. First, (1.2) is approximated by (2.17) min u i

φ(Di u) +

μ Ku − f 2

2

,

where φ : R2 → R is a 2D Huber-like function deﬁned as φ(t) = β 2

t

2 1 2β

t −

if t ≤ 1/β, otherwise

for β 0. Then a similar additive construction in two dimensions, instead of in one dimension, gives Q(t, s) = β t − s 2 and ψ(s) = s , where s, t ∈ R2 , which satisfy the 2D version of 2 (2.15). Clearly, s and β t − s 2 are nothing but the ﬁrst and second terms in our model 2 (1.3) for s = wi and t = Di u. Therefore, (1.3) turns out to be a long overdue extension to Nikolova and Ng’s half-quadratic model [34].

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

256

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

Since the work of Geman and Reynolds [17] and Geman and Yang [18], the convergence of alternating minimization algorithms for various half-quadratic models has been extensively studied, including that for multiplicative formulations in [24, 2, 34, 12] and for additive formulations in [24, 2, 34]. However, these analyses were done under some strict convexity assumptions on objective functions; see [24, 34, 2, 12]. Moreover, the convergence rate results obtained were r-linear at best; see [34, 2], for example. Unlike these previous results, our convergence results, presented in the next section, require neither nonsingularity of the operator K nor strict convexity of the objective function as a whole (notice that φ = · 2 is not strictly convex). Moreover, we establish a ﬁnite convergence property for some variables and strong q-linear convergence for the rest. 3. Convergence analysis. The convergence of the quadratic penalty method as the penalty parameter goes to inﬁnity is well known (see Theorem 17.1 in [35], for example). That is, as β → ∞, the solution of (1.3) converges to that of (1.2). However, it is generally diﬃcult to determine theoretically how large a β value must be to attain a given accuracy. We will later choose β based on numerical experiments. In this section, we analyze the convergence properties of Algorithm 1 for a ﬁxed β > 0. First, we prove that the sequence {(wk , uk )} generated by Algorithm 1 from any initial point converges to a solution of (1.3). Then, using a ﬁnite convergence property for {wk }, we establish q-linear convergence rates for Algorithm 1. In what follows, we omit β for notational simplicity. For a ∈ R2 , the 2D shrinkage operator s : R2 → R2 is deﬁned as (3.1) s(a) max a − 1 ,0 β a , a

where the convention 0 · (0/0) = 0 is followed. It is easy to see that s(a) = a − P(a), where P(·) PB (·) : R2 → R2 is the projection onto the closed disc B For vectors u, v ∈ RN , N ≥ 1, we deﬁne S(u; v) : R2N → R2N by S(u; v) (s(a1 ); . . . ; s(aN )) , where ai = ui vi ; {a ∈ R2 | a ≤ 1/β}.

i.e., S applies 2D shrinkage to each pair (ui , vi ) ∈ R2 , for i = 1, 2, . . . , N . The ﬁrst convergence result is Theorem 3.4, and we take a few steps to prove it. The ﬁrst step is to prove the nonexpansiveness of the shrinkage operator s. The following proposition is known to hold for projections P onto any closed convex set. We include a proof for the sake of completeness. Proposition 3.1. For any a, b ∈ R2 , it holds that s(a) − s(b)

2

≤ a−b

2

− P(a) − P(b) 2 .

Furthermore, if s(a) − s(b) = a − b , then s(a) − s(b) = a − b.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

257

Proof. Since B is convex, P(·) satisﬁes (a − P(a))T (P(a) − P(b)) ≥ 0 ∀ a, b ∈ R2 .

Exchanging a and b gives us −(b − P(b))T (P(a) − P(b)) ≥ 0. Adding the two yields (3.2) (a − b)T (P(a) − P(b)) ≥ P(a) − P(b) 2 .

Therefore, from s(a) = a − P(a), we get s(a) − s(b)

2

= a − b − (P(a) − P(b)) = a−b ≤ a−b

2 2

2 2

− 2(a − b)T (P(a) − P(b)) + P(a) − P(b) − P(a) − P(b) 2 ,

where the last inequality follows from (3.2). Moreover, if s(a) − s(b) = a − b , then P(a) = P(b) and s(a) − s(b) = a − b − (P(a) − P(b)) = a − b. We will make use of the following assumption in our convergence analysis. It is a mild assumption that has been used in most, if not all, previous analysis of a similar type. Assumption 1. N (K) ∩ N (D) = {0}, where N (·) represents the null space of a matrix. The following two symmetric positive deﬁnite matrices will be used in our analysis: (3.3) M = DT D + μ T K K and T = DM −1 DT . β

Assumption 1 ensures the nonsingularity of M ; thus T is well deﬁned. It is worth noting that 2 2 ρ(T ) ≤ 1. We also deﬁne a linear operator h : R2n → R2n as h(w) = (h(1) (w); h(2) (w)), where μ h(j) (w) = D(j) M −1 DT w + K T f , j = 1, 2. β Using the deﬁnitions of S and h, we can rewrite the two iterative steps in Algorithm 1 as (3.4) (3.5) wk+1 = S(D(1) uk ; D(2) uk ) = S ◦ h(wk ), μ uk+1 = M −1 DT wk+1 + K T f . β

Since the objective function in (1.3) is convex, bounded below, and coercive (i.e., its value goes to inﬁnity as (w, u) → ∞), (1.3) has at least one minimizer pair (w∗ , u∗ ) that cannot be decreased by our alternating minimization scheme and thus must satisfy (3.6) (3.7) w∗ = S(D(1) u∗ ; D(2) u∗ ) = S ◦ h(w∗ ), μ u∗ = M −1 DT w∗ + K T f . β

Particularly, (3.6) means that w∗ is a ﬁxed point of S ◦ h. We need to show that h is nonexpansive. 2 Proposition 3.2. For any w = w in R2n , it holds that ˜ h(w) − h(w) ≤ w − w , ˜ ˜

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

258

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

and the equality holds if and only if h(w) − h(w) = w − w. ˜ ˜ Proof. From the deﬁnitions of M and T , it can be shown that ρ(T ) = ρ(DM −1 DT ) ≤ 1 and h(w) − h(w) = T (w − w) ≤ ρ(T ) w − w ≤ w − w . ˜ ˜ ˜ ˜ Let T = QT ΛQ be the eigendecomposition of T , where Q is an orthogonal matrix and Λ is a diagonal matrix with elements 0 ≤ λi ≤ 1, for all i. If h(w) − h(w) = w − w , namely, ˜ ˜ QT ΛQ(w − w) = w − w , then ΛQ(w − w) = Q(w − w) . Since Λ is diagonal and ˜ ˜ ˜ ˜ 0 ≤ λi ≤ 1, ΛQ(w − w) = Q(w − w) holds true. Multiplying both sides by QT completes the ˜ ˜ proof. The following lemma gives a useful property for ﬁxed points of the operator S ◦ h. Lemma 3.3. Let w be any ﬁxed point of S ◦ h. For any w, we have S ◦ h(w) − S ◦ h(w) < ˆ ˆ w − w unless w is a ﬁxed point of S ◦ h. ˆ Proof. From Propositions 3.1 and 3.2 and the deﬁnition of S, it holds that S ◦ h(w) − S ◦ h(w) ≤ w − w . If equality holds, then, again from Propositions 3.1 and 3.2, we get ˆ ˆ S ◦ h(w) − S ◦ h(w) = h(w) − h(w) = w − w. ˆ ˆ ˆ Since w = S ◦ h(w), we get w = S ◦ h(w). ˆ ˆ Now we are ready to prove the convergence of Algorithm 1. Theorem 3.4. For any ﬁxed β > 0, the sequence {(wk , uk )} generated by Algorithm 1 from any starting point (w0 , u0 ) converges to a solution (w∗ , u∗ ) of (1.3). Proof. We prove the convergence of {wk } in three steps. First, from (3.4) and the nonexpansiveness of S and h, it is easy to show that the sequence {wk } lies in a compact region and thus has at least one limit point, say w∗ = limj→∞ wkj . Letting w be any ﬁxed point of ˆ S ◦ h, i.e., w = S ◦ h(w), we get ˆ ˆ ˆ ˆ ˆ wk − w = S ◦ h(wk−1 ) − S ◦ h(w) ≤ wk−1 − w . Therefore, the following limit exists: (3.8) k→∞ ˆ ˆ ˆ lim wk − w = lim wkj − w = w∗ − w , j→∞ which means that all limit points of {wk }, if more than one, have an equal distance to w. By ˆ the continuity of S ◦ h, we have S ◦ h(w∗ ) = lim S ◦ h(wkj ) = lim wkj +1 . j→∞ j→∞

Hence, S ◦ h(w∗ ) is also a limit point of {wk } that must have the same distance to w as w∗ ˆ does; i.e., w∗ − w = S ◦ h(w∗ ) − w = S ◦ h(w∗ ) − S ◦ h(w) . ˆ ˆ ˆ It follows from Lemma 3.3 that w∗ = S ◦h(w∗ ). Since w is any ﬁxed point of S ◦h, by replacing ˆ w with w∗ in (3.8), we establish the convergence of {wk }: limk→∞ wk = w∗ . The convergence ˆ of {uk } to some u∗ follows directly from (3.5). Therefore, (w∗ , u∗ ) satisﬁes (3.6)–(3.7) and is a solution of (1.3).

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

259

We note that Theorem 3.4 does not require the objective function in (1.3) to be strictly convex in (u, w) (even though it is so with respect to each individual variable). As such, (1.3) may have multiple solutions, and the limit (w∗ , u∗ ) of the sequence generated by Algorithm 1 depends on starting point (w0 , u0 ). Next we develop a ﬁnite convergence property for the auxiliary variable w. Let hi (w) = (hi (w); hi (w)) ∈ R2 ,

(1) (2)

i = 1, . . . , n2 ;

namely, hi (w) is a vector formed by stacking the ith components of h(1) (w) and h(2) (w). We will make use of the following two index sets: (3.9) L= i : Di u∗ = hi (w∗ ) < 1 β and E = {1, . . . , n2 } \ L.

Theorem 3.5 (ﬁnite convergence). The sequence {(wk , uk )} generated by Algorithm 1 from k ∗ any starting point (w0 , u0 ) satisﬁes wi = wi = 0, for all i ∈ L, for all but a ﬁnite number of 0 − w ∗ 2 /ω 2 , where iterations that does not exceed w (3.10) ω min i∈L 1 − hi (w∗ ) β

> 0.

Proof. By the nonexpansiveness of s(·), for all i, (3.11) k+1 ∗ wi − wi = s ◦ hi (wk ) − s ◦ hi (w∗ ) ≤ hi (wk ) − hi (w∗ ) .

k+1 Suppose that at iteration k there exists at least one index i ∈ L such that wi = s ◦ hi (wk ) = ∗ ) < 1/β, h (w k ) ≥ 1/β, and w∗ = s ◦ h (w ∗ ) = 0. Therefore, 0. Then, hi (w i i i

(3.12)

k+1 ∗ wi − wi

2

= s ◦ hi (wk )

2

= ( hi (wk ) − 1/β)2

2 2

≤ { hi (wk ) − hi (w∗ ) − (1/β − hi (w∗ ) )}2 ≤ hi (wk ) − hi (w∗ ) ≤ hi (wk ) − hi (w∗ ) − (1/β − hi (w∗ ) )2 − ω2,

where the ﬁrst inequality is a triangular inequality, the second follows from the fact that hi (wk ) − hi (w∗ ) ≥ 1/β − hi (w∗ ) > 0, and the last uses the deﬁnition of ω in (3.10). Combining (3.11), (3.12), and the nonexpansiveness of h(·), we obtain wk+1 − w∗

2

= i k+1 ∗ wi − wi

2

≤ i hi (wk ) − hi (w∗ )

2

2

− ω2

= h(wk ) − h(w∗ )

2

− ω 2 ≤ wk − w∗

− ω2.

k+1 Therefore, the number of iterations k with wi = 0 does not exceed w0 − w∗ 2 /ω 2 . k = w∗ for i ∈ L, we next show the linear convergence of Given the ﬁnite convergence of wi i k and wk for i ∈ E. Denote w = ((w ) ; (w ) ), where (w ) and (w ) are the subvectors u 1 E 2 E 1 E 2 E E i of w1 and w2 , respectively, with components (w1 )i , (w2 )i , i ∈ E. Theorem 3.6 (q-linear convergence). Let M and T be deﬁned as in (3.3) and let TEE = [Ti,j ]i,j∈E∪(n2 +E) . Under Assumption 1, the sequence {(wk , uk )} generated by Algorithm 1 satisﬁes

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

260

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

k+1 ∗ k ∗ 1. wE − wE ≤ ρ((T 2 )EE ) wE − wE ; 2. D(uk+1 − u∗ ) ≤ ρ((T 2 )EE ) D(uk − u∗ ) ; 3. uk+1 − u∗ M ≤ ρ(TEE ) uk − u∗ M for all k suﬃciently large. Proof. From (3.4)–(3.7) and the nonexpansiveness of S, we get

(3.13) and (3.14) wk+1 − w∗

2

uk+1 − u∗ = M −1 DT (wk+1 − w∗ )

= S(D(1) uk ; D(2) uk ) − S(D(1) u∗ ; D(2) u∗ )

2

≤ D(uk − u∗ ) 2 .

Combining the recursion (3.13), (3.14) and the deﬁnition of T , it holds that wk+1 − w∗

2

≤ T (wk − w∗ ) 2 .

Since we are interested in the asymptotic behavior of Algorithm 1, without loss of generality, k k k ∗ hereafter we assume that wL ((w1 )L ; (w2 )L ) = wL = (0; 0). Therefore, the above inequality becomes k+1 ∗ wE − wE 2 k ∗ k ∗ k ∗ ≤ (wE − wE )T (T 2 )EE (wE − wE ) ≤ ρ((T 2 )EE ) wE − wE 2

.

k ∗ Multiplying D on both sides of (3.13), from wL = wL = 0 and (3.14), we get

D(uk+1 − u∗ )

2

≤ ρ((T 2 )EE ) wk+1 − w∗

2

≤ ρ((T 2 )EE ) D(uk − u∗ ) 2 .

k The above two inequalities imply assertions 1 and 2 of this theorem. From (3.13) and wL = ∗ wL = 0, we have

uk+1 − u∗

2 M

= (wk+1 − w∗ )T T (wk+1 − w∗ ) ≤ ρ(TEE ) wk+1 − w∗ 2 .

Considering (3.14) and the deﬁnition of M , the above inequality becomes uk+1 − u∗

M

≤

ρ(TEE ) D(uk − u∗ ) ≤

ρ(TEE ) uk − u∗

M;

namely, the last conclusion holds, and thus uk converges to u∗ q-linearly. Theorem 3.6 states that Algorithm 1 converges q-linearly with the convergence rate depending on the spectral radii of the submatrices TEE and (T 2 )EE rather than on that of the whole matrix T . Since ρ(T ) ≤ 1 and TEE is a minor of T , it holds that ρ(TEE ) ≤ ρ(T ) ≤ 1. Similarly, ρ((T 2 )EE ) ≤ ρ(T 2 ) ≤ 1. In the next section, we will demonstrate that in practice ρ(TEE ) can be much smaller than ρ(T ). 4. A continuation scheme. In this section, we develop a continuation scheme on the penalty parameter β based on our convergence results and demonstrate its eﬀectiveness. It is easy to see from (3.3) that a smaller β generally gives a smaller ρ(T ). As indicated in (3.9), a smaller β should also give a smaller set E and hence a smaller ρ(TEE ) (this fact will be veriﬁed by our numerical experiments). These observations strongly suggest the use of a continuation scheme on β in which β is initially small and gradually increases to the ﬁnal

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

261

value. In this scheme, earlier subproblems corresponding to smaller β values can be solved quickly, and the later subproblems can also be solved relatively quickly by warm starting from the previous solutions. To demonstrate the need for continuation and its eﬀectiveness, we carried out a set of experiments on the 128 × 128 image Checkerboard, a blocky image available in MATLAB. Our experiments were done under MATLAB using its Image Processing Toolbox. A motion blurring kernel of motion distance “len = 5” and angle “theta = 135” was applied to the image with additive Gaussian noise of zero mean and standard deviation 10−3 through MATLAB functions “fspecial,” “imfilter,” and “imnoise.” In our experiments, we set μ = 5 × 104 . In the ﬁrst experiment, we compared the quantities ρ(T ) and ρ(TEE ) for diﬀerent β values, both related to q-convergence rates. In addition, we also compared the theoretical upper bound, N0 w0 − w∗ /ω 2 , for ﬁnite convergence (see Theorem 3.5) with the observed k ˆ ˆ iteration number, N0 , for ﬁnite convergence (i.e., after N0 iterations wi , i ∈ L, remained zero until termination). To do these calculations, we ran Algorithm 1 for three β values: β = 24 , 29 , and 214 . For ∗ each β, an “exact” solution (u∗ , wβ ) of (1.3) was obtained using an extremely tight stopping β tolerance = 10−15 in (2.10). Afterwards, we reran the algorithm for each β and stopped ˆ once uk − u∗ / u∗ < 1.2 × 10−6 . The calculated values of ρ(T ), ρ(TEE ), N0 , and N0 β β for the three β values are given in Table 1.

Table 1 Comparison of theoretical convergence rates for diﬀerent β. β ρ(T ) ρ(TEE ) N0 ˆ N0 24 0.9999 0.5188 5.6850e+006 5 29 0.9999 0.9955 1.1990e+017 139 214 0.9999 0.9999 1.2103e+018 3658

Data in Table 1 show that on the tested image (i) ρ(TEE ) can be much smaller than ρ(T ) when β is small and (ii) the theoretical worst-case upper bound for ﬁnite convergence can be vastly overconservative in practice, especially when β is small. In Figure 1, we plot the histories of the relative error, uk − u∗ / u∗ , and the observed β β convergence factor, uk+1 − u∗ M / uk − u∗ M , generated by Algorithm 1 for the three tested β β β values. The plots clearly illustrate that the convergence speed of Algorithm 1 is much faster for smaller β values. Speciﬁcally, when β increased from 24 to 214 , the required iteration number for the given accuracy jumped from less than 10 to over 3000. Similarly, the observed q-convergence factor increased from around 0.2 to close to 1. Our theoretical and numerical results strongly suggest the use of a continuation scheme on β. How to do this continuation most eﬃciently should be an interesting research issue in its own right, but will not be studied in depth in the current paper. Instead, we implement a very basic version, likely still far from optimal, to test the viability of our proposed framework. We call the resulting algorithm “fast total variation deconvolution” or FTVd, which is given below.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

262

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

Relative error

10

0

Observed convergence rate β=24 1 0.9 0.8 0.7 0.6 0.5 β=24 β=29 β=214 β=2

9

10

−1

β=214

10

−2

10

−3

10

−4

0.4 0.3 0.2

10

−5

10

−6

10

0

10

1

10

2

10

3

10

4

0.1 0 10

10

1

10

2

10

3

10

4

Figure 1. Convergence behavior for β = 24 , 29 , and 214 . The left plot is the relative error ek (β) = uk −u∗ / u∗ , and the right one is the observed q-linear convergence factor qk (β) = uk+1 −u∗ M / uk −u∗ M . β β β β In each plot, the horizontal axis represents the number of iterations.

Algorithm 2 (FTVd). Input f , K, μ > 0, β0 > 0, and βmax > β0 . Initialize u = f , up = 0, β = β0 , and > 0. While β ≤ βmax , Do (1) Run Algorithm 1 until (2.10) is met. (2) β ← 2 ∗ β. End Do One can modify the above framework to make it more ﬂexible, though as it is the above basic implementation already works surprisingly well. For example, one could adaptively increase β and choose from one outer iteration to another (or use another stopping criterion for Algorithm 1). We tested the eﬀectiveness of FTVd for β0 = 1 and βmax = 214 and in Figure 2 present the results in comparison to the results without doing continuation. From this comparison, it is abundantly clear that the algorithm would not be practically eﬀective without continuation. 5. Numerical experiments. In this section, we present detailed numerical results comparing FTVd to the LD algorithm [41, 10], a state-of-the-art method for solving the isotropic TV deblurring model (1.2), as well as to some non-TV deblurring algorithms. 5.1. Overall assessment of several algorithms. We tested several popular algorithms for solving the TV deblurring model (1.2). Based on our experience, we have reached the following overall impression. The artiﬁcial time-marching algorithm used by Rudin and Osher in [37] and by Rudin, Osher, and Fatemi in [38] is easy to implement, but usually takes a large number of iterations, each with a small step size governed by the CFL condition, to reach a modest accuracy. On the other hand, the second-order cone programming (SOCP) approach [19] recovers solutions with a high accuracy in a small number of iterations (typically 30–50), but takes much more running time and memory in each iteration even on medium-sized images. A comprehensive comparison between the artiﬁcial time-marching algorithm and the SOCP approach is presented in [19]. We did not test the interior-point primal-dual implicit quadratic methods

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

263

Relative error

10

0

Observed convergence rate

1 0.9 0.8 0.7 0.6

with continuation without continuation 10

−1

10

−2

0.5 0.4

10

−3

0.3 0.2 0.1 with continuation without continuation 10

1

10

−4

0

500

1000

1500

2000

2500

3000

0 0 10

10

2

10

3

10

4

Figure 2. Continuation vs. no continuation: u∗ is an “exact” solution corresponding to β = 214 . The left plot is the relative error ek = uk − u∗ / u∗ , and the right one is the observed q-linear convergence factor qk = uk+1 − u∗ M / uk − u∗ M . In each plot, the horizontal axis represents the number of iterations.

presented in [9, 3, 8, 30], but expect their performance to be similar to SOCP because they also require solving large systems of linear equations at each iteration. We also tested the recently proposed iteratively reweighted least squares algorithm [43] and obtained speed and accuracy consistent with that reported by the authors. Speciﬁcally, it was slower than the LD method devised by Vogel and Oman [41, 10] for solving the TV/L2 model (1.2). Therefore, among all tested algorithms, LD is the most eﬃcient for solving (1.2). Clearly, there are other deblurring algorithms not included in our discussion, but a more exhaustive comparison is beyond the scope of this work. 5.2. Test images and test platforms. We used two images, Lena and Man, in our experiments; see Figure 3. Image Lena is a good test image because it has a nice mixture of detail, ﬂat regions, shading area, and texture. Image Man also consists of complex components in diﬀerent scales, with diﬀerent patterns and under inhomogeneous illuminations. Both images are suited for our experiments. We tested several kinds of blurring kernels including Gaussian, average, and motion, and found that the running times of FTVd and LD remained essentially constant for diﬀerent blurring kernels as long as other conditions (e.g., kernel size or image size) were the same. Therefore, without loss of generality, we present only results using a “Gaussian” blurring kernel with diﬀerent kernel sizes. In all tests, the additive noise used was Gaussian noise with zero mean and standard deviation 10−3 . This level of noise is substantial in the context of deblurring. Detailed information of our experiment set-up is summarized in Table 2. FTVd was implemented in MATLAB. All blurring eﬀects were generated using the MATLAB function “imfilter” with periodic boundary conditions. For LD, we used an eﬃcient C++ code NUMIPAD [36]. NUMIPAD uses Dirichlet boundary conditions for TV calculation (i.e., areas outside of the image boundary are treated as 0) instead of periodic boundary conditions. This diﬀerence caused the two codes to generate images with minor diﬀerences near boundaries, but should not noticeably aﬀect their running times.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

264

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

Figure 3. Test images: Lena (left) and Man (right).

Table 2 Information on test images and blurring setting. Test no. 1. 2. Image Lena Man Size 512 × 512 1024 × 1024 Blurring type Gaussian Gaussian Blurring kernel parameters hsize={3, 5, . . . , 19, 21}, σ = 10 hsize={3, 5, . . . , 19, 21}, σ = 10

All comparisons between FTVd and LD were performed under GNU/Linux Release 2.6.955.0.2 and MATLAB v7.2 (R2006a) running on a Dell Optiplex GX620 with Dual Intel Pentium D CPU 3.20GHz (only one processor was used by MATLAB) and 3 GB of memory. The C++ code of LD was a part of the library package NUMIPAD [36] and was compiled with gcc v3.4.6. Since the code of ForWaRD had compiling problems on our Linux workstation, we chose to compare FTVd with ForWaRD, and also MATLAB functions “deconvwnr” and “deconvreg” under Windows XP and MATLAB v7.5 (R2007b) running on a Lenovo laptop with an Intel Core 2 Duo CPU at 2 GHz and 2 GB of memory. 5.3. Parameter values. In model (1.2), the parameter μ controls the amount of penalty applied to the squared L2 -distance between Ku and f . According to (1.1), an appropriate μ should give a solution u of (1.2) satisfying Ku − f ≈ ω . In our experiments, we used μ = 0.05/σ 2 , where σ is the standard deviation of the additive Gaussian noise ω in (1.1). This formula is based on the observation that μ should be inversely proportional to the noise variance, while the constant 0.05 was determined empirically so that the restored images had reasonable signal-to-noise ratios (to be deﬁned later) and relative errors. More practical techniques exist for choosing μ, often by testing on a small window of image. The reader interested in this topic is referred to [39]. The issue of how to optimally select μ is important, but is outside of the scope of this work. As is usually done, we measure the quality of restoration by the signal-to-noise ratio

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

265

23 22 21 20

SNR (dB)

19 18 17 16 15 14 13 0 2 4 6 8 10 12 14 Lena Man

log β

2

16

18

Figure 4. SNRs of images recovered from (1.3) for diﬀerent β.

(SNR), deﬁned by SNR 10 ∗ log10 u0 − u ˜ 0−u u

2 2

,

where u0 is the original image and u is the mean intensity value of u0 . Generally, the quality ˜ of the restored image is expected to increase as β increases (hence, (1.3) becomes a closer approximation to (1.2)). In practice, we observed that the SNR value of recovered images from (1.3) stabilized once β reached a reasonably large value. To see this, we plot the SNR values of restored images corresponding to β = 20 , 21 , . . . , 218 in Figure 4. In this experiment, “Gaussian” out-of-focus blurring was applied to Lena with blurring size hsize = 11 and standard deviation σ = 5, and the “motion” blurring was applied to Man with motion distance “len = 21” and “theta = 135.” As can be seen from Figure 4, the SNR values on both images essentially remain constant for β ≥ 27 . This suggests that β need not be excessively large from a practical point of view. In our numerical experiments comparing FTVd with LD, we set β0 = 1 and βmax = 27 in Algorithm 1. For each β, the inner iteration was stopped once (2.10) was satisﬁed with = 0.05. 5.4. Comparison of FTVd and LD. In this subsection, we compare FTVd with LD in terms of both convergence speed and quality of recovery. As mentioned in subsection 5.2, the two methods use diﬀerent boundary conditions, which causes minor diﬀerences. The blurry images with the circular boundary condition usually have slightly higher SNRs than those blurred with the Neumann boundary condition, as do the restored images. To have a fairer comparison, we measure the restoration quality by the so-called improvement SNR (ISNR) deﬁned as f − u0 2 , ISNR 10 ∗ log10 u − u0 2

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

266

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

9

10

4

8.5 10 8

3

FTVd: Test No.1 LD: Test No.1 FTVd: Test No.2 LD: Test No.2

Running times (s)

FTVd: Test No.1 LD: Test No.1 FTVd: Test No.2 LD: Test No.2 4 6 8 10 12 14 16 18 20

ISNR (dB)

7.5

10

2

7

10

1

6.5

6

10

0

4

6

8

10

12

14

16

18

20

hsize

hsize

Figure 5. ISNR (left) and CPU time (right) comparisons for FTVd and LD methods. In both ﬁgures, the horizontal axis represents the size of the blurring kernel.

where u0 , u, and f are the original, restored, and observed images, respectively. This quantity measures the quality of a restored image u relative to the blurry and noisy observation f . The ISNRs of the recovered images computed by the FTVd and LD methods for the two test images are given in the left-hand plot of Figure 5. We observe that FTVd produces slightly lower ISNRs when hsize is small and slightly higher ISNRs when hsize is large. However, since both FTVd and LD solve the same model (boundary conditions aside), the restored images have few visible diﬀerences. Therefore, we will not display the restored images here. The right-hand plot in Figure 5 shows the comparison of running time between FTVd and LD. While generating the restored images with similar qualities, our MATLAB code FTVd is much faster than the C++ code for LD. The running time of FTVd essentially remains constant for diﬀerent values of hsize and increases only moderately when image size is doubled. However, the running time for LD increases dramatically with the increase in both the image size and the blurring kernel size. As hsize becomes relatively large, the CPU times consumed by FTVd and LD become orders of magnitude apart. The reason is clear: Larger hsize makes the matrix K denser and more ill-conditioned; hence, solving a linear system involving K becomes more diﬃcult for the preconditioned conjugate gradient (CG) method used in NUMIPAD [36]. On the other hand, the increase of kernel size does not noticeably increase the cost of doing FFT on the function represented by K. 5.5. Observations on FTVd’s behavior and a note. As has been mentioned, the periteration computation of FTVd is dominated by three FFTs, each at the cost of O(n2 log(n)). The question is how many iterations are needed in general for FTVd to attain a required accuracy for images of diﬀerent sizes. In our experiments, the number of total inner iterations taken by the basic version of FTVd, with the given “default” parameters, is almost always around 12. These total inner iteration numbers are quite reasonable, given that there were 8 outer iterations corresponding to the β-sequence {20 , 21 , . . . , 27 }. For each β > 1, since FTVd has already obtained a good starting point, on average it takes only one or two inner

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

267

iterations to reach the prescribed accuracy. Given that three FFTs are required at each inner iteration, the total number of FFTs taken by FTVd is around 40 for all tests. Upon examining detailed timing information, we observed that over two-thirds of the CPU time was spent on the FFT and related calculations in (2.4). Finally, FTVd is numerically stable and insensitive to ill-conditioning, apparently because it does not require any matrix operations. After we ﬁnished the writing of this paper, we found a newly released paper by Huang, Ng and Wen [23] which solves the following approximation to the TV/L2 model (after necessary notational changes for consistency): (5.1) min u,v i

Di v + α v − u

2

+

μ Ku − f 2

2

,

by alternating minimization. For a ﬁxed v, the minimization with respect to u involves two FFTs, one less than what is required by the corresponding subproblem in our method. However, for a ﬁxed u, the minimization with respect to v is nothing but a TV denoising problem which does not have a closed-form solution with linear complexity as our corresponding subproblem has (see formula (2.2)). In [23], the TV denoising problem is solved iteratively by Chambolle’s projection algorithm [6]. While the per-iteration computational complexity of our method is dominated by three FFTs, that of [23] is dominated by the cost of solving a TV denoising problem in addition to two FFTs. According to the reported numerical results in [23], their algorithm appears to require at least as many outer iterations as ours. 5.6. Comparison of FTVd with other methods. We also conducted tests to show how FTVd compares with some other state-of-the-art, non-TV methods. A hybrid Fourier-wavelet regularized deconvolution (ForWaRD) algorithm by Neelanmani, Choi, and Baraniuk [32] uses shrinkage in both Fourier and wavelet domains to recover images while preserving edges and removing noise. Their code was written in the C programming language with a MATLAB interface.1 We compared the quality of images recovered by FTVd, ForWaRD, and two MATLAB deblurring functions “deconvwnr” (Wiener ﬁlter) and “deconvreg” (Laplacian lowpass ﬁlter). These methods are all based on distinct models. In this experiment, we used the image Lena, which was blurred with a “Gaussian” kernel of hsize = 21 and standard deviation σ = 11. We ran FTVd twice, once with β = 25 and tolerance = 0.05 in (2.10), and again with β = 27 and = 0.002. The blurry and noisy and the restored images are presented in Figure 6. SNRs of the restored images and running times of diﬀerent algorithms are also given. As can be seen from Figure 6, FTVd is comparable to ForWaRD in speed and attains a better image quality. The MATLAB functions are faster as they solve simpler models, but they generate images of lower quality. Speciﬁcally, the artifacts of multiple rings or ripples are clearly more visible in the restored images by the three methods than in the ones by FTVd. By visually comparing the restored images, we give preference to the results of FTVd. 5.7. Summary of numerical results. First, FTVd can be orders of magnitude faster than the LD method. The performance gap between them will continue to grow as the image and blurring kernel sizes increase. Since LD compared favorably to many other existing algorithms

1

The code ForWaRD can be downloaded at http://www.dsp.rice.edu/software/ward.shtml.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

268

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

Blurry&Noisy. SNR: 5.19dB ForWaRD: SNR: 12.46dB, t = 4.88s

FTVd: β = 25, SNR: 12.58dB, t = 1.83s

FTVd: β = 2 , SNR: 13.11dB, t = 14.10s

7

deconvwnr: SNR: 11.51dB, t = 0.05s

deconvreg: SNR: 11.20dB, t = 0.34s

Figure 6. The ﬁrst row contains the blurry and noisy image and the image recovered by ForWaRD; the second row contains the results of FTVd (left: β = 25 , = 5 × 10−2 ; right: β = 27 , = 2 × 10−3 ); and the third row is recovered by deconvwnr (left) and deconvreg (right).

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

269

for solving (1.2) with isotropic TV, the lead of FTVd obviously extends to those algorithms. Second, compared to some other non-TV deblurring algorithms such as ForWaRD and the MATLAB deblurring functions, FTVd generally produces images of higher quality in our experiments on numerous natural images. Under the “default” parameter setting, the speed of FTVd is not as fast as the MATLAB functions “deconvwnr” and “deconvreg.” However, with a less restrictive stopping tolerance and a smaller βmax , the running time of FTVd can be faster than that of “deconvreg” on large images and close to that of “deconvwnr” while still producing images of higher quality. 6. Extensions. There are some obvious extensions in which one can apply the variablesplitting and alternating-minimization methodology to solving other problems. Here we mention a few. Color or other multichannel image reconstruction with TV regularization. This extension is practically important, as straightforward as it may be, because it enables multichannel TV image reconstruction problems to be solved at a speed not seen before. We are currently working on this extension and will report our ﬁndings in a forthcoming paper. In addition, we can also extend this methodology to models where the ﬁdelity term Ku − f is not measured in the 2-norm but in the 1-norm in order to handle impulsive noise (such as “salt and pepper”) other than white noise. The use of higher-order ﬁnite diﬀerence operators. In this case, we may have models in the form of (6.1) min u i

Di u +

(p)

μ Ku − f 2

2

,

for some p > 2, which can be approximated by (6.2) min w,u i

wi +

β (p) wi − Di u 2

(p)

2

+

μ Ku − f 2

2

.

As long as the diﬀerence operators Di are linear and shift-invariant, minimization with respect to u is as simple as a few FFTs for a ﬁxed w, while minimization with respect to w can be done by shrinkage in an appropriate dimension (depending on the order of the ﬁnite diﬀerence involved). Higher-order derivatives of u are used in [27, 44], for example, to reduce staircasing eﬀects. Operator K is a partial Fourier transform matrix. In this case, Ku = P F(u), where 2 P ∈ Rm×n is a selection matrix consisting of m < n2 rows of the identity matrix. This case occurs, for example, in compressive sensing [4, 16] where a subset of noisy Fourier coeﬃcients, f , is used to recover a signal u0 with a sparse gradient ﬁeld via solving the model (6.3) min u i

Di u +

μ P F(u) − f 2

2

.

It is easy to verify that Algorithm 1 can be applied to the above problem with minor modiﬁ2 cations to formula (2.4), including replacing F(K) by a section vector p ∈ Rn whose entries are 1 corresponding to rows selected by P and 0 otherwise. We have applied this approach to

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

270

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

solving simulation problems in compressive magnetic resonance (MR) imaging with promising results. Multiple regularization terms. Such examples arise in both image processing (e.g., [28, 8]) and compressive MR imaging (e.g., [22, 26, 29]). For example, the model (6.4) min u i

Di u + α Φu

1

+

μ P F(u) − f 2

2

,

where Φ is an orthonormal matrix, can be approximated by (6.5) min wi + i u,v,w

β wi − Di u 2

2

+α

v

1

+

γ v − Φu 2

2

+

μ P F(u) − f 2

2

.

In an alternating minimization setting, w can be solved by (2.2) as before, v can be solved by 1D shrinkage, and u can be solved by four FFTs (two for w1 and w2 , one for ΦT v noting that v − Φu = ΦT v − u , and one inverse FFT). In general, one needs to introduce an auxiliary variable for each regularization term and manages to obtain easily solvable subproblems. Finally, we mention that even in the case where the quadratic part in (1.3) cannot be minimized with respect to u via fast transforms as in formula (2.4), model (1.3) may still lead to an eﬃcient algorithm if the quadratic in u can be eﬀectively decreased by a few iterations of an iterative method such as a preconditioned CG method. This topic is of interest for further investigation. 7. Concluding remarks. We proposed, analyzed, and tested a new algorithm FTVd for solving the TV/L2 model (1.2). Our approximation model (1.3) was derived from the classic quadratic penalty method ﬁrst proposed by Courant in 1943, but also falls into the more recent half-quadratic class for image processing after our extensions. In fact, it can be regarded as a long missing, isotropic TV model in the half-quadratic class whose quadratic portion has a block circulant Hessian, thus allowing fast transforms to be utilized in quadratic minimization. We established strong convergence results for our algorithm and validated a continuation scheme. Our numerical results convincingly show that FTVd can be orders of magnitude faster than the LD algorithm. In particular, our results demonstrate that the TV/L2 model can be solved just as quickly as some other more widely used models, putting it on an equal footing with others in terms of solution speed. The results of this paper can be extended to a number of models involving TV regularization. More studies are underway to further explore the potential of the proposed approach. Acknowledgments. We are grateful to Paul Rodr´ iguez and Brendt Wohlberg at the Los Alamos National Lab for valuable discussions and for making their algorithm library NUMIPAD [36] available to us.

REFERENCES [1] R. Acar and C. R. Vogel, Analysis of total variation penalty methods, Inverse Problems, 10 (1994), pp. 1217–1229. [2] M. Allain, J. Idier, and Y. Goussard, On global and local convergence of half-quadratic algorithms, IEEE Trans. Image Process., 15 (2006), pp. 1130–1142.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

271

[3] P. Blomgren, T. F. Chan, and C. K. Wong, Total variation image restoration: Numerical methods and extensions, in Proceedings of the IEEE International Conference on Image Processing, Vol. 3, 1997, pp. 384–387. ` [4] E. Candes, J. Romberg, and T. Tao, Stable signal recovery from incomplete and inaccurate measurements, Comm. Pure Appl. Math., 59 (2006), pp. 1207–1223. [5] A. S. Carasso, Linear and nonlinear image deblurring: A documented study, SIAM J. Numer. Anal., 36 (1999), pp. 1659–1689. [6] A. Chambolle, An algorithm for total variation minimization and applications, J. Math. Imaging Vision, 20 (2004), pp. 89–97. [7] A. Chambolle and P. L. Lions, Image recovery via total variation minimization and related problems, Numer. Math., 76 (1997), pp. 167–188. [8] T. F. Chan, S. Esedoglu, F. Park, and A. Yip, Total variation image restoration: Overview and recent developments, in Handbook of Mathematical Models in Computer Vision, edited by N. Paragios, Y. Chen, and O. Faugeras, Springer-Verlag, New York, 2006, pp. 17–31. [9] T. F. Chan, G. H. Golub, and P. Mulet, A nonlinear primal-dual method for total variation-based image restoration, SIAM J. Sci. Comput., 20 (1999), pp. 1964–1977. [10] T. F. Chan and P. Mulet, On the convergence of the lagged diﬀusivity ﬁxed point method in total variation image restoration, SIAM J. Numer. Anal., 36 (1999), pp. 354–367. [11] T. F. Chan and J. Shen, Theory and Computation of Variational Image Deblurring, IMS Lecture Notes, Institute of Mathematical Statistics, Beachwood, OH, 2006. ´ [12] P. Charbonnier, L. Blanc-Feraud, G. Aubert, and M. Barlaud, Deterministic edge preserving regularization in computed imaging, IEEE Trans. Image Process., 6 (1997), pp. 298–311. [13] R. Courant, Variational methods for the solution of problems with equilibrium and vibration, Bull. Amer. Math. Soc., 49 (1943), pp. 1–23. [14] G. Demoment, Image reconstruction and restoration: Overview of common estimation structures and problems, IEEE Trans. Acoust. Speech Signal Process., 37 (1989), pp. 2024–2036. [15] D. C. Dobson and F. Santosa, Recovery of blocky images from noisy and blurred data, SIAM J. Appl. Math., 56 (1996), pp. 1181–1198. [16] D. Donoho, Compressed sensing, IEEE Trans. Inform. Theory, 52 (2006), pp. 1289–1306. [17] D. Geman and G. Reynolds, Constrained restoration and the recovery of discontinuities, IEEE Trans. Pattern Anal. Mach. Intell., 14 (1992), pp. 367–383. [18] D. Geman and C. Yang, Nonlinear image recovery with half-quadratic regularization, IEEE Trans. Image Process., 4 (1995), pp. 932–946. [19] D. Goldforb and W. Yin, Second-order cone programming methods for total variation-based image restoration, SIAM J. Sci. Comput., 27 (2005), pp. 622–645. [20] T. Goldstein and S. Osher, The Split Bregman Algorithm for L1 Regularized Problems, UCLA CAM Report 08-29, UCLA, Los Angeles, 2008. [21] R. Gonzalez and R. Woods, Digital Image Processing, Addison-Wesley, Reading, MA, 1992. [22] L. He, T.-C. Chang, S. Osher, T. Fang, and P. Speier, MR Image Reconstruction by Using the Iterative Reﬁnement Method and Nonlinear Inverse Scale Space Methods, UCLA CAM Report 06-35, UCLA, Los Angeles, 2006. [23] Y. Huang, M. K. Ng, and Y.-W. Wen, A fast total variation minimization method for image restoration, Multiscale Model. Simul., 7 (2008), pp. 774–795. [24] J. Idier, Convex half-quadratic criteria and interacting auxiliary variables for image restoration, IEEE Trans. Image Process., 10 (2001), pp. 1001–1009. [25] A. K. Katsaggelos, ed., Digital Image Restoration, Springer-Verlag, Berlin, 1991. [26] M. Lustig, D. Donoho, and J. Pauly, Sparse MRI: The application of compressed sensing for rapid MR imaging, Magn. Reson. Med., 58 (2007), pp. 1182–1195. [27] M. Lysaker, A. Lundervold, and X. C. Tai, Noise removal using fourth-order partial diﬀerential equations with applications to medical magnetic resonance images in space and time, IEEE Trans. Image Process., 12 (2002), pp. 1579–1590. [28] O. M. Lysaker and X.-C. Tai, Iterative image restoration combining total variation minimization and a second-order functional, Internat. J. Comput. Vision, 66 (2006), pp. 5–18.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

272

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

[29] S. Ma, W. Yin, Y. Zhang, and A. Chakraborty, Compressed MR imaging based on wavelets and total variation, in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’08), 2008. [30] A. Marquina and S. Osher, Explicit algorithms for a new time dependent model based on level set motion for nonlinear deblurring and noise removal, SIAM J. Sci. Comput., 22 (2000), pp. 387–405. [31] R. Neelamani, H. Choi, and R. G. Baraniuk, Wavelet-based deconvolution for ill-conditioned systems, in Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing Vol. 6, 1999, pp. 3241–3244. [32] R. Neelamani, H. Choi, and R. G. Baraniuk, ForWaRD: Fourier-wavelet regularized deconvolution for ill-conditioned systems, IEEE Trans. Signal Process., 52 (2004), pp. 418–433. [33] M. K. Ng, R. H. Chan, and W.-C. Tang, A fast algorithm for deblurring models with Neumann boundary conditions, SIAM J. Sci. Comput., 21 (1999), pp. 851–866. [34] M. Nikolova and M. K. Ng, Analysis of half-quadratic minimization methods for signal and image recovery, SIAM J. Sci. Comput., 27 (2005), pp. 937–966. [35] J. Nocedal and S. J. Wright, Numerical Optimization, Springer-Verlag, New York, 1999. [36] P. Rodr´ iguez and B. Wohlberg, Numerical Methods for Inverse Problems and Adaptive Decomposition (NUMIPAD), http://numipad.sourceforge.net/. [37] L. Rudin and S. Osher, Total variation based image restoration with free local constraints, in Proceedings of the 1st IEEE International Conference on Image Processing, Vol. 1, 1994, pp. 31–35. [38] L. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation based noise removal algorithms, Phys. D, (1992), pp. 259–268. [39] D. Strong and T. F. Chan, Edge-preserving and scale-dependent properties of total variation regularization, Inverse Problems, 19 (2003), pp. 165–187. [40] A. Tikhonov and V. Arsenin, Solution of Ill-Posed Problems, Winston, Washington, DC, 1977. [41] C. R. Vogel and M. E. Oman, Iterative methods for total variation denoising, SIAM J. Sci. Comput., 17 (1996), pp. 227–238. [42] Y. Wang, W. Yin, and Y. Zhang, A Fast Algorithm for Image Deblurring with Total Variation Regularization, CAAM Technical Report 07-10, Rice University, Houston, TX, 2007. [43] B. Wohlberg and P. Rodr´ iguez, An iteratively reweighted norm algorithm for minimization of total variation functionals, IEEE Signal Process. Lett., 14 (2007), pp. 948–951. [44] Y. L. You and M. Kaveh, Fourth-order partial diﬀerential equations for noise removal, IEEE Trans. Image Process., 9 (2000), pp. 1723–1730.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.…...

Premium Essay

...Since the dawn of time, men and women have aspired to have the best looking image. This aspiration is often measured by the appearances of others – the appearances that appeal to the majority. People strive to achieve this by giving in to advertisements which eloquently sell materials that will help them reach this goal. The media is the main contributor to self-image and how the majority thinks they should look. Today’s society has become too attached to this idea of self-image, and the media has become a crutch on people’s outlook on beauty. This single image of beauty created by the media comes with a great deal of problems for society. Mass media has been able to shape popular culture and influence public opinion, and when abused, the power of the media can harm the general population. Certain tactics and strategies have been used to create the picturesque image in today’s world. This image is seen in every sort of corporate material or entertainment source. Corporations have gone as far as to promote this body image – the supposed “perfect body” – in the toys they sell to children. For instance, in her article “Drugs, Sports, Body Image and G.I. Joe,” Natalie Angier looks at the “G.I. Joe [which has gotten] more muscular and sharply defined, or ‘cut,’ than the model before” (Shea 486). This doll which exploits the “perfect” body image is one of the most popular toys for male children. Magazines are another source that influences the people on what is and is not......

Words: 638 - Pages: 3

Free Essay

...INTRODUCTION The term image refer to two dimensional light intensity function f(x, y), where x and y denote spatial coordinates and the value of at any point (x, y) is proportional to the brightness (or gray level) of the image at that point. A digital image is an image f(x, y) that has been discretized both in spatial coordinate and brightness. A digital image can be considered a matrix whose row and column indices identify a point in the image and the corresponding matrix element value identifies the gray level at that point. The elements of such a digital array are called image elements, picture elements. The term digital image processing refers to processing of a two dimensional picture by a digital computer. An image given in the form of a transparency, slide, photograph or chart is first digitized and stored as matrix of binary digits in computer memory. This digitized image can be processed and/or displayed on a high resolution monitor. Segmentation subdivides an image into its constituent regions or objects .The level to which the segmentation is carried depends on the problem being solved. Segmentation is carried until we are able to distinguish the object from its backgrounds. Image segmentation is based on one of the two basic properties of intensity values: discontinuity and similarity. Thresholding is a similarity approach of segmenting an image. The simplest thresholding technique is to partition the......

Words: 560 - Pages: 3

Premium Essay

...Introduction Brand image is that the current outlook of the shoppers about a brand. It is characterized as a exclusive bundle of associations at gaps the minds of target customers. It signifies what the emblem presently stands for. it is a assembly of convictions order about a exact brand. In short, it is not anything but the consumers’ perception about the product. it is the kind inside which a particular emblem is positioned in the market. Emblem likeness expresses emotional worth and not just a mental representation. Emblem image is not anything but aide degree organization’s feature. it is aide degree accumulation of contact and fact by individuals external to an association. It ought to focus aide degree organization’s operation and dream to any or all. The major parts of positive emblem image are- exclusive emblem reflective organization’s likeness, saying recounting organization’s business in short and emblem identifier carrying the key values. Brand image is that the overall impression in consumers’ mind that's formed from all sources. shoppers develop various associations with the brand. supported these associations, they form brand image. an image is made about the brand on the idea of subjective perceptions of associations bundle that the shoppers have about the brand. Volvo is associated with safety. Toyota is associated with reliableness Image of a brand has a great impact on consumer’s preferences, because it’s the current view of a product according to......

Words: 1557 - Pages: 7

Free Essay

...IMAGE OF THE CITY KEVIN LYNCH Lynch influenced the field of city planning through his work on the theory of city form, and studies relating to human perceptions of the city y y , g p p y on the perception of the city environment and its consequences for city design. Lynch says "Looking at cities can give a special pleasure, however y y g g p p , commonplace the sight may be. Like a piece of architecture, the city is a construction in space, but of a vast scale, . . . perceived only in the course of long spans of time . . . At every instant, there is more than the eye can see, more than the ear can hear, a setting or view waiting to be explored. , , g g p Nothing is experienced by itself, but always in relation to its surroundings, the sequences of events leading up to it, the memory of past experiences . . . Every citizen has had long associations with some part of his city, and his image is soaked in memories and meanings . . . “ g g Theory of Kevin L h K i Lynch IMAGE OF THE CITY KEVIN LYNCH Image of the city 1. 2. 3. 4. 5. The Image of the Environment Three cities The city Image and Its elements City Form A new Scale Book contents Appendices 1. 1 Some references to orientation 2. The Use of the methodology 3. Two examples of analysis IMAGE OF THE CITY KEVIN LYNCH Methodology of working: Make visual plan Analyze the existing form and public image of the area. Understand the critical problems, opportunities and image elements and use......

Words: 1627 - Pages: 7

Free Essay

...Image of God Chart part 11 1 Grand Canyon University: HTH-359 July 24, 2013 Image of God Chart part 11 2 Terry Harris is an African American Black woman who resides in Kent, Ohio and appears to work hard for what she has. Terry’s temperament is sweet for she is a kind and gentle individual but, she struggles because of her crack, alcohol, and marijuana addictions. Although she attends Church on a regular basis she still needs a support system/group to help or rather show her how to get past her addictions. Angela (my wife) was her support system but, since we have relocated Terry has down spiraled. Angela picks up the telephone weekly to encourage her by having Bible study and reading scriptures with her. She also has invited her to come to stay with us as a vacation to get her away from that scene. Terry's is in her sixties and Angela told her "You would not look good wearing an orange jail jumpsuit, it that what you want because, if you do not stop you will end up in jail or even worse the graveyard." It is pretty much Terry's daughter Kim who constantly keeps bringing the drugs into her mother's home. Angela told Terry "Kim is grown and needs her own place so put her out." Spiritually she needs all the help she can get, her health is okay but she is also the type to get mad at others when they attempt to step in and her help. She will tell you "I don't get mad I get even." What she does not realize if she does not get help soon......

Words: 1128 - Pages: 5

Free Essay

...Analyzing Image Essay A picture is worth a thousand words, yet even some pictures of delectable items are not so appealing to the eye, and capturing the essence of food is quite the job. The portrayal of scrumptious or attractive décor has a price, and in some cases can be deceiving. Though there are not any words in these advertisements, a lot is certainly depicted in each. Not everything that glitters is always gold. The drive to have and consume the extravagant is compelling. Seen the first photo advertisement for Clementine Granita is an example of sweet delectable bites. This is certainly one of the better ads that appeal to the average foodie. The photographer used simple yet effective lighting and background, perhaps to instill simplicity. It conveys a measure of healthy due to the white background and clear glass Pyrex. Pure ingredients are silently conveyed in this advertisement for Clementine desserts. Clearly, the objectivity is to get the viewer to buy and eat this fruit. This one-sided framed shot suggests the audience too want what appears to be already eaten or is not present. The top illumination sings of bright, fresh, outdoor, color, which insists that one’s mouth salivates while reading the recipe. Consider though the other side too this advertisement, which is what they do not tell the consumer. Consider how much sugar is in this dish. The calories in each little beautiful morsel contain enough to put someone in a diabetic comma if consumed. The......

Words: 657 - Pages: 3

Free Essay

...Analysis of SVD based Image Fusion: A review xxxxxxxxxxx xxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxx ABSTRACT: In this paper, we have reviewed different types of Image fusion techniques based on Singular Value Decomposition (SVD) technique. Basically, Image fusion can be described as a technique which is used to generate a single good quality image from one or more images. Image fusion can be applied at many levels viz. pixel level, feature level, signal level and decision level. Image fusion can be applied in many areas like recognition of patterns, to enhance visual features, detection of objects, area surveillance etc.[2] The techniques that are used mostly are Intensity-Hue-Saturation (IHS), high pass filtering, principal component analysis (PCA), different arithmetic combinations, multi resolution analysis based methods (pyramid algorithm and wavelet transform), Artificial Neural Networks(ANN), Singular Value Decomposition (SVD) etc.[4] Nowadays, SVD is becoming very popular technique for image fusion due to many factors like conceptuality, stability and it is also a robust and reliable orthogonal decomposition technique. A huge advantage of SVD is that it can also adjust the variations that are present in the local statistics of an image[2]. In this paper, we have compared and reviewed different types of modifications that can be added to the basic SVD technique. Keywords: Singular Value Decomposition, Tensors, Image fusion. 1. INTRODUCTION ......

Words: 2109 - Pages: 9

Free Essay

...Images Durch die Entwicklung der Gesellschaft zu einer Mediengesellschaft, wird der Radius des eigenen Erlebens vom unmittelbaren Umfeld auf die Welt erweitert. Daraus ergibt sich, dass wir viele Erfahrungen nicht mehr selbst machen sondern wir sie von Medien zugänglich gemacht bekommen. Bsp.: Jeder weis über das angebliche Sexualverhalten eines amerikanischen Präsidenten mit seiner Sekretärin bescheid, hat ihn aber nie persönlich kennengelernt. Um die Flut an Informationen zu bewältigen, müssen sie gebündelt und vereinfacht werden. Die Mediengesellschaft kreiert deshalb eine Art Stellvertreter. Ein solcher Stellvertreter ist das Image Erste Definition Ein Image ist die „Gesamtheit der Voststellungen, Einstellungen und Gefühle, die eine Person im Hinblick auf ein Obkekt (Person, Organisation, Produkt, Idee) besitzt“ (Fuchs et al. 1978, S.330) 1. Eigenheiten von Images nur selten wahrhaftig -werden interessenbezogen u. zielgruppenorientiert konstruiert ---> schliesst Objektivität oder Wahrhaftigkeit aus -müssen nicht wahrhaftig sein sondern nur für zutreffend gehalten werden -Goffman: Image als Image einer Person -jeder Mensch ist bestrebt gegenüber der Öffentlichkeit einen guten Eindruck zu machen, also ein Image aufzubauen und dieses aufrecht zu erhalten Bsp: mehr Schein als Sein: man versucht beispielesweise selbstbewusst, cool und ruhig zu wirken... ist aber unsicher und aufgeregt --->das nach außen......

Words: 1191 - Pages: 5

Free Essay

...“Body Image” Have you see a fashion show on the television or read a magazine? If you have then you have been exposed to the media revealing models as beautiful. What you don’t know is that a decent percentage of these models are suffering from eating disorders. I trust the media is to blame for our country’s epidemic of eating disorders because, not only do magazines and television portray skinny to be in, but also songs in our nation deliver the attitude. The burdens to encounter the world’s demands to reaching self-satisfaction with one’s body image emotional have influenced the impression of eating disorders. As people are exposed to countless forms of media not only just on television or in a magazine but on the radio, internet & part of our everyday life media being a huge method of communication has been characterized as the leading reason of why people suffer from eating disorders. The role of media profoundly adds to the increase of abnormal eating behaviors within an individual. Without society and the media generating a fabricated image of attractiveness, the calculation of women suffering from eating disorders would decline extremely, women would discontinue trying to grasp a body weight that is unnatural also practically unmanageable to accomplish. Body image has been an ongoing issue in the world and eating disorders are frequently developed, especially in females, in an effort to keep their bodies “fit” and “in shape.” The physical appearance of a person...

Words: 2834 - Pages: 12

Premium Essay

...Image Cafe Summary: Clarence Wooten, Jr. is the founder and CEO of Image Cafe which the World's first online superstore of pre-fabricated websites for small businesses. When he was a child, his childhood dream was become a rich man. He did a lot of things about his hobbies and interests. He was fascinated video games, then he discovered the game cartridges were too expensive for him to buying. Therefore, he found the computers did not require the cartridges to play. He learned the computer games could be transferred. When he was 12, Wooten became so immersed the computer world. He parents banned him using computer. However, he became an athletic when his high school career but Wooten never lost interest in the computer. Both of his parents become self-employed when Wooten was a teenager. The income of his family was depending on the success of his parent own business but it was just small business. He moved to transfer to different school systems but eight times. Transitioning is hard for him. However, he quickly solved this problem. He want his life can be succeed. Wooten and his group members used the computers to crack anticopying features enabling electronic games to be duplicated. Wooten had become well known as part of competitive and cream of the crop computer underworld. He started his interest which is his love of video games. He became obsessed with acquiring more and more games quickly. The group was a team. They developed the software was distributed to pirates...

Words: 1512 - Pages: 7

Premium Essay

...Images Project Paper Week Three HR587 Managing Organizational Change Professor Maxine Walker Mercedric Golden Keller Graduate School of Management Devry University The organization I decided to do my change analysis research paper is my Army Reserve unit located in Grand Prairie, TX. I was assigned to the unit after coming off active duty with the Army in September of 2009. The unit is a battalion sized training unit with ninety percent of its members being male Soldiers. The battalion mission is conduct training readiness oversight and mobilization of designated active and reserve component forces in the western are of responsibility in order to provide trained and ready forces to regional combatant commanders. The battalion supports pre-mobilization training for reserve component forces in accordance with our Higher Headquarters, First Army, Division West located at Fort Hood, TX. Some of the specific tasks of the unit is to assess and report pre-mobilization readiness for reserve component forces; conduct mobilization and demobilization operations; conduct counter-improvised explosive device, counter insurgency and escalation of force training; provide command and control over assigned and mobilized forces; and provide operational force protection. Most of these training tasks and activities have traditionally been performed by all male Soldiers since it has long been considered a male’s job to perform any type of combat related duty or......

Words: 984 - Pages: 4

Free Essay

...1 Introduction to Image Processing and Compression Using Matlab/Octave Table of content Image Processing and Compression Using Matlab ........................................................................... 1 1.0 Basic image operation.............................................................................................................. 1 1.1 Image display ....................................................................................................................... 1 1.2: Color Image Display........................................................................................................... 2 1.3 Image Array Indexing .......................................................................................................... 3 1.4 Converting image array to image file .................................................................................. 8 1.5 Image Manipulation ............................................................................................................. 8 1.5.1 Resize image ..................................................................................................................... 9 1.6 Image Analysis Using Histogram ...................................................................................... 10 1.6.1 Generating feature vector for image using greyscale histogram .................................... 14 1.6.2 Comparing Image Descriptor ...............................................................

Words: 3272 - Pages: 14

Premium Essay

...and also, both share the similar style of being a mid ankle rise shoe. Image one is from the 1995 era, and image two is from the 1985 area. Comparing the two advertisements shows the effective type of advertising that Nike used in two different time periods. These advertisements for Nike’s basketball shoes are used to attract customers to buying their shoes by using athletes to endorse them, the images surroundings and background to give this you a personality, and the text on the advertisement to support the reasoning behind the shoe. In image one, we see the shoe up close with a spotlight shining down to enhance the detail with a solid black background. It is a jet black shoe with an air cushion on the sole. It has a combination of mesh and solid material that reflects the light. To the left of the shoe, the text contains both white and gray color. It says and large, bold text, "unbreakable. Lightweight explosiveness. The LeBron X” (Nike. The LeBron X.). The Nike logo also called the swoosh is above the text. From the type of swoosh, the advertisement is of the 1995 to present era. You can call this advertisement by the name of the shoe, LeBron advertisement. However, the second image is very different. This image is more focused on Michael Jordan as he is holding the shoe in his right hand at chest level. Jordan is pointing at the shoe and smiling while wearing his sunglasses. The text of this Image are not as large as the first, but they are noticeable and can attract......

Words: 1279 - Pages: 6

Free Essay

...Bradley Emery i think one cultural change that would be benifical to the mental health of Americans would be the image that society portrays as the "perfect female". i think that we as a society stress that attractive women need to be skinny and clear skin that it puts a lot of pressure on younger/ teenage girls to try to look like that. willing to do anything to be considered attractive or perfect young girls put their body through some situations that their bodies are not meant to handle such as anorexia, bulemia, and binge eating. In recent studies (Smolak, 2011) 40-60 percent of girls starting at the age of are worried about becoming overweight and continue to worry about it for the span of their life. The detrimental part about this is the fact that it is causing health problems to young innocent boys and girls that are under what they think is a microscope of their peers. There is so much pressure on the people that are worried about what they look like and worried about what others think about them that they turn to a horrible option in eating disorders. Not only do eating disorders cause muscle loss, reduction of bone density, lack of nutrients, and fainting but it also can result in death. A change in this i think would be more advertisments showing that you do not have to be super skinny to be considered attractive or to be a super model. commercials and magzines could put girls that are bigger than the normal cover girl or bigger than the thin super model......

Words: 503 - Pages: 3

Premium Essay

...management process in order to make the require changes happen. When an organization change its processes and working style for adding more value to the organization and its operations known as the change management initiative. The change management initiative can be done in various ways which depends on the manager of the organization. Each and every change management process has its own way of changing the processes and operations of an organization. The images of managing change are a set of six points which has the various ways of change management. Each and every image for managing changes has its own particular kind of pattern and formula of implementing the change which is different from one to another. There are six images for the managing change which are director, navigator, caretaker, coach and nurturer. These six images provide various style of implementing changes to the managers which may be different according to the nature and size of the organization. The managers need to choose the appropriate image out of six images which can be beneficial for their organization. Application Analysis The change management case study is related to the company Nike where thy company has made the major changes in its product range in order to reach maximum customers and to cover more market segments and customers. The NIKE was started through the production of sportswear which has changed now and the company is dealing in various other products as well. The change......

Words: 634 - Pages: 3