On the image inpainting problem from the viewpoint of a nonlocal Cahn-Hilliard type equation.

Motivated by the fact that the fractional Laplacean generates a wider choice of the interpolation curves than the Laplacean or bi-Laplacean, we propose a new non-local partial differential equation inspired by the Cahn-Hilliard model for recovering damaged parts of an image. We also note that our model is linear and that the computational costs are lower than those for the standard Cahn-Hilliard equation, while the inpainting results remain of high quality. We develop a numerical scheme for solving the resulting equations and provide an example of inpainting showing the potential of our method.

Introduction

Digital image inpainting is the problem of modifying parts of an image such that the resulting changes are not trivially detectable by an ordinary observer. It is used to recover the missing or damaged regions of an image based on the data from the known regions. It represents an ill-posed problem because the missing or damaged regions can never be recovered correctly with absolute certainty unless the initial image is completely known.

In this paper we are concerned with the following problem. Let be a square image domain and an open region with smooth boundary such that . Let f be the original image, known only on , and let . For a constant , we aim to solve the following problemwhere is the interpolation of the original image f. In the special case when and (so called double-well potential), Eq. (1) is the famous stationary Cahn-Hilliard equation (CHE) (see the original paper [10]), a well-known macroscopic field model for the phase separation of a binary alloy at a fixed temperature. It is derived from the Helmholtz free energywhere u typically denotes the concentration, is the range of intermolecular forces, is the free energy density and the last term is a contribution to the free energy originating from the spatial fluctuations of u.

Note that it is almost a rule that nonlinear PDEs (like Perona-Malik or the CHE mentioned above; see also [35]) often capture the most interesting phenomena. This increases the computational costs and makes the numerical procedure more complicated.

In this contribution we assume H that to be quadratic, which yields a linear equation, but instead of integer order derivatives we deal with a fractional order equation. The motivation for this comes from a simple observation from fractional calculus. Namely, recall that for given boundary values linear diffusion yields only linear solutions. On the other hand, has a much wider set of solutions and due to this, it is reasonable to expect that the image inpainting using fractional equations produces images that seem more natural (see Fig. 1).

Fig. 1

Example of 1D inpainting problem on the . Biharmonic equation (green), integer order CHTE (magenta) and fractional order CHTE (red).

The aim of this paper is to study the application of the fractional generalization of the Cahn-Hilliard type equations (CHTE) given in (1) to the image inpainting problem and to propose a fast algorithm for obtaining its numerical solutions. Through several examples, we are going to show that fractional PDEs produce superior results over integer order PDEs.

To this end, we derive a fast algorithm based on the matrix decomposition that solves (1) (formulated as (27)) in the local as well as in the non-local case. In both cases, the idea is to use appropriate arrangements of the discrete equations obtained by the finite difference method (see [12]) so that the computed matrix of the linear system exhibits a sparse structure with block symmetry (see (35)). This structure enables us to derive the recursive relations for the computation of the decomposition that, by using simple backward and forward substitutions, yields the solution. We also carry out a comparison of this approach with the standard algorithms for numerical solutions of the sparse linear system.

We would like to emphasize that the discretized fractional order partial differential equation (PDE) under the consideration serves as a motivation for the construction of a fast and efficient inpainting algorithms rather than as a problem from a purely mathematical point of view that will be submitted to the rigorous numerical analysis.

The rest of the paper is structured as follows. In the next section, we give a short overview of the previous approaches, motivations and ideas underlying the inpainting problem. In Section “Numerical method”, we introduce the notion of the discrete Laplacean and its fractional powers with the applications to the equation under consideration. Together with that, we derive an algorithm based on matrix decompositions for both local and nonlocal case for the purposes of fast image inpainting. In Section “Results” we present the application of the introduced ideas on several testing images, comparing it with well known linear methods. Finally in Section “Conclusions and further work”, we finish with a short discussion and ideas for the future work.

A short overview of previous results

The literature regarding the PDEs with the applications to the image inpainting problems is extensive. The terminology of digital image inpainting first appeared in the paper of Bertalmio in [4], based on the discretization of the transport-like PDE modelwhich is, for stabilization purposes, coupled with the anisotropic diffusionwhere is a smooth cut-off function that forces the equation to act only . Furthermore, represents the perpendicular gradient of the image and this is the term that controls the speed of the transport. In this model is the curvature along the isophotes (curves on a surface that connect points of equal brightness), is a measure of image smoothness and is the propagation direction, i.e., the direction of smallest spatial change. The idea was to extend the image intensity in the direction of the isophotes arriving to the subset , where is the inpainting domain. It can be shown that the steady state equation of (4) is the equation satisfied by the steady state inviscid flow in the two dimensional incompressible Navier-Stokes equation [5]. In this concept we can identify image intensity as a stream function for which the Laplacian of the image intensity models the vorticity that results in an algorithm that continues the isophotes while matching gradient vectors at the boundary of the inpaiting domain.

Given the subjective nature of the image inapainting problem it is reasonable to argue that the brain interpolates broken missing edges using elastica-type curves. More precisely, if we slightly extend the inpainting domain in and denote it by , one can extrapolate the isophotes of an image u by a collection of curves with no mutual crossing, which coincide with the isophotes of u on and minimize the energywhere is the intensity span and is the curvature. For some parameters and , depending on the specific application, this energy penalizes a generalized Euler’s elastica energy. In [33] the authors proposed two novel inpainting models based on the seminal Mumford-Shah image model [24] and its high order correction, called Mumford-Shah-Euler image model. The second one improves the first model by replacing the embedded straight-line curve model with Euler’s elastica, first introduced by Mumford in the context of curve modeling. This approach is not very computationally efficient, and attempts have been made to create more effective schemes, and various extensions involving augmented Lagrangians were considered (see [34], [37], [38]).

More recently, modified CH and Allen–Cahn equations for the inpainting of binary images have been analyzed [6], [7], [9], [15], [18]. In the integer order case, besides the standard double well-potential, researchers have investigated the nonsmooth double obstacle potential [8] (considered in this contribution) with the applications to the grayscale images, as well as the logarithmic potential [13]. They successfully demonstrated that a simpler class of models (comparing to [33]) can be modified to achieve fast inpainting of simple binary shapes, text reparation, road interpolation, and image upscaling. This was the motivation for the development of advanced numerical methods based on finite difference methods [11], finite element methods with preconditioning [8], spectral methods [9], and operator splitting [20] in order to make the practical computation faster and more efficient. For a mathematical analysis of the fourth order models we refer to [35].

Several years later the investigation of fractional models started in signal and image processing, a tool already widely used in physics [30]. Fractional derivatives with respect to space have been used in the attempts to describe more accurately the anomalous diffusion or dispersion, where a particle plume spreads at a rate inconsistent with the well known Brownian model of motion, and the plume may be asymmetric. The application of fractional calculus resulted in superior algorithms for the edge detection [21], filters for texture improvement [25], noise removal [3], [39], etc. Roughly speaking, the idea is to solve the following equationon either the entire or a part of the image domain where is an appropriately chosen function, depending on the specific application. For a numerical treatment of such equations see [12], [22], and for applications in image inpainting cf. [1].

In general, fractional differential equations are characterized by nonlocal and spatially heterogeneous properties in which classical models fail to provide the adequate results. Regarding image inpainting problems it has been shown that they improve the image quality as well as the peak signal to noise ratio (PSNR) [1], [19], [41]. For a review of the field, starting from simple harmonic inpainting to the state of the art methods in PDE based inpainting see [31], [32].

Physical reasoning in the integer case

Even though ad hoc adjustments of the governing equations have been known to produce impressive results (for example the Perona-Malik equation [36]) it is important to keep in mind the underlying thermodynamic theory for the construction of this class of tools, at least when dealing with integer order equations. Letbe the free energy functional, where is the local free energy per molecule. Next, we expand in a Taylor series around and neglect higher order termsImposing that the free energy is invariant under all rotations and reflections i.e.we get the local free energyAfter integration over the domain and integration by parts we obtain the total free energywhere .

Let us consider the mixture of two miscible phases, where and are relative concentrations of the components such that and . In the general case, the corresponding flux is given bywhere is the (generalized) diffusivity, and are the chemical potentials of the components. By Fick’s first law, the gradient of two chemical potentials can be calculated as a variation of a corresponding free energy potentialCombining (14), (15) and assuming a general anisotropic situation, as is often the case in image processing, one obtainsNow if we assume that the mass is conserved we obtain a class of equations depending on the choice of energy functionalIn typical situations, while constructing PDE interpolation or a filter, one can choose either the specific diffusivity or the free energy functional in an effort to process the image. In the simplest case where we neglect terms with derivatives and take to be quadratic, one obtains the linear diffusion equation (Gaussian filter), the very first PDE model for harmonic inpainting and image processing. Note that harmonic inpainting is a linear extension scheme, and, because of this, images obtained by employing such a technique do not produce very convincing results. In practice, one replaces the Laplace operator with the biharmonic one to define cubic inpainting by taking the total free energy in the form (13) to get

As for the application in image processing, we assume that (or we can replace by a non-isotropic constant coefficient elliptic operator), and we take the stationary case, i.e. we arrive at (1) with quadratic H and .

Extension to the fractional case

Continuing in the same direction as explained in the introduction, it is natural to go beyond integer derivatives in order to increase the variety of curves which can be used during the inpainting procedure (we recall that in the harmonic case it is a linear curve and in the bi-harmonic case it is a third order polynomial; see Fig. 1 and Fig. 2).

Fig. 2

Graphical representation of the image inpainting problem of the Runge function using integer order and fractional order equations with on .

This intention can be supported by the following arguments. The system is expected to evolve so that Helmholtz energy (3) decreases in time and approaches a minimum. For a Hilbert space V, Eq. (1) with can be viewed as a gradient flowwhere D is some positive constant and the gradient of at a point is defined as followsLet denote standard scalar product, now by ignoring boundary terms for the moment we can obtain thatNow we can identify two distinct situations. In the first one, we can choose to obtain the gradientand the associated gradient flow isthat gives already mentioned Allen-Cahn equation. It is well known that this equation does not preserve mass. Alternatively, one can take with the scalar product . In this case the associated gradient flow is Eq. (1) with that does preserves mass.

Next natural step would be to explore the gradient flow in the space where with the appropriate scalar product . In this case, a gradient is given byIn conclusion, one can consider a gradient flow in the Sobolev space , for where the choice one recovers Allen-Cahn equation and yields CHE. At this point, further generalizations could be obtained by allowing the other Laplacian in (25) to be of fractional order. For more details please see [2].

Numerical method

In order to introduce a discretization procedure, we shall first rewrite (1), (2) within a more suitable way for the numerical treatment. Denote by a positive constant (in the applications below, we take ), by for the indicator function of the . Namely, we will investigate the following equivalent variant of (1), (2)

Fractional derivatives can be defined in several, essentially equivalent, ways (see e.g. the classical books [26], [29]). However, depending on the situation, certain variants of the definition of fractional derivatives provide better operational aspects.

The fractional power of the discrete Laplace operator can be found in [12, Theorem 1.1].

Let and . Furthermore, let be such thatThenwhere the discrete kernel is given by

Discretization of the integer order problem

In this section, motivated by applicability of the algorithm and methodical reasons, we want to lay out the main ideas of the numerical scheme in the integer order case, i.e. that will be extended to the non-integer case.

We discretize at grid points in the square domain which are at with and , with , where n represents the resolution of the image. Let us abbreviate and . By using the standard three-point discretization to approximate and , and applying it to , one can derive the following:Using this approximation we can write the first term on the left hand side of (27) in the following formwhere the truncation error is bounded by . For defined by (40), Eq. (40) yields a set of linear equations in unknowns as follows

Let us now rewrite the equations given by (33) as a single matrix equation by making the following arrangement of the unknowns . We will use the linear isomorphism that reshapes an matrix into a vector of elements, in the following wayThis approach leads to a symmetric block pentadiagonal matrixabbreviated by , where the matrices and C are defined as followsThus we arrive at the problem of solving the sparse symmetric linear system . If we consider the columns of S to be vectors in , we can easily conclude that they are linearly independent, so S is a regular matrix. Hence we can conclude that the linear system has a unique solution. Before turning to the question of solving it, we make a small note, relevant for the practical implementation. Namely, for implementation purposes, the matrix S can be constructed easily using the Kronecker productwhere I is identity matrix, , and , where is the standard Kronecker symbol.

Formulation of the linear system

Now we shift our focus to the discrete form of (27). Taking into account the finite difference equations from the previous section we want to solvewhere Z is defined as , where I is the identity matrix, is the original image and is the characteristic function of the set defined as in (27). Note that, in general, Z is not a symmetric matrix.

Solutions of the linear system

We aim to apply a suitable factorization to the symmetric matrix so that , where L is a lower-diagonal matrix. To this end we suppose that L (and ) is a lower (upper) triangular matrix of the following formLet us note that this is a very general form and that for practical purposes the matrix might be even more sparse, depending of the size of the inpainting domain. We also observe that the special form of the matrix Z enables us to perform this symmetrisation in only operations. Furthermore, keeping in mind that each and are matrices we obtain the following recursive scheme for determining the elements of the lower-diagonal matrix L by using only the definition of the matrix :where to to and are given by

Note that in each step one needs to compute the inverse of and this can be done in operations. Now we proceed in two steps: a) solve , b) use the computed Y to solve and obtain the solution X. After the computation of the lower-diagonal matrix L we have the following recursionwhere and have already been determined, beginning withBecause the matrices have already been computed (during the decomposition step) we can perform a forward substitutionin order to obtain Y. The solution X is finally computed through the backward substitution given byHere, the boundary cases are defined in the following wayBecause the inverses of are known, it is easy to see that both, forward and backward substitutions can be done in operations.

In total, this yields operations.

Discretization of fractional differential equations

Now, we are ready to deal with the fractional variant of (1), (2). In our simulations we have fixed because numerical experiments have indicated that this could be the compromise between the quality of the inpainting results and keeping the numerical scheme relatively simple. In this way, with the appropriate selection of the parameters a and , Eq. (1) can be viewed as a fractional inpainting with a corrective local term of high order. Using Theorem 1 (and notations from there) and assuming , we haveHere we have neglected terms containing the constant , with because the constants rapidly tend to zero for increasing m and because taking more terms does not seem to influence the subjective assessment of the inpainted image. Thus, it has a negligible influence in minimizing the relative error (see Section “Results”).

Next, we proceed similarly as in the integer case by defining.

, where the matrices and G are given byNow we define the matrix Z as in (40) and follow the same steps as in the integer case.

Results

In Fig. 2 we set and compare applications of different equations on a 1D inapainting problem on the domain . Let us note that the simple integer order inpainting using the Laplace equation resulted in a relative -error of 0.59974, the linear integer order CHTE () (27) produced a relative -error of 0.38415, and on the other hand the (linear) fractional order CHTE () produced a relative -error of 0.05859. Besides the lower relative error, fractional order equations seem to preserve the image features (also see Fig. 1) and thereby produces images that look more natural. Clearly, the fractional order CHTE delivers superior results compared to the linear integer order equations.

Furthermore, we have performed the experiment on 100 1D test images with 3 different inpainting domains and different fractional orders in (69) and compared it to the results obtained by the integer order equations. By choosing the result of the fractional order inpainting with the least -relative error (with respect to the undamaged image) and comparing it to the error produced by integer order inpainting, we conclude that the fractional order inpainting approach yields on average lower -relative error. The selection of the parameters was done by exhaustive enumeration method. Selected images from this experiment are presented in Fig. 1. Further details are given in Table 1.

Table 1

Values of the parameters for the inpainting results shown on Fig. 1.

Inpainting method	∊^2	min	max	Lrel2
Image (a)	–	0.00	0.53	–
Integer order biharmonic eq.	1	0.06	0.53	0.223
Integer order CHTE	10	0.14	0.53	0.746
Fractional order CHTE μ=0.30	1	0.00	0.53	0.014
Image (b)	–	0.00	1.00	–
Integer order biharmonic eq.	1	0.60	1.15	0.878
Integer order CHTE	10	0.63	0.96	0.757
Fractional order CHTE μ=0.80	1	−0.18	1.17	0.320
Image (c)	–	0.14	0.75	–
Integer order biharmonic eq.	1	0.14	0.66	0.118
Integer order CHTE	100	0.14	0.63	0.163
Fractional order CHTE μ=0.40	1	0.14	0.74	0.010

Moreover, for different dimensions n of the system and different numerical methods we have performed 6 independent measurements of the running times required for solving the inpainting problem. The average computational times are presented in Fig. 3, where one can compare the computational performance of the proposed approach as compared to the classical numerical methods for solving such a system. Note that the running time of the algorithm under the consideration in the case was only 1.87 s whereas other methods were not able to yield a solution within 100 s. This experiment was performed on the standard desktop computer.

Fig. 3

Comparison of different solvers for the system. PM denotes the approach proposed in this paper.

In Fig. 4, the test example consists of a gray scale image that contains a wide damaged area in the shape of a rectangle in the middle of the image. For three different parameters the proposed method was tested against the MATLAB inpainting function called inpaintn, transport equation of Bertalmío [5], Laplace and biharmonic inpainting as well as integer order CHTE and two total variation (TV and high order TV denoted by TV4) inpainting methods. The MATLAB files for Laplace, transport and total variation methods are available on the Web.1 2 All simulations are run with standard parameters that were determined by the exhaustive enumeration approach (Fig. 4, results (d)–(i)).

Fig. 4

Inpainted gray scale stripes image using different PDE based inpainting methods. (a) Original image, (b) Image with the inpainting domain, (c) Matlab inpaintn function, (d) Transport equation of Bertalmío, (e) Local Laplace inpainting, (f) Local biharmonic inpainting, (g) Integer order CHTE, (h) TV inpainting, (i) TV4 inpainting, and fractional order CHTE with (j) , k) , (l) .

For further details on values of the parameters please see Table 2.

Table 2

Results for the gray scale stripes inpainting using different models, presented on Fig. 4.

Inpainting method	a	∊^2	PSNR	SSIM	Lrel2
Image with crack	–	–	9.8941	0.78550	0.46369
Matlab inpaintn function	–	–	26.266	0.56987	0.064359
Transport eq.	–	–	18.284	0.80683	0.17605
Integer order Laplace eq.	0	1	24.257	0.80114	0.088172
Integer order biharmonic eq.	1	0	27.307	0.81967	0.061465
Integer order CHTE	1	10	25.742	0.81167	0.073854
TV	–	–	19.762	0.78139	0.14869
TV4	–	–	24.460	0.78756	0.085689
Fractional order CHTE μ=0.7	10	1	37.580	0.89977	0.019116
Fractional order CHTE μ=0.8	10	1	31.774	0.82172	0.037167
Fractional order CHTE μ=0.9	10	1	28.810	0.80908	0.052209

In the Fig. 5 we have applied the same methods as in the Fig. 4, but this time to a real and more complex image – the famous Lena with a 2.5× zoom covering details of the nose and right eye region. For further details on values of the parameters please see Table 3.

Fig. 5

Inpainted gray scale Lena using different PDE based inpainting methods with 2.5× zoom over the nose and right eye region. (a) Original image, (b) Image with the inpainting domain in blue, (c) Matlab inpaintn function, (d) Transport equation of Bertalmío, (e) Local Laplace inpainting, (f) Local biharmonic inpainting, (g) Integer order CHTE, (h) TV inpainting, (i) TV4 inpainting, and fractional order CHTE with (j) , (k) , (l) .

Table 3

Results for gray scale Lena image inpainting using different models, presented on Fig. 5.

Inpainting method	a	∊^2	PSNR	SSIM	Lrel2
Image with crack	–	–	20.612	0.92782	0.146790
Matlab inpaintn function	–	–	39.661	0.98699	0.008313
Transport eq.	–	–	35.339	0.97638	0.022761
Integer order Laplace eq.	0	1	38.639	0.98532	0.009822
Integer order biharmonic eq.	1	0	40.080	0.98733	0.009100
Integer order CHTE	1	10	38.982	0.98602	0.009420
TV	–	–	30.468	0.45713	0.019342
TV4	–	–	30.746	0.54271	0.010896
Fractional order CHTE μ=0.7	10	1	40.892	0.98958	0.007518
Fractional order CHTE μ=0.8	10	1	39.643	0.98802	0.008540
Fractional order CHTE μ=0.9	10	1	39.265	0.98700	0.009382

Moreover, we have performed thousands of experiments with above mentioned approaches (shown in Fig. 4 and Fig. 5, results (d)–(i)), however, we do not exclude the possibility that even better results for the other approaches could be obtained by a fine tuning of the parameters. Based on the tests performed for the fractional order CHTE, it seems that the best inpainting results are obtained for the values of close to , although this is probably subjected to the features of the image under the consideration.

In addition, in Fig. 6, we have applied the proposed method on the RGB image where each color channel was treated separately. We see that the difference between the original image and the inpainted one is almost not detectable.

Fig. 6

Example of the image produced by the proposed equation. (a) Original image, (b) Image with inpainting domain, (c) Result of the inpainting.

Conclusions and further work

The success of the inpainting depends on the choice of curves which can be used to interpolate damaged parts of the image. If we have only linear curves or third order polynomials as in the case of the harmonic or biharmonic inpainting approach, we cannot obtain satisfactory results. One way to overcome this limitation and still remain in the framework of analysis of harmonic and biharmonic PDEs is to add nonlinear terms (this is the case with the CHE), but such an approach decreases computational efficiency and usually requires a non-standard numerical treatment.

In the current contribution, we extended the choice of possible interpolating curves not by adding a nonlinear (correcting) terms, but by replacing integer by fractional order derivatives, staying at the same time in the linear setting. This significantly simplifies the numerical treatment of the problem and decreases computational costs. On the other hand, we find the obtained results at least equally convincing as the ones obtained using the integer order CHTE.

In future work, we shall try to extend this approach by introducing equations with nonlinear coefficients and derivatives of variable order [14] or derivatives of complex order [29]. We shall also continue in the direction of rigorously proving the convergence of the scheme and optimizing the order of the equation used for the inpainting.

Compliance with Ethics Requirements

This article does not contain any studies with human or animal subjects.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.