[Insight-users] bug in BSplineInterpolationWeightFunction, BSplineFloor()

[Tom Vercauteren]

Hi Stefan,

> Until there is consensus on the best implementation of the fast
> BSplineFloor(), wouldn't it be best to opt for the safe choice, vcl_floor()?

I would vote for such a fix to.

By the way, there's another place in ITK where a similar function is used.

Regards,
Tom


>
> It's quite a nasty bug, like it is now.
>
> Kind regards,
> Stefan.
>
>
>
> Tom Vercauteren wrote:
>>
>> Hi Stefan,
>>
>> You'll find more information about that bug here;
>> http://www.itk.org/Bug/view.php?id=2078
>> http://www.itk.org/Bug/view.php?id=5692
>>
>> http://sourceforge.net/mailarchive/forum.php?thread_name=28392e8b0808050138r31574c90wddd19c0072a7da3c%40mail.gmail.com&forum_name=vxl-maintainers
>>
>> Unfortunately, a consensus doesn't seem to exist right now.
>>
>> Regards,
>> Tom Vercauteren
>>
>> On Mon, Sep 8, 2008 at 12:42 PM, Stefan Klein <s.klein at erasmusmc.nl>
>> wrote:
>>>
>>> Hi,
>>>
>>> It seems there is a bug in the BSplineFloor() function defined in
>>> itkBSplineInterpolationWeightFunction.txx.
>>>
>>> It floors 97.999999 to 98...
>>>
>>> The code in the attachment reproduces the problem and illustrates how it
>>> may
>>> lead to severe crashes when using the BSplineDeformableTransform (that's
>>> how
>>> I found it).
>>>
>>> Can anybody reproduce the problem? I'm using Visual C++ 2008, Windows XP
>>> 32bit, ITK 3.8.
>>>
>>> It can be easily fixed by changing in
>>> itkBSplineInterpolationWeightFunction.txx the following lines:
>>>
>>>  union { unsigned int hilo[2]; double d; } u;
>>>  u.d = x + 103079215104.0;
>>>  return (int)((u.hilo[1]<<16)|(u.hilo[0]>>16));
>>>
>>> into:
>>>
>>>  return static_cast<int>(vcl_floor(x));
>>>
>>> But that might give a performance loss, as the comment explains:
>>>
>>> // The 'floor' function on x86 and mips is many times slower than these
>>> // and is used a lot in this code, optimize for different CPU
>>> architectures
>>>
>>> Does anybody have a better solution?
>>>
>>> Kind regards,
>>> Stefan
>>>
>>>
>>>
>>>
>>> //#define TESTBSPLINEFLOOR
>>>
>>> #include "itkBSplineDeformableTransform.h"
>>> #include "itkNumericTraits.h"
>>> #include <iostream>
>>>
>>> int main( int argc, char ** argv )
>>> {
>>>  const unsigned int Dimension = 3;
>>>  typedef double CoordRepType;
>>>
>>>  typedef itk::BSplineDeformableTransform<CoordRepType, Dimension, 3>
>>>   TransformType;
>>>  typedef TransformType::ParametersType ParametersType;
>>>  typedef TransformType::SizeType SizeType;
>>>  typedef TransformType::RegionType RegionType;
>>>  typedef TransformType::OriginType OriginType;
>>>  typedef TransformType::SpacingType SpacingType;
>>>  typedef TransformType::InputPointType PointType;
>>>  typedef PointType::VectorType VectorType;
>>>
>>>  TransformType::Pointer transform = TransformType::New();
>>>  SizeType size;
>>>  OriginType origin;
>>>  SpacingType spacing;
>>>  RegionType region;
>>>  ParametersType parameters;
>>>  PointType inpoint;
>>>  PointType outpoint;
>>>  VectorType vec;
>>>
>>>  for (unsigned int i =0; i < Dimension ; ++i)
>>>  {
>>>   size[i]=100;
>>>   origin[i]=0.0;
>>>   spacing[i]=1.0;
>>>  }
>>>  region.SetSize( size);
>>>  transform->SetGridRegion( region );
>>>  transform->SetGridOrigin( origin );
>>>  transform->SetGridSpacing( spacing );
>>>  unsigned int N = transform->GetNumberOfParameters();
>>>  parameters.SetSize(N);
>>>  parameters.Fill(0.0);
>>>  for ( unsigned int p = 0; p < N; ++p )
>>>  {
>>>   parameters[p]=static_cast<double>(p);
>>>  }
>>>
>>>  transform->SetParameters( parameters );
>>>
>>>  // the following explains the problem
>>>  int testint = ::BSplineFloor(97.999999);
>>>  std::cerr << "::BSplineFloor(97.999999) = " << testint << std::endl;
>>>
>>>  // consequently, this loop crashes at i=1
>>>  inpoint[0]=70.0;
>>>  inpoint[1]=50.0;
>>>  inpoint[2]=100.0;
>>>  for ( unsigned int i = 0; i <= 100; i++)
>>>  {
>>>   inpoint[2] = 97.900000 + 0.099999*static_cast<double>(i);
>>>   std::cerr << "i =" << i << " inpoint " << inpoint[2] << " -> ";
>>>   outpoint = transform->TransformPoint(inpoint);
>>>   vec=outpoint-inpoint;
>>>   std::cerr << outpoint[2] << "\t" << "vec=" << vec << std::endl;
>>>  }
>>>
>>>  return 0;
>>>
>>> } // end main
>>>
>>>
>>>
>>> CMAKE_MINIMUM_REQUIRED(VERSION 2.6)
>>>
>>> # This project is intended to be built outside the Insight source tree
>>> PROJECT(BugTest)
>>>
>>> # Find ITK.
>>> FIND_PACKAGE(ITK)
>>> IF(ITK_FOUND)
>>>  INCLUDE(${ITK_USE_FILE})
>>> ELSE(ITK_FOUND)
>>>  MESSAGE(FATAL_ERROR
>>>         "Cannot build without ITK.  Please set ITK_DIR.")
>>> ENDIF(ITK_FOUND)
>>>
>>> ADD_EXECUTABLE(bugtest bugtest.cxx)
>>>
>>> TARGET_LINK_LIBRARIES(bugtest
>>> ITKBasicFilters
>>> ITKNumerics
>>> ITKIO
>>> ITKCommon )
>>>
>>> INSTALL_TARGETS(/. bugtest)
>>>
>>> _______________________________________________
>>> Insight-users mailing list
>>> Insight-users at itk.org
>>> http://www.itk.org/mailman/listinfo/insight-users
>>>
>>>
>