计算几何模板 （更新中）

1 判断线段是否非规范相交

下面模板是根据跨立实验得到的结果

http://download.csdn.net/download/welon123/3813929

#define eps 1e-8
struct point
{
double x;
double y;
};
double multi(point p0, point p1, point p2)//j计算差乘
{   return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);}
bool is_cross(point s1,point e1,point s2,point e2)//判断线段是否相交（非规范相交）
{
   return max(s1.x,e1.x) >= min(s2.x,e2.x)&&
          max(s2.x,e2.x) >= min(s1.x,e1.x)&&
          max(s1.y,e1.y) >= min(s2.y,e2.y)&&
          max(s2.y,e2.y) >= min(s1.y,e1.y)&&
          multi(s1,e1,s2)*multi(s1,e1,e2) <= 0&&
          multi(s2,e2,s1)*multi(s2,e2,e1) <= 0;
}

2 判断直线与线段相交

bool across(point &a,point &b,point &c,point &d)//直线ab和线段cd是否相交
{
double p=xmulit(a,b,c),p1=xmulit(a,b,d);
if( fabs(p1) <= eps || fabs(p) <= eps ) return true;
if( p*p1 <= 1e-8 )
return true;
return false;
}
bool is_equal(point &a,point &b)//判断点a和点b是否相等
{
return (fabs(a.x-b.x) <= eps) && (fabs(a.y-b.y) <=eps);
}

3 判断结果是否为0

bool zero(double a)//判断结果是否为0
{
return fabs(a) <= eps;
}

4 判断两直线是否平行

bool parallel(line &u,line &v)//判断直线u和直线v是否平行
{
return zero((u.a.x-u.b.x)*(v.a.y-v.b.y) - (v.a.x-v.b.x)*(u.a.y-u.b.y));
}

bool parallel(point &u1,point &u2,point &v1,point &v2)//判断直线u1 u2 and v1 v2 是否平行
{
return zero((u1.x-u2.x)*(v1.y-v2.y)-(v1.x-v2.x)*(u1.y-u2.y));
}

5 判断点是否在直线上

bool dot_in_line(point &a,point &b,point &c)
{
return parallel(b,a,a,c);
}

6 求两直线的交点 注意是直线不是线段

point intersection(line &u,line &v)
{
point ret=u.a;
double t=((u.a.x-v.a.x)*(v.a.y-v.b.y) - (u.a.y-v.a.y)*(v.a.x-v.b.x))/((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));
ret.x+=(u.b.x-u.a.x)*t;
ret.y+=(u.b.y-u.a.y)*t;
return ret;
}

判断线段是否相交

bool on_segment(point pi,point pj,point pk)//判断点pk是否在线段pi, pj上
{
    if(xmulit(pi, pj, pk)==0)
    {
        if(pk.x>=min(pi.x,pj.x)&&pk.x<=max(pi.x,pj.x)&&pk.y>=min(pi.y,pj.y)&&pk.y<=max(pi.y,pj.y))
            return true;
    }
    return false;
}
bool segments_intersect(point p1,point p2,point p3,point p4)//判断线段是否相交
{
    double d1=xmulit(p3,p4,p1);
    double d2=xmulit(p3,p4,p2);
    double d3=xmulit(p1,p2,p3);
    double d4=xmulit(p1,p2,p4);
    if(d1*d2<0&&d3*d4<0)
        return true;
    else if(d1==0&&on_segment(p3,p4,p1))
        return true;
    else if(d2==0&&on_segment(p3,p4,p2))
        return true;
    else if(d3==0&&on_segment(p1,p2,p3))
        return true;
    else if(d4==0&&on_segment(p1,p2,p4))
        return true;
    return false;
}

7 求凸包的graham扫描法

输入n表示有多少个点

接下来n个点的坐标

example：

输入：

9

200 400

300 400

300 300

400 300

400 400

500 400

500 200

350 200

200 200

输出：

500.000000 200.000000
500.000000 400.000000
400.000000 400.000000
300.000000 400.000000
200.000000 400.000000
200.000000 200.000000

#include <iostream>
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <cmath>
using namespace std;
#define eps 1e-6
#define PI 3.14159265
struct point
{
double x;
double y;
}po[1500],temp;
int n,pos;
bool zero(double a)
{
return fabs(a) < eps;
}
double dis(point &a,point &b)//返回两点之间距离的平方
{
return (a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y);
}
double across(point &a,point &b,point &c)//求a b and a c 的X积
{
return (b.x-a.x)*(c.y-a.y) - (b.y-a.y)*(c.x-a.x);
}
int cmp(const void *a,const void *b)
{
return across(po[0],*(point*)a,*(point*)b) > 1e-8 ? -1 : 1;
}
int select()
{
int i,j,k=1;
for(i=2;i<n;i++)
{
if(zero(across(po[0],po[k],po[i])))
{
if(dis(po[0],po[k]) < dis(po[0],po[i]))
po[k]=po[i];
}
else
po[++k]=po[i];
}
return k+1;
}
int graham(int num)
{
int i,j,k=2;
//
po[num]=po[0];//fangbian 
num++;
for(i=3;i<num;i++)
{
while(across(po[k-1],po[k],po[i]) < -eps)
{k--;}
po[++k]=po[i];//就这个循环结束，不需要了！
}
for(i=0;i<k;i++)
printf("%lf %lf\n",po[i].x,po[i].y);
return 0;
}
int main()
{
int i,j,k;
point my_temp;
while(scanf("%d",&n)!=EOF)
{
scanf("%lf%lf",&po[0].x,&po[0].y);
temp=po[0];
pos=0;
for(i=1;i<n;i++)
{
scanf("%lf%lf",&po[i].x,&po[i].y);
if(po[i].y < temp.y)
temp=po[i],pos=i;
}
my_temp=po[0];
po[0]=po[pos];
po[pos]=my_temp;
qsort(po+1,n-1,sizeof(po[0]),cmp);
graham(select());
}
return 0;
}

最近平面点对 见其他博客

给定一个多边形判断是否为凸包

要求逆时针给定多边形的各个顶点

放在po数组中

struct point
{
int x;
int y;
}po[1000];
int across(point &a,point &b,point &c)
{
return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x);
}
int is_convex(int num,point *p)
{
int i;
p[num]=p[0];
num++;
p[num]=po[1];
num++;
int count=0;
for(i=2;i < num;i++){
if(across(p[i-2],p[i-1],p[i]) < 0)
return 0;
}
return 1;
}

已知正方形的一条边的两点坐标，求另外两点的坐标

已知： (x1,y1)  (x2,y2)
则：   x3=x1+(y1-y2)   y3= y1-(x1-x2)
x4=x2+(y1-y2)   y4= y2-(x1-x2)
或
x3=x1-(y1-y2)   y3= y1+(x1-x2)
x4=x2-(y1-y2)   y4= y2+(x1-x2)

扩展KMP 模板

#include<iostream>
#include<string>
using namespace std;
const int MM=100005;
int next[MM],extand[MM];
char S[MM],T[MM];
void GetNext(const char *T){
     int len=strlen(T),a=0;
     next[0]=len;
     while(a<len-1 && T[a]==T[a+1]) a++;
     next[1]=a;
     a=1;
     for(int k=2;k<len;k++){
         int p=a+next[a]-1,L=next[k-a];
         if( (k-1)+L >= p){
             int j = (p-k+1)>0 ? (p-k+1) : 0;
             while(k+j<len && T[k+j]==T[j]) j++;
             next[k]=j;
             a=k; 
         } 
         else
             next[k]=L; 
     } 
} 
void GetExtand(const char *S,const char *T){
     GetNext(T);
     int slen=strlen(S),tlen=strlen(T),a=0; 
     int MinLen = slen < tlen ? slen : tlen;
     while(a<MinLen && S[a]==T[a]) a++;
     extand[0]=a;
     a=0;
     for(int k=1;k<slen;k++){
         int p=a+extand[a]-1, L=next[k-a];
         if( (k-1)+L >= p){
             int j= (p-k+1) > 0 ? (p-k+1) : 0;
             while(k+j<slen && j<tlen && S[k+j]==T[j]) j++;
             extand[k]=j;
             a=k; 
         }
         else 
             extand[k]=L; 
     } 
} 
int main(){
    while(scanf("%s%s",S,T)==2){
         GetExtand(S,T);
         for(int i=0;i<strlen(T);i++)
             printf("%d ",next[i]);
         puts("");
         for(int i=0;i<strlen(S);i++)
             printf("%d ",extand[i]);
         puts(""); 
    } 
    return 0;
}

欧拉四面体公式 ：HDU 1411 

注意输入点边的长度依次是：AB，AC，AD，BC，BD，CD

#include <iostream>
#include <string.h>
#include <cmath>
#include <stdio.h>
#include <algorithm>
using namespace std;
double l,m,n,p,q,r;
double area()
{
    return sqrt((p*p*(q*q*r*r-((q*q+r*r-l*l)/2)*((q*q+r*r-l*l)/2))
    +((p*p+r*r-m*m)/2)*(((p*p+q*q-n*n)/2)*((q*q+r*r-l*l)/2)-((p*p+r*r-m*m)/2)*q*q)
    -((p*p+q*q-n*n)/2)*(((p*p+q*q-n*n)/2)*r*r-((p*p+r*r-m*m)/2)*(q*q+r*r-l*l)/2))/36);
}
int main()
{
    while(scanf("%lf%lf%lf%lf%lf%lf",&n,&m,&p,&l,&q,&r)!=EOF)
    {
      printf("%.4lf\n",area());
    }
    return 0;
}

求两个圆相交部分的面积：HDU 3264

函数area

 a：第一个圆心 r1 第一个圆的半径

 b：第二个圆心 r2 第二个圆的半径             

struct point
{
 double x;
 double y;
};
double area(point a,double r1,point b,double r2)
{
 double d=sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));//圆心距
 if(r1>r2)
 {
 double temp=r1;
 r1=r2;
 r2=temp;
 }//r1取小
 if(r1+r2<=d)
 return 0;//相离
 else if(r2-r1>=d)
 return pi*r1*r1;//内含
 else
 {
 double a1=acos((r1*r1+d*d-r2*r2)/(2.0*r1*d));
 double a2=acos((r2*r2+d*d-r1*r1)/(2.0*r2*d));
 return (a1*r1*r1+a2*r2*r2-r1*d*sin(a1));
 }//相交
}

多边形重心模板适合任意多边形，但是多边形的边要顺序给出！

point gravity(point *p, int n)
{
double area = 0;
point center;
center.x = 0;
center.y = 0;
for (int i = 0; i < n-1; i++)
{
   area += (p[i].x*p[i+1].y - p[i+1].x*p[i].y)/2;
   center.x += (p[i].x*p[i+1].y - p[i+1].x*p[i].y) * (p[i].x + p[i+1].x);
   center.y += (p[i].x*p[i+1].y - p[i+1].x*p[i].y) * (p[i].y + p[i+1].y);
}
area += (p[n-1].x*p[0].y - p[0].x*p[n-1].y)/2;
center.x += (p[n-1].x*p[0].y - p[0].x*p[n-1].y) * (p[n-1].x + p[0].x);
center.y += (p[n-1].x*p[0].y - p[0].x*p[n-1].y) * (p[n-1].y + p[0].y);
center.x /= 6*area;
center.y /= 6*area;
return center;
}

点到线段之间距离 线段到线段之间距离

double point_to_seg(point a,point b,point c)//c 到直线a-b的距离  
{  
    point ab,ac;  
    ab.x=b.x-a.x;  
    ab.y=b.y-a.y;  
    ac.x=c.x-a.x;  
    ac.y=c.y-a.y;  
    double f=dot(ab,ac);  
    if(f<0)return dis(a,c);  
    double f1=dot(ab,ab);  
    if(f>f1)return dis(b,c);  
    f=f/f1;  
    point d;  
    d.x=a.x+ab.x*f;  
    d.y=a.y+ab.y*f;  
    return dis(d,c);  
}  
double seg_to_seg(point a1,point b1,point a2,point b2)  // a1-b1    a2-b2
{  
    return min(min(point_to_seg(a1,b1,a2),point_to_seg(a1,b1,b2)),min(point_to_seg(a2,b2,a1),point_to_seg(a2,b2,b1)));  
}  

判断点是否在多边形内部a表示要判断的那个点，po是多边形（按一定顺序）点集，n多边形点的个数

int inpoto(point a,point *po,int n)//判断点是否在多边形的内部  
{  
    int i;  
    point b,c,d;  
    b.y=a.y;  
    b.x=1e15;//定义射线  
    int flag=0;  
    int count=0;  
    for(i=0;i<n;i++)  
    {  
        c = po[i];  
        d = po[i + 1];  
        if(on_segment(c,d,a))//该点在多边形的一条边上  
            return 1;  
        if(abs(c.y-d.y)<eps)  
            continue;  
        if(on_segment(a,b,c))//和顶点相交的情况，如果y值较大则取  
        {  
            if(c.y>d.y)  
                count++;  
        }  
        else if(on_segment(a,b,d))//和顶点相交的情况，如果y值较大则取  
        {  
            if(d.y>c.y)  
                count++;  
        }  
        else if(segments_intersect(a,b,c,d))//和边相交  
            count++;  
    }  
    return count%2;//当L和多边形的交点数目C是奇数的时候，P在多边形内，是偶数的话P在多边形外。  
}  

求任意多边形与圆的相交部分的面积

a b 多边形上所有边的两个顶点 c圆心 r 半径

#define eps 1e-8
struct Point
{
    double x,y;
};
double x_mult(Point sp, Point ep, Point op)
{
    return (sp.x-op.x)*(ep.y-op.y)-(sp.y-op.y)*(ep.x-op.x);
}
double cross(Point a,Point b,Point c)
{
    return (a.x-c.x)*(b.x-c.x)+(a.y-c.y)*(b.y-c.y);
}
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 cal_area(Point a,Point b,Point c,double r)//a b多边形上所有边的两个端点 c圆心 r半径
{
    double A,B,C,x,y,tS;
    A=dist(b,c);
    B=dist(a,c);
    C=dist(b,a);
    if(A<r&&B<r)
    return x_mult(a,b,c)/2;
    else if(A<r&&B>=r)
    {
        x=(cross(a,c,b)+sqrt(r*r*C*C-x_mult(a,c,b)*x_mult(a,c,b)))/C;
        tS=x_mult(a,b,c)/2;
        return asin(tS*(1-x/C)*2/r/B*(1-eps))*r*r/2+tS*x/C;
    }
    else if(A>=r&&B<r)
    {
        y=(cross(b,c,a)+sqrt(r*r*C*C-x_mult(b,c,a)*x_mult(b,c,a)))/C;
        tS=x_mult(a,b,c)/2;
        return asin(tS*(1-y/C)*2/r/A*(1-eps))*r*r/2+tS*y/C;
    }
    else if(fabs(x_mult(a,b,c))>=r*C||cross(b,c,a)<=0||cross(a,c,b)<=0)
    {
        if(cross(a,b,c)<0)
            if(x_mult(a,b,c)<0)
                return (-acos(-1.0)-asin(x_mult(a,b,c)/A/B*(1-eps)))*r*r/2;
            else return (acos(-1.0)-asin(x_mult(a,b,c)/A/B*(1-eps)))*r*r/2;
        else return asin(x_mult(a,b,c)/A/B*(1-eps))*r*r/2;
    }
    else
    {
        x=(cross(a,c,b)+sqrt(r*r*C*C-x_mult(a,c,b)*x_mult(a,c,b)))/C;
        y=(cross(b,c,a)+sqrt(r*r*C*C-x_mult(b,c,a)*x_mult(b,c,a)))/C;
        tS=x_mult(a,b,c)/2;
        return (asin(tS*(1-x/C)*2/r/B*(1-eps))+asin(tS*(1-y/C)*2/r/A*(1-eps)))*r*r/2+tS*((y+x)/C-1);
    }
}

·

坐标绕原点旋转公式新坐标公式(x,y)->(x1,y1)

x1=cos(angle)*x-sin(angle)*y;

y1=sin(angle)*x+cos(angle)*y;

点集的最小圆覆盖

/*http://blog.sina.com.cn/s/blog_5caa94a001015dr3.html*/
#include<stdio.h>
#include<math.h>
#include<string.h>
#define MAXN 20
struct pointset{
    double x, y;
};
const double precison=1.0e-8;
pointset maxcic, point[MAXN];
double radius;
int curset[MAXN], posset[3];
int set_cnt, pos_cnt;
inline double dis_2(pointset &from, pointset& to){
    return ((from.x-to.x)*(from.x-to.x)+(from.y-to.y)*(from.y-to.y));
}
int in_cic(int pt){
    if(sqrt(dis_2(maxcic, point[pt]))<radius+precison) return 1;
    return 0;
}
int cal_mincic(){
    if(pos_cnt==1 || pos_cnt==0)
        return 0;
    else if(pos_cnt==3){
        double A1, B1, C1, A2, B2, C2;
        int t0=posset[0], t1=posset[1], t2=posset[2];
        A1=2*(point[t1].x-point[t0].x);
        B1=2*(point[t1].y-point[t0].y);
        C1=point[t1].x*point[t1].x-point[t0].x*point[t0].x+
            point[t1].y*point[t1].y-point[t0].y*point[t0].y;
        A2=2*(point[t2].x-point[t0].x);
        B2=2*(point[t2].y-point[t0].y);
        C2=point[t2].x*point[t2].x-point[t0].x*point[t0].x+
            point[t2].y*point[t2].y-point[t0].y*point[t0].y;
        maxcic.y=(C1*A2-C2*A1)/(A2*B1-A1*B2);
        maxcic.x=(C1*B2-C2*B1)/(A1*B2-A2*B1);
        radius=sqrt(dis_2(maxcic, point[t0]));
    }
    else if(pos_cnt==2){
        maxcic.x=(point[posset[0]].x+point[posset[1]].x)/2;
        maxcic.y=(point[posset[0]].y+point[posset[1]].y)/2;
        radius=sqrt(dis_2(point[posset[0]], point[posset[1]]))/2;
    }
    return 1;
}
int mindisk(){
    if(set_cnt==0 || pos_cnt==3){
        return cal_mincic();
    }
    int tt=curset[--set_cnt];
    int res=mindisk();
    set_cnt++;
    if(!res || !in_cic(tt)){
        set_cnt--;
        posset[pos_cnt++]=curset[set_cnt];
        res=mindisk();
        pos_cnt--;
        curset[set_cnt++]=curset[0];
        curset[0]=tt;
    }
    return res;
}
int main(){
    int n;
    while(scanf("%d", &n)!=EOF){
        if(n==0) break;
        int i;
        for(i=0; i<n; i++)
            scanf("%lf %lf", &point[i].x, &point[i].y);
            if(n==1){
                maxcic.x=point[0].x;
                maxcic.y=point[0].y;
                radius=0;
                printf("%.2lf %.2lf %.2lf\n", maxcic.x, maxcic.y, radius);
                continue;
            }
            set_cnt=n; pos_cnt=0;
            for(i=0 ;i<n ;i++)  curset[i]=i;
            mindisk();
            printf("%.2lf %.2lf %.2lf\n", maxcic.x, maxcic.y, radius);
    }

    return 0;
}

求ab 和ac的夹角

这个求出来的结果只是两个向量之间夹角小的那个，如果想求顺时针的可以判断ab和cd的时针关系，如果ab在cd的逆时针方向，用PI*2去减

double angle(point a, point b, point c) {
  double ux = b.x - a.x, uy = b.y - a.y;
  double vx = c.x - a.x, vy = c.y - a.y;
  return acos((ux*vx + uy*vy) /
              sqrt((ux*ux + uy*uy) * (vx*vx + vy*vy)));
}

#include <string.h>
#include <stdio.h>
#include <algorithm>
#include <iostream>
#include <cmath>
using namespace std;
const double PI=acos(-1.0);
struct point{
    double x,y;
};
double angle(point a, point b, point c) {
  double ux = b.x - a.x, uy = b.y - a.y;
  double vx = c.x - a.x, vy = c.y - a.y;
  return acos((ux*vx + uy*vy) /
              sqrt((ux*ux + uy*uy) * (vx*vx + vy*vy)));
}
double cross(point &a,point &b,point &c){
    return (b.x-a.x)*(c.y-a.y)-(c.x-a.x)*(b.y-a.y);
}
int main(){
    int i,j,k;
    point a,b,c;
    c.x=c.y=0;
    double my_angle,ans;
    while(scanf("%lf %lf %lf %lf",&a.x,&a.y,&b.x,&b.y)!=EOF){
        my_angle=angle(c,a,b);
        if(cross(c,a,b)<0){
            printf("YES\n");
            my_angle=PI*2-my_angle;
        }
        ans=180*my_angle/PI;
        printf("%lf\n",ans);
    }
    return 0;
}

POJ 2954 PICK定理

/*整点多边形面积=多边形内含点数+边上整点数/2-1
外边界上的整点个数可以根据GCD来求，具体看代码
题意：给你一个整点三角的三个顶点让你求三角形内部整点个数
*/
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <iostream>
using namespace std;
struct point{
    int x,y;
}a,b,c;
int gcd(int a,int b){
    if(a<b) a^=b,b^=a,a^=b;
    if(b==0) return a;
    return gcd(b,a%b);
}
int cross(point &a,point &b,point &c){
    return (b.x-a.x)*(c.y-a.y)-(c.x-a.x)*(b.y-a.y);
}
int main(){
    int i,j,k,num;
    while(scanf("%d%d%d%d%d%d",&a.x,&a.y,&b.x,&b.y,&c.x,&c.y)){
        if(a.x==0 && b.x==0 && c.x==0 && a.y==0 && b.y==0 && c.y==0)
        return 0;
        num=gcd(abs(a.x-b.x),abs(a.y-b.y))+gcd(abs(a.x-c.x),abs(a.y-c.y))+gcd(abs(c.y-b.y),abs(c.x-b.x));
        k=cross(a,b,c);
        printf("%d\n",(abs(k)-num+2)/2);
    }
    return 0;
}