Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update halfplane-intersection.md #1367

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Update halfplane-intersection.md #1367

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
90 changes: 53 additions & 37 deletions src/geometry/halfplane-intersection.md
View file
Open in desktop
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 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);
}
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
// Actual algorithm
// Actual Algorithm
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 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<Point>();
// 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<Point>();

// 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;
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
Loading