Cris’ Image Analysis Blog

Last month I posted a tutorial on snakes. I explained that with snakes, the segmentation solution is obtained by minimizing an energy functional. The gradient descent minimization then yields an update function for the snake that can be seen as a set of internal and external forces acting on the snake. I was curious to see where we would get if we were to define the snake simply as a set of balls, connected with each other through springs, moving over the image with momentum, friction, and all the rest. For this, we will completely abandon the notions of gradient descent, function minimization, and so forth, and write a simple Newtonian physics simulation.

We start, as with last month’s post, defining a simple image a (with noise and an edge, requires DIPimage) and an initial position [px,py] of our snake (or rather the set of connected balls):

a = noise(50+50*gaussf(rr>85,2),'gaussian',30)
px = 126+50*cos(0:0.5:2*pi)';
py = 132+55*sin(0:0.5:2*pi)';
line([px;px(1)],[py;py(1)],'color',[0,0.8,0],'marker','o');

Next, we set the initial speed [vx,vy] for each point to 0:

vx = zeros(size(px));
vy = vx;

We now have the initial conditions for our simulation, and must decide what forces will act upon our system. The most important force is the one corresponding to the external force in the snake. This is the one that will pull the balls towards the edges in the input image. The gradient magnitude is large at edges, and the gradient of this gradient magnitude always points towards those edges. Thus, our first force [Fx,Fy] acting on the balls is given by (this part also requires DIPimage):

f = gradient(gradmag(a,10));
Fx = get_subpixel(f{1},[px,py],'linear');
Fy = get_subpixel(f{2},[px,py],'linear');

Traditional snakes have two components to their internal energy, one termed elasticity and one termed rigidity (also known as the membrane and thin plate energies). We will include two forces with a similar effect. Our elastic force is the one caused by the springs in between the balls, and forces the balls to be distributed equally along the contour. It pushes balls apart when they are close together, and pulls them together when they are far apart. To avoid instabilities, I chose to have each ball be attracted towards the neighbor that is furthest; the spring on the other side does not do anything. This is computed in the following way. We first calculate vectors to the two neighbors and their lengths. Then we select the longer of the two vectors as a force. The parameter alpha controls the strength of this force.

d1x = px([2:end,1])-px;
d1y = py([2:end,1])-py;
d1 = sqrt(d1x.^2+d1y.^2);
d2x = px([end,1:end-1])-px;
d2y = py([end,1:end-1])-py;
d2 = sqrt(d2x.^2+d2y.^2);
% Elasticity
f1 = d1-d2;
f1(f1<0) = 0;
f2 = d2-d1;
f2(f2<0) = 0;
Fx = Fx + alpha * ( f1.*d1x + f2.*d2x );
Fy = Fy + alpha * ( f1.*d1y + f2.*d2y );

(I know that this is not the most efficient programming, it is written for understandability.)

Our rigidity is also caused by the springs between the balls, but this time they both pull equally hard. What this does is move each point towards the middle point between its two neighbors, straightening the contour. This force is controlled by the parameter beta.

% Rigidity
Fx = Fx + beta * (d1x+d2x)/2;
Fy = Fy + beta * (d1y+d2y)/2;

Given an initially round snake, both these forces are directed inwards, causing the snake to shrink to a point. To counteract this effect, we can add an additional force that is perpendicular to the boundary. This would be equivalent to the balloon force in snakes. I call this the pressure force, and compute it by rotating the vector from the left to the right neighbor by 90 degrees:

% Pressure
dx = px([2:end,1])-px([end,1:end-1]);
dy = py([2:end,1])-py([end,1:end-1]);
Fx = Fx + kappa * dy;
Fy = Fy - kappa * dx;

The parameter kappa controls the strength of this force.

Leaving the friction out for a moment, all we need to compute now is the acceleration, , the new speed, , and the new positions, . Given these equations, we can compute a good value for , the step size, so that the maximum distance traveled by the balls is half a pixel. As long as we keep the time step short enough, this simple simulation will be OK. The simplest way to obtain this step size is by computing the maximum acceleration in the set, the maximum current speed in the set, and solve the quadratic equation . For this we simply set the mass to 1:

ax = Fx;
ay = Fy;
v = max(sqrt( vx.^2 + vy.^2 ));
a = max(sqrt( ax.^2 + ay.^2 ));
stepsz = (sqrt(v^2+2*a)-v)/(2*a);

Now that we know how large our time step is, we can compute the friction as well. Friction in general is proportional to the square of the velocity, and the deceleration is therefore proportional to . But if is small in comparison to the time step , we could run into problems by reversing the speed; deceleration should bring the speed towards 0, but never reverse it. Hence the following, slightly convoluted code that adds a check to make sure the speed never becomes negative:

% Friction
tmp = (1-(stepsz^2/mass)*sqrt( vx.^2 + vy.^2 ));
tmp(tmp<0) = 0;
vx = vx.*tmp;
vy = vy.*tmp;

This system of updating the speed before adding the acceleration is only to avoid this problem. Otherwise we could have simply computed a friction force and added it to [Fx,Fy]. Now we are ready to update the speed and position vectors:

vx = vx + stepsz*ax;
vy = vy + stepsz*ay;
px = px + stepsz*vx;
py = py + stepsz*vy;

Repeating this set of steps a few hundred times and plotting the results shows the balls moving around the image. It is important to finely tune the parameters. I have chosen them as follows so that, given a 0 external force, the contour defined by the balls does not change:

alpha = 0.005;
beta = 0.001;
kappa = 0.0002;
mass = 1;

I’ve put all of the code above into a function called snakemove. It contains some additional checks not discussed here, as well as code to plot the contour every couple of iterations. [Note: this function was written for DIPimage 2, and will not work unchanged with DIPimage 3.]

Let’s run 600 iterations of this code with the parameters and input as set above:

a
snakemove([px,py],f,alpha,beta,kappa,mass,600);

We can see how some of the balls are attracted to the contour first. Because of the rigidity and elasticity parameters, these balls pull the other balls with them, and eventually they all end up on the edge. If we add an initial velocity, we can see the effect of the mass of the balls:

a
vx = 20*cos(0:0.5:2*pi)';
vy = 20*sin(0:0.5:2*pi)';
snakemove([px,py,vx,vy],f,alpha,beta,kappa,mass,120);

Now the balls move outwards due to their initial velocity, only very slightly slowed by friction, but slowed more strongly by the external force once they overshot the edge. The external force slows them, and eventually sends them in the opposite direction, back towards the edge. The static situation is reached in much fewer steps.

We can do the same thing with many more balls. I’ve kept all the other parameters the same, and edited snakemove to not draw a circle for each ball. As you can see, the effect is approximately the same, but the “snake” is much more detailed.

a
px = 126+50*cos(0:0.1:2*pi)';
py = 132+55*sin(0:0.1:2*pi)';
snakemove([px,py],f,alpha,beta,kappa,mass,600);

a
vx = 20*cos(0:0.1:2*pi)';
xy = 20*sin(0:0.1:2*pi)';
snakemove([px,py,vx,vy],f,alpha,beta,kappa,mass,120);