basegfx/source/workbench/convexhull.cxx

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