## A simple implementation of snakes

In this post I want to show how easy it is to implement snakes using MATLAB and DIPimage. The next release of DIPimage will have some functionality implementing snakes, and this post is based on code written for that. The first implementation of snakes I made for DIPimage used a custom class, `dip_snake`

, to simplify usage. We decided eventually to do without the class. The code discussed here uses that first implementation, which you can download here. The most important bits of the code are identical to that in the next release of DIPimage, but the wrapping is not. Adding the class is a nice excuse to show how to modify DIPimage to understand a new parameter type. I will discuss that in a future post.

### Snake theory

A snake (a.k.a. active contour) is a flexible 2D line, which is moved around the image to minimize an energy functional. The line is parametrized by through a variable *p* that goes from 0 to 1. Two functions of this variable, *x*(*p*) and *y*(*p*), define the coordinates of the points along the line. For simplicity we use the vector * s*(

*p*) = (

*x*(

*p*),

*y*(

*p*))

^{T}. More often than not, the snake is implemented as a closed contour, because it is most often used to segment an object. This is enforced setting

*(0) =*

__s__*(1).*

__s__The energy functional that we want to minimize is defined as:

E = ∫ E_{int}(* s*(

*p*)) + E

_{ext}(

*(*

__s__*p*)) d

*p*

E_{int} is an internal energy, meant to force the snake to be small and smooth. It avoids obviously wrong solutions. E_{ext} is an external energy, the one responsible for the snake finding the edges of an object in the image. E_{ext} is therefore defined on the image. A common external energy is the inverse of the gradient magnitude: low energies at the location of the edges, higher energies everywhere else. This energy definition comes from the original paper by Kass, Witkin and Terzopoulos (1988).

Let’s discuss the internal energy first. It is defined as

E_{int}(* s*(

*p*)) = ½ {

*α*(

*p*)|

*′(*

__s__*p*)|

^{2}+

*β*(

*p*)|

*″(*

__s__*p*)|

^{2}}

This looks rather complicated, but it need not be. For starters, *α* and *β* are usually constant (the method has enough parameters as it is, there is no point in making two of the parameters different for every point along the snake!). So what we have left is the magnitude of the first derivative, which is larger for longer snakes, and the magnitude of the second derivative, which is larger for sharper bends. The first part keeps the snake short, the second part keeps it straight. The two parameters, *α* and *β*, define the relative importance of these two terms.

To minimize this energy function we use the Euler-Lagrange equation. After quite some work we end up with the following equation:

*α** s*″(

*p*) –

*β*

*″″(*

__s__*p*) – ∇E

_{ext}(

*(*

__s__*p*)) = 0

(By the way, when you derive this is when you understand why the original definition of the internal energy has that ½ in it!) ∇E_{ext}, the gradient of the external energy, is a force. We call this the external force, __F___{ext}, and this is where most of the improvements to snakes were made. I won’t discuss them here.

Solving the equation above is accomplished with the gradient descent method: you convert the snake * s* into a function of time, and replace the 0 with the partial derivative of

*to time. You can see the reason for this if you assume that, after a long time when the snake has converged to a minimum, it’s derivative to time will be zero, and so the equation above is satisfied.*

__s__∂* s*(

*p*) / ∂t =

*α*

*″(*

__s__*p*,

*t*) –

*β*

*″″(*

__s__*p*,

*t*) +

__F__

_{ext}(

*(*

__s__*p*,

*t*))

Next we need to discretize the parameter *p* so we can actually implement this numerically. So we set *p* = *i*/*N*, with *i* going from 0 to *N*-1. Because the snake is periodic, * s*[

*N*] =

*[*

__s__*0*],

*[*

__s__*N+1*] =

*[*

__s__*1*], etc. (Note I’m using the square brackets to denote the discrete indexing, as opposed to the round brackets for continuous-domain functions.) This is also the time when we split the vector

*into its two components*

__s__*x*and

*y*. The two components of the external force

__F__

_{ext}are

*f*and

_{x}*f*. We use the finite difference approximation to the spatial derivatives:

_{y}*x*″[*i*] = *x*[*i*-1]-2*x*[*i*]+*x*[*i*+1]

and

*x*″″[*i*] = *x*[*i*-2]-4*x*[*i*-1]+6*x*[*i*]-4*x*[*i*+1]+*x*[*i*+2].

This results in the following set of 2*N* equations:

∂*x _{t}*[

*i*] / ∂t =

*α*(

*x*[

_{t}*i*-1]-2

*x*[

_{t}*i*]+

*x*[

_{t}*i*+1]) +

*β*(-

*x*[

_{t}*i*-2]+4

*x*[

_{t}*i*-1]-6

*x*[

_{t}*i*]+4

*x*[

_{t}*i*+1]-

*x*[

_{t}*i*+2]) +

*f*(

_{x}*x*[

_{t}*i*],

*y*[

_{t}*i*])

∂*y _{t}*[

*i*] / ∂t =

*α*(

*y*[

_{t}*i*-1]-2

*y*[

_{t}*i*]+

*y*[

_{t}*i*+1]) +

*β*(-

*y*[

_{t}*i*-2]+4

*y*[

_{t}*i*-1]-6

*y*[

_{t}*i*]+4

*y*[

_{t}*i*+1]-

*y*[

_{t}*i*+2]) +

*f*(

_{y}*x*[

_{t}*i*],

*y*[

_{t}*i*])

We can write the whole thing, except the external force, as a *N*-by-*N* matrix *A*:

-2α-6β |
α+4β |
–β |
0 | 0 | 0 | ⋯ | –β |
α+4β |

α+4β |
-2α-6β |
α+4β |
–β |
0 | 0 | ⋯ | 0 | –β |

–β |
α+4β |
-2α-6β |
α+4β |
–β |
0 | ⋯ | 0 | 0 |

0 | –β |
α+4β |
-2α-6β |
α+4β |
–β |
⋯ | 0 | 0 |

0 | 0 | –β |
α+4β |
-2α-6β |
α+4β |
⋯ | 0 | 0 |

⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋱ | ⋮ | ⋮ |

–β |
0 | 0 | 0 | 0 | 0 | ⋯ | -2α-6β |
α+4β |

α+4β |
–β |
0 | 0 | 0 | 0 | ⋯ | α+4β |
-2α-6β |

such that:

∂*X _{t}* / ∂t =

*AX*+

_{t}*f*(

_{x}*X*,

_{t}*Y*)

_{t}∂*Y _{t}* / ∂t =

*AY*+

_{t}*f*(

_{y}*X*,

_{t}*Y*)

_{t}Note the simple structure of this matrix. Only 5 diagonals have values in them, the rest of the matrix is 0, no matter how large *N* is. This is called a cyclic pentadiagonal banded matrix, and can easily be inverted using Cholesky decomposition. If the snake were not circular (the two ends not connected) this matrix would be even simpler!

The last step is to discretize the time parameter, adding a step size *γ*:

∂*X _{t}* / ∂t = (

*X*–

_{t}*X*) /

_{t-1}*γ*.

To make this easy to solve, we need to do a little trick: we will assume that the movement of the snake is small enough that the external force doesn’t change much from one time step to the next. We can then substitute in the set of equations above the external force at time *t*-1, yielding a simple to solve set of equations:

*X _{t}* = (

*I*–

*γ*

*A*)

^{-1}{

*X*

_{t-1}+

*γ*

*f*(

_{x}*X*

_{t-1},

*Y*

_{t-1})}

*Y _{t}* = (

*I*–

*γ*

*A*)

^{-1}{

*Y*

_{t-1}+

*γ*

*f*(

_{y}*X*

_{t-1},

*Y*

_{t-1})}

### Snake implementation

Finally! We’re ready to start writing code…

Let’s start with creating an image to test our algorithm on:

img = noise(50+100*gaussf(rr>85,2),'gaussian',20)

We will also need an initial snake to deform:

x = 120+50*cos(0:0.1:2*pi)'; y = 140+60*sin(0:0.1:2*pi)'; hold on, plot([x;x(1)],[y;y(1)],'g')

Note that I duplicated the first point at the end when plotting, so that the snake would be a closed contour. I’ll define the parameters now. This is a tricky part, I played around with them until I was satisfied. I haven’t yet seen a clever way of determining these. Because there is no weight for the external force, the *γ* parameter needs to be increased if the external force is too weak. To compensate, both *α* and *β* should be decreased by the same factor. Decreasing *β* makes the snake less smooth. Increasing *α* increases the force that wants to shrink the snake; because we started inside the boundary, we want it very small.

alpha = 0.001; beta = 0.4; gamma = 100; iterations = 50;

Next we will compute the matrix inversion. We need to create the matrix *I* – *γ**A* and invert it. This is the most complex part of the whole algorithm, but in MATLAB it’s a breeze:

N = length(x); a = gamma*(2*alpha+6*beta)+1; b = gamma*(-alpha-4*beta); c = gamma*beta; P = diag(repmat(a,1,N)); P = P + diag(repmat(b,1,N-1), 1) + diag( b, -N+1); P = P + diag(repmat(b,1,N-1),-1) + diag( b, N-1); P = P + diag(repmat(c,1,N-2), 2) + diag([c,c],-N+2); P = P + diag(repmat(c,1,N-2),-2) + diag([c,c], N-2); P = inv(P);

`a`

, `b`

and `c`

are the elements along the 5 diagonals that are not 0. I used the function `diag`

to populate the matrix. `inv`

takes care of selecting an appropriate inversion algorithm. Nothing to worry about here. Note that we need to do this inversion only once, not at every iteration. The matrix is independent of the position of the snake.

Next up is the external force. We need to compute the gradient of the image representing the external energy. This is typically the gradient magnitude. Using DIPimage this is again a breeze! We’ll use `gradmag`

and `gradient`

:

f = gradient(gradmag(img,30));

I’ve used a very large value for the parameter *σ* here. This is because the snake starts off far away from the edge we’re looking for, and the external force is the only thing that is able to pull it towards these edges. Many of the improvements to the snake so far have been therefore to increase the “pull” range of edges by modifying the external force. Since I don’t want to go into that, I’ve opted to use a large *σ*, which spreads out the gradient over a large area. Another thing to note here is that `f`

is a *vector image*, that is, each pixel has two values, one for the gradient along the x-axis, and one along the y-axis. `f{1}`

extracts a scalar image containing only the x-component of the gradient, whereas `f(100,100)`

extracts a vector corresponding to the pixel at (100,100).

The last thing to do is to get the value (at sub-pixel locations by interpolation) of the external force at the snake’s control points, and update the position of these control points according to the gradient descent equation we derived. We’ll do this in a loop, `iterations`

times, and plot the resulting snake every 5 iterations. I’m plotting the final result in red:

for ii = 1:iterations % Calculate external force coords = [x,y]; fex = get_subpixel(f{1},coords,'linear'); fey = get_subpixel(f{2},coords,'linear'); % Move control points x = P*(x+gamma*fex); y = P*(y+gamma*fey); if mod(ii,5)==0 plot([x;x(1)],[y;y(1)],'b') end end plot([x;x(1)],[y;y(1)],'r')

Ok. So now we’ve got a snake, and got it moving towards an edge! Too bad that the snake seems to converge to a point inside the edge. This is because the Gaussian gradient used to compute the external energy had a very large parameter *σ*. In a Gaussian smoothed image, edges always move inwards. We could run the loop above again, but using an external force computed with a smaller *σ*, since the snake is closer to the edge now. However, there are better alternatives to the external force that do not require this type of multi-scale approach. See for example gradient vector flow or vector field convolution. I won’t go into these, but the next release of DIPimage will have functions to compute these external forces.

Yet another modification to the snake is the so called balloon force. It is an internal force that moves the snake outwards (or inwards if you make it’s weight negative). It can very easily be added to the loop above, as you can see in the function `snakeminimize`

. Go ahead, download and install it into your DIPimage directory. There’s a README file with instructions inside. That function uses a class `dip_snake`

, which I will discuss in a future post.

Excellent tutorial!! Easily one of the best tutorials on snakes I have come across.

Thanks! Great tutorial! Now I undestand the theory of the snake and its computational implementation.

I’ve one question. What does the parameter σ means in gradmag?? How does it calculate the gradient of a region?? I didn’t find any documentation about this function…

Thanks again!

Lucas,

gradmagis the gradient magnitude. It uses Gaussian gradients, σ is the scale parameter. See this blog post for more information regarding Gaussian gradients.Hello Sir,

This tutorial is very useful and informative.

I am a PG student and I am doing my project in Medical Image Processing.(Lung Cancer Diagnosis)

So in order to extract the Lung Nodules(tissues) which is my ROI, Can i use this Snakes.

How should i approach my problem sir.

Regards,

Karthik.

Karthik,

It is difficult to give directions on how to accomplish such a task without seeing the images in question. In general for snakes, you need to initialise your snake close to the final solution, meaning you need some way of detecting the nodules

beforeyou apply the snake!Sir,

Thank you for this wonderful post. It has helped me to begin my coding of segmentation.

I faced an issue while implementing your algorithm.

In some cases, for a given iteration,the formulas

x = P*(x+gamma*fex); and

y = P*(y+gamma*fey);

give negative values. What could be the reason and how do we overcome it?

Hi Nathanael,

Yes, sometimes the snake “walks” out of the image. You need to check x and y after each iteration to make sure they’re still inside the image domain, and if not, move them back in. For example,

`x(x<0) = 0;`

x(x>max_x) = max_x;

y(y<0) = 0;

y(y>max_y) = max_y;

Hi first of all we want to say u thanks for ur article, it’s the first concrete implementation of snake that we have found on the net.

We are working at a projet for create a snake filter in silverlight, but we still have a doubt even after reading your article, and we would be gratefull to you if you can help us.

At the end of your article u said :

Xt = ((I – γA)^-1){Xt-1 + γ fx(Xt,Yt)}

Yt = ((I – γA)^-1){Yt-1 + γ fy(Xt,Yt)}

so if I is the matrix that contain image luminosity (suppose size n*m)

A is the matrix u explained, that we can expand to have the same size of the I.

Now we got that the other term of expression should have also n*m size because we have to multiply 2 matrix.

So what exactly Xt will contain?

suppose at the time 0 we have a snake with 3 point

(1,5) (3,9) (30,9)

and the size of the image is 30 x 30 pixel

what will be our X0?

Thank you

Hi,

The matrix

Iis the identity matrix (all elements are zero except the main diagonal, that is all ones).fand_{x}fare the two images that you need to create based on the input image: one is the x-derivative and the other the y-derivative of the external energy image._{y}In your example with a snake with 3 points (1,5) (3,9) (30,9),

X= (1,3,30)_{0}^{T}andY= (5,9,9)_{0}^{T}encode the initial coordinates for the snake.Hope this clarifies things!

hi I read ur very interesting article, and i would be grateful to u if u can answer my question.

Suppose we implement the snake algorithm in C when we have to calculate fx and fy as the x derivative and y derivative of the intensity of the picture how can we obtain a vector if we use the kernel matrix that need not only a point. Here is an example: we have the point (30,30) of

intensity 1 and we have to calculate fx so we use the matrix Gx =((1,0,-1),(2,0,-2),(1,0,-1)). As we need to do Gx*A to calculate fx we would need not only our point (30,30) but at least 3 point (29,30),(30,30),(31,30), in this case the result will not be a scalar but a vector. So in the end fx will be a matrix and not a vector as it will be a vector of vector and not a vector of scalar. Is there another way to proceed?

Thank u

Dott.Marika

Hi Dott.Marika,

In your example, fx = Gx*A is an image, same size as A. You need to take the value of this function at (30,30), that is the derivative at that point. Of course you use more pixels in the calculation, but the end result is a single value.

Thanks for the tutorial Prof. Luengo. I was wondering how matrix A would change if we had an open contour instead of a closed one. Also, periodicity of the snake as described above cannot be assumed in such a case. Will that affect the implementation in anyway ?

In response to Nathanael’s comment: If the pixel intensities in the image you are using varies from 0 to 255, for the alpha,beta and gamma values assumed in this tutorial, the initial contour goes beyond the image boundaries in the first iteration itself. Converting it to [0,1] solves the problem.

Dear Vikram,

The matrix

Awould have zeros also in the top right and bottom left corners. That is, it would be a simple band matrix (pentadiagonal matrix). This matrix is even easier to invert. The algorithm doesn’t change at all. However, you might need to include some constraint on the end points, so that the snake doesn’t converge to a single point. I have never implemented such a snake, but I suspect that the first and last rows of the matrix should be changed so that these points are not pulled inwards.Just wanting to clarify, for

fex = get_subpixel(f{1},coords,’linear’);

Is this just returning an intensity at the specified coordinates of the x-grad matrix? When you do,

x = P*(x+gamma*fex);

y = P*(y+gamma*fey);

It looks like you are reassigning the coordinate of the snake to be a matrix (P) multiplied by a constant (x + gamma*fcx), which will mean your coordinates of your matrix on the next pass through will then be matrices? What am I missing?

Thanks

James,

`x`

is a vector, as is`fex`

. Thus,`P*(x+gamma*fex)`

is a matrix times a vector, which results in a vector again. In effect, this statement computes the new coordinates for all snake control points in one operation.Hi Cris, great tutorial. I was just trying to implement this without the toolbox, I’ve replaced

f = gradient(gradmag(img,30));

with,

[FX, FY] = gradient(double(img));

H = fspecial(‘gaussian’, 30);

FY = imfilter(FY, H);

FX = imfilter(FX, H);

And it works fine, I tried replacing

fex = get_subpixel(f{1},coords,’linear’);

fey = get_subpixel(f{2},coords,’linear’);

with,

F = griddedInterpolant(FX, ‘linear’);

for i = 1:length(x)

fex=F(x,y);

end

G = griddedInterpolant(FY, ‘linear’);

for i = 1:length(y)

fey=G(x,y);

end

etc, and it does something very different.. Is my inderstanding of the sub_pixel function incorrect?

Thanks!

Hi Sam,

I think the difference might be in the indexing: for

`dip_image`

objects, the first index goes horizontally, and start at 0. You probably need to switch the x and y coordinates, and might need to add 1 (depending on how you defined the snake earlier).Dear Cris

It’s a long time I’m searching to find a tutorial on snake implementation, here I find it perfect.

No doubt your weblog is the best tutorial and analysis web I have ever seen.

I surely wish you best…

Hello,

Your explanations is clear.

I wrote the algorithm in ITK. I share my implementation if someone finds it useful.

I only skipped the sub-pixel sampling.

One can invoke like this:

Best wishes for you!

[code in comment moved to these files: CMakeLists.txt and Snakes.cxx. That they’re hosted here does not mean I’ve tested or vetted them. — Cris]

I don’t understand how the snake is getting larger if it’s an energy minimizing problem. Does this have to do with the parameters you chose? I thought a growing snake was the balloon snake method as opposed to the original snake method which starts outside the boundary.

Thanks!

Michelle

Hi Michelle,

If the larger snake has a lower energy, then the snake will grow. The energy depends on image features and snake shape. One of the internal energy terms is smaller for a shorter boundary, but the other internal term is smaller for a less curved boundary. A shorter boundary also needs a stronger curvature, and therefore, these two internal energy terms oppose each other in some way. What then makes the snake grow or shrink is the external energy (i.e. the location of the boundary) or the balloon term. In the example above, it is the boundary “pulling” the snake towards itself.

hi, thank you so much for your tutorial.

I have a silly question about the derivative of the external force. If you can read the paper “Real-time Non-Rigid Surface Detection” by Julien Pilet, Vincent Lepetit, Pascal Fua, in this paper, they used the snake model to solve the energy minimization problem for surface detection. the formula were quite similar. however, the external force they used was not the usual image gradient but the sum of the square distances between corresponding feature points in two images.

In this case, I ‘m not sure how to calculate the fx(Xt,Yt) or fy(Xt,Yt) in that paper, because I thought the f in the paper is not a vector but a single value, so do you have any thought?

Sorry about my English…(not mother tongue…)

Thank you!

Hi Xiaofei,

I have not read the paper you refer to, but I presume that the “sum of square distances between corresponding feature points in two images” is the external energy (let’s call it E). To get the external force, take the derivative along x and y: fx = dx(E), fy = dy(E).

Dear Prof Chris,

Good day!

Thank you very much for the algorithm. It is very helpful for an amateur like me.

I am perusing my Masters degree in Image processing and using your algo to help me understand how snake algorithm works. Could you please help me in solving this error?

I am getting an error at

>> fex = subpixel(f(1),coords,’linear’);

??? Undefined function or method ‘subpixel’ for input arguments of type ‘double’.

Error in ==> snakemethod at 34

fex = subpixel(f(1),coords,’linear’);

Hi Swapnil,

Did you install DIPimage?

Dear Chris,

Good day!

Thank you very much. It was very much valuable and helpful. I have installed DIPimage toolbox 2.1 and its working fine.

I may require your valuable resource further as well.

Thanks & Regards,

swapnil R kulkarni

Pune, India.

Thank you for a very readable explanation on snakes! I have one question though. The last equation assumes that the “force does not change between iterations” so solving an implicit scheme can be avoided. But shouldn’t the last parts of the formulae then be something like f_x(X_{t-1}, Y_{t-1}), f_y similar, otherwise the implicitness of the equations remains?

Chris,

Yes, indeed, you are right. My mistake. I will fix it right away.

Thank you!

Hi,

I was trying to implement snakes using open contours, after talk with Cris, I make it works using the same algorithm but changing the 2 first and 2 lasts rows of A:

That means we are using partial derivatives on the boundaries where second and fourth derivative are 0 at the extreme points. At the second and second to last points only the fourth is 0.

Hi Jordi,

setting P(0,0) to 1 and the rest of the row to 0 means that the first and last point of the open contour is only affected by the external forces and not the internal forces.

Depending on the external forces this might be a bad choice.

Stefan,

This is how I see it: If the end points of an open contour are affected by internal forces, this probably means that the contour will be shortened until it is a single point. Forces need to be in equilibrium after a number of iterations. If an end point is being pulled or pushed from one side, but not from the other, then there can never be an equilibrium.

Hello sir, thanks for your tutorial, i’m PhD student and my project on “improve gvf snake algorithm to restore broken characters” , i have a code about gvf snake but i want help to improve it and i will grateful for that . B.R.

hi sir, my inquiry about what is the value of “rr” when want to create an image.

Qusay,

`rr`

is a function in the DIPimage toolbox, it generates an image where each pixel’s value is its distance to the centre of the image.I’m sorry, but I won’t have time to help you develop algorithms or code.

Hello sir,

when I’m trying to execute the code, it says that a file named “dip_gradientmagnitude.mexw32” is missing. But I check the foler, the file is indeed here. Do you know how to solve this problem ? I’m a beginner in snake algorithms. Thanks a lot.

Hi Michael,

It is probable that the MEX-file cannot be loaded because one of the required libraries cannot be found. Did you follow the installation instructions in the user manual? Did you call

`dip_initialise`

? If so, did that call produce any error message?Sorry, my mistake. I didn’t install it correctly. Thanks for sharing your cleverness.

Hi,

thanks for this great article. Your blog is really interesting! I was a bit confused by the fact that you are apparently calculating the gradient out of the “gradient magnitude” and not of the image. Is this really what is happening?

Thanks

Hi Felix,

Yes indeed, that is what I’m doing in this case. The snake uses the gradient of the “external energy image”, which you can make however you want it. This image needs to have low values at the locations where you want the snake to move to. Thus I use the inverse of the gradient magnitude, which has low values at the edges of the object.

[…] It’s a bit older from 2009 but by far the clearest simple explanation with implementation… http://www.crisluengo.net/index.php/archives/217 […]

Great Tutorial! Thank you.

I am having a problem using the gradmag function. I am getting the following error when I use the “f = gradient(gradmag(img,30))” command:

Error using gradmag

Too many input arguments.

Do you know what is causing this?

Aneeq,

Did you install DIPimage? There’s a link on the right column of this page.

If yes, check that you don’t have another function with the same on the path. Type

`which gradmag`

.Hi Chris,

On your code you use gradmag(img, 30) but in all the papers and references they are using gradmag(img, 30)^2. Am I missing something?

Jordi,

Indeed, the original paper by Kass et al. proposes the square of the magnitude (I just checked). I guess it’s cheaper to compute, as you leave out a square root computation. Other than that, it increases the difference between sharp and smooth edges, I’m not sure if this is desirable.

Note that the original paper also proposes a lot of other external energy functions, and what you use as external energy depends on the application.

Hello ! I’d like to thank you for the really nice post about snakes .. i finally got around understanding

some more complicated aspects of the work…

even so, I was wondering if you could develop a bit more about the transition from:

AXt + fx(Xt,Yt)

to

(I – γA)-^-=1{Xt-1 + γ fx(Xt-1,Yt-1)}

… i don’t quite get the need of the identity matrix there..

sorry if it is a dumb question, but i’m kinda helpless with matrix operations

regards,

Wagner

Wagner,

If you move the

Xin_{t}(Xto one side, and multiply by γ, you get_{t}-X_{t-1})/γ = AX_{t}+ f_{x}(X. The left-hand side then simplifies to_{t}– γAX_{t}) = X_{t-1}+ γf_{x}(I – γA) Xbecause_{t}X. You cannot write_{t}= I X_{t}(1 – γA) X, that is not the same thing._{t}Great tutorial!

Request you to explain bit more about final converging criterion.

Thanks

Madhuri,

Good question. There is no guarantee that the snake will converge. This is why there’s a set number of iterations to take. But it is very possible that in one step the points don’t move much. You can detect that situation and terminate, assuming convergence.

Absolutely fantastic tutorial! Helped me greatly with implementing another snake algorithm. Great to see someone using a common sence approach to teaching this method by providing an implementation along with mathematics (which can often be overwhelming)

Hi Sir,

I found your tutorial very helpful. I just have a quick question. I’d like to know whether we can use this method for 3D images, like ultrasound cardiographs, or not?

Your time is greatly appreciated.

Regards,

Ara

Hello Sir,

Could you kindly tell me how will the matrix code change if I want to use this for open contours ?

Thank you

Ara,

There are 3D implementations of snakes, but they are not nearly as trivial as the 2D version discussed here. You might be better off using level sets instead.

Thanuja,

Please see the earlier comments for discussions on open contours.

Hi

what is rr in the first line?

`rr`

is part of the DIPimage toolbox (see link on the sidebar). It creates an image where each pixel’s value is the distance to the origin (centre of image by default). Thresholding it gives a disk.Very nice tutorial! Of course, you know this given that it looks like you’ve been getting these types of comments for 5 years now 🙂

In case anyone’s interested, I’ve implemented the above in Python here: https://github.com/chestervonwinchester/snakes

Hi Cris

Thanks for your nice explanation, I was wondering if I can use any other external force with snakes, for example: texture, distance measures, gray level, shapes, etc… and how can I incorporate this new external force into the equation?

Leiner,

I have not seen any definition of an external force that uses texture directly, but you could use some texture filter (e.g. Gabor) that converts the texture into gray levels.

Distance measures can be interpreted as gray levels also.

To use gray levels to define an external force, see the blog post. This is what snakes have traditionally used (either edges, i.e. change of gray values, or lines).

Hello Cris,

I know that this was blog was published a while ago, but I am trying to implement an open contour snakes algorithm and I am having difficulties.

I tried to implement the suggested changes to matrix A in the first 2 and last 2 rows, but I get very strange results that I do not understand. This is the code that I am trying to use (http://pastebin.com/icfM7H1k), as you will see the first snake iteration is all over the place. Any thoughts would be greatly appreciated.

Thanks!

Sean,

If your first iteration is all over the place, maybe you need to set gamma to a smaller value, so that each iteration makes a smaller step. That usually works in keeping things under control.

Thanks so much. This gives me the clearest understanding of how to apply Snake model on image contour tracking!!!

Hello Sir,

I am a pg student working on to detect the nails shape ,will this snake algorithm can be used to detect the nails of a human body.

Your time is greatly appreciated.Thanks alot!

Shilpa,

With snakes you always need to know first where the target object is. If you have a nail detection algorithm, you can then apply the snake to refine the outline, or to put pieces of outline together into a single shape.

Good luck!

need the code source for snakes please

Hiba,

I have updated the link to the code in the article. It got broken when I moved the blog to this new address. Thanks for letting me know!

Cris, thanks for a nice post! On my way to understanding. It really helps.

Have a question regarding matrix A:

Can you please explain why do you add these terms to your diagonals:

diag( b, -N+1);

diag( b, N-1);

diag([c,c],-N+2);

diag([c,c], N-2);

Sanzhar,

Those are the top-right and bottom-left corners of the matrix A. They represent the effect of points on the start of the snake on those at the end, and vice versa. If you leave these six values out, then the snake is not circular, it’s a cord with two loose ends.

Hi Cris,

I am trying to understand how you computed the matrix A whose values are represented in the combinations of α and β. For instance, how did you arrive at the A[0,0] value -2α-6β? I would highly appreciate your help. Thank you for your time!

Surya,

The matrix elements come from the two equations that are written just above it. The matrix elements at (0,0), (1,1), … (i,i) (i.e. the main diagonal) correspond to the component of the equation that multiplies

x(i). The diagonals one to the right and left of the main diagonal correspond to thex(i+1) andx(i-1) components, etc.Hi Cris,

You’ve got: AXt + fx = (Xt–Xt-1) / γ.

But in the original paper in equation (17)

They’ve got: AXt + fx= -γ(Xt–Xt-1).

Could you please clarify me why the original paper have minus for the gamma value?

Teera,

Honestly, I don’t remember why this difference, or who’s definition I was following when I wrote this post. But given that γ is a constant, you can set γ

_{mine}= -1/γ_{theirs}, and get the same equation back.I highly recommend that you write out the derivation for yourself, you’ll learn a lot from it, and you might discover typos in Kass’ paper or in my blog post. 🙂

Hi Cris,

You really make it very easy to understand.I am trying to extracting roads using ribbon snake. However I can’t reduplicate the code for ribbon snake. Can you write annother blog about ribbon snake? Thank you very much.

Hanting,

I’m glad this was useful to you. I have never heard of ribbon snakes before. I suggest you read whatever paper that proposes it a few more times, understanding new work is not always trivial. And if you get really stuck, contact its authors; they might be willing to share their code.

Hi Cris,

Very helpful tutorial!! I wrote my own program in C++ based on Kass’ paper, but my results are quite surprising.

The tricky task is the tuning of the weights alpha, beta, gamma. Imagine we do not consider the image energy at first time, only continuity and smoothness. Is the circle supposed to keep the same shape and size throughout the iterations? Mine is growing endlessly…

Is there a way to ensure a stable internal force?

Thank you very much,

Tristan,

Yes, it’s a balancing act. The smoothness constraint wants to minimize curvature, and a larger circle has a lower curvature than a smaller circle. You can add a balloon force (there’s a link to the paper towards the bottom of the blog post) to counteract that tendency, but if you initialize your snake close enough to the edges this should not be necessary.

Hey, great tutorial! I have one question: in the finite differences part, normally you would divide by h^2 for the second derivative and h^4 for the fourth, right? Where did you leave those, or what did you assume to discard them (or something else)? 🙂

Kate,

Yes, well spotted! Indeed, you would normally divide by

h^{2}andh^{4}for the second and fourth derivative, assuming the distances between samples are allh. In reality, the distances between samples (the points that form the snake) are all different, so we should divide by whatever distances are really there. But computing these distances every iteration would make the method significantly more expensive, as the values in the matrix we invert would keep changing every iteration, meaning we would have to invert that (large!) matrix every iteration.So to keep things simple, in the derivation of the snake it is assumed that these distances are all more or less the same. These factors

hare then folded into the constants α and β. If the points move closer together or further apart, these two constants effectively change, meaning that the speed of the snake changes too. It is good to resample the snake every certain number of iterations to avoid the distance between points to deviate too much. This resampling also avoids other issues that happen when points get too close together, and prevents the snake from folding in on itself. I think I didn’t discuss the resampling in this post, but it is in the linked implementation.Do you perhaps have a source for this specific discretization?

Kate,

I’m not are what you mean. This is a finite difference approximation to the derivative. Additionally we approximate by assuming *h* doesn’t change, even though it does.