Lattice-Based High-Dimensional Gaussian Filtering and the Permutohedral Lattice

Abstract

High-dimensional Gaussian filtering is a popular technique in image processing, geometry processing and computer graphics for smoothing data while preserving important features. For instance, the bilateral filter, cross bilateral filter and non-local means filter fall under the broad umbrella of high-dimensional Gaussian filters. Recent algorithmic advances therein have demonstrated that by relying on a sampled representation of the underlying space, one can obtain speed-ups of orders of magnitude over the naïve approach. The simplest such sampled representation is a lattice, and it has been used successfully in the bilateral grid and the permutohedral lattice algorithms. In this paper, we analyze these lattice-based algorithms, developing a general theory of lattice-based high-dimensional Gaussian filtering. We consider the set of criteria for an optimal lattice for filtering, as it offers a good tradeoff of quality for computational efficiency, and evaluate the existing lattices under the criteria. In particular, we give a rigorous exposition of the properties of the permutohedral lattice and argue that it is the optimal lattice for Gaussian filtering. Lastly, we explore further uses of the permutohedral-lattice-based Gaussian filtering framework, showing that it can be easily adapted to perform mean shift filtering and yield improvement over the traditional approach based on a Cartesian grid.

Appendix

1.1 A.1 Useful Lemmas

This section outlines a number of useful lemmas required for computing the SNR terms for K _x and Z _x. Most of them concern the permutohedral lattice—more machinery is required for the derivations on $A^{*}_{d}$. They are included here so that Sect. A.2 and A.3 can remain succinct and parallel.

Lemma A.1

The power of a d-dimensional Gaussian kernel with covariance θI is (4πθ)^−d/2.

Proof

as desired. □

Lemma A.2

Let $\mathcal{S}$ be a d-dimensional simplex defined by vertices $\mathbf {s_{0}},\ldots,\mathbf {s_{d}}$, and let B ₀(x),B ₁(x),…,B _d(x) be the barycentric coordinates of x for the vertices in the given order. Then, ∀k,j∈{0,…,d},

$$\frac{\int_{\mathcal{S}}B_{k}(\mathbf {x})B_j(\mathbf {x})d\mathbf {x}}{\int_{\mathcal{S}}d\mathbf {x}} =\begin{cases} \frac{2}{(d+1)(d+2)},&k=j,\\ \noalign {\vspace {3pt}} \frac{1}{(d+1)(d+2)},&k\ne j. \end{cases} $$

Proof

Let P be the face defined by all vertices except $\mathbf {s_{k}}$. Then $\mathcal{S}$ is a pyramid with base P and apex $\mathbf {s_{k}}$. Let h be the height of the altitude and let y be the foot of the altitude. Consider slicing $\mathcal{S}$ along a hyperplane parallel to P that contains $(1-t)\cdot \mathbf {s_{k}}+t\cdot \mathbf {y}$, for some t∈[0,1]. The slice will be similar to P with a scale factor of t, and by definition of barycentric coordinates, B _k(⋅) will evaluate to a constant, namely 1−t, across the whole slice. Therefore,

(A.1)

Note that the same integral without the weight term (1−t)² will simply yield the volume of $\mathcal{S}$:

(A.2)

Dividing (A.1) by (A.2) yields the first of the two cases. When k≠j, we consider applying the same decomposition to P: let P′ be the face of P spanned by all vertices of $\mathcal{S}$ except $\mathbf {s_{k}}$ and $\mathbf {s_{j}}$. Then P is again a pyramid with base P′ and height h′, and y′ the foot of the altitude. Each slice of P parallel to P′ cutting the altitude at $(1-r)\cdot \mathbf {s_{j}}+r\cdot \mathbf {y'}$ is again a copy of P′ scaled by r, and the barycentric weight B _j(⋅) is uniformly 1−r. Therefore,

(A.3)

Without the barycentric weights,

(A.4)

Dividing (A.3) by (A.4) yields the second of the two cases, as desired. □

Lemma A.3

Let $\mathbf {x},\mathbf {y}\in A^{*}_{d}$. Then the number of Delaunay cells containing both x and y is given by

$$\begin{cases} (d+1-k)!k!,&\parbox[t]{2.15in}{if $\mathbf {x}-\mathbf {y}$ is some permutation of the remainder-$k$ vertex of the canonical simplex,} \\ \noalign {\vspace {3pt}} 0,&\textrm{otherwise.} \end{cases} $$

Proof

Since every Delaunay cell is given by permuting the coordinates of the canonical simplex and translating the vertices, if x−y cannot be written as a permutation of a vertex of the canonical simplex, then they must not belong to any common Delaunay cell. So we can assume that x−y is a permutation of some vertex of the canonical simplex.

Without loss of generality, let y=0. Then any simplex containing y has a corresponding permutation in S _d+1. Among those, the number of simplices also containing x is the number of permutations that sends the remainder-k vertex of the canonical simplex to x.

Because components of x have two possible values, namely k and k−(d+1), and the number of these two values are d+1−k and k respectively, the number of possible permutations is (d+1−k)!k!. □

The next pair of lemmas illustrate how to integrate arbitrary functions over the canonical simplex with particular weighting schemes.

Lemma A.4

Let X be the union of Delaunay cells containing the origin in $A^{*}_{d}$, and let f:X→ℝ be a function. Then,

Proof

This is a direct consequence of Proposition 5.4. In short, weighted sampling corresponding to K _s can be achieved by uniformly sampling in the unit cube and projecting the sample onto H _d. □

Lemma A.5

Let $f:A^{*}_{d}\rightarrow\mathbb{R}$ be a function on $A^{*}_{d}$. Let $b:A^{*}_{d}\rightarrow\mathbb{R}$ and the vectors $\mathbf {w_{k}}$ be as in Sect. 5.2. Then,

Proof

The blur kernel b(x) is given by the convolution of $K^{0}_{b},\ldots,K^{d}_{b}$. Therefore, it is the accumulation of the weights corresponding to offsets in the d+1 axes that sum to x. The weight equals $\frac{1}{2}$ if the offset t _k in a given axis is zero, and $\frac{1}{4}$ if it is t _k=±1. Hence we can write it succinctly as $\frac{1}{2^{1+|t_{k}|}}$. The claim follows. □

1.2 A.2 Evaluation of SNR for ℤ^d

The kernel Z _x repeats periodically because of the translational symmetry of ℤ^d. In order to evaluate the expected value of an expression involving Z _x, it suffices to consider x∈[0,1]^d. In other words,

$$\eta={\frac{\int_{\mathbb{R}^d} N(\mathbf {y};5/6)^2d\mathbf {y}}{\int_{[0,1]^d} \int_{\mathbb{R}^d}\left(Z_{\mathbf {x}}(\mathbf {y}-\mathbf {x})-N(\mathbf {y}-\mathbf {x};5/6)\right)^2d\mathbf {y}d\mathbf {x}}}. $$

Expanding the quadratic term in the integrand, we obtain $\eta={\frac{(A.6)}{(A.5)+(A.6)-2(A.7)}}$, containing the following three expressions:

(A.5)

(A.6)

(A.7)

First of all, Lemma A.1 tells us that (A.6) is (10π/3)^−d/2.

(A.5) is similarly separable:

Note that B _p(x) and B _p′(x) are both nonzero only if p=p′ or |p−p′|=1. In the former case, the integral of B _p(x)B _p′(x) with respect to x over ℝ is 2/3; in the latter case, 1/6. Similarly, B _q(y) and B _q′(y) are both nonzero over y∈[0,1] iff q,q′=0 or 1. Setting q=q′ yields 1/3 for the value of the integral, and q≠q′ yields 1/6. Summing over all possible values of p′,q,q′ with nonzero summands, we obtain that (A.5) equals

The cross term (A.7) remains. We can similarly expand the expression and collect terms into the product of sums, and exploit the fact that only few values of p,q yield nonzero summands. Unfortunately, because of the Gaussian term, there is no succinct analytic expression.

Corollary A.6

The SNR term for ℤ^d is

$$\eta={\frac{\sqrt{3/(10\pi)}^d}{(29/96)^d+\sqrt{3/(10\pi)}^d-2\cdot 0.30456^d}}. $$

1.3 A.3 Evaluation of SNR for $A^{*}_{d}$

The error analysis for $A^{*}_{d}$ is analogous to that of ℤ^d, with a few subtle changes. First, K _x does not sum to 1 over its domain and should be normalized first before being compared to a Gaussian.

Lemma A.7

$$\int_{H_d}K_{\mathbf {x}}(\mathbf {y}-\mathbf {x})d\mathbf {y}=(d+1)^{d-1/2}. $$

Proof

The last line follows as the sum of splatting weights for x over all lattice points is also 1. □

Since K _x is translation-invariant (for translations by a lattice point) and symmetric with respect to permutation of the axes, it suffices to iterate x over the canonical simplex C ^∗, in order to compute the expected value of the L ₂ distance. Recall from Proposition 4.11 that the volume of C ^∗ is given by (d+1)^d−1/2/d!. Therefore,

$$\eta=\frac{(d+1)^{d-1/2}/d!\int_{H_d} N(\mathbf {y};2(d+1)^2/3)^2d\mathbf {y}}{\int_{C^*}\int_{H_d}[\frac{K_{\mathbf {x}}(\mathbf {y}-\mathbf {x})}{(d+1)^{d-1/2}} -N(\mathbf {y}-\mathbf {x};2(d+1)^2/3)]^2d\mathbf {y}d\mathbf {x}}. $$

η can be written as ${\frac{(A.9)}{(A.8)+(A.9)-2(A.10)}}$, where the components are as follows:

(A.8)

(A.9)

(A.10)

By Lemma A.1, (A.9) equals $(\frac{8(d+1)^{2}\pi}{3})^{-d/2}$.

Let us compute (A.8) without the constant term. The integrand may be expanded as follows:

Our expression is a weighted sum in which the weights are $b(\mathbf {q_{1}})b(\mathbf {q_{2}})$. By Lemma A.5, it is equivalent to

$$ \sum_{\substack{\mathbf {p_1}\in C^*\\\mathbf {p_2}\in C^*}} \sum_{\substack{\mathbf {t}\in\{-1,0,1\}^{d+1}\\\mathbf {s}\in\{-1,0,1\}^{d+1}}}[\cdots][\cdots] \prod_k\frac{1}{2^{1+|t_k|}}\frac{1}{2^{1+|s_k|}}, $$

(A.11)

where $\mathbf {q_{1}}=\sum t_{k}\cdot \mathbf {w_{k}}$ and $\mathbf {q_{2}}=\sum s_{k}\cdot \mathbf {w_{k}}$.

Note that the first bracketed term is effectively given by Lemma A.2. The second term is similar, but the integral is over H _d, so we must additively consider each simplex whose vertex contains both $\mathbf {p_{1}}-\mathbf {q_{1}}$ and $\mathbf {p_{2}}-\mathbf {q_{2}}$. The number of such simplices is given by Lemma A.3.

(A.11) is more attractive than its original form because its summands are now simple numerical quantities. There is no succinct analytic expression for the sum, but its numerical computation is straightforward. The number of summands is (d+1)²⋅3^2(d+1).

That leaves us with (A.10). Once again, the integral is expanded and manipulated:

Using Lemma A.5, we obtain

where $\mathbf {q}=\sum t_{k}\cdot \mathbf {w_{k}}$. By symmetry, we can swap out C ^∗ with the union of Delaunay cells containing the origin. Then we have an integral weighted by B _p(x) and B _p+q(y), which can be computed by Monte-Carlo sampling the (d+1)-dimensional unit cube, by Lemma A.4.

In summary, we iterate over t∈{−1,0,1}^d+1 and p∈C ^∗, and for each pair of vectors, sample x and y via Lemma A.4, and evaluate the Gaussian of y−x, weighted by $\prod_{k}\frac{1}{2^{1+|t_{k}|}}$. Multiplying the result by an appropriate constant yields (A.11).