Re: [boost] Geometry and spatial indexes, my opinion

Subject: Re: [boost] Geometry and spatial indexes, my opinion
From: Simonson, Lucanus J (lucanus.j.simonson_at_[hidden])
Date: 2008-10-09 17:33:20

Brandon wrote:
>Well, the jury is still out as to how effective the *real* solution
will >be. 
>I start by adding the ability to perform type dispatch based on the
number 
>type of the coordinate. I manage that by again using a traits class
(for >the 
>coordinates):

Interesting, I have coordinate_traits and register the coordinate type
for tag dispatching purposes.  My coordinate traits register related
data types such as area data type.  For 32 bit integer coordinates I use
a 64 bit integer to compute and return area values, and I register it
through the 32 bit integer coordinate traits.

>Then I can do a dispatch on the traits. The only issue I've tackled so
far 
>is for comparing segments on a sweepline. I did this by using 
>boost::rational ( so overflow may still be a problem in some cases ). 

I use mpq_class because some cases is way to many cases for me.  I
compute the y as a mpq_class at a given x, which allows me exact,
overflow proof comparison, but I don't store the mpq_class value in the
sweepline data structure.  I compute it on the fly as needed.

>The 
>other issue I have yet to tackle is intersection between two segments
for 
>integers. I have an implementation that uses O'Rourke's method from 
>Computational Geometry in C. In his method the segments are first
converted 
>into parametric coordinates and then compared to see if the resulting
>system 
>of equations results in a point which is contained in each segment. As
you 
>well know you can't do this with segments in integral coordinates as
the 
>intersection point often has fractional coordinates. I suppose the only

>thing you can do with that is resort to rationals. 

Which is what I'm doing.  I simplified the algebra so that there is only
one equation for x and y respectively.

x = (x11 * dy1 * dx2 - x21 * dy2 * dx1 + y21 * dx1 * dx2 - y11 * dx1 *
dx2) / (dy1 * dx2 - dy2 * dx1);

transpose x and y and it gives you the y value as well.  When evaluated
with gmp rational it gives the exact intersection.

I also have a predicate that determines whether two line segments
intersect, but does not compute the intersection point.  I do that by
comparing slope cross products looking for both points of one segment
being on the same side of the other line segment.  It does not need more
than 64 bit integer to deal with 32 bit coordinate segments.

>This will be necessary as 
>you must track all event points during the sweep (including the 
>intersections), and reporting of intersections requires that all
segments 
>intersecting in the same point be associated with that point in the 
>sweepline. This leads me to think that the easiest solution is to make
all 
>the coordinates rational in the event structure (for sorting 
>lexicographically) and then when actually reporting, making the
rationals 
>available to the visitor class (which is how the reporting is done)
along 
>with the offending segments. This way users who are able to use the 
>fractional intersections can still use them.. and those who cannot
still 
>have the segments which intersect (and I assume they can glean whatever

>information they need from them.) In cases like yours where you
constraint 
>the slopes to be 0,1 or infinity, you can no doubt do a lot better on
the 
>intersection, but a general rational solution for intersection would 
>probably still work. I still need to think about it more to see how I
can 
>mitigate issues with overflow.

Don't forget -1.  Actually, I've been extending the library to general
polygons recently.

You have pretty much no other option than to use gmp to avoid overflow,
as far as I can tell.  In general you need 96 bits of integer precision
to compute the numerator of the x value in my expression above.  Other
than not including the gmp header file in boost licensed source, it is
no problem to depend on it because it is open source and more than
portable in that it explicitly targets almost every platform.

My arbitrary line segment intersection algorithm snaps the intersection
points to the integer grid and deals with the problems that causes when
line segments move and collide with nearby verticies.  This would be
quite a bit harder to do when snapping to the floating point grid, I
would imagine.

The output of my integer line intersection algorithm is the input
segments broken into smaller segments at the intersection points and
snapped to the integer grid.  The invariant it guarantees is that no
pair of output segments intersect at a non-vertex position.  Once you
have that you can scan it quite simply and build a graph or do Booleans
or whatever else you need to.

I'm hoping to get the algorithms completed some time this quarter and
released to boost.  The intention is to have 100% robust integer
linescan and boolean operations.  I really haven't thought about how to
accomplish the same with floating point coordinates.

Luke