2D Misc

Define

#include <bits/stdc++.h>
using namespace std;
// this can be safely changed to long double
typedef double ld;

const ld INF = 1e100;
const ld EPS = 1e-12;
const ld PI = acos(-1.0);

// -1: less than 0
// 1: greater than 0
inline int dcmp(ld x) {
    return fabs(x) < EPS ? 0 : (x > 0 ? 1 : -1);
}

struct PT {
    ld x, y;
    PT() { x = 0; y = 0; }
    PT(ld x, ld y): x(x), y(y) {}
    PT(const PT &p) : x(p.x), y(p.y) {}
    PT operator + (const PT &p) const { return PT(x + p.x, y + p.y); }
    PT operator - (const PT &p) const { return PT(x - p.x, y - p.y); }
    PT operator * (ld c) const { return PT(x * c, y * c); }
    PT operator / (ld c) const { return PT(x / c, y / c); }
    inline int operator < (const PT& p) const {
        return dcmp(x - p.x) < 0 || (dcmp(x - p.x) == 0 && dcmp(y - p.y) < 0);
    }
    inline bool operator == (const PT& p) const { return dcmp(x - p.x) == 0 && dcmp(y - p.y) == 0; }
    inline bool operator != (const PT& p) const { return !((*this) == p); }
};

ld dot(PT p, PT q) { return p.x * q.x + p.y * q.y; }
ld dist2(PT p, PT q) { return dot(p - q, p - q); }
ld cross(PT p, PT q) { return p.x * q.y - p.y * q.x; }
ostream &operator << (ostream &os, const PT &p) {
    return (os << "(" << p.x << "," << p.y << ")");
}
istream &operator >> (istream &is, PT &p) {
    return (is >> p.x >> p.y);
}

Basics

Rotate a point CCW or CW around the origin

PT RotateCCW90(PT p) { return PT(-p.y, p.x); }
PT RotateCW90(PT p) { return PT(p.y, -p.x); }
PT RotateCCW(PT p, ld t) {
    return PT(p.x * cos(t) - p.y * sin(t), p.x * sin(t) + p.y * cos(t));
}

Project point $c$ onto line through $a$ and $b$, assuming $a \ne b$

PT ProjectPointLine(PT a, PT b, PT c) {
    return a + (b - a) * dot(c - a, b - a) / dot(b - a, b - a);
}

Project point $c$ onto line segment through $a$ and $b$

PT ProjectPointSegment(PT a, PT b, PT c) {
    ld r = dot(b - a, b - a);
    if (dcmp(r) == 0) return a;
    r = dot(c - a, b - a) / r;
    if (r < 0) return a;
    if (r > 1) return b;
    return a + (b - a) * r;
}

Is point $p$ on line through $a$ and $b$?

inline bool PointOnLine(PT a, PT b, PT p) { return dcmp(cross(a - p, b - p)) == 0; }

Is point $p$ on segment from $a$ to $b$?

inline bool PointOnSegment(PT a, PT b, PT p) { return PointOnLine(a, b, p) &&
                                                dcmp(dot(a - p, b - p)) <= 0; }

ld DistanceBetweenPoints(PT a, PT b) { return sqrt(dist2(a, b)); }
ld DistancePointSegment(PT a, PT b, PT c) { return sqrt(dist2(c, ProjectPointSegment(a, b, c))); }
// 3D: compute distance between point(x, y, z) and plane ax+by+cz=d
ld DistancePointPlane(ld x, ld y, ld z, ld a, ld b, ld c, ld d) {
    return fabs(a * x + b * y + c * z - d) / sqrt(a * a + b * b + c * c);
}

Determine if lines from $a$ to $b$ and $c$ to $d$ and parallel or collinear

bool LinesParallel(PT a, PT b, PT c, PT d) { return dcmp(cross(b - a, c - d)) == 0; }
bool LinesCollinear(PT a, PT b, PT c, PT d) {
    return LinesParallel(a, b, c, d) &&
           dcmp(cross(a - b, a - c)) == 0 && dcmp(cross(c - d, c - a)) == 0;
}

Determine if line segment from $a$ to $b$ intersects with line segment from $c$ to $d$ (without points)

bool SegmentsProperIntersect(PT a, PT b, PT c, PT d) {
    ld c1 = cross(b - a, c - a), c2 = cross(b - a, d - a),
           c3 = cross(d - c, a - c), c4 = cross(d - c, b - c);
    return dcmp(c1) * dcmp(c2) < 0 && dcmp(c3) * dcmp(c4) < 0;
}

Determine if line segment from $a$ to $b$ intersects with line segment from $c$ to $d$

bool SegmentsIntersect(PT a, PT b, PT c, PT d) {
    return PointOnSegment(c, d, a) || PointOnSegment(c, d, b) ||
           PointOnSegment(a, b, c) || PointOnSegment(a, b, d) ||
           SegmentsProperIntersect(a, b, c, d);
}

Compute intersection of line passing through $a$ and $b$ with line passing through $c$ and $d$, assuming that unique intersection exists; for segment intersection, check if semgents intersect first

PT ComputeLineIntersection(PT a, PT b, PT c, PT d) {
    b = b - a; d = c - d; c = c - a;
    assert (dcmp(dot(b, b)) > 0 && dcmp(dot(d, d)) > 0);
    return a + b * cross(c, d) / cross(b, d);
}

Compute center of circle given three points

PT ComputeCircleCenter(PT a, PT b, PT c) {
    b = (a + b) / 2;
    c = (a + c) / 2;
    return ComputeLineIntersection(b, b + RotateCW90(a - b), c , c + RotateCW90(a - c));
}

Determine if point is in a possibly non-convex polygon (by William Randolph Franklin); returns 1 for strictly interior points, 0 for strictly exterior points, and 0 or 1 for the remaining points. Note that it is possible to convert this into an exact test using integer arithmetic by taking care of the division appropriately (making sure to deal with signs properly) and then by writing exact tests for checking point on polygon boundary

bool PointInPolygon(const vector<PT> &p, PT q) {
    bool c = 0;
    for (unsigned i = 0; i < p.size(); ++i) {
        unsigned j = (i + 1) % p.size();
        if (((dcmp(p[i].y - q.y) <= 0 && dcmp(q.y - p[j].y) < 0) ||
             (dcmp(p[j].y - q.y) <= 0 && dcmp(q.y - p[i].y) < 0)) &&
            dcmp(q.x - (p[i].x + (p[j].x - p[i].x) * (q.y - p[i].y) / (p[j].y - p[i].y))) < 0)
            c = !c;
    }
    return c;
}

Determine if point is on the boundary of a polygon

bool PointOnPolygon (const vector<PT> &p, PT q){
    for (unsigned i = 0; i < p.size(); ++i)
        if (PointOnSegment(p[i], p[(i + 1) % p.size()], q))
            return true;
  return false;
}

Compute intersection of line through points $a$ and $b$ with circle centered at $c$ with radius $r > 0$

vector<PT> CircleLineIntersection(PT a, PT b, PT c, ld r) {
    vector<PT> ret;
    b = b - a;
    a = a - c;
    ld A = dot(b, b), B = dot(a, b), C = dot(a, a) - r * r;
    ld D = B * B - A * C;
    if (dcmp(D) < 0) return ret;
    ret.push_back(c + a + b * (-B + sqrt(D + EPS)) / A);
    if (dcmp(D) > 0)
        ret.push_back(c + a + b * (-B - sqrt(D)) / A);
    return ret;
}

Compute intersection of circle centered at $a$ with radius $r$ with circle centered at $b$ with radius $R$.

If they have the same center, nothing will be returned.

vector<PT> CircleCircleIntersection(PT a, PT b, ld r, ld R) {
    vector<PT> ret;
    ld d = sqrt(dist2(a, b));
    if (dcmp(d) == 0 || dcmp(d - (r + R)) > 0 || dcmp(d + min(r, R) - max(r, R)) < 0) return ret;
    ld x = (d * d - R * R + r * r) / (2 * d);
    ld y = sqrt(r * r - x * x);
    PT v = (b - a) / d;
    ret.push_back(a + v * x + RotateCCW90(v) * y);
    if (dcmp(y) > 0)
        ret.push_back(a + v * x - RotateCCW90(v) * y);
    return ret;
}

Counterclockwise fashion. Note that the centroid is often known as the "center of gravity" or "center of mass".

ld ComputeSignedArea(const vector<PT> &p) {
    ld area = 0;
    for (unsigned i = 0; i < p.size(); i++) {
        unsigned j = (i + 1) % p.size();
        area += p[i].x * p[j].y - p[j].x * p[i].y;
    }
    return area / 2.0;
}

ld ComputeArea(const vector<PT> &p) {
    return fabs(ComputeSignedArea(p));
}

PT ComputeCentroid(const vector<PT> &p) {
    PT c(0, 0);
    ld scale = 6.0 * ComputeSignedArea(p);
    for (unsigned i = 0; i < p.size(); i++) {
        unsigned j = (i + 1) % p.size();
        c = c + (p[i] + p[j]) * (p[i].x * p[j].y - p[j].x * p[i].y);
    }
    return c / scale;
}

Tests whether or not a given polygon (in CW or CCW order) is simple

bool IsSimple(const vector <PT> &p) {
    for (unsigned i = 0; i < p.size(); i++) {
        for (unsigned k = i + 1; k < p.size(); k++) {
            unsigned j = (i + 1) % p.size();
            unsigned l = (k + 1) % p.size();
            if (i == l || j == k) continue;
            if (SegmentsIntersect(p[i], p[j], p[k], p[l]))
                return false;
        }
    }
    return true;
}

Test Example

int main() {
    // expected: (-5,2)
    cerr << RotateCCW90(PT(2,5)) << endl;
    // expected: (5,-2)
    cerr << RotateCW90(PT(2,5)) << endl;
    // expected: (-5,2)
    cerr << RotateCCW(PT(2,5),M_PI/2) << endl;
    // expected: (5,2)
    cerr << ProjectPointLine(PT(-5,-2), PT(10,4), PT(3,7)) << endl;
    // expected: (5,2) (7.5,3) (2.5,1)
    cerr << ProjectPointSegment(PT(-5,-2), PT(10,4), PT(3,7)) << " "
        << ProjectPointSegment(PT(7.5,3), PT(10,4), PT(3,7)) << " "
        << ProjectPointSegment(PT(-5,-2), PT(2.5,1), PT(3,7)) << endl;
    // expected: 6.78903
    cerr << DistancePointPlane(4,-4,3,2,-2,5,-8) << endl;
    // expected: 1 0 1
    cerr << LinesParallel(PT(1,1), PT(3,5), PT(2,1), PT(4,5)) << " "
        << LinesParallel(PT(1,1), PT(3,5), PT(2,0), PT(4,5)) << " "
        << LinesParallel(PT(1,1), PT(3,5), PT(5,9), PT(7,13)) << endl;
    // expected: 0 0 1
    cerr << LinesCollinear(PT(1,1), PT(3,5), PT(2,1), PT(4,5)) << " "
        << LinesCollinear(PT(1,1), PT(3,5), PT(2,0), PT(4,5)) << " "
        << LinesCollinear(PT(1,1), PT(3,5), PT(5,9), PT(7,13)) << endl;
    // expected: 1 1 1 0
    cerr << SegmentsIntersect(PT(0,0), PT(2,4), PT(3,1), PT(-1,3)) << " "
        << SegmentsIntersect(PT(0,0), PT(2,4), PT(4,3), PT(0,5)) << " "
        << SegmentsIntersect(PT(0,0), PT(2,4), PT(2,-1), PT(-2,1)) << " "
        << SegmentsIntersect(PT(0,0), PT(2,4), PT(5,5), PT(1,7)) << endl;
    // expected: (1,2)
    cerr << ComputeLineIntersection(PT(0,0), PT(2,4), PT(3,1), PT(-1,3)) << endl;
    // expected: (1,1)
    cerr << ComputeCircleCenter(PT(-3,4), PT(6,1), PT(4,5)) << endl;
    vector<PT> v;
    v.push_back(PT(0,0));
    v.push_back(PT(5,0));
    v.push_back(PT(5,5));
    v.push_back(PT(0,5));
    // expected: 1 1 1 0 0
    cerr << PointInPolygon(v, PT(2,2)) << " "
        << PointInPolygon(v, PT(2,0)) << " "
        << PointInPolygon(v, PT(0,2)) << " "
        << PointInPolygon(v, PT(5,2)) << " "
        << PointInPolygon(v, PT(2,5)) << endl;
    // expected: 0 1 1 1 1
    cerr << PointOnPolygon(v, PT(2,2)) << " "
        << PointOnPolygon(v, PT(2,0)) << " "
        << PointOnPolygon(v, PT(0,2)) << " "
        << PointOnPolygon(v, PT(5,2)) << " "
        << PointOnPolygon(v, PT(2,5)) << endl;
    // expected: (1,6)
    // (5,4) (4,5)
    // blank line
    // (4,5) (5,4)
    // blank line
    // (4,5) (5,4)
    vector<PT> u = CircleLineIntersection(PT(0,6), PT(2,6), PT(1,1), 5);
    for (unsigned i = 0; i < u.size(); i++) cerr << u[i] << " ";
    cerr << endl;
    u = CircleLineIntersection(PT(0,9), PT(9,0), PT(1,1), 5);
    for (unsigned i = 0; i < u.size(); i++) cerr << u[i] << " ";
    cerr << endl;
    u = CircleCircleIntersection(PT(1,1), PT(10,10), 5, 5);
    for (unsigned i = 0; i < u.size(); i++) cerr << u[i] << " ";
    cerr << endl;
    u = CircleCircleIntersection(PT(1,1), PT(8,8), 5, 5);
    for (unsigned i = 0; i < u.size(); i++) cerr << u[i] << " ";
    cerr << endl;
    u = CircleCircleIntersection(PT(1,1), PT(4.5,4.5), 10, sqrt(2.0)/2.0);
    for (unsigned i = 0; i < u.size(); i++) cerr << u[i] << " ";
    cerr << endl;
    u = CircleCircleIntersection(PT(1,1), PT(4.5,4.5), 5, sqrt(2.0)/2.0);
    for (unsigned i = 0; i < u.size(); i++) cerr << u[i] << " ";
    cerr << endl;
    // area should be 5.0
    // centroid should be (1.1666666, 1.166666)
    PT pa[] = { PT(0,0), PT(5,0), PT(1,1), PT(0,5) };
    vector<PT> p(pa, pa+4);
    PT c = ComputeCentroid(p);
    cerr << "Area: " << ComputeArea(p) << endl;
    cerr << "Centroid: " << c << endl;
}