Half-Plane Intersection

A half-plane is one side of a directed line; intersecting many describes a convex feasible region. The standard angle-sorted deque algorithm discards boundaries that no longer contribute, giving O(n log n) processing. Unbounded and lower-dimensional outcomes need explicit semantics. This implementation computes the intersection with a user-supplied square [−B,B]² and reports its positive area, or EMPTY OR DEGENERATE; the box is part of the problem, not a substitute for infinity.

Prerequisites: Lines and Intersections; convex polygons; deque; floating-point conditioning.

Build the idea

1. Define the allowed side and numerical contract

   A directed boundary A→C permits points P with cross(C−A,P−A)≥0, its closed left side. Each direction is normalized to unit length so cross(direction,offset) measures signed distance. Inputs use distinct integer endpoints within ±10^4 and 1≤B≤10^4. The numerical model uses distance tolerance 10^−10 and parallel-direction tolerance 10^−14. Regions at or below the area threshold are deliberately grouped with empty or degenerate results.

2. Sort boundaries by direction

   Sort atan2 angles with ordinary strict comparisons, never an epsilon comparator inside std::sort. Processing directions cyclically follows the candidate convex boundary. Same-direction parallel constraints keep only the stricter allowed side. Opposing parallel constraints that become adjacent after pruning cannot support a positive-area bounded polygon in this algorithm; return the non-positive-area category. Nearly parallel configurations remain numerically sensitive.

   $$\mathrm{allowed}\left(P\right)\iff \mathrm{cross}(\mathrm{direction},P-\mathrm{boundaryPoint})\ge 0$$

3. Maintain the deque invariant

   The deque holds candidate boundary lines in angle order, and adjacent intersections are corners of the current candidate region. Before adding a half-plane, remove back corners outside it; do the same at the front. Each removal means that corner cannot belong to the final feasible region, so the associated obsolete boundary can be discarded. The code never intentionally intersects a parallel adjacent pair.

4. Close the cyclic boundary

   After all insertions, prune the back against the front constraint and the front against the back constraint. At least three nonparallel cyclic boundaries are needed for a positive-area polygon. Intersect adjacent lines, including last with first, and use shoelace area. The explicit square guarantees boundedness of the requested clipped region. A half-plane alone may be unbounded, but its intersection with this square is not.

5. Trace the sample and state limitations

   With B=2, boundary (0,1)→(0,−1) keeps x≥0. Boundary (0,0)→(1,0) keeps y≥0. Their clipped feasible region is [0,2]×[0,2], area 4. A contradictory pair or a feasible line/point produces EMPTY OR DEGENERATE. The program does not distinguish those cases or infer whether the original unclipped intersection was bounded. Exact degeneracy classification needs stronger predicates and a richer result representation.

This is a complete program. Compile with C++17 and enter the input described below.

#include<iostream>
#include<vector>
#include<deque>
#include<algorithm>
#include<cmath>
#include<iomanip>
struct Point{
    long double x,y;
    Point operator+(Point p)const{
        return {x+p.x,y+p.y};
    }
    Point operator-(Point p)const{
        return {x-p.x,y-p.y};
    }
    Point operator*(long double t)const{
        return {x*t,y*t};
    }
};
long double cross(Point a,Point b){
    return a.x*b.y-a.y*b.x;
}
long double dot(Point a,Point b){
    return a.x*b.x+a.y*b.y;
}
const long double EPS=1e-10L,PAR=1e-14L;
struct HalfPlane{
    Point p,d;
    long double angle;
    HalfPlane(Point a,Point b):p(a),d(b-a){
        long double len=std::hypotl(d.x,d.y);
        d=d*(1/len);
        angle=std::atan2(d.y,d.x);
    }
    bool outside(Point q)const{
        return cross(d,q-p)<-EPS;
    }
};
bool parallel(HalfPlane a,HalfPlane b){
    return std::fabs(cross(a.d,b.d))<PAR;
}
Point intersection(HalfPlane a,HalfPlane b){
    return a.p+a.d*(cross(b.p-a.p,b.d)/cross(a.d,b.d));
}
std::vector<Point>intersect(std::vector<HalfPlane>lines){
    // Exact angle comparator: never insert an epsilon into std::sort's ordering.
    std::sort(lines.begin(),lines.end(),[](HalfPlane a,HalfPlane b){
        return a.angle<b.angle;
    });
    std::deque<HalfPlane>q;
    for(HalfPlane h:lines){
        while(q.size()>1&&h.outside(intersection(q[q.size()-2],q.back()))){
            q.pop_back();
        }
        while(q.size()>1&&h.outside(intersection(q[0],q[1]))){
            q.pop_front();
        }
        if(!q.empty()&&parallel(q.back(),h)){
            if(dot(q.back().d,h.d)<0){
                return {}; // No positive-area bounded polygon can remain.
            }
            if(h.outside(q.back().p)){
                q.pop_back();
            }
            else{
                continue; // Keep the stricter same-direction constraint.
            }
        }
        q.push_back(h);
    }
    while(q.size()>2&&q.front().outside(intersection(q[q.size()-2],q.back()))){
        q.pop_back();
    }
    while(q.size()>2&&q.back().outside(intersection(q[0],q[1]))){
        q.pop_front();
    }
    if(q.size()<3||parallel(q.back(),q.front())){
        return {};
    }
    std::vector<Point>polygon;
    for(int i=0;i<(int)q.size();++i){
        polygon.push_back(intersection(q[i],q[(i+1)%q.size()]));
    }
    return polygon;
}
int main(){
    int n;
    long double b;
    std::cin>>n>>b;
    std::vector<HalfPlane>lines;
    for(int i=0;i<n;++i){
        Point a,c;
        std::cin>>a.x>>a.y>>c.x>>c.y;
        if(a.x==c.x&&a.y==c.y){
            std::cout<<"INVALID LINE\n";
            return 0;
        }
        lines.emplace_back(a,c);
    }
    // This is an explicit part of the requested region, NOT an infinity approximation.
    lines.emplace_back(Point{-b,-b},Point{b,-b});
    lines.emplace_back(Point{b,-b},Point{b,b});
    lines.emplace_back(Point{b,b},Point{-b,b});
    lines.emplace_back(Point{-b,b},Point{-b,-b});
    auto polygon=intersect(lines);
    long double twice=0;
    for(int i=0;i<(int)polygon.size();++i){
        twice+=cross(polygon[i],polygon[(i+1)%polygon.size()]);
    }
    long double area=std::fabs(twice)/2;
    if(area<=EPS){
        std::cout<<"EMPTY OR DEGENERATE\n";
    }
    else{
        std::cout<<std::fixed<<std::setprecision(6)<<area<<'\n';
    }
}

Complexity analysis

O(n log n) time from angle sorting; each line enters and leaves the deque at most once, so pruning is O(n). O(n) lines/deque/vertex storage. The guarantee is for the stated floating predicate model; almost parallel or tolerance-scale regions may need robust exact alternatives.

Example input and output

Read n B, 0≤n≤100000, integer B in 1..10000. Then n directed lines ax ay bx by with distinct integer endpoints within ±10^4. Compute left-half-plane intersection clipped to [−B,B]². Print approximate positive area, EMPTY OR DEGENERATE for area≤10^−10, or INVALID LINE.

Input

2 2
0 1 0 -1
0 0 1 0

Output

4.000000

The two allowed sides are x≥0 and y≥0; clipping to the square leaves its first-quadrant quarter.

Practical applications

Use for bounded linear constraints, visibility kernels and convex-region intersections. Repeated polygon clipping is simpler for small inputs but can cost quadratic time. General unbounded feasibility and exact point/segment outputs require additional handling beyond this bounded-area implementation.

Check your understanding

What does the clipping square mean in this implementation?

Put it into practice

Compare the deque area with incremental clipping on small random half-planes.