[Insight-users] bug in BSplineInterpolationWeighFunction?

[Stefan Klein]

Hi,

In addition to my previous mail: in the BSplineInterpolateImageFunction the 
startIndex is computed differently than in the 
BSplineInterpolationWeightFunction.

Method: DetermineRegionOfSupport

The following code does the work:


   long indx;

// compute the interpolation indexes
   for (unsigned int n = 0; n< ImageDimension; n++)
     {
     if (splineOrder & 1)     // Use this index calculation for odd splineOrder
       {
       indx = (long)vcl_floor((float)x[n]) - splineOrder / 2;
       for (unsigned int k = 0; k <= splineOrder; k++)
         {
         evaluateIndex[n][k] = indx++;
         }
       }
     else                       // Use this index calculation for even 
splineOrder
       {
       indx = (long)vcl_floor((float)(x[n] + 0.5)) - splineOrder / 2;
       for (unsigned int k = 0; k <= splineOrder; k++)
         {
         evaluateIndex[n][k] = indx++;
         }
       }
     }


I think this behaves exactly like the solution I proposed in my previous 
mail, but distinguishes between the cases of odd and even spline order (so 
might be slightly less efficient?; especially since the function is often 
invoked many times).

Let me know what you think about it.

Regards!
Stefan.



At 15:01 12/10/06, Stefan Klein wrote:
>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.
>
>_______________________________________________
>Insight-users mailing list
>Insight-users at itk.org
>http://www.itk.org/mailman/listinfo/insight-users

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20061012/02d34bef/attachment-0001.html