CGAL users discussion list

[Benjamin Kehlet <>]

2014-07-02 11:11 GMT+02:00 Andreas Fabri 
<>:
> Hi Benjamin,
>
> thank you for your bug report and even the analysis.
> We will add a floor function for the CGAL number types
> in order to avoid passing to floating point arithmetic.

Sounds good! Thanks!

>
> Can you work around the problem for the time being?

Yes, it is not beatiful, but it seems to work for now :)

> Note that even without the bug, segments may end up
> in the same pixel.

Yes, I'm aware of that. Thanks.

Regards

Benjamin

>
> Best,
>
> Andreas
>
>
> On 28/06/2014 20:23, Benjamin Kehlet wrote:
>>
>>
>> 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
>>
>>
>
> --
> Andreas Fabri, PhD
> Chief Officer, GeometryFactory
> Editor, The CGAL Project
>
> phone: +33.492.954.912    skype: andreas.fabri
>
> --
> You are currently subscribed to cgal-discuss.
> To unsubscribe or access the archives, go to
> https://sympa.inria.fr/sympa/info/cgal-discuss
>
>