Package gov.nih.mipav.model.algorithms

Class AlgorithmLevelSet

java.lang.Object

java.lang.Thread

gov.nih.mipav.model.algorithms.AlgorithmBase

gov.nih.mipav.model.algorithms.AlgorithmLevelSet

All Implemented Interfaces:

AlgorithmInterface, ActionListener, WindowListener, Runnable, EventListener

public class AlgorithmLevelSet
extends AlgorithmBase
implements AlgorithmInterface

This algorithm iteratively expands or contracts one or more contours to a boundary.

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

Nested Class Summary

Nested classes/interfaces inherited from class java.lang.Thread

Thread.Builder, Thread.State, Thread.UncaughtExceptionHandler

Field Summary

Fields

Modifier and Type	Field	Description
private static final int	CONTRACT	DOCUMENT ME!
private float	deltaT	Time step The defalut is 0.05 for 2D and 0.025 for 3D.
private float	edgeAttract	edgeAttract is 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.
private float	epsilon	epsilon is used to smooth high curvature regions. epsilon is between 0 and 1 with larger epsilon doing more smoothing
private static final int	EXPAND	DOCUMENT ME!
private static final int	EXPAND_CONTRACT	DOCUMENT ME!
private float[]	GxData	Storage location of the first derivative of the Gaussian in the X direction.
private float[]	GyData	Storage location of the first derivative of the Gaussian in the Y direction.
private float[]	GzData	Storage location of the first derivative of the Gaussian in the Z direction.
private boolean	image25D	If image25D is true in 3D images, process each slice separately
private int	iterations	Number of iterations of the diffusion.
private int[]	kExtents	Dimensionality of the kernel.
private int	movement	EXPAND, EXPAND_CONTRACT, or CONTRACT.
private float[]	outputBuffer	Stores output of AlgorithmConvolver
private float[]	sigmas	Standard deviations of the gaussian used to calculate the kernels.
private int	testIters	If testIters > 0, check if the boundary is unchanged every testIters iterations.

Fields inherited from class gov.nih.mipav.model.algorithms.AlgorithmBase

destFlag, destImage, mask, maxProgressValue, minProgressValue, multiThreadingEnabled, nthreads, progress, progressModulus, progressStep, runningInSeparateThread, separable, srcImage, threadStopped

Fields inherited from class java.lang.Thread

MAX_PRIORITY, MIN_PRIORITY, NORM_PRIORITY

Constructor Summary

Constructors

Constructor	Description
AlgorithmLevelSet(ModelImage srcImg, float[] sigmas, int movement, int iter, float deltaT, float epsilon, float edgeAttract, int testIters, boolean image25D)	LevelSet.

Method Summary

Modifier and Type	Method	Description
void	algorithmPerformed(AlgorithmBase algorithm)	Called after an algorithm this listener is registered to exits (maybe successfully, maybe not).
private void	calc25D()	calc25D - Calculates level set from contours, propagates level set, and obtains new contour voi from level set.
private void	calc2D()	calc2D - Calculates level set from contours, propagates level set, and obtains new contour voi from level set.
private void	calc3D()	calc3D - Calculates level set from contours, propagates level set, and obtains new contour voi from level set.
void	cleanUp()	DOCUMENT ME!
void	finalize()	finalize - sets class storages arrays to null so that System.gc() can free the memory.
private void	makeKernels2D()	makeKernals2D - creates the derivative kernels used to calculate the gradient magnitude and kernel for the diffusion process.
private void	makeKernels3D()	makeKernals3D - creates the derivative kernels used to calculate the gradient magnitude and kernel for the diffusion process.
void	runAlgorithm()	starts the algorithm.

Methods inherited from class gov.nih.mipav.model.algorithms.AlgorithmBase

actionPerformed, addListener, addProgressChangeListener, calculateImageSize, calculatePrincipleAxis, computeElapsedTime, computeElapsedTime, convertIntoFloat, delinkProgressToAlgorithm, delinkProgressToAlgorithmMulti, displayError, errorCleanUp, fireProgressStateChanged, fireProgressStateChanged, fireProgressStateChanged, fireProgressStateChanged, fireProgressStateChanged, generateProgressValues, getDestImage, getElapsedTime, getMask, getMaxProgressValue, getMinProgressValue, getNumberOfThreads, getProgress, getProgressChangeListener, getProgressChangeListeners, getProgressModulus, getProgressStep, getProgressValues, getSrcImage, isCompleted, isImage25D, isMultiThreadingEnabled, isRunningInSeparateThread, isThreadStopped, linkProgressToAlgorithm, linkProgressToAlgorithm, makeProgress, notifyListeners, removeListener, removeProgressChangeListener, run, setCompleted, setImage25D, setMask, setMaxProgressValue, setMinProgressValue, setMultiThreadingEnabled, setNumberOfThreads, setProgress, setProgressModulus, setProgressStep, setProgressValues, setProgressValues, setRunningInSeparateThread, setSrcImage, setStartTime, setThreadStopped, startMethod, windowActivated, windowClosed, windowClosing, windowDeactivated, windowDeiconified, windowIconified, windowOpened

Methods inherited from class java.lang.Thread

activeCount, checkAccess, clone, countStackFrames, currentThread, dumpStack, enumerate, getAllStackTraces, getContextClassLoader, getDefaultUncaughtExceptionHandler, getId, getName, getPriority, getStackTrace, getState, getThreadGroup, getUncaughtExceptionHandler, holdsLock, interrupt, interrupted, isAlive, isDaemon, isInterrupted, isVirtual, join, join, join, join, ofPlatform, ofVirtual, onSpinWait, resume, setContextClassLoader, setDaemon, setDefaultUncaughtExceptionHandler, setName, setPriority, setUncaughtExceptionHandler, sleep, sleep, sleep, start, startVirtualThread, stop, suspend, threadId, toString, yield

Methods inherited from class java.lang.Object

equals, getClass, hashCode, notify, notifyAll, wait, wait, wait

Field Details

EXPAND

private static final int EXPAND

DOCUMENT ME!

See Also:

• Constant Field Values

EXPAND_CONTRACT

private static final int EXPAND_CONTRACT

DOCUMENT ME!

See Also:

• Constant Field Values

CONTRACT

private static final int CONTRACT

DOCUMENT ME!

See Also:

• Constant Field Values

deltaT

private float deltaT

Time step The defalut is 0.05 for 2D and 0.025 for 3D.

edgeAttract

private float edgeAttract

edgeAttract is 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.

epsilon

private float epsilon

epsilon is used to smooth high curvature regions. epsilon is between 0 and 1 with larger epsilon doing more smoothing

GxData

private float[] GxData

Storage location of the first derivative of the Gaussian in the X direction.

GyData

private float[] GyData

Storage location of the first derivative of the Gaussian in the Y direction.

GzData

private float[] GzData

Storage location of the first derivative of the Gaussian in the Z direction.

iterations

private int iterations

Number of iterations of the diffusion.

kExtents

private int[] kExtents

Dimensionality of the kernel.

movement

private int movement

EXPAND, EXPAND_CONTRACT, or CONTRACT.

sigmas

private float[] sigmas

Standard deviations of the gaussian used to calculate the kernels.

testIters

private int testIters

If testIters > 0, check if the boundary is unchanged every testIters iterations.

image25D

private boolean image25D

If image25D is true in 3D images, process each slice separately

outputBuffer

private float[] outputBuffer

Stores output of AlgorithmConvolver

Constructor Details

AlgorithmLevelSet

public AlgorithmLevelSet(ModelImage srcImg,
 float[] sigmas,
 int movement,
 int iter,
 float deltaT,
 float epsilon,
 float edgeAttract,
 int testIters,
 boolean image25D)

LevelSet.

Parameters:

srcImg - reference to the source image

sigmas - sigmas used to describe the gaussian that is used in the calculation of the gradient magnitude

movement - EXPAND, EXPAND_CONTRACT, or CONTRACT

iter - number of iterations (t) of the diffusion equation

deltaT - time step

epsilon - smoothing factor between 0 and 1 with higher values giving greater smoothing

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

testIters - If testIters > 0, check if the boundary is unchanged every testIters iterations

image25D - If truein 3D images, process each slice separately

Method Details

cleanUp

public void cleanUp()

DOCUMENT ME!

finalize

public void finalize()

finalize - sets class storages arrays to null so that System.gc() can free the memory.

Overrides:

finalize in class AlgorithmBase

runAlgorithm

public void runAlgorithm()

starts the algorithm.

Specified by:

runAlgorithm in class AlgorithmBase

calc2D

private void calc2D()

calc2D - Calculates level set from contours, propagates level set, and obtains new contour voi from level set.

calc25D

private void calc25D()

calc25D - Calculates level set from contours, propagates level set, and obtains new contour voi from level set. Calculates one slice at a time in a multislice image

calc3D

private void calc3D()

calc3D - Calculates level set from contours, propagates level set, and obtains new contour voi from level set.

makeKernels2D

private void makeKernels2D()

makeKernals2D - creates the derivative kernels used to calculate the gradient magnitude and kernel for the diffusion process.

makeKernels3D

private void makeKernels3D()

makeKernals3D - creates the derivative kernels used to calculate the gradient magnitude and kernel for the diffusion process.

algorithmPerformed

public void algorithmPerformed(AlgorithmBase algorithm)

Description copied from interface: AlgorithmInterface

Called after an algorithm this listener is registered to exits (maybe successfully, maybe not). If the algorithm is run in a separate thread, this call will be made within that thread. If not, this call will be made from that same, shared thread.

Specified by:

algorithmPerformed in interface AlgorithmInterface

Parameters:

algorithm - the algorithm which has just completed