Fast Asymmetric Fronts Propagation for Image Segmentation

The remaining of this paper is organized as follows: In Section 2, we introduce the geodesic distance map associated to a general Finsler metric, the Voronoi regions and the relevant numerical tool. Section 3 presents the construction of the asymmetric Finsler metrics associated to different image segmentation applications. The numerical considerations for the Finsler metrics-based fronts propagation are introduced in Section 4. The experimental results and the conclusion are respectively presented in Sections 5 and 6.

In this section, we consider a more general case for which a family of source point sets are provided. Suppose that these source point sets are denoted by $\mathfrak{s}_{k}$ and are indexed by $k\in\{1,2,\cdots,n\}$ with $n$ the total number of source point sets. For the sake of simplicity, we note $\mathfrak{s}=\cup_{k=1}^{n}\mathfrak{s}_{k}$.

For a given geodesic metric $\mathcal{F}$, we can compute the respective geodesic distance map $\mathcal{U}_{k}$ associated to each source point set $\mathfrak{s}_{k}$ by Eq. (6). A Voronoi index map is a function $\mathcal{L}:\Omega\to\{1,2,\cdots,n\}$ such that for any source point $\mathbf{x}\in\mathfrak{s}_{k}$

and for any domain point $\mathbf{x}\in\Omega\backslash\mathfrak{s}$, the Voronoi index $\mathcal{L}(\mathbf{x})$ is identical to the index of the closest source point set in the sense of the geodesic distance arbelaez2004energy ; benmansour2009fast . One can construct a Voronoi index map $\mathcal{L}$ in terms of the geodesic distance maps $\mathcal{U}_{k}$ ($1\leq k\leq n$) by

By the Voronoi index map $\mathcal{L}$, the whole domain $\Omega$ can be partitioned into $n$ Voronoi regions $\mathcal{V}_{k}\subset\Omega$

The common boundary $\Gamma_{i,j}:=\partial\mathcal{V}_{i}\cap\partial\mathcal{V}_{j}$ of two adjacent Voronoi regions $\mathcal{V}_{i}$ and $\mathcal{V}_{j}$ is comprised of a collection of equidistant points to the source point sets $\mathfrak{s}_{i}$ and $\mathfrak{s}_{j}$, i.e.,

Finally, we consider a geodesic distance map $\mathcal{U}_{\mathfrak{s}}$ associated to the set $\mathfrak{s}=\cup_{k}\mathfrak{s}_{k}$ which can be expressed by

where $\mathcal{U}_{k}$ is the geodesic distance map with respect to the source point set $\mathfrak{s}_{k}$ indexed by $k$.

Many approaches bornemann2006finite ; weber2008parallel ; yatziv2006n ; rouy1992viscosity can be used to estimate the geodesic distance map $\mathcal{U}_{\mathfrak{s}}$. Among them, the fast marching method sethian1999fast ; tsitsiklis1995efficient solves the Eikonal equation in a very efficient way. It has a similar distance estimation scheme with Dijkstra’s graph-based shortest path algorithm dijkstra1959note . One crucial ingredient of the fast marching method is the construction of the stencil map $\Lambda$, where $\Lambda(\mathbf{x})$ defines the neighbourhood of a grid point $\mathbf{x}$. The isotropic fast marching methods sethian1999fast ; tsitsiklis1995efficient are established on an orthogonal 4-connectivity neighbourhood system, which suffers some difficulties for the distance computation associated to asymmetric Finsler metrics mirebeau2014efficient . In order to find accurate solutions to the Finsler Eikonal equation, more complicated neighbourhood systems are taken into account mirebeau2014anisotropic ; mirebeau2014efficient ; sethian2003ordered ; mirebeau2017fast. These neighbourhood systems or stencils are usually constructed depending on the geodesic metrics. In this paper, we adopt the Finsler variant of the fast marching method proposed by Mirebeau mirebeau2014efficient . It invokes a geometry tool of anisotropic stencil refinement and leads to a highly accurate solution, but requires a low computation complexity.

In this section, we present the construction method of the Finsler metric which is suitable for fronts propagation and image segmentation. Suppose that a vector field $\mathfrak{g}:\Omega\to\mathbb{R}^{2}$ has been provided such that $\mathfrak{g}(\mathbf{x})$ points to the object edges at least when $\mathbf{x}$ is nearby them. In this case, the orthogonal vector field $\mathfrak{g}^{\perp}$ indicates the tangents of the edges.

Basically, the Eikonal equation-based fronts propagation models malladi1998real perform the segmentation scheme through a geodesic distance map. In order to find a good solution for image segmentation, the used geodesic metric should be able to reduce the risk of front leaking problem. For this purpose, we search for a direction-dependent metric $\mathcal{F}_{\mathfrak{g}}$ satisfying the following inequality

Recall that for an edge point $\mathbf{x}$, both the feature vectors $\mathfrak{g}^{\perp}(\mathbf{x})$ or $-\mathfrak{g}^{\perp}(\mathbf{x})$ are propositional to the tangent of the edge at $\mathbf{x}$. When the fast marching front arrives at the vicinity of image edges, it prefers to travel along the edge feature vectors $\mathfrak{g}^{\perp}(\mathbf{x})$ and $-\mathfrak{g}^{\perp}(\mathbf{x})$, instead of passing through the edges, i.e., prefers to travel along the direction $-\mathfrak{g}(\mathbf{x})$.

The inequality (20) requires the geodesic metric $\mathcal{F}_{\mathfrak{g}}$ to be anisotropic and asymmetric with respect to its second argument. Thus, we consider a Finsler metric with a Randers form randers1941asymmetrical involving a symmetric quadratic term and a linear asymmetric term for any $\mathbf{x}\in\mathbb{R}^{2}$ and any vector $\vec{u}\in\mathbb{R}^{2}$

where $\mathcal{M}_{\mathfrak{g}}:\Omega\to S^{+}_{2}$ is a positive symmetric definite tensor field and $\vec{\omega}_{\mathfrak{g}}:\Omega\to\mathbb{R}^{2}$ is a vector field that is sufficiently small. The function $\mathfrak{C}:\Omega\to\mathbb{R}^{+}$ is a positive scalar-valued potential which gets small values in the homogeneous regions and large values around the image edges. It can be derived from the image data such as the coherence measurements of the image features, which will be discussed in detail in Section 4.3.

The tensor field $\mathcal{M}_{\mathfrak{g}}$ and the vector field $\vec{\omega}_{\mathfrak{g}}$ should satisfy the constraint

in order to guarantee the positiveness mirebeau2014efficient of the Randers metric $\mathcal{F}$.

We reformulate the Randers metric $\mathcal{F}_{\mathfrak{g}}$ in Eq. (21) as

where $\mathcal{G}_{\mathfrak{g}}:\Omega\times\mathbb{R}^{2}\to[0,\infty]$ is still a Randers metric which can be formulated by

The remaining part of this section will be devoted to the construction of the Randers metric $\mathcal{G}_{\mathfrak{g}}$ in terms of the vector field $\mathfrak{g}$ which is able to characterize the directions orthogonal to the image edges.

Let us define a new vector field $\bar{\mathfrak{g}}:\Omega\to\mathbb{R}^{2}$ by

The tensor field $\mathcal{M}_{\mathfrak{g}}$ used in Eq. (21) can be constructed dependently on two scalar-valued coefficient functions $\eta_{1}$ and $\eta_{2}$ such that

where $\bar{\mathfrak{g}}^{\perp}(\mathbf{x})$ is the orthogonal vector of $\bar{\mathfrak{g}}(\mathbf{x})$ and $\otimes$ denotes the tensor product, i.e., $\mathbf{u}\otimes\mathbf{u}=\mathbf{u}\mathbf{u}^{T}$. Note that the eigenvalues of $\mathcal{M}_{\mathfrak{g}}(\mathbf{x})$ are $\eta_{1}^{2}(\mathbf{x})/\|\mathfrak{g}(\mathbf{x})\|^{2}$ and $\eta_{2}(\mathbf{x})/\|\mathfrak{g}(\mathbf{x})\|^{2}$, respectively corresponding to the eigenvectors $\mathfrak{g}(\mathbf{x})/\|\mathfrak{g}(\mathbf{x})\|$ and $\mathfrak{g}^{\perp}(\mathbf{x})/\|\mathfrak{g}(\mathbf{x})\|$.

The vector $\vec{\omega}_{\mathfrak{g}}(\mathbf{x})$ is positively collinear to field $\mathfrak{g}(\mathbf{x})$ for all $\mathbf{x}\in\Omega$

where $\tau:\Omega\to\mathbb{R}$ is a scalar-valued coefficient function.

We estimate the coefficient functions $\eta_{1}$, $\eta_{2}$ and $\tau$ through two cost functions $\psi_{\rm f},\,\psi_{\rm b}:\Omega\to(1,\,+\infty)$, which assign the cost values $\psi_{\rm f}(\mathbf{x})$, $\psi_{\rm b}(\mathbf{x})$ and $1$ to the Randers metric $\mathcal{G}_{\mathfrak{g}}$ respectively along the directions $\mathfrak{g}(\mathbf{x})$, $-\mathfrak{g}(\mathbf{x})$ and $\mathfrak{g}^{\perp}(\mathbf{x})$ for any point $\mathbf{x}\in\Omega$ such that

Combining Eqs. (25) and (27) yields that

for any $\mathbf{x}\in\Omega$.

The positive symmetric definite tensor field $\mathcal{M}_{\mathfrak{g}}$ and the vector field $\vec{\omega}_{\mathfrak{g}}$ thus can be respectively expressed in terms of the cost functions $\psi_{\rm f}$ and $\psi_{\rm b}$ by

and

Based on the tensor field $\mathcal{M}_{\mathfrak{g}}$ and the vector field $\vec{\omega}_{\mathfrak{g}}$ respectively formulated in Eqs. (3.1) and (30), the positiveness constraint (22) is satisfied due to the assumption that $\psi_{\rm f}(\mathbf{x})>1$ and $\psi_{\rm b}(\mathbf{x})>1,\,\forall\mathbf{x}\in\Omega$. The cost functions $\psi_{\rm f}$ and $\psi_{\rm b}$ can be derived from the image edge information such as the image gradients, which will be discussed in Section 4.

Note that if we set $\psi_{\rm f}\equiv\psi_{\rm b}$, the vector field $\vec{\omega}_{\mathfrak{g}}$ will vanish, i.e., $\vec{\omega}_{\mathfrak{g}}\equiv\mathbf{0}$ (see Eq. (30)). In this case, one has $\langle\vec{\omega}_{\mathfrak{g}}(\mathbf{x}),\mathbf{u}\rangle=0$ for any point $\mathbf{x}\in\Omega$ and any vector $\mathbf{u}\in\mathbb{R}^{2}$, leading to a special form of the Randers metric $\mathcal{G}_{\mathfrak{g}}$. This special form is a symmetric (potentially anisotropic) Riemannian metric $\mathcal{R}(\mathbf{x},\mathbf{u})=\|\mathbf{u}\|_{\mathcal{M}_{\mathfrak{g}}(\mathbf{x})}$ which depends only on the tensor field $\mathcal{M}_{\mathfrak{g}}$.

A basic tool for studying and visualizing the geometry distortion induced from a geodesic metric is the Tissot’s indicatrix defined as the collection of control sets in the tangent space chen2017global . For an arbitrary geodesic metric $\mathcal{F}:\Omega\times\mathbb{R}^{2}\to[0,\infty]$, the control set $\mathcal{B}(\mathbf{x})$ for any point $\mathbf{x}\in\Omega$ is defined as the unit ball centered at $\mathbf{x}$ such that

We demonstrate the control sets $\mathcal{B}(\mathbf{q})$ in Fig. 3 for the Randers metric $\mathcal{G}_{\mathfrak{g}}(\mathbf{q},\cdot)$ with different values of $\psi_{\rm f}(\mathbf{q})$ and $\psi_{\rm b}(\mathbf{q})$ at a point $\mathbf{q}\in\Omega$. The vector $\mathfrak{g}(\mathbf{q})$ is set as

In Fig. 3a, we show the control sets for the Randers metric $\mathcal{G}_{\mathfrak{g}}$ with respect to different values of $\psi_{\rm b}(\mathbf{q})$ and a fixed value $\psi_{\rm f}(\mathbf{q})=5$. One can point out that the common origin of these control sets have shifted from the original center of the ellipses²²2These ellipses are the boundaries of the control sets. due to the asymmetric property as formulated in Eq. (3). In Fig. 3b, the control sets for the Randers metric $\mathcal{G}_{\mathfrak{g}}$ associated to $\psi_{\rm f}(\mathbf{q})=\psi_{\rm b}(\mathbf{q})>1$ are demonstrated, where $\mathcal{G}_{\mathfrak{g}}(\mathbf{q},\cdot)$ gets to be anisotropic and symmetric on its second argument. When $\psi_{\rm f}(\mathbf{q})=\psi_{\rm b}(\mathbf{q})=1$, the values of the Randers metric $\mathcal{G}_{\mathfrak{g}}(\mathbf{q},\vec{u})$ turn to be invariant with respect to $\vec{u}$ as shown in Fig. 3c. In this case, the tensor $\mathcal{M}_{\mathfrak{g}}(\mathbf{q})$ is propositional to the identity matrix.

In Fig. 4, we show the geodesic distance maps associated to $\mathcal{G}_{\mathfrak{g}}$ with different values of the cost functions $\psi_{\rm f}$ and $\psi_{\rm b}$. The vector field $\mathfrak{g}$ is set to

In Figs. 4a and  4c, we can see that the geodesic distance maps have strongly asymmetric and anisotropic appearance. In Fig. 4b, we observe that the geodesic distance map appears to be symmetric and strongly anisotropic. This is because the respective propagation speed of the fast marching fronts along the directions $(\cos(\pi/4),\sin(\pi/4))^{T}$ and $-(\cos(\pi/4),\sin(\pi/4))^{T}$ is identical to each other.

In this section, we present the computation method for the cost functions $\psi_{\rm f}$ and $\psi_{\rm b}$ using the image edge information, based on which the image data-driven Randers metric $\mathcal{G}_{\mathfrak{g}}$ can be obtained.

Let $\mathbf{I}=(I_{1},I_{2},I_{3}):\Omega\to\mathbb{R}^{3}$ be a vector-valued image in the chosen color space and let $G_{\sigma}$ be a Gaussian kernel with variance $\sigma$ (We set $\sigma=1$ through all the experiments of this paper). The gradient of the image $\mathbf{I}$ at each point $\mathbf{x}=(x,y)$ is a $2\times 3$ Jacobian matrix

involving the Gaussian-smoothed first-order derivatives along the axis directions $x$ and $y$

Let $\rho:\Omega\to\mathbb{R}$ be an edge saliency map. It has high values in the vicinity of image edges and low values inside the flatten regions. For each domain point $\mathbf{x}\in\Omega$, the value of $\rho(\mathbf{x})$ can be computed by the Frobenius norm of the Jacobian matrix $\nabla\mathbf{I}_{\sigma}(\mathbf{x})$

For a gray level image $I:\Omega\to\mathbb{R}$, the edge saliency map $\rho$ can be simply computed by the norm of the Euclidean gradient of the image $I$ such that

The vector field $\mathfrak{g}$ for the construction of the Randers metric $\mathcal{G}_{\mathfrak{g}}$ can be obtained by normalizing the vector field $\vec{h}$

The cost functions $\psi_{\rm f}$ and $\psi_{\rm b}$ used in Eq. (27) for the foreground and background segmentation application can be expressed for any $\mathbf{x}\in\Omega$ by

where $\alpha_{\rm f}$ and $\alpha_{\rm b}$ are nonnegative constants dominating how anisotropic and asymmetric the Randers metric $\mathcal{G}_{\mathfrak{g}}$ is.

Once the cost functions $\psi_{\rm f}$ and $\psi_{\rm b}$ have been computed by Eqs. (40) and (41), we can construct the tensor filed $\mathcal{M}_{\mathfrak{g}}$ and the vector field $\vec{\omega}_{\mathfrak{g}}$ respectively via Eqs. (3.1) and (30). Indeed, one has $\psi_{\rm f}(\mathbf{x})\approx\psi_{\rm b}(\mathbf{x})\approx 1$ for the points $\mathbf{x}$ located in the homogeneous region of the image $\mathbf{I}$ where $\rho(\mathbf{x})\approx 0$. In this case, the data-driven Randers metric $\mathcal{G}_{\mathfrak{g}}(\mathbf{x},\cdot)$ expressed in Eq. (24) approximates to be an isotropic Riemannian metric. For each point $\mathbf{x}$ around the image edges such that the value of $\rho(\mathbf{x})$ is large, the Randers metric $\mathcal{G}_{\mathfrak{g}}(\mathbf{x},\cdot)$ will appear to be strongly anisotropic and asymmetric.

In this section, we take into account a feature vector field $\mathbf{p}:\Omega\to\mathbb{R}^{2}$ which extracts local orientation features from the image. A feature vector $\mathbf{p}(\mathbf{x})$ characterizes the orientation that a tubular structure should have at a point $\mathbf{x}$ belonging to this structure. The feature vector field $\mathbf{p}$ can be estimated by the steerable filters jacob2004design; law2008three; frangi1998multiscale; freeman1991design.

Based on the cost functions $\psi_{\rm f}$ and $\psi_{\rm b}$ in Eqs. (39) and (40), we are able to build a positive symmetric definite tensor field $\mathcal{M}_{\mathfrak{g}}$ for the Randers metric $\mathcal{G}_{\mathfrak{g}}$. For each point $\mathbf{x}$ inside the tubular structures where the gray levels vary slowly, the value of the gradient norm $\rho(\mathbf{x})$ nearly vanishes, leading to $\psi_{\rm b}(\mathbf{x})\approx\psi_{\rm f}(\mathbf{x})\approx 1$. In this case, the Randers metric $\mathcal{G}_{\mathfrak{g}}(\mathbf{x},\mathbf{u})$ is approximately independent of the directions $\mathbf{u}$ (see Section 4.1). However, with respect to tubular structure segmentation application, we expect that inside the tubular structure the fast marching fronts travel faster along the directions $\mathbf{p}(\cdot)$ or $-\mathbf{p}(\cdot)$ than along the directions $\mathbf{p}^{\perp}(\cdot)$ or $-\mathbf{p}^{\perp}(\cdot)$ in order to reduce the risk of the front leakages. For this purpose, we make use of a new tensor field $\tilde{\mathcal{M}}_{\mathfrak{g}}:\Omega\to S_{2}^{+}$ based on $\mathcal{M}_{\mathfrak{g}}$ in Eq. (3.1) such that

where $\mu\in\mathbb{R}^{+}$ is a constant. It is relevant to the anisotropy property of the tensor field $\tilde{\mathcal{M}}_{\mathfrak{g}}$. In the numerical experiments of this paper, we set $\mu=\|\psi_{\rm f}\|^{2}_{\infty}$.

In practice, we observe that inside the tubular structures, the feature vector field $\mathbf{p}$ can be well approximated by the vector field $\mathfrak{g}^{\perp}$ which is the orthogonal vector field of $\mathfrak{g}$ derived from Eqs. (36) and (39). In order to reduce the computation time, we construct the tensor field $\tilde{\mathcal{M}}_{\mathfrak{g}}$ by the vector field $\mathfrak{g}$ such that

In this case, we obtain a new Randers metric $\tilde{\mathcal{G}}_{\mathfrak{g}}$ for the application of tubularity segmentation, which depends on the tensor field $\tilde{\mathcal{M}}_{\mathfrak{g}}$ and can be formulated as

Note that we build $\tilde{\mathcal{G}}_{\mathfrak{g}}$ by using the same vector field $\vec{\omega}_{\mathfrak{g}}$ with the Randers metric $\mathcal{G}_{\mathfrak{g}}$. Based on the new Randers metric $\tilde{\mathcal{G}}_{\mathfrak{g}}$, inside the tubular structures the fast marching fronts will travel fast along the directions $\mathbf{p}(\cdot)$ or $-\mathbf{p}(\cdot)$, but slow when the fronts arrive at the boundaries and slower when the fronts tend to pass through them.

We present the computation methods for the potential function $\mathfrak{C}$ used by the data-driven Randers metric that is formulated in Eq. (23). Basically, the potential function $\mathfrak{C}$ should have small values in the homogeneous regions and large values in the vicinity of the image edges.