Simultaneous Inpainting and Denoising by Directional Global Three-part Decomposition: Connecting Variational and Fourier Domain Based Image Processing
SSimultaneous Inpainting and Denoising byDirectional Global Three-part Decomposition:Connecting Variational and Fourier DomainBased Image Processing
D.H. Thai ∗ and C. Gottschlich † Abstract
We consider the very challenging task of restoring images (i) whichhave a large number of missing pixels, (ii) whose existing pixels are cor-rupted by noise and (iii) the ideal image to be restored contains bothcartoon and texture elements. The combination of these three propertiesmakes this inverse problem a very difficult one. The solution proposed inthis manuscript is based on directional global three-part decomposition(DG3PD) [1] with directional total variation norm, directional G-normand (cid:96) ∞ -norm in curvelet domain as key ingredients of the model. Imagedecomposition by DG3PD enables a decoupled inpainting and denoisingof the cartoon and texture components. A comparison to existing ap-proaches for inpainting and denoising shows the advantages of the pro-posed method. Moreover, we regard the image restoration problem fromthe viewpoint of a Bayesian framework and we discuss the connectionsbetween the proposed solution by function space and related image rep-resentation by harmonic analysis and pyramid decomposition. Keywords
Image decomposition, variational calculus, inverse problems, image inpainting,image denoising, cartoon image, texture image, noise, residual image, featureextraction. ∗ Institute for Mathematical Stochastics, University of Goettingen, Goldschmidtstr. 7,37077 G¨ottingen, Germany, and Statistical and Applied Mathematical Science Institute(SAMSI), USA. Email: [email protected] † Institute for Mathematical Stochastics, University of Goettingen, Goldschmidtstr. 7,37077 G¨ottingen, Germany. Email: [email protected] a r X i v : . [ c s . C V ] J un Introduction and Related Work
Image enhancement and image restoration are two superordinate concepts inimage processing which encompass a plethora of methods to solve a multitudeof important real-world problems [2, 3]. Image enhancement has the goal ofimproving an input image for a specific application, e.g. in areas such as med-ical image processing, biometric recognition, computer vision, optical charac-ter recognition, texture recognition or machine inspection of surfaces [4, 5, 6].Methods for image enhancement can be grouped by the domain in which theyperform their operations: Images are processed in the spatial domain or Fourierdomain, or modified e.g. in the wavelet or curvelet domain [7]. Types of en-hancement methods include contextual filtering, e.g. for fingerprint image en-hancement [8, 9, 10], contrast enhancement, e.g. by histogram equalization [11],and image superresolution [12]. Image restoration is connected to the notionthat a given input image suffers from degradation and the goal is restore anideal version of it. Degradations are caused by various types of noise, missingpixels or blurring and their countermeasures are denoising, inpainting and de-blurring. In general, one has to solve a linear or nonlinear inverse problem toreconstruct the ideal image from its given degraded version. Denoising aims toremove noise from an image and denoising methods include total variation min-imization based approaches [13], the application of nonlocal means (NL-Means)[14] or other dictionaries of image patches for smoothing, and adaptive thresh-olding in the Wavelet domain [15]. Inpainting [16] is the filling-in of missingpixels from the available information in the image and it is applied for scratchremoval from scanned photographs, for occlusion filling, for removing objects orpersons from images (in image forgery [17] or for special effects), for filling-in ofpixels which were lost during the transmission of an image or left out on purposefor image compression [18]. Deblurring [19] addresses the removal of blurringartifacts and is not in the focus of this paper.Rudin, Osher, and Fatemi [20] pioneered two-part image decomposition bytotal variation (TV) regularization for denoising. Shen and Chan [21] appliedTV regularization to image inpainting, called TV inpainting model. and theyalso suggested image inpainting by curvature-driven diffusions (CDD), see [22].Starck et. al. [23] defined a model for two-part decomposition based on dic-tionary approach. Then, Elad et al. [24] applied this decomposition idea forimage inpainting by introducing the indicator function in the (cid:96) norm of theresidual, see Eq. (6) in [24]. Esedoglu and Shen [25] introduced two inpaintingmodels based on the Mumford-Shah model [26] and its higher order correction- the Mumford-Shah-Euler image model. They also presented numerical com-putation based on the Γ-convergence approximations [27, 28]. Shen et al. [29]proposed image inpainting based on bounded variation and elastica models fornon-textured images.Image inpainting can be an easy or difficult problem depending on the amountof missing pixels [22], the complexity of the image content and whether prior2nowledge about the image content is available. Methods have been proposedwhich perform only cartoon inpainting (also referred to as structure inpainting)[21, 29, 30] or only texture inpainting [31]. Images which consist of both cartoon(structure) and texture components are more challenging to inpaint. Bertalmio et al. [32], Elad et al. [24] and Cai et al. [33] have proposed methods forinpainting which can handle images with both cartoon (structure) and texturecomponents.In this paper, we tackle an even more challenging problem. Consider an inputimage f which has the following three properties:(i) a large percentage of pixels in f are missing and shall be inpainted.(ii) the known pixels in f are corrupted by noise.(iii) f contains both cartoon and texture elements.The co-occurrence of noise and missing pixels in an image with cartoon andtexture components increases the difficulty of both the inpainting problem andthe denoising problem. A multitude of methods has been proposed for inpaintingand denoising. Existing inpainting methods in the literature typically assumethat the non-missing pixels in a given image contain only a small amount of noiseor are noise-free, and existing methods for denoising typically assume all pixelsof the noisy image are known. The proposed method for solving this challengingproblem is inspired by the works of Efros and Leung [31], Bertalmio et al. [32],Vese and Osher [34], Aujol and Chambolle [35], Buades et al. [14] and Elad et al. [24], and it is based on the directional global three-part decomposition (DG3PD)[1]. The DG3PD method decomposes an image into three parts: a cartoonimage, a texture image and a residual image. Advantages of the DG3PD modellie in the properties which are enforced on the cartoon and texture images. Thegeometric objects in the cartoon image have a very smooth surface and sharpedges. The texture image yields oscillating patterns on a defined scale which isboth smooth and sparse. Recently, the texture images have been applied as avery useful feature for fingerprint segmentation [1, 36, 37].We address the challenging task of simultaneous inpainting and denoising inthe following way: The advanced DG3PD model introduced in the next sectiondecomposes a noisy input image f (with missing regions D ) into cartoon u ,texture v and residual (cid:15) components. At the same time, the missing regions D of the cartoon component u are interpolated and the available regions of u are denoised by the advantage of multi-directional bounded variation. Thiseffect benefits from the help of the indicator function in the measurement ofthe residual, i.e. (cid:13)(cid:13) C{ χ cD · × (cid:15) } (cid:13)(cid:13) (cid:96) ∞ in (1). However, texture v is not interpolateddue to the ”cancelling” effect of this supremum norm for residual in unknownregions. Therefore, the texture component v is inpainted and denoised by adictionary based approach instead. The DG3PD decomposition drives noise intothe residual component (cid:15) which is discarded. The reconstruction of the ideal3igure 1: Overview over the DG3PD image inpainting and denoising process.version of f is obtained by summation of the inpainted and denoised cartoonand texture components (see Figure 1 for an visual overview).Moreover, we uncover the link between the calculus of variations [38, 39, 40]and filtering in the Fourier domain [41] by analyzing the solution of the convexminimization in Eq. (1), i.e. roughly speaking the solution of the DG3PDinpainting model which can be understood as the response of the lowpass filter (cid:99) LP( ω ), highpass filter (cid:99) HP( ω ) and bandpass filter (cid:99) BP( ω ) and the unity conditionis satisfied, i.e. (cid:99) LP( ω ) + (cid:99) BP( ω ) + (cid:99) HP( ω ) = 1 , where ω ∈ [ − π , π ] is a coordinator in the Fourier domain. We observe thatthis decomposition is similar to wavelet or pyramidal decomposition scheme[42, 43, 44]. However, the basis elements obtaining the decomposition, i.e. scal-ing function and frame (or wavelet-like) function, are constructed by discretedifferential operators (due to the discrete setting in minimizing (1)) which arereferred to as wavelet-like operators in [45]. In particular, • the scaling function and wavelet-like function for the cartoon u are fromthe effect of the multi-directional TV norm, • the scaling function and wavelet-like function to extract the texture v arereconstructed by the effect of the multi-directional G-norm, • the effect of (cid:96) ∞ norm (cid:13)(cid:13) C{ χ cD · × (cid:15) } (cid:13)(cid:13) (cid:96) ∞ is to remove the remaining signalin the known regions of the residual (cid:15) (due to the duality property of (cid:96) ∞ ).We also describe flowcharts to show that the method of variational calculus(or the DG3PD inpainting) is a closed loop pyramidal decomposition which is4ifferent from an open loop one, e.g. wavelet [46], curvelet [47], see Figure 12. Bynumerics, we observe that the closed loop filter design by the calculus of variationwill result in lowpass, highpass and bandpass filters which are “unique” fordifferent images, see Figure 15. We also analyze the DG3PD inpainting modelfrom a view of the Bayesian framework and then define a discrete innovationmodel for this inverse problem.The setup of paper is as follows. In Section 2, we describe the DG3PD modelfor image inpainting and denoising. In Section 3, we show how to computethe solution of the convex minimization in the DG3PD inpainting problem byaugmented Lagrangian method. In Section 4, we describe the proposed methodfor texture inpainting and denoising. In Section 5, we compare the proposedmethod to existing ones (TVL2 inpainting, Telea [48] and Navier Stokes [49]).In Section 6, we consider our inverse problem from a statistical point of view,i.e. the Bayesian framework, to describe how to select priors for cartoon u andtexture v . We analyze the relation between the calculus of variations and thetraditional pyramid decomposition scheme, e.g. Gaussian pyramid, in Section7. We conclude the study with Section 8. For more detailed notation andmathematical preliminaries, we refer the reader to [36, 50]. We define a method for a restoration of the original noisy/non-noisy image f witha set of known missing regions D . The proposed model is a generalized versionof DG3PD [1] for the inpainting and denoising problem. This modification forour DG3PD-inpainting is inspired by [24].In particular, the discrete image f (size m × n ) with missing regions D issimultaneously decomposed into cartoon u , texture v and noise (cid:15) , especially,cartoon u is interpolated on the missing regions due to the indicator function χ D in the modified model.A set of missing regions D on a bounded domain Ω is defined by the indicatorfunction χ D whose complement is χ cD [ k ] = (cid:40) , k ∈ Ω \ D , k ∈ D .
Instead of putting χ cD on (cid:96) norm of the residual (see Eq. (4) and (6) in[24]), we introduce χ cD for the residual in the DG3PD model with the point-wise multiplication operator · × , i.e. (cid:13)(cid:13)(cid:13) C (cid:8) χ cD · × (cid:15) (cid:9)(cid:13)(cid:13)(cid:13) (cid:96) ∞ . We propose the DG3PD-5npainting as follows( u ∗ , v ∗ , (cid:15) ∗ , g ∗ ) = argmin ( u , v , (cid:15) , g ) ∈ X S (cid:26) L − (cid:88) l =0 (cid:13)(cid:13)(cid:13)(cid:13) cos (cid:0) πlL (cid:1) uD T n + sin (cid:0) πlL (cid:1) D m u (cid:13)(cid:13)(cid:13)(cid:13) (cid:96) (cid:124) (cid:123)(cid:122) (cid:125) := (cid:107) ∇ + L u (cid:107) (cid:96) + µ S − (cid:88) s =0 (cid:107) g s (cid:107) (cid:96) (cid:124) (cid:123)(cid:122) (cid:125) := (cid:107) v (cid:107) GS + µ (cid:107) v (cid:107) (cid:96) s.t. (cid:13)(cid:13)(cid:13) C (cid:8) χ cD · × (cid:15) (cid:9)(cid:13)(cid:13)(cid:13) (cid:96) ∞ ≤ ν , v = S − (cid:88) s =0 (cid:104) cos (cid:0) πsS (cid:1) g s D T n + sin (cid:0) πsS (cid:1) D m g s (cid:105)(cid:124) (cid:123)(cid:122) (cid:125) =div + S g , f = u + v + (cid:15) (cid:27) . (1)Note that if there are no missing regions, i.e. χ cD [ k ] = 1 , ∀ k ∈ Ω, the mini-mization problem (1) becomes the DG3PD model [1]. Next, we discuss how tosolve the DG3PD-inpainting model (1).
In this section we present a numerical algorithm for obtaining the solution ofthe DG3PD-inpainting model stated in (1). To simplify the calculation, weintroduce three new variables: r b = cos (cid:0) πbL (cid:1) uD T n + sin (cid:0) πbL (cid:1) D m u , b = 0 , . . . , L − , w a = g a , a = 0 , . . . , S − , e = χ cD · × (cid:15) . and denote G ∗ (cid:0) e ν (cid:1) as the indicator function on the feasible convex set A ( ν ) of(1), i.e. A ( ν ) = (cid:110) e ∈ X : (cid:13)(cid:13)(cid:13) C (cid:8) e (cid:9)(cid:13)(cid:13)(cid:13) (cid:96) ∞ ≤ ν (cid:111) and G ∗ (cid:0) e ν (cid:1) = (cid:40) , (cid:15) ∈ A ( ν )+ ∞ , else . Due to the constrained minimization problem in (1), the augmented Lagrangianmethod is applied to turn it into an unconstrained one asmin (cid:0) u , v , (cid:15) , e , (cid:2) r l (cid:3) L − l =0 , (cid:2) w s (cid:3) S − s =0 , (cid:2) g s (cid:3) S − s =0 (cid:1) ∈ X L +2 S +4 L (cid:16) · ; (cid:2) λ l (cid:3) L − l =0 , (cid:2) λ s (cid:3) S − s =0 , λ , λ , λ (cid:17) , (2)6here the Lagrange function is L ( · ; · ) = L − (cid:88) l =0 (cid:107) r l (cid:107) (cid:96) + µ S − (cid:88) s =0 (cid:107) w s (cid:107) (cid:96) + µ (cid:107) v (cid:107) (cid:96) + G ∗ (cid:0) e ν (cid:1) + β L − (cid:88) l =0 (cid:13)(cid:13)(cid:13)(cid:13) r l − cos (cid:0) πlL (cid:1) uD T n − sin (cid:0) πlL (cid:1) D m u + λ l β (cid:13)(cid:13)(cid:13)(cid:13) (cid:96) + β S − (cid:88) s =0 (cid:13)(cid:13)(cid:13)(cid:13) w s − g s + λ s β (cid:13)(cid:13)(cid:13)(cid:13) (cid:96) + β (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) v − S − (cid:88) s =0 (cid:104) cos (cid:0) πsS (cid:1) g s D T n + sin (cid:0) πsS (cid:1) D m g s (cid:105) + λ β (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:96) + β (cid:13)(cid:13)(cid:13)(cid:13) f − u − v − (cid:15) + λ β (cid:13)(cid:13)(cid:13)(cid:13) (cid:96) + β (cid:13)(cid:13)(cid:13)(cid:13) e − χ cD · × (cid:15) + λ β (cid:13)(cid:13)(cid:13)(cid:13) (cid:96) (cid:41) . Similar to [1], the alternating directional method of multipliers is applied tosolve the unconstrained minimization problem (2). Its minimizer is numericallycomputed through iterations t = 1 , , . . . (cid:16) u ( t ) , v ( t ) , (cid:15) ( t ) , e ( t ) , (cid:2) r ( t ) l (cid:3) L − l =0 , (cid:2) w ( t ) s (cid:3) S − s =0 , (cid:2) g ( t ) s (cid:3) S − s =0 (cid:17) =argmin L (cid:16) u , v , (cid:15) , e , (cid:2) r l (cid:3) L − l =0 , (cid:2) w s (cid:3) S − s =0 , (cid:2) g s (cid:3) S − s =0 ; (cid:2) λ ( t − l (cid:3) L − l =0 , (cid:2) λ ( t − s (cid:3) S − s =0 , λ ( t − , λ ( t − , λ ( t − (cid:17) (3)and the Lagrange multipliers are updated after every step t . We also initial-ize u (0) = f , v (0) = (cid:15) (0) = e (0) = (cid:2) r (0) l (cid:3) L − l =0 = (cid:2) w (0) s (cid:3) S − s =0 = (cid:2) g (0) s (cid:3) S − s =0 = (cid:2) λ (0) l (cid:3) L − l =0 = (cid:2) λ (0) a (cid:3) S − a =0 = λ (0) = λ (0) = λ (0) = , where is a m × n zeromatrix.In each iteration, we first solve the seven subproblems in the listed order:“ (cid:2) r l (cid:3) L − l =0 -problem”, “ (cid:2) w s (cid:3) S − s =0 -problem”, “ (cid:2) g s (cid:3) S − s =0 -problem”, “ v -problem”, “ u -problem”, “ e -problem”, “ (cid:15) -problem”, and then we update the five Lagrangemultipliers, namely (cid:2) λ (cid:3) L − l =0 , (cid:2) λ (cid:3) S − a =0 , λ , λ , λ , see Algorithm 1.In this section, we only present the solution of the “ e -problem” and “ (cid:15) -problem” as well as the updated Lagrange multiplier λ . We refer the read-ers to [1] for the detailed explanation of the other subproblems and the otherLagrange multipliers. The e-problem:
Fix u , v , (cid:15) , (cid:2) r l (cid:3) L − l =0 , (cid:2) w s (cid:3) S − s =0 , (cid:2) g s (cid:3) S − s =0 andmin e ∈ X (cid:110) G ∗ (cid:0) e ν (cid:1) + β (cid:13)(cid:13)(cid:13)(cid:13) e − (cid:16) χ cD · × (cid:15) − λ β (cid:17)(cid:13)(cid:13)(cid:13)(cid:13) (cid:96) (cid:111) (4)The solution of (4) is defined as a projection operator on a convex set A ( ν ), i.e.7 A ( ν ) (cid:16) χ cD · × (cid:15) − λ β (cid:17) : e ∗ = (cid:16) χ cD · × (cid:15) − λ β (cid:17) − CST (cid:16) χ cD · × (cid:15) − λ β , ν (cid:17) (5) The (cid:15) -problem:
Fix u , v , e , (cid:2) r l (cid:3) L − l =0 , (cid:2) w s (cid:3) S − s =0 , (cid:2) g s (cid:3) S − s =0 andmin (cid:15) ∈ X (cid:110) β (cid:13)(cid:13)(cid:13)(cid:13) f − u − v − (cid:15) + λ β (cid:13)(cid:13)(cid:13)(cid:13) (cid:96) + β (cid:13)(cid:13)(cid:13)(cid:13) e − χ cD · × (cid:15) + λ β (cid:13)(cid:13)(cid:13)(cid:13) (cid:96) (cid:111) (6)The solution of (6) with the point-wise division operator / . is (cid:15) ∗ = (cid:104) β (cid:0) f − u − v + λ β (cid:1) + β χ cD · × (cid:0) e + λ β (cid:1)(cid:105) / . (cid:104) β mn + β χ cD (cid:105) (7)and mn is a m × n matrix of ones. Updated Lagrange Multiplier λ ∈ X : λ = λ + β (cid:16) e − χ cD · × (cid:15) (cid:17) Choice of Parameters
Due to an adaptability to specific images and a balance between the smoothingterms and updated terms for the solution of the above subproblems, the selectionof parameters ( µ , µ , β , β , β ) is described in [1] and β is defined as β = c β . The term (cid:13)(cid:13) C{ χ cD · × (cid:15) } (cid:13)(cid:13) (cid:96) ∞ in (1) serves as a good measurement for the amountof noise due to the substraction in Eq. (5), see [51, 36, 1] (thanks to themultiscale and multidirection properties of the curvelet transform and its du-ality). However, due to the indicator function χ cD in the noise measurement (cid:13)(cid:13) C{ χ cD · × (cid:15) } (cid:13)(cid:13) (cid:96) ∞ , the interpolated texture by DG3PD inpainting is almost inthe residual, see Figure 3 (a), and its smoothed version by curvelet shrinkage(b). Because of the substraction operator in Eq. (5), these interpolated tex-tures are canceled out in e , see (c). Figure 3 (d) illustrates the estimated texturebefore the substraction in Eq. (5) as v texture = ( v + e · × χ D ) · × ROI , where ROI is the region of interest obtained by the morphological operator fortexture image in [37].In order to analyze the effects of the proposed model in terms of interpola-tion (for u ) and decomposition, we consider a 1-dimensional signal (which is8igure 2: An ideal image (a) with cartoon and texture components. Missingpixels (b) are shown in white. Image (c) combines missing pixels and noise. (cid:15) [ k ] ∼ N (0 , σ = 100 ) , k ∈ Ω.Figure 3: The “canceling effect” is caused by the noise measurement (cid:13)(cid:13) C{ χ cD · × (cid:15) } (cid:13)(cid:13) (cid:96) ∞ .extracted along the red line in Figure 4 (a)-(d)). By the DG3PD inpainting,the mean values of the cartoon u in (e) “nearly” remain unchanged. The ho-mogeneous regions have “sharp” edges and a “stair-case” effect does not occur(thanks to the directional TV). Texture v in (f) is extracted in areas whichcontain repeated pattern in (b). Moreover, small scales objects, e.g. noise, areremoved from v .However, the term (cid:13)(cid:13) C{ χ cD · × (cid:15) } (cid:13)(cid:13) (cid:96) ∞ causes a “canceling effect” which removesthe interpolated texture in unknown area during the decomposition process.Therefore, texture inpainting and denoised is tackled separately as described inthe next section by a generalized version of the dictionary learning for texturesynthesis [31]. 9igure 4: Interpolation for 1-dimensional signal after DG3PD inpainting (with-out texture inpainting): (c) is an 1-dimensional signal extracted along the redline of the original image (a). (d) is its corrupted signal by Gaussian noise anda missing regions D , see (b). Thanks to the directional TV norm, the DG3PDinpainting produces “sharp” edges without “stair-case” effect in a cartoon u (e)and the mean values of u “nearly” remains. Due to noise and the shrinkagesoft-thresholding operator by (cid:107) v (cid:107) (cid:96) in (1), texture v is reconstructed with a“sharp” transition on areas where texture presents in the original signal (c).However, due to a “heavy” noise (cid:15) on the known domains in (e), the DG3PDfails to reconstruct a “full” texture. In general, this is a challenging problemto recover an oscillating (or weak) signal corrupted by heavy noise. (h) is areconstructed image (without step of inpainting texture) and (i) is a quantizedsignal, i.e. truncate to [0, 255] and quantize.10 lgorithm 1 The Discrete DG3PD Inpainting Model
Initialization: u (0) = f , v (0) = (cid:15) (0) = (cid:2) r (0) l (cid:3) L − l =0 = (cid:2) w (0) s (cid:3) S − s =0 = (cid:2) g (0) s (cid:3) S − s =0 = (cid:2) λ (0) l (cid:3) L − l =0 = (cid:2) λ (0) a (cid:3) S − a =0 = λ (0) = λ (0) = . for t = 1 , . . . , T do1. Compute (cid:16)(cid:2) r ( t ) b (cid:3) L − b =0 , (cid:2) w ( t ) a (cid:3) S − a =0 , (cid:2) g ( t ) a (cid:3) S − a =0 , v ( t ) , u ( t ) , e ( t ) , (cid:15) ( t ) (cid:17) ∈ X L +2 S +4 : r ( t ) b = Shrink (cid:16) cos (cid:0) πbL (cid:1) u ( t − D T n + sin (cid:0) πbL (cid:1) D m u ( t − − λ ( t − b β , β (cid:17) , b = 0 , . . . , L − w ( t ) a = Shrink (cid:16) t w a := g ( t − a − λ ( t − a β , µ β (cid:17) , a = 0 , . . . , S − g ( t ) a = Re (cid:20) F − (cid:110) B ( t ) a ( z ) A ( t ) a ( z ) (cid:111)(cid:21) , a = 0 , . . . , S − v ( t ) = Shrink (cid:18) t v := β β + β (cid:18) S − (cid:88) s =0 (cid:104) cos (cid:0) πsS (cid:1) g ( t ) s D T n + sin (cid:0) πsS (cid:1) D m g ( t ) s (cid:105) − λ ( t − β (cid:19) + β β + β (cid:18) f − u ( t − − (cid:15) ( t − + λ ( t − β (cid:19) , µ β + β (cid:19) u ( t ) = Re (cid:20) F − (cid:110) Y ( t ) ( z ) X ( t ) ( z ) (cid:111)(cid:21) e ( t ) = (cid:16) χ cD · × (cid:15) ( t − − λ ( t − β (cid:17) − CST (cid:16) χ cD · × (cid:15) ( t − − λ ( t − β , ν (cid:17) , ν = c ν (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) C{ χ cD · × (cid:15) ( t − − λ ( t − β } (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:96) ∞ (cid:15) ( t ) = (cid:104) β (cid:0) f − u ( t ) − v ( t ) + λ ( t − β (cid:1) + β χ cD · × (cid:0) e ( t ) + λ ( t − β (cid:1)(cid:105) / . (cid:104) β mn + β χ cD (cid:105)
2. Update (cid:16)(cid:2) λ ( t ) b (cid:3) L − b =0 , (cid:2) λ ( t ) a (cid:3) S − a =0 , λ ( t ) , λ ( t ) , λ ( t ) (cid:17) ∈ X L + S +3 : λ ( t ) b = λ ( t − b + γβ (cid:16) r ( t ) b − cos (cid:0) πbL (cid:1) u ( t ) D T n − sin (cid:0) πbL (cid:1) D m u ( t ) (cid:17) , b = 0 , . . . , L − λ ( t ) a = λ ( t − a + γβ (cid:16) w ( t ) a − g ( t ) a (cid:17) , a = 0 , . . . , S − λ ( t ) = λ ( t − + γβ (cid:16) v ( t ) − S − (cid:88) s =0 (cid:2) cos (cid:0) πsS (cid:1) g ( t ) s D T n + sin (cid:0) πsS (cid:1) D m g ( t ) s (cid:3)(cid:17) λ ( t ) = λ ( t − + γβ (cid:0) f − u ( t ) − v ( t ) − (cid:15) ( t ) (cid:1) λ ( t ) = λ ( t − + β (cid:16) e ( t ) − χ cD · × (cid:15) ( t ) (cid:17) end for a ( z ) = β mn + β (cid:12)(cid:12)(cid:12)(cid:12) sin πaS ( z −
1) + cos πaS ( z − (cid:12)(cid:12)(cid:12)(cid:12) , B a ( z ) = β (cid:104) W a ( z ) + Λ a ( z ) β (cid:105) + β (cid:104) sin (cid:0) πaS (cid:1) ( z − −
1) + cos (cid:0) πaS (cid:1) ( z − − (cid:105) × (cid:20) V ( z ) − (cid:88) s =[0 ,S − \{ a } (cid:104) cos (cid:0) πsS (cid:1) ( z −
1) + sin (cid:0) πsS (cid:1) ( z − (cid:105) G s ( z ) + Λ ( z ) β (cid:21) , X ( z ) = β mn + β L − (cid:88) l =0 (cid:12)(cid:12)(cid:12)(cid:12) sin (cid:0) πlL (cid:1) ( z −
1) + cos (cid:0) πlL (cid:1) ( z − (cid:12)(cid:12)(cid:12)(cid:12) , Y ( z ) = β (cid:104) F ( z ) − V ( z ) − E ( z ) + Λ ( z ) β (cid:105) + β L − (cid:88) l =0 (cid:104) sin (cid:0) πlL (cid:1) ( z − −
1) + cos (cid:0) πlL (cid:1) ( z − − (cid:105)(cid:104) R l ( z ) + Λ l ( z ) β (cid:105) . Choice of Parameters µ = c µ β · max k ∈ Ω (cid:0)(cid:12)(cid:12) t w a [ k ] (cid:12)(cid:12) (cid:1) , µ = c µ ( β + β ) · max k ∈ Ω (cid:0)(cid:12)(cid:12) t v [ k ] (cid:12)(cid:12) (cid:1) , β = c β , β = θ − θ β , β = c β , β = c β . The proposed method for reconstructing the texture component combines thefollowing ideas in five subsequent steps. • Texture extraction by DG3PD [1]. • Morphological processing [36, 37]. • Texture inpainting inspired by the work of Efros and Leung [31]. • Texture denoising by nonlocal means [14]. • Image synthesis by summation with the cartoon component u .The DG3PD model enforces sparsity and smoothness on the obtained texturecomponent v . Sparsity means that the vast majority of coefficients in v areequal to zero (the percentage of zero coefficients depends on how much texturea specific image contains). We make good use of this property to answer thequestion: which pixels of the missing shall be inpainted with texture? First,the texture component v is segmented into zero (gray pixels in Figure 5(a))and nonzero coefficients (black and white pixels in Figure 5(a)). Morphologicalprocessing (as described in [36, 37] for fingerprint segmentation) obtains textureregions (see Figure 5(b)). The mask of missing pixels D shown in (c) is dilated12igure 5: All coefficients of the obtained texture component v equal to zero arevisualized by gray pixels in (a), positive coefficients are indicated by white pixelsand negative coefficients by black pixels. (b) depicts the resulting segmentation.The mask of missing pixels D shown in (c) is dilated and then intersected with(b), leading to the (white) pixels which are to be inpainted with texture in (d).(we used a circular structure element with 5 pixels radius in order to avoidborder effects on the margin between existing and missing pixels) and then, itis intersected with the texture region segmentation to obtain the mask I shownin Figure 5(d) which are the pixels to be inpainted with texture.The proposed texture inpainting proceeds in two phases: First, we build adictionary of texture patches and second, we inpaint all pixels of mask I in aspecific order. For the dictionary, we select all patches of size s × s pixels (we used s = 15 in our experiments) from texture component v which meet the followingtwo criteria: At least p percent of the pixels are known (and additional, thecentral pixel of the patch has to exist) and at least p percent of the coefficientsare nonzero (we used p = 70% and p = 60% in our experiments). The firstcriterion excludes patches which contain too many missing (unknown) pixelsand the second criterion excludes patches without texture from the dictionary.Next, we iterate over all pixels of mask I and consider the pixel to be inpaintedas the central pixel of an image patch with size s × s . We count the numberof known pixels inside the patch and compute the percentage of known pixels.Pixels are inpainted in the following order. We start with a threshold percentage t = 90% which is decreased in steps of 5% per iteration and in each iteration, weinpaint all pixels of mask I for which at least t percent are known. The rationalebehind this ordering is that the more pixels are known (or already inpainted) inthe neighborhood around a missing pixels, the better it can inpainted, becausethis additional information improves the chances of finding a good match in thetexture dictionary. A third constraint ensures that the overlap of known pixelsbetween a dictionary patch and the patch around the missing pixels containsat least p percent of known pixels (we used p = 30% in our experiments).Inpainting a pixel means that we find the best fitting patch in our texturedictionary and set the pixel value to the central pixel of that patch. Best fittingis defined as minimum sum of squared differences per pixel, divided by the13igure 6: Images (a) to (e) visualize the inpainting progress (added to thecartoon component) from 0 to 100% in steps of 25%. Images (f) to (h) displaythe denoising results after one, five and 20 iterations, respectively.number of pixels which overlap.After all missing texture pixels have been inpainted the texture region is de-noised n iterations of nonlocal means. In each iteration, we first construct adictionary considering all patches in the texture region with known and pre-viously inpainted pixels. Next, for each pixel to be denoised, we find the top k fitting patches in the dictionary using the same distance function (sum ofsquared pixelwise differences) and set the denoised pixel to the average value ofthe top k central pixel values. In our experiments, we used k = 5 and we haveobserved that after about n = 10 iterations, the image has reached a steadystate. See Figure 6 for a visualization of the status at several intermediate stepsduring the inpainting and denoising process.Finally, the full image is reconstructed by summation of the inpainted car-toon component u with the inpainted and denoised texture component v , seeFigure 9. 14 Comparison of DG3PD to Further InpaintingMethods
Figure 7 shows inpainting results for the considered challenging problem (seeFigure 2 (c)) obtained by well-known inpainting methods from the literature: anapproach based on Navier Stokes equations from fluid dynamics which has beenproposed by Bertalmio et al. [49], an inpainting method suggested by Telea [48]and a well-known TVL2 inpainting approach with its synthesis image shown inFigure 7 (b) and (c).The image shown in Figure 2(c) has the three properties discussed in Section 1:(i) a large percentage of pixels in f are missing and shall be inpainted.(ii) the known pixels in f are corrupted by noise.(iii) f contains both cartoon and texture elements.Furthermore, a noise-free, ideal version f depicted in Figure 2(a) is availablewhich serves as ground truth for evaluating the quality of inpainted and denoisedimages by different methods. The availability of a noise-free image is an advan-tage for comparing and evaluating different inpainting and denoising methods.In contrast, images taken by a digital camera contain a certain amount of noiseand for them, an ideal version is not available, see e.g. the Barbara image weused in previous comparisons [1]. The choice of the ideal image (shown in Figure2) enables to clearly see and understand the advantages and limitations of thecompared methods (see Figure 7).Due to the fidelity term, the L -norm, in the TVL2 inpainting model, noise isreduced and the homogeneous areas are well interpolated. However, the methodis known to cause the “stair case” effect on the homogeneous regions and texture,while small scale objects tend to be eliminated. If the parameter β is chosensmaller (e.g. β = 0 . In this section, we describe the DG3PD inpainting model from a view of theBayesian framework through maximum a posterior (MAP) estimator. We as-sume that u and v are independent random variables (r.v.).15igure 7: Comparison of inpainting results obtained by (a) Navier Stokes [49],by (b) Telea [48], by (c,d)TVL2 and by (e) DG3PD.16 Prior of cartoon image u:
A cartoon u = (cid:2) u [ k ] (cid:3) k ∈ Ω consists of elementpixel u [ k ] which is considered as an independent r.v. with a prior p U k .Given k ∈ Ω, denote r [ k ] = ∇ + L u [ k ] = (cid:104) r l [ k ] (cid:105) L − l =0 . To describe the distri-bution of a r.v. U k , we firstly define the distribution of L -dimensional r.v. R k = [ R ,k , . . . , R L − ,k ]. We assume that r.v. R k has a multi-dimensionalLaplace distribution which is a part of multi-dimensional power exponen-tial distribution [52, 53], i.e. R k ∼ PE L (0 , Σ , ) = Laplace L (0 , Σ) withΣ = (cid:18) γ γ (cid:19) : p R k ( r [ k ] | , Σ ,
12 ) = 18 πγ exp (cid:110) − (cid:16) r [ k ]... r L − [ k ] T γ (cid:18) γ γ (cid:19) r [ k ]... r L − [ k ] (cid:17) (cid:111) = 18 πγ exp (cid:110) − γ (cid:104) L − (cid:88) l =0 (cid:12)(cid:12) r l [ k ] (cid:12)(cid:12) (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) := | r [ k ] | = | ∇ + L u [ k ] | (cid:111) . Thus, the distribution of a r.v. U k is a multi-dimensional Laplace distri-bution with operator ∇ + L : p U k ( u [ k ]) = 18 πγ exp (cid:110) − γ (cid:12)(cid:12)(cid:12) ∇ + L u [ k ] (cid:12)(cid:12)(cid:12) (cid:111) . The joint probability density function (p.d.f.) of a prior u is p U ( u ) = (cid:89) k ∈ Ω p U k ( u [ k ]) = 1(8 πγ ) | Ω | exp (cid:110) − γ (cid:88) k ∈ Ω (cid:12)(cid:12)(cid:12) ∇ + L u [ k ] (cid:12)(cid:12)(cid:12) (cid:111) = 1(8 πγ ) | Ω | e − γ (cid:107) ∇ + L u (cid:107) (cid:96) . We choose γ = . The potential function of u in a matrix form isΦ U ( u ) = − log p U ( u ) = (cid:13)(cid:13)(cid:13) ∇ + L u (cid:13)(cid:13)(cid:13) (cid:96) + | Ω | log(2 π ) . Note that the original u is not a Laplace distribution, but a transform of u under an operator ∇ + L has an independent multi-dimensional Laplacedistribution. • Prior of texture image v:
Under a transform operator T , the textureimage v is decomposed in different orientations. Operator T is suitablychosen to capture the texture components in the original image f . Asdefinition of discrete multi-dimensional G-norm [1], the definition of thetransform T (in a direction s , s = 0 , . . . , S −
1) is defined by its inversion(note that
T T − = Id): T { v } = (cid:2) T s { v } (cid:124) (cid:123)(cid:122) (cid:125) := g s (cid:3) S − s =0 = g ⇔ v = T − (cid:110)(cid:2) g s (cid:3) S − s =0 (cid:111) = div + S g .
17e assume that a texture v is sparse in a transform domain T and sparsein the spatial domain (nonzero coefficients of v are only due to texture).To satisfy these conditions, we assume that a r.v. V k has a mixture ofa Laplace distribution in spatial domain with a p.d.f. p V k , ( v [ k ]) (w.r.t.sparse in spatial domain) and multi-dimensional Laplace distribution (inan operator T ) with a p.d.f. p V k , ( v [ k ]) (w.r.t. sparse in transform domain T ). Thus, we have a p.d.f. of r.v. V k as follows p V k ( v [ k ]) = p V k ( v [ k ]) · p V k ( v [ k ]) . (8)Since V k ∼ Laplace(0 , γ ), we have p V k ( v [ k ]) = 14 γ exp (cid:110) − γ (cid:12)(cid:12) v [ k ] (cid:12)(cid:12) (cid:111) . (9)To define the distribution of a r.v. V k , we define the distribution of S -dimension r.v. G k = [ G ,k , . . . , G L − ,k ]. Similar to U k , we assumethat r.v. G k has a multi-dimensional Laplace distribution, i.e. G k ∼ Laplace S (0 , Σ) with Σ = (cid:18) γ γ (cid:19) . It is easy to see that (with g [ k ] = T { v } [ k ] , k ∈ Ω) p G k ( g [ k ]) = 18 πγ exp (cid:110) − γ (cid:12)(cid:12) g [ k ] (cid:12)(cid:12) (cid:111) or p V k ( v [ k ]) = 18 πγ exp (cid:110) − γ (cid:12)(cid:12) T { v } [ k ] (cid:12)(cid:12) (cid:111) . (10)Put (9) and (10) to (8), we have the p.d.f. of r.v. V k as p V k ( v [ k ]) = 132 πγ γ exp (cid:110) − γ (cid:12)(cid:12) T { v } [ k ] (cid:12)(cid:12) − γ (cid:12)(cid:12) v [ k ] (cid:12)(cid:12) (cid:111) . The joint p.d.f. of a texture image v = (cid:2) v [ k ] (cid:3) k ∈ Ω is p V ( v ) = (cid:89) k ∈ Ω p V k ( v [ k ]) = 1(32 πγ γ ) | Ω | exp (cid:110) − γ (cid:13)(cid:13) T { v } (cid:13)(cid:13) (cid:96) − γ (cid:107) v (cid:107) (cid:96) (cid:111) . we choose µ = γ and µ = γ . The potential function of v with ananisotropic version (cid:0)(cid:13)(cid:13) T { v } (cid:13)(cid:13) (cid:96) = (cid:80) S − s =0 (cid:13)(cid:13) T s { v } (cid:13)(cid:13) (cid:96) = (cid:80) S − s =0 (cid:107) g s (cid:107) (cid:96) (cid:1) in amatrix form is defined asΦ V ( v ) = µ S − (cid:88) s =0 (cid:107) g s (cid:107) (cid:96) + µ (cid:107) v (cid:107) (cid:96) + | Ω | log 4 πµ µ . • Likelihood (or the joint p.d.f. of p F | U,V ( f | u , v ) ): If we assume thatthe residual (cid:15) in an image f has a power exponential distribution [52, 53],18.e. E k ∼ PE( µ, σ , ξ ), e.g. Gaussian ( ξ = 1) or Laplacian ( ξ = ), itsdensity function at k ∈ Ω is p E k ( (cid:15) [ k ] | µ ) = 1Γ(1 + ξ )2 ξ σ exp (cid:110) − (cid:12)(cid:12) (cid:15) [ k ] − µ (cid:12)(cid:12) ξ σ ξ (cid:111) . Due to the inpainting problem with a missing region D , the likelihood isdefined on a known domain Ω \ D as p F | U,V ( f | u , v ) = (cid:89) k ∈ Ω p F k | U k ,V k ( f [ k ] | u [ k ] , v [ k ])= 1 (cid:16) Γ(1 + ξ )2 ξ σ (cid:17) | Ω | exp (cid:110) − σ ξ (cid:88) k ∈ Ω (cid:12)(cid:12) χ cD [ k ]( f [ k ] − u [ k ] − v [ k ]) (cid:12)(cid:12) ξ (cid:124) (cid:123)(cid:122) (cid:125) := (cid:107) χ cD · × ( f − u − v ) (cid:107) ξ ξ (cid:111) . • A posterior:
Since u and v are independent to each other, a posterioris written as p U,V | F ( u , v | f ) = p F | U,V ( f | u , v ) p U ( u ) p V ( v ) p F ( f ) ∝ p F | U,V ( f | u , v ) p U ( u ) p V ( v ) . Let (cid:15) = f − u − v , the MAP estimator, i.e. max ( u , v ) ∈ X p U,V | F ( u , v | f ), isdefined asmin ( u , v ) ∈ X (cid:110) − log p U ( u ) (cid:124) (cid:123)(cid:122) (cid:125) =Φ U ( u ) − log p V ( v ) (cid:124) (cid:123)(cid:122) (cid:125) =Φ V ( v ) − log p U,V | F ( u , v | f ) (cid:111) = min ( u , v ) ∈ X (cid:110)(cid:13)(cid:13)(cid:13) ∇ + L u (cid:13)(cid:13)(cid:13) (cid:96) + µ S − (cid:88) s =0 (cid:107) g s (cid:107) (cid:96) + µ (cid:107) v (cid:107) (cid:96) + 12 σ ξ (cid:13)(cid:13) χ cD · × (cid:15) (cid:13)(cid:13) ξ ξ , s.t. v = S − (cid:88) s =0 (cid:2) cos( πsS ) g s D T n + sin( πsS ) D m g s (cid:3) , f = u + v + (cid:15) (cid:111) . However, in practice, we do not know which types of noise are observed in asignal. Instead of the (cid:96) ξ norm of the residual (cid:15) on a known domain Ω \ D tocharacterize the properties of noise, we control the smoothness of the solutionby the maximum of curvelet coefficient of (cid:15) on a domain Ω \ D by a constant ν ,i.e. (cid:13)(cid:13) C{ χ cD · × (cid:15) } (cid:13)(cid:13) (cid:96) ∞ ≤ ν .In [54], if noise in image f is Gaussian, as in extreme value theory, the r.v. (cid:13)(cid:13) C{ χ cD · × (cid:15) } (cid:13)(cid:13) (cid:96) ∞ has a Gumbel distribution. Thus, there is a condition to choose ν . However, in practice, if noise cannot be identified, ν is chosen by α -quantile.Note that the condition (cid:13)(cid:13) C{ χ cD · × (cid:15) } (cid:13)(cid:13) (cid:96) ∞ ≤ ν is similar to the Dantzig selector[55]. 19igure 8 illustrates the empirical density functions of the solution of the min-imization (3). The statistical properties of the solution can be characterized inthe Bayesian framework with the priors, likelihood and posterior as mentionedabove.Note that the DG3PD inpainting model (1) can be described as an inverseproblem in image restoration (e.g. image inpainting and denoising), see Figure9 for the discrete innovation model. In this section, we discuss connections and similarities between the proposedDG3PD inpainting model and existing work in the area of pyramid represen-tations, multiresolution analysis, and scale-space representation. Historically,Gaussian pyramids and Laplacian pyramids have been developed for applica-tions like texture synthesis, image compression and denoising. Early works inthis area include a paper by Burt and Adelson in 1983 [42] applying Laplacianpyramids for image compression. Pyramidal decomposition schemes [42, 43, 44]can be grouped into two categories. Their main property is either lowpass orbandpass filtering. Very related are transforms from multiresolution analysis,including wavelet, curvelet [47, 51], contourlet [57], and steerable wavelet [58]. Adifference of the proposed DG3PD inpainting model is that its basis elements forobtaining the decomposition are constructed by discrete differential operatorsand the decomposition is solved in a non-linear way which adapts to each image(enabled by update iterations, corresponding to the loop in Figure 12). Thediscussion of commonalities and differences is made more precise and detailedin the following rest of the section.Let denote F ←→ a discrete Fourier transform pair. Given z = e jω and z = e jω and the Dirac delta function δ ( · ), the impulse response of the discretedirectional derivative operators ( l = 0 , . . . , L −
1) and their spectra are • the forward operator: ∂ + l δ [ k ] = (cid:104) cos( πlL ) ∂ + x +sin( πlL ) ∂ + y (cid:105) δ [ k ] F ←→ cos( πlL )( z − πlL )( z − , • the backward operator: ∂ − l δ [ k ] = (cid:104) cos( πlL ) ∂ − x +sin( πlL ) ∂ − y (cid:105) δ [ k ] F ←→ − cos( πlL )( z − − − sin( πlL )( z − − . In order to describe the variational analysis for the u -problem from the view offiltering in the Fourier domain for a pyramid decomposition scheme, we define20igure 8: The figure depicts the empirical density functions of the solutionin (3) by ADMM: (a), (b), (c), (d), (e), (f) and (g) are the empirical densityfunctions of u , v , (cid:15) , r , g , w and the QQ-plot of the residual (cid:15) , respectively.We assume that the prior of ∇ + L u has a multi-Laplacian distribution, i.e. achanged variable r = (cid:2) r l (cid:3) L − l =0 = ∇ + L u has a multi-Laplace distribution, or r l has an univariate Laplace distribution. This effect causes the sparseness of r l due to the shrinkage operator, see (d). The sparsity for ∇ + L u causes thesmoothness of the objects in the cartoon u , see (a). The smooth and sparsetexture v (due to the expansion of a whole range in intensity value of its non-zero coefficients caused by (cid:107) v (cid:107) (cid:96) in (1)) is illustrated in (b). Because of anadditive white Gaussian noise in a simulation, the distribution of (cid:15) in Ω \ D is approximately normal, see plot (c) and its QQ-plot (g). The reasons of thedifferences near the boundary in (g) may cause by: (1) a simulation of a Gaussiannoise, i.e. the tail of a Gaussian noise does not go to infinity in a simulation. (2)the remaining texture in (cid:15) (due to a selection of ν ). Changing variable w = g causes a non-sparsity in g and a sparsity in w (which is shrinked from g ), c.f.(e) and (f). 21igure 9: Overview over the discrete innovation model for the DG3PD basedimage inpainting and denoising. Note that the L − does not exist in general,but the reconstructed image is obtained by the function space of the inverseproblem.the discrete inverse Fourier transform of X ( z ) and Y ( z ) in the “ u -problem” in[1] as X ( z ) F − ←→ X [ k ] = (cid:104) β − β L − (cid:88) l =0 ∂ − l ∂ + l (cid:105) δ [ k ] , Y ( z ) F − ←→ Y [ k ] = β (cid:0) f [ k ] − v [ k ] − (cid:15) [ k ] + λ [ k ] β (cid:17) − β L − (cid:88) l =0 ∂ − l (cid:110) Shrink (cid:16) ∂ + l u [ k ] − λ l [ k ] β , β (cid:17)(cid:124) (cid:123)(cid:122) (cid:125) := r l [ k ] + λ l [ k ] β (cid:111) . Note that the function X ( z ) satisfies0 < β ≤ X ( z ) < + ∞ , ∀ ω ∈ [ − π , π ] and it is similar to the autocorrelation function in [56] which follows the condi-tion of the Riesz basis. We rewrite cartoon u in the Fourier and spatial domainin the scheme of filter design as U ( z ) = X − ( z ) · Y ( z ) = Φ( z ) (cid:16) F ( z ) − V ( z ) − E ( z ) + Λ ( z ) β (cid:17) + L − (cid:88) l =0 ˜Ψ( z ) (cid:104) R l ( z ) + Λ l ( z ) β (cid:105) F − ←→ u [ k ] = (cid:16) φ ∗ ( f − v − (cid:15) + λ β ) (cid:17) [ k ] (cid:124) (cid:123)(cid:122) (cid:125) := u update [ k ] + L − (cid:88) l =0 (cid:16) ˜ ψ l ∗ (cid:104) Shrink (cid:16) ψ l ∗ u − λ l β , β (cid:17) + λ l β (cid:105)(cid:17) [ k ] (cid:124) (cid:123)(cid:122) (cid:125) := FST ψ, ˜ ψ,c ( u , λ lβ , β ) := u smooth [ k ] (frame soft thresholding) (11)22iven β = c β , the spectra and impulse responses of the scaling and framefunctions in (11) ( l = 0 , . . . , L −
1) areΦ( z ) = β X − ( z ) F ←→ φ [ k ] = (cid:104) − c L − (cid:88) l =0 ∂ − l ∂ + l (cid:105) − δ [ k ] , ˜Ψ l ( z ) = β X − ( z ) (cid:104) ( z − −
1) sin πlL + ( z − −
1) cos πlL (cid:105) F ←→ ˜ ψ l [ k ] = − c (cid:104) − c L − (cid:88) l =0 ∂ − l ∂ + l (cid:105) − ∂ − l δ [ k ] , Ψ l ( z ) = ( z −
1) sin πlL + ( z −
1) cos πlL F ←→ ψ l [ k ] = ∂ + l δ [ k ] . The equation (11) is somewhat similar to the projection operator in the pyramiddecomposition scheme (see curvelet [47, 51], contourlet [57], steerable wavelet[58], etc.), with: • scaling function ϕ ( · ) and its dual ˜ ϕ ( · ) with the interpolant φ = ϕ ∗ ˜ ϕ , • frame function ψ l ( · ) and its dual ˜ ψ l ( · ), • proximity operator for frame elements, i.e. Shrink( · , · ).Note that the frame soft thresholding FST( · , · , · ) in (11) consists of frame element ψ l , its dual ˜ ψ l and shrinkage operator Shrink( · , · ) together with the updated La-grange multiplier. Although FST( · , · , · ) is similar to wavelet/curvelet shrinkageoperator, this operator is obtained from the calculus of variation.It is easy to obtain the unity property of frame elements ( φ , ψ l , ˜ ψ l ) in theFourier domain as Φ( z ) + L − (cid:88) l =0 ˜Ψ l ( z )Ψ l ( z ) = 1 , Φ( e j ) = 1 and ˜Ψ l ( e j ) = Ψ l ( e j ) = 0 , which satisfies the condition of the perfect reconstruction in the scheme ofwavelet pyramid, see Figure 10 for illustration with L = 4 directions. When thenumber of directions L increases, the bandwidth of the scaling function Φ( z ) inthe spectral domain reduces, see Figure 17. Due to the effect of this lowpassfilter, the cartoon u is smoother. These frame elements ( φ, ψ l , ˜ ψ l ) which consistof partial differential operators can be considered as a kind of a wavelet-likeoperator [45]. However, the procedure to obtain these elements is from thecalculus of variation. (cid:2) g s (cid:3) S − s =0 -problem”: In analogy to the u -problem in the variational framework, we define the frameoperators to extract the directional features of the texture v in pyramid decom-23osition scheme. The discrete inverse Fourier transforms of A ( z ) and B ( z ) inthe solution of the “ v -problem” [1] are defined as A a ( z ) F − ←→ A a [ k ] = (cid:104) β − β ∂ − a ∂ + a (cid:105) δ [ k ] , B a ( z ) F − ←→ B a [ k ] = β (cid:104) Shrink (cid:16) g a [ k ] − λ a [ k ] β , µ β (cid:17) + λ a [ k ] β (cid:105) − β ∂ − a (cid:110) v [ k ] − (cid:88) s =[0 ,S − \{ a } ∂ + s g s [ k ] (cid:124) (cid:123)(cid:122) (cid:125) := ∂ + a g a [ k ] + λ [ k ] β (cid:111) Similar to X ( z ), the “autocorrelation” function A a ( z ) is bounded as in thecondition of the Riesz basis [56], i.e.0 < β ≤ A a ( z ) < + ∞ , ∀ ω ∈ [ − π , π ] . Texture v is rewritten in the Fourier and spatial domain with a = 0 , . . . , S − G a ( z ) = A − a ( z ) · B a ( z ) = Ξ a ( z ) (cid:104) W a ( z ) + Λ a ( z ) β (cid:105) + Ψ (cid:48) a ( z ) (cid:104) Ψ a ( z ) G a ( z ) + Λ ( z ) β (cid:105) F − ←→ g a [ k ] = (cid:18) ξ a ∗ (cid:104) Shrink (cid:16) g a − λ a β , µ β (cid:17) + λ a β (cid:124) (cid:123)(cid:122) (cid:125) := ST (cid:0) g a , λ a β , µ β (cid:1) (cid:105)(cid:19) [ k ] + (cid:18) ψ (cid:48) a ∗ (cid:104) ψ a ∗ g a + λ β (cid:105)(cid:19) [ k ] . (12)Given β = c β , the spectra and impulse responses of the frame functions in(12) areΞ a ( z ) = β A − ( z ) F − ←→ ξ a [ k ] = c (cid:0) c − ∂ − a ∂ + a (cid:1) − δ [ k ] , Ψ (cid:48) a ( z ) = β A − ( z ) (cid:104) ( z − −
1) sin πaS + ( z − −
1) cos πaS (cid:105) F − ←→ ψ (cid:48) a [ k ] = − (cid:0) c − ∂ − a ∂ + a (cid:1) − ∂ − a δ [ k ]Ψ a ( z ) = ( z −
1) cos πaS + ( z −
1) sin πaS F − ←→ ψ a [ k ] = ∂ + a δ [ k ] . Similar to the u -problem, it is easy to obtain the unity property of frameelements ( ξ , ψ a , ψ (cid:48) a ) with a ∈ [0 , S −
1] in Fourier domain:Ξ a ( z ) + Ψ a ( z )Ψ (cid:48) a ( z ) = 1 , Ξ a ( e j ) = 1 and Ψ a ( e j ) = Ψ (cid:48) a ( e j ) = 0which satisfies the condition of the perfect reconstruction, see Figure 10 with L = 4 directions and Figures 17 and 18 with L = 8 directions.24 .3 The “v-problem”: The solution of the v -problem is rewritten as v ∗ = Shrink (cid:16) t v , c µ · max k ∈ Ω (cid:0)(cid:12)(cid:12) t v [ k ] (cid:12)(cid:12) (cid:1)(cid:17) with t v = θ (cid:104) S − (cid:88) s =0 ψ s ∗ g s − λ β (cid:124) (cid:123)(cid:122) (cid:125) t v smooth (cid:105) + (1 − θ ) (cid:104) v + λ β (cid:124) (cid:123)(cid:122) (cid:125) t v update (cid:105) . Note that there is a shrinkage operator with two terms in texture v , namely asmoothing term and a updated term, balanced by a parameter θ . The noise term (cid:13)(cid:13) C{ χ cD · × (cid:15) } (cid:13)(cid:13) (cid:96) ∞ in Eq. (1) is suitable to capture high oscillatingpatterns, especially small-scale types of noise, due to the advantage of the multi-orientation and multi-scale in a curvelet domain and the supremum norm. Fora convenient calculation, we introduce a variable e as a residual (cid:15) in a knowndomain Ω \ D . The explanation for using the supremum norm of (cid:15) in the curveletdomain can be found in (5). In particular, there are two terms in (5), includinga residual in Ω \ D with an updated Lagrange multiplier λ , i.e. (cid:16) χ cD · × (cid:15) − λ β (cid:17) ,and its curvelet smoothing term at a level ν , i.e. CST (cid:16) χ cD · × (cid:15) − λ β , ν (cid:17) . Assumeat an iteration t , there are some remaining signals, e.g. texture, in residual (cid:15) inΩ \ D . The curvelet soft-thresholding operator CST( · , · ) reduces noise (or smallscale objects) in (cid:0) χ cD · × (cid:15) − λ β (cid:1) at a level ν in different scales and orientations.By a subtraction operator in (5), e at a iteration t contains mainly noise orsmall scale objects, see Figure 11. (cid:15) -problem”: At iteration t , we denote the updated residual by the unity condition as (cid:15) ( t )unity = f − u ( t ) − v ( t ) + λ ( t − β , and rewrite (7) with the indicator functions on unknown domain D and knowndomain Ω \ D as (cid:15) ( t ) = (cid:104) β β + β (cid:15) ( t )unity + β β + β (cid:0) e ( t ) + λ ( t − β (cid:1)(cid:105)(cid:124) (cid:123)(cid:122) (cid:125) := (cid:15) ( t )known · × χ cD + β β + β (cid:15) ( t )unity (cid:124) (cid:123)(cid:122) (cid:125) := (cid:15) ( t )unknown · × (1 − χ cD ) .
25e see that there are two terms in the residual (cid:15) ( t ) , including the updated term (cid:15) ( t )known for Ω \ D and (cid:15) ( t )unknown for D . And there are other two terms for up-dating (cid:15) ( t )known in Ω \ D , namely (cid:15) ( t )unity (from an unity condition) and e ( t ) (froma curvelet smoothing operator). In a contrast to an open loop in a pyramidaldecomposition (e.g. Laplacian pyramid, wavelet, curvelet, etc), the solution ofDG3PD inpainting can be considered as a closed loop of the pyramid schemewith the true solution ( u ∗ , v ∗ , (cid:15) ∗ , g ∗ ) in the minimization problem. The loopof the updated Lagrange multipliers can be seen as a refinement. If we removethis loop from the scheme in Figure 12 (a), the minimization in DG3PD be-comes the quadratic penalty method. The removal would cause an imperfectreconstruction for the constraints, e.g. u + v + (cid:15) (cid:54) = f , see Figure 7 (p) and (q)in [1]. We note that the lowpass, bandpass, and highpass filters for cartoon,texture and residual components are ’custom-made’ for each image, resultingin different filters for different images, see Figure 16 and 15. The directionalproperties of texture and their directional filter banks are illustrated in Figure13 and 14. In this paper, we have addressed the very challenging task of restoring imageswith a large number of missing pixels, whose existing pixels are corrupted bynoise and importantly, the image to be restored contains both cartoon and tex-ture elements. The proposed DG3PD model for inpainting and denoising cancope with this threefold difficult problem. The task of simultaneous inpaintingand denoising for cartoon and texture components is solved by DG3PD decom-position, followed by inpainting and denoising both components separately andfinally, image restoration by synthesis of the restored components. More specif-ically, the DG3PD inpainting model is based on a regularization in the Banachspace in a discrete setting which is solved by ALM and ADMM.In summary, the decomposition step of the proposed method has two majoradvantages for tackling this challenging problem. Firstly, the decomposition stepdenoises simultaneously the cartoon and texture component. Secondly, it allowsto handle cartoon and texture in different ways. The cartoon image is inpaintedby DG3PD as described in Section 2 and 3. Separately, the texture componentis inpainted followed by further denoising as described in Section 4. Therefore,the proposed DG3PD decomposition can also be understood as breaking the fullproblem down into ’smaller’ subproblems (e.g. texture only inpainting) whichare easier to solve.Image restoration (or image denoising and inpainting) can be understood asan inverse problem and it can described by a discrete innovation model (seeFigure 9). The assumption in this procedure is that signal has a sparse repre-sentation in suitable transform domains. It is known from the probability theorythat the sparsity (by (cid:96) ) is connected to the Laplace distribution of the signal26hich results in the heavy tail distribution. Note that the Laplace distributionwith the (cid:96) norm is one of the best approximation of the sparsity by (cid:96) norm.Moreover, by choosing the priors and posterior according to the Bayesianframework and a maximum of a posterior (MAP), we can understand the se-lection of parameters (see Figure 8 for the results of the selection of the heavy-tailed distribution, e.g. the Laplace distribution). Note that MAP requires anassumption on the noise, e.g. Gaussian or Laplacian, to establish a minimiza-tion. In our proposed model, there is no requirement for an assumption onnoise distribution, e.g. independent identical distributed (i.i.d), or correlatedand Gaussian or non-Gaussian (due to the measurement (cid:13)(cid:13) C{·} (cid:13)(cid:13) (cid:96) ∞ , which issimilar to the Dantzig selector [55]). Acknowledgements
D.H. Thai is supported by the National Science Foundation under Grant DMS-1127914 to the Statistical and Applied Mathematical Sciences Institute.C. Gottschlich gratefully acknowledges the support of the Felix-Bernstein-Institutefor Mathematical Statistics in the Biosciences and the Niedersachsen Vorab ofthe Volkswagen Foundation.
References [1] D.H. Thai and C. Gottschlich. Directional global three-part image decom-position.
EURASIP Journal on Image and Video Processing , 2016(12):1–20, March 2016.[2] R.C. Gonzalez and R.E. Woods.
Digital Image Processing . Prentice Hall,Upper Saddle River, NJ, USA, 2002.[3] R. Szeliski.
Computer Vision: Algorithms and Applications . Springer,London, United Kingdom, 2011.[4] M. Sonka, V. Hlavac, and R. Boyle.
Image Processing, Analysis, and Ma-chine Vision . Thomson, Toronto, Canada, 2008.[5] C. Steger, M. Ulrich, and C. Wiedemann.
Machine Vision Algorithms andApplications . Wiley, Weinheim, Germany, 2008.[6] E.R. Davies.
Machine Vision: Theory, Algorithms, Practicalities . Aca-demic Press, Waltham, MA, USA, 2012.[7] J. Ma and G. Plonka. The curvelet transform.
IEEE Signal ProcessingMagazin , 27(2):118–133, March 2010.278] C. Gottschlich. Curved-region-based ridge frequency estimation and curvedGabor filters for fingerprint image enhancement.
IEEE Transactions onImage Processing , 21(4):2220–2227, April 2012.[9] C. Gottschlich and C.-B. Sch¨onlieb. Oriented diffusion filtering for enhanc-ing low-quality fingerprint images.
IET Biometrics , 1(2):105–113, June2012.[10] J.S. Bart˚unˇek, M. Nilsson, B. S¨allberg, and I. Claesson. Adaptive finger-print image enhancement with emphasis on preprocessing of data.
IEEETransactions on Image Processing , 22(2):644–656, February 2013.[11] S. Marukatat. Image enhancement using local intensity distribution equal-ization.
EURASIP Journal on Image and Video Processing , 2015(31),September 2015.[12] J. Yang, J. Wright, T.S. Huang, and Y. Ma. Image super-resolutionvia sparse representation.
IEEE Transactions on Image Processing ,19(11):2861–2873, November 2010.[13] L.A. Vese and S. Osher. Image denoising and decomposition with totalvariation minimization and oscillatory functions.
Journal of MathematicalImaging and Vision , 20(1-2):7–18, January 2004.[14] A. Buades, B. Coll, and J.M. Morel. A review of image denoising algo-rithms, with a new one.
Multiscale Modeling and Simulation , 4(2):490–530,July 2005.[15] S.G. Chang, B. Yu, and M. Vetterli. Adaptive wavelet thresholding forimage denoising and compression.
IEEE Transactions on Image Processing ,9(9):1532–1546, September 2000.[16] C.-B. Schoenlieb.
Modern PDE Techniques for Image Inpainting . PhDthesis, University of Cambridge, Cambridge, UK, June 2009.[17] F. Battisti, M. Carli, and A. Neri. Image forgery detection by using no-reference quality metrics. In
Proc. MWSF , pages 1–9, Burlingame, CA,USA, January 2012.[18] D. Salomon.
Data Compression . Springer, London, UK, fourth editionedition, 2007.[19] N. Joshi, S.B. Kang, C.L. Zitnick, and R. Szeliski. Image deblurring us-ing inertial measurement sensors. In
Proc. SIGGRAPH , pages 1–9, LosAngeles, CA, USA, July 2010.[20] L. Rudin, S. Osher, and E. Fatemi. Nonlinear total variation based noiseremoval algorithms.
Physica D , 60(1-4):259–268, November 1992.[21] J. Shen and T.F. Chan. Mathematical models for local nontexture inpaint-ings.
SIAM Journal on Applied Mathematics , 62(3):1019–1043, 2002.2822] T.F. Chan and J. Shen. Non-texture inpainting by curvature-driven diffu-sions (CDD).
Journal of Visual Communication and Image Representation ,12(4):436–449, December 2001.[23] J.-L. Starck, M. Elad, and D.L. Donoho. Image decomposition via thecombination of sparse representations and a variational approach.
IEEETransactions on Image Processing , 14(10):1570–1582, October 2005.[24] M. Elad, J.-L. Starck, P. Querre, and D.L. Donoho. Simultaneous car-toon and texture image inpainting using morphological component analysis(MCA).
Appl. Comput. Harmon. Anal. , 19:340–358, 2005.[25] S. Esedoglu and J. Shen. Digital inpainting based on the Mumford-Shah-Euler image model.
European J. Appl. Math , 13:353–370, 2002.[26] D. Mumford and J. Shah. Optimal approximations by piecewise smoothfunctions and associated variational problems.
Communications on Pureand Applied Mathematics , 42(5):577–685, July 1989.[27] L. Ambrosio and V.M. Tortorelli. Approximation of functional depend-ing on jumps by elliptic functional via Γ-convergence.
Comm. Pure Appl.Math. , 43(8):999–1036, December 1990.[28] E.D. Giorgi. Some remarks on Γ-convergence and least squares method.
Composite Media and Homogenization Theory , 5:135–142, 1991.[29] J. Shen, S.H. Kang, and T.F. Chan. Euler’s elastica and curvature-basedinpainting.
SIAM Journal on Applied Mathematics , 63(2):564–592, 2002.[30] C. Ballester, M. Bertalmio, V. Caselles, G. Sapiro, and J. Verdera. Filling-in by joint interpolation of vector fields and grey levels.
IEEE Transactionson Image Processing , 10(8):1200–1211, August 2001.[31] A.A. Efros and T.K. Leung. Texture synthesis by nonparametric sampling.In
Proc. ICCV , pages 1033–1038, Corfu, Greece, September 1999.[32] M. Bertalmio, L. Vese, G. Sapiro, and S. Osher. Simultaneous structureand texture image inpainting.
IEEE Transactions on Image Processing ,12(8):882–889, August 2003.[33] J.-F. Cai, R.H. Chan, and Z. Shen. Simultaneous cartoon and textureinpainting.
Inverse Problems and Imaging , 4(3):379–395, August 2010.[34] L.A. Vese and S. Osher. Modeling textures with total variation minimiza-tion and oscillatory patterns in image processing.
Journal of ScientificComputing , 19(1-3):553–572, December 2003.[35] J.-F. Aujol and A. Chambolle. Dual norms and image decomposition mod-els.
International Journal of Computer Vision , 63(1):85–104, June 2005.2936] D.H. Thai and C. Gottschlich. Global variational method for fingerprintsegmentation by three-part decomposition.
IET Biometrics , 5(2):120–130,June 2016.[37] D.H. Thai, S. Huckemann, and C. Gottschlich. Filter design and per-formance evaluation for fingerprint image segmentation.
PLoS ONE ,11(5):e0154160, May 2016.[38] G. Aubert and P. Kornprobst.
Mathematical Problems in Image Processing .Springer, New York, NY, USA, 2002.[39] O. Scherzer, M. Grasmair, H. Grossauer, M. Haltmeier, and F. Lenzen.
Variational methods in imaging . Springer, New York, NY, USA, 2009.[40] O. Scherzer, editor.
Handbook of Mathematical Methods in Imaging .Springer, New York, NY, USA, 2011.[41] P. Prandoni and M. Vetterli.
Signal Processing for Communications . CRCPress, Boca Raton, FL, USA, 2008.[42] P.J. Burt and E.H. Adelson. The laplacian pyramid as a compact imagecode.
IEEE Transactions on Communications , 31(4):532–540, April 1983.[43] E.P. Simoncelli and W.T. Freeman. The steerable pyramid: a flexiblearchitecture for multi-scale derivative computation. In
Proc. ICIP , pages444–447, Washington, DC, October 1995.[44] M. Unser, N. Chenouard, and D. Van De Ville. Steerable pyramids andtight wavelet frames in L ( R d ). IEEE Transactions on Image Processing ,20(10):2705–2721, October 2011.[45] I. Khalidov and M. Unser. From differential equations to the contruc-tion of new wavelet-like bases.
IEEE Transactions on Signal Processing ,54(4):1256–1267, April 2006.[46] S. Mallat.
A Wavelet Tour of Signal Processing . Academic Press, SanDiego, CA, USA, 2008.[47] E. Cand`es and D. Donoho. New tight frames of curvelets and optimalrepresentations of objects with piecewise singularities.
Communications onPure and Applied Mathematics , 57(2):219–266, February 2004.[48] A. Telea. An image inpainting technique based on the fast marchingmethod.
Journal of Graphics Tools , 9(1):25–36, 2004.[49] M. Bertalmio, A.L. Bertozzi, and G. Sapiro. Navier-stokes, fluid dynamics,and image and video inpainting. In
Proc. CVPR , pages 355–362, Kauai,HI, USA, December 2001. 3050] D.H. Thai.
Fourier and Variational Based Approaches for Fingerprint Seg-mentation . PhD thesis, University of Goettingen, Goettingen, Germany,January 2015.[51] E. Cand`es, L. Demanet, D. Donoho, and L. Ying. Fast discrete curvelettransforms.
Multiscale Model. Simul. , 5(3):861–899, September 2006.[52] S. Kotz, T.J. Kozubowski, and K. Podg´orski.
The Laplace Distribution andGeneralizations . Springer, 2001.[53] S. Kotz and S. Nadarajah.
Multivariate T-Distributions and Their Appli-cations . Cambridge University Press, 2004.[54] M. Haltmeier and A. Munk. Extreme value analysis of empirical framecoefficients and implications for denoising by soft-thresholding.
Appliedand Computational Harmonic Analysis , 36(3):434–460, May 2014.[55] E. Cand`es and T. Tao. The Dantzig selector: Statistical estimation when p is much larger than n . The Annals of Statistics , 35(6):2313–2351, December2007.[56] M. Unser. Sampling - 50 years after Shannon.
Proceedings of the IEEE ,88(4):569–587, April 2000.[57] M.N. Do and M. Vetterli. The contourlet transform: An efficient direc-tional multiresolution image representation.
IEEE Transactions on ImageProcessing , 14(12):2091–2106, December 2005.[58] M. Unser and D. Van De Ville. Wavelet steerability and the higher-orderRiesz transform.