Using the Fast-Marching for segmentation

Originally introduced by Malladi and Sethian [111], the can be used as a fast initialization algorithm for image segmentation. Despite the fact that it requires low computational cost, it has several drawbacks, in particular that the speed $F$ in equation (4.20) is always positive or always negative. Therefore, the front evolves and propagates in the whole image domain, without being stop by any edge, and it cannot be used for cases with curvature-dependent speed functions. The stopping criterion still relies on a user perception of the segmentation, like visual inspection of the domain visited by the algorithm during propagation, or the choice for an appropriate stopping time, as done in [111]. But the nice segmentation properties of the algorithm shall be taken into account, as shown in figure 4.10, where we segment a brain on the basis of the method described in [111].

Still, the stopping criterion in this case was to visually control the segmentation process, because the automatic detection of the crossing of the interesting edges is difficult to implement.

Noticing that is used to model the continuous formulation of the watershed transform, we have developed an interactive segmentation tool based on the flooding process of this morphological algorithm introduced in [180]. It is based on the same principle that segmentation boundaries are defined by pixels where several evolving fronts collide.

Comparison with the watershed algorithm

The classical morphological segmentation paradigm [156] is based on the watershed transform [181]. The method is very simple:

• the gradient image $\Vert\nabla I\Vert$ is computed (and enhanced);
• for each object of interest, an inside particle is detected (either interactively or automatically);
• flooding waves are propagated from the set of markers and flood the topographic surface $\Vert\nabla I\Vert$.

The points where the flooding waves meet each other form the segmentation boundaries. Different regions between boundaries are called catchment basins. Usually markers are the regional minima of the gradient image. An example of the watershed transform is shown in figure 5.9.

Figure: Watershed transform of a 2D X Ray image of the heart (LV): The seed points are regional minima of the gradient of image 5.9-left; the flooding was implemented in the continuous framework of , and the regional minima were filtered in order to have a smaller set of markers.

Very often, the minima are extremely numerous, leading to an over-segmentation, therefore for practical cases, the watershed transform will take as seed points a smaller set of markers, identified through pre-processing techniques.

The flooding waves are propagated according to the topographic map $\Vert\nabla I\Vert$, where the pixels belonging to a catchment basin are nearer to its corresponding regional minimum than to any other minima. In this framework, the construction of catchment basins can be seen as a shortest path problem between a marker and the image points. It corresponds to compute the gray-weighted distance transform (GWDT) of the image to the set of markers. There are several ways to compute this : as shown in [118], viewing the topographic image domain as a refractive index, as explained in chapter 1, the distance between a marker and any point in the image is the line integral of the penalty $\Vert\nabla I\Vert$ along the optical path length. We reformulate 4.22 with

$$\Vert\nabla T(\mathbf{x}) \Vert = \Vert\nabla I(\mathbf{x})\Vert$$ (5.7)

and the can be computed on the whole image domain, using the markers as sources pixels $p$ where $T(\mathbf{p}) = 0$.

Maragos has shown in [118] that the continuous segmentation approach based on the outperforms the discrete segmentation results. However several problems remain important:

• the final segmentation relies on the number of the markers: too many markers lead to an over segmentation of the image, and merging algorithms are needed as post-processing;
• it relies also on the quality of the edge detector applied to the real image: if there is a gap in the edge information, the nearest marker will flood the adjacent region.

We have devised a very simple method based on user interactivity, and different front speeds, based on the formulation of the flooding of the whole image. Each front, initialized by a user-defined seed point, will propagate according to its own propagation speed, derived from local image information (and not just gradients), and will stop when it meet other propagating fronts. This method has been implemented to provide a quick and dirty initialization for the method.

Interactive segmentation with the Fast-Marching

Flooding algorithm

In the following, we assume that we have an undefined number of seed points for front propagation, that can be labeled with a known number of labels. We detail below the flooding algorithm for multi-label front propagation.

Definition

• a set of starting point $\mathbf{p}_1,\ldots,\mathbf{p}_n$, located inside the image domain;
• each starting point $\mathbf{p}_i$ is assigned a label $l_i$ $i\in[1;n]$ (not necessarily different);
• a label map ${\cal L}$, for controlling collision;
• a penalty image ${\cal P}$ which will drive the front propagation, and which is a function of the position and the label of the front;
• We initialize the usual set of data-structures for front propagation, including an action map ${\cal A}$ and (only) one min-heap structure ${\cal H}$;

Initialization

• we initialize a classical front propagation method, setting ${\cal A}(p_i) = 0$ $\forall i\in[1;n]$ and storing all seed points in the min-heap structure;
• ${\cal L}(\mathbf{p}_i) = l_i$, ${\cal L} = -1$ elsewhere;

Loop: at any iteration

• Let $\mathbf{x}_{min}$ be the Trial point with the smallest action ${\cal A}$;
• Move it from the Trial to the Alive set (i.e. ${\cal A}_{\mathbf{x}_{min}}$ is frozen);
• Update ${\cal P}({\cal L}(\mathbf{x}_{min}))$
• For each neighbor $\mathbf{x}$ (6-connexity in 3D) of $\mathbf{x}_{min}$:
  • If $\mathbf{x}$ is Far
    1. add it to the Trial set;
    2. ${\cal L}(\mathbf{x}) = {\cal L}(\mathbf{x}_{min})$
    3. compute ${\cal A}$ according to equation
       

       $$\Vert \nabla {\cal A} \Vert = {\cal P}\left({\cal L}(\mathbf{x}_{min})\right)$$ (5.8)

       
       
       using the up-wind discretization scheme, involving only the neighbors of $\mathbf{x}$ with label ${\cal L}(\mathbf{x}_{min})$;
  • If $(\mathbf{x})$ is Trial
    1. recompute the action ${\cal A}(\mathbf{x})$ with equation( 5.8), involving only the neighbors of $\mathbf{x}$ with label ${\cal L}(\mathbf{x}_{min})$;
    2. if the new value is smaller, then ${\cal L}(\mathbf{x}) = {\cal L}(\mathbf{x}_{min})$

Termination

• stop when all pixels on the image domain are visited;
• the label map ${\cal L}$ represents the final segmentation into $k$ regions $R_k$.

In this algorithm, we have deliberately left blank the update procedure of the potential. Using a very simple potential based on the grey level information. Having the $n$ seed points $\mathbf{p}_1,\ldots,\mathbf{p}_n$, we can set at initialization

$$\hbox{ if } {\cal L}(\mathbf{x}) = {\cal L}(\mathbf{p}_i) \... ...}(\mathbf{x}) = \max( I(\mathbf{x}) - I(\mathbf{p}_i) , 0) + w$$ (5.9)

Then the speed is inversely proportional to the difference between the grey-level of the starting seed point and other points in the image. Result of this strategy can be observed in figure 5.10 where different seed points were manually located on the image 5.10-left, leading to the segmentation of figure 5.10-right.

Figure 5.10: Front competition on a 2D X Ray image of the heart (LV): The seed points were manually located on the dataset on the left image; they are visible as white dots on the right image.

The resulting is not convincing, but the use of this kind of information to differentiate the different front speeds has a bigger potential than the classical watersheds formulation. If there is a gap in the gradients between two basins, one of the two fronts could flood completely into the adjacent basin with the watersheds, whereas the speed in the adjacent region could still be discriminating the unwanted flooding, in the case of our competitive fronts algorithm. This is mainly the difference between a difference of a plateau and a barrier in the penalty (where plateau is a region with small variations in the intensity).

Thus, different strategies can be implemented, as many as there are different possible potentials.

Setting the penalty information

We can devise a very simple algorithm that do not use the information coming from the seed points, but for the $m$ first visited voxels for each label. This reduces greatly the influence of the user interaction in setting the markers position. Therefore for each label,

• At iteration $i$ in the Fast-Marching we examine the point ${\mathbf x}$ in the min-heap whose action is minimal;
• We remove it from the heap and we examine its label $l = {\cal L(\mathbf{x})}$ and its grey level $I(\mathbf{x})$;
• For the front points with label $l$, we know $N_l$ the number of points considered and $\mu_l$ the mean grey level value of those points;
• if $N_l < m$, we update this mean grey level value for the label $l$
  

  $$\mu_l = ( N_l * \mu_l + I(\mathbf{x}) )/ ( N_l + 1 ) \hbox{ , } N_l++$$ (5.10)

  
  

• considering that the image grey level of the pixel in the desired region are a random distribution of mean $\mu_l$ we can also update the variance with Koenig-Huyghens formula $\sigma_l^2[X] = \mu_l[X^2] - \mu_l[X]^2$ for the random variable $X$;
• we use this new mean values $\mu_l$ and $\sigma_l$ in the computations of the action at each neighbors of $\mathbf{x}$.

In figure 5.11 is shown an example of using this strategy for setting the penalty term.

Figure 5.11: Front competition on a 2D scanner slice of the lungs: The seed points were manually located on the dataset on the left image; they are visible as white dots on the right image.

The speed function in the image shown is computed with $m = 10$ for all labels. The speed function for each label $l$ is computed with the penalty ${\cal P}(\mathbf{x}) = e^{-\vert I(\mathbf{x}) - \mu_l\vert /2 \sigma_l^2}$.

Unfortunately, this algorithmic trick makes the penalty a time-dependent function, and thus the minimality principle is violated, and discretization cannot be applied in theory during the iterations where ${\cal P} = {\cal P}(\mathbf{x},t)$. Therefore, this heuristic must be seen as a region growing algorithm to operate a statistical study at the beginning of the segmentation process before visiting the whole image domain. It can also be replaced by a similar statistical study in a sphere around the seed points.

Figure 5.12 displays the result of the same strategy on the brain dataset used for test in figure 4.10.

Figure 5.12: Front competition result on the brain: Simultaneous segmentation of the white and grey matter of the brain, plus the skull and the skin; each row corresponds to a region, and the intersection between the regions and the datasets is represented in pink, in the left column.

This result was obtained by setting different seed points in the different objects in the brain MR image which is represented on the left column of figure 5.12. Each one of the regions was initialized with one marker, and they are represented when all fronts collide super-imposed on three orthogonal views of the dataset in the left column of figure 5.12. Inconvenient of this kind of segmentation is that it assumed that the correct number of classes is known a priori and that the user can clearly indicate to the algorithm the location of each connected components of this class. But target here is to provide a quick and dirty initialization to a more complex and time consuming method (like ).

The same test was done on a lung image, where we try to separate the airways from the soft tissue and vessels, in figure 5.13. Results

Figure 5.13: Front competition result on the lungs: separating air and soft tissues in a multi-slice CT dataset of the Lungs; seed points are manually located in the different parts of the image (on in the airways, on in the tissues, and on outside the object).

The underlying dataset is a volume-of-interest extracted from a multi-slice CT scanner of the lungs. The dataset, almost isotropic, contains a huge number of pixels, and even a fast algorithm like ours is not able to achieve this result in ``interactive time''. But we have shown the ability of this simple algorithm to achieve difficult segmentation tasks.

Further, in this algorithm, there is no need to use several data-structures (like min-heap) to store each different front. Due to the minimality principle of the construction, the point with minimum energy, independently from its label, will still be at the root of the min-heap data-structure. And each point of the image will be visited only one time, leading to approximately the same computing cost than the propagation of a single front in the image domain. It seems that a similar version of the multiple front propagation algorithm has been developed in [103,167], for initializing color and texture segmentation with the methods.