diff options
Diffstat (limited to '')
-rw-r--r-- | src/3rdparty/adaptagrams/libcola/convex_hull.cpp | 129 |
1 files changed, 129 insertions, 0 deletions
diff --git a/src/3rdparty/adaptagrams/libcola/convex_hull.cpp b/src/3rdparty/adaptagrams/libcola/convex_hull.cpp new file mode 100644 index 0000000..3ecbfea --- /dev/null +++ b/src/3rdparty/adaptagrams/libcola/convex_hull.cpp @@ -0,0 +1,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 |