[Insight-users] Segmenting Visible Human Data : RGB ConfidenceConnected.

Wed, 17 Dec 2003 14:54:50 -0500

This is a multi-part message in MIME format.
--------------040506090806040102030906
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit

Hi Luis,

I have generated the Vector ThresholdSegmentationLevelSetImageFilter. 
The files are attached. You can add them to the toolkit.

I wanted to write a test too, but I did not figure out how I could do 
this. Maybe I could do this by modifying the tests for the 
ThresholdSegmentationLevelSetImageFilter, but I don't know how. If you 
have a good starting point I will try to write this, too.

And there is another question regarding the MahalanobisDistance. In the 
filter I used the square root of the Evaluate method of  the 
MahalanobisDistanceMembershipFunction. Is this reasonable?
Are there any resources available on the Internet, where example values 
are given so that it is possible to make estimations. I think I got some 
good values, but I am interested how it is working. I did not find much 
useful by using google.

Thank you
stefan

Luis Ibanez wrote:

>
> Hi Stefan,
>
> You are right,
> Region growing filters that only based on intensity values
> are prone to producing leaks.
>
> You may reduce this tendency by first applying a smoothing
> filter like the VectorGradientAnisotropic smoothing. This
> may help, but still it is not possible to guarranty that it
> will prevent the leaks.
>
> In practice a common approach is to use the RegionGrowing
> methods for producing a very sketchy representation of the
> object, then solidify this representation using mathematical
> morphology filters (like dilation-erosion sequences). Then
> use this as an initialization for level set filters that
> have better capabilities for dealing with leaks.
>
> The ThresholdLevelSet filter is certainly one of the best
> first options to try out.  Unfortunately this filter is
> not yet extended to RGB data.  In make sense, as you
> suggested, to use the Mahalanobis distance in this case
> for controling the Threholding values.
>
> The good news for you is that is should be relatively easy
> to get this RGB level set filter done. You simply need to
> create a class
>
>
>    itkVectorThresholdSegmentationLevelSetFunction
>
> based on the code of the current
>
>    itkThresholdSegmentationLevelSetFunction
>
>
> whose only goal is to compute the Speed image used by the
> level set.
>
> Once you have this new class, you can easily create the
>
>    VectorThresholdSegmentationLevelSetImageFilter
>
> by copy/pasting code from the current class:
>
>         ThresholdSegmentationLevelSetImageFilter
>
>
>
> Notice that the level set code only sees this speed image,
> not the original RGB data.
>
> The only trick here is that you will have to compute the
> mean and covariance of a "sample region" in order to
> feed the Mahalanobis distance function.  You may also want
> to look at the speed function before you start using the
> level set method.  A quick look at the speed functions
> will give you a feeling on the chances of segmenting the
> region using the level set method. A high quality speed
> image is a fundamental requirement for getting good results
> from level set methods.
>
>
> Please let us know if you have any problems in writing
> this class.  We will also be very interested in adding it
> to the toolkit  :-)
>
>
> Thanks
>
>
>
>    Luis
>
>
>
> -------------------------------
> Stefan Lindenau wrote:
>
>> Hi Luis,
>>
>> ok I have read the parts of the software guide that you mentioned again.
>>
>> Now I want to realize the segmentation of the Visible Human Data by 
>> using the VectorConfidenceConnectedImageFilter to get the mean vector 
>> and the covariant matrix of my tissue. I cannot use the segmentation 
>> of this filter directly because it is leaking.
>> With this data I want to initialize a Levelset filter that is almost 
>> similar to the ThresholdLevelset filter, but it should use the 
>> Mahalanobis distance for generating the speed image.
>>
>> I think that I have to write this LevelsetFilter by myself or is 
>> there a implementation for such a problem available?
>>
>>
>> Thanks
>> Stefan
>>
>> Luis Ibanez wrote:
>>
>>> Hi Stefan,
>>>
>>> When you use ConfidenceConnected you only need to provide the 
>>> multiplier
>>> for the variance.  The range of intensities is computed by the filter
>>> based on the mean and the variance of intensities around the seed
>>> points.
>>>
>>> The range is simply:
>>>
>>>      lower limit =   mean   -   standardDeviation * multiplier
>>>      upper limit =   mean   +   standardDeviation * multiplier
>>>
>>> The mean and standardDeviation are computed by the filter.
>>> You only need to tune the value of the multiplier, and
>>> experiement with the number of iterations.
>>>
>>> This holds for RGB confidence connected, where instead of a scalar mean
>>> you have a mean vector of three components (RGB components), and 
>>> instead
>>> of standardDeviation you have a covariance matrix, intead of lower and
>>> upper limits the filter computes the Mahalanobis distance in RGB space.
>>> Therefore you only need to provide the value for the multiplier.
>>>
>>> You may want to read again the description of this method in the
>>> SoftwareGuide.
>>>
>>>         http://www . itk . org/ItkSoftwareGuide . pdf
>>>
>>> It is in Section 9.1.3.
>>> In particular look at equation 9.2 in pdf-page 348.
>>>
>>> We used the RGB Confidence connected filter for producing most of the
>>> segmentations shown in the cover of the SoftwareGuide printed version.
>>>
>>> The code we used for creating the cover is available in
>>>
>>>      InsightDocuments/SoftwareGuide/Cover/Source
>>>
>>>
>>>
>>> Regards,
>>>
>>>
>>>   Luis
>>>
>>>
>>> ------------------------
>>> Stefan Lindenau wrote:
>>>
>>>> Hi Luis,
>>>>
>>>> I tried to get the example of Josh working but I failed on VC6 and 
>>>> Cygwin to compile it. At the moment I want to give your suggestion 
>>>> with the ConfidenceConnected and the ThresholdConnected filter a try.
>>>> I read the Software Guide and I think that I am now knowing how 
>>>> these filters are working. The only thing that I do not understand 
>>>> is how I can get the intensity range values from the 
>>>> ConfidenceConnected filter. I can get/set the multiplier, but I see 
>>>> no access method to these values.
>>>>
>>>> Maybe I could get them by comparing the input image of the 
>>>> ConfidenceConnectedFilter and the output Image, but this seems a 
>>>> bit to complicated to me. Is there a more elegant solution? Did I 
>>>> miss a method?
>>>>
>>>> Thank you
>>>> Stefan
>>>>
>>>> P.S.: as I have progressed with my work I have seen that the data I 
>>>> need can be reduced to 500MB (unsigned char RGB).
>>>>
>>>> Luis Ibanez wrote:
>>>>
>>>>>
>>>>> Hi Stefan,
>>>>>
>>>>>
>>>>> The reason for postprocessing the joint regions is that
>>>>> if you take two contiguous pieces of the image and run
>>>>> level sets on each one, the level sets will evolve in
>>>>> a different way at each side of the boundary, and it
>>>>> is likely that if you try to put the two level sets
>>>>> together just by joining the two blocks of data, the
>>>>> zero set surface will not be contiguous from one block
>>>>> to the next.
>>>>>
>>>>> I would anticipate that some smoothing will be needed
>>>>> for ironing out any discontinuity in the connection.
>>>>> taking the joint region (a region around the boundary
>>>>> of the two block and running some more iterations of
>>>>> the level set there may help to smooth out the transition
>>>>> between the blocks.
>>>>>
>>>>> You could certainly attempt this post-processing-smoothing
>>>>> with other methods. For example, a simple median filter
>>>>> has proved to be powerful enough for smoothing out
>>>>> transitions and it will be a much faster approach too.
>>>>>
>>>>> You may want to start by trying Josh's suggestion since
>>>>> he and his group are the ones who experimented more
>>>>> deeply into this issue.
>>>>>
>>>>>
>>>>> Please let of know of your findings,
>>>>>
>>>>>
>>>>>  Thanks
>>>>>
>>>>>
>>>>>   Luis
>>>>>
>>>>>
>>>>> -----------------------
>>>>> Stefan Lindenau wrote:
>>>>>
>>>>>> Hi Luis,
>>>>>>
>>>>>> thank you for your quick and comprehensive answer. I will just 
>>>>>> have to cut the image in pieces.
>>>>>>
>>>>>> Only one thing I still do not understand:
>>>>>>
>>>>>>> If you use level sets, you could post process
>>>>>>> the joint regions between your contiguous pieces
>>>>>>> in order to smooth out the potential differences
>>>>>>> between the level set obtained in one region and
>>>>>>> the level set obtained in the facing region.
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> Why is it dependend on the level sets to postprocess the the 
>>>>>> joint region. In my comprehension I will just cut the data in big 
>>>>>> pieces,process it and put it together just after the processing. 
>>>>>> Then such a postprocessing should be possible with any of the 
>>>>>> methods. Or did I ignore some facts?
>>>>>>
>>>>>> Maybe I can get it working with the streaming example for 
>>>>>> watershed algorithms as Joshua proposed. I will just have to test 
>>>>>> it out.
>>>>>>
>>>>>>
>>>>>> thanks
>>>>>> Stefan
>>>>>>
>>>>>> _______________________________________________
>>>>>> 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
>>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> 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
>>
>
>
>
>

--------------040506090806040102030906
Content-Type: text/plain;
 name="itkVectorThresholdSegmentationLevelSetFunction.txx"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="itkVectorThresholdSegmentationLevelSetFunction.txx"

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkThresholdSegmentationLevelSetFunction.txx,v $
  Language:  C++
  Date:      $Date: 2003/09/10 14:28:39 $
  Version:   $Revision: 1.5 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www . itk . org/HTML/Copyright . htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef __itkVectorThresholdSegmentationLevelSetFunction_txx_
#define __itkVectorThresholdSegmentationLevelSetFunction_txx_

#include "itkVectorThresholdSegmentationLevelSetFunction.h"
#include "itkImageRegionIterator.h"
#include "itkGradientAnisotropicDiffusionImageFilter.h"
#include "itkLaplacianImageFilter.h"
#include "itkImageFileWriter.h"

namespace itk {

template <class TImageType, class TFeatureImageType>
void VectorThresholdSegmentationLevelSetFunction<TImageType, TFeatureImageType>
::CalculateSpeedImage()
{

  ImageRegionConstIterator<FeatureImageType>
    fit(this->GetFeatureImage(), this->GetFeatureImage()->GetRequestedRegion());
  ImageRegionIterator<ImageType>
    sit(this->GetSpeedImage(), this->GetFeatureImage()->GetRequestedRegion());

  ScalarValueType threshold;
  for ( fit.GoToBegin(), sit.GoToBegin(); ! fit.IsAtEnd(); ++sit, ++fit)
    {
	  threshold = m_Threshold - sqrt(m_Mahalanobis->Evaluate(fit.Get()));
      sit.Set( static_cast<ScalarValueType>(threshold) );

    }

}

} // end namespace itk

#endif

--------------040506090806040102030906
Content-Type: text/plain;
 name="itkVectorThresholdSegmentationLevelSetImageFilter.h"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="itkVectorThresholdSegmentationLevelSetImageFilter.h"

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkVectorThresholdSegmentationLevelSetImageFilter.h,v $
  Language:  C++
  Date:      $Date: 2003/12/05 14:28:39 $
  Version:   $Revision: 1.00 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www . itk . org/HTML/Copyright . htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef __itkVectorThresholdSegmentationLevelSetImageFilter_h_
#define __itkVectorThresholdSegmentationLevelSetImageFilter_h_

#include "itkSegmentationLevelSetImageFilter.h"
#include "itkVectorThresholdSegmentationLevelSetFunction.h"

namespace itk {

/** \class ThresholdSegmentationLevelSetImageFilter
 *    \brief Segments structures in images based on intensity values.
 *
 *   \par IMPORTANT
 *   The SegmentationLevelSetImageFilter class and the
 *   VectorThresholdSegmentationLevelSetFunction class contain additional information necessary
 *   to the full understanding of how to use this filter.
 *
 *    \par OVERVIEW
 *    This class is a level set method segmentation filter.  It constructs a
 *    speed function which is close to zero where the Mahalabonian Distance
 *    exceeds a certain threshold, effectively locking the propagating front onto those
 *    edges.  Elsewhere, the front will propagate quickly.
 *
 *    \par INPUTS
 *    This filter requires two inputs.  The first input is a seed
 *    image.  This seed image must contain an isosurface that you want to use as the
 *    seed for your segmentation.  It can be a binary, graylevel, or floating
 *    point image.  The only requirement is that it contain a closed isosurface
 *    that you will identify as the seed by setting the IsosurfaceValue parameter
 *    of the filter.  For a binary image you will want to set your isosurface
 *    value halfway between your on and off values (i.e. for 0's and 1's, use an
 *    isosurface value of 0.5).
 *
 *    \par
 *    The second input is the feature image.  This is the image from which the
 *    speed function will be calculated the feature image has to be a Vector Image.
 *    For most applications, this is the
 *    image that you want to segment. The desired isosurface in your seed image
 *    should lie within the region of your feature image that you are trying to
 *    segment. Note that this filter does no preprocessing of the feature image
 *    before thresholding.
 *
 *    \par
 *    See SegmentationLevelSetImageFilter for more information on Inputs.
 *
 *    \par OUTPUTS
 *    The filter outputs a single, scalar, real-valued image.
 *    Positive values in the output image are inside the segmentated region
 *    and negative values in the image are outside of the inside region.  The
 *    zero crossings of the image correspond to the position of the level set
 *    front.
 *
 *   \par
 *   See SparseFieldLevelSetImageFilter and
 *   SegmentationLevelSetImageFilter for more information.
 *
 *   \par PARAMETERS
 *   In addition to parameters described in SegmentationLevelSetImageFilter,
 *   this filter adds the Threshold, the Mean and the Covariance.  See
 *   VectorThresholdSegmentationLevelSetFunction for a description of how this value
 *   affect the segmentation.
 *
 *   \sa SegmentationLevelSetImageFilter
 *   \sa ThresholdSegmentationLevelSetFunction,
 *   \sa SparseFieldLevelSetImageFilter */
template <class TInputImage,
          class TFeatureImage,
          class TOutputPixelType = float >
class ITK_EXPORT VectorThresholdSegmentationLevelSetImageFilter
  : public SegmentationLevelSetImageFilter<TInputImage, TFeatureImage, TOutputPixelType, Image<TOutputPixelType, ::itk::GetImageDimension<TInputImage>::ImageDimension> >
{
public:
  /** Standard class typedefs */
  typedef VectorThresholdSegmentationLevelSetImageFilter Self;
  typedef  SegmentationLevelSetImageFilter<TInputImage, TFeatureImage, TOutputPixelType, Image<TOutputPixelType, ::itk::GetImageDimension<TInputImage>::ImageDimension> > Superclass;
  typedef SmartPointer<Self>  Pointer;
  typedef SmartPointer<const Self>  ConstPointer;

  /** Inherited typedef from the superclass. */
  typedef typename Superclass::ValueType ValueType;
  typedef typename Superclass::OutputImageType OutputImageType;
  typedef typename Superclass::FeatureImageType FeatureImageType;

  /** Type of the segmentation function */
  typedef typename VectorThresholdSegmentationLevelSetFunction<OutputImageType,FeatureImageType> ThresholdFunctionType;
  typedef typename ThresholdFunctionType::Pointer ThresholdFunctionPointer;
  typedef typename ThresholdFunctionType::MeanVectorType MeanVectorType;
  typedef typename ThresholdFunctionType::CovarianceMatrixType CovarianceMatrixType;
  typedef typename ThresholdFunctionType::ScalarValueType ScalarValueType;

  /** Run-time type information (and related methods). */
  itkTypeMacro(VectorThresholdSegmentationLevelSetImageFilter, SegmentationLevelSetImageFilter);

  /** Method for creation through the object factory */
  itkNewMacro(Self);

  /** Set/Get mean and covariance that will be used to calculate the speed function */
  void SetMean(const MeanVectorType &mean) 
  {	
	  m_ThresholdFunction->SetMean(mean);
	  this->Modified();
  }
  const MeanVectorType & GetMean() const 
  {	
	  return m_ThresholdFunction->GetMean();
  }

  void SetCovariance(const CovarianceMatrixType &cov) 
  { 
	  m_ThresholdFunction->SetCovariance(cov);
	  this->Modified();
  }
  const CovarianceMatrixType & GetCovariance() const
  { 
	  return m_ThresholdFunction->GetCovariance();
  }

  /** Set/Get the threshold for the Mahanalobis Distance */
  void SetThreshold(ScalarValueType thr) 
  {
	m_ThresholdFunction->SetThreshold(thr);
	this->Modified();
  }
  ScalarValueType GetThreshold ()
  {
	return m_ThresholdFunction->GetThreshold();
  }

protected:
  ~VectorThresholdSegmentationLevelSetImageFilter() {}
  VectorThresholdSegmentationLevelSetImageFilter();

  virtual void PrintSelf(std::ostream &os, Indent indent) const; 

  ThresholdSegmentationLevelSetImageFilter(const Self &); // purposely not impl.
  void operator=(const Self&); //purposely not implemented
private:
  ThresholdFunctionPointer m_ThresholdFunction;
};

} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkVectorThresholdSegmentationLevelSetImageFilter.txx"
#endif

#endif

--------------040506090806040102030906
Content-Type: text/plain;
 name="itkVectorThresholdSegmentationLevelSetImageFilter.txx"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="itkVectorThresholdSegmentationLevelSetImageFilter.txx"

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkVectorThresholdSegmentationLevelSetImageFilter.txx,v $
  Language:  C++
  Date:      $Date: 2003/09/10 14:28:39 $
  Version:   $Revision: 1.11 $

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef __itkVectorThresholdSegmentationLevelSetImageFilter_txx_
#define __itkVectorThresholdSegmentationLevelSetImageFilter_txx_

#include "itkVectorThresholdSegmentationLevelSetImageFilter.h"

namespace itk {

template <class TInputImage, class TFeatureImage, class TOutputType>
VectorThresholdSegmentationLevelSetImageFilter<TInputImage, TFeatureImage, TOutputType>
::VectorThresholdSegmentationLevelSetImageFilter(void)
{
  m_ThresholdFunction = ThresholdFunctionType::New();
  m_ThresholdFunction->SetThreshold(0);

  this->SetSegmentationFunction(m_ThresholdFunction);
}

template <class TInputImage, class TFeatureImage, class TOutputType>
void
VectorThresholdSegmentationLevelSetImageFilter<TInputImage, TFeatureImage, TOutputType>
::PrintSelf(std::ostream &os, Indent indent) const
{
  Superclass::PrintSelf(os, indent);
  os << "ThresholdFunction: " << m_ThresholdFunction;
}

}// end namespace itk

#endif

--------------040506090806040102030906
Content-Type: text/plain;
 name="itkVectorThresholdSegmentationLevelSetFunction.h"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="itkVectorThresholdSegmentationLevelSetFunction.h"

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkVectorThresholdSegmentationLevelSetFunction.h,v $
  Language:  C++
  Date:      $Date: 2003/09/10 14:28:39 $
  Version:   $Revision: 1.12 $

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef __itkVectorThresholdSegmentationLevelSetFunction_h_
#define __itkVectorThresholdSegmentationLevelSetFunction_h_

#include "itkSegmentationLevelSetFunction.h"
#include "itkNumericTraits.h"
#include "itkMahalanobisDistanceMembershipFunction.h"
namespace itk {

/** \class ThresholdSegmentationLevelSetFunction
 *    
 * \brief This function is used in VectorThresholdSegmentationLevelSetImageFilter to
 * segment structures in images based on the Mahalanobis distance.
 *
 * \par  SegmentationLevelSetFunction is a subclass of the generic LevelSetFunction.
 * It useful for segmentations based on intensity values in an image.  It works
 * by constructing a speed term (feature image) with positive values inside an
 * intensity window (between a low and high threshold) and negative values
 * outside that intensity window.  The evolving level set front will lock onto
 * regions that are at the edges of the intensity window.
 *
 *
 *  \par
 *  Image $f$ is thresholded pixel by pixel using threshold $T$
 *  according to the following formula.
 *
 * \par
 *  \f$  f(x) = T - MahalanobisDistance(x)
 *
 * \sa SegmentationLevelSetImageFunction
 *  \sa ThresholdSegmentationLevelSetImageFilter
 *  \sa MahalanobisDistanceMembershipFunction*/
template <class TImageType, class TFeatureImageType>
class ITK_EXPORT VectorThresholdSegmentationLevelSetFunction
  : public SegmentationLevelSetFunction<TImageType, TFeatureImageType>
{
public:
  /** Standard class typedefs. */
  typedef VectorThresholdSegmentationLevelSetFunction Self;
  typedef SegmentationLevelSetFunction<TImageType, TFeatureImageType> Superclass;
  typedef SmartPointer<Self> Pointer;
  typedef SmartPointer<const Self> ConstPointer;
  typedef TFeatureImageType FeatureImageType;

  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** Run-time type information (and related methods) */
  itkTypeMacro( VectorThresholdSegmentationLevelSetFunction, SegmentationLevelSetFunction );

  /** Extract some parameters from the superclass. */
  typedef typename Superclass::ImageType ImageType;
  typedef typename Superclass::ScalarValueType ScalarValueType;
  typedef typename Superclass::FeatureScalarType FeatureScalarType;
  typedef typename Superclass::RadiusType RadiusType;

  /** Extract some parameters from the superclass. */
  itkStaticConstMacro(ImageDimension, unsigned int,
                      Superclass::ImageDimension);

  typedef Statistics::MahalanobisDistanceMembershipFunction<FeatureScalarType> MahalanobisFunctionType;
  typedef typename MahalanobisFunctionType::Pointer MahalanobisFunctionPointer;
  typedef typename MahalanobisFunctionType::MeanVectorType MeanVectorType;
  typedef typename MahalanobisFunctionType::CovarianceMatrixType CovarianceMatrixType;

  /** Set/Get mean and covariance */
  void SetMean(const MeanVectorType &mean) 
  {	m_Mahalanobis->SetMean(mean); }
  const MeanVectorType & GetMean() const 
  {	return m_Mahalanobis->GetMean(); }

  void SetCovariance(const CovarianceMatrixType &cov) 
  { m_Mahalanobis->SetCovariance(cov); }
  const CovarianceMatrixType & GetCovariance() const
  { return m_Mahalanobis->GetCovariance(); }

  /** Set/Get the threshold value for the MahanalobisDistance */
  void SetThreshold(ScalarValueType thr) 
  {
	m_Threshold = thr;

  }
  ScalarValueType GetThreshold() 
  {
	return m_Threshold;
  }

  virtual void CalculateSpeedImage();

  virtual void Initialize(const RadiusType &r)
  {
    Superclass::Initialize(r);

    this->SetAdvectionWeight( NumericTraits<ScalarValueType>::Zero);
    this->SetPropagationWeight(-1.0 * NumericTraits<ScalarValueType>::One);
    this->SetCurvatureWeight(NumericTraits<ScalarValueType>::One);
  }

protected:
  VectorThresholdSegmentationLevelSetFunction()
  {
	  MeanVectorType mean(FeatureImageType::GetImageDimension());
	  CovarianceMatrixType covariance(FeatureImageType::GetImageDimension(),
									 FeatureImageType::GetImageDimension());

	  mean.fill(NumericTraits<FeatureScalarType::ValueType>::Zero);
	  covariance.fill(NumericTraits<FeatureScalarType::ValueType>::Zero);

	  m_Mahalanobis = MahalanobisFunctionType::New();
	  m_Mahalanobis->SetMean(mean);
	  m_Mahalanobis->SetCovariance(covariance);

    this->SetAdvectionWeight(0.0);
    this->SetPropagationWeight(1.0);
    this->SetThreshold(1.8);
  }
  virtual ~VectorThresholdSegmentationLevelSetFunction(){}

  ThresholdSegmentationLevelSetFunction(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemented

  void PrintSelf(std::ostream& os, Indent indent) const
  {
    Superclass::PrintSelf(os, indent );
    os << indent << "MahalanobisFunction: " << m_Mahalanobis << std::endl;
	os << indent << "ThresholdValue: " << m_Threshold << std::endl;
  }

  MahalanobisFunctionPointer m_Mahalanobis;
  ScalarValueType	m_Threshold;

};

} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkVectorThresholdSegmentationLevelSetFunction.txx"
#endif

#endif

--------------040506090806040102030906--