[Insight-users] Fwd: problem with ParallelSparseFieldLevelSetImageFilter - Fix

[Luca Antiga]

Hey,
  so I looked into the issue with the  
ParallelSparseFieldLevelSetImageFilter

Luis, I need your wisdom here (sorry for the long message ahead).

As the documentation for itk::LevelSetFunction says:
  * \warning You MUST call Initialize() in the constructor of  
subclasses of this
  * object to set it up properly to do level-set Calculations.  The  
argument
  * that you pass Initialize is the radius of the neighborhood needed  
to perform
  * the calculations.  If your subclass does not do any additional  
neighborhood
  * processing, then the default radius should be 1 in each direction.

Now, this is pretty high up in the hierarchy, so the message is very  
likely not to be read by users that are doing what Kun rightfully  
wants to do.

To me, this call
>       ThresholdFunctionType::RadiusType    r;
>       r[0] = r[1] = r[2] =1;
>       m_LevelSetFunction->Initialize(r);

shouldn't be required when using SparseFieldImageFilter directly, as  
in Kun's case.

In usual level set filters, the Initialize(r) method is called in
itk::SegmentationLevelSetImageFilter::SetSegmentationFunction,
which does the following

   virtual void SetSegmentationFunction(SegmentationFunctionType *s)
   {
     m_SegmentationFunction = s;

     typename SegmentationFunctionType::RadiusType r;
     r.Fill( 1 );

     m_SegmentationFunction->Initialize(r);
     this->SetDifferenceFunction(m_SegmentationFunction);
     this->Modified();
   }

One drawback is that initialize is called with radius of 1, and it a  
function subclass needs a radius of 2 it has to call Initialize(r)  
again but forcibly after the SetSegmentationFunction call, which it's  
not clean (and not documented).

Now, to me a concrete level set function should know what  
neighborhood radius it uses by default and should be able to set it  
up independently from the image filter, it shouldn't be up to the  
image filter to set it.
In fact, the documentation for LevelSetFunction says that Initialize  
should be called in the constructor of subclasses.

One option could be to call Initialize(r) in the constructor of all  
concrete subclasses and remove it from SetSegmentationFunction
(We can't call it once for all in the constructor of  
itk::SegmentationLevelSetFunction because Initialize is virtual and  
it's being overridden in the subclasses).

This is a backward incompatible change, because people who coded  
custom level set functions should also add the call to Initialize(r)  
in the constructor.

There is a backward compatible option:
- to introduce a Initialize() method (without the radius as an  
argument) that uses the m_Radius member in FiniteDifferenceFunction.
- set m_Radius to some proper value in the constructors of subclasses  
(or default it to 1 in the constructor of  
itk::SegmentationLevelSetFunction)
- modify SetSegmentationFunction in this way:
  	if m_Radius is zero in all directions (the way it's set in the  
constructor of FiniteDifferenceFunction), use the Initialize(r)  
signature with radius 1 the way it is now;
	else, use Initialize(), which uses m_Radius;

This way we keep backward compatibility and give level set functions  
a chance to have a radius != 1.

In this case, the user is required to call function->Initialize()  
before using it, but at least he/she doesn't deal with radius.

Last, we could make itk::LevelSetFunction throw an exception if it's  
being used uninitialized (which would require the addition of a  
initialized state flag), in order to avoid plain crashes like the one  
Kun reported.

I hope the above is intelligible... it's easier to do it than to  
write about it for sure :-)


Luca



On Jul 30, 2008, at 1:10 AM, Kun wrote:

>  Thanks, Luca,
>
>     I have found the problem. The LevelSetFunction should be  
> initialized first,otherwise the radius will set 0, then such an  
> error happens. Also, it seems that the SpeedImage should be  
> calculated before sending to the parallelfilter.
>
>     The code as follows:
>
>       m_LevelSetFunction    =     ThresholdFunctionType::New();
>       m_LevelSetFunction    ->    SetUpperThreshold(150);
>       m_LevelSetFunction    ->    SetLowerThreshold(100);
>       m_LevelSetFunction    ->    SetFeatureImage(SourceImage);
>
>       ThresholdFunctionType::RadiusType    r;
>       r[0] = r[1] = r[2] =1;
>       m_LevelSetFunction->Initialize(r);
>
>       m_LevelSetFunction->AllocateSpeedImage();
>       m_LevelSetFunction->CalculateSpeedImage();
>
>       parallelfilter    =     ParallelSparseFieldFilterType::New();
>       parallelfilter    ->    SetInput(InitImage);
>       parallelfilter    ->    SetNumberOfLayers(3);
>       parallelfilter    ->    SetIsoSurfaceValue(0.0);
>       parallelfilter    ->    SetDifferenceFunction 
> (m_LevelSetFunction);
>       parallelfilter    ->    Update();
>
>
>     Appreciate and Thanks all.
>
> Kun
>
> From: Luca Antiga <luca.antiga at gmail.com>
> Date: Tue, Jul 29, 2008 at 10:59 AM
> Subject: Re: [Insight-users] problem with  
> ParallelSparseFieldLevelSetImageFilter
> To: ITK Users <insight-users at itk.org>
>
>
> Dear Kun,
>  I'll try to reproduce the problem and get back to you.
> Keep in touch
>
> Luca
>
>
>
>
>
>
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users