Content-Length: 6168 | pFad | http://github.com/cp-algorithms/cp-algorithms/pull/1367.diff

thub.com diff --git a/src/geometry/halfplane-intersection.md b/src/geometry/halfplane-intersection.md index c0334fd3a..17f90186e 100644 --- a/src/geometry/halfplane-intersection.md +++ b/src/geometry/halfplane-intersection.md @@ -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 double eps = 1e-9, inf = 1e9; +const long double eps = 1e-8; +const long double inf = 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); } @@ -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) { @@ -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() {} @@ -131,10 +127,10 @@ 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) { @@ -142,12 +138,31 @@ struct Halfplane { return s.p + (s.pq * alpha); } }; + +// Function to calculate a dynamic bounding box based on input points. +Point calculate_bounding_box(const std::vector& 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 -// Actual algorithm +// Actual Algorithm vector hp_intersect(vector& H) { Point box[4] = { // Bounding box in CCW order @@ -157,7 +172,7 @@ vector hp_intersect(vector& 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); } @@ -166,40 +181,41 @@ vector hp_intersect(vector& H) { sort(H.begin(), H.end()); deque 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 while the last 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 while the first 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(); - - // Same direction half-plane: keep only 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(); + + // Keep only 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 + + // Add the new half-plane dq.push_back(H[i]); ++len; } - // Final cleanup: Check half-planes at the front against the back and vice-versa + // Final cleanup: Check half-planes at the front and back while (len > 2 && dq[0].out(inter(dq[len-1], dq[len-2]))) { dq.pop_back(); --len; @@ -210,12 +226,12 @@ vector hp_intersect(vector& H) { --len; } - // Report empty intersection if necessary + // If less than 3 planes, report empty intersection if (len < 3) return vector(); - // Reconstruct the convex polygon from the remaining half-planes. + // Reconstruct the convex polygon from the remaining half-planes vector 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]);

pFad - (p)hone/(F)rame/(a)nonymizer/(d)eclutterfier! Saves Data!