A fast algorithm for computation of discrete Euclidean distance transform in three or more dimensions on vector processing architectures

Abstract

In this note, we introduce a function for calculating Euclidean distance transform in large binary images of dimension three or higher in Matlab. This function uses transparent and fast line-scan algorithm that can be efficiently implemented on vector processing architectures such as Matlab and significantly outperforms the Matlab’s standard distance transform function “bwdist” both in terms of the computation time and the possible data sizes. The described function also can be used to calculate the distance transform of the data with anisotropic voxel aspect ratios. These advantages make this function especially useful for high-performance scientific and engineering applications that require distance transform calculations for large multidimensional and/or anisotropic datasets in Matlab. The described function is publicly available from the Matlab Central website under the name “bwdistsc”, “Euclidean Distance Transform for Variable Data Aspect Ratio”.

Appendix: Proof of Lemma 1

Here, we present a formal proof of Lemma 1, namely that for a family of one-dimensional parabolas $\mathcal F =\{F_{y}(x)=|x-y|^2 + g(y) \}$ and an arbitrary parabola $f(x)=x^2 + b$ there can exist at most a single interval $(x_1,x_2)$ such that $f(x)<\mathcal G (x)$, where $\mathcal G (x)$ is the lower envelop of the parabolas in the family $\mathcal F$.

We begin by noting that, in general, three possibilities can exist above: $f(x)>\mathcal G (x)$ everywhere on $x$; $f(x)<\mathcal G (x)$ everywhere on $x$; and $f(x)<\mathcal G (x)$ at some points and $f(x)>\mathcal G (x)$ at other points. The first two possibilities are obviously accommodated by assuming $(x_1,x_2)=\emptyset$ and $(x_1,x_2)=(-\infty ,\infty )$, respectively.

For the third case, taking into account continuity of $\mathcal G (x)$ and$f(x)$, there should exist at least one crossing point such that $\mathcal G (x)=f(x)$. Let $x_1$ be such a crossing point, $\mathcal G (x_1)=f(x_1)$, and let $F_{y_1}(x)$ be the parabola in $\mathcal F$ such that $\mathcal G (x)=F_{y}(x)$ at $x=x_1$. We shall dismiss the trivial case $F_{y_1}(x)\equiv f(x)$, in which case $\mathcal G (x) = f(x)$ for the entire interval where $\mathcal G (x)=F_{y_1}(x)$ and, thus, $x_1$ cannot be a crossing point since $\mathcal G (x) > f(x)$ nowhere in that entire interval. For a real crossing point, therefore, $F_{y_1}(x) \not \equiv f(x)$ should hold. Then, by an obvious property of the intersection of the parabolas of the form given above, $f(x)>F_{y_1}(x)$ necessarily on a single semi-interval $(-\infty ,x_1)$ or $(x_1,\infty )$, and therefore $f(x)>F_{y_1}(x)\ge \mathcal G (x)$ on the same semi-interval. Without loss of generality we shall assume now that $f(x)>\mathcal G (x)$ on the semi-interval $(-\infty ,x_1)$.

Let us now assume that $x_2$ is another crossing point $\mathcal G (x_2)=f(x_2)$, and $F_{y_2}(x)$ is the parabola in $\mathcal F$ such that $\mathcal G (x)=F_{y}(x)$ at $x=x_2$. In that case, by the same obvious property of the intersection of the parabolas in $\mathcal F$ and $f(x)$, $f(x)>\mathcal G (x)$ necessarily on a single semi-interval $(-\infty ,x_2)$ or $(x_2,\infty )$. Since we already have that $f(x)>\mathcal G (x)$ in $(-\infty ,x_1)$ while $f(x_1)=\mathcal G (x_1)$, the only self-consistent possibility for $x_2$ to be a crossing point is $f(x)>\mathcal G (x)$ on $(x_2,\infty )$ while $x_1<x_2$.

Note that it is impossible to self-consistently add a third crossing point $x_3$, since that would imply that $f(x)>\mathcal G (x)$ on either $(-\infty ,x_3)$ or $(x_3,\infty )$, and that would necessarily violate either $f(x_1)=\mathcal G (x_1)$ or $f(x_2)=\mathcal G (x_2)$, or $f(x_3)=\mathcal G (x_3)$.

Therefore, there can be at most two crossing points for $f(x)$ and $\mathcal G (x)$, $x_1$ and $x_2$, and $f(x)<\mathcal G (x)$ can be achieved on at most a single interval $(x_1,x_2)$, respectively.