Skip to content

Instantly share code, notes, and snippets.

@zhenghaoz

zhenghaoz/geometry.hpp

Last active June 9, 2018 12:48
Show Gist options
  • Save zhenghaoz/8abaea616d74aa29c5f9fdc9ec9c4817 to your computer and use it in GitHub Desktop.
Save zhenghaoz/8abaea616d74aa29c5f9fdc9ec9c4817 to your computer and use it in GitHub Desktop.
Download ZIP
Geometry
Raw
geometry.hpp
#include <vector>
#include <algorithm>
#include <iostream>
#include <utility>
#include <cmath>
#include <exception>
// verbose output
/*#define VERBOSE*/
#ifdef VERBOSE
#define VLOG(text) do { std::cout << text << std::endl; } while (0)
#else
#define VLOG(text) do {} while (0);
#endif
// point
// point data structure
struct Point
{
int x, y;
Point(): Point(0, 0) {}
Point(int x, int y): x(x), y(y) {}
Point operator+(const Point &p2) const { return Point(x+p2.x, y+p2.y); }
Point operator-(const Point &p2) const { return Point(x-p2.x, y-p2.y); }
bool operator==(const Point &p2) const { return x==p2.x && y==p2.y; }
bool operator!=(const Point &p2) const { return !((*this)==p2); }
};
std::ostream &operator<<(std::ostream &out, const Point &point)
{
out << '(' << point.x << ',' << point.y << ')';
return out;
}
#define INNER_PRODUCT(p1,p2) ((p1).x*(p2).y-(p1).y*(p2).x)
// segment
// segment data structure
struct Segment
{
Point start, end;
Segment() = default;
Segment(const Point &start, const Point &end): start(start), end(end) {}
bool operator==(const Segment &seg2) const { return start == seg2.start && end == seg2.end; }
bool operator!=(const Segment &seg2) const { return !((*this) == seg2); }
};
std::ostream &operator<<(std::ostream &out, const Segment &seg)
{
out << seg.start << "->" << seg.end;
return out;
}
// return positive number if p1->p3 is on the clockwise direction of p1->p2, negative number if anticlockwise, zero elsewise
#define DIRECTION(p1,p2,p3) (INNER_PRODUCT((p3)-(p1),(p2)-(p1)))
#define MIN(a,b) ((a)<(b)?(a):(b))
#define MAX(a,b) ((a)>(b)?(a):(b))
// when p and s are on the same line, return true if p on s
#define ON_SEGMENT(p,s) (MIN((s).start.x, (s).end.x) <= (p).x && (p).x <= MAX((s).start.x, (s).end.x)\
&& MIN((s).start.y, (s).end.y) <= (p).y && (p).y <= MAX((s).start.y, (s).end.y))
bool intersect(const Segment &seg1, const Segment &seg2)
{
int d1 = DIRECTION(seg1.start, seg1.end, seg2.start);
int d2 = DIRECTION(seg1.start, seg1.end, seg2.end);
int d3 = DIRECTION(seg2.start, seg2.end, seg1.start);
int d4 = DIRECTION(seg2.start, seg2.end, seg1.end);
if (((d1 < 0 && d2 > 0) || (d1 > 0 && d2 < 0))
&& ((d3 < 0 && d4 > 0) || (d3 > 0 && d4 < 0)))
return true;
if (d1 == 0 && ON_SEGMENT(seg2.start, seg1)) return true;
if (d2 == 0 && ON_SEGMENT(seg2.end, seg1)) return true;
if (d3 == 0 && ON_SEGMENT(seg1.start, seg2)) return true;
if (d4 == 0 && ON_SEGMENT(seg1.end, seg2)) return true;
return false;
}
bool intersect(const std::vector<Segment> &segs)
{
// intern class
struct Event
{
bool end;
Point point;
Segment segment;
Event() = default;
Event(bool end, const Point &point, const Segment &segment): end(end), point(point), segment(segment) {}
};
// construct event list
std::vector<Event> events(segs.size()*2);
for (int i = 0; i < segs.size(); i++) {
Point left = segs[i].start;
Point right = segs[i].end;
if (left.x > right.x)
std::swap(left, right);
events[2*i] = Event(false, left, segs[i]);
events[2*i+1] = Event(true, right, segs[i]);
}
std::sort(events.begin(), events.end(), [](const Event &lhs, const Event &rhs)->bool{
if (lhs.point.x != rhs.point.x)
return lhs.point.x < rhs.point.x;
if (lhs.end != lhs.end)
return lhs.end < rhs.end;
return lhs.point.y < rhs.point.y;
});
// check intersect
int x;
auto compare = [&x](const Segment &lhs, const Segment &rhs)->bool{
int dx1 = lhs.end.x - lhs.start.x;
int dx2 = rhs.end.x - rhs.start.x;
int dy1 = lhs.end.y - lhs.start.y;
int dy2 = rhs.end.y - rhs.start.y;
return (dy1*(x-lhs.start.x)+dx1*lhs.start.y)*dx2
< (dy2*(x-rhs.start.x)+dx2*rhs.start.y)*dx1;
};
VLOG("checking intersect...");
std::set<Segment, decltype(compare)> segSet(compare);
for (const Event &event : events) {
x = event.point.x;
if (event.end == false) {
auto ret = segSet.insert(segSet.begin(), event.segment);
auto precursor = std::prev(ret);
auto successor = std::next(ret);
#ifdef VERBOSE
VLOG("inserted "<<event.segment);
if (precursor != segSet.end() && precursor != ret)
VLOG("precursor "<<*precursor);
if (successor != segSet.end())
VLOG("successor "<<*successor);
#endif
if (precursor != segSet.end() && precursor != ret && intersect(*precursor, event.segment))
return true;
if (successor != segSet.end() && intersect(*successor, event.segment))
return true;
} else {
auto ret = segSet.find(event.segment);
auto precursor = std::prev(ret);
auto successor = std::next(ret);
#ifdef VERBOSE
VLOG("removing... "<<event.segment);
if (precursor != segSet.end() && precursor != ret)
VLOG("precursor "<<*precursor);
if (successor != segSet.end())
VLOG("successor "<<*successor);
#endif
if (precursor != segSet.end() && precursor != ret &&
successor != segSet.end() && intersect(*precursor, *successor))
return true;
if (ret != segSet.end())
segSet.erase(ret);
}
}
return false;
}
// convex hull
// assert there is no duplicate points
std::vector<Point> graham(const std::vector<Point> &points)
{
// require more than 2 points
if (points.size() < 3) return std::vector<Point>();
// find the downmost point (choose the leftmost in case of a tie)
int startPointIndex = 0;
for (int i = 1; i < points.size(); i++)
if (points[i].y < points[startPointIndex].y
|| (points[i].y == points[startPointIndex].y && points[i].x < points[startPointIndex].x))
startPointIndex = i;
Point startPoint = points[startPointIndex];
// sort rest points by polar angle (compare distance in case of a tie)
std::vector<Point> restPoints(points.size()-1);
for (int i = 0, offset = 0; i < points.size(); i++)
if (i != startPointIndex)
restPoints[offset++] = points[i];
std::sort(restPoints.begin(), restPoints.end(), [&startPoint](const Point &lhs, const Point &rhs)->bool{
int d = DIRECTION(startPoint, rhs, lhs);
return d == 0 ? abs(lhs.x-startPoint.x) < abs(rhs.x-startPoint.x) : d > 0;
});
VLOG("sorted rest points:");
for (const Point &point : restPoints)
VLOG(point);
// construct convex hull
std::vector<Point> hull = {startPoint};
Point top = restPoints[0];
for (int i = 1; i < restPoints.size(); i++) {
if (DIRECTION(hull[hull.size()-1], top, restPoints[i]) < 0) // turn left
hull.push_back(top);
top = restPoints[i];
}
hull.push_back(top);
return hull;
}
// return true if the polar angle of p1->p3 is greater than p1->p2 (compare polar radius in case of a tie)
bool comparePolar(const Point &p1, const Point &p2, const Point &p3)
{
int d2 = p2.y - p1.y;
int d3 = p3.y - p1.y;
if ((d2 > 0 && d3 > 0) || (d2 < 0 && d3 < 0)) {
int d = DIRECTION(p1, p2, p3);
return (d == 0) ? (abs(p2.x-p1.x) > abs(p3.x-p1.x)) : (d < 0);
}
if (d2 > 0 && d3 < 0)
return false;
if (d2 < 0 && d3 > 0)
return true;
if (d2 == 0 && d3 == 0)
if (p2.x-p1.x > 0 && p3.x-p1.x < 0)
return true;
else if (p2.x-p1.x < 0 && p3.x-p1.x > 0)
return false;
else return abs(p2.x-p1.x) > abs(p3.x-p1.x); // in case of a tie
return DIRECTION(p1, p2, p3) < 0;
}
// assert there is no duplicate points
std::vector<Point> jarvis(const std::vector<Point> &points)
{
// require more than 2 points
if (points.size() < 3) return std::vector<Point>();
// find the downmost point (choose the leftmost in case of a tie)
Point startPoint = points[0];
for (int i = 1; i < points.size(); i++)
if (points[i].y < startPoint.y
|| (points[i].y == startPoint.y && points[i].x < startPoint.x))
startPoint = points[i];
// construct convex hull
std::vector<Point> hull = {startPoint};
Point currentPoint = startPoint;
VLOG("vertexes:");
while (true) {
// find the point with the minimum polar angle (choose the farthest point in case of a tie)
bool first = true;
Point nextPoint;
for (const Point &point : points)
if (point != currentPoint)
if (first) {
first = false;
nextPoint = point;
} else if (comparePolar(currentPoint, point, nextPoint)) {
nextPoint = point;
}
VLOG(nextPoint);
if (nextPoint == startPoint) break;
currentPoint = nextPoint;
hull.push_back(currentPoint);
}
return hull;
}
// closest pair
#define SQUARE(a) ((a)*(a))
#define DISTANCE(a,b) (SQUARE((a).x-(b).x)+SQUARE((a).y-(b).y))
std::pair<int, int> closestPair(const std::vector<Point> &points, std::vector<bool> &bitmap,
const std::vector<int> &x, const std::vector<int> &y)
{
if (x.size() == 2) {
return std::make_pair(x[0], x[1]);
} else if (x.size() == 3) {
int mind = DISTANCE(points[x[0]], points[x[1]]);
int d02 = DISTANCE(points[x[0]], points[x[2]]);
int d12 = DISTANCE(points[x[1]], points[x[2]]);
std::pair<int, int> minp = std::make_pair(x[0], x[1]);
if (d02 < mind) {
mind = d02;
minp = std::make_pair(x[0], x[2]);
}
if (d12 < mind)
minp = std::make_pair(x[1], x[2]);
return minp;
} else {
// split x
int mid = x.size() / 2;
std::vector<int> lx(mid), rx(x.size()-mid);
for (int i = 0; i < mid; i++) {
lx[i] = x[i];
bitmap[x[i]] = false;
}
for (int i = mid; i < x.size(); i++) {
rx[i-mid] = x[i];
bitmap[x[i]] = true;
}
// split y
int lyoffset = 0, ryoffset = 0;
std::vector<int> ly(mid), ry(y.size()-mid);
for (int i = 0; i < y.size(); i++)
if (bitmap[y[i]] == false)
ly[lyoffset++] = y[i];
else
ry[ryoffset++] = y[i];
// get subanswer
std::pair<int, int> p1 = closestPair(points, bitmap, lx, ly);
std::pair<int, int> p2 = closestPair(points, bitmap, rx, ry);
int d1 = DISTANCE(points[p1.first], points[p1.second]);
int d2 = DISTANCE(points[p2.first], points[p2.second]);
int mind = d1;
std::pair<int, int> minp = p1;
if (d2 < d1) {
mind = d2;
minp = p2;
}
// find candidate pairs
std::vector<int> cy;
for (int i = 0; i < x.size(); i++)
if (SQUARE(points[x[i]].x-points[x[mid]].x) <= mind)
cy.push_back(x[i]);
std::sort(cy.begin(), cy.end(), [&points](int lhs, int rhs)->bool{ return points[lhs].y < points[rhs].y; });
for (int i = 0; i < cy.size(); i++)
for (int j = i+1; j < cy.size() && j-i < 8; j++) {
int d = DISTANCE(points[cy[i]], points[cy[j]]);
if (d < mind) {
mind = d;
minp = std::make_pair(cy[i], cy[j]);
}
}
return minp;
}
}
std::pair<int, int> closestPair(const std::vector<Point> &points)
{
// require more than one points
if (points.size() < 2) throw std::invalid_argument("closestPair: require more than one points");
// construct p, x, y, bitmap
std::vector<int> x(points.size());
std::vector<bool> bitmap(points.size());
for (int i = 0; i < x.size(); i++)
x[i] = i;
std::vector<int> y(x);
std::sort(x.begin(), x.end(), [&points](int lhs, int rhs)->bool{ return points[lhs].x < points[rhs].x; });
std::sort(y.begin(), y.end(), [&points](int lhs, int rhs)->bool{ return points[lhs].y < points[rhs].y; });
return closestPair(points, bitmap, x, y);
}
Raw
main.cpp
#define CATCH_CONFIG_MAIN
#include "catch.hpp"
#include "geometry.hpp"
#include <vector>
#include <utility>
TEST_CASE("direction of segment", "[direction]") {
REQUIRE(DIRECTION(Point(), Point(1,2), Point(2,1)) > 0);
REQUIRE(DIRECTION(Point(), Point(2,1), Point(1,2)) < 0);
REQUIRE(DIRECTION(Point(), Point(1,2), Point(1,2)) == 0);
}
TEST_CASE("intersect of two segments", "[intersect]") {
REQUIRE(intersect(Segment(Point(), Point(6,4)), Segment(Point(0,5), Point(5,2))));
REQUIRE(intersect(Segment(Point(), Point(6,4)), Segment(Point(5,2), Point(0,5))));
REQUIRE(intersect(Segment(Point(6,4), Point()), Segment(Point(0,5), Point(5,2))));
REQUIRE(intersect(Segment(Point(6,4), Point()), Segment(Point(5,2), Point(0,5))));
REQUIRE(intersect(Segment(Point(), Point(6,4)), Segment(Point(0,5), Point(3,2))));
REQUIRE(intersect(Segment(Point(), Point(6,4)), Segment(Point(3,2), Point(0,5))));
REQUIRE(intersect(Segment(Point(), Point(3,2)), Segment(Point(0,5), Point(6,-1))));
REQUIRE(intersect(Segment(Point(3,2), Point()), Segment(Point(0,5), Point(6,-1))));
REQUIRE(!intersect(Segment(Point(), Point(6,4)), Segment(Point(0,5), Point(1,1))));
}
TEST_CASE("graham algorithm", "[graham]") {
// {}
std::vector<Point> points, hull;
REQUIRE(graham(points).size() == 0);
// {(0,0)}
points.push_back(Point());
REQUIRE(graham(points).size() == 0);
// {(0,0),(3.0)}
points.emplace_back(3,0);
REQUIRE(graham(points).size() == 0);
// {(0,0),(3.0),(2,2)}
points.emplace_back(2,2);
hull = {Point(0,0), Point(3,0), Point(2,2)};
REQUIRE(graham(points) == hull);
// {(0,0),(3.0),(2,2),(1,1)}
points.emplace_back(1,1);
hull = {Point(0,0), Point(3,0), Point(2,2)};
REQUIRE(graham(points) == hull);
// {(0,0),(3.0),(2,2),(1,1),(4,4)}
points.emplace_back(4,4);
hull = {Point(0,0), Point(3,0), Point(4,4)};
REQUIRE(graham(points) == hull);
// {(0,0),(3.0),(2,2),(1,1),(4,4),(5,6)}
points.emplace_back(5,6);
hull = {Point(0,0), Point(3,0), Point(5,6)};
REQUIRE(graham(points) == hull);
// {(0,0),(3.0),(2,2),(1,1),(4,4),(5,6),(7,12)}
points.emplace_back(7,12);
hull = {Point(0,0), Point(3,0), Point(7,12)};
REQUIRE(graham(points) == hull);
}
TEST_CASE("jarvis algorithm", "[jarvis]") {
// {}
std::vector<Point> points, hull;
REQUIRE(jarvis(points).size() == 0);
// {(0,0)}
points.push_back(Point());
REQUIRE(jarvis(points).size() == 0);
// {(0,0),(3.0)}
points.emplace_back(3,0);
REQUIRE(jarvis(points).size() == 0);
// {(0,0),(3.0),(2,2)}
points.emplace_back(2,2);
hull = {Point(0,0), Point(3,0), Point(2,2)};
REQUIRE(jarvis(points) == hull);
// {(0,0),(3.0),(2,2),(1,1)}
points.emplace_back(1,1);
hull = {Point(0,0), Point(3,0), Point(2,2)};
REQUIRE(jarvis(points) == hull);
// {(0,0),(3.0),(2,2),(1,1),(4,4)}
points.emplace_back(4,4);
hull = {Point(0,0), Point(3,0), Point(4,4)};
REQUIRE(jarvis(points) == hull);
// {(0,0),(3.0),(2,2),(1,1),(4,4),(5,6)}
points.emplace_back(5,6);
hull = {Point(0,0), Point(3,0), Point(5,6)};
REQUIRE(jarvis(points) == hull);
// {(0,0),(3.0),(2,2),(1,1),(4,4),(5,6),(7,12)}
points.emplace_back(7,12);
hull = {Point(0,0), Point(3,0), Point(7,12)};
REQUIRE(jarvis(points) == hull);
}
bool samePair(const std::pair<int, int> &p, int a, int b)
{
return p.first == a && p.second == b || p.first == b && p.second == a;
}
TEST_CASE("closest pair", "[closestPair]") {
// {(0,0),(1,4)}
std::vector<Point> points{Point(), Point(1,4)};
REQUIRE(samePair(closestPair(points), 0, 1));
// {(0,0),(1,4),(3,0)}
points.emplace_back(3,0);
REQUIRE(samePair(closestPair(points), 0, 2));
// {(0,0),(1,4),(3,0),(5,1)}
points.emplace_back(5,1);
REQUIRE(samePair(closestPair(points), 2, 3));
// {(0,0),(1,4),(3,0),(5,1),(0,2)}
points.emplace_back(0,2);
REQUIRE(samePair(closestPair(points), 0, 4));
// {(0,0),(1,4),(3,0),(5,1),(0,2),(1,-1)}
points.emplace_back(1,-1);
REQUIRE(samePair(closestPair(points), 0, 5));
// {(0,0),(1,4),(3,0),(5,1),(0,2),(1,-1),(2,-1)}
points.emplace_back(2,-1);
REQUIRE(samePair(closestPair(points), 5, 6));
// {(0,0),(1,4),(3,0),(5,1),(0,2),(1,-1),(2,-1),(4,4),(3,4)}
points.emplace_back(3,4);
REQUIRE(samePair(closestPair(points), 5, 6));
// {(0,0),(1,4),(3,0),(5,1),(0,2),(1,-1),(2,-1),(4,4),(3,4),(-1,-1)}
points.emplace_back(-1,-1);
REQUIRE(samePair(closestPair(points), 5, 6));
}
TEST_CASE("intersect of segments", "[intersect]") {
// {}
std::vector<Segment> segs;
REQUIRE(intersect(segs) == false);
// {(0,0)->(2,0),(0,3)->(3,3)}
segs.emplace_back(Point(), Point(0,2));
segs.emplace_back(Point(0,3), Point(3,3));
REQUIRE(intersect(segs) == false);
// {(0,0)->(2,0),(0,3)->(3,3)}
segs.clear();
segs.emplace_back(Point(0,0), Point(4,4));
segs.emplace_back(Point(0,4), Point(4,0));
REQUIRE(intersect(segs) == true);
// {(0,0)->(2,0),(0,3)->(3,3),(0,1)->(1,1)}
segs.emplace_back(Point(0,2), Point(1,2));
REQUIRE(intersect(segs) == true);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment