【计算几何 01】叉积

这几天闲来无事去学习了一下计算几何，发现其实不（sang）是（xin）太（bing）难（kuang）😛
今天就重点介绍一下简单的叉积及其简单的运用（毕竟作为蒟蒻，难的搞不来啊）

什么是计算几何？

“对几何外形信息的计算机表示、分析和综合”——福雷斯特

其实所谓计算几何，就是用计算机去解决不同维度上的几何问题，而我们要做的，就是设计一套算法，让计算机去分析和解决一系列的几何问题。

一些基本知识

坐标的存储

一般来说，计算几何都涉及到了小数坐标，这时开一个含两个double类型的结构体就比较稳健，我一般的定义如下：

//定义一个点
struct node { double x, y;} a[11000];

线段的存储

聪明绝顶的你肯定已经想到了，线段就是包含两个点的一个结构体，所以我一般的定义如下

//定义一条线
struct line { node p1, p2;} list[11000];

叉积

终于到今天的主角—叉积了！😝
那么什么是叉积呢？叉积其实是两个矢量模的乘积再乘夹角正弦，经过推导可以发现两个向量\(A(x1,y1),B(x2,y2)\)的叉积为：

\[x1*y2-x2*y1 \]

先丢出代码：

//计算叉积，p1，p2，p0都为点
double multi(node p1, node p2, node p0) {
    double x1, y1, x2, y2;
    x1 = p1.x - p0.x;
    y1 = p1.y - p0.y;
    x2 = p2.x - p0.x;
    y2 = p2.y - p0.y;
    return x1 * y2 - x2 * y1;
}

这就是传说中的叉积计算了！是不是很短呢，那就让我们来理解一下这个multi函数吧。

由于叉积是向量间的运算，大家都知道向量可以用末坐标减去首坐标得到，那么其实multi函数计算的就是向量A(x1,y1)与向量B(x2,y2)的叉积。

对于multi函数的意义，首先multi的正负是有特殊含义的：以p0为参考点，如果multi大于0，则p2在p1的逆时针方向，反正，如果multi小于0，则p2在p1的顺时针方向，特殊的，当multi等于0，p1、p2、p0三点共线。

其次，multi的值也是有含义的！通过后面的学习你就会发现，multi的绝对值就是以p0、p1、p2三点构成的三角形的面积的两倍！利用这个规律，可以完成很多与计算几何有关的题目。

一些简单的运用

【caioj 1212 判断线段相交（跨立实验）】

题目描述
【题意】
有n条线段（编号为1_n），按1n的顺序放在二维坐标系上（就是先放1号，再放2号……），
要求输出最上面的那些线段的编号（就是没有其他线段压在它上面的那些线段）
【输入格式】
第一行第一个数n（ 1 <= n <= 10000）表示这组数据有n条线段。
下来n行，每行两个坐标，表示第i条线段的两个端点。
【输出格式】
一行。输出最上面的线段的编号（从小到大）。相邻两个编号用空格隔开，最后一个编号没有空格。
【样例1输入】
5
1 1 4 2
2 3 3 1
1 -2.0 8 4
1 4 8 2
3 3 6 -2.0
【样例1输出】
2 4 5
【样例2输入】
3
0 0 1 1
1 0 2 1
2 0 3 1
【样例2输出】
1 2 3

这道题很明显要让我们判断任意两条线段是否相交，这就要用到刚刚讲到的叉积的第一种用法了!

很明显，相交有两种情况，相交且穿过和相交但不穿过（如下图）

我们首先先讨论第一种情况，如何判断两条直线相交且穿过？

如上图，先假设两条直线为L1和L2，每条直线的两个端点分别为p1和p2，那么如果以
L1.p1为一个参考点，如果L1.p2在L2.p1与L2.p2中间；再以L2.p1为参考点，如果L2.p2在L1.p1与L1.p2中间，那么就可以判断L1与L2相交且穿过（如下图）

这种情况非常好理解吧，下一种情况就更加简单了，仅需讨论哪三个点在一条“直线”上（毕竟不穿过，所以必有三点共线），然后判断一下
不穿过的那个点是否在一条“线段”上（如下图）

具体代码如下

double check(line l1, line l2) { //判断是否相交
    if (multi(l2.p1, l1.p2, l1.p1)*multi(l1.p2, l2.p2, l1.p1) > 0) //判断以L1.p1为基准，L1.p2是否在L2.p1与L2.p2中间
        if (multi(l1.p1, l2.p2, l2.p1)*multi(l2.p2, l1.p2, l2.p1) > 0) //判断以L2.p1为基准,L2.p2是否在L1.p1与L2.p2中间
            return true;//大部分情况下，能判断两直线相交
    if (multi(l1.p1, l1.p2, l2.p1) == 0)if (mymin(l1.p1.x, l1.p2.x) <= l2.p1.x)if (mymax(l1.p1.x, l1.p2.x) >= l2.p1.x)if (mymin(l1.p1.y, l1.p2.y) <= l2.p1.y)if (mymax(l1.p1.y, l1.p2.y) >= l2.p1.y)return true;
    if (multi(l1.p1, l1.p2, l2.p2) == 0)if (mymin(l1.p1.x, l1.p2.x) <= l2.p2.x)if (mymax(l1.p1.x, l1.p2.x) >= l2.p2.x)if (mymin(l1.p1.y, l1.p2.y) <= l2.p2.y)if (mymax(l1.p1.y, l1.p2.y) >= l2.p2.y)return true;
    if (multi(l2.p1, l2.p2, l1.p1) == 0)if (mymin(l2.p1.x, l2.p2.x) <= l1.p1.x)if (mymax(l2.p1.x, l2.p2.x) >= l1.p1.x)if (mymin(l2.p1.y, l2.p2.y) <= l1.p1.y)if (mymax(l2.p1.y, l2.p2.y) >= l1.p1.y)return true;
    if (multi(l2.p1, l2.p2, l1.p2) == 0)if (mymin(l2.p1.x, l2.p2.x) <= l1.p2.x)if (mymax(l2.p1.x, l2.p2.x) >= l1.p2.x)if (mymin(l2.p1.y, l2.p2.y) <= l1.p2.y)if (mymax(l2.p1.y, l2.p2.y) >= l1.p2.y)return true;
    //判断两条直线相交且一条直线端点在另一条直线上
    return false;
}

在之前的版本里，我在判断两条直线相交且一条直线端点在另一条直线上的情况中忽略了两条直线互相一条与x轴垂直，一条与y轴垂直，且某一条延伸出去恰好交于另一条直线交点的情况，于是在一次比赛中WA了十多发QwQ，下面放上错误代码以作为警示：

//判断是否相交 
double check(line l1,line l2){
    if(multi(l2.p1,l1.p2,l1.p1)*multi(l1.p2,l2.p2,l1.p1)>0)//判断以L1.p1为基准，L1.p2是否在L2.p1与L2.p2中间
    if(multi(l1.p1,l2.p2,l2.p1)*multi(l2.p2,l1.p2,l2.p1)>0)//判断以L2.p1为基准,L2.p2是否在L1.p1与L2.p2中间
    return true;//大部分情况下，能判断两直线相交   if(multi(l2.p1,l1.p1,l1.p2)==0)if(mymin(l1.p1.x,l1.p2.x<=l2.p1.x)if(mymax(l1.p1.x,l1.p2.x)>=l2.p1.x)return true;
    if(multi(l2.p2,l1.p1,l1.p2)==0)if(mymin(l1.p1.x,l1.p2.x<=l2.p2.x)if(mymax(l1.p1.x,l1.p2.x)>=l2.p2.x)return true;
    if(multi(l1.p1,l2.p1,l2.p2)==0)if(mymin(l2.p1.x,l2.p2.x<=l1.p1.x)if(mymax(l2.p1.x,l2.p2.x)>=l1.p1.x)return true;
    if(multi(l1.p2,l2.p1,l2.p2)==0)if(mymin(l2.p1.x,l2.p2.x<=l1.p2.x)if(mymax(l2.p1.x,l2.p2.x)>=l1.p2.x)return true;
    //判断两条直线相交且一条直线端点在另一条直线上 
    return false;
}

主程序就很简单了，一条一条边读入，读入第i条边时判断他与前i-1条边是否相交，若相交，这之前的边j被盖住，判断数组bo【j】=1。最后统计bo【i】==0的边就好啦
附上完整代码：

#include<cmath>
#include<cstdlib>
#include<cstdio>
#include<cstring>
using namespace std;
//定义一个点
struct node { double x, y;};
//定义一条线
struct line {  node p1, p2;} list[11000];
int a[11000] = {0};
bool bo[11000] = {0};
int i, j, k, m, n, o, p, js, jl, jk, len;

//================================
//求max
double mymax(double x, double y) { if (x > y)return x; else return y;}
//求min
double mymin(double x, double y) { if (x < y)return x; else return y;}
//=================================

double multi(node p1, node p2, node p0) { //计算叉积：而且值为三点构成三角形的面积的两倍
    double x1, y1, x2, y2;
    x1 = p1.x - p0.x;
    y1 = p1.y - p0.y;
    x2 = p2.x - p0.x;
    y2 = p2.y - p0.y;
    return x1 * y2 - x2 * y1;
}

double check(line l1, line l2) { //判断是否相交
    if (multi(l2.p1, l1.p2, l1.p1)*multi(l1.p2, l2.p2, l1.p1) > 0) //判断以L1.p1为基准，L1.p2是否在L2.p1与L2.p2中间
        if (multi(l1.p1, l2.p2, l2.p1)*multi(l2.p2, l1.p2, l2.p1) > 0) //判断以L2.p1为基准,L2.p2是否在L1.p1与L2.p2中间
            return true;//大部分情况下，能判断两直线相交
    if (multi(l1.p1, l1.p2, l2.p1) == 0)if (mymin(l1.p1.x, l1.p2.x) <= l2.p1.x)if (mymax(l1.p1.x, l1.p2.x) >= l2.p1.x)if (mymin(l1.p1.y, l1.p2.y) <= l2.p1.y)if (mymax(l1.p1.y, l1.p2.y) >= l2.p1.y)return true;
    if (multi(l1.p1, l1.p2, l2.p2) == 0)if (mymin(l1.p1.x, l1.p2.x) <= l2.p2.x)if (mymax(l1.p1.x, l1.p2.x) >= l2.p2.x)if (mymin(l1.p1.y, l1.p2.y) <= l2.p2.y)if (mymax(l1.p1.y, l1.p2.y) >= l2.p2.y)return true;
    if (multi(l2.p1, l2.p2, l1.p1) == 0)if (mymin(l2.p1.x, l2.p2.x) <= l1.p1.x)if (mymax(l2.p1.x, l2.p2.x) >= l1.p1.x)if (mymin(l2.p1.y, l2.p2.y) <= l1.p1.y)if (mymax(l2.p1.y, l2.p2.y) >= l1.p1.y)return true;
    if (multi(l2.p1, l2.p2, l1.p2) == 0)if (mymin(l2.p1.x, l2.p2.x) <= l1.p2.x)if (mymax(l2.p1.x, l2.p2.x) >= l1.p2.x)if (mymin(l2.p1.y, l2.p2.y) <= l1.p2.y)if (mymax(l2.p1.y, l2.p2.y) >= l1.p2.y)return true;
    //判断两条直线相交且一条直线端点在另一条直线上
    return false;
}
//==================================
int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; i++)scanf("%lf%lf%lf%lf", &list[i].p1.x, &list[i].p1.y, &list[i].p2.x, &list[i].p2.y);
    for (int i = 1; i < n; i++)
        for (int j = i + 1; j <= n; j++)
            if (check(list[i], list[j])) { //判断相交
                bo[i] = 1;
                break;
            }
    len = 0;
    for (int i = 1; i <= n; i++)if (bo[i] == 0)a[++len] = i;
    for (int i = 1; i < len; i++)printf("%d ", a[i]);
    printf("%d\n", a[len]);
    return 0;
}

【caioj 1213 面积】

题目描述
【题意】
在一个平面坐标系上随意画一条有n个点的封闭折线（按画线的顺序给出点的坐标），保证封闭折线的任意两条边都不相交。最后要计算这条路线包围的面积。

【输入格式】
第一行整数 n (3 <= n <= 1000),表示有n个点。
下来n行，每行两个整数x（横坐标）和y（纵坐标），表示点坐标（-10000<x,y<=10000）。
【输出格式】
一行一个实数，即封闭折线所包围的面积（保留4位小数）。
【样例1输入】
4
2 1
5 1
5 5
2 5
【样例1输出】
12.0000
【样例2输入】
5
2 1
5 1
3 2
5 3
2 3
【样例2输出】
4.0000

这道题对比上一道题就简单很多了，这道题要用到叉积的第二种用法——计算面积！对于任意一个多边形，只要以一号点p1为参考点，枚举任意两个相邻的点pi和pi-1（i>=3）,带符号计算p1和pi、pi-1所构成的三角形的面积（multi（pi，pi-1，p1）/2），累加取绝对值就是答案了。

至于为什么要带符号运算，请看我用下面这幅图演示

对于这个多边形，首先先加上multi(p3,p2,p1)/2的值（由于p3在p2的逆时针方向，所以为正）即为加上了红色部分

然后再加上multi(p4,p3,p1)/2的值（由于p4在p3的顺时针方向，所以为负）即为减去了绿色部分

然后再加上multi(p5,p4,p1)/2的值（由于p5在p4的逆时针方向，所以为正)即为加上了黄色部分，就是答案了

最后附上代码

#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;
struct node {
    double x, y;
} a[11000];

int i, j, k, m, n, o, p, js, jl;
double ans;

double multi(node p1, node p2, node p0) {
    double x1 = p1.x - p0.x;
    double y1 = p1.y - p0.y;
    double x2 = p2.x - p0.x;
    double y2 = p2.y - p0.y;
    return (x1 * y2 - x2 * y1);
}

int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; i++)scanf("%lf%lf", &a[i].x, &a[i].y);

    for (int i = 2; i <= n; i++) ans = ans + multi(a[i], a[i - 1], a[1]);
    ans = ans / 2.0;
    if (ans < 0)ans = -ans;

    printf("%0.4lf", ans);
    return 0;
}

结语

通过上面的学习和例题，是不是已经初步了解了计算几何的一点点精（皮）髓（毛）呢，希望你能有所收获，在计算几何的路上越走越远吧（蒟蒻的我可能是走不下去了）
最后Orz各位大佬，写的不好请各位神犇提出建议和意见，感谢你的阅读:neckbeard: