计算几何

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

好久没写博客了= =

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

图形存储

• 点：我们可以直接存储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}+tvec{v_1},vec{p_2}+tvec{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_1cdot x_2+y_1cdot y_2=|vec{a}|cdot|vec{b}|cdot cosleft langlevec{a},vec{b} ight angle)

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_1cdot y_2-x_2cdot y_1=|vec{a}|cdot|vec{b}|cdot sinleft langlevec{a},vec{b} ight angle)

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

• 向量旋转(alpha)弧度：(x'=xcosalpha-ysinalpha,y'=xsinalpha+ycosalpha)

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_1vec{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，有空以后补上