BUG: itkPolyLineParametricPath infinite loop for degenerate paths

itkPolyLineParametricPath inherits EvalateDerivative and IncrementInput from its
superclass: itkParametricPath. IncrementInput uses a default
timestep and a derivative estimated from that default timestep to implement a
tooBig/tooSmall search strategy to find the next pixel along the path.

This strategy will result in oscillation between tooBig and tooSmall for
degenerate itkPolyLineParametricPath paths, such as those with a lot of
verticies at the beginning less than one image index away from each other
followed by a single vertex at the end that is multiple image indices away.

This patch overrides the inherited IncrementInputs and EvaluateDerivative
function and takes advantage of the fact that it is possible to calculate an
instantaneous derivative which allows to calculate the exact timestep to effect
a change in one pixel.

A test is also added for one of the degenerate paths mentioned in paragraph 2.

The ExtractOrthogonalSwatch2DImageFilterTest baseline image was updated due to
the old tooSmall/tooLarge search strategy not always honoring the 8-connected
neighborhood.

Change-Id: Ibe2d5b8fb913e6311ee0f3b977fdadf1321f0653

Modules/Filtering/Path/include/itkPolyLineParametricPath.h

@@ -119,6 +119,17 @@ class ITK_EXPORT PolyLineParametricPath:public
   /** Return the container of Vertices as a const object. */
   itkGetConstObjectMacro(VertexList, VertexListType);
 
+  /** This function overrides the superclass IncrementInput and calculates
+   *  the next pixel along the path to visit by using the instantaneous
+   *  partial derivatives to calculate the timestep needed to move along the
+   *  path by one pixel */
+  virtual OffsetType IncrementInput(InputType & input) const;
+
+  /** This function overrides the superclass EvaluateDerivative and instead
+   *  calculates the instantaneous derivative of input by taking the index
+   *  of the previous and next integral timepoints and subtracting them */
+  virtual VectorType EvaluateDerivative(const InputType & input) const;
+
 protected:
   PolyLineParametricPath();
   ~PolyLineParametricPath(){}

Modules/Filtering/Path/include/itkPolyLineParametricPath.hxx

@@ -36,7 +36,7 @@ PolyLineParametricPath< VDimension >
     }
 
   const VertexType vertex0 = static_cast<const VertexListType*>(this->m_VertexList)->ElementAt( int(input) );
-  const VertexType vertex1 = static_cast<const VertexListType*>(this->m_VertexList)->ElementAt( int(input + 1.0) );
+  const VertexType vertex1 = static_cast<const VertexListType*>(this->m_VertexList)->ElementAt( int(input) + 1 );
 
   const double fractionOfLineSegment = input - int(input);
 
@@ -53,12 +53,108 @@ PolyLineParametricPath< VDimension >
   return output;
 }
 
-//template<unsigned int VDimension>
-//typename PolyLineParametricPath<VDimension>::VectorType
-//PolyLineParametricPath<VDimension>
-//::EvaluateDerivative(const InputType & input) const
-//{
-//}
+template<unsigned int VDimension>
+typename PolyLineParametricPath<VDimension>::VectorType
+PolyLineParametricPath<VDimension>
+::EvaluateDerivative(const InputType & input) const
+{
+  //Get next integral time-point
+  const InputType nextTimepoint = std::min(std::floor(input + 1.0), this->EndOfInput());
+
+  //Get previous integral time-point
+  const InputType previousTimepoint = nextTimepoint - 1.0;
+
+  //Calculate the continuous index for both points
+  const ContinuousIndexType nextIndex = this->Evaluate(nextTimepoint);
+  const ContinuousIndexType prevIndex = this->Evaluate(previousTimepoint);
+
+  //For some reason, there's no way to convert ContinuousIndexType to VectorType
+  VectorType partialDerivatives;
+  for (unsigned int i = 0; i < VDimension; ++i)
+    {
+    partialDerivatives[i] = nextIndex[i] - prevIndex[i];
+    }
+
+  return partialDerivatives;
+}
+
+template<unsigned int VDimension>
+typename PolyLineParametricPath<VDimension>::OffsetType
+PolyLineParametricPath<VDimension>
+::IncrementInput(InputType & input) const
+{
+  //Save this input index since we will use it to calculate the offset
+  const OutputType originalIndex = this->EvaluateToIndex(input);
+
+  InputType potentialTimestep = itk::NumericTraits< InputType >::ZeroValue();
+  bool timeStepSmallEnough = false;
+  while (!timeStepSmallEnough)
+    {
+    if (input == this->EndOfInput())
+      {
+      return this->GetZeroOffset();
+      }
+
+    //Check to make sure we aren't already at a place with an offset of 1 pixel
+    const OutputType potentialIndex = this->EvaluateToIndex(input);
+    //For some reason, there's no way to convert OutputType to OffsetType
+    OffsetType offset;
+    for (unsigned int i = 0; i < VDimension; ++i)
+      {
+      offset[i] = potentialIndex[i] - originalIndex[i];
+      }
+
+    if (offset != this->GetZeroOffset())
+      {
+      return offset;
+      }
+    //Get the derivative at the current input
+    VectorType partialDerivatives = this->EvaluateDerivative(input);
+
+    //Find the largest of all the partial derivatives
+    unsigned int maxPartialDerivativeIndex = 0;
+    for (unsigned int i = 1; i < VDimension; ++i)
+      {
+      if (std::abs(partialDerivatives[i]) > std::abs(partialDerivatives[maxPartialDerivativeIndex]))
+        {
+        maxPartialDerivativeIndex = i;
+        }
+      }
+
+    //Calculate the timestep required to effect a 1 pixel change
+    potentialTimestep = 1.0/std::abs(partialDerivatives[maxPartialDerivativeIndex]);
+
+    //Check to make sure the timestep doesn't put the input past the next integral timestep
+    //(since the derivatives can change)
+    if (input + potentialTimestep > std::floor(input + 1.0))
+      {
+      input = std::floor(input + 1.0); //Set the input to the next integral time-step
+      }
+    else
+      {
+      timeStepSmallEnough = true;
+      }
+    }
+
+  //Finalize the potential timestep into THE timestep
+  const InputType timestep = potentialTimestep;
+
+  //Get the index at the next timestep so we can calculate the offset
+  const OutputType nextIndex = this->EvaluateToIndex(input + timestep);
+
+  //For some reason, there's no way to convert OutputType to OffsetType
+  OffsetType offset;
+  for (unsigned int i = 0; i < VDimension; ++i)
+    {
+    offset[i] = nextIndex[i] -  originalIndex[i];
+    }
+
+  //Update input timestep
+  input += timestep;
+
+  //Return the offset
+  return offset;
+}
 
 template< unsigned int VDimension >
 PolyLineParametricPath< VDimension >

Modules/Filtering/Path/test/itkPolyLineParametricPathTest.cxx

@@ -84,6 +84,35 @@ int itkPolyLineParametricPathTest(int, char* [])
   offset = path->IncrementInput( input );
   std::cout << "Incrementing the input from 0.5 to " << input << ":  " << offset << std::endl;
 
+  //Test a degenerate path
+  std::cout << "Generating degenerate path" << std::endl;
+  PathType::Pointer path2 = PathType::New();
+
+  //Add a bunch of points closely spaced together
+  for (double k = 0; k < 10; k+=0.1)
+    {
+      v.Fill(k);
+      path2->AddVertex(v);
+    }
+
+  //Add a point that is very far away from the first points
+  v.Fill(100);
+  path2->AddVertex(v);
+  PathType::InputType path2Input = path2->StartOfInput();
+
+  PathType::OffsetType zeroOffset;
+  zeroOffset.Fill(0);
+
+  std::cout << "Starting degenerate path test" << std::endl;
+  while (true)
+    {
+    offset = path2->IncrementInput(path2Input);
+    if (offset == zeroOffset)
+      {
+      break;
+      }
+    }
+
   if (passed)
     {
     std::cout << "PolyLineParametricPath tests passed" << std::endl;

Testing/Data/Baseline/BasicFilters/ExtractOrthogonalSwath2DImageFilterTest.png.md5

@@ -1 +1 @@
-e28ea77794e99f1da041ac67fa55e7ed
+b79c0a553750ed6fce228b23df79f7d8