(0002729)
Nicolas Savoire (reporter)
2005-07-27 11:55

A similarly incorrect function is used in BilateralImageFilter.

(0002771)
Lydia Ng (developer)
2005-08-05 15:02

Triaged 8/5/05 tcon.
Jim will look at it.

(0002810)
Nicolas Savoire (reporter)
2005-08-08 03:45

I forgot to add that the bug affects only the i386 platform.

Replacing
  unsigned int hilo[2];
  *((double *)hilo) = x + 103079215104.0;
  return (int)((hilo[1]<<16)|(hilo[0]>>16));
by
  union { unsigned int hilo[2]; double d; } u;
  u.d = x + 103079215104.0;
  return (int)((u.hilo[1]<<16)|(u.hilo[0]>>16));
seems to fix the problem.

(0002917)
Nicolas Savoire (reporter)
2005-08-18 08:12

Another thing: the bug seems to appear only with gcc 3.4 and flag -O2/-O3

(0003285)
user390
2005-11-20 23:34

Committed to the Main branch.
http://www.itk.org/cgi-bin/viewcvs.cgi/Code/Common/itkBSplineInterpolationWeightFunction.txx?rev=1.10&root=Insight&view=log [^]
http://www.itk.org/cgi-bin/viewcvs.cgi/Code/BasicFilters/itkBilateralImageFilter.txx?rev=1.25&root=Insight&view=log [^]

(0009206)
Luis Ibanez (manager)
2007-09-21 15:56

correcting 'resolved' status to 'closed'

(0009269)
Nicolas Savoire (reporter)
2007-09-24 02:44

It seems that the returned value is wrong for inputs really close but smaller than a given integer:
BSplineFloor(0.999999) returns 1 instead of 0.

I think the issue is the addition of a large constant value which triggers a loss of precision.

(0009270)
Luis Ibanez (manager)
2007-09-24 08:34

Nicos,

Can you please provide a test program that illustrate this problem ?
We need a test in order to track the issue.

Also, is it still the case that this error is only observed
in gcc 3.4 with flags -O2/O3 ?

For the long term maintenance it seems better to use vcl_floor().

Whether this is a performance issue or not, is something that we
have to verify by using a profiling test. Are we saving 1% of the
computation time ?
is it 10% ?
is it 50% ?

For 1% of computation time savings, we probably shouldn't sacrifice
portability and maintainability...

Thanks

Luis

(0009272)
Nicolas Savoire (reporter)
2007-09-24 08:46

The problem appears on i386 (linux and OSX) platform. Compilation flags and gcc version have no influence on the result.

Test case:
#include <assert.h>

inline int BSplineFloor(double x)
{
#if defined mips || defined sparc || defined __ppc__
  return (int)((unsigned int)(x + 2147483648.0) - 2147483648U);
#elif defined i386 || defined _M_IX86
  union { unsigned int hilo[2]; double d; } u;
  u.d = x + 103079215104.0;
  return (int)((u.hilo[1]<<16)|(u.hilo[0]>>16));
#else
  return int(floor(x));
#endif
}

int main ()
{
  assert(BSplineFloor(0.999999)==0);
}

(0009279)
dbrussak (reporter)
2007-09-24 20:35

Hi,
  This particular BSplineFloor() precision issue is the one I reported here in issue 0005692:
http://www.itk.org/Bug/view.php?id=5692 [^]

I'm happy to delete/rescind my issue and just wait for resolution of this one. Is that possible?

  Daniel R.

(0009942)
Aviv Hurvitz (reporter)
2007-12-16 09:34

I've also run into the problem where BSplineFloor returns an incorrect result. This caused BSPlineDeformableTransform ::TransformPoint() to SEGFAULT.

Tested on a Pentium 4. Compiled with:
- gcc version 4.1.3 20070629 (prerelease) (Debian 4.1.2-13)
- MS Visual C++ 2005

A sample program:

>>>>>>>>>>>>>>>>>>>>>>>>>

#include <iostream>
#include <math.h>

inline int BSplineFloor(double x)
{
#if defined mips || defined sparc || defined __ppc__
    return (int)((unsigned int)(x + 2147483648.0) - 2147483648U);
#elif defined i386 || defined _M_IX86
    union { unsigned int hilo[2]; double d; } u;
    u.d = x + 103079215104.0;
    return (int)((u.hilo[1]<<16)|(u.hilo[0]>>16));
#else
    return int(floor(x));
#endif
}

int main(int argc, char *argv[])
{
    double x = 16.999995;

    int bsplineFloorResult = BSplineFloor(x);
    int standardFloor = static_cast<int>( floor(x) );

    std::cout.precision(20);
    std::cout << "x = " << x << std::endl;
    std::cout << "BSspline_floor = " << bsplineFloorResult << std::endl;
    std::cout << "Standard floor = " << standardFloor << std::endl;

    /* Prints:
       x = 16.999994999999998413
       BSspline_floor = 17
       Standard floor = 16
    */
    return 0;
}
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

(0016727)
Tom Vercauteren (developer)
2009-06-16 04:04

Fixed in HEAD:
http://www.itk.org/cgi-bin/viewcvs.cgi/Code/Common/itkBSplineInterpolationWeightFunction.txx?root=Insight&r1=1.15&r2=1.16&sortby=date [^]
http://www.itk.org/cgi-bin/viewcvs.cgi/Code/BasicFilters/itkBilateralImageFilter.txx?root=Insight&r1=1.29&r2=1.30&sortby=date [^]

Issue History
Date Modified	Username	Field	Change
2007-09-21 15:56	Luis Ibanez	Status	resolved => closed
2007-09-21 15:56	Luis Ibanez	Note Added: 0009206
2007-09-24 02:44	Nicolas Savoire	Status	closed => feedback
2007-09-24 02:44	Nicolas Savoire	Resolution	fixed => reopened
2007-09-24 02:44	Nicolas Savoire	Note Added: 0009269
2007-09-24 08:34	Luis Ibanez	Note Added: 0009270
2007-09-24 08:46	Nicolas Savoire	Note Added: 0009272
2007-09-24 20:35	dbrussak	Note Added: 0009279
2007-10-08 08:21	Luis Ibanez	Relationship added	related to 0005692
2007-12-16 09:34	Aviv Hurvitz	Note Added: 0009942
2009-06-16 04:04	Tom Vercauteren	Note Added: 0016727
2009-06-16 04:04	Tom Vercauteren	Relationship replaced	duplicate of 0005692
2009-06-16 04:04	Tom Vercauteren	Duplicate ID	0 => 5692
2009-06-16 04:04	Tom Vercauteren	Status	feedback => resolved
2009-06-16 04:04	Tom Vercauteren	Resolution	reopened => fixed
2009-07-07 05:23	Tom Vercauteren	Status	resolved => closed

Notes
(0002729) Nicolas Savoire (reporter) 2005-07-27 11:55	A similarly incorrect function is used in BilateralImageFilter.

(0002771) Lydia Ng (developer) 2005-08-05 15:02	Triaged 8/5/05 tcon. Jim will look at it.

(0002810) Nicolas Savoire (reporter) 2005-08-08 03:45	I forgot to add that the bug affects only the i386 platform. Replacing unsigned int hilo[2]; ((double )hilo) = x + 103079215104.0; return (int)((hilo[1]<<16)\|(hilo[0]>>16)); by union { unsigned int hilo[2]; double d; } u; u.d = x + 103079215104.0; return (int)((u.hilo[1]<<16)\|(u.hilo[0]>>16)); seems to fix the problem.

(0002917) Nicolas Savoire (reporter) 2005-08-18 08:12	Another thing: the bug seems to appear only with gcc 3.4 and flag -O2/-O3

(0003285) user390 2005-11-20 23:34	Committed to the Main branch. http://www.itk.org/cgi-bin/viewcvs.cgi/Code/Common/itkBSplineInterpolationWeightFunction.txx?rev=1.10&root=Insight&view=log [^] http://www.itk.org/cgi-bin/viewcvs.cgi/Code/BasicFilters/itkBilateralImageFilter.txx?rev=1.25&root=Insight&view=log [^]

(0009206) Luis Ibanez (manager) 2007-09-21 15:56	correcting 'resolved' status to 'closed'

(0009269) Nicolas Savoire (reporter) 2007-09-24 02:44	It seems that the returned value is wrong for inputs really close but smaller than a given integer: BSplineFloor(0.999999) returns 1 instead of 0. I think the issue is the addition of a large constant value which triggers a loss of precision.

(0009270) Luis Ibanez (manager) 2007-09-24 08:34	Nicos, Can you please provide a test program that illustrate this problem ? We need a test in order to track the issue. Also, is it still the case that this error is only observed in gcc 3.4 with flags -O2/O3 ? For the long term maintenance it seems better to use vcl_floor(). Whether this is a performance issue or not, is something that we have to verify by using a profiling test. Are we saving 1% of the computation time ? is it 10% ? is it 50% ? For 1% of computation time savings, we probably shouldn't sacrifice portability and maintainability... Thanks Luis

(0009272) Nicolas Savoire (reporter) 2007-09-24 08:46	The problem appears on i386 (linux and OSX) platform. Compilation flags and gcc version have no influence on the result. Test case: #include <assert.h> inline int BSplineFloor(double x) { #if defined mips \|\| defined sparc \|\| defined __ppc__ return (int)((unsigned int)(x + 2147483648.0) - 2147483648U); #elif defined i386 \|\| defined _M_IX86 union { unsigned int hilo[2]; double d; } u; u.d = x + 103079215104.0; return (int)((u.hilo[1]<<16)\|(u.hilo[0]>>16)); #else return int(floor(x)); #endif } int main () { assert(BSplineFloor(0.999999)==0); }

(0009279) dbrussak (reporter) 2007-09-24 20:35	Hi, This particular BSplineFloor() precision issue is the one I reported here in issue 0005692: http://www.itk.org/Bug/view.php?id=5692 [^] I'm happy to delete/rescind my issue and just wait for resolution of this one. Is that possible? Daniel R.

(0009942) Aviv Hurvitz (reporter) 2007-12-16 09:34	I've also run into the problem where BSplineFloor returns an incorrect result. This caused BSPlineDeformableTransform ::TransformPoint() to SEGFAULT. Tested on a Pentium 4. Compiled with: - gcc version 4.1.3 20070629 (prerelease) (Debian 4.1.2-13) - MS Visual C++ 2005 A sample program: >>>>>>>>>>>>>>>>>>>>>>>>> #include <iostream> #include <math.h> inline int BSplineFloor(double x) { #if defined mips \|\| defined sparc \|\| defined __ppc__ return (int)((unsigned int)(x + 2147483648.0) - 2147483648U); #elif defined i386 \|\| defined _M_IX86 union { unsigned int hilo[2]; double d; } u; u.d = x + 103079215104.0; return (int)((u.hilo[1]<<16)\|(u.hilo[0]>>16)); #else return int(floor(x)); #endif } int main(int argc, char argv[]) { double x = 16.999995; int bsplineFloorResult = BSplineFloor(x); int standardFloor = static_cast<int>( floor(x) ); std::cout.precision(20); std::cout << "x = " << x << std::endl; std::cout << "BSspline_floor = " << bsplineFloorResult << std::endl; std::cout << "Standard floor = " << standardFloor << std::endl; / Prints: x = 16.999994999999998413 BSspline_floor = 17 Standard floor = 16 */ return 0; } >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

(0016727) Tom Vercauteren (developer) 2009-06-16 04:04	Fixed in HEAD: http://www.itk.org/cgi-bin/viewcvs.cgi/Code/Common/itkBSplineInterpolationWeightFunction.txx?root=Insight&r1=1.15&r2=1.16&sortby=date [^] http://www.itk.org/cgi-bin/viewcvs.cgi/Code/BasicFilters/itkBilateralImageFilter.txx?root=Insight&r1=1.29&r2=1.30&sortby=date [^]