Oct 15, 2024 · Apr 16, 2025 · Apr 17, 2025
diff --git a/src/geometry/halfplane-intersection.md b/src/geometry/halfplane-intersection.md
 Simple point/vector and half-plane structs:

 ```cpp
 // Redefine epsilon and infinity as necessary. Be mindful of precision errors.
 const long doubleeps =1e-9, inf = 1e9;
 const long double eps = 1e-8;
 const long doubleinf =1e6; // Modify as needed based on input scale.

 // Basic point/vector struct.
 struct Point {

 struct Point {
    long double x, y;
    explicit Point(long double x = 0, long double y = 0) : x(x), y(y) {}

    // Addition, substraction, multiply by constant, dot product, cross product.

    // Addition, subtraction, multiply by constant, dot product, cross product.
    friend Point operator + (const Point& p, const Point& q) {
        return Point(p.x + q.x, p.y + q.y);
    }

    friend Point operator * (const Point& p, const long double& k) {
        return Point(p.x * k, p.y * k);
    }

    }

    friend long double dot(const Point& p, const Point& q) {
 return p.x * q.x + p.y * q.y;
 return p.x * q.x + p.y * q.y;
    }

    friend long double cross(const Point& p, const Point& q) {
 };

 // Basic half-plane struct.
 struct Halfplane {

    // 'p' is a passing point of the line and 'pq' is the direction vector of the line.
    Point p, pq;
 struct Halfplane {
    Point p, pq; // 'p' is a passing point of the line, 'pq' is the direction vector of the line.
    long double angle;

    Halfplane() {}
        return cross(pq, r - p) < -eps;
    }

    // Comparator for sorting.
    // Comparator for sorting.
    bool operator < (const Halfplane& e) const {
        return angle < e.angle;
    }
    }

    // Intersection point of the lines of two half-planes. It is assumed they're never parallel.
    friend Point inter(const Halfplane& s, const Halfplane& t) {
        long double alpha = cross((t.p - s.p), t.pq) / cross(s.pq, t.pq);
        return s.p + (s.pq * alpha);
    }
 };

 // Function to calculate a dynamic bounding box based on input points.
 Point calculate_bounding_box(const std::vector<Point>& points) {
    long double maxX = -inf, maxY = -inf;
    long double minX = inf, minY = inf;

    // Traverse all points to find the bounding box
    for (const Point& pt : points) {
        if (pt.x > maxX) maxX = pt.x;
        if (pt.x < minX) minX = pt.x;
        if (pt.y > maxY) maxY = pt.y;
        if (pt.y < minY) minY = pt.y;
    }

    // Use the most extreme points to define the bounding box
    Point bbox((maxX - minX) * 2, (maxY - minY) * 2);  // Making the box larger for stability
    return bbox;
 }

 ```

 Algorithm:

 ```cpp
 // Actualalgorithm
 // ActualAlgorithm
 vector<Point> hp_intersect(vector<Halfplane>& H) {

    Point box[4] = {  // Bounding box in CCW order
        Point(inf, -inf)
    };

    for(int i = 0; i<4; i++) { // Add bounding box half-planes.
    for(int i = 0; i <4; i++) { // Add bounding box half-planes.
        Halfplane aux(box[i], box[(i+1) % 4]);
        H.push_back(aux);
    }
    sort(H.begin(), H.end());
    deque<Halfplane> dq;
    int len = 0;

    for(int i = 0; i < int(H.size()); i++) {

        // Remove from the back of the deque while last half-plane is redundant
        // Remove from the back of the deque whilethelast half-plane is redundant
        while (len > 1 && H[i].out(inter(dq[len-1], dq[len-2]))) {
            dq.pop_back();
            --len;
        }

        // Remove from the front of the deque while first half-plane is redundant
        // Remove from the front of the deque whilethefirst half-plane is redundant
        while (len > 1 && H[i].out(inter(dq[0], dq[1]))) {
            dq.pop_front();
            --len;
        }

        // Special case check: Parallel half-planes

        // Special case: Check for parallel half-planes
        if (len > 0 && fabsl(cross(H[i].pq, dq[len-1].pq)) < eps) {
 // Opposite parallel half-planes that ended up checked against each other.
 if (dot(H[i].pq, dq[len-1].pq) < 0.0)
 return vector<Point>();

 // Same direction half-plane: keeponly the leftmost half-plane.
 if (H[i].out(dq[len-1].p)) {
 dq.pop_back();
 --len;
 }
 else continue;
 // Opposite parallel half-planes
 if (dot(H[i].pq, dq[len-1].pq) < 0.0)
 return vector<Point>();

    // Keeponly the leftmost half-plane for same direction
 if (H[i].out(dq[len-1].p)) {
 dq.pop_back();
 --len;
 }
 else continue;
        }

        // Add new half-plane

        // Addthenew half-plane
        dq.push_back(H[i]);
        ++len;
    }

    // Final cleanup: Check half-planes at the frontagainst the backandvice-versa
    // Final cleanup: Check half-planes at the front andback
    while (len > 2 && dq[0].out(inter(dq[len-1], dq[len-2]))) {
        dq.pop_back();
        --len;
        --len;
    }

    //Report empty intersection if necessary
    //If less than 3 planes, report empty intersection
    if (len < 3) return vector<Point>();

    // Reconstruct the convex polygon from the remaining half-planes.
    // Reconstruct the convex polygon from the remaining half-planes
    vector<Point> ret(len);
    for(int i = 0; i+1 < len; i++) {
    for(int i = 0; i+1 < len; i++) {
        ret[i] = inter(dq[i], dq[i+1]);
    }
    ret.back() = inter(dq[len-1], dq[0]);
Original file line number	Diff line number	Diff line change
Expand Up		@@ -81,17 +81,15 @@ Here is a sample, direct implementation of the algorithm, with comments explaini
		Simple point/vector and half-plane structs:

		```cpp
		// Redefine epsilon and infinity as necessary. Be mindful of precision errors.
		const long doubleeps =1e-9, inf = 1e9;
		const long double eps = 1e-8;
		const long doubleinf =1e6; // Modify as needed based on input scale.

		// Basic point/vector struct.
		struct Point {

		struct Point {
		long double x, y;
		explicit Point(long double x = 0, long double y = 0) : x(x), y(y) {}

		// Addition, substraction, multiply by constant, dot product, cross product.

		// Addition, subtraction, multiply by constant, dot product, cross product.
		friend Point operator + (const Point& p, const Point& q) {
		return Point(p.x + q.x, p.y + q.y);
		}
Expand All		@@ -102,10 +100,10 @@ struct Point {

		friend Point operator * (const Point& p, const long double& k) {
		return Point(p.x * k, p.y * k);
		}

		}

		friend long double dot(const Point& p, const Point& q) {
		return p.x * q.x + p.y * q.y;
		return p.x * q.x + p.y * q.y;
		}

		friend long double cross(const Point& p, const Point& q) {
Expand All		@@ -114,10 +112,8 @@ struct Point {
		};

		// Basic half-plane struct.
		struct Halfplane {

		// 'p' is a passing point of the line and 'pq' is the direction vector of the line.
		Point p, pq;
		struct Halfplane {
		Point p, pq; // 'p' is a passing point of the line, 'pq' is the direction vector of the line.
		long double angle;

		Halfplane() {}
Expand All		@@ -131,23 +127,42 @@ struct Halfplane {
		return cross(pq, r - p) < -eps;
		}

		// Comparator for sorting.
		// Comparator for sorting.
		bool operator < (const Halfplane& e) const {
		return angle < e.angle;
		}
		}

		// Intersection point of the lines of two half-planes. It is assumed they're never parallel.
		friend Point inter(const Halfplane& s, const Halfplane& t) {
		long double alpha = cross((t.p - s.p), t.pq) / cross(s.pq, t.pq);
		return s.p + (s.pq * alpha);
		}
		};

		// Function to calculate a dynamic bounding box based on input points.
		Point calculate_bounding_box(const std::vector<Point>& points) {
		long double maxX = -inf, maxY = -inf;
		long double minX = inf, minY = inf;

		// Traverse all points to find the bounding box
		for (const Point& pt : points) {
		if (pt.x > maxX) maxX = pt.x;
		if (pt.x < minX) minX = pt.x;
		if (pt.y > maxY) maxY = pt.y;
		if (pt.y < minY) minY = pt.y;
		}

		// Use the most extreme points to define the bounding box
		Point bbox((maxX - minX) * 2, (maxY - minY) * 2); // Making the box larger for stability
		return bbox;
		}

		```

		Algorithm:

		```cpp
		// Actualalgorithm
		// ActualAlgorithm
		vector<Point> hp_intersect(vector<Halfplane>& H) {

		Point box[4] = { // Bounding box in CCW order
Expand All		@@ -157,7 +172,7 @@ vector<Point> hp_intersect(vector<Halfplane>& H) {
		Point(inf, -inf)
		};

		for(int i = 0; i<4; i++) { // Add bounding box half-planes.
		for(int i = 0; i <4; i++) { // Add bounding box half-planes.
		Halfplane aux(box[i], box[(i+1) % 4]);
		H.push_back(aux);
		}
Expand All		@@ -166,40 +181,41 @@ vector<Point> hp_intersect(vector<Halfplane>& H) {
		sort(H.begin(), H.end());
		deque<Halfplane> dq;
		int len = 0;

		for(int i = 0; i < int(H.size()); i++) {

		// Remove from the back of the deque while last half-plane is redundant
		// Remove from the back of the deque whilethelast half-plane is redundant
		while (len > 1 && H[i].out(inter(dq[len-1], dq[len-2]))) {
		dq.pop_back();
		--len;
		}

		// Remove from the front of the deque while first half-plane is redundant
		// Remove from the front of the deque whilethefirst half-plane is redundant
		while (len > 1 && H[i].out(inter(dq[0], dq[1]))) {
		dq.pop_front();
		--len;
		}

		// Special case check: Parallel half-planes

		// Special case: Check for parallel half-planes
		if (len > 0 && fabsl(cross(H[i].pq, dq[len-1].pq)) < eps) {
		// Opposite parallel half-planes that ended up checked against each other.
		if (dot(H[i].pq, dq[len-1].pq) < 0.0)
		return vector<Point>();

		// Same direction half-plane: keeponly the leftmost half-plane.
		if (H[i].out(dq[len-1].p)) {
		dq.pop_back();
		--len;
		}
		else continue;
		// Opposite parallel half-planes
		if (dot(H[i].pq, dq[len-1].pq) < 0.0)
		return vector<Point>();

		// Keeponly the leftmost half-plane for same direction
		if (H[i].out(dq[len-1].p)) {
		dq.pop_back();
		--len;
		}
		else continue;
		}

		// Add new half-plane

		// Addthenew half-plane
		dq.push_back(H[i]);
		++len;
		}

		// Final cleanup: Check half-planes at the frontagainst the backandvice-versa
		// Final cleanup: Check half-planes at the front andback
		while (len > 2 && dq[0].out(inter(dq[len-1], dq[len-2]))) {
		dq.pop_back();
		--len;
Expand All		@@ -210,12 +226,12 @@ vector<Point> hp_intersect(vector<Halfplane>& H) {
		--len;
		}

		//Report empty intersection if necessary
		//If less than 3 planes, report empty intersection
		if (len < 3) return vector<Point>();

		// Reconstruct the convex polygon from the remaining half-planes.
		// Reconstruct the convex polygon from the remaining half-planes
		vector<Point> ret(len);
		for(int i = 0; i+1 < len; i++) {
		for(int i = 0; i+1 < len; i++) {
		ret[i] = inter(dq[i], dq[i+1]);
		}
		ret.back() = inter(dq[len-1], dq[0]);
Expand Down