[Insight-users] BSplineInterpolateImageFunction and Offset

[maquette]

hi, I have a problem using BSplineInterpolateImageFunction on an Image 
where the Index
of the buffered region is NOT (0,0). I'm talking about the start Index.

This causes ApplyMirrorBoundaryConditions to compute evaluateIndex the 
wrong way. Then it
contains indices outside the buffered region and crashes in silence.

In my opinion the class can not work as long as it only uses 
InputImage->GetBufferedRegion->GetSize() and 
BSplineInterpolateImageFunction is ignored.

My next thought was, that I did not understand ContinuousIndex right. Is 
it relative to the image origin, or relative to the buffered regions offset?

So I used the "frontend" functions provided by itk::Image for this 
example (failes) and see that the continuous index is approximately the 
same as physical point
and the integer index:

-----SNIP---------

   const unsigned int NDimensions = 2;

   typedef unsigned char                       PixelT;
   typedef itk::Image<PixelT, NDimensions>     ImageT;
   typedef ImageT::IndexType                   ImageIndexT;
   typedef ImageT::Pointer                     ImagePointerT;
   typedef ImageT::RegionType                  ImageRegionT;
   typedef ImageT::SizeType                    ImageSizeT;
   typedef ImageT::OffsetType                  ImageOffsetT;
   typedef ImageSizeT::SizeValueType           SizeValueT;
   typedef ImageIndexT::IndexValueType         IndexValueT;
   typedef ImageOffsetT::OffsetValueType       OffsetValueT;

   // test the bspline
   typedef itk::BSplineInterpolateImageFunction<ImageT> BSplineFunctionT;
   BSplineFunctionT::Pointer F = BSplineFunctionT::New();
   F->SetSplineOrder(1);

   // use a region "far out" but with small size, image is very 
elongated in dimension 0
   SizeValueT chunksize( 1024*1024 ); // 1MiB
   // about the end of signed 64bit range:
   const SizeValueT largeIntegerOffset( 4503599627369471 );// double 
mantissa can still resolve points between image gridpoints
   const OffsetValueT chunklocation = largeIntegerOffset - chunksize; // 
test starts below

   ImageOffsetT  testOffset = {{chunksize-100, 0}}; // offset within 
chunk (near the end)
   ImageIndexT   zero = {{ testlocation, 0}};        // start of chunk 
in image coordinates
   ImageSizeT    size = {{ chunksize, 1 }};         // size of chunk
   ImageIndexT   testLocation = zero + testOffset;   // start of test 
sequence in image coordinates
   ImageRegionT  largeRegion(zero,size);

   // chose a bigger largest possible region starting from 0,0 size 
largeIntegerOffset,1
   ImageIndexT   LPRzero = {{ 0, 0 }};
   ImageSizeT    LPRsize = {{ largeIntegerOffset, 1}};
   ImageRegionT  LPR( LPRzero, LPRsize );

   ImageT::Pointer large = ImageT::New();
   large->SetRegions(largeRegion);
   large->SetLargestPossibleRegion(LPR);
   large->Allocate();

   // large as input to the interpolator
   F->SetInputImage(large);

   BSplineFunctionT::ContinuousIndexType cidx;
   ImageT::PointType point;

   // OK use official functions to translate from Index to Point to 
continuous Index
   // double is used as default realtype in point and cidx
   large->TransformIndexToPhysicalPoint( testLocation, point );
   // debugger: point[0] == 4503599627369371.0
   large->TransformPhysicalPointToContinuousIndex( point, cidx );
   cidx[0] += 0.5;
   // debugger:  cidx[0] == 4503599627369371.5
   PixelT interpolate = F->EvaluateAtContinuousIndex(cidx); // CRASH!

// in BSplineInterpolateImageFunction::EvaluateAtContinuousIndex(x)
//   ...
// ApplyMirrorBoundaryConditions() seems to mess up EvaluateIndex, 
negative index occurs
//   ...
//   interpolated += w * m_Coefficients->GetPixel(coefficientIndex);

--------------------------

So to clear things up:
1) how is ContinuousIndex supposed to work? relative to what?
2) if its relative to image origin (ImageBase::TransformXToY functions 
seem to indicate this):
    can BSplineInterpolateImageFunction be fixed by
     storing the buffered region INDEX of the input image AND
     ( in EvaluateAtContinuousIndex: subtract index from x, and add it 
somewhere (where?) else XOR
       use index only within ApplyMirrorBoundaryConditions() )
3) I dont get it how m_Coefficients is to be used...