Re: [boost] GGL Review

Subject: Re: [boost] GGL Review
From: Brandon Kohn (blkohn_at_[hidden])
Date: 2009-11-16 13:18:24

Barend Gehrels wrote:
>> There are many buggy floating point boolean op libraries out there. 
>> Why do we need another? 
> Can you explain this? I've done several benchmarks comparing 7 libraries 
> and I think they all reported the correct results. Do you refer to cgal, 
> geos, gpc, terralib, gtl (b.p.), wykobi or ggl? Or do you use them in 
> other domains than we've tested them?

Hi Barend,

I'm referring to OpenGL's CSG boolean ops and PolyBoolean, both of which 
routinely fail in an application with which I work. (It requires 
partitioning space by iteratively subtracting isovists from a triangle 
partition of a polygon with holes.) I have read the CGAL also has issues 
for similar reasons.

>> And are we talking about numerical floating point issues, or
>> algorithmic bugs?  I don't consider those to be the same thing.
> 
> I totally agree with this. Algorithmic bugs should not be there, not in 
> any library, of course. Numerical floating point issues can be addressed 
> using arbitrary arithmetic types, and besides that other measures can be 
> taken using the default floating point types. Only the very last thing 
> is not done in GGL.

I think this goes without saying. I was of course talking about issues 
where an algorithm fails because the underlying numeric type cannot be 
used correctly in the algorithm.

As for GGL's working with arbitrary arithmetic types, I would have to 
take your word for it. I have found a few bits in the code that seem to 
be converting to double.. and so I'm not so sure about your claim.

//excerpt from relate_cartesian_segments::relate
double ra = double(da) / double(d);
double rb = double(db) / double(d);

Also when I perform:

test_all<ggl::point_xy<int> >(); //to test correct

It uses the following predicate (presumably for sorting coordinates):

std::less<double> predicate

...

//From correct_polygon:
correct_ring<ring_type, std::greater<double> >::apply(*it);

//* From get_intersection_points.hpp
double eps = 1.0e-10;
if (dir.how == 'i'
     && (dir.ra < eps
     || dir.rb < eps
     || 1.0 - dir.ra < eps
     || 1.0 - dir.rb <  eps

So these are at least four places where double/ fp issues could creep 
into an calculations with arbitrary precision types.

Regards,

Brandon