[Insight-users] Nasty bug in Canny level set + another nasty bug??

[Joshua Cates]

Hi Zachary,

Thanks for your work and analysis on this issue.  Rather than modifying
"ReverseExpansionDirection", my suggestion is to deprecate its use
entirely and (eventually) force users to explicitly reverse the
propagation and/or advection signs where appropriate for their sign
conventions and methods.  The bottom line is that this method adds
unneccessary complexity and confusion to the level-set function. This is
not the first time it has been an issue.  I think, for example, that when
it was first introduced, "ReverseExpansionDirection" existed in the form
you suggest, i.e.  without the flip of the advection sign.

Note that some groups currently do make use of
"ReverseExpansionDirection", so we will have to be careful in phasing this
method out or in making any changes in the short term.

What do others think about this?

Josh.

______________________________
 Josh Cates			
 Scientific Computing and Imaging Institute
 University of Utah
 http://www.sci.utah.edu/~cates


On Wed, 16 Feb 2005, Zachary Pincus wrote:

> Hello,
> 
> I think I've discovered another problematic bug in the ITK level sets  
> implementation. In answering Hiraki Hideaki's question about whether  
> the advection direction problems were the result of whether inside  
> values were assumed to be positive or negative (it's not; see below), I  
> think I found a problem with ReverseExpansionDirection().
> 
> Now, the advection direction is uninfluenced by whether the "inside" of  
> the level set is assumed to be the positive or negative values of the  
> phi function. One can verify this by referring to the figures I made (  
> http://www.stanford.edu/~zpincus/GeodesicLevelSet.pdf and  
> http://www.stanford.edu/~zpincus/CannyLevelSet.pdf ) and noting that  
> nowhere does the question of "which side is inside" come into play. The  
> Geodesic case works and the Canny case does not regardless of your  
> convention of positive and negative.
> 
> You can also verify this by inverting the input level sets in the test  
> examples I provided and noting that they have the exact same behavior  
> to the un-inverted version. When advection is considered alone, apart  
> from propagation, the question of inside versus outside isn't relevant:  
> the only matter of importance is which way the zero-level set moves  
> with respect to image edges.
> 
> However, calling ReverseExpansionDirectionOn() on the level set filter  
> "fixed" the problem, which puzzled me, since in my test examples I had  
> set the propagation coefficient to zero, so which way the expansion  
> direction is set should have no influence. I investigated, and found  
> that this function causes the advection term to be negated.
> 
> I believe this behavior is entirely incorrect. As above, the advection  
> direction is entirely decoupled from whether the inside is assumed  
> positive or negative. The speed image (and hence the "expansion"  
> direction) is *not* decoupled from this assumption. Thus, to reverse  
> the expansion direction, it is correct to negate the propagation  
> coefficient, but it is *not* correct to negate the advection  
> coefficient. Negating the advection coefficient ONLY serves to drive  
> the zero-level set away from edges (provided the advection image is  
> properly defined!)
> 
> I thus strongly believe that this is a second nasty bug in the ITK  
> level sets implementation. ReverseExpansionDirectionOn() called on the  
> level set segmentation filter (which causes ReverseExpansionDireciton()  
> of the level set function object to be called) should not cause the  
> advection term to be negated.
> 
> As such, lines 24 - 30 of itkSegmentationLevelSetFunction.txx should be  
> changed from:
> 
> template <class TImageType, class TFeatureImageType>
> void SegmentationLevelSetFunction<TImageType, TFeatureImageType>
> ::ReverseExpansionDirection()
> {
>    this->SetPropagationWeight( -1.0 * this->GetPropagationWeight() );
>    this->SetAdvectionWeight( -1.0 * this->GetAdvectionWeight() );
> }
> 
> to
> 
> template <class TImageType, class TFeatureImageType>
> void SegmentationLevelSetFunction<TImageType, TFeatureImageType>
> ::ReverseExpansionDirection()
> {
>    this->SetPropagationWeight( -1.0 * this->GetPropagationWeight() );
> }
> 
> Zach
> 
> 
> On Feb 16, 2005, at 1:59 AM, HIRAKI Hideaki wrote:
> 
> > Hello Zach,
> >
> > Thank you for providing good illustration and test code.
> > I think the advection direction is reversed because your
> > InputLevelsets are following the negative inside convention.
> >
> > You might be interested in our discussion on "Level Set methods":
> > http://public.kitware.com/pipermail/insight-users/2003-June/ 
> > thread.html#4018
> > I hope you would clear the complication.
> >
> > Regards,
> >
> > Hideaki Hiraki
> >
> >
> > At Tue, 15 Feb 2005 22:27:41 -0800, Zachary Pincus wrote:
> >> I think I've found a problem in the Canny level set segmentation
> >> classes that cause the level set to actually *flee* from the Canny
> >> edges rather than be attracted to them.
> >>
> >> I have made some simple test cases which illustrate this problem.  
> >> These
> >> tests can be made to work properly by implementing a simple fix in the
> >> code (see below), or just setting the "AdvectionWeight" term to a
> >> negative value instead of a positive value. The test cases (and movies
> >> of correct and incorrect level set evolution) are at:
> >> http://www.stanford.edu/~zpincus/CannyLevelSetTester.zip [2.44 MB]
> >>
> >> Basically, the issue is that the advection term calculated by
> >> itkCannySegmentationLevelSetFunction.txx is the negative of what it
> >> needs to be to work properly. This advection term is made by taking  
> >> the
> >> distance map from the Canny edges of an image and multiplying it by  
> >> its
> >> gradient. If you think through the implications, this creates an
> >> advection image where the vectors point away from the Canny edges.  
> >> This
> >> causes the level set to be drawn away from these images.
> >>
> >> This all took me a while to think through, so I made some diagrams to
> >> help clarify the issue. Here are PDFs illustrating the ITK
> >> implementation of the GeodesicLevelSet advection term (correct) and  
> >> the
> >> Canny advection term (incorrect, I think):
> >> http://www.stanford.edu/~zpincus/GeodesicLevelSet.pdf
> >> http://www.stanford.edu/~zpincus/CannyLevelSet.pdf
> >>
> >> Stepping through the diagram for the Geodesic case:
> >> (1) shows a 1-dimensional "image" with a step edge.
> >> (2) shows a typical speed image (a bounded reciprocal or sigmoid of a
> >> smoothed derivative of the image)
> >> (3) shows the gradient/derivative (same thing in one dimension) of the
> >> speed image.
> >> (4) shows the negative gradient of the speed image. This is what is
> >> used by the GeodesicActiveContourLevelSetFunction class as the
> >> advection image.
> >> (5) shows the advection term of the ITK level set equation (eq. 9.3
> >> from the ITK software guide)
> >> (6) shows the A(x) term with a phi function (blue dashed line)
> >> superimposed. Note that the gradient of phi is positive, as is A(x)
> >> over the "narrow band", so assuming a positive alpha (which is the
> >> usual case), d/dt(phi) is therefore negative. This means that after  
> >> the
> >> next time step, the phi function is translated downward, moving the
> >> zero-level set toward the edge, as illustrated in (7). It is left as a
> >> trivial exercise to the reader to determine that for other  
> >> orientations
> >> and positions of phi, the zero-level set is always translated toward
> >> the edge.
> >>
> >> The corresponding diagram for the Canny case should make it quite  
> >> clear
> >> what the problem is:
> >> (1) shows the 1-d image
> >> (2) shows an impulse canny edge
> >> (3) shows the distance map to the edge
> >> (4) shows the gradient of this distance map
> >> (5) shows the distance map times the gradient: this is A(x) for the
> >> Canny level set filter.
> >> (6) is the advection term from the level set equation
> >> (7) shows A(x) with a phi function. As above, d/dt(phi) can be
> >> calculated; in this case it is positive, which draws the zero-level  
> >> set
> >> *away* from the edge, as shown in (8). This is the case for any
> >> position/orientation of phi.
> >>
> >> Fortunately, the problem can be fixed by changing line 74 of
> >> itkCannySegmentationLevelSetFunction.txx from:
> >>      it.Set(it_a.Get());
> >> to
> >>      it.Set(it_a.Get() * -1.0);
> >>
> >> In the test cases I constructed after identifying this problem, I  
> >> could
> >> see the level set pushed away from the canny edges. When I made the
> >> above modification (or set the advection weight alpha to a negative
> >> value, which has the same effect), the level set was clearly drawn to
> >> the Canny edges. Again, the test cases can be found here:
> >> http://www.stanford.edu/~zpincus/CannyLevelSetTester.zip
> >>
> >> Zach Pincus
> >>
> >> Department of Biochemistry and Program in Biomedical Informatics
> >> Stanford University School of Medicine
> >>
> > _______________________________________________
> > Insight-users mailing list
> > Insight-users at itk.org
> > http://www.itk.org/mailman/listinfo/insight-users
> >
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>