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 */ 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: 58 ThetaCompare( const PointType& rRefPoint ) : maRefPoint( rRefPoint ) {} 59 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 66 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 */ 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 */ 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