On surface completion and image inpainting by biharmonic functions: Numerical aspects

Numerical experiments with smooth surface extension and image inpainting using harmonic and biharmonic functions are carried out. The boundary data used for constructing biharmonic functions are the values of the Laplacian and normal derivatives of the functions on the boundary. Finite difference schemes for solving these harmonic functions are discussed in detail.


Introduction
The smooth function extension problem is a classical problem and has been studied extensively in the literature from various view points. Some of the well-known results include the Uryshons Lemma, the Tietze Extension, and the Whitneys Extension Theorem (see, e.g, [20,21,4,10,11,14,13,18,23]).
Inpainting was first introduced in [1] and then have been studied extensively by several authors (see, e.g., [2], [8]). Although smooth image inpainting is a smooth function extension problem, the most common approach in image inpainting so far is to use the solution to some PDE which are obtained from minimum-energy models as the recovered image. The most commonly used density function for these energy models is total derivation.
In [8] by considering smooth inpainting as a smooth surface extension problem, the author studied methods for linear inpainting and cubic inpainting. Error bounds for these inpainting methods are derived in [8]. In [6], several surface completion methods have been studied. An optimal bound for the errors of the cubic inpainting method in [8] is given. Applications to smooth inpainting have also been discussed in [6]. There, error bounds of completion methods are derived in terms of the radius of the domains on which the functions are completed. In one of the methods in [6], the author proposed to use p-harmonic functions for smooth surface completion and smooth surface inpainting. Later in [7], p-harmonic functions are also studied for smooth surface completion and inpaiting. The differences of the method using pharmonic functions in [7] and [6] are: the method in [6] uses △ i u| S , i = 0, 1, .., p − 1 as boundary data while the method in [7] uses ∂ i ∂N i u| S , i = 0, 1, ..., p−1 as boundary data to solve for a p-harmonic function. Here u is the function to be inpainted and S is the boundary of the inpainted region.
The goals of this paper are to implement and compare the performance of the two surface completion schemes in [6] and [7]. In particular, we focus our study on smooth surface completion and smooth surface inpainting by biharmonic functions.

Surface completion by biharmonic functions
Let D be a simply connected region in R 2 with C 1 -boundary S = ∂D and d be the diameter of D. Let u 0 be a smooth function on any region containing D. Assume that u 0 is known on a neighbourhood outside D. The surface completion problem consists of finding a function u on a region containing D such that There are several ways to construct the function u so that (2.1) holds. For smooth surface completion, one often is interested in finding a sufficiently smooth function u satisfying (2.1). An application of smooth surface completion is in smooth image inpainting. In smooth image inpainting, one has a smooth image u 0 which is known in a neighbourhood outside of a region D while the data inside D is missing. The goal of image inpainting is to extend the function u over the region D in such a way that the extension over the missing region is not noticeable with human eyes.
In image inpainting, an inpainting scheme is said to be of linear order, or simply linear inpainting, if for any smooth test image u 0 , as the diameter d of the inpainting region D shrinks to 0, one has where u is the image obtained from the inpainting scheme. Here · D denotes the L ∞ (D) norm. Here and throughout f = O(g) if |f | |g| is bounded uniformly by some constant C > 0.
Note that harmonic inpainting, i.e., the extension found from the equation is a linear inpainting scheme ( [8]). In [8] the following result for cubic inpainting is proved Theorem 2.1 (Cubic Inpainting, Theorem 6.5 [8]). Let u 1 be the harmonic inpainting of u 0 . Let u ℓ be a linear inpainting of △u 0 on D (not necessarily being harmonic), and let u(x) be defined by where u 2 solves the Poisson's equation Then u defines a cubic inpainting of u 0 , i.e., Remark 2.2. If u ℓ is the harmonic inpainting of △u 0 in D, i.e., u ℓ solves the equation Then the element u defined by (2.2) is a biharmonic function which solves the following problem In [6], this result is generalized to a multi-resolution approximation extension scheme in which the Laplacian is replaced by more general lagged diffusivity anisotropic operators. It is proved in [6] that if u solves the equation A sharper error bound than (2.7) is obtained in [6]. Equation (2.6) can be written as a system of Poisson's equations as follows Thus, the problem of solving (2.6) is reduced to the problem of solving Poisson's equations of the form Numerical methods for solving equation (2.9) have been extensively studied in the literature.
Note that the normal derivatives ∂ i ∂N i u 0 | S are not presented in equation (2.6). Thus, the extension using (2.6) may be not differentiable across the boundary S. To improve the smoothness of the extension across S, it is suggested to find u from the equation It is proved in [7] that if u is the solution to equation (2.10), then (2.7) also holds. Equation (2.10) cannot be reduced to a system of Poisson's equations as (2.6). In fact, to solve (2.10) one often uses a finite difference approach which consists of finding discrete approximations to the operators △ n and ∂ i ∂N i , i = 1, ..., n − 1. For 'large' n, it is quite complicated, although possible, to obtain these approximations.
As we can see from the above discussions, equation (2.6) is easier to solve numerically than equation (2.10). However, scheme (2.6) does not use any information about the normal derivatives of the surface on S. Thus, the extension surface, obtained from scheme (2.6), may not be smooth cross the boundary S. On the contrary, equation (2.10) uses normal derivatives as boundary data and, therefore, is expected to yield better results than scheme (2.6) does.
In the next section we will implement and compare the two surface completion schemes using equations (2.6) and (2.10). In particular, we focus our study on biharmonic functions which are solutions to (2.6) and (2.10) for n = 2.

Implementation
Let us discuss a numerical method for solving the equation To solve this equation, one often defines v = △u and solves for u from the following system Thus, the problem of solving (3.2) is reduced to the problem of solving the following Poisson's equation To solve equation (3.3), we use a 5-points finite difference scheme to approximate the Laplacian operator. This 5-points scheme is based on the following well-known formula Here, h is the discretization step-size. This scheme is well-defined at points P , whose nearest neighbours are interior points of D. If P has a neighbour Q ∈ ∂D, then we use a stencil of the form In the above formula Q is the nearest neighbour to the left of P . Similar formulae when Q is the nearest neighbour on the top, on the bottom, and to the right of P can be obtained easily.
In our experiments, we choose D as a square and the solution u on the computation grid is presented as a vector. Using the 5-points finite difference scheme above equation (3.3) is reduced to the following algebraic system where A is the 5-points finite difference approximation to the Laplacian andg 1 is a vector containing boundary values of u on S at suitable entries. The matrix A is is a tridiagonal matrix, that is, all non-zero elements of A lie on the main diagonal, and the first diagonals above and below the main diagonal. The matrix −h 2 A can be obtained by the function delsq available in Matlab. Let us discuss a numerical method for solving the equation For a discrete approximation to the bilaplacian we use a 13-points finite difference scheme which is based on the following formula (see [3]) This stencil is well-defined for a grid point P if all its nearest neighbours are in the interior of the domain. If P has a neighbour Q ∈ ∂D and Q is on the left of P , then we use the following formula Using the above finite difference scheme equation (3.7) is reduced to a linear algebraic system of the form Au = b where A is a five-diagonal matrix. Numerical solutions to u on the grid can be obtained by solving this linear algebraic system.
Before we proceed with numerical experiments we need: 3.1. Quantative comparisons. It is constructive to provide quantitative correlations between original and processed images and in particular code to compare figures such as those below. In order to calculate these required correlations (and many are provided), we refer the reader to a free access code for our method in a sckit-image processing in python unit competely at the disposal of the reader which readily provides quantative correlations-in particular for the figures below. The sckit-image-image processing package is at: http://scikit-image.org/ and is a collection of open access algorithms for image processing with peer-reviewed code.
Moreover, for the benefit of the reader, many comparison methods with associate code are given in this unit, see below for some, but see full package for a longer list with references.

4.1.
Smooth surface completion. Let us first do some numerical experiments with smooth surface completion. In our experiments, we compare numerical solutions from the three surface completion methods: the method by harmonic functions, the method by biharmonic functions in [8] and [6], and the method with biharmonic functions in [7]. The function u to be completed in our first experiment is u(x, y) = xy + x 2 (y + 1), x, y ∈ D = [−1, 1].
Note that this function is a biharmonic function. The domain D is discretized by a grid of size (n + 1) × (n + 1) points. In the first experiment we used n = 50.
In our numerical experiments, we denote by u H the extension by a harmonic function, by u L the biharmonic extension from [6], and by u N the biharmonic extension from [7]. Figure 1 plots the original function and the error of reconstruction by a harmonic function. Figure 2 plots the errors of surface reconstructions by biharmonic functions from [6] and [7].
From Figures 1 and 2, one can see that the biharmonic reconstructions from [6] and [7] are much better than the reconstruction by harmonic functions. The method in [6] in this experiment yields numerical results with accuracy a bit higher than the method in [7]. However, this doesn't imply that the method in [6] is better in term of accuracy than the method in [7]. The condition number of A, the finite difference approximation to the bilaplacian in this experiment is larger than that of the finite difference approximation to the Laplacian. Due to these condition numbers, the algorithm using the method in [6] yields results with higher accuracy then the algorithm using the method in [7]. As we can see in later experiments, the method in [7] often gives better results than the method in [6]. The conclusion from this example is that both the methods in [6] and [7] yield numerical solutions at very high accuracy. The harmonic reconstruction in this experiment is not very good. This comes from the fact that the function to be reconstructed is not harmonic.
In the next experiment the function to be reconstructed is chosen by This function u(x, y) is not a biharmonic function. Figure 3 and 4 plot the errors of harmonic and biharmonic reconstructions. From these figures it is clear that the method [7] yields the best approximation. The harmonic reconstruction is the worse amongst the 3 methods in this experiment.

4.2.
Image inpainting. Let us do some numerical experiments with image inpainting. Figure 5 plots a damaged image and a reconstructed image by harmonic functions. Figure 6 plots restored images by biharmonic functions following the methods from [6] and [7].
It can be seen from Figures 5 and 6, that the biharmonic extension method from [7] yields the best reconstruction. Although the biharmonic extension method from [6] is better than harmonic extension in our experiments with smooth surface completion, it is not as good as the harmonic extension in this experiment. This is understandable since our image contains edges and is not a smooth function. It can be seen from the restored image by the method in [6] in Figure 6 that the reconstruction may not be smooth, or even not differentiable, across the boundary.   Table 1, one can see that the harmonic reconstruction has an order of accuracy of 2 while the biharmonic reconstruction methods have an order of accuracy of 4. This agrees with the theoretical estimates in [8], [6], and [7].     [7] yields the best restoration while the biharmonic reconstruction in [6] yields the worst reconstruction.     Damaged Image Harmonic inpainting