CGAL users discussion list

[Benjamin Kehlet <>]

Hello!

This seems to be a bug in include/CGAL/Snap_rounding_traits_2.h:67-68
in Snap_2::operator(). The "pixel number" is computed as floor(x /
pixel_size). This floor operation is carried out in double precision,
which is not safe.

I tried to replace the call to std::floor with the simplest possible
implementation and it solves the problem for my case.

My very silly floor() implementation essentially just counts the
number of integers smaller than x). Here is the diff:

diff --git a/3rdparty/CGAL-4.4/include/CGAL/Snap_rounding_traits_2.h
b/3rdparty/CGAL-4.4/include/CGAL/Snap_rounding_traits_2.h
index 2c252e5..5538e0d 100644
--- a/3rdparty/CGAL-4.4/include/CGAL/Snap_rounding_traits_2.h
+++ b/3rdparty/CGAL-4.4/include/CGAL/Snap_rounding_traits_2.h
@@ -26,6 +26,36 @@

 #include <CGAL/Arr_segment_traits_2.h>

+
+namespace
+{
+template<typename NT>
+int floor_for_testing_only(NT x)
+{
+  const bool negative = x < 0;
+  if (negative)
+    x = -x;
+
+  int f = 0;
+  while (x >= 1)
+  {
+    f++;
+    x -= 1;
+  }
+
+  if (negative)
+  {
+    if (x != 0)
+      f++;
+
+    f = -f;
+  }
+
+  return f;
+}
+}
+
 namespace CGAL {

 template<class Base_kernel>
@@ -64,10 +94,8 @@ public:
       NT x_tmp = p.x() / pixel_size;
       NT y_tmp = p.y() / pixel_size;

-      double x_floor = std::floor(CGAL::to_double(x_tmp));
-      double y_floor = std::floor(CGAL::to_double(y_tmp));
-      x = NT(x_floor) * pixel_size + pixel_size / NT(2.0);
-      y = NT(y_floor) * pixel_size + pixel_size / NT(2.0);
+      x = floor_for_testing_only(x_tmp) * pixel_size + pixel_size / NT(2.0);
+      y = floor_for_testing_only(y_tmp) * pixel_size + pixel_size / NT(2.0);
     }
   };

I don't know how this should be implemented in CGAL in a both
efficient and number type agnostic way. (The only way I could think of
was to use std::floor() with double precision and then check and
adjust the result afterwards, but that sounds suboptimal), but I do
hope that a fix will find its way into CGAL. Let me know if there is
anything I can do in this regard.

Best regards

Benjamin Kehlet


2014-06-27 22:11 GMT+02:00 Benjamin Kehlet 
<>:
> Hello again!
>
> Here is a simpler program that (if I understand the snap rounding
> package correctly) reveals a bug
>
> ---
> #include <CGAL/Cartesian.h>
> #include <CGAL/Quotient.h>
> #include <CGAL/MP_Float.h>
> #include <CGAL/Snap_rounding_traits_2.h>
> #include <CGAL/Snap_rounding_2.h>
>
> typedef CGAL::Quotient<CGAL::MP_Float>           Number_type;
> typedef CGAL::Cartesian<Number_type>             Kernel;
> typedef CGAL::Snap_rounding_traits_2<Kernel>     Traits;
> typedef Kernel::Segment_2                        Segment_2;
> typedef Kernel::Point_2                          Point_2;
>
> int main(int argc, char** argv)
> {
>   Point_2 p1(1, 0);
>   Point_2 p2(0, 0);
>   Number_type pixel_width(0.1);
>
>   CGAL::Hot_pixel<Traits> hp(p1, pixel_width);
>   CGAL::Segment_data<Traits> sd(p1, p2);
>
>   CGAL::SEG_Direction dir;
>   sd.determine_direction(dir);
>
>   if (!hp.intersect(sd, dir))
>     std::cout << "Error! Segment did not intersect hot pixel of
> endpoint" << std::endl;
> }
> ---
>
> This hits the bug, but only with a debug build. However, my guess is
> that this is not related to code that is removed when doing a release
> build. I think it is more likely that something in the creation of the
> hot pixel or in the intersection test is somehow prone to roundoff
> errors.
>
> I will try to investigate more. I would be nice if someone who know
> the package better than me could confirm that I'm not doing something
> completely wrong.
>
> Best regards
>
> Benjamin Kehlet
>
> 2014-06-27 1:58 GMT+02:00 Benjamin Kehlet 
> <>:
>> Hello!
>>
>> I am attempting to use the 2D snap rounding package to safely get from
>> exact numbers to float in our mesh generation tool in FEniCS. However,
>> I get some (to me) unexpected results. Here is a simple example:
>> When running this simple modification of
>> examples/Snap_rounding_2/snap_rounding.cpp:
>> ---
>> #include <CGAL/Cartesian.h>
>> #include <CGAL/Quotient.h>
>> #include <CGAL/MP_Float.h>
>> #include <CGAL/Snap_rounding_traits_2.h>
>> #include <CGAL/Snap_rounding_2.h>
>>
>> typedef CGAL::Quotient<CGAL::MP_Float>           Number_type;
>> typedef CGAL::Cartesian<Number_type>             Kernel;
>> typedef CGAL::Snap_rounding_traits_2<Kernel>     Traits;
>> typedef Kernel::Segment_2                        Segment_2;
>> typedef Kernel::Point_2                          Point_2;
>> typedef std::list<Segment_2>                     Segment_list_2;
>> typedef std::list<Point_2>                       Polyline_2;
>> typedef std::list<Polyline_2>                    Polyline_list_2;
>>
>> int main()
>> {
>>   Segment_list_2 seg_list;
>>   Polyline_list_2 output_list;
>>
>>   seg_list.push_back(Segment_2(Point_2(1, 0), Point_2(0, 0)));
>>
>>   // Generate an iterated snap-rounding representation, where the centers 
>> of
>>   // the hot pixels bear their original coordinates, using 5 kd trees:
>>   
>> CGAL::snap_rounding_2<Traits,Segment_list_2::const_iterator,Polyline_list_2>
>>     (seg_list.begin(), seg_list.end(), output_list, .1, true, false, 5);
>>
>>   int counter = 0;
>>   Polyline_list_2::const_iterator iter1;
>>   for (iter1 = output_list.begin(); iter1 != output_list.end(); ++iter1) {
>>     std::cout << "Polyline number " << ++counter << ":\n";
>>     Polyline_2::const_iterator iter2;
>>     for (iter2 = iter1->begin(); iter2 != iter1->end(); ++iter2)
>>       std::cout << "    (" << iter2->x() << ":" << iter2->y() << ")\n";
>>   }
>>
>>   return(0);
>> }
>> ---
>>
>> I get this result:
>> benjamik@benjamik-ThinkPad-X300:~/software/CGAL-4.4/examples/Snap_rounding_2/build$
>> ./snap_rounding
>> Polyline number 1:
>>     (0.1/2:0.1/2)
>> benjamik@benjamik-ThinkPad-X300:~/software/CGAL-4.4/examples/Snap_rounding_2/build$
>>
>> So the one segment in the input seems to be "reduced" to a single
>> point. Is this expected? If the pixel size argument in the call to
>> CGAL::snap_rounding_2() is changed to "Number_type(1, 10)", then I get
>> the output that I would expect:
>>
>> benjamik@benjamik-ThinkPad-X300:~/software/CGAL-4.4/examples/Snap_rounding_2/build$
>> ./snap_rounding
>> Polyline number 1:
>>     (210/200:10/200)
>>     (10/200:10/200)
>> benjamik@benjamik-ThinkPad-X300:~/software/CGAL-4.4/examples/Snap_rounding_2/build$
>>
>> This may be due to my limited knowledge of the
>> Quotient<CGAL::MP_float> number type, but from the reference manual it
>> seems that it should be constructable from a double.
>>
>> This occurs with a debug build of CGAL 4.4. I get the equal results on
>> a 32-bit Ubuntu 13.10 box with gcc 4.8.1 and with a 64-bit debian
>> Jessie box with gcc 4.9.0.
>>
>> Best regards
>>
>> Benjamin Kehlet