Reference: Level Set Methods and Fast Marching Methods by J.A. Sethian, Cambridge University Press, 1999

First the distance phi of every point from a contour is determined. If the point is inside a contour, the distance is negative. If the point is outside all contours the distance is positive. The routine slicePointToContour(int slice, int x, int y) in VOI.java is used to obtain phi for 2D. In 3D pointToContour(int x, int y, int z) is used.

Equation 17.8 on page 222 of the reference is used: d(phi)/dt + gI(x,y)(1 - (epsilon)(kappa))|grad(phi)| - (beta)(grad(P(x,y))) dot grad(phi) = 0 where dot represents the dot product between (beta)(grad(P(x,y))) and grad(phi). Note that this equation is a development of the most basic form of the level set equation 1.6 on page 7: d(phi)/dt + F|grad phi| = 0

Equation 17.7 on page 220 gives: gI(x,y) = 1/(1 + |grad(Gaussian sigma * I(x,y))|) = 1/(1 - P(x,y)) An equation on page 222 gives: P(x,y) = -|grad(Gaussian sigma * I(x,y))| Note that in this code the |grad(Gaussian sigma * I(x,y))| is normalized to vary from 0.0 to 100.0. Equation 6.35 on page 69 gives the 2D curvature: The curvature kappa = the divergence of the normalized grad(phi) = ((phixx)(phiy)**2 - 2(phiy)(phix)(phixy) + (phiyy)(phix)**2)/ ((phix)**2 + (phiy)**2)**(3/2) Derivatives based on the central difference approximation give: phix = ((phi)i+1 - (phi)i-1)/(2(deltaX)) phiy = ((phi)j+1 - (phi)j-1)/(2(deltaY)) phixx = ((phi)i+1 - 2(phi)i + (phi)i-1)/(deltaX**2) phiyy = ((phi)j+1 - 2(phi)j + (phi)j-1)/(deltaY**2) phixy = ((phi)i+1,j+1 - (phi)i+1,j-1 - (phi)i-1,j+1 + (phi)i-1,j-1)/(4(deltaX)(deltaY))

The solution to this equation is based on solution 6.45 on page 74 to equation 6.44 on page 73. Equation 6.44 is similar to equation 17.8. Note that the equation is currently set up to perform an expansion. The discussion on page 220 says use 1 - (epsilon)(kappa) for expansion and -1 - (epsilon)(kappa) for contraction. (phi)i,j,n+1 = (phi)i,j + deltaT[c1 + c2 + c3] where: c1 = -[max(gI(i,j),0)gradPlus + min(gI(i,j),0)gradMinus] for expansion c1 = -[max(-gI(i,j),0)gradPlus + min(-gI(i,j),0)gradMinus] for contraction Since gI(i,j) is always between 0 and 1: c1 expansion = -gI(i,j)gradPlus for movement = EXPAND or EXPAND_CONTRACT c1 contraction = gI(i,j)gradMinus for movement = CONTRACT gradPlus is given by equation 6.18 on page 65: gradPlus = [max(Dminusx,0)**2 + min(Dplusx,0)**2 + max(Dminusy,0)**2 + min(Dplusy,0)**2]**0.5 gradMinus is given by equation 6.19 on page 65: gradMinus = [max(Dplusx,0)**2 + min(Dminusx,0)**2 + max(Dplusy,0)**2 + min(Dminusy,0)**2]**0.5 Dminusx = ((phi)i,j - (phi)i-1,j)/deltaX Dplusx = ((phi)i+1,j - (phi)1,j)/deltaX Dminusy = ((phi)i,j - (phi)i,j-1)/deltaY Dplusy = ((phi)i,j+1 - (phi)i,j)/deltaY c2 = -[max(u(i,j),0)Dminusx + min(u(i,j),0)Dplusx + max(v(i,j),0)Dminusy + min(v(i,j),0)Dplusy] where u(i,j) is the first vector component of -beta(grad(P(x,y))) and v(i,j) is the second vector component of -beta(grad(P(x,y))). c3 = epsilon(kappa(i,j))(gI(i,j))(((Dzerox)**2 + (Dzeroy)**2)**(1/2)) Dzerox = phix = ((phi)i+1,j - (phi)i-1,j)/(2(deltaX)) Dzeroy = phiy = ((phi)i,j+1 - (phi)i,j-1)/(2(deltaY))

The user supplies 7 constants needed for the solution of the equation: 1.) epsilon 2.) sigma for the gaussian derivative convolution 3.) The deltaT time step value 4.) The number of iterations 5.) movement = EXPAND, EXPAND_CONTRACT, or CONTRACT EXPAND_CONTRACT is the EXPAND mechanism except that it does not prevent a value of phi from transiting from a negative value to value >= 0.0. movement == EXPAND prevents a boundary from ever retreating - phi is prevented from ever changing from negative to >= 0. movement == CONTRACT prevents a boundary from ever expanding - phi is prevented from ever changing from >= 0.0 to negative. 6.) edgeAttract the ratio of the maximum value of the absolute value of the edge attractive force to the maximum value of the absolute value of the sum of the propagation and curvature forces. 7.) testIters if testIters > 0, check if the boundary is unchanged every testIters iterations. On page 221 the author used time step = 0.001 and 391 iterations on a 64 by 64 grid. In the paper A Fast Level Set Based Algorithm for Topology-Independent Shape Modeling by Malladi, Sethian, and Vemuri a similar scheme was used. epsilons of 0.05, 0.25, and 0.75 were used. As epsilon was increased, the level set attained a smoother final configuration. The (epsilon)(kappa) term smoothes out the high curvature regions and has the same regularizing effect as the internal deformation energy term in thin-plate-membrane splines. Time steps ranging from 0.00025 to 0.001 were used. One run on a 64 by 64 grid using 0.001 time steps took 391 iterations to complete and one run on a 128 by 128 grid using 0.0005 time steps took 575 iterations to complete. a sigma = 3.25 was used. In Shape Modeling with Front Propagation: A Level Set Approach by Malladi, Sethian, and Vemuri a 128 by 128 image with 0.00025 time steps took 1,180 iterations.

The term grad(P) dot grad(phi) denotes the projection of an (attractive) force vector on the surface normal. This force results from gradient of a potential field P(x,y) = -|grad(Gaussian sigma * I(x,y)| that attracts the contours to the image edges. Beta controls the strength of this attraction.

When the iterations are complete, the area inside the curve has negative phi values and the area outside the curve has positive phi values. Bitset a mask from the pixels with negative phi values and perform a mask to VOI conversion to obtain the output curve.

Extensions to 3D are very straightforward except possibly for the 3D mean curvature. Equation 6.36 on page 70 is used: mean curvature = ((phiyy + phizz)phix**2 + (phixx + phizz)phiy**2 + (phixx + phiyy)phiz**2 - 2(phix)(phiy)(phixy) - 2(phix)(phiz)(phixz) -2(phiy)(phiz)(phiyz))/ (phix**2 + phiy**2 + phiz**2)**1.5