Sparse Field Active Contours

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:
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

Median Filter and Morphological Dilation in Python

Python is a very nice programming language. Fast. Simple. Free. I recently spent some time learning it for a class on computer vision. I was using the PIL and numpy packages to make Python feel more like my old friend Matlab.

The two functions that I couldn’t find, and missed the most (especially when writing hack-y code for class projects) were median filtering and morphological dilation. So, in hopes of sparing other the pain of writing them… here they are! The function has both functions.

The medfilt() function uses the PIL filtering code. The dilate() function was written from scratch with NumPy.

Active Contour Matlab Code Demo

My new post: Sparse Field Active Contours
implements quicker, more accurate active contours.

Today, I added demo code for the Hybrid Segmentation project. This segmentation algorithm (in the publications section) can be used to find the boundary of objects in images. This approach uses localized statistics and sometimes gets better results than classic methods. For an example, see the video below: The contour begins as a rectangle, but deforms over time so that it finally forms the outline of the monkey.

This can be used to segment many different classes of image. To try it out, download the demo below and run >>localized_seg_demo

This code is based on a standard level set segmentation; it just optimizes a different energy. I’ve also made a demo which implements the well-known Chan-Vese segmentation algorithm. This technique is similar to the one above, but it looks at global statistics. This makes it more robust to initialization, but it also means that more constraints are placed on the image. Download it and see what you think! Again, unzip the file and run >>region_seg_demo (New! Described Here) (old and slow)

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!

A Short MEX Tutorial and Demo

Matlab is a great programming language/environment because of its ease of use, great visualization, and rapid prototyping abilities. Raw speed is not one of its strong suits. MEX (Matlab Executables) are the answer. These functions allow you to program in C or C++ (ultra fast languages), but be able to call and use them from Matlab programs. This post is a short intro to mex files which should get you up and running.

What This Post Teaches

In this post, I show how to create a mex file, how to set up inputs and outputs, how to get access to Matlab objects, and how to manipulate them. I also give a skeleton mex program that might be helpful. There is a lot more to learn, and I’d refer you to the mex manual regardless.

Continue reading “A Short MEX Tutorial and Demo”

GrowCut Segmentation In Matlab

I came across a cute segmentation idea called “Grow Cut” [pdf]. This paper by Vladimir Vezhnevets and Vadim Konouchine presents a very simple idea that has very nice results. I always feel that the simplest ideas are the best! Below I give a brief description of the algorithm and link to the Matlab/C/mex code.

GrowCut Region Growing Algorithm

This algorithm is presented as an alternative to graph-cuts. The operation is very simple, and can be thought of with a biological metaphor: Imagine each image pixel is a “cell” of a certain type. These cells can be foreground, background, undefined, or others. As the algorithm proceeds, these cells compete to dominate the image domain. The ability of the cells to spread is related to the image pixel intensity.

The authors give some pseudocode that very concisely describes the algorithm.

//for every cell p
for all p in image
  //copy previous state
  labels_new = labels;
  strength_new = strength;
  // all neighbors q of p attack
  for all q neighbors
      labels_new(p) = labels(q)
      strength(p) = strength_new(q)
    end if
  end for
end for

Segmentation Results

Once implemented, this is a nice way to get segmentations. It is quite fast, and the initialization is very intuitive. Consider this picture of a lotus flower:

growcut image

I made an initialization by clicking 20 points in the flower and 30 points outside. I then made a “label map” where unlabeled pixels are 0 (gray), foreground pixels are 1 (white) and background pixels are -1 (black).

growcut seeds

Based on this simple initialization, we obtain a very decent segmentation:

growcut output

As you can see, it isn’t perfect, but it is quite good. Its possible to interactively refine the seed points to improve the segmentation, but I didn’t do that here.

Matlab Code Downloads

I implemented this code in Matlab (using mex files due to the extensive use of for loops). You can download this below with compiled binaries for mac, linux, and windows. Unzip the file and run >>growcut_test for a demo.

UPDATE: I’ve fixed some bugs thanks to reader, Lin. The code works much better now!

Source & Compiled Binaries (96k) [zip]
“GrowCut” Paper [pdf]

Please let me know if you find this useful, and if you make improvements! Also, check out these related segmentation posts:

Related Segmentation Posts

Quick and simple derivatives in Matlab

Ever try to compute the directional derivatives on a matrix? Google turn up menacing formulas for Taylor expansions, grid spacing, and boundary conditions?

A quick and simple way of computing derivatives is to perform arithmetic on shifted versions of the matrix, and vectorized indexing help Matlab speed things up. For example, use the following for the (central) x-derivative:

dx = (M(:,[2:end end]) - M(:,[1 1:end-1]))/2

You can even define inline functions to perform the shift and derivative operations. Below are definitions for the shift operations and first and second order derivatives.

% shift operations
shiftD = @(M) M([1 1:end-1],:);
shiftL = @(M) M(:,[2:end end]);
shiftR = @(M) M(:,[1 1:end-1]);
shiftU = @(M) M([2:end end],:);

% derivatives
Dx = @(M) (shiftL(M) - shiftR(M))/2;
Dy = @(M) (shiftU(M) - shiftD(M))/2;
Dxx = @(M) (shiftL(M) - 2*M + shiftR(M));
Dyy = @(M) (shiftU(M) - 2*M + shiftD(M));
Dxy = @(M) (shiftU(M) - shiftD(M) + shiftL(M) - shiftR(M))/4;

And use them like this:
>> Ax = Dx(A)
Ax =
-7.0000 -6.5000 5.5000 5.0000
3.0000 2.5000 -1.5000 -1.0000
-1.0000 -1.5000 2.5000 3.0000
5.0000 5.5000 -6.5000 -7.0000

>> Ayy = Dyy(A)
Ayy =
-11 9 7 -5
15 -13 -11 9
-9 11 13 -15
5 -7 -9 11

Mean Shift Segmentation in Matlab


Recently I have decided to explore tracking from 3D point clouds extracted from stereo vision cameras. Step 1: Extract 3D point cloud from stereo vision cameras. So right now I’m implementing Segment-Based Stereo Matching Using Belief Propogation and Self-Adapting Dissimilarity Measure” by Klaus, Sormann, and Karner. This paper is defined by the source on stereo vision to be the best one around. This paper has two parts. Part 1: Segment the image. Part 2: Compute disparity (and depth) from the segments. Well, today I finished Part 1.

Stereo Cameras

First Try

The authors refer to a mean-shift segmentation algorithm presented in Mean Shift: A Robust Approach Toward Feature Space Analysis” [pdf] by Comaniciu and Meer to do the image segmentation. This paper (unlike some of my own previous work) leans towards oversegmentation of an image. Meaning that you prefer to get lots of little bits rather than the “right object” after the algorithm has run.

Well, after looking over the paper and getting a grasp for the mathematics, I took a crack at implementing it. Easily done… HOWEVER, my first attempt, written in Matlab, was painfully slow. (For a simple image it took 6 hours to run!) So, I got on the internet and came up with a better solution!

The Solution

Some great guys at Rutgers University implemented this paper in C++ and made the code available to the public under the name EDISON. (there’s also a nice GUI that goes along with this if you want to just play to see if these codes will work for you). Okay, so I had C++ codes that worked well (only 2 sec to do an image rather than 6 hours). The next step was to bring the code into Matlab.

Mean Shift Segmentation Results
These were the type of results I was trying for

I cracked my knuckles and got ready to write a MEX wrapper for this EDISON code. Then I said to myself, “Self, maybe you should check the ‘net first.” Turns out I had a good point. I found the website of Shai Bagon. Mr. Bagon had already made the MEX wrapper! Awesome.

I downloaded the codes and put them together. Mr. Bagon’s stuff worked right out of the box, although it would have saved me about an hour if I would have had this information (alternative readme.txt for Matlab Interface for EDISON). I also wrote my own wrapper-wrapper so that I could process grayscale images, and do simpler calls to accomplish what I wanted. If you’d like the code, download my wrapper-wrapper here (msseg.m).


Here is a sample of the output of this algorithm. The first image is a regular photo of some posed objects. The second image is the segmented version. Notice how the regions of the image are much, much more constant. This image has been broken into “tiles” of constant color.

Left Image
The original image (part of a standard pair of test images).
segmented image
The segmented image (ready to be processed in step 2)


Don’t re-invent the wheel. Taking a first crack at the implementation was good, and it helped me understand the algorithm. However, there was no need for me to spend a week tweaking it to be super-fast or two days getting the Matlab interface working. These things had already been done! It feels nice to knock out a task that you thought was going to take a week in a few hours : ) Stay tuned for the stereo part of this paper coming soon. Then maybe people will be writing about my page!

Ellipse Selection in Matlab

Much of the work on image segmentation that I do requires an initial guess of the answer. This initialization gives the algorithm something to improve. Hopefully, it improves all the way to the correct answer. You can use any manner of contour as an initialization. I like to select slightly different contours while I’m testing to make sure I’m not ‘over-tuning’ my system for some specific set of initial conditions.

Sometimes, it is preferable to use a very close initial outline like this:

Close initilization

(by the way, if you are interested in initializations like this see this code)

However, when you are trying to demonstrate the robustness of your technique (and/or you don’t want to waste time drawing a complicated contour) it is nice to use geometric shapes. Squares are popular, and can be captured using this simple Matlab code.

>>[rect] = getrect(1);

Many things in life are shaped more like ellipses, though. For this case, there are no built-in Matlab functions to get this type of initialization. Fear not! I wrote some simple code that allows you to capture graphically (or define in terms of parameters) an ellipse very simply. The params are major and minor radius (a, b), center location (x0, y0) and angle of rotation (rho).

>>[mask, a, b, x0, y0, rho] = get_ellipse(I, a, b, x0, y0, rho);

Here is an example of the type of initialization you can get (shown in white). Also in the image below you can see the final segmentation in green.

Ellipse Initialization

Download this .m file and try it out yourself.