View Issue Details [ Jump to Notes ] [ Print ]
IDProjectCategoryView StatusDate SubmittedLast Update
0004021ITKpublic2006-11-06 02:282010-10-21 10:37
ReporterStefan Klein 
Assigned ToLuis Ibanez 
PriorityhighSeveritymajorReproducibilityalways
StatusclosedResolutionfixed 
PlatformOSOS Version
Product Version 
Target VersionFixed in Version 
Summary0004021: BSplineInterpolationWeightFunction computes startIndex incorrectly
DescriptionIn 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.

See also:

http://public.kitware.com/pipermail/insight-users/2006-October/019744.html [^]

http://public.kitware.com/pipermail/insight-users/2006-October/019747.html [^]
TagsNo tags attached.
Resolution Date
Sprint
Sprint Statusbacklog
Attached Files

 Relationships

  Notes
(0010852)
Tomas Kazmar (reporter)
2008-03-19 08:10

 I think support indices are computed correctly, so we should not change them. The error is in incorrect, as Stefan reports, in incorrect positions used for B-spline kernel computation which uses positions outside support. In case of second order B-spline, current implementation uses -1..2 but us should use -1.5..1.5 as the support is symmetric around zero.

The patch is this simple (itkBSplineInterpolationWeightFunction.txx):
@@ -154,7 +154,8 @@ void BSplineInterpolationWeightFunction<
   Matrix<double,SpaceDimension,SplineOrder + 1> weights1D;
   for ( j = 0; j < SpaceDimension; j++ )
     {
- double x = index[j] - static_cast<double>( startIndex[j] );
+ const double shiftEven = 0.5*(SplineOrder%2-1);
+ double x = index[j] - static_cast<double>( startIndex[j] ) + shiftEven;

     for( k = 0; k <= SplineOrder; k++ )
       {
(0010876)
Luis Ibanez (manager)
2008-03-20 20:44

 Index: itkBSplineInterpolationWeightFunction.txx
===================================================================
RCS file: /cvsroot/Insight/Insight/Code/Common/itkBSplineInterpolationWeightFunction.txx,v
retrieving revision 1.13
diff -r1.13 itkBSplineInterpolationWeightFunction.txx
150c150
< BSplineFloor( index[j] - static_cast<double>( SplineOrder / 2 ) ) );
---
> BSplineFloor( index[j] - static_cast<double>( SplineOrder - 1 ) / 2.0 ) );

The fix turned out to be needed not in the computation of the kernel argument but in the computation of the startIndex.
(0010877)
Luis Ibanez (manager)
2008-03-21 10:57

 A Fix has been committed to the CVS repository:

http://www.itk.org/cgi-bin/viewcvs.cgi/Code/Common/itkBSplineInterpolationWeightFunction.txx?root=Insight&r1=1.13&r2=1.14&sortby=date [^]

Tests were added for 1D resampling using Spline orders 1,2 and 3:
http://www.itk.org/cgi-bin/viewcvs.cgi/Testing/Code/Common/itkBSplineDeformableTransformTest2.cxx?root=Insight&sortby=date&view=log [^]

A test of kernel symmetry has been added:
http://www.itk.org/cgi-bin/viewcvs.cgi/Testing/Code/Common/itkBSplineInterpolationWeightFunctionTest.cxx?root=Insight&r1=1.4&r2=1.6&sortby=date [^]
(0010878)
Luis Ibanez (manager)
2008-03-21 10:59

 A fix has been committed along with additional tests.

We are waiting for Stefan to confirm that the fix work for his reported problem. Then we will close the bug.
(0010890)
Luis Ibanez (manager)
2008-03-24 13:45

 Reopen the bug, due the following email from Tomas;

IMHO, the patch you decided to use is not correct.
The tests shadow the error as they define the grid, so that it works.
The grid is defined the same irrespective of the spline order.
There are five nodes on image and three outside the image
and therefore you need to correct it changing the starting index.

 const unsigned int numberOfGridNodesOutsideTheImageSupport = 3;
 const unsigned int numberOfGridNodesInsideTheImageSupport = 5;
 const unsigned int numberOfGridNodes = numberOfGridNodesInsideTheImageSupport +
                       numberOfGridNodesOutsideTheImageSupport;
                       As a user I would expect that I setup the grid with n nodes outside the image
for a spline of order n.

 const unsigned int numberOfGridNodesOutsideTheImageSupport = VSplineOrder;
 const unsigned int numberOfGridNodesInsideTheImageSupport = 5;
 const unsigned int numberOfGridNodes = numberOfGridNodesInsideTheImageSupport +
                       numberOfGridNodesOutsideTheImageSupport;


Maybe this way has some implementation advantages, say in implementation
of B-spline derivatives, I don't know. So the question: is the
BSplineDeformableTransform meant to be used always with the grid with three
nodes outside the image? Consider this, as after inclusion in ITK, there will
be no way to correct it afterwards because of the backwards compatibility.


Regards,
Tomas
(0010891)
Luis Ibanez (manager)
2008-03-24 13:52

 Original email from Tomas in the ITK developers list:
http://www.itk.org/mailman/private/insight-developers/2008-March/010079.html [^]
(0010902)
Tomas Kazmar (reporter)
2008-03-25 06:25

 Thank you for reopening, now with the new test I see it works well and it successfully warps my images for both odd and even order splines.

Sorry for the confusion.
(0022557)
Hans Johnson (developer)
2010-10-21 10:37

 already fixed

 Issue History
Date Modified Username Field Change
2008-03-19 08:10 Tomas Kazmar Note Added: 0010852
2008-03-19 13:24 Luis Ibanez Assigned To user390 => Luis Ibanez
2008-03-20 20:44 Luis Ibanez Note Added: 0010876
2008-03-21 10:57 Luis Ibanez Note Added: 0010877
2008-03-21 10:59 Luis Ibanez Status assigned => resolved
2008-03-21 10:59 Luis Ibanez Resolution open => fixed
2008-03-21 10:59 Luis Ibanez Note Added: 0010878
2008-03-24 13:45 Luis Ibanez Note Added: 0010890
2008-03-24 13:45 Luis Ibanez Status resolved => confirmed
2008-03-24 13:52 Luis Ibanez Note Added: 0010891
2008-03-25 06:25 Tomas Kazmar Note Added: 0010902
2010-10-21 10:37 Hans Johnson Sprint Status => backlog
2010-10-21 10:37 Hans Johnson Note Added: 0022557
2010-10-21 10:37 Hans Johnson Status confirmed => closed

Copyright © 2000 - 2018 MantisBT Team