itkRegionBasedLevelSetFunction.h

Go to the documentation of this file.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkRegionBasedLevelSetFunction.h,v $
  Language:  C++
  Date:      $Date: 2010-02-24 22:42:16 $
  Version:   $Revision: 1.19 $

  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 __itkRegionBasedLevelSetFunction_h
#define __itkRegionBasedLevelSetFunction_h

#include "itkFiniteDifferenceFunction.h"
#include "itkRegularizedHeavisideStepFunction.h"
#include "vnl/vnl_matrix_fixed.h"

namespace itk {

template < class TInput, // LevelSetImageType
  class TFeature, // FeatureImageType
  class TSharedData >
class ITK_EXPORT RegionBasedLevelSetFunction: public
FiniteDifferenceFunction< TInput >
{
public:
  typedef RegionBasedLevelSetFunction                Self;
  typedef FiniteDifferenceFunction< TInput >         Superclass;
  typedef SmartPointer<Self>                         Pointer;
  typedef SmartPointer<const Self>                   ConstPointer;

  itkStaticConstMacro(ImageDimension, unsigned int, Superclass::ImageDimension);

  // itkNewMacro() is purposely not provided since this is an abstract class.

  itkTypeMacro( RegionBasedLevelSetFunction, FiniteDifferenceFunction );

  typedef double                                          TimeStepType;
  typedef typename Superclass::ImageType                  ImageType;
  typedef typename Superclass::PixelType                  PixelType;
  typedef PixelType                                       ScalarValueType;
  typedef typename Superclass::RadiusType                 RadiusType;
  typedef typename Superclass::NeighborhoodType           NeighborhoodType;
  typedef typename Superclass::NeighborhoodScalesType     NeighborhoodScalesType;
  typedef typename Superclass::FloatOffsetType            FloatOffsetType;
  typedef FixedArray< ScalarValueType, itkGetStaticConstMacro(ImageDimension) >
                                                          VectorType;

  /* This structure is derived from LevelSetFunction and stores intermediate
  values for computing time step sizes */
  struct GlobalDataStruct
    {
    GlobalDataStruct()
      {
      ScalarValueType null_value = NumericTraits<ScalarValueType>::Zero;

      m_MaxCurvatureChange   = null_value;
      m_MaxAdvectionChange   = null_value;
      m_MaxGlobalChange      = null_value;
      }

    ~GlobalDataStruct() {}

    vnl_matrix_fixed<ScalarValueType,
      itkGetStaticConstMacro(ImageDimension),
      itkGetStaticConstMacro(ImageDimension)> m_dxy;

    ScalarValueType m_dx[itkGetStaticConstMacro(ImageDimension)];

    ScalarValueType m_dx_forward[ itkGetStaticConstMacro( ImageDimension ) ];
    ScalarValueType m_dx_backward[ itkGetStaticConstMacro( ImageDimension ) ];

    ScalarValueType m_GradMagSqr;
    ScalarValueType m_GradMag;

    ScalarValueType m_MaxCurvatureChange;
    ScalarValueType m_MaxAdvectionChange;
    ScalarValueType m_MaxGlobalChange;
    };


  typedef TInput                                    InputImageType;
  typedef typename InputImageType::ConstPointer     InputImageConstPointer;
  typedef typename InputImageType::Pointer          InputImagePointer;
  typedef typename InputImageType::PixelType        InputPixelType;
  typedef typename InputImageType::IndexType        InputIndexType;
  typedef typename InputImageType::IndexValueType   InputIndexValueType;
  typedef typename InputImageType::SizeType         InputSizeType;
  typedef typename InputImageType::SizeValueType    InputSizeValueType;
  typedef typename InputImageType::RegionType       InputRegionType;
  typedef typename InputImageType::PointType        InputPointType;

  typedef TFeature                                  FeatureImageType;
  typedef typename FeatureImageType::ConstPointer   FeatureImageConstPointer;
  typedef typename FeatureImageType::PixelType      FeaturePixelType;
  typedef typename FeatureImageType::IndexType      FeatureIndexType;
  typedef typename FeatureImageType::SpacingType    FeatureSpacingType;
  typedef typename FeatureImageType::OffsetType     FeatureOffsetType;

  typedef TSharedData                               SharedDataType;
  typedef typename SharedDataType::Pointer          SharedDataPointer;

  typedef HeavisideStepFunctionBase< InputPixelType, InputPixelType > HeavisideFunctionType;
  typedef typename HeavisideFunctionType::ConstPointer HeavisideFunctionConstPointer;

  void SetDomainFunction( const HeavisideFunctionType * f )
    {
    this->m_DomainFunction = f;
    }

  virtual void Initialize(const RadiusType &r)
    {
    this->SetRadius(r);

    // Dummy neighborhood.
    NeighborhoodType it;
    it.SetRadius( r );

    // Find the center index of the neighborhood.
    m_Center =  it.Size() / 2;

    // Get the stride length for each axis.
    for(unsigned int i = 0; i < ImageDimension; i++)
      {
      m_xStride[i] = it.GetStride(i);
      }
    }


  void SetSharedData( SharedDataPointer sharedDataIn )
    {
    this->m_SharedData = sharedDataIn;
    }

  void UpdateSharedData( bool forceUpdate );

  void *GetGlobalDataPointer() const
    {
    return new GlobalDataStruct;
    }

  TimeStepType ComputeGlobalTimeStep(void *GlobalData) const;

  virtual PixelType ComputeUpdate(const NeighborhoodType &neighborhood,
    void *globalData, const FloatOffsetType& = FloatOffsetType(0.0));

  void SetInitialImage(InputImageType *f)
    {
    m_InitialImage = f;
    }

  virtual const FeatureImageType *GetFeatureImage() const
    { return m_FeatureImage.GetPointer(); }
  virtual void SetFeatureImage(const FeatureImageType *f)
    {
    m_FeatureImage = f;

    FeatureSpacingType spacing = m_FeatureImage->GetSpacing();
    for(unsigned int i = 0; i < ImageDimension; i++)
      {
      this->m_InvSpacing[i] = 1/spacing[i];
      }
    }

  virtual VectorType AdvectionField(const NeighborhoodType &,
                                    const FloatOffsetType &, GlobalDataStruct * = 0)  const
    { return this->m_ZeroVectorConstant; }

  void SetAreaWeight( const ScalarValueType& nu)
    { this->m_AreaWeight = nu; }
  ScalarValueType GetAreaWeight() const
    { return this->m_AreaWeight; }

  void SetLambda1( const ScalarValueType& lambda1 )
    { this->m_Lambda1 = lambda1; }
  ScalarValueType GetLambda1() const
    { return this->m_Lambda1; }

  void SetLambda2( const ScalarValueType& lambda2 )
    { this->m_Lambda2 = lambda2; }
  ScalarValueType GetLambda2() const
    { return this->m_Lambda2; }

  void SetOverlapPenaltyWeight( const ScalarValueType& gamma )
    { this->m_OverlapPenaltyWeight = gamma; }
  ScalarValueType GetOverlapPenaltyWeight() const
    { return this->m_OverlapPenaltyWeight; }

  virtual void SetCurvatureWeight(const ScalarValueType c)
    { m_CurvatureWeight = c; }
  ScalarValueType GetCurvatureWeight() const
    { return m_CurvatureWeight; }

  void SetAdvectionWeight( const ScalarValueType& iA)
    { this->m_AdvectionWeight = iA; }
  ScalarValueType GetAdvectionWeight() const
    { return this->m_AdvectionWeight; }

  void SetReinitializationSmoothingWeight(const ScalarValueType c)
    { m_ReinitializationSmoothingWeight = c; }
  ScalarValueType GetReinitializationSmoothingWeight() const
    { return m_ReinitializationSmoothingWeight; }

  void SetVolumeMatchingWeight( const ScalarValueType& tau )
    { this->m_VolumeMatchingWeight = tau; }
  ScalarValueType GetVolumeMatchingWeight() const
    { return this->m_VolumeMatchingWeight; }

  void SetVolume( const ScalarValueType& volume )
    { this->m_Volume = volume; }
  ScalarValueType GetVolume() const
    { return this->m_Volume; }

  void SetFunctionId( const unsigned int& iFid )
    { this->m_FunctionId = iFid; }

  virtual void ReleaseGlobalDataPointer(void *GlobalData) const
  { delete (GlobalDataStruct *) GlobalData; }

  virtual ScalarValueType ComputeCurvature(const NeighborhoodType &,
    const FloatOffsetType &, GlobalDataStruct *gd );

  virtual ScalarValueType LaplacianSmoothingSpeed(
    const NeighborhoodType &,
    const FloatOffsetType &, GlobalDataStruct * = 0) const
    { return NumericTraits<ScalarValueType>::One; }

  virtual ScalarValueType CurvatureSpeed(const NeighborhoodType &,
                                         const FloatOffsetType &, GlobalDataStruct * = 0
                                         ) const
    { return NumericTraits<ScalarValueType>::One; }

  virtual void CalculateAdvectionImage() {}

protected:

  RegionBasedLevelSetFunction();
  virtual ~RegionBasedLevelSetFunction() {}

  InputImageConstPointer m_InitialImage;

  FeatureImageConstPointer m_FeatureImage;

  SharedDataPointer                 m_SharedData;
  HeavisideFunctionConstPointer     m_DomainFunction;

  ScalarValueType           m_AreaWeight;

  ScalarValueType           m_Lambda1;

  ScalarValueType           m_Lambda2;

  ScalarValueType           m_OverlapPenaltyWeight;

  ScalarValueType           m_VolumeMatchingWeight;

  ScalarValueType           m_Volume;

  ScalarValueType           m_CurvatureWeight;

  ScalarValueType           m_AdvectionWeight;

  ScalarValueType           m_ReinitializationSmoothingWeight;

  unsigned int              m_FunctionId;

  std::slice x_slice[itkGetStaticConstMacro(ImageDimension)];
  ::size_t m_Center;
  ::size_t m_xStride[itkGetStaticConstMacro(ImageDimension)];
  double m_InvSpacing[itkGetStaticConstMacro(ImageDimension)];

  static double m_WaveDT;
  static double m_DT;

  void ComputeHImage();

  ScalarValueType ComputeGlobalTerm(
    const ScalarValueType& imagePixel,
    const InputIndexType& inputIndex );

  virtual ScalarValueType ComputeInternalTerm(const FeaturePixelType& iValue,
    const FeatureIndexType& iIdx ) = 0;

  virtual ScalarValueType ComputeExternalTerm(const FeaturePixelType& iValue,
    const FeatureIndexType& iIdx ) = 0;

  virtual ScalarValueType ComputeOverlapParameters( const FeatureIndexType& featIndex,
    ScalarValueType& pr ) = 0;

  ScalarValueType ComputeVolumeRegularizationTerm( );

  ScalarValueType ComputeLaplacian( GlobalDataStruct *gd );

  void ComputeHessian( const NeighborhoodType &it,
    GlobalDataStruct *globalData );

  virtual void ComputeParameters() = 0;

  virtual void UpdateSharedDataParameters() = 0;

  bool m_UpdateC;

  static VectorType InitializeZeroVectorConstant();

  static VectorType m_ZeroVectorConstant;

private:
  RegionBasedLevelSetFunction(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemented
};

} // end namespace itk

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

#endif