CGAL users discussion list

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

Can you work around the problem for the time being?
Note that even without the bug, segments may end up
in the same pixel.

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