Half-plane Intersection

struct Line {
    PT p, v;
    double ang;

    Line() {}
    Line(PT from, PT to) : p(from), v(to - from) { ang = atan2(v.y, v.x); }
    friend bool operator<(Line a, Line b) {
        return a.ang < b.ang;
    }
};

bool OnLeft(Line L, PT p) {
    return dcmp(cross(L.v, p - L.p)) >= 0;
}

PT GetIntersection(Line a, Line b) {
    PT u = a.p - b.p;
    ld t = cross(b.v, u) / cross(a.v, b.v);
    return a.p + a.v * t;
}

vector<PT> HalfplaneIntersection(vector<Line>& L) {
    int n = L.size();
    sort(L.begin(), L.end());

    int first, last;
    vector<PT> p(n);
    vector<Line> q(n);
    q[first = last = 0] = L[0];
    for (int i = 1; i < n; i++) {
        while (first < last && !OnLeft(L[i], p[last - 1])) last--;
        while (first < last && !OnLeft(L[i], p[first])) first++;
        q[++last] = L[i];
        if (dcmp(cross(q[last].v, q[last - 1].v)) == 0) {
            last--;
            if (OnLeft(q[last], L[i].p)) q[last] = L[i];
        }
        if (first < last) p[last - 1] = GetIntersection(q[last - 1], q[last]);
    }
    while (first < last && !OnLeft(q[first], p[last - 1])) last--;
    if (last - first <= 1) return vector<PT>();
    p[last] = GetIntersection(q[last], q[first]);

    return vector<PT>(p.begin() + first, p.begin() + last + 1);
}

vector<PT> convexIntersection(const vector<PT> &v1, const vector<PT> &v2) {
    vector<Line> h;
    int n = v1.size(), m = v2.size();
    for (int i = 0; i < n; ++i)
        h.push_back(Line(v1[i], v1[(i + 1) % n]));
    for (int i = 0; i < m; ++i)
        h.push_back(Line(v2[i], v2[(i + 1) % m]));
    return HalfplaneIntersection(h);
}

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));
}