MantisBT - ITK
View Issue Details
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 
Resolution Date
Sprint
Sprint Statusbacklog
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 [^]
Steps To Reproduce
Additional Information
TagsNo tags attached.
Relationships
Attached Files
Issue History
Date ModifiedUsernameFieldChange
2008-03-19 08:10Tomas KazmarNote Added: 0010852
2008-03-19 13:24Luis IbanezAssigned Touser390 => Luis Ibanez
2008-03-20 20:44Luis IbanezNote Added: 0010876
2008-03-21 10:57Luis IbanezNote Added: 0010877
2008-03-21 10:59Luis IbanezStatusassigned => resolved
2008-03-21 10:59Luis IbanezResolutionopen => fixed
2008-03-21 10:59Luis IbanezNote Added: 0010878
2008-03-24 13:45Luis IbanezNote Added: 0010890
2008-03-24 13:45Luis IbanezStatusresolved => confirmed
2008-03-24 13:52Luis IbanezNote Added: 0010891
2008-03-25 06:25Tomas KazmarNote Added: 0010902
2010-10-21 10:37Hans JohnsonSprint Status => backlog
2010-10-21 10:37Hans JohnsonNote Added: 0022557
2010-10-21 10:37Hans JohnsonStatusconfirmed => closed

Notes
(0010852)
Tomas Kazmar   
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   
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   
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   
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   
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   
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   
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   
2010-10-21 10:37   
already fixed