实用算法实现-第 29 篇 计算几何学

29.1 最近点对

“最近”是指通常意义上的欧几里德距离。《算法导论》中对这个问题进行了介绍。[i]问对于经典的分治法求最近点对有深入的介绍和详尽的伪代码。

29.1.1实例

PKU JudgeOnline, 3714, Raid.

29.1.2问题描述

给定n个岗位和n个战士的坐标(0 ≤ X ≤1000000000, 0 ≤ Y ≤ 1000000000)。问这些战士当中，离某个岗位最近的距离是多少？

29.1.3输入

2

4

00

01

10

11

22

23

32

33

4

00

00

00

00

00

00

00

0 0

29.1.4输出

1.414

0.000

29.1.5分析

这也是一个最近距离的问题，只不过是两个集合中的元素的最近距离的问题。采用的方法不同于经典的最近距离点对的分治算法，但是思想上也有类似的地方。

对于每个战士(x0, y0)，首先找到一个点，这个点是x最接近x0的所有点中y最接近y0的点，求出(x, y)、(x0, y0)两点距离dis。若minDis > dis，令minDis = dis。在[x - minDis, x + minDis]的横坐标范围内移动x，并对每个x，寻找y最接近y0的点(x, y)，若minDis > dis，令minDis = dis。

上述描述中，“最接近x0的x”或者“最接近y0的y”都通过对有序数组的二分查找法来找出。

29.1.6程序

#include<stdio.h> #include<string.h> #include <stdlib.h> #include <math.h> #define ONLINE_JUDGE 0 #define maxNum 100001 #define MAX 0x7FFFFFFF #define __int64MAX 0x7FFFFFFFFFFFFFFF struct position { int x; int y; }; int cmp(const void*p1, const void*p2) { int y1; int y2; int x1; int x2; y1 = (*(position *)p1).y; y2 = (*(position *)p2).y; x1 = (*(position *)p1).x; x2 = (*(position *)p2).x; if(x1 >x2) { return1; }else if(x1 < x2){ return-1; }else{ if(y2< y1){ return1; }else{ return-1; } } } positionstation[maxNum]; positionagent[maxNum]; __int64 minDis; struct coordinate{ int pos; int x; }; int findMinY(int begin, int end, int t) { int y; int low,high, mid; low = begin; high = end - 1; while(low< high){ mid = (low + high) / 2; if(t> station[mid].y){ low = mid + 1; }else{ high = mid; } } //现在有关系: t< station[high].y && t > station[high - 1].y y = station[high].y - t; //如果直接check[high+ 1].pos可能会越界，因为high有可能等于top if(high> begin && y > t - station[high - 1].y) { y = t - station[high - 1].y; } return y; } coordinatecheck[maxNum]; int top; int N; void cal(position now) { int low,high, mid; int t; int x; int y; int temp; int begin; int end; int pos; __int64bigTemp; t = now.x; low = 0; high = top; while(low< high){ mid = (low + high) / 2; if(t> check[mid].x){ low = mid + 1; }else{ high = mid; } } //现在有关系: t< check[high].x && t > check[high - 1].x x = check[high].x - t; begin = check[high].pos; end = high + 1; pos = high; //如果直接check[high+ 1].pos可能会越界，因为high有可能等于top if(high> 0 && x > t - check[high - 1].x) { pos = high - 1; x = t - check[pos].x; begin = check[pos].pos; end = high; } end = check[end].pos; //这下即使士兵的横坐标大于所有电站的横坐标，check也不可能越界了， //因为只要top不为（即电站数目不为，这是可以保证的）， //就存在一个电站，使得电站与士兵的横坐标之差小于MAX - t的距离。 y = findMinY(begin, end, now.y); bigTemp = ((__int64)x)* ((__int64)x)+ ((__int64)y)* ((__int64)y); if(minDis> bigTemp){ minDis = bigTemp; } temp = pos - 1; while(temp>= 0){ x = now.x - check[temp].x; bigTemp = ((__int64)x)* ((__int64)x); if(bigTemp> minDis) { break; } begin = check[temp].pos; end = check[temp + 1].pos; y = findMinY(begin, end, now.y); bigTemp += ((__int64)y)* ((__int64)y); if(minDis> bigTemp){ minDis = bigTemp; } temp--; } temp = pos + 1; while(temp< top){ x = check[temp].x - now.x; bigTemp = ((__int64)x)* ((__int64)x); if(bigTemp> minDis) { break; } begin = check[temp].pos; end = check[temp + 1].pos; y = findMinY(begin, end, now.y); bigTemp += ((__int64)y)* ((__int64)y); if(minDis> bigTemp){ minDis = bigTemp; } temp++; } } int main() { #ifndef ONLINE_JUDGE FILE *fin; fin = freopen("test.txt", "r", stdin); if( !fin ) { printf("reopen in file failed...\n"); while(1){} return 0; } freopen("out.txt", "w", stdout); #endif int cases; int i; int temp; scanf("%d",&cases); for(; cases> 0; cases--){ scanf("%d",&N); for(i =0; i < N; i++){ scanf("%d%d",&station[i].x, &station[i].y); } for(i =0; i < N; i++){ scanf("%d%d",&agent[i].x, &agent[i].y); } qsort(station, N, sizeof(station[0]), cmp); top = 0; temp = -MAX; for(i =0; i < N; i++){ if(temp!= station[i].x) { temp = station[i].x; check[top].x = temp; check[top].pos = i; top++; } } check[top].x = MAX; check[top].pos = N; //check设立的最后这个元素，有诸多好处 minDis = __int64MAX; for(i =0; i < N; i++){ cal(agent[i]); } //printf("%d\n",minDis); printf("%.3f\n",sqrt((double)minDis)); } #ifndef ONLINE_JUDGE // fclose(stdout); fclose(stdin); #endif return 1; } 29.2 线段相交

关于线段相交，《算法艺术与信息学竞赛》介绍得比较细致。

29.2.1实例

PKU JudgeOnline, 1039, Pipe.

29.2.2问题描述

如上图所示，有一条管子，拐角处的两点坐标为(x, y)和(x, y + 1)。有一束光从入口射入，遇到管壁则被阻挡。求光所能达到的最大的x。

29.2.3输入

4

01

22

41

64

6

01

2-0.6

5-4.45

7-5.57

12-10.8

17-16.55

0

29.2.4输出

4.67

Through all the pipe.

29.2.5分析

《算法艺术与信息学竞赛》对这个题有结题报告。
29.2.6程序

#include <stdio.h> #include <math.h> typedef struct { double x,y; }point; int dblcmp (double d) //判断d的符号 { if (fabs(d) < 1e-6) return 0; return (d> 0) ? 1 : -1; } void getline (point a, point b, double&k, double &t, double&c) { //根据点a，b求直线kx+ty+c=0 k = a.y - b.y; t = b.x - a.x; c = a.x*b.y - b.x*a.y; return ; } int cal (double k, doublet, double c, point a) {//计算点a在直线kx+ty+c=0的上方还是下方还是直线上 returndblcmp (k*a.x + t*a.y + c); } int n; const int size = 21; pointup[size], down[size]; double max; void updata (double k, double t, double c,point a, point b) {//更新最大值 double k2,t2, c2; getline (a, b, k2, t2, c2); double xp =(t*c2 - t2*c) / (k*t2 - t*k2); if (xp >max) { max = xp; } return ; } bool judge (point a, point b) { // 判断点a，b构成的直线是否能通过整条管道 double k,t, c; getline (a, b, k, t, c); //判断光线能够从第一个通道口射入 if (cal(k,t, c, up[0]) * cal(k, t, c, down[0]) > 0) returnfalse; int i; //寻找第一个拐角，在这个拐角前，光线已穿越管壁 for (i = 1;i < n; i++) if(cal(k, t, c, up[i]) * cal(k, t, c, down[i]) > 0) break; if (i == n)return true; if (cal(k,t, c, up[i]) * cal(k, t, c, up[i-1]) <= 0) updata (k, t, c, up[i], up[i-1]); if (cal(k,t, c, down[i]) * cal(k, t, c, down[i-1]) <= 0) updata (k, t, c, down[i],down[i-1]); return false; } bool solve () { for (int i = 0; i < n; i++) for (int j = 0; j < n; j++){ if(i == j) continue; if(judge (up[i], down[j])) return true; } return false; } int main () { // freopen("in.txt", "r", stdin); while(scanf ("%d", &n), n){ for (int i = 0; i < n; i++){ scanf ("%lf%lf",&up[i].x, &up[i].y); down[i].x = up[i].x; down[i].y =up[i].y - 1; } max = up[0].x; if(solve ()) printf ("Through all thepipe.\n"); elseprintf ("%.2f\n", max); } }
29.3 Melkman在线凸包算法

《算法艺术与信息学竞赛》比较全面、生动地介绍了凸包算法及其应用。

Melkman凸包算法继承Graham-Scan算法的主要思想，并更近一步地采用双端队列，动态地在队列两头进行增删操作，维护“凸性”。

一个Melkman凸包算法的实例如下：

[ii]中描述的Melkman凸包算法如下所示：

1. t ← 1; b ← 0;//这里似乎Melkman弄反了，应该改成：t ← 0;b ←1;

vl← input; v2 ← input; v3 ← input;

if (Vl, v2, v3) > 0

then begin push vl; push v2; end;

else begin push v2; push v1; end;

pushv3 ; insert v3;

2. v ← input;

until (v, db, db+l) < 0 or (dt-1, dt, v) < 0

do v ← inputend;

3. until (dt-l,dt, v) > 0 do pop dt end;

pushv;

4. until(v, db, db+l) > 0 do remove db end;

insertv;

goto 2.

其中定义：

(x, y, z)为根据z在x到y的向量右边、同一条直线上、或左边，而返回1, 0, -1的函数。

push操作就是使得t = t + 1; dt = v。

pop操作就是使得 t = t - 1。

insert操作就是使得b = b - 1; db = v。

remove操作就是使得b = b + 1。

29.3.1实例

PKU JudgeOnline, 1113, Wall.

29.3.2问题描述

已知N个点的坐标，要将这些点用围墙围起来，并使得围墙距离点的距离不超过L。求围墙的周长。输出四舍五入。

先输入N和L，然后是N个点的坐标。

1.3.3输入

9100

200400

300400

300300

400300

400400

500400

500200

350200

200 200

29.3.4输出

1628

29.3.5分析

先求这些点的凸包。

可以很容易地证明：周长最短的围墙是由由于凸包的每条边，垂直向外平移L的距离，再将这些线段用半径为L的圆弧连接起来而构成的。

同时又很容易证明：这些圆弧的角度之和为360度。

故此这个问题直接转化为求凸包的周长问题。

凸包的求法中三点共线是十分需要注意的。Melkman算法描述中略去了对三点共线的许多细节的处理。

29.3.6程序

#include<stdio.h> #include<string.h> #include<math.h> #include<iostream> using namespace std; #define maxNum 2010 #define BEGIN 1005 #define PI 3.1415926 struct Vector{ int x; int y; }; Vectord[maxNum]; inline int crossProduct(Vector p1,Vector p2) { return p1.x* p2.y - p1.y * p2.x; } inline int fun(Vector v1,Vectorv2,Vector v3) { Vector tempv1; Vector tempv2; tempv1.x = v2.x - v1.x; tempv1.y = v2.y - v1.y; tempv2.x = v3.x - v1.x; tempv2.y = v3.y - v1.y; returncrossProduct(tempv1, tempv2); } inline double dis(Vector p1, Vectorp2) { double x; double y; x = p2.x - p1.x; y = p2.y - p1.y; x = x * x; y = y * y; returnsqrt(x + y); } inline int between(Vector p0,Vector p1, Vector p2) { int x1; int x2; if(p1.x> p2.x){ x1 = p1.x; x2 = p2.x; }else{ x1 = p2.x; x2 = p1.x; } if(p0.x< x1 && p0.x > x2){ return1; }else{ return0; } } #define push(v) t++;d[t]= v; #define pop() t--; #define insert(v) b--;d[b] =v; #define remove() b++; int N; void Melkman() { int t; int b; int i; double sum; int temp; Vector v, v1, v2, v3; int L; cin >> L; cin >> v1.x >> v1.y; cin >> v2.x >> v2.y; cin >> v3.x >> v3.y; t = 0; b = 1; t += BEGIN; b += BEGIN; if(fun(v1,v2, v3) > 0){ push(v1); push(v2); }else{ push(v2); push(v1); } push(v3); insert(v3); for(i = 3;i < N; i++){ cin >> v.x >> v.y; /* //如果有坐标相同的点，需要加入如下的判断 for(j = b; j< t; j++){ if(v.x ==d[j].x && v.y == d[j].y) { //int*a = (int *)0; //a =0; break; } } if(j < t) { continue; } */ //if(!(fun(v,d[b], d[b + 1]) < 0 || fun(d[t - 1], d[t], v) < 0)) //把共点的连同共线的都过滤了，过滤点的条件太松 if((fun(v,d[b], d[b + 1]) > 0 && fun(d[t - 1], d[t], v) > 0)) //把共点的连同共线的都没过滤，过滤点的条件太紧 { //这里必须如此，不能是上面的，否则共线的出错 continue; } if(fun(v,d[b], d[b + 1]) == 0 && fun(d[t - 1], d[t], v) == 0) { if(between(v,d[b], d[b + 1]) || between(v, d[t - 1], d[t])) { //这里进行的是将共线的，但是处在已得凸包边上的点滤除 continue; } } while(t- b >= 2 && fun(d[t - 1], d[t], v) <= 0){ //只要共线的都删除重新加了，直到只剩下两个点为止 pop(); // cout<< "poped"<< endl; } push(v); while(t- b >= 2 && fun(v, d[b], d[b + 1]) <= 0){ //只要共线的都删除重新加了，直到只剩下两个点为止 remove(); // cout<< "removed"<<endl; } insert(v); } sum = 0; for(i = b;i < t; i++){ sum += dis(d[i], d[i + 1]); } sum += 2 * PI * L; temp = (int)sum; sum -= temp; if(sum>= 0.5) { temp++; } cout << temp << endl; } int main() { while(cin>> N){ Melkman(); } }
29.4 实例

29.4.1线段相交的实例

PKU JudgeOnline, 1039, Pipe.

29.4.2凸包的实例

PKU JudgeOnline, 1113, Wall.

本文章欢迎转载，请保留原始博客链接http://blog.csdn.net/fsdev/article

---

[i] Discrete Mathematics, Sixth Edition. Richard Johnsonbaugh.

[ii] On-line Construction of the Convex Hull of a Simple Polygon. MelkmanA. Information Processing Letters 25, p11-12 , (1987).