itkFastMarchingImageFilter.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 itkFastMarchingImageFilter_h
#define itkFastMarchingImageFilter_h

#include "itkImageToImageFilter.h"
#include "itkImageRegionConstIteratorWithIndex.h"
#include "itkLevelSet.h"
#include "itkMath.h"

#include <functional>
#include <queue>
#include "itkMath.h"

namespace itk
{
template<
  typename TLevelSet,
  typename TSpeedImage = Image< float,  TLevelSet ::ImageDimension > >
class ITK_TEMPLATE_EXPORT FastMarchingImageFilter:
  public ImageToImageFilter< TSpeedImage, TLevelSet >
{
public:
  ITK_DISALLOW_COPY_AND_ASSIGN(FastMarchingImageFilter);

  using Self = FastMarchingImageFilter;
  using Superclass = ImageToImageFilter< TSpeedImage, TLevelSet >;
  using Pointer = SmartPointer< Self >;
  using ConstPointer = SmartPointer< const Self >;

  itkNewMacro(Self);

  itkTypeMacro(FastMarchingImageFilter, ImageToImageFilter);

  using LevelSetType = LevelSetTypeDefault< TLevelSet >;
  using LevelSetImageType = typename LevelSetType::LevelSetImageType;
  using LevelSetPointer = typename LevelSetType::LevelSetPointer;
  using PixelType = typename LevelSetType::PixelType;
  using NodeType = typename LevelSetType::NodeType;
  using NodeIndexType = typename NodeType::IndexType;
  using NodeContainer = typename LevelSetType::NodeContainer;
  using NodeContainerPointer = typename LevelSetType::NodeContainerPointer;
  using OutputSizeType = typename LevelSetImageType::SizeType;
  using OutputRegionType = typename LevelSetImageType::RegionType;
  using OutputSpacingType = typename LevelSetImageType::SpacingType;
  using OutputDirectionType = typename LevelSetImageType::DirectionType;
  using OutputPointType = typename LevelSetImageType::PointType;

  class AxisNodeType:public NodeType
  {
public:
    AxisNodeType()  {}
    int GetAxis() const { return m_Axis; }
    void SetAxis(int axis) { m_Axis = axis; }
    const AxisNodeType & operator=(const NodeType & node)
    { this->NodeType::operator=(node); return *this; }

private:
    int m_Axis{0};
  };

  using SpeedImageType = TSpeedImage;

  using SpeedImagePointer = typename SpeedImageType::Pointer;
  using SpeedImageConstPointer = typename SpeedImageType::ConstPointer;

  static constexpr unsigned int SetDimension = LevelSetType::SetDimension;
  static constexpr unsigned int SpeedImageDimension = SpeedImageType::ImageDimension;

  using IndexType = Index< Self::SetDimension >;

  enum LabelType { FarPoint = 0, AlivePoint,
                   TrialPoint, InitialTrialPoint, OutsidePoint };

  using LabelImageType = Image< unsigned char, Self::SetDimension >;

  using LabelImagePointer = typename LabelImageType::Pointer;

  template< typename TPixel >
  void SetBinaryMask( Image< TPixel, SetDimension >* iImage )
    {
    using InternalImageType = Image< TPixel, SetDimension >;
    using InternalRegionIterator = ImageRegionConstIteratorWithIndex<InternalImageType>;
    InternalRegionIterator b_it( iImage, iImage->GetLargestPossibleRegion() );
    b_it.GoToBegin();

    TPixel zero_value = NumericTraits< TPixel >::ZeroValue();
    typename NodeContainer::ElementIdentifier NumberOfPoints = 0;

    NodeType node;
    node.SetValue( 0. );

    while( !b_it.IsAtEnd() )
      {
      if( Math::ExactlyEquals(b_it.Get(), zero_value) )
        {
        if( NumberOfPoints == 0 )
          {
          m_OutsidePoints = NodeContainer::New();
          }
        node.SetIndex( b_it.GetIndex() );
        m_OutsidePoints->InsertElement( NumberOfPoints++, node );

        }
      ++b_it;
      }
    this->Modified();
    }

  void SetOutsidePoints(NodeContainer *points)
  {
    m_OutsidePoints = points;
    this->Modified();
  }

  void SetAlivePoints(NodeContainer *points)
  {
    m_AlivePoints = points;
    this->Modified();
  }

  NodeContainerPointer GetAlivePoints()
  {
    return m_AlivePoints;
  }

  void SetTrialPoints(NodeContainer *points)
  {
    m_TrialPoints = points;
    this->Modified();
  }

  NodeContainerPointer GetTrialPoints()
  {
    return m_TrialPoints;
  }

  LabelImagePointer GetLabelImage() const
  {
    return m_LabelImage;
  }

  void SetSpeedConstant(double value)
  {
    m_SpeedConstant = value;
    m_InverseSpeed = -1.0 * itk::Math::sqr(1.0 / m_SpeedConstant);
    this->Modified();
  }

  itkGetConstReferenceMacro(SpeedConstant, double);

  itkSetMacro(NormalizationFactor, double);
  itkGetConstMacro(NormalizationFactor, double);

  itkSetMacro(StoppingValue, double);

  itkGetConstReferenceMacro(StoppingValue, double);

  itkSetMacro(CollectPoints, bool);

  itkGetConstReferenceMacro(CollectPoints, bool);
  itkBooleanMacro(CollectPoints);

  NodeContainerPointer GetProcessedPoints() const
  {
    return m_ProcessedPoints;
  }

  virtual void SetOutputSize(const OutputSizeType & size)
  { m_OutputRegion = size; }
  virtual OutputSizeType GetOutputSize() const
  { return m_OutputRegion.GetSize(); }
  itkSetMacro(OutputRegion, OutputRegionType);
  itkGetConstReferenceMacro(OutputRegion, OutputRegionType);
  itkSetMacro(OutputSpacing, OutputSpacingType);
  itkGetConstReferenceMacro(OutputSpacing, OutputSpacingType);
  itkSetMacro(OutputDirection, OutputDirectionType);
  itkGetConstReferenceMacro(OutputDirection, OutputDirectionType);
  itkSetMacro(OutputOrigin, OutputPointType);
  itkGetConstReferenceMacro(OutputOrigin, OutputPointType);
  itkSetMacro(OverrideOutputInformation, bool);
  itkGetConstReferenceMacro(OverrideOutputInformation, bool);
  itkBooleanMacro(OverrideOutputInformation);

#ifdef ITK_USE_CONCEPT_CHECKING
  // Begin concept checking
  itkConceptMacro( SameDimensionCheck,
                   ( Concept::SameDimension< SetDimension, SpeedImageDimension > ) );
  itkConceptMacro( SpeedConvertibleToDoubleCheck,
                   ( Concept::Convertible< typename TSpeedImage::PixelType, double > ) );
  itkConceptMacro( DoubleConvertibleToLevelSetCheck,
                   ( Concept::Convertible< double, PixelType > ) );
  itkConceptMacro( LevelSetOStreamWritableCheck,
                   ( Concept::OStreamWritable< PixelType > ) );
  // End concept checking
#endif

protected:
  FastMarchingImageFilter();
  ~FastMarchingImageFilter() override = default;
  void PrintSelf(std::ostream & os, Indent indent) const override;

  virtual void Initialize(LevelSetImageType *);

  virtual void UpdateNeighbors(const IndexType & index,
                               const SpeedImageType *, LevelSetImageType *);

  virtual double UpdateValue(const IndexType & index,
                             const SpeedImageType *, LevelSetImageType *);

  const AxisNodeType & GetNodeUsedInCalculation(unsigned int idx) const
  { return m_NodesUsed[idx]; }

  void GenerateData() override;

  void GenerateOutputInformation() override;

  void EnlargeOutputRequestedRegion(DataObject *output) override;

  itkGetConstReferenceMacro(LargeValue, PixelType);

  OutputRegionType m_BufferedRegion;
  using LevelSetIndexType = typename LevelSetImageType::IndexType;
  LevelSetIndexType m_StartIndex;
  LevelSetIndexType m_LastIndex;

  itkGetConstReferenceMacro(StartIndex, LevelSetIndexType);
  itkGetConstReferenceMacro(LastIndex, LevelSetIndexType);

private:
  NodeContainerPointer m_AlivePoints;
  NodeContainerPointer m_TrialPoints;
  NodeContainerPointer m_OutsidePoints;

  LabelImagePointer m_LabelImage;

  double m_SpeedConstant;
  double m_InverseSpeed;
  double m_StoppingValue;

  bool                 m_CollectPoints;
  NodeContainerPointer m_ProcessedPoints;

  OutputRegionType    m_OutputRegion;
  OutputPointType     m_OutputOrigin;
  OutputSpacingType   m_OutputSpacing;
  OutputDirectionType m_OutputDirection;
  bool                m_OverrideOutputInformation;

  typename LevelSetImageType::PixelType m_LargeValue;
  AxisNodeType m_NodesUsed[SetDimension];

  using HeapContainer = std::vector< AxisNodeType >;
  using NodeComparer = std::greater< AxisNodeType >;
  using HeapType =
      std::priority_queue< AxisNodeType, HeapContainer, NodeComparer >;

  HeapType m_TrialHeap;

  double m_NormalizationFactor;
};
} // namespace itk

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

#endif