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

[Stefan Klein]

Hey Tom,

i found one other location:
BasicFilters/itkBilateralImageFilter.txx

by grepping on 'union'.
was that the one you have in mind?

Stefan

Tom Vercauteren wrote:
> 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
>>>>
>>>>