位置和方向的世界，计算几何的基本问题

ACM算法日常

发布于 2020-06-29 15:28:41

8710

发布于 2020-06-29 15:28:41

文章被收录于专栏：ACM算法日常

缘起

本文从最基本的线段相交问题出发，从解析几何进入计算几何，介绍点积和叉积这个最基本的计算几何工具，引入计算几何这个关于位置和方向的大航海世界~

分析

本文要讲清楚的两个基本问题是:

如何判断线段相交
进一步地，如果存在唯一交点，试求出相交的交点坐标

判断线段相交

考虑以下基本问题:

判断平面上两条线段是否相交

输入：4个点，分别表示第一条线段的两个端点和第二条线段的两个端点.

输出：Yes/No

线段相交，分为两种

规范相交，即两条线段交点恰有一个，而且该交点不是线段的任何一个端点. 例如

非规范相交，也就是不是"规范相交"的相交. 例如两条线段有重合部分或者唯一交点恰好是某条线段的一个端点. 例如（让我想起了GTA里面警察的警棍~）

回到本题，学习过大学解析几何（或者高中立体几何）的我们，立马就想到了写出直线方程，联立求解即可知道是否相交. 但是这种做法有诸多弊端. 至少有以下两点

情形繁多，尤其是特殊情形. 例如就拿直线方程而言，如果你要写斜截式的话，就要考虑斜率是否存在的问题.
可能涉及到除法、甚至三角函数. 而我们知道计算机处理浮点数可能丢失精度. 所以我们首先应该树立一个基本观念就是——能避免除法或三角函数的使用就尽量避免.

可是，问题本身仅仅对相交与否感兴趣而已（虽然后续的计算几何的问题会涉及到求交点坐标）, 于是，我们希望发展更为简洁高效的算法来解决这个问题.

首先，两条线段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

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-06-27，如有侵权请联系 cloudcommunity@tencent.com 删除

bash

本文分享自 ACM算法日常微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

bash

登录后参与评论

0 条评论

热度