3498

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 24, NO. 11, NOVEMBER 2015

Annihilating Filter-Based Low-Rank Hankel Matrix Approach for Image Inpainting Kyong Hwan Jin and Jong Chul Ye, Senior Member, IEEE

Abstract— In this paper, we propose a patch-based image inpainting method using a low-rank Hankel structured matrix completion approach. The proposed method exploits the annihilation property between a shift-invariant filter and image data observed in many existing inpainting algorithms. In particular, by exploiting the commutative property of the convolution, the annihilation property results in a low-rank block Hankel structure data matrix, and the image inpainting problem becomes a low-rank structured matrix completion problem. The block Hankel structured matrices are obtained patch-by-patch to adapt to the local changes in the image statistics. To solve the structured low-rank matrix completion problem, we employ an alternating direction method of multipliers with factorization matrix initialization using the low-rank matrix fitting algorithm. As a side product of the matrix factorization, locally adaptive dictionaries can be also easily constructed. Despite the simplicity of the algorithm, the experimental results using irregularly subsampled images as well as various images with globally missing patterns showed that the proposed method outperforms existing state-of-the-art image inpainting methods. Index Terms— Annihilating filter, low rank structured matrix completion, image inpainting, block Hankel matrix, Markov random field, partial differential equation, ADMM.

I. I NTRODUCTION

A

N INPAINTING technique is a method to estimate the missing pixels in an image by exploiting redundant information from adjacent pixels. Applications of inpainting can be found in many scientific and engineering contexts, e.g., for the removal of metal artefacts in x-ray CT [1], restoration of deteriorated images [2], [3], and synthesis of natural texture patterns [4]. In MR applications, the number of k-space samples is directly proportional to the acquisition time. Hence, in order to achieve a faster temporal resolution, k-space down-sampling is often necessary. However, to avoid aliasing artifacts from down-sampling, a high performance interpolation scheme for the missing k-space data is required [5], [6]. Inpainting is also closely related to sampling theory from irregular, non-uniform or scattered samples. Aldroubi and Gröchenig provided a unified theory for

Manuscript received October 27, 2014; revised March 2, 2015 and May 18, 2015; accepted June 8, 2015. Date of publication June 17, 2015; date of current version July 7, 2015. This work was supported by the Korea Science and Engineering Foundation under Grant NRF2014R1A2A1A11052491. The associate editor coordinating the review of this manuscript and approving it for publication was Dr. Lei Zhang. The authors are with the Department of Bio and Brain Engineering, Korea Advanced Institute of Science and Technology, Daejeon 305-338, Korea (e-mail: [email protected]; [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TIP.2015.2446943

uniform and non-uniform sampling in shift-invariant spaces [7]. Another approach is to define new basis functions that are better suited to non-uniform sampling structures. For example, an irregularly sampled function can be interpolated using radial basis functions and thin plate splines [8]. Kybic et al. [9], [10] generalized this approach using a variational framework. Another important recent research trend in inpainting comes from data driven approaches. For example, the kernel regression method [11] has been developed to adapt the local nature of images using data-driven kernels. The authors further reported that there exists an equivalent kernel that produces the cubic-B spline interpolation formulation, and asserted that their kernel regression method is more powerful because it enables a rich class of data-driven kernels to be used. On the other hand, sparsity driven interpolation methods have been extensively studied [12]–[18], among which K-SVD [12], [13] is the most well-known one. Here, the main idea is to find the data-adaptive local basis and fit the data based on the sparsity of the regression coefficients. The way to learn the local basis is also driven by the sparsity of the coefficients, which is now classified as sparse basis learning methods. By extending the idea of K-SVD, patch-based inpainting techniques have been extensively studied [17], [19]–[22]. Patch-based inpainting algorithms try to match the local features of an incomplete patch to a reference’s patches by controlling the scale and orientation [19], [21] or by grouping similar patches into one family to exploit the low-rankness of the similarity group [17], [20], [22]. Yet another important class of classical inpainting methods is based on the statistical and deterministic structure of images. For example, the Gaussian Markov random field (GMRF) model [23]–[26] is often used to model the textures of images in a stochastic sense, which is again used to estimate the missing data [27]. Deterministic generalization has also been extensively studied using a variational framework and partial differential equation (PDE) for a (non)-linear diffusion process [28]–[31]. The main theme in most of the aforementioned inpainting methods is to find an optimized interpolation filter or kernel that can be used to obtain the missing pixels. However, one of the main contributions of this study is to show that this step is not necessary if the interpolation filter is locally spatially invariant. Specifically, we show that the fixed point solution of many inpainting algorithms such as GMRF or anisotropic diffusion provides an annihilating filter [32], [33] relationship at each image patch. Moreover, this phenomenon has an important Fourier domain interpretation in that the annihilating

1057-7149 © 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

JIN AND YE: ANNIHILATING FILTER-BASED LOW-RANK HANKEL MATRIX APPROACH

filter relationship indeed originates from the spectrum sparsity of the local patch. Then, by exploiting the commutative property of convolution, a rank deficient block Hankel structured data matrix can be obtained at each patch from the corresponding annihilating filter relationship. Therefore, the proposed algorithm uses a low rank matrix completion algorithm for the Hankel structured matrix [34], [35], from which the missing samples in the patch can be easily obtained. To adapt to the changes in the local image statistics and avoid any blocking artefacts, interpolation steps are applied on overlapping patches to obtain the missing pixel values, and the estimated values are averaged. Note that our approach is fundamentally different from the existing patch-based approaches. Even though some of the state-of-the-art patch-based inpainting approaches [17], [18], [20] also impose low rankness or sparsity to represent patches with missing pixels, groups of similar patches need to be collected from a neighborhood, from which the local manifold structure is learned. On the other hand, in the proposed method, the manifold structure is estimated internally without using the dictionary obtained from neighborhood patches thanks to the Hankel structured matrix formation. Hence, sophisticated local adaptation schemes based on image statistics [17], [18], [20] to choose the proper dictionary are not required. Despite its simplicity and regular structure without sophisticated adaptation steps, the new algorithm is either comparable or outperforms the state-of-the-art inpainting algorithms both qualitatively and quantitatively. We are aware of an earlier work on texture synthesis and inpainting by a state space model [36], which results in similar low rank Hankel matrix structures. However, this approach uses the Hankel matrix constructed from all images rather than patches due to the difficulty in dealing with the boundary condition. On the other hand, our method is more general because it exploits the annihilation property for each overlapping image patch and can be easily adapted to the local change in image statistics. Another important contribution is that the proposed method not only provides a unified inpainting framework for GMRF and PDE-based inpainting approaches but also for data-driven inpainting algorithms using sparse dictionary learning like the K-SVD. We provide an explicit algorithm that constructs a local dictionary at each patch. The remainder of this paper is as follows. Section II discusses the annihilation structures that are commonly observed in image inpainting algorithms. In Section III, a detailed description of the proposed method is given. Section IV then explains the link between the proposed method and the inpainting approaches using dictionary learning. Experimental results are provided in Section V which is followed by the conclusion in Section VI. II. A NNIHILATING F ILTER S TRUCTURE OF I NPAINTING A. Existing Examples In this sub-section, we review existing inpainting algorithms such as GMRF and PDE-based approaches that

3499

naturally provide an annihilating filter structure within local patches. 1) Inpainting Using Gaussian Markov Random Field Model: A 2D GMRF is described as  gs−r X r + Es , (1) Xs = r∈∂s

where ∂s is the set of nearest neighbors in a 2D grid defined as ∂s = {r = (xr , yr ) : r = s and |xr − x s | ≤ P and |yr − ys | ≤ P},

(2)

where P denotes the size of the neighborhood, and gs denotes the coefficients for a non-causal prediction filter with g0 = 0. The prediction error Es is, in general, not white. Specifically, using the Markovian property E{Et X s } = σ 2 δ(s − t), its correlation matrix Rs,t can be calculated using σ 2 δ(s − t) = E {Et X s }     = E Et gs−r X r + Es =



(3)

r∈∂s 2

σ gs−r δ(r − t) + Rs,t

r∈∂s 2

= σ gs−t + Rs,t .

(4)

Moreover, if we let x = V EC(X), e = V EC(E) where V EC(·) denotes the column major vectorization operator of a matrix, then the probability density function for E becomes   1 p(e) = exp − eT R −1 e /(2π) N/2 |R|1/2 , 2 and using the coordinate transform from e to x shown in (1), we have   1 1 T 1/2 p(X = x) = |B| exp − x Bx , (5) (2π) N/2 2 where the inverse covariance matrix, B, has the elements 1 (δi− j − gi− j ). (6) σ2 Therefore, at the fixed point, the maximum likelihood estimator [27] for a given GMRF should satisfy the following annihilation property: B(i, j ) =

∂ log p(x) =0 ⇒ Bx = 0, (7) ∂x with the Toeplitz matrix B constructed by the shift-invariant annihilating filter kernel (δi− j − gi− j )/σ 2 . 2) Inpainting Using Anisotropic Diffusion: Inpainting approaches using anisotropic diffusion [28] assume that the image is a discretization version of the continuum u(r), r ∈ R2 that satisfies the following heat equation: ∂u = div(g(|∇u|)∇u), u(0) = u 0 , (8) ∂t where g(·) denotes the edge preserving the nonlinear function and u 0 is the initial image. Now, let us define an inner product v, u = v(r)u(r)dr D

3500

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 24, NO. 11, NOVEMBER 2015

where D denotes the image domain. Using the variational formulation, for any test function v(r) that satisfies the Dirichlet boundary condition, we have L(u, v) ≡ v, ∂u/∂t = −g(|∇u|)∇v, ∇u

(9) (10)

where we use divT = −∇ and Green’s theorem. Let u(r) be represented as linear combinations of the basis set {u j (r)} Nj=1 in D: u(r) =

N 

x j u j (r).

j =1

Suppose, furthermore, that the anisotropic diffusion coefficient g(|∇u|) is relatively constant, i.e. g(|∇u|) g0 within a patch domain D (this is another reason why our algorithm uses patch processing). Then, we have L(u, v) = −

N 

x j g0 ∇v, ∇u j .

(11)

At a fixed point, there exists no more diffusion and we have L(u, v i ) = 0 for any test function v i , i = 1, · · · , M. Accordingly, we have T

where x := x 1 · · · x N and ⎡ ∇v 1 , ∇u 1  ⎢ .. B := g0 ⎣ .

∇v M , ∇u 1 

··· .. . ···

patch in Fig. 1(c), where the spectral components are mainly concentrated on the fundamental frequencies of the patterns. In general, if a patch x(r) ∈ R2 has sparse spectral components, we can show that there exists a corresponding annihilating filter in the image domain. More specifically, if the spectrum of an image patch is described by F {x(r)} = x(ω) ˆ =

k−1 

c j δ(ω − ω j ), ω = (ωx , ω y )

(13)

j =0

j =1

Bx = 0

Fig. 1. Spectral contents for (a) a smoothly varying patch, (b) a patch with abrupt changes, and (c) a textured patch.

(12) ⎤ ∇v 1 , ∇u N  ⎥ .. ⎦. . ∇v M , ∇u N 

Among the various choices of v i , we choose test functions such that ∇v i , ∇u j  = h(i − j ). For example, the choice of v i = u i for non-overlapping u i satisfies the condition. Then, B becomes a Toeplitz matrix, and we have the annihilation property that can be converted to a convolution with a shift invariant filter. B. Fourier Domain Interpretation The existence of annihilating filters in GMRF or PDE based inpainting approaches is not just coincidental but has a very important fundamental Fourier domain interpretation related to the sampling theory of the finite rate of innovations (FRI) [32], [33]. For example, an image patch that can be very well represented using GMRF includes the smooth part, textures, etc [23]–[26]. Moreover, it is well known that a smooth varying image patch can be modelled as homogeneous heat diffusion in (8) [28]. In fact, one of the important observations is that an image patch which has smooth transitions or textures is usually sparse in their Fourier domain in Fig. 1(a)-(c). For example, as shown in Fig. 1(a), a smoothly varying patch usually has spectrum content in the low frequency regions, and the other frequency regions have very little spectral components. For the case of an abrupt transition along the edge as shown in Fig. 1(b), the spectral components are mostly localized along the ωx axis. A similar spectrum sparsity can be observed in a texture

where k denotes the number of non-zero spectral components and F (·) denotes the Fourier transform. Then, it is easy to find ˆ an annihilating function h(ω) in the spectrum domain that does not overlap with the support of x(ω), ˆ i.e. ˆ h(ω) x(ω) ˆ = 0.

(14)

This implies the existence of the annihilating filter ˆ in the image domain: h(r) := F −1 {h(ω)} h(r) ∗ x(r) = 0.

(15)

This annihilation relation can be reformulated as matrix-vector formulation: Bx = 0,

(16)

where B is the convolution matrix comprised of filter h, and x is the vectorized patch signal. This reproduces the annihilating filter relationships in (7) and (12). III. PATCH -BASED L OW R ANK H ANKEL M ATRIX C OMPLETION FOR I NPAINTING A. Patch-Based Formulation The annihilation property in (7), (12) and (16) is an extremely useful property, so we are interested in exploiting this for image inpainting. However, image structure and statistics can vary across images, along with the resulting annihilation property. Therefore, we are interested in partitioning images into overlapping patches, within which the same annihilation property holds. Moreover, we can find a dual representation that is more suitable for algorithmic implementation. Specifically, thanks to the Toeplitz structure of B, if we let X ∈ R M×N be an image patch, then we can find a 2D annihilating filter H = [h−n , · · · , hn ] ∈ R(2m+1)×(2n+1) such that the annihilation property holds by representing B as the convolution matrix for H. Note that if H is an annihilating filter, then its convolution with any finite tap filter is also annihilating filter. Hence, we choose the filter tap size for H that is sufficiently larger. However, care needs

JIN AND YE: ANNIHILATING FILTER-BASED LOW-RANK HANKEL MATRIX APPROACH

3501

where patch-based block Hankel structured matrix H {X} is constructed as ⎤ ⎡ H{x1 } H{x2 } · · · H{x2n+1 } ⎢ H{x } H{x3 } · · · H{x2n+2 } ⎥ 2 ⎥ ⎢ ⎥ H {X} = ⎢ .. .. .. ⎥ ⎢ .. . ⎦ ⎣ . . . H{x N } H{x N−2n } H{x N−2n+1 } · · · ∈ R(M−2m)(N−2n)×(2m+1)(2n+1) .

(23)

B. Low-Rank Hankel Matrix Completion Fig. 2.

Region within a patch where the annihilation property holds.

to be taken when constructing the convolution matrix due to the effect from the boundary of a patch in both column and row directions. Accordingly, the annihilation relation (7), (12) and (16) should hold only inside of the patch as shown in Fig. 2 where the adjusted matrix B˜ is defined as: B˜ = C {H}

(17)

where C {H} denotes the 2D chopped convolution matrix: ⎡ ⎤ 0 C{hn } · · · C{h−n } ⎢ ⎥ .. .. .. ⎢ ⎥ . . . ⎢ ⎥ C {H} = ⎢ ⎥ .. .. .. ⎢ ⎥ . . . ⎣ ⎦ 0 C{hn } · · · C{h−n } ∈ R(M−2m)(N−2n)×M N

(18)

and a chopped 1D convolution matrix C{hi } with the filter hi ∈ R2m+1 is defined by ⎤ ⎡ 0 hi (m) · · · hi (−m) ⎥ ⎢ .. .. .. ⎥ ⎢ . . . ⎥ ⎢ C{hi } = ⎢ ⎥ .. .. .. ⎥ ⎢ . . . ⎦ ⎣ 0 hi (m) · · · hi (−m) ∈ R(M−2m)×M .

(19)

(20)

where the overbar denotes the order reversal operator and H{xi } is a Hankel matrix given by ⎡ ⎤ xi (2) · · · xi (2m + 1) xi (1) ⎢ xi (2) xi (3) · · · xi (2m + 2) ⎥ ⎢ ⎥ ⎢ ⎥ .. .. .. .. ⎣ ⎦ . . . . xi (M − 2m) ∈R

xi (M − 2m (M−2m)×(2m+1)

+ 1)

···

xi (M)

.

(21)

Using Eq. (20), the annihilation property can be represented as C {H}VEC(X) = H {X}VEC(H) = 0,

R ANK (H {X}) = r.

(22)

(24)

This relationship leads us to an algorithm that estimates the missing elements in H {X}. More specifically, the missing element estimation problem becomes a low rank structured matrix completion problem: min R ANK (H {X}) subject to X(i, j ) = M(i, j ), ∀(i, j ) ∈ ,

(25) (26)

where M is the measurement matrix with measured values at positions (i, j ) ∈  and 0’s elsewhere. A variety of low rank matrix completion algorithms have been developed [35], [37]–[42] among which convex relaxation using nuclear norm [37], [38] is most well-known. However, the original nuclear norm minimization algorithm does not exploit the Hankel structure of the matrix, and there exist many variations that exploit the low rank structured matrix [34], [43], [44]. Among these, we employ the SVD-free structured rank minimization algorithm [34] with an initialization using the low-rank factorization model (LMaFit) algorithm [35]. These algorithms do not utilize the singular value decomposition (SVD), so the computational complexity can be significantly reduced. Specifically, the algorithm is based on the following observation [45]:

A ∗ =

Now we have the following commutative property: C{h j }xi = H{xi }h j ,

The annihilation structure in (22) implies that there exists null space in H {X}. As will be discussed later, if the patch image X lives in a low dimensional manifold (say, r -dimensional manifold where r < (2m + 1)(2n + 1), we can show that the block Hankel matrix is also low dimensional and its rank is given by

min

U,V:A=UV H

U 2F + V 2F

(27)

Here A ∗ denotes the nuclear norm of the matrix A [37], [38]. Accordingly, we are interested in the nuclear norm minimization under the matrix factorization constraint: min U 2F + V 2F X

subject to H {X} = UV H X(i, j ) = m i j , ∀(i, j ) ∈ ,

(28) (29)

The first constraint (28) can be handled using alternating direction method of multipliers (ADMM) [46] which will be explained soon. The second constraint (29) is a set constraint imposing that X lives on the set C: C = {X ∈ R M×N | X(i, j ) = m i j , ∀(i, j ) ∈ }.

(30)

3502

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 24, NO. 11, NOVEMBER 2015

Therefore, by combining the two constraints, we have the following cost function for ADMM:  1 L(U, V, X, ) := ιC (X) +

U 2F + V 2F 2 μ + H {X} − UV H +  2F 2

(38) (31)

X(k+1) = arg min ιC (X) X

V(k+1)

(k+1)

(32)

(33)

(34) (35)

It is easy to show that the first step can be simply reduced to   X(k+1) = Pc H † U(k) V(k)H − (k) + P M, (36) where P is a projection mapping on the set  and H † corresponds to the Moore-Penrose pseudo-inverse mapping from our block Hankel structure to a patch, which is calculated as  −1 ∗ H . H † = H ∗H

V

(k+1)

H  −1  = μ H {X(k+1) }+ (k) U(k+1) I + μU(k+1)H U(k+1) . (39)

One of the advantages of the ADMM formulation is that each subproblem is simply obtained from (31). More specifically, we have

U(k+1)

U(k+1)

  −1  = μ H {X(k+1) } + (k) V(k) I + μV(k)H V(k)

where ιC (X) denotes the indicator function:  0, if X ∈ C ιC (X) = ∞, otherwise

μ + H {X} − U(k) V(k)H + (k) 2F 2 1 = arg min U 2F U 2 μ + H {X(k+1) } − UV(k)H + (k) 2F 2 1 = arg min V 2F V 2 μ + H {X(k+1) } − U(k+1) V H + (k) 2F 2 = H {X(k+1) } − U(k+1) V(k+1)H + (k)

the closed-form update equations for U and V are given by

Even though the original Hankel matrix H {X} has large dimension, it is important to note that our algorithm using (38) and (39) only require the matrix inversion of k ×k matrix, where k denotes the estimated rank of the Hankel matrix. This significantly reduces the overall computational complexity. Now, the remaining issue is how to initialize U and V. For this, we employed another SVD-free algorithm called low-rank factorization model (LMaFit) [35]. More specifically, for a low-rank matrix Z, LMaFit solves the following optimization problem: 1

UV H − Z 2F 2 subject to Z(i, j ) = M(i, j ), ∀(i, j ) ∈ , min

U,V,Z

(40)

and Z is initialized with H {X}. LMaFit solves a linear equation with respect to U and V to find their updates and relaxes the updates by taking the average between the previous iteration. Moreover, the rank update can be done automatically by detecting the abrupt changes of diagonal elements of QR factorization [35]. Even though the problem (40) is non-convex due to the multiplication of U and V, the convergence of LMaFit to a stationary point was analyzed in detail [35]. However, the LMaFit alone cannot recover the block Hankel structure, which is the reason we use ADMM later to impose the structure. C. Image Partitioning

(37)

Note that the adjoint operator H ∗ {A} adds multiple elements of A and put it back to the patch coordinate, and (H ∗ H )−1 denotes the division by the number of multiple correspondence; hence, the role of the pseudo-inverse is taking the average value and put it back to the patch coordinate. Next, the subproblems for U and V can be easily calculated by taking the derivative with respect to each matrix. For example, the derivative of cost function for U is given by   ∂ 1 ∂L μ 2 H 2 =

U F + H {X} − UV +  F ∂U ∂U 2 2   = U − μ H {X} − UV H +  V   = U I + μV H V − μ (H {X} + ) V, and the closed-form solution of the subproblem for U is obtained by setting ∂ L/∂U = 0. In the similar way, the derivative with respect to V can be obtained. Accordingly,

To adapt the local image statistics, we partition an image into M × N squared patches and exploiting the annihilation property at each patch. As shown in Fig. 2, only the central (M − 2m) × (N − 2n) regions of an image patch X is not effected by the boundary, so when we partition the image, each image patch should be overlapped such that there exist no areas where the annihilation property is not utilised. Moreover, it is better to allow smooth transition of the image structure between patches, so we partition the image into overlapping patches such that the half of the central (M−2m)×(N−2n) are overlapped between two adjacent image patches. The overall flowchart of the patch based processing is illustrated in Fig. 3. IV. D ICTIONARY /M ANIFOLD L EARNING T HROUGH S TRUCTURED L OW-R ANK PATCHES As described in the Introduction, another important class of modern inpainting algorithms is based on dictionary learning that constructs local dictionaries to fit the patch data by exploiting the sparsity of the regression coefficients [12]. In fact, this sparse dictionary learning is closely related to

JIN AND YE: ANNIHILATING FILTER-BASED LOW-RANK HANKEL MATRIX APPROACH

Fig. 3.

3503

The reconstruction flowchart of the proposed annihilating filter based low rank Hankel approach for image inpainting.

localized manifold learning [47] where the local image patches are assumed to be in a low dimensional manifold [48]. In conventional approaches, the manifold structure is learned from similar patches collected from the neighborhood. Because there exists no explicit step to construct dictionaries, the proposed algorithm appears different from the conventional dictionary/manifold learning algorithms; however, the goal of this section is to show that the proposed algorithm indeed learns the low dimensional manifold, and we can explicitly construct the local dictionaries. More specifically, in system identification literature [34], [43], [44], it is well known that a linear time invariant dynamic system can be estimated from the Hankel matrix, and the rank of the Hankel matrix corresponds to the model order. To adapt this for a 2D image patch, one could use a separable dynamic state space model for each coordinate and identify the 2D image manifold. However, because the image signal is usually non-causal and separable manifold expression often restricts its richness in dictionary representation, here we provide an alternative way of constructing dictionary atoms. More specifically, H {X} = UV H and we have  r   X = H † {UV H } = H † ui viH (41) i=1

where ui and vi denotes the i -th vector of U and V, respectively. Since the Moore-Penrose pseudo-inverse is linear,

Fig. 4. Representative dictionary atoms from 45 × 45 patches of marked regions. In this figure, the values of each atom image are normalized for high contrast visualization.

we have X=

r 

  H † ui viH .

(42)

i=1

This implies that X  can be represented  as  a linear combination of {H † ui viH }ri=1 ; and H † ui viH becomes an (un-normalized) atom of the corresponding dictionary. Note that the proposed dictionary construction method has many unique features. Compared to conventional dictionary

3504

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 24, NO. 11, NOVEMBER 2015

Fig. 5. Rotation invariancy comparison results. First row: data images and the block location. Singular value distribution with respect to a rotation angle for the case of the image patches (second row) and their block Hankel structured matrix (third row). σ 2 is a variance of compression ratios at numerical ranks and μ is an averaged compression ratio at numerical ranks. The black lines correspond to the 5% singular value location. TABLE I C OMPARISON OF N UMERICAL R ANK AND THE C OMPRESSION R ATIO B ETWEEN THE 37 × 37 I MAGE PATCH AND I TS

B LOCK H ANKEL S TRUCTURED M ATRIX OF 625 × 169 S IZE

learning algorithms that use a fixed number of dictionaries, the size of the dictionary in the proposed method can easily be determined patch by patch during the low-rank structured matrix completion. Moreover, our dictionary estimation provides the number of significant atoms that reflect the local image statistics. For example, it is usually small for a flat background or repetitive structure, whereas for more complicated unorganized regions, the estimate is usually big, and the resulting dictionary atoms increase as well. More specifically, Fig. 4 shows the normalized atoms of the dictionaries for a hat, Lena’s eye region and decorations, respectively. In this example, the significant atoms were determined by choosing atoms whose singular values were greater than 5% of the maximum singular value. We can confirm that the flat regions have a much smaller number of atoms. Furthermore, the lower order atom retains the important features of an original patch, whereas the many higher order atoms usually encode the texture structure. We believe that this local adaptive nature of our dictionary learning is another important advantage of the proposed method. As a side remark, note that {H † {ui viH }} are not the basis because the atoms are not necessarily independent. This is

because the pseudo-inverse operator is just a left inverse that has a non-trivial null space. Instead, the estimated dictionary is closely related to the frame representation. V. E XPERIMENTAL R ESULTS A. Low Rank Property Currently, many applications for the low rank matrix completion in image processing are based on the low rankness of an image or image patch itself. However, we assert that there exists much more useful intrinsic low rank representation rather than the image or patch itself, and we claim that it indeed comes from the proposed Hankel structured matrix. To confirm our claim, we analyzed the low rankness of various textures and images. For example, the first three columns in Fig. 5 are from the Brodatz texture set [49], and the last two columns in Fig. 5 are from the USC-SIPI image database. Then, we calculated the SVD of the original image patch X and its block Hankel structure H {X}. The numerical rank is calculated as the number of singular values above 5% of the maximum singular value. In Fig. 5, the black lines correspond to the 5% singular value location. As shown in Fig. 5 and Table I, even though the image patch itself

JIN AND YE: ANNIHILATING FILTER-BASED LOW-RANK HANKEL MATRIX APPROACH

3505

TABLE II H YPER -PARAMETERS U SED IN THE P ROPOSED A LGORITHM . T HE S IZE OF I NTERPOLATION F ILTER I S 13 × 13 AND I TERATION OF ADMM I S 50. F OR LMaF IT, I NITIAL R ANK W ITH I NCREASING S TRATEGY I S O NE AND T OLERANCE I S 1.25 × 10−4

Fig. 6.

Part of reconstructed Barbara images by various methods from 80% missing samples. Missing pixels are in black.

has low rank properties, the compression ratio of the block Hankel structure is much more significant, and the singular values drop much more quickly. Here, the compression ratio is defined as numerical rank . compression ratio := number of columns For example, the block Hankel structured matrix was 625 × 169, which was much larger than the original 37 × 37 patch size, but the numerical rank estimation was still about the same, so the compression ratio was significantly improved. Another technical limitation of imposing a low rank constraint to the image patch itself is that the low rank structure varies depending on the patch orientation. This can be easily confirmed from the results in the second row of Fig. 5, where the singular value distribution of a patch X rapidly changes by varying the rotation angle from 0 to 90 degrees. This implies that the numerical rank of the patch itself is not an intrinsic property of the patches, but it is more determined by a geometric transform. To address this problem, Zhang et al. [50] proposed the so-called TILT (transform invariant low rank texture) algorithm that simultaneously estimates the transform to provide the low rank structure. However, we assert that the proposed block Hankel structure is better in retaining the intrinsic low rank structure of the image patch. Specifically, as shown in the third row of Fig. 5, the singular value distribution

of H {X} does not vary with respect to the rotation and retains the nearly the same values. To quantify this, we calculated the mean (μ) and variance (σ 2 ) of the compression ratios along the rotation in Fig. 5. Here, the numerical rank was again defined as the number of singular values bigger than 5% of the maximum singular value. Both the mean and variance of the compression ratio of patch X are consistently larger than those from H {X}, which clearly indicates that the singular values from the block Hankel structured data matrix is rotationally less sensitive. This rotational invariance property can easily be expected from the derivation of H {X}, which originates from the convolution property. When the image patch is rotated, the filter kernel rotates as well so that the H {X} retains the same annihilation structure. Therefore, we believe that the low rankness of our block Hankel structure matrix is in fact an intrinsic property. B. Image Inpainting Experiments To evaluate the image inpainting performance of the proposed algorithm, we first performed inpainting experiments by randomly removing 50% (or 80%) of the pixels in images. The test sets were Baboon, Barbara, Boat, Cameraman, House, Lena and Peppers images. All test images were rescaled to have values between 0 and 1. The image sizes are given in Table II. The sampling masks were generated once using a

3506

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 24, NO. 11, NOVEMBER 2015

Fig. 7.

Part of reconstructed Lena images by various methods from 80% missing samples. Missing pixels are in black.

uniform random distribution, and the same sampling mask was used for all test images. For comparison, a basic interpolation method (triangular based linear mesh interpolation, MATLAB built-in function ‘griddata’, indicated as Mesh in the figures) was used as the simplest reference algorithm, and state-ofthe-art algorithms such as the Kernel regression [11], IPPO (image processing by patch-ordering) [20], K-SVD [12], and C-SALSA [51] were also used. The IPPO [20] permutes a target image into a smooth 1D signal which was compared with similarity patches. Thanks to the permutation, the patch with the missing pixels can be interpolated with a 1D cubic spline interpolator very effectively using similar patches in the permuted dictionary. The parameters for the comparative algorithms were optimized to have the best performances. More specifically, the Delaunay linear triangular based mesh interpolation has no regularization parameters, so we just initialized the sampled positions from the sampling mask and interpolated scattered data along a uniformly distributed regular grid. In the K-SVD reconstruction, in the case of 50% (resp. 80%) missing samples, the size of the patch was 4×4 (8×8); the number of atoms was set to 256 (512); the number of seeds was set to the total number of partitioned patches, and the sparsity was set to 8 (32) for dictionary learning from the training set. In the IPPO, the same parameters were used as in the inpainting experiments in [20]. In the C-SALSA, μ1 and μ2 were optimally chosen from the range of [0.01, 0.1] and [1, 2], respectively. The tolerance was set to 10−4 . For the kernel regression using a steering kernel, h (global smoothing parameter), the size of the kernel, the size of the local analysis window, λ (regularization for the elongation parameter), and α (structure sensitive parameter) were optimally chosen from the ranges of [0.4, 2], [13, 17], [5, 13], [0.1, 0.5], and [0.1, 0.5], respectively. In the proposed method, the parameters in Table II were used.

Fig. 8. Magnified views of the results in Fig. 6 and 7. (First column) original images, (second column) IPPO results, and (third column) results by the proposed method.

For the evaluation, we used two quantitative metrics: PSNR (peak signal-to-noise ratio) and SSIM (structural similarity) index [52]. When the reference signal (y) is given, the PSNR of the reconstructed image (x) is calculated as  

y ∞ √ , x, y ∈ R N . PSNR(x) = 20 log10 1/ N × y − x 2 On the other hand, the SSIM index in [52] (when α = β = γ = 1 and C3 = C2 /2) was defined as    N 2μxi μyi + C1 2σxi yi + C2 1    , SSIM(x) = N 2 2 σx2i + σy2i + C2 i=1 μxi + μyi + C1 C1 = (K 1 L)2 , C2 = (K 2 L)2 ,

JIN AND YE: ANNIHILATING FILTER-BASED LOW-RANK HANKEL MATRIX APPROACH

3507

TABLE III R ECONSTRUCTION I MAGE PSNR BY VARIOUS I MAGE I NPAINTING A LGORITHMS F ROM 50% AND 80% M ISSING S AMPLES . T HE H IGHEST PSNR IN E ACH I MAGE A RE H IGHLIGHTED W ITH B OLDFACE

TABLE IV R ECONSTRUCTION I MAGE SSIM BY VARIOUS I MAGE I NPAINTING A LGORITHMS F ROM 50% AND 80% M ISSING S AMPLES . T HE H IGHEST SSIM IN E ACH I MAGE A RE H IGHLIGHTED W ITH B OLDFACE

where {K 1 , K 2 }  1 are constants; L is the dynamic range of the pixel values; μ(·) is the mean value of a patch region; σ(·) is the standard deviation of a patch region, and σab is the square root value of the cross covariance between both patch regions with different centers. PSNR and SSIM become higher when the quality of image is enhanced. We summarized the results in terms of these quantitative metrics in Tables III and IV, and some of the reconstructed images and their magnified images for the case of 80% missing data are shown in Figs. 6-8. The proposed method outperforms most of the state-of-the-art algorithms in PSNR and SSIM except in the cameraman and pepper images, where the IPPO had the best PSNR and SSIM values. The cameraman and pepper images have more smoothly-varying image content, in which case the permutation based interpolation scheme in IPPO appears to effectively remove the interpolation noises. However, for general images with textures or patterns, the proposed algorithm mostly outperformed IPPO. Barbara (Fig. 6) and Lena (Fig. 7) were typical examples that

showed noticeable enhancement by the proposed method in both visual quality and quantitative measures. In Figs. 6 and 7, the reconstruction quality by IPPO and the proposed method appeared similar at first glance. However, in the magnified views in Fig. 8, the texture of the weaves and the region of the eye in Barbara were much better reconstructed by the proposed method. Similarly, the texture of the hat in Lena was blurred in the IPPO result, whereas the details of the texture were well reconstructed in the proposed algorithm. To show the usefulness of the proposed method in various inpainting problems, we also performed inpainting experiments with a text inlayed mask (Fig. 9), line shape scratches (Fig. 10), and object removal (Fig. 11). As a basic reference, a triangular-based linear mesh interpolation was used. In addition, state-of-the-art algorithms such as PatchMatch [19] and IPPO [20] were also used for comparative studies. The PatchMatch algorithm tries to match a target patch with a missing region to a near neighborhood patch from the rest of the image through a randomized nearest neighbor algorithm. After the searching process, the incomplete patch

3508

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 24, NO. 11, NOVEMBER 2015

TABLE V H YPER -PARAMETERS U SED IN THE P ROPOSED A LGORITHM FOR T EXT I NLAYED M ASK , L INE S HAPE S CRATCHES AND O BJECT R EMOVAL . T HE μ AND I TERATION N UUMBER OF ADMM I S 100 AND 200, R ESPECTIVELY. F OR LMaF IT, I NITIAL R ANK W ITH I NCREASING S TRATEGY I S O NE AND T OLERANCE I S 1.25 × 10−4

Fig. 10. Reconstructed Lizard images by various methods from line shape scratch damages. Missing pixels are in white.

Fig. 9. Reconstructed Rope images using various methods from text inlayed samples. Missing pixels are in white.

is completed using the matched patch. Therefore, PatchMatch could not be used for the previous 50% and 80% irregular sampled experiments because there is no nearest neighborhood without missing pixels. However, in the following inpainting problems with global missing patterns, PatchMatch is effective because it can find the neighborhood from the surroundings. The test sets were chosen from the Berkeley segmentation dataset [53]. To implement the experiments, PatchMatch was conducted with Adobe Photoshop (built-in operation ‘Content-Aware Fill’), and for the case of object removal, the size of the patch in IPPO had to be increased due to the large size of the removed mask. In the proposed method, the parameters in Table V were used. For the text inlayed experiment, English and Korean characters meaning “rope” were inlayed shown in Fig. 9. The inpainted regions of the proposed method correctly recovered the original texture of the rope whereas the other algorithms produced blurry or unnatural textures. Next, several line scratches were added to the Lizard image shown in Fig. 10. The proposed method reconstructed the textures of the lizard and pebbles accurately, whereas the other algorithms generated blurred patterns noticeable to the naked eye. In the case of the object removal experiment in Fig. 11, due to the large missing area, the size of the annihilation

Fig. 11. Reconstructed Bear images after object removal. Missing pixels are in white.

filter in the proposed method and the size of the patch in the IPPO had to be increased. The result from the Mesh interpolation was significantly distorted. The inpainted image from PatchMatch was somewhat better, but the missing pattern on the background could not be recovered. In this example, the IPPO failed due to the large missing area whose missing information may be difficult to restore using patch reordering. On the other hand, the proposed method could recover the background grass and provide the most natural looking images, even though a relatively small blurring artifact still remained. In Fig. 12, we considered various types of inpainting scenarios: object removal with a textured background, multiple object removals, and color inpainting. The test images were chosen from the Berkeley segmentation dataset [53], and other inpainting related papers [21], [54]. The results clearly show that the proposed method in practice can be used for various applications.

JIN AND YE: ANNIHILATING FILTER-BASED LOW-RANK HANKEL MATRIX APPROACH

3509

Fig. 12. Reconstruction results using the proposed method for various inpainting scenario. The first column represents original images, the second column shows image with missing area and the third column shows the reconstructed images from the proposed method. Missing pixels are in black or white.

VI. C ONCLUSION In this paper, we proposed a low-rank patch-based block Hankel structured matrix approach for image inpainting. The proposed method exploits the annihilation property that is commonly observed in existing inpainting algorithms. We further provided a Fourier domain interpretation of the

annihilating filter relationship. We showed that the annihilating filter relationship implied the existence of a low-rank block Hankel structured matrix that is constructed for each image patch. By demonstrating the rotationally invariant low-rankness of the block Hankel matrix, we confirmed that the proposed low-rank property is really an intrinsic

3510

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 24, NO. 11, NOVEMBER 2015

property of image patches. To exploit the low-rankness for the estimation of missing elements in the structured matrix, we employed the ADMM algorithm that alternatively estimates the factorization matrices, which were initialized using the LMaFit algorithm. Because our algorithm is SVD free, the computational speed was much faster. The proposed method showed performance enhancement in both quantitative and qualitative aspects compared with existing inpainting algorithms. We also demonstrated that a locally adaptive dictionary can be easily obtained from the resulting matrix factorization. One of the drawbacks of the proposed method is the computational complexity due to the large number of matrix inversions. Even though the computational complexity in the ADMM step is determined by the inversion of the matrix whose size is determined by the estimated rank of the Hankel matrix, we need to do the low rank matrix completion for each overlapping patch, which increases the computational complexity. For the case of the object removal, the Hankel matrix size was 3249 × 625 for each patch, and the resulting reconstruction times are about 3, 000s using Matlab implementation on an Intel i7-4770k (3.5GHz) platform. However, with smart optimization of the algorithm (for example, using multi-scale approach to reduce the missing object size at low resolution and refine the results at a finer resolution level using a smaller patch size and GPU implementation), the computational time could be significantly reduced. This direction of research is important which we will address in subsequent publications. Thanks to the unique advantages of the proposed algorithm both theoretically and performance-wise, we believe that our algorithm has significant potential in image inpainting and irregular sampling problems in general. ACKNOWLEDGMENT The authors would like to thank Prof. Michael Unser at EPFL for helpful discussions. R EFERENCES [1] X. Duan, L. Zhang, Y. Xiao, J. Cheng, Z. Chen, and Y. Xing, “Metal artifact reduction in CT images by sinogram TV inpainting,” in Proc. IEEE Nucl. Sci. Symp. Conf. Rec., Oct. 2008, pp. 4175–4177. [2] A. C. Kokaram, R. D. Morris, W. J. Fitzgerald, and P. J. W. Rayner, “Detection of missing data in image sequences,” IEEE Trans. Image Process., vol. 4, no. 11, pp. 1496–1508, Nov. 1995. [3] A. C. Kokaram, R. D. Morris, W. J. Fitzgerald, and P. J. W. Rayner, “Interpolation of missing data in image sequences,” IEEE Trans. Image Process., vol. 4, no. 11, pp. 1509–1519, Nov. 1995. [4] M. Ashikhmin, “Synthesizing natural textures,” in Proc. ACM Symp. Interact. 3D Graph., 2001, pp. 217–226. [5] P. J. Shin et al., “Calibrationless parallel imaging reconstruction based on structured low-rank matrix completion,” Magn. Reson. Med., vol. 72, no. 4, pp. 959–970, Oct. 2014. [6] M. A. Griswold et al., “Generalized autocalibrating partially parallel acquisitions (GRAPPA),” Magn. Reson. Med., vol. 47, no. 6, pp. 1202–1210, Jun. 2002. [7] A. Aldroubi and K. Gröchenig, “Nonuniform sampling and reconstruction in shift-invariant spaces,” SIAM Rev., vol. 43, no. 4, pp. 585–620, 2001. [8] F. L. Bookstein, “Principal warps: Thin-plate splines and the decomposition of deformations,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 11, no. 6, pp. 567–585, Jun. 1989. [9] J. Kybic, T. Blu, and M. Unser, “Generalized sampling: A variational approach—Part I. Theory,” IEEE Trans. Signal Process., vol. 50, no. 8, pp. 1965–1976, Aug. 2002.

[10] J. Kybic, T. Blu, and M. Unser, “Generalized sampling: A variational approach—Part II. Applications,” IEEE Trans. Signal Process., vol. 50, no. 8, pp. 1977–1985, Aug. 2002. [11] H. Takeda, S. Farsiu, and P. Milanfar, “Kernel regression for image processing and reconstruction,” IEEE Trans. Image Process., vol. 16, no. 2, pp. 349–366, Feb. 2007. [12] M. Aharon, M. Elad, and A. Bruckstein, “K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation,” IEEE Trans. Signal Process., vol. 54, no. 11, pp. 4311–4322, Nov. 2006. [13] M. Elad and M. Aharon, “Image denoising via sparse and redundant representations over learned dictionaries,” IEEE Trans. Image Process., vol. 15, no. 12, pp. 3736–3745, Dec. 2006. [14] J. Mairal, M. Elad, and G. Sapiro, “Sparse representation for color image restoration,” IEEE Trans. Image Process., vol. 17, no. 1, pp. 53–69, Jan. 2008. [15] J.-L. Starck, M. Elad, and D. L. Donoho, “Image decomposition via the combination of sparse representations and a variational approach,” IEEE Trans. Image Process., vol. 14, no. 10, pp. 1570–1582, Oct. 2005. [16] M. Elad, J.-L. Starck, P. Querre, and D. L. Donoho, “Simultaneous cartoon and texture image inpainting using morphological component analysis (MCA),” Appl. Comput. Harmon. Anal., vol. 19, no. 3, pp. 340–358, Nov. 2005. [17] Z. Xu and J. Sun, “Image inpainting by patch propagation using patch sparsity,” IEEE Trans. Image Process., vol. 19, no. 5, pp. 1153–1165, May 2010. [18] G. Yu, G. Sapiro, and S. Mallat, “Solving inverse problems with piecewise linear estimators: From Gaussian mixture models to structured sparsity,” IEEE Trans. Image Process., vol. 21, no. 5, pp. 2481–2499, May 2012. [19] C. Barnes, E. Shechtman, A. Finkelstein, and D. B. Goldman, “PatchMatch: A randomized correspondence algorithm for structural image editing,” ACM Trans. Graph., vol. 28, no. 3, p. 24, Aug. 2009. [20] I. Ram, M. Elad, and I. Cohen, “Image processing using smooth ordering of its patches,” IEEE Trans. Image Process., vol. 22, no. 7, pp. 2764–2774, Jul. 2013. [21] A. Criminisi, P. Pérez, and K. Toyama, “Region filling and object removal by exemplar-based image inpainting,” IEEE Trans. Image Process., vol. 13, no. 9, pp. 1200–1212, Sep. 2004. [22] A. Wong and J. Orchard, “A nonlocal-means approach to exemplarbased inpainting,” in Proc. 15th IEEE Int. Conf. Image Process. (ICIP), Oct. 2008, pp. 2600–2603. [23] G. R. Cross and A. K. Jain, “Markov random field texture models,” IEEE Trans. Pattern Anal. Mach. Intell., vol. PAMI-5, no. 1, pp. 25–39, Jan. 1983. [24] R. Chellappa and S. Chatterjee, “Classification of textures using Gaussian Markov random fields,” IEEE Trans. Acoust., Speech, Signal Process., vol. 33, no. 4, pp. 959–963, Aug. 1985. [25] F. S. Cohen, Z. Fan, and M. A. Patel, “Classification of rotated and scaled textured images using Gaussian Markov random field models,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 13, no. 2, pp. 192–202, Feb. 1991. [26] B. S. Manjunath and R. Chellappa, “Unsupervised texture segmentation using Markov random field models,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 13, no. 5, pp. 478–482, May 1991. [27] M. Li and T. Q. Nguyen, “Markov random field model-based edgedirected image interpolation,” IEEE Trans. Image Process., vol. 17, no. 7, pp. 1121–1128, Jul. 2008. [28] T. F. Chan and J. Shen, “Nontexture inpainting by curvature-driven diffusions,” J. Vis. Commun. Image Represent., vol. 12, no. 4, pp. 436–449, Dec. 2001. [29] F. Catté, P.-L. Lions, J.-M. Morel, and T. Coll, “Image selective smoothing and edge detection by nonlinear diffusion,” SIAM J. Numer. Anal., vol. 29, no. 1, pp. 182–193, 1992. [30] L. Alvarez, P.-L. Lions, and J.-M. Morel, “Image selective smoothing and edge detection by nonlinear diffusion. II,” SIAM J. Numer. Anal., vol. 29, no. 3, pp. 845–866, 1992. [31] M. J. Black, G. Sapiro, D. H. Marimont, and D. Heeger, “Robust anisotropic diffusion,” IEEE Trans. Image Process., vol. 7, no. 3, pp. 421–432, Mar. 1998. [32] M. Vetterli, P. Marziliano, and T. Blu, “Sampling signals with finite rate of innovation,” IEEE Trans. Signal Process., vol. 50, no. 6, pp. 1417–1428, Jun. 2002. [33] P. Shukla and P. L. Dragotti, “Sampling schemes for multidimensional signals with finite rate of innovation,” IEEE Trans. Signal Process., vol. 55, no. 7, pp. 3670–3686, Jul. 2007.

JIN AND YE: ANNIHILATING FILTER-BASED LOW-RANK HANKEL MATRIX APPROACH

[34] M. Signoretto, V. Cevher, and J. A. Suykens, “An SVD-free approach to a class of structured low rank matrix optimization problems with application to system identification,” ESAT-SISTA, KU Leuven, Leuven, Belgium, Tech. Rep. 13-44, 2013. [35] Z. Wen, W. Yin, and Y. Zhang, “Solving a low-rank factorization model for matrix completion by a nonlinear successive over-relaxation algorithm,” Math. Program. Comput., vol. 4, no. 4, pp. 333–361, Dec. 2012. [36] M. Sznaier and O. Camps, “A Hankel operator approach to texture modelling and inpainting,” in Proc. 4th Int. Workshop Texture Anal. Synth., Oct. 2005, pp. 125–130. [37] E. J. Candès and B. Recht, “Exact matrix completion via convex optimization,” Found. Comput. Math., vol. 9, no. 6, pp. 717–772, Dec. 2009. [38] B. Recht, M. Fazel, and P. A. Parrilo, “Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization,” SIAM Rev., vol. 52, no. 3, pp. 471–501, 2010. [39] J. Wright, A. Ganesh, S. Rao, Y. Peng, and Y. Ma, “Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization,” in Proc. Neural Inf. Process. Syst. (NIPS), 2009, pp. 2080–2088. [40] J.-F. Cai, E. J. Candès, and Z. Shen, “A singular value thresholding algorithm for matrix completion,” SIAM J. Optim., vol. 20, no. 4, pp. 1956–1982, 2010. [41] R. H. Keshavan, A. Montanari, and S. Oh, “Matrix completion from a few entries,” IEEE Trans. Inf. Theory, vol. 56, no. 6, pp. 2980–2998, Jun. 2010. [42] K. Lee and Y. Bresler, “ADMiRA: Atomic decomposition for minimum rank approximation,” IEEE Trans. Inf. Theory, vol. 56, no. 9, pp. 4402–4416, Sep. 2010. [43] Z. Liu and L. Vandenberghe, “Interior-point method for nuclear norm approximation with application to system identification,” SIAM J. Matrix Anal. Appl., vol. 31, no. 3, pp. 1235–1256, 2009. [44] M. Fazel, T. K. Pong, D. Sun, and P. Tseng, “Hankel matrix rank minimization with applications to system identification and realization,” SIAM J. Matrix Anal. Appl., vol. 34, no. 3, pp. 946–977, 2013. [45] N. Srebro, “Learning with matrix factorizations,” Ph.D. dissertation, Dept. Elect. Eng. Comput. Sci., Massachusetts Inst. Technol., Cambridge, MA, USA, 2004. [46] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Found. Trends Mach. Learn., vol. 3, no. 1, pp. 1–122, Jan. 2011. [47] L. K. Saul and S. T. Roweis, “Think globally, fit locally: Unsupervised learning of low dimensional manifolds,” J. Mach. Learn. Res., vol. 4, pp. 119–155, Jan. 2003. [48] G. Peyré, “Manifold models for signals and images,” Comput. Vis. Image Understand., vol. 113, no. 2, pp. 249–260, Feb. 2009. [49] P. Brodatz, Textures: A Photographic Album for Artists and Designers. New York, NY, USA: Dover, 1966. [50] Z. Zhang, A. Ganesh, X. Liang, and Y. Ma, “TILT: Transform invariant low-rank textures,” Int. J. Comput. Vis., vol. 99, no. 1, pp. 1–24, 2012.

3511

[51] M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo, “An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems,” IEEE Trans. Image Process., vol. 20, no. 3, pp. 681–695, Mar. 2011. [52] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: From error visibility to structural similarity,” IEEE Trans. Image Process., vol. 13, no. 4, pp. 600–612, Apr. 2004. [53] P. Arbelaez, M. Maire, C. Fowlkes, and J. Malik, “Contour detection and hierarchical image segmentation,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 33, no. 5, pp. 898–916, May 2011. [54] J. Jia and C.-K. Tang, “Image repairing: Robust image synthesis by adaptive ND tensor voting,” in Proc. Conf. Comput. Vis. Pattern Recognit., vol. 1. Jun. 2003, pp. I-643–I-650.

Kyong Hwan Jin received the B.S., integrated M.S., and Ph.D. degrees from the Department of Bio and Brain Engineering, Korea Advanced Institute of Science and Technology (KAIST), Daejeon, Korea, in 2008 and 2015, respectively. He is currently a Post-Doctoral Scholar at KAIST. His research interests include low rank matrix completion, sparsity promoted signal recovery, sampling theory, biomedical imaging, and image processing in various applications.

Jong Chul Ye received the B.Sc. and M.Sc. degrees from Seoul National University, Korea, and the Ph.D. degree from Purdue University, West Lafayette. He joined the Korea Advanced Institute of Science and Technology (KAIST), Daejeon, Korea, in 2004. He worked at Philips Research, and GE Global Research, both in New York. He is currently a Full Professor with the Department of Bio/Brain Engineering, KAIST. His current research interests include compressed sensing and statistical signal processing for various imaging modalities, such as MRI, NIRS, CT, PET, and optics. He was the winner of the 2009 Recon Challenge at the ISMRM Workshop with his k-t FOCUSS dynamic MR imaging algorithm, and a coauthor of Student’s Best Paper Award at the 2013 IEEE Symposium on Biomedical Imaging. He is an Associate Editor of the IEEE T RANSACTIONS ON I MAGE P ROCESSING, and the IEEE T RANSACTIONS ON C OMPUTATIONAL I MAGING, and an Editorial Board Member of M AGNETIC R ESONANCE IN M EDICINE.