计算几何¶

点和线¶

精度控制¶

const double eps = 1e-10;
double dcmp(double x) { //判断 正/负/0
    if(fabs(x) < eps) return 0;
    else return x > 0 ? 1 : -1;
}

点和向量存储¶

struct Vector{
    double x, y;
    Point(double x = 0, y = 0):x(x), y(y){}
    Vector rotate_90() {return Vector(y, -x);} //顺时针90度
    Vector rotate_90_c() {return Vector(-y, x);} //逆时针90度
};
typedef Vector Point;

向量的四则运算¶

//这种返回方式似乎有时候会出问题
//向量+向量=向量  点+向量=点
Vector operator+(Vector a, Vector b){ 
    return Vector(a.x+b.x, a.y+b.y); 
}
//点-点=向量
Vector operator-(Vector a, Vector b){ 
    return Vector(a.x-b.x, a.y-b.y); 
}
//向量*数=向量
Vector operator*(Vector a, double p){ 
    return Vector(a.x*p, a.y*p); 
}
//向量/数=向量
Vector operator/(Vector a, double p){ 
    return Vector(a.x/p. a.y/p); 
}

点积叉积及应用¶

\(\vec{a} \cdot \vec{b} = |a||b|cos\theta\)

\(|\vec{a} \times \vec{b}| = |a||b|sin\theta\)，正负由两向量的相对关系决定

//点积
double Dot(Vector a, Vector b){ 
    return a.x*b.x+a.y*b.y; 
}
//叉积
double Cross(Vector a, Vector b){ 
    return a.x*b.y-a.y*b.x; 
}

//应用
//求模长
double Length(Vector a){ 
    return sqrt(Dot(a, a)); 
}
//求角度
double Angle(Vector a, Vector b){ 
    return acos(Dot(a, b)/Length(a)/Length(b));
}

判断点是否在线段上¶

bool OnSegment(Point p1, Point p2, Point q){  //p1-p2是线段
    return Cross(p1-q, p2-q) == 0 && Dot(p1-q, p2-q) <= 0;
}

求两直线的交点¶

由 \(q_2,q_1,p_2,p_1\)确定的两条直线，假设交点为 \(p_1+t(p_2-p_1)\)

有 \((q_2-q_1) \times (p_1 + t(p_2-p_1)-q_1)=0\)

得 \(t = \frac{(q_2-q_1)\times(q_1-p_1)}{(q_2-q_1)\times(p_2-p_1)}\)

/*Point GetLineIntersection(Point p1, Point p2, Point q1, Point q2){
    return p1+(p2-p1)*Cross(q2-q1, q1-p1)/Cross(q2-q1, p2-p1)
}*/
Point GetLineIntersection(Point p, Vector v1, Point q, Vector v2){
    return p+v1*Cross(v2, q-p)/Cross(v2, v1);
} //平行特判

线段相交判定¶

跨立实验：判断q1,q2在线段p2-p1两侧且 p1,p2在线段q2-q1两侧

bool SegmentProperIntersection(Point p1, Point p2, Point q1, Point q2){
    double a1 = Cross(p2-p1, q2-p1), a2 = Cross(p2-p1, q1-p1);
    double a3 = Cross(q2-q1, p1-q1), a4 = Cross(q2-q1, p2-q1);
    return dcmp(a1)*dcmp(a2)<0 && dcmp(a3)*dcmp(a4)<0;
} 
/*特判端点在另一条线段的情况
if(OnSegment(p[i],q[i],p[j]) 
|| OnSegment(p[i],q[i],q[j]) 
|| OnSegment(p[j],q[j],p[i])
|| OnSegment(p[j],q[j],q[i])
)  G[i][j] = G[j][i] = true;
*/

线段和直线相交判定¶

bool fun(Point a, Point b, Point c, Point d) { //a b 确定直线     c d 确定线段
    return Cross(a - b, a - c) * Cross(a - b, a - d) <= eps
}

向量旋转与单位化¶

逆时针：

\[ \left[ \matrix{ x^{'} \\ y^{'} } \right] = \left[ \matrix{ x \\ y } \right] \left[ \matrix{ cos\theta & sin\theta\\ -sin\theta & cos\theta } \right] \]

顺时针：

\[ \left[ \matrix{ x^{'} \\ y^{'} } \right] = \left[ \matrix{ x \\ y } \right] \left[ \matrix{ cos\theta & -sin\theta\\ sin\theta & cos\theta } \right] \]

//rad 是弧度, 逆时针
Vector Rotate(Vector a, double rad) {
    return Vector(a.x*cos(rad)-a.y*sin(rad), a.x*sin(rad)+a.y*soc(rad));
}
//单位化
Vector Normal(Vector a) {
    double L = length(a);
    return Vector(a.x/L, a.y/L);
}

点到直线的距离¶

double DistaceToLine(Point p, Point a, Point b) {
    Vector v1 = b-a, v2 = p-a;
    return fabs(Cross(v1, v2))/Length(v1); //平行四边形面积 除/底=高
}

圆¶

圆的存储¶

struct Circle{
    Point c;
    double r;
    Circle(Point a, double b = 0):c(a), r(b){}
}

三角形外接圆¶

//三点不共线
Circle circumcircle(Point p1, Point p2, Point p3) {
    Point p = GetLineIntersection((p1+p2)/2, (p2-p1).rotate_90(), (p1+p3)/2, (p3-p1).rotate_90());
    return Circle(p, Length(p1-p));
}

最小圆覆盖¶

随机增量法 \(O(n)\)

//洛谷 P1742
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;

const int maxn = 1e5+5;
typedef struct Vector{
    double x, y;
    Vector(double a = 0, double b = 0):x(a),y(b){}
    Vector rotate_90(){return Vector(y, -x); }
}Point;

struct Circle{
    Point c;
    double r;
    Circle(Point a, double b):c(a),r(b){}
};

Point operator-(const Point& a, const Point& b){ return Vector(a.x-b.x, a.y-b.y); }
Point operator+(const Point& a, const Point& b){ return Vector(a.x+b.x, a.y+b.y); }
Point operator/(const Point& a, const double& b){ return Vector(a.x/b, a.y/b); }
Point operator*(const Point& a, const double& b){ return Vector(a.x*b, a.y*b); }
double Dot(Vector a, Vector b){ return a.x*b.x+a.y*b.y; }
double Cross(Vector a, Vector b){ return a.x*b.y-a.y*b.x; }
double Length(Vector a){ return sqrt(Dot(a, a)); }

Point GetLineIntersection(Point p, Vector v1, Point q, Vector v2){
    return p+v1*Cross(v2, q-p)/Cross(v2, v1);
}

Circle circumcircle(Point p1, Point p2, Point p3){
    Point p = GetLineIntersection((p1+p2)/2, (p2-p1).rotate_90(), (p1+p3)/2, (p3-p1).rotate_90());
    return Circle(p, Length(p1-p));
}

Circle min_circle(Point *p, int n)  { //最小圆覆盖
    random_shuffle(p+1, p+1+n);
    Point o;
    double r = 0;
    for(int i = 1; i <= n; i++)
        if(Length(p[i]-o) > r)
        {
            o = p[i], r = 0;
            for(int j = 1; j < i; j++)
                if(Length(p[j]-o) > r)
                {
                    o = (p[i]+p[j])/2, r = Length(p[i]-p[j])/2;
                    for(int k = 1; k < j; k++)
                        if(Length(o-p[k]) > r)
                        {
                            Circle tmp = circumcircle(p[i], p[j], p[k]);
                            o = tmp.c, r = tmp.r;
                        }
                }
        }
    return Circle(o, r);
} 

int main(){
    int n;
    Point p[maxn];
    cin >> n;
    for(int i = 1; i <= n; i++) cin >> p[i].x >> p[i].y;
    Circle res =  min_circle(p, n);
    printf("%.10lf\n%.10lf %.10lf\n", res.r, res.c.x, res.c.y);
    return 0;
}

凸包¶

Graham扫描法算凸包¶

按x，y坐标升序排列，确定凸包的2的顶点，然后再确定上下两个凸壳

排序 \(O(nlogn)\) 其他 \(O(n)\)

bool cmp(Point &p, Point &q){
    if(p.x != p.x) return p.x < q.x;
    else return p.y < q.y;
}

vector<Point> ConvexHull(Point *ps, int n){
    sort(ps, ps+n, cmp);
    int k = 0;  //凸包顶点数
    vector<Point> qs(n*2);
    //下凸壳
    for(int i = 0; i < n; i++){  //while要注意
        while(k > 1 && Cross(qs[k-1]-qs[k-2], ps[i]-qs[k-1] <= 0)) k--;
        qs[k++] = ps[i];
    }
    //上凸壳
    for(int i = n-2, t = k; i >= 0; i--){
        while(k > t && Cross(qs[k-1]-qs[k-2], ps[i]-qs[k-1] <= 0)) k--;
        qs[k++] = ps[i];
    }
    qs.resize(k-1);
    return qs;
}