1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
/*
* vim: ts=4 sw=4 et tw=0 wm=0
*
* libcola - A library providing force-directed network layout using the
* stress-majorization method subject to separation constraints.
*
* Copyright (C) 2006-2008 Monash University
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* See the file LICENSE.LGPL distributed with the library.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
*
*/
#include <valarray>
#include <cfloat>
#include <algorithm>
#include "libvpsc/assertions.h"
#include "libcola/convex_hull.h"
namespace hull {
using namespace std;
/**
* CrossProduct of three points: If the result is 0, the points are collinear;
* if it is positive, the three points (in order) constitute a "left turn",
* otherwise a "right turn".
*/
inline double crossProduct(
double x0, double y0,
double x1, double y1,
double x2, double y2) {
return (x1-x0)*(y2-y0)-(x2-x0)*(y1-y0);
}
struct CounterClockwiseOrder {
CounterClockwiseOrder(
const unsigned p,
std::valarray<double> const & X, std::valarray<double> const & Y)
:px(X[p]), py(Y[p]), X(X), Y(Y) {};
bool operator() (unsigned i, unsigned j) {
// o=crossProduct(px,py,X[i],Y[i],X[j],Y[j]);
double xi=X[i]-px;
double xj=X[j]-px;
// since py is the min y pos, yi and yj must be positive
double yi=Y[i]-py;
double yj=Y[j]-py;
// use cross product rule
double o = xi*yj-xj*yi;
if(o!=0) {
return o>0;
}
// in case of ties choose point farthest from p
return (xi*xi+yi*yi) < (xj*xj+yj*yj);
}
const double px;
const double py;
std::valarray<double> const & X;
std::valarray<double> const & Y;
};
void convex(const unsigned n, const double* X, const double* Y, vector<unsigned> & h) {
const valarray<double> XA(X,n);
const valarray<double> YA(Y,n);
convex(XA,YA,h);
}
/**
* Implementation of Graham's scan convex hull finding algorithm.
* X and Y give the horizontal and vertical positions of the pointset.
* The result is returned in hull as a list of indices referencing points in X and Y.
*/
void convex(valarray<double> const & X, valarray<double> const & Y, vector<unsigned> & h) {
unsigned n=X.size();
COLA_ASSERT(n==Y.size());
unsigned p0=0;
// find point p0 with min Y position, choose leftmost in case of tie.
// This is our "pivot" point
double minY=DBL_MAX,minX=DBL_MAX;
for(unsigned i=0;i<n;i++) {
if ( (Y[i] < minY) || ((Y[i] == minY) && (X[i] < minX)) )
{
p0=i;
minY=Y[i];
minX=X[i];
}
}
// sort remaining points by the angle line p0-p1 (p1 in points) makes
// with x-axis
vector<unsigned> points;
for(unsigned i=0;i<n;i++) {
if(i!=p0) points.push_back(i);
}
CounterClockwiseOrder order(p0,X,Y);
sort(points.begin(),points.end(),order);
// now we maintain a stack in h, adding points while each successive
// point is a "left turn", backtracking if we make a right turn.
h.clear();
h.push_back(p0);
h.push_back(points[0]);
for(unsigned i=1;i<points.size();i++) {
double o=crossProduct(
X[h[h.size()-2]],Y[h[h.size()-2]],
X[h[h.size()-1]],Y[h[h.size()-1]],
X[points[i]],Y[points[i]]);
if(o==0) {
h.pop_back();
h.push_back(points[i]);
} else if(o>0) {
h.push_back(points[i]);
} else {
while(o<=0 && h.size()>2) {
h.pop_back();
o=crossProduct(
X[h[h.size()-2]],Y[h[h.size()-2]],
X[h[h.size()-1]],Y[h[h.size()-1]],
X[points[i]],Y[points[i]]);
}
h.push_back(points[i]);
}
}
}
} // namespace hull
|
