[Insight-users] Re: Level Set methods
Luis Ibanez
luis . ibanez at kitware . com
Sun, 15 Jun 2003 13:41:46 -0400
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
>