[Insight-users] Contributions for a few enhancements

[lucantiga at softhome.net]

Hi Louis,
 long message ahead!
 I've got a few enhancements I'd like to contribute to ITK. I describe
them below. All the changes are backward-compatible. If you're interested, 
I can provide you with the code.
 
1. GeodesicActiveContourLevelSetFunction::CalculateAdvectionImage now
 uses GradientRecursiveGaussianImageFilter to perform derivatives. This
 is fine, but sometimes, when you have tiny objects like small vessels, 
 advection-driven level sets collapse because of the Gaussian filtering
 introduced in performing derivatives. The enhacement in this class is
 to use GradientImageFilter to perform derivatives when DerivativeSigma
 is set to zero. 
 
2. I've extended FastMarchingImageFilter to do two additional things:
 a. to accept a set of TargetPoints in input. After one or all
 (depending on user's choice) the target points are reached, the
 evolution stops. A variable TargetOffset controls how much time after
 targets are reached the evolution must stop (so to ensure a smooth level 
 set corresponding to the time the target is reached).
 b. to optionally compute the upwind gradient vector each time a Trial 
 point becomes Accepted and to store it in a GradientImage. This
 basically generates a vector field of reciprocal of velocities (i.e. the
 argument of the modulus in the Eikonal equation). It's not the same
 as computing the gradient with a GradientImageFilter since here we
 use upwind finite differences and we only do it for accepted points.

3. I wrote a new class which I find very useful for level sets
 initialization. It's called CollidingFrontsImageFilter. Given two sets
 of seeds, the filter propagates two fronts with the
 FastMarchingImageFilter extension described above with the input
 image as speed, computing upwind gradients and stopping when the other set is met. 
 A scalar product between the two upwind gradient images is then taken.
 The region where such scalar product is negative (and optionally flood-fill
 connected to the seeds) is the region where the two fronts collide one 
 against the other. In ascii-art: .)))) ((((. that is, the region
 "between" the two points. The regions where the scalar product is
 positive can be ignored, in ascii-art: ((((.   .))))
 This is really handy, for example, if you have a bright blood image of a
 branching vessel. Placing two seeds in two points of the vascular
 tree, the filter will extract the region comprised between the
 seeds ignoring side branches (because both fronts will seep into side
 branches in the same direction, i.e. positive scalar product).
 Morevoer, if snr is good, both fronts will exit the vessel boundary
 in the same direction (since the front there is slow), so one doesn't
 even need to set propagation speed to zero outside the vessel.
 
Just let me know if you're interested. If not, I'll just release these
classes with my library (and maybe move them to itk at a later time)

Last, a question: 
 the code for CurvesLevelSetFunction is basically the same 
 as GeodesicActiveContourLevelSetFunction. However, Curves is the name of Lorigo's 
 method to compute codimension-two geodesic active contours, which is
 not what the class does. Couriously, the txx file includes
 "vnl/algo/vnl_symmetric_eigensystem.h" which is what Lorigo's method
 needs for codimension-two geodesic active contours.
 Was CurvesLevelSetFunction originally supposed to implement
 codimension-two geodesic active contour?

Sorry for the lengthy message!

Luca


-- 
Luca Antiga, PhD
----------------------------------------------
Biomedical Technologies Laboratory
Bioengineering Department, 
Mario Negri Institute
Villa Camozzi, 24020, Ranica (BG), Italy
----------------------------------------------
phone: +39 035 4535-381
email: antiga at marionegri.it
web:   http://villacamozzi.marionegri.it/~luca
----------------------------------------------