Skip to main content

Advanced Geometry

Sweep line algorithm, closest pair of points, line segment intersection, polygon operations, and 3D geometry basics

Explore advanced geometric algorithms for complex spatial problems and computational geometry.

Sweep Line Algorithm​

Event-Driven Sweep Line​

struct Event {
    double x;
    int type; // 0: start, 1: end
    int id;
    
    Event(double x, int type, int id) : x(x), type(type), id(id) {}
    
    bool operator<(const Event& other) const {
        if (abs(x - other.x) < 1e-9) {
            return type < other.type; // Start events before end events
        }
        return x < other.x;
    }
};

// Sweep line for line segment intersection
vector<Point> sweepLineIntersection(vector<LineSegment>& segments) {
    vector<Event> events;
    vector<Point> intersections;
    
    // Create events for each segment
    for (int i = 0; i < segments.size(); i++) {
        double x1 = segments[i].p1.x;
        double x2 = segments[i].p2.x;
        
        if (x1 > x2) {
            swap(x1, x2);
        }
        
        events.push_back(Event(x1, 0, i)); // Start event
        events.push_back(Event(x2, 1, i)); // End event
    }
    
    sort(events.begin(), events.end());
    
    set<int> activeSegments;
    
    for (const Event& event : events) {
        if (event.type == 0) {
            // Start event - add segment to active set
            activeSegments.insert(event.id);
            
            // Check intersections with other active segments
            for (int id : activeSegments) {
                if (id != event.id) {
                    Point intersection = lineSegmentIntersection(segments[event.id], segments[id]);
                    if (intersection.x != INFINITY) {
                        intersections.push_back(intersection);
                    }
                }
            }
        } else {
            // End event - remove segment from active set
            activeSegments.erase(event.id);
        }
    }
    
    return intersections;
}

Sweep Line for Rectangle Union​

struct Rectangle {
    double x1, y1, x2, y2;
    
    Rectangle(double x1, double y1, double x2, double y2) 
        : x1(x1), y1(y1), x2(x2), y2(y2) {}
};

// Calculate union area of rectangles using sweep line
double rectangleUnionArea(vector<Rectangle>& rectangles) {
    vector<Event> events;
    
    // Create events for each rectangle
    for (int i = 0; i < rectangles.size(); i++) {
        events.push_back(Event(rectangles[i].x1, 0, i)); // Start
        events.push_back(Event(rectangles[i].x2, 1, i)); // End
    }
    
    sort(events.begin(), events.end());
    
    double totalArea = 0;
    double prevX = 0;
    set<int> activeRectangles;
    
    for (const Event& event : events) {
        if (!activeRectangles.empty()) {
            double width = event.x - prevX;
            double height = calculateActiveHeight(rectangles, activeRectangles);
            totalArea += width * height;
        }
        
        if (event.type == 0) {
            activeRectangles.insert(event.id);
        } else {
            activeRectangles.erase(event.id);
        }
        
        prevX = event.x;
    }
    
    return totalArea;
}

double calculateActiveHeight(vector<Rectangle>& rectangles, set<int>& active) {
    vector<pair<double, int>> yEvents;
    
    for (int id : active) {
        yEvents.push_back({rectangles[id].y1, 1}); // Start
        yEvents.push_back({rectangles[id].y2, -1}); // End
    }
    
    sort(yEvents.begin(), yEvents.end());
    
    double height = 0;
    int count = 0;
    double prevY = 0;
    
    for (const auto& event : yEvents) {
        if (count > 0) {
            height += event.first - prevY;
        }
        count += event.second;
        prevY = event.first;
    }
    
    return height;
}

Closest Pair of Points​

Divide and Conquer Approach​

// Find closest pair of points using divide and conquer
pair<Point, Point> closestPair(vector<Point>& points) {
    int n = points.size();
    if (n < 2) return {Point(), Point()};
    
    // Sort points by x-coordinate
    sort(points.begin(), points.end(), [](const Point& a, const Point& b) {
        return a.x < b.x;
    });
    
    return closestPairRecursive(points, 0, n - 1);
}

pair<Point, Point> closestPairRecursive(vector<Point>& points, int left, int right) {
    if (right - left <= 3) {
        // Base case: brute force for small number of points
        return bruteForceClosestPair(points, left, right);
    }
    
    int mid = (left + right) / 2;
    Point midPoint = points[mid];
    
    // Recursively find closest pair in left and right halves
    pair<Point, Point> leftPair = closestPairRecursive(points, left, mid);
    pair<Point, Point> rightPair = closestPairRecursive(points, mid + 1, right);
    
    double leftDist = leftPair.first.distanceTo(leftPair.second);
    double rightDist = rightPair.first.distanceTo(rightPair.second);
    
    double minDist = min(leftDist, rightDist);
    pair<Point, Point> minPair = (leftDist < rightDist) ? leftPair : rightPair;
    
    // Check for closer pairs across the dividing line
    vector<Point> strip;
    for (int i = left; i <= right; i++) {
        if (abs(points[i].x - midPoint.x) < minDist) {
            strip.push_back(points[i]);
        }
    }
    
    // Sort strip by y-coordinate
    sort(strip.begin(), strip.end(), [](const Point& a, const Point& b) {
        return a.y < b.y;
    });
    
    // Check for closer pairs in strip
    for (int i = 0; i < strip.size(); i++) {
        for (int j = i + 1; j < strip.size() && (strip[j].y - strip[i].y) < minDist; j++) {
            double dist = strip[i].distanceTo(strip[j]);
            if (dist < minDist) {
                minDist = dist;
                minPair = {strip[i], strip[j]};
            }
        }
    }
    
    return minPair;
}

pair<Point, Point> bruteForceClosestPair(vector<Point>& points, int left, int right) {
    double minDist = INFINITY;
    pair<Point, Point> minPair;
    
    for (int i = left; i <= right; i++) {
        for (int j = i + 1; j <= right; j++) {
            double dist = points[i].distanceTo(points[j]);
            if (dist < minDist) {
                minDist = dist;
                minPair = {points[i], points[j]};
            }
        }
    }
    
    return minPair;
}

Line Segment Intersection​

Bentley-Ottmann Algorithm​

struct Segment {
    Point p1, p2;
    int id;
    
    Segment(const Point& p1, const Point& p2, int id) : p1(p1), p2(p2), id(id) {}
    
    double getYAtX(double x) const {
        if (abs(p2.x - p1.x) < 1e-9) return p1.y;
        return p1.y + (p2.y - p1.y) * (x - p1.x) / (p2.x - p1.x);
    }
};

// Bentley-Ottmann algorithm for line segment intersection
vector<Point> bentleyOttmannIntersection(vector<Segment>& segments) {
    vector<Event> events;
    vector<Point> intersections;
    
    // Create events for each segment
    for (int i = 0; i < segments.size(); i++) {
        double x1 = segments[i].p1.x;
        double x2 = segments[i].p2.x;
        
        if (x1 > x2) {
            swap(x1, x2);
        }
        
        events.push_back(Event(x1, 0, i)); // Start
        events.push_back(Event(x2, 1, i)); // End
    }
    
    sort(events.begin(), events.end());
    
    set<Segment> activeSegments;
    
    for (const Event& event : events) {
        if (event.type == 0) {
            // Start event
            Segment current = segments[event.id];
            activeSegments.insert(current);
            
            // Check intersections with adjacent segments
            auto it = activeSegments.find(current);
            if (it != activeSegments.begin()) {
                auto prev = it;
                prev--;
                Point intersection = lineSegmentIntersection(current, *prev);
                if (intersection.x != INFINITY) {
                    intersections.push_back(intersection);
                }
            }
            
            auto next = it;
            next++;
            if (next != activeSegments.end()) {
                Point intersection = lineSegmentIntersection(current, *next);
                if (intersection.x != INFINITY) {
                    intersections.push_back(intersection);
                }
            }
        } else {
            // End event
            Segment current = segments[event.id];
            auto it = activeSegments.find(current);
            
            if (it != activeSegments.begin() && it != activeSegments.end()) {
                auto prev = it;
                prev--;
                auto next = it;
                next++;
                
                if (next != activeSegments.end()) {
                    Point intersection = lineSegmentIntersection(*prev, *next);
                    if (intersection.x != INFINITY) {
                        intersections.push_back(intersection);
                    }
                }
            }
            
            activeSegments.erase(current);
        }
    }
    
    return intersections;
}

Polygon Operations​

Point in Polygon​

// Check if point is inside polygon using ray casting
bool pointInPolygon(const Point& point, const vector<Point>& polygon) {
    int n = polygon.size();
    bool inside = false;
    
    for (int i = 0, j = n - 1; i < n; j = i++) {
        if (((polygon[i].y > point.y) != (polygon[j].y > point.y)) &&
            (point.x < (polygon[j].x - polygon[i].x) * (point.y - polygon[i].y) / 
             (polygon[j].y - polygon[i].y) + polygon[i].x)) {
            inside = !inside;
        }
    }
    
    return inside;
}

// Check if point is on polygon boundary
bool pointOnPolygonBoundary(const Point& point, const vector<Point>& polygon) {
    int n = polygon.size();
    
    for (int i = 0; i < n; i++) {
        int j = (i + 1) % n;
        LineSegment edge(polygon[i], polygon[j]);
        
        if (edge.containsPoint(point)) {
            return true;
        }
    }
    
    return false;
}

Polygon Intersection​

// Check if two polygons intersect
bool polygonsIntersect(const vector<Point>& poly1, const vector<Point>& poly2) {
    // Check if any vertex of poly1 is inside poly2
    for (const Point& p : poly1) {
        if (pointInPolygon(p, poly2)) {
            return true;
        }
    }
    
    // Check if any vertex of poly2 is inside poly1
    for (const Point& p : poly2) {
        if (pointInPolygon(p, poly1)) {
            return true;
        }
    }
    
    // Check if any edges intersect
    for (int i = 0; i < poly1.size(); i++) {
        int j = (i + 1) % poly1.size();
        LineSegment edge1(poly1[i], poly1[j]);
        
        for (int k = 0; k < poly2.size(); k++) {
            int l = (k + 1) % poly2.size();
            LineSegment edge2(poly2[k], poly2[l]);
            
            if (edge1.intersects(edge2)) {
                return true;
            }
        }
    }
    
    return false;
}

Polygon Triangulation​

// Triangulate polygon using ear clipping
vector<vector<Point>> triangulatePolygon(vector<Point>& polygon) {
    vector<vector<Point>> triangles;
    
    while (polygon.size() > 3) {
        bool earFound = false;
        
        for (int i = 0; i < polygon.size(); i++) {
            int prev = (i - 1 + polygon.size()) % polygon.size();
            int next = (i + 1) % polygon.size();
            
            if (isEar(polygon, prev, i, next)) {
                // Add triangle
                triangles.push_back({polygon[prev], polygon[i], polygon[next]});
                
                // Remove ear vertex
                polygon.erase(polygon.begin() + i);
                earFound = true;
                break;
            }
        }
        
        if (!earFound) {
            break; // No ear found, polygon might be degenerate
        }
    }
    
    // Add final triangle
    if (polygon.size() == 3) {
        triangles.push_back(polygon);
    }
    
    return triangles;
}

bool isEar(const vector<Point>& polygon, int prev, int curr, int next) {
    Point p1 = polygon[prev];
    Point p2 = polygon[curr];
    Point p3 = polygon[next];
    
    // Check if triangle is convex
    if ((p2 - p1).cross(p3 - p1) <= 0) {
        return false;
    }
    
    // Check if any other vertex is inside the triangle
    for (int i = 0; i < polygon.size(); i++) {
        if (i == prev || i == curr || i == next) continue;
        
        Point p = polygon[i];
        if (pointInTriangle(p, p1, p2, p3)) {
            return false;
        }
    }
    
    return true;
}

bool pointInTriangle(const Point& p, const Point& a, const Point& b, const Point& c) {
    double denom = (b.y - c.y) * (a.x - c.x) + (c.x - b.x) * (a.y - c.y);
    double alpha = ((b.y - c.y) * (p.x - c.x) + (c.x - b.x) * (p.y - c.y)) / denom;
    double beta = ((c.y - a.y) * (p.x - c.x) + (a.x - c.x) * (p.y - c.y)) / denom;
    double gamma = 1 - alpha - beta;
    
    return alpha >= 0 && beta >= 0 && gamma >= 0;
}

3D Geometry Basics​

3D Point Operations​

struct Point3D {
    double x, y, z;
    
    Point3D(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}
    
    Point3D operator+(const Point3D& other) const {
        return Point3D(x + other.x, y + other.y, z + other.z);
    }
    
    Point3D operator-(const Point3D& other) const {
        return Point3D(x - other.x, y - other.y, z - other.z);
    }
    
    Point3D operator*(double scalar) const {
        return Point3D(x * scalar, y * scalar, z * scalar);
    }
    
    // Dot product
    double dot(const Point3D& other) const {
        return x * other.x + y * other.y + z * other.z;
    }
    
    // Cross product
    Point3D cross(const Point3D& other) const {
        return Point3D(
            y * other.z - z * other.y,
            z * other.x - x * other.z,
            x * other.y - y * other.x
        );
    }
    
    // Distance to another point
    double distanceTo(const Point3D& other) const {
        return sqrt((x - other.x) * (x - other.x) + 
                   (y - other.y) * (y - other.y) + 
                   (z - other.z) * (z - other.z));
    }
    
    // Magnitude
    double magnitude() const {
        return sqrt(x * x + y * y + z * z);
    }
};

3D Plane Operations​

struct Plane {
    double a, b, c, d; // ax + by + cz + d = 0
    
    Plane(double a = 0, double b = 0, double c = 0, double d = 0) 
        : a(a), b(b), c(c), d(d) {}
    
    // Plane from three points
    Plane(const Point3D& p1, const Point3D& p2, const Point3D& p3) {
        Point3D v1 = p2 - p1;
        Point3D v2 = p3 - p1;
        Point3D normal = v1.cross(v2);
        
        a = normal.x;
        b = normal.y;
        c = normal.z;
        d = -(a * p1.x + b * p1.y + c * p1.z);
    }
    
    // Distance from point to plane
    double distanceToPoint(const Point3D& p) const {
        return abs(a * p.x + b * p.y + c * p.z + d) / sqrt(a * a + b * b + c * c);
    }
    
    // Check if point is on plane
    bool containsPoint(const Point3D& p) const {
        return abs(a * p.x + b * p.y + c * p.z + d) < 1e-9;
    }
};

Performance Analysis​

Time Complexity​

  • Sweep line intersection: O(n log n + k) where k is number of intersections
  • Closest pair: O(n log n)
  • Point in polygon: O(n)
  • Polygon intersection: O(n × m)
  • 3D operations: O(1) for basic operations

Space Complexity​

  • Sweep line: O(n)
  • Closest pair: O(n)
  • Polygon operations: O(n)
  • 3D operations: O(1)

Common Patterns​

  1. Sweep line for intersection problems
  2. Divide and conquer for optimization problems
  3. Event-driven processing for geometric algorithms
  4. Cross product for orientation and area
  5. Ray casting for point-in-polygon tests

Applications​

  • Computer graphics: Rendering and collision detection
  • GIS: Geographic information systems
  • Robotics: Path planning and obstacle avoidance
  • Game development: Physics and collision detection
  • Computational geometry: Algorithm design and analysis