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