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

Subject: Re: [boost] [Review] GGL review starts today, November 5th
From: Simonson, Lucanus J (lucanus.j.simonson_at_[hidden])
Date: 2009-11-12 18:58:20

Barend Gehrels wrote:
> Hi Luke,
> 
> 
> Simonson, Lucanus J wrote:
>> 
>> I still don't know what the find_intersections algorithm's
>> complexity is or how it works.  What data structure does it
>> currently use for indexing?  Based on the documentation I'm guessing
>> you use a 1D spatial index for red blue overlap detection between
>> monotonic sections.  Is that correct?    
> 
> As in the documentation: monotonic sections, and the algorithm is
> named "get_intersection_points". How it works: I worked quite a day
> on that 
> page and I hope it is clear enough for the detail needed for a review.
> If not, I'm sorry to hear that but I don't know if I can still improve
> this before the review period ends.

I have to read the code to learn the algorithm...

>From get_intersection_points.hpp
00390         for (typename boost::range_const_iterator<sections1_type>::type
00391                     it1 = sec1.begin();
00392             it1 != sec1.end();
00393             ++it1)
00394         {
00395             for (typename boost::range_const_iterator<sections2_type>::type
00396                         it2 = sec2.begin();
00397                 it2 != sec2.end();
00398                 ++it2)
00399             {
00400                 if (! ggl::detail::disjoint::disjoint_box_box(
00401                                 it1->bounding_box, it2->bounding_box))
00402                 {
00403                     get_ips_in_sections
00404                     <
00405                         Geometry1,
00406                         Geometry2,
00407                         typename boost::range_value<sections1_type>::type,
00408                         typename boost::range_value<sections2_type>::type,
00409                         IntersectionPoints
00410                     >::apply(
00411                             source_id1, geometry1, *it1,
00412                             source_id2, geometry2, *it2,
00413                             false,
00414                             intersection_points, trivial);
00415                 }
00416             }
00417         }

Your algorithm is: for all red montonic sections, inspect all blue monotonic sections for bounding box overlap and if overlap exists inspect segments of monotonic section pair for line intersection and report said intersections.  It has O(n^2) runtime complexity wrt. monotonic sections.  Monotonic sections are in the worst case equivilent to the number of edges.

>> [...] Do you realize that pathological cases of large number of 45
>> degree edges with no intersection between them with result in O(n^2)
>> work even with a good 2d spatial index based algorithm for finding
>> overlaps because 2d spatial indexing works on bounding boxes. [...] 
>> 
> A large number of edges will result in multiple monotonic sections, as
> described in the new doc-page, so also then it will not be quadratic,
> in 
> most cases. You're referring (I think) to a very large number of
> diagonal parallel edges and indeed all those bounding boxes will
> overlap. And I understand most spatial indexes will not help here (did
> you scan all variants?). As this seems to be problematic for both our
> libraries we can discuss this further (but at this moment it is not so
> convenient for me).

My library is not relevant to the discussion.  I don't understand what you mean my not convenient.  Do you mean it is not convenient to discuss the algorithmic complexity and worst case inputs for your algorithm during the review, or just that you are preoccupied with other matters?
 
> Anyway, I think our figures are clear enough, there is no library
> found 
> which even comes close, within several factors. So I would not worry
> too 
> much about its current performance. If you find a real test-case where
> it is too slow, I'm interested, we're happy to (try to) improve it.

My concern about performance is that your benchmarking may not expose performance problems.  Your claims about performance are based upon one benchmark (with sever minor variations of one of the two input geometries) from which you generalize about performance.  The algorithm has a fatal weakness to 

>> Holes are obviously associated to their enclosing outer rings by
>> your algorithm because that is the semantics of the output data
>> structure.  My question is, how is it identified by your algorithm
>> that a hole is enclosed by a polygon?  This is something that is
>> prone to pathological performance in the cases where the number of
>> holes is very large, so it is an important part of the algorithm.   
>> 
> I agree, and it's currently sorted on descending area to get obvious
> possibilities. There are alternatives here such as bboxes, spatial
> indexes or multiple-point-in-polygon, there is room for research and
> possible improvement.

This is not a complete description of an algorithm.  Are you saying you sort outer rings and inner rings by descending area and then do all pairs point in polygon check between pair of outer and inner ring that could possibly have containment relationship by area (and I'd guess bounding box) criteria to find which outer shells contain which holes?  That algorithm is O(n^2) wrt edges in the pathological case, which is exactly what I feared.  This is really the fundamental weakness of the Atherton-Wieler algorithm: that it forces handling of hole containment to post processing.

The algorithm is O(n^2) wrt. both monotonic sections and holes.  This means that for all inputs with large numbers of holes or low average length of monotonic sections but large total number of edges the algorithm with suffer alogorithmic (quadratic) slowness.  These weaknesses are not exposed in your benchmark testing.

>>> - interior_type is there a polygon may have holes. If library users
>>> use polygons without holes, they can take a ring (or, stated more
>>> precisely: a type fulfilling the Ring Concept)
>>> - interior_type is expected to be a Boost.Range (for a read-only
>>> polygon) and to be STL container-conformant (for a mutable polygon).
>>> 
>> 
>> I'm not too keen on this requirement.  It turns out that my own
>> polygon with holes has a private data member which is an stl
>> container that contains the holes, but I don't expose that in the
>> public interface.  I don't believe I could adapt my
>> polygon_with_holes_data to your polygon concept without changing its
>> implementation, however, I am certain that I could adapt your
>> polygon data type to my polygon_with_holes_concept.  I doubt there
>> are many legacy polygon types that can be successfully adapted to
>> your mutable polygon concept.  This calls into question to the
>> utility of the polygon concept.  I would rather see it done the same
>> way as linestring where you can either use stl style pushback and
>> clear or supply your own clear and append traits.           
>> 
> Also here we're happy to provide it, it is not satisfying your needs,
> especially if this enhances interoperability we're certainly willing
> to 
> adapt and finetune things. Actually I'm also glad to hear positive
> feedback about the way we designed it for linestring.

I think fixing the polygon traits to be like the linestring/ring traits would improve both consistency and genericity of the interfaces and should be required before the library is accepted into boost.

Regards,
Luke