itkMultiphaseFiniteDifferenceImageFilter.h

Go to the documentation of this file.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkMultiphaseFiniteDifferenceImageFilter.h,v $
  Language:  C++
  Date:      $Date: 2009-07-29 15:14:03 $
  Version:   $Revision: 1.7 $

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

#include "itkArray.h"
#include "itkInPlaceImageFilter.h"
#include "itkFiniteDifferenceFunction.h"
#include <vnl/vnl_vector.h>
#include "itkImageRegionIterator.h"

#include "itkVector.h"
#include "itkListSample.h"
#include "itkKdTree.h"
#include "itkKdTreeGenerator.h"

namespace itk {

 template < class TInputImage,
  class TFeatureImage,
  class TOutputImage,
  class TFiniteDifferenceFunction = FiniteDifferenceFunction<TOutputImage>,
  typename TIdCell = unsigned int >
class ITK_EXPORT 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 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 0;
      }
    }

  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, 0 );

    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, 0 );

    // 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 = 0;
    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.txx"
#endif

#endif