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