Simultaneous inpainting and denoising by directional global three-part decomposition: connecting variational and Fourier domain-based image processing

1. Introduction and related work

Image enhancement and image restoration are two superordinate concepts in image processing which encompass a plethora of methods to solve a multitude of important real-world problems [1,2]. Image enhancement has the goal of improving an input image for a specific application, e.g. in areas such as medical image processing, biometric recognition, computer vision, optical character recognition, texture recognition or machine inspection of surfaces [3–5]. Methods for image enhancement can be grouped by the domain in which they perform their operations: images are processed in the spatial domain or Fourier domain, or modified, e.g. in the wavelet or curvelet domain [6]. The types of enhancement methods include contextual filtering, e.g. for fingerprint image enhancement [7–9], contrast enhancement, e.g. by histogram equalization [10], and image super-resolution [11]. Image restoration is connected to the notion that a given input image suffers from degradation and the goal is restore an ideal version of it. Degradations are caused by various types of noise, missing pixels or blurring and their countermeasures are denoising, inpainting and deblurring. In general, one has to solve a linear or nonlinear inverse problem to reconstruct the ideal image from its given degraded version. Denoising aims to remove noise from an image and denoising methods include total variation (TV) minimization-based approaches [12], the application of non-local means (NL means) [13] or other dictionaries of image patches for smoothing, and adaptive thresholding in the wavelet domain [14]. Inpainting [15] is the filling-in of missing pixels from the available information in the image, and it is applied for scratch removal from scanned photographs, for occlusion filling, for removing objects or persons from images (in image forgery [16] or for special effects), and for filling-in of pixels which were lost during the transmission of an image or left out on purpose for image compression [17]. Deblurring [18] addresses the removal of blurring artefacts and is not the focus of this paper.

Rudin et al. [19] pioneered two-part image decomposition by TV regularization for denoising. Shen & Chan [20] applied TV regularization to image inpainting, called the TV inpainting model, and they also suggested image inpainting by curvature-driven diffusions [21]. Starck et al. [22] defined a model for two-part decomposition based on the dictionary approach. Then, Elad et al. [23] applied this decomposition idea for image inpainting by introducing the indicator function in the ℓ₂ norm of the residual; see eqn (6) in [23]. Esedoglu & Shen [24] introduced two inpainting models based on the Mumford–Shah model [25] and its higher order correction—the Mumford–Shah–Euler image model. They also presented numerical computation based on the Γ-convergence approximations [26,27]. Shen et al. [28] proposed image inpainting based on bounded variation and elastica models for non-textured images.

Image inpainting can be an easy or difficult problem depending on the amount of missing pixels [21], the complexity of the image content and whether prior knowledge about the image content is available. Methods have been proposed which perform only cartoon inpainting (also referred to as structure inpainting) [20,28,29] or only texture inpainting [30]. Images which consist of both cartoon (structure) and texture components are more challenging to inpaint. Bertalmio et al. [31], Elad et al. [23] and Cai et al. [32] have proposed methods for inpainting which can handle images with both cartoon (structure) and texture components.

In this paper, we tackle an even more challenging problem. Consider an input image 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 and texture components increases the difficulty of both the inpainting problem and the denoising problem. A multitude of methods have been proposed for inpainting and denoising. Existing inpainting methods in the literature typically assume that the non-missing pixels in a given image contain only a small amount of noise or are noise free, and existing methods for denoising typically assume that all pixels of the noisy image are known. The proposed method for solving this challenging problem is inspired by the works of Efros & Leung [30], Bertalmio et al. [31], Vese & Osher [33], Aujol & Chambolle [34], Buades et al. [13] and Elad et al. [23], and is based on the directional global three-part decomposition (DG3PD) [35]. The DG3PD method decomposes an image into three parts: a cartoon image, a texture image and a residual image. Advantages of the DG3PD model lie in the properties which are enforced on the cartoon and texture images. The geometric objects in the cartoon image have a very smooth surface and sharp edges. The texture image yields oscillating patterns on a defined scale which is both smooth and sparse. Recently, the texture images have been applied as a very useful feature for fingerprint segmentation [35–37].

We address the challenging task of simultaneous inpainting and denoising in the following way. The advanced DG3PD model introduced in the next section decomposes a noisy input image f (with missing regions D) into cartoon u, texture v and residual ϵ 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. This effect benefits from the help of the indicator function in the measurement of the residual, i.e. $\parallel \mathcal{C}\{{\chi }_D^c{\cdot }^{\times }\mathit{ϵ}\}{\parallel }_{{\ell }_{\mathrm{\infty }}}$ in (2.1). However, texture v is not interpolated due to the ‘cancelling’ effect of this supremum norm for residual in unknown regions. Therefore, the texture component v is inpainted and denoised by a dictionary-based approach instead. The DG3PD decomposition drives noise into the residual component ϵ, which is discarded. The reconstruction of the ideal version of f is obtained by summation of the inpainted and denoised cartoon and texture components (see figure 1 for a visual overview).

• Download figure
• Open in new tab
• Download PowerPoint

Moreover, we uncover the link between the calculus of variations [38–40] and filtering in the Fourier domain [41] by analysing the solution of the convex minimization in equation (2.1), i.e. roughly speaking the solution of the DG3PD inpainting model which can be understood as the response of the lowpass filter $\widehat{\mathrm{LP}}(\mathit{\omega })$, highpass filter $\widehat{\mathrm{HP}}(\mathit{\omega })$ and bandpass filter $\widehat{\mathrm{BP}}(\mathit{\omega }),$ and the unity condition is satisfied, i.e.

where ω∈[ − π, π]² is a coordinator in the Fourier domain. We observe that this decomposition is similar to the wavelet or pyramidal decomposition scheme [42–44]. However, the basis elements obtaining the decomposition, i.e. the scaling function and frame (or wavelet-like) function, are constructed by discrete differential operators (due to the discrete setting in minimizing (2.1)), which are referred to as wavelet-like operators in [45]. In particular,

• — the scaling function and wavelet-like function for the cartoon u are from the effect of the multi-directional TV norm,

• — the scaling function and wavelet-like function to extract the texture v are reconstructed by the effect of the multi-directional G-norm,

• — the effect of the ℓ_∞ norm $\parallel \mathcal{C}\{{\chi }_D^c{\cdot }^{\times }\mathit{ϵ}\}{\parallel }_{{\ell }_{\mathrm{\infty }}}$ is to remove the remaining signal in the known regions of the residual ϵ (due to the duality property of ℓ_∞).

We also describe flowcharts to show that the method of variational calculus (or the DG3PD inpainting) is a closed loop pyramidal decomposition, which is different from an open loop one, e.g. wavelet [46], curvelet [47], see §7. By numerics, we observe that the closed loop filter design by the calculus of variation will result in lowpass, highpass and bandpass filters which are ‘unique’ for different images (figure 17). We also analyse the DG3PD inpainting model from the perspective of a Bayesian framework and then define a discrete innovation model for this inverse problem.

This paper is organized as follows. In §2, we describe the DG3PD model for image inpainting and denoising. In §3, we show how to compute the solution of the convex minimization in the DG3PD inpainting problem by the augmented Lagrangian method. In §4, we describe the proposed method for texture inpainting and denoising. In §5, we compare the proposed method with existing ones (TVL2 inpainting [48,49]). Moreover, in order to enhance our evaluation in terms of simultaneous inpainting and denoising effects, we also perform a comparison with the NL-means filter for denoising, namely block-matching and three-dimensional filtering (BM3D), in [50]. In §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 and texture v. We analyse the relation between the calculus of variations and the traditional pyramid decomposition scheme, e.g. Gaussian pyramid, in §7. We conclude the study with §8. For more detailed notation and mathematical preliminaries, we refer the reader to [36,51].

2. Inpainting by DG3PD

We define a method for restoration of the original noisy/non-noisy image f with a set of known missing regions D. The proposed model is a generalized version of DG3PD [35] for the inpainting and denoising problem. This modification for our DG3PD inpainting is inspired by Elad et al. [23].

In particular, the discrete image f (size m × n) with missing regions D is simultaneously decomposed into cartoon u, texture v and noise ϵ; in particular, 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 indicator function χ_D whose complement is

Instead of putting χ^c_D on the ℓ₂ norm of the residual (see eqns (4) and (6) in [23]), we introduce χ^c_D for the residual in the DG3PD model with the point-wise multiplication operator · × , i.e. $\parallel \mathcal{C}\{{\chi }_D^c{\cdot }^{\times }\mathit{ϵ}\}{\parallel }_{{\ell }_{\mathrm{\infty }}}$. We propose the DG3PD inpainting as follows:

where $\mathcal{C}$ is a curvelet transform [47], g_s is an unknown variable in the G-norm [52] and a finite difference matrix is and ${\mathbf{D}}_{\mathbf{n}}\in {\mathbb{R}}^{n\times n}$ is similar.

Note that if there are no missing regions, i.e. χ^c_D[k] = 1, ∀k∈Ω, the minimization problem (2.1) becomes the DG3PD model [35]. Next, we discuss how to solve the DG3PD-inpainting model (2.1).

3. Solution of the DG3PD inpainting model

In this section, we present a numerical algorithm for obtaining the solution of the DG3PD-inpainting model stated in (2.1). To simplify the calculation, we introduce three new variables

and denote G*(e/ν) as the indicator function on the feasible convex set A(ν) of (2.1), i.e. Owing to the constrained minimization problem in (2.1), the augmented Lagrangian method is applied to turn it into an unconstrained one as where the Lagrange function is Similar to [35], the alternating directional method of multipliers is applied to solve the unconstrained minimization problem (3.1). Its minimizer is numerically computed through iterations, t = 1, 2, …,

and the Lagrange multipliers are updated after every step t. We also initialize u⁽⁰⁾ = f, v⁽⁰⁾ = ϵ⁽⁰⁾ = e⁽⁰⁾ = [r⁽⁰⁾_l]^L−1_l=0 = [w⁽⁰⁾_s]^S−1_s=0 = [g⁽⁰⁾_s]^S−1_s=0 = [λ⁽⁰⁾_1l]^L−1_l=0 = [λ⁽⁰⁾_2a]^S−1_a=0 = λ⁽⁰⁾₃ = λ⁽⁰⁾₄ = λ⁽⁰⁾₅ = 0, where 0 is an m × n zero matrix.

In each iteration, we first solve the seven subproblems in the listed order: ‘[r_l]^L−1_l=0-problem’, ‘[w_s]^S−1_s=0-problem’, ‘[g_s]^S−1_s=0-problem’, ‘v-problem’, ‘u-problem’, ‘e-problem’, ‘ϵ-problem’, and then we update the five Lagrange multipliers, namely [λ_1l]^L−1_l=0, [λ_2a]^S−1_a=0, λ₃, λ₄, λ₅ (see algorithm 1).

In this section, we only present the solution of the ‘e-problem’ and ‘ϵ-problem’ as well as the updated Lagrange multiplier λ₅. We refer readers to [35] for a detailed explanation of the other subproblems and the other Lagrange multipliers.

The e-problem: Fix u, v, ϵ, [r_l]^L−1_l=0, [w_s]^S−1_s=0, [g_s]^S−1_s=0 and

The solution of (3.3) is defined as a projection operator on a convex set A(ν), i.e. ${\mathbb{P}}_{A(\nu )}({\chi }_D^c{\cdot }^{\times }\mathit{ϵ}-{\mathbf{\lambda }}_{\mathbf{5}}/{\beta }_5)$:

The ϵ-problem: Fix u, v, e, [r_l]^L−1_l=0, [w_s]^S−1_s=0, [g_s]^S−1_s=0 and

The solution of (3.5) with the point-wise division operator /_. is and 1_mn is an m × n matrix of ones.

Updated Lagrange multiplierλ₅∈X:

Choice of parameters

Owing to an adaptability to specific images and a balance between the smoothing terms and updated terms for the solution of the above subproblems, the selection of parameters (μ₁, μ₂, β₁, β₂, β₃) is described in [35] and β₅ is defined as

The term ${\parallel \mathcal{C}\{{\chi }_D^c{\cdot }^{\times }\mathit{ϵ}\}\parallel }_{{\ell }_{\mathrm{\infty }}}$ in (2.1) serves as a good measurement for the amount of noise due to the substraction in equation (3.4) [35,36,53] (owing to the multi-scale and multi-direction properties of the curvelet transform and its duality). However, because of the indicator function χ^c_D in the noise measurement ${\parallel \mathcal{C}\{{\chi }_D^c{\cdot }^{\times }\mathit{ϵ}\}\parallel }_{{\ell }_{\mathrm{\infty }}}$, the interpolated texture by DG3PD inpainting is almost in the residual; see figure 2a, and its smoothed version by curvelet shrinkage (figure 2b). Because of the substraction operator in equation (3.4), these interpolated textures are cancelled out in e (figure 2c). Figure 2d illustrates the estimated texture before the substraction in equation (3.4) as

where ROI is the region of interest obtained by the morphological operator for the texture image in [37].

• Download figure
• Open in new tab
• Download PowerPoint

In order to analyse the effects of the proposed model in terms of interpolation (for u) and decomposition, we consider a one-dimensional signal (which is extracted along the red line in figure 3a–d). By the DG3PD inpainting, the mean values of the cartoon u in (figure 3e) remain almost unchanged. The homogeneous regions have ‘sharp’ edges and a ‘staircase’ effect does not occur (owing to the directional TV). Texture v in (figure 3f) is extracted in areas which contain a repeated pattern (figure 3b). Moreover, small-scale objects, e.g. noise, are removed from v.

• Download figure
• Open in new tab
• Download PowerPoint

However, the term ${\parallel \mathcal{C}\{{\chi }_D^c{\cdot }^{\times }\mathit{ϵ}\}\parallel }_{{\ell }_{\mathrm{\infty }}}$ causes a cancelling effect which removes the interpolated texture in an unknown area during the decomposition process. Therefore, texture inpainting and denoising are tackled separately, as described in the next section, by a generalized version of the dictionary learning for texture synthesis [30].

4. Texture inpainting and denoising

The proposed method for reconstructing the texture component combines the following ideas in five subsequent steps:

• — texture extraction by DG3PD [35],

• — morphological processing [36,37],

• — texture inpainting inspired by the work of Efros & Leung [30] or [54],

• — texture denoising by NL means [13],

• — image synthesis by summation with the cartoon component u.

The DG3PD model enforces sparsity and smoothness on the obtained texture component v. Sparsity means that the vast majority of coefficients in v are equal to zero (the percentage of zero coefficients depends on how much texture a specific image contains). We make good use of this property to answer the following question: which pixels of those that are missing shall be inpainted with texture? First, the texture component v is segmented into zero (grey pixels in figure 4a) and non-zero coefficients (black and white pixels in figure 4a). Morphological processing (as described in [36,37] for fingerprint segmentation) obtains texture regions (figure 4b). The mask of missing pixels D shown in figure 4c is dilated (we used a circular structure element with a radius of 5 pixels to avoid border effects on the margin between existing and missing pixels) and then it is intersected with texture region segmentation to obtain the mask I shown in figure 4d, which shows the pixels to be inpainted with texture.

• Download figure
• Open in new tab
• Download PowerPoint

The proposed texture inpainting proceeds in two phases. First, we build a dictionary of texture patches and, second, we inpaint all pixels of mask I in a specific 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 following two criteria: (i) at least p₁ per cent of the pixels are known (also the central pixel of the patch has to exist) and (ii) at least p₂ per cent of the coefficients are non-zero (we used p₁ = 70% and p₂ = 60% in our experiments). The first criterion excludes patches which contain too many missing (unknown) pixels and 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 inpainted as the central pixel of an image patch with size s × s. We count the number of 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%; this is decreased in steps of 5% per iteration and, in each iteration, we inpaint all pixels of mask I for which at least t per cent are known. The rationale behind this ordering is that the more pixels are known (or already inpainted) in the neighbourhood around a missing pixel, the better it can inpainted, because this additional information improves the chances of finding a good match in the texture dictionary. A third constraint ensures that the overlap of known pixels between a dictionary patch and the patch around the missing pixels contains at least p₃ per cent of known pixels (we used p₃ = 30% in our experiments). Inpainting a pixel means that we find the best fitting patch in our texture dictionary and set the pixel value to the central pixel of that patch. Best fitting is defined as the minimum sum of squared differences per pixel, divided by the number of pixels which overlap.

After all missing texture pixels have been inpainted, the texture region is denoised by n iterations of NL means. In each iteration, we first construct a dictionary considering all patches in the texture region with known and previously 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 of squared pixel-wise differences) and set the denoised pixel to the average value of the top k central pixel values. In our experiments, we used k = 5 and we have observed that, after about n = 10 iterations, the image has reached a steady state. See figure 5 for a visualization of the status at several intermediate steps during the inpainting and denoising process.

• Download figure
• Open in new tab
• Download PowerPoint

Finally, the full image is reconstructed by summation of the inpainted cartoon component u with the inpainted and denoised texture component v (figure 10).

5. Comparison of DG3PD with further inpainting and denoising methods

Figure 6 shows inpainting results for the challenging problem under consideration (figure 7c) obtained by well-known inpainting methods from the literature: an approach based on Navier–Stokes equations from fluid dynamics which was proposed by Bertalmio et al. [49], an inpainting method suggested by Telea [48] and a well-known TVL2 inpainting approach with its synthesis image shown in figure 6b,c.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

The image shown in figure 7c has the three properties discussed in §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 7a is available which serves as the ground truth for evaluating the quality of inpainted and denoised images by different methods. The availability of a noise-free image is an advantage for comparing and evaluating different inpainting and denoising methods. By contrast, images taken by a digital camera contain a certain amount of noise and, for them, an ideal version is not available, see, for example, the Barbara image we used in previous comparisons [35]. The choice of the ideal image (shown in figure 7) enables us to clearly see and understand the advantages and limitations of the compared methods (figure 6).

Owing to the fidelity term, the L₂-norm, in the TVL2 inpainting model, noise is reduced and the homogeneous areas are well interpolated. However, the method is known to cause a ‘staircase’ effect on the homogeneous regions and texture, while small-scale objects tend to be eliminated. If the parameter β is chosen smaller (e.g. β₅ = 0.005), the resulting inpainted image is smoother, i.e. the ‘staircase’ effect is reduced, but also more texture parts are removed.

DG3PD inpainting produces a smooth result without the staircase effect. In comparison with the other approaches, it is the only method which can reconstruct the texture regions in a satisfactory way.

In order to evaluate the proposed model in terms of image denoising, we compared it with BM3D in [50]. This NL-mean filter enhances the sparse representation of an image in a transform domain by grouping similar two-dimensional image blocks into three-dimensional ‘groups’. Collaborative filtering is applied for these three-dimensional groups by using three successive steps: three-dimensional transformation of a group, shrinkage of the transform spectrum, and inverse three-dimensional transformation. Finally, aggregation (i.e. an averaging procedure) is applied for different estimates of each pixel due to several overlapped blocks.

BM3D is a state-of-the-art method for image denoising. However, it is not designed for inpainting and, in order to be able to run the BM3D code available at http://www.cs.tut.fi/foi/GCF-BM3D/, we have to provide input values for the missing pixels; therefore, we

• — set missing pixels in a corrupted image in figure 7c as uniform random values in [0, 255] (figure 8a),

• — let standard deviation of noise in known regions σ = 100, and

• — linearly rescale an input image to [0, 1].

• Download figure
• Open in new tab
• Download PowerPoint

In the end, we realized that, while BM3D can separate noise reasonably well from a corrupted image in figure 8a, the resultant image illustrated in figure 8b does not have a good performance. In particular,

• — from an inpainting aspect, BM3D cannot interpolate mixing areas on homogeneous or textured regions and

• — from a denoising aspect, reconstructed areas by BM3D are blurred without sharp edges and with a loss of contrast as well as with several artefacts on homogeneous regions.

6. From the Bayesian framework to the DG3PD inpainting model

Following [55], in this section, we describe the DG3PD inpainting model from the perspective of a Bayesian framework through a maximum a posteriori (MAP) estimator. We assume that u and v are independent random variables (r.v.).

• — Prior of cartoon image u: A cartoon u = [u[k]]_k∈Ω consists of an element pixel u[k] which is considered as an independent r.v. with a prior p_Uk.

  Given k∈Ω, denote r[k] = ∇⁺_Lu[k] = [r_l[k]]^L−1_l=0. To describe the distribution of a r.v. U_k, we firstly define the distribution of an L-variate r.v. R_k = [R_0,k, …, R_L−1,k]. We assume that r.v. R_k has a multi-variate Laplace distribution which is a part of multi-variate power exponential distribution [56–58], i.e. ${\mathbf{R}}_k\sim {\mathrm{PE}}_L(0,\mathit{\Sigma },\frac{1}{2})={\mathrm{Laplace}}_L(0,\mathit{\Sigma })$ with $\mathit{\Sigma }=\left(\begin{array}{ccc}{\gamma }^2& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & {\gamma }^2\end{array}\right)\in {\mathbb{R}}_{+}^{L\times L}\mathrm{\setminus }\{0\}$ and Γ( · ) function:

  $$\begin{array}{l}{p}_{R_k}\left(r[k]\mid 0,\Sigma ,\frac{1}{2}\right)\\ =\frac{\Gamma (L/2)}{{\pi }^{L/2}\Gamma (L)2^{1+L}{\gamma }^L}\exp \left\{-\frac{1}{2}{\left({\left(\begin{array}{c}{r}_0[k]\\ ⋮\\ r_{L-1}[k]\end{array}\right)}^{\text{T}}\frac{1}{{\gamma }^2}\left(\begin{array}{ccc}1& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & 1\end{array}\right)\left(\begin{array}{c}{r}_0[k]\\ ⋮\\ r_{L-1}[k]\end{array}\right)\right)}^{1/2}\right\}\\ =\frac{\Gamma (L/2)}{{\pi }^{L/2}\Gamma (L)2^{1+L}{\gamma }^L}\underset{:=|r[k]|=|{\nabla }_L^{+}u[k]|}{\underbrace{{\left[{\displaystyle \sum _{l=0}^{L-1}|}{r}_l[k]{|}^2\right]}^{1/2}}}.\end{array}$$
  Thus, the distribution of a r.v. U_k is a multi-dimensional Laplace distribution with operator ∇⁺_L:
  $$p_{U_k}(u[\mathit{k}])=\frac{\mathit{\Gamma }(L/2)}{{\pi }^{L/2}\mathit{\Gamma }(L)2^{1+L}{\gamma }^L}\exp \{-\frac{1}{2\gamma }|{\mathrm{\nabla }}_L^{+}u[\mathit{k}]|\}.$$
  The joint probability density function (p.d.f.) of a prior u is
  $$p_U(\mathbf{\text{u}})=\prod _{\mathit{k}\in \mathit{\Omega }}{p}_{U_k}(u[\mathit{k}])={\left[\frac{\mathit{\Gamma }(L/2)}{{\pi }^{L/2}\mathit{\Gamma }(L)2^{1+L}{\gamma }^L}\right]}^{|\mathit{\Omega }|}{e}^{-(1/2\gamma )}{\parallel {\mathrm{\nabla }}_L^{+}\mathbf{\text{u}}\parallel }_{{\ell }_1}.$$
  We choose $\gamma =\frac{1}{2}$. The potential function of u in matrix form is
  $${\mathit{\Phi }}_U(\mathbf{\text{u}})=-\log p_U(\mathbf{\text{u}})={\parallel {\mathrm{\nabla }}_L^{+}\mathbf{\text{u}}\parallel }_{{\ell }_1}+|\mathit{\Omega }|\log \frac{2{\pi }^{L/2}\mathit{\Gamma }(L)}{\mathit{\Gamma }(L/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 Laplace distribution.

• — Prior of texture image v: Under a transform operator $\mathcal{T}$, the texture image v is decomposed in different orientations. Operator $\mathcal{T}$ is suitably chosen to capture the texture components in the original image f. As for the definition of a discrete multi-dimensional G-norm [35], the transform $\mathcal{T}$ (in a direction s, s = 0, …, S − 1) is defined by its inversion (note that $\mathcal{T}{\mathcal{T}}^{-1}=\mathrm{Id}$):

  $$\mathcal{T}\{\mathbf{v}\}={[\underset{:={\mathbf{g}}_s}{\underbrace{{\mathcal{T}}_s\{\mathbf{v}\}}}]}_{s=0}^{S-1}=\mathbf{g}\iff \mathbf{v}={\mathcal{T}}^{-1}\{{[{\mathbf{g}}_s]}_{s=0}^{S-1}\}={\mathrm{div}}_S^{+}\mathbf{g}.$$
  We assume that a texture v is sparse in a transform domain $\mathcal{T}$ and sparse in the spatial domain (non-zero coefficients of v are only due to texture). To satisfy these conditions, we assume that a r.v. V_k has a mixture of a Laplace distribution in a spatial domain with a p.d.f. p_Vk,2(v[k]) (w.r.t. sparse in the spatial domain) and a multi-variate Laplace distribution (in an operator $\mathcal{T}$) with a p.d.f. p_Vk,1(v[k]) (w.r.t. sparse in the transform domain $\mathcal{T}$). Thus, we have a p.d.f. of r.v. V_k as follows:
  $$p_{V_k}(v[\mathit{k}])=p_{V_{1k}}(v[\mathit{k}])\cdot p_{V_{2k}}(v[\mathit{k}]).$$6.1
  Since V_2k∼Laplace(0, γ²₁), we have
  $$p_{V_{2k}}(v[\mathit{k}])=\frac{1}{4{\gamma }_2}\exp \{-\frac{1}{2{\gamma }_2}\left|v[\mathit{k}]\right|\}.$$6.2
  To define the distribution of a r.v. V_1k, we define the distribution of the S-dimension r.v. G_k = [G_0,k, …, G_S−1,k]. Similar to U_k, we assume that r.v. G_k has a multi-variate Laplace distribution, i.e. G_k∼Laplace_S(0, Σ) with $\mathit{\Sigma }=\left(\begin{array}{ccc}{\gamma }_1^2& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & {\gamma }_1^2\end{array}\right)\in {\mathbb{R}}_{+}^{S\times S}\mathrm{\setminus }\{0\}$. It is easy to see that (with $g[\mathit{k}]=\mathcal{T}\{\mathbf{v}\}[\mathit{k}],\mathit{k}\in \mathit{\Omega }$)
  $$p_{{\mathbf{\text{G}}}_{\mathbf{\text{k}}}}(g[\mathit{k}])=\frac{\mathit{\Gamma }(S/2)}{{\pi }^{S/2}\mathit{\Gamma }(S)2^{1+S}{\gamma }_1^S}\exp \{-\frac{1}{2{\gamma }_1}\left|g[\mathit{k}]\right|\}$$
  or
  $$p_{V_{1k}}(v[\mathit{k}])=\frac{\mathit{\Gamma }(S/2)}{{\pi }^{S/2}\mathit{\Gamma }(S)2^{1+S}{\gamma }_1^S}\exp \{-\frac{1}{2{\gamma }_1}\left|\mathcal{T}\{\mathbf{\text{v}}\}[\mathit{k}]\right|\}.$$6.3
  Putting (6.2) and (6.3) into (6.1), we have the p.d.f. of r.v. V_k as
  $$p_{V_k}(v[\mathit{k}])=\frac{\mathit{\Gamma }(S/2)}{{\pi }^{S/2}\mathit{\Gamma }(S)2^{3+S}{\gamma }_1^S{\gamma }_2}\exp \{-\frac{1}{2{\gamma }_1}\left|\mathcal{T}\{\mathbf{\text{v}}\}[\mathit{k}]\right|-\frac{1}{2{\gamma }_2}\left|v[\mathit{k}]\right|\}.$$
  The joint p.d.f. of a texture image v = [v[k]]_k∈Ω is
  $$p_V(\mathbf{\text{v}})=\prod _{\mathit{k}\in \mathit{\Omega }}{p}_{V_k}(v[\mathit{k}])={\left[\frac{\mathit{\Gamma }(S/2)}{{\pi }^{S/2}\mathit{\Gamma }(S)2^{3+S}{\gamma }_1^S{\gamma }_2}\right]}^{|\mathit{\Omega }|}\exp \{-\frac{1}{2{\gamma }_1}{\parallel \mathcal{T}\{\mathbf{\text{v}}\}\parallel }_{{\ell }_1}-\frac{1}{2{\gamma }_2}{\parallel \mathbf{\text{v}}\parallel }_{{\ell }_1}\},$$
  and we choose μ₂ = 1/2γ₂ and μ₁ = 1/2γ₁. The potential function of v with an anisotropic version $(\parallel \mathcal{T}\{\mathbf{v}\}{\parallel }_{{\ell }_1}=\sum _{s=0}^{S-1}\parallel {\mathcal{T}}_s\{\mathbf{v}\}{\parallel }_{{\ell }_1}=\sum _{s=0}^{S-1}\parallel {\mathbf{g}}_{\mathbf{s}}{\parallel }_{{\ell }_1})$ in matrix form is defined as
  $${\mathit{\Phi }}_V(\mathbf{\text{v}})={\mu }_1\sum _{s=0}^{S-1}{\parallel {\mathbf{\text{g}}}_{\mathbf{\text{s}}}\parallel }_{{\ell }_1}+{\mu }_2\parallel \mathbf{\text{v}}{\parallel }_{{\ell }_1}+|\mathit{\Omega }|\log \frac{4{\pi }^{S/2}\mathit{\Gamma }(S)}{\mathit{\Gamma }(S/2){\mu }_1^S{\mu }_2}.$$

• — Likelihood (or the joint p.d.f. of p_F|U,V (f|u, v)): If we assume that the residual ϵ in an image f has a power exponential distribution [56,57], i.e. E_k∼PE(μ, σ², ξ), e.g. Gaussian (ξ = 1) or Laplacian $(\xi =\frac{1}{2})$, its density function at k∈Ω is

  $$p_{E_k}(ϵ[\mathit{k}]\mid \mu )=\frac{1}{\mathit{\Gamma }(1+1/2\xi )2^{1+1/2\xi }\sigma }\exp \{-\frac{1}{2}\frac{{|ϵ[\mathit{k}]-\mu |}^{2\xi }}{{\sigma }^{2\xi }}\}.$$
  Owing to the inpainting problem with a missing region D, the likelihood is defined on a known domain Ω\D as
  $$\begin{array}{c}{p}_{F\mid U,V}(\mathbf{f}\mid \mathbf{u},\mathbf{v})=\prod _{\mathbf{k}\in \mathit{\Omega }}{p}_{F_k\mid U_k,V_k}(f[\mathbf{k}]\mid u[\mathbf{k}],v[\mathbf{k}])\\ =\frac{1}{{(\mathit{\Gamma }(1+1/2\xi )2^{1+1/2\xi }\sigma )}^{|\mathit{\Omega }|}}\exp \left\{-\frac{1}{2{\sigma }^{2\xi }}\underset{:{\parallel ={\chi }_D^c{\cdot }^{\times }(\mathbf{f}-\mathbf{u}-\mathbf{v})\parallel }_{2\xi }^{2\xi }}{\underbrace{\sum _{\mathbf{k}\in \mathit{\Omega }}{\left|{\chi }_D^c[\mathbf{k}](f[\mathbf{k}]-u[\mathbf{k}]-v[\mathbf{k}])\right|}^{2\xi }}}\right\}.\end{array}$$

• — A posteriori: Since u and v are independent of each other, a posteriori is written as

  $$p_{U,V\mid F}(\mathbf{u},\mathbf{v}\mid \mathbf{f})=\frac{p_{F\mid U,V}(\mathbf{f}\mid \mathbf{u},\mathbf{v})p_U(\mathbf{u})p_V(\mathbf{v})}{p_F(\mathbf{f})}\propto p_{F\mid U,V}(\mathbf{f}\mid \mathbf{u},\mathbf{v})p_U(\mathbf{u})p_V(\mathbf{v}).$$

Let ϵ = f − u − v, and the MAP estimator, i.e. $\underset{(\mathbf{u},\mathbf{v})\in X^2}{max}{p}_{U,V\mid F}(\mathbf{u},\mathbf{v}\mid \mathbf{f})$, is defined as

However, in practice, we do not know which types of noise are observed in a signal. Instead of the ℓ_2ξ norm of the residual ϵ on a known domain Ω\D to characterize the properties of noise, we control the smoothness of the solution by the maximum of the curvelet coefficient of ϵ on a domain Ω\D by a constant ν, i.e. ${\parallel \mathcal{C}\{{\chi }_D^c{\cdot }^{\times }\mathit{ϵ}\}\parallel }_{{\ell }_{\mathrm{\infty }}}\le \nu$.

In [59], if noise in image f is Gaussian, as in extreme value theory, the r.v. ${\parallel \mathcal{C}\{{\chi }_D^c{\cdot }^{\times }\mathit{ϵ}\}\parallel }_{{\ell }_{\mathrm{\infty }}}$ 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 ${\parallel \mathcal{C}\{{\chi }_D^c{\cdot }^{\times }\mathit{ϵ}\}\parallel }_{{\ell }_{\mathrm{\infty }}}\le \nu$ is similar to the Dantzig selector [60].

Figure 9 illustrates the empirical density functions of the solution of the minimization (3.2). The statistical properties of the solution can be characterized in the Bayesian framework with the priors, likelihood and posterior as mentioned above.

• Download figure
• Open in new tab
• Download PowerPoint

Note that the DG3PD inpainting model (2.1) can be described as an inverse problem in image restoration (e.g. image inpainting and denoising); see figure 10 for the discrete innovation model.

• Download figure
• Open in new tab
• Download PowerPoint

7. Relation between variational analysis and pyramid decomposition

In this section, we discuss connections and similarities between the proposed DG3PD inpainting model and existing work in the area of pyramid representations, multi-resolution analysis and scale-space representation. Historically, Gaussian pyramids and Laplacian pyramids have been developed for applications such as texture synthesis, image compression and denoising. Early works in this area include a paper by Burt & Adelson [42] applying Laplacian pyramids for image compression. Pyramidal decomposition schemes [42–44] can be grouped into two categories. Their main property is either lowpass or bandpass filtering. Closely related are transforms from multi-resolution analysis, including wavelet, curvelet [47,53], contourlet [61] and steerable wavelet [62]. A difference in the proposed DG3PD inpainting model is that its basis elements for obtaining the decomposition are constructed by discrete differential operators and the decomposition is solved in a nonlinear way which adapts to each image (enabled by update iterations, corresponding to the loop in figure 15). The discussion of commonalities and differences is made more precise and detailed in the rest of the section.

Let us denote $\stackrel{\mathcal{F}}{⟷}$ as a discrete Fourier transform pair. Given z₁ = e^jω₁ and z₂ = e^jω₂ and the Dirac delta function δ( · ), the impulse responses of the discrete directional derivative operators (l = 0, …, L − 1) and their spectra are

• — the forward operator:

  $$\begin{array}{rl}{\mathrm{\partial }}_l^{+}\delta [\mathit{k}]& =[\cos \left(\frac{\pi l}{L}\right){\mathrm{\partial }}_x^{+}+\sin \left(\frac{\pi l}{L}\right){\mathrm{\partial }}_y^{+}]\delta [\mathit{k}]\\ & \stackrel{\mathcal{F}}{⟷}\cos \left(\frac{\pi l}{L}\right)(z_2-1)+\sin \left(\frac{\pi l}{L}\right)(z_1-1),\end{array}$$

• — the backward operator:

  $$\begin{array}{rl}{\mathrm{\partial }}_l^{-}\delta [\mathit{k}]& =[\cos \left(\frac{\pi l}{L}\right){\mathrm{\partial }}_x^{-}+\sin \left(\frac{\pi l}{L}\right){\mathrm{\partial }}_y^{-}]\delta [\mathit{k}]\\ & \stackrel{\mathcal{F}}{⟷}-\cos \left(\frac{\pi l}{L}\right)(z_2^{-1}-1)-\sin \left(\frac{\pi l}{L}\right)(z_1^{-1}-1).\end{array}$$

7.1. The ‘u-problem’

In order to describe the variational analysis for the u-problem from the point of view of filtering in the Fourier domain for a pyramid decomposition scheme, we define the discrete inverse Fourier transform of $\mathcal{X}(\mathit{z})$ and $\mathcal{Y}(\mathit{z})$ in the ‘u-problem’ in [35] as

Note that the function $\mathcal{X}(\mathit{z})$ satisfies and is similar to the autocorrelation function in [63], which follows the condition of the Riesz basis. We rewrite cartoon u in the Fourier and spatial domain in the scheme of the filter design as

Given β₁ = c₁β₄, the spectra and impulse responses of the scaling and frame functions in (7.1) (l = 0, …, L − 1) are

Equation (7.1) is somewhat similar to the projection operator in the pyramid decomposition scheme (see curvelet [47,53], contourlet [61], steerable wavelet [62], etc.) with:

• — scaling function φ( · ) and its dual $\tilde{\phi }(\cdot )$ with the interpolant $\varphi =\phi \ast \tilde{\phi }$,

• — frame function ψ_l( · ) and its dual ${\tilde{\psi }}_l(\cdot )$,

• — proximity operator for frame elements, i.e. Shrink( · , · ).

Note that the frame soft-thresholding FST( · , · , · ) in (7.1) consists of frame element ψ_l, its dual ${\tilde{\psi }}_l$ and shrinkage operator Shrink( · , · ) together with the updated Lagrange multiplier. Although FST( · , · , · ) is similar to the wavelet/curvelet shrinkage operator, this operator is obtained from the calculus of variation.

It is easy to obtain the unity property of frame elements $(\varphi ,{\psi }_l,{\tilde{\psi }}_l)$ in the Fourier domain as

which satisfies the condition of the perfect reconstruction in the scheme of the wavelet pyramid; see figure 11 for illustration with L = 4 directions. When the number of directions L increases, the bandwidth of the scaling function Φ(z) in the spectral domain reduces (figure 12). Owing to the effect of this lowpass filter, the cartoon u is smoother. These frame elements $(\varphi ,{\psi }_l,{\tilde{\psi }}_l)$, which consist of partial differential operators, can be considered as a kind of a wavelet-like operator [45]. However, the procedure to obtain these elements is from the calculus of variation.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

7.2. The ‘[g_s]^S−1_s=0-problem’

In analogy to the u-problem in the variational framework, we define the frame operators to extract the directional features of the texture v in the pyramid decomposition scheme. The discrete inverse Fourier transforms of $\mathcal{A}(\mathit{z})$ and $\mathcal{B}(\mathit{z})$ in the solution of the ‘v-problem’ [35] are defined as

Similar to $\mathcal{X}(\mathit{z})$, the ‘autocorrelation’ function ${\mathcal{A}}_a(\mathit{z})$ is bounded as in the condition of the Riesz basis [63], i.e. Texture v is rewritten in the Fourier and spatial domains with a = 0, …, S − 1, as Given β₂ = c₂β₃, the spectra and impulse responses of the frame functions in (7.2) are

Similar to the u-problem, it is easy to obtain the unity property of frame elements (ξ, ψ_a, ψ′_a) with a∈[0, S − 1] in the Fourier domain:

which satisfies the condition of the perfect reconstruction; see figure 11 with L = 4 directions and figures 12 and 13 with L = 8 directions.

7.3. The ‘v-problem’

The solution to the v-problem is rewritten as

with Note that there is a shrinkage operator with two terms in texture v, namely a smoothing term and an updated term, balanced by a parameter θ.

7.4. The ‘e-problem’

The noise term ${\parallel \mathcal{C}\{{\chi }_D^c{\cdot }^{\times }\mathit{ϵ}\}\parallel }_{{\ell }_{\mathrm{\infty }}}$ in equation (2.1) is suitable to capture high oscillating patterns, especially small-scale types of noise, because of the advantages of the multi-orientation and multi-scale in a curvelet domain and the supremum norm. For a convenient calculation, we introduce a variable e as a residual ϵ in a known domain Ω\D. The explanation for using the supremum norm of ϵ in the curvelet domain can be found in (3.4). In particular, there are two terms in (3.4), including a residual in Ω\D with an updated Lagrange multiplier λ₅, i.e. (χ^c_D · × ϵ − λ₅/β₅), and its curvelet smoothing term at a level ν, i.e. CST(χ^c_D · × ϵ − λ₅/β₅, ν). Assume that, at iteration t, there are some remaining signals, e.g. texture, in residual ϵ in Ω\D. The curvelet soft-thresholding operator CST( · , · ) reduces noise (or small-scale objects) in (χ^c_D · × ϵ − λ₅/β₅) at a level ν in different scales and orientations. By a subtraction operator in (3.4), e at iteration t contains mainly noise or small-scale objects (figure 14).

• Download figure
• Open in new tab
• Download PowerPoint

7.5. The ‘ϵ-problem’

At iteration t, we denote the updated residual by the unity condition as

and rewrite (3.6) with the indicator functions on unknown domain D and known domain Ω\D as We see that there are two terms in the residual ϵ^(t), including the updated term ϵ^(t)_known for Ω\D and ϵ^(t)_unknown for D. Also, there are two other terms for updating ϵ^(t)_known in Ω\D, namely ϵ^(t)_unity (from a unity condition) and e^(t) (from a curvelet smoothing operator). In contrast to an open loop in a pyramidal decomposition (Laplacian pyramid, wavelet, curvelet, etc.), the solution of DG3PD inpainting can be considered as a closed loop of the pyramid scheme with the true solution (u*, v*, ϵ*, g*) in the minimization problem. The loop of the updated Lagrange multipliers can be seen as a refinement. If we remove this loop from the scheme in figure 15a, the minimization in DG3PD becomes the quadratic penalty method. The removal would cause an imperfect reconstruction for the constraints, e.g. u + v + ϵ≠f; see fig. 7(p) and (q) in [35]. We note that the lowpass, bandpass and highpass filters for the cartoon, texture and residual components are ‘custom-made’ for each image, resulting in different filters for different images (figures 16 and 17). The directional properties of texture and their directional filter banks are illustrated in figures 18 and 19.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

8. Conclusion

In this paper, we have addressed the very challenging task of restoring images with a large number of missing pixels, whose existing pixels are corrupted by noise and, importantly, the image to be restored contains both cartoon and texture elements. The proposed DG3PD model for inpainting and denoising can cope with this threefold difficult problem. The task of simultaneous inpainting and denoising for cartoon and texture components is solved by DG3PD decomposition, followed by inpainting and denoising both components separately and, finally, image restoration by synthesis of the restored components. More specifically, the DG3PD inpainting model is based on a regularization in the Banach space in a discrete setting which is solved by the augmented Lagrangian method and the alternating direction method of multipliers.

In summary, the decomposition step of the proposed method has two major advantages for tackling this challenging problem. Firstly, the decomposition step denoises simultaneously the cartoon and texture component. Secondly, it allows the cartoon and texture to be handled in different ways. The cartoon image is inpainted by DG3PD, as described in §§2 and 3. Separately, the texture component is inpainted followed by further denoising as described in §4. Therefore, the proposed DG3PD decomposition can also be understood as breaking the full problem down into ‘smaller’ subproblems (e.g. texture only inpainting) which are easier to solve.

Image restoration (or image denoising and inpainting) can be understood as an inverse problem and it can be described by a discrete innovation model (figure 10). The assumption in this procedure is that signal has a sparse representation in suitable transform domains. It is known from probability theory that the sparsity (by ℓ₁) is connected to the Laplace distribution of the signal, which results in the heavy tail distribution. Note that the Laplace distribution with the ℓ₁ norm is one of the best approximations of the sparsity by the ℓ₀ norm.

Moreover, by choosing the priors and posterior according to the Bayesian framework and an MAP, we can understand the selection of parameters (see figure 9 for the results of the selection of the heavy-tailed distribution, e.g. the Laplace distribution). Note that the MAP requires an assumption on the noise, e.g. Gaussian or Laplacian, to establish a minimization. In our proposed model, there is no requirement for an assumption on noise distribution, e.g. independent identical distributed, or correlated and Gaussian or non-Gaussian (due to the measurement ${\parallel \mathcal{C}\{\cdot \}\parallel }_{{\ell }_{\mathrm{\infty }}}$, which is similar to the Dantzig selector [60]).

Data accessibility

Data are available from https://figshare.com/s/76d30e70e6e131da2a08.

Author's contributions

D.H.T. and C.G. wrote the manuscript and conducted the experiments. Both authors gave their final approval for publication.

Competing interests

We declare we have no competing interests.

Funding

D.H.T. is supported by the National Science Foundation under grant DMS-1127914 to the Statistical and Applied Mathematical Sciences Institute.

Footnotes