# On Surface Completion and Image Inpainting by Biharmonic Functions: Numerical Aspects

Authors: S. B. Damelin; N. S. Hoang
Journal: International Journal of Mathematics and Mathematical Sciences
 (2018)

---

## 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.


---

## Body

## 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 p-harmonic functions for smooth surface completion and smooth surface inpainting. Later in [14], p-harmonic functions are also studied for smooth surface completion and inpainting. The differences of the method using p-harmonic functions in [13, 14] are as follows: the method in [13] uses ΔiuS, i=0,1,…,p-1, as boundary data while the method in [14] uses ∂i/∂NiuS, 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 [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

LetD be a simply connected region in R2 with C1-boundary S=∂D and d be the diameter of D. Let u0 be a smooth function on any region containing D. Assume that u0 is known on a neighbourhood outside D. The surface completion problem consists of finding a function u on a region containing D such that (1)u=u0outsideD.There are several ways to construct the function u so that (1) holds. For smooth surface completion, one is often interested in finding a sufficiently smooth function u satisfying (1).An application of smooth surface completion is in smooth image inpainting. In smooth image inpainting, one has a smooth imageu0 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 imageu0, as the diameter d of the inpainting region D shrinks to 0, one has (2)u-u0D=Od2,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, that is, the extension found from the equation(3)Δu=0inD,uS=u0S,is a linear inpainting scheme [12].In [12], the following result for cubic inpainting is proved.Theorem 1 (cubic inpainting, Theorem  6.5 [12]).
Letu1 be the harmonic inpainting of u0. Let ul be a linear inpainting of Δu0 on D (not necessarily harmonic), and let u(x) be defined by(4)ux=u1x+u2x,x∈D,where u2 solves Poisson’s equation(5)Δu2=ul,inD,u2S=0.Then, u defines a cubic inpainting of u0; that is,(6)u-u0D=Od4.Remark 2.
Iful is the harmonic inpainting of Δu0 in D, that is, ul solves the equation (7)Δul=0inD,ulS=Δu0S,then the element u defined by (4) is a biharmonic function which solves the following problem: (8)Δ2u=0,ΔuS=Δu0S,uS=u0S.
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 u solves the equation (9)Δnu=0inD,ΔiuS=Δiu0S,i=0,1,…,n-1,then (10)u-u0D=Od2n.A sharper error bound than (10) is obtained in [13].
Equation (9) can be written as a system of Poisson’s equations as follows: (11)v0=0,inD,Δvi=vi-1inD,viS=Δn-iu0S,i=1,2,…,nu=vn.Thus, the problem of solving (9) is reduced to the problem of solving Poisson’s equations of the form (12)Δu=finD,uS=g.Numerical methods for solving (12) have been extensively studied in the literature.
Note that the normal derivatives∂i/∂Niu0S are not presented in (9). Thus, extension using (9) may not be differentiable across the boundary S. To improve the smoothness of the extension across S, it is suggested to find u from the equation (13)Δnu=0inD,∂i∂Niu=∂i∂Niu0,i=0,1,…,n-1.It is proved in [14] that if u 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 Δn and ∂i/∂Ni, i=1,…,n-1. For “large” n, 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 S. Thus, the extension surface, obtained from scheme (9), may not be smooth across the boundary S. 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 n=2.


## 3. Implementation

Let us discuss a numerical method for solving the equation(14)Δ2u=0inD,ΔuS=f,uS=g.To solve this equation, one often defines v=Δu and solves for u from the following system: (15)Δv=0,vS=f,Δu=v,uS=g.Thus, the problem of solving (15) is reduced to the problem of solving the following Poisson’s equation: (16)Δu=f1inD,uS=g1.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:(17)Δ5u≔1h211-411u=Δu+Oh2.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(18)Δ5uP≔1h21-411uP+uQh2=ΔuP+Oh2.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, in the bottom, and to the right of P can be obtained easily.In our experiments, we chooseD as a square and the solution u on the computation grid is presented as a vector. Using the 5-point finite difference scheme, (16) is reduced to the following algebraic system:(19)Au=f1-g~1h2,where A is the 5-point finite difference approximation to the Laplacian and g~1 is a vector containing boundary values of u on S at suitable entries. Matrix A is a tridiagonal matrix; that is, all nonzero elements of A lie on the main diagonal and the first diagonals above and below the main diagonal. Matrix -h2A can be obtained by the function delsq available in Matlab.Let us discuss a numerical method for solving the equation(20)Δ2u=0,uNS=f,uS=g.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]): (21)Δ132u=1h412-821-820-812-821u=Δ2u+Oh2.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:(22)Δ132uP≔1h412-82-820-812-821uP+uNQh3+uQh4=Δ2u+Oh-1.Using the above finite difference scheme, (20) 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 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 athttp://scikit-image.org/ and is a collection of open access algorithms for image processing with peer-reviewed code.For our method, seehttp://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.


## 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 athttp://scikit-image.org/ and is a collection of open access algorithms for image processing with peer-reviewed code.For our method, seehttp://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 functionu to be completed in our first experiment is (23)ux,y=xy+x2y+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 byuH the extension by a harmonic function, by uL the biharmonic extension from [13], and by uN the biharmonic extension from [14].Figure1 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].Figure 1Figure 2From Figures1 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 ofA, 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(24)ux,y=1+cos⁡x1+cos⁡y4,x,y∈-1,1.This function u(x,y) is not a biharmonic function.Figures3 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.Figure 3Figure 4


### 4.2. Image Inpainting

Let us do some numerical experiments with image inpainting.Figure5 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].Figure 5
Damaged image and restored image by harmonic functions.Figure 6
Restored image by biharmonic functions.It can be seen from Figures5 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.


## 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 functionu to be completed in our first experiment is (23)ux,y=xy+x2y+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 byuH the extension by a harmonic function, by uL the biharmonic extension from [13], and by uN the biharmonic extension from [14].Figure1 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].Figure 1Figure 2From Figures1 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 ofA, 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(24)ux,y=1+cos⁡x1+cos⁡y4,x,y∈-1,1.This function u(x,y) is not a biharmonic function.Figures3 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.Figure 3Figure 4


## 4.2. Image Inpainting

Let us do some numerical experiments with image inpainting.Figure5 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].Figure 5
Damaged image and restored image by harmonic functions.Figure 6
Restored image by biharmonic functions.It can be seen from Figures5 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.



---

*Source: 3950312-2018-02-27.xml*