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