Minimal Paths by Sum of Distance Transforms

Abstract

Minimal path extraction is a frequently used tool in image processing with applications in for example centerline extraction and edge tracking. This paper presents results and methods for (i) extracting minimal paths and for (ii) utilizing local direction in the minimal path computation. Minimal path extraction is based on sum of distance transforms resulting in a stable method, without need for local decisions, which is needed for methods based on backtracking. Local direction utilization is discrete derivative based concepts such as the structure tensor, the Hessian and the Beltrami framework. The combination of minimal path extraction and local direction utilization gives a strong framework for minimal path extraction.

These keywords were added by machine and not by the authors. This process is experimental and the keywords may be updated as the learning algorithm improves.

1 Introduction

Many image processing and analysis applications use distance transforms [1, 3, 23, 25, 30, 31]. By taking the pixel intensity values into consideration, the homogeneous regions can be extracted [5, 6, 25, 26]. Not only pixel intensity values, but also other high-level features, based on for example pointwise texture features or local direction, can be utilized. Minimal path is a frequently used tool for, e.g., centerline or boundary extraction. One approach for minimal path extraction is based on backtracking the distance transform values from a given point to the closest seed point by steepest descent [15, 33]. A limitation with this approach is that it is based on local decisions, which, in general, does not give unique solutions. This manuscript presents an alternative approach based on sum of distance transforms for finding minimal paths. This approach has been used in mathematical morphology [27, 28], level-set based distance computation [13] and for a specific family of geodesic distance functions [10]. Here, we show that this approach is valid for a family of distance functions, including for example the fuzzy distance [6, 25, 27] and the geodesic distance [6, 32]. In Sect. 4, the approach presented in Sect. 3 is applied to minimal path formulations that utilize prior knowledge on the prefered direction at each point.

2 Notation

Given some cost function $C:\pi _{p,q} \mapsto \mathbb {R}$, where p, q are points in a digital image and $\pi _{p,q}$ is a path between p and q, the minimal path between two points p and q is $\arg \min _{\pi \in \varPi _{p,q}} C(\pi )$, where $\varPi _{p,q}$ is the set of all paths between p and q. The cost of the minimal cost path between p and q is ${C}_{\text {opt}}(p,q)=\min _{\pi \in \varPi _{p,q}} C(\pi )$.

Definition 1

A path $\pi =\langle p=p_0, p_1, \dots , p_n=q \rangle$ from p to q is a minimal path if there is no other path between p and q with lower cost/distance.

A function d such that, for all x, y, z in its domain is a metric if

• 1. $d(x,x)= 0$ (identity)

• 2. $d(x,y) > 0$ for all $x\ne y$ (positivity)

• 3. $d(x,y)=d(y,x)$ (symmetry)

• 4. $d(x,z)\le d(x,y)+d(y,z)$ (triangle inequality)

A distance transform from a single source point p is denoted $DT_p$ and

$$ DT_p(q)=d(p,q). $$

Intensity-weighted distance functions can be used to model some intuitive, expected distance in the image, such that the obtained minimal paths are optimal on dark or bright regions in the image, or at strong edges in the image and higher-level features can be used to extract minimal paths on a centerline or boundary of an object. Distance transforms utilizing higher-order features can for example be obtained by pre-filtering the image by, e.g., an edge detector or a “vesselness” filter [19], before using the distance functions below. Other approaches, where the prior information is incorporated in the distance function, will be given in Sect. 4.

3 Minimal Paths as Minima in the Sum of Distance Transforms

Minimal cost paths can be computed by an assumption that the minimal cost path between two points p and q is in the set of minima in $DT_{pq}(x)= DT_p(x)+DT_q(x)$ [10, 13, 27, 28].

The results in this section are valid for distance functions such that all metricity properties fulfilled, except that the strict positivity is replaced by non-negativity, and the following property holds.

Property 1

For any minimal path $\langle p=p_0 , p_1 , \dots , p_n=q \rangle$, it holds that

$$ d(p, q)=d(p,p_i)+d(p_i,q) \text { for all } i:0\le i \le n. $$

Examples of distance functions that fulfill Property 1 are

• The fuzzy distance [6, 25, 27]. The fuzzy distance is the minimal cost path when the cost between two adjacent points is

  $$c(p,q)= \frac{f(p)+f(q)}{2}\cdot \Vert p-q \Vert ,$$

  where f(p) is the intensity value at pixel p.

• The geodesic distance [6, 32]. The geodesic distance is the minimal cost path when the cost between two adjacent points is

  $$c(p,q)= \omega \left| f(p)-f(q)\right| + \Vert p-q \Vert . $$

  The parameter $\omega$ affects trade-off between the fuzzy membership values and the Euclidean distance.

In Theorem 1, we prove that the minimal path between two points is in the set of minimum of sum distance transforms from the points. For the proof we will use Lemmas 1 and 2, both related to the property that line segments between two given points contain points of the shortest path between the points [12, 20].

Lemma 1

Let $M=\{ x : d(p,q)=d(p,x)+d(x,q) \}$ and $R=\left\{ x:DT_{pq}(x)= \min _x DT_{pq}(x) \right\}$. Then $M=R$.

Proof

This is a direct consequence of the triangle inequality: For any point r, $d(p,q)\le d(p,r)+d(r,q)$. If $r\notin M$, then $DT_{pq}(r)=DT_p(r)+DT_q(r)=d(p,r)+d(r,q)>d(p,q) = \min _x DT_{pq}(x)$. If $r\in M$, then $DT_{pq}(r)=DT_p(r)+DT_q(r)=d(p,r)+d(r,q)=d(p,q) = \min _x DT_{pq}(x)$.   $\square$

Lemma 2

If $x\in M=\{ x : d(p,q)=d(p,x)+d(x,q) \}$, then x is in a minimal path between two points p and q.

Proof

If $d(p,q)=d(p,x)+d(x,q)$, then the path $\pi =\langle p=p_0, p_1, \dots , p_m =x \rangle \cdot \langle x=q_0, q_1, \dots , q_n =q \rangle$ is minimal, since its cost is $d(p,x)+d(x,q)=d(p,q)$, which is optimal by the triangular inequality.    $\square$

Theorem 1

Given a (pseudo)-metric such that Property 1 is fulfilled, the points in minimal paths between p and q are contained in the set

$$\begin{aligned} R=\left\{ x:DT_{pq}(x)=\min _x DT_{pq}(x) \right\} \end{aligned}$$

(1)

where $DT_{pq}(x)=DT_p(x)+DT_q(x)$.

Proof

The theorem follows directly from Lemmas 1 and 2.   $\square$

Theorem 1 can be applied by Algorithm 1. Examples of realizations of Algorithm 1 are shown in Fig. 1.

4 Minimal Paths with Relative Position Prior

A common approach for extracting vessels and tubular structures is to use a so-called vesselness filter [4, 7], which is based on the Hessian matrix, i.e., second derivatives of the intensity values in the image data.

Fig. 1.

Sets containing the minimal paths (in dark grey) between the upper left point to the lower right point obtained by Algorithm 1. (a): Input image with intensity values. (b–d): Fuzzy distance, (e–g): Geodesic distance with $\omega =1$.

4.1 Minimal Paths by Prefiltering the Input Image

Given a window, or smoothing, function W, the structure tensor S is obtained by integrating the outer product of the gradient $\nabla f = [f_x f_y]$ with itself over a region given by a window function:

$$S=\nabla f \nabla f^T * W=\left[ \begin{array}{cc}f_x f_x\ &{}f_x f_y\\ f_y f_x&{} f_y f_y\end{array} \right] * W, $$

where $*$ is convolution and W is some window function or a Gaussian. Note that the derivatives should be derived using an appropriate discrete operator [9, 34]. The structure tensor, which captures information on the predominant directions of the gradient vector field in a specified neighborhood of a point, and its variations, has been widely used for enhancing local structure and direction [2, 14, 18]. Let $\lambda _1$ be the largest eigenvalue of S and $\lambda _2$ the smallest eigenvalue. Then

• The derivative is close to zero in all directions in constant areas ($|\lambda _1| \approx |\lambda _2| \approx ~0$),

• there is a dominant direction with a strong gradient magnitude at edges ($|\lambda _1| \gg |\lambda _2| \approx 0$) and

• the derivative is strong also in the orthogonal directon at corners, so $|\lambda _1| \approx |\lambda _2| \gg 0$.

The structure tensor is based on the first derivative and is therefore not ideal for tubular structures since both $f_x$ and $f_y$ are close to zero at the center of a tubular structure. For tubular structures, the second derivative has properties that can be utilized since we can expect low intensity variation along the tubular structure and a local max/min in the perpendicular direction, i.e., a high second derivative.

To enhance tubular structures, the Hessian

$$H=\nabla ^2 f =\left[ \begin{array}{cc}f_{xx} &{} f_{xy}\\ f_{yx} &{} f_{yy}\end{array} \right] , $$

which is based on second derivatives, is a good choice. Again, let $\lambda _1$ be the largest eigenvalue and $\lambda _2$ the smallest eigenvalue. Now,

• $\lambda _1 \approx \lambda _2 \approx 0$ in constant areas,

• $|\lambda _1| \gg |\lambda _2| \approx 0$ at tubular structures and

• $\lambda _1 \approx \lambda _2$ and $|\lambda _2| \gg 0$ for blob-like structures.

Based on this, different measures of “vesselness” and “blobness” have been proposed [8, 16, 17].

4.2 Minimal Paths by Prior Knowledge on Prefered Directions

The local “direction” in the image can be extracted in many ways, for example by the Hessian or structure tensor as in Sect. 4.1 and the prefered local direction can be based on measures on local features corresponding to, e.g., vessels, blobs, edges or corners. The resulting preferred local directions can be represented by a vector field, which, in turn can be used to guide the path cost function. In this section, we will explain how the local direction can be used for guiding the minimal path extraction by including information based on the gradient in the definition of the grey-weighted distance function.

Minimal Paths Guided by the Gradient

Metrics can be defined on the topographic representation of the image intensity function f, or rather, on the graph of f defined as $\varphi (x,y)\rightarrow (x,y, f(x,y))$. If the metric is a (pseudo-) Riemannian metric expressed as a metric tensor g on the graph of f, then the metric can be pulled back to the image plane by the Beltrami framework, [11]. The pullback metric is given by $J^T g J$, where J is the Jacobian of $\varphi$, i.e.,

$$J=\left[ \begin{array}{cc}1 &{} 0\\ 0 &{} 1\\ f_x &{} f_y\end{array} \right] .$$

Following the Beltrami framework, the distance (pullback metric) between two points p and $q=p+u$ is given by

$$\begin{aligned} c(p,q)=\sqrt{ u^T B u }, \end{aligned}$$

(2)

where

$$B=J^T g J = \left[ \begin{array}{cc}1 + f_x^2 &{} f_x f_y\\ f_y f_x &{} 1 + f_x^2\end{array} \right] = \left[ \begin{array}{cc}1 &{} 0\\ 0 &{} 1\end{array} \right] + \nabla f \nabla f^T. $$

If we set

$$g=\left[ \begin{array}{ccc}1 &{} 0 &{} 0\\ 0 &{} 1 &{} 0\\ 0 &{} 0 &{} 1\end{array} \right] ,$$

we get the geodesic path length on the surface. Similar to the geodesic distance, a topographic representation of the image is considered, where the intensity values correspond to the height.

Also similar to the geodesic distance, we can set a parameter that gives the trade-off between the intensity range and the spatial distance between x and y;

$$g=\left[ \begin{array}{ccc}1 &{} 0 &{} 0\\ 0 &{} 1 &{} 0\\ 0 &{} 0 &{} \alpha ^2\end{array} \right] ,$$

corresponds to scaling of f with $\alpha$, leading to the following following effect on B:

$$B= \left[ \begin{array}{cc}1 &{} 0\\ 0 &{} 1\end{array} \right] + \alpha ^2 \nabla f \nabla f^T, $$

This idea has also been generalized by using the structure tensor instead of the outer product of gradients, i.e., exchanging B for

$$S= \left[ \begin{array}{cc}1 &{} 0\\ 0 &{} 1\end{array} \right] + \alpha ^2 \nabla f \nabla f^T * G_{\sigma }, $$

where $G_{\sigma }$ is a Gaussian function with standard deviation $\sigma$ [18]. This is a generalization since the standard Beltrami pullback metric is achieved by setting $\sigma =0$ and the generalization extends the Beltrami framework to penalize both image gradients (edges) and ridges (lines) by the structure tensor [18].

Note that this general framework allows other (pseudo-) metrics to be used and pulled back to the image plane. For example,

$$g=\left[ \begin{array}{ccc}0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 1\end{array} \right] $$

gives a pseudo-metric that does not take spatial distance in the image plane into account. This idea has been used for edge detection in color images [24].

4.3 General Metric Tensor Minimal Path-Cost

By combining path costs given by the intensity values in the image, as given in Eq. 2, with the cost given by the local direction, a more general path cost framework is obtained by using a path cost functions on the form

$$ c(p,p+u)= E\left( c_{\text {intensity}}(p,p+u), c_{\text {tensor}}(p,p+u) \right) , $$

where E is an appropriate combination of a grey-weighted path cost function $c_{\text {intensity}}$ and a tensor-based path cost function $c_{\text {tensor}}$ that takes local direction into account.

The tensor-based path cost function $c_{\text {tensor}}$ is defined by a metric tensor M(p) i.e.,

$$\begin{aligned} c_{\text {tensor}}(p,p+u)=\sqrt{ u^T M(p) u }, \end{aligned}$$

(3)

where M(p) typically is restricted to symmetric and positive definite matrices. The positive definiteness is sometimes relaxed to non-negativity and sometimes M(p) is restricted to diagonal matrices [22].

Example 1. Circular paths, for example, annual rings in log ends, can be extracted by utilizing the polar geometry of the paths [21, 29]. By matrix diagonalization, we can define a tensor field where eigenvalues and eigenvectors are given by P and D:

$$\begin{aligned} M(p)=P D P^{-1}, \end{aligned}$$

(4)

where

$$\begin{aligned} P=\left[ \begin{array}{rc}- \sin \theta &{} \cos \theta \\ \cos \theta &{} \sin \theta \end{array} \right] \text { and } D=\left[ \begin{array}{cc}1/r^2 &{} 0\\ 0 &{} 1\end{array} \right] \end{aligned}$$

(5)

and r and $\theta$ are polar coordinate representation of p.

A path cost function based on multiplicative weighting of image intensity differences and spatial distance in polar coordinates can be defined by

$$c(p,p+v)= \frac{f(p)+f(p+v)}{2} \cdot \sqrt{ v^T M(p) v }, $$

with M(p) as given in Eq. 4. An example of distance transforms and a minimal path using this path-cost function is shown in Fig. 2. By applying a standard Taylor expansion of $(p-(p+v))' P D P^{-1} (p-(p+v))$, P and D are as in Eq. 5 and

$$p=\left[ \begin{array}{c} r \cos \theta \\ r \sin \theta \end{array}\right] \text { and } p+v=\left[ \begin{array}{c} (r+\varDelta r) \cos (\theta + \varDelta \theta )\\ (r+ \varDelta r) \sin (\theta + \varDelta \theta )\end{array}\right] ,$$

we get

$$ c(p,p+v)\approx \frac{f(p)+f(p+v)}{2} \cdot \sqrt{ \varDelta \theta ^2 + \varDelta r^2 + R_2(\varDelta r,\varDelta \theta ) }, $$

where $R_2(\varDelta r,\varDelta \theta )$ is the remainder term (a weighted sum of monomials of degrees higher than 2). This leads to a special case of the path cost function used in [21] for extracting circular paths:

$$c(p,p+v)= \frac{f(p)+f(p+v)}{2} \cdot \sqrt{ (\omega _r \cdot \varDelta r(p,v))^2 + (\omega _\theta \cdot \varDelta \theta (p,v))^2}, $$

where $\varDelta r(p,v)=|r_p - r_{p+v} |$ and $\varDelta \theta (p,v)=\min \left( |\theta _p - \theta _{p+v}| , 2\pi -|\theta _p - \theta _{p+v}| \right)$.

Fig. 2.

Distance transform in polar coordinates on a constant intensity input image by the general metric tensor and minimal path extraction. (a): Visualization of the tensor field. (b): The sum of distance transforms from two source points on equal distance from the center of the image using polar weights. (c): The set containing the minimal path by Algorithm 1.

5 Conclusion

Results on minimal path extraction for path-based distance functions have been presented together with results on path cost functions that utilizes local direction. This combination gives a general and practically useful framework for minimal path extraction with many potential applications based on for example, centerline extraction and edge tracking.