计算几何-HPI

 

This article is made by Jason-Cow.
Welcome to reprint.
But please post the article's address.

 

 

 

在线笛卡尔坐标系绘图网站

核心思想

　　1.积角排序

　　2.双端队列维护半平面交

木有啦～～

 1 int HPI(L*l,int n,D*ans){
 2   int head,tail,m=0;D*P=new D[n];L*q=new L[n];
 3   sort(l+1,l+n+1),q[head=tail=0]=l[1];
 4   for(int i=2;i<=n;i++){
 5     while(head<tail && !Left(l[i],P[tail-1]))tail--;
 6     while(head<tail && !Left(l[i],P[head]))  head++;
 7     q[++tail]=l[i];//双端队列<(￣3￣)>
 8     if(fabs(Cross(q[tail].v,q[tail-1].v))<eps){
 9       tail--;//判断Cross(q[tail].v,q[tail-1].v)==0♪(^∇^*)
10       if(Left(q[tail],l[i].P))q[tail]=l[i];
11     }
12     if(head<tail)P[tail-1]=Intersect(q[tail-1],q[tail]);
13   }
14   while(head<tail && !Left(q[head],P[tail-1]))tail--;
15   if(tail-head<=1)return 0;
16   P[tail]=Intersect(q[tail],q[head]);
17   for(int i=head;i<=tail;i++)ans[++m]=P[i];
18   return m;
19 }

多边形相关

db Area(D*R,int n){
    db S=0.0;
    for(int i=1;i<n;i++)S+=Cross(R[i]-R[1],R[i+1]-R[1]);
    return S/2;
}

db Length(D*R,int n){
    db C=0.0;
    for(int i=2;i<=n;i++)C+=Dis(R[i],R[i-1]);
    return C+Dis(R[n],R[1]);
}

 

 

主程序

 1 L l[500];D A[500];
 2 int main(){
 3   for(int n;scanf("%d",&n)&&n;){
 4     for(int i=1;i<=n;i++){
 5       db a,b,c,d;
 6       scanf("%lf%lf%lf%lf",&a,&b,&c,&d);
 7       D A(a,b),B(c,d);
 8       l[i]=L(A,B-A);
 9     }
10     int m=HPI(l,n,A);
11     cout<<"m="<<m<<endl;
12     for(int i=1;i<=m;i++)printf("(%.4lf,%.4lf)\n",A[i].x,A[i].y);
13     printf("S=%.2lf\n",Area(A,m));
14     printf("C=%.2lf\n",Length(A,m));
15   }
16   return 0;
17 }

 

 

读入，如图（橙、红、蓝、绿）

4
2 1 1 2
1 2 0 1
1 0 2 1
0 1 1 0

输出

m=4
(0.0000,1.0000)
(1.0000,0.0000)
(2.0000,1.0000)
(1.0000,2.0000)
S=2.00
C=5.66

读入

 

6
1 2 0 1
0 1.8 0.8 0.2
0 1 1 0
1 0 2 1
3 0 1 2
0.6 1.6 0 0.7

 

输出

m=6
(0.6000,1.6000)
(0.3143,1.1714)
(0.8000,0.2000)
(1.0000,0.0000)
(2.0000,1.0000)
(1.0000,2.0000)
S=1.76
C=5.28

完整代码

  1 #include <algorithm>
  2 #include <iostream>
  3 #include <cstring>
  4 #include <cstdlib>
  5 #include <cstdio>
  6 #include <vector>
  7 #include <cmath>
  8 #include <queue>
  9 #include <map>
 10 #include <set>
 11 using namespace std;
 12 #define sqr(x) ((x)*(x))
 13 #define RG register
 14 #define op operator
 15 #define IL inline
 16 typedef double db;
 17 typedef bool bl;
 18 const db pi=acos(-1.0),eps=1e-10;
 19 struct D{
 20   db x,y;
 21   D(db x=0.0,db y=0.0):x(x),y(y){}
 22 };
 23 typedef D V;//不知道什么鬼，emacs 将 operator -> op 会萎掉  对！不！齐！
 24 bl operator<(D A,D B){return A.x<B.x||(A.x==B.x&&A.y<B.y);}
 25 V operator+(V A,V B){return V(A.x+B.x,A.y+B.y);}
 26 V operator-(V A,V B){return V(A.x-B.x,A.y-B.y);}
 27 V operator*(V A,db N){return V(A.x*N,A.y*N);}
 28 V operator/(V A,db N){return V(A.x/N,A.y/N);}
 29 
 30 db Ang(db x){return(x*180.0/pi);}
 31 db Rad(db x){return(x*pi/180.0);}
 32 V Rotate(V A,db a){return V(A.x*cos(a)-A.y*sin(a),A.x*sin(a)+A.y*cos(a));}
 33 db Dis(D A,D B){return sqrt(sqr(A.x-B.x)+sqr(A.y-B.y));}
 34 db Cross(V A,V B){return A.x*B.y-A.y*B.x;}
 35 
 36 db Area(D*R,int n){
 37     db S=0.0;
 38     for(int i=1;i<n;i++)S+=Cross(R[i]-R[1],R[i+1]-R[1]);
 39     return S/2;
 40 }
 41 
 42 db Length(D*R,int n){
 43     db C=0.0;
 44     for(int i=2;i<=n;i++)C+=Dis(R[i],R[i-1]);
 45     return C+Dis(R[n],R[1]);
 46 }
 47 
 48 struct L{
 49   D P,v;
 50   db a;
 51   L(){}
 52   L(D P,V v):P(P),v(v){a=atan2(v.y,v.x);}
 53   bool operator<(const L x)const{return a<x.a;}
 54 };
 55 
 56 D Intersect(L a,L b){
 57   V u=a.P-b.P;
 58   return a.P+a.v*(Cross(b.v,u)/Cross(a.v,b.v));
 59 }
 60 
 61 bool Left(L l,D A){
 62   return Cross(l.v,A-l.P)>0;
 63 }
 64 
 65 int HPI(L*l,int n,D*ans){
 66   int head,tail,m=0;D*P=new D[n];L*q=new L[n];
 67   sort(l+1,l+n+1),q[head=tail=0]=l[1];
 68   for(int i=2;i<=n;i++){
 69     while(head<tail && !Left(l[i],P[tail-1]))tail--;
 70     while(head<tail && !Left(l[i],P[head]))  head++;
 71     q[++tail]=l[i];//双端队列<(￣3￣)>
 72     if(fabs(Cross(q[tail].v,q[tail-1].v))<eps){
 73       tail--;//判断Cross(q[tail].v,q[tail-1].v)==0♪(^∇^*)
 74       if(Left(q[tail],l[i].P))q[tail]=l[i];
 75     }
 76     if(head<tail)P[tail-1]=Intersect(q[tail-1],q[tail]);
 77   }
 78   while(head<tail && !Left(q[head],P[tail-1]))tail--;
 79   if(tail-head<=1)return 0;
 80   P[tail]=Intersect(q[tail],q[head]);
 81   for(int i=head;i<=tail;i++)ans[++m]=P[i];
 82   return m;
 83 }
 84 
 85 L l[500];
 86 D A[500];
 87 int main(){
 88   cout<<4*sqrt(2)<<endl;
 89   for(int n;scanf("%d",&n)&&n;){
 90     for(int i=1;i<=n;i++){
 91       db a,b,c,d;
 92       scanf("%lf%lf%lf%lf",&a,&b,&c,&d);
 93       D A(a,b),B(c,d);
 94       l[i]=L(A,B-A);
 95     }
 96     int m=HPI(l,n,A);
 97     cout<<"m="<<m<<endl;for(int i=1;i<=m;i++)printf("(%.4lf,%.4lf)\n",A[i].x,A[i].y);
 98     printf("S=%.2lf\n",Area(A,m));
 99     printf("C=%.2lf\n",Length(A,m));
100   }
101   return 0;
102 }
HPI

 

posted @ 2017-03-25 17:26  墨鳌  阅读(178)  评论(0)  编辑  收藏  举报