Snakes: Active Contour Models

Implements snakes or active contour models for image segmentation.

Description

This demo implements the Active Contour Models as proposed by Kass et al.

To run it with GUI
   1. Type guide on the matlab prompt.
   2. Click on "Go to Existing GUI"
   3. Select the snk.fig file in the same directory as this file
   4. Click the green arrow at the top to launch the GUI

   Once the GUI has been launched, you can use snakes by
   1. Click on "New Image" and load an input image. Samples image are provided.
   2. Set the smoothing parameter "sigma" or leave it at its default value and click "Filter". This will smooth the image.
   3. As soon as you click "Filter", cross hairs would appear and using them and left click of you mouse you can pick initial contour location on the image. A red circle would appear everywhere you click and in most cases you should click all the way around the object you want to segment. The last point must be picked using a right-click in order to stop matlab for asking for more points.
   4. Set the various snake parameters (relative weights of energy terms in the snake objective function) or leave them with their default value and click "Iterate" button. The snake would appear and move as it converges to its low energy state.

Copyright (c) Ritwik Kumar, Harvard University 2010
              www.seas.harvard.edu/~rkkumar

This code implements “Snakes: Active Contour Models” by Kass, Witkin and Terzopolous incorporating Eline, Eedge and Eterm energy factors. See the included report and the paper to learn more.

If you find this useful, also look at Radon-Like Features based segmentation in the following paper:
Ritwik Kumar, Amelio V. Reina & Hanspeter Pfister, Radon-Like Features and their Application to Connectomics”, IEEE Computer Society Workshop on Mathematical Methods in Biomedical Image Analysis (MMBIA) 2010
 http://seas.harvard.edu/~rkkumar
Its code is also available on MATLAB Central

• Required Products: Image Processing Toolbox
• MATLAB release: MATLAB 6.5 (R13)

[Kirsty]

Quick point - you can speed up interate.m significantly if you take the expression for eterm out of the loop... i.e. replace with single calculation:

eterm = (cyy.*cx.*cx -2.*cxy.*cx.*cy + cxx.*cy.*cy)./((1+cx.*cx + cy.*cy).^1.5);

[Keith]

I received an error loading an image in the GUI using Mac OSX.

To fix this I removed the second argument of the strcat function on line 547 of snk.m - the '\' character.

[Fan Lin]

sometimes it work,but sometimes something bad happen:

??? Subscript indices must either be real positive integers or logicals.

Error in ==> interp2>linear at 336
    F = ( arg3(ndx).*(onemt) + arg3(ndx+1).*t ).*(1-s) + ...

Error in ==> interp2 at 215
    zi = linear(ExtrapVal,x,y,z,xi,yi);

Error in ==> interate at 124
    ssx = gamma*xs - kappa*interp2(fx,xs,ys);

Error in ==> snk>pushbutton2_Callback at 531
interate(handles.smth, handles.xs, handles.ys, alpha_val, beta_val, gamma_val, kappa_val,
weline_val, weedge_val, weterm_val, inter_val);

Error in ==> gui_mainfcn at 96
        feval(varargin{:});

Error in ==> snk at 85
    gui_mainfcn(gui_State, varargin{:});

Error in ==>
@(hObject,eventdata)snk('pushbutton2_Callback',hObject,eventdata,guidata(hObject))

 
??? Error while evaluating uicontrol Callback

[Junghun]

 

[omar]

good work man, but I have a problem in using snakes algorithm.

what I want is to get the contour of a mountain given only the top of the mountain.
I tried it using your implementation and gave it the surrounding indices of the top pixel but it was keeping smaller and smaller till vanishing. I know that is the idea of the algorithm, minimizing the area of the snake.

I want to give it the surrounding pixels of the top of the mountain and the region grows up till the stoppage point.

Is what I want feasible?

[Sean Lawson]

I have the same problem as Fan Lin.

Also, in the eterm part in interation.m, I notice that the original paper gives the formula (discarded the i,j here): eterm = (cyy.*cx.*cx -2.*cxy.*cx.*cy + cxx.*cy.*cy)./((cx.*cx + cy.*cy).^1.5); while you are using: eterm = (cyy.*cx.*cx -2.*cxy.*cx.*cy + cxx.*cy.*cy)./((1 + cx.*cx + cy.*cy).^1.5), could you please explain a little bit why you have a "1+" part in the denominator? Thank you

[Jeremy]

@ omar. The type of segmentation algorithm that you are describing is possible, but this program was not designed for it. I suggest you read about the gradient vector flow (GVF) on page 15 of "Biomedical Image Analysis: Segmentation" by Acton and Ray http://www.morganclaypool.com/doi/abs/10.2200/S00133ED1V01Y200807IVM009?prevSearch=allfield%253A%2528segmentation%2529&searchHistoryKey=

Additionally, I don't believe that the algorithm used in this program tries to minimize the area of the contour. Notice how the snake will expand to fill an area bounded by an edge.

[Adam]

 Good day
 I am a student of Biomedical Engineering.

 Should I request to you, I do not understand the function used Interam. The
 specific location:
 
 %masks for taking various derivatives
 m1 = [-1 1];
 m2 = [-1;1];
 m3 = [1 -2 1];
 m4 = [1;-2;1];
 m5 = [1 -1;-1 1];
 
 cx = conv2(smth,m1,'same');
 cy = conv2(smth,m2,'same');
 cxx = conv2(smth,m3,'same');
 cyy = conv2(smth,m4,'same');
 cxy = conv2(smth,m5,'same');
 
 and

 %populating the penta diagonal matrix
 A = zeros(m,m);
 b = [(2*alpha + 6 *beta) -(alpha + 4*beta) beta];
 brow = zeros(1,m);
 brow(1,1:3) = brow(1,1:3) + b;
 brow(1,m-1:m) = brow(1,m-1:m) + [beta -(alpha + 4*beta)]; % populating a
 template row
 for i=1:m
     A(i,: ) = brow;
     brow = circshift(brow',1)'; % Template row being rotated to egenrate
 different rows in pentadiagonal matrix
 end
 
 [L U] = lu(A + gamma .* eye(m,m));
 Ainv = inv(U) * inv(L); % Computing Ainv using LU factorization
 
 

 Many times I thank