MantisBT - ITK
View Issue Details
IDProjectCategoryView StatusDate SubmittedLast Update
0008311ITKpublic2008-12-17 20:022011-05-18 17:33
ReporterDaniel Betz 
Assigned Tokentwilliams 
PrioritynormalSeveritymajorReproducibilityalways
StatusfeedbackResolutionreopened 
PlatformOSOS Version
Product VersionITK-3-10 
Target VersionFixed in Version 
Resolution Date
Sprint
Sprint Statusbacklog
Summary0008311: PolyLineParametricPath converges too slow
DescriptionThe ParametricPath has a methode IncrementInput() which is used by the PathIterator to go along a path until the next pixel is reached. The implementation converges far too slow (possibly sometimes never) and the allowed maximum of 10000 iterations get reached so that an exception is thrown. Some example pathes that fail with this or the related implementation http://www.seas.upenn.edu/~jtduda/code/itkSlowPolyLineParametricPath.txx [^]
is attached.

This bug is related to:
http://www.itk.org/Bug/bug_view_page.php?bug_id=5371&history=1 [^]

and this posting from the itk user mailinglist:
http://www.itk.org/pipermail/insight-users/2005-July/013922.html [^]
Steps To Reproduce
Additional InformationThe solution is to use interval bisection, once an upper boundary was found through doubling the input step size.

A new implemetation of IncrementInput could look like this (This never failed for me until now):

template<unsigned int VDimension>
typename FastPolyLineParametricPath<VDimension>::OffsetType
FastPolyLineParametricPath<VDimension>
::IncrementInput(InputType & input) const
{
  unsigned iterationCount;
  bool tooSmall;
  bool tooBig;
  InputType inputStepSize;
  InputType tooSmallInputStepSize;
  InputType tooBigInputStepSize;
  InputType finalInputValue;
  OffsetType offset;
  IndexType currentImageIndex;
  IndexType nextImageIndex;
  IndexType finalImageIndex;

  iterationCount = 0;
  inputStepSize = this->GetDefaultInputStepSize();
  tooSmallInputStepSize = 0;
  tooBigInputStepSize = 0;

  // Are we already at (or past) the end of the input ?
  finalInputValue = this->EndOfInput();
  currentImageIndex = this->EvaluateToIndex( input );
  finalImageIndex = this->EvaluateToIndex( finalInputValue );
  offset = finalImageIndex - currentImageIndex;

  if( ( offset == this->GetZeroOffset() && input != this->StartOfInput() ) ||
      ( input >= finalInputValue ) )
    return this->GetZeroOffset();

  do
  {
    if( iterationCount++ > 10000 ) { itkExceptionMacro(<<"Too many iterations"); }

    nextImageIndex = this->EvaluateToIndex( input + inputStepSize );
    offset = nextImageIndex - currentImageIndex;
    tooSmall = ( offset == this->GetZeroOffset() );
    tooBig = false;

    if( tooSmall )
    {
      // increase the input step size, but don't go past the end of the input
      if( (input + inputStepSize) >= finalInputValue )
      {
        inputStepSize = finalInputValue - input;
      }
      else
      {
        tooSmallInputStepSize = inputStepSize;
        if( tooBigInputStepSize == 0 )
          inputStepSize *= 2;
        else
          inputStepSize += ( tooBigInputStepSize - tooSmallInputStepSize ) / 2;
      }
    }
    else
    {
      // Search for an offset dimension that is too big
      for( unsigned i = 0; i < VDimension && ! tooBig; ++i )
      {
        tooBig = ( offset[i] >= 2 || offset[i] <= -2 );
      }
      if( tooBig )
      {
        tooBigInputStepSize = inputStepSize;
        inputStepSize -= (tooBigInputStepSize - tooSmallInputStepSize) / 2;
      }
    }
  }
  while( tooSmall || tooBig );

  input += inputStepSize;
  return offset;
}


TagsNo tags attached.
Relationships
Attached Filescxx SomePolyLines.cxx (2,606) 2008-12-17 20:02
https://public.kitware.com/Bug/file/1930/SomePolyLines.cxx
Issue History
Date ModifiedUsernameFieldChange
2008-12-17 20:02Daniel BetzNew Issue
2008-12-17 20:02Daniel BetzFile Added: SomePolyLines.cxx
2009-01-18 11:40Daniel BetzNote Added: 0014601
2010-11-02 13:36Hans JohnsonStatusnew => assigned
2010-11-02 13:36Hans JohnsonAssigned To => kentwilliams
2010-11-22 17:43kentwilliamsNote Added: 0023500
2010-12-09 07:01Hans JohnsonSprint Status => backlog
2010-12-09 07:01Hans JohnsonNote Added: 0023835
2010-12-09 07:01Hans JohnsonStatusassigned => closed
2010-12-09 07:01Hans JohnsonResolutionopen => unable to reproduce
2010-12-26 12:29Daniel BetzNote Added: 0024328
2010-12-26 12:29Daniel BetzStatusclosed => feedback
2010-12-26 12:29Daniel BetzResolutionunable to reproduce => reopened
2011-05-18 17:33kentwilliamsNote Added: 0026545

Notes
(0014601)
Daniel Betz   
2009-01-18 11:40   
if( ( offset == this->GetZeroOffset() && input != this->StartOfInput() ) ||
       ( input >= finalInputValue ) )
     return this->GetZeroOffset();

should be replaced with

  if( offset == this->GetZeroOffset() || ( input >= finalInputValue ) )
     return this->GetZeroOffset();
(0023500)
kentwilliams   
2010-11-22 17:43   
Mr. Betz -- the version in ITK4 no longer resembles the code you posted here; I don't know if it's been refactored or if I'm looking in the wrong place for the code. At any rate, here's the current implementation.


template< unsigned int VDimension >
typename PolyLineParametricPath< VDimension >::OutputType
PolyLineParametricPath< VDimension >
::Evaluate(const InputType & input) const
{
  OutputType output;
  PointType outputPoint;
  VertexType vertex0;
  VertexType vertex1;
  double fractionOfLineSegment;

  // Handle the endpoint carefully, since there is no following vertex
  if ( input >= InputType(m_VertexList->Size() - 1) )
    {
    return m_VertexList->ElementAt(m_VertexList->Size() - 1); // the last vertex
    }

  vertex0 = m_VertexList->ElementAt( int(input) );
  vertex1 = m_VertexList->ElementAt( int(input + 1.0) );

  fractionOfLineSegment = input - int(input);

  outputPoint = vertex0 + ( vertex1 - vertex0 ) * fractionOfLineSegment;

  // For some stupid reason, there is no easy way to cast from a point to a
  // continuous index.
  for ( unsigned int i = 0; i < VDimension; i++ )
    {
    output[i] = outputPoint[i];
    }

  return output;
}
(0023835)
Hans Johnson   
2010-12-09 07:01   
Attempted to review old bug, but code looked much different than that described in the bug. It may have already been fixed, so closing this bug because no feed back has been given on how to reproduce it.
(0024328)
Daniel Betz   
2010-12-26 12:29   
Steps to reproduce:
See http://review.source.kitware.com/629 [^]

Solution:
See http://review.source.kitware.com/630 [^]
(0026545)
kentwilliams   
2011-05-18 17:33   
Those two patches are abandoned. Daniel, will you come up with a new one? This has been hanging fire for a very long time!