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