Modeling variational inpainting methods with splines

Mathematical methods of image inpainting involve the discretization of given continuous models. We present a method that avoids the standard pointwise discretization by modeling known variational approaches, in particular total variation (TV), using a finite dimensional spline space. Besides the analysis of the resulting model, we present a numerical implementation based on the alternating method of multipliers. We compare the results numerically with classical TV inpainting and give examples of applications.

1. Introduction. In this paper we investigate the modeling of a continuous inpainting method for digital images. Overviews of mathematical models/methods in image processing including inpainting are given in [4,9]. Due to the well-known work of Rudin and Osher [13], and Rudin, Osher, and Fatemi [14], functions of bounded variation are often considered in the continuous model. For a numerical solution and its implementation, the models need to be discretized. For digital images, given by a set of uniformly distributed pixels, those functions and their derivatives are typically discretized using difference schemes for the pixel values. In contrast to that, we use a finite dimensional function space here, specifically, the space of tensor product spline functions, see e.g., [7,15], and thus avoid a pointwise discretization.
In general, inpainting or filling in pursues the goal to restore parts Ω of an image Ω , Ω ⊂ Ω ⊂ R d , where the information has been removed, damaged, or is missing by using the information from the remaining and well-maintained part. There exists a variety of models for this purpose, see e.g., [1,2,8] and also [4,9]. We will incorporate a variational approach by, roughly speaking, minimizing some functional over an extended inpainting area U (Ω) ⊃ Ω subject to side conditions which ensure the reproduction of the well-maintained part of the image. Typically, the functional involves the image function u and some of its derivatives. The side condition, on the other hand, is chosen over some neighborhood B ⊆ Ω \ Ω of the inpainting area. We will model the class of functions for the variational method with tensor product B-splines together with a focus on TV inpainting, [5], that is |∇u(x x x)| dx x x, the effective, often used ROF-functional, see [14], which results in level curves of minimal length.
We begin by giving some basic properties and notation for the tensor product splines used in this paper in Section 2. Afterwards, in Section 3 we show how to model the inpainting problem with those splines and analyze its properties. Numerical results are given in Section 4. We present an implementation using the alternating method of multipliers (ADMM) [3], compare the results of our method to standard TV inpainting and show some examples of applications.

Preliminaries and notation.
In the sequel, we will use the following notations and basic facts about B-splines and their derivatives. Even if images are usually only considered in 2D, we will present the theory in an arbitrary number of variables; applications in higher dimensions would, for example include inpainting in voxel data provided by computerized tomography. The tensor product B-splines for inpainting which we will use in Section 3 are defined over a rectangle R = ⊗ d j=1 [a j , b j ], d ∈ N. They are based on a tensor product knot grid which is given as (1) T := ⊗ d j=1 {τ j,1 , . . . , τ j,mj +2nj } with τ j,i < τ j,i+nj for n j ≤ i ≤ m j + n j . At the boundary, we request multiple knots (2) τ j,1 = · · · = τ j,nj = a j and τ j,mj +nj +1 = · · · = τ j,mj +2nj = b j , for n j , m j ∈ N, 1 ≤ j ≤ d. We set n n n := (n 1 , . . . , n d ) and m m m := (m 1 , . . . , m d ) for the dimensions and degrees in the individual coordinates, respectively. The grid width h h h := (h 1 , . . . , h d ), defined as h j := max k |τ j,k − τ j,k+1 |, is known to influence the approximation properties of the spline space. A grid cell is denoted by Z k k k : The tensor product B-splines of order n n n with respect to the grid T are denoted by i.e., α α α ∈ I m m m,n n n := ⊗ d j=1 {1, . . . , m j + n j }. They have the support respectively. The spline space S n n n (T, Ω) of order n n n restricted to a domain Ω ⊂ R d is spanned by all B-splines B n n n α α α with S n n n α α α ∩ Ω = ∅. The index set of all B-splines relevant for Ω is defined as I Ω := {α α α ∈ I m m m,n n n |S n n n α α α ∩ Ω = ∅} and their number as #(I Ω ). If Ω = R, then I R = I m m m,n n n and #(I R ) = d j=1 m j + n j . A tensor product spline of order n n n with respect to T over a domain Ω is given by with coefficients f α α α ∈ R. Let f = f α α α α α α∈IΩ and B(x x x) T := B n n n α α α (x x x) α α α∈IΩ denote the vector of the coefficients and B-splines, respectively, then the spline s is given in matrix notation by With respect to inpainting methods an important property of the spline functions is that their derivatives can be expressed in terms of the coefficients of the spline function itself, that is denotes the derivative of the B-splines. Here, ε j is the j-th unit vector. The p -norm of vectors will be denoted by · p and the L p -norm of functions by · p,Ω .
3. Modeling of inpainting problem. Assume we want to reproduce a rectangular picture Ω ⊂ R d by a tensor product spline. This can easily be done by choosing a tensor product grid with multiple knots on the boundary ∂Ω such that Ω is the domain of definition. This method guarantees that the B-spline basis is stable.
, d ∈ N and any spline order n n n = (n 1 , . . . , n d ) ∈ N d , the B-spline basis {B n n n α α α (x x x)} α α α∈Im m m,n n n with respect to knots defined as in (1) and (2) is stable, that is, there exist constants c, C > 0 such that The constants are only depending on n n n and d.
Due to the tensor product structure and the multiple knots on the boundary, this result can be proven analogously to the classical univariate results for intervals, see, e.g., [6]. For further information about stability of tensor product B-spline bases see for example [11]. In addition to stability, multiple knots on the boundary of R help us to avoid artifacts at the boundary since they result in an interpolation at the boundary. For arbitrary Ω ⊂ R d , the inpainting problem considered here can now be described in the following way: Given two bounded domains Ω and Ω such that that Ω ⊂ R ⊂ Ω and ∂Ω ∩ ∂R = ∅, reconstruct a function or picture g on Ω that is known only on Ω * := Ω \ Ω by using tensor product splines defined over R. The situation is illustrated in Figure 1. We require that the spline fulfills a side condition on R * := R \ Ω ⊂ Ω * and minimizes a functional F over some neighborhood U (Ω). The general variational inpainting model using splines can now be formulated as follows.
SPLINE INPAINTING MODEL 1. Let U n n n h h h (Ω) be some neighborhood of Ω only depending on h h h and n n n with Ω ⊆ U n n n h h h (Ω) ⊆ R. The spline inpainting model is given by: Determine s ∈ S n n n (T, R) by REMARK 3.2. Since hard constraints can be problematic, especially in the presence of noise, the minimization problem is often relaxed to This relaxed formulation will be discussed briefly in Section 4.
For the explicit modeling or the discretization of the Spline Inpainting model 1, the concrete functional in the minimization is, of course, crucial. As already mentioned in Section 2, we will follow the most frequently used approach and focus on TV inpainting [5] and, therefore, on the ROF-functional [14]. It seems worthwhile, however, to mention that the spline approach also works with other functionals.
SPLINE INPAINTING MODEL 2. Let U n n n h h h (Ω) be some neighborhood of Ω only depending on h h h and n n n with Ω ⊆ U n n n h h h (Ω) ⊆ R. The TV spline inpainting model is given by: Determine s ∈ S n n n (T, R) by where In the next two subsections we show that the continuous Spline Inpainting Models 1-2 applied to discrete problems (e.g., digital images) is already a discrete model, except for the discretization of the integral. Clearly, there are different possible interpretations of a discrete image. In our situation, where we evaluate functions at points in R, those interpretations influence the value of the image that we assume to be at that points. Since it does not change the method or modeling, we interpret, for the sake of simplicity, the discrete values as piecewise constant functions with constant values over rectangles. Given an image g consisting of µ 1 × · · · × µ d pixels, those pixel rectangles are of the form P β β β : The value over P β β β is chosen corresponding to the value at the center c(P β β β ) := a j + (β j − 1 2 )p j d j=1 ; the j-th coordinate will be denoted by c(P β β β ) j . 3.1. Side condition. In the reproduction of digital images on R * , the side condition (7) corresponds to the interpolation of the pixel values in R * . In the space of tensor product splines an interpolation problem is solvable if the number of interpolation points is not greater than the degrees of freedom #(I R ) and if there is a relationship between the interpolation sites and the knots, known as the Schoenberg-Whitney condition. Since we also need some degrees of freedom for the minimization problem over U (Ω), this forces us to choose a spline space of sufficiently large dimension. DEFINITION 3.3. Given a tensor product grid T over R and a set of discrete points For functions defined over R and grids T satisfying (1) and (2) this interpolation problem is uniquely solvable if the Schoenberg-Whitney condition (see, e.g., [10]) is satisfied. Since R is a rectangle, this can be guaranteed by choosing the tensor product Greville abscissae as interpolation sites, that is, by setting It should be noted that the position of the Greville abscissae depends on the position of the knots of T , hence, the interpolation points depend on T . Therefore, the tensor product knot grid must be chosen according to the interpolation points. In case of digital images, where each rectangle P β β β represents a pixel of the image, we have the following result for the reproduction of all known image pixel values.
of known pixel values and an order n n n of a spline space, there exist a knot grid T and a set of interpolation points Ξ * such that for s ∈ S n n n (T, R).
Proof. Due to the tensor product structure, it suffices to consider a single coordinate direction j for the construction of the grid and set of interpolation points.
In doing so, we have to distinguish between even and odd values of n j . We set m j = µ j − 1 for n j odd and m j = µ j otherwise, and choose τ j,nj +kj := a j + k j · p j n j odd, τ j,nj +kj := a j + (k j − 1 2 ) · p j n j even for 1 ≤ k j ≤ m j and p j := bj −aj µj . The boundary knots are determined with proper multiplicity according to (2). Recalling (8), the j-th coordinates of the associated Greville abscissae for odd n j are given by Replacing the terms of the sums by − 1/2, results in the Greville abscissae for even n j , respectively. Therefore, we have for odd n j and for even n j , respectively. Now ξ γ γ γ,j coincides, for n j ≤ γ j ≤ m j + 1, with c(P β β β ) j for some 1 ≤ β j ≤ µ j ; more precisely, with centers c(P β β β ) j such that The Greville abscissae of boundary B-splines, that is, of B-splines with some knots at a j or b j , are not necessarily placed at a center, but they can be replaced by the center that is closest to them without loss of the Schoenberg-Whitney condition B nj γj (ξ j ) > 0, due to property (3). We now claim: For every β j ∈ {1, . . . , µ j } \I, there exists γ j such that There are four cases to be considered: small and large β j (or γ j ) and even or odd order. We prove the claim for small β j or γ j and odd order, the other cases can be verified in the same way. For γ j ∈ {2, . . . , n j − 1} and odd n j we have Therefore, we need to show that for every β j ∈ {1, . . . , (n j − 1)/2} there exists a γ j ∈ {2, . . . , n j − 1} such that nj −1 ≤ 1 and, thus, the claim is valid. Having this property at hand, we replace those ξ γ γ γ,j by c(P β ) j that are closest to the center, to guarantee that all centers are contained in the set of interpolation points. For β ∈ {1, . . . , µ j } \I set I := γ j ∈ {1, . . . , m j + n j } : ξ γ γ γ,j = miñ γj {argmin ξγ γ γ,j |ξγ γ γ,j − c(P β ) j |} andξ γ γ γ,j := argmin c(P β )j |ξγ γ γ,j − c(P β ) j | for γ j ∈Ĩ. The set Ξ j := ξ γ γ γ,j : γ j ∈ {1, . . . , m j + n j } \Ĩ ∪ ξ γ γ γ,j : γ j ∈Ĩ gives us the j-th coordinate of Ξ, the set of possible interpolation points. Finally, the set Ξ * of interpolation points is determined by restriction of Ξ to R * , that is ξ ξ ξ γ γ γ ∈ Ξ * := Ξ ∩ R * . By construction, c(P β β β ) ∈ Ξ * for β β β ∈ J R * and there exists a spline interpolant s ∈ S n n n (T, R) for g at Ξ * .

Minimization.
Depending on Ω, there are still some degrees of freedom left for minimization; their number is given by #(Ξ \ Ξ * ). We define U n n n T (Ω) as the union of all supports of B-splines corresponding to ξ α α α ∈ Ξ \ Ξ * , that is (10) U n n n T (Ω) := α α α∈I R \I Ξ * S n n n α α α .
The index set of all cells in U n n n T (Ω) is denoted by K Ω := {k k k ∈ K| Z k k k ⊂ U n n n T (Ω)}. PROPOSITION 3.5. Given T and Ω ⊂ R. For U n n n T (Ω) defined as in (10), the minimization problem (6) reduces to REMARK 3.6. It is worthwhile to point out that the derivatives ∂ j B n n n α α α (x x x) are just weighted differences of B-splines of lower order, the explicit formula being given in Section 2. This fact is of crucial importance for the efficient implementation of the method as it leads to sparse difference matrices.
Proof of Proposition 3.5. The result follows by applying (4) and (10). The first one directly gives us ∇s( with unknown coefficients f . The area of integration U n n n T (Ω) consists of a union of grid cells and, thus, the integral can be split up into those grid cells Z k k k for k k k ∈ K Ω .
Typically, for digital images the integration and the derivatives need to be discretized, as already mentioned in Section 1. In our case, we only need to discretize the integral and the way the problem is modeled helps with that, too: the discretization of the integrals can be done simply and efficiently by using a Gaussian quadrature formula for the grid cells Z k k k for k k k ∈ K Ω . The tensor product structure allows us to use the univariate Gauss-Legendre quadrature formula in each coordinate direction. Let Θ ⊂ U n n n T (Ω) be the set of the nodes of the Gauss-Legendre quadrature and w θ θ θ be the corresponding weights. We get where · 2 denotes the Euclidean norm. DISCRETE SPLINE INPAINTING MODEL 1. Using Proposition 3.5 and 3.4 (or (9)) the discrete form of TV Spline Inpainting Model 2 is given by: Determine s ∈ S n n n (T, R) by . This approach is also suitable for other types of Spline Inpainting Models 1, as long as the functionals depend on the function and its derivatives and it can be applied as soon as the explicit functional is given.

Numerical implementation and experiments. Define the convex operators
∞ otherwise, and the linear operator K : R #(IΩ) → R 2×#(Θ) with K(f ) = (w θ θ θ B θ θ θ f ) θ θ θ∈Θ . Using this, the Discrete Spline Inpainting Model 1 can be reformulated as an optimization problem in the following way: OPTIMIZATION PROBLEM 1. Our Spline Inpainting Model 1 is equivalent to the unconstrained problem This formulation is suitable for the alternating method of multipliers (ADMM) [3]. Therefore, the λ-proximity operators prox λF * , prox λG need to be calculated. Note that the convex conjugate function F * : R d×#(Θ) → R is given by Both, F * and G are indicator functions and, thus, the proximity operators are equivalent to the projections: prox λF * (y 1 , . . . , y #(Θ) ) = y 1 max (1, y 1 2 ) , . . . , where B + Ξ * is the pseudoinverse of B Ξ * . We remark that both operators are independent of λ. In our experiments, we use the MATLAB implementation of ADMM due to Gabriel Peyre [12] and apply the tests to several cartoon-like images and natural images in different sizes; some are shown in Figure 2.
We consider two different scenarios for the inpainting region Ω. In one scenario, the image is damaged by one or several "scratches" of variable width; in the second case we consider randomly missing pixels similar to "salt-and-pepper" noise. Figure 3 gives an example for both types (in comparison to Figure 2 the contrast is changed to emphasize the inpainting area).
In the next subsection we discuss the optimal spline order for the inpainting problem as well as reasonable strategies to guess a starting value for ADMM. We compare our results with the standard TV inpainting using, again, an implementation by Gabriel Peyre [12]. Afterwards, we demonstrate the benefits of our method by two examples: Text removal and salt-and-pepper denoising.

Numerical evaluation.
4.1.1. Starting guess. Both implementations, spline inpainting using ADMM and standard TV inpainting, are iterative methods. Thus, they rely on a suitable initial value to return a good solution after a reasonable number of iterations. We tried several strategies among which two stood out and will be detailed in what follows. In a first strategy we choose random uniformly distributed values in [0, 255] (for standard grayscale images). Note that  standard TV inpainting directly iterates on the pixel values of the inpainting area while our spline approach works on the spline coefficients as its variables. Thus, although both algorithms use uniformly distributed values, the starting guess strategies differ slightly. For our second strategy, we calculate the mean value ω mean of the image outside the inpainting area. In contrast to that, standard TV sets the pixels inside the inpainting area to ω mean . Our spline approach uses the starting values f = prox λG (ω mean 1), where 1 is a vector of ones. Since splines form a partition of unity, the coefficient vector ω mean 1 generates a constant image with value ω mean . Using prox λG , this vector is projected onto the interpolation space. Figure 4 illustrates the mean signal-to-noise ratio (SNR) for both starting strategies using 100 experiments on several images (128 × 128 pixels) and inpainting areas with 100 iterations of ADMM. We use two types of inpainting areas: on the one hand 3% of the pixels are randomly set to zero (left), on the other hand 3 scratches with a 4 pixel width are used (right). The mean SNR is calculated for cartoon-like images (top) and natural images (bottom) separately.
We see that the mean value starting guess performs much better on natural images and in case of cartoon like images the SNR values are at least of the same order as for a random starting guess. Only for standard TV on scratch inpainting areas the reconstruction quality of the random method outperforms the mean value. Thus, we suggest to use the mean value starting guess for spline inpainting, while for standard TV a random guess on cartoon-like images and the mean value on natural images should be chosen. This will also be the strategies used in the following experiments. Figure 5 shows the reconstruction of the "Lena" image shown in Figure 3 with randomly missing pixels. We use splines of order 2; first with random starting guess and a second time with mean value. The obtained SNR values are 20.00 and 26.47, respectively. Figure 4 already makes clear that a higher spline order not necessarily yields a better reconstruction. Indeed, the optimal order seems to depend on the image type as well as on the inpainting area. Higher order splines are preferable for natural images, which usually have a more complex structure, as well as for inpainting regions that produce larger gaps, e.g., scratches. For more simple structures and smaller inpainting regions lower order splines return comparable or even better results. This will be analyzed in more details in the next experiment.

Spline order.
Again, we calculate the mean SNR over 100 runs using 128 × 128 pixel images. This time, we plot the SNR against the percentage of unknown pixels (random inpainting area) or scratch width (scratch inpainting area) for both image types. Figure 6 shows the obtained SNR for different spline orders and standard TV as comparison. For both algorithms, we use the optimal starting guess just discussed.
Obviously, natural images are better reconstructed using splines of order 3 or higher. Especially, when the inpainting area is a connected scratch, higher orders can improve the result. When the image geometry becomes more simple, as e.g., for cartoon-like images, splines of order 2 − 3 perform best. For this reason, we recommend to use splines of order 2 or 3 for cartoon-like images and order 3 to 4 for natural images, depending on the inpainting region. Note that standard TV inpainting performs comparable to spline inpainting with order 2 for natural images. For cartoon-like images and randomly missing pixels it even performs better than all splines of order ≥ 3, but is outperformed by spline order 2. Figure 7 gives an example how the reconstruction quality can increase if a higher order spline is used on natural images. The scratched image was reconstructed using splines of order 2 and 3 obtaining a SNR of 25.73 and 30.67. As a comparison, the SNR in case of standard TV is 25.72. 4.2. Application 1: text removal. An application of inpainting methods is the removal of unwanted objects in images, such as text. In Figure 8 we illustrate how an image might be covered by an advert or other text. Here, we use the "Lena" image of size 256 × 256 pixels. We compare the TV reconstruction with our spline approach of order 3, 4 and 5. Since we use a natural image and the text inpainting area creates large gaps in the data, the reconstruction quality strongly profits from a higher spline order. The reconstruction has an SNR of only 23.2 for spline order 3, but increases to 28.45 for order 4 and even 31.25 for order 5. Standard TV inpainting results in an SNR of 23.58.

Application 2:
Salt-and-pepper noise. As a last example, we consider an image that is corrupted by salt-and-pepper noise, i.e., some pixels are randomly set to the minimal or maximal value. Since this noise deletes all information about the original pixel value, we can model the reconstruction as inpainting problem where the inpainting region is given by all pixels with minimal or maximal value. This coincides with our random inpainting area model. So far, we discussed noiseless data outside the inpainting domain. Therefore, the operator G was chosen such that only interpolating solutions were allowed. However, when dealing with noisy images it is not unlikely that the image is corrupted by several types of noise. Hence, we now use the relaxed model (5). We can still use the ADMM, but now with the operator G(f ) = ε 2 B Ξ * f − g Ξ * 2 2 and its proximity operator where I is the identity matrix. Note that we cannot use the relaxed model (5) to denoise the image far away from the inpainting domain since the TV minimization is only performed in a surrounding of the unknown pixels. Hence, for a good reconstruction we need to choose a higher spline order such that the spline support enlarges. Moreover, a large constant ε should be chosen, otherwise the solution of the relaxed model (5) is nearly constant. A simultaneous denoising may be obtained by expanding the TV term in model (5) to the complete image. However, this will drastically increase the number of Gauss-Legendre points used for the quadrature formula and, thus, increase the computational complexity. Figure 9 shows the result of spline inpainting with order 4 using the relaxed model on an image that was corrupted by Gaussian and salt-and-pepper noise. The parameter ε is chosen such that the best SNR is obtained. This is illustrated in Figure 10 where the SNR for different choices of ε is shown. We notice that for small ε the solution is nearly constant while a too large parameter can reproduce more of the noise in the reconstruction. Figure 9 shows the reconstruction for ε = 50 which leads to an SNR of 16.7. As a comparison, the SNR of the noised image is 4.04 in case of Gaussian and salt-and-pepper noise, and 15.69 if only Gaussian noise is considered.

Conclusion.
We presented a new approach to model the discrete inpainting problem using TV regularization and splines. The spline order can be chosen adapted to the underlying image and inpainting area. The advantages of this method were demonstrated in numerical experiments. Especially, when the images are complex and/or the inpainting domains are large, the reconstruction quality can highly profit from the new concepts by choosing a higher spline order. But also for cartoon-like images and low spline orders the method returns Although the method increases the reconstruction quality, it is still based on TV minimization and, thus, cannot overcome the known issues of TV inpainting completely. Therefore, considering other functionals is future work. Moreover, we want to extend the method on more general domains. Depending on the image, it might not be necessary to interpolate every pixel which could be exploited by the use of non-uniform knots (and interpolation points) and/or the combination with hierarchical splines.