Modeling variational inpainting methods with splines

TL;DR

This paper introduces a novel approach to image inpainting by modeling variational methods, especially total variation, using finite-dimensional spline spaces, providing an alternative to traditional discretization techniques.

Contribution

The paper develops a spline-based variational inpainting model and analyzes its properties, offering a new discretization method that improves upon classical pointwise approaches.

Findings

Spline-based inpainting yields comparable or better results than classical TV inpainting.

The proposed method effectively handles image inpainting tasks with smooth and textured regions.

Numerical experiments demonstrate the advantages of the spline approach over traditional discretization.

Abstract

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.

Figures

Equations

$${\int }_{U(\Omega )}|\nabla u(\mathbf{x})|d\mathbf{x},$$

$$T:={\otimes }_{j=1}^d\{{\tau }_{j,1},\dots ,{\tau }_{j,m_j+2n_j}\}$$

$${\tau }_{j,1}=\cdots ={\tau }_{j,n_j}=a_j\text{ and }{\tau }_{j,m_j+n_j+1}=\cdots ={\tau }_{j,m_j+2n_j}=b_j,$$

$$B_{\mathbf{\alpha }}^{\mathbf{n}}(\mathbf{x})=\prod _{j=1}^d{B}_{{\alpha }_j}^{n_j}(x_j),\mathbf{\alpha }\le \mathbf{m}+\mathbf{n},$$

$$S_{\mathbf{\alpha }}^{\mathbf{n}}:={\otimes }_{j=1}^d{S}_{{\alpha }_j}^{n_j},S_{{\alpha }_j}^{n_j}:=[{\tau }_{j,{\alpha }_j},{\tau }_{j,{\alpha }_j+n_j}],$$

$$s(\mathbf{x})=\sum _{\mathbf{\alpha }\in {\mathrm{I}}_{\Omega }}{f}_{\mathbf{\alpha }}{B}_{\mathbf{\alpha }}^{\mathbf{n}}(\mathbf{x}),\mathbf{x}\in \Omega ,$$

$${\partial }_{j}s(\mathbf{x}):=\frac{\partial }{\partial x_j}s(\mathbf{x})=\sum _{\mathbf{\alpha }\in {\mathrm{I}}_{\Omega }}{f}_{\mathbf{\alpha }}{\partial }_j{B}_{\mathbf{\alpha }}^{\mathbf{n}}(\mathbf{x})$$

$${\partial }_j{B}_{\mathbf{\alpha }}^{\mathbf{n}}(\mathbf{x}):=\frac{\partial }{\partial x_j}{B}_{\mathbf{\alpha }}^{\mathbf{n}}(\mathbf{x})=(n_j-1)\left(\frac{B_{\mathbf{\alpha }}^{\mathbf{n}-{\epsilon }_j}(\mathbf{x})}{|S_{{\alpha }_j}^{n_j-1}|}-\frac{B_{\mathbf{\alpha }+{\epsilon }_j}^{\mathbf{n}-{\epsilon }_j}(\mathbf{x})}{|S_{{\alpha }_j+1}^{n_j-1}|}\right)$$

$$c\Vert \mathbf{f}{\Vert }_p\le \Vert \mathbf{Bf}{\Vert }_{p,\Omega }\le C\Vert \mathbf{f}{\Vert }_p,\mathbf{f}\in {\mathbb{R}}^{I_{\Omega }},$$

$$F(s,\nabla s,\dots )(\mathbf{x})$$

$$s(\mathbf{x})=g(\mathbf{x})$$

$$\text{minimize}F(s,\nabla s,\dots ){|}_{U_{\mathbf{h}}^{\mathbf{n}}(\Omega )}+\frac{\epsilon }{2}\Vert s-g{\Vert }_{2,R^{\ast }}^2\text{for some }\epsilon >0.$$

$${\int }_{U_{\mathbf{h}}^{\mathbf{n}}(\Omega )}|\nabla s(\mathbf{x})|d\mathbf{x}$$

$$s(\mathbf{x})=g(\mathbf{x})$$

$${\lambda }_{\mathbf{\beta }}(s)=\sum _{\mathbf{\alpha }\in {\mathrm{I}}_R}{f}_{\mathbf{\alpha }}{B}_{\mathbf{\alpha }}^{\mathbf{n}}({\mathbf{\xi }}_{\mathbf{\beta }})={\lambda }_{\mathbf{\beta }}(f),{\mathbf{\xi }}_{\mathbf{\beta }}\in \Xi .$$

$$\Xi_{Gr}:=\bigg{\{}\boldsymbol{\xi}_{\boldsymbol{\gamma}}\in\mathbb{R}^{d}|\xi_{\boldsymbol{\gamma},j}:=\frac{1}{n_{j}-1}\sum_{i=1}^{n_{j}-1}\tau_{j,\gamma_{j}+i},\,\boldsymbol{\gamma}\in\mathrm{I}_{R}\bigg{\}}.$$

$$\{c(P_{\mathbf{\beta }}){\}}_{\mathbf{\beta }\in {\mathrm{J}}_{R^{\ast }}}\subset {\Xi }^{\ast }\text{and}{\lambda }_{\mathbf{\gamma }}(s)={\lambda }_{\mathbf{\gamma }}(g),{\mathbf{\xi }}_{\mathbf{\gamma }}\in {\Xi }^{\ast },$$

$$\{\begin{array}{ll}{\tau }_{j,n_j+k_j}:=a_j+k_j\cdot p_j& n_j\text{ odd},\\ {\tau }_{j,n_j+k_j}:=a_j+(k_j-\frac{1}{2})\cdot p_j& n_j\text{ even}\end{array}$$

$$\xi_{\boldsymbol{\gamma},j}=\begin{cases}a_{j}&\text{for }\gamma_{j}=1,\\ a_{j}+\frac{1}{n_{j}-1}\big{(}\sum_{\ell=1}^{\gamma_{j}-1}\ell\big{)}\,p_{j}&\text{for }2\leq\gamma_{j}\leq n_{j}-1,\\ a_{j}+\frac{1}{n_{j}-1}\big{(}\sum_{\ell=\gamma_{j}+1-n_{j}}^{\gamma_{j}-1}\ell\big{)}\,p_{j}&\text{for }n_{j}\leq\gamma_{j}\leq m_{j}+1,\\ b_{j}-\frac{1}{n_{j}-1}\big{(}\sum_{\ell=1}^{m_{j}+n_{j}-\gamma_{j}}\ell\big{)}\,p_{j}\phantom{p.}&\text{for }m_{j}+2\leq\gamma_{j}\leq m_{j}+n_{j}-1,\\ b_{j}&\text{for }\gamma_{j}=m_{j}+n_{j}.\end{cases}$$

$$\xi_{\boldsymbol{\gamma},j}=\begin{cases}a_{j}+\frac{(\gamma_{j}-1)\gamma_{j}}{2(n_{j}-1)}\,p_{j}&\text{for }1\leq\gamma_{j}\leq n_{j}-1,\\ a_{j}+\big{(}\gamma_{j}-\frac{n_{j}}{2}\big{)}\,p_{j}&\text{for }n_{j}\leq\gamma_{j}\leq m_{j}+1,\\ b_{j}-\frac{(m_{j}+n_{j}-\gamma_{j}+1)(m_{j}+n_{j}-\gamma_{j})}{2(n_{j}-1)}\,p_{j}&\text{for }m_{j}+2\leq\gamma_{j}\leq m_{j}+n_{j}.\end{cases}$$

$$\xi_{\boldsymbol{\gamma},j}=\begin{cases}a_{j}&\text{for }\gamma_{j}=1,\\ a_{j}+\frac{(\gamma_{j}-1)^{2}}{2(n_{j}-1)}\,p_{j}&\text{for }2\leq\gamma_{j}\leq n_{j}-1,\\ a_{j}+\big{(}\gamma_{j}-\frac{n_{j}+1}{2}\big{)}\,p_{j}&\text{for }n_{j}\leq\gamma_{j}\leq m_{j}+1,\\ b_{j}-\frac{(m_{j}+n_{j}-\gamma_{j})^{2}}{2(n_{j}-1)}\,p_{j}&\text{for }m_{j}+2\leq\gamma_{j}\leq m_{j}+n_{j}-1,\\ b_{j}&\text{for }\gamma=m_{j}+n_{j}.\end{cases}$$

$${\beta }_j\in \mathring{\mathrm{I}}:=\{\begin{array}{ll}\{\frac{n_j+1}{2},\dots ,{\mu }_j-\frac{n_j-1}{2}\}& \text{for odd }{n}_j\\ \{\frac{n_j}{2},\dots ,{\mu }_j+1-\frac{n_j}{2}\}& \text{for even }{n}_j\end{array}.$$

$$-\frac{1}{2}{p}_j<{\xi }_{\mathbf{\gamma },j}-c(P_{\beta }{)}_j\le \frac{1}{2}{p}_j.$$

$$\xi_{\boldsymbol{\gamma},j}-c(P_{\beta})_{j}=\Big{(}\frac{(\gamma_{j}-1)\gamma_{j}}{2(n_{j}-1)}-\beta_{j}+\frac{1}{2}\Big{)}\cdot p_{j}.$$

$${\beta }_j-1<\frac{({\gamma }_j-1){\gamma }_j}{2(n_j-1)}\le {\beta }_j.$$

$$\tilde{\mathrm{I}}:=\Big{\{}\gamma_{j}\in\{1,\dots,m_{j}+n_{j}\}:\xi_{\boldsymbol{\gamma},j}=\min_{\tilde{\gamma}_{j}}\{\operatorname{argmin}_{\xi_{\tilde{\boldsymbol{\gamma}},j}}|\xi_{\tilde{\boldsymbol{\gamma}},j}-c(P_{\beta})_{j}|\}\Big{\}}$$

$${\tilde{\xi }}_{\mathbf{\gamma },j}:={\mathrm{argmin}}_{c(P_{\beta }{)}_j}|{\xi }_{\tilde{\mathbf{\gamma }},j}-c(P_{\beta }{)}_j|\text{for }{\gamma }_j\in \tilde{\mathrm{I}}.$$

$$\Xi_{j}:=\big{\{}\xi_{\boldsymbol{\gamma},j}:\,\gamma_{j}\in\{1,\dots,m_{j}+n_{j}\}\setminus\tilde{\mathrm{I}}\big{\}}\cup\big{\{}\tilde{\xi}_{\boldsymbol{\gamma},j}:\,\gamma_{j}\in\tilde{\mathrm{I}}\big{\}}$$

$${\mathbf{s}}_{{\Xi }^{\ast }}={\mathbf{B}}_{{\Xi }^{\ast }}\mathbf{f}={\mathbf{g}}_{{\Xi }^{\ast }}$$

$$\mathbf{g}_{\Xi^{\ast}}=\big{(}g(\xi_{\boldsymbol{\beta}})\big{)}_{\boldsymbol{\beta}\in{\Xi^{\ast}}}\in\mathbb{R}^{\#({\Xi^{\ast}})},\qquad\mathbf{B}_{\Xi^{\ast}}=\bigg{[}B_{\boldsymbol{\alpha}}^{\boldsymbol{n}}(\xi_{\boldsymbol{\beta}}):\begin{array}[]{l}\xi_{\boldsymbol{\beta}}\in{\Xi^{\ast}}\\ \boldsymbol{\alpha}\in\mathrm{I}_{R}\end{array}\bigg{]}\in\mathbb{R}^{\#({\Xi^{\ast}})\times\#(\mathrm{I}_{R})}.$$

$$U_T^{\mathbf{n}}(\Omega ):=\bigcup _{\mathbf{\alpha }\in {\mathrm{I}}_R\setminus {\mathrm{I}}_{{\Xi }^{\ast }}}{S}_{\mathbf{\alpha }}^{\mathbf{n}}.$$

Full text

{AMS}

41A15, 41A63, 49M25, 65D05, 65D07, 65K10

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 $\Omega$ of an image $\Omega^{\prime}$, $\Omega\subset\Omega^{\prime}\subset\mathbb{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(\Omega)\supset\Omega$ 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\subseteq\Omega^{\prime}\setminus\Omega$ 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

[TABLE]

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.

2 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=\otimes_{j=1}^{d}[a_{j},b_{j}]$, $d\in\mathbb{N}$. They are based on a tensor product knot grid which is given as

[TABLE]

with $\tau_{j,i}<\tau_{j,i+n_{j}}$ for $n_{j}\leq i\leq m_{j}+n_{j}$. At the boundary, we request multiple knots

[TABLE]

for $n_{j},\,m_{j}\in\mathbb{N}$, $1\leq j\leq d$. We set $\boldsymbol{n}:=(n_{1},\dots,n_{d})$ and $\boldsymbol{m}:=(m_{1},\dots,m_{d})$ for the dimensions and degrees in the individual coordinates, respectively. The grid width $\boldsymbol{h}:=(h_{1},\dots,h_{d})$, defined as $h_{j}:=\max_{k}|\tau_{j,k}-\tau_{j,k+1}|$, is known to influence the approximation properties of the spline space. A grid cell is denoted by $Z_{\boldsymbol{k}}:=\otimes_{j=1}^{d}[\tau_{j,k_{j}},\tau_{j,k_{j}+1}]$ for $\boldsymbol{k}\in\mathrm{K}:=\otimes_{j=1}^{d}\{1,\dots,m_{j}+2n_{j}-1\}$.

The tensor product B-splines of order $\boldsymbol{n}$ with respect to the grid $T$ are denoted by

[TABLE]

i.e., $\boldsymbol{\alpha}\in\mathrm{I}_{\boldsymbol{m},\boldsymbol{n}}:=\otimes_{j=1}^{d}\{1,\dots,m_{j}+n_{j}\}$. They have the support

[TABLE]

respectively. The spline space $\mathscr{S}^{\boldsymbol{n}}(T,\Omega)$ of order $\boldsymbol{n}$ restricted to a domain $\Omega\subset\mathbb{R}^{d}$ is spanned by all B-splines $B_{\boldsymbol{\alpha}}^{\boldsymbol{n}}$ with $S_{\boldsymbol{\alpha}}^{\boldsymbol{n}}\cap\Omega\neq\emptyset$. The index set of all B-splines relevant for $\Omega$ is defined as $\mathrm{I}_{\Omega}:=\{\boldsymbol{\alpha}\in\mathrm{I}_{\boldsymbol{m},\boldsymbol{n}}|S_{\boldsymbol{\alpha}}^{\boldsymbol{n}}\cap\Omega\neq\emptyset\}$ and their number as $\#(\mathrm{I}_{\Omega})$. If $\Omega=R$, then $\mathrm{I}_{R}=\mathrm{I}_{\boldsymbol{m},\boldsymbol{n}}$ and $\#(\mathrm{I}_{R})=\prod_{j=1}^{d}m_{j}+n_{j}$.

A tensor product spline of order $\boldsymbol{n}$ with respect to $T$ over a domain $\Omega$ is given by

[TABLE]

with coefficients $f_{\boldsymbol{\alpha}}\in\mathbb{R}$. Let \mathbf{f}=\big{(}f_{\boldsymbol{\alpha}}\big{)}_{\boldsymbol{\alpha}\in\mathrm{I}_{\Omega}} and \mathbf{B}(\boldsymbol{x})^{T}:=\big{(}B_{\boldsymbol{\alpha}}^{\boldsymbol{n}}(\boldsymbol{x})\big{)}_{\boldsymbol{\alpha}\in\mathrm{I}_{\Omega}} denote the vector of the coefficients and B-splines, respectively, then the spline $s$ is given in matrix notation by $s(\boldsymbol{x})=\mathbf{B}(\boldsymbol{x})\mathbf{f}$.

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

[TABLE]

where

[TABLE]

denotes the derivative of the B-splines. Here, $\varepsilon_{j}$ is the $j$-th unit vector.

The $\ell^{p}$-norm of vectors will be denoted by $\|\cdot\|_{p}$ and the $L^{p}$-norm of functions by $\|\cdot\|_{p,\Omega}$.

3 Modeling of inpainting problem

Assume we want to reproduce a rectangular picture $\Omega^{\prime}\subset\mathbb{R}^{d}$ by a tensor product spline. This can easily be done by choosing a tensor product grid with multiple knots on the boundary $\partial\Omega^{\prime}$ such that $\Omega^{\prime}$ is the domain of definition. This method guarantees that the B-spline basis is stable.

Theorem 3.1**.**

For any domain $R=\otimes_{j=1}^{d}[a_{j},b_{j}]$, $d\in\mathbb{N}$ and any spline order $\boldsymbol{n}=(n_{1},\dots,n_{d})\in\mathbb{N}^{d}$, the B-spline basis $\{B_{\boldsymbol{\alpha}}^{\boldsymbol{n}}(\boldsymbol{x})\}_{\boldsymbol{\alpha}\in\mathrm{I}_{\boldsymbol{m},\boldsymbol{n}}}$ with respect to knots defined as in (1) and (2) is stable, that is, there exist constants $c,C>0$ such that

[TABLE]

*The constants are only depending on $\boldsymbol{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 $\Omega^{\prime}\subset\mathbb{R}^{d}$, the inpainting problem considered here can now be described in the following way: Given two bounded domains $\Omega$ and $\Omega^{\prime}$ such that

$\Omega\subset\Omega^{\prime}\subset\mathbb{R}^{d}$, 2. 2.

there exists $R=\otimes_{j=1}^{d}[a_{j},b_{j}]$ such that $\Omega\subset R\subset\Omega^{\prime}$ and $\partial\Omega\cap\partial R=\emptyset$,

reconstruct a function or picture $g$ on $\Omega$ that is known only on $\Omega^{\ast}:=\Omega^{\prime}\setminus\Omega$ 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^{\ast}:=R\setminus\Omega\subset\Omega^{\ast}$ and minimizes a functional $F$ over some neighborhood $U(\Omega)$. The general variational inpainting model using splines can now be formulated as follows.

Spline Inpainting Model 1**.**

Let $U_{\boldsymbol{h}}^{\boldsymbol{n}}(\Omega)$ be some neighborhood of $\Omega$ only depending on $\boldsymbol{h}$ and $\boldsymbol{n}$ with $\Omega\subseteq U_{\boldsymbol{h}}^{\boldsymbol{n}}(\Omega)\subseteq R$. The spline inpainting model is given by: Determine $s\in\mathscr{S}^{\boldsymbol{n}}(T,R)$ by

[TABLE]

Remark 3.2**.**

Since hard constraints can be problematic, especially in the presence of noise, the minimization problem is often relaxed to

[TABLE]

*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_{\boldsymbol{h}}^{\boldsymbol{n}}(\Omega)$ be some neighborhood of $\Omega$ only depending on $\boldsymbol{h}$ and $\boldsymbol{n}$ with $\Omega\subseteq U_{\boldsymbol{h}}^{\boldsymbol{n}}(\Omega)\subseteq R$. The TV spline inpainting model is given by: Determine $s\in\mathscr{S}^{\boldsymbol{n}}(T,R)$ by

[TABLE]

*where |\nabla s(\boldsymbol{x})|:=\sqrt{\sum_{i=1}^{d}\big{(}\partial_{i}s(\boldsymbol{x})\big{)}^{2}}. *

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 $\mu_{1}\times\cdots\times\mu_{d}$ pixels, those pixel rectangles are of the form $P_{\boldsymbol{\beta}}:=\otimes_{j=1}^{d}[a_{j}+(\beta_{j}-1)p_{j},a_{j}+\beta_{j}p_{j})$ for $p_{j}:=(b_{j}-a_{j})/\mu_{j}$ and $1\leq\beta_{j}\leq\mu_{j}$. The value over $P_{\boldsymbol{\beta}}$ is chosen corresponding to the value at the center c(P_{\boldsymbol{\beta}}):=\big{(}a_{j}+(\beta_{j}-\frac{1}{2})p_{j}\big{)}_{j=1}^{d}; the $j$-th coordinate will be denoted by $c(P_{\boldsymbol{\beta}})_{j}$.

3.1 Side condition

In the reproduction of digital images on $R^{\ast}$, the side condition (7) corresponds to the interpolation of the pixel values in $R^{\ast}$. 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 $\#(\mathrm{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(\Omega)$, 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 $\Xi:=\{\boldsymbol{\xi}_{\boldsymbol{\beta}}\in R|\,\boldsymbol{\beta}\in\mathbb{Z}^{d};\#(\boldsymbol{\beta})\leq\#(\mathrm{I}_{R})\}$, let $\lambda_{\boldsymbol{\beta}}$ be the point evaluation functionals $\lambda_{\boldsymbol{\beta}}(f):=f(\boldsymbol{\xi}_{\boldsymbol{\beta}})$, $f:\Xi\to\mathbb{R}$. A spline $s\in\mathscr{S}^{\boldsymbol{n}}(T,R)$ is the spline interpolant of a function $f$ at $\Xi$ if

[TABLE]

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

[TABLE]

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_{\boldsymbol{\beta}}$ represents a pixel of the image, we have the following result for the reproduction of all known image pixel values.

Theorem 3.4**.**

Given an image $g$ of size $\otimes_{j=1}^{d}\mu_{j}$, an area $R^{\ast}:=\cup_{\boldsymbol{\beta}\in\mathrm{J}_{R^{\ast}}}P_{\boldsymbol{\beta}}$, $\mathrm{J}_{R^{\ast}}:=\{\boldsymbol{\beta}\in\mathbb{Z}^{d}:1\leq\beta_{j}\leq\mu_{j},\,1\leq j\leq d\text{ and }c(P_{\boldsymbol{\beta}})\in R\setminus\Omega\}$ of known pixel values and an order $\boldsymbol{n}$ of a spline space, there exist a knot grid $T$ and a set of interpolation points $\Xi^{\ast}$ such that

[TABLE]

*for $s\in\mathscr{S}^{\boldsymbol{n}}(T,R)$. *

Proof 3.5**.**

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}=\mu_{j}-1$ for $n_{j}$ odd and $m_{j}=\mu_{j}$ otherwise, and choose

[TABLE]

for $1\leq k_{j}\leq m_{j}$ and $p_{j}:=\frac{b_{j}-a_{j}}{\mu—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

[TABLE]

Replacing the terms $\ell$ of the sums by $\ell-1/2$, results in the Greville abscissae for even $n_{j}$, respectively. Therefore, we have

[TABLE]

for odd $n_{j}$ and

[TABLE]

for even $n_{j}$, respectively. Now $\xi_{\boldsymbol{\gamma},j}$ coincides, for $n_{j}\leq\gamma_{j}\leq m_{j}+1$, with $c(P_{\boldsymbol{\beta}})_{j}$ for some $1\leq\beta_{j}\leq\mu_{j}$; more precisely, with centers $c(P_{\boldsymbol{\beta}})_{j}$ such that

[TABLE]

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_{\gamma_{j}}^{n_{j}}(\xi_{j})>0$, due to property (3).

We now claim: For every $\beta_{j}\in\{1,\dots,\mu_{j}\}\setminus\mathring{\mathrm{I}}$, there exists $\gamma_{j}$ such that

[TABLE]

There are four cases to be considered: small and large $\beta_{j}$ (or $\gamma_{j}$) and even or odd order. We prove the claim for small $\beta_{j}$ or $\gamma_{j}$ and odd order, the other cases can be verified in the same way. For $\gamma_{j}\in\{2,\dots,n_{j}-1\}$ and odd $n_{j}$ we have

[TABLE]

Therefore, we need to show that for every $\beta_{j}\in\{1,\dots,(n_{j}-1)/2\}$ there exists a $\gamma_{j}\in\{2,\dots,n_{j}-1\}$ such that

[TABLE]

Set $f(\gamma_{j}):=\frac{(\gamma_{j}-1)\gamma_{j}}{2(n_{j}-1)}$. We have $0<f(2)\leq 1$ and $\frac{n_{j}-1}{2}-1<f(n_{j}-1)\leq\frac{n_{j}-1}{2}$. As a consequence, $0<f(\gamma_{j})\leq\frac{n_{j}-1}{2}$ for all $\gamma_{j}\in\{2,\dots,n_{j}-1\}$. Clearly, $f(\gamma_{j})$ is monotonically increasing and $\#\gamma_{j}=n_{j}-2$ distributed over $(n_{j}-1)/2$ intervals. Further, $f(\gamma_{j})-f(\gamma_{j}-1)=\frac{\gamma_{j}-1}{n_{j}-1}\leq 1$ and, thus, the claim is valid.

Having this property at hand, we replace those $\xi_{\boldsymbol{\gamma},j}$ by $c(P_{\beta})_{j}$ that are closest to the center, to guarantee that all centers are contained in the set of interpolation points. For $\beta\in\{1,\dots,\mu_{j}\}\setminus\mathring{\mathrm{I}}$ set

[TABLE]

and

[TABLE]

The set

[TABLE]

*gives us the $j$-th coordinate of $\Xi$, the set of possible interpolation points. Finally, the set $\Xi^{\ast}$ of interpolation points is determined by restriction of ${\Xi}$ to $R^{\ast}$, that is $\boldsymbol{\xi}_{\boldsymbol{\gamma}}\in\Xi^{\ast}:={\Xi}\cap R^{\ast}$. By construction, $c(P_{\boldsymbol{\beta}})\in\Xi^{\ast}$ for $\boldsymbol{\beta}\in\mathrm{J}_{R^{\ast}}$ and there exists a spline interpolant $s\in\mathscr{S}^{\boldsymbol{n}}(T,R)$ for $g$ at $\Xi^{\ast}$. *

In the sequel, $\Xi$ and $\Xi^{\ast}$ always will refer to the sets of interpolation sites determined in the way described in the proof of Theorem 3.4. For applications we also want to describe the side conditions in matrix notation. For the coefficients \mathbf{f}=\big{(}f_{\boldsymbol{\alpha}}\big{)}_{\boldsymbol{\alpha}\in\mathrm{I}_{R}}\in\mathbb{R}^{\#(\mathrm{I}_{R})}, we get

[TABLE]

with

[TABLE]

3.2 Minimization

Depending on $\Omega$, there are still some degrees of freedom left for minimization; their number is given by $\#(\Xi\setminus\Xi^{\ast})$. We define $U_{T}^{\boldsymbol{n}}(\Omega)$ as the union of all supports of B-splines corresponding to $\xi_{\boldsymbol{\alpha}}\in\Xi\setminus\Xi^{\ast}$, that is

[TABLE]

The index set of all cells in $U_{T}^{\boldsymbol{n}}(\Omega)$ is denoted by $\mathrm{K}_{\Omega}:=\{\boldsymbol{k}\in\mathrm{K}|\,Z_{\boldsymbol{k}}\subset U_{T}^{\boldsymbol{n}}(\Omega)\}$.

Proposition 3.6**.**

Given $T$ and $\Omega\subset R$. For $U_{T}^{\boldsymbol{n}}(\Omega)$ defined as in (10), the minimization problem (6) reduces to

[TABLE]

*where \mathbf{f}=\big{(}f_{\boldsymbol{\alpha}}\big{)}_{\boldsymbol{\alpha}\in\mathrm{I}_{\Omega}}, $f_{\boldsymbol{\alpha}}\in\mathbb{R}$, and [\mathbf{B}^{j}(\boldsymbol{x})]^{T}:=\big{(}\partial_{j}B_{\boldsymbol{\alpha}}^{\boldsymbol{n}}(\boldsymbol{x})\big{)}_{\boldsymbol{\alpha}\in\mathrm{I}_{R}}\in\mathbb{R}^{\#(\mathrm{I}_{R})}. *

Remark 3.7**.**

*It is worthwhile to point out that the derivatives $\partial_{j}B_{\boldsymbol{\alpha}}^{\boldsymbol{n}}(\boldsymbol{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 3.8** (Proof of Proposition 3.6).**

*The result follows by applying (4) and (10). The first one directly gives us \nabla s(x)=\big{(}\mathbf{B}^{j}(\boldsymbol{x})\mathbf{f}\big{)}_{j=1}^{d} with unknown coefficients $\mathbf{f}$. The area of integration $U_{T}^{\boldsymbol{n}}(\Omega)$ consists of a union of grid cells and, thus, the integral can be split up into those grid cells $Z_{\boldsymbol{k}}$ for $\boldsymbol{k}\in\mathrm{K}_{\Omega}$. *

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_{\boldsymbol{k}}$ for $\boldsymbol{k}\in\mathrm{K}_{\Omega}$. The tensor product structure allows us to use the univariate Gauss-Legendre quadrature formula in each coordinate direction. Let $\Theta\subset U_{T}^{\boldsymbol{n}}(\Omega)$ be the set of the nodes of the Gauss-Legendre quadrature and $w_{\boldsymbol{\theta}}$ be the corresponding weights. We get

[TABLE]

where $\|\cdot\|_{2}$ denotes the Euclidean norm.

Discrete Spline Inpainting Model 1**.**

Using Proposition 3.6 and 3.4 (or (9)) the discrete form of TV Spline Inpainting Model 2 is given by: Determine $s\in\mathscr{S}^{\boldsymbol{n}}(T,R)$ by

[TABLE]

*where \mathbf{B}_{\boldsymbol{\theta}}:=\big{(}\mathbf{B}^{j}(\boldsymbol{\theta})\big{)}_{j=1}^{d}\in\mathbb{R}^{d\times\#(\mathrm{I}_{R})}. *

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.

4 Numerical implementation and experiments

Define the convex operators $F:\mathbb{R}^{d\times\#(\Theta)}\rightarrow\mathbb{R}$ and $G:\mathbb{R}^{\#(\mathrm{I}_{\Omega})}\rightarrow\mathbb{R}$ with

[TABLE]

and the linear operator $K:\mathbb{R}^{\#(\mathrm{I}_{\Omega})}\rightarrow\mathbb{R}^{2\times\#(\Theta)}$ with $K(\mathbf{f})=(w_{\boldsymbol{\theta}}\mathbf{B}_{\boldsymbol{\theta}}\mathbf{f})_{\boldsymbol{\theta}\in\Theta}$. 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

[TABLE]

This formulation is suitable for the alternating method of multipliers (ADMM) [3]. Therefore, the $\lambda$-proximity operators $\mathop{prox}_{\lambda F^{*}}$, $\mathop{prox}_{\lambda G}$ need to be calculated. Note that the convex conjugate function $F^{*}:\mathbb{R}^{d\times\#(\Theta)}\rightarrow\mathbb{R}$ is given by

[TABLE]

Both, $F^{*}$ and $G$ are indicator functions and, thus, the proximity operators are equivalent to the projections:

[TABLE]

where $\mathbf{B}_{\Xi^{\ast}}^{+}$ is the pseudoinverse of $\mathbf{B}_{\Xi^{\ast}}$. We remark that both operators are independent of $\lambda$. 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 $\Omega$. 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.

4.1 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 $\omega_{\text{mean}}$ of the image outside the inpainting area. In contrast to that, standard TV sets the pixels inside the inpainting area to $\omega_{\text{mean}}$. Our spline approach uses the starting values $\mathbf{f}=\mathop{prox}_{\lambda G}(\omega_{\text{mean}}{\mathbf{1}})$, where $\mathbf{1}$ is a vector of ones. Since splines form a partition of unity, the coefficient vector $\omega_{\text{mean}}{\mathbf{1}}$ generates a constant image with value $\omega_{\text{mean}}$. Using $\mathop{prox}_{\lambda 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\times 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.

4.1.2 Spline order

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.

Again, we calculate the mean SNR over $100$ runs using $128\times 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 $\geq 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\times 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$.

4.3 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(\mathbf{f})=\frac{\varepsilon}{2}\|\mathbf{B}_{\Xi^{\ast}}\mathbf{f}-\mathbf{g}_{\Xi^{\ast}}\|_{2}^{2}$ and its proximity operator

[TABLE]

where $\mathop{\mathbf{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 $\varepsilon$ 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 $\varepsilon$ is chosen such that the best SNR is obtained. This is illustrated in Figure 10 where the SNR for different choices of $\varepsilon$ is shown. We notice that for small $\varepsilon$ 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 $\varepsilon=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.

5 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 reasonable results. As a slight disadvantage, the new methods requires precalculation of the spline basis, Greville abscissae and Gauss-Legendre points what increases the runtime of the algorithm.

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.