实用算法实践-第 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.youkuaiyun.com/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).