计算几何学(Computational Geometry)

分类：算法学习 ｜ 作者：酷~行天下 ｜ 发表于2011/11/27 2条评论 3,299 views

看名字就知道是研究几何问题的算法。

一、线段算法基础

首先要讨论的是三个基础问题：

1）已知两条有向线段p0p1，和p0p2，相对于它们的公共端点p0来说，p0p1是否在p0p2的顺时针方向上？

2）已知两条线段p0p1和p1p2，如果先通过p0p1再通过p1p2，在点p1处是不是要向左旋？

3）线段p1p2和p3p4是否相交？

1）叉积(Cross Product)

叉积是线段算法的中心，图 a) 两个向量 p1，p2，可以把 p1 ×p2看作是由点（0, 0），p1，p2 和 p1 + p2 = (x1 + x2，y1 + y2) 所形成平行四边形的面积，另一种等价定义是把叉积定义为一个矩阵的行列式：

如果 p1×p2 为正数，则相对于原点来说，p1 在 p2 的顺时针方向上；如果 p1×p2 为负数，则 p1 在 p2 的逆时针方向上。如图b) 浅阴影区域包含了 p 的顺时针方向上的向量。深阴影区包含了 p 的逆时针方向上的向量。

对于上边提到的第一个问题，解决方法是，将p0当做原点即可，计算叉积：

如果叉积为政，则 p0p1 在 p0p2 的顺时针方向上，如果为负，则 p0p1 在 p0p2 的逆时针方向上。

2）确定线段旋转方向

第二个问题，如果知道叉积，那再简单不过了，相对于有向线段 p0p1，检查有向线段 p0p2 是在其顺时针方向还是逆时针方向即可，只需计算叉积（p2 – p0）× （p1 – p0）,如果叉积为负，则 p0p2 在 p0p1 逆时针方向，因此左旋，反之，同理。

3）判断两个线段是否相交

第三个问题，以下两个条件，至少满足一个，即可断定两线段相交：

1）每个线段都跨越（straddle）包含了另一线段的直线（即最普通的相交）

2）一个线段的某一个端点位于另一线段上（边界情况，包含重合）

SEGMENTS-INTERSECT 用于判断是否相交，DIRECTION利用叉积方法，计算线段方位；ON-SEGMENT确定一个与某一线段共线的点是否位于该线段上。伪代码：

SEGMENTS-INTERSECT(p1, p2, p3, p4)
1  d1 ← DIRECTION(p3, p4, p1)
2  d2 ← DIRECTION(p3, p4, p2)
3  d3 ← DIRECTION(p1, p2, p3)
4  d4 ← DIRECTION(p1, p2, p4)
5  if ((d1 > 0 and d2 < 0) or (d1 < 0 and d2 > 0)) and
((d3 > 0 and d4 < 0) or (d3 < 0 and d4 > 0))
6     then return TRUE
7  elseif d1 = 0 and ON-SEGMENT(p3, p4, p1)
8     then return TRUE
9  elseif d2 = 0 and ON-SEGMENT(p3, p4, p2)
10     then return TRUE
11  elseif d3 = 0 and ON-SEGMENT(p1, p2, p3)
12     then return TRUE
13  elseif d4 = 0 and ON-SEGMENT(p1, p2, p4)
14     then return TRUE

15  else return FALSE

DIRECTION(pi, pj, pk)

1  return (pk – pi) × (pj – pi)

ON-SEGMENT(pi, pj, pk)
1  if min(xi, xj) ≤ xk ≤ max(xi, xj) and min(yi, yj) ≤ yk ≤ max(yi, yj)
2     then return TRUE

3     else return FALSE

伪代码很容易看懂，判断过程，有四种情况需要选择，a）相交，叉积异号；b）不想交，叉积同号；c）p3和 p1，p2 共线，4）p3 和 p1,p2共线，但是不相交。（最容易忽略就是第四种情况）

代码实现：

+ expand source 帮助
+ expand source
01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
//新写一个，这个代码可以顺便AC掉 HDU 1086
#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
  
struct Point {
     double x, y;
};
  
double Direction(Point p1, Point p2, Point p0 ) {
     return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);
}
  
bool OnSegment( Point p1, Point p2, Point p0) {
     if (min(p1.x, p2.x) <= p0.x && p0.x <= max(p1.x, p2.x) &&
         min(p1.y, p2.y) <= p0.y && p0.y <= max(p1.y, p2.y)) {
             return true ;
     } else {
         return false ;
     }
}
  
bool SegmentsIntersert(Point p1, Point p2, Point p3, Point p4) {
     double d1 = Direction(p3, p4, p1);
     double d2 = Direction(p3, p4, p2);
     double d3 = Direction(p1, p2, p3);
     double d4 = Direction(p1, p2, p4);
  
     if (d1 * d2 < 0 && d3 * d4 < 0) {
         return true ;
     } else if (d1 == 0 && OnSegment (p3, p4, p1)) {
         return true ;
     } else if (d2 == 0 && OnSegment (p3, p4, p2)) {
         return true ;
     } else if (d3 == 0 && OnSegment (p1, p2, p3)) {
         return true ;
     } else if (d4 == 0 && OnSegment (p1, p2, p4)) {
         return false ;
     }
     return 0;
}
  
int main() {
     int n;
     Point a[105], b[105];
     while (~ scanf ( "%d" , &n), n) {
         int cnt = 0;
         for ( int i = 0; i < n; i++) {
             scanf ( "%lf %lf %lf %lf" , &a[i].x, &a[i].y, &b[i].x, &b[i].y);
         }
         for ( int i = 0; i < n; i++) {
             for ( int j = i + 1; j < n; j++) {
                     if (SegmentsIntersert (a[i], b[i], a[j], b[j])) {
                         cnt++;
                     }
             }
         }
         printf ( "%d\n" , cnt);
     }
     return 0;
}

二、寻找凸包(Convex Hull)

点集Q的凸包是一个最小的凸多边形P，可以把Q中的每个点都想象成是露在一块板外的铁钉，那么凸包就是包围了所有这些铁钉的一条拉紧了的橡皮绳所构成的形状。如图点集Q和凸包CH（Q）。

寻找凸包，可以使用Graham扫描法（Graham’s scan）。借助一个堆栈S，输入集合Q中的每个点都被压入栈一次，非 CH(Q) 中顶点的点最终将被弹出堆栈，算法结束时，堆栈S中仅包含CH（Q）中的顶点，其顺序为各点在边界上出现的逆时针方向排列的顺序。

伪代码：

GRAHAM-SCAN(Q)
1  let p0 be the point in Q with the minimum y-coordinate,
or the leftmost such point in case of a tie
2  let 〈p1, p2, …, pm〉 be the remaining points in Q,
sorted by polar angle in counterclockwise order around p0
(if more than one point has the same angle, remove all but
the one that is farthest from p0)
3  PUSH(p0, S)
4  PUSH(p1, S)
5  PUSH(p2, S)
6  for i ← 3 to m
7       do while the angle formed by points NEXT-TO-TOP(S), TOP(S),
and pi makes a nonleft turn
8              do POP(S)
9          PUSH(pi, S)

10  return S

第一行找到 y 坐标最小的 p0，第二行根据相对于p0的极角排序（利用叉积），3~5 行堆栈初始化，第 6~9行的迭代对一个点执行一次，意图是：在对点pi进行处理后，在栈S中，按由底到顶得顺序，包含CH（{p0，p1, ……，pi}）中按逆时针方向排列的各个顶点。第 7~8 行由 while 循环把发现不是凸包中的顶点从堆栈中移去。沿逆时针方向通过凸包时，在每个点处应该左转。因此while循环每次发现在一个顶点处没有向左转时，就把该顶点从堆栈中弹出。图示：

如图 a)→b) 过程中，原本p2在栈S中，到 b) 时，根据p1p3，p1p2的叉积计算，得知 p1p3 在 p1p2 的顺时针方向行，所以 p2 不在CH（Q）中，所以 p2 出栈，p3入栈，整个算法都循环这个过程，后续完整图示：

代码实现：

+ expand source 帮助
+ expand source
01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
//可以顺便AC掉  HDU 1392，去年的代码略加修改
#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
  
typedef struct {
     double x;
     double y;
}Point;
Point p[110], stack[110];
  
int N, top;
  
double Direction(Point p1, Point p2, Point p3) {
     return (p2.x - p1.x) * (p3.y - p1.y) - (p2.y - p1.y) * (p3.x - p1.x);
}
  
double Dis(Point a, Point b) {
     return sqrt ((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
  
int cmp( const void *a, const void *b) {
     Point c = *(Point *)a;
     Point d = *(Point *)b;
     double k = Direction(p[0], c, d);
     if (k < 0 || (!k && Dis(c, p[0]) > Dis(d, p[0]))) {
         return 1;
     }
     return -1;
}
  
void Convex() {
     for ( int i = 1; i < N; i++) {
         Point temp;
         if (p[i].y < p[0].y || (p[i].y == p[0].y && p[i].x < p[0].x)) {
             temp = p[i];
             p[i] = p[0];
             p[0] = temp;
         }
     }
     qsort (p + 1, N - 1, sizeof (p[0]), cmp);
     stack[0] = p[0];        
     stack[1] = p[1];
     stack[2] = p[2];
     top =2;
     for ( int i = 3; i < N; i++) {
         while (top >= 2 && Direction(stack[top - 1], stack[top], p[i]) <= 0) {
             top--;
         }
         top++;
         stack[top] = p[i];
     }
}
  
int main() {
     while (~ scanf ( "%d" , &N), N){
         for ( int i = 0; i < N; i++) {
             scanf ( "%lf %lf" , &p[i].x, &p[i].y);
         }
         cout.setf(ios::fixed);
         cout.precision(2);
  
         if (N == 1) {
             puts ( "0.00" );
             continue ;
         } else if (N == 2) {
             printf ( "%.2f\n" , Dis(p[0], p[1]));
             continue ;
         }
  
         Convex();
         double ans = 0;
         for ( int i = 0; i < top; i++) {
             ans += Dis(stack[i],stack[i + 1]);
         }
  
         ans += Dis(stack[top], stack[0]);
         printf ( "%.2f\n" , ans);
     }
     return 0;
}
        

三、寻找最近点对

有 n 个点的集合 Q，找出最近的两点，这一问题可以应用于交通控制等系统。在空中或海洋交通控制系统中，需要发现两个距离最近的交通工具，以便检测出可能发生的相撞事故。

暴力必然可以，但是时间复杂度太高，因为要查看 Cn2 = O（n2）个点对，算法导论中介绍的是一种分治策略，输入 P，X 和 Y 的递归调用首先检查是否 |P| <= 3。如果是，则用暴力枚举C|p|2个点对，并返回最近点对，如果 |P| > 3，则分治递归。

找出一条垂直线，近似平分，左边的为 PL，右边的为 PR，类似地X数组划分为两个 XL 和 XR，分别包含 PL和 PR 中的店，按 x 坐标递增排序。类似的，Y 数组划分为连个 YL 和 YR，分别包含 PL 和 PR 中的店，按 y 坐标递增排序。

把 P 分为 PL 和 PR 后，再两次递归调用，一次找出 PL 中的最近点对，另一次找出 PR 中的最近点对。第一次调用输入为 PL，XL，YL，第二次调用输入为 PR，XR，YR。设分别返回的最近距离为 δL 和 δR，则δ = min（δL，δR）

算出左右两边独自的最小距离后，还要计算跨域左右两边的最小距离，即，一个点在PL中，另外一个在PR中。比δ小，且在分布在两边，必定是以下情况：

即，必定在（-δ，+δ）内，如图a)，另外，最多有 8 个点可能处于 δ × 2δ 矩形区域内，原因是：比如左半边，因为 PL 中的所有点之间的距离至少为 δ 单位，所以至多有 4 个点可能位于该正方形内，PR 同理。这个结论的用处是，在暴力枚举可能范围内的点时，每个点只需考虑紧随其后的 7 个点即可。部分伪代码：

1  length[YL] ← length[YR] ← 0
2  for i ← 1 to length[Y]
3       do if Y[i] ∈ PL
4             then length[YL] ← length[YL] + 1
5                  YL[length[YL]] ← Y[i]
6             else length[YR] ← length[YR] + 1

7                  YR[length[YR]] ← Y[i]

代码实现：

+ expand source 帮助
+ expand source
01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cfloat>
using namespace std;
const int maxn = 10000;
  
struct Point {
     double x, y;
     int index;
}X[maxn], Y[maxn];
double ans;
  
bool Cmpx( const Point &a, const Point &b) {
     return a.x < b.x;
}
  
bool Cmpy( const Point &a, const Point &b) {
     return a.y < b.y;
}
  
double Dis( const Point &a, const Point &b) {
     return sqrt ((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
  
//暴力枚举
double ForceMinDis( int a, int b) {
     for ( int i = a; i < b; i++) {
         for ( int j = i + 1; j <= b; j++) {
             double temp = Dis(X[i], X[j]);
             if (ans > temp) {
                 ans = temp;
             }
         }
     }
     return ans;
}
  
double MinDis( int a, int b, Point *Y) {
     if (b - a <= 3) {
         return ForceMinDis(a, b);
     }
  
     int mid = (a + b) >> 1;
     Point *Yleft = new Point[mid - a + 1];
     Point *Yright = new Point[b - mid];
  
     for ( int i = 0, j = 0, k = 0; i <= b - a; i++) {
         if (Y[i].index <= mid) {
             Yleft[j++] = Y[i];
         } else {
             Yright[k++] = Y[i];
         }
     }
  
     double left_min = MinDis(a, mid, Yleft);
     double right_min = MinDis(mid + 1, b, Yright);
     ans = min(left_min, right_min);
  
     double line = X[mid].x;
     int ind;
     for ( int i = 0, ind = 0; i <= b - a; i++) {
         if ( fabs (Y[i].x - line) < ans) {
             Y[ind++] = Y[i];
         }
     }
     for ( int i = 0; i < ind - 1; i++) {
         for ( int j = i + 1; j <= i + 7 && j < ind; j++) {
             double temp = Dis(Y[i], Y[j]);
             if (ans > temp) {
                 ans = temp;
             }
         }
     }
     delete Yleft;
     delete Yright;
     return ans;
}
  
int main() {
     int n, i;
     while ( scanf ( "%d" , &n), n) {
         ans = DBL_MAX;
         for (i = 0; i < n; i++) {
             scanf ( "%lf %lf" , &X[i].x, &X[i].y);
         }
         sort(X, X + n, Cmpx);
         for (i = 0; i < n; i++) {
             X[i].index = i;
         }
         memcpy (Y, X, n * sizeof (Point));
         sort(Y, Y + n, Cmpy);
         printf ( "%lf\n" , MinDis(0, n - 1, Y));
     }
     return 0;
}

这里有个网页《计算几何算法概览》，太全面了

计算几何学(Computational Geometry)相关推荐

  1. 算法导论之计算几何学

    计算几何学是计算机科学的一个分支,专门研究集合问题的解决的算法.计算几何学的问题一般输入关于一组集合对象的描述,如一组点.一组线段:输出是对问题的回答,如直线是否相交.三维空间和高维空间很难视觉化,这 ...

  2. 计算几何学 | 知识点及C++代码实现汇总

    博客上看到的非常优质的计算几何学c++代码实现. 代码可直接访问github仓库 计算几何学 计算几何学 | 几何对象的基本元素与表现 | Basic 计算几何基础知识整理大全 | Vector | ...

  3. 计算几何学 | 逆时针方向 | Counter-Clockwise | C/C++实现

    问题描述 对于三个点p0.p1.p2,请按照下列情况进行输出: p0.p1.p2成逆时针方向 COUNTER_CLOCKWISE p0.p1.p2成顺时针方向 CLOCKWISE p2.p0.p1依次 ...

  4. 计算几何学 | 圆与直线的交点 | Cross Points of a Circle and a Line | C/C++实现

    问题描述 求圆c与直线 l l l的交点. 输入: 输入按照下述格式给出: c x cx cx c y cy cy r r r q q q L i n e 1 Line_1 Line1​ L i n ...

  5. 计算思维Computational Thinking

    计算思维 Jeannette M. Wing ( 周以真 )  ( 翻译:徐韵文,王飞跃 , 校对:王飞跃) 它代表着一种普遍的认识和一类普适的技能,每一个人,不仅仅是计算机科学家,都应热心于它的学习 ...

  6. Counter-Clockwise 三点位置关系 计算几何学

    题目链接:https://vjudge.net/problem/Aizu-CGL_1_C #include<bits/stdc++.h> using namespace std; #def ...

  7. Binary Search(二分搜索)

    转载请注明出处 leonchen1024.com/2018/08/14/- 二分搜索(binary search),也叫做 折半搜索(half-interval search),对数搜索(logari ...

  8. 卡内基·梅隆大学计算机科学系主任周以真的父母是中国人吗,计算思维（Computational Thinking）...

    最近两天的讲座中连续两天听到两位老师不约而同的提到一个他们认为非常重要的概念"计算思维". 资料:2006年3月,美国卡内基·梅隆大学计算机科学系主任周以真(Jeannette M ...

  9. PostGIS 距离计算规范 - 投影 与 球 坐标系, geometry 与 geography 类型

    平面坐标,球面坐标 SRID=2369  投影坐标 SRID=4326  球面坐标 ST_Transform(ST_GeomFromText('POINT(120.08 30.96)', 4326), ...

最新文章

  1. 招聘|追觅科技VSLAM​、CV算法实习生招聘
  2. python 开发版-MicroPython开发之物联网快速开发板
  3. 代码大全阅读笔记02
  4. ASP.NET AJAX深入浅出系列课程
  5. 记录本地Hexo博客部署到服务器上
  6. linux系统编程：IO读写过程的原子性操作实验
  7. Windows系统常用技巧总结
  8. 牛牛牛！干翻Sci-hub!
  9. Java——volatile关键字详解
  10. cx_Oracle.DatabaseError: DPI-1047: 64-bit Oracle Client library cannot be loaded 解决方法
  11. 教你win10系统无法识别语音识别的解决方法
  12. 开源工业物联网数据库 Apache IoTDB 毕业成为 Apache 顶级项目！
  13. 南大被骂到上热搜！Nature杂志回应南京大学拟花120万发校庆特刊！
  14. 《自抗扰控制技术》——第二遍（仿真）
  15. 【中科院】分子生物学-朱玉贤第四版-笔记-第10讲 分子生物学操作技术
  16. 安卓pdf阅读器_【软件分享】自用的一款PDF阅读器——悦书PDF阅读器，支持护眼模式、注释涂鸦、PDF转换，功能齐全，界面简洁美观。...
  17. 《逐梦旅程 WINDOWS游戏编程之从零开始》笔记6——四大变换光照与材质
  18. antd select.option选项加入额外属性
  19. oracle ebs form查询,Oracle EBS FORM 更改记录状态
  20. 如何从CRAN上下载R语言程序包

热门文章

  1. 红外光学雨量传感器的工作原理
  2. 【CTF】加密1——滴答~滴+聪明的小羊
  3. linux系统kate,Manjaro Linux 19.0 系统正式发布 对文本编辑器Kate进行补充
  4. 将安卓项目部署云服务器,将app项目部署到云服务器上
  5. 解决ie6、ie7下float为right换行的情况
  6. 攻克论文写作系列之2：阅读文献的三个进阶和文献笔记法
  7. 永洪科技王桐：人生要么是一段大胆的冒险，要么什么都不是
  8. linux PS命令
  9. JData数据处理及高潜用户购买意向预测
  10. JData大数据竞赛18年赛题-如期而至-用户购买时间预测