[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
   // want to test more than 32bit range
   const SizeValueT largeIntegerOffset( 4503599627369471 );// double mantissa can still resolve points between image gridpoints
   const OffsetValueT chunklocation = largeIntegerOffset - chunksize; // test starts belowlargeIntegerOffset

   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\
   // debugger:testLocation[0] ==  4503599627369371

   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...