1 /**************************************************************
2  *
3  * Licensed to the Apache Software Foundation (ASF) under one
4  * or more contributor license agreements.  See the NOTICE file
5  * distributed with this work for additional information
6  * regarding copyright ownership.  The ASF licenses this file
7  * to you under the Apache License, Version 2.0 (the
8  * "License"); you may not use this file except in compliance
9  * with the License.  You may obtain a copy of the License at
10  *
11  *   http://www.apache.org/licenses/LICENSE-2.0
12  *
13  * Unless required by applicable law or agreed to in writing,
14  * software distributed under the License is distributed on an
15  * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
16  * KIND, either express or implied.  See the License for the
17  * specific language governing permissions and limitations
18  * under the License.
19  *
20  *************************************************************/
21 
22 
23 
24 // MARKER(update_precomp.py): autogen include statement, do not remove
25 #include "precompiled_basegfx.hxx"
26 
27 #include <algorithm>
28 #include <functional>
29 #include <vector>
30 
31 #include "bezierclip.hxx"
32 
33 // -----------------------------------------------------------------------------
34 
35 /* Implements the theta function from Sedgewick: Algorithms in XXX, chapter 24 */
theta(const PointType & p1,const PointType & p2)36 template <class PointType> double theta( const PointType& p1, const PointType& p2 )
37 {
38     typename PointType::value_type dx, dy, ax, ay;
39     double t;
40 
41     dx = p2.x - p1.x; ax = absval( dx );
42     dy = p2.y - p1.y; ay = absval( dy );
43     t = (ax+ay == 0) ? 0 : (double) dy/(ax+ay);
44     if( dx < 0 )
45         t = 2-t;
46     else if( dy < 0 )
47         t = 4+t;
48 
49     return t*90.0;
50 }
51 
52 /* Model of LessThanComparable for theta sort.
53  * Uses the theta function from Sedgewick: Algorithms in XXX, chapter 24
54  */
55 template <class PointType> class ThetaCompare : public ::std::binary_function< const PointType&, const PointType&, bool >
56 {
57 public:
ThetaCompare(const PointType & rRefPoint)58     ThetaCompare( const PointType& rRefPoint ) : maRefPoint( rRefPoint ) {}
59 
operator ()(const PointType & p1,const PointType & p2)60     bool operator() ( const PointType& p1, const PointType& p2 )
61     {
62         // return true, if p1 precedes p2 in the angle relative to maRefPoint
63         return theta(maRefPoint, p1) < theta(maRefPoint, p2);
64     }
65 
operator ()(const PointType & p) const66     double operator() ( const PointType& p ) const
67     {
68         return theta(maRefPoint, p);
69     }
70 
71 private:
72     PointType maRefPoint;
73 };
74 
75 /* Implementation of the predicate 'counter-clockwise' for three points, from Sedgewick: Algorithms in XXX, chapter 24 */
ccw(const PointType & p0,const PointType & p1,const PointType & p2,CmpFunctor & thetaCmp)76 template <class PointType, class CmpFunctor> typename PointType::value_type ccw( const PointType& p0, const PointType& p1, const PointType& p2, CmpFunctor& thetaCmp )
77 {
78     typename PointType::value_type dx1, dx2, dy1, dy2;
79     typename PointType::value_type theta0( thetaCmp(p0) );
80     typename PointType::value_type theta1( thetaCmp(p1) );
81     typename PointType::value_type theta2( thetaCmp(p2) );
82 
83 #if 0
84     if( theta0 == theta1 ||
85         theta0 == theta2 ||
86         theta1 == theta2 )
87     {
88         // cannot reliably compare, as at least two points are
89         // theta-equal. See bug description below
90         return 0;
91     }
92 #endif
93 
94     dx1 = p1.x - p0.x; dy1 = p1.y - p0.y;
95     dx2 = p2.x - p0.x; dy2 = p2.y - p0.y;
96 
97     if( dx1*dy2 > dy1*dx2 )
98         return +1;
99 
100     if( dx1*dy2 < dy1*dx2 )
101         return -1;
102 
103     if( (dx1*dx2 < 0) || (dy1*dy2 < 0) )
104         return -1;
105 
106     if( (dx1*dx1 + dy1*dy1) < (dx2*dx2 + dy2*dy2) )
107         return +1;
108 
109     return 0;
110 }
111 
112 /*
113  Bug
114  ===
115 
116  Sometimes, the resulting polygon is not the convex hull (see below
117  for an edge configuration to reproduce that problem)
118 
119  Problem analysis:
120  =================
121 
122  The root cause of this bug is the fact that the second part of
123  the algorithm (the 'wrapping' of the point set) relies on the
124  previous theta sorting. Namely, it is required that the
125  generated point ordering, when interpreted as a polygon, is not
126  self-intersecting. This property, although, cannot be
127  guaranteed due to limited floating point accuracy. For example,
128  for two points very close together, and at the same time very
129  far away from the theta reference point p1, can appear on the
130  same theta value (because floating point accuracy is limited),
131  although on different rays to p1 when inspected locally.
132 
133  Example:
134 
135                     /
136              P3    /
137                |\ /
138                | /
139                |/ \
140              P2    \
141                     \
142       ...____________\
143      P1
144 
145  Here, P2 and P3 are theta-equal relative to P1, but the local
146  ccw measure always says 'left turn'. Thus, the convex hull is
147  wrong at this place.
148 
149  Solution:
150  =========
151 
152  If two points are theta-equal and checked via ccw, ccw must
153  also classify them as 'equal'. Thus, the second stage of the
154  convex hull algorithm sorts the first one out, effectively
155  reducing a cluster of theta-equal points to only one. This
156  single point can then be treated correctly.
157 */
158 
159 
160 /* Implementation of Graham's convex hull algorithm, see Sedgewick: Algorithms in XXX, chapter 25 */
convexHull(const Polygon2D & rPoly)161 Polygon2D convexHull( const Polygon2D& rPoly )
162 {
163     const Polygon2D::size_type N( rPoly.size() );
164     Polygon2D result( N + 1 );
165     ::std::copy(rPoly.begin(), rPoly.end(), result.begin()+1 );
166     Polygon2D::size_type min, i;
167 
168     // determine safe point on hull (smallest y value)
169     for( min=1, i=2; i<=N; ++i )
170     {
171         if( result[i].y < result[min].y )
172             min = i;
173     }
174 
175     // determine safe point on hull (largest x value)
176     for( i=1; i<=N; ++i )
177     {
178         if( result[i].y == result[min].y &&
179             result[i].x > result[min].x )
180         {
181             min = i;
182         }
183     }
184 
185     // TODO: add inner elimination optimization from Sedgewick: Algorithms in XXX, chapter 25
186     // TODO: use radix sort instead of ::std::sort(), calc theta only once and store
187 
188     // setup first point and sort
189     ::std::swap( result[1], result[min] );
190     ThetaCompare<Point2D> cmpFunc(result[1]);
191     ::std::sort( result.begin()+1, result.end(), cmpFunc );
192 
193     // setup sentinel
194     result[0] = result[N];
195 
196     // generate convex hull
197     Polygon2D::size_type M;
198     for( M=3, i=4; i<=N; ++i )
199     {
200         while( ccw(result[M], result[M-1], result[i], cmpFunc) >= 0 )
201             --M;
202 
203         ++M;
204         ::std::swap( result[M], result[i] );
205     }
206 
207     // copy range [1,M] to output
208     return Polygon2D( result.begin()+1, result.begin()+M+1 );
209 }
210