计算几何学(Computational Geometry)

2024-05-30 19:18
文章标签 计算 geometry 几何学 computational

本文主要是介绍计算几何学(Computational Geometry)，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

计算几何学(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)的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/1016758。

相关文章

poj 1113 凸包+简单几何计算

poj 1113 凸包+简单几何计算

题意： 给N个平面上的点，现在要在离点外L米处建城墙，使得城墙把所有点都包含进去且城墙的长度最短。 解析： 韬哥出的某次训练赛上A出的第一道计算几何，算是大水题吧。 用convexhull算法把凸包求出来，然后加加减减就A了。 计算见下图： 好久没玩画图了啊好开心。 代码： #include <iostream>#include <cstdio>#inclu
阅读更多...
uva 1342 欧拉定理（计算几何模板）

uva 1342 欧拉定理（计算几何模板）

题意： 给几个点，把这几个点用直线连起来，求这些直线把平面分成了几个。 解析： 欧拉定理： 顶点数 + 面数 - 边数= 2。 代码： #include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#include <cstring>#include <cmath>#inc
阅读更多...
uva 11178 计算集合模板题

uva 11178 计算集合模板题

题意： 求三角形行三个角三等分点射线交出的内三角形坐标。 代码： #include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#include <cstring>#include <cmath>#include <stack>#include <vector>#include <
阅读更多...
XTU 1237 计算几何

XTU 1237 计算几何

题面： Magic Triangle Problem Description: Huangriq is a respectful acmer in ACM team of XTU because he brought the best place in regional contest in history of XTU. Huangriq works in a big compa
阅读更多...
音视频入门基础：WAV专题（10）——FFmpeg源码中计算WAV音频文件每个packet的pts、dts的实现

音视频入门基础：WAV专题（10）——FFmpeg源码中计算WAV音频文件每个packet的pts、dts的实现

一、引言 从文章《音视频入门基础：WAV专题（6）——通过FFprobe显示WAV音频文件每个数据包的信息》中我们可以知道，通过FFprobe命令可以打印WAV音频文件每个packet（也称为数据包或多媒体包）的信息，这些信息包含该packet的pts、dts： 打印出来的“pts”实际是AVPacket结构体中的成员变量pts，是以AVStream->time_base为单位的显
阅读更多...
计算数组的斜率，偏移，R2

计算数组的斜率，偏移，R2

模拟Excel中的R2的计算。         public bool fnCheckRear_R2(List<double[]> lRear, int iMinRear, int iMaxRear, ref double dR2)         {             bool bResult = true;             int n = 0;             dou
阅读更多...
GPU 计算 CMPS224 2021 学习笔记 02

GPU 计算 CMPS224 2021 学习笔记 02

并行类型 （1）任务并行 （2）数据并行 CPU & GPU CPU和GPU拥有相互独立的内存空间，需要在两者之间相互传输数据。 （1）分配GPU内存 （2）将CPU上的数据复制到GPU上 （3）在GPU上对数据进行计算操作 （4）将计算结果从GPU复制到CPU上 （5）释放GPU内存 CUDA内存管理API （1）分配内存 cudaErro
阅读更多...
Java - BigDecimal 计算分位(百分位)

Java - BigDecimal 计算分位(百分位)

日常开发中，如果使用数据库来直接查询一组数据的分位数，就比较简单，直接使用对应的函数就可以了，例如：         PERCENT_RANK() OVER(PARTITION BY 分组列名 ORDER BY 目标列名) AS 目标列名_分位数         如果是需要在代码逻辑部分进行分位数的计算，就需要我们自己写一个工具类来支持计算了 import static ja
阅读更多...
OpenStack离线Train版安装系列—2计算节点-环境准备

OpenStack离线Train版安装系列—2计算节点-环境准备

本系列文章包含从OpenStack离线源制作到完成OpenStack安装的全部过程。 在本系列教程中使用的OpenStack的安装版本为第20个版本Train（简称T版本），2020年5月13日，OpenStack社区发布了第21个版本Ussuri（简称U版本）。 OpenStack部署系列文章 OpenStack Victoria版 安装部署系列教程 OpenStack Ussuri版
阅读更多...
新一代车载(E/E)架构下的中央计算载体---HPC软件架构简介

新一代车载(E/E)架构下的中央计算载体---HPC软件架构简介

老规矩，分享一段喜欢的文字，避免自己成为高知识低文化的工程师： 屏蔽力是信息过载时代一个人的特殊竞争力，任何消耗你的人和事，多看一眼都是你的不对。非必要不费力证明自己，无利益不试图说服别人，是精神上的节能减排。 无人问津也好,技不如人也罢,你都要试着安静下来,去做自己该做的事.而不是让内心的烦躁、焦虑、毁掉你本就不多的热情和定力。 时间不知不觉中，快要来到夏末秋初。一年又过去了一大半，成
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2