[Insight-users] Re: Level Set methods

[hhiraki@lab.nig.ac.jp]

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