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

[Zachary Pincus]

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
>