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 s(0) = s(1).

The energy functional that we want to minimize is defined as: E = ∫ Eint(s(p)) + Eext(s(p)) dp Eint is an internal energy, meant to force the snake to be small and smooth. It avoids obviously wrong solutions. Eext is an external energy, the one responsible for the snake finding the edges of an object in the image. Eext 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 Eint(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) – ∇Eext(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!) ∇Eext, the gradient of the external energy, is a force. We call this the external force, Fext, 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 s 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(p) / ∂t = αs″(p,t) – βs″″(p,t) + Fext(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 s into its two components x and y. The two components of the external force Fext are fx and fy. We use the finite difference approximation to the spatial derivatives: x″[i] = x[i-1]-2x[i]+x[i+1] and x″″[i] = x[i-2]-4x[i-1]+6x[i]-4x[i+1]+x[i+2]. This results in the following set of 2N equations: xt[i] / ∂t = α(xt[i-1]-2xt[i]+xt[i+1]) + β(-xt[i-2]+4xt[i-1]-6xt[i]+4xt[i+1]-xt[i+2]) + fx(xt[i],yt[i])yt[i] / ∂t = α(yt[i-1]-2yt[i]+yt[i+1]) + β(-yt[i-2]+4yt[i-1]-6yt[i]+4yt[i+1]-yt[i+2]) + fy(xt[i],yt[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: Xt / ∂t = AXt + fx(Xt,Yt)Yt / ∂t = AYt + fy(Xt,Yt) 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 γ: ∂Xt / ∂t = (XtXt-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: Xt = (IγA)-1{Xt-1 + γ fx(Xt-1,Yt-1)}Yt = (IγA)-1{Yt-1 + γ fy(Xt-1,Yt-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

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. That package also shows how to integrate such a new parameter type into the DIPimage interface, which I will discuss in another future post.

74 Responses to “A simple implementation of snakes”

  1. On March 14th, 2010, at 12:32, Rahul said:

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

  2. On September 6th, 2012, at 15:37, Lucas said:

    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!

  3. On September 6th, 2012, at 15:58, Cris Luengo said:


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

  4. On September 20th, 2012, at 12:46, karthik said:

    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.



  5. On September 21st, 2012, at 8:42, Cris Luengo said:


    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 before you apply the snake!

  6. On October 8th, 2012, at 12:20, Nathanael said:


    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?

  7. On October 8th, 2012, at 15:58, Cris Luengo said:

    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;

  8. On October 17th, 2012, at 18:29, Students Of University of Catania said:

    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

  9. On October 18th, 2012, at 7:29, Cris Luengo said:


    The matrix I is the identity matrix (all elements are zero except the main diagonal, that is all ones). fx and fy are 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.

    In your example with a snake with 3 points (1,5) (3,9) (30,9), X0 = (1,3,30)T and Y0 = (5,9,9)T encode the initial coordinates for the snake.

    Hope this clarifies things!

  10. On October 23rd, 2012, at 15:32, Dott.Marika said:

    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


  11. On October 23rd, 2012, at 19:59, Cris Luengo said:

    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.

  12. On February 7th, 2013, at 11:59, Vikram said:

    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 ?

  13. On February 7th, 2013, at 12:05, Vikram said:

    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.

  14. On February 7th, 2013, at 13:28, Cris Luengo said:

    Dear Vikram,

    The matrix A would 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.

  15. On April 28th, 2013, at 11:18, James said:

    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?


  16. On April 28th, 2013, at 12:03, Cris Luengo said:


    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.

  17. On May 2nd, 2013, at 6:30, Sam said:

    Hi Cris, great tutorial. I was just trying to implement this without the toolbox, I’ve replaced

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


    [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’);


    F = griddedInterpolant(FX, ‘linear’);

    for i = 1:length(x)



    G = griddedInterpolant(FY, ‘linear’);

    for i = 1:length(y)



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


  18. On May 2nd, 2013, at 7:47, Cris Luengo said:

    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).

  19. On May 25th, 2013, at 14:26, M.Jalali said:

    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…

  20. On July 10th, 2013, at 14:12, Edgar Delgado said:


    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:

       Snakes.exe 100 0.01 0.4 0.07 4 6000

    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]

  21. On October 15th, 2013, at 23:11, Michelle said:

    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.



  22. On October 16th, 2013, at 9:23, Cris Luengo said:

    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.

  23. On October 16th, 2013, at 15:04, xiaofei said:

    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!

  24. On October 16th, 2013, at 23:32, Cris Luengo said:

    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).

  25. On October 18th, 2013, at 11:23, swapnil kulkarni said:

    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’);

  26. On October 21st, 2013, at 14:07, Cris Luengo said:

    Hi Swapnil,

    Did you install DIPimage?

  27. On October 25th, 2013, at 7:38, swapnil kulkarni said:

    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.

  28. On November 2nd, 2013, at 18:19, Chris Zeinstra said:

    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?

  29. On November 3rd, 2013, at 14:00, Cris Luengo said:


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

    Thank you!

  30. On November 13th, 2013, at 10:44, Jordi Villar said:


    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:

    1	0	0	0	0	0	⋯	0	0	0
    α	-2α	α	0	0	0	⋯	0	0	0
    -β	α+4β	-2α-6β	α+4β	-β	0	⋯	0	0	0
    0	-β	α+4β	-2α-6β	α+4β	-β	⋯	0	0	0
    0	0	-β	α+4β	-2α-6β	α+4β	⋯	0	0	0
    ⋮	⋮	⋮	⋮	⋮	⋮	⋱	⋮	⋮	⋮
    0	0	0	0	0	0	⋯	α	-2α	α
    0	0	0	0	0	0	⋯	0	0	1

    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.

  31. On January 10th, 2014, at 18:07, Stefan said:

    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.

  32. On January 10th, 2014, at 22:40, Cris Luengo said:


    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.

  33. On February 2nd, 2014, at 12:59, Qusay Mosa said:

    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.

  34. On February 2nd, 2014, at 13:43, Qusay Mosa said:

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

  35. On February 2nd, 2014, at 16:52, Cris Luengo said:


    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.

  36. On February 6th, 2014, at 10:29, Michael B said:

    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.

  37. On February 6th, 2014, at 13:56, Cris Luengo said:

    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?

  38. On February 7th, 2014, at 10:24, Michael B said:

    Sorry, my mistake. I didn’t install it correctly. Thanks for sharing your cleverness.

  39. On February 16th, 2014, at 23:29, felix said:


    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?


  40. On February 23rd, 2014, at 17:34, Cris Luengo said:

    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.

  41. On March 8th, 2014, at 18:06, Safar said:

    Hello, i’m a beginner with matlab and this family of algortihms, so i found some problems when executing the program below, i need some help please (i’m using matalb 7)sorry for my poor english.

  42. On March 10th, 2014, at 16:30, Cris Luengo said:


    I have no way of possibly helping you if you can’t explain what problem you encountered.

  43. On March 12th, 2014, at 23:29, Awesome Simple Tutorial on Snakes / Active Contours by Cris Luengo! | S1nch@i said:

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

  44. On April 8th, 2014, at 9:49, Aneeq said:

    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?

  45. On April 8th, 2014, at 15:01, Cris Luengo said:


    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.

  46. On April 15th, 2014, at 15:26, Jordi Villar said:

    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?

  47. On April 15th, 2014, at 15:49, Cris Luengo said:


    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.

  48. On May 29th, 2014, at 1:01, Wagner Rampon said:

    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)


    (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



  49. On June 2nd, 2014, at 10:01, Cris Luengo said:


    If you move the Xt in (Xt-Xt-1)/γ = AXt + fx to one side, and multiply by γ, you get (Xt – γAXt) = Xt + γfx. The left-hand side then simplifies to (I – γA) Xt because Xt = I Xt. You cannot write (1 – γA) Xt, that is not the same thing.

  50. On June 20th, 2014, at 13:27, Madhuri said:

    Great tutorial!

    Request you to explain bit more about final converging criterion.


  51. On June 23rd, 2014, at 14:14, Cris Luengo said:


    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.

  52. On August 12th, 2014, at 12:19, Dan Shayler said:

    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)

  53. On April 16th, 2015, at 0:40, Ara said:

    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.



  54. On April 18th, 2015, at 7:51, thanuja said:

    Hello Sir,

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

    Thank you

  55. On April 21st, 2015, at 10:29, Cris Luengo said:


    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.

  56. On April 21st, 2015, at 10:29, Cris Luengo said:


    Please see the earlier comments for discussions on open contours.

  57. On May 24th, 2015, at 14:01, vahide said:


    what is rr in the first line?

  58. On May 25th, 2015, at 12:27, Cris Luengo said:

    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.

  59. On May 29th, 2015, at 1:34, chestervonwinchester said:

    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

  60. On July 24th, 2015, at 20:33, Leiner said:

    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?

  61. On July 29th, 2015, at 14:03, Cris Luengo said:


    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).

  62. On April 6th, 2016, at 19:02, Sean said:

    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.


  63. On April 7th, 2016, at 5:10, Cris said:


    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.

  64. On April 19th, 2016, at 22:02, Johnson said:

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

  65. On August 17th, 2016, at 5:52, Shilpa Rajpal said:

    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!

  66. On August 17th, 2016, at 21:22, Cris said:


    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!

  67. On September 7th, 2016, at 16:38, Hiba said:

    need the code source for snakes please

  68. On September 7th, 2016, at 23:46, Cris said:


    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!

  69. On November 16th, 2016, at 6:04, Sanzhar said:

    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);

  70. On November 16th, 2016, at 8:56, Cris said:


    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.

  71. On April 15th, 2017, at 17:35, Surya said:

    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!

  72. On April 15th, 2017, at 22:31, Cris said:


    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 the x(i+1) and x(i-1) components, etc.

  73. On May 24th, 2017, at 9:29, Teera said:

    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?

  74. On May 25th, 2017, at 23:11, Cris said:


    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. 🙂

Leave a Reply

You can use these HTML tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

Note: I moderate all comments. Comments without a clear relation to the text above will not be published.