- Email: [email protected]
An approximate sparsity model for inpainting
Appl. Comput. Harmon. Anal. 37 (2014) 171–184
Contents lists available at ScienceDirect
Applied and Computational Harmonic Analysis www.elsevier.com/locate/acha
Letter to the Editor
An approximate sparsity model for inpainting ✩ Lixin Shen a,b,∗ , Yuesheng Xu b,a , Na Zhang b,c a
Department of Mathematics, Syracuse University, Syracuse, NY 13244, USA Guangdong Province Key Laboratory of Computational Science, School of Mathematics and Computational Sciences, Sun Yat-sen University, Guangzhou 510275, PR China c Center for Computer Vision, School of Mathematics and Computational Sciences, Sun Yat-sen University, Guangzhou 510275, PR China b
a r t i c l e
i n f o
Article history: Received 14 May 2013 Received in revised form 18 October 2013 Accepted 7 November 2013 Available online 12 November 2013 Communicated by Charles K. Chui Keywords: Tight framelet Sparsity Inpainting
a b s t r a c t Existing sparse inpainting models often suffer from their over-constraints on the sparsity of the transformed recovered images. Due to the fact that a transformed image of a wavelet or framelet transform is not truly sparse, but approximately sparse, we introduce an approximate sparsity model for inpainting. We formulate the model as minimizing the number of nonzero components of the soft-thresholding operator applied to the transformed image. The key difference of the proposed model from the existing ones is the use of a soft-thresholding operator which shrinkages the components of the transformed image. To efficiently solve the resulting nonconvex optimization problem, we rewrite the 0 norm, which counts the number of nonzero components, as a weighted 1 norm with a nonlinear discontinuous weight function, which is then approximated by a continuous weight function. We overcome the nonlinearity in the weight function by an iteration which leads to a numerical scheme for solving the nonconvex optimization problem. In each iteration, we solve a weighted 1 convex optimization problem. We then focus on understanding the existence of solutions of the weighted 1 convex optimization problem and characterizing them as fixed-points of a nonlinear mapping. The fixed-point formulation allows us to employ efficient iterative algorithms to find the fixed-points. Numerical experiments are shown to demonstrate improvement in performance of the proposed model over the existing models for image inpainting. © 2013 Elsevier Inc. All rights reserved.
1. Introduction In many applications, we often have an image with some pixels missing due to various reasons including impulsive noise caused by malfunctioning pixels in camera sensors, a scratch on text or signature in a picture. It is highly desirable to recover its original image from a damaged image. Image inpainting, first introduced to digital image processing in [1], is a mean to fulfill this goal. Specifically, image inpainting refers to the recovery of missing pixels in a digital image based on pixels available in the observed image. Many of its applications can be found in image processing, such as old films restoration [1], text or scratch removal [2,7], and cloud removal from remotely sensed images [28]. A variety of successful inpainting methods were proposed in the last decade. A group of classical image inpainting algorithms were developed based on numerical solutions of partial differential equations [1,2,13,14]. Algorithms of this kind propagate available information from an observed region to missing regions and experiments show that they perform well
✩ Supported in part by the US National Science Foundation under grant DMS-1115523, by Guangdong Provincial Government of China through the “Computational Science Innovative Research Team” program, by the Natural Science Foundation of China under grants 11071286 and 91130009. Corresponding author at: Department of Mathematics, Syracuse University, Syracuse, NY 13244, USA. Fax: +1 315 443 1475. E-mail address: [email protected] (L. Shen).
*
1063-5203/$ – see front matter © 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.acha.2013.11.002
172
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
for piecewise smooth images with sharp edges. Another popular group of methods are path-based [4,22,34]. Algorithms in this category fill in missing pixels by exploring repetitions of local information and they perform well for missing regions with texture and of large sizes. Recently, approaches based on sparse representation were applied to image inpainting [5,6,8,26]. More precisely, due to the assumption that images are piecewise smooth, these approaches assume that the images can be sparsely represented by some redundant systems, which are generated by a set of transforms such as wavelets, the discrete cosine transform, framelets and curvelets. In this setting, the inpainting problem is modeled as an optimization problem, whose objective function involves the 0 norm of the transformed coefficients of the underlying image with a redundancy system. However, the 0 norm is nonconvex and the associated optimization problem becomes a combinatorial problem which is NP-complete. A common alternative is to replace the 0 norm that appears in the cost functional by the 1 norm. As a result, it relaxes the nonconvex optimization problem to a convex one and therefore it allows us to develop fast algorithms for solving the convex optimization problem. The use of the 1 norm in sparse signal recovery has been well studied in the literature (see e.g. [9,10,21]). Existing sparsity models often suffer from their over-constraints on the sparsity of the recovered images in the transformed domain. Outputs of a wavelet or framelet transform are often not truly sparse and instead, they are approximately sparse in the sense that most of the components of the vectors obtained from such a transform are not zero but have small absolute values compared to the remaining components. It is well known that a wavelet transform, a DCT transform or a framelet transform makes the energy of a transformed image concentrated on certain small regions. In other words, the transformed image is approximately sparse even though it is not truly sparse. This demands a model suitable for approximately sparse signals. The purpose of this paper is to introduce an approximate sparsity inpainting model which takes this important feature of such a transform into account. The model is expressed as the minimization of an objective function which is the 0 norm of the soft-thresholding operator applied to the transformed image over the constrained set described by the region having missing pixels. The key difference of the proposed model from the existing ones is the use of the soft-thresholding operator to truncate the components of the transformed image which are smaller than a threshold. The resulting optimization problem contains the 0 norm in its objective function and thus it is nonconvex. To solve the nonconvex optimization problem efficiently, we rewrite the 0 norm as a weighted 1 norm with a nonlinear discontinuous weight function. We then approximate the discontinuous weight function by a continuous one and use it to design an iteration scheme for solving the nonconvex optimization problem. At each step of the iteration scheme, we solve a weighted 1 convex optimization problem, which is the core minimization problem. Solutions of the core minimization problem are characterized as a system of fixed-point equations, for which efficient iterative algorithms were developed in the last few years [23–25,30,31]. Numerical experiments show that the proposed model outperforms the FBIA method in [8]. This paper is organized in four sections. In Section 2, we propose an approximate sparsity inpainting model and introduce a heuristic method to solve the resulting nonconvex optimization problem that describes the model. The method leads to an iteration scheme for solving the nonconvex problem, each step of whose iteration is to solve a core minimization problem. Section 3 is devoted to a study of the solvability of the core minimization problem and the fixed-point formulation of its solutions. In Section 4, we consider in details the algorithm for solving the approximate sparsity inpainting model. In Appendix A, numerical experiments are presented to demonstrate the accuracy of the proposed model in comparison with well-known existing sparse models. 2. Inpainting model based on approximate sparsity In this section, we propose an approximate sparsity inpainting model and develop an efficient numerical algorithm for solving it. We first present a generic form of an observation model for inpainting. Let f ∈ Rn be the original image defined on the domain Ω := {1, 2, . . . , n}. Let Λ be a proper subset of Ω and Λc be the complementary set of Λ in Ω . The observation model for an image g to be inpainted on Λc may be described as
gi =
fi,
if i ∈ Λ;
arbitrary, otherwise.
(1)
This observation model says that the image g to be inpainted is identical to the original image f on set Λ. Therefore, an inpainted image for the given observation g and set Λ should be sought in the set
C := f : f ∈ Rn and P Λ f = P Λ g ,
(2)
where P Λ is an n × n diagonal matrix whose i-th diagonal entry is 1 if i ∈ Λ and 0 otherwise. The set C describes the constraints that the inpainted image must satisfy. This inpainting problem of inferring unavailable information on set Λc from available ones on set Λ is ill-posed. Therefore, it should be regularized based upon prior knowledge on the original image f . An image is often assumed to be piecewise smooth. It can, therefore, have a sparse representation under a frame transform which has vanishing moments of a high order [29]. Furthermore, tight framelets are often employed since tight framelet transforms can be applied and inverted efficiently [15,17]. A discrete tight framelet is a matrix B such that B B = I ,
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
173
where B is the transpose of B and I is the identity matrix. In particular, we consider in this paper that the matrix B has the form
B=
B0 , W
(3)
where B 0 is an n × n matrix formed by a low-pass filter while W is an m × n matrix formed by high-pass filters associated with the underlying tight framelet system. In general, an image, especially photograph-like one, is not exactly sparse under any tight framelet systems. Instead, it is approximately sparse (or compressibly sparse) [18] in the sense that an accurate approximation of the image can be constructed by using a small number of vectors from the columns of B . These vectors such that the magnitudes of their include all those formed from the columns of B 0 and those from the columns of W inner product with the underlying image are greater than pre-determined thresholds. Since B 0 is formed by a low-pass filter, we do not expect the inner product of its columns with the underlying image to have a small magnitude. One can only expect the inner product of some columns of W with the underlying image to have small absolute values. From this view point, the following inpainting model
min D W f 1
(4)
f ∈C
was studied in [8,26]. Here D is a diagonal matrix whose diagonal entries are positive. The diagonal components of D impose multiparameters for regularization [27]. Many algorithms including the split Bregman iteration [8,26], the first-order primal–dual algorithm [12], and proximity based algorithms [24] are proposed to solve model (4). Model (4) has two drawbacks. First, it looks for an image whose coefficients under the redundant system W are sparse. As a result of this over-constrained on the sparsity, subtle features of the underlying image are missed in the inpainted image. In addition, model (4) uses the 1 norm to impose the sparsity requirement. Intrinsically, for the sparsity requirement, we want our solution to have as few nonzero wavelet (or framelet) coefficients as possible. A natural norm for this should be the 0 norm. We propose below a new inpainting model to address the drawbacks of model (4) by employing the property of the approximate sparsity of the underlying image and using the 0 norm directly. In other words, we need a model which preserves crucial subtle features of the underlying image and at the same time its truncated transformed signal is as sparse as possible. To this end, for a given vector γ ∈ Rd , recall that the hard- and soft-thresholding operators Hγ : Rd → Rd and Sγ : Rd → Rd at u are defined, respectively, as follows:
u i , if |u i | > γi ; Hγ (u ) i := 0, otherwise
Sγ (u ) i := sign(u i ) max |u i | − γi , 0 ,
and
where i ∈ Nd := {1, 2, . . . , d}. By using the hard- and soft-thresholding operators, we expect the vectors Hγ ( W f ) and Sγ ( W f ) are sparse and their sparsity can be naturally counted by the corresponding nonzero components. That is, it can be measured by the 0 norm. Although both hard- and soft-thresholding operators can serve our purpose, we, however, choose to use only the soft-thresholding operator in this paper due to the numerical convenience that will be explained at the end of this section. With this in mind, we introduce the approximate sparsity inpainting model (ASIM):
min Sγ ( W f ) 0 , f ∈C
(5)
where γ ∈ Rm has positive components and · 0 denotes the 0 norm. In the ASIM, the sparsity of the vector Sγ ( W f ) as prior knowledge on the underlying image f is explored in finding a desirable inpainted image from the set C . On the other hand, the objective function of the ASIM is nonconvex and nonsmooth. Thus, it is challenging to develop computationally efficient and mathematically sound algorithms for solving the optimization problem (5). Specifically, the proposed model (5) contains the 0 norm which makes the optimization problem nonconvex. One may replace the 0 norm by the 1 norm to relax the nonconvex optimization problem to a convex one. This is not the best choice since the 1 norm is not as effective to catch the sparsity as the 0 norm. Our idea is to follow the 0 norm as faithfully as possible. To this end, we shall rewrite the 0 norm as a weighted 1 norm with a nonlinear discontinuous weight function. We then approximate the discontinuous weight function by a continuous one and use an iteration to treat the nonlinearity of the weight function. This naturally leads to an iteration scheme for solving the nonconvex optimization problem. Although the idea of approximating the 0 norm with the weighted 1 norm has been reported in [11], we are dealing with the composition of the 0 norm with a soft-thresholding operator, and, only the 0 norm itself was considered in [11]. Furthermore, as we will see in Section 4, how we weight the 1 norm is different from the one suggested in [11]. We first deal with the 0 norm in model (5). We define F : R → R at x ∈ R as
F (x) :=
1
|x| ,
if x = 0;
0,
otherwise.
174
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
Then, for u ∈ Rd we have
u 0 =
d
F (u i )|u i |,
(6)
i =1
that is, the 0 norm of the vector u can be expressed as its “weighted” 1 norm with the weight F (u i ) of its i-th component. The function F has a discontinuity at x = 0 and it presents difficulties for computation. We approximate the discontinuous weight function F by a sequence of positive even continuous functions F on R parameterized by , which have the property that
lim F (x) = F (x)
→0+
for all x ∈ R \ {0}.
A useful example of the function sequence will be given in Section 4. Immediately, for u ∈ Rd we have that
u 0 = lim
d
→0+
F (u i )|u i |.
i =1
Namely, u 0 is approximated by a sequence of weighted 1 norms with the weights depending on the parameter Accordingly, model (5) is approximated by the model depending on the parameter as follows:
min f ∈C
m
F
i =1
Sγ ( W f ) i Sγ ( W f ) i .
and u.
(7)
That is, model (5) is approximated by the weighted 1 model with a nonlinear weight. The nonlinearity of the weight in model (7) may be treated numerically by iteration. This leads to a heuristic method for solving model (7): 1. Initialization: choose an initial guess f (0) . 2. Updating f (k+1) :
f (k+1) ∈ arg min f ∈C
m i =1
F
Sγ W f (k) i Sγ ( W f ) i .
(8)
3. Set k ← k + 1; Repeat step 2 until some stopping criteria are satisfied. We see that the main step of the above procedure is to solve the convex optimization problem (8). We call it the core optimization problem for inpainting. This core optimization problem fits into the general setting studied recently in [12,24]. We remark that if the soft-thresholding operator in (8) is replaced by the hard-thresholding operator, then the optimization problem (8) becomes
arg min f ∈C
m i =1
F
Hγ W f (k) i Hγ ( W f ) i ,
whose objective function is nonconvex. It turns out that the corresponding optimization problem is still difficult to solve. This is the reason why the soft-thresholding is adopted in (5). 3. Solutions of the core optimization problem In this section, we investigate the existence of solutions of the core optimization problem and describe a fixed-point based algorithm to solve it. The existence of the solution for the core optimization problem (8) is essentially determined by the wavelet matrix W appeared in (3). We begin with a description of the construction of matrix W from a set of finite-length filters {τ : = 0, 1, . . . , L } that associate with a tight framelet system. The Fourier series of each filter τ is a trigonometric polynomial and τ . We assume that τ0 is low-pass, that is, τ0 (0) = 1, and the remaining filters are high-pass. For each ω ∈ R, is denoted by these filters satisfy the conditions
τ0 (ω)τ0 (ω + ρπ ) +
L
τ (ω)τ (ω + ρπ ) = δρ ,0 ,
=1
where
ρ = 0, 1. In what follows, we view all filters τ as row vectors. We define recursively
τ,1 = τ ,
τ,k := τ,k−1 ⊗ [1 0], k 2.
(9)
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
175
These vectors τ,k will be used to form the multi-level framelet decomposition matrix B in (3). For simplicity, we consider only the one-dimensional case, which can be extended to the two-dimensional case straightforwardly. Suppose that the dimension of the vector τ,k is smaller than n. For each k, we generate an n × n circulant matrix T ,k whose first row is the vector τ,k with augmentation by zeros from its right end. For a K -level framelet decomposition, the matrices B 0 and W mentioned in (3) are defined by
⎡
B0 =
K −1
T 0, K −k
k =0
T 1, K
K −1 k =1
T 0, K −k
⎤
.. ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ K −. 1 ⎥ ⎢ T L,K T k=1 0, K −k ⎥ ⎢ ⎥ ⎢ . and W = ⎢ .. ⎥. ⎥ ⎢ ⎥ ⎢ T 1,1 ⎥ ⎢ ⎥ ⎢ .. ⎦ ⎣ .
(10)
T L ,1 The next technical lemma states a crucial property of the matrix W for establishing the existence of a solution of the core optimization problem. Lemma 3.1. Let τ0 be the low-pass filter and τ , 1 L, be the high-pass filters associated with a tight framelet system. If the | τ0 (ω)| = 1 occurs only at ω = 0 on the interval [0, π ], then the matrix W after deleting any one of its columns has rank n − 1. Proof. Set H := [ T 1,1 · · · T L,1 ] . It suffices to show that the matrix H after deleting any one of its columns has rank n − 1. Property (9) implies that
T 0,1 T 0,1 + H H = I
(11)
and the eigenvalues of T 0,1 are all less than or equal to 1. Moreover, 1 is the simple eigenvalue of both T 0,1 and T 0,1 T 0,1 due to the assumption that | τ0 (ω)| = 1 occurs only at ω = 0 on the interval [0, π ]. We denote by 1 the all ones n-dimensional vector. Then the vector 1 is the eigenvector of both T 0,1 and T 0,1 T 0,1 corresponding to eigenvalue 1 because τ0 is the low-pass filter. By (11), the rank of the matrix H is at most n − 1. Next, we show that the rank of the matrix H after deleting any one of its columns is exactly n − 1. Suppose this statement is false and the last column of H is removed. Then there exists a nonzero vector v ∈ Rn−1 such that H [ v 0] = 0. This together with identity (11) ensures that the vector [ v 0] is the eigenvector of T 0,1 T 0,1 corresponding to eigenvalue 1.
This contradicts the assumption that 1 is the simple eigenvalue of T 0,1 T 0,1 .
2
We remark that the conditions in Lemma 3.1 imposed on the low-pass filter τ0 are satisfied by many well-known tight framelets. Examples of this type include the compactly supported orthogonal Daubechies wavelets of vanishing moment order q (see [16]),
τ0 (ω) = cos2q (ω/2)
q −1 2q − 1 j =0
j
sin2 j (ω/2) cos2(q−1− j ) (ω/2)
and the spline framelets (see [32]),
τ0 (ω) = e i jω/2 cosq (ω/2),
where j = 0 when q is even , j = 1 when q is odd.
It is straightforward to check that 0 | τ0 (ω)| 1 for all ω ∈ [0, π ] and | τ0 (ω)| equals to 1 only at examples can be found in [15,17]. Essentially, the core optimization problem (8) has the form
min Ξ Sγ ( W f ) 1 , f ∈C
ω = 0. Many other
(12)
where Ξ is a diagonal matrix with positive diagonal entries. Let ξi be the i-th diagonal entry of Ξ . In particular, when ξi is set to F ((Sγ ( W f (k) ))i ), then problem (12) reduces to problem (8). Note that the objective function of problem (12) is · 1 ◦ Ξ ◦ Sγ ◦ W . To establish the existence of the solutions to problem (12), we need the following technical lemma. Lemma 3.2. For an m × n matrix Q , an m × m diagonal matrix Ξ with positive diagonal entries, and a vector γ in Rm with all positive components, if · 1 ◦ Q is coercive, then so is · 1 ◦ Ξ ◦ Sγ ◦ Q .
176
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
Proof. For all f ∈ Rn the following inequality holds:
Ξ Sγ ( Q f ) Ξ Q f 1 − Ξ γ 1 ξmin Q f 1 − Ξ γ 1 , 1 where ξmin is the smallest diagonal entry of Ξ . Since · 1 ◦ Q is coercive, that is, Q f 1 → +∞ as f 2 → +∞, hence, by the above inequality Ξ Sγ ( Q f )1 goes to infinity as well. 2 The constraint set C for the inpainted image that appears in the constrained optimization problem (12) may present computational difficulties in finding its solution. We convert the constrained optimization problem (12) to an unconstrained problem with the help of the indicator function. Specifically, we let ιC denote the indicator function of the closed convex set C , which is defined at f as
ιC :=
if f ∈ C ;
0,
+∞, else.
Then, the constrained problem (12) is written as an equivalent unconstrained form
min ιC ( f ) + Ξ Sγ ( W f ) 1 : f ∈ Rn .
(13)
Next, we establish the existence of a solution of the optimization problem (13). Proposition 3.3. Let τ0 be the low-pass filter and τ , 1 L, be the high-pass filters associated with a tight framelet system, and let C be defined by (2). If | τ0 (ω)| = 1 occurs only at ω = 0 on the interval [0, π ] and the matrix W is constructed by (10), then for any m × m diagonal matrix Ξ with positive diagonal elements, the optimization problem (13) has a minimizer. Proof. It suffices to show that the objective function of the optimization problem (13) is convex and coercive. First, notice that the function max{| · | − a, 0} for any fixed non-negative a is convex on R. Then m Ξ Sγ ( W f ) = ξ j max ( W f ) j − γ j , 0 1 j =1
is a convex function of f on Rn . Next, we show that the objective function is coercive. By Lemma 3.2 we only need to prove that · 1 ◦ W + ιC is coercive on Rn . Clearly, if f is not in C , then W f 1 + ιC ( f ) = +∞. Therefore, we consider only the case f ∈ C . In this case,
W f 1 + ιC ( f ) = W f 1 . By W Λ and W Λc we denote the submatrices extracted from W with column indices in Λ and Λc , respectively. By f Λ and f Λc we denote the subvectors extracted from f with component indices in Λ and Λc , respectively. Hence, f Λ is a constant vector and W Λ f Λ 1 is a constant. Notice that
f 22 = f Λ 22 + f Λc 22 and
W f 1 W Λc f Λc 1 − W Λ f Λ 1 . Further, f 2 → +∞ implies that f Λc 2 → +∞. By Lemma 3.1, matrix W Λc is of full rank. Hence, W Λc f Λc 1 approaches infinity as f Λc 2 goes to infinity. Therefore, W f 1 approaches infinity as f 2 goes to infinity. We conclude that · 1 ◦ W + ιC is coercive on Rn . This completes the proof of this proposition. 2 Proposition 3.3 with the property of the function F ensures that for each k the core optimization problem (8) has a solution. In the remaining part of this section, we characterize a solution of problem (8) in terms of a system of fixed-point equations based on which we describe a fixed-point based algorithm for solving problem (8). We consider a somewhat general minimization problem. By Γ0 (Rd ) we denote the class of all lower semicontinuous convex functions h : Rd → (−∞, +∞] such that dom(h) := {u ∈ Rd : h(u ) < +∞} = ∅. Suppose both ϕ and ψ are in Γ0 and we consider the convex minimization problem:
min
ϕ ( f ) + ψ( W f ): f ∈ Rn .
(14)
Apparently, with ϕ := ιC and ψ(·) := Ξ Sγ (·)1 the optimization problem (8) (or its equivalent problem (13)) can be viewed as a special case of optimization problem (14). In the following discussion, we assume that problem (14) has a solution. Many algorithms were developed for solving problem (14) from different perspectives in the literature. For example, the alternating direction method of multipliers (ADMM) [20] based approaches make use of variable splitting to
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
177
treat the term ψ( W f ). Other techniques, namely Bregman and split Bregman methods [8,21,35], and Douglas–Rachford splitting [3,19] are closely related to ADMM. Fixed-point based algorithms for (14) were developed recently in [24,30,31]. We next describe a fixed-point based algorithm for solving problem (14). For this purpose, we present a characterization of the solutions of optimization problem (14) in terms of a system of fixed-point equations via the proximity operators of the functions ϕ and ψ . Recall that the proximity operator of h, denoted by proxh , is defined by
proxh (x) := arg min
1 2
u − x22 + h(u ): u ∈ Rd
(15)
for all x ∈ Rd . The conjugate of h ∈ Γ0 (Rd ) is the function h∗ ∈ Γ0 (Rd ) defined at y ∈ Rd by h∗ ( y ) := sup{ x, y − h(x): x ∈ Rd }. In the next theorem, we characterize solutions of optimization problem (14). Theorem 3.4. Let ϕ ∈ Γ0 (Rn ), ψ ∈ Γ0 (Rm ), and W be an m × n matrix. If f ∈ Rn is a solution of problem (14), then for any λ, β > 0 there exists a vector y ∈ Rm such that
f = proxλϕ f − λ W y ,
(16)
y = proxβψ ∗ ( y + β W f ).
(17)
Conversely, if there exist λ, β > 0, f ∈ R , and y ∈ R satisfying Eqs. (16) and (17), then f is a solution of problem (14). n
m
A proof of the above theorem may be found in [24,25]. Our purpose is to design an algorithm for solving problem (14) based upon the characterization of solutions to the problem presented in Theorem 3.4. With x = ( f , y ) ∈ Rn × Rm , define
Ψ (x) := λϕ ( f ) + βψ ∗ ( y ) and
E :=
−λ W
In βW
Im
(18)
(19)
.
Note that
proxΨ (x) = proxλϕ ( f ), proxβψ ∗ ( y ) ∈ Rn × Rm . Then Eqs. (16) and (17) can be refreshed into a single equation as follows:
x = proxΨ ( Ex).
(20)
We further note that E 2 , the spectral norm of E (i.e., the largest singular value of E) is strictly greater than 1 for any nonzero matrix W . As a consequence, the operator proxΨ ◦ E may not be nonexpansive. Hence, we cannot seek a solution to problem (14) via a simple iteration
xk+1 = proxΨ Exk
for any initial estimate x0 ∈ Rn+m because convergence of the resulting sequence {xk : k ∈ N} cannot be guaranteed. To overcome the difficulty caused by the expansiveness of the matrix E, we split E as the sum of M and E − M for a proper (m + n) × (m + n) matrix M. An appropriate splitting of matrix E may lead to a convergent iteration scheme
xk+1 = proxΨ Mxk+1 + ( E − M )xk ,
(21)
for an initial estimate x0 ∈ Rn+m . The corresponding limit is a solution of the fixed-point equation (20). There are many ways to choose the matrix M. In general, iteration (21) is an implicit scheme. Systematic study of the implicit scheme (21) in a general context will be investigated in a different occasion. Here, we give two choices of the matrix M which leads to explicit iterations. Specifically, when we choose
0 M= 2β W
0 , 0
the iterative scheme (21) becomes
f k+1 = proxλϕ f k − λ W yk ,
yk+1 = proxβψ ∗ yk + β W 2 f k+1 − f k .
(22)
178
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
When we choose
M=
0 −2λ W , 0 0
the iterative scheme (21) becomes
yk+1 = proxβψ ∗ yk + β W f k ,
f k+1 = proxλϕ f k − λ W 2 yk+1 − yk .
(23)
These two iterations in fact coincide with the primal–dual algorithms which were originally developed from a different point of view [12] applied to problem (14). It was known (see, e.g., [12,25,36]) that if λβ < 1/ W 22 , then the sequence { f k : k ∈ N} generated from (22) or from (23) converges to a solution of problem (14). 4. An algorithm for solving the ASIM In this section, we discuss crucial issues for implementation of the heuristic method for the ASIM (5) in the context of inpainting. First of all, for the inpainting model to be efficient we need to choose an appropriate transform W . Moreover, the main step of the heuristic method is to find a solution of the optimization problem (8), which can be solved by the proposed fixed-point equation based algorithms via proximity operators. They are the iterative scheme (22) or (23). The implementation of these algorithms requires that the proximity operators of both λϕ and βψ ∗ are specifically given. We shall present specific formula for computing the proximity operators of convex functions λϕ and βψ ∗ . We discuss the choice of the transform W . The piecewise cubic tight framelet system (see [32]) will be adopted in the our model. The same system is also used in the models which will be used for comparison. The filters associated with the system are
⎧ ⎪ τ0 = ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ τ1 = ⎪ ⎪ ⎪ ⎪ ⎨ τ2 = ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ τ3 = ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩τ = 4
1 [1, 4, 6, 4, 1], 16 1 [1, 2, 0, −2, −1], 8
√
6
16 1
[−1, 0, 2, 0, −1],
[−1, 2, 0, −2, 1],
8 1
16
[1, −4, 6, −4, 1].
The transform matrix W in (5) is then built from these filters. The approximate sparsity model involves the function F which approximates the function F . We now choose
1
F (x) =
,
if 0 |x| ;
1 |x| ,
if |x| > .
Clearly, F is continuous and it converges to F pointwise as goes to zero. We remark that the function F (x) = 1/(|x| + ) was used in [11]. For the inpainting model, we are required to compute the proximity operator of the two functions
ϕ := ιC and ψ(·) := · 1 ◦ Ξ ◦ Sγ .
(24)
Here C is the set in which the inpainted image lies, Ξ is a diagonal matrix having positive diagonal entries determined by the approximate function F and γ is a vector with non-negative components. In the heuristic method, the function ϕ is kept unchanged while the function ψ changes iteratively via updating the matrix Ξ . One can easily check that for f ∈ Rn
proxλϕ ( f ) = P Λ g + ( I − P Λ ) f . Note that the right-hand side of the formula above is independent of λ due to the fact of proxλϕ = proxϕ . Recall that P Λ is an n × n diagonal matrix whose i-th diagonal entry is 1 if i ∈ Λ and 0 otherwise. We now derive the proximity operator of the function ψ ∗ . For any y ∈ Rm , one has that
ψ( y ) =
m i =1
ξi · max | y i | − γi , 0 .
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
179
Fig. 1. Plots of (a) ξ, , (b) ∗ξ, , and (c) prox∗ . ξ,
For given non-negative numbers ξ and , we define a function ξ, at x ∈ R by
ξ, (x) := ξ · max |x| − , 0 . We then have that ψ( y ) =
ψ ∗ ( y) =
m
m
i =1 ξi ,γi ( y i )
which is separable. Therefore,
∗ξi ,γi ( y i )
i =1
and
proxβψ ∗ ( y ) = proxβ∗
ξ1 ,γ1
( y 1 ), proxβ∗ξ
2 ,γ2
( y 2 ), . . . , proxβ∗ξ
m ,γm
( ym ) .
Thus, the proximity operator of βψ ∗ can be obtained with the help of the following lemma. Lemma 4.1. For two non-negative numbers ξ and , we have at x ∈ R that
∗ξ, (x) =
|x|,
if |x| ξ ;
(25)
+∞, otherwise.
Further, for any positive number β , we have that
⎧ if |x| > β + ξ ; ⎨ sign(x)ξ, proxβ∗ (x) = x − sign(x)β , if β |x| β + ξ ; ξ, ⎩ 0, otherwise.
(26)
Proof. By the definition of the conjugate function, we have that ∗ξ, (x) = sup{xy − ξ, ( y ): y ∈ R} at x ∈ R. From the definition of ξ, we observe that
xy − ξ, ( y ) =
xy ,
if | y | ;
xy − ξ(| y | − ),
otherwise.
Hence, for the fixed x, the supremum of xy − ξ, ( y ) is |x| for | y | ; or +∞ if | y | > and |x| > ξ ; or |x| if | y | > and |x| ξ . Therefore, Eq. (25) holds. Now, we prove the formula (26). For x ∈ R, we write y = proxβ∗ (x). We have that x − y ∈ ∂(β∗ξ, )( y ) by using Propoξ,
sition 2.6 in [30], where ∂ denotes the subdifferential. Since ∂(β∗ξ, ) = β∂∗ξ, , we get (x − y )/β ∈ ∂∗ξ, ( y ). On the other hand, for any y ∈ dom(∗ξ, ) = [−ξ, ξ ] we have that (it is also clear from Fig. 1(b))
⎧ (−∞, − ], ⎪ ⎪ ⎪ ⎪ − ⎪ ⎨ , ∗ ∂ξ, ( y ) = [− , ], ⎪ ⎪ ⎪ , ⎪ ⎪ ⎩ [ , ∞),
y = −ξ ;
−ξ < y < 0; y = 0; 0 < y < ξ; y = ξ,
which, together with the inclusion (x − y )/β ∈ ∂∗ξ, ( y ), yields Eq. (26).
2
180
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
Algorithm 1 (Iterative scheme for the ASIM (5)). Input: the m × n matrix W ; positive numbers λ, β , τ , and α with λβ < 1 and 0 < α < 1; positive integers kmax and j max ; 0 0 Initialization: choose 0 > 0, γ ∈ Rm + , f , y , diagonal matrix Ξ with positive diagonal entries, and ψ = · 1 ◦ Ξ ◦ Sγ ; set k = 0. while k < kmax do Step 1: Compute f k+1 : Initialize x0 = f k , u 0 = yk , and j = 0 while j < j max and x j − x j −1 2 > τ x j 2 do
⎧ j +1 ← P Λ g + ( I − P Λ )(x j − λ W u j ) ⎪ ⎨x j +1 u ← proxβψ ∗ (u j + β W (2x j +1 − x j )) ⎪ ⎩ j← j+1 end while Step 2: Update f k+1 ← x j , yk+1 ← u j ,
k+1 ← α k , and ψ ← · 1 ◦ Ξ ◦ Sγ , where Ξ = diag( F k+1 ( w 1 ), F k+1 ( w 2 ), . . . , F k+1 ( w m )) with w = Sγ ( W f k+1 ) and w i being the i-th component of w for i = 1, 2, . . . , m. Step 3: k ← k + 1. end while
We make two remarks on the function ξ, . First, the function ξ, (see Fig. 1(a)) is the -insensitive loss function scaled by ξ that is widely used in support vector machines and regression problems [33]. In viewing model (14) with ϕ and ψ given by (24), the i-th framelet coefficient ( W f )i has no cost if it lies in the γi -neighborhood of the origin; otherwise, its contribution will be weighted by the parameter ξi in the objective function of model (14). Second, in the proximity operator of β∗ξ, at x ∈ R, the number β acts as a threshold. The operator proxβ∗ (see Fig. 1(c)) performs as ξ,
a thresholding operator with β as its threshold on the interval [−(β + ξ ), (β + ξ )]; when x is outside this interval, the proxβ∗ (x) saturates at ξ if x is positive; and −ξ otherwise. ξ,
To develop a convergent algorithm for solving problem (8), we also need to know the 2 norm of matrix W . The size of W is m × n with m = 24nK , where n is the number of pixels of the underlying image and K is the decomposition level. The 2 norm of W is always equal to 1 for any K . Thus, in iterative scheme (22) or (23), we choose the parameters λ and β such that λβ < 1. In the experiments to be presented in Appendix A, we specifically choose K = 1. Now, we are ready to present a complete and implementable algorithm, as described in Algorithm 1, for solving the ASIM (5). There are two “while” loops in the algorithm. The outer-loop is to update ψ and rescale by a factor α < 1 in m each iteration. Changing the value of parameter is mainly to enforce i =1 F (u i )|u i | to be close enough to u 0 . The maximum number of outer-loops is kmax . The inner-loop is to solve the core optimization problem (8) by using the iterative scheme (22). The maximum number of inner-loops is set to be j max . The algorithm for solving the ASIM by using the iterative scheme (23) can be given accordingly. Appendix A. Examples In this appendix, we present numerical results to demonstrate the performance of our proposed Algorithm 1 for image inpainting and compare it with Algorithm FBIA [8] since other existing algorithms either have a numerical behavior similar to that of FBIA or are not comparable to FBIA. The objective quality of restored images is evaluated quantitatively by the peak signal-to-noise ratio (PSNR). All algorithms are implemented in Matlab and executed on MATLAB R2010a running on a PC equipped with an Intel Core 2 Quad CPU at 2.2.00 GHz and 2G RAM memory. The grayscale images of “Cameraman” and “House” of size 256 × 256 and the color images of “MonaLisa” of size 330 × 450 and “Madonna&Child” of size of 880 × 1200 shown in Fig. 2 are used as our test images. The grayscale images are superimposed with “Text 1” (see the first column of Fig. 3) and “Text 2” (see the first column of Fig. 5) while the color images with scratches are displayed in the first column of Fig. 7. The inpainting problems that we consider here are to restore these images overlying with texts or having scratches by the proposed Algorithm 1 and the FBIA [8]. For Algorithm 1, we fix λ = 0.5, α = 0.8, γi = 0.25 for i = 1, 2, . . . , 24n, and τ = 0.0005 for all the experiments. To ensure the convergence of Algorithm 1, we choose β = 0.999/λ. The initial points for all the experiments are chosen as f 0 = g, y 0 = 0, 0 = 1, and Ξ = diag(0n×n , 20I 24n×24n ). Moreover, for all experiments on both images of “Cameraman” and “House”, we choose ( jmax , kmax ) = (30, 4) for experiments on “Text 1” and ( jmax , kmax ) = (60, 7) for those on “Text 2”. For a color image, we apply Algorithm 1 for its red, green, and blue channels, separately. For the image of “MonaLisa”, we set jmax = 30 for all its red, green, and blue channels, while kmax = 3, 5, 7 for its red, green, and blue channels, respectively. For the image of “Madonna&Child”, we set ( jmax , kmax ) to be (20, 6), (70, 3), and (12, 8) for its red, green, and blue channels respectively. For FBIA, the parameters are determined experimentally for the restored image to achieve the best possible PSNR-value. We present the numerical results in one table and six figures. Table 1 lists the PSNR values of the inpainted images from model (5) by using Algorithm 1 and from model (4) by using Algorithm FBIA. It clearly indicates that model (5) with Algorithm 1 outperforms model (4) with Algorithm FBIA. Fig. 3 shows the visual comparison of the inpainted images from models (5) and (4) for the images of “Cameraman” and “House” overlaid by “Text 1”. Although it can be observed that both models inpaint the images well, differences in visual quality of the restored images can be observed from the zoom-in images shown in Fig. 4. One can see that model (5) with Algorithm 1 can efficiently remove the superposed words over the regions having cartoon components. Although the words in “Text 2” are denser and thicker than those in “Text 1”, from
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
181
Fig. 2. The original images.
Fig. 3. The results of recovering images overlying with “Text 1”: (a) images superposed by “Text 1”; (b) restored images by model (4) with FBIA; (c) restored images by model (5) with Algorithm 1. The zoomed-in part of each image is marked by a square in the corresponding image that will be highlighted in Fig. 4.
Fig. 4. Column 1: Zoom-in portions of the original test images; Column 2: Zoom-in portions of the original test images with “Text 1”; Column 3: inpainted images by model (4) with FBIA; Column 4: inpainted images by model (5) with Algorithm 1.
182
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
Table 1 PSNR (dB) of restored images by model (5) with Algorithm 1 and model (4) with FBIA. Methods
Algorithm 1 FBIA
“Text 1”
“Text 2”
“Scratch”
Cameraman
House
Cameraman
House
MonaLisa
Madonna&Child
30.91 30.29
38.29 37.37
27.32 26.83
33.21 31.51
46.60 45.45
46.62 46.24
Fig. 5. The results of recovering images overlying with “Text 2”: (a) images superposed by “Text 2”; (b) restored images by model (4) with FBIA; (c) restored images by model (5) with Algorithm 1. The zoomed-in part of each image is marked by a square in the corresponding image that will be highlighted in Fig. 6.
Fig. 6. Column 1: Zoom-in portions of the original test images; Column 2: Zoom-in portions of the original test images with “Text 2”; Column 3: inpainted images by model (4) with FBIA; Column 4: inpainted images by model (5) with Algorithm 1.
Figs. 5 and 6 we can draw the same conclusion for the images having “Text 2”. The restored color images are shown in Figs. 7 and 8. Again, these figures show that model (5) with Algorithm 1 outperforms model (4) with Algorithm FBIA. In the spirit of reproducible research, the Matlab source codes for reproducing the results in this section can be obtained by sending an email to [email protected].
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
183
Fig. 7. The results of recovering images overlying with scratch: (a) images superposed by scratch; (b) restored images by model (4) with FBIA; (c) restored images by model (5) with Algorithm 1. The zoomed-in part of each image is marked by a square in the corresponding image that will be highlighted in Fig. 8.
Fig. 8. Column 1: Zoom-in portions of the original test images; Column 2: Zoom-in portions of the original test images with “Text 3”; Column 3: inpainted images by model (4) with FBIA; Column 4: inpainted images by model (5) with Algorithm 1.
References [1] M. Bertalmio, G. Sapiro, V. Caselles, C. Ballester, Image inpainting, in: Proceedings of SIGGRAPH, New Orleans, LA, 2000, pp. 417–424. [2] M. Bertalmio, L. Vese, G. Sapiro, S. Osher, Simultaneous structure and texture image inpainting, IEEE Trans. Image Process. 12 (2003) 882–889. [3] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, Distributed optimization and statistical learning via alternating direction method of multipliers, Found. Trends Mach. Learn. 3 (1) (2010) 1–122. [4] A. Bugeau, M. Bertalmio, V. Caselles, A comprehensive framework for image inpainting, IEEE Trans. Image Process. (2010) 2634–2645. [5] J.-F. Cai, R. Chan, L. Shen, Z. Shen, Simultaneously inpainting in image and transformed domains, Numer. Math. 112 (2009) 509–533. [6] J.-F. Cai, R. Chan, Z. Shen, A framelet-based image inpainting algorithm, Appl. Comput. Harmon. Anal. 24 (2007) 131–149. [7] J.-F. Cai, R.H. Chan, L. Shen, Z. Shen, Convergence analysis of tight framelet approach for missing data recovery, Adv. Comput. Math. 31 (2009) 87–113. [8] J.-F. Cai, S. Osher, Z. Shen, Split Bregman methods and frame based image restoration, Multiscale Model. Simul. 2 (2009) 337–369.
184
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
[9] E. Candes, J. Romberg, T. Tao, Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information, IEEE Trans. Inform. Theory 52 (2006) 489–509. [10] E. Candes, J. Romberg, T. Tao, Stable signal recovery from incomplete and inaccurate measurements, Comm. Pure Appl. Math. 59 (2006) 1207–1223. [11] E. Candes, M.B. Wakin, S.P. Boyd, Enhancing sparsity by reweighted 1 minimization, J. Fourier Anal. Appl. 14 (2008) 877–905. [12] A. Chambolle, T. Pock, A first-order primal-dual algorithm for convex problems with applications to imaging, J. Math. Imaging Vision 40 (2011) 120–145. [13] T. Chan, S.-H. Kang, J. Shen, Euler’s elastica and curvature-based image inpainting, SIAM J. Appl. Math. 63 (2002) 564–592. [14] T. Chan, J. Shen, Mathematical models for local nontexture inpaintings, SIAM J. Appl. Math. 62 (2002) 1019–1043. [15] C.K. Chui, W. He, J. Stockler, Compactly supported tight and sibling frames with maximum vanishing moments, Appl. Comput. Harmon. Anal. 13 (2002) 224–262. [16] I. Daubechies, Ten Lectures on Wavelets, CBMS-NSF Regional Conf. Ser. in Appl. Math., vol. 61, SIAM, Philadelphia, 1992. [17] I. Daubechies, B. Han, A. Ron, Z. Shen, Framelets: MRA-based constructions of wavelet frames, Appl. Comput. Harmon. Anal. 14 (2003) 1–46. [18] D. Donoho, M. Vetterli, R. DeVore, I. Daubechies, Data compression and harmonic analysis, IEEE Trans. Inform. Theory 44 (1998) 2435–2476. [19] J. Eckstein, D.P. Bertsekas, On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators, Math. Program. 55 (1992) 293–318. [20] D. Gabay, B. Mercier, A dual algorithm for the solution of nonlinear variational problems via finite-element approximations, Comput. Math. Appl. 2 (1976) 17–40. [21] T. Goldstein, S. Osher, The split Bregman method for 1 regularization problems, SIAM J. Imaging Sci. 2 (2009) 323–343. [22] N. Komodakis, G. Tziritas, Image completion using efficient belief propagation via priority scheduling and dynamic pruning, IEEE Trans. Image Process. 16 (2007) 2649–2661. [23] A. Krol, S. Li, L. Shen, Y. Xu, Preconditioned alternating projection algorithms for maximum a posteriori ECT reconstruction, Inverse Problems 28 (2012) 115005, 34pp. [24] Q. Li, C.A. Micchelli, L. Shen, Y. Xu, A proximity algorithm accelerated by Gauss–Seidel iterations for L1/TV denoising models, Inverse Problems 28 (2012) 095003, 20 pp. [25] Q. Li, L. Shen, Y. Xu, N. Zhang, Multi-step fixed-point proximity algorithms for solving a class of regularization problems arising from image processing, UCLA CAM Report 12-65, 2012. [26] Q. Li, L. Shen, L. Yang, Split-Bregman iteration for framelet based image inpainting, Appl. Comput. Harmon. Anal. 32 (2012) 145–154. [27] Y. Lu, L. Shen, Y. Xu, Multi-parameter regularization methods for high-resolution image reconstruction with displacement errors, IEEE Trans. Circuits Syst. I, Fundam. Theory Appl. 54 (2007) 1788–1799. [28] A. Maalouf, P. Carre, B. Augereau, C. Fernandez-Maloigne, A bandelet-based inpainting technique for clouds removal from remotely sensed images, IEEE Trans. Geosci. Remote Sens. 47 (2009) 2363–2371. [29] S. Mallat, A Wavelet Tour of Signal Processing, 2nd ed., Academic Press, New York, 1999. [30] C.A. Micchelli, L. Shen, Y. Xu, Proximity algorithms for image models: Denoising, Inverse Problems 27 (2011), 045009 (30 pp). [31] C.A. Micchelli, L. Shen, Y. Xu, X. Zeng, Proximity algorithms for the L1/TV image denoising model, Adv. Comput. Math. 38 (2013) 401–426. [32] A. Ron, Z. Shen, Affine system in L 2 ( R d ): the analysis of the analysis operator, J. Funct. Anal. 148 (1997) 408–447. [33] V. Vapnik, The Nature of Statistical Learning Theory (Information Science and Statistics), 2nd ed., Springer-Verlag, New York, 1999. [34] Y. Wexler, E. Shechtman, M. Irani, Space–time completion of video, IEEE Trans. Pattern Anal. Mach. Intell. 29 (2007) 463–476. [35] W. Yin, S. Osher, D. Goldfarb, J. Darbon, Bregman iterative algorithms for 1 minimization with applications to compressed sensing, SIAM J. Imaging Sci. 1 (2008) 143–168. [36] X. Zhang, M. Burger, S. Osher, A unified primal–dual algorithm framework based on Bregman iteration, J. Sci. Comput. 46 (2011) 20–46.
Contents lists available at ScienceDirect
Applied and Computational Harmonic Analysis www.elsevier.com/locate/acha
Letter to the Editor
An approximate sparsity model for inpainting ✩ Lixin Shen a,b,∗ , Yuesheng Xu b,a , Na Zhang b,c a
Department of Mathematics, Syracuse University, Syracuse, NY 13244, USA Guangdong Province Key Laboratory of Computational Science, School of Mathematics and Computational Sciences, Sun Yat-sen University, Guangzhou 510275, PR China c Center for Computer Vision, School of Mathematics and Computational Sciences, Sun Yat-sen University, Guangzhou 510275, PR China b
a r t i c l e
i n f o
Article history: Received 14 May 2013 Received in revised form 18 October 2013 Accepted 7 November 2013 Available online 12 November 2013 Communicated by Charles K. Chui Keywords: Tight framelet Sparsity Inpainting
a b s t r a c t Existing sparse inpainting models often suffer from their over-constraints on the sparsity of the transformed recovered images. Due to the fact that a transformed image of a wavelet or framelet transform is not truly sparse, but approximately sparse, we introduce an approximate sparsity model for inpainting. We formulate the model as minimizing the number of nonzero components of the soft-thresholding operator applied to the transformed image. The key difference of the proposed model from the existing ones is the use of a soft-thresholding operator which shrinkages the components of the transformed image. To efficiently solve the resulting nonconvex optimization problem, we rewrite the 0 norm, which counts the number of nonzero components, as a weighted 1 norm with a nonlinear discontinuous weight function, which is then approximated by a continuous weight function. We overcome the nonlinearity in the weight function by an iteration which leads to a numerical scheme for solving the nonconvex optimization problem. In each iteration, we solve a weighted 1 convex optimization problem. We then focus on understanding the existence of solutions of the weighted 1 convex optimization problem and characterizing them as fixed-points of a nonlinear mapping. The fixed-point formulation allows us to employ efficient iterative algorithms to find the fixed-points. Numerical experiments are shown to demonstrate improvement in performance of the proposed model over the existing models for image inpainting. © 2013 Elsevier Inc. All rights reserved.
1. Introduction In many applications, we often have an image with some pixels missing due to various reasons including impulsive noise caused by malfunctioning pixels in camera sensors, a scratch on text or signature in a picture. It is highly desirable to recover its original image from a damaged image. Image inpainting, first introduced to digital image processing in [1], is a mean to fulfill this goal. Specifically, image inpainting refers to the recovery of missing pixels in a digital image based on pixels available in the observed image. Many of its applications can be found in image processing, such as old films restoration [1], text or scratch removal [2,7], and cloud removal from remotely sensed images [28]. A variety of successful inpainting methods were proposed in the last decade. A group of classical image inpainting algorithms were developed based on numerical solutions of partial differential equations [1,2,13,14]. Algorithms of this kind propagate available information from an observed region to missing regions and experiments show that they perform well
✩ Supported in part by the US National Science Foundation under grant DMS-1115523, by Guangdong Provincial Government of China through the “Computational Science Innovative Research Team” program, by the Natural Science Foundation of China under grants 11071286 and 91130009. Corresponding author at: Department of Mathematics, Syracuse University, Syracuse, NY 13244, USA. Fax: +1 315 443 1475. E-mail address: [email protected] (L. Shen).
*
1063-5203/$ – see front matter © 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.acha.2013.11.002
172
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
for piecewise smooth images with sharp edges. Another popular group of methods are path-based [4,22,34]. Algorithms in this category fill in missing pixels by exploring repetitions of local information and they perform well for missing regions with texture and of large sizes. Recently, approaches based on sparse representation were applied to image inpainting [5,6,8,26]. More precisely, due to the assumption that images are piecewise smooth, these approaches assume that the images can be sparsely represented by some redundant systems, which are generated by a set of transforms such as wavelets, the discrete cosine transform, framelets and curvelets. In this setting, the inpainting problem is modeled as an optimization problem, whose objective function involves the 0 norm of the transformed coefficients of the underlying image with a redundancy system. However, the 0 norm is nonconvex and the associated optimization problem becomes a combinatorial problem which is NP-complete. A common alternative is to replace the 0 norm that appears in the cost functional by the 1 norm. As a result, it relaxes the nonconvex optimization problem to a convex one and therefore it allows us to develop fast algorithms for solving the convex optimization problem. The use of the 1 norm in sparse signal recovery has been well studied in the literature (see e.g. [9,10,21]). Existing sparsity models often suffer from their over-constraints on the sparsity of the recovered images in the transformed domain. Outputs of a wavelet or framelet transform are often not truly sparse and instead, they are approximately sparse in the sense that most of the components of the vectors obtained from such a transform are not zero but have small absolute values compared to the remaining components. It is well known that a wavelet transform, a DCT transform or a framelet transform makes the energy of a transformed image concentrated on certain small regions. In other words, the transformed image is approximately sparse even though it is not truly sparse. This demands a model suitable for approximately sparse signals. The purpose of this paper is to introduce an approximate sparsity inpainting model which takes this important feature of such a transform into account. The model is expressed as the minimization of an objective function which is the 0 norm of the soft-thresholding operator applied to the transformed image over the constrained set described by the region having missing pixels. The key difference of the proposed model from the existing ones is the use of the soft-thresholding operator to truncate the components of the transformed image which are smaller than a threshold. The resulting optimization problem contains the 0 norm in its objective function and thus it is nonconvex. To solve the nonconvex optimization problem efficiently, we rewrite the 0 norm as a weighted 1 norm with a nonlinear discontinuous weight function. We then approximate the discontinuous weight function by a continuous one and use it to design an iteration scheme for solving the nonconvex optimization problem. At each step of the iteration scheme, we solve a weighted 1 convex optimization problem, which is the core minimization problem. Solutions of the core minimization problem are characterized as a system of fixed-point equations, for which efficient iterative algorithms were developed in the last few years [23–25,30,31]. Numerical experiments show that the proposed model outperforms the FBIA method in [8]. This paper is organized in four sections. In Section 2, we propose an approximate sparsity inpainting model and introduce a heuristic method to solve the resulting nonconvex optimization problem that describes the model. The method leads to an iteration scheme for solving the nonconvex problem, each step of whose iteration is to solve a core minimization problem. Section 3 is devoted to a study of the solvability of the core minimization problem and the fixed-point formulation of its solutions. In Section 4, we consider in details the algorithm for solving the approximate sparsity inpainting model. In Appendix A, numerical experiments are presented to demonstrate the accuracy of the proposed model in comparison with well-known existing sparse models. 2. Inpainting model based on approximate sparsity In this section, we propose an approximate sparsity inpainting model and develop an efficient numerical algorithm for solving it. We first present a generic form of an observation model for inpainting. Let f ∈ Rn be the original image defined on the domain Ω := {1, 2, . . . , n}. Let Λ be a proper subset of Ω and Λc be the complementary set of Λ in Ω . The observation model for an image g to be inpainted on Λc may be described as
gi =
fi,
if i ∈ Λ;
arbitrary, otherwise.
(1)
This observation model says that the image g to be inpainted is identical to the original image f on set Λ. Therefore, an inpainted image for the given observation g and set Λ should be sought in the set
C := f : f ∈ Rn and P Λ f = P Λ g ,
(2)
where P Λ is an n × n diagonal matrix whose i-th diagonal entry is 1 if i ∈ Λ and 0 otherwise. The set C describes the constraints that the inpainted image must satisfy. This inpainting problem of inferring unavailable information on set Λc from available ones on set Λ is ill-posed. Therefore, it should be regularized based upon prior knowledge on the original image f . An image is often assumed to be piecewise smooth. It can, therefore, have a sparse representation under a frame transform which has vanishing moments of a high order [29]. Furthermore, tight framelets are often employed since tight framelet transforms can be applied and inverted efficiently [15,17]. A discrete tight framelet is a matrix B such that B B = I ,
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
173
where B is the transpose of B and I is the identity matrix. In particular, we consider in this paper that the matrix B has the form
B=
B0 , W
(3)
where B 0 is an n × n matrix formed by a low-pass filter while W is an m × n matrix formed by high-pass filters associated with the underlying tight framelet system. In general, an image, especially photograph-like one, is not exactly sparse under any tight framelet systems. Instead, it is approximately sparse (or compressibly sparse) [18] in the sense that an accurate approximation of the image can be constructed by using a small number of vectors from the columns of B . These vectors such that the magnitudes of their include all those formed from the columns of B 0 and those from the columns of W inner product with the underlying image are greater than pre-determined thresholds. Since B 0 is formed by a low-pass filter, we do not expect the inner product of its columns with the underlying image to have a small magnitude. One can only expect the inner product of some columns of W with the underlying image to have small absolute values. From this view point, the following inpainting model
min D W f 1
(4)
f ∈C
was studied in [8,26]. Here D is a diagonal matrix whose diagonal entries are positive. The diagonal components of D impose multiparameters for regularization [27]. Many algorithms including the split Bregman iteration [8,26], the first-order primal–dual algorithm [12], and proximity based algorithms [24] are proposed to solve model (4). Model (4) has two drawbacks. First, it looks for an image whose coefficients under the redundant system W are sparse. As a result of this over-constrained on the sparsity, subtle features of the underlying image are missed in the inpainted image. In addition, model (4) uses the 1 norm to impose the sparsity requirement. Intrinsically, for the sparsity requirement, we want our solution to have as few nonzero wavelet (or framelet) coefficients as possible. A natural norm for this should be the 0 norm. We propose below a new inpainting model to address the drawbacks of model (4) by employing the property of the approximate sparsity of the underlying image and using the 0 norm directly. In other words, we need a model which preserves crucial subtle features of the underlying image and at the same time its truncated transformed signal is as sparse as possible. To this end, for a given vector γ ∈ Rd , recall that the hard- and soft-thresholding operators Hγ : Rd → Rd and Sγ : Rd → Rd at u are defined, respectively, as follows:
u i , if |u i | > γi ; Hγ (u ) i := 0, otherwise
Sγ (u ) i := sign(u i ) max |u i | − γi , 0 ,
and
where i ∈ Nd := {1, 2, . . . , d}. By using the hard- and soft-thresholding operators, we expect the vectors Hγ ( W f ) and Sγ ( W f ) are sparse and their sparsity can be naturally counted by the corresponding nonzero components. That is, it can be measured by the 0 norm. Although both hard- and soft-thresholding operators can serve our purpose, we, however, choose to use only the soft-thresholding operator in this paper due to the numerical convenience that will be explained at the end of this section. With this in mind, we introduce the approximate sparsity inpainting model (ASIM):
min Sγ ( W f ) 0 , f ∈C
(5)
where γ ∈ Rm has positive components and · 0 denotes the 0 norm. In the ASIM, the sparsity of the vector Sγ ( W f ) as prior knowledge on the underlying image f is explored in finding a desirable inpainted image from the set C . On the other hand, the objective function of the ASIM is nonconvex and nonsmooth. Thus, it is challenging to develop computationally efficient and mathematically sound algorithms for solving the optimization problem (5). Specifically, the proposed model (5) contains the 0 norm which makes the optimization problem nonconvex. One may replace the 0 norm by the 1 norm to relax the nonconvex optimization problem to a convex one. This is not the best choice since the 1 norm is not as effective to catch the sparsity as the 0 norm. Our idea is to follow the 0 norm as faithfully as possible. To this end, we shall rewrite the 0 norm as a weighted 1 norm with a nonlinear discontinuous weight function. We then approximate the discontinuous weight function by a continuous one and use an iteration to treat the nonlinearity of the weight function. This naturally leads to an iteration scheme for solving the nonconvex optimization problem. Although the idea of approximating the 0 norm with the weighted 1 norm has been reported in [11], we are dealing with the composition of the 0 norm with a soft-thresholding operator, and, only the 0 norm itself was considered in [11]. Furthermore, as we will see in Section 4, how we weight the 1 norm is different from the one suggested in [11]. We first deal with the 0 norm in model (5). We define F : R → R at x ∈ R as
F (x) :=
1
|x| ,
if x = 0;
0,
otherwise.
174
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
Then, for u ∈ Rd we have
u 0 =
d
F (u i )|u i |,
(6)
i =1
that is, the 0 norm of the vector u can be expressed as its “weighted” 1 norm with the weight F (u i ) of its i-th component. The function F has a discontinuity at x = 0 and it presents difficulties for computation. We approximate the discontinuous weight function F by a sequence of positive even continuous functions F on R parameterized by , which have the property that
lim F (x) = F (x)
→0+
for all x ∈ R \ {0}.
A useful example of the function sequence will be given in Section 4. Immediately, for u ∈ Rd we have that
u 0 = lim
d
→0+
F (u i )|u i |.
i =1
Namely, u 0 is approximated by a sequence of weighted 1 norms with the weights depending on the parameter Accordingly, model (5) is approximated by the model depending on the parameter as follows:
min f ∈C
m
F
i =1
Sγ ( W f ) i Sγ ( W f ) i .
and u.
(7)
That is, model (5) is approximated by the weighted 1 model with a nonlinear weight. The nonlinearity of the weight in model (7) may be treated numerically by iteration. This leads to a heuristic method for solving model (7): 1. Initialization: choose an initial guess f (0) . 2. Updating f (k+1) :
f (k+1) ∈ arg min f ∈C
m i =1
F
Sγ W f (k) i Sγ ( W f ) i .
(8)
3. Set k ← k + 1; Repeat step 2 until some stopping criteria are satisfied. We see that the main step of the above procedure is to solve the convex optimization problem (8). We call it the core optimization problem for inpainting. This core optimization problem fits into the general setting studied recently in [12,24]. We remark that if the soft-thresholding operator in (8) is replaced by the hard-thresholding operator, then the optimization problem (8) becomes
arg min f ∈C
m i =1
F
Hγ W f (k) i Hγ ( W f ) i ,
whose objective function is nonconvex. It turns out that the corresponding optimization problem is still difficult to solve. This is the reason why the soft-thresholding is adopted in (5). 3. Solutions of the core optimization problem In this section, we investigate the existence of solutions of the core optimization problem and describe a fixed-point based algorithm to solve it. The existence of the solution for the core optimization problem (8) is essentially determined by the wavelet matrix W appeared in (3). We begin with a description of the construction of matrix W from a set of finite-length filters {τ : = 0, 1, . . . , L } that associate with a tight framelet system. The Fourier series of each filter τ is a trigonometric polynomial and τ . We assume that τ0 is low-pass, that is, τ0 (0) = 1, and the remaining filters are high-pass. For each ω ∈ R, is denoted by these filters satisfy the conditions
τ0 (ω)τ0 (ω + ρπ ) +
L
τ (ω)τ (ω + ρπ ) = δρ ,0 ,
=1
where
ρ = 0, 1. In what follows, we view all filters τ as row vectors. We define recursively
τ,1 = τ ,
τ,k := τ,k−1 ⊗ [1 0], k 2.
(9)
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
175
These vectors τ,k will be used to form the multi-level framelet decomposition matrix B in (3). For simplicity, we consider only the one-dimensional case, which can be extended to the two-dimensional case straightforwardly. Suppose that the dimension of the vector τ,k is smaller than n. For each k, we generate an n × n circulant matrix T ,k whose first row is the vector τ,k with augmentation by zeros from its right end. For a K -level framelet decomposition, the matrices B 0 and W mentioned in (3) are defined by
⎡
B0 =
K −1
T 0, K −k
k =0
T 1, K
K −1 k =1
T 0, K −k
⎤
.. ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ K −. 1 ⎥ ⎢ T L,K T k=1 0, K −k ⎥ ⎢ ⎥ ⎢ . and W = ⎢ .. ⎥. ⎥ ⎢ ⎥ ⎢ T 1,1 ⎥ ⎢ ⎥ ⎢ .. ⎦ ⎣ .
(10)
T L ,1 The next technical lemma states a crucial property of the matrix W for establishing the existence of a solution of the core optimization problem. Lemma 3.1. Let τ0 be the low-pass filter and τ , 1 L, be the high-pass filters associated with a tight framelet system. If the | τ0 (ω)| = 1 occurs only at ω = 0 on the interval [0, π ], then the matrix W after deleting any one of its columns has rank n − 1. Proof. Set H := [ T 1,1 · · · T L,1 ] . It suffices to show that the matrix H after deleting any one of its columns has rank n − 1. Property (9) implies that
T 0,1 T 0,1 + H H = I
(11)
and the eigenvalues of T 0,1 are all less than or equal to 1. Moreover, 1 is the simple eigenvalue of both T 0,1 and T 0,1 T 0,1 due to the assumption that | τ0 (ω)| = 1 occurs only at ω = 0 on the interval [0, π ]. We denote by 1 the all ones n-dimensional vector. Then the vector 1 is the eigenvector of both T 0,1 and T 0,1 T 0,1 corresponding to eigenvalue 1 because τ0 is the low-pass filter. By (11), the rank of the matrix H is at most n − 1. Next, we show that the rank of the matrix H after deleting any one of its columns is exactly n − 1. Suppose this statement is false and the last column of H is removed. Then there exists a nonzero vector v ∈ Rn−1 such that H [ v 0] = 0. This together with identity (11) ensures that the vector [ v 0] is the eigenvector of T 0,1 T 0,1 corresponding to eigenvalue 1.
This contradicts the assumption that 1 is the simple eigenvalue of T 0,1 T 0,1 .
2
We remark that the conditions in Lemma 3.1 imposed on the low-pass filter τ0 are satisfied by many well-known tight framelets. Examples of this type include the compactly supported orthogonal Daubechies wavelets of vanishing moment order q (see [16]),
τ0 (ω) = cos2q (ω/2)
q −1 2q − 1 j =0
j
sin2 j (ω/2) cos2(q−1− j ) (ω/2)
and the spline framelets (see [32]),
τ0 (ω) = e i jω/2 cosq (ω/2),
where j = 0 when q is even , j = 1 when q is odd.
It is straightforward to check that 0 | τ0 (ω)| 1 for all ω ∈ [0, π ] and | τ0 (ω)| equals to 1 only at examples can be found in [15,17]. Essentially, the core optimization problem (8) has the form
min Ξ Sγ ( W f ) 1 , f ∈C
ω = 0. Many other
(12)
where Ξ is a diagonal matrix with positive diagonal entries. Let ξi be the i-th diagonal entry of Ξ . In particular, when ξi is set to F ((Sγ ( W f (k) ))i ), then problem (12) reduces to problem (8). Note that the objective function of problem (12) is · 1 ◦ Ξ ◦ Sγ ◦ W . To establish the existence of the solutions to problem (12), we need the following technical lemma. Lemma 3.2. For an m × n matrix Q , an m × m diagonal matrix Ξ with positive diagonal entries, and a vector γ in Rm with all positive components, if · 1 ◦ Q is coercive, then so is · 1 ◦ Ξ ◦ Sγ ◦ Q .
176
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
Proof. For all f ∈ Rn the following inequality holds:
Ξ Sγ ( Q f ) Ξ Q f 1 − Ξ γ 1 ξmin Q f 1 − Ξ γ 1 , 1 where ξmin is the smallest diagonal entry of Ξ . Since · 1 ◦ Q is coercive, that is, Q f 1 → +∞ as f 2 → +∞, hence, by the above inequality Ξ Sγ ( Q f )1 goes to infinity as well. 2 The constraint set C for the inpainted image that appears in the constrained optimization problem (12) may present computational difficulties in finding its solution. We convert the constrained optimization problem (12) to an unconstrained problem with the help of the indicator function. Specifically, we let ιC denote the indicator function of the closed convex set C , which is defined at f as
ιC :=
if f ∈ C ;
0,
+∞, else.
Then, the constrained problem (12) is written as an equivalent unconstrained form
min ιC ( f ) + Ξ Sγ ( W f ) 1 : f ∈ Rn .
(13)
Next, we establish the existence of a solution of the optimization problem (13). Proposition 3.3. Let τ0 be the low-pass filter and τ , 1 L, be the high-pass filters associated with a tight framelet system, and let C be defined by (2). If | τ0 (ω)| = 1 occurs only at ω = 0 on the interval [0, π ] and the matrix W is constructed by (10), then for any m × m diagonal matrix Ξ with positive diagonal elements, the optimization problem (13) has a minimizer. Proof. It suffices to show that the objective function of the optimization problem (13) is convex and coercive. First, notice that the function max{| · | − a, 0} for any fixed non-negative a is convex on R. Then m Ξ Sγ ( W f ) = ξ j max ( W f ) j − γ j , 0 1 j =1
is a convex function of f on Rn . Next, we show that the objective function is coercive. By Lemma 3.2 we only need to prove that · 1 ◦ W + ιC is coercive on Rn . Clearly, if f is not in C , then W f 1 + ιC ( f ) = +∞. Therefore, we consider only the case f ∈ C . In this case,
W f 1 + ιC ( f ) = W f 1 . By W Λ and W Λc we denote the submatrices extracted from W with column indices in Λ and Λc , respectively. By f Λ and f Λc we denote the subvectors extracted from f with component indices in Λ and Λc , respectively. Hence, f Λ is a constant vector and W Λ f Λ 1 is a constant. Notice that
f 22 = f Λ 22 + f Λc 22 and
W f 1 W Λc f Λc 1 − W Λ f Λ 1 . Further, f 2 → +∞ implies that f Λc 2 → +∞. By Lemma 3.1, matrix W Λc is of full rank. Hence, W Λc f Λc 1 approaches infinity as f Λc 2 goes to infinity. Therefore, W f 1 approaches infinity as f 2 goes to infinity. We conclude that · 1 ◦ W + ιC is coercive on Rn . This completes the proof of this proposition. 2 Proposition 3.3 with the property of the function F ensures that for each k the core optimization problem (8) has a solution. In the remaining part of this section, we characterize a solution of problem (8) in terms of a system of fixed-point equations based on which we describe a fixed-point based algorithm for solving problem (8). We consider a somewhat general minimization problem. By Γ0 (Rd ) we denote the class of all lower semicontinuous convex functions h : Rd → (−∞, +∞] such that dom(h) := {u ∈ Rd : h(u ) < +∞} = ∅. Suppose both ϕ and ψ are in Γ0 and we consider the convex minimization problem:
min
ϕ ( f ) + ψ( W f ): f ∈ Rn .
(14)
Apparently, with ϕ := ιC and ψ(·) := Ξ Sγ (·)1 the optimization problem (8) (or its equivalent problem (13)) can be viewed as a special case of optimization problem (14). In the following discussion, we assume that problem (14) has a solution. Many algorithms were developed for solving problem (14) from different perspectives in the literature. For example, the alternating direction method of multipliers (ADMM) [20] based approaches make use of variable splitting to
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
177
treat the term ψ( W f ). Other techniques, namely Bregman and split Bregman methods [8,21,35], and Douglas–Rachford splitting [3,19] are closely related to ADMM. Fixed-point based algorithms for (14) were developed recently in [24,30,31]. We next describe a fixed-point based algorithm for solving problem (14). For this purpose, we present a characterization of the solutions of optimization problem (14) in terms of a system of fixed-point equations via the proximity operators of the functions ϕ and ψ . Recall that the proximity operator of h, denoted by proxh , is defined by
proxh (x) := arg min
1 2
u − x22 + h(u ): u ∈ Rd
(15)
for all x ∈ Rd . The conjugate of h ∈ Γ0 (Rd ) is the function h∗ ∈ Γ0 (Rd ) defined at y ∈ Rd by h∗ ( y ) := sup{ x, y − h(x): x ∈ Rd }. In the next theorem, we characterize solutions of optimization problem (14). Theorem 3.4. Let ϕ ∈ Γ0 (Rn ), ψ ∈ Γ0 (Rm ), and W be an m × n matrix. If f ∈ Rn is a solution of problem (14), then for any λ, β > 0 there exists a vector y ∈ Rm such that
f = proxλϕ f − λ W y ,
(16)
y = proxβψ ∗ ( y + β W f ).
(17)
Conversely, if there exist λ, β > 0, f ∈ R , and y ∈ R satisfying Eqs. (16) and (17), then f is a solution of problem (14). n
m
A proof of the above theorem may be found in [24,25]. Our purpose is to design an algorithm for solving problem (14) based upon the characterization of solutions to the problem presented in Theorem 3.4. With x = ( f , y ) ∈ Rn × Rm , define
Ψ (x) := λϕ ( f ) + βψ ∗ ( y ) and
E :=
−λ W
In βW
Im
(18)
(19)
.
Note that
proxΨ (x) = proxλϕ ( f ), proxβψ ∗ ( y ) ∈ Rn × Rm . Then Eqs. (16) and (17) can be refreshed into a single equation as follows:
x = proxΨ ( Ex).
(20)
We further note that E 2 , the spectral norm of E (i.e., the largest singular value of E) is strictly greater than 1 for any nonzero matrix W . As a consequence, the operator proxΨ ◦ E may not be nonexpansive. Hence, we cannot seek a solution to problem (14) via a simple iteration
xk+1 = proxΨ Exk
for any initial estimate x0 ∈ Rn+m because convergence of the resulting sequence {xk : k ∈ N} cannot be guaranteed. To overcome the difficulty caused by the expansiveness of the matrix E, we split E as the sum of M and E − M for a proper (m + n) × (m + n) matrix M. An appropriate splitting of matrix E may lead to a convergent iteration scheme
xk+1 = proxΨ Mxk+1 + ( E − M )xk ,
(21)
for an initial estimate x0 ∈ Rn+m . The corresponding limit is a solution of the fixed-point equation (20). There are many ways to choose the matrix M. In general, iteration (21) is an implicit scheme. Systematic study of the implicit scheme (21) in a general context will be investigated in a different occasion. Here, we give two choices of the matrix M which leads to explicit iterations. Specifically, when we choose
0 M= 2β W
0 , 0
the iterative scheme (21) becomes
f k+1 = proxλϕ f k − λ W yk ,
yk+1 = proxβψ ∗ yk + β W 2 f k+1 − f k .
(22)
178
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
When we choose
M=
0 −2λ W , 0 0
the iterative scheme (21) becomes
yk+1 = proxβψ ∗ yk + β W f k ,
f k+1 = proxλϕ f k − λ W 2 yk+1 − yk .
(23)
These two iterations in fact coincide with the primal–dual algorithms which were originally developed from a different point of view [12] applied to problem (14). It was known (see, e.g., [12,25,36]) that if λβ < 1/ W 22 , then the sequence { f k : k ∈ N} generated from (22) or from (23) converges to a solution of problem (14). 4. An algorithm for solving the ASIM In this section, we discuss crucial issues for implementation of the heuristic method for the ASIM (5) in the context of inpainting. First of all, for the inpainting model to be efficient we need to choose an appropriate transform W . Moreover, the main step of the heuristic method is to find a solution of the optimization problem (8), which can be solved by the proposed fixed-point equation based algorithms via proximity operators. They are the iterative scheme (22) or (23). The implementation of these algorithms requires that the proximity operators of both λϕ and βψ ∗ are specifically given. We shall present specific formula for computing the proximity operators of convex functions λϕ and βψ ∗ . We discuss the choice of the transform W . The piecewise cubic tight framelet system (see [32]) will be adopted in the our model. The same system is also used in the models which will be used for comparison. The filters associated with the system are
⎧ ⎪ τ0 = ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ τ1 = ⎪ ⎪ ⎪ ⎪ ⎨ τ2 = ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ τ3 = ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩τ = 4
1 [1, 4, 6, 4, 1], 16 1 [1, 2, 0, −2, −1], 8
√
6
16 1
[−1, 0, 2, 0, −1],
[−1, 2, 0, −2, 1],
8 1
16
[1, −4, 6, −4, 1].
The transform matrix W in (5) is then built from these filters. The approximate sparsity model involves the function F which approximates the function F . We now choose
1
F (x) =
,
if 0 |x| ;
1 |x| ,
if |x| > .
Clearly, F is continuous and it converges to F pointwise as goes to zero. We remark that the function F (x) = 1/(|x| + ) was used in [11]. For the inpainting model, we are required to compute the proximity operator of the two functions
ϕ := ιC and ψ(·) := · 1 ◦ Ξ ◦ Sγ .
(24)
Here C is the set in which the inpainted image lies, Ξ is a diagonal matrix having positive diagonal entries determined by the approximate function F and γ is a vector with non-negative components. In the heuristic method, the function ϕ is kept unchanged while the function ψ changes iteratively via updating the matrix Ξ . One can easily check that for f ∈ Rn
proxλϕ ( f ) = P Λ g + ( I − P Λ ) f . Note that the right-hand side of the formula above is independent of λ due to the fact of proxλϕ = proxϕ . Recall that P Λ is an n × n diagonal matrix whose i-th diagonal entry is 1 if i ∈ Λ and 0 otherwise. We now derive the proximity operator of the function ψ ∗ . For any y ∈ Rm , one has that
ψ( y ) =
m i =1
ξi · max | y i | − γi , 0 .
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
179
Fig. 1. Plots of (a) ξ, , (b) ∗ξ, , and (c) prox∗ . ξ,
For given non-negative numbers ξ and , we define a function ξ, at x ∈ R by
ξ, (x) := ξ · max |x| − , 0 . We then have that ψ( y ) =
ψ ∗ ( y) =
m
m
i =1 ξi ,γi ( y i )
which is separable. Therefore,
∗ξi ,γi ( y i )
i =1
and
proxβψ ∗ ( y ) = proxβ∗
ξ1 ,γ1
( y 1 ), proxβ∗ξ
2 ,γ2
( y 2 ), . . . , proxβ∗ξ
m ,γm
( ym ) .
Thus, the proximity operator of βψ ∗ can be obtained with the help of the following lemma. Lemma 4.1. For two non-negative numbers ξ and , we have at x ∈ R that
∗ξ, (x) =
|x|,
if |x| ξ ;
(25)
+∞, otherwise.
Further, for any positive number β , we have that
⎧ if |x| > β + ξ ; ⎨ sign(x)ξ, proxβ∗ (x) = x − sign(x)β , if β |x| β + ξ ; ξ, ⎩ 0, otherwise.
(26)
Proof. By the definition of the conjugate function, we have that ∗ξ, (x) = sup{xy − ξ, ( y ): y ∈ R} at x ∈ R. From the definition of ξ, we observe that
xy − ξ, ( y ) =
xy ,
if | y | ;
xy − ξ(| y | − ),
otherwise.
Hence, for the fixed x, the supremum of xy − ξ, ( y ) is |x| for | y | ; or +∞ if | y | > and |x| > ξ ; or |x| if | y | > and |x| ξ . Therefore, Eq. (25) holds. Now, we prove the formula (26). For x ∈ R, we write y = proxβ∗ (x). We have that x − y ∈ ∂(β∗ξ, )( y ) by using Propoξ,
sition 2.6 in [30], where ∂ denotes the subdifferential. Since ∂(β∗ξ, ) = β∂∗ξ, , we get (x − y )/β ∈ ∂∗ξ, ( y ). On the other hand, for any y ∈ dom(∗ξ, ) = [−ξ, ξ ] we have that (it is also clear from Fig. 1(b))
⎧ (−∞, − ], ⎪ ⎪ ⎪ ⎪ − ⎪ ⎨ , ∗ ∂ξ, ( y ) = [− , ], ⎪ ⎪ ⎪ , ⎪ ⎪ ⎩ [ , ∞),
y = −ξ ;
−ξ < y < 0; y = 0; 0 < y < ξ; y = ξ,
which, together with the inclusion (x − y )/β ∈ ∂∗ξ, ( y ), yields Eq. (26).
2
180
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
Algorithm 1 (Iterative scheme for the ASIM (5)). Input: the m × n matrix W ; positive numbers λ, β , τ , and α with λβ < 1 and 0 < α < 1; positive integers kmax and j max ; 0 0 Initialization: choose 0 > 0, γ ∈ Rm + , f , y , diagonal matrix Ξ with positive diagonal entries, and ψ = · 1 ◦ Ξ ◦ Sγ ; set k = 0. while k < kmax do Step 1: Compute f k+1 : Initialize x0 = f k , u 0 = yk , and j = 0 while j < j max and x j − x j −1 2 > τ x j 2 do
⎧ j +1 ← P Λ g + ( I − P Λ )(x j − λ W u j ) ⎪ ⎨x j +1 u ← proxβψ ∗ (u j + β W (2x j +1 − x j )) ⎪ ⎩ j← j+1 end while Step 2: Update f k+1 ← x j , yk+1 ← u j ,
k+1 ← α k , and ψ ← · 1 ◦ Ξ ◦ Sγ , where Ξ = diag( F k+1 ( w 1 ), F k+1 ( w 2 ), . . . , F k+1 ( w m )) with w = Sγ ( W f k+1 ) and w i being the i-th component of w for i = 1, 2, . . . , m. Step 3: k ← k + 1. end while
We make two remarks on the function ξ, . First, the function ξ, (see Fig. 1(a)) is the -insensitive loss function scaled by ξ that is widely used in support vector machines and regression problems [33]. In viewing model (14) with ϕ and ψ given by (24), the i-th framelet coefficient ( W f )i has no cost if it lies in the γi -neighborhood of the origin; otherwise, its contribution will be weighted by the parameter ξi in the objective function of model (14). Second, in the proximity operator of β∗ξ, at x ∈ R, the number β acts as a threshold. The operator proxβ∗ (see Fig. 1(c)) performs as ξ,
a thresholding operator with β as its threshold on the interval [−(β + ξ ), (β + ξ )]; when x is outside this interval, the proxβ∗ (x) saturates at ξ if x is positive; and −ξ otherwise. ξ,
To develop a convergent algorithm for solving problem (8), we also need to know the 2 norm of matrix W . The size of W is m × n with m = 24nK , where n is the number of pixels of the underlying image and K is the decomposition level. The 2 norm of W is always equal to 1 for any K . Thus, in iterative scheme (22) or (23), we choose the parameters λ and β such that λβ < 1. In the experiments to be presented in Appendix A, we specifically choose K = 1. Now, we are ready to present a complete and implementable algorithm, as described in Algorithm 1, for solving the ASIM (5). There are two “while” loops in the algorithm. The outer-loop is to update ψ and rescale by a factor α < 1 in m each iteration. Changing the value of parameter is mainly to enforce i =1 F (u i )|u i | to be close enough to u 0 . The maximum number of outer-loops is kmax . The inner-loop is to solve the core optimization problem (8) by using the iterative scheme (22). The maximum number of inner-loops is set to be j max . The algorithm for solving the ASIM by using the iterative scheme (23) can be given accordingly. Appendix A. Examples In this appendix, we present numerical results to demonstrate the performance of our proposed Algorithm 1 for image inpainting and compare it with Algorithm FBIA [8] since other existing algorithms either have a numerical behavior similar to that of FBIA or are not comparable to FBIA. The objective quality of restored images is evaluated quantitatively by the peak signal-to-noise ratio (PSNR). All algorithms are implemented in Matlab and executed on MATLAB R2010a running on a PC equipped with an Intel Core 2 Quad CPU at 2.2.00 GHz and 2G RAM memory. The grayscale images of “Cameraman” and “House” of size 256 × 256 and the color images of “MonaLisa” of size 330 × 450 and “Madonna&Child” of size of 880 × 1200 shown in Fig. 2 are used as our test images. The grayscale images are superimposed with “Text 1” (see the first column of Fig. 3) and “Text 2” (see the first column of Fig. 5) while the color images with scratches are displayed in the first column of Fig. 7. The inpainting problems that we consider here are to restore these images overlying with texts or having scratches by the proposed Algorithm 1 and the FBIA [8]. For Algorithm 1, we fix λ = 0.5, α = 0.8, γi = 0.25 for i = 1, 2, . . . , 24n, and τ = 0.0005 for all the experiments. To ensure the convergence of Algorithm 1, we choose β = 0.999/λ. The initial points for all the experiments are chosen as f 0 = g, y 0 = 0, 0 = 1, and Ξ = diag(0n×n , 20I 24n×24n ). Moreover, for all experiments on both images of “Cameraman” and “House”, we choose ( jmax , kmax ) = (30, 4) for experiments on “Text 1” and ( jmax , kmax ) = (60, 7) for those on “Text 2”. For a color image, we apply Algorithm 1 for its red, green, and blue channels, separately. For the image of “MonaLisa”, we set jmax = 30 for all its red, green, and blue channels, while kmax = 3, 5, 7 for its red, green, and blue channels, respectively. For the image of “Madonna&Child”, we set ( jmax , kmax ) to be (20, 6), (70, 3), and (12, 8) for its red, green, and blue channels respectively. For FBIA, the parameters are determined experimentally for the restored image to achieve the best possible PSNR-value. We present the numerical results in one table and six figures. Table 1 lists the PSNR values of the inpainted images from model (5) by using Algorithm 1 and from model (4) by using Algorithm FBIA. It clearly indicates that model (5) with Algorithm 1 outperforms model (4) with Algorithm FBIA. Fig. 3 shows the visual comparison of the inpainted images from models (5) and (4) for the images of “Cameraman” and “House” overlaid by “Text 1”. Although it can be observed that both models inpaint the images well, differences in visual quality of the restored images can be observed from the zoom-in images shown in Fig. 4. One can see that model (5) with Algorithm 1 can efficiently remove the superposed words over the regions having cartoon components. Although the words in “Text 2” are denser and thicker than those in “Text 1”, from
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
181
Fig. 2. The original images.
Fig. 3. The results of recovering images overlying with “Text 1”: (a) images superposed by “Text 1”; (b) restored images by model (4) with FBIA; (c) restored images by model (5) with Algorithm 1. The zoomed-in part of each image is marked by a square in the corresponding image that will be highlighted in Fig. 4.
Fig. 4. Column 1: Zoom-in portions of the original test images; Column 2: Zoom-in portions of the original test images with “Text 1”; Column 3: inpainted images by model (4) with FBIA; Column 4: inpainted images by model (5) with Algorithm 1.
182
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
Table 1 PSNR (dB) of restored images by model (5) with Algorithm 1 and model (4) with FBIA. Methods
Algorithm 1 FBIA
“Text 1”
“Text 2”
“Scratch”
Cameraman
House
Cameraman
House
MonaLisa
Madonna&Child
30.91 30.29
38.29 37.37
27.32 26.83
33.21 31.51
46.60 45.45
46.62 46.24
Fig. 5. The results of recovering images overlying with “Text 2”: (a) images superposed by “Text 2”; (b) restored images by model (4) with FBIA; (c) restored images by model (5) with Algorithm 1. The zoomed-in part of each image is marked by a square in the corresponding image that will be highlighted in Fig. 6.
Fig. 6. Column 1: Zoom-in portions of the original test images; Column 2: Zoom-in portions of the original test images with “Text 2”; Column 3: inpainted images by model (4) with FBIA; Column 4: inpainted images by model (5) with Algorithm 1.
Figs. 5 and 6 we can draw the same conclusion for the images having “Text 2”. The restored color images are shown in Figs. 7 and 8. Again, these figures show that model (5) with Algorithm 1 outperforms model (4) with Algorithm FBIA. In the spirit of reproducible research, the Matlab source codes for reproducing the results in this section can be obtained by sending an email to [email protected].
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
183
Fig. 7. The results of recovering images overlying with scratch: (a) images superposed by scratch; (b) restored images by model (4) with FBIA; (c) restored images by model (5) with Algorithm 1. The zoomed-in part of each image is marked by a square in the corresponding image that will be highlighted in Fig. 8.
Fig. 8. Column 1: Zoom-in portions of the original test images; Column 2: Zoom-in portions of the original test images with “Text 3”; Column 3: inpainted images by model (4) with FBIA; Column 4: inpainted images by model (5) with Algorithm 1.
References [1] M. Bertalmio, G. Sapiro, V. Caselles, C. Ballester, Image inpainting, in: Proceedings of SIGGRAPH, New Orleans, LA, 2000, pp. 417–424. [2] M. Bertalmio, L. Vese, G. Sapiro, S. Osher, Simultaneous structure and texture image inpainting, IEEE Trans. Image Process. 12 (2003) 882–889. [3] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, Distributed optimization and statistical learning via alternating direction method of multipliers, Found. Trends Mach. Learn. 3 (1) (2010) 1–122. [4] A. Bugeau, M. Bertalmio, V. Caselles, A comprehensive framework for image inpainting, IEEE Trans. Image Process. (2010) 2634–2645. [5] J.-F. Cai, R. Chan, L. Shen, Z. Shen, Simultaneously inpainting in image and transformed domains, Numer. Math. 112 (2009) 509–533. [6] J.-F. Cai, R. Chan, Z. Shen, A framelet-based image inpainting algorithm, Appl. Comput. Harmon. Anal. 24 (2007) 131–149. [7] J.-F. Cai, R.H. Chan, L. Shen, Z. Shen, Convergence analysis of tight framelet approach for missing data recovery, Adv. Comput. Math. 31 (2009) 87–113. [8] J.-F. Cai, S. Osher, Z. Shen, Split Bregman methods and frame based image restoration, Multiscale Model. Simul. 2 (2009) 337–369.
184
L. Shen et al. / Appl. Comput. Harmon. Anal. 37 (2014) 171–184
[9] E. Candes, J. Romberg, T. Tao, Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information, IEEE Trans. Inform. Theory 52 (2006) 489–509. [10] E. Candes, J. Romberg, T. Tao, Stable signal recovery from incomplete and inaccurate measurements, Comm. Pure Appl. Math. 59 (2006) 1207–1223. [11] E. Candes, M.B. Wakin, S.P. Boyd, Enhancing sparsity by reweighted 1 minimization, J. Fourier Anal. Appl. 14 (2008) 877–905. [12] A. Chambolle, T. Pock, A first-order primal-dual algorithm for convex problems with applications to imaging, J. Math. Imaging Vision 40 (2011) 120–145. [13] T. Chan, S.-H. Kang, J. Shen, Euler’s elastica and curvature-based image inpainting, SIAM J. Appl. Math. 63 (2002) 564–592. [14] T. Chan, J. Shen, Mathematical models for local nontexture inpaintings, SIAM J. Appl. Math. 62 (2002) 1019–1043. [15] C.K. Chui, W. He, J. Stockler, Compactly supported tight and sibling frames with maximum vanishing moments, Appl. Comput. Harmon. Anal. 13 (2002) 224–262. [16] I. Daubechies, Ten Lectures on Wavelets, CBMS-NSF Regional Conf. Ser. in Appl. Math., vol. 61, SIAM, Philadelphia, 1992. [17] I. Daubechies, B. Han, A. Ron, Z. Shen, Framelets: MRA-based constructions of wavelet frames, Appl. Comput. Harmon. Anal. 14 (2003) 1–46. [18] D. Donoho, M. Vetterli, R. DeVore, I. Daubechies, Data compression and harmonic analysis, IEEE Trans. Inform. Theory 44 (1998) 2435–2476. [19] J. Eckstein, D.P. Bertsekas, On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators, Math. Program. 55 (1992) 293–318. [20] D. Gabay, B. Mercier, A dual algorithm for the solution of nonlinear variational problems via finite-element approximations, Comput. Math. Appl. 2 (1976) 17–40. [21] T. Goldstein, S. Osher, The split Bregman method for 1 regularization problems, SIAM J. Imaging Sci. 2 (2009) 323–343. [22] N. Komodakis, G. Tziritas, Image completion using efficient belief propagation via priority scheduling and dynamic pruning, IEEE Trans. Image Process. 16 (2007) 2649–2661. [23] A. Krol, S. Li, L. Shen, Y. Xu, Preconditioned alternating projection algorithms for maximum a posteriori ECT reconstruction, Inverse Problems 28 (2012) 115005, 34pp. [24] Q. Li, C.A. Micchelli, L. Shen, Y. Xu, A proximity algorithm accelerated by Gauss–Seidel iterations for L1/TV denoising models, Inverse Problems 28 (2012) 095003, 20 pp. [25] Q. Li, L. Shen, Y. Xu, N. Zhang, Multi-step fixed-point proximity algorithms for solving a class of regularization problems arising from image processing, UCLA CAM Report 12-65, 2012. [26] Q. Li, L. Shen, L. Yang, Split-Bregman iteration for framelet based image inpainting, Appl. Comput. Harmon. Anal. 32 (2012) 145–154. [27] Y. Lu, L. Shen, Y. Xu, Multi-parameter regularization methods for high-resolution image reconstruction with displacement errors, IEEE Trans. Circuits Syst. I, Fundam. Theory Appl. 54 (2007) 1788–1799. [28] A. Maalouf, P. Carre, B. Augereau, C. Fernandez-Maloigne, A bandelet-based inpainting technique for clouds removal from remotely sensed images, IEEE Trans. Geosci. Remote Sens. 47 (2009) 2363–2371. [29] S. Mallat, A Wavelet Tour of Signal Processing, 2nd ed., Academic Press, New York, 1999. [30] C.A. Micchelli, L. Shen, Y. Xu, Proximity algorithms for image models: Denoising, Inverse Problems 27 (2011), 045009 (30 pp). [31] C.A. Micchelli, L. Shen, Y. Xu, X. Zeng, Proximity algorithms for the L1/TV image denoising model, Adv. Comput. Math. 38 (2013) 401–426. [32] A. Ron, Z. Shen, Affine system in L 2 ( R d ): the analysis of the analysis operator, J. Funct. Anal. 148 (1997) 408–447. [33] V. Vapnik, The Nature of Statistical Learning Theory (Information Science and Statistics), 2nd ed., Springer-Verlag, New York, 1999. [34] Y. Wexler, E. Shechtman, M. Irani, Space–time completion of video, IEEE Trans. Pattern Anal. Mach. Intell. 29 (2007) 463–476. [35] W. Yin, S. Osher, D. Goldfarb, J. Darbon, Bregman iterative algorithms for 1 minimization with applications to compressed sensing, SIAM J. Imaging Sci. 1 (2008) 143–168. [36] X. Zhang, M. Burger, S. Osher, A unified primal–dual algorithm framework based on Bregman iteration, J. Sci. Comput. 46 (2011) 20–46.