[Insight-users] problems with PolyLineParametricPath and PathIterator

[Jeffrey Duda]

I've been having some problems using the PathIterator.  I'm using the 
PolyLineParametricPath, but the problem occurs in IncrementInput() which 
is defined in the ParametricPath.  The problem occurs when there are a 
large number of vertices in a row, that are all in the same voxel.  When 
the inputStepSize doesn't result in a new voxel, it is doubled.  This 
quickly leads to a large inputStepSize.  Its been my experience that 
this has a tendency to then cause the inputStepSize to bounce back and 
forth between 2 values until the max number of iterations is reached and 
an exception is thrown.  While this problem isn't necessarily that 
common, I'm working with sets that contain thousands of paths (dti fiber 
tractography) and its almost a guarantee that this problem occurs in at 
least 1 path per set.

The IncrementInput() function is virtual so it would be possible to 
write a new class that replaces it, but that seems like a bit of a waste 
since the PolyLineParametricPath is exactly what I need, and it would be 
almost a complete replication of code.  Things work just fine if 
IncrementInput is modified to take an even slower but steady approach by 
making the following modification:

"inputStepSize *= 2.0;" becomes "inputStepSize += m_DefaultInputStepSize;"
and
"inputStepSize /= 1.5;" becomes "inputStepSize -= 
m_DefaultInputStepSize/0.5;"

In summary:
1) Why was the decision made to double the step size when it is too small?
2) Would it be better to more slowly increment it, giving a more stable 
iterator?
3) Should I stop complaining and just write my own class?

thanks,
-jeffrey duda

Here is some simple code that breaks the PathIterator

typedef itk::Image<int, 3> ImageType;
typedef itk::PolyLineParametricPath<3> PathType;
typedef itk::PathIterator< ImageType, PathType > IteratorType;
typedef PathType::ContinuousIndexType VertexType;

ImageType::Pointer img = ImageType::New();
ImageType::SizeType imgSize;
imgSize[0] = imgSize[1] = 10;
imgSize[2] = 2;
img->SetRegions( imgSize );
img->Allocate();

PathType::Pointer path = PathType::New();
VertexType pts[12];

pts[0][0] = 3.5396; pts[0][1] = 4.5323; pts[0][2] = 0.42581;
pts[1][0] = 4.0457; pts[1][1] = 5.1173; pts[1][2] = 0.47410;
pts[2][0] = 4.4586; pts[2][1] = 5.6914; pts[2][2] = 0.52771;
pts[3][0] = 4.6246; pts[3][1] = 5.9129; pts[3][2] = 0.60156;
pts[4][0] = 4.6535; pts[4][1] = 5.9860; pts[4][2] = 0.67940;
pts[5][0] = 4.6541; pts[5][1] = 6.0289; pts[5][2] = 0.75745;
pts[6][0] = 4.6705; pts[6][1] = 6.0968; pts[6][2] = 0.83538;
pts[7][0] = 4.7229; pts[7][1] = 6.2189; pts[7][2] = 0.91279;
pts[8][0] = 4.8351; pts[8][1] = 6.4354; pts[8][2] = 0.98848;
pts[9][0] = 5.0322; pts[9][1] = 6.7869; pts[9][2] = 1.0598;
pts[10][0] = 5.1955; pts[10][1] = 7.2585; pts[10][2] = 1.1254;
pts[11][0] = 4.9012; pts[11][1] = 7.5336; pts[11][2] = 1.1968;

for (unsigned int i=0; i<12; i++)
   path->AddVertex( pts[i] );

// Create iterator following path through image
IteratorType iter( img, path );

while (!iter.IsAtEnd())
   {
     std::cout << iter.GetPathPosition() << " -> " << std::flush;
     try
       {
	++iter;
       }
     catch( itk::ExceptionObject & err )
       {
	std::cout << "ExceptionObject caught !" << std::endl;
	std::cout << err << std::endl;
	return -1;
        }

   }