Re: [boost] [Review] GGL review starts today, November 5th

Subject: Re: [boost] [Review] GGL review starts today, November 5th
From: Barend Gehrels (barend_at_[hidden])
Date: 2009-11-12 14:39:52

Hi Luke,

Herewith answers on questions regarding FP-precision.

> 7. When intersection points are snapped to the floating point grid at the output of the polygon pair intersection algorithm edges may move laterally some small distance which may introduce self intersection artifacts.  How is this handled by the polygon intersection algorithm?
>   

This is not handled explicitly. On segment-segment intersection in the 
epsilon-regions, the intersection will fall between the endpoints of 
both segments. Even if it is rounded, it will not surpass those 
endpoints. Therefore, self-intersections will be unlikely in most cases. 
I'll not say that it is impossible, but it is unlikely.

Our implementation is not yet checked on 100% numerical robustness. If 
using floating point (IEEE-single), there might be errors in the 1e-3x 
regions. If using double precision, there will be less errors. That is 
to say: in close-to-100% of the use-cases: no errors. But in the regions 
of 1e-3xx differences, errors might be expected, and in some cases they 
might even cause self-intersections, I'll not deny that now, we've to do 
more research here.

But you can go higher then, if you want that precision, use long double 
(on MSVC). Or use the high-precision arithmetic support that we're 
providing, which is slower but more exact. It's explained below in some 
more detail.

> 11. When a polygon is intersected against a bounding box the same problem can happen that self intersections are introduced in the output polygon(s) when intersection points are snapped to the floating point grid.  These artifacts cannot be addressed by a linear time algorithm.  Are the output polygons of the polygon intersects bounding box algorithm implemented by this library potentially self intersecting?
>   

I cannot say much more here than I answered on your question 7. In other 
words: if you really want IEE-single, you can expect this effect in rare 
cases. Using IEEE-double, it will be very unlikely. If you're using 
coordinates with differences in the regions of 1e-3XX, I would advice a 
high-precision number data type.

The additional mechanism I already mentioned in previous answer, on the 
integer/floating point issue: even if you're using a point type based 
on  IEEE-single coordinates, you can still do the calculations in 
IEEE-double or long double. Or high precision.

In all my career (18 years GIS) I didn't encounter problems, using 
IEEE-double (using single, I did). Therefore, and because of recent 
tests, I mentioned repeatedly that it is rare or unlikely. But it is not 
impossible, and therefore we're willing to support users who want to go 
beyond certain limits. We decided to go first for high precision 
arithmetic numbers. See below the answer on your next question. That is 
exact (or quite exact) but slow. We're also willing to tweak things to 
be more robust in the (long) double regions, provided that it'll not 
affect performance. But as I'm not an expert on numerical robustness 
(did read the introductions though), we'll probably need support for 
that and Fernando did offer his help here on this list.

> 17. Can you tell us more about the support and usage of high precision numerical datatypes?
>   

Sure. We worked on a wrapper for GMP and CLN and tested with that. That 
wrapper (numeric_adaptor) might be redundant because there is also 
support in Boost.Math bindings. We were somehow not aware of that but, 
however, it was a useful exercise because it works in a similar way, and 
we could adapt our library to using high-precision scenario. Not all 
algorithms are adapted, currently, but we did some to research and to 
prove that it works well.

The basic example is in x02 (x02__numeric__adaptor__example.cpp), on the 
website and in the example folder. You can use a point which is based on 
GMP or CLN (or other arbitrary precision types) coordinate types. All 
coordinates are stored in that type then, and calculations are done in 
that type. Besides that, you can store a point in double or long double, 
and do calculations in high-precision. Like a mentioned above and in 
previous answer.

Hope this makes this more clear.

Regards, Barend