计算几何

把草稿纸上的集合体转变成代码其中还是有很多细节问题需要注意的，《算法导论》主要介绍了叉积在几何学中的应用。向量的两个运算点积与叉积,点积：ab = |a||b|cos 叉积：a x b = |a||b|sin,由于叉积的形式与面积公式相似，也常用于计算多边形的面积，同时如果a x b > 0则a 向量在b向量的顺时针方向上，利用这个性质，《算法导论》介绍了凸包的求解方法。

HDU2036就是利用叉积求多边形面积的例子：

 1 #include<stdio.h>
 2 #include<math.h>
 3 
 4 #define MAXN 105
 5 
 6 struct POINT
 7 {
 8     double x,y;
 9 }p[MAXN];
10 
11 double polygon_area(POINT p[],int n);
12 
13 int main()
14 {
15     int n;
16 
17     //freopen("Sample Input.txt","r",stdin);
18 
19     while(scanf("%d",&n) != EOF && n)
20     {
21         for(int i = 0;i < n;i++)
22         {
23             scanf("%lf%lf",&p[i].x,&p[i].y);
24         }
25 
26         printf("%.1f\n",polygon_area(p,n));
27     }
28 }
29 
30 double polygon_area(POINT p[],int n)
31 {
32     double s = 0;
33 
34     for(int i = 0;i < n;i++)
35     {
36         s += (p[i].x * p[(i + 1) % n].y - p[i].y * p[(i + 1) % n].x);
37     }
38 
39     return fabs(s / 2.0);
40 }

 

double polygon_area(POINT p[],int n)就是将多边形以原点为定点分成小三角形利用叉积求解的例子，其中输入点对为逆时针方向，如果杂乱则可用下面的方法排序。

HDU1392介绍了利用叉积求解凸包的方法：

  1 #include<stdio.h>
  2 #include<algorithm>
  3 #include<math.h>
  4 
  5 #define MAXN 105
  6 
  7 
  8 struct POINT
  9 {
 10     int x,y;
 11 }trees[MAXN];
 12 
 13 int cmp(const void * a,const void * b);
 14 double get_dis(POINT a,POINT b);
 15 int cross_product(POINT a,POINT b,POINT c);
 16 int graham_scan(POINT trees[],POINT ch[],int n);
 17 
 18 int main()
 19 {
 20     int n;
 21 
 22     //freopen("Sample Input.txt","r",stdin);
 23 
 24     while(scanf("%d",&n) && n)
 25     {
 26         double dis = 0;
 27 
 28         for(int i = 0;i < n;i++)
 29         {
 30             scanf("%d%d",&trees[i].x,&trees[i].y);
 31         }
 32         
 33         if(n == 2)                                        
 34         {
 35             dis = get_dis(trees[0],trees[1]);
 36         }
 37         
 38         if(n >= 3)
 39         {
 40             POINT ch[MAXN];
 41             int cnt = graham_scan(trees,ch,n);
 42             for(int i = 0;i < cnt - 1;i++)
 43             {
 44                 dis += get_dis(ch[i],ch[i + 1]);
 45             }
 46             dis += get_dis(ch[cnt - 1],ch[0]);
 47         }
 48 
 49         printf("%.2f\n",dis);
 50     }
 51 }
 52 
 53 int cross_product(POINT a,POINT b,POINT c)
 54 {
 55     return (a.x - c.x) * (b.y - c.y) - (a.y - c.y) * (b.x - c.x);
 56 }
 57 
 58 int cmp(const void * a,const void * b)
 59 {
 60     double k = cross_product(*(POINT *)a,*(POINT *)b,trees[0]);
 61     if(k < 0 || (!k && get_dis(*(POINT *)a,trees[0]) > get_dis(*(POINT *)b,trees[0])))return 1;                    //不容易删除的情况下折中的办法
 62     return -1;
 63 }
 64 
 65 double get_dis(POINT a,POINT b)
 66 {
 67     return sqrt(pow(double(a.x - b.x),2) + pow(double(a.y - b.y),2));
 68 }
 69 
 70 void swap(POINT * a,POINT * b)
 71 {
 72     POINT t = * a;
 73     * a = * b;
 74     * b = t;
 75 }
 76 
 77 int graham_scan(POINT trees[],POINT ch[],int n)
 78 {
 79     int min = 0;
 80     for(int i = 0;i < n;i++)
 81     {
 82         if(trees[i].y < trees[min].y || (trees[i].y == trees[min].y && trees[i].x < trees[min].x))
 83         {
 84             min = i;
 85         }
 86     }
 87     swap(&trees[0],&trees[min]);
 88 
 89     qsort(trees + 1,n - 1,sizeof(trees[0]),cmp);
 90 
 91     ch[0] = trees[0];
 92     ch[1] = trees[1];
 93     ch[2] = trees[2];
 94     int top = 2;
 95 
 96     for(int i = 3;i < n;i++)
 97     {
 98         while(cross_product(trees[i],ch[top],ch[top - 1]) > 0 && top >= 2)
 99         {
100             top--;
101         }
102         ch[++top]=trees[i];
103     }
104     
105     return top + 1;
106 }

 

这里利用扫描法解决问题，先讲点按逆时针排序，然后扫描新点并检查以前的点最终的到结果，每一个状态都能保证暂时是凸包。


计算几何学博大精深啊，以后还会选修《计算机图形学》，里边应该会有很多有趣的东西，期待。。。

 

 posted on 2013-03-17 13:49  莫扎特的代码  阅读(195)  评论(0)  编辑  收藏  举报