Active contour methods for image segmentation allow a contour to deform iteratively to partition an image into regions. Active contours are often implemented with level sets. The primary drawback, however, is that they are slow to compute. This post presents a technical report describing, in detail, the sparse field method (SFM) proposed by Ross Whitaker [pdf], which allows one to implement level set active contours very efficiently. The algorithm is described in detail, specific notes are given about implementation, and source code is provided.

### Fast Level Sets Demo

The links below point to the technical report and a demo written in C++/MEX that can be run directly in MATLAB. The demo implements the Chan-Vese segmentation energy, but many energies can be minimized using the provided framework.

Sparse Field Method – Technical Report [pdf]

Sparse Field Method – Matlab Demo [zip]

To run the MATLAB demo, simply unzip the file and run:

`>>sfm_chanvese_demo`

at the command line. On the first run, this will compile the MEX code on your machine and then run the demo. If the MEX compile fails, please check your MEX setup. The demo is for a 2D image, but the codes work for 3D images as well.

My hope is that other researchers wishing to quickly implement Whitaker’s method can use this information to easily understand the intricacies of the algorithm which, in my opinion, were not presented clearly in Whitaker’s original paper. Personally, these codes have SUBSTANTIALLY sped up my segmentations, and are allowing me to make much faster progress towards completing my PhD!

Thanks to Ernst Schwartz and Andy for helping to find small bugs in the codes and documentation. (they’re fixed now!)

This code can be used according to the MIT license. As long as this work is appropriately cited and attributed, and not being used for proprietary or commercial purposes, I’m fully supportive of you using it. Please drop me a line if it helps you!

For more information regarding active contour, segmentation, and computer vision, check here: Computer Vision Posts

Hi Shawn,

Could you please tell me the reason why you did not normalize the curvature force like you did for the external region force?

I am referring to this line btw.

F = F/(max(abs(F))+eps) – 1.0*kappa(phi, C);

Eran

Dear Shawn,

It’s easy and useful of you provided code. It is great work. I will study it careful and do some applications. I will cite the paper of yours. It seems that, you are a very kind man.

Thanks and Regards.

XinHK

I can’t get it to run. I’m using R2009b, and I have the path set correctly. I get:

>> sfm_chanvese_demo

??? Undefined function or variable ‘sfm_chanvese_demo’.

So I tried:

>> sfm_local_demo

??? Undefined function or method ‘fat_contour’ for input arguments of type

‘double’.

Error in ==> visualize_phi at 4

[h1 h2] = fat_contour(phi);

Error in ==> sfm_local_chanvese at 79

[phi c] = sfm_local_chanvese_mex(img,mask,iterations,display,lam,rad);

Error in ==> sfm_local_demo at 22

[seg] = sfm_local_chanvese(img,mask,iterations,lambda,rad);

Any idea what I’m doing wrong?

Okay, got it. I had to

1. Download fat_contour.m

2. Comment out the delete(h1) and delete(h2)

3. make the change described above with the ~

4. clear all

And then it works!

you have a bug in your code provided at mathworks (http://www.mathworks.co.jp/matlabcentral/fileexchange/23847-sparse-field-methods-for-active-contours) in line 132, please check the signs.

the code here is already patched, maybe you want to update the files at mathworks

regards

@Travis

The error is still there. Just downloaded the code from the homepage and “Undefined function or method ‘fat_contour'”

How to run this this code for a 3D image( a stack of slices)

Hey Shawn:

So apparently the technical write up is broken or something can’t really open it. Can you upload another copy? Thanks.

@Matias Perez

Working on getting a non-corrupted copy… I wonder how this happened!?

–Fixed!

Dear Shawn,

Thank you for your thorough work on the sparse field level set algorithm. It is helpful. I would like to ask if you have any idea about how to deal with anisotropic 2D or 3D image spacing. Apart from the correct calculation of the gradients and the forces, the sparse field algorithm supposes that the neighbors of, say, the zero list member, are all either layer-1 or layer-minus-1 members, whereas when the z-factor (i.e. intra slices spacing over inter slice spacing) is larger than 1.5 this is impossible. Do I miss anything?

Thanks a lot

Hadar

@Hadar

As you said, correct computation of the gradients (i.e., based on the true inter-pixel spacing) is important.

However, as long as the spacing is on a regular grid (even if x, y, and z spacings are different), the SFM should continue to work because of the floating point “level set” structure that’s maintained. That structure allows you to interpolate the “zero level set” at sub-pixel resolution regardless of pixel spacing.

You can also consider having more adjacent layers if you’re spacing gets really huge! Whittaker describes how to do this, and it’s a trival extension of my published algorithm.

Finally, a colleague of mine wrote codes to perform SFM on very non-uniform data (e.g., point clouds). I could try to track those codes down if people would find them helpful.

hello,

First thanks for your program. really thank you.

i have 2 questions.

first about CFL condition i don’t understand why this condition and also why you multiply it by dphidt in order to evolve contour.

Second: why the expression of curvature is different from whitaker(page 10 “A Level-Set Approach to 3D Reconstruction From Range Data”).

thanks an advance

hi every body,

i test the region_seg_demo. it works very well. but i have some questions about the equations used in this code.

F = (I(idx)-u).^2-(I(idx)-v).^2; % force from image information

curvature = get_curvature(phi,idx); % force from curvature penalty

dphidt = F./max((F)) + alpha*curvature; % gradient descent to minimize energy

%– maintain the CFL condition

dt = 3/(max(dphidt)+eps);

%– evolve the curve

phi(idx) = phi(idx) + dt.*dphidt;

%– Keep SDF smooth

phi = sussman(phi, .5);

could some one clarify my these equations.

thanks an advance.

i have error in sfm_local_chanvese_mex.cpp file.

Error is as follow:

error C2660: ‘en_chanvese_compute’ : function does not take 6 arguments

Please suggest solution.

Is the algorithm in your technical report up to date?

Since I don’t think that it preserves the distance invariant, which leads to inconsistencies in the distance function, duplicates in the lists, and ends in a mess.

@Zohar The algorithm is up-to-date as far as I know… Not sure what you observed, but I had good results. When using non-uniform pixel spacing, you need to pass those as arguments to the functions.

Thanks for replying.

So it probably my mistake, since I’m not sure what is the invariant the algorithm keeps in each level. For example I assumed the algorithm tries to preserve for each level its data range. For example Lp2 should always have values between (1.5,2.5]. But according to line 7 at procedure 2, p would still have a Lp1 value and it would be transferred to Lp2. I find it hard to follow all the cases, and see if this is crucial or not, but currently I admit that the algorithm behaves well, and my previous impression was due to a bug in my implementation. Maybe the worst that could happen are orphan cells (as implied by [Whitaker98], or wrong approximation of the distance function at higher levels. I would appreciate your opinion on the matter.

Thanks,

Zohar

Also I was wondering if you noticed any precision degradation due to reinitialization of the level set in each iteration. To quote [Sethian99] ‘Motion by Intrinsic Laplacian of Curvature’: “The use of reinitialization provides the extra stability required to complete the numerical method, however it comes at a price. Reinitialization can also be seen as a form of artificial diffusion. The amount of artificial diffusion added to the problem and the amount of the error introduced in slight movement of the interface varies between reinitialization methods. In any case we use reinitialization as little as possible without compromising the stability of the method”.

Which means according to Sethian that your method should be very stable (you also preserve the CFL condition by normalizing the force, and limiting it to 0.5), but at a price of precision.

Did you try a fourth order flow such as the Laplacian curvature of Willmore flows?

Concerning my remark from September 19th, 2011. There is a basic thing that I don’t understand here, and I’ll appreciate an explanation.

Why don’t we simply in each iteration calculate F for the zero level set, then reconstruct from scratch all other levels? More specifically:

1. Calculate F for Lz.

2. Go over all the other lists Ln[], Lp[], and delete them, while setting the cells they were pointing to to the highest label (+/-3 in your case).

3. Iterate Lz:

3.1 Update the cell value according to F.

3.2 If the value is still in Lz interval, add the neighbors with label > 0 to Lp1 with a value of the Lz cell + 1, and neighbors with label i

to Lp(i+1). Update the value of all neighbors with label i+1 to Lpi’s value + 1.

It would replace your 2,3 procedures, won’t need the intermediate Sx lists, would be faster, and it would be cleaner, in the sense that it would keep the distance invariant.

@Zohar

Thanks for the comments!

The approach you lay out may be faster in some situations (e.g., where there are few voxels in Lz). In the case where the number of voxels gets large (and only a sub-set of those voxels cross the zero-level-set on each iteration) the method presented here will be MUCH faster.

In terms of invariant distances, it’s important to ensure floating-point accuracy on the phi function to ensure sub-voxel accuracy of the surface.

I urge you to please consider my point again. I think that my method is faster: In your method you go over all the items in all the lists, update the values (we update the same way with floating point accuracy), and then move the necessary elements to Sx lists, including new discovered cells. Then you go over the Sx lists, and move the items to their corresponding Lx. I also iterate the *same* amount of items in the lists, when I rediscover them from scratch, and instead of the Sx lists, I deleted the items beforehand. You can argue that allocating new cells for the list would take time, but you can hold an unoccupied cell pointer list, and just move pointers from list to list without any allocation. Again the time you spend to update an element according to its neighbors, would be the same amount of time I spend to rediscover the item from scratch.

Let’s check my claim on an extreme case, where your method should have a better advantage over mine. The new F would be zero for all Lz. I would iterate once over Ln and Lx clearing all cells value. Then I would iterate Lz again, and rediscover all the neighbors. It would mean that if we have N cells in the the active, I would have 2N accesses. You on the hand would have only N access.

Considering the other extreme case that the whole Lz moved, you would need to relocate all your items through Sx. So you’ll spend 3N, while I still spend 2N. Then we should consider the amortized average case, although if you normalize the force as you do, the level set should change quite often.

Well, I admit my examples were a bit sketchy, but I hope that you get the idea, that we access the memory, and move the list points quite in the same amount, and I don’t think it would matter much in terms of performance.

Why my solution is cleaner: The code is shorter and simpler, and it’s really easy to track the changes, and prove that the distance invariant is kept (cell in level i has it value in the allowed distance range), as opposed to the one example I gave you that your invariant can fail (I noticed it also in my experiments).

Shawn Lankton

Hi, thank you for sharing this code, it’s works for what I want to do. but it hard for me to understand all of the mathematical equations on your paper “Localized Region Based Active Contours”(code:Chan&Vese code); I wondering if you have any link where can I find the explication or more details?. if you have something, this is my @mail: kasmireda@gmail.com.

thanks and regards

Reda

Hai,

Ur code is really helpful. May I know how the 3D SFM works?

Hi Shawn, thanks for sharing code.

Please, can you explain me why you use

seg = phi<=0

to create the final label map result instead of

seg = phi<=0.5

since the range for the zero level set in phi is [-0.5,0.5]?

thanks in advance

Hi Shawn, thanks for sharing code.

Please, can you explain me why you use

seg = phi<=0

to create the final binary map instead of

seg = phi<=0.5

if the range for the zero level set in phi is [-0.5,0.5]?

thanks in advance

Hi Shawn, thanks for sharing code.

Please, can you explain me why I can’t get it to run. I’m using R2007b. I get while running demo code:

Compiling sfm_chanvese_mex…

Error: Could not detect a compiler on local system

which can compile the specified input file(s)

C:\PROGRA~1\MATLAB\R2007B\BIN\MEX.PL: Error: No compiler options file could be found to compile source code. Please run “mex -setup” to rectify.

??? Error using ==> mex at 208

Unable to complete successfully.

Error in ==> compile_sfm_local_chanvese at 6

mex sfm_local_chanvese_mex.cpp llist.cpp energy3c.cpp lsops3c.cpp

Error in ==> sfm_local_chanvese at 30

compile_sfm_local_chanvese

Error in ==> sfm_local_demo at 22

[seg] = sfm_local_chanvese(img,mask,iterations,lambda,rad);

thanks an advance.

Annet

Download code from following link and check :

http://www.mathworks.com/matlabcentral/fileexchange/23847

@annet

Is there any one can tell me how to decide the which part of the mask equals to 0?

“mask(86:218,109:238) = 0”. Or how is the number 86:218 gotten?

Thank you very much.

bonjour peut ton donne un methode pour ecrier le code de contour actif en VHDL

bonjour sawan peut ton donne un methode pour ecrier le code de contour actif en VHDL

bonjour Dear Shawn

peut ton donne un methode pour ecrier le code de contour actif en VHDL

hi sir,

first and foremost tnx 4 d code. when im trying to run it im getting following errors. plz anybody help me out. im using it for some medical applications. plz help me fast as soon as possible.

??? Error: File: sfm_local_chanvese_mex.m Line: 6 Column: 21

Unexpected MATLAB operator.

Error in ==> sfm_local_chanvese at 79

[phi c] = sfm_local_chanvese_mex(img,mask,iterations,display,lam,rad);

Error in ==> sfm_local_demo at 23

[seg] = sfm_local_chanvese(img,mask,iterations,lambda,rad);

Hi,first thanks for sharing your code.

I have a problem, you use double link sparse field method to evlove the curve. can I change it into a single link? If I change it into a single link, is there any problem?

Nothing comes to mind as to why it wouldn’t work… I’d say try it and see how it goes! If it works and it’s faster let me know and I’ll link to your code.

@Shawn Lankton

yes, I have changed the code to single link(refer your code and J Malcolm’s code), it can work. I test the code on my data set, The running time is not declining, but the accuracy is increased. It is strange. I don’t know why it is happend.

another question, do you have any ideas to reduce running time? @Shawn Lankton

@Shawn Lankton hey can you please tell me which papers or which articles you used for this script. I just want to know on which parameters i can play this algorithm.

@yogesh Check the bibliography of the technical report.

Sorry, I make an error in the program. Single link is more faster than double link and the accuracy is about the same….Thanks! @Dily

@Dily That’s awesome – great work! If you want to share a link I’ll post it here so folks can find your code.

Thanks, I’d like to share the code. Because the code has changed a lot(for personal reasons), I will need some time before it is finished.@Shawn Lankton

Hi All,

I am looking for code for level set based shape prior segmentation. Does anybody know whether such a code is available?

Thanks.

Hello sir,

I need a code for detecting lip, lip motion and lip segmentation. Plz help me to find the code for this.

Regards

J.T.

How i can use the same code for 3d array

Lets say i have an array img = 512*512*100;

what will be mask ?

@Dily: Do you finish single link? Please share link