[Insight-users] bug in BSplineInterpolationWeighFunction?

[Stefan Klein]

Hi,

I think i found a bug in the implementation of the 
BSplineInterpolationWeightFunction.

In the Evaluate function, at lines 146-151, the startIndex of the support 
region is computed, given a continuous input index (denoted by 'index'):

   // Find the starting index of the support region
   for ( j = 0; j < SpaceDimension; j++ )
     {
     startIndex[j] = static_cast<typename IndexType::IndexValueType>(
       BSplineFloor( index[j] - static_cast<double>( SplineOrder / 2 ) ) );
     }


If I understand well, for a zero order bspline, this formula should result 
in a 'round' operation. This is not the case with the formula above:

suppose the continuous input index = 0.9 and the SplineOrder = 0.
then:

startindex = floor( 0.9 - 0/2 ) = floor (0.9) = 0.

Consequently, in line 161, where the bsplinekernel weight is evaluated, the 
result is zero.
BsplineZeroOrder( inputindex - startindex ) = BSplineZeroOrder( 0.9 ) = 0.0.
Which must be wrong, since for a zeroth order bspline the weight value 
should be nonzero everywhere in its support region.

For splineorder=1 the formula works. For splineorder=2 it goes wrong again 
and for splineorder=3 it works. So, it's wrong for even spline orders.

I suggest changing it to:

      startIndex = floor( inputIndex + bsplineoffset )

  where

     bsplineOffset =  0.5 - static_cast<double>(SplineOrder) / 2.0 ;

which results in bsplineOffset = 1/2, 0, -1/2, or -1 respectively for order 
0, 1, 2, and 3.

In the case of zero order bspline this obviously gives the desired round 
operation (floor(input+0.5)).
In general: for odd spline orders it works the same as the current formula, 
but for even orders the result is different.

Please let me know if I made a mistake, or if it's really a bug.

Kind regards,
Stefan.