Consolidate interpolator

examples/CMakeLists.txt

@@ -16,3 +16,6 @@ include(${ITK_USE_FILE})
 
 add_executable(MinimalPathExamples example.cxx)
 target_link_libraries(MinimalPathExamples ${ITK_LIBRARIES})
+
+add_executable(interptest interpolator.cxx)
+target_link_libraries(interptest ${ITK_LIBRARIES})

include/itkArrivalFunctionToPathFilter.h

@@ -72,7 +72,8 @@ class ITK_TEMPLATE_EXPORT ArrivalFunctionToPathCommand : public itk::Command
   void
   Execute(const itk::Object * object, const itk::EventObject & event) override
   {
-    if (!itk::IterationEvent().CheckEvent(&event))
+    if (!itk::IterationEvent().CheckEvent(&event) &&
+	!itk::EndEvent().CheckEvent(&event))
     {
       return;
     }
@@ -249,6 +250,11 @@ class ITK_TEMPLATE_EXPORT ArrivalFunctionToPathFilter : public ImageToPathFilter
   typename OptimizerType::MeasureType m_TerminationValue;
   std::vector<PointsContainerType>    m_PointList;
   unsigned int                        m_CurrentOutput;
+
+  // variable for communicating state of a multisegment path
+  // Sometimes an optimizer may terminate during one segment of the
+  // path, and needs to be restarted
+  bool m_TrackingPath;
 };
 
 } // end namespace itk

include/itkArrivalFunctionToPathFilter.hxx

@@ -158,7 +158,8 @@ ArrivalFunctionToPathFilter<TInputImage, TOutputPath>::GenerateData()
   typename CommandType::Pointer observer = CommandType::New();
   observer->SetFilter(this);
 
-  unsigned long observerTag = m_Optimizer->AddObserver(itk::IterationEvent(), observer);
+  unsigned long observerTagIter = m_Optimizer->AddObserver(itk::IterationEvent(), observer);
+  unsigned long observerTagEnd  = m_Optimizer->AddObserver(itk::EndEvent(), observer);
 
   // Do for each output
   for (unsigned int n = 0; n < numberOfOutputs; n++)
@@ -194,11 +195,15 @@ ArrivalFunctionToPathFilter<TInputImage, TOutputPath>::GenerateData()
     m_Optimizer->SetInitialPosition(end);
 
     // Use optimizer to back propagate from end point
-    m_Optimizer->StartOptimization();
+    m_TrackingPath = false;
+    do {
+      m_Optimizer->StartOptimization();
+    } while(m_TrackingPath);
   }
 
   // Clean up by removing observer
-  m_Optimizer->RemoveObserver(observerTag);
+  m_Optimizer->RemoveObserver(observerTagIter);
+  m_Optimizer->RemoveObserver(observerTagEnd);
 }
 
 template <typename TInputImage, typename TOutputPath>

include/itkLinearInterpolateSelectedNeighborsImageFunction.h

@@ -0,0 +1,182 @@
+/*=========================================================================
+ *
+ *  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 itkLinearInterpolateSelectedNeighborsImageFunction_h
+#define itkLinearInterpolateSelectedNeighborsImageFunction_h
+
+#include "itkInterpolateImageFunction.h"
+#include "itkVariableLengthVector.h"
+
+namespace itk
+{
+/** \class LinearInterpolateSelectedNeighborsImageFunction
+ * \brief Linearly interpolate an image at specified positions.
+ *
+ * LinearInterpolateSelectedNeighborsImageFunction linearly interpolates image intensity at
+ * a non-integer pixel position. This class is templated
+ * over the input image type and the coordinate representation type
+ * (e.g. float or double).
+ *
+ * This function works for N-dimensional images.
+ *
+ * This function uses a functor to verify whether neighbours are valid.
+ * The intention is to enable reliable estimations of pixel values in
+ * the presence of indicator pixels, as generated by fast marching methods.
+ *
+ * It is not optimised for specific dimensions.
+ *
+ * \ingroup ImageFunctions ImageInterpolators
+ * \ingroup ITKImageFunction
+ *
+ */
+
+  namespace Functor
+{
+
+template< typename TInput >
+class AllNeighbors
+{
+public:
+  AllNeighbors() = default;
+  ~AllNeighbors() = default;
+  bool operator!=(const AllNeighbors &) const
+  {
+    return false;
+  }
+
+  bool operator==(const AllNeighbors & other) const
+  {
+    return !( *this != other );
+  }
+
+  inline bool operator()(const TInput & A) const
+  { return true; }
+};
+}
+
+
+  template< typename TInputImage, typename TCoordRep = double, typename TNeighborCheckFunction = Functor::AllNeighbors<typename TInputImage::PixelType> >
+class ITK_TEMPLATE_EXPORT LinearInterpolateSelectedNeighborsImageFunction:
+  public InterpolateImageFunction< TInputImage, TCoordRep >
+{
+public:
+  ITK_DISALLOW_COPY_AND_ASSIGN(LinearInterpolateSelectedNeighborsImageFunction);
+
+  /** Standard class type aliases. */
+  using Self = LinearInterpolateSelectedNeighborsImageFunction;
+  using Superclass = InterpolateImageFunction< TInputImage, TCoordRep >;
+  using Pointer = SmartPointer< Self >;
+  using ConstPointer = SmartPointer< const Self >;
+
+  /** Run-time type information (and related methods). */
+  itkTypeMacro(LinearInterpolateSelectedNeighborsImageFunction, InterpolateImageFunction);
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  /** OutputType type alias support */
+  using OutputType = typename Superclass::OutputType;
+
+  /** InputImageType type alias support */
+  using InputImageType = typename Superclass::InputImageType;
+
+  /** InputPixelType type alias support */
+  using InputPixelType = typename Superclass::InputPixelType;
+
+  /** RealType type alias support */
+  using RealType = typename Superclass::RealType;
+
+  /** Size type alias support */
+  using SizeType = typename Superclass::SizeType;
+
+  /** Dimension underlying input image. */
+  static constexpr unsigned int ImageDimension = Superclass::ImageDimension;
+
+  /** Index type alias support */
+  using IndexType = typename Superclass::IndexType;
+
+  /** ContinuousIndex type alias support */
+  using ContinuousIndexType = typename Superclass::ContinuousIndexType;
+  using InternalComputationType = typename ContinuousIndexType::ValueType;
+
+  using NeighborCheckFunctionType = TNeighborCheckFunction;
+
+  /** Evaluate the function at a ContinuousIndex position
+   *
+   * Returns the linearly interpolated image intensity at a
+   * specified point position. No bounds checking is done.
+   * The point is assume to lie within the image buffer.
+   *
+   * ImageFunction::IsInsideBuffer() can be used to check bounds before
+   * calling the method. */
+  OutputType EvaluateAtContinuousIndex(const
+                                                ContinuousIndexType &
+                                                index) const override
+  {
+    return this->EvaluateUnoptimized(index);
+  }
+
+  SizeType
+  GetRadius() const override
+  {
+    return SizeType::Filled(1);
+  }
+
+protected:
+  LinearInterpolateSelectedNeighborsImageFunction();
+  ~LinearInterpolateSelectedNeighborsImageFunction() override;
+  void PrintSelf(std::ostream & os, Indent indent) const override;
+
+private:
+  NeighborCheckFunctionType m_NeighborCheck = NeighborCheckFunctionType();
+  /** Evaluate interpolator at image index position. */
+  virtual inline OutputType EvaluateUnoptimized(
+    const ContinuousIndexType & index) const;
+
+  /** \brief A method to generically set all components to zero
+   */
+  template<typename RealTypeScalarRealType>
+    void
+    MakeZeroInitializer(const TInputImage * const inputImagePtr,
+      VariableLengthVector<RealTypeScalarRealType> & tempZeros) const
+      {
+      // Variable length vector version to get the size of the pixel correct.
+      typename TInputImage::IndexType idx;
+      idx.Fill(0);
+      const typename TInputImage::PixelType & tempPixel = inputImagePtr->GetPixel(idx);
+      const unsigned int sizeOfVarLengthVector = tempPixel.GetSize();
+      tempZeros.SetSize(sizeOfVarLengthVector);
+      tempZeros.Fill(NumericTraits< RealTypeScalarRealType >::ZeroValue());
+      }
+
+  template<typename RealTypeScalarRealType>
+    void
+    MakeZeroInitializer(const TInputImage * const itkNotUsed( inputImagePtr ),
+      RealTypeScalarRealType & tempZeros) const
+      {
+      // All other cases
+      tempZeros = NumericTraits< RealTypeScalarRealType >::ZeroValue();
+      }
+
+};
+} // end namespace itk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "itkLinearInterpolateSelectedNeighborsImageFunction.hxx"
+#endif
+
+#endif

include/itkLinearInterpolateSelectedNeighborsImageFunction.hxx

@@ -0,0 +1,134 @@
+/*=========================================================================
+ *
+ *  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 itkLinearInterpolateSelectedNeighborsImageFunction_hxx
+#define itkLinearInterpolateSelectedNeighborsImageFunction_hxx
+
+#include "itkConceptChecking.h"
+#include "itkLinearInterpolateSelectedNeighborsImageFunction.h"
+
+#include "itkMath.h"
+
+namespace itk
+{
+  template< typename TInputImage, typename TCoordRep, typename TNeighborCheckFunction >
+  LinearInterpolateSelectedNeighborsImageFunction< TInputImage, TCoordRep, TNeighborCheckFunction >
+::LinearInterpolateSelectedNeighborsImageFunction() = default;
+
+  template< typename TInputImage, typename TCoordRep, typename TNeighborCheckFunction >
+  LinearInterpolateSelectedNeighborsImageFunction< TInputImage, TCoordRep, TNeighborCheckFunction >
+::~LinearInterpolateSelectedNeighborsImageFunction() = default;
+
+  template< typename TInputImage, typename TCoordRep, typename TNeighborCheckFunction >
+  typename LinearInterpolateSelectedNeighborsImageFunction< TInputImage, TCoordRep, TNeighborCheckFunction >
+::OutputType
+  LinearInterpolateSelectedNeighborsImageFunction< TInputImage, TCoordRep, TNeighborCheckFunction >
+::EvaluateUnoptimized(const ContinuousIndexType & index) const
+{
+  // Avoid the smartpointer de-reference in the loop for
+  // "return m_InputImage.GetPointer()"
+  const TInputImage * const inputImagePtr = this->GetInputImage();
+
+  // Compute base index = closest index below point
+  // Compute distance from point to base index
+
+  IndexType baseIndex;
+  InternalComputationType    distance[ImageDimension];
+  for ( unsigned int dim = 0; dim < ImageDimension; ++dim )
+    {
+    baseIndex[dim] = Math::Floor< IndexValueType >(index[dim]);
+    distance[dim] = index[dim] - static_cast< InternalComputationType >( baseIndex[dim] );
+    }
+  // The interpolated value is the weighted sum of each of the surrounding
+  // neighbors. The weight for each neighbor is the fraction overlap
+  // of the neighbor pixel with respect to a pixel centered on point.
+
+  // When RealType is VariableLengthVector, 'value' will be resized properly
+  // below when it's assigned again.
+  Concept::Detail::UniqueType< typename NumericTraits< RealType >::ScalarRealType >();
+
+  RealType value;
+  // Initialize variable "value" with overloaded function so that
+  // in the case of variable length vectors the "value" is initialized
+  // to all zeros of length equal to the InputImagePtr first pixel length.
+  this->MakeZeroInitializer( inputImagePtr, value );
+
+  Concept::Detail::UniqueType< typename NumericTraits< InputPixelType >::ScalarRealType >();
+
+  // Number of neighbors used in the interpolation
+  constexpr unsigned long numberOfNeighbors = 1 << TInputImage::ImageDimension;
+  InternalComputationType TotalOverlap = 0.0;
+
+  for ( unsigned int counter = 0; counter < numberOfNeighbors; ++counter )
+    {
+    InternalComputationType overlap = 1.0;    // Fraction overlap
+    unsigned int upper = counter;  // Each bit indicates upper/lower neighbour
+    IndexType    neighIndex( baseIndex );
+
+    // Get neighbor index and overlap fraction
+    for ( unsigned int dim = 0; dim < ImageDimension; ++dim )
+      {
+      if ( upper & 1 )
+        {
+        ++(neighIndex[dim]);
+        // Take care of the case where the pixel is just
+        // in the outer upper boundary of the image grid.
+        if ( neighIndex[dim] > this->m_EndIndex[dim] )
+          {
+          neighIndex[dim] = this->m_EndIndex[dim];
+          }
+        overlap *= distance[dim];
+        }
+      else
+        {
+        // Take care of the case where the pixel is just
+        // in the outer lower boundary of the image grid.
+        if ( neighIndex[dim] < this->m_StartIndex[dim] )
+          {
+          neighIndex[dim] = this->m_StartIndex[dim];
+          }
+        overlap *= 1.0 - distance[dim];
+        }
+
+      upper >>= 1;
+      }
+    InputPixelType nval = inputImagePtr->GetPixel(neighIndex);
+
+    if (m_NeighborCheck(nval))
+      {
+      TotalOverlap += overlap;
+      value += static_cast< RealType >( nval ) * overlap;
+      }
+    }
+  if (TotalOverlap == 0)
+    {
+      itkWarningMacro("LinearInterpolateSelectedNeighbors - no valid neighbors");
+      return(static_cast< OutputType >(inputImagePtr->GetPixel(baseIndex)));
+    }
+  return ( static_cast< OutputType >( value/TotalOverlap ) );
+}
+
+  template< typename TInputImage, typename TCoordRep, typename TNeighborCheckFunction >
+void
+  LinearInterpolateSelectedNeighborsImageFunction< TInputImage, TCoordRep, TNeighborCheckFunction >
+::PrintSelf(std::ostream & os, Indent indent) const
+{
+  this->Superclass::PrintSelf(os, indent);
+}
+} // end namespace itk
+
+#endif

include/itkPhysicalCentralDifferenceImageFunction.h

@@ -78,7 +78,13 @@ class ITK_TEMPLATE_EXPORT PhysicalCentralDifferenceImageFunction
   using PointType = typename Superclass::PointType;
 
   /** Linear interpolate function type alias. */
-  using InterpolateImageFunctionType = LinearInterpolateImageFunction<TInputImage, TCoordRep>;
+  /** Type of the Interpolator class */
+  using InterpolatorType = InterpolateImageFunction< InputImageType, TCoordRep >;
+  using DefaultInterpolatorType = LinearInterpolateImageFunction< InputImageType, TCoordRep >;
+
+  /** Get/set the Interpolator. */
+  itkSetObjectMacro( Interpolator, InterpolatorType );
+  itkGetConstObjectMacro( Interpolator, InterpolatorType );
 
   /** Set the input image.
    * \warning this method caches BufferedRegion information.
@@ -137,7 +143,8 @@ class ITK_TEMPLATE_EXPORT PhysicalCentralDifferenceImageFunction
   PrintSelf(std::ostream & os, Indent indent) const override;
 
 private:
-  typename InterpolateImageFunctionType::Pointer m_Interpolator;
+  typename InterpolatorType::Pointer m_Interpolator;
+
 };
 
 } // end namespace itk