cam11-32 Minimization for Wavelet Frame Based Image Restoration

�0 Minimization for Wavelet Frame Based Image Restoration Yong Zhang ∗ Bin Dong † Zhaosong Lu‡ May 12, 2011 Abstract The theory of (tight) wavelet frames has been extensively studied in the past twenty years and they are currently widely used for image restoration and other image processing and analysis problems. The success of wavelet frame based models, including balanced approach [20, 8] and analysis based approach [13, 32, 49], is due to their capability of sparsely approximating piecewise smooth functions like images. Motivated by the balanced approach and analysis based approach, we shall propose a wavelet frame based �0 minimization model, where the �0 of the frame coefficients are penalized. We adapt the penalty decomposition (PD) method of [40] to solve the proposed optimization problem. Numerical results showed that the proposed model solved by the PD method can generate images with better quality than those obtained by either analysis based approach or balanced approach in terms of restoring sharp features as well as maintaining smoothness of the recovered images. Some convergence analysis of the PD method will also be provided. Key words: �0 minimization, wavelet frame, image restoration. 1 Introduction Mathematics has been playing an important role in the modern developments of image processing and analysis. Image restoration, including image denoising, deblurring, inpainting, tomography, etc., is one of the most important areas in image processing and analysis. Its major purpose is to enhance the quality of a given image that is corrupted in various ways during the process of imaging, acquisition and communication, and enable us to see crucial but subtle objects residing in the image. Therefore, image restoration is an important step to take towards accurate interpretations of the physical world and making optimal decisions. 1.1 Image Restoration Image restoration is often formulated as a linear inverse problem. For the simplicity of the notations, we denote the images as vectors in Rn with n equals to the total number of pixels. A typical image restoration problem is formulated as f = Au + η, (1.1) where f ∈ Rd is the observed image (or measurements), η denotes white Gaussian noise with variance σ2, and A ∈ Rd×n is some linear operator. The objective is to find the unknown true image u ∈ Rn from the observed image f . Typically, the linear operator in (1.1) is a convolution operator for image deconvolution problems, a projection operator for image inpainting and partial Radon transform for computed tomography. ∗Department of Mathematics, Simon Fraser University, Burnaby, BC, V5A 1S6, Canada. (Email: yza30@sfu.ca). †Department of Mathematics, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA. (Email: dongbin725@gmail.com) ‡Department of Mathematics, Simon Fraser University, Burnaby, BC, V5A 1S6, Canada. (Email: zhaosong@sfu.ca). This author was supported in part by NSERC Discovery Grant. 1 To solve u from (1.1), one of the most natural choices is the following least square problem min u∈Rn ‖Au− f‖22, where ‖ · ‖2 denotes the �2-norm. This is, however, not a good idea in general. Taking image deconvolution problem as an example, since the matrix A is ill-conditioned, the noise η possessed by f will be amplified after solving the above least squares problem. Therefore, in order to suppress the effect of noise and also preserve key features of the image, e.g., edges, various regularization based optimization models were proposed in the literature. Among all regularization based models for image restoration, variational methods and wavelet frames based approaches are widely adopted and have been proven successful. The trend of variational methods and partial differential equation (PDE) based image processing started with the refined Rudin-Osher-Fatemi (ROF) model [44] which penalizes the total variation (TV) of u. Many of the current PDE based methods for image denoising and decomposition utilize TV regularization for its beneficial edge preserving property (see e.g., [41, 45, 42]). The ROF model is especially effective on restoring images that are piecewise constant, e.g., binary images. Other types of variational models were also proposed after the ROF model. We refer the interested readers to [41, 42, 2, 24] and the references therein for more details. Wavelet frame based approaches are relatively new and came from a different path. The basic idea for wavelet frame based approaches is that images can be sparsely approximated by properly designed wavelet frames, and hence, the regularization used for wavelet frame based models is the �1-norm of frame coefficients. Although wavelet frame based approaches take similar forms as variational methods, they were generally considered as different approaches than variational methods because, among many other reasons, wavelet frame based approaches is defined for discrete data, while variational methods assume all variables are functions. This impression was changed by the recent paper [11], where the authors established a rigorous connection between one of the wavelet frame based approaches, namely the analysis based approach, and variational models. It was shown in [11] that the analysis based approach can be regarded as a finite difference approximation of a certain type of general variational model, and such approximation will be exact when image resolution goes to infinity. Furthermore, the solutions of the analysis based approach also approximate, in some proper sense, the solutions of corresponding variational model. Such connections not only grant geometric interpretation to wavelet frame based approaches, but also lead to even wider applications of them, e.g., image segmentation [28] and 3D surface reconstruction from unorganized point sets [30]. On the other hand, the discretization provided by wavelet frames was shown, in e.g., [20, 22, 12, 13, 11, 29], to be superior than the standard discretizations for some of the variational models, due to the multiresolution structure and redundancy of wavelet frames which enable wavelet frame based models to adaptively choose a proper differential operators in different regions of a given image according to the order of the singularity of the underlying solutions. For these reasons, as well as the fact that digital images are always discrete, we use wavelet frames as the tool for image restoration in this paper. 1.2 Wavelet Frame Based Approaches We now briefly introduce the concept of tight frames and tight wavelet frame, and then recall some of the frame based image restoration models. Interesting readers should consult [43, 25, 26] for theories of frames and wavelet frames, [46] for a short survey on theory and applications of frames, and [29] for a more detailed survey. A countable set X ⊂ L2(R) is called a tight frame of L2(R) if f = ∑ h∈X 〈f, h〉h ∀f ∈ L2(R), where 〈·, ·〉 is the inner product of L2(R). The tight frame X is called a tight wavelet frame if the elements of X are generated by dilations and translations of finitely many functions called framelets. The construction of framelets can be obtained by the unitary extension principle (UEP) of [43]. In our implementations, we will use the piecewise linear B-spline framelets constructed by [43]. Given a 1-dimensional framelet system 2 for L2(R), the s-dimensional tight wavelet frame system for L2(Rs) can be easily constructed by using tensor products of 1-dimensional framelets (see e.g., [25, 29]). In the discrete setting, we will use W ∈ Rm×n with m ≥ n to denote fast tensor product framelet decomposition and use W� to denote the fast reconstruction. Then by the unitary extension principle [43], we have W�W = I, i.e., u = W�Wu for any image u. We will further denote an L-level framelet decomposition of u as Wu = (. . . ,Wl,ju, . . .) � for 0 ≤ l ≤ L− 1, j ∈ I, where I denotes the index set of all framelet bands and Wl,ju ∈ Rn. Under such notation, we have m = L× |I| × n. We will also use α ∈ Rm to denote the frame coefficients, i.e., α = Wu, where α = (. . . , αl,j , . . .) � , with αl,j = Wl,ju. More details on discrete algorithms of framelet transforms can be found in [29]. Since tight wavelet frame systems are redundant systems (i.e., m > n), the mapping from the image u to its coefficients is not one-to-one, i.e., the representation of u in the frame domain is not unique. Therefore, there are mainly three formulations utilizing the sparseness of the frame coefficients, namely, analysis based approach, synthesis based approach, and balanced approach. Detailed and integrated descriptions of these three methods can be found in [29]. The wavelet frame based image processing started from [20, 21] for high-resolution image reconstructions, where the proposed algorithm was later analyzed in [8]. These work lead to the following balanced approach [10, 14] min α∈Rm 1 2 ‖AW�α− f‖2D + κ 2 ‖(I −WW�)α‖22 + ∥∥∥∥∥∥∥ L−1∑ l=0 ⎛ ⎝∑ j∈I λl,j |αl,j |p ⎞ ⎠ 1/p ∥∥∥∥∥∥∥ 1 , (1.2) where p = 1 or 2, 0 ≤ κ ≤ ∞, λl,j ≥ 0 is a scalar parameter, and ‖ · ‖D denotes the weighted l2-norm with D positive definite. This formulation is referred to as the balanced approach because it balances the sparsity of the frame coefficient and the smoothness of the image. The balanced approach (1.2) was applied to various applications in [7, 9, 19, 23, 47, 39]. When κ = 0, only the sparsity of the frame coefficient is penalized. This is called the synthesis based approach, as the image is synthesized by the sparsest coefficient (see e.g., [27, 33, 34, 35, 36]). When κ = +∞, only the sparsity of canonical wavelet frame coefficients, which corresponds to the smoothness of the underlying image, is penalized. For this case, problem (1.2) can be rewritten as min u∈Rn 1 2 ‖Au− f‖2D + ∥∥∥∥∥∥∥ L−1∑ l=0 ⎛ ⎝∑ j∈I λl,j |Wl,ju|p ⎞ ⎠ 1/p ∥∥∥∥∥∥∥ 1 . (1.3) This is called the analysis based approach, as the coefficient is in range of the analysis operator (see, for example, [13, 32, 49]). Note that if we take p = 1 for the last term of (1.2) and (1.3), it is known as the anisotropic �1-norm of the frame coefficients, which is the case used for earlier frame based image restoration models. The case p = 2, called isotropic �1-norm of the frame coefficients, was proposed in [11] and was shown to be superior than anisotropic �1-norm. Therefore, we will choose p = 2 for our simulations. 1.3 Motivations and Contributions For most of the variational models and wavelet frame based approaches, the choice of norm for the regularization term is the �1-norm. Taking wavelet frame based approaches for example, the attempt of minimizing the �1-norm of the frame coefficients is to increase their sparsity, which is the right thing to do since piecewise smooth functions like images can be sparsely approximated by tight wavelet frames. Although the �1-norm of a vector does not directly correspond to its cardinality in contrast to �0-norm, it can be 3 regarded as a convex approximation to �0-norm. Such approximation is also an excellent approximation for many cases. It was shown by [15], which generalizes the exciting results of compressed sensing [16, 18, 17, 31], that for a given wavelet frame, if the operator A satisfies certain conditions, and if the unknown true image can be sparsely approximated by the given wavelet frame, one can robustly recover the unknown image by penalizing the �1-norm of the frame coefficients. For image restoration, however, the conditions on A as required by [15] are not generally satisfied, which means penalizing �0 “norm” and �1 norm may produce different solutions. Although both the balanced approach (1.2) and analysis based approach (1.3) can generate restored images with very high quality, one natural question is whether using �0 “norm” instead of �1-norm can further improve the results. On the other hand, it was observed, in e.g., [29] (also see Figure 2 and Figure 3), that balanced approach (1.2) generally generates images with sharper features like edges than the analysis based approach (1.3), because balanced approach emphasizes more on the sparsity of the frame coefficients. However, the recov- ered images from balanced approach usually contains more artifact (e.g., oscillations) than analysis based approach, because the regularization term of the analysis based approach has a direct link to the regularity of u (as proven by [11]) comparing to balanced approach. Although such trade-off can be controlled by the parameter κ in the balanced approach (1.2), it is not very easy to do in practice. Furthermore, when a large κ is chosen, some of the numerical algorithms solving (1.2) will converge slower than choosing a smaller κ (see e.g., [47, 29]). Since penalizing �1-norm of Wu ensures smoothness while not as much sparsity as balanced approach, we propose to penalize �0 “norm” of Wu instead. Intuitively, this should provide us a balance between sharpness of the features and smoothness for the recovered images. The difficulty here is that �0 minimization problems are generally hard to solve. Recently, penalty decomposition (PD) methods were proposed by [40] for a general �0 minimization problem that can be used to solve our proposed model due to its generality. The computational results for compressed sensing problems, sparse logistic regression and sparse inverse covariance selection problems demonstrated that their methods generally outperform the existing methods in terms of quality of solutions and/or computational efficiency. This motivates us to adapt one of their PD methods to solve our proposed �0 minimizations. The PD method uses the block coordinate descent (BCD) method to solve each penalty subproblem. In this paper, we will prove that any fixed point of the BCD method is a local minimizer of the corresponding subproblem. Under certain conditions, we further show that the sequence generated by BCD method converges and the limit is a local minimizer to the corresponding subproblem. These convergence analysis was missing from [40]. Furthermore, our numerical simulations show that the PD method generates recovered images with better quality than the results of both balanced approach (1.2) and analysis based approach (1.3). We now leave the details of the model and algorithm to Section 2 and details of simulations to Section 3. 2 Model and Algorithm We start by introducing some simple notations. The space of symmetric n× n matrices will be denoted by Sn. If X ∈ Sn is positive definite, we write X � 0. We denote by I the identity matrix, whose dimension should be clear from the context. Given an index set J ⊆ {1, . . . , n}, xJ denotes the sub-vector formed by the entries of x indexed by J . For any real vector, ‖ · ‖0 and ‖ · ‖2 denote the cardinality (i.e., the number of nonzero entries) and the Euclidean norm of the vector, respectively. In addition, ‖x‖D denotes the weighted l2-norm defined by ‖x‖D = √ x�Dx with D � 0. 2.1 Model We now propose the following optimization model for image restoration, min u∈Y,α=Wu 1 2 ‖Au− f‖2D + ∑ i λi‖αi‖0, (2.1) where Y is some convex subset of Rn. Here we are using the multi-index i and denote αi (similarly for λi) the value of α at a given pixel location within a certain level and band of wavelet frame transform. Comparing 4 to the analysis based model, we are now penalizing the number of nonzero elements of Wu. As mentioned earlier that if we emphasize too much on the sparsity of the frame coefficients as in the balanced approach or synthesis based approach, the recovered image will contain artifacts, although features like edges will be sharp; if we emphasize too much on the regularity of u like in analysis based approach, features in the recovered images will be slightly blurred, although artifacts and noise will be nicely suppressed. Therefore, by penalizing the �0 “norm” of Wu as in (2.1), we can indeed achieve a better balance between sharpness of features and smoothness of the recovered images. Given that the �0 “norm” is an integer-valued, discontinuous and nonconvex function, problem (2.1) is generally hard to solve. Some algorithms proposed in the literature, e.g., iterative hard thresholding algorithms [5, 6, 38], cannot be directly applied to the proposed model (2.1) unless W = I. Recently, Lu and Zhang [40] proposed a penalty decomposition (PD) method to solve the following general �0 minimization problem: min x∈X f(x) + ν‖xJ‖0 (2.2) for some ν > 0 controlling the sparsity of the solution, where X is a closed convex set in Rn, f : Rn → R is a continuously differentiable function, and ‖xJ‖0 denotes the cardinality of the subvector formed by the entries of x indexed by J . In view of [40], we adapt their PD method to tackle problem (2.1) directly in this paper. In addition, we apply the non-monotone gradient projection method proposed in [4] to solve one of the subproblem in the PD method. 2.2 Algorithm for Problem (2.1) In this section, we discuss how the PD method proposed in [40] solving (2.2) can be adapted to solve problem (2.1). Letting x = (u, α), J = {n + 1, . . . , n + m}, J¯ = {1, . . . , n}, f(x) = 12‖AxJ¯ − f‖2D andX = {x ∈ Rn+m : xJ = WxJ¯ and xJ¯ ∈ Y}, we can clearly see that the problem (2.1) takes the same form as (2.2). In addition, there obviously exists a feasible point (ufeas, αfeas) for problem (2.1) when Y �= ∅. Thus, the PD method studied in [40] can now be used to solve (2.1). We now discuss the implementation details of the PD method when solving the proposed wavelet frame based model (2.1). Given a penalty parameter > 0, the associated quadratic penalty function for (2.1) is defined as p�(u, α) := 1 2 ‖Au− f‖2D + ∑ i λi‖αi‖0 + 2‖Wu− α‖ 2 2. (2.3) Then we have the following PD method for problem (2.1) where each penalty subproblem is approximately solved by a block coordinate descent (BCD) method (see [40] for details). Penalty Decomposition (PD) Method for (2.1): Let 0 > 0, δ > 1 be given. Choose an arbitrary α0,0 ∈ Rm and a constant Υ such that Υ ≥ max{ 12‖Aufeas− f‖2D + ∑ i λi‖αfeasi ‖0,minu∈Y p�0(u, α0,0)}. Set k = 0. 1) Set q = 0 and apply the BCD method to find an approximate solution (uk, αk) ∈ Y × Rm for the penalty subproblem min{p�k(u, α) : u ∈ Y, α ∈ Rm} (2.4) by performing steps 1a)-1c): 1a) Solve uk,q+1 ∈ Argmin u∈Y p�k(u, α k,q). 1b) Solve αk,q+1 ∈ Arg min α∈Rn p�k(u k,q+1, α). 1c) Set (uk, αk) := (uk,q+1, αk,q+1). 2) Set k+1 := δ k. 3) If min u∈Y p�k+1(u, α k) > Υ, set αk+1,0 := ufeas. Otherwise, set αk+1,0 := αk. 5 4) Set k ← k + 1 and go to step 1). end Remark. In the practical implementation, we terminate the BCD method based on the relative progress of p�k(u k,q, αk,q) which can be described as follows: |p�k(uk,q, αk,q)− p�k(uk,q+1, αk,q+1)| max(|p�k(uk,q+1, αk,q+1)|, 1) ≤ �I . (2.5) Moreover, we terminate the outer iterations of the PD method once ‖Wuk − αk‖2 max(|p�k(uk, αk)|, 1) ≤ �O. (2.6) Next we discuss how to solve two subproblems arising in step a) and b) of the BCD method. 2.2.1 The BCD subproblem in step 1a) The BCD subproblem in step 1a) is in the form of min u∈Y 1 2 〈u,Qu〉 − 〈c, u〉 (2.7) for some Q � 0 and c ∈ Rn. Obviously, when Y = Rn, problem (2.7) is an unconstrained quadratic programming problem that can be solved by the conjugate gradient method. Nevertheless, the pixel values of an image are usually bounded. For example, the pixel values of a CT image should be always greater than or equal to zero and the pixel values of a grayscale image is between [0, 255]. Then the corresponding Y of these two examples are Y = {x ∈ Rn : xi ≥ lb ∀i = 1, . . . , n} with lb = 0 and Y = {x ∈ Rn : lb ≤ xi ≤ ub ∀i = 1, . . . , n} with lb = 0 and ub = 255. To solve these types of the constrained quadratic programming problems, we apply the non-monotone projected gradient method proposed in [4] and terminate it using the duality gap and dual feasibility conditions (if necessary). For Y = {x ∈ Rn : xi ≥ lb ∀i =