A TWO-STEP MIXED INPAINTING METHOD WITH CURVATURE-BASED ANISOTROPY AND SPATIAL ADAPTIVITY

. The image inpainting problem consists of restoring an image from a (possibly noisy) observation, in which data from one or more regions are missing. Several inpainting models to perform this task have been developed, and although some of them perform reasonably well in certain types of images, quite a few issues are yet to be sorted out. For instance, if the image is expected to be smooth, the inpainting can be made with very good results by means of a Bayesian approach and a maximum a posteriori computation [2]. For non-smooth images, however, such an approach is far from being satisfactory. Even though the introduction of anisotropy by prior smooth gradient inpainting to the latter methodology is known to produce satisfactory results for slim missing regions [2], the quality of the restoration decays as the occluded regions widen. On the other hand, Total Variation (TV) inpainting models based on high order PDE diﬀusion equations can be used whenever edge restoration is a priority. [4], has allowed inpainted images with well deﬁned edges and enhanced object connectivity. The CDD approach, nonetheless, is not quite suitable wherever the image is smooth, as it tends to produce piecewise constant restorations. In this work we present a two-step inpainting process. The ﬁrst step consists of using a CDD inpainting to build a pilot image from which to infer a-priori structural information on the image gradient. The second step is inpainting the image by minimizing a mixed spatially variant anisotropic functional, whose weight and penalization directions are based upon the aforementioned pilot image. Results are presented along with comparison measures in order to illustrate the performance of this inpainting method.


(Communicated by Weihong Guo)
Abstract. The image inpainting problem consists of restoring an image from a (possibly noisy) observation, in which data from one or more regions are missing. Several inpainting models to perform this task have been developed, and although some of them perform reasonably well in certain types of images, quite a few issues are yet to be sorted out. For instance, if the image is expected to be smooth, the inpainting can be made with very good results by means of a Bayesian approach and a maximum a posteriori computation [2]. For nonsmooth images, however, such an approach is far from being satisfactory. Even though the introduction of anisotropy by prior smooth gradient inpainting to the latter methodology is known to produce satisfactory results for slim missing regions [2], the quality of the restoration decays as the occluded regions widen. On the other hand, Total Variation (TV) inpainting models based on high order PDE diffusion equations can be used whenever edge restoration is a priority. More recently, the introduction of spatially variant conductivity coefficients on these models, such as in the case of Curvature-Driven Diffusion (CDD) [4], has allowed inpainted images with well defined edges and enhanced object connectivity. The CDD approach, nonetheless, is not quite suitable wherever the image is smooth, as it tends to produce piecewise constant restorations.
In this work we present a two-step inpainting process. The first step consists of using a CDD inpainting to build a pilot image from which to infer a-priori structural information on the image gradient. The second step is inpainting the image by minimizing a mixed spatially variant anisotropic functional, whose weight and penalization directions are based upon the aforementioned pilot image. Results are presented along with comparison measures in order to illustrate the performance of this inpainting method.
1. Introduction. Let us consider an image from which data from one or more regions are missing. The inpainting problem can be described as that of using the available information from the occluded image for filling up the missing areas in a way that looks natural and intuitive to the human eye. Although it is clear then that defining a "good" inpainting method involves a high degree of subjectivity, some properties are to be expected. For instance, if an occluded region runs over a smooth area, we would expect the inpainting to be smooth, while if it goes through an edge we would expect the two parts set apart by the occlusion to result "appropriately" connected. These desired properties are often hard to include in a mathematical model, even more so, simultaneously. In this article we shall combine the Curvature Driven Diffusion (CDD) model proposed in [4] with the ideas of including structural prior information [2] into regularization methods to develop a two-step inpainting process which adequately complies with the aforementioned criteria.
We shall begin by fixing some notation. Our (grayscale) image will be associated to a function u : Ω ⊂ R 2 → [0, 1], where u(x, y) represents the light intensity of the point (x, y), corresponding u = 1 to white, and u = 0 to black. With D we shall denote the occlusion and with v . = u| Ω\D the known part of u. Before we proceed to build our two-step method, it is timely to describe three traditional inpainting methods that are closely bonded to it.
1.1. Tikhonov-Phillips inpainting. If the image to be inpainted is expected to be smooth, then there are traditional methods that turn out to be quite satisfactory. With the Tikhonov-Phillips regularization method of order 1 (T1), the inpainted image is defined as where T : L 2 (Ω) → L 2 (Ω \ D) is the D-occlusion operator, defined by T u = u | Ω\D , and λ > 0 is a regularization parameter. It is well known, however, that if the original image has edges or borders set apart by an occlusion, then such an inpainting method fails to extend them inside the occluded region (see Figure 1(b)). Edge preservation, however, can be fairly enhanced by appropriately modifying the penalizing term in (1). One way of doing so is to introduce an anisotropyinducing matrix field A [2], which attenuates penalization on the maximum gradient direction. In this case, the inpainted image is defined as The matrix field A makes the penalizer take into account not only the norm of the gradient, but also its direction. The properties of A and a way to build it will be presented later, but it is timely to point out here that its construction requires apriori information of the expected image. Figure 1(c) depicts an order 1 Tikhonov-Phillips anisotropic inpainting for which A was build upon an a-priori order 2 Tikhonov-Phillips inpainting of the gradient field, as suggested in [2]. Although we can see that this method works significantly better than the isotropic one for "thin" occluded regions, it becomes inappropriate as those regions widen (see Figure 1(c)). We then describe an inpainting method known to possess better edge preserving properties. 1.2. Total variation inpainting. One of the most commonly used edge preserving methods is the so called total variation (TV) regularization [1], [10]. With this method, the inpainted image is given by By using the corresponding Euler-Lagrange equation, the minimizer in (3) can be thought of as the steady-state solution of the following diffusion PDE (see [7]) TV inpainting tends to produce piecewise constant solutions with sharp edges inside the occlusion (this can be seen in Figure 2). Even though this inpainting process seems to fairly comply with the edge preservation requirement, it presents a peculiar property. In fact, if we look at the black vertical bar in Figure 2(a), we notice that it is wider than the occluded region, but whenever the occluded region is wider than the bar (Figure 3(a)), the TV solution is very different (see Figure  3(b)). Although there is no formal reason to sustain that the inpainted image in Figure 3(b) is not a "good" solution, it is a fact that most humans will prefer the image in Figure 2(b) over the one in Figure 3   Once again, this puts into evidence that any inpainting process can entail a high level of subjectivity. That being said, a "good" inpainting method could be conceived as one that would most frequently emulate what most humans would do. The rest of the article is strongly aligned with this belief.
1.3. Curvature driven diffusion inpainting. An isophote of an image represented by u can be roughly thought of as any level line of u that separates regions of different light intensities. It is highly desired for an inpainting method to be able to connect isophotes inside the occluded areas, since this will result in that the edges "broken" by the occlusion will be reconnected. In fact, in Figure 3(b) we observe that the original isophotes (the edges of the black bar) do not result connected as we would hope. Note also that the isophotes of the inpainted image have corners, meaning that their curvature κ at those points is ±∞, in contrast with the zero curvature in the isophotes of the "expected" inpainted image ( Figure 2(b)). This observation led Chan and Shen [4] to include the curvature into the diffusion model (4). By lettingD = 1 |∇u| , equation (4) restricted to the occlusion reads The ideas in [4] consist essentially in appropriately modifying the diffusion coefficient D in order to take curvature into account. This can be done by redefiningD = g(|κ|) |∇u| , where g : R + 0 → R + 0 is an increasing function such that g(0) = 0 and g(∞) = ∞. In this way, diffusion is strong where curvature is large, while it is weak where curvature is small. This modification leads to the so called Curvature-Driven Diffusion (CDD) equation Now, since the curvature of the level lines is ±∞ at the corners of the TV inpainted image in Figure 3(b), clearly such an image cannot be the steady-state of equation (5). It turns out that in this case the method strongly favors connectivity on the steady state, as it can be seen in Figure 4. Nonetheless, it occurs that this method tends to produce piecewise constant restorations, and therefore it results unsuitable for inpainting over smooth regions (see Figure 4). Before we go any further, we should mention that the presence of noise is common in this kind of problems. The CDD inpainting model (5) can take noise into account, for instance, by using TV regularization outside the occluded region. Thus, the resulting diffusion PDE takes the form Note that (6) reduces to equation (5) inside the occlusion, and to (4) outside of it.
The main objective of this article is to develop an inpainting model that complies with the three aforementioned properties. Namely, it must preserve edges while enhancing object connectivity, and allowing for smooth inpainting wherever it is deemed appropriate.
2. Inpainting via combined CDD and T1-TV methods. Regularization methods can be combined by using spatially-varying weighted averages of two or more penalizers [8] to encompass the characteristic properties of each one. We shall consider the case of mixed T1 and TV regularization, taking advantage of the fact that the first one tends to produce smooth restorations, while the latter is better suited for edge preservation. The inpainted image will then be defined as the minimizer of is an anisotropy matrix field and θ = θ(x, y) ∈ [0, 1] is a spatially-varying function weighting both penalizers at each point. Note that θ = 0 corresponds to pure anisotropic T1 regularization, while θ = 1 leads to pure anisotropic TV regularization.
The matrix field A : Ω → R 2×2 is used for introducing anisotropy to the regularization model. One way of constructing A can be found in [2], where it is built upon the gradient field of an a-priori estimation u p (x, y) of u(x, y). The first step for doing this is to define a continuous and decreasing function h : R + 0 → (0, 1], satisfying h(0) = 1 and lim t→∞ h(t) = 0, that will serve to determine the eigenvalues of A(x, y), constructed as follows: As a consequence, A has the following properties: With the required properties on h, the eigenvalues σ 1 (x, y) turn out to be small at points (x, y) where the norm of the estimated gradient ∇u p (x, y) is large, while σ 2 (x, y) remains constant. This is desirable since it will later translate into a decay on penalization in the expected gradient direction while keeping it unchanged in its normal direction. The function h can be chosen, for instance, as h(t) = 1/(1 + (t/τ ) k ), where τ, k > 0 are control parameters that could be roughly thought of as a lower threshold for the values of |∇u p | starting from which we infer the image has an edge and the width of the transition region, respectively. In [2], the gradient field estimation ∇u p was built by inpainting the components of the gradient using a second order Tikhonov-Phillips method. However, since this method produces smooth restorations, edge preservation remains compromised for wide occluded regions, as it can be seen in Figure 1(c). To overcome this disadvantage, we shall build the matrix field A based upon the gradients of an a-priori CDD inpainting of the image, which as previously noted, favors object connectivity.
In the same fashion, the weighting function θ, which also requires of a-priori information about the local regularity of u, can be built from a "pilot image" u p , obtained by means of a CDD-inpainting process. We take θ(x, y) = w(|∇u p (x, y)|), where w is an appropriately chosen increasing function, with 0 ≤ w(t) ≤ 1 ∀t ∈ R + 0 . Due to the way in which the pilot image is built, if used direclty as obtained by the CDD-inpainting process, the staircasing effect on u p will have a negative impact on the effect of the weighting and anisotropy functions, which in turn might lead to sharp artificial edges appearing in the final inpainted image. To overcome this, the pilot image is smoothed out with a small-variance Gaussian kernel.
Our full inpainting process can then be stated as follows: Step 1: CDD inpainting.: Perform a CDD inpainting on the occluded image to obtain a first estimation. That is, computeũ p as the steady state of the equation Step 2: Pilot image smoothing.: Smooth outũ p by means of a low-pass filter by computing u p . = G * ũ p , where G is a low-variance Gaussian kernel, and " * " denotes convolution. This step will result in the good definition of the gradient field ∇u p , while attenuating the staircasing effect that may appear in the final restoration.
Step 3: Construction of the anisotropy matrix field.: Use u p to construct A as Step 4: Construction of the weighting function.: Use u p to build θ as θ(x, y) = |∇u p (x, y)| max (w,z)∈Ω |∇u p (w, z)| .
Step 5: Final inpainting.: Use the previously computed A and θ to build the mixed weighted T1-TV anisotropic functional and compute the restored image as (9)û = arg min Having stated our new inpainting method, we shall proceed to briefly describe an appropriate numerical implementation.
3. Numerical implementation. We start the implementation of our inpainting method by performing a discretization over the image domain. We assume that the grayscale image domain is Ω = [0, 1] × [0, 1] and we discretize it to obtain an M -by-M pixel grid and an M -by-M matrix U , consisting its entries of the values of the function u at the centerpoints of the pixels. Next, we stack the columns of the matrix U to get a vector u ∈ R M 2 so that u M (l−1)+m = U m,l ∀l, m = 1, 2, . . . , M . For better understanding we will often identify u M (l−1)+m with u(x, y).
= g(|κ|) |∇u| ∇u, whose divergence we need to approximate. We do it by computing where h = 1/M is the pixel-width. Now, in order to compute (10) at the midpoints between adjacent pixels, we need to estimate both the gradient of u and the curvature κ at those points. Firstly, for the points of the form ∇u(x + h/2, y), we compute the components of ∇u as .
At the points of the form (x, y + h/2) the construction is analogous. As for the isophote curvature κ, we compute it explicitly as a function of the gradient of u: At points of the form (x + h/2, y), we approximate it with The computation of κ at points of the form (x, y + h/2) is analogous. To avoid division by zero in the last equation (when ∇u = 0), we replace |s| by √ s 2 + 2 , with > 0 chosen sufficiently small.
We are then ready to state the algorithm for the CDD inpainting process. It is timely to point out here that due to numerical stability issues, instead of using Euler's method to solve the IVP (as stated in [4]), we use an Adams-Moulton Adams-Bashford predictor-corrector method, which showed better performance. Our modified iterative algorithm reads as follows: Step 1: Initializing.: Define an arbitrary initial estimation u (0) , coinciding with the data outside the occlusion, and let n = 0.
Step 2: Updating.: For m : 1 . . . M , define and compute u (n+1) m as follows: where for each m, κ (n) and ∇u (n) are computed from u (n) as previously described. A second order Runge-Kutta method is used for the first two steps.
Step 3: Stopping.: If an appropriate stopping criterion (defined upon the decay of the sum of the curvature at each pixel) is reached, the process stops and the inpainted image is defined as u (n+1) . Otherwise, n is increased (n = n + 1) and the algorithm continues from Step 2. Next, we show how the mixed weighted anisotropic T1-TV regularization can be numerically implemented.

3.2.
Mixed anisotropic T1-TV implementation. To find the minimizer of the T1-TV inpainting functional given in (7), we consider the discretized version where T ∈ R M 2 ×M 2 is the diagonal matrix associated to the occlusion operator T , θ m and A m are the evaluations of the weighting function and of the anisotropy matrix field at the centerpoint of the m th pixel, and M denotes the set of indices corresponding to the interior pixels, on which the gradient is estimated, that is Although there are several efficient ways to estimate the minimizer of (13), we have chosen a half-quadratic approach, which is described below, for it has proven to cope very well with the great dimensions that this problem can present (details can be found in [5]).
First, we approximate J by a differentiable functional and make use of a duality relation to iteratively approach its minimizer. We begin by replacing the last term in (13) by a differentiable approximation in order to find the minimizer by considering the first order necessary condition. We do so by replacing, for w ∈ R 2 , the value of w 2 by φ( w 2 ), where φ : R → R is given by φ(t) . = t 2 + η 2 −η, for η sufficiently small. With this choice of φ, it can be shown ( [9]) that there exists a function ψ satisfying the following duality relation and therefore Let L x and L y be the M 2 -by-M 2 first order finite difference approximating matrices for the components of the gradient, and let R 1 and R 2 be the M 2 -by-M 2 matrices defined as R 1 . = A 1,1 L x + A 1,2 L y and R 2 . = A 2,1 L x + A 2,2 L y . Finally, let I be the M 2 -by-M 2 identity matrix, and define the functional It can be shown ( [5]) that where J φ is the functional obtained from J by using approximation (15). Hence, our problem turns out to be equivalent to minimizing K with respect to both u and s simultaneously. Note that the first order necessary condition on K with respect to u can be written as (20) In order to minimize K θ,φ (u, s) with respect to s, we define where t m,1 and t m,2 are defined as in (16) and (17), respectively, and resort to (14) to deduce that b m must satisfy Although details on de derivation of equation (21) can be found in [6], we give here a brief sketch of the proof. Let f : R 2 → R be defined as f (s, t) .
Step 3: Updating b.: Update b j using equation (21): Step 4: Updating u.: Update u j by solving the linear system where B j is the M 2 -by-M 2 diagonal matrix with elements b j m for m ∈ M and 0 otherwise. It is worth noticing that this linear system is well posed ( [8]) since the matrices appearing as a consequence of the penalization terms have strictly positive eigenvalues, and hence solving it entails no difficulties.
Step 5: Convergence.: if a previously defined convergence criterion is satisfied, the algorithm ends and our inpainted imageû is defined as u j (in our case, this criterion was defined upon the norm u j − u j−1 ). Otherwise, the algorithm repeats from step 2. In the next section we present some examples of the performance of the previously described inpainting method. 4. Inpainting application results. We begin by comparing the performance of our new approach with the isotropic order-one Tikhonov-Phillips and Curvature-Driven Diffusion methods on the previously used test image, occluded over both smooth and piecewise constant regions. A 1% Gaussian white noise was added to the grayscale image. Figure 5 depicts the occluded noisy image together with the results obtained with three different inpainting methods. As it can be clearly seen, the mixed anisotropic T1-TV inpainting outperforms the isotropic T1 method in terms of edge preservation, while it also works better than CDD inpainting on the smooth region of the image. Although this is a somewhat subjective analysis, these conclusions are supported by the corresponding peak signal-to-noise ratios (PSNR) shown in Table 1. The PSNR is defined as follows: PSNR(û) = 10 log 10 1 MSE(û, u 0 ) , whereû, u 0 ∈ R M 2 are the vectors associated to the inpainted image and to the original image (unknown in real problems), respectively, and Isotropic T1 CDD Anisotropic T1-TV PSNR 20.140 35.497 36.329 Table 1. PSNR values for the test image ( Figure 5).
Up next, we show the performance of our new method on a more realistic example. A 300 × 300 pixel grayscale image was occluded and contaminated with 2% Gaussian white noise. Figure 6 shows the resulting occluded noisy image, along with the isotropic T1, CDD and anisotropic T1-TV inpainted images. Once again we observe an improvement on the quality of the inpainting with the two-step method, which is reflected on the PSNR values shown on Table 2.  Table 2. PSNR values for grayscale image ( Figure 6).
Next we show an example of the performance of our method on a color image. In this case, the inpainting process was performed separately over the red, green and blue layers of a 300 × 300 pixel image. Here again, a 2% Gaussian white noise was added to the occluded image. Figure 7 shows the occluded noisy image, along with the CDD, isotropic T1 and the anisotropic T1-TV inpainted images. A close look at Figure 7(c) shows that CDD inpainting tends to produce artificial edges, which are not produced by the combined T1-TV method. Note also that although no artificial borders are generated with the latter method, significantly better edge preservation properties are achieved with respect to the pure T1 method. Table 3 shows the obtained PSNR values for the three inpainting methods. Now we compare the performance of our two-step method with that of the CKS inpainting model ( [3]), which is based upon Euler's elastica and curvature. This  Table 3. PSNR values for color image (Figure 7). model allows for smooth inpainting, and so it provides a way to avoid the artifacts which appear using CDD. Nonetheless, this smoothing might entail loosing sharpness of inpainted edges and borders. To illustrate this, inpainting results for the same images used before are displayed on Figure 8, obtained with a 200 × 200 image, previously contaminated with a 2% Gaussian additive noise. It can be observed that our two-step method has better performance, as the spatial adaptivity provides a better compromise between edge preservation and smooth inpainting. Table 4 shows the PSNR values for the displayed restorations. Finally, we present an example of object removal via masking-inpainting. If one wishes to remove an object from a certain image in an "unnoticeable" way, an artificial occlusion can be made over the object, and then proceed with the inpainting. Figure 9 shows an example of the results of such a procedure over a  Table 4. PSNR values for color image (Figure 7). 300 × 300 color image. The compromise between edge preservation and smooth inpainting of the mixed anisotropic T1-TV model can be once again appreciated here. Note that computing the PSNR values does not make sense in this case. Computational cost: Although certainly there is room for improvement in the implemented algorithms, after several trials we can assure that the computational cost for the mixed inpainting approach does not increase severely with respect to the pure CDD method. Our programs were run on a PC with an Intel Core i7-3930K processor, CPU@3.20 GHz x 12. For instance, for the 300 × 300 grayscale image displayed on Figure 6, the CPU-time for the CDD process was 2482.4 seconds, while the CPU-time for the T1-TV model (having previously computed the pilot image u p and the parameters) was 56.8 seconds. The CPU-time for the full algorithm, including building the pilot image, parameter estimations and auxiliar computations was 3572.4 seconds. We mention however that optimizing the algorithm was not a main goal of our work. Experimentation suggests, for instance, that a variable time-step on the diffusion process and refining parameter estimation on the T1-TV inpainting might lead to a significantly faster performance.

5.
Conclusions. In this article we developed a two-step method for image inpainting. The first step consists of implementing a Curvature-Driven Diffusion process, which serves to obtain a good approximation of the gradient of the image inside the occlusion. Such an approximation is then used to construct an appropriate weighting function and an anisotropy-inducing matrix field. The second step consists of using these functions to build a mixed Tikhonov-Total Variation spatially varying anisotropic functional, whose global minimizer defines the final inpainting. The use of a CDD inpainting process in the construction of the weighting function and the anisotropy matrix field results in the fact that the well known objectconnectivity property of this method is incorporated into the mixed T1-TV model. This, along with the spatial adaptivity and anisotropy features of the mixed T1-TV method results in an inpainting model having both good object-connectivity and edge preservation properties as well as high-quality inpainting performance over smooth regions.
The performance of the method was shown through several examples, in which it produced better results than the CDD, T1 and CKS inpainting models. Furthermore, a thorough analysis of the color examples shows that while the CDD inpainting seems to sometimes produce little artifacts near edges inside the occlusions, such artifacts do not appear when using our two-step method.
Finally, it is appropriate to mention that there is room for further improvements. For instance, other means of constructing the weighting function θ and the matrix field A can be explored. Also, if there is adequate available a-priori information about the "expected" image, that information could be embebed into the model through the functions θ and A.