本文从最基本的线段相交问题出发,从解析几何进入计算几何,介绍点积和叉积这个最基本的计算几何工具,引入计算几何这个关于位置和方向的大航海世界~
本文要讲清楚的两个基本问题是:
考虑以下基本问题:
判断平面上两条线段是否相交
输入:4个点,分别表示第一条线段的两个端点和第二条线段的两个端点.
输出:Yes/No
线段相交,分为两种
回到本题,学习过大学解析几何(或者高中立体几何)的我们,立马就想到了写出直线方程,联立求解即可知道是否相交. 但是这种做法有诸多弊端. 至少有以下两点
可是,问题本身仅仅对相交与否感兴趣而已(虽然后续的计算几何的问题会涉及到求交点坐标), 于是,我们希望发展更为简洁高效的算法来解决这个问题.
首先,两条线段AB 和 CD相交等价于
A、B分别在 直线CD的两侧且C、D分别在直线 AB 的两侧. 形象的说,如果将A、B想象成人的两只脚,则就是人跨立在直线CD的两侧,类似的,C、D跨立在直线 AB 两侧. 而如何刻画两点跨立在直线两侧呢? 这就不得不引入一个极为有用的工具——叉积(或者叫外积).
令 , , 那么 和 的叉积定义如下
其中 det 是行列式.
叉积为什么这么重要呢? 因为它能有效刻画 和 之间的方位关系, 确切讲, 在 的逆时针旋转方向的充要条件为
其实这也就是高中时学过的右手定则.
那么 A、B两点位于直线 CD 两侧就很容易使用叉积等价刻画了. 图1中,A、B跨立在直线CD 两侧的充要条件就是
类似的,C、D跨立在直线 AB 两侧的充要条件是
上面两个不等式被形象的称为跨立实验(cross test)
跨立实验能帮助我们知道两条线段是否规范相交,那么非规范相交怎么处理呢?
非规范相交有以下两种情况
情况1: 线段重合, 线段有无穷多个交点.
情况2: 线段有唯一交点,但是此交点恰好是其中一条线段的一个端点.
无论哪种情况,上面的跨立实验中涉及到的4个外积中至少会有一个为0. 但是仅仅为0就能判定非规范相交了吗? 非也非也,例如
上图中A、B、C、D四点共线. 这种情况会使得跨立实验中涉及的四个叉积全部为0,但是线段AB 和线段CD 依旧不相交. 所以还必须做快速排斥实验. 具体见下面的伪代码.
给出判断线段相交的伪代码如下
#define cp const Point
struct Point
{
double x, y;
Point(double x = 0, double y = 0):x(x), y(y){}
Point operator - (cp &o) const
{
return Point(x - o.x, y - o.y);
}
double operator / (cp &o) const
{
return x * o.y - y * o.x;
}
};
ild direction(cp &p1, cp &p2, cp &p3)
{
return (p2 - p1) / (p3 - p1);
}
ili onsegment(cp &p1, cp &p2, cp &p3) // 快速排斥
{
return min(p1.x, p2.x) <= p3.x && p3.x <= max(p1.x, p2.x) && min(p1.y, p2.y) <= p3.y && p3.y <= max(p1.y, p2.y);
}
ili intersect(cp &p1, cp &p2, cp &p3, cp &p4) // 线段p1p2 和 线段p3p4 是否相交, 相交返回1, 否则返回0
{
double d1 = direction(p1, p2, p3);
double d2 = direction(p1, p2, p4);
double d3 = direction(p3, p4, p1);
double d4 = direction(p3, p4, p2);
if (d1 * d2 < 0 && d3 * d4 < 0) return 1; // 跨立实验
if (fabs(d1) < eps && onsegment(p1, p2, p3)) return 1;
if (fabs(d2) < eps && onsegment(p1, p2, p4)) return 1;
if (fabs(d3) < eps && onsegment(p3, p4, p1)) return 1;
if (fabs(d4) < eps && onsegment(p3, p4, p2)) return 1;
return 0;
}
特别的,如果我们想要判断的是直线 p1p2 和 线段 p3p4是否相交的话,则伪代码精简为
ili intersect(cp &p1, cp &p2, cp &p3, cp &p4) // 直线p1p2 和 线段p3p4 是否相交, 相交返回1, 否则返回0
{
double d1 = direction(p1, p2, p3);
double d2 = direction(p1, p2, p4);
if (d1 * d2 < 0) return 1;
if (fabs(d1) < eps) return 1;
if (fabs(d2) < eps) return 1;
return 0;
}
比较一下上述计算几何解法和我们伊始想用的解析几何解法,我们会发现计算几何的一个巨大的好处——不涉及三角函数以及除法,仅仅涉及加减法和乘法,性能高,精度好,而且不会遗漏任何特殊情况.
好了,讲清楚了判断线段相交的问题,进一步的问题就是计算交点坐标.
已知平面上两直线 L1(P, u), L2(Q, v) 相交,且恰有一个交点, 试计算该交点坐标.
其中 u 是 L1 的方向向量, v 是 L2 的方向向量.
首先说答案, 真的是一个十分精简的公式哟~
至于为什么,可以使用叉积为0得到以下方程
然后利用坐标展开,使用线性代数中的 Crammer 法则就可以得到上面的公式.
计算相交线段(或者直线)交点的坐标的伪代码如下
Point getintersect(cp &p, cp &u, cp &q, cp &v)
{
return p + sc(u, ((q - p) / v) / (u / v)); // sc 是 scale 的意思, 即向量数乘
}
最后,我们来看一道经典的例题 Poj 1039 Pipe
有宽度为1的折线管道,如下图所示. 上管道各个顶点的坐标为 (x1,y1),...,(xn, yn)
下管道相应的顶点为 (x1,y1-1),...,(xn, yn-1),管道的材料是不透明也不能反射光线的.
光线从左边入口的(x1,y1),(x1,y1-1)之间射入,各种方向进行直线传播,
问光线最远能射到哪里(求出x坐标)或者能穿透整个管道.
输入
多样例. 每个样例开始于 n (2<=n<=20), 然后 n行,每行两个浮点数 xi yi, 输入以 n = 0 结束.
输出
每个样例输出最远点的 x 坐标(保留两位小数), 如果光线能穿透整个管道, 打印
Through all the pipe.
样例输入
4
0 1
2 2
4 1
6 4
6
0 1
2 -0.6
5 -4.45
7 -5.57
12 -10.8
17 -16.55
0
样例输出
4.67
Through all the pipe.
限制
1s 10MB
来源
CEOI 1995
此题咋一看,无从下手——因为我们总不能枚举入射光线的角度吧?
但是 n 很小,所以本题各种暴力枚举都是允许的. 因为两点确定一条直线,所以每次枚举2个点,那么这两个点是哪两个点呢? 这就需要我们大力 YY 一下最后射的最远的光线的长相——一定是擦过一个管道的上顶点并且擦过一个管道的下顶点的光线. 于是我们就知道了,每次只需要枚举一个管道的上顶点和枚举一个管道的下顶点,这样就将光线确定下来了. 然后再去验证这条直线是否和线段 相交. 如果相交的话,则按照 的顺序去验证光线是否和垂直线段 相交. 假设 是第一个使得光线和 不相交的 , 则这说明光线不是和
通往 的上管道部分相交,就是和下管道部分相交. 就像下图一样. 当然,如果这个最小的 的话,则光线可以穿透所有管道.
于是可以写下如下代码
//#include "stdafx.h"
//#define LOCAL
#pragma GCC optimize(2)
#pragma G++ optimize(2)
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include <stdio.h>
#include <iostream>
#include <iomanip>
#include <string>
#include <ctype.h>
#include <string.h>
#include <fstream>
#include <sstream>
#include <math.h>
#include <map>
//#include <unordered_map>
#include <algorithm>
#include <vector>
#include <stack>
#include <queue>
#include <set>
#include <time.h>
#include <stdlib.h>
#include <bitset>
using namespace std;
//#define int unsigned long long
//#define int long long
#define re register int
#define ci const int
#define ui unsigned int
typedef pair<int, int> P;
#define FE(cur) for(re h = head[cur], to; ~h; h = g[h].nxt)
#define ilv inline void
#define ili inline int
#define ilc inline char
#define ild inline double
#define ilp inline P
#define LEN(cur) (hjt[cur].r - hjt[cur].l)
#define MID(cur) (hjt[cur].l + hjt[cur].r >> 1)
#define SQUARE(x) ((x) * (x))
typedef vector<int>::iterator vit;
typedef set<int>::iterator sit;
typedef map<int, int>::iterator mit;
const int inf = ~0u>>1;
const double PI = acos(-1.0), eps = 1e-8;
namespace fastio
{
const int BUF = 1 << 21;
char fr[BUF], fw[BUF], *pr1 = fr, *pr2 = fr;int pw;
ilc gc() { return pr1 == pr2 && (pr2 = (pr1 = fr) + fread(fr, 1, BUF, stdin), *pr2 = 0, pr1 == pr2) ? EOF : *pr1++; }
ilv flush() { fwrite(fw, 1, pw, stdout); pw = 0; }
ilv pc(char c) { if (pw >= BUF) flush(); fw[pw++] = c; }
ili read(int &x)
{
x = 0; int f = 1; char c = gc(); if (!~c) return EOF;
while(!isdigit(c)) { if (c == '-') f = -1; c = gc(); }
while(isdigit(c)) x = (x << 3) + (x << 1) + (c ^ 48), c = gc();
x *= f; return 1;
}
ili read(double &x)
{
int xx = 0; double f = 1.0, fraction = 1.0; char c = gc(); if (!~c) return EOF;
while (!isdigit(c)) { if (c == '-') f = -1.0; c = gc(); }
while (isdigit(c)) { xx = (xx << 3) + (xx << 1) + (c ^ 48), c = gc(); }
x = xx;
if (c ^ '.') { x = f * xx; return 1; }
c = gc();
while (isdigit(c)) x += (c ^ 48) * (fraction /= 10), c = gc();
x *= f; return 1;
}
ilv write(int x) { if (x < 0) pc('-'), x = -x; if (x > 9) write(x / 10); pc(x % 10 + 48); }
ilv writeln(int x) { write(x);pc(10); }
ili read(char *x)
{
char c = gc(); if (!~c) return EOF;
while(!isalpha(c) && !isdigit(c)) c = gc();
while (isalpha(c) || isdigit(c)) *x++ = c, c = gc();
*x = 0; return 1;
}
ili readln(char *x)
{
char c = gc(); if (!~c) return EOF;
while(c == 10) c = gc();
while(c >= 32 && c <= 126) *x++ = c, c = gc();
*x = 0; return 1;
}
ilv write(char *x) { while(*x) pc(*x++); }
ilv writeln(char *x) { write(x); pc(10); }
ilv write(char c) { pc(c); }
ilv writeln(char c) { write(c); pc(10); }
} using namespace fastio;
#define cp const Point
const int maxn = 105;
int n;
struct Point
{
double x, y;
Point(double x = 0, double y = 0):x(x), y(y){}
Point operator + (cp &o) const
{
return Point(x + o.x, y + o.y);
}
Point operator - (cp & o) const
{
return Point(x - o.x, y - o.y);
}
double operator / (cp &o) const
{
return x * o.y - y * o.x;
}
Point sc(double a) const
{
return Point(a * x, a * y);
}
};
struct Line
{
Point s, e;
} ls[maxn];
ilv mmax(double &a, double b)
{
if (a < b)
{
a = b;
}
}
ild direction(cp &p1, cp &p2, cp &p3)
{
return (p2 - p1) / (p3 - p1);
}
ili intersect(cp &p1, cp &p2, cp &p3, cp &p4)
{
double d1 = direction(p1, p2, p3);
double d2 = direction(p1, p2, p4);
if (d1 * d2 < 0) return 1;
if (fabs(d1) < eps) return 1;
if (fabs(d2) < eps) return 1;
return 0;
}
Point getintersect(cp &p, cp &u, cp &q, cp &v)
{
return p + u.sc(((q - p) / v) / (u / v));
}
ild tt(cp &s, cp &e)
{
int i = 1;
while (intersect(s, e, ls[i].s, ls[i].e) && i <= n) ++i;
if (i == 1) return -inf;
else if (i == n + 1) return inf;
else return max(getintersect(s, e - s, ls[i].s, ls[i].s - ls[i - 1].s).x, getintersect(s, e - s, ls[i].e, ls[i].e - ls[i - 1].e).x);
}
ild kk()
{
double ans = ls[1].s.x;
for (re i = 1; i < n; i++)
{
for (re j = i + 1; j <= n; j++)
{
mmax(ans, tt(ls[i].s, ls[j].e));
mmax(ans, tt(ls[i].e, ls[j].s));
}
}
return ans;
}
signed main()
{
#ifdef LOCAL
FILE *ALLIN = freopen("d:\\data.in", "r", stdin);
// freopen("d:\\my.out", "w", stdout);
#endif
int cnt = 0;
while (read(n), n)
{
for (re i = 1; i <= n; i++)
{
read(ls[i].s.x), read(ls[i].s.y);
ls[i].e.x = ls[i].s.x, ls[i].e.y = ls[i].s.y - 1;
}
double ans = kk();
fabs(ans - inf) < eps ? puts("Through all the pipe.") : printf("%.2lf\n", ans);
}
flush();
#ifdef LOCAL
fclose(ALLIN);
#endif
return 0;
}
ac情况
C++ Accepted 2208K 16MS