彩色图像分割

Color image decomposition and restoration Jean-Franc¸ois Aujol a and Sung Ha Kang b,1,2 aDepartment of Mathematics, UCLA, 405 Hilgard Avenue, Los Angeles, California 90095-1555 bDepartment of Mathematics, 715 Patterson Office Tower, University of Kentucky, Lexington, KY 40506 Abstract Y. Meyer has recently introduced an image decomposition model to split an image into two components: a geometrical component and a texture (oscillatory) compo- nent. Inspired by his work, numerical models have been developed to carry out the decomposition of gray scale images. In this paper, we propose a decomposition al- gorithm for color images. We introduce a generalization of Meyer’s G norm to RGB vectorial color images, and use Chromaticity and Brightness color model with total variation minimization. We illustrate our approach with numerical examples. Key words: Total variation, structure, texture, color, image decomposition, image restoration. 1 Introduction In [1], Meyer introduced an image decomposition model based on the Rudin- Osher-Fatemi’s total variation minimization (TV) model [2]. A given image f is separated into f = u + v by minimizing the following functional, inf (u,v)∈BV×G/f=u+v ∫ Ω |∇u|+ α‖v‖G . (1) 1 We acknowledge supported by grants from the NSF under contracts DMS- 0312223, DMS-9973341, ACI-0072112, INT-0072863, the ONR under contract N00014-03-1-0888, the NIH under contract P20 MH65166, and the NIH Roadmap Initiative for Bioinformatics and Computational Biology U54 RR021813 funded by the NCRR, NCBC, and NIGMS. 2 Email :aujol@math.ucla.edu (Aujol), skang@ms.uky.edu (Kang) Preprint submitted to Elsevier Science 7 December 2004 The first term is a TV minimization which reduces u as the bounded variation (BV) component of the original image f . It is well known that BV is well suited to model the structure component of f [2], which means the edges of f are in the BV component u. The second term gives the v component containing the oscillatory part of the image, which is textures and noise. The Banach space G contains such signals with large oscillations. A distribution v belongs to G if v can be written as v = ∂1g1 + ∂2g2 = Div(g) g1, g2 ∈ L ∞. The G-norm ‖v‖G in (1) is defined as the infimum of all ‖g‖L∞ = supx∈Ω |g(x)|, where v = Div(g) and |g(x)| = √ |g1|2 + |g2|2(x). Any function belonging to the space G can present strong oscillations, nonetheless have a small norm [1,3]. Meyer’s model (1) was first successfully implemented by Vese and Osher [4,5], by considering the space Gp(Ω) = {v = Div(g) | g1, g2 ∈ L p(Ω)} with ‖v‖Gp = inf ‖ √ g21 + g 2 2‖p. The authors minimized the energy with respect to u, g1 and g2, and they solved the associated Euler-Lagrange equations. In [6], the authors used the H−1 norm instead of Meyer’s G norm. A different approach has been proposed in [7,8] by minimizing a convex functional which depends on the two variables u and v: inf u∈BV, ‖v‖G≤µ ∫ Ω |∇u|+ 1 2λ ‖f − u− v‖2L2 . (2) In this paper, we follow the image decomposition model of [7]. We review this method in section 2.1. More literature on image decomposition models can be found in [3,9–13]. These decomposition models are mostly devoted to gray- scale images. In this paper, we propose a decomposition algorithm for color images. There are various ways to deal with color images. For example, it can be treated as 3-dimensional vectorial functions [14], as tensor products of dif- ferent color components such as Chromaticity and Brightness (CB) or HSV nonlinear color models. Many related literature on color image models can be found at [15–21]. In this paper, we use 3D vectorial TV model [14], as well as Chromaticity and Brightness model. In [15], the authors showed that using Chromaticity and Brightness (CB) model gives a better color control and detail recovery for color image denoising, compared to channel by channel denoising or vectorial denoising. This CB model also provides a better color recovery compared to denoising HSV color system separately. Typically color images are represented by RGB (Red, Green, and Blue) color system, u : Ω → R3+ = {(r, g, b) : r, g, b > 0}. In the CB model, u is separated into the Brightness component ub = ‖u‖, and 2 the Chromaticity component uc = u/‖u‖ = u/ub. The Brightness component ub can be treated as a gray-scale image, while the Chromaticity component uc stores the color information which takes values on the unit sphere S2. The main contribution of the paper is to propose a color image decomposition model using TV minimization for color images. This paper is organized as fol- lows. In section 2, we generalize Meyer’s definition to color textures through the theory of convex analysis, following the review of [7] in section 2.1. We introduce a functional whose minimizers correspond to the color image de- composition, and we conclude this section by presenting some mathematical results. In section 3, we illustrate the details of numerical computations, and present numerical examples. In section 4, we conclude the paper with some final remarks. 2 Color Decomposition Model Meyer has introduced the G norm to capture textures in a noise free image [1]. The idea is that, while BV is a good space to model piecewise constant images, a space close to the dual of BV is well suited to model oscillating patterns. However, the dual of BV is not a separable space, and in [7], the authors considered the polar semi-norm associated to the total variation semi-norm for such purpose. We first review this approach. 2.1 Review of an image decomposition model for gray-scale images In [7], the authors introduced the following image decomposition model: inf u∈BV, ‖v‖G≤µ { J1D(u) + 1 2λ ‖f − u− v‖2L2 } , where J1D(u) = ∫ Ω |∇u| is the 1D total variation of u. The parameter λ controls the L2-norm of the residual f − u − v. The smaller λ is, the smaller the L2-norm of the residual gets. The bound µ controls the G norm of the oscillating component v. It is shown in [7] that solving (2) is equivalent to computing Meyer’s image decomposition (1). Let us denote by: µBG = {v ∈ G such that ‖v‖G ≤ µ} . We recall that the Legendre-Fenchel transform of F is given by F ∗(v) = supu (〈u, v〉L2 − F (u)), where 〈., .〉L2 stands for the L 2 inner product [22,23]. 3 It is shown in [7] that: J∗1D ( v µ ) = χµBG(v) =   0 if v ∈ µBG +∞ otherwise . Then, one can rewrite problem (2) above as: inf (u,v) { J1D(u) + J ∗ 1D ( v µ ) + 1 2λ ‖f − u− v‖2L2 } . (3) This functional (3) can be minimized with respect to the two variables u and v alternatively. First, fix u and solve for v which is the solution of infv∈µBG ‖f − u − v‖ 2, then fix v and solve for u which is the solution of infu ( J1D(u) + 1 2λ ‖f − u− v‖2 ) . In [7], the authors use Chambolle’s projec- tion algorithm [22] to compute the solution of each minimization problem. 2.2 Proposed model for color images Let us first define J as the total variation for 3D vector: J(u) = ∫ Ω √ |∇ur|2 + |∇ug|2 + |∇ub|2 , where r, g, and b stands for RGB channels. Let us denote by J ∗ the Legendre- Fenchel transform of J [24]. Then, since J is 1-homogeneous (that is J(λu) = λJ(u) for every u and λ > 0), it is a standard fact in convex analysis [22,23] that J∗ is the indicator function of a closed convex set K. We have: J∗(v) = χK(v) =   0 if v ∈ K +∞ otherwise . (4) We define the ~G norm by setting: ‖v‖ ~G = inf {µ > 0 | v ∈ µK} = inf { µ > 0 ∣∣∣ J∗ ( v µ ) = 0 } . We use ~G notation for 3-dimensional G norm. Notice that in one dimensional, this definition is exactly the same as Meyer’s original G norm [1], and this new definition is a natural extension of Meyer’s to the color case. Here, K is quite a complicated set; however, the simplest characterization is that K = {v | J∗(v) = 0}. We now propose a functional to split a color image f into a bounded variation 4 component u and a texture component v: inf u+v=f {J(u) + α‖v‖ ~G} . (5) In order to derive a partical numerical scheme, we slightly modify this func- tional by adding a L2 residual: inf u,v { J(u) + 1 2λ ‖f − u− v‖2 + J∗ ( v µ )} . (6) We show the details of the relation between equation (5) and equation (6) in subsection 2.3. Here, λ is to be small so that the residual f − u − v is negligeable, and µ controls the ‖.‖ ~G norm of v. Since the functional (6) is convex, a natural way to handle the problem is to minimize the functional with respect to each of the variable u and v alternatively, i.e. • First v being fixed, we search for u as a solution of: inf u ( J(u) + 1 2λ ‖f − u− v‖2 ) . (7) • Then, u being fixed, we search for v as a solution of: inf v∈µK ‖f − u− v‖2. (8) To solve these two minimization problems, we use the dual approach to the one used in [7]: we consider the direct total variation minimization approach. This will allow us to use total variation minimization algorithms devoted to color images [15]. We cannot use the numerical approach of [7], since the projection algorithm of [22] only works for gray-scale images. 2.3 Mathematical analysis of our color image decomposition model In this section, we only consider the discrete setting (for the sake of clarity), and present some mathematical results of the proposed functional. First, in the following proposition, we show that (8) can be solved by using direct TV minimization. Proposition 1 v˜ is the solution of (8), if and only if, w˜ = f − u − v˜ is the solution of: inf w ( J(w) + 1 2µ ‖f − u− w‖2 ) . (9) 5 Proof : This is a classical convex analysis result [3,22]. Let’s denote ∂H as the subdifferential of H (see [23,24]), and we recall that: w ∈ ∂H(u) ⇐⇒ H(v) ≥ H(u) + 〈w, v − u〉L2 , for all v in L 2. Here 〈., .〉L2 stands for the L 2 inner product. First, we remark that v˜ is the solution of (8) if and only if it minimizes: inf v { ‖f − u− v‖2 + J∗ ( v µ )} , (10) since J∗ is defined by (4). Let v˜ be the solution of (10). Then, as in [22,23], this is equivalent to v˜ + u − f ∈ ∂J∗ ( v˜ µ ) , which means v˜ µ ∈ ∂J (v˜ + u− f). Since w˜ is defined as w˜ = f − u− v˜, we get 0 ∈ ∂J (w˜) + 1 µ (−f + u + w˜), i.e. w˜ is a solution of (8). � We have just shown that equation (8) can be solved by direct computation of TV minimization (9). All the following lemma and propositions can be proved by straightforward generalization of results in [7]; therefore, we refer readers to [7] for more details, and we omit the proofs in this paper. Lemma 2 Problem (6) admits a unique solution (uˆ, vˆ). Outline of the proof : The existence of a solution comes from the convexity and the coercivity of the functional [24]. For the uniqueness, we first remark that (6) is strictly convex on BV ×µK, except in the direction (u,−u) . Then, with simple computations, it can be shown that if (uˆ, vˆ) is a minimizer of (6), then, for t 6= 0, (uˆ + tuˆ, vˆ − tuˆ) is not a minimizer of (6). � From lemma 2, we know that problem (6) has a unique solution. To compute it, we consider alternatively equations (7) and (8). This means that we consider the following sequence (un, vn): we set u0 = v0 = 0, define un+1 as the solution of infu ( J(u) + 1 2λ ‖f − u− vn‖ 2 ) , and vn+1 as the solution of infv∈µK ‖f − un+1 − v‖ 2. As a consequence of lemma 2, we get the convergence of (un, vn) to the unique solution (uˆ, vˆ) of (6). Proposition 3 The sequence (un, vn) converges to (uˆ, vˆ), the unique solution of problem (6), when n → +∞. From this result, we see that solving iteratively (7) and (8) amounts to solving (6). This justify the algorithm that we will propose in section 3.1. In section 2.2, we claimed that solving (6) is a way to solve (5). To explain this equivalence, we first introduce the following problem: inf u+v=f (J(u) + J∗(v/µ)) . (11) 6 The next result states the link between (5) and (11). Proposition 4 For a fixed α > 0, let (uˆ, vˆ) be a solution of problem (5). Then, if µ = ‖vˆ‖ ~G in (11), • (uˆ, vˆ) is also a solution of problem (11). • Conversely, any solution (u˜, v˜) of (11) (with µ = ‖vˆ‖ ~G) is a solution of (5). This proposition says that (5) and (11) are equivalent. To close the link be- tween (5) and (6), we check what happens when λ goes to zero in problem (6). This is explained in the followong result. Proposition 5 For a fixed α > 0 in (5), let α = ‖vˆ‖ ~G in (6) and (11). Let (uλn , vλn) be the solution of problem (6) with λ = λn. Then, when λn goes to 0, any cluster point of (uλn , vλn) is a solution of problem (11). All these results show that solving(7) and (8) iteratively is a way to solve problem (5): this is the theoretical justification of the decomposition algorithm that we will propose in section 3.1. In the following section, we detail the algorithm we use, and we show numerical examples to illustrate its efficiency. 3 Numerical Experiments For numerical computation for color image decomposition (5), we minimize the two functionals (7) and (9) alternatively. For the minimization, we compute the associated Euler-Lagrange equations. Notice that the two functionals (7) and (9) are almost exactly the same as the classical TV minimizing functional. Therefore, we can utilize all the benefits of well-known TV minimization tech- niques. Since we deal with color images, we use the results in [15]. In [15], the authors showed that Chromaticity and Brightness (CB) model gives the best denoising results, compared to denosing RGB channel by channel separately, denoising HSV channel by channel separately, or even vectorial color TV [14]. For solving (7), we use the CB model for the BV component u, and we use color TV for the texture component v by solving equation (9). This u component is the 3D RGB vector which is the structure component of the image. For u, we separate it into Chromaticity uc and Brightness ub components, and we denoise them with TV minimization separately, i.e. u = uc×ub. We believe this is the best way to keep the edges sharp, and get a good BV component of the image, as in [15]. As for the texture component v, we keep it as 3D RGB color vector and use color TV for denoising [14]. We could use the CB model for this v component as well; however, we kept it as one vector for the following three reasons. First of all, if color texture is one component it is easier to control. 7 Since the BV part u is already well kept by the CB model, the rest will also represents the texture well, and having one vector for v will be good enough and easy to handle. Secondly, by using color TV instead of the CB model, only one iteration is needed unlike two iterations in the CB model. Thirdly, which is the main reason for our choice, it is better to introduce some relations (coupled information) between the Chromaticity and Brightness components of u, and this is precisely what we do in the following algorithm. 3.1 Algorithm (1) Initially we set, f = fo, u = fo and v = 0 (fo is the original given image). (2) iterate m times: (a) Separate u to Chromaticity uc and Brightness ub components. (Also, separate f = fc × fb and v = vb × vb to chromaticity and brightness component respectively.) (b) For the Chromaticity component uc, solve the Euler-Lagrange equa- tion of (7) with vc and fc, and iterate n times: un+1 − un ∆t = ∇ · ∇un |∇un| + 1 λ (fc − un − vc). (12) (c) For the Brightness component ub, solve the 1D version of (12) with vb and fb, and iterate n times. (d) With updated uc and ub, let new u = uc × ub and w = f − u − v. Solve the Euler-Lagrange equation of (9) for w, and iterate n times; wn+1 − wn ∆t = ∇ · ∇wn |∇wn| + 1 µ (f − u− wn). (e) Update v = f − u− w. (3) Stopping test: we stop if max(|un+1 − un|, |vn+1 − vn|) ≤ � This algorithm decomposes a color image into f = u + v, where u is the structure component of the image, and v is the color texture component of the image. As a numerical computation, we used digital TV filter [25] type computation for non-flat TV denoising [15,26] as well as color TV [14]. For numerical computational details, we refer readers to [15,25]. 3.2 Numerical Experiments For typical experiments, image intensity was between 0 and 1. We used total iteration m = 5 and subiteration n = 30. For λ, for chromaticity we used 8 (a) f (b) u (c) v Fig. 1. (a) original image f , (b) BV component u, (c) texture v component (v + 0.5 plotted). Some of of the thicker branches are in the BV part u, while the thin and narrow branches in the bottom middle are in the v component. u as well as v are both color images. (a) f (b) u (c) v Fig. 2. (a) original image f , (b) BV component u, (c) texture v component (v + 0.5 plotted). All the details of image are in v, while the BV component is well kept in u. λc = 0.04, for brightness λb = 0.01, and we used µ = 0.1. The first two figures, Fig. 1 and Fig. 2, are image decomposition examples. From the original image f in (a), f is separated into two components u in image (b) and v in image (c). The texture part of the image, v clearly shows the color texture and the details of the images, while the u component captures the BV part of the image. The second example is applying image decomposition model to image denois- ing problems. In all our restoration examples, we have used a white Gaussian noise with standard deviation σ = 0.8, image values in each channel rank from 0 to 1 (this is equivalent to σ = 204 for image intensity ranging from 0 to 255). In the following three experiments, Fig. 3, Fig. 4 and Fig. 5, we 9 (a) f (b) u (c) CB model true image (b) v (c) CB residual Fig. 3. (a) noisy original image f . (b) u, (c) CB denoising model [15]. In the second row : (b) v component (v + 0.5 plotted), and (c) is the residual of CB denoising model (f − CB result + 0.5 plotted). compared image decomposition model with CB denoising model [15] which is one of the best color image denoising models using TV minimization. The denoising results in the top raw are similar; however, in the residual of the results (the v component), we see that image decomposition model have less edge information compared to the residual from CB denoising model. In all the experiments presented, the parameters have been tuned so that we show the best numerical results for both our color decomposition model as well as the CB denoising model. Notice that the CB denoising results displaid here are of the same quality as the ones presented in the original paper [15]. The final example, Fig. 6 is comparing different numerical implementations fo our color image decomposition model (5). When we solve the two coupled equations (7) and (9), we have many options. For example, we can treat both u and v as 3D vectors and use color TV model (two iterations: one for u and one for v), as in image (a), or we can treat both u and v with CB model and have four iterations (two sets of iterations for u and v each) as in image (c). In Fig. 6, we consider the noisy image displaid in Fig. 3. It shows the comparison between (a) using both color TV, (b) our model (CB model for 10 (a) f (b) u (c) CB model true image (b) v (c) CB residual Fig. 4. (a) noisy original image f . (b) u, (c) CB denoising model [15]. In the second row : (b) v component (v + 0.5 plotted), and (c) is the residual of CB denoising model (f − CB result + 0.5 plotted). u and vectorial TV for v), and (c) CB model for both u and v. The denoised results u are quite similar to each other; therefore, we only plot the noise v components (Top rows). Comparing (a) and (b), using CB model results in better color and detail control as in[15], and v component looks better in (b). Comparing (b) and (c), the results are almost similar (top row (b) and (c)); nevertheless, there are some difference. We separated v components to chromaticity component vc = v ‖v‖ and brightness component vb = ‖v‖ in second and third rows for better comparison (for image (c), vc and vb are given from the algorithm). Interesting