A Real Time Morphological Snakes Algorithm

This workshop may eventually be submitted and published.

Overview
References
Algorithm
Implementation details
Examples
Online Demo
Source Code

Luis Alvarez lalvarez@dis.ulpgc.es
Luis Baumela lbaumela@fi.upm.es
Pablo Marquez-Neila pmarquezneila@gmail.com
Pedro Henriquez phenriquez@ctim.es

Edited by Agustin Salgado asalgado@dis.ulpgc.es

Overview

Abstract

Active contours, or snakes, are computer-generated curves that move within images to find out salient image structures like object boundaries. Energy based formulations using a level set approach have been successfully used to model the snake evolution. The Euler-Lagrange equation associated to such energies yields to partial differential equations (PDE) which are usually solved using level set methods which involve contour distance function estimation and standard methods to discretize the PDE. Recently in reference [1] we have proposed a morphological approach to snake evolution. First, we observe that the differential operators used in the standard PDE snake models can be approached using morphological operations. By combining the morphological operators associated to the PDE components we achieve a new morphological approach to the PDE snakes evolution. This new approach is based on numerical methods which are very simple and fast. Moreover, since the level set is just a binary piecewise constant function, this approach does not require to estimate a contour distance function to define the level set. In this paper we present an algorithm to compute morphological snake evolution in real time.

Short presentation

Morphological snakes are based on a morphological discretization of the PDE model proposed in reference [2] and given by the equation:

$\frac{\partial u}{\partial t}= g(I)|\nabla u|\left( \mathrm{div} \left( \frac{\nabla u}{|\nabla u|} \right)+ \nu \right)+ \nabla g(I)\nabla u \qquad \qquad \qquad (1)$

We will approach the solution of the above equation by decomposing it in 3 terms : (i) the ballon force term, (ii) the smoothing term and (iii) the attraction force term. Next we will analyse these terms separately

The ballon force term

the balloon force operator term is given by the equation

$\frac{\partial u}{\partial t} = g(I) \nu |\nabla{}u|. \qquad \qquad \qquad \qquad \qquad (2)$

The above equation is strongly related to the classical erosion and dilation morphological operators. In fact, according with the sign of $\nu$ the solution of the above equation is equivalent to an erosion or dilation on the level set images and can be approached using the numerical scheme

$u^{n+1}(x)= \begin{cases} D_{d}u^{n}(x) & \text{if $g(I)(x) \geq t_1 $ and $ \nu >0$} \\ E_{d} u^{n}(x) & \text{if $g(I)(x) \geq t_1 $ and $ \nu <0$} \\ u^{n}(x) & \text{otherwise} \end{cases} \qquad (3)$

where D_d and E_d are the discrete versions of dilation and erosion. A disk with radius 1 and formed with eight neighbours and the pixel, is used as a structure element. Erosions and dilations are implemented in a very simple way by iterations of 8 or 5 neighborhood maxima (or minima) computation. In the demo is used the 5 neighbors version, and at the boundaries we use an homogeneous Neumann type boundary condition. Another useful option to make evolves this ballon type operator is to use an image interval value, that is :

$u^{n+1}(x)= \begin{cases} D_{d}u^{n}(x) & \text{if $I(x) \in [I_0,I_1] $ and $\nu >0$} \\ E_{d} u^{n}(x) & \text{if $I(x) \in [I_0,I_1]$ and $\nu <0$} \\ u^{n}(x) & \text{otherwise} \end{cases} \qquad (4)$

The smoothing term

Next we consider the regularization term given by the equation

$\frac{\partial u}{\partial t}= g(I)|\nabla u|\left( \mathrm{div} \left( \frac{\nabla u}{|\nabla u|} \right) \right) \qquad \qquad \qquad (5)$

The solution of the above PDE can be approached using line morphological operators. Let B the set of all line segments of length 2 centered at the origin . We define the morphological line operators as

$(SI_d u)(x) = \sup_{S \in B} \inf_{y \in x+hS} u(y)$

$(IS_d u)(x) = \inf_{S \in B} \sup_{y \in x+hS} u(y)$

Using these line morphological operators the solution of the above PDE can be approached using the numerical scheme (approximates mean curvature motion)

$u^{n+1}(x)= \begin{cases} (SI_{d} \circ IS_{d}u^{n})(x) & \text{if} \quad g(I)(x) \geq t_2 \\ u^{n}(x) & \text{otherwise} \end{cases} \qquad (6)$

where $SI_{d} \text{ and } IS_{d}$ are discrete versions of the above morphological continuous line operators. Both $SI_{d} \text{ and } IS_{d}$ have their own version of the set $\mathcal{B}, \mathcal{P}$ , which is a collection of four discretized segments centered at the origin:

$\mathcal{P}=\left( \begin{array} [c]{c}% P_{0}=\{(0,0),(1,0),(-1,0)\}\\ P_{1}=\{(0,0),(0,1),(0,-1)\}\\ P_{2}=\{(0,0),(1,1),(-1,-1)\}\\ P_{3}=\{(0,0),(1,-1),(-1,1)\} \end{array} \right)$

Below there are some examples of the effect of the $SI_{d} \text{ and } IS_{d}$ operator on individual pixels of binary images. In those cases where a straight line is found (marked in red), the central pixel remains active (a) and (b). When the central pixel does not belong to a straight line of active pixels, it is made inactive (c) and (d). For exemplification purposes, we assume the pixels on the borders are not affected by the operator.

Examples of the $SI_{h}$ operator

Example of the $SI_{h} \circ IS_{h}$ operator iterated until convergence.

The attraction force term

The attraction force PDE term is given by the equation

$\frac{\partial u}{\partial t}=\nabla g(I)\nabla u \qquad \qquad \qquad \qquad (7)$

The morphological discretization of this PDE term is very simple and given by

$u^{n+1}(x)= \begin{cases} 1 & \text{if } \nabla u^{n}(x)\nabla g(I)(x)>0 \text{ and } g(I)(x) \geq t_3 \\ 0 & \text{if } \nabla u^{n}(x)\nabla g(I)(x)<0 \text{ and } g(I)(x) \geq t_3 \\ u^{n}(x) & \text{otherwise}% \end{cases} \qquad (8)$ Gradient is implemented using basic finite difference schemes, we use central finite difference scheme. We observe that the image level set u(x) we use is a binary function taking values in {0,1} and the snake is given by the level set contour.

References

M. Aleman-Flores, L. Alvarez and P. Henriquez. AmiLabContours: A tool for image structure segmentation, To appear in IADIS Multi Conference on Computer Science and Information Systems 2010. http://ami.dis.ulpgc.es/ami/amilabcontours.html
L. Alvarez, L. Baumela, P. Marquez-Neila and P. Henriquez. Morphological Snakes. CVPR 2010 http://www.sciweavers.org/publications/morphological-snakes
V. Caselles, R. Kimmel, and G. Sapiro. Geodesic active contours. International Journal of Computer Vision, 22(1):61–79, 1997. http://www.springerlink.com/content/w63mp82m702nk816/
F. Catté, F. Dibos, and G. Koepfler. A morphological scheme for mean curvature motion and applications to anisotropic diffusion and motion of level sets. SIAM Journal on Numerical Analysis, 32(6):1895–1909, 1995. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.49.8784&rep=rep1&type=pdf
F. Guichard, J.Morel, and R.Ryan. Contrast invariant image analysis and PDE’s. http://mw.cmla.ens-cachan.fr/~morel/JMMBookOct04.pdf.

Algorithm

To design the algorithm for the snake evolution we take into account that the solution of the snake evolution PDE can be approached by composition of the solution of the 3 terms the original equation has been decomposed, so the algorithm is very simple, it is just iterations of the composition of the numerical schemes (3) or (4),(6) and (8). One important observation is that these 3 numerical schemes are morphological in the sense that they do not create extra level set values, that is if the input level set is a binary image, the output level sel is also a binary image.

Parameter selection

1.- Snake initialization

The snake is initiallized using a polygon provided by the user. From this polygon we build a binary level set image u(x) assigning the value 1 inside the contour and 0 outside the contour.

2.- Attraction force g(I). According with the application requirements we can choose between

$\begin{array}{ll} g(I)=\dfrac{1}{\sqrt{10^{-20}+|\nabla G_{s }\ast I|}} & \text{Edge border attractor} \medskip \\ g(I)=|G_{s}\ast I| & \text{Dark center line attractor} \medskip \\ g(I)=1 & \text{No attractor defined}% \end{array}$

where G_s is the Gaussian function with standard deviation s.

3.- Choice of t₁, t₂ and t₃ in expressions (3),(5) and (7).

To simplify the parameter choice t₂ and t₃ are fitted to 0, that is regularization and snake attractor are always applied. t₁ (when applied) is estimated using a percentage of g(I) histogram values. That is, the user provided a number $0\leq p \leq 1$ and we take t₁ such that

$\frac{\left\vert \left\{ x\in \Omega :g(I)(x)\geq t_{1}\right\} \right\vert }{\left\vert \Omega \right\vert }=p \qquad \qquad (9)$

where $\Omega$ is the image domain. In the case ballon option (4) is used, we use as parameter the radius of interval [I₀,I₁]. The interval center is obtained using the median of the image value in a contour neighborhood region.

4.- Maximun Number of iterations

Since we deal with an iterative procedure user also provides the maximun number of iterations

Algorithm Output

1.- Contour collection

The main algorithm output is a contour collection with the contour points for all algorithm iterations

2.- Illustrative image

To illustrate the results we draw in the original image the contour for some algorithm iterations

Implementation details

As explained above, the snake evolution is performed by iterations of the composition of the numerical schemes (3) or (4), (6) and (8). To speed up the algorithm we store in memory the binary level set snake representation and a contour representation where only the snake contour points are stored. Using this representation numerical scheme iterations are performed in the following way :

step 1 : Initialitation of contour points and binary snake level set from the contour provided by the user
step 2 : In each iteration of operators (3) (4), (6) or (8):
- a) For each point of the contour we compute the operator in a 3x3 neigborhood of the point to update the binary snake level set.
- b) We update the contour using the new updated binary snake level set obtained above

Using this hybrid scheme where we manage the contour representation of the snake and the associated level set we obtain that at each iteration, only a 3x3 neigborhood of the contour is considered for computation.

We observe that the snake evolution is always performed in the level set representation, never in the contour representation. On the other hand contour update from snake level set function is very easy : x is part of the contour line if u(x) = 1 and at least one of its 4 neighbors is equal to 0. For snake evolution purpose we do not need any special contour organization. Contour is a non-ordered collection of 2D points.

The linear differential operators we use, as the gradient, are implemented using simple centered finite difference scheme

We use always homogeneous Neumann type boundary conditions in operator implementation.

Examples

We present some examples to illustrate the perfomance of snake algorithm.

First, we present 2 medical images. The initial snake is in blue inside the object. Using different colours we draw different iterations of snake algorithm.

Next we present the results of snake algorithm in the case of a electronic microscopy image

Online Demo

An on-line demo of the algorithm is available. The demo allows:

(i) Upload any image and to segment an image structure.

(ii) Fit the main snake evolution parameters which are (we strongly simplify the parameter choice) :

The interest structure type. We have 2 options: (a) HIGHCONTRASTREGION : a contrasted region in the image or (b) DARKCENTERLINES : centerlines of dark band structures.
Snake Ballon force type. We have 3 options: (a) NOBALLONFORCE : No Ballon force is used, (b) EXPANSIONBALLONFORCE an expansion Ballon force or (c) CONTRACTIONBALLONFORCE : a contraction Ballon force.
The number of iterations.
The standard deviation of the Gaussian to fit the scale to look for interest structures (used in g(I) snake operator function). The mask size is showed with a circle to give the user more information about the changes produced by this parameter.

(iii) Define the initial snake by clicking in the image.

(iv) Run the snake algorithm. The output of snake algorithm is an ASCII file with the final contour and an illustrative image with some iterations of the snake algorithm.

(v) Run the snake algorithm using a previous output as new input.

Source Code

Source Code is avalaible in MorphologicalSnakes 20111228.zip

Source Code documentation is avalaible in Readme.pdf

A Basic Source Code where only .bmp image format is used is avalaible in MorphologicalSnakes basic 20111228.zip, and its documentation is available online.

image credits

courtesy Servicio de Radiodiagnóstico del Hospital Universitario Insular de Gran Canaria

courtesy Center for Biomedical Technology, Universidad Politécnica de Madrid