PageRenderTime 16ms CodeModel.GetById 8ms app.highlight 3ms RepoModel.GetById 1ms app.codeStats 1ms

C++ Header | 214 lines | 3 code | 5 blank | 206 comment | 0 complexity | 04720261abd659af4e1a8e8419d16a56 MD5 | raw file
  1// Boost.Geometry (aka GGL, Generic Geometry Library)
  3// Copyright (c) 2007-2011 Barend Gehrels, Amsterdam, the Netherlands.
  4// Use, modification and distribution is subject to the Boost Software License,
  5// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  8#ifndef _DOXYGEN_SETS_HPP
  9#define _DOXYGEN_SETS_HPP
 15\page sets Spatial set-theoretic operations (union, intersection, difference)
 20\section sets_par1 Introduction
 22The GGL implementation of the algorithms intersection, union, difference and symmetric difference is based on set theory (a.o.
 23This theory is applied for spatial sets.
 25Intersection and union are so-called set-theoretic operations. Those operations
 26work on sets, and geometries (especially
 27polygons and multi-polygons) can be seen as sets, sets of points.
 30The first section will repeat a small, relevant part of the algebra of sets, also describing the notation used in this page. The next section will extend this algebra of sets for spatial sets (polygons).
 32\section sets_par2 Algebra of sets in a spatial context
 34- A ∩ B : the intersection of two sets A and B is the set that contains all elements of A that also belong to B (aka AND)
 35- A ∪ B : the union of two sets A and B is the set that contains all elements that belong to A or B (aka OR)
 36- A<small><sup>c</sup></small> :    the complement of set A is the set of elements which do not belong to A
 37- A \ B :    the difference of two sets A and B is the set of elements which belong to A but not to B
 38- A &#8710; B : the symmetric difference of two sets A and B is the set of elements which belong to either A or to B, but not to A and B (aka XOR)
 40(Source of most definitions:
 42There are several laws on sets and we will not discuss them all here. The most important for this page are:
 43- B \ A = A<small><sup>c</sup></small> &#8745; B and, vice versa, A \ B = B<small><sup>c</sup></small> &#8745; A
 44- A &#8710; B = (B \ A) &#8746; (A \ B)  (
 47\section sets_par3 Polygons
 50Polygons are sets of points, and, therefore polygons follow all definitions and laws for sets. For pragmatic reasons and implementations in computer programs, polygons have an orientation, clockwise or counter clockwise. Orientation is not part of most set theory descriptions, but is an important aspect for the appliance of sets to polygon operations.
 52If a polygon is (arbitrarily) defined as having its vertices in clockwise direction:
 53- then its interior lies on the right side of the edges []
 54- and its exterior lies, therefore, on the left side of the edges
 56This definition is important for the spatial interpretation sets.
 58- If set A describes the interior of a polygon, then A<small><sup>c</sup></small>, its complement, describes the exterior of a polygon.
 59- Stated differently, set A is a polygon, all points belonging to A are inside the polygon. Its complement, A<small><sup>c</sup></small>, contains all points not belonging to A.
 60- If A is a polygon with its vertices oriented clockwise, A<small><sup>c</sup></small> is a polygon with the same vertices, but in reverse order, so counter clockwise. Both sets have their points belonging to them at the right side of their edges
 62\image html set_a_ac.png
 64The last observation is helpful in calculating the difference and the symmetric difference:
 65- the difference B \ A is defined above by the law B \ A = A<small><sup>c</sup></small> &#8745; B.
 66    In polygon terms it is therefore the intersection of the "reverse of A" with B.
 67    To calculate it, it is enough to have polygon A input in reverse order, and intersect this with polygon B
 68- the symmetric difference A &#8710; B is defined above by the law (B \ A) &#8746; (A \ B), which leads to
 69    (A<small><sup>c</sup></small> &#8745; B) &#8746; (B<small><sup>c</sup></small> &#8745; A).
 70    To calculate the symmetric difference, it is enough to have polygon A input in reverse order,
 71    intersect this with polygon B, store the result; then have polygon B input in reverse order,
 72    intersect this with polygon A, add this to the result and this is the symmetric difference.
 73    The combination of both sub-results does not have to be intersected as it is only touching
 74    on vertices and do not have overlaps.
 77\section sets_par4 Implementation of intersection and union
 78All spatial set-theoretic operations are implemented in shared code. There is hardly any difference in code
 79between the calculation of an intersection or a union. The only difference is that at each intersection point,
 80for an intersection the right turn should be taken. For union the left turn should be taken.
 82\image html set_int_right_union_left.png
 85This is implemented as such in GGL. The turn to be taken is a variable.
 87There is an alternative to calculate union as well:
 88- the union A &#8746; B equals to the complement of the intersection of the complements of the inputs,
 89    (A<small><sup>c</sup></small> &#8745; B<small><sup>c</sup></small>)<small><sup>c</sup></small> (De Morgan law,
 92There is an additional difference in the handling of disjoint holes (holes which are not intersected). This is also
 93implemented in the same generic way (the implementation will still be tweaked a little to have it even more generic).
 95For a counter clockwise polygon, the behaviour is the reverse: for intersection take the left path, for union
 96take the right path. This is a trivial thing to implement, but it still has to be done (as the orientation was introduced
 97in a later phase in GGL).
100\section sets_par5 Iterating forward or reverse
101As explained above, for a difference, the vertices of the first polygon should be iterated by a forward iterator, but
102the vertices of the second polygon should be iterated by a reverse iterator (or vice versa). This (trivial) implementation
103still has to be done. It will <b>not</b> be implemented by creating a copy, reversing it, and presenting it as input to the
104set operation (as outlined above). That is easy and will work but has a performance penalty. Instead a reversible iterator will used,
105extended from Boost.Range iterators, and decorating a Boost.Range iterator at the same time, which can travel forward or backward.
107It is currently named \b reversible_view and usage looks like:
111template <int Direction, typename Range>
112void walk(Range const & range)
114    typedef reversible_view<Range, Direction> view_type;
115    view_type view(range);
117    typename boost::range_const_iterator<view_type>::type it;
118    for (it = boost::begin(view); it != boost::end(view); ++it)
119    {
120        // do something
121    }
124walk<1>(range); // forward
125walk<-1>(range); // backward
130\section sets_par6 Characteristics of the intersection algorithm
131The algorithm is a modern variant of the graph traversal algorithm, after Weiler-Atherton
134It has the following characteristics (part of these points are deviations from Weiler-Atherton)
135- No copy is necessary (the original Weiler-Atherton, and its descendants, insert intersections in (linked) lists,
136    which require to make first copies of both input geometries).
137- All turning points (which are intersection points with side/turn information) are unaware of the context, so we have (and need) no information about
138    if, before the intersection point, a segment was inside or outside of the other geometry
139- It can do intersections, unions, symmetric differences, differences
140- It can handle polygons with holes, non-convex polygons, polygons-with-polygons, polygons-with-boxes (clip), rings, multi-polygons
141- It can handle one polygon at the time (resolve self-intersections), two polygons (the normal use case), or more polygons (applicable
142    for intersections and unions)
143- It can handle clockwise and counter-clockwise geometries
146\section sets_par7 Outline of the intersection algorithm
147The actual implementation consists of the next phases.
149\b 1 the input geometries are indexed (if necessary). Currently we use monotonic sections for the index. It is done
150by the algorithm \b sectionalize. Sections are created is done on the fly, so no preparation is required before (though this would improve
151performance - it is possible that there will be an alternative variant where prepared sections or other indexes are part of the input).
152For box-polygon this phase is not necessary and skipped. Sectionalizing is done in linear time.
154\b 2, intersection points are calculated. Segments of polygon A are compared with segments of polygon B. Segment intersection is only
155done for segments in overlapping sections. Intersection points are not inserted into the original polygon or in a copy. A linked list is therefore not necessary.
156This phase is called \b get_intersection_points. This function can currently be used for one or two input geometries, for self-intersection
157or for intersection.
158Because found intersections are provided with intersection-information, including a reference to their source, it is possible (but currently not
159implemented) to have more than two geometry inputs.
161The complexity of getting the intersections is (much) less than quadratic (n*m) because of the monotonic sections. The exact complexity depends on
162the number of sections, of how the input polygons look like. In a worst case scenario, there are only two monotonic sections per polygon and both
163overlap. The complexity is then quadratic. However, the sectionalize algorithm has a maximum number of segments per section, so for large polygons
164there are always more monotonic sections and in those cases they do not overlap by the definition of "monotonic".
165For boxes, the complexity is linear time.
167To give another idea of how sections and indexes work:
168For a test processing 3918 polygons (but not processing those of which envelopes do not overlap):
169- brute force (O(n<small><sup>2</sup></small>)): 11856331 combinations
170- monotonic sections: 213732  combinations (55 times less)
171- a spatial index: 34877 combinations (1/6 of using monotonic sections)
172So there can still be gained some performance by another index. However the current spatial index implementation (in an extension,
173not in Formal Review) will still be revisited, so it is currently not used.
175<i>In "normal" cases 84% of the time is spent on finding intersection points. These divisions in %'s refers to the performance test described elsewhere</i>
177One piece of information per intersection points is if it is \b trivial. It is trivial if the intersection is not located at segment end points.
179\b 3, the found intersection points are merged (\b merge_intersection_points), and some intersections can be deleted (e.g. in case of
180collinearities). This merge process consists of sorting the intersection points in X (major) and Y (minor) direction, and merging intersections with a common location together. Intersections with common
181locations do occur as soon as segments are collinear or meet at their end points.
182This phase is skipped if all intersection points are trivial.
184<i>About 6% is spent on merging.</i>
186\b 4, some turns need to be adapted. If segments intersect in their interiors, this is never necessary. However, if segments intersect on their
187end points, it is sometimes necessary to change "side" information to "turn" information. This phase is called \b adapt_turns.
189The image below gives one example when adapting turns is necessary. There is side information, both segments have sides \b left and \b right, there
190is also \b collinear.
191However, for an intersection no turn should be taken at all, so no right turn. For a union, both polygons have to be travelled.
192In this case the side information is adapted to turn information, both turns will be \b left. This phase is skipped if all intersection points are trivial.
194\image html set_adapt_turns.png
197\b 5, the merged intersection points are enriched (\b enrich_intersection_points) with information about a.o. the next intersection point (travel information).
199<i>About 3% is spent on enrichment.</i>
201\b 6, polygons are traversed (\b traverse) using the intersection points, enriched with travel information. The input polygons are traversed and
202at all intersection poitns a direction is taken, left for union, right for intersection point (for counter clockwise polygons this is the other way
203round). In some cases separate rings are produced. In some cases new holes are formed.
205<i>About 6% is spent on traversal.</i>
207\b 7, the created rings are assembled (\b assemble) into polygon(s) with exterior rings and interior rings. Even if there are no intersection points found, this process can be important to find containment and coverage.
209<i>Timing of this phase is not yet available, as the comparison program work on rings.</i>
214#endif // _DOXYGEN_SETS_HPP