[Insight-users] Re: Level Set methods

Joshua Cates cates at sci . utah . edu
Fri, 20 Jun 2003 10:56:17 -0600 (MDT)


Hi all,

I've been out of town for several days so I've missed much of this
discussion, but as I understand it, there is some confusion about the sign
of the terms of the level-set equation.  Note that it is really the
*relative* sign of the terms with respect to each other that is important.  
The advection and propagation terms should both be positive (negative)  
and the curvature term should be negative (positive).  My intention
with the "UseNegativeFeatures" parameter ( yes, perhaps it was poorly
named ) was to reverse the signs of the terms and therefore effectively
reverse the positive-inside condition to negative-inside.  If anyone has 
experienced problems with this please let me know ASAP.

When analyzing the numerical implementation, take note of the direction in
which finite difference derivatives are computed.  This affects the
ultimate sign of the terms and may be another source of confusion.

All of our experience thus far indicates that the numerical implementation
of the basic level set equation is correct.  Major instabilities should
result if the signs of the terms are not correct. (Try flipping one and
you'll see what I mean.) The sign of the curvature term in the
documentation, however, was in error.  Thanks to Luis for correcting that!


Josh.


______________________________
 Josh Cates			
 School of Computer Science	
 University of Utah
 Email: cates at sci . utah . edu
 Phone: (801) 587-7697
 URL:   http://www . sci . utah . edu/~cates


On Fri, 20 Jun 2003 hhiraki at lab . nig . ac . jp wrote:

> Hi Luis,
> 
> Thanks for the fix. I'm sorry my comments were hasty.
> Still there may be misunderstanding.
> 
> > 1) I agree with you about equation 8.3.
> > 
> >     It has been fixed now.
> >     The negative signs were added to
> >     the advection and propagation terms.
> 
> The fixed equation 8.3 seems to follow the "negative inside" 
> convention. I think, if in the "positive inside" way, the 
> advection term to be negative, the propagation term to be 
> positive, the curvature term to be negative. 
> 
> 
> > 2) You are right also about the two filters
> > 
> >     ShapeDetection
> >     GeodesicActiveContours
> 
> As we have agreed, FastMarching, ShapeDetection and 
> GeodesicActiveContours are following the convention of 
> negative inside. The example of ThresholdLevelSet is also 
> following the convention of negative inside. It is shown 
> in two lines in the file 
> Examples/Segmentation/ThresholdSegmentationLevelSetImageFilter.cxx:
>  227:  thresholdSegmentation->SetUseNegativeFeaturesOn();
>  242:  thresholdSegmentation->SetInput( fastMarching->GetOutput() );
> 
> The examples left are of Canny-edgeLevelSet and LaplacianLevelSet. 
> These are actually following the convention of positive 
> inside. I assume the reason why these are working correct 
> as follows.
> 
> On the curvature term, as the "curvature_term" in itkLevelSetFuntion.txx 
> is calculated from the derivatives of the levelset function, 
> this term changes its sign automatically on switching the 
> "positive/negative inside" conventions. If "kappa" in equation 
> 8.3 means the curvature of the segmented region, the 
> "curvature_term" in ComputeUpdate() is "- kappa" when following 
> the convention of positive inside. The calculations 
> in the file itkLevelSetFunction.txx seems to be assuming 
> the convention of negative inside.
> 
> On the propagation term, the sign change seems to be 
> realized in GeneratateData() in itkSegmentationLevelSetImageFilter.txx.
> If m_UseNegativeFeatures == false (this is the default 
> "positive inside"), the propagation weight "beta" is 
> forced to be negative. Then as the propagation term in 
> the fixed eq. 8.3 is positive (where P(x)>0), the region 
> is expanding.
> 
> On the advection term, the sign is forced to be negative 
> in itkSegmentationLevelSetImageFilter.txx, too. This seems 
> to be wrong. But it doesn't affect on LaplacianLevelSet as 
> this term is not used. I haven't understood Canny-edgeLevelSet 
> well. The advection image seems to be the weighted gradient 
> of distance from the Canny-edges (equation 8.5). As the 
> gradient points away from the edges and is the reverse of 
> the advection image in GeodesicActiveContours, the fault 
> in itkSegmentationLevelSetImageFilter.txx may be saved.
> 
> Could you try GeodesicActiveContour with high advection 
> weight following the convention of positive inside? I'm 
> afraid the contour might be expeled from the edges.
> 
> 
> IMHO, the name "UseNegativeFeatures" is misleading. One 
> possible implementation is using the signs of the three 
> weights as is the user set without the UseNegativeFeatures 
> trick and the required consistency of the signs is 
> stated in the documentations.
> 
> 
> > 3) The ReinitializeLevelSetImageFilter
> >     is not embedded in the LevelSet framework.
> >     It is rather a helper class.
> 
> > 4) Here is an example on how the Reinitialize
> >     filter may use used. (Note that it expects
> >     the input to be a binary image with values
> >     -0.5 and 0.5.)
> <snip>
> >    typedef typename DistanceFilterType::NodeContainerPointer
> >                                                NodeContainerPointer;
> > 
> >    NodeContainerPointer nodes =  m_DistanceFilter->GetOutputNarrowBand();
> 
> How should I set this "nodes" to the LevelSet framework? 
> 
> 
> Thanks,
> 
> Hideaki Hiraki
> 
> 
> From: Luis Ibanez <luis . ibanez at kitware . com>
> Subject: Re: [Insight-users] Re: Level Set methods
> Date: Thu, 19 Jun 2003 09:55:42 -0400
> > Hi Hideaki,
> > 
> > Thanks for your detailed comments,
> > 
> > 1) I agree with you about equation 8.3.
> > 
> >     It has been fixed now.
> >     The negative signs were added to
> >     the advection and propagation terms.
> > 
> > 
> > 2) You are right also about the two filters
> > 
> >     ShapeDetection
> >     GeodesicActiveContours
> > 
> >     The are following the convention of
> >     negative inside, as opposed to the
> >     other level set filters. This also
> >     changes the way in which equation
> >     8.3 may be interpreted. We could
> >     remove the line
> >     this->SetUseNegativeFeatures( true );
> >     from the constructors, and rather
> >     make the point in the documentation...
> > 
> > 3) The ReinitializeLevelSetImageFilter
> >     is not embedded in the LevelSet framework.
> >     It is rather a helper class.
> > 
> >     As you pointed out, it is a Distance
> >     transform implemented usin FastMarching.
> >     The advantage of this filter with respect
> >     to the DanielssonMap is that it produces
> >     a signed distance to the edge of a binary
> >     image.  The DanielssonMap will only produce
> >     a positive distance to the binary object
> >     itself.
> > 
> >     You may want to use the Reinitialize filter for
> >     removing computation artifacts from the level
> >     set every so many iterations. (but you will have
> >     to do it at the application level, not as part
> >     of the normal ITK pipeline).
> > 
> >     If you look at any of the level set images,
> >     after a couple hundred iterations, you will
> >     notice a good number of rectilinear artifacs,
> >     in particular in the regions far from the zero
> >     set. These artifacts cumulate over time and
> >     degrade the underlying assumption that the
> >     \phi function provides a linear approximation
> >     in the neighbohood of the zero set.
> > 
> > 
> > 4) Here is an example on how the Reinitialize
> >     filter may use used. (Note that it expects
> >     the input to be a binary image with values
> >     -0.5 and 0.5.)
> > 
> >    typedef typename itk::ReinitializeLevelSetImageFilter<
> >                                              RealImageType >
> >                                                      DistanceFilterType;
> >    m_DistanceFilter = DistanceFilterType::New();
> >    m_DistanceFilter->SetInput( inputBinaryImage );
> >    m_DistanceFilter->NarrowBandingOn();
> >    m_DistanceFilter->SetNarrowBandwidth( m_BandWidth );
> > 
> >    m_DistanceFilter->Update();
> > 
> >    typedef typename DistanceFilterType::NodeContainerPointer
> >                                                NodeContainerPointer;
> > 
> >    NodeContainerPointer nodes =  m_DistanceFilter->GetOutputNarrowBand();
> > 
> > 
> > 
> > Regards,
> > 
> > 
> >     Luis
> > 
> > 
>