itkMultiphaseFiniteDifferenceImageFilter.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 __itkMultiphaseFiniteDifferenceImageFilter_h
#define __itkMultiphaseFiniteDifferenceImageFilter_h

#include "itkInPlaceImageFilter.h"
#include "itkFiniteDifferenceFunction.h"
#include "vnl/vnl_vector.h"
#include "itkImageRegionIterator.h"

#include "itkListSample.h"
#include "itkKdTreeGenerator.h"

namespace itk
{
template< typename TInputImage,
          typename TFeatureImage,
          typename TOutputImage,
          typename TFiniteDifferenceFunction = FiniteDifferenceFunction< TOutputImage >,
          typename TIdCell = unsigned int >
class MultiphaseFiniteDifferenceImageFilter:
  public InPlaceImageFilter< TFeatureImage, TOutputImage >
{
public:

  typedef MultiphaseFiniteDifferenceImageFilter             Self;
  typedef InPlaceImageFilter< TFeatureImage, TOutputImage > Superclass;
  typedef SmartPointer< Self >                              Pointer;
  typedef SmartPointer< const Self >                        ConstPointer;

  itkTypeMacro(MultiphaseFiniteDifferenceImageFilter, InPlaceImageFilter);

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

  typedef TInputImage                               InputImageType;
  typedef typename InputImageType::Pointer          InputImagePointer;
  typedef typename InputImageType::PointType        InputPointType;
  typedef typename InputPointType::CoordRepType     InputCoordRepType;
  typedef typename InputImageType::IndexType        InputIndexType;
  typedef typename InputIndexType::IndexValueType   InputIndexValueType;
  typedef typename InputImageType::SizeType         InputSizeType;
  typedef typename InputSizeType::SizeValueType     InputSizeValueType;
  typedef typename InputImageType::RegionType       InputRegionType;
  typedef typename InputImageType::PixelType        InputPixelType;
  typedef typename InputImageType::SpacingType      InputSpacingType;
  typedef typename InputImageType::OffsetValueType  InputOffsetValueType;

  typedef TFeatureImage                          FeatureImageType;
  typedef typename FeatureImageType::Pointer     FeatureImagePointer;
  typedef typename FeatureImageType::RegionType  FeatureRegionType;
  typedef typename FeatureImageType::SizeType    FeatureSizeType;
  typedef typename FeatureImageType::SpacingType FeatureSpacingType;
  typedef typename FeatureImageType::PointType   FeaturePointType;
  typedef typename FeatureImageType::PixelType   FeaturePixelType;

  typedef TOutputImage                             OutputImageType;
  typedef typename OutputImageType::Pointer        OutputImagePointer;
  typedef typename OutputImageType::PixelType      OutputPixelType;
  typedef typename OutputImageType::RegionType     OutputRegionType;
  typedef typename OutputImageType::SizeType       OutputSizeType;
  typedef typename OutputImageType::SizeValueType  OutputSizeValueType;
  typedef typename OutputImageType::IndexType      OutputIndexType;
  typedef typename OutputImageType::IndexValueType OutputIndexValueType;

  typedef TIdCell                   IdCellType;
  typedef std::vector< IdCellType > VectorIdCellType;

  typedef TFiniteDifferenceFunction                           FiniteDifferenceFunctionType;
  typedef typename FiniteDifferenceFunctionType::Pointer      FiniteDifferenceFunctionPointer;
  typedef typename FiniteDifferenceFunctionType::TimeStepType TimeStepType;
  typedef typename std::vector< TimeStepType >                TimeStepVectorType;
  typedef typename FiniteDifferenceFunctionType::RadiusType   RadiusType;

  typedef Vector< float, itkGetStaticConstMacro(ImageDimension) >
  CentroidVectorType;
  typedef Statistics::ListSample< CentroidVectorType > SampleType;
  typedef Statistics::KdTreeGenerator< SampleType >    KdTreeGeneratorType;
  typedef typename KdTreeGeneratorType::Pointer        KdTreeGeneratorPointer;
  typedef typename KdTreeGeneratorType::KdTreeType     KdTreeType;
  typedef typename KdTreeType::Pointer                 KdTreePointer;

  virtual const FiniteDifferenceFunctionPointer GetDifferenceFunction(
    const IdCellType & functionIndex) const
  {
    if ( functionIndex < m_FunctionCount )
      {
      return ( this->m_DifferenceFunctions[functionIndex] );
      }
    else
      {
      return ITK_NULLPTR;
      }
  }

  virtual void SetDifferenceFunction(const IdCellType & functionIndex,
                                     FiniteDifferenceFunctionPointer function)
  {
    if ( functionIndex < m_FunctionCount )
      {
      this->m_DifferenceFunctions[functionIndex] = function;
      }
  }

  itkSetMacro(NumberOfIterations, unsigned int);
  itkGetConstReferenceMacro(NumberOfIterations, unsigned int);

  itkSetMacro(UseImageSpacing, bool);
  itkBooleanMacro(UseImageSpacing);
  itkGetConstReferenceMacro(UseImageSpacing, bool);

  itkSetMacro(MaximumRMSError, double);
  itkGetConstReferenceMacro(MaximumRMSError, double);

  itkSetMacro(RMSChange, double);
  itkGetConstReferenceMacro(RMSChange, double);

  itkSetMacro(InitializedState, bool);
  itkGetConstReferenceMacro(InitializedState, bool);
  itkBooleanMacro(InitializedState);

  itkSetMacro(ManualReinitialization, bool);
  itkGetConstReferenceMacro(ManualReinitialization, bool);
  itkBooleanMacro(ManualReinitialization);

  itkSetMacro(ElapsedIterations, unsigned int);

  itkGetConstReferenceMacro(ElapsedIterations, unsigned int);

  void SetLevelSet(const IdCellType & i, const InputImageType *levelSet)
  {
    m_LevelSet[i] = InputImageType::New();
    m_LevelSet[i]->SetRequestedRegion( levelSet->GetRequestedRegion() );
    m_LevelSet[i]->SetBufferedRegion( levelSet->GetBufferedRegion() );
    m_LevelSet[i]->SetLargestPossibleRegion( levelSet->GetLargestPossibleRegion() );
    m_LevelSet[i]->Allocate();
    m_LevelSet[i]->CopyInformation(levelSet);

    ImageRegionConstIterator< InputImageType > in ( levelSet, levelSet->GetBufferedRegion() );
    ImageRegionIterator< InputImageType >      cp ( m_LevelSet[i], levelSet->GetBufferedRegion() );

    in.GoToBegin();
    cp.GoToBegin();

    while ( !in.IsAtEnd() )
      {
      cp.Set( in.Get() );
      ++in;
      ++cp;
      }
  }

  InputImagePointer GetLevelSet(const IdCellType & i)
  {
    if ( i >= m_FunctionCount )
      {
      itkExceptionMacro("Request for level set #" << i
                                                  << " but there are only " << m_FunctionCount);
      }
    else
      {
      return m_LevelSet[i];
      }
  }

  void SetLookup(VectorIdCellType lookup)
  {
    this->m_Lookup = lookup;
  }

  void SetKdTree(KdTreeType *kdtree)
  {
    this->m_KdTree = kdtree;
  }

  void SetFunctionCount(const IdCellType & n)
  {
    m_FunctionCount = n;

    m_DifferenceFunctions.resize(m_FunctionCount, ITK_NULLPTR);

    RadiusType radius;
    radius.Fill(1);

    for ( unsigned int i = 0; i < this->m_FunctionCount; i++ )
      {
      this->m_DifferenceFunctions[i] = FiniteDifferenceFunctionType::New();
      this->m_DifferenceFunctions[i]->Initialize(radius);
      }

    // Initialize the images
    m_LevelSet.resize(m_FunctionCount, ITK_NULLPTR);

    // Initialize the lookup table
    this->m_Lookup.resize(m_FunctionCount);

    IdCellType k = 1;

    typedef typename std::vector< IdCellType >::iterator VectorIteratorType;

    VectorIteratorType it = this->m_Lookup.begin();

    while ( it != this->m_Lookup.end() )
      {
      *it = k;
      ++it;
      ++k;
      }
  }

protected:
  MultiphaseFiniteDifferenceImageFilter()
  {
    this->m_KdTree = ITK_NULLPTR;
    this->m_ElapsedIterations = 0;
    this->m_MaximumRMSError = vnl_math::eps;
    this->m_RMSChange = NumericTraits< double >::max();
    this->m_UseImageSpacing = true;
    this->m_ManualReinitialization = false;
    this->m_InitializedState = false;
    this->m_NumberOfIterations = NumericTraits< unsigned int >::max();
    this->m_FunctionCount = 0;
    this->InPlaceOff();
  }

  ~MultiphaseFiniteDifferenceImageFilter(){}

  IdCellType                       m_FunctionCount;
  std::vector< InputImagePointer > m_LevelSet;
  VectorIdCellType                 m_Lookup;
  KdTreePointer                    m_KdTree;

  unsigned int m_ElapsedIterations;
  double       m_MaximumRMSError;
  double       m_RMSChange;
  unsigned int m_NumberOfIterations;

  std::vector< FiniteDifferenceFunctionPointer > m_DifferenceFunctions;

  bool m_UseImageSpacing;

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

  virtual void AllocateUpdateBuffer() = 0;

  virtual void ApplyUpdate(TimeStepType dt) = 0;

  virtual TimeStepType CalculateChange() = 0;

  virtual void CopyInputToOutput() = 0;

  virtual void GenerateData();

  virtual void GenerateInputRequestedRegion();

  virtual bool Halt();

  virtual bool ThreadedHalt( void *itkNotUsed(threadInfo) )
  {
    return this->Halt();
  }

  virtual void Initialize() {}

  virtual void InitializeIteration()
  {
    for ( IdCellType i = 0; i < this->m_FunctionCount; i++ )
      {
      this->m_DifferenceFunctions[i]->InitializeIteration();
      }
  }

  inline TimeStepType ResolveTimeStep(const TimeStepVectorType & timeStepList,
                                      const std::vector< bool > & valid);

  virtual void PostProcessOutput() {}

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

  bool m_ManualReinitialization;

  bool m_InitializedState;
};
} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkMultiphaseFiniteDifferenceImageFilter.hxx"
#endif

#endif