我的计算几何之路

            计算几何之路

 

计算几何 Part.1---点，线，面，形基本关系，点积叉积的理解

计算几何 Part.2---凸包问题

计算几何 Part.3---面积公式矩形切割

计算几何 Part.4---半平面交

计算几何 Part.5---计算几何背景，实际上解题的关键是其他问题（数据结构、组合数学，或者是枚举思想）若干道经典的离散化＋扫描线的题目

计算几何 Part.6---解析几何

计算几何 Part.7---旋转卡壳

计算几何 Part.8---极角排序

 

 

END：推荐BLOG

 

 

计算几何题的特点与做题要领：
1.大部分不会很难，少部分题目思路很巧妙
2.做计算几何题目，模板很重要，模板必须高度可靠。
3.要注意代码的组织，因为计算几何的题目很容易上两百行代码，里面大部分是模板。如果代码一片混乱，那么会严重影响做题正确率。
4.注意精度控制。
5.能用整数的地方尽量用整数，要想到扩大数据的方法（扩大一倍，或扩大sqrt2）。因为整数不用考虑浮点误差，而且运算比浮点快。

 

 

 

 

 

计算几何 Part.1---点，线，面，形基本关系，点积叉积的理解

1.

链接：http://poj.org/problem?id=2318

 

题意：给了m个点，落在n+1个区域中，问各个区域有多少个点。

就是利用叉积去判断点在线段的哪一侧，可以二分去做，比较快。

 

代码：

#include <iostream>

#include <stdio.h>

#include <string.h>

#include <algorithm>

#include <queue>

#include <map>

#include <vector>

#include <set>

#include <string>

#include <math.h>

 

using namespace std;

struct Point

{

   int x,y;

   Point(){}

   Point(int _x,int _y)

    {

       x = _x;y = _y;

    }

   Point operator -(const Point &b)const

    {

       return Point(x - b.x,y - b.y);

    }

   int operator *(const Point &b)const

    {

       return x*b.x + y*b.y;

    }

   int operator ^(const Point &b)const

    {

       return x*b.y - y*b.x;

    }

};

struct Line

{

   Point s,e;

   Line(){}

   Line(Point _s,Point _e)   

    {

       s = _s;e = _e;

    }

};

 

int xmult(Point p0,Point p1,Point p2) //计算p0p1 X p0p2

{

   return (p1-p0)^(p2-p0);

}

const int MAXN = 5050;

Line line[MAXN];

int ans[MAXN];

int main()

{

   //freopen("in.txt","r",stdin);

   //freopen("out.txt","w",stdout);

   int n,m,x1,y1,x2,y2;

   bool first = true;

   while(scanf("%d",&n) == 1 && n)

    {

       if(first)first = false;

       else printf("\n");

       scanf("%d%d%d%d%d",&m,&x1,&y1,&x2,&y2);

       int Ui,Li;

       for(int i = 0;i < n;i++)

       {

           scanf("%d%d",&Ui,&Li);

           line[i] = Line(Point(Ui,y1),Point(Li,y2));

       }

       line[n] = Line(Point(x2,y1),Point(x2,y2));

       int x,y;

       Point p;

       memset(ans,0,sizeof(ans));

       while( m-- )

       {

           scanf("%d%d",&x,&y);

           p = Point(x,y);

           int l = 0,r = n;

           int tmp;

           while( l <= r)

           {

                int mid = (l + r)/2;

               if(xmult(p,line[mid].s,line[mid].e) < 0)

                {

                    tmp = mid;

                    r = mid - 1;

                }

                else l = mid + 1;

           }

           ans[tmp]++;

       }

       for(int i = 0; i <= n;i++)

           printf("%d: %d\n",i,ans[i]);

    }

   return 0;

}

 

2. http://poj.org/problem?id=1269

 

先判断两条直线是不是同线，不是的话再判断是否平行，再不是的话就只能是相交的，求出交点。

如何判断是否同线？由叉积的原理知道如果p1，p2，p3共线的话那么(p2-p1)X(p3-p1)=0。因此如果p1，p2，p3共线，p1，p2，p4共线，那么两条直线共线。direction（）求叉积，叉积为0说明共线。

如何判断是否平行？由向量可以判断出两直线是否平行。如果两直线平行，那么向量p1p2、p3p4也是平等的。即((p1.x-p2.x)*(p3.y-p4.y)-(p1.y-p2.y)*(p3.x-p4.x))==0说明向量平等。

如何求出交点？这里也用到叉积的原理。假设交点为p0(x0,y0)。则有：

(p1-p0)X(p2-p0)=0

(p3-p0)X(p2-p0)=0

展开后即是

(y1-y2)x0+(x2-x1)y0+x1y2-x2y1=0

(y3-y4)x0+(x4-x3)y0+x3y4-x4y3=0

将x0,y0作为变量求解二元一次方程组。

假设有二元一次方程组

a1x+b1y+c1=0;

a2x+b2y+c2=0

那么

x=(c1*b2-c2*b1)/(a2*b1-a1*b2);

y=(a2*c1-a1*c2)/(a1*b2-a2*b1);

因为此处两直线不会平行，所以分母不会为0。

 

#include<iostream>

#include<stdio.h>

using namespace std;

const double epx=1e-10;

struct Point

{

    double x;

   double y;

};

//求解二元一次方程

Point solve(double a1,double b1,doublec1,double a2,double b2,double c2)

{

   Point p;

   p.x=(c1*b2-c2*b1)/(a2*b1-a1*b2);

   p.y=(a2*c1-a1*c2)/(a1*b2-a2*b1);

   return p;

}

//p1p3,p1p2的叉积

double direction(Point p1,Point p2,Pointp3)

{

   return (p3.x-p1.x)*(p2.y-p1.y)-(p2.x-p1.x)*(p3.y-p1.y);

}

//判断两直线的关系

/*

在这里有三种关系：1共线 2平行 3相交

1 共线可通过叉积来判断

2 平行通过向量来判断

3 通了上面2种情况的其他情况

  求交点可通过叉积及解二元一次方程来求解

*/

int N;

Point p1,p2,p3,p4;

Point p0;//交点

double a1,b1,c1,a2,b2,c2;

int main()

{

   scanf("%d",&N);

   int i;

   printf("INTERSECTING LINES OUTPUT\n");

   for(i=0;i<N;++i)

    {

       scanf("%lf%lf%lf%lf%lf%lf%lf%lf", &p1.x, &p1.y,&p2.x, &p2.y, &p3.x, &p3.y, &p4.x, &p4.y);

 

       if(direction(p3,p4,p1)==0 && direction(p3,p4,p2)==0)//共线

           printf("LINE\n");

       else

       {

           if( ((p1.x-p2.x)*(p3.y-p4.y)-(p1.y-p2.y)*(p3.x-p4.x))==0 )//平行

                printf("NONE\n");

           else

           {

                a1=p1.y-p2.y;b1=p2.x-p1.x;c1=p1.x*p2.y-p2.x*p1.y;

               a2=p3.y-p4.y;b2=p4.x-p3.x;c2=p3.x*p4.y-p4.x*p3.y;

                p0=solve(a1,b1,c1,a2,b2,c2);

                printf("POINT %.2f %.2f\n",p0.x,p0.y);

           }

       }

    }

   printf("END OF OUTPUT\n");

   return 0;

}

 

3. http://poj.org/problem?id=3304

 

题目大意：给出n条线段两个端点的坐标，问所有线段投影到一条直线上，如果这些所有投影至少相交于一点就输出Yes！，否则输出No！。
解题思路：如果有存在这样的直线，过投影相交区域作直线的垂线，该垂线必定与每条线段相交，问题转化为问是否存在一条线和所有线段相交

 

直线肯定经过两个端点。

枚举端点，判断直线和线段是否相交。

 

细节要注意，判断重合点。

 

还有就是加入只有一条线段的话，刚好直线是过同一条直线的。

所以保险的做法是枚举所有的两个端点，包括同一条直线的。

 

 

#include<cstdio>

#include<string>

#include<cstring>

#include<iostream>

#include<cmath>

#include<algorithm>

#include<vector>

using namespace std;

 

#define all(x) (x).begin(), (x).end()

#define for0(a, n) for (int (a) = 0; (a)< (n); (a)++)

#define for1(a, n) for (int (a) = 1; (a)<= (n); (a)++)

#define mes(a,x,s)  memset(a,x,(s)*sizeof a[0])

#define mem(a,x)  memset(a,x,sizeof a)

#define ysk(x)  (1<<(x))

typedef long long ll;

typedef pair<int, int> pii;

const int INF =0x3f3f3f3f;

const int maxn=  100  ;

int n;

struct Point

{

    double x,y;

   Point(double x=0,double y=0):x(x),y(y) {};

};

typedef Point Vector;

Vector operator +(Vector A,Vector B){return Vector(A.x+B.x,A.y+B.y); }

Vector operator -(Vector A,Vector B){return Vector(A.x-B.x,A.y-B.y); }

Vector operator *(Vector A,double p){return Vector(A.x*p,A.y*p); }

Vector operator /(Vector A,double p){return Vector(A.x/p,A.y/p); }

Vector operator -(Vector A)  {return Vector(-A.x,-A.y);}

 

const double eps=1e-8;

 

int dcmp(double x)

{

   if(fabs(x)<eps)  return 0;

   else return x<0?-1:1;

}

bool operator ==(Point A,Point B) {returndcmp(A.x-B.x)==0&&dcmp(A.y-B.y)==0;}

 

double Cross(Vector A,Vector B)//叉乘

{

   return A.x*B.y-A.y*B.x;

}

 

struct Line

{

   Point p,p2;

   Vector v;

   Line(){}

   Line(Point p,Vector v):p(p),v(v){}//点线式

   void twoPointIntial(Point p,Point p2)//两点式

    {

       this->p=p;

       this->p2=p2;

       v=  p2-p;

    }

 

}seg[maxn+5];

 

typedef Line Seg;

 

 

bool intersect(Line a,Seg b)//判断直线和线段是否相交

{

   Point A=b.p;

   Point B=b.p2;

   Vector v1=a.p-A;

   Vector v2=a.p2-A;

   Vector v3=a.p-B;

   Vector v4=a.p2-B;

   return dcmp(Cross(v1,v2))*dcmp(Cross(v3,v4))<=0; //如果线段的端点在直线上,肯定相交

}

bool solve(Point A,Point B)

{

    Line L;L.twoPointIntial(A,B);

    if( A==B )  return false;

    for(int i=0;i<n;i++)

    {

        if(!intersect(L,seg[i])    )   return false;

    }

    return true;

 

}

int main()

{

  std::ios::sync_with_stdio(false);

  int T;cin>>T;

  while(T--)

   {

     cin>>n;

     for0(i,n)

     {

         cin>>seg[i].p.x>>seg[i].p.y>>seg[i].p2.x>>seg[i].p2.y;

     }

     if(n==1) {puts("Yes!");continue;}

     bool find=false;

     for(int i=0;i<n&&!find;i++)

     {

         for(int j=i+1;j<n;j++)

         {

              find|=  solve(seg[i].p ,seg[j].p);

             find|=  solve(seg[i].p ,seg[j].p2);

              find|=  solve(seg[i].p2 ,seg[j].p);

              find|=  solve(seg[i].p2 ,seg[j].p2);

              if(find)  break;

         }

     }

 

     puts(find?"Yes!":"No!");

   }

  return 0;

}

 

3. http://poj.org/problem?id=2653

题目大意：有n根木条，一根一根的往一个坐标系上丢(给出木条两点的坐标)，问最后不被覆盖的木条有哪些，即丢的木条如果和前面丢的木条交叉的话，就会覆盖前面那根木条
解题思路：用到线段交叉的判断，我用队列模拟是否被覆盖，没被覆盖的入队，并且将当前火柴和前面没被覆盖的木条比较如果交叉则前面那根出队，这样的暴力做法在时间上不知道怎么优化了，不过比较了下，我的内存还是很低的

#include <iostream>

#include <queue>

#include <stdio.h>

using namespace std;

#define min(a,b) (a>b?b:a)

#define max(a,b) (a<b?b:a)

 

typedef struct

{

       intcnt;

       doublex1,y1,x2,y2;

}Line;

 

int n;

 

double direction(double x,double y,doublex1,double y1,double x2,double y2)

{     //叉积

       doublea1=x1-x;

       doubleb1=y1-y;

       doublea2=x2-x;

       doubleb2=y2-y;

       returna1*b2-a2*b1;

}

 

int on_segment(double x1,double y1,doublex2,double y2,double x,double y)

{            //判断共线

       if((min(x1,x2)<=x&& x<=max(x1,x2)) && (min(y1,y2)<=y &&y<=max(y1,y2)))

              return1;

       return0;

}

 

int cross(Line v,Line t)

{            //是否交叉，根据算法导论写的

       doubled1,d2,d3,d4;

       d1=direction(t.x1,t.y1,t.x2,t.y2,v.x1,v.y1);   //算叉积

       d2=direction(t.x1,t.y1,t.x2,t.y2,v.x2,v.y2);

       d3=direction(v.x1,v.y1,v.x2,v.y2,t.x1,t.y1);

       d4=direction(v.x1,v.y1,v.x2,v.y2,t.x2,t.y2);

       if(d1*d2<0&& d3*d4<0) return 1;             //直接和0比较的话时间是625MS

       if(!d1&& on_segment(t.x1,t.y1,t.x2,t.y2,v.x1,v.y1)) return 1;

       if(!d2&& on_segment(t.x1,t.y1,t.x2,t.y2,v.x2,v.y2)) return 1;

       if(!d3&& on_segment(v.x1,v.y1,v.x2,v.y2,t.x1,t.y1)) return 1;

       if(!d4&& on_segment(v.x1,v.y1,v.x2,v.y2,t.x2,t.y2)) return 1;

       return0;

}

 

void input()

{

       queue<Line>q;

       Linev,t;

       inti;

       scanf("%lf%lf%lf%lf",&v.x1,&v.y1,&v.x2,&v.y2);

       v.cnt=1;

       q.push(v);

 

       for(i=2;i<=n;i++)

  {

              scanf("%lf%lf%lf%lf",&t.x1,&t.y1,&t.x2,&t.y2);

              t.cnt=i;

              q.push(t);              //用当前的作为队尾

              while(!q.empty())

    {

                     v=q.front();q.pop();              //一个一个出队，直到队尾

                     if(t.cnt==v.cnt)

                     {

                            q.push(t);

                            break;

                     }

                     if(!cross(v,t))q.push(v);        //如果不交叉继续入队

              }

       }

 

       v=q.front();q.pop();

       printf("Topsticks: %d",v.cnt);

       while(!q.empty())

  {

              v=q.front();

   q.pop();

              printf(",%d",v.cnt);

       }

       printf(".\n");

 

}

 

int main()

{

       while(scanf("%d",&n)&& n)

  {

              input();

       }

       return0;

}

 

2.

http://acm.hdu.edu.cn/showproblem.php?pid=5721

 

问题描述

为了寻找失去的爱人Cupid，Psyche需要完成Venus的最后一项任务：前往冥界，收集一盒冥界皇后Prosperina的美貌。

 

冥界有nn个神殿，可以表示为平面上的nn个整点。Psyche想要找到这nn座神殿中，最近的两座神殿之间的距离。传说那就是通往主殿的密码。

 

但是冥界神秘莫测，在不同的时刻，这nn座神殿中的某一座会消失。

 

Psyche想要知道，对于nn座神殿中的任意一座消失的情况，最近的两座神殿之间的距离。你只需要输出它们的和。

 

为避免精度误差，定义两点(x_1, y_1), (x_2, y_2)(x​1​​,y​1​​),(x​2​​,y​2​​)间的距离为d = (x_1 - x_2) ^ 2 + (y_1 - y_2) ^ 2d=(x​1​​−x​2​​)​2​​+(y​1​​−y​2​​)​2​​。

输入描述

第一行，一个整数TT(1 \le T \le 5)(1≤T≤5)，代表数据组数。

 

对于每组数据，第一行，一个整数nn(3 \le n \le 10 ^ 5)(3≤n≤10​5​​)，代表神殿个数。

 

下面nn行，每行两个整数x, yx,y(-10 ^ 5 \le x,y \le 10 ^ 5)(−10​5​​≤x,y≤10​5​​)，代表神殿的位置在(x, y)(x,y)。

 

注意可能存在两座神殿坐落在同一位置。

输出描述

输出TT行，对于每组数据，输出nn座神殿中的任意一座消失的情况，最近两座神殿之间的距离的和。

输入样例

1

3

0 0

1 1

2 2

输出样例

12

Hint

神殿(0,0)(0,0)消失时，d = (1-2) ^ 2 + (1 - 2) ^ 2 = 2d=(1−2)​2​​+(1−2)​2​​=2；

 

神殿(1,1)(1,1)消失时，d = (0-2) ^ 2 + (0 - 2) ^ 2 = 8d=(0−2)​2​​+(0−2)​2​​=8；

 

神殿(2,2)(2,2)消失时，d = (0-1) ^ 2 + (0-1) ^ 2 = 2d=(0−1)​2​​+(0−1)​2​​=2；

 

故答案为2 + 8 + 2 = 122+8+2=12。

 

解题思路：

感觉就是很裸的平面最近点对的模板题= =题解上的说的也很明确了，只有三种情况

一种情况是删的点不属于最近点对里面的点

剩下两种情况各删的最近点对里面的其中一点

 

 

#include <cstdio>

#include <cstring>

#include <algorithm>

using namespace std;

 

#define INF 1e18

#define Min(a, b) ((a) <= (b) ? (a) :(b))

 

typedef long long LL;

typedef struct node{

   LL x, y;

   LL index;

}Coord;

 

const int maxn = 1e5 + 5;

 

LL p1, p2, ans;

Coord p[maxn], ppp[maxn], tmp[maxn];

 

inline LL f(LL x){

   return (x >= 0 ? x : -x);

}

inline bool cmp(Coord a, Coord b){

   return a.x < b.x;

}

inline bool cmp2(Coord a, Coord b){

   return a.y < b.y;

}

inline LL Dist(Coord a, Coord b){

   return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);

}

LL Cal_Closest(LL left, LL right, intstep){

   LL d = INF;

   if(left == right){

       return d;

    }

   if(left + 1 == right){

       LL tmp2 = Dist(p[left], p[right]);

       if(step == 0){

           if(step == 0 && tmp2< ans){

                p1 = left;

                p2 = right;

                ans = tmp2;

           }

       }

       return tmp2;

    }

   

   LL mid = (left + right) >> 1;

   LL d1 = Cal_Closest(left, mid, step);

   LL d2 = Cal_Closest(mid + 1, right, step);

    d= Min(d1, d2);

   

   LL k = 0;

   for(LL i = left; i <= right; ++i){

       if(f(p[mid].x - p[i].x) <= d){

           tmp[k++] = p[i];

       }

    }

   sort(tmp, tmp + k, cmp2);

   

   for(LL i = 0; i < k; ++i){

       for(LL j = i + 1; j < k && tmp[j].y - tmp[i].y < d; ++j){

           if(d > Dist(tmp[i], tmp[j])){

                d = Dist(tmp[i], tmp[j]);

                if(step == 0 && d <ans){

                    p1 = tmp[i].index;

                   p2 = tmp[j].index;

                    ans = d;

                }

           }

       }

    }

   return d;

}

 

int main()

{

   LL t, n;

   scanf("%lld", &t);

   while(t--){

       scanf("%lld", &n);

       memset(p, 0, sizeof(p));

       memset(ppp, 0, sizeof(ppp));

       for(LL i = 0; i < n; ++i){

           scanf("%lld%lld", &p[i].x, &p[i].y);

       }

       

       sort(p, p + n, cmp);

       for(LL i = 0; i < n; ++i) p[i].index = i;

       for(LL i = 0; i < n; ++i) ppp[i] = p[i];

       p1 = 0; p2 = 0;

       ans = INF;

       

       ans = Cal_Closest(0, n - 1, 0);

       ans = ans * (n - 2);

       

       LL t = 0;

       for(LL i = 0; i < n; ++i){

           if(ppp[i].index == p1) continue;

           p[t++] = ppp[i];

       }

       sort(p, p + t, cmp);

       ans += Cal_Closest(0, t - 1, 1);

       

       t = 0;

       for(LL i = 0; i < n; ++i){

           if(ppp[i].index == p2) continue;

           p[t++] = ppp[i];

       }

       sort(p, p + t, cmp);

        ans += Cal_Closest(0, t - 1, 1);

       printf("%lld\n", ans);

    }

   return 0;

}

计算几何 Part.2---凸包问题

关于凸包的严格定义，这里不打算写出来，大家可以自行Google或者百度，因为严格的数学定义反而不太好理解，用最通俗的话来解释凸包：给定二维平面上的点集，凸包就是将最外层的点连接起来构成的凸多边型，它能包含点集中所有的点。

 

1.       

acm.hdu.edu.cn/showproblem.php?pid=3847

 

/*

求出多边形最窄的地段长度

枚举边，求出所有点中到边的距离最大的值

这些值中最小的就是答案

*/

#include <cstdio>

#include <cstring>

#include <cmath>

#include <cstdlib>

const int MAXN = 109;

const double eps = 1e-4;

struct point{

    double x,y;

}p[MAXN],h[MAXN];

inline double distance(const point &p1,const point &p2){

    return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));

}//   求两点之间的距离

inline double multiply(const point &sp,const point &ep,const point &op){

      return ((sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y));

}//判断sp,ep，op是否满足左转

int cmp(const void *a,const void *b){//按极角排序

    point *p1 = (point *)a;

    point *p2 = (point *)b;

    double t = multiply(*p2,*p1,p[0]);

    if(t>eps) return 1;

    else if(fabs(t)<=eps) 

    {

    if(distance(p[0],*p1)>distance(p[0],*p2)) return 1;

    else return -1;

    }

    else return -1;

}

void anglesort(point p[],int n){//找到最左下方的点

    int i,k=0;

    point temp;

    for(i=1;i<n;i++)

        if(p[i].x<p[k].x ||( fabs(p[i].x-p[k].x)<eps && (p[i].y<p[k].y)))

            k=i;

    temp=p[0],p[0]=p[k],p[k]=temp;

    qsort(p+1,n-1,sizeof(point),cmp);

}

void Graham_scan(point p[],point ch[],int n,int &len){//建立凸包

    int i,top=2;

    anglesort(p,n);

    if(n<3){

        for(i=0,len=n;i<n;i++) ch[i]=p[i];

        return;

    }

    ch[0]=p[0],ch[1]=p[1],ch[2]=p[2];

    p[n]=p[0];

    for(i=3;i<n;i++){

        while(multiply(p[i],ch[top],ch[top-1])>=eps) top--;

        ch[++top]=p[i];

    }

    len=top+1;

}

double judge(point _x,point _y,int len)

{

    double res=0;

    for(int i=0;i<len;i++)

    {

        double tmp=fabs(multiply(h[i],_x,_y))/distance(_x,_y);//面积法求出距离

        if(tmp>res)

        res=tmp;

    }

    return res;

}

double solve(int len)

{

    h[len]=h[0];

    double ans=1<<30;

    for(int i=0;i<len;i++)

    {

        double tmp=judge(h[i],h[i+1],len);//枚举边

        if(tmp<ans)

        ans=tmp;

    }

    return ans;

}

int main(){

    int n,len;

    int d=0;

    while(scanf("%d",&n),n)

    {

        for(int i=0;i<n;i++) scanf("%lf %lf",&p[i].x,&p[i].y);

        Graham_scan(p,h,n,len);

        printf("Case %d: %.2lf\n",++d,solve(len)+0.005);

    }

    return 0;

}

 

 

 

 

 

2.

 

http://poj.org/problem?id=2187

 

大致题意：

给定平面上的一些散点集，求最远两点距离的平方值。

 

解题思路：

别想着暴力枚举任意亮点距离找最大，行不通，想想三点共线吧！

 

平面上的散点集的最远的两点距离必然在这个散点集的凸包的某两个顶点上出现。

那么先求凸包，再枚举顶点距离就OK了。

 

别看是3000ms就想用简单的卷包裹，这题数据规模极大，卷包裹铁超（我一开始就是这么做的。。。） 万般无奈不得不用GrahamScanAlgorithm。。。。O(nlogn)用来做这题还是相当可观的。

 

GrahamScan理解是不困难的，同学们百度搜搜就知道基本思想了：数据结构：栈。而且每个点最多入栈出栈一次。

 

问题是用栈构造凸包之前要对散点集进行排序，水平序我不说了，我没看懂，听说很简单。我用的是极坐标排序。

利用叉积的旋向 配合 比较排序 （如插入排序，冒泡，快排等）可以对极坐标排序，推荐用qsort，不要用冒泡之类的，太慢了，用Graham都是想快而已。

 

qsort（）难点在于 比较规则，（我的程序为cmp函数），必须把“qsort对cmp返回值的处理、cmp返回值、叉积旋向返回值”三者有机结合，否则一定排序出错，详见我的程序，有详细解释。

Cmp比较规则为“按极角大小逆时针排序，极角相等则近极坐标的点优先”。在网上看到有些童鞋在“极角相同”的情况下，利用第二关键字“距离”计算排序时，是通过计算两点的横坐标到 原点横坐标之差 作为判断“距离”依据。乍一看好像没错，也能AC，其实那是POJ数据库不强大，试试多点与原点的连线垂直于x轴的情况吧！

 

最后注意的：在使用qsort之前，必须先找到 散点集中 最左下角的点，把它放在散点集数组的最后一个元素位，作为 极坐标原点（快排的参考点），否则排序也会出错。

 

其实只要qsort的cmp函数能处理好了，基本这题就过了，难度不大。

 

 

#include<iostream>

#include<algorithm>

using namespace std;

 

const int inf=10001;

 

typedef class location

{

   public:

       int x,y;

}node;

 

/*点A到点B的距离的平方*/

 

int distsquare(node A,node B)

{

   return (B.x-A.x)*(B.x-A.x)+(B.y-A.y)*(B.y-A.y);

}

 

/*计算叉积*/

 

int det(int x1,int y1,int x2,int y2)

{

   return x1*y2-x2*y1;

}

 

int cross(node A,node B,node C,node D)

{

   return det(B.x-A.x , B.y-A.y , D.x-C.x , D.y-C.y);

}

 

/*qsort比较规则*/

 

node* p; //极坐标原点

int cmp(const void* a,const void* b)

{

   node* s=(node*)a;

   node* t=(node*)b;

 

   int temp=cross(*p,*s,*p,*t);  //叉乘ps X pt

   if(temp>0)    //说明pt向量的极角小于 ps向量的极角

       return -1;  //从逆时针排序角度，t点位置在s点前面，即t<s ，根据qsort规则返回-1

 

   else if(temp==0) //说明pt向量的极角 等于 ps向量的极角

       return distsquare(*p,*t)-distsquare(*p,*s);  //距离原点近的点优先排序，用减法能够实现3出口：- 0 +

                        //注意，网上有些程序在这里不是比较 距离之差，而是比较横坐标绝对值 之差

                        //这是欠缺考虑  多点与原点连线垂直于x轴，不完善，之所以能AC是因为POJ的数据库不足够大而已

   else

       return 1; //pt向量的极角 大于 ps向量的极角

}

 

int main()

{

   int n,i,j;

   while(cin>>n)

    {

       node* farm=new node[n+1];

       int* conbag=new int[n+2]; //conbag[]顺序记录输入的点中构造成凸包的顶点集的各点在farm[]中的下标

 

       /*Input & search the first vertex*/

 

       int min_x=inf;

       int fi=0;

       for(i=1;i<=n;i++)

       {

           cin>>farm[i].x>>farm[i].y;

 

           if(min_x > farm[i].x)

           {

                min_x = farm[i].x;   //先以横坐标为第一关键字搜索最左下角的点

                fi=i;

           }

           else if(min_x == farm[i].x)

           {

                if(farm[fi].y >farm[i].y)   //若第一关键字相同，则以纵坐标作为第二关键搜索

                    fi=i;

           }

       }

 

       /*Quicksort Point Sets*/

 

       farm[0]=farm[n];   //这三步非常关键，是使用qsort前的准备工作

       farm[n]=farm[fi];  //目的是把前面找到的最左下角的点作为 极坐标原点

       farm[fi]=farm[0];  //即把第fi个点移到farm[]最后的位置,qsort则会把它作为排序的参考点

 

       p=&farm[n]; //极坐标原点传参

 

       qsort(farm+1,n,sizeof(node),cmp); //farm[]散点集排序

 

       /*Graham Scan Algorithm*/

 

       int pc=1;  //conbag[]指针

       conbag[1]=n;  //(约定)极坐标原点为凸包第1个顶点

       conbag[++pc]=1; //(在前面基础上,)有序点集farm[]的第一个点 (必)为凸包第2个顶点

       conbag[0]=2;   //凸包顶点计数器

 

       for(i=2;i<=n;)

           if(cross(farm[ conbag[pc-1] ],farm[ conbag[pc] ],farm[ conbag[pc]],farm[i]) >= 0)

           {   //检查向量(前一点pc-1,当前点pc) 与 向量(当前点pc,待入栈点i) 的旋转关系

                conbag[++pc]=i++;  //入栈

                conbag[0]++;

           }

           else

           {

                pc--;         //当前点pc出栈

                conbag[0]--;

           }

 

        /*Search The Max distant*/

 

       int maxdist=0;

       for(i=1;i<=conbag[0]-2;i++)         //散点集的两点最大距离必定出现在该散点集的凸包的某两个顶点之间

           for(j=i+1;j<=conbag[0]-1;j++)

           {

                int temp=distsquare(farm[conbag[i] ],farm[ conbag[j] ]);

                if(maxdist < temp)

                    maxdist = temp;

           }

 

       /*Output*/

 

       cout<<maxdist<<endl;

 

 

       delete farm;

       delete conbag;

    }

   return 0;

}

 

 

计算几何 Part.3---面积公式矩形切割

1.      http://acm.hdu.edu.cn/showproblem.php?pid=1700

题意：以原点为圆心，给出圆上的一点，要求两位两点，是的这三个点的距离和最大，很容易想到这是一个等边三角形，而事实上，经过对题目给出样例的测试也证明了这确实是一个等边三角形

思路：几何水题

我们可以得到方程组

x^2+y^2 = r^2

(a-x)^2+(b-y)^2=3r^2

解方程组得到的两点即为三角形的另外两点

#include <stdio.h>

#include <math.h>

 

int main()

{

   int t;

   double x,y,x2,y2,r;

   double ax,ay,bx,by,k,m,l,A,B,C;

   scanf("%d",&t);

   while(t--)

    {

       scanf("%lf%lf",&x,&y);

       r = x*x+y*y;

       A = r;

       B = y*r;

        C = r*r/4-r*x*x;

       ay = (-B-sqrt(B*B-4*A*C))/(2*A);

       by = (-B+sqrt(B*B-4*A*C))/(2*A);

       if(fabs(x-0)<1e-7)//防止除数出现0的情况

       {

           ax=-sqrt(r-ay*ay);

           bx=sqrt(r-by*by);

       }

       else

       {

           ax=(-r/2-ay*y)/x;//由于ay必定小于by，所以ax也必定小于bx，所以无需进行大小判定

           bx=(-r/2-by*y)/x;

       }

       printf("%.3lf %.3lf %.3lf %.3lf\n",ax,ay,bx,by);

    }

   return 0;

}

 

 

计算几何 Part.4---半平面交

1.

题目链接： http://poj.org/problem?id=1755

 

题目：铁人三项，每个人在某一项中有确定的速度，裁判可以决定某一项比赛的路程为多少，问对于某个人，是否存在一种安排能使他拿到第一，而且不能是并列。

我们假设三项的路程分别人X,Y,Z。

比较其中的两个人。A的时间为X / U1+Y / V1+Z /W1     B的时间为X / U2 +Y / V2 +Z / W2

如果A想要获胜妈，X / U1+Y / V1+Z / W1- X / U2 +Y / V2 +Z / W2 < 0

由于我写的是顺时针的，所以把不等式变个号。这里一定要注意，小于0大于0和顺时针逆时针的对应关系

这个不等式还是有3个未知数。显然用三维就太麻烦了。由于我们最终不需要求出X,Y,Z具体为多少，而且Z>0，所以把不等式两边同时除以Z，则把X/Z看成一个未知量，Y/Z看成另外一个。

问题转化成一系列的不等式是否为解，而且注意条件X>0Y>0

那么可以设立一个初始范围，(0,0)(0,inf)(inf,inf)(inf,0)

然后通过两个人的参数，求出AX+BY+C>0，通过半平面交解决，最终判断面积是否为0

有几个地方需要注意：

这题要求的精度很高，大家一味的说需要1e-16，其实不然，1e-8也过了。主要是中间的处理细节，对于1/U1-1/U2，普通的处理是需要两次除法，精度严重受损，可以改成(U2-U1)/(U1*U2)。

另外就是特判，题目要求是不能并列，所以最终结果是大于0才行。而且如果遇到A==0&&B==0&&C<=0说明不等式无解，直接返回

#include <stdio.h>

#include <string.h>

#include <iostream>

#include <algorithm>

#include <vector>

#include <queue>

#include <set>

#include <map>

#include <string>

#include <math.h>

#include <stdlib.h>

#include <time.h>

using namespace std;

const double eps = 1e-18;

int sgn(double x)

{

    if(fabs(x) < eps)return 0;

    if(x < 0)return -1;

    elsereturn 1;

}

struct Point

{

    double x,y;

    Point(){}

    Point(double _x,double _y)

    {

        x = _x; y = _y;

    }

    Point operator -(const Point &b)const

    {

        return Point(x - b.x, y - b.y);

    }

    double operator ^(const Point &b)const

    {

        return x*b.y - y*b.x;

    }

    double operator *(const Point &b)const

    {

        return x*b.x + y*b.y;

    }

};

//计算多边形面积

double CalcArea(Point p[],intn)

{

    double res = 0;

    for(int i = 0;i < n;i++)

        res += (p[i]^p[(i+1)%n]);

    return fabs(res/2);

}

//通过两点，确定直线方程

void Get_equation(Pointp1,Point p2,double &a,double &b,double &c)

{

    a = p2.y - p1.y;

    b = p1.x - p2.x;

    c = p2.x*p1.y - p1.x*p2.y;

}

//求交点

Point Intersection(Pointp1,Point p2,double a,double b,double c)

{

    double u = fabs(a*p1.x + b*p1.y + c);

    double v = fabs(a*p2.x + b*p2.y + c);

    Point t;

    t.x = (p1.x*v + p2.x*u)/(u+v);

    t.y = (p1.y*v + p2.y*u)/(u+v);

    return t;

}

Point tp[110];

void Cut(double a,doubleb,double c,Point p[],int &cnt)

{

    int tmp = 0;

    for(int i = 1;i <= cnt;i++)

    {

        //当前点在左侧，逆时针的点

        if(a*p[i].x + b*p[i].y + c <eps)tp[++tmp] = p[i];

        else

        {

            if(a*p[i-1].x + b*p[i-1].y + c <-eps)

                tp[++tmp] =Intersection(p[i-1],p[i],a,b,c);

            if(a*p[i+1].x + b*p[i+1].y + c <-eps)

                tp[++tmp] =Intersection(p[i],p[i+1],a,b,c);

        }

    }

    for(int i = 1;i <= tmp;i++)

        p[i] = tp[i];

    p[0] = p[tmp];

    p[tmp+1] = p[1];

    cnt = tmp;

}

double V[110],U[110],W[110];

int n;

const double INF = 100000000000.0;

Point p[110];

bool solve(int id)

{

    p[1] = Point(0,0);

    p[2] = Point(INF,0);

    p[3] = Point(INF,INF);

    p[4] = Point(0,INF);

    p[0] = p[4];

    p[5] = p[1];

    int cnt = 4;

    for(int i = 0;i < n;i++)

        if(i != id)

        {

            double a = (V[i] -V[id])/(V[i]*V[id]);

            double b = (U[i] -U[id])/(U[i]*U[id]);

            double c = (W[i] -W[id])/(W[i]*W[id]);

            if(sgn(a) == 0 && sgn(b) ==0)

            {

                if(sgn(c) >= 0)return false;

                else continue;

            }

            Cut(a,b,c,p,cnt);

        }

    if(sgn(CalcArea(p,cnt)) == 0)return false;

    else return true;

}

int main()

{

   //freopen("in.txt","r",stdin);

   //freopen("out.txt","w",stdout);

    while(scanf("%d",&n) == 1)

    {

        for(int i = 0;i < n;i++)

           scanf("%lf%lf%lf",&V[i],&U[i],&W[i]);

        for(int i = 0;i < n;i++)

        {

           if(solve(i))printf("Yes\n");

            else printf("No\n");

        }

    }

    return 0;

}

 

 

2.

http://poj.org/problem?id=1279

对于多边形核不懂的，自行查资料了。

第一道半平面交，只会写N^2。

这里默认是顺时针的顺序，否则就要调整一下。

将每条边化作一个不等式，ax+by+c>0，所以要固定顺序，方便求解。

半平面交其实就是对一系列的不等式组进行求解可行解。

如果某点在直线右侧，说明那个点在区域内，否则出现在左边，就可能会有交点，将交点求出加入。

#include<iostream>

#include<fstream>

#include<iomanip>

#include<cstdio>

#include<cstring>

#include<algorithm>

#include<cstdlib>

#include<cmath>

#include<set>

#include<map>

#include<queue>

#include<stack>

#include<string>

#include<vector>

#include<sstream>

#include<cassert>

#define LL long long

#define eps 1e-7

#define inf 1<<30

using namespace std;

struct Point{

    double x,y;

}p[1505],tp[1505],q[1505];

//叉积

double xmul(Point p0,Point p1,Point p2){

    return(p1.x-p0.x)*(p2.y-p0.y)-(p1.y-p0.y)*(p2.x-p0.x);

}

//通过两点，确定直线方程

double Get_equation(Point p1,Point p2,double &a,double&b,double &c){

    a=p2.y-p1.y;

    b=p1.x-p2.x;

    c=p2.x*p1.y-p1.x*p2.y;

}

//求交点

Point Intersection(Point p1,Point p2,double a,double b,double c){

    doubleu=fabs(a*p1.x+b*p1.y+c);

    doublev=fabs(a*p2.x+b*p2.y+c);

    Point t;

   t.x=(p1.x*v+p2.x*u)/(u+v);t.y=(p1.y*v+p2.y*u)/(u+v);

    return t;

}

//求面积，正为顺时针，和叉积写法有关

double Get_area(Point p[],int n){

    double area=0;

    for(int i=2;i<n;i++)

       area+=xmul(p[1],p[i],p[i+1]);

    return -area/2.0;

}

//改变顺序

double Change_dir(Point p[],int n){

    for(int i=1;i<=n/2;i++)

        swap(p[i],p[n+1-i]);

}

//加入一条边，切割

void Cut(double a,double b,double c,Point p[],int &cnt){

    int tmp=0;

    for(inti=1;i<=cnt;i++){

        //当前点就在右侧

       if(a*p[i].x+b*p[i].y+c>-eps) tp[++tmp]=p[i];

        else{

            //前一个点在右侧，产生交点

           if(a*p[i-1].x+b*p[i-1].y+c>eps)

                tp[++tmp]=Intersection(p[i-1],p[i],a,b,c);

            //同理

           if(a*p[i+1].x+b*p[i+1].y+c>eps)

               tp[++tmp]=Intersection(p[i],p[i+1],a,b,c);

        }

    }

    for(int i=1;i<=tmp;i++)

        p[i]=tp[i];

    p[0]=p[tmp];p[tmp+1]=p[1];

    cnt=tmp;

}

void slove(Point p[],int n){

    //默认顺时针，通过面积判断一下

    if(Get_area(p,n)<eps)Change_dir(p,n);

    p[0]=p[n];p[n+1]=p[1];

    //原来的点要备份一遍，查了好久

    for(int i=0;i<=n+1;i++)q[i]=p[i];

    int cnt=n;

    for(int i=1;i<=n;i++){

        double a,b,c;

       Get_equation(q[i],q[i+1],a,b,c);

        Cut(a,b,c,p,cnt);

    }

    printf("%.2f\n",fabs(Get_area(p,cnt)));

}

int main(){

    int t,n;

   scanf("%d",&t);

    while(t--){

       scanf("%d",&n);

        for(inti=1;i<=n;i++)

            scanf("%lf%lf",&p[i].x,&p[i].y);

        slove(p,n);

    }

    return 0;

}

 

 

 

3.

题目链接:http://poj.org/problem?id=2451

题意：在（0，10000）*（0，10000）的坐标系上，给定n个半平面，求出它们围成的图形的面积

每个半平面由两点（x1，y1）（x2，y2）确定的直线确定，规定半平面为直线的左边，即存在一点（x，y） 

使得(x – x1) *(y – y2) – (x – x2) * (y – y1) = (x1 – x) * (y2 – y) – (x2 – x) * (y1 – y)>=0成立

//模板题直接copy网上模板

#include <iostream>

#include <cstdio>

#include <cmath>

#include <algorithm>

using namespace std;

 

const double eps = 1e-8;

int n, pn, dq[20005], top, bot;

struct Point {

    double x, y;

} p[20005];

 

struct Line {

    Point a, b;

    double angle;

    Line& operator= (Linel) {

       a.x = l.a.x; a.y =l.a.y;

       b.x = l.b.x; b.y =l.b.y;

       angle = l.angle;

       return *this;

    }

} l[20005];

 

int dblcmp(double k) {

    if (fabs(k) < eps)return 0;

    return k > 0 ? 1 : -1;

}

 

double multi(Point p0, Point p1, Point p2) {

    return(p1.x-p0.x)*(p2.y-p0.y)-(p1.y-p0.y)*(p2.x-p0.x);

}

 

bool cmp(const Line& l1, const Line& l2) {

    int d =dblcmp(l1.angle-l2.angle);

    if (!d) returndblcmp(multi(l1.a, l2.a, l2.b)) > 0; //极角相同时，将更靠半平面里面的放在前面

    return d < 0;

}

 

void addLine(Line& l, double x1, double y1, double x2, doubley2) {

    l.a.x = x1; l.a.y = y1;

    l.b.x = x2; l.b.y = y2;

   /*atan2(y,x)返回值的取值范围为-PI到PI，实际上就是根据向量（x,y）确定极角，

      若向量在1，2象限，则值大于0；3，4象限，则值小于0

     所以atan2（1，1）与atan（-1，-1）值不同

    */

    l.angle = atan2(y2-y1,x2-x1);

}

 

//求交点，写得略复杂

 void getIntersect(Line l1,Line l2, Point& p) {

    double A1 = l1.b.y -l1.a.y;

    double B1 = l1.a.x -l1.b.x;

    double C1 = (l1.b.x -l1.a.x) * l1.a.y - (l1.b.y - l1.a.y) * l1.a.x;

    double A2 = l2.b.y -l2.a.y;

    double B2 = l2.a.x -l2.b.x;

    double C2 = (l2.b.x -l2.a.x) * l2.a.y - (l2.b.y - l2.a.y) * l2.a.x;

    p.x = (C2 * B1 - C1 * B2)/ (A1 * B2 - A2 * B1);

    p.y = (C1 * A2 - C2 * A1) /(A1 * B2 - A2 * B1);

}

 

bool judge(Line l0, Line l1, Line l2) {

     Point p;

     getIntersect(l1, l2, p);

     return dblcmp(multi(p,l0.a, l0.b)) < 0;

}

 

void HalfPlaneIntersect( ){

    int i, j;

   /*排序是在满足所有半平面A*x+B*y+C>0或（<,<=,>=）,

      也就是所有半平面的符号均相同的情况下对极角进行排序，

      此题中题目对半平面的定义就等价于这个条件

 */

    sort(l, l+n, cmp);

    for (i = 0, j = 0; i <n; i++)

        if(dblcmp(l[i].angle-l[j].angle) > 0) //极角相同时，只保留最靠里面的那条

            l[++j] = l[i];

    n = j + 1;

    dq[0] = 0; //双端队列

    dq[1] = 1;

    top = 1;   //顶部和底部

    bot = 0;

    for (i = 2; i < n; i++){

       //当栈顶的两条直线交点在当前半平面外部时，弹栈

        while (top > bot&& judge(l[i], l[dq[top]], l[dq[top-1]])) top--;

       /*由于求的是一个凸多边形，所以当半平面转过接近一圈时，某个半平面满足上一个while的条件后，

          它又会影响到底部的两条直线，当底部的两条直线的交点，在当前的半平面外部时，底部弹栈

        */

        while (top > bot&& judge(l[i], l[dq[bot]], l[dq[bot+1]])) bot++;

        dq[++top] = i;//当前半平面入栈

    }

   //当最顶部的两条直线的交点不在最底部的半平面内时，顶部的那个半平面是多余的，顶部弹栈

    while (top > bot&& judge(l[dq[bot]], l[dq[top]], l[dq[top-1]])) top--;

   //当最底部的两条直线的交点不在最顶部的半平面内时，底部的那个半平面是多余的，底部弹栈

    while (top > bot&& judge(l[dq[top]], l[dq[bot]], l[dq[bot+1]])) bot++;

    dq[++top] = dq[bot];//将最底部的半平面放到最顶部来，方便下面求顶点

    for (pn = 0, i = bot; i< top; i++, pn++)

        getIntersect(l[dq[i+1]],l[dq[i]], p[pn]);

}

 

double getArea() { //叉积求面积

    if (pn < 3) return 0;

    double area = 0;

    for (int i = 1; i <pn-1; i++)

        area += multi(p[0],p[i], p[i+1]);

    if (area < 0) area =-area;

    return area/2;

}

 

int main()

{

    int i;

    double x1, y1, x2, y2;

 

    while (scanf("%d", &n) != EOF) {

       //输入半平面，由一条线段确定，半平面为在线段的哪边题目已经给出

        for (i = 0; i < n;i++) {

             scanf("%lf%lf%lf%lf", &x1, &y1, &x2, &y2);

             addLine(l[i], x1,y1, x2, y2);

        }

       //这四个半平面的添加要满足题目给出的半平面的定义

        addLine(l[n++], 0, 0,10000, 0);

        addLine(l[n++], 10000,0, 10000, 10000);

        addLine(l[n++], 10000,10000, 0, 10000);

        addLine(l[n++], 0,10000, 0, 0);

        HalfPlaneIntersect();

        printf("%.1lf\n", getArea());

    }

    return 0;

}

计算几何 Part.5---计算几何背景，实际上解题的关键是其他问题（数据结构、组合数学，或者是枚举思想）若干道经典的离散化＋扫描线的题目

1.

http://poj.org/problem?id=1151

 

 

离散化: 将所有的x轴坐标存在一个数组里..排序.当进入一条线段时..通过二分的方式确定其左右点对应的离散值...

扫描线..可以看成一根平行于x轴的直线..至y=0开始往上扫..直到扫出最后一条平行于x轴的边..但是真正在做的时候..不需要完全模拟这个过程..扫描线的做法是从最下面的边开始扫到最上面的边.

线段树: 本题用于动态维护扫描线在往上走时..x哪些区域是有合法面积的..

 

/*

POJ 1151 Atlantis

求矩形并的面积（线段树+离散化）

*/

#include<stdio.h>

#include<iostream>

#include<algorithm>

using namespace std;

#define MAXN 201

struct Node

{

   int l,r;//线段树的左右整点

   int c;//c用来记录重叠情况

   double cnt,lf,rf;//

   //cnt用来计算实在的长度，rf,lf分别是对应的左右真实的浮点数端点

}segTree[MAXN*3];    

struct Line

{

   double x,y1,y2;

   int f;

}line[MAXN];

//把一段段平行于y轴的线段表示成数组 ，

//x是线段的x坐标，y1,y2线段对应的下端点和上端点的坐标

//一个矩形，左边的那条边f为1，右边的为-1，

//用来记录重叠情况，可以根据这个来计算，nod节点中的c

 

bool cmp(Line a,Line b)//sort排序的函数

{

   return a.x < b.x;

}    

 

double y[MAXN];//记录y坐标的数组

void Build(int t,int l,int r)//构造线段树

{

   segTree[t].l=l;segTree[t].r=r;

   segTree[t].cnt=segTree[t].c=0;

   segTree[t].lf=y[l];

   segTree[t].rf=y[r];

   if(l+1==r)  return;

   int mid=(l+r)>>1;

   Build(t<<1,l,mid);

   Build(t<<1|1,mid,r);//递归构造

}    

void calen(int t)//计算长度

{

   if(segTree[t].c>0)

    {

       segTree[t].cnt=segTree[t].rf-segTree[t].lf;

       return;

   }   

   if(segTree[t].l+1==segTree[t].r) segTree[t].cnt=0;

   else segTree[t].cnt=segTree[t<<1].cnt+segTree[t<<1|1].cnt;

}    

void update(int t,Line e)//加入线段e，后更新线段树

{

   if(e.y1==segTree[t].lf&&e.y2==segTree[t].rf)

    {

       segTree[t].c+=e.f;

       calen(t);

       return;

   }   

   if(e.y2<=segTree[t<<1].rf) update(t<<1,e);

   else if(e.y1>=segTree[t<<1|1].lf)  update(t<<1|1,e);

   else

    {

       Line tmp=e;

       tmp.y2=segTree[t<<1].rf;

       update(t<<1,tmp);

       tmp=e;

       tmp.y1=segTree[t<<1|1].lf;

       update(t<<1|1,tmp);

   }   

   calen(t);

}   

int main()

{

   int i,n,t,iCase=0;

   double x1,y1,x2,y2;

   while(scanf("%d",&n),n)

    {

       iCase++;

       t=1;

       for(i=1;i<=n;i++)

       {

           scanf("%lf%lf%lf%lf",&x1,&y1,&x2,&y2);

            line[t].x=x1;

           line[t].y1=y1;

           line[t].y2=y2;

           line[t].f=1;

           y[t]=y1;

           t++;

           line[t].x=x2;

           line[t].y1=y1;

           line[t].y2=y2;

           line[t].f=-1;

           y[t]=y2;

           t++;

       }

       sort(line+1,line+t,cmp);

       sort(y+1,y+t);  

       Build(1,1,t-1);

       update(1,line[1]);

       double res=0;

       for(i=2;i<t;i++)

       {

           res+=segTree[1].cnt*(line[i].x-line[i-1].x);

           update(1,line[i]);

       }   

       printf("Test case #%d\nTotal explored area: %.2f\n\n",iCase,res);

       //看来POJ上%.2f可以过，%.2lf却不行了

   }   

   return 0;

}

 

计算几何 Part.6---解析几何

1.

http://acm.sgu.ru/problem.php?contest=0&problem=110

题目大意
空间探测器在星球M上发现了巨大的地牢，地牢被明亮的球充满，探测器发现光线能按自然规律被球表面反射（入射角等于反射角，入射光线、反射光线、法线在同一平面）。古老的传说说如果光按一定顺序被球表面反射，房间的门就会打开。
你不需要去猜这个顺序；你的任务更简单一些。你会知道球的位置和半径、激光发射的位置及光传播的方向。你要找出光被球反射的顺序。
输入

n(1<=n<=50)球的数量 下面n行读入球坐标，半径xi, yi, zi, ri(integer范围内)。
最后一行包含6个数。
前三个是激光发射的坐标(发射点严格的在任何球体外)。
后三个告诉你激光发射的方向(这个坐标在光线上)。
输出

光被球反射的顺序(球在输入中从1开始依次编号)。
如果球被反射10次以上，输出前十次，然后用一个空格和'etc.'隔开（不含引号）。
注意:如果光的轨迹是某个球的切线,认为光被球反射。

对于点(xs,ys,zs)和光线向量(xl,yl,zl)，在这条光线上的点(x,y,z)满足：
x=xs+k*xl;
y=ys+k*yl;
z=zs+k*zl;
其中，k∈R

对于球(x0,y0,z0)和半径r，在这个球上的点(x,y,z)满足：
(x0-x)^2+(y0-y)^2+(z0-z)^2=r^2;

联立上述四个式子，可以得到一个关于k的一元二次方程。
若该方程无解，则证明两者不相交。

在得到交点之后，问题就是反射了(请先行了解向量加减法的几何意义和光的反射定律)
法线向量v1为从球心指向交点的向量。
入射的反向量v2为从交点指向出发点的向量。
我们可以利用公式求出v1在v2上的投影向量v3
然后将(v3-v1)*2+v1得到反射向量(无法理解的话可以在二维平面上模拟画图帮助理解)。
反射出发点变为交点。

注意事项：1.一元二次方程的解可能有两个，我们要取的是离出发点近的那个交点。
          2.光线是一条射线而不是直线，所以一元二次方程的解应该>0。
          3.在变换向量时入射向量不等于光线向量(长度不同)，不能混用。
          4.开始时记得计算光线向量(输入所给并非光线向量，而是这上面的一个点)。

#include <stdio.h>

#include <stdlib.h>

#include <stdbool.h>

#include <math.h>

struct b

{

 double x,y,z;

 double r;

 bool pass;

}ball[51],now,line,next,zero;

int n;

double dis(struct b a1,struct b a2)

{

 return(sqrt((a1.x-a2.x)*(a1.x-a2.x)+(a1.y-a2.y)*(a1.y-a2.y)+(a1.z-a2.z)*(a1.z-a2.z)));

}

struct b equation(double a,doubleb,double c) //解方程ak^2+bk+c=0和对应的交点

{

 double k1,k2;

 struct b x1,x2;

 x1.pass=false;

 k1=(-b+sqrt(b*b-4*a*c))/(2*a);

 x1.x=k1*line.x+now.x;

 x1.y=k1*line.y+now.y;

 x1.z=k1*line.z+now.z;

 k2=(-b-sqrt(b*b-4*a*c))/(2*a);

 x2.x=k2*line.x+now.x;

 x2.y=k2*line.y+now.y;

 x2.z=k2*line.z+now.z;

 if ((k1==k2 || dis(now,x1)<dis(now,x2)) && k1>0)

   {

   x1.pass=true;

   return x1;

   }

 else if (k2>0)

   {

   x2.pass=true;

   return x2;

   }

 return x1;

}

struct b find(int k) //判断交点

{

 double a,b,c;

 struct b re;

 if (k==2)

   k=2;

 re.pass=false;

 a=line.x*line.x+line.y*line.y+line.z*line.z;

 b=2*line.x*(now.x-ball[k].x)+2*line.y*(now.y-ball[k].y)+2*line.z*(now.z-ball[k].z);

 c=(now.x-ball[k].x)*(now.x-ball[k].x)+(now.y-ball[k].y)*(now.y-ball[k].y)+(now.z-ball[k].z)*(now.z-ball[k].z)-ball[k].r*ball[k].r;

 if (a!=0 && b*b>=4*a*c)

   re=equation(a,b,c);

 else if (a==0 && b!=0 && (-c)/b>0)

   {

      re.x=(-c)/b*line.x+now.x;

   re.y=(-c)/b*line.y+now.y;

   re.z=(-c)/b*line.z+now.z;

   re.pass=true;

      }

 return re;

}

double point(struct b a1,struct b a2)

{

 return a1.x*a2.x+a1.y*a2.y+a1.z*a2.z;

}

struct b makev(struct b a1,struct ba2)

{

 struct b re;

 re.x=a2.x-a1.x;

 re.y=a2.y-a1.y;

 re.z=a2.z-a1.z;

 return re;

}

struct b change(struct b a1,double i)

{

 struct b re;

 re.x=a1.x*i;

 re.y=a1.y*i;

 re.z=a1.z*i;

 return re;

}

struct b add(struct b a1,struct b a2)

{

 struct b re;

 re.x=a1.x+a2.x;

 re.y=a1.y+a2.y;

 re.z=a1.z+a2.z;

  returnre;

}

void turn(int k)

{

 struct b v1,v2,v3,v4,v5;

 double i;

 v1=makev(next,now);  //从交点向起始点指一个向量

 v2=makev(ball[k],next);  //从球心向交点指一个向量

 i=point(v1,v2)/(dis(zero,v2)*dis(zero,v2));

 v3=change(v2,i);  //以上两步求出投影向量v3

 v4=makev(v1,v3);  //从v1向v3指一个向量(实际就是减法)

 v5=change(v4,2);  //将上面的向量延伸

 line=add(v1,v5);  //计算反射向量

 return ;

}

void work()

{

 int t,i,ballk,lastball=0;

 struct b k;

 t=1;

 while (t<=11)

   {

      next.pass=false;

      for(i=1;i<=n;i++)

        if (i!=lastball) //因为起点在球上，特判

          {

         k=find(i); //对第i个球进行判断

          if (k.pass && (!next.pass ||dis(k,now)<dis(next,now))) //判断交点更新

            {

             next=k;

             next.pass=true;

             ballk=i;

             }

        }

   if (!next.pass)

     break;

   turn(ballk); //光线反射

   now=next;  //更改入射点

   if (t<=10)

     if (t==1)

        printf("%d",ballk);

     else

        printf(" %d",ballk);

   else

     printf(" etc.");

   t++;

   lastball=ballk; //记录上一个到达的球

      }

 printf("\n");

 return ;

}

void init()

{

 int i;

 scanf("%d",&n);

 for (i=1;i<=n;i++)

   scanf("%lf%lf%lf%lf",&ball[i].x,&ball[i].y,&ball[i].z,&ball[i].r);

 scanf("%lf%lf%lf%lf%lf%lf",&now.x,&now.y,&now.z,&line.x,&line.y,&line.z);

 line.x-=now.x;

 line.y-=now.y;

 line.z-=now.z; //计算向量

 work();

 return ;

}

int main()

{

 init();

 return 0;

}

 

 

 

计算几何 Part.7---旋转卡壳

1.

 http://poj.org/problem?id=3608

 

凸多边形间最小距离

给定两个非连接（比如不相交）的凸多边形 P 和 Q，目标是找到拥有最小距离的点对 (p,q) （p 属于 P 且 q 属于Q）。 

事实上，多边形非连接十分重要，因为我们所说的多边形包含其内部。如果多边形相交，那么最小距离就变得没有意义了。然而，这个问题的另一个版本，凸多边形顶点对间最小距离对于相交和非相交的情况都有解存在。 

回到我们的主问题：直观的，确定最小距离的点不可能包含在多边形的内部。与最大距离问题相似，我们有如下结论： 

两个凸多边形 P 和 Q 之间的最小距离由多边形间的对踵点对确立。存在凸多边形间的三种多边形间的对踵点对，因此就有三种可能存在的最小距离模式：

1. “顶点-顶点”的情况
2. “顶点-边”的情况
3. “边-边”的情况

换句话说，确定最小距离的点对不一定必须是顶点。
   
给定结果，一个基于旋转卡壳的算法自然而然的产生了: 
考虑如下的算法，算法的输入是两个分别有 m和 n个顺时针给定顶点的凸多边形 P和 Q。

1. 计算 P 上 y 坐标值最小的顶点（称为 yminP ） 和 Q 上 y 坐标值最大的顶点（称为 ymaxQ）。
2. 为多边形在 yminP 和 ymaxQ 处构造两条切线 LP 和 LQ 使得他们对应的多边形位于他们的右侧。 此时 LP和 LQ 拥有不同的方向， 并且 yminP 和 ymaxQ 成为了多边形间的一个对踵点对。
3. 计算距离(yminP,ymaxQ) 并且将其维护为当前最小值。
4. 顺时针同时旋转平行线直到其中一个与其所在的多边形的边重合。
5. 如果只有一条线与边重合， 那么只需要计算“顶点-边”对踵点对和“顶点-顶点”对踵点对距离。 都将他们与当前最小值比较， 如果小于当前最小值则进行替换更新。 如果两条切线都与边重合， 那么情况就更加复杂了。 如果边“交叠”， 也就是可以构造一条与两条边都相交的公垂线（但不是在顶点处相交）， 那么就计算“边-边”距离。 否则计算三个新的“顶点-顶点”对踵点对距离。 所有的这些距离都与当前最小值进行比较， 若小于当前最小值则更新替换。
6. 重复执行步骤4和步骤5， 直到新的点对为(yminP,ymaxQ)。
7. 输出最大距离。

旋转卡壳模式保证了所有的对踵点对（和所有可能的子情况）都被考虑到。此外，整个算法拥有现行的时间复杂度，因为（除了初始化），只有与顶点数同数量级的操作步数需要执行。 

最小距离和最大距离的问题表明了旋转卡壳模型可以用在不同的条件下（与先前的直径和宽度问题比较）。这个模型可以应用于两个多边形的问题中。 
“最小盒子”问题（最小面积外接矩形）通过同一多边形上两个正交切线集合展示了另一种条件下旋转卡壳的应用

 

#include <iostream>

#include <string.h>

#include <algorithm>

#include <stdio.h>

#include <math.h>

 

using namespace std;

 

const int N=50000;

const double eps=1e-9;

const double INF=1e99;

 

struct Point

{

   double x,y;

};

 

Point P[N],Q[N];

 

double cross(Point A,Point B,Point C)

{

   return (B.x-A.x)*(C.y-A.y)-(B.y-A.y)*(C.x-A.x);

}

 

double dist(Point A,Point B)

{

   return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y));

}

 

double multi(Point A,Point B,Point C)

{

   return (B.x-A.x)*(C.x-A.x)+(B.y-A.y)*(C.y-A.y);

}

 

//顺时针排序

void anticlockwise(Point p[],int n)

{

   for(int i=0;i<n-2;i++)

    {

       double tmp=cross(p[i],p[i+1],p[i+2]);

       if(tmp>eps) return;

       else if(tmp<-eps)

       {

           reverse(p,p+n);

           return;

       }

    }

}

 

//计算C点到直线AB的最短距离

double Getdist(Point A,Point B,Point C)

{

   if(dist(A,B)<eps) return dist(B,C);

   if(multi(A,B,C)<-eps) return dist(A,C);

   if(multi(B,A,C)<-eps) return dist(B,C);

   return fabs(cross(A,B,C)/dist(A,B));

}

 

//求一条直线的两端点到另外一条直线的距离，反过来一样，共4种情况

double MinDist(Point A,Point B,PointC,Point D)

{

   return min(min(Getdist(A,B,C),Getdist(A,B,D)),min(Getdist(C,D,A),Getdist(C,D,B)));

}

 

double Solve(Point P[],Point Q[],int n,intm)

{

   int yminP=0,ymaxQ=0;

   for(int i=0;i<n;i++)

      if(P[i].y<P[yminP].y)

         yminP=i;

   for(int i=0;i<m;i++)

      if(Q[i].y>Q[ymaxQ].y)

         ymaxQ=i;

    P[n]=P[0];

   Q[m]=Q[0];

   double tmp,ans=INF;

   for(int i=0;i<n;i++)

    {

       while(tmp=cross(P[yminP+1],Q[ymaxQ+1],P[yminP])-cross(P[yminP+1],Q[ymaxQ],P[yminP])>eps)

           ymaxQ=(ymaxQ+1)%m;

       if(tmp<-eps) ans=min(ans,Getdist(P[yminP],P[yminP+1],Q[ymaxQ]));

       else        ans=min(ans,MinDist(P[yminP],P[yminP+1],Q[ymaxQ],Q[ymaxQ+1]));

       yminP=(yminP+1)%n;

    }

   return ans;

}

 

int main()

{

   int n,m;

   while(cin>>n>>m)

    {

       if(n==0&&m==0) break;

       for(int i=0;i<n;i++)

          cin>>P[i].x>>P[i].y;

       for(int i=0;i<m;i++)

          cin>>Q[i].x>>Q[i].y;

       anticlockwise(P,n);

       anticlockwise(Q,m);

       printf("%.5lf\n",min(Solve(P,Q,n,m),Solve(Q,P,m,n)));

    }

   return 0;

}

 

 

计算几何 Part.8---极角排序

1.

http://poj.org/problem?id=1106

 

题意：给一个半圆的半径和圆心坐标，再给平面上的n个点，这个半圆可以绕圆心任意旋转，问半圆最多能覆盖多少个点？

 

分析：本题可以这样做，我们只需要考虑到圆心距离小于或者等于半径的那些点，因为大于半径的点一定不能覆盖到，那么这

样我们把符合条件的点全部存入一个数组p[]，那么然后就二重循环枚举每一个点与圆心所连的直线的的左侧有多少个点，记录

最大即可。

 

 

#include <iostream>

#include <string.h>

#include <algorithm>

#include <stdio.h>

#include <math.h>

 

using namespace std;

 

const int N=50000;

const double eps=1e-9;

 

struct Point

{

   double x,y;

};

 

Point p[N];

 

double cross(Point A,Point B,Point C)

{

   return (B.x-A.x)*(C.y-A.y)-(B.y-A.y)*(C.x-A.x);

}

 

double dist(Point A,Point B)

{

   return (A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y);

}

 

int main()

{

   int k,n;

   Point t,cir;

   double d,r;

   while(cin>>cir.x>>cir.y>>r)

    {

        k=0;

       if(r<0) break;

       cin>>n;

       for(int i=0;i<n;i++)

       {

           cin>>t.x>>t.y;

           d=dist(t,cir);

           if(d<r*r||fabs(d-r*r)<eps) //表示点在圆内或者在圆上

           {

                p[k].x=t.x;

                p[k].y=t.y;

                k++;

           }

       }

       int maxval=0;

       for(int i=0;i<k;i++)

       {

           int count=1;

           for(int j=0;j<k;j++)

           {

               if(i!=j&&cross(cir,p[i],p[j])>=0)

                    count++;

           }

           if(count>maxval) maxval=count;

       }

       cout<<maxval<<endl;

    }

   return 0;

}

 

 

 

END

 

推荐博客：

1. http://blog.csdn.net/zsc09_leaf/article/details/6331809

2. http://blog.csdn.net/yew1eb/article/details/20706115

3. http://blog.csdn.net/clasky/article/details/9990235