计算几何

超级全的计算几何全家桶~

好久没写博客了= =

正好刚听学长讲了计算几何，就收录下一些相关操作吧

图形存储

• 点：我们可以直接存储x,y坐标或x,y,z坐标

• 向量：起点在原点的向量用x,y坐标或x,y,z坐标表示，任意一个向量可以用两个向量相减的形式表示

struct node //或Vector
{
    double x,y;
};

• 直线：用两个点或一个点和一个向量的形式

struct line  //直线
{
    node p,v;
};

• 多边形：按顺时针或逆时针把点存起来，相邻点有连线

struct polygon  //多边形
{
    int n;  //边数
    node p[N + 5];  //点
};

误差避免

• 计算机存储浮点数会有些许误差，所以一般我们设一个很小的数\(eps=10^{-6}\)来避免误差

const double eps = 1e-6;
int com(double a,double b)  //比较a和b的大小
{
    double x = a - b;
    if (fabs(x) < eps)
        return 0;
    else
        if (x > 0)
            return 1;
        else
            return -1;
}

基础计算

后面的\(\vec{a},\vec{b}\)均为\(\vec{a}=(x_1,y_1),\vec{b}=(x_2,y_2)\)，直线\(l_1,l_2\)均为\(\vec{p_1}+t\vec{v_1},\vec{p_2}+t\vec{v_2}\)

• 两点之间距离：\(dist=\sqrt{(x_1-x_2)^2+(y_1-y_2)^2}\)

double dist(node a,node b)  //两点之间的距离
{
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}

• 向量加、减、数乘、数除：x,y对应操作就行了

node operator +(const node &a,const node &b)
{
    return (node){a.x + b.x,a.y + b.y};
}
node operator -(const node &a,const node &b)
{
    return (node){a.x - b.x,a.y - b.y};
}
node operator *(const node &a,const double &b)
{
    return (node){a.x * b,a.y * b};
}
node operator /(const node &a,const double &b)
{
    return (node){a.x / b,a.y / b};
}

• 点乘：\(dot(\vec{a},\vec{b})=x_1\cdot x_2+y_1\cdot y_2=|\vec{a}|\cdot|\vec{b}|\cdot \cos\left \langle\vec{a},\vec{b}\right \rangle\)

double dot(node a,node b)   //点乘
{
    return a.x * b.x + a.y * b.y;
}

• 模长：\(|\vec{a}|=\sqrt{|dot(\vec{a},\vec{a})|}\)

double molen(node a)  //模长
{
    return sqrt(fabs(dot(a,a)));
}

• 叉积：\(cross(\vec{a},\vec{b})=x_1\cdot y_2-x_2\cdot y_1=|\vec{a}|\cdot|\vec{b}|\cdot \sin\left \langle\vec{a},\vec{b}\right \rangle\)

double cross(node a,node b)  //叉积
{
    return a.x * b.y - a.y * b.x;
}

• 向量旋转\(\alpha\)弧度：\(x'=x\cos\alpha-y\sin\alpha,y'=x\sin\alpha+y\cos\alpha\)

node rot(node x,double a)  //旋转
{
    return (node){x.x * cos(a) - x.y * sin(a),x.x * sin(a) + x.y * cos(a)};
}

• 直线交点：由方程算出交点在任何一条直线上的\(t\)，再算交点，设\(\vec{u}=\vec{p_1}-\vec{p_2},t_1=\frac{cross(\vec{v_1},\vec{u})}{cross(\vec{v_1},\vec{v_2})}\)，交点即为\(\vec{p_1}+t_1\vec{v_1}\)

node linepoint(line a,line b)  //直线交点
{
    node u = a.p - b.p;
    double t = cross(a.v,u) / cross(a.v,b.v);
    return a.p + a.v * t;
}

• 点到直线的距离：在直线任取两个点用等积法求高

double pointlinedist(node a,line l)  //点到直线的距离
{
    node x = l.p + l.v * 2.33,y = l.p + l.v * 6.66;
    return fabs(cross(a - x,y - x)) / molen(x - y);
}

• 点到线段的距离：要分两种情况讨论

1. 过点的垂线不过线段：距离为点和线段端点的连线长度
2. 过点的垂线过线段：距离就是垂线长度

至于怎么判断：利用点乘的性质(夹角小于\(90^\circ\)为正否则为负)，设线段两端点为\(A,B\)，线段外的点为\(C\)，若\(dot(\vec{AC},\vec{AB})<0\)或\(dot(\vec{AB},\vec{BC})>0\)（反向），则为情况1，否则为情况2

double pointline_Seg(node c,node a,node b)  //点到线段的距离
{
    node l1 = b - a,l2 = c - a,l3 = c - b;
    if (com(dot(l1,l2),0) < 0)
        return molen(l2);
    else
        if (com(dot(l1,l3),0) > 0)
            return molen(l3);
        else
            return fabs(cross(l1,l2)) / molen(l1);
}

• 多边形面积计算：将一个有\(n\)个边的划分成\(n-2\)个三角形，然后算面积（也就是叉积）的和，注意叉积不能取绝对值，要最后取绝对值，因为对于凹多边形会有多算的面积

double polygon_S(polygon a)
{
    double s = 0;
    for (int i = 1;i < a.n - 1;i++)
        s += cross(a.p[i + 1] - a.p[1],a.p[i + 2] - a.p[1]) / 2;
    return fabs(s);
}

先写这么多吧QAQ，有空以后补上