[Insight-users] Re: Level Set methods
hhiraki@lab.nig.ac.jp
hhiraki at lab . nig . ac . jp
Mon, 16 Jun 2003 16:19:54 +0900
Hi Luis,
Thank you for the detailed answers for my wordy questions. All
my assumptions were from browsing SoftwareGuide without playing
with the codes. I'll try some trivial tests and report the result
if I find unexpected behaviour. Before that, I have two comments.
> Some ITK methods are to be interpreted as
> generating \phi() negative inside. For example,
> FastMarching
This seems to be a source of confusion. What are the filters
actually follow the positive inside convention? ShapeDetection in
section 8.3.2, for example, uses FastMarching's output as its input
without ReinitializeLevelSetImageFilter. This means the negative
inside there, doesn't it?
> Note that these methods do not use millimeters.
> Image spacing is not taken into account in the
> PDE solver. If you have strongly anisotropic
> pixels you may want to resample your data before
> using the level set filters.
Thanks for pointing this out. It counters my intuition using
properties like curvature without taking into account spacing.
I hope those might help more clarifying the documents.
Regards,
Hideaki Hiraki
Luis Ibanez <luis . ibanez at kitware . com> wrote:
>
> Hi Hideaki,
>
> Thanks for your detailed questions regarding the ITK
> level set methods. Please find comments below.
>
>
> 1) Yes, equation 8.3 in the SoftwareGuide is to
> be integrated over time ("t").
>
> \phi(X,t) is the explicit function in which
> the zero set (\phi(X,t)==0) is embebded as an
> implicit function.
>
>
> 2) The definition on whether the LevelSet is
> positive inside and negative outside has been
> a source of (friendly) discussion for a while.
>
> The decision is totally arbitrary, so in some
> way it is a matter of tossing a coin and sticking
> to the outcome. Using positive inside is convenient
> when you use these methods in combination with
> visualization systems, since you can easily
> perform volume rendering on the data set.
>
> Some ITK methods are to be interpreted as
> generating \phi() negative inside. For example,
> FastMarching which is actually generating a
> time-crossing map from a set of seeds. If we
> assume the seeds to be inside the object to
> be segmented, then the outcome of FastMarching
> will be a \phi(x,t), negative inside the object.
>
> However you may object that the time-crossing
> map shouldn't be used as a level set directly,
> and instead it must be thresholded and passed
> through the ReinitializeLevelSetImageFilter
> http://www . itk . org/Insight/Doxygen/html/classitk_1_1ReinitializeLevelSetImageFilter . html
> in order to generate a formal level set with
> \phi positive inside.
>
> 3) The values of \phi are not necesarity equal to
> the distance to the zero set. This is only
> the case when you reinitialize the level set
> with the ReinitilizeLevelSetImageFilter.
> Otherewise, after a couple of iterations of
> the PDE solver, \phi will just approximate
> such distance, close to the zero set.
>
> Note that these methods do not use millimeters.
> Image spacing is not taken into account in the
> PDE solver. If you have strongly anisotropic
> pixels you may want to resample your data before
> using the level set filters.
>
> Take a look at the ResamplingImageFilter in
> http://www . itk . org/Insight/Doxygen/html/classitk_1_1ResampleImageFilter . html
> Section 5.7.1, pdf-pag 142 of the SoftwareGuide
> http://www . itk . org/ItkSoftwareGuide . pdf
>
> 4) Positive values of d phy / dt, at a point means
> that the zero set is expanding at that point.
> (Again, under the assumption of 'positive inside').
>
>
> 5) nabla phi is the normal vector to the zero set.
> Note that this is not necessarily a unit vector,
> and, as you said, points towards the inside of the
> object (due to the 'positive inside').
>
> 6) abs(nabla phi) is 1 on the boundaries only after
> using the Reinitialize filter. After that this
> condition is only an approximation, since the
> PDE solver iteration will not necesarily produce
> pure advection. Specially in regions of high
> curvature.
>
> 7) kappa is the curvature of phi. Negative in
> concave regions, positive in convex regions
> (assuimg the convention of phi = 'positive
> inside')
>
> 8) alpha, beta and gamma in equation 8.3 are
> respectively the weights for the terms
>
> - advection
> - propagation
> - curvature
>
> FastMarching is like
> alpha == 0
> gamma == 1
> so, only propagation term
>
> ShapeDetection is like
> alpha == 0
> so: propagation and curvature terms
>
> GeodesicAcitveContours
> uses the three terms:
> propagation, curvature and advection.
>
>
> 9) These terms are initialized to 1.0
> at the level of the base class, the
> itkSegmentationLevelSetImageFilter
> http://www . itk . org/Insight/Doxygen/html/classitk_1_1SegmentationLevelSetImageFilter . html
>
> alpha = advection scaling = 1.0
> beta = propagation scaling = 1.0
> gamma = curvature scaling = 1.0
>
>
> 10) Yeap, P(x) is set as the feature image.
> It is computed as the sigmoid of the
> gradient magnitude. But this is only
> one possibility among many others.
> Using an exponential filter is another
> common option. The advantage of the
> sigmoid is that it is easy to control
> in order to obtain 1.0 in homogeneous
> regions and close to 0.0 in strong
> edges.
>
>
> 11) The vector map A(x), also called the
> "Advection image" is computed inside
> the GeodesicActiveContours filter using
> the GradientRecursiveGaussian image filter
> over the feature image. The Advection image
> is simply the negative of the gradient of
> the feature map. In this way, the vectors
> of this map point towards the edges of the
> original image. You can see this computation
> in the file.
>
> /Insight/Code/Algorithms/
> itkGeodesicActiveContourLevelSetFunction.txx
>
>
> 12) You can verify the whole computation of
> equation 8.3 (SoftwareGuide.pdf) in the
> class itk::LevelSetFunction, in particular
> in the file
>
> Insight/Code/Common/itkLevelSetFunction.txx
>
> in the method ComputeUpdate().
>
> Note that the contributions of the
> advection term are verified in lines 224,225
>
> The final goal is to move the zero set toward
> the edges. If you find that this equation is
> doing something different, please let us know
> so it can be fixed.
>
>
> 13) When you say that the curvature terms seems
> to have the sign inversed, do you meant it
> in equation 8.3 of the SoftwareGuide ?
>
> or in the computations in
>
> Insight/Code/Common/itkLevelSetFunction.txx ?
>
>
>
> It may be worth to run some tests in trivial cases
> where we could determine if all these terms are
> behaving correctly.
>
> In any case, so far we have use the GeodesicActive
> contour filter to succesfully segment brain structures
> from 2D images, as well as Abdominal Aortic sections
> from contrast enhanced CT in 3D.
>
> You may want to play with the examples available in
> the Insight/Examples/Segmentation directory, by running
> them over synthetic images. That may prove useful to
> evaluate any hypothesis of wrong signs you may have.
> For example, a simple case is to try segmenting a
> binary triangle. Playing with the curvature scaling
> you should be able to regulate how deep the zero set
> will go into the corners.
>
>
>
> Please let us know if you find any case in which those
> terms are not behaving as expected, we will be happy to
> fix them.
>
>
>
> Thanks
>
>
> Luis
>
>
>
> ----------------------------------------
> hhiraki at lab . nig . ac . jp wrote:
> > Hello Luis,
> >
> > Thank you very much for the good SoftwareGuide.
> >
> > Trying to realize the Geodesic Active Contour Segmentation in ITK, I am
> > in confusion now. Would you please correct my misunderstandings in the
> > following assumptions?
> >
> >
> > The equation 8.3 in the section 8.3 of the SoftwareGuide is numerically
> > integrated on "t". (By the way, the equations in doxygen generated
> manual
> > http://www . itk . org/Doxygen/html/classitk_1_1LevelSetFunction . html are
> > corrupted.) Here "phi(x,t)" is the levelset funtion that is positive
> > inside the segmentated region (contrary to some literatures). The
> > absolute value of "phi(x,t)" is the distance (in millimeter) from the
> > region boundaries. Positive value of the equation means the region is
> > expanding.
> > "nabla phi" is the normal vector of the contours of "phi" and directs
> > toward inside the region. "abs(nable phi)" is 1 on the boundaries.
> > "kappa" is the mean curvature of "phi" at the time. On the boundaries,
> > "kappa" is negative if the region is concave, and positive if convex.
> >
> > The three term in the right hand side is linearly combined with the
> > weights "alpha", "beta" and "gamma". If "alpha"=0, the result is exactly
> > same as the Shape Detection Segmentation. And if "alpha"="gamma"=0, the
> > result is same as the Fast Marching Segmentation except for some errors
> > from floating point computation. The default weights in Geodesic Active
> > Contour Segmentation are "alpha"="beta"="gamma"=1.
> >
> > The scalar map "P(x)" is set as FeatureImage that is gradient magnitude
> > of denoised input image transformed by sigmoid filter. At the boundaries
> > of regions, the gradient magnitude is large and the FeatureImage is
> about
> > zero. Inside regions, the gradient magnitude is small and the
> FeatureImage
> > is about one.
> > The vector map "A(x)" is negative gradient of the FeatureImage. Near the
> > boundaries of regions, the vectors direct toward the boundary because
> the
> > FeatureImage is changing from one to zero.
> > The scalar map "Z(x)" is constant 1 (not used for now).
> >
> > The first term of the r. h. s. is advection term which works to attract
> > the boundaries to the edges (where the FeatureImage is small). But where
> > the edge is just outside of the boundary, as the inner product of "A(x)"
> > and "nabla phi" is negative, the region shrinks. This term seems to work
> > to expel the boundaries from the edges???
> >
> > The second term of the r. h. s. is expansion term which works to inflate
> > the regions constantly and to slow down at the edges.
> >
> > The third term of the r. h. s. is smoothing term which works to minimize
> > the length (or area in 3D) of the boundaries. The sign of this term also
> > seems to be inverse???
> >
> >
> > I am very sorry for my poor English. I would appreciate any
> clarification.
> >
> > Regards,
> >
> > Hideaki Hiraki
> >