Understanding image inpainting with the help of the Helmholtz equation

Partial differential equations have recently been used for image compression purposes. One of the most successful frameworks solves the Laplace equation using a weighting scheme to determine the importance of individual pixels. We provide a physical interpretation of this approach in terms of the Helmholtz equation which explains its superiority. For better reconstruction quality, we subsequently formulate an optimisation task for the corresponding finite difference discretisation to maximise the influence of the physical traits of the Helmholtz equation. Our findings show that sharper contrasts and lower errors in the reconstruction are possible.


Introduction
The reconstruction of an image from a sparse subset of all pixels is known as inpainting [1]. For the application of image compression, the selection of the data has to be optimised. Let us emphasise that this is by no means a simple task. Selecting 5% of the pixels from a 256 Â 256 pixel image offers more than 10 5000 possible choices. In [2][3][4][5][6][7], corresponding strategies are devised via partial differential equations (PDEs). Related models for the inpainting step using the Allen-Cahn model have also been considered in [8], whereas the authors of [9] analysed the Cahn-Hilliard equation. Finally, a broader discussion on fluid dynamics for image reconstruction tasks has been discussed in [10]. The results from [2,6] motivated the authors of [11] to suggest an optimal control-based approach with a relaxed formulation of the Laplace equation given for known data f by and additional boundary conditions. The support of the function c indicates that the data locations used for reconstruction should be minimised while preserving a good reconstruction quality. In [11], a local contrast enhancing effect was also observed if c maps to values outside of [0, 1]. Based on [6], this finding was reinforced in [12] where an equivalence with a tuning of the data f was proven. A concrete explanation for this behaviour was not given. However, the understanding of the influence of c on the reconstruction u is crucial for the understanding and improvement of current and future approaches to optimise the inpainting data. Our Contributions. We show that (1) is related to the Helmholtz equation with a non-constant refraction index if cðxÞ [ 1 and deduce from this interpretation the benefits of non-binary-valued functions c and thus the observations in [11].
In addition, we provide details on how to maximise the benefits gained from this insight. As discussed in [13], it is important to assert that solutions of (1) exist and are unique for each feasible choice of c, and 8 / 7 is stated as its upper bound. To improve this finding, we formulate the finite difference discretisation as an optimisation task. This allows us to specify larger feasible ranges for the values of & Laurent Hoeltgen hoeltgen@b-tu.de 1 Chair for Applied Mathematics, Brandenburg University of Technology Cottbus-Senftenberg, Platz der Deutschen Einheit 1, 03046 Cottbus, Germany c than in [5,13]. We obtain more accurate reconstructions and a stronger contrast enhancing effect.

Inpainting and the Helmholtz equation
Let us briefly recall the mechanisms of inpainting with homogeneous diffusion (IHD). Let f : X ! R be a smooth function on some bounded and open domain X & R 2 with a sufficiently regular boundary oX. Moreover, let us assume that there exists a closed non-empty set of known data X K $X that we interpolate. IHD considers the following PDE with mixed boundary conditions: and where o n u denotes the derivative of u in the outer normal direction. We refer to [14] for an extensive study on the existence and uniqueness of solutions. Following [6], we introduce the confidence function c : X ! R indicating the presence of data. We set cðxÞ to 1 for x 2 X K and 0 otherwise. Lets us rewrite (2) as with Neumann boundary conditions along oX n oX K . As shown in [2,6], the choice of c has a substantial influence on the solution. Interestingly, (3) also makes sense when c is not binary-valued but takes values in R. This has been exploited in [11], where (3) is complemented by a convex energy to obtain a sparse and optimal support for c.
Let us now combine the idea of a non-binary-valued confidence function c with the mixed boundary value problem given in (2). We consider: with u ¼ f on X K ; and o n u x ð Þ ¼ 0 on oX n oX K : ð4Þ for non-binary-valued c and equivalent to (2) for binary-valued c with range 0; 1 f g.
The major difference between (4) and the previous formulations lies in the distinction between c ¼ 1 and c 6 ¼ 1. We now proceed to the first important finding of this paper. ð Þ [ 1, we can subdivide X n X K into disjoint regions, where c\1 and where c [ 1. These regions are separated by X K on which the solutions u are enforced to coincide with f. Thus, the problem decouples and allows us to discuss these two cases independently. The case 0 6 c x ð Þ\1 has already been discussed in [13] and will not be investigated further in this paper. Inside those regions, where c [ 1, we can divide (4) by 1 À c x ð Þ and obtain the following formulation: where (5) is the inhomogeneous Helmholtz equation with a refraction g and mixed boundary conditions [15]. h Theorem 1 gives us valuable insight into the properties of our inpainting equation. The Dirichlet data in X K in (5) can be interpreted in physical terms as a radiation source. The solutions u model the spread of this radiation inside X n X K and the superposition of radiated waves causes the contrast enhancement. Thus, our observations provide a physically motivated explanation for the equivalence shown in [12] and the usage of mask values outside of the range 0; 1 ½ . Furthermore, the largest refraction numbers g are obtained for values of c slightly above 1 (g ! 1 for c & 1). If c ! 1, then g ! 1, and the sharpening effect vanishes. This explains why values for c significantly larger than 1 are rarely observed in practice.
To get the best possible results for image reconstructions, it is essential to maximise the admissible range of g. At the same time, it should be asserted that the discrete setup is well posed. In [13], the author analysed (3) and provided bounds that guarantee that the discretised PDE in (3) has a unique solution. For standard finite difference schemes, these bounds are given by 0 6 c i 6 8=7 for all stencil points c i and 0\c i for at least one c i .

Improved schemes for the inpainting equation
We now follow the philosophy to optimise key features of our discrete operator for fixed grid parameters and maximise the feasible range for the refraction directly within the design process. For simplicity, we restrict ourselves to 3 Â 3 stencils.
To pursue our goals, we need the 2D Taylor expansion. For a sufficiently smooth function f : R 2 ! R, the Taylor expansion of order n around x 0 is given by x À x 0 ð Þ a , where we employ multi-index notation on a 2 N 2 and where D a f x 0 ð Þ denotes the partial derivatives of f evaluated at x 0 . Given a regular 2D grid with constant step sizes h in each direction, we use the Taylor approximations in all 8 neighbouring pixels. Each Taylor expansion uses x 0 :¼ x 0 ; y 0 ð Þ as centre. This yields the 9 combinations f x 0 þ kh; y 0 þ jh ð Þ % P a j j6n D a f x 0 ;y 0 ð Þ a! kh; jh ð Þ a with k and j 2 À1; 0; 1 f geach. In a next step, we express the desired derivative D a f x 0 ; y 0 ð Þ as a linear combination of all these positions in our stencil: Inserting corresponding Taylor expansions and performing a comparison of coefficients leads to a linear system of equations. Its solutions represent the coefficient vectors k for the stencil. Omitting, for clarity, the common argument ðx 0 ; y 0 Þ, we have up to Oðh 3 Þ. Unfortunately, this approach does not lead, in general, to a square system matrix. The number of unknowns coincides with the number of stencil positions.
On the other hand, if we perform a Taylor expansion up to order n, we obtain nðn þ 1Þ=2 equations. These numbers rarely match. Nevertheless, unless the equations contradict, it is still possible to determine a particular solution k ¼ i¼1 as well as a basis m 1 ; m 2 ; . . . f gfor the nullspace of the matrix and express all solutions as k þ P j b j m j . By computing a particular solution of the linear system for an approximation of the second-order derivatives f xx , we arrive at the parametric representations presented in Fig. 1. The stencil for f yy is obtained analogously and corresponds to transposing the stencil of o xx f . The stencil for the Laplacian as in Fig. 1 is obtained by adding the stencils for f xx and f yy : Let us now compute the stencil entries for our inpainting task from (3) resp. (4). As discussed in [11], these equations can be discretised and written as follows: where c, u, and f are the discretised variants of c, u, and f respectively. The matrix L is the discrete analogue of the Laplace operator with incorporated boundary conditions, whereas I is the identity matrix. We assume that L is discretised with the stencil from Fig. 1. The system matrix of (8) is large and banded. A straightforward computation for non-boundary pixels i shows that the stencil is given as in Fig. 2. For b ¼ À1=2, the corresponding stencils for the Laplacian as well as for the inpainting matrix perform an undesirable odd-even decoupling. As a remedy, values for b should be chosen in the range À1; À1=2 ½ Þ . If b ¼ À1, we obtain the well-known standard five point stencil for the Laplacian.
Following [5,13], we can use the stencil entries to obtain estimates on the c i for which invertibility of the system matrix is asserted below. This finding extends results in [5] with feasible range 0; 1 f g and from [13] with feasible range 0; 8=7 ½ .

Theorem 2
The inpainting matrix A c ð Þ from (8) corresponding to the stencil from Fig. 2 is invertible when all c i lie in the range 0; 4=3 ½ for h ¼ 1 and b ¼ À1=2. This is the largest possible range for the Laplacian from Fig. 1. Furthermore, the lower bound c i > 0 for the mask entries c i can only be asserted when b 2 À1; À1=2 ½ .
Proof We follow the ideas from [5,13]. By applying Geršgorin's circle theorem at the rows of the matrix A c ð Þ, we obtain pointwise estimates for the position of its eigenvalues. We note that all non-zero entries in any of the Fig. 1 Differential operators and corresponding stencils. The free parameter b stems from the fact that the nullspace of the matrix is one dimension For h ¼ 1, a lengthy but simple computation shows that 0\c i \8b=ð1 þ 8bÞ must hold whenever b 6 À1=2. Other values for b yield ranges which do not include the interval (0, 1] and thus are of no interest to us. Maximising the upper bound is possible for b ¼ À1=2 and yields the range (0, 4 / 3]. An identical analysis for the cases where a pixel is placed along a boundary or in a corner does not yield further restrictions. h Let us now demonstrate the benefits of larger feasible ranges for the mask values c. Our objective in this is to improve image reconstructions by IHD. Thus, we consider for quantitative evaluation a small synthetic image and measure the reconstruction error in function of the mask values. Our test image consists of 25 Â 25 pixels representing a sampled version of the function x 2 þ y 2 over 0; 1 ½ Â 0; 1 ½ . Our mask contains 5 non-zero entries. In each corner, we fix the mask value to be 1. The remaining non-zero entry is placed at the centre of the image and we measure the mean-squared error (MSE) of the reconstruction as a function of the mask value at this position. We study this setup for different values b, as shown in Fig. 3.
In a second experiment, we show that our discretisations allow a contrast enhancing effect. To this end, we select a small 10 Â 12 grey-scale image patch with pixel value 0.25 on the left half and 0.75 on the right half. Our masks consists of a small strip with 2 pixels width along the image edge. We set all non-zero mask values to the same value 8b=ð1 þ 8bÞ, and let b vary in the admissible range À1; À1=2 ½ . For each value of b, we consider the reconstruction at the two neighbouring pixels from each side of the image. Table 1 confirms the desired contrast improvements.

Conclusions
This paper shows the relation between the classic Helmholtz equation and the inpainting problem. Thus, we relate two up to now unconnected fields from science. In the future, we hope to develop better performing inpainting models based on this insight.