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