Sparse Non-Negative Stencils for Anisotropic Diffusion ^††thanks: This work was partly supported by ANR grant MESANGE ANR-08-BLAN-0198.

Our numerical scheme, Anisotropic Diffusion using Lattice Basis Reduction (AD-LBR), involves the construction of stencils whose geometry is tailored after the local diffusion tensor. Its essential feature is non-negativity: the discrete energy ${\mathcal{E}}_{h}(u)$ is written as a sum (3) of squared differences of values of $u$, with non-negative weights $\gamma_{z}(e)\geq 0$. This discretization is consistent if for each $z\in\Omega_{h}$, and any smooth $u$,

$$h^{d}\|\nabla u(z)\|_{\operatorname{\mathbf{D}}(z)}^{2}=h^{d-2}\sum_{e\in V(z)}\gamma_{z}(e)\,\langle\nabla u(z),he\rangle^{2}.$$

(4)

Indeed, the left hand side approximates the contribution of the “voxel” $z+[-h/2,h/2]^{d}$ to the integral (1), while the right hand side is obtained by inserting the first order approximation $u(z+he)\simeq u(z)+\langle\nabla u(z),he\rangle$ in (3). The identity (4) is in turn equivalent to

$$\operatorname{\mathbf{D}}(z)=\sum_{e\in V(z)}\gamma_{z}(e)\,ee^{\mathrm{T}}.$$

(5)

The next lemma shows how to obtain such a decomposition in 2D. We denote by $u^{\perp}:=(-b,a)$ the rotation of a vector $u=(a,b)\in{\mathbb{R}}^{2}$ by $\pi/2$, in such way that for all $v\in{\mathbb{R}}^{2}$:

$$\langle u^{\perp},v\rangle=\det(u,v).$$

The diffusion tensor $\operatorname{\mathbf{D}}$ is meant to measure gradients, as in (1). In order to measure angles between vectors, we introduce a Riemannian metric ³³3The Laplace Beltrami operator associated to $\operatorname{\mathbf{M}}$ does not coincide with $\operatorname{div}(\operatorname{\mathbf{D}}\nabla\cdot)$, unless $\operatorname{\mathbf{D}}$ is identically of determinant $1$. This is not an issue for our application. $\operatorname{\mathbf{M}}$ on the domain $\Omega$, which is proportional to the inverse of $\operatorname{\mathbf{D}}$: for all $z\in\Omega$

$$\operatorname{\mathbf{M}}(z):=\operatorname{\mathbf{d}}(z)\operatorname{\mathbf{D}}(z)^{-1},\,\text{ where}\,\operatorname{\mathbf{d}}(z):=\det(\operatorname{\mathbf{D}}(z))^{\frac{1}{d}}.$$

(7)

The normalizing factor $\operatorname{\mathbf{d}}(z)$ was chosen so as to normalize the metric determinant: $\det(\operatorname{\mathbf{M}}(z))=1$. This normalization reflects the fact that the construction of our stencil $V(z)$ depends on the preferred direction of diffusion, and on the amount of anisotropy, whereas the absolute strength of diffusion is irrelevant. In dimension $d=2$, one easily checks that for any $z\in\Omega$ and any $e,f\in{\mathbb{R}}^{2}$, one has

$$\langle e^{\perp},\operatorname{\mathbf{D}}(z)f^{\perp}\rangle=\operatorname{\mathbf{d}}(z)\langle e,\operatorname{\mathbf{M}}(z)f\rangle.$$

(8)

The AD-LBR is based on decompositions (5), given by the previous lemma, with a family of vectors $(e_{i})_{i=0}^{2}$ chosen so that the scalar products appearing in (6) are non-positive. The adequate concept is that of $M$-obtuse superbase of ${\mathbb{Z}}^{d}$ [4].

In dimension $d\leq 3$, there exists for each $M\in S_{d}^{+}$ at least one $M$-obtuse superbase of ${\mathbb{Z}}^{d}$ [4]. The practical construction of such superbases is discussed in §2, and based on lattice basis reduction algorithms described in [11, 23, 17] (hence the name of our numerical scheme). This construction has a logarithmic numerical cost ${\mathcal{O}}(\ln\kappa(M))$ in the anisotropy ratio of the matrix $M$:

$$\kappa(M):=\max_{|u|=|v|=1}\frac{\|u\|_{M}}{\|v\|_{M}}=\sqrt{\|M\|\|M^{-1}\|}.$$

(9)

The AD-LBR energy ${\mathcal{E}}_{h}:L^{2}(\Omega_{h})\to{\mathbb{R}}_{+}$, see (3), is in two dimensions written in terms of the following stencils and coefficients. Let $z\in\Omega$, and let $e_{0},e_{1},e_{2}$ be an $\operatorname{\mathbf{M}}(z)$-obtuse superbase of ${\mathbb{Z}}^{2}$. We set

$$V(z):=\{e_{0},e_{1},e_{2},\,-e_{0},-e_{1},-e_{2}\},\\$$

(10)

and for $0\leq i\leq 2$, with the convention $e_{i+3}:=e_{i}$,

$$\gamma_{z}(\pm e_{i}):=-\frac{1}{2}\langle e^{\perp}_{i+1},\operatorname{\mathbf{D}}(z)e^{\perp}_{i+2}\rangle.$$

(11)

Lemma 1 implies the announced decomposition (5), and the weights $\gamma_{z}$ are non-negative in view of (8). These weights $\gamma_{z}:{\mathbb{Z}}^{2}\to{\mathbb{R}}_{+}$, extended by $0$ outside $V(z)$, do not depend on the choice of $\operatorname{\mathbf{M}}(z)$-obtuse superbase $(e_{0},e_{1},e_{2})$, see Lemma 11. Stencils of the three dimensional AD-LBR are described in §2, and involve a construction of Selling⁴⁴4The authors would like to thank Professor P. Q. Nguyen for pointing out this $12$ points 3D stencil, which is simpler and sparser than the $14$ points stencil proposed by the authors in an earlier version of the manuscript. [22]. The above description of the stencils $V(z)$ is suitable for periodic, reflected, and Dirichlet boundary conditions (extending $u$ by zero outside $\Omega_{h}$ in the latter case). In the case of Neumann boundary conditions, a slight modification is in order:

$$V(z;\,h):=\{e\in V(z);\,z+he\in\Omega_{h}\}.$$

We have so far established three strongpoints of the AD-LBR:

Non-negativity.

Off diagonal coefficients of the symmetric semi-definite $N\times N$ matrix, $N=\#(\Omega_{h})$, associated to the energy ${\mathcal{E}}_{h}$ are non-positive, while diagonal coefficients are positive.

Sparsity.

Stencil cardinality is uniformly bounded, without restriction on the anisotropy ratio $\kappa(\operatorname{\mathbf{D}}(z))$ of the diffusion tensor.

Complexity.

The construction of the stencil $V(z)$, and of the associated coefficients $\gamma_{z}$, has a logarithmic cost ${\mathcal{O}}(\ln\kappa(\operatorname{\mathbf{D}}(z)))$ in the anisotropy ratio of the diffusion tensor.

The next result, Theorem 1, restricted to the two dimensional case, establishes that AD-LBR is asymptotically equivalent to a more involved and computationally intensive procedure: a finite element discretization of the energy (1), on an Anisotropic Delaunay Triangulation (ADT, see [10] and below) of the domain $\Omega$. Under the assumptions of Theorem 1, AD-LBR benefits from two additional guarantees, that we state informally and without proof.

No chessboard artifacts.

Some numerical schemes for anisotropic diffusion suffer from chessboard artifacts, in the sense that periodic artifacts develop at the pixel level. Such artifacts cannot develop in finite element discretizations, since they would lead to high frequency oscillations of the finite element interpolant, and therefore to an increase of the energy (12). The asymptotic equivalence of the AD-LBR with a finite element discretization also rules out these defects.

Spectral correctness.

The $n$-th smallest eigenvalue $\lambda_{n}(h)$ of the symmetric matrix associated to $h^{-d}{\mathcal{E}}_{h}$ (3), converges as $h\to 0$ towards the $n$-th smallest eigenvalue $\lambda_{n}$ of the continuous operator $-\operatorname{div}(\operatorname{\mathbf{D}}\nabla)$, for any given integer $n\geq 0$. See Figure 8 page 8 for an illustration. This follows from a similar property of the finite element energy ${\mathcal{E}}^{\prime}_{h}$ (12), and from the asymptotic equivalence (13).

Our convergence result, Theorem 1 below, is specialized to the case of a square periodic domain, which covers reflecting boundary conditions frequently used in image processing. Since the grid discretization must be compatible with the boundary conditions, any scale parameter $h$ appearing in the rest of the paper is assumed to be the inverse of a positive integer:

$$h\in\{1/n;\,n\geq 1\}.$$

Let us mention that the finite element discretization on an ADT is a more general procedure than AD-LBR, since it does not require the domain $\Omega$ to be sampled on a grid. This flexibility can be used to locally increase the density of vertices, in places where solution $u$ is expected to be less regular, or to insert vertices exactly on $\partial\Omega$ for a better discretization of boundary conditions. (Such refinements are however generally incompatible with image processing since the unknowns, the pixel values, lie by construction on a fixed and given cartesian grid.) Here and as often, the performance of AD-LBR is at the cost of its specialization.

The proof of Theorem 1 is postponed to §3, but for the sake of concreteness, we describe here the concept of Anisotropic Delaunay Triangulation (ADT) [10]. In the rest of this introduction, and in §3, we assume as in Theorem 1 that the diffusion tensor $\operatorname{\mathbf{D}}$ is defined on the square $[0,1]^{2}$ and satisfies periodic boundary conditions. We extend it, as well as the metric $\operatorname{\mathbf{M}}$, to the whole plane ${\mathbb{R}}^{2}$ by periodicity.

We specialize the concept of ADT [10], to the domain ${\mathbb{R}}^{2}$ and the collection of vertices $h{\mathbb{Z}}^{2}$. For that purpose, we introduce some notations. For all $p,q\in{\mathbb{R}}^{2}$, we denote by $\delta_{p}(q)$ the distance from $p$ to $q$, as measured by the metric at the point $p$:

$$\delta_{p}(q):=\|q-p\|_{\operatorname{\mathbf{M}}(p)}.$$

(14)

We denote by $\Delta_{h}(q)$ the least distance from a point $q\in{\mathbb{R}}^{2}$, to the grid $h{\mathbb{Z}}^{2}$:

$$\Delta_{h}(q):=\min_{p\in h{\mathbb{Z}}^{2}}\delta_{p}(q).$$

(15)

We introduce the Voronoi cell $\operatorname{Vor}_{h}(p)$ of a grid point $p\in h{\mathbb{Z}}^{2}$, which is the collection of points $q\in{\mathbb{R}}^{2}$ closer to $p$ than to any other grid point:

$$\operatorname{Vor}_{h}(p):=\{q\in{\mathbb{R}}^{2};\,\delta_{p}(q)=\Delta_{h}(q)\}.$$

(16)

The collection of Voronoi cells is referred to as the Voronoi diagram, see Figure 2. A Voronoi vertex is a point $q\in{\mathbb{R}}^{2}$ at which at least three distinct Voronoi regions intersect: $(\operatorname{Vor}_{h}(p_{i}))_{i=1}^{k}$, $k\geq 3$, $p_{i}\in h{\mathbb{Z}}^{2}$. We attach to $q$ a dual Voronoi cell $T_{q}$, defined as the convex hull of the points $(p_{i})_{i=1}^{k}$, see Figure 3.

The geometric dual ${\mathcal{Q}}_{h}$, of the Voronoi diagram, is defined as the collection of all dual Voronoi cells $T_{q}$. Note that, generically on the metric $\operatorname{\mathbf{M}}$, no more than three Voronoi regions can intersect at any point in ${\mathbb{R}}^{2}$, thus the elements of ${\mathcal{Q}}_{h}$ are generically triangles. If $h$ is small enough, we show in §3 (using the Dual Triangulation Theorem in [10]) that $T_{q}$ is a strictly convex polygon, of vertices $(p_{i})_{i=1}^{k}$ with the above notations, and that ${\mathcal{Q}}_{h}$ is a polygonization (generically a triangulation) of ${\mathbb{R}}^{2}$, with vertices $h{\mathbb{Z}}^{2}$.

Since the metric $\operatorname{\mathbf{M}}$ and the vertices $h{\mathbb{Z}}^{2}$ are periodic (recall that $h=1/n$ for some integer $n\geq 1$), arbitrarily triangulating the elements of ${\mathcal{Q}}_{h}$, respecting periodicity, yields a periodic triangulation ${\mathcal{T}}_{h}$.

We establish in §3.1 the existence of the ADT ${\mathcal{T}}_{h}$. Incidentally, we show in Lemma 7 (iii) page 7 that the angles of the elements of ${\mathcal{T}}_{h}$, measured with respect to the local metric $\operatorname{\mathbf{M}}$, are asymptotically acute. This geometrical property (which holds thanks to our special choice of triangulation vertices, on a grid) is linked to the non-negativity of AD-LBR: indeed, it is known that finite elements discretizations such as (12) yield non-negative numerical schemes, and the discrete maximum principle, if the mesh satisfies a non-obtuse angle condition, see Lemma 3.1 in [9].

Subsection §3.2 is devoted to the study of $M$-obtuse superbases of ${\mathbb{Z}}^{2}$, and their cousins $M$-reduced bases of ${\mathbb{Z}}^{2}$, on which the AD-LBR relies: we discuss their characterization, uniqueness and stability properties. We study in §3.3 the finite element stencils, defined for $p\in h{\mathbb{Z}}^{2}$ by

$$V_{h}(p):=\{e\in{\mathbb{Z}}^{2};[p,p+he]\text{ is an edge of }{\mathcal{T}}_{h}\},$$

(17)

see Figure 3 (right). We show that $V_{h}(p)$ coincides with the AD-LBR stencil $V(p)$, unless the lattice ${\mathbb{Z}}^{2}$ admits a basis almost orthogonal with respect to the scalar product associated to $\operatorname{\mathbf{M}}(p)$, see Lemma 13. This is tied to the fact that orthogonal grids admit several (usual) Delaunay triangulations. Overcoming this technical difficulty, we conclude the proof of Theorem 1.

Note that the construction of the ADT ${\mathcal{T}}_{h}$ is not easy to parallelize, in particular when anisotropy is pronounced since the Voronoi regions of far away points interact. The construction of ${\mathcal{T}}_{h}$ also involves solving polynomial equations of degree four, because Voronoi regions boundaries are conics, and Voronoi vertices must be identified at their intersections. In contrast, the AD-LBR stencils are independent of each other, and the numerical cost of their construction only grows logarithmically with the metric anisotropy.

Algorithms for the construction of privileged bases of lattices, consisting of short and almost orthogonal vectors, have attracted an important research effort from the mathematical community, over a long period of time. The first such algorithm dates back to Lagrange [11], and is restricted to two dimensional lattices. Methods for high dimensional lattices, such as the LLL algorithm [12], are of key importance for integer programming and cryptography [18]. AD-LBR is based on the original algorithm of Lagrange [11], and on its recent extension to three dimensional lattices [23, 17]. These methods output a basis of ${\mathbb{Z}}^{d}$ reduced in the sense of Minkowski, which in dimension $d\leq 4$ is equivalent to the following definition.

We denote by $e_{1}{\mathbb{Z}}+\cdots+e_{k}{\mathbb{Z}}$ the sub-lattice of ${\mathbb{Z}}^{d}$ generated by vectors $e_{1},\cdots,e_{k}\in{\mathbb{Z}}^{d}$. This sub-lattice equals $\{0\}$ by convention if $k=0$.

For each $d\leq 4$, and each $M\in S_{d}^{+}$, there exists at least one $M$-reduced basis [17]. In contrast, there exists $M\in S_{5}^{+}$ for which no basis of ${\mathbb{Z}}^{d}$ satisfies (18). The norms of the elements $(e_{i})_{i=1}^{d}$ of an $M$-reduced basis,

$$\lambda_{i}(M):=\|e_{i}\|_{M},$$

(19)

are called the Minkowski minima, and are independent of the choice of $M$-reduced basis. In particular, $e_{1}$ is the shortest vector of ${\mathbb{Z}}^{d}$, with respect to the norm $\|\cdot\|_{M}$, and $e_{2}$ is the shortest linearly independent vector.

We emphasize that obtaining an $M$-reduced basis, i.e. solving the minimization problems (18), is both simple and cheap numerically. In dimension $d=2$, this is the object of Lagrange’s algorithm [11] (later rediscovered by Gauss and often erroneously called Gauss’s algorithm, see [17]): initialize $(e,f)$ as the canonical basis of ${\mathbb{Z}}^{2}$, and

		$\displaystyle\text{\bf Do }(e,\ f):=(f,\ e-\mathrm{Round}(\langle e,Mf\rangle/\|f\|_{M}^{2})\,f),$		(20)
		$\displaystyle\text{\bf while }\|e\|_{M}>\|f\|_{M}.$

This algorithm can be regarded as a two dimensional geometrical generalization of greatest common divisor computation. It can be extended to higher dimension and, in dimension up to four, outputs an $M$-reduced basis after at most ${\mathcal{O}}(\ln\kappa(M))$ iterations [17], each consisting of ${\mathcal{O}}(1)$ operations among reals.

The elements of an $M$-reduced basis are heuristically never very far from being orthogonal, as illustrated by the following lemma.

The practical construction of the two dimensional AD-LBR stencil at a point $z\in\Omega$ amounts to (i) compute an $\operatorname{\mathbf{M}}(z)$-reduced basis $(e,f)$ using Lagrange’s algorithm (20), (ii) replace $f$ with $-f$, if necessary, so that $\langle e,Mf\rangle\leq 0$, and (iii) define the stencil $V(z)$ and the weights $\gamma_{z}$ in terms of the $M$-obtuse superbase $(e,f,g)$ of ${\mathbb{Z}}^{2}$, where $g=-e-f$, as described in (10) and (11).

The rest of this section is devoted to the description of the three dimensional AD-LBR stencils. In constrast with the two dimensional case, the construction of the 3D stencil $V(z)$ at a point $z\in\Omega$, involves a $\operatorname{\mathbf{D}}(z)$-obtuse basis, instead of an $\operatorname{\mathbf{M}}(z)$-obtuse basis.

In view of the previous Proposition, obtaining a $D$-obtuse superbase of ${\mathbb{Z}}^{3}$ has numerical cost ${\mathcal{O}}(\ln\kappa(D))$. Indeed a $D$-reduced basis needs to be computed in a preliminary step, after what Proposition 1 is applied for a negligible ${\mathcal{O}}(1)$ cost. An alternative method for the construction of $D$-obtuse superbases of ${\mathbb{Z}}^{3}$ is presented in [4] and in appendix B of [3], but its numerical complexity is not known to the authors.

The three dimensional AD-LBR is defined by the following stencils and coefficients. Let $z\in\Omega$, let $D:=\operatorname{\mathbf{D}}(z)$, and let $(e_{i})_{i=0}^{3}$ be a $D$-obtuse superbase of ${\mathbb{Z}}^{3}$. We set

$$V(z):=\{e_{k}\times e_{l};\ k,l\in\{0,1,2,3\},\ k\neq l\},$$

and if $\{i,j,k,l\}=\{0,1,2,3\}$, $i\neq j$ and $k\neq l$, then

$$\gamma_{z}(e_{k}\times e_{l}):=-\frac{1}{2}\langle e_{i},De_{j}\rangle.$$

As announced, $\#(V(z))=12$, and the weights $\gamma_{z}$ are non-negative. The proof of the scheme consistency (5), due to Selling [22], is reproduced in the next lemma for completeness. A generalization, appearing in Appendix B of [3], allows in arbitrary dimension to build a non-negative decomposition of the form (25) from a $D$-obtuse superbase of ${\mathbb{Z}}^{d}$. However the non existence of such a superbase, for some matrices $D\in S_{4}^{+}$, forbids a straightforward extension of AD-LBR to higher dimension.

This section is devoted to the proof of Theorem 1: the asymptotic equivalence of AD-LBR with a finite element discretization on an Anisotropic Delaunay Triangulation (ADT). We use the notations of §1. The existence of the ADT ${\mathcal{T}}_{h}$ is established in the first subsection, for $h$ sufficiently small, as well as a few of its properties. The second subsection is devoted to the study of $M$-reduced bases. Theorem 1 is proved in the third subsection, by comparing the stencils of the AD-LBR and of the finite element discretization.

We denote by ${\boldsymbol{\kappa}}$ the maximum anisotropy ratio (9) of the diffusion tensor

$${\boldsymbol{\kappa}}:=\max_{z\in\Omega}\kappa(\operatorname{\mathbf{D}}(z)).$$

(26)

Observing that $\kappa(\operatorname{\mathbf{D}}(z))=\kappa(\operatorname{\mathbf{M}}(z))$, and recalling that $\det(\operatorname{\mathbf{M}}(z))=1$, one easily checks that

$${\boldsymbol{\kappa}}^{-\frac{1}{2}}\|e\|\leq\|e\|_{\operatorname{\mathbf{M}}(z)}\leq{\boldsymbol{\kappa}}^{\frac{1}{2}}\|e\|,$$

(27)

for all $z\in\Omega$ and all $e\in{\mathbb{R}}^{2}$.