Package gov.nih.mipav.model.algorithms

Class AlgorithmLevelSet

  • All Implemented Interfaces:
    AlgorithmInterface, java.awt.event.ActionListener, java.awt.event.WindowListener, java.lang.Runnable, java.util.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

    • Field Detail

      • 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 Detail

      • 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 Detail

      • 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

      • 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