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

[Stefan Klein]

Hi Tom,

Thanks, I should have checked the bugtracker!

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()?

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
>>
>>