1. Introduction. Image inpainting aims at restoring partially damaged image or missing parts of an image in a visually appealing manner. It has a wide number of practical applications such as art restoration [40] , object removal [23] , image compression [32, 63] , and even image denoising [2] . The term "inpainting" itself was introduced for digital images by Bertalmío et al. in [7] , but similar concepts were already explored in earlier work under different names such as image restoration, interpolation, or disocclusion [16, 30, 38, 52, 57] .
Any inpainting model needs to assume some kind of relation between known and unknown data. As there is a variety of plausible assumptions for such relations, many solutions to an inpainting problem exist.
Based on the underlying assumptions, the inpainting methods from the literature can be grouped into certain main categories [6, 34] . One class is based on variational models and partial differential equations (PDEs) [64] , comprising e.g. Euler's elastica [13, 17, 18, 52, 56] , transport-like equations [7, 11] , anisotropic diffusion processes [12, 32, 75] , harmonic and biharmonic inpainting [19, 32] , total variation restoration [19] , the Mumford-Shah functional [19] , and the Cahn-Hilliard equation [9, 15] . Exemplar-based approaches emerged from texture synthesis and exploit the notion of patch similarity [23, 25, 28] . Other techniques rely on overcomplete dictionaries and the concept of sparsity [26, 27, 51] , and more recently also deep learning concepts have been proposed [39, 58, 72] . Each of these strategies has its advantages and disadvantages depending on the type of image it is applied to. For example, exemplarbased techniques perform fairly well on highly textured images, PDE-based methods are more suited for geometrical structures, and deep learning approaches can capture high-level semantics from images. This has led to the development of hybrid approaches which combine the strengths of different methods [4, 8, 62, 69] .
A subclass of inpainting problems deals with the recovery of a whole image from a small amount of sparsely distributed data [1, 12, 28, 35] . These kind of problems are encountered particularly in the context of compression [21, 32, 60, 63] . The sparsity of available data makes it feasibly to consider scattered data interpolation, e.g. by radial basis functions [22, 43, 48, 71, 76] or by Shepard interpolation [1, 42, 60, 67] as inpainting technique.
A key observation for applications in compression is that the data can be chosen freely from the image. Thus, a careful selection of the sparse set of pixels to store such that it fits a chosen inpainting method is essential for a good performance [21, 32, 35, 41, 50, 63] . Astonishingly, simple linear methods such as homogeneous diffusion inpainting show remarkable quality if combined with optimally chosen data [5, 10, 21, 32, 36, 61] and can even compete with the widely used JPEG [59] and JPEG2000 [70] standards [37, 49, 61] . Anisotropic diffusion approaches perform even better [32, 35, 63] and can outperform JPEG and JPEG2000 for high compression ratios.
1.1. Goals and Contributions. The goal of our paper is to show that a hitherto hardly explored class of scattered data interpolation methods based on Smoothed Particle Hydrodynamics (SPH) can provide excellent results on sparse inpainting problems, if one improves them with a number of refined concepts.
SPH was originally introduced to solve astrophysical problems [47] , but has also been applied to problems that deal with large deformations [14] , computational fluid mechanics [54] , and soil mechanics [53] . In SPH, the solution to a given problem is represented by a set of particles and functions. Derivatives and integrals are approximated using those particles.
Di Blasi et al. [24] have introduced the SPH method for sparse image interpolation problems, and applications to non-sparse inpainting are studied in [3] . We have not found more work on SPH-based image inpainting. One reason for this lack of popularity might lie in the fact that in the naive formulation, not even constant functions are interpolated correctly. However, in our paper we show that one can come up with highly competitive approaches by integrating more sophisticated concepts. Our key contributions are the following: 1. We restore particle consistency and perform inpainting with a new method that adapts its consistency order to the local approximation error. 2. To define a measure for the area of influence of a given particle (mask point), we combine a Voronoi tessellation with the Euclidean distance transform. 3. We use this tessellation to propose a novel strategy for spatial data optimization. 4. We optimize not only the data locations, but also their values (tonal optimization). 5. To incorporate anisotropy in the process, we use anisotropic kernels. 6. We assess the performance of different smoothing kernels and compare to some of the best sparse inpainting methods for optimised data.
Structure. This paper is organized as follows: In section 2, we give a brief summary of the ideas behind Smoothed Particle Hydrodynamics. This includes its origin from an integral approximation, techniques to restore consistency in a discrete setting, and a brief overview on common smoothing kernels used for our experiments. Section 3 explains how SPH can be used for inpainting. Here we introduce Voronoi tessellation to determine parameters of the method and show how performance can be enhanced by combining results from methods of different consistency order. For further improvement, we explain our data optimization strategies both with respect to data locations (spatial optimization) as well as data values (tonal optimization) in section 4. We proceed by comparing results from our method to results from other techniques in subsection 4.5 and draw conclusions in section 5.
2. SPH in a Nutshell.
2.1. Essential Ideas. We are interested in approximating a function f on the domain Ω in R 2 . The point of departure for SPH is the idea to replace the value of f at a point q by a weighted average of the function, i.e.
Equation (2.1) is also known as the kernel approximation of the function with the smoothing kernel W (·, h) and its smoothing length h, which represents the effective width of W . The kernel should be a monotonically decreasing positive mollifier, i.e. it should have the following properties:
• monotonicity: W is monotonically decreasing function w.r.t q − p . Here, · denotes the Euclidean distance. Positivity is not strictly necessary, but desired in order for the approximated function values to have physical meaning. Allowing the kernel to take negative values in parts of the domain can lead to unnatural approximated values and corrupt the entire computation [45] . The same holds for monotonicity, which is connected with the usual behavior of physical forces to decrease with increasing distance.
Discretizing the integral of (2.1) yields the particle approximation of f given by
Here, in order to approximate the function value at point q we sum over its nearest neighbors p j , where N (q) denotes the index set of the nearest neighbors and we assume that 1 ≤ j ≤ M with M the total number of particles under consideration. Each of these particles is related to a specific area of influence (or weight) V j in this quadrature rule. The particle approximation is interpolating at the particles p j if the kernel W satisfies
with the Kronecker delta δ k,j . This requirement is separate from the desired properties in the continuous setting and, in general, not satisfied by kernels with said properties. Some modifications which achieve interpolation at the particles are discussed in subsection 2.2. In order to determine the nearest neighbors, we use the so-called scatter approach. Here, the neighbors of a point q, are the particles p j that include q in the support domain of the kernels centered at these particles p j , cf. Figure 1 . The scatter approach is preferable for our purpose as it is less dependent on the particle distribution than its alternatives [44] and allows to consider different smoothing lengths h j for each particle.
Consistency in SPH is defined in the sense of (local) polynomial reproduction [76] . While the restriction of kernels to positive mollifiers guarantees that linear polynomials are reconstructed in the continuous formulation (2.1), this is no longer the case in the discrete formulation (2.2); a phenomenon known as particle inconsistency [55] . Thus, (2.2) needs to be modified in order to restore consistency in the discrete setting.
For the particle approximation to satisfy zero order consistency, it needs to be able to reproduce constants. This requirement is satisfied by the well known Shepard interpolation formula [20, 67, 76] 
The fact that zero order consistency already requires to modify (2.2) shows that, in general, using the particle approximation directly cannot even reproduce constant functions, i.e. (2.2) has no consistency. One way to interpret Shepard interpolation is that the kernel W is replaced by a modified kernel W such that
and the original interpolation formula is used with the modified kernel, i.e.,
By comparison with (2.4), we see that b 0 (q) is given by
Thus, b 0 (q) is a function of q, but is constant with respect to the particle p j for a fixed q. In other words, zero order consistency can be restored by multiplying the kernel by a constant. This motivates the attempt to restore first order consistency by multiplying the kernel with a function which is linear in the difference q − p j , i.e.
where x q and y q denote the coordinates of q and analogous for p j . For a fixed q, this means that we have to determine three coefficients such that we can reproduce linear polynomials in R 2 . This yields the system of equations
Equation (2.9) allows to multiply the right-hand sides of (2.10) and (2.11) by the left-hand side of (2.9) without changing the equations, such that this linear system can be recast as
The matrix D(q) can be expressed as a sum of matrices of rank 1 in the form
It is positive semidefinite as for any z ∈ R 3 holds (2.15)
since the smoothing kernel W (·) is nonnegative and V j > 0. However, D(q) is singular unless q has at least three nearest neighbors which are not collinear. In order to achieve an SPH interpolation with first order consistency at a pixel q, it is necessary to solve the 3 × 3 linear system (2.12). To inpaint a whole image from first order SPH interpolation, (2.12) needs to be solved for each unknown pixel. With the solution b(q) of (2.12), a first order consistent SPH approximation can be written as
The particular method used here to restore first order consistency was derived in a longer way in [78] , whereas other methods that modify the kernel to restore first order consistency can be found in [44, 46] . The method described here has the advantage that it does not involve derivatives of f to restore first order consistency.
The benefit of a higher order consistency does not come without a prize. The zero order consistent Shepard interpolation (2.4) only modifies the kernel such that it satisfies a discrete partition of unity property. With this modified kernel, Shepard interpolation produces the value of u at q as a convex combination of the values of f at the neighboring particles p j . Thus, it prevents over-and undershoots. If we want the first order consistent method (2.16) to prevent over-and undershoots, we have to put a restriction on the positions of particles, since we have to satisfy (2.10) and (2.11). These equations can be written in a compact way as
In an ideal situation, u at q would be a convex combination of the values of f at the neighboring particles p j with weights given by W (q − p j , h) V j . However, (2.17) along with (2.9) implies that this can only be the case if the position q is also a convex combination of the particle positions p j with the same weights. In most cases, this condition on the positions of particles is violated. In order to achieve first order consistency regardless of the spatial distribution of particles, the modified kernel W violates some of the properties defined in subsection 2.1. In particular violation of the positivity requirement results in visible artifacts as can be seen in the two bottom rows of Figure 5 .
In SPH, most smoothing kernels W incorporate the smoothing length as a scaling parameter, such that they can be expressed in the form
Further, it is common to choose radial kernels such that they can be written in the form
Here, ρ is a normalization factor to satisfy the continuous unity property. Probably the most common kernel is of Gaussian type:
However, as the Gaussian does not have a compact support, it is truncated at η = 1. Thus, the parameter should be chosen in a way that values of the resulting kernel W (η) for η > 1 can be safely neglected. For the value of the Gaussian at η = 1, we obtain
which inspires the condition
Together with the continuous unity property, this allows us to determine both parameters ρ and , such that we use the Gaussian in the form
Here, we have expressed ρ as a function of . An alternative to the Gaussian are Matérn kernels [29] . Contrary to the Gaussian which is arbitrarily often continuously differentiable, Matérn kernels differ in smoothness. The C 0 -Matérn kernel, which is not differentiable but just continuous at η = 0, is given by
for which we chose = 6.52. A higher regularity at η = 0 can be achieved with the C 2 -Matérn kernel
which we use in our experiments with = 8.04. Although the truncated Gaussian is a common choice, the original SPH paper [47] already introduced a kernel with a compact support, namely
which we will call Lucy kernel. Other commonly used kernels with compact support are the cubic spline [45] (2.27)
and the Wendland C 4 kernel [76] (2.28)
3. SPH Inpainting. As inpainting problem, we consider the reconstruction of an image f from a sparse set of values at scattered pixel locations. These locations, called mask points, take the role of particles p j for our SPH-inspired inpainting procedure. The set of all mask points is the inpainting mask c.
Lengths. In order to use (2.6) for inpainting, whether with the original particle approximation, Shepard interpolation, or the first order consistent method, we still have to determine an area of influence V j for each given mask point p j . Further, we want to enhance the adaptivity of the method by allowing for different smoothing lengths h j of the kernels centered at the individual particles. This adaptivity is motivated by the results in [24] .
A reasonable idea is to assume that the area of influence of a given mask point p j is the set of all points which are closer to p j than to any other mask point in c. This idea leads to a Voronoi tessellation of the domain Ω with seeds given by the mask points. Voronoi cells have been used before in the context of SPH [33, 68] with promising results.
As we are working in a discrete setting where the smallest unit of area is a pixel, a method which determines approximate Voronoi diagrams based on the squared Euclidean distance transform is our tool of choice for this task. This is a rather natural approach since the Voronoi cell Ω j associated to the mask point p j is defined as
Given a binary image g which only takes the values 0 and ∞ throughout a domain Ω, the distance transforms assigns to each pixel q in Ω its squared distance to the nearest pixel p with g(p) = 0. In our case, g takes the value 0 at the mask points and ∞ everywhere else. For practical applications, ∞ can be replaced by a sufficiently large number. In a two-dimensional domain, the squared distance transform is given by
such that it can be computed by two consecutive squared distance transforms in one dimension. We used the algorithm from [31, 73] to compute the distance transform, which is shown to have a complexity of O(n x n y ), where n x and n y are the number of pixels in the image domain Ω in x-and y-direction, respectively. For visualization purposes each and every Voronoi cell is depicted with a different color per Voronoi cell in Figure 2c .
After the Voronoi tessellation, each mask point p j is assigned as area of influence the area of its corresponding Voronoi cell. In our discrete setting, this area is defined as the sum of pixels that belong to that cell, i.e., pixels are assumed to have area equal to 1.
It seems natural to determine the smoothing length h j in relation to the volume V j , e.g. half the diameter of the Voronoi cell associated to p j . However, this particular choice is prone to result in pixels for which the requirements of the minimal necessary number of nearest neighbors are not satisfiable. If a pixel does not lie within the support of any kernel, it cannot be inpainted. A straightforward remedy would be to multiply the diameter of each Voronoi cell with a constant factor chosen such that each pixel has at least the desired minimum number of nearest neighbors. Unfortunately, this would result in oversmoothing and blurring as the resulting kernel supports would be rather large. Instead, we follow the adaptive, iterative approach of [24] for the choice of smoothing lengths, which also enforces that any pixel is inpainted with at least a specified minimum number of neighbors. The scheme starts by assigning each mask point p j an initial smoothing length h j,init . Using the corresponding kernels, unknown pixels are inpainted, but only if they lie within the support of a least a fixed number of kernels. All pixels which do not satisfy this requirement are not inpainted. Afterwards we check whether there are still pixels with no assigned value left. If so, we increase all smoothing lengths according to a certain rule and try again to inpaint those pixels which are not yet assigned a value. This procedure is repeated iteratively until each pixel is inpainted. As growing strategy for the smoothing lengths, we increase them linearly with the number of iterations.
The original method in [24] assigned initial smoothing lengths which are connected to the choice of V j as made there. However, our experiments showed that it is beneficial if each kernel starts with a minimal smoothing length of h j,init = 1. Thus, mask points can initially only be recognized as neighbors within a 3 × 3-patch around them such that the process starts with the smallest sensible isotropic support for each kernel. The smoothing length h j is in each step equal to the number of iterations.
The last parameter that we have to set is the required minimal number of nearest neighbors. For Shepard interpolation, we need at least one neighbor to perform an inpainting, whereas for the method with first order consistency, any pixel must be contained in the support of at least three kernels. For all methods, it is reasonable to choose a slightly larger necessary minimal number of nearest neighbors as this improves results. In particular for the method of first order consistency that reduces the chance to encounter cases in which all mask points closest to an unknown pixel are collinear. For our experiments we require a minimal number of five nearest neighbors.
3.2. Results for Random Inpainting Masks. We follow here a didactic approach and consider the test image "trui" of size 256 × 256 pixels as an example. Results and comparisons for further images can be found in subsection 4.5 and in the supplementary material. To get a first impression, we equip the test image with a mask c for which we randomly selected 5% of the pixels. Both, test image and mask, are shown in Figure 3 . We notice that the performance gains with first order consistency compared to zero order consistency are not significant whereas for kernels with compact support (Lucy, cubic spline, and C 4 -Wendland) the results are even worse with the first order consistency method. This seemingly paradox result is exactly the consequence of having a distribution of mask points such that (2.17) cannot be satisfied in a way that gives a convex combination on the righthand side for every unknown pixel q. In other words, the modified kernels at such a pixel do no longer obey the positivity requirements, resulting in over-and undershoots which can become quite severe. Moreover, the 3 × 3-system that needs to be solved at each pixel for each modified kernel may become almost singular, leading to further instabilities. To achieve better results we introduce and analyze a new approach in the following section.
Method. The goal of mixed order consistency is to combine zero order and first order consistency to get the best possible result. For this purpose, two inpaintings are done, one with zero order consistency and one with first order consistency. Both results are compared, and the one with the better reconstruction error is kept. The new For the same 5% random mask as before with a minimum of five neighbors we get the results depicted in Figure 6 , which are clearly superior to the results achievable with either the zero or first order consistent method. 
Input : Original image f , mask c Output : Reconstruction u Initialize: Perform Voronoi tessellation. Assign areas of influence V j to mask points.
Assign initial smoothing lengths h j = h j,init = 1 to mask points. Set k = 1. 1 while not all pixels q have been inpainted do 2 for each pixel q do 3 Detect neighboring mask points p j , j ∈ N (q) of q. 4 if number of neighbors is larger than or equal to required minimum and neighbors are not collinear then 5 Inpaint q with zero order consistency. 6 Inpaint q with zero order consistency. 7 if error of zero order consistency is less than error of first order consistency then 
Our spatial optimization relies on a novel densification strategy. Instead of a probabilistic approach as in [37] , we base our method on Voronoi tessellation. The algorithm starts with an empty mask and, as an initial step, inserts the minimum number of neighbors required at random positions. After this step, an initial inpainting takes place and a Voronoi tessellation is performed with the initial mask points as "seeds". Once this is done, we detect the Voronoi cell with the highest error and insert a new mask point at the pixel with the highest error within the cell. The error of the reconstruction at a pixel q j,k in the Voronoi cell Ω j is defined as
with f being the original image and u the reconstruction, whereas the error for the Voronoi cell Ω j is given by the sum of the reconstruction errors at all pixels in the cell, i,e.,
A new inpainting as well as a new Voronoi tessellation are then computed with the new mask and the process continues in the same manner until the required mask density is achieved. The densification algorithm is described in Algorithm 2. We remark that it is possible to insert more than one new mask point in each step to speed up the procedure. However, inserting too many mask points at once deteriorates the quality of the final mask. For the sake of completeness, we also mention that a densification approach using the L 1 -error within Voronoi cells has been used in [65] for nearest-neighbor and piecewise constant interpolation.
In our experiments, using the L 1 -error always yielded inferior results.
Input : Original image f , minimum number of neighbors, number of mask points to add per iteration, required density Output : Mask c, reconstruction u Initialize: Insert minimum number of neighbors at random positions. Perform Voronoi tessellation and initial inpainting. 1 while mask density < required mask density do 2 Find Voronoi cell(s) Ω j with highest error E Ω j . 3 Find pixel(s) q j,k in cell(s) Ω j with highest error(s) E q j,k . 
Apart from spatial optimization, we also incorporate a gray value optimization of the mask points for a fixed mask c. The goal of this process is to find the optimal gray values g such that the mean square error of the reconstructed image is minimal. Indeed, for a fixed mask c, the inpainting is given by (2.6). However, due to the adaptive nature of the smoothing length h j which is not only determined in dependence of the mask point p j , but also in dependence of the pixel q currently under consideration for inpainting, we have to perform one inpainting for our final mask first to determine all necessary smoothing lengths. Once this is done, (2.6) can be written as a matrix vector multiplication of the form u = A f where u is a vector containing the values of the reconstruction at every pixel, A is a matrix containing the values of all modified kernels at all pixels multiplied with their area of influence V j , and f is a vector containing the values of the original image at all mask points p j . Thus, as A is fixed now, u can be interpreted as the solution to an interpolation problem at the mask points. With this interpretation, it is straightforward to consider the corresponding least-squares problem. With the above form for u, it can be written as
where f denotes the vector with the original image values at all pixels and g is a vector with gray values at the mask points, which can be determined by solving the normal equations
Although the least squares problem can be solved directly, we prefer an iterative solver instead, specifically the conjugate gradients (CG) method. We always initialize the iterates g k with the zero vector. As stopping criterion we use a threshold on the relative residual defined such that
The above procedure is clear for the zero and first order consistency method as they use the same kind of modified kernel in every pixel. It stays valid for the mixed order consistency method as in this approach, the kernel that is used at each pixel q is of the same type for all mask points p j contributing to the inpainting at q. Thus, for the mixed consistency method, the type of the modified kernel changes with the rows in A, but stays the same within each row over all columns. We can still write the whole inpainting process as a matrix-vectormultiplication and solve the associated least-square problem to perform tonal optimization.
Optimized Data with Isotropic Kernels. As an example for inpainting on an optimized mask with zero order consistency, we present the results produced with a Gaussian kernel in Figure 7 . Even without tonal optimization, the MSE improves by roughly a factor 6.5. With tonal optimization, the MSE improves by a factor 10 with respect to the random mask and by almost 35% with respect to the result on the spatially optimized mask without tonal optimization. As far as spatial optimization is concerned, the densification process prefers to capture the geometry of the image, by adding more mask points near edges compared to rather homogeneous regions of the image. Changing the SPH inpainting method from zero order consistency to the mixed order consistency method improves the result even further as can be seen in Figure 8 . Using the same isotropic Gaussian kernel as before, spatial optimization improves the MSE by almost a factor 8 compared to the the random mask result in Figure 6a . A comparison with respect to the inpainting method instead of with respect to the mask shows improvement of almost a factor 2 with respect to the MSE compared to the results in Figure 7 . Unfortunately, there seems to be no structure in the distribution of pixels for which a first order consistency and for which a zero order consistency method performs better, respectively, as can be seen in Figure 8b . Table 1 summarizes MSEs of inpainting results for "trui" with the other isotropic kernels used in Figures 4 and 5 if these kernels are equipped with optimized masks containing 5% of all pixels and tonal optimization is performed. Once again, the benefits of mixed order Table 1 : MSE comparison between zero and mixed order consistency optimized inpainting with isotropic kernels on "trui" for 5% masks.
consistency are quite substantial since a significant decrease of the MSE has been achieved in all cases compared to zero order consistency. Even the best performing kernel for the zero order consistency method, the C 2 -Matérn kernel, has an MSE which is approximately 50% larger than the MSE of the worst performing kernel in the mixed consistency setting and almost double of the MSE of the best performing kernel in the mixed consistency setting which is the C 2 -Matérn kernel as well. Further, it appears that compactly supported kernels perform worse in this setting for the zero order consistency SPH inpainting than truncated kernels, but are competitive in the mixed order consistency method.
Inpainting with Anisotropic Kernels. The observation that optimized mask points tend to cluster around edges and the fact that edges are clearly oriented structures suggest to adapt the support of kernels to account for this by incorporating anisotropy. This is further supported by results which show that incorporating anisotropy in other inpainting strategies can improve reconstruction quality compared to the related isotropic method [63] .
For SPH, anisotropic kernels have been used in the so-called Adaptive Smoothed Particle Hydrodynamics (ASPH) formulation [66] to better account for the actual distribution of particles. Here, we replace the smoothing length h by a symmetric positive definite tensor G ∈ R 2×2 for a two-dimensional problem and redefine η as
G has units of inverse length and in the isotropic case it is given by a diagonal matrix with each diagonal element equal to 1 h . This observation makes it clear that we also have to adapt the normalization of our kernels from a factor ρ h 2 in (2.19) to a factor ρ det(G). For SPH inpainting, we determine the anisotropy from the distribution of mask points. For this purpose, we follow the approach of [77] by constructing a weighted local covariance matrix C within a fixed predetermined window around each mask point. For a known mask point p j , the covariance matrix is given by
Here, numbers the mask points within a neighborhood of p j . It is necessary to restrict the set of mask points under consideration to such a neighborhood to catch the locally prevalent direction of structures in the image. Next, we perform a singular value decomposition (SVD). As C j is symmetric and positive semidefinite by construction, this is the same as the eigenvalue decomposition
with a rotation matrix Q and a matrix D with nonnegative eigenvalues along the diagonal in decreasing order. As C j is constructed from the positions of mask points, its eigenvalues can be assigned a unit of length. The eigenvectors in Q correspond to the directions of major and minor axis of an ellipse whose orientation is in line with the locally prevalent orientation in the distribution of mask points. Hence, the tensor G is given by
such that it has units of inverse length as desired.
In our experiments, we incorporate anisotropy after spatially optimizing mask points for isotropic kernels. We fix the window size for construction of covariance matrices to 25 × 25 pixels and demand a minimum number of 15 mask points within that window. If this minimum number of mask points is not satisfied, the corresponding kernel stays isotropic. This behavior is desirable since the densification process results in masks where the majority of mask points are placed near discontinuities rather than in homogeneous areas of the image. Thus, a low local density of mask points implies homogeneous areas of the image. The results achieved with mixed order consistency and an anisotropic Gaussian kernel are depicted in Figure 9 . Compared to the result in Figure 8 , we observe an improvement in MSE of only 1.18, which is roughly 10%, whereas a large amount of mask points is now equipped with anisotropic kernels. This relatively moderate improvement may be explained by the fact that our method of determining anisotropy relies on the local spatial distribution of mask points whereas many important structures in the given image live on a mesoscale. This behavior cannot be captured by increasing the size of the search window as the covariance matrix becomes prone to incorporating the orientations of neighboring structures, resulting in a more isotropic behavior instead of a better orientation along the mesoscale structures. Table 2 summarizes MSEs of inpainting results for "trui" with the other anisotropic kernels used if these kernels are equipped with optimized masks containing 5% of all pixels and tonal optimization is performed.
Kernel MSE Gaussian 10.41 C 0 -Matérn 9.49 C 2 -Matérn 9.73 Lucy 11.63 cubic spline 11.28 C 4 -Wendland 11.67 Table 2 : MSE with mixed order consistency optimized inpainting with anisotropic kernels on "trui" for 5% masks.
Methods. In order to assess the performance of our inpainting method we compare it with other existing methods in the literature. As a first example, we consider results for homogeneous diffusion inpainting on "trui" for an optimized 5% mask from [50] . There, spatial optimization was done with a probabilistic sparsification and further improved with a Nonlocal Pixel Exchange (NLPE). With optimally chosen mask points and gray values, homogeneous diffusion inpainting shows an impressive quality in reconstructing the original image. We compare these results with the ones we got for a mixed order consistency SPH inpainting with C 2 -Matérn kernels. To be fair, we only consider isotropic kernels since homogeneous diffusion has no way to incorporate anisotropy. The results are summarized in Table 3 [50] and results from our method with an isotropic C 2 -Matérn kernel with mixed order consistency.
form homogeneous diffusion, both without and with tonal optimization. We also compared to further results in [50] and always achieved MSEs which were 20% to 40% smaller than the ones reported for homogeneous diffusion inpainting.
To evaluate SPH inpainting with anisotropic kernels, we consider the results from [35] achieved with edge-enhancing diffusion (EED) inpainting for "trui" with a mask of density 4% that is constructed by probabilistic sparsification. The authors report the MSE of inpainting on this mask without tonal optimization and improve the location of mask points further with NLPE before considering tonal optimization. As competitor, we used a mixed order consistency SPH inpainting with anisotropic C 0 -Matérn kernels on an optimized 4% mask. We also tried to improve our inpaintings with NLPE. However, the reduction on MSE was negligible. Results for both inpainting methods are summarized in Table 4 Table 4 : MSE for inpaintings of "trui" on optimized 4% masks. Compared are results from EED inpainting [35] and our method with an anisotropic C 0 -Matérn kernel with mixed order consistency.
outperform EED if we only incorporate spatial optimization, but no tonal optimization. By construction, EED should perform better in preserving edges [74] . Thus, we conjecture that probabilistic sparsification, which only relies on pointwise errors, is inferior to our Voronoibased densification method as long as the former is not improved by a consecutive NLPE.
To evaluate the performance of our method for images rich in texture, we consider the exemplar-based inpainting technique from [41] and the results given there for an inpainting of a gray value version of the "baboon" image. In this setting, the authors report an MSE of 518.52 on a mask constructed with "densification by dithering" and a consecutive NLPE. As tonal optimization is not considered in [41] , we compare to the result our method could achieve for mixed order consistency inpainting with isotropic Gaussians on an optimized 10% mask without tonal optimization in Figure 10 . Here, the MSE is 290.64, which means we outperform the exemplar-based inpainting method by almost 44%.
(a) Image "baboon" (b) Optimized 10% mask (c) MSE = 290.64 Figure 10 : Inpainting of "baboon" (rescaled to 256 × 256 pixels as in [41] ) with a spatially optimized mask for a mixed order consistency method with an isotropic Gaussian kernel. From left to right: Original image (a), optimized 10% mask for mixed order consistency method (b), and mixed order consistency inpainting result with given mask without tonal optimization (c).
We have shown that smoothed particle hydrodynamics is a highly competitive method for the challenging problem of sparse data inpainting. It can produce results on par or even better than other, better explored PDE-or exemplar-based inpainting strategies.
The success of SPH for sparse inpainting relies on several novel modifications. With regard to the interpolation procedure, we presented a way to combine the strength of first and zero order methods into a mixed order method. Moreover, we presented a better approach to choose method parameters based on Voronoi tesselations.
The main ingredient to reveal the potential of SPH is the use of optimally chosen data. We proposed a new densification process based on Voronoi tesselations, which lead naturally to a strategy based on a regional error instead of a completely local pointwise error. Thus, a larger amount of data is considered in the optimization procedure, yielding better suitable inpainting masks. Furthermore, we introduced a so far unexplored formulation that allows to use SPH for least-square approximation, i.e., the optimization of data not only in the spatial, but also in the tonal domain.
We hope that our work will help to give SPH-based inpainting the attention that it deserves. Moreover, we believe that some of our novel concepts, e.g. the Voronoi-based densification for data optimization, will also be useful in applications beyond SPH-based inpainting. Figure SM5 : Inpainting of "baboon" with a 10% spatially and tonally optimized mask with a mixed order consistency method and anisotropic kernels.