2D Misc

From EOJ Wiki
Revision as of 00:50, 11 March 2018 by Ultmaster (talk | contribs) (Created page with "== Define == <syntaxhighlight lang='cpp'> #include <bits/stdc++.h> using namespace std; // this can be safely changed to long double typedef double ld; const ld INF = 1e100;...")
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

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$. 

vector<PT> CircleCircleIntersection(PT a, PT b, ld r, ld R) {
    vector<PT> ret;
    ld d = sqrt(dist2(a, b));
    if (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;
}
Retrieved from "http://acm.ecnu.edu.cn/wiki/index.php?title=2D_Misc&oldid=66"