[Insight-developers] Nasty bug in Canny level set segmentation (?)

[Zachary Pincus]

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