Fast probabilistic snake algorithm

Abstract

Few people use the probability theory in order to achieve image segmentation with snake models. In this article¹¹1International Conference on Image Processing (ICIP), Barcelona, Spain, September 2003, we are presenting an active contour algorithm based on a probability approach inspired by A. Blake work and P. Réfrégier’s team research in France. Our algorithm, both very fast and highly accurate as far as contour description is concerned, is easily adaptable to any specific application.

Keywords: snake, energy minimization, probability, regularization.

1 Introduction

In the last fifteen years the active contours have been successfully applied in many different ways. Snake can be modelled according to several mathematical formulation. Each one of them entails its own drawback but provides its own advantages at the same time. In 1988, Kass et al. [1] proposed an energy-based formulation,

\begin{split}E(C)=\int_{0}^{1}\alpha(s)\left\|\dfrac{\partial C(s)}{\partial s}\right\|^{2}+\beta(s)\left\|\dfrac{\partial^{2}C(s)}{\partial s^{2}}\right\|^{2}\\ -\left\|\nabla I(C(s))\right\|^{2}ds\end{split}

(1)

in which $C$ represented the curve, $I$ the image, $\alpha$ and $\beta$ the parameters used for the control of the curve properties. Though it works, Kass’s model is limited in its sensitivity to initialization. What is more, it has no natural adaptability to changes in topology. As for the improvements proposed (see [2],[3]), they add to the complexity of the implementation.

In 1997 geodesic active contours were proposed as a geometric alternative to snake (see [4],[5]). They included as well an energy-based term. Yet they managed to overcome the deficiencies of the classical snake formulation. Thanks to geodesic active region [6] and to levelset implementation, geodesic active contours stopped being sensitive to initialization and at the same time topology changes became implicit.

Lately the question has been examined from an innovative angle by A. Black et al. [7] who resort to probabilities and B-Spline curves. And as for him C. Chesnaud [8] proposes a region-based criterion so as to minimize,

J(s,\vec{C})=N_{a}(\vec{C})H(\widehat{\theta}_{a})+N_{b}(\vec{C})H(\widehat{\theta}_{b})+K_{l}

(2)

In which $N_{i}(\vec{C})$ is the number of pixels inside the curve and $H_{i}(\widehat{\theta_{i}})$ represents a measure of probability wich depends on the a priori model chosen to figure the pixels distribution, see [8] for more details (we shall call this model CASP).

CASP is actually interesting but its main drawback is that it cannot be used with open curves. That is why we have endeavoured to develop a new formulation of CASP so as to adapt it equally to open or closed curves. The following section gives the formulation of our model, the third one exposes the implementation of our algorithm and the last one states our results and compares them with other models.

2 Statistical snake model

Let us assume that:

•

$I$ is the image
•

$C(s)$ represents the curve wich corresponds to $N$ knots ( $N_{i}$ )
•

$B$ represents the edges in the image
•

$p(C|B)$ is the probability density that the snake is on the edge
•

$p(B|C)$ is the probability density that there is an edge on the current position of the snake

According to the Bayes’ rule we can write,

p(C|B)=\dfrac{p(B|C)p(C)}{p(B)}

(3)

Since $p(B)$ is constant, the final snake position is the curve $C$ which minimizes the likelihood $p(C|B)$ , or
$p(B|C)p(C)$ .
If we consider that each knot is independent we can write,

p(B|C)=\prod_{i\in[0..N-1]}p(N_{i}|I)

(4)

The choice of the density $p(N_{i}|I)$ depends on the problem we want to solve. For example, if we want to detect the edges in the image, we can calculate the variance inside a window centered in $N_{i}$ (see next section).
The density $p(C)$ represents the constraints on the curve such as [8] where,

p(C)=A\exp\{-\frac{1}{2\varphi^{2}}\sum_{i=0}^{N-1}d_{i}^{2}\}

(5)

in which $d_{i}$ is the distance between $N_{i}$ and the center of segment $N_{i-1}N_{i+1}$ . It limits the irregularity of the curve. Another example of density will be given in the next section.
The upper model has a very simple formulation and we can easily adapt it to different problems. A drawback is that the topology changes are not included in the formulation.

3 Implementation of the statistical snake

In this section, we will give some practical pieces of information to implement our model. As in the case of geodesic snake and, as in [7], we consider only the movement along the normal of each knot of the curve. Along this normal, we fix a research depth (named $L$ ), in order to find the new position of the current knot. On each point of this segment we calculate a measure of variance inside a window so as to estimate $p(N_{i}|I)$ . Thus an edge corresponds to a maximum of the variance (see Fig.1).

Refer to caption — Figure 1: Ideal edge image and variance along the normal of $N_{i}$

If we intend to regularize the snake, we can use a function which directly acts on the density function $p(N_{i}|I)$ . For example if we want to give top priority to the nearest contour to the current position, we can filter $p(N_{i}|I)$ with a “low-pass filter” considering that “low frequencies” corresponds to the positions close to the current knot and “high frequencies” to the distant positions (see Fig.2).

With open curves, the normal chosen at end points is orthogonal to the last segment.

So as to improve the convergence of the snake, we process three steps:

•

Optimizing the curve with a few points and high research depth.
•

Resampling the snake to have high resolution contour.
•

Re-optimizing the curve with all points but with a little research depth.

4 Experimental results

Here we propose to compare first our model with the classical snake model and the CASP model of [8] in the case of closed snakes. Then we will compare our model with the classical one in the case of open snakes.

4.1 Closed curves

Figure 3 shows the results in the closed curve case. In (a) we have the initialization, in (b) we can see classical snake result with the parameters: $\alpha=0.1$ , $\beta=0.1$ , $\lambda=0.5$ where $\lambda$ is the regularisation parameter in the evolution equation (see [9]). In (c) we can see the output of the CASP algorithm when we choose max $\_$ deviation $=5$ , Nb $\_$ iteration $=3000$ , regularization $=0.2$ and a gaussian density for each region. The last, (d), shows the result with our model with:

•

First pass: research depth = 25, window analysis size of 7 pixels, regularisation = 0
•

Resampling with a maximum of 4 pixels between two consecutive knots
•

Second pass: research depth = 5, same window size and regularisation = 1

Our model gives better resolution than the classical snakes. It is as efficient as CASP but our parameters can be found more easily than any.

Tab.1 gives the different computing times for the two algorithms. The tests were made on a bi-P3-866 with 1Gb of RAM. We do not give the computing time for the CASP model because we haven’t implemented the fastest algorithm described in [8].

Algorithm	Time
Classical snake	566 $\mu$ s
Our model	9.95 $\mu$ s

Table 1: Time convergence for closed curves

Fig.4 gives another example of segmentation in infrared image for military application. We use the model with only 14 knots for the first step, resampling with a maximum of 4 pixels between two consecutive knots in the second step and a new optimizing pass in the third step. As we can see, the tank is extracted with high accuracy.

4.2 Open curves

Fig.5 highlights the results obtained on open curves. Initialization is shown in (a). (b) gives the results of the classical snake using $\alpha=0.1$ , $\beta=0.5$ , $\lambda=0.99$ . (c) depicts our model results using the following steps:

•

First pass: research depth = 20, window analysis size of 7 pixels, regularisation = 0
•

Resampling with a maximum of 4 pixels between two consecutives knots
•

Second pass: research depth = 5, same window size and regularisation = 0

The parameters that control the classical snake were found with difficulty. Tab.2 gives the different computing times showing that our algorithm if 200 times faster than the classical snake model.

Algorithm	Time
Classical snake	565 $\mu$ s
Our model	2.2 $\mu$ s

Table 2: Time convergence for open curve

5 Conclusion

We have proposed a very fast snake algorithm based on a probability approach. This model accuracy is quite high when using two steps that speed the convergence using few knots and then resample the curve. Futhermore this model uses only three parameters, one of them (the window size) being set to a fixed value for most applications.
We are now working on its adaptability to fit multiple objects detection (topology changes) and on the regularization process to increase its stability.

References

[1] M.Kass, A.Witkin, and D.Terzopoulos, “Snakes: actives contours models,” International Journal of Computer Vision, vol. 1, pp. 312–331, 1988.
[2] L.D.Cohen, “On active contour models and balloons,” CVGIP: Image understanding, vol. 53, no. 2, pp. 211–218, 1991.
[3] Y.Elomary and J-M.Chassery, “Cooperation between active contours and multiresolution for edge-based segmentation,” Proceedings of the IEEE IMDSP Conference, Cannes, France, pp. 206–207, 1993.
[4] V.Caselles, R.Kimmel, and G.Sapiro, “Geodesic active contours,” International Journal of Computer Vision, vol. 22, pp. 61–79, 1997.
[5] S.Kichenassamy, A.Kumar, P.Olver, A.Tannenbaum, and A.Yezzi, “Gradient flows and geometric active contour models,” IEEE International Conference on Computer Vision, vol. 1, pp. 810–815, 1995.
[6] N.K.Paragios, Geodesic active regions and level set methods: contributions and applications in artificial vision, Ph.D. thesis, University of Nice Sophia Antipolis, France, 2000.
[7] A.Blake and M.Isard, Active Contours, Springer, 1999.
[8] C.Chesnaud, Techniques statistiques de segmentation par contour actif et mise en œuvre rapide, Ph.D. thesis, ENSPM, University of Aix-Marseille, France, 2000.
[9] Y.Elomary, Modèles déformables et multirésolution pour la détection de contours en traitement d’images, Ph.D. thesis, University of Joseph Fourier, Grenoble, France, 1994.