itkRegionBasedLevelSetFunction.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef __itkRegionBasedLevelSetFunction_h
#define __itkRegionBasedLevelSetFunction_h

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

namespace itk
{
template< typename TInput,   // LevelSetImageType
          typename TFeature, // FeatureImageType
          typename TSharedData >
class 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);
      }
  }

#if !defined( CABLE_CONFIGURATION )
  void SetSharedData(SharedDataPointer sharedDataIn)
  {
    this->m_SharedData = sharedDataIn;
  }
#endif

  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)];
  OffsetValueType m_Center;
  OffsetValueType 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.hxx"
#endif

#endif