Abstract
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.
1. Introduction
The smooth function extension problem is a classical problem that has been studied extensively in the literature from various viewpoints. Some of the well-known results include Urysohn’s Lemma, the Tietze Extension, and Whitney’s Extension Theorem (see, e.g., [1–9]).
Inpainting was first introduced in [10] and then has been studied extensively by several authors (see, e.g., [11, 12]). 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 [12], 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 [12]. In [13], several surface completion methods have been studied. An optimal bound for the errors of the cubic inpainting method in [12] is given. Applications to smooth inpainting have also been discussed in [13]. 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 [13], the author proposed to use -harmonic functions for smooth surface completion and smooth surface inpainting. Later in [14], -harmonic functions are also studied for smooth surface completion and inpainting. The differences of the method using -harmonic functions in [13, 14] are as follows: the method in [13] uses , , as boundary data while the method in [14] uses , , as boundary data to solve for a -harmonic function. Here, is the function to be inpainted and 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 [13, 14]. In particular, we focus our study on smooth surface completion and smooth surface inpainting by biharmonic functions.
2. Surface Completion by Biharmonic Functions
Let be a simply connected region in with -boundary and be the diameter of . Let be a smooth function on any region containing . Assume that is known on a neighbourhood outside . The surface completion problem consists of finding a function on a region containing such that There are several ways to construct the function so that (1) holds. For smooth surface completion, one is often interested in finding a sufficiently smooth function satisfying (1).
An application of smooth surface completion is in smooth image inpainting. In smooth image inpainting, one has a smooth image which is known in a neighbourhood outside of a region while the data inside is missing. The goal of image inpainting is to extend the function over the region 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 , as the diameter of the inpainting region shrinks to 0, one has where is the image obtained from the inpainting scheme. Here, denotes the norm. Here and throughout, if is bounded uniformly by some constant .
Note that harmonic inpainting, that is, the extension found from the equation is a linear inpainting scheme [12].
In [12], the following result for cubic inpainting is proved.
Theorem 1 (cubic inpainting, Theorem 6.5 [12]). Let be the harmonic inpainting of . Let be a linear inpainting of on (not necessarily harmonic), and let be defined bywhere solves Poisson’s equationThen, defines a cubic inpainting of ; that is,
Remark 2. If is the harmonic inpainting of in , that is, solves the equation then the element defined by (4) is a biharmonic function which solves the following problem: In [13], this result is generalized to a multiresolution approximation extension scheme in which the Laplacian is replaced by more general lagged diffusivity anisotropic operators. It is proved in [13] that if solves the equation then A sharper error bound than (10) is obtained in [13].
Equation (9) can be written as a system of Poisson’s equations as follows: Thus, the problem of solving (9) is reduced to the problem of solving Poisson’s equations of the form Numerical methods for solving (12) have been extensively studied in the literature.
Note that the normal derivatives are not presented in (9). Thus, extension using (9) may not be differentiable across the boundary . To improve the smoothness of the extension across , it is suggested to find from the equation It is proved in [14] that if is the solution to (13), then (10) also holds.
Equation (13) cannot be reduced to a system of Poisson’s equations as (9). In fact, to solve (13), one often uses a finite difference approach which consists of finding discrete approximations to the operators and , . For “large” , it is quite complicated, though possible, to obtain these approximations.
As we can see from the above discussions, (9) is easier to solve numerically than (13). However, scheme (9) does not use any information about the normal derivatives of the surface on . Thus, the extension surface, obtained from scheme (9), may not be smooth across the boundary . On the contrary, (13) uses normal derivatives as boundary data and, therefore, is expected to yield better results than scheme (9) does.
In the next section, we will implement and compare the two surface completion schemes using (9) and (13). In particular, we focus our study on biharmonic functions which are solutions to (9) and (13) for .
3. Implementation
Let us discuss a numerical method for solving the equation To solve this equation, one often defines and solves for from the following system: Thus, the problem of solving (15) is reduced to the problem of solving the following Poisson’s equation: To solve (16), we use a 5-point finite difference scheme to approximate the Laplacian operator. This 5-point scheme is based on the following well-known formula:Here, is the discretization step size. This scheme is well defined at points , whose nearest neighbours are interior points of . If has a neighbour , then we use a stencil of the formIn the above formula, is the nearest neighbour to the left of . Similar formulae when is the nearest neighbour on the top, in the bottom, and to the right of can be obtained easily.
In our experiments, we choose as a square and the solution on the computation grid is presented as a vector. Using the 5-point finite difference scheme, (16) is reduced to the following algebraic system:where is the 5-point finite difference approximation to the Laplacian and is a vector containing boundary values of on at suitable entries. Matrix is a tridiagonal matrix; that is, all nonzero elements of lie on the main diagonal and the first diagonals above and below the main diagonal. Matrix 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 bi-Laplacian, we use a 13-point finite difference scheme which is based on the following formula (see [15]): This stencil is well defined for a grid point if all its nearest neighbours are in the interior of the domain. If has a neighbour and is on the left of , then we use the following formula:
Using the above finite difference scheme, (20) is reduced to a linear algebraic system of the form , where is a five-diagonal matrix. Numerical solutions to on the grid can be obtained by solving this linear algebraic system.
Before we proceed with numerical experiments, we need the following.
3.1. Quantitative 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 scikit-image processing in python unit completely at the disposal of the reader which readily provides quantitative correlations—in particular for the figures. The scikit-image processing package is at http://scikit-image.org/ and is a collection of open access algorithms for image processing with peer-reviewed code.
For our method, see http://scikit-image.org/docs/stable/api/skimage.restoration.html?highlight=biharmonicskimage.restoration.inpaint-biharmonic.
Moreover, for the benefit of the reader, many comparison methods with associated code are given in this unit; see some below, but see the full package for a longer list with references.(i)Denoise-bilateral (see [16]).(ii)Denoise-nl-means (see [17–19]).(iii)Denoise-tv-bregman, denoise-tv-chambolle, denoise-wavelet, estimate-sigma, inpaint-biharmonic (the present paper), nl-means-denoising, richardson-lucy, unsupervised-wiener, unwrap-phase, Wiener-Hunt deconvolution.
4. Numerical Experiments
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 with harmonic functions, the method with biharmonic functions in [12, 13], and the method with biharmonic functions in [14].
The function to be completed in our first experiment is Note that this function is a biharmonic function. The domain is discretized by a grid of size points. In the first experiment, we used .
In our numerical experiments, we denote by the extension by a harmonic function, by the biharmonic extension from [13], and by the biharmonic extension from [14].
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 [13, 14].
From Figures 1 and 2, one can see that the biharmonic reconstructions from [13, 14] are much better than the reconstruction by harmonic functions. The method in [13] in this experiment yields numerical results with accuracy a bit higher than the method in [14]. However, this does not imply that the method in [13] is better in terms of accuracy than the method in [14]. In the condition number of A, the finite difference approximation to the bi-Laplacian 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 [13] yields results with higher accuracy than the algorithm using the method in [14]. As we can see in later experiments, the method in [14] often gives better results than the method in [13]. The conclusion from this example is that both methods in [13, 14] yield numerical solutions at very high accuracy. The harmonic reconstruction in this experiment is not very good. This stems 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 is not a biharmonic function.
Figures 3 and 4 plot the errors of harmonic and biharmonic reconstructions. From these figures, it is clear that the method in [14] yields the best approximation. The harmonic reconstruction is the worst 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 [13, 14].
It can be seen from Figures 5 and 6 that the biharmonic extension method from [14] yields the best reconstruction. Although the biharmonic extension method from [13] 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 [13] in Figure 6 that the reconstruction may not be smooth, or even differentiable, across the boundary.
Appendix
Table 1 presents numerical results for the function defined by (24) and on the domain , . The diameter of is . From Table 1, one can see that the harmonic reconstruction has an order of accuracy of while the biharmonic reconstruction methods have an order of accuracy of 4. This agrees with the theoretical estimates in [12–14].
Figures 7 and 8 plot a damaged picture of peppers and reconstructed images by the 3 methods. From these figures, one gets the same conclusion as in the previous experiment. The biharmonic reconstruction in [14] yields the best restoration while the biharmonic reconstruction in [13] yields the worst reconstruction.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.