basegfx/source/workbench/convexhull.cxx

   1 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
   2 /*
   3  * This file is part of the LibreOffice project.
   4  *
   5  * This Source Code Form is subject to the terms of the Mozilla Public
   6  * License, v. 2.0. If a copy of the MPL was not distributed with this
   7  * file, You can obtain one at http://mozilla.org/MPL/2.0/.
   8  *
   9  * This file incorporates work covered by the following license notice:
  10  *
  11  *   Licensed to the Apache Software Foundation (ASF) under one or more
  12  *   contributor license agreements. See the NOTICE file distributed
  13  *   with this work for additional information regarding copyright
  14  *   ownership. The ASF licenses this file to you under the Apache
  15  *   License, Version 2.0 (the "License"); you may not use this file
  16  *   except in compliance with the License. You may obtain a copy of
  17  *   the License at http://www.apache.org/licenses/LICENSE-2.0 .
  18  */
  19
  20 #include <algorithm>
  21 #include <functional>
  22 #include <vector>
  23
  24 #include "bezierclip.hxx"
  25
  26 // -----------------------------------------------------------------------------
  27
  28 /* Implements the theta function from Sedgewick: Algorithms in XXX, chapter 24 */
  29 template <class PointType> double theta( const PointType& p1, const PointType& p2 )
  30 {
  31     typename PointType::value_type dx, dy, ax, ay;
  32     double t;
  33
  34     dx = p2.x - p1.x; ax = absval( dx );
  35     dy = p2.y - p1.y; ay = absval( dy );
  36     t = (ax+ay == 0) ? 0 : (double) dy/(ax+ay);
  37     if( dx < 0 )
  38         t = 2-t;
  39     else if( dy < 0 )
  40         t = 4+t;
  41
  42     return t*90.0;
  43 }
  44
  45 /* Model of LessThanComparable for theta sort.
  46  * Uses the theta function from Sedgewick: Algorithms in XXX, chapter 24
  47  */
  48 template <class PointType> class ThetaCompare : public ::std::binary_function< const PointType&, const PointType&, bool >
  49 {
  50 public:
  51     ThetaCompare( const PointType& rRefPoint ) : maRefPoint( rRefPoint ) {}
  52
  53     bool operator() ( const PointType& p1, const PointType& p2 )
  54     {
  55         // return true, if p1 precedes p2 in the angle relative to maRefPoint
  56         return theta(maRefPoint, p1) < theta(maRefPoint, p2);
  57     }
  58
  59     double operator() ( const PointType& p ) const
  60     {
  61         return theta(maRefPoint, p);
  62     }
  63
  64 private:
  65     PointType maRefPoint;
  66 };
  67
  68 /* Implementation of the predicate 'counter-clockwise' for three points, from Sedgewick: Algorithms in XXX, chapter 24 */
  69 template <class PointType, class CmpFunctor> typename PointType::value_type ccw( const PointType& p0, const PointType& p1, const PointType& p2, CmpFunctor& thetaCmp )
  70 {
  71     typename PointType::value_type dx1, dx2, dy1, dy2;
  72     typename PointType::value_type theta0( thetaCmp(p0) );
  73     typename PointType::value_type theta1( thetaCmp(p1) );
  74     typename PointType::value_type theta2( thetaCmp(p2) );
  75
  76 #if 0
  77     if( theta0 == theta1 ||
  78         theta0 == theta2 ||
  79         theta1 == theta2 )
  80     {
  81         // cannot reliably compare, as at least two points are
  82         // theta-equal. See bug description below
  83         return 0;
  84     }
  85 #endif
  86
  87     dx1 = p1.x - p0.x; dy1 = p1.y - p0.y;
  88     dx2 = p2.x - p0.x; dy2 = p2.y - p0.y;
  89
  90     if( dx1*dy2 > dy1*dx2 )
  91         return +1;
  92
  93     if( dx1*dy2 < dy1*dx2 )
  94         return -1;
  95
  96     if( (dx1*dx2 < 0) || (dy1*dy2 < 0) )
  97         return -1;
  98
  99     if( (dx1*dx1 + dy1*dy1) < (dx2*dx2 + dy2*dy2) )
 100         return +1;
 101
 102     return 0;
 103 }
 104
 105 /*
 106  Bug
 107  ===
 108
 109  Sometimes, the resulting polygon is not the convex hull (see below
 110  for an edge configuration to reproduce that problem)
 111
 112  Problem analysis:
 113  =================
 114
 115  The root cause of this bug is the fact that the second part of
 116  the algorithm (the 'wrapping' of the point set) relies on the
 117  previous theta sorting. Namely, it is required that the
 118  generated point ordering, when interpreted as a polygon, is not
 119  self-intersecting. This property, although, cannot be
 120  guaranteed due to limited floating point accuracy. For example,
 121  for two points very close together, and at the same time very
 122  far away from the theta reference point p1, can appear on the
 123  same theta value (because floating point accuracy is limited),
 124  although on different rays to p1 when inspected locally.
 125
 126  Example:
 127
 128                     /
 129              P3    /
 130                |\ /
 131                | /
 132                |/ \
 133              P2    \
 134                     \
 135       ...____________\
 136      P1
 137
 138  Here, P2 and P3 are theta-equal relative to P1, but the local
 139  ccw measure always says 'left turn'. Thus, the convex hull is
 140  wrong at this place.
 141
 142  Solution:
 143  =========
 144
 145  If two points are theta-equal and checked via ccw, ccw must
 146  also classify them as 'equal'. Thus, the second stage of the
 147  convex hull algorithm sorts the first one out, effectively
 148  reducing a cluster of theta-equal points to only one. This
 149  single point can then be treated correctly.
 150 */
 151
 152
 153 /* Implementation of Graham's convex hull algorithm, see Sedgewick: Algorithms in XXX, chapter 25 */
 154 Polygon2D convexHull( const Polygon2D& rPoly )
 155 {
 156     const Polygon2D::size_type N( rPoly.size() );
 157     Polygon2D result( N + 1 );
 158     ::std::copy(rPoly.begin(), rPoly.end(), result.begin()+1 );
 159     Polygon2D::size_type min, i;
 160
 161     // determine safe point on hull (smallest y value)
 162     for( min=1, i=2; i<=N; ++i )
 163     {
 164         if( result[i].y < result[min].y )
 165             min = i;
 166     }
 167
 168     // determine safe point on hull (largest x value)
 169     for( i=1; i<=N; ++i )
 170     {
 171         if( result[i].y == result[min].y &&
 172             result[i].x > result[min].x )
 173         {
 174             min = i;
 175         }
 176     }
 177
 178     // TODO: add inner elimination optimization from Sedgewick: Algorithms in XXX, chapter 25
 179     // TODO: use radix sort instead of ::std::sort(), calc theta only once and store
 180
 181     // setup first point and sort
 182     ::std::swap( result[1], result[min] );
 183     ThetaCompare<Point2D> cmpFunc(result[1]);
 184     ::std::sort( result.begin()+1, result.end(), cmpFunc );
 185
 186     // setup sentinel
 187     result[0] = result[N];
 188
 189     // generate convex hull
 190     Polygon2D::size_type M;
 191     for( M=3, i=4; i<=N; ++i )
 192     {
 193         while( ccw(result[M], result[M-1], result[i], cmpFunc) >= 0 )
 194             --M;
 195
 196         ++M;
 197         ::std::swap( result[M], result[i] );
 198     }
 199
 200     // copy range [1,M] to output
 201     return Polygon2D( result.begin()+1, result.begin()+M+1 );
 202 }
 203
 204 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */