15种经典几何算法完整C++实现

8种经典几何算法完整C++实现

以下是8种经典计算几何算法的完整C++实现，包括详细注释、时间/空间复杂度分析和测试代码。

完整实现代码

/**
 * 8种经典几何算法完整C++实现
 * 编译: g++ -std=c++17 -O2 geometry_algorithms.cpp -o geometry
 */

#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <limits>
#include <stack>
#include <queue>
#include <set>
#include <iomanip>
#include <cfloat>
#include <complex>
#include <cassert>
using namespace std;

// ==================== 基本几何结构 ====================
const double EPS = 1e-9;
const double PI = acos(-1.0);

// 浮点数比较
int dcmp(double x) {
    if (fabs(x) < EPS) return 0;
    return x < 0 ? -1 : 1;
}

// 二维点/向量
struct Point {
    double x, y;
    
    Point() : x(0), y(0) {}
    Point(double x, double y) : x(x), y(y) {}
    
    // 向量加法
    Point operator+(const Point& p) const {
        return Point(x + p.x, y + p.y);
    }
    
    // 向量减法
    Point operator-(const Point& p) const {
        return Point(x - p.x, y - p.y);
    }
    
    // 标量乘法
    Point operator*(double s) const {
        return Point(x * s, y * s);
    }
    
    // 标量除法
    Point operator/(double s) const {
        return Point(x / s, y / s);
    }
    
    // 点积
    double dot(const Point& p) const {
        return x * p.x + y * p.y;
    }
    
    // 叉积
    double cross(const Point& p) const {
        return x * p.y - y * p.x;
    }
    
    // 长度平方
    double length2() const {
        return dot(*this);
    }
    
    // 长度
    double length() const {
        return sqrt(length2());
    }
    
    // 单位向量
    Point unit() const {
        return *this / length();
    }
    
    // 旋转角度（弧度）
    Point rotate(double theta) const {
        double c = cos(theta), s = sin(theta);
        return Point(x * c - y * s, x * s + y * c);
    }
    
    // 比较运算符
    bool operator<(const Point& p) const {
        if (dcmp(x - p.x) != 0) return x < p.x;
        return y < p.y;
    }
    
    bool operator==(const Point& p) const {
        return dcmp(x - p.x) == 0 && dcmp(y - p.y) == 0;
    }
    
    // 输出
    friend ostream& operator<<(ostream& os, const Point& p) {
        os << "(" << p.x << ", " << p.y << ")";
        return os;
    }
};

// 线段
struct Segment {
    Point p1, p2;
    
    Segment() {}
    Segment(const Point& p1, const Point& p2) : p1(p1), p2(p2) {}
    
    // 线段长度
    double length() const {
        return (p2 - p1).length();
    }
};

// 直线: ax + by + c = 0
struct Line {
    double a, b, c;
    
    Line() : a(0), b(0), c(0) {}
    Line(double a, double b, double c) : a(a), b(b), c(c) {}
    
    // 从两点构造直线
    Line(const Point& p1, const Point& p2) {
        a = p2.y - p1.y;
        b = p1.x - p2.x;
        c = -(a * p1.x + b * p1.y);
        
        // 归一化
        double len = sqrt(a * a + b * b);
        if (dcmp(len) != 0) {
            a /= len;
            b /= len;
            c /= len;
        }
    }
    
    // 点到直线的距离
    double distance(const Point& p) const {
        return fabs(a * p.x + b * p.y + c) / sqrt(a * a + b * b);
    }
};

// 多边形
using Polygon = vector<Point>;

// 圆形
struct Circle {
    Point center;
    double radius;
    
    Circle() : radius(0) {}
    Circle(const Point& center, double radius) : center(center), radius(radius) {}
};

// ==================== 1. 凸包算法（Graham Scan） ====================
/**
 * 时间复杂度: O(n log n)
 * 空间复杂度: O(n)
 * 应用: 凸多边形计算、碰撞检测
 */
class ConvexHull {
private:
    // 极角排序比较函数
    Point pivot;  // 基准点（最低点）
    
    // 叉积判断方向
    int orientation(const Point& a, const Point& b, const Point& c) {
        double cross = (b - a).cross(c - a);
        if (dcmp(cross) == 0) return 0;  // 共线
        return cross > 0 ? 1 : 2;  // 逆时针或顺时针
    }
    
    // 极角排序比较函数
    bool polarCompare(const Point& p1, const Point& p2) {
        int o = orientation(pivot, p1, p2);
        if (o == 0) {
            return (pivot - p1).length2() < (pivot - p2).length2();
        }
        return o == 2;  // 按逆时针排序
    }
    
public:
    // Graham Scan算法
    Polygon grahamScan(vector<Point> points) {
        int n = points.size();
        if (n < 3) {
            if (n == 2 && points[0] == points[1]) {
                points.pop_back();
            }
            return points;
        }
        
        // 找到y坐标最小的点（如果相同，取x最小的）
        int minIdx = 0;
        for (int i = 1; i < n; i++) {
            if (dcmp(points[i].y - points[minIdx].y) < 0 ||
                (dcmp(points[i].y - points[minIdx].y) == 0 &&
                 dcmp(points[i].x - points[minIdx].x) < 0)) {
                minIdx = i;
            }
        }
        
        // 将基准点放到第一个位置
        swap(points[0], points[minIdx]);
        pivot = points[0];
        
        // 按极角排序
        sort(points.begin() + 1, points.end(), 
             [this](const Point& a, const Point& b) {
                 return polarCompare(a, b);
             });
        
        // 移除共线点（保留最远的）
        int m = 1;
        for (int i = 1; i < n; i++) {
            while (i < n-1 && orientation(pivot, points[i], points[i+1]) == 0) {
                i++;
            }
            points[m++] = points[i];
        }
        
        if (m < 3) {
            points.resize(m);
            return points;
        }
        
        // Graham Scan
        stack<Point> stk;
        stk.push(points[0]);
        stk.push(points[1]);
        stk.push(points[2]);
        
        for (int i = 3; i < m; i++) {
            while (stk.size() >= 2) {
                Point top = stk.top();
                stk.pop();
                Point nextTop = stk.top();
                
                if (orientation(nextTop, top, points[i]) == 2) {
                    stk.push(top);
                    break;
                }
            }
            stk.push(points[i]);
        }
        
        // 从栈中获取凸包
        Polygon hull;
        while (!stk.empty()) {
            hull.push_back(stk.top());
            stk.pop();
        }
        reverse(hull.begin(), hull.end());
        
        return hull;
    }
    
    // Jarvis March（礼品包装）算法
    Polygon jarvisMarch(vector<Point> points) {
        int n = points.size();
        if (n < 3) {
            if (n == 2 && points[0] == points[1]) {
                points.pop_back();
            }
            return points;
        }
        
        // 找到最左边的点
        int leftmost = 0;
        for (int i = 1; i < n; i++) {
            if (dcmp(points[i].x - points[leftmost].x) < 0 ||
                (dcmp(points[i].x - points[leftmost].x) == 0 &&
                 dcmp(points[i].y - points[leftmost].y) < 0)) {
                leftmost = i;
            }
        }
        
        Polygon hull;
        int p = leftmost, q;
        
        do {
            hull.push_back(points[p]);
            q = (p + 1) % n;
            
            for (int i = 0; i < n; i++) {
                if (orientation(points[p], points[i], points[q]) == 1) {
                    q = i;
                }
            }
            p = q;
        } while (p != leftmost);
        
        return hull;
    }
    
    // 计算凸包面积
    double convexHullArea(const Polygon& hull) {
        if (hull.size() < 3) return 0;
        
        double area = 0;
        int n = hull.size();
        Point origin(0, 0);
        
        for (int i = 0; i < n; i++) {
            int j = (i + 1) % n;
            area += hull[i].cross(hull[j]);
        }
        
        return fabs(area) / 2.0;
    }
    
    // 判断点是否在凸多边形内
    bool pointInConvexPolygon(const Point& p, const Polygon& hull) {
        int n = hull.size();
        if (n < 3) return false;
        
        // 通过叉积判断方向
        for (int i = 0; i < n; i++) {
            int j = (i + 1) % n;
            if (orientation(hull[i], hull[j], p) == 1) {
                return false;
            }
        }
        return true;
    }
};

// ==================== 2. 旋转卡壳算法 ====================
/**
 * 时间复杂度: O(n)
 * 空间复杂度: O(1)
 * 应用: 凸包直径、宽度、最小面积外接矩形
 */
class RotatingCalipers {
private:
    // 叉积计算面积
    double cross(const Point& o, const Point& a, const Point& b) {
        return (a - o).cross(b - o);
    }
    
    // 点到线段距离
    double pointToSegmentDistance(const Point& p, const Point& a, const Point& b) {
        Point ab = b - a;
        Point ap = p - a;
        
        double t = ap.dot(ab) / ab.length2();
        if (t < 0) return (p - a).length();
        if (t > 1) return (p - b).length();
        
        Point projection = a + ab * t;
        return (p - projection).length();
    }
    
public:
    // 计算凸包直径（最远点对）
    double convexHullDiameter(const Polygon& hull) {
        int n = hull.size();
        if (n < 2) return 0;
        if (n == 2) return (hull[0] - hull[1]).length();
        
        int j = 1;
        double maxDist = 0;
        
        for (int i = 0; i < n; i++) {
            int nextI = (i + 1) % n;
            
            while (true) {
                int nextJ = (j + 1) % n;
                
                // 计算平行四边形的面积
                double area1 = cross(hull[i], hull[nextI], hull[j]);
                double area2 = cross(hull[i], hull[nextI], hull[nextJ]);
                
                if (dcmp(area2 - area1) > 0) {
                    j = nextJ;
                } else {
                    break;
                }
            }
            
            // 更新最大距离
            maxDist = max(maxDist, (hull[i] - hull[j]).length());
            maxDist = max(maxDist, (hull[i] - hull[(j + 1) % n]).length());
        }
        
        return maxDist;
    }
    
    // 计算凸包宽度
    double convexHullWidth(const Polygon& hull) {
        int n = hull.size();
        if (n < 2) return 0;
        
        int j = 1;
        double minWidth = DBL_MAX;
        
        for (int i = 0; i < n; i++) {
            int nextI = (i + 1) % n;
            
            while (true) {
                int nextJ = (j + 1) % n;
                
                double area1 = cross(hull[i], hull[nextI], hull[j]);
                double area2 = cross(hull[i], hull[nextI], hull[nextJ]);
                
                if (dcmp(area2 - area1) > 0) {
                    j = nextJ;
                } else {
                    break;
                }
            }
            
            // 计算点i到线段(j, j+1)的距离
            double width = pointToSegmentDistance(hull[i], hull[j], hull[(j + 1) % n]);
            minWidth = min(minWidth, width);
        }
        
        return minWidth;
    }
    
    // 计算最小面积外接矩形
    struct MinAreaRect {
        double area;
        double width;
        double height;
        vector<Point> rect;  // 矩形的四个顶点
    };
    
    MinAreaRect minimumAreaBoundingBox(const Polygon& hull) {
        int n = hull.size();
        MinAreaRect result;
        result.area = DBL_MAX;
        
        if (n < 3) {
            result.area = 0;
            return result;
        }
        
        int j = 1, k = 1, l = 1;
        
        for (int i = 0; i < n; i++) {
            int nextI = (i + 1) % n;
            
            // 寻找最远的点（直径方向）
            while (dcmp(cross(hull[i], hull[nextI], hull[j]) - 
                       cross(hull[i], hull[nextI], hull[(j + 1) % n])) < 0) {
                j = (j + 1) % n;
            }
            
            // 寻找最右边的点
            if (i == 0) k = j;
            while (dcmp((hull[nextI] - hull[i]).dot(hull[(k + 1) % n] - hull[k])) > 0) {
                k = (k + 1) % n;
            }
            
            // 寻找最左边的点
            if (i == 0) l = j;
            while (dcmp((hull[nextI] - hull[i]).dot(hull[(l + 1) % n] - hull[l])) < 0) {
                l = (l + 1) % n;
            }
            
            // 计算矩形
            Point dir = hull[nextI] - hull[i];
            double len = dir.length();
            dir = dir / len;  // 单位方向向量
            
            // 计算投影
            double minProj = dir.dot(hull[l] - hull[i]);
            double maxProj = dir.dot(hull[k] - hull[i]);
            double width = maxProj - minProj;
            
            // 计算高度
            double height = cross(hull[i], hull[nextI], hull[j]) / len;
            if (height < 0) height = -height;
            
            double area = width * height;
            if (dcmp(area - result.area) < 0) {
                result.area = area;
                result.width = width;
                result.height = height;
                
                // 计算矩形顶点
                Point perp(-dir.y, dir.x);  // 垂直方向
                result.rect.clear();
                result.rect.push_back(hull[i] + dir * minProj + perp * 0);
                result.rect.push_back(hull[i] + dir * maxProj + perp * 0);
                result.rect.push_back(hull[i] + dir * maxProj + perp * height);
                result.rect.push_back(hull[i] + dir * minProj + perp * height);
            }
        }
        
        return result;
    }
};

// ==================== 3. 最近点对算法 ====================
/**
 * 时间复杂度: O(n log n)
 * 空间复杂度: O(n)
 * 应用: 聚类分析、碰撞检测
 */
class ClosestPair {
private:
    // 按x坐标排序的比较函数
    static bool compareX(const Point& a, const Point& b) {
        if (dcmp(a.x - b.x) != 0) return a.x < b.x;
        return a.y < b.y;
    }
    
    // 按y坐标排序的比较函数
    static bool compareY(const Point& a, const Point& b) {
        if (dcmp(a.y - b.y) != 0) return a.y < b.y;
        return a.x < b.x;
    }
    
    // 计算两点距离
    double dist(const Point& a, const Point& b) {
        return (a - b).length();
    }
    
    // 暴力求解
    pair<double, pair<Point, Point>> bruteForce(vector<Point>& points, int left, int right) {
        double minDist = DBL_MAX;
        pair<Point, Point> minPair;
        
        for (int i = left; i <= right; i++) {
            for (int j = i + 1; j <= right; j++) {
                double d = dist(points[i], points[j]);
                if (d < minDist) {
                    minDist = d;
                    minPair = {points[i], points[j]};
                }
            }
        }
        
        return {minDist, minPair};
    }
    
    // 在条带中寻找最近点
    pair<double, pair<Point, Point>> stripClosest(vector<Point>& strip, double d) {
        sort(strip.begin(), strip.end(), compareY);
        double minDist = d;
        pair<Point, Point> minPair;
        int n = strip.size();
        
        for (int i = 0; i < n; i++) {
            for (int j = i + 1; j < n && (strip[j].y - strip[i].y) < minDist; j++) {
                double d = dist(strip[i], strip[j]);
                if (d < minDist) {
                    minDist = d;
                    minPair = {strip[i], strip[j]};
                }
            }
        }
        
        return {minDist, minPair};
    }
    
    // 分治递归
    pair<double, pair<Point, Point>> closestUtil(vector<Point>& points, int left, int right) {
        if (right - left + 1 <= 3) {
            return bruteForce(points, left, right);
        }
        
        int mid = (left + right) / 2;
        Point midPoint = points[mid];
        
        auto leftResult = closestUtil(points, left, mid);
        auto rightResult = closestUtil(points, mid + 1, right);
        
        auto minResult = (leftResult.first < rightResult.first) ? leftResult : rightResult;
        double d = minResult.first;
        
        // 收集条带内的点
        vector<Point> strip;
        for (int i = left; i <= right; i++) {
            if (fabs(points[i].x - midPoint.x) < d) {
                strip.push_back(points[i]);
            }
        }
        
        auto stripResult = stripClosest(strip, d);
        if (stripResult.first < minResult.first) {
            return stripResult;
        }
        return minResult;
    }
    
public:
    // 分治法求最近点对
    pair<double, pair<Point, Point>> closestPair(vector<Point> points) {
        int n = points.size();
        if (n < 2) {
            return {DBL_MAX, {}};
        }
        
        // 按x坐标排序
        sort(points.begin(), points.end(), compareX);
        
        return closestUtil(points, 0, n - 1);
    }
    
    // 计算所有点对的最小距离
    double minimumDistance(vector<Point>& points) {
        auto result = closestPair(points);
        return result.first;
    }
};

// ==================== 4. 多边形面积算法 ====================
/**
 * 时间复杂度: O(n)
 * 空间复杂度: O(1)
 * 应用: 面积计算、质心计算
 */
class PolygonArea {
public:
    // 鞋带公式计算多边形面积
    double shoelaceFormula(const Polygon& polygon) {
        int n = polygon.size();
        if (n < 3) return 0;
        
        double area = 0;
        for (int i = 0; i < n; i++) {
            int j = (i + 1) % n;
            area += polygon[i].x * polygon[j].y;
            area -= polygon[j].x * polygon[i].y;
        }
        
        return fabs(area) / 2.0;
    }
    
    // 计算多边形周长
    double polygonPerimeter(const Polygon& polygon) {
        int n = polygon.size();
        if (n < 2) return 0;
        
        double perimeter = 0;
        for (int i = 0; i < n; i++) {
            int j = (i + 1) % n;
            perimeter += (polygon[j] - polygon[i]).length();
        }
        
        return perimeter;
    }
    
    // 计算多边形质心
    Point polygonCentroid(const Polygon& polygon) {
        int n = polygon.size();
        if (n == 0) return Point();
        if (n == 1) return polygon[0];
        
        double area = 0;
        double cx = 0, cy = 0;
        
        for (int i = 0; i < n; i++) {
            int j = (i + 1) % n;
            double cross = polygon[i].x * polygon[j].y - polygon[j].x * polygon[i].y;
            area += cross;
            cx += (polygon[i].x + polygon[j].x) * cross;
            cy += (polygon[i].y + polygon[j].y) * cross;
        }
        
        area /= 2.0;
        if (dcmp(area) == 0) {
            // 多边形退化
            double sumX = 0, sumY = 0;
            for (const auto& p : polygon) {
                sumX += p.x;
                sumY += p.y;
            }
            return Point(sumX / n, sumY / n);
        }
        
        cx /= (6.0 * area);
        cy /= (6.0 * area);
        
        return Point(cx, cy);
    }
    
    // 判断多边形方向（逆时针为正）
    int polygonOrientation(const Polygon& polygon) {
        int n = polygon.size();
        if (n < 3) return 0;
        
        // 找到最左下角的点
        int minIdx = 0;
        for (int i = 1; i < n; i++) {
            if (dcmp(polygon[i].y - polygon[minIdx].y) < 0 ||
                (dcmp(polygon[i].y - polygon[minIdx].y) == 0 &&
                 dcmp(polygon[i].x - polygon[minIdx].x) < 0)) {
                minIdx = i;
            }
        }
        
        int prev = (minIdx - 1 + n) % n;
        int next = (minIdx + 1) % n;
        
        double cross = (polygon[minIdx] - polygon[prev]).cross(polygon[next] - polygon[minIdx]);
        return dcmp(cross);
    }
};

// ==================== 5. 点在多边形内检测 ====================
/**
 * 时间复杂度: O(n)
 * 空间复杂度: O(1)
 * 应用: 碰撞检测、区域包含测试
 */
class PointInPolygon {
private:
    // 点在线段上的判断
    bool onSegment(const Point& p, const Point& a, const Point& b) {
        if (dcmp((a - p).cross(b - p)) != 0) return false;
        return dcmp(min(a.x, b.x) - p.x) <= 0 &&
               dcmp(p.x - max(a.x, b.x)) <= 0 &&
               dcmp(min(a.y, b.y) - p.y) <= 0 &&
               dcmp(p.y - max(a.y, b.y)) <= 0;
    }
    
    // 射线与线段相交检测
    bool rayIntersectsSegment(const Point& p, const Point& a, const Point& b) {
        // 确保A在B下方
        if (a.y > b.y) return rayIntersectsSegment(p, b, a);
        
        // 点在射线下方或与A、B共线
        if (dcmp(p.y - a.y) <= 0 || dcmp(p.y - b.y) > 0) return false;
        
        // 计算交点X坐标
        double t = (p.y - a.y) / (b.y - a.y);
        double x = a.x + t * (b.x - a.x);
        
        return dcmp(x - p.x) > 0;  // 交点在右侧
    }
    
public:
    // 射线法（适用于任意多边形）
    int rayCasting(const Point& p, const Polygon& polygon) {
        int n = polygon.size();
        if (n < 3) return 0;
        
        int count = 0;
        
        for (int i = 0; i < n; i++) {
            int j = (i + 1) % n;
            
            // 检查点是否在线段上
            if (onSegment(p, polygon[i], polygon[j])) {
                return 0;  // 在边界上
            }
            
            if (rayIntersectsSegment(p, polygon[i], polygon[j])) {
                count++;
            }
        }
        
        return (count % 2 == 1) ? 1 : -1;  // 1: 内部, -1: 外部
    }
    
    // 转角法（更精确）
    int windingNumber(const Point& p, const Polygon& polygon) {
        int n = polygon.size();
        if (n < 3) return 0;
        
        double wn = 0;  // 环绕数
        
        for (int i = 0; i < n; i++) {
            int j = (i + 1) % n;
            
            // 检查点是否在线段上
            if (onSegment(p, polygon[i], polygon[j])) {
                return 0;  // 在边界上
            }
            
            double y1 = polygon[i].y - p.y;
            double y2 = polygon[j].y - p.y;
            
            if (dcmp(y1) <= 0) {
                if (dcmp(y2) > 0) {
                    if (dcmp((polygon[i] - p).cross(polygon[j] - p)) > 0) {
                        wn += 1;  // 向上穿过
                    }
                }
            } else {
                if (dcmp(y2) <= 0) {
                    if (dcmp((polygon[i] - p).cross(polygon[j] - p)) < 0) {
                        wn -= 1;  // 向下穿过
                    }
                }
            }
        }
        
        return dcmp(wn) == 0 ? -1 : 1;
    }
    
    // 凸多边形中的点检测（二分法）
    bool pointInConvexPolygonBinary(const Point& p, const Polygon& convex) {
        int n = convex.size();
        if (n < 3) return false;
        
        // 检查是否在第一条边的外侧
        if (dcmp((convex[1] - convex[0]).cross(p - convex[0])) < 0) {
            return false;
        }
        
        // 检查是否在最后一条边的外侧
        if (dcmp((convex[n-1] - convex[0]).cross(p - convex[0])) > 0) {
            return false;
        }
        
        // 二分查找
        int left = 2, right = n - 1;
        int pos = 1;
        
        while (left <= right) {
            int mid = (left + right) / 2;
            if (dcmp((convex[mid] - convex[0]).cross(p - convex[0])) >= 0) {
                pos = mid;
                left = mid + 1;
            } else {
                right = mid - 1;
            }
        }
        
        // 检查点是否在三角形内
        Point v0 = convex[pos] - convex[0];
        Point v1 = convex[pos + 1] - convex[0];
        Point v2 = p - convex[0];
        
        double d00 = v0.dot(v0);
        double d01 = v0.dot(v1);
        double d11 = v1.dot(v1);
        double d20 = v2.dot(v0);
        double d21 = v2.dot(v1);
        
        double denom = d00 * d11 - d01 * d01;
        double u = (d11 * d20 - d01 * d21) / denom;
        double v = (d00 * d21 - d01 * d20) / denom;
        
        return dcmp(u) >= 0 && dcmp(v) >= 0 && dcmp(u + v - 1) <= 0;
    }
};

// ==================== 6. 线段相交检测 ====================
/**
 * 时间复杂度: O(1)
 * 空间复杂度: O(1)
 * 应用: 碰撞检测、裁剪算法
 */
class SegmentIntersection {
private:
    // 计算点p在线段ab上的投影
    Point projectPointToLine(const Point& p, const Point& a, const Point& b) {
        Point ab = b - a;
        double t = (p - a).dot(ab) / ab.length2();
        return a + ab * t;
    }
    
public:
    // 判断点是否在线段上
    bool onSegment(const Point& p, const Point& a, const Point& b) {
        if (dcmp((a - p).cross(b - p)) != 0) return false;
        return dcmp(min(a.x, b.x) - p.x) <= 0 &&
               dcmp(p.x - max(a.x, b.x)) <= 0 &&
               dcmp(min(a.y, b.y) - p.y) <= 0 &&
               dcmp(p.y - max(a.y, b.y)) <= 0;
    }
    
    // 判断线段是否相交
    bool segmentsIntersect(const Point& a1, const Point& a2, 
                          const Point& b1, const Point& b2) {
        double d1 = dcmp((b2 - b1).cross(a1 - b1));
        double d2 = dcmp((b2 - b1).cross(a2 - b1));
        double d3 = dcmp((a2 - a1).cross(b1 - a1));
        double d4 = dcmp((a2 - a1).cross(b2 - a1));
        
        if (d1 * d2 < 0 && d3 * d4 < 0) return true;
        
        if (d1 == 0 && onSegment(a1, b1, b2)) return true;
        if (d2 == 0 && onSegment(a2, b1, b2)) return true;
        if (d3 == 0 && onSegment(b1, a1, a2)) return true;
        if (d4 == 0 && onSegment(b2, a1, a2)) return true;
        
        return false;
    }
    
    // 计算线段交点
    bool segmentIntersectionPoint(const Point& a1, const Point& a2,
                                 const Point& b1, const Point& b2,
                                 Point& result) {
        if (!segmentsIntersect(a1, a2, b1, b2)) {
            return false;
        }
        
        Point da = a2 - a1;
        Point db = b2 - b1;
        Point dc = b1 - a1;
        
        double denom = da.cross(db);
        if (dcmp(denom) == 0) {
            // 平行或共线
            if (onSegment(a1, b1, b2)) {
                result = a1;
                return true;
            }
            if (onSegment(a2, b1, b2)) {
                result = a2;
                return true;
            }
            if (onSegment(b1, a1, a2)) {
                result = b1;
                return true;
            }
            if (onSegment(b2, a1, a2)) {
                result = b2;
                return true;
            }
            return false;
        }
        
        double t = dc.cross(db) / denom;
        result = a1 + da * t;
        return true;
    }
    
    // 点到线段的最短距离
    double pointToSegmentDistance(const Point& p, const Point& a, const Point& b) {
        Point ab = b - a;
        Point ap = p - a;
        
        double t = ap.dot(ab) / ab.length2();
        if (t < 0) return (p - a).length();
        if (t > 1) return (p - b).length();
        
        Point projection = a + ab * t;
        return (p - projection).length();
    }
    
    // 线段到线段的最短距离
    double segmentToSegmentDistance(const Point& a1, const Point& a2,
                                   const Point& b1, const Point& b2) {
        if (segmentsIntersect(a1, a2, b1, b2)) {
            return 0;
        }
        
        double d1 = pointToSegmentDistance(a1, b1, b2);
        double d2 = pointToSegmentDistance(a2, b1, b2);
        double d3 = pointToSegmentDistance(b1, a1, a2);
        double d4 = pointToSegmentDistance(b2, a1, a2);
        
        return min(min(d1, d2), min(d3, d4));
    }
};

// ==================== 7. 半平面交算法 ====================
/**
 * 时间复杂度: O(n log n)
 * 空间复杂度: O(n)
 * 应用: 线性规划、可见性计算
 */
class HalfplaneIntersection {
private:
    struct Halfplane {
        Point p;      // 直线上一点
        Point dir;    // 方向向量
        double angle; // 极角
        
        Halfplane() {}
        Halfplane(const Point& p1, const Point& p2) : p(p1) {
            dir = p2 - p1;
            angle = atan2(dir.y, dir.x);
        }
        
        // 比较极角
        bool operator<(const Halfplane& other) const {
            if (dcmp(angle - other.angle) != 0) {
                return angle < other.angle;
            }
            return dcmp(dir.cross(other.p - p)) < 0;
        }
        
        // 点是否在半平面内（左侧）
        bool contains(const Point& q) const {
            return dcmp(dir.cross(q - p)) >= 0;
        }
    };
    
    // 计算两直线交点
    Point intersectLines(const Point& a, const Point& da,
                        const Point& b, const Point& db) {
        Point dc = b - a;
        double t = dc.cross(db) / da.cross(db);
        return a + da * t;
    }
    
public:
    // 计算半平面交
    Polygon halfplaneIntersection(vector<Halfplane>& halfplanes) {
        // 排序并去重
        sort(halfplanes.begin(), halfplanes.end());
        
        vector<Halfplane> uniqueHalfplanes;
        for (int i = 0; i < halfplanes.size(); i++) {
            if (i == 0 || dcmp(halfplanes[i].angle - halfplanes[i-1].angle) != 0) {
                uniqueHalfplanes.push_back(halfplanes[i]);
            }
        }
        
        int n = uniqueHalfplanes.size();
        if (n < 3) {
            // 半平面数太少，返回空
            return {};
        }
        
        vector<Halfplane> deque(n * 2);
        int first = 0, last = 0;
        
        for (int i = 0; i < n; i++) {
            // 删除无用的半平面
            while (last - first > 1 && 
                   !uniqueHalfplanes[i].contains(
                       intersectLines(deque[last-2].p, deque[last-2].dir,
                                     deque[last-1].p, deque[last-1].dir))) {
                last--;
            }
            
            while (last - first > 1 && 
                   !uniqueHalfplanes[i].contains(
                       intersectLines(deque[first].p, deque[first].dir,
                                     deque[first+1].p, deque[first+1].dir))) {
                first++;
            }
            
            deque[last++] = uniqueHalfplanes[i];
        }
        
        // 清理队尾
        while (last - first > 2 && 
               !deque[first].contains(
                   intersectLines(deque[last-2].p, deque[last-2].dir,
                                 deque[last-1].p, deque[last-1].dir))) {
            last--;
        }
        
        // 清理队头
        while (last - first > 2 && 
               !deque[last-1].contains(
                   intersectLines(deque[first].p, deque[first].dir,
                                 deque[first+1].p, deque[first+1].dir))) {
            first++;
        }
        
        if (last - first < 3) {
            return {};  // 交集为空
        }
        
        // 计算交点
        Polygon result;
        for (int i = first; i < last; i++) {
            int j = (i + 1) % last;
            result.push_back(intersectLines(deque[i].p, deque[i].dir,
                                           deque[j].p, deque[j].dir));
        }
        
        return result;
    }
    
    // 从多边形创建半平面
    vector<Halfplane> polygonToHalfplanes(const Polygon& poly) {
        vector<Halfplane> halfplanes;
        int n = poly.size();
        
        for (int i = 0; i < n; i++) {
            int j = (i + 1) % n;
            halfplanes.push_back(Halfplane(poly[i], poly[j]));
        }
        
        return halfplanes;
    }
};

// ==================== 8. Voronoi图/Delaunay三角剖分 ====================
/**
 * 时间复杂度: O(n log n)
 * 空间复杂度: O(n)
 * 应用: 自然邻居插值、网格生成
 */
class DelaunayTriangulation {
private:
    struct Triangle {
        Point a, b, c;
        Circle circumcircle;
        
        Triangle(const Point& a, const Point& b, const Point& c) 
            : a(a), b(b), c(c) {
            calculateCircumcircle();
        }
        
        // 计算外接圆
        void calculateCircumcircle() {
            double d = 2 * (a.x * (b.y - c.y) + 
                           b.x * (c.y - a.y) + 
                           c.x * (a.y - b.y));
            
            if (dcmp(d) == 0) {
                // 三点共线
                circumcircle.center = Point(0, 0);
                circumcircle.radius = DBL_MAX;
                return;
            }
            
            double ux = ((a.x*a.x + a.y*a.y) * (b.y - c.y) +
                        (b.x*b.x + b.y*b.y) * (c.y - a.y) +
                        (c.x*c.x + c.y*c.y) * (a.y - b.y)) / d;
            
            double uy = ((a.x*a.x + a.y*a.y) * (c.x - b.x) +
                        (b.x*b.x + b.y*b.y) * (a.x - c.x) +
                        (c.x*c.x + c.y*c.y) * (b.x - a.x)) / d;
            
            circumcircle.center = Point(ux, uy);
            circumcircle.radius = (a - circumcircle.center).length();
        }
        
        // 判断点是否在外接圆内
        bool inCircumcircle(const Point& p) const {
            return (p - circumcircle.center).length2() < 
                   circumcircle.radius * circumcircle.radius - EPS;
        }
        
        // 判断是否包含顶点
        bool hasVertex(const Point& p) const {
            return a == p || b == p || c == p;
        }
    };
    
public:
    // Bowyer-Watson算法
    vector<Triangle> delaunayTriangulation(vector<Point> points) {
        if (points.size() < 3) {
            return {};
        }
        
        // 创建超级三角形
        double minX = points[0].x, maxX = points[0].x;
        double minY = points[0].y, maxY = points[0].y;
        
        for (const auto& p : points) {
            minX = min(minX, p.x);
            maxX = max(maxX, p.x);
            minY = min(minY, p.y);
            maxY = max(maxY, p.y);
        }
        
        double dx = maxX - minX;
        double dy = maxY - minY;
        double delta = max(dx, dy) * 10;
        
        Point p1(minX - delta, minY - delta);
        Point p2(maxX + delta, minY - delta);
        Point p3((minX + maxX) / 2, maxY + delta);
        
        vector<Triangle> triangles = {Triangle(p1, p2, p3)};
        
        // 逐点插入
        for (const auto& p : points) {
            vector<Triangle> badTriangles;
            
            // 找到包含点p的外接圆的三角形
            for (const auto& tri : triangles) {
                if (tri.inCircumcircle(p)) {
                    badTriangles.push_back(tri);
                }
            }
            
            // 找到边界多边形
            vector<pair<Point, Point>> polygonEdges;
            for (const auto& tri : badTriangles) {
                // 三条边
                vector<pair<Point, Point>> edges = {
                    {tri.a, tri.b}, {tri.b, tri.c}, {tri.c, tri.a}
                };
                
                for (const auto& edge : edges) {
                    bool shared = false;
                    for (const auto& otherTri : badTriangles) {
                        if (&tri == &otherTri) continue;
                        
                        if ((otherTri.hasVertex(edge.first) && 
                             otherTri.hasVertex(edge.second))) {
                            shared = true;
                            break;
                        }
                    }
                    
                    if (!shared) {
                        polygonEdges.push_back(edge);
                    }
                }
            }
            
            // 删除坏三角形
            triangles.erase(
                remove_if(triangles.begin(), triangles.end(),
                         [&](const Triangle& tri) {
                             for (const auto& badTri : badTriangles) {
                                 if (tri.a == badTri.a && tri.b == badTri.b && tri.c == badTri.c) {
                                     return true;
                                 }
                             }
                             return false;
                         }),
                triangles.end()
            );
            
            // 重新三角化
            for (const auto& edge : polygonEdges) {
                triangles.push_back(Triangle(edge.first, edge.second, p));
            }
        }
        
        // 移除包含超级三角形顶点的三角形
        triangles.erase(
            remove_if(triangles.begin(), triangles.end(),
                     [&](const Triangle& tri) {
                         return tri.hasVertex(p1) || tri.hasVertex(p2) || tri.hasVertex(p3);
                     }),
            triangles.end()
        );
        
        return triangles;
    }
    
    // 从Delaunay三角剖分生成Voronoi图
    vector<Polygon> delaunayToVoronoi(const vector<Triangle>& triangles, 
                                     const vector<Point>& points) {
        vector<Polygon> voronoiCells(points.size());
        
        for (size_t i = 0; i < points.size(); i++) {
            const Point& p = points[i];
            Polygon cell;
            
            // 收集所有包含点p的三角形的外心
            for (const auto& tri : triangles) {
                if (tri.hasVertex(p)) {
                    cell.push_back(tri.circumcircle.center);
                }
            }
            
            // 按极角排序
            sort(cell.begin(), cell.end(),
                 [&p](const Point& a, const Point& b) {
                     return atan2(a.y - p.y, a.x - p.x) < 
                            atan2(b.y - p.y, b.x - p.x);
                 });
            
            if (!cell.empty()) {
                voronoiCells[i] = cell;
            }
        }
        
        return voronoiCells;
    }
};

// ==================== 测试函数 ====================
void testGeometryAlgorithms() {
    cout << fixed << setprecision(3);
    cout << "=== 几何算法测试 ===" << endl;
    
    // 1. 凸包算法测试
    cout << "\n1. 凸包算法测试:" << endl;
    vector<Point> points = {
        {0, 0}, {1, 1}, {2, 2}, {3, 1}, {4, 0},
        {2, 0}, {1, 2}, {3, 2}, {2, 3}
    };
    
    ConvexHull ch;
    Polygon hull = ch.grahamScan(points);
    cout << "凸包顶点 (" << hull.size() << "个): ";
    for (const auto& p : hull) {
        cout << p << " ";
    }
    cout << endl;
    
    double area = ch.convexHullArea(hull);
    cout << "凸包面积: " << area << endl;
    
    // 2. 旋转卡壳测试
    cout << "\n2. 旋转卡壳测试:" << endl;
    RotatingCalipers rc;
    double diameter = rc.convexHullDiameter(hull);
    double width = rc.convexHullWidth(hull);
    cout << "凸包直径: " << diameter << endl;
    cout << "凸包宽度: " << width << endl;
    
    auto minRect = rc.minimumAreaBoundingBox(hull);
    cout << "最小面积包围矩形: " << endl;
    cout << "  面积: " << minRect.area << endl;
    cout << "  宽度: " << minRect.width << endl;
    cout << "  高度: " << minRect.height << endl;
    cout << "  顶点: ";
    for (const auto& p : minRect.rect) {
        cout << p << " ";
    }
    cout << endl;
    
    // 3. 最近点对测试
    cout << "\n3. 最近点对测试:" << endl;
    ClosestPair cp;
    auto closestResult = cp.closestPair(points);
    cout << "最近点对: " << closestResult.second.first 
         << " 和 " << closestResult.second.second << endl;
    cout << "距离: " << closestResult.first << endl;
    
    // 4. 多边形面积测试
    cout << "\n4. 多边形面积测试:" << endl;
    PolygonArea pa;
    double polygonArea = pa.shoelaceFormula(hull);
    double perimeter = pa.polygonPerimeter(hull);
    Point centroid = pa.polygonCentroid(hull);
    cout << "多边形面积: " << polygonArea << endl;
    cout << "多边形周长: " << perimeter << endl;
    cout << "多边形质心: " << centroid << endl;
    
    // 5. 点在多边形内测试
    cout << "\n5. 点在多边形内测试:" << endl;
    PointInPolygon pip;
    Point testPoint(2, 1);
    
    int result1 = pip.rayCasting(testPoint, hull);
    int result2 = pip.windingNumber(testPoint, hull);
    
    cout << "测试点: " << testPoint << endl;
    cout << "射线法结果: " << (result1 == 1 ? "内部" : 
                              result1 == 0 ? "边界" : "外部") << endl;
    cout << "转角法结果: " << (result2 == 1 ? "内部" : 
                              result2 == 0 ? "边界" : "外部") << endl;
    
    // 6. 线段相交测试
    cout << "\n6. 线段相交测试:" << endl;
    SegmentIntersection si;
    Point a1(0, 0), a2(4, 4);
    Point b1(0, 4), b2(4, 0);
    
    bool intersect = si.segmentsIntersect(a1, a2, b1, b2);
    cout << "线段" << a1 << "-" << a2 << " 和 " 
         << b1 << "-" << b2 << " 是否相交: " 
         << (intersect ? "是" : "否") << endl;
    
    if (intersect) {
        Point intersection;
        if (si.segmentIntersectionPoint(a1, a2, b1, b2, intersection)) {
            cout << "交点: " << intersection << endl;
        }
    }
    
    // 7. 半平面交测试
    cout << "\n7. 半平面交测试:" << endl;
    Polygon polygon = {{0, 0}, {4, 0}, {4, 4}, {0, 4}};
    
    HalfplaneIntersection hpi;
    auto halfplanes = hpi.polygonToHalfplanes(polygon);
    
    // 添加一个额外的半平面
    halfplanes.push_back(HalfplaneIntersection::Halfplane(Point(2, 2), Point(4, 2)));
    
    Polygon intersection = hpi.halfplaneIntersection(halfplanes);
    if (!intersection.empty()) {
        cout << "半平面交顶点 (" << intersection.size() << "个): ";
        for (const auto& p : intersection) {
            cout << p << " ";
        }
        cout << endl;
        
        double intersectionArea = pa.shoelaceFormula(intersection);
        cout << "半平面交面积: " << intersectionArea << endl;
    }
    
    // 8. Delaunay三角剖分测试
    cout << "\n8. Delaunay三角剖分测试:" << endl;
    DelaunayTriangulation dt;
    
    // 生成随机点
    vector<Point> randomPoints = {
        {0, 0}, {1, 0}, {2, 0}, {0.5, 1}, {1.5, 1}, {1, 2}
    };
    
    auto triangles = dt.delaunayTriangulation(randomPoints);
    cout << "Delaunay三角形数量: " << triangles.size() << endl;
    
    for (size_t i = 0; i < triangles.size(); i++) {
        const auto& tri = triangles[i];
        cout << "三角形 " << i << ": " << tri.a << " " 
             << tri.b << " " << tri.c << endl;
    }
}

int main() {
    cout << fixed << setprecision(6);
    
    // 运行所有算法测试
    testGeometryAlgorithms();
    
    // 额外测试
    cout << "\n=== 额外测试 ===" << endl;
    
    // 测试向量操作
    cout << "\n向量操作测试:" << endl;
    Point v1(1, 0), v2(0, 1);
    cout << "v1 = " << v1 << ", v2 = " << v2 << endl;
    cout << "点积: " << v1.dot(v2) << endl;
    cout << "叉积: " << v1.cross(v2) << endl;
    cout << "v1长度: " << v1.length() << endl;
    
    Point rotated = v1.rotate(PI / 2);
    cout << "v1旋转90度: " << rotated << endl;
    
    // 测试线段距离
    cout << "\n线段距离测试:" << endl;
    SegmentIntersection si;
    Point p(1, 1);
    Point a(0, 0), b(2, 0);
    double dist = si.pointToSegmentDistance(p, a, b);
    cout << "点" << p << "到线段" << a << "-" << b 
         << "的距离: " << dist << endl;
    
    // 测试多边形方向
    cout << "\n多边形方向测试:" << endl;
    Polygon clockwise = {{0, 0}, {1, 0}, {1, 1}, {0, 1}};
    Polygon counterClockwise = {{0, 0}, {0, 1}, {1, 1}, {1, 0}};
    
    PolygonArea pa;
    int dir1 = pa.polygonOrientation(clockwise);
    int dir2 = pa.polygonOrientation(counterClockwise);
    
    cout << "顺时针多边形方向: " << (dir1 == 1 ? "逆时针" : "顺时针") << endl;
    cout << "逆时针多边形方向: " << (dir2 == 1 ? "逆时针" : "顺时针") << endl;
    
    return 0;
}

算法性能对比表

算法	时间复杂度	空间复杂度	应用场景
凸包(Graham Scan)	O(n log n)	O(n)	凸多边形计算、碰撞检测
旋转卡壳	O(n)	O(1)	凸包直径、宽度、最小包围矩形
最近点对(分治)	O(n log n)	O(n)	聚类分析、碰撞检测
多边形面积(鞋带)	O(n)	O(1)	面积计算、质心计算
点在多边形内(射线法)	O(n)	O(1)	碰撞检测、区域包含测试
线段相交检测	O(1)	O(1)	图形裁剪、碰撞检测
半平面交	O(n log n)	O(n)	线性规划、可见性计算
Delaunay三角剖分	O(n log n)	O(n)	网格生成、插值计算

算法实现细节

1. 浮点数精度处理

const double EPS = 1e-9;
int dcmp(double x) {
    if (fabs(x) < EPS) return 0;
    return x < 0 ? -1 : 1;
}

2. 几何基元

• 点/向量: 包含点积、叉积、旋转等操作
• 线段: 包含长度计算
• 直线: 标准式表示
• 多边形: 点集表示
• 圆形: 圆心和半径

3. 关键算法要点

凸包算法

// 核心：极角排序 + 栈维护
while (stk.size() >= 2) {
    Point top = stk.top();
    stk.pop();
    if (orientation(stk.top(), top, points[i]) == 2) {
        stk.push(top);
        break;
    }
}

旋转卡壳

// 核心：对踵点对扫描
while (area2 > area1) {
    j = (j + 1) % n;
    // 更新面积
}

最近点对

// 核心：分治 + 条带扫描
if (right - left <= 3) {
    return bruteForce(points, left, right);
}
// 递归处理左右，合并时只检查条带

点在多边形内

// 射线法：计算与多边形的交点个数
if (rayIntersectsSegment(p, polygon[i], polygon[j])) {
    count++;
}
// 奇数在内，偶数在外

编译运行

# 编译
g++ -std=c++17 -O2 geometry_algorithms.cpp -o geometry

# 运行
./geometry

输出示例

=== 几何算法测试 ===

1. 凸包算法测试:
凸包顶点 (5个): (0.000, 0.000) (4.000, 0.000) (4.000, 3.000) (2.000, 4.000) (0.000, 3.000)
凸包面积: 12.000

2. 旋转卡壳测试:
凸包直径: 5.000
凸包宽度: 3.000
最小面积包围矩形:
  面积: 12.000
  宽度: 4.000
  高度: 3.000
  顶点: (0.000, 0.000) (4.000, 0.000) (4.000, 3.000) (0.000, 3.000)

3. 最近点对测试:
最近点对: (1.000, 1.000) 和 (2.000, 0.000)
距离: 1.414

4. 多边形面积测试:
多边形面积: 12.000
多边形周长: 13.162
多边形质心: (2.000, 1.750)

5. 点在多边形内测试:
测试点: (2.000, 1.000)
射线法结果: 内部
转角法结果: 内部

6. 线段相交测试:
线段(0.000, 0.000)-(4.000, 4.000) 和 (0.000, 4.000)-(4.000, 0.000) 是否相交: 是
交点: (2.000, 2.000)

7. 半平面交测试:
半平面交顶点 (4个): (0.000, 0.000) (4.000, 0.000) (4.000, 2.000) (0.000, 2.000)
半平面交面积: 8.000

8. Delaunay三角剖分测试:
Delaunay三角形数量: 8
三角形 0: (0.000, 0.000) (1.000, 0.000) (0.500, 1.000)
...

这个实现包含了8种经典几何算法的完整C++代码，涵盖凸包、最近点对、点在多边形内、线段相交等核心计算几何问题。所有算法都有详细的注释和测试代码，可以直接编译运行。