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 ITK_TEMPLATE_EXPORT MultiphaseFiniteDifferenceImageFilter:
  public InPlaceImageFilter< TFeatureImage, TOutputImage >
{
public:
  ITK_DISALLOW_COPY_AND_ASSIGN(MultiphaseFiniteDifferenceImageFilter);

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

  itkTypeMacro(MultiphaseFiniteDifferenceImageFilter, InPlaceImageFilter);

  static constexpr unsigned int ImageDimension = TOutputImage::ImageDimension;

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

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

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

  using IdCellType = TIdCell;
  using VectorIdCellType = std::vector< IdCellType >;

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

  using CentroidVectorType =
      Vector< float, Self::ImageDimension >;
  using SampleType = Statistics::ListSample< CentroidVectorType >;
  using KdTreeGeneratorType = Statistics::KdTreeGenerator< SampleType >;
  using KdTreeGeneratorPointer = typename KdTreeGeneratorType::Pointer;
  using KdTreeType = typename KdTreeGeneratorType::KdTreeType;
  using KdTreePointer = typename KdTreeType::Pointer;

  virtual const FiniteDifferenceFunctionPointer GetDifferenceFunction(
    const IdCellType & functionIndex) const
  {
    if ( functionIndex < m_FunctionCount )
      {
      return ( this->m_DifferenceFunctions[functionIndex] );
      }
    else
      {
      return 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, 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, nullptr);

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

    IdCellType k = 1;

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

protected:
  MultiphaseFiniteDifferenceImageFilter()
  {
    this->m_KdTree = nullptr;
    this->m_ElapsedIterations = 0;
    this->m_MaximumRMSError = itk::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() override{}

  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 override;

  virtual void AllocateUpdateBuffer() = 0;

  virtual void ApplyUpdate(TimeStepType dt) = 0;

  virtual TimeStepType CalculateChange() = 0;

  virtual void CopyInputToOutput() = 0;

  void GenerateData() override;

  void GenerateInputRequestedRegion() override;

  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:
  bool m_ManualReinitialization;

  bool m_InitializedState;
};
} // end namespace itk

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

#endif