[原创]计算几何个人模板
都是自己整理的一些简单的计算几何模板，如果里面有问题请指出，毕竟我自己学的也不怎么好，由于计算几何部分东西很多，所以此文章不够全面，但是有新内容我会随时更新此文章的！
1. 前序
1. 注意舍入方式(0.5的舍入方向);防止输出-0.
2. 几何题注意多测试不对称数据.
3. 整数几何注意xmult和dmult是否会出界;
符点几何注意eps的使用.
4. 避免使用斜率;注意除数是否会为0.
5. 公式一定要化简后再代入.
6. 判断同一个2*PI域内两角度差应该是
abs(a1-a2)<beta||abs(a1-a2)>pi+pi-beta;
相等应该是
abs(a1-a2)<eps||abs(a1-a2)>pi+pi-eps;
7. 需要的话尽量使用atan2,注意:atan2(0,0)=0,
atan2(1,0)=pi/2,atan2(-1,0)=-pi/2,atan2(0,1)=0,atan2(0,-1)=pi.
8. cross product = |u|*|v|*sin(a)
dot product = |u|*|v|*cos(a)
9. (P1-P0)x(P2-P0)结果的意义:
正: <P0,P1>在<P0,P2>顺时针(0,pi)内
负: <P0,P1>在<P0,P2>逆时针(0,pi)内
0 : <P0,P1>,<P0,P2>共线,夹角为0或pi
10. 误差限缺省使用1e-？
几何公式：
三角形:
1. 半周长 P=(a+b+c)/2
2. 面积 S=aHa/2=absin(C)/2=sqrt(P(P-a)(P-b)(P-c))
3. 中线 Ma=sqrt(2(b^2+c^2)-a^2)/2=sqrt(b^2+c^2+2bccos(A))/2
4. 角平分线 Ta=sqrt(bc((b+c)^2-a^2))/(b+c)=2bccos(A/2)/(b+c)
5. 高线 Ha=bsin(C)=csin(B)=sqrt(b^2-((a^2+b^2-c^2)/(2a))^2)
6. 内切圆半径 r=S/P=asin(B/2)sin(C/2)/sin((B+C)/2)
=4Rsin(A/2)sin(B/2)sin(C/2)=sqrt((P-a)(P-b)(P-c)/P)
=Ptan(A/2)tan(B/2)tan(C/2)
7. 外接圆半径 R=abc/(4S)=a/(2sin(A))=b/(2sin(B))=c/(2sin(C))
四边形:
D1,D2为对角线,M对角线中点连线,A为对角线夹角
1. a^2+b^2+c^2+d^2=D1^2+D2^2+4M^2
2. S=D1D2sin(A)/2
(以下对圆的内接四边形)
3. ac+bd=D1D2
4. S=sqrt((P-a)(P-b)(P-c)(P-d)),P为半周长
正n边形:
R为外接圆半径,r为内切圆半径
1. 中心角 A=2PI/n
2. 内角 C=(n-2)PI/n
3. 边长 a=2sqrt(R^2-r^2)=2Rsin(A/2)=2rtan(A/2)
4. 面积 S=nar/2=nr^2tan(A/2)=nR^2sin(A)/2=na^2/(4tan(A/2))
圆:
1. 弧长 l=rA
2. 弦长 a=2sqrt(2hr-h^2)=2rsin(A/2)
3. 弓形高 h=r-sqrt(r^2-a^2/4)=r(1-cos(A/2))=atan(A/4)/2
4. 扇形面积 S1=rl/2=r^2A/2
5. 弓形面积 S2=(rl-a(r-h))/2=r^2(A-sin(A))/2
棱柱:
1. 体积 V=Ah,A为底面积,h为高
2. 侧面积 S=lp,l为棱长,p为直截面周长
3. 全面积 T=S+2A
棱锥:
1. 体积 V=Ah/3,A为底面积,h为高
(以下对正棱锥)
2. 侧面积 S=lp/2,l为斜高,p为底面周长
3. 全面积 T=S+A
棱台:
1. 体积 V=(A1+A2+sqrt(A1A2))h/3,A1.A2为上下底面积,h为高
(以下为正棱台)
2. 侧面积 S=(p1+p2)l/2,p1.p2为上下底面周长,l为斜高
3. 全面积 T=S+A1+A2
圆柱:
1. 侧面积 S=2PIrh
2. 全面积 T=2PIr(h+r)
3. 体积 V=PIr^2h
圆锥:
1. 母线 l=sqrt(h^2+r^2)
2. 侧面积 S=PIrl
3. 全面积 T=PIr(l+r)
4. 体积 V=PIr^2h/3
圆台:
1. 母线 l=sqrt(h^2+(r1-r2)^2)
2. 侧面积 S=PI(r1+r2)l
3. 全面积 T=PIr1(l+r1)+PIr2(l+r2)
4. 体积 V=PI(r1^2+r2^2+r1r2)h/3
球:
1. 全面积 T=4PIr^2
2. 体积 V=4PIr^3/3
球台:
1. 侧面积 S=2PIrh
2. 全面积 T=PI(2rh+r1^2+r2^2)
3. 体积 V=PIh(3(r1^2+r2^2)+h^2)/6
球扇形:
1. 全面积 T=PIr(2h+r0),h为球冠高,r0为球冠底面半径
2. 体积 V=2PIr^2h/3
2. 线段问题
A.给你两个线段（八个点），求这两个线段的位置关系，如果相交，求出交点。
代码：
#include<iostream>
using namespace std;
struct point
{
double x;
double y;
};
//判断是否共线叉乘
double gongxian(point a,point b,point c)
{
return((a.x-c.x)*(b.y-c.y)-(b.x-c.x)*(a.y-c.y));
}
//判断是否平行
int judge(point a,point b,point c,point d)
{
if((gongxian(a,b,c))==0&&(gongxian(b,c,d))==0)
return 0;//说明共线
else if( ((a.x-b.x)*(c.y-d.y)-(a.y-b.y)*(c.x-d.x))==0 )
return -1;//平行
else
return 1;//相交
}
//求两线段的交点利用叉乘
void jiaodian(point a,point b,point c,point d)
{
double a1,a2,b1,b2,c1,c2;
a1=a.y-b.y;b1=b.x-a.x;c1=a.x*b.y-b.x*a.y;
a2=c.y-d.y;b2=d.x-c.x;c2=c.x*d.y-d.x*c.y;
point p;
p.x=(c1*b2-c2*b1)/(a2*b1-a1*b2);
p.y=(a2*c1-a1*c2)/(a1*b2-a2*b1);
printf("两线段相交于( %.2f,%.2f)\n",p.x,p.y);
}
int main()
{
int n;
point p0,p1,p2,p3;
scanf("%d",&n);
while(n--)
{
//输出八个点确定两条线段
scanf("%lf%lf%lf%lf",&p0.x,&p0.y,&p1.x,&p1.y);
scanf("%lf%lf%lf%lf",&p2.x,&p2.y,&p3.x,&p3.y);
double t=judge(p0,p1,p2,p3);
if(t==0)
printf("两线段共线\n");
if(t==-1)
printf("两线段平行\n");
if(t==1)
jiaodian(p0,p1,p2,p3);
}
return 0;
}
B.点到线段的距离
#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
const double eps = 1e-6;
double xmult(double x0,double y0,double x1,double y1,double x2,double y2)
{
return(x0-x2)*(y1-y2)-(x1-x2)*(y0-y2);
}
double distance(double x1,double y1,double x2,double y2)
{
return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}
double disptoseg(double x0,double y0,double x1,double y1,double x2,double y2)//x0,y0是点
{
double x3=x0,y3=y0;
x3+=y1-y2,y3+=x2-x1;
if(xmult(x1,y1,x3,y3,x0,y0)*xmult(x2,y2,x3,y3,x0,y0)>eps)
return distance(x0,y0,x1,y1)<distance(x0,y0,x2,y2)?distance(x0,y0,x1,y1):distance(x0,y0,x2,y2);
return fabs(xmult(x0,y0,x1,y1,x2,y2))/distance(x1,y1,x2,y2);//如果单独使用这个就是点到直线的距离
}
int main()
{
double x0,y0,x1,y1,x2,y2;
while(~scanf("%lf%lf",&x0,&y0))
{
cout<<"请输入线段两点"<<endl;
scanf("%lf%lf%lf%lf",&x1,&y1,&x2,&y2);
double t= disptoseg(x0,y0,x1,y1,x2,y2)+eps;
printf("%.3f\n",t);
}
return 0;
}
3. 匹克定理
----在一个平面坐标系中计算图形内的点（都是整点）的个数----
//匹克定理S = A + E/2 - 1
//S表示多边形面积
//A表示多边形内部的点的个数
//E表示在多边形上的点的个数
//A = S - E/2 + 1
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
//最大公约数
int gcd(int a,int b)
{
if(!b)
return a;
else
return gcd(b, a%b);
}
//叉乘等于三角形面积的两倍
int area(int x1,int y1,int x2,int y2,int x3,int y3)
{
return abs((x1-x2) * (y2-y3) - (x2-x3) * (y1-y2));
}
//计算边上的点的个数
int edga(int x1,int y1,int x2,int y2,int x3,int y3)
{
int sum= 0;
sum+= gcd(abs(x1-x2), abs(y1-y2));
sum+= gcd(abs(x2-x3), abs(y2-y3));
sum+= gcd(abs(x3-x1), abs(y3-y1));
return sum;
}
int main()
{
int x1,y1,x2,y2,x3,y3,i,j,s,e,a;
while(scanf("%d%d%d%d%d%d", &x1, &y1, &x2, &y2, &x3, &y3))
{
if(x1==0&&y1==0 && x2==0 && y2==0&& x3==0 && y3==0)
{
break;
}
s=area(x1,y1,x2,y2,x3,y3)/2;
e=edga(x1,y1,x2,y2,x3,y3);
a=s-e/2+1;
printf("%d\n",a);
}
return 0;
}
4. 多边形面积
给你n条边数，再给出每条边对应的点（前提必须按顺序，如不按顺序就得凸包）；
#include<iostream>
#include<stdio.h>
#include<math.h>
using namespace std;
int a[1000005],b[1000005];
__int64 sum;
int main()
{
int tcase,i,j,len;
scanf("%d",&tcase);
while(tcase--)
{
cout<<"输入多边形n=：";
cin>>len;
for(i=0;i<len;i++)
{
scanf("%d%d",&a[i],&b[i]);
}
sum=0;a[len]=a[0];b[len]=b[0];
for(i=0;i<len;i++)
{
sum+=a[i]*b[i+1]-a[i+1]*b[i];
}
if(sum<0)
sum=-sum;
if(sum%2==0)
printf("%I64d\n",sum/2);
else
printf("%I64d.5\n",sum/2);
}
return 0;
}
求三角形面积可以用叉乘还可以用海伦公式（只适用于三角形）：
#include<iostream>
#include<stdio.h>
#include<math.h>
using namespace std;
//叉乘等于三角形面积的两倍
double area(double x1,double y1,double x2,double y2,double x3,double y3)
{
return (abs)((x1-x2) * (y2-y3) - (x2-x3) * (y1-y2));
}
int main()
{
int tcase;
double x1,y1,x2,y2,x3,y3,a,b,c;
scanf("%d",&tcase);
while(tcase--)
{
cout<<"输入三角形三边坐标分别为："<<endl;
cin>>x1>>y1>>x2>>y2>>x3>>y3;
a=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
b=sqrt((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2));
c=sqrt((x1-x3)*(x1-x3)+(y1-y3)*(y1-y3));
double p=(a+b+c)/2;//半周长
double s=sqrt(p*(p-a)*(p-b)*(p-c));
printf("叉积求面积得：%.2f\n",area(x1,y1,x2,y2,x3,y3)/2);
printf("海伦公式求面积得：%.2f\n",s);
}
return 0;
}
5. 多边形重心（质心）
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include <algorithm>
#include <iostream>
#define eps 1e-8
using namespace std;
struct point
{
double x , y ;
};
double Area( point p0 , point p1 ,point p2 )
{
double area = 0 ;
area = p0.x * p1.y + p1.x * p2.y + p2.x * p0.y - p1.x * p0.y - p2.x * p1.y - p0.x * p2.y;
return area / 2 ;
}
int main ()
{
point p0 , p1 , p2 ;
int n;
double sum_x , sum_y , sum_area , area;
while ( ~ scanf ("%d",&n ) )
{
sum_x = sum_y = sum_area = 0;
scanf ("%lf%lf",&p0.x,&p0.y);
scanf ("%lf%lf",&p1.x,&p1.y);
for ( int i = 2 ; i < n ; ++ i )
{
scanf ( "%lf%lf" , &p2.x , &p2.y ) ;
area = Area(p0,p1,p2) ;
sum_area += area ;
sum_x += (p0.x + p1.x + p2.x) * area ;
sum_y += (p0.y + p1.y + p2.y) * area ;
p1 = p2 ;
}
sum_x=sum_x / (sum_area * 3);
sum_y=sum_y / (sum_area * 3) ;
printf ( "%.3f %.3f\n" , sum_x +eps, sum_y+eps);
}
return 0 ;
}
6. 三角形各种心
#include <iostream>
#include <cmath>
#include <algorithm>
#include<stdio.h>
using namespace std;
const double EPS = 1e-8;
//叉乘等于三角形面积的两倍
double area(double x1,double y1,double x2,double y2,double x3,double y3)
{
return (abs)((x1-x2) * (y2-y3) - (x2-x3) * (y1-y2));
}
//两点间距离
double Distance( double x1, double y1, double x2, double y2)
{
return sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
}
//外心
void waixin(double x1,double y1,double x2,double y2,double x3, double y3)
{
double s = 0.5*fabs((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1)); //三角形面积
double a = Distance(x1, y1, x2, y2);
double b = Distance(x2, y2, x3, y3);
double c = Distance(x3, y3, x1, y1);
double r = a*b*c/(4*s);//外接圆半径
double C1 = (x1*x1+y1*y1-x2*x2-y2*y2)/2.0;
double C2 = (x1*x1+y1*y1-x3*x3-y3*y3)/2.0;
double x0 = (C1*(y1-y3)-C2*(y1-y2))/((x1-x2)*(y1-y3)-(x1-x3)*(y1-y2))+EPS;//外接圆圆心x轴坐标
double y0 = (C1*(x1-x3)-C2*(x1-x2))/((y1-y2)*(x1-x3)-(y1-y3)*(x1-x2))+EPS;//外接圆圆心y轴坐标
printf("该三角形的外接圆半径为： %.3f\n该三角形的外心为 （%.3f,%.3f）\n",r,x0,y0);
}
double CrossProduct(double x1, double y1, double x2, double y2)
{
return (x1*y2-x2*y1);
}
//垂心
void chuixin(double x1,double y1,double x2,double y2,double x3, double y3)
{
double B1 = y2 - y1;
double A1 = x2 - x1;
double C1 = y1*y3 - y2*y3 + x1*x3 - x2*x3;
double B2 = y3 - y1;
double A2 = x3 - x1;
double C2 = y1*y2 - y2*y3 + x1*x2 - x2*x3;
double x0 = CrossProduct(C2, B2, C1, B1)/CrossProduct(A1, B1, A2, B2)+EPS;//调节精度垂心x轴坐标
double y0 = CrossProduct(C1, A1, C2, A2)/CrossProduct(A1, B1, A2, B2)+EPS;//垂心y轴坐标
printf("该三角形的垂心为（%.3f ,%.3f）\n", x0, y0);
}
//重心
void zhongxin(double x1,double y1,double x2,double y2,double x3,double y3)
{
double A1=(x1+x2)/2;
double B1=(y1+y2)/2;
double C1=(x1+x3)/2;
double D1=(y1+y3)/2;
double t=((A1-C1)*(D1-y2)-(B1-D1)*(C1-x2))/((A1-x3)*(D1-y2)-(B1-y3)*(C1-x2));
double x0=A1+(x3-A1)*t+EPS;
double y0=B1+(y3-B1)*t+EPS;
printf("该三角形的重心为（%.3f ，%.3f）\n",x0,y0);
}
//内心
void neixin(double x1,double y1,double x2,double y2,double x3,double y3)
{
double m1=atan2(y2-y1,x2-x1);
double n1=atan2(y3-y1,x3-x1);
double A1=x1+cos((m1+n1)/2);
double B1=y1+sin((m1+n1)/2);
double m2=atan2(y1-y2,x1-x2);
double n2=atan2(y3-y2,x3-x2);
double C1=x2+cos((m2+n2)/2);
double D1=y2+sin((m2+n2)/2);
double t=((x1-x2)*(y2-D1)-(y1-y2)*(x2-C1))/((x1-A1)*(y2-D1)-(y1-B1)*(x2-C1));
double x0=(A1-x1)*t+EPS;
double y0=(B1-y1)*t+EPS;
double r=area(x1,y1,x2,y2,x0,y0)/Distance(x1,y1,x2,y2);
printf("该三角形的内接圆半径为： %.3f\n该三角形的内心为（%.3f ， %.3f）\n",r,x0,y0);
}
//费马点到三角形三顶点距离之和最小的点
void fermentpoint(double x1,double y1,double x2,double y2,double x3, double y3)
{
double step = fabs(x1) + fabs(y1) + fabs(x2) + fabs(y2) + fabs(x3) + fabs(y3);
int i, j, k;
double x0 = (x1+x2+x3)/3;
double y0 = (y1+y2+y3)/3;
while (step > EPS)
{
for (k = 0; k < 10; step /= 2, k ++)
{
for (i = -1; i <= 1; i ++)
{
for (j =- 1; j <= 1; j ++)
{
double A1 = x0 + step * i;
double B1 = y0 + step * j;
if (Distance(x0,y0,x1,y1) + Distance(x0,y0,x2,y2) + Distance(x0,y0,x3,y3) > Distance(A1,B1,x1,y1) + Distance(A1,B1,x2,y2) + Distance(A1,B1,x3,y3))
{
x0=A1;y0=B1;
}
}
}
}
}
printf("该三角形的费马点为（%.3f ， %.3f）\n",x0,y0);
}
int main()
{
double x1, y1, x2, y2, x3, y3;
while (scanf("%lf%lf%lf%lf%lf%lf", &x1, &y1, &x2, &y2, &x3, &y3) != EOF)
{
printf("该三角形面积为： %.3f\n",area(x1,y1,x2,y2,x3,y3)/2) ;
waixin(x1,y1,x2,y2,x3,y3);
chuixin(x1,y1,x2,y2,x3,y3);
zhongxin(x1,y1,x2,y2,x3,y3);
neixin(x1,y1,x2,y2,x3,y3);
fermentpoint(x1,y1,x2,y2,x3,y3);
}
return 0;
}
7. 多边形的费马点（凸包+费马点）
给你一个多边形不是按顺序输入的点求其费马点（卡随机淬火算法）
代码：
#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <ctime>
using namespace std;
typedef struct pnode
{
double x,y,d;
}point;
point Pn,P[105];
double dist1( point a, point b )
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
double dist( point p, int N )
{
double sum = 0.0;
for ( int i = 0 ; i <= N ; ++ i )
sum += dist1( p, P[i] );
return sum;
}
double crossproduct( point a, point b, point c )
{
return (b.x-a.x)*(c.y-a.y)-(c.x-a.x)*(b.y-a.y);
}
bool cmp1( point a, point b )
{
if ( a.x == b.x ) return a.y < b.y;
else return a.x < b.x;
}
bool cmp2( point a, point b )
{
return crossproduct( P[0], a, b ) > 0;
}
bool cmp3( point a, point b )
{
double cp = crossproduct( P[0], a, b );
if ( cp == 0 ) {
if ( crossproduct( P[0], a, Pn ) == 0 )
return a.d > b.d;
else return a.d < b.d;
}else return cp > 0;
}
double Graham( int N )
{
sort( P+0, P+N, cmp1 );
sort( P+1, P+N, cmp2 );
for ( int i = 1 ; i < N ; ++ i )
P[i].d = dist1( P[0], P[i] );
Pn = P[N-1];
sort( P+1, P+N, cmp3 );
int top = N;
if ( N > 2 )
{
top = 2;
for ( int i = 3 ; i < N ; ++ i )
{
while ( crossproduct( P[top-1], P[top], P[i] ) < 0 ) -- top;
P[++ top] = P[i];
}
//删掉共线的点
int now = 1;
for ( int i = 2 ; i <= top ; ++ i )
{
if ( crossproduct( P[now-1], P[now], P[i] ) == 0 )
P[now] = P[i];
else P[++ now] = P[i];
}
top = now;
}
//随机增量逼近法
point t,s = P[0];
double sp = 10000.0,esp = 0.01;
double temp,min = dist( s, top );
while ( sp > esp )
{
int flag = 0;
for ( int i = 0 ; i < 10 ; ++ i )
{
t = s;
int R = rand()%361;
t.x += sp*cos(double (R));
t.y += sp*sin(double (R));
temp = dist( t, top );
if ( min > temp )
{
min = temp;
s = t;
flag = 1;
}
}
if ( !flag ) sp /= 2.0;
}
printf("%.3f %.3f\n",s.x,s.y);
return min;
}
int main()
{
srand( time(NULL));
int N;
while (~scanf("%d",&N))
{
for ( int i = 0 ; i < N ; ++ i )
scanf("%lf%lf",&P[i].x,&P[i].y);
printf("%.3f\n",Graham( N ));
}
return 0;
}
8. 凸包问题
a. 平面上给你一些让你求围成的最大凸包，
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
const int N = 10005;
const double eps = 1e-8;
struct gPoint
{
double x, y;
}p[N];
int n;
gPoint convex[N<<1];
int vCount;
int dbcmp(double x)
{
if(x < -eps)
return -1;
else if(x > eps)
return 1;
else
return 0;
}
bool cmp(const gPoint &a, const gPoint &b)
{
if(dbcmp(a.x-b.x) == 0)
return dbcmp(a.y-b.y)<=0;
else
return dbcmp(a.x-b.x)<=0;
}
double Dis(gPoint a, gPoint b)
{
double tx=b.x-a.x;
double ty=b.y-a.y;
return sqrt(tx*tx+ty*ty);
}
bool rCheck(const gPoint &a, const gPoint &b, const gPoint &c)
{
double x1 = b.x - a.x;
double y1 = b.y - a.y;
double x2 = c.x - a.x;
double y2 = c.y - a.y;
double tmp = x1*y2-x2*y1;
return dbcmp(tmp) < 0;
}
void grahamScan()
{
vCount = 2;
convex[0] = p[0];
if(n == 1) return;
convex[1] = p[1];
for(int i = 2; i < n; i++)
{
while(vCount>=2 && rCheck(convex[vCount-2],convex[vCount-1], p[i]))
vCount--;
convex[vCount++] = p[i];
}
for(int i = n-2; i >= 0; i--)
{
while(vCount>=2 && rCheck(convex[vCount-2],convex[vCount-1], p[i]))
vCount--;
convex[vCount++] = p[i];
}
}
int main()
{
scanf("%d", &n);
for(int i = 0; i < n; i++) scanf("%lf%lf", &p[i].x, &p[i].y);
sort(p, p+n, cmp);
grahamScan();
double ans = 0.0;
cout<<"凸包各边点为："<<endl;
for(int i = 1; i < vCount; i++)
{
cout<<convex[i].x<<" "<<convex[i].y<<endl;
ans+=Dis(convex[i-1], convex[i])+eps;
}
printf("凸包周长为 ： %.3f\n", ans);
return 0;
}
b. 求凸包中两点最大距离
#include <cmath>
#include <algorithm>
#include <iostream>
using namespace std;
#define MAXN 50005
struct Point
{
int x, y;
}pset[MAXN],ch[MAXN];
bool cmp(Point a,Point b)
{
return b.y<a.y||(a.y==b.y&&b.x<a.x);
}
int cross(Point a,Point b,Point c)
{
return (a.x - c.x) * (b.y - c.y) - (b.x - c.x) * (a.y - c.y);
}
void convex_hull(Point *p,Point *ch,int n,int &len)
{
sort(p, p+n,cmp);
ch[0]=p[0];
ch[1]=p[1];
int top=1;
for(int i=2;i<n;i++)
{
while(top>0&&cross(ch[top],p[i],ch[top-1])<=0)
top--;
ch[++top]=p[i];
}
int tmp=top;
for(int i=n-2;i>=0;i--)
{
while(top>tmp&&cross(ch[top],p[i],ch[top-1])<=0)
top--;
ch[++top]=p[i];
}
len=top;
}
int dist2(Point a,Point b)
{
return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
//前面的可以作为整数凸包模板
int rotating_calipers(Point *ch,int n)
{
int q=1,ans=0;
ch[n]=ch[0];
for(int p=0;p<n;p++)
{
while(cross(ch[p+1],ch[q+1],ch[p])>cross(ch[p+1],ch[q],ch[p]))
q=(q+1)%n;
ans=max(ans,max(dist2(ch[p],ch[q]),dist2(ch[p+1],ch[q+1])));
}
return ans;
}
int main()
{
int n, len;
while(scanf("%d", &n)!=EOF)
{
for(int i = 0;i < n;i++)
{
scanf("%d %d",&pset[i].x,&pset[i].y);
}
convex_hull(pset,ch,n,len);
for(int i=0;i<len;i++)
cout<<ch[i].x<<" "<<ch[i].y<<endl;
printf("%d\n",rotating_calipers(ch,len));
}
return 0;
}
c. 求两凸包点对距离最近
#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,Point C,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,int m)
{
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("%.5f\n",min(Solve(P,Q,n,m),Solve(Q,P,m,n)));
}
return 0;
}
9. 点集最小圆覆盖(随机增量法)
给你n个点，求一个最小的圆把所有点覆盖了，这些点可以在圆上；
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <cstdlib>
using namespace std;
#define N 505
struct POINT
{
double x, y;
} p[N];
int n;
inline double dist(const POINT &a, const POINT &b)
{
return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y));
}
POINT circumcenter(POINT &a, POINT &b, POINT &c)
{
POINT ret;
double a1=b.x-a.x, b1=b.y-a.y, c1=(a1*a1+b1*b1)/2;
double a2=c.x-a.x, b2=c.y-a.y, c2=(a2*a2+b2*b2)/2;
double d = a1*b2 - a2*b1;
ret.x = a.x + (c1*b2-c2*b1)/d;
ret.y = a.y + (a1*c2-a2*c1)/d;
return ret;
}
void solve()
{
random_shuffle(p, p+n); //随机化序列，std里面的随机函数
POINT c;
double r = 0;
for (int i=1; i<n; i++)
{
if (dist(p[i], c) <= r) continue;
c = p[i];
r = 0;
for (int j=0; j<i; j++)
{
if (dist(p[j], c) <= r) continue;
c.x = (p[i].x+p[j].x)/2;
c.y = (p[i].y+p[j].y)/2;
r = dist(p[j], c);
for (int k=0; k<j; k++)
{
if (dist(p[k], c) <= r) continue;
c = circumcenter(p[i], p[j], p[k]);
r = dist(p[i], c);
}
}
}
printf("%.2f %.2f %.2f\n", c.x, c.y, r);
}
int main()
{
while (scanf(" %d", &n) == 1 && n)
{
for (int i=0; i<n; i++)
scanf(" %lf %lf", &p[i].x, &p[i].y);
solve();
}
return 0;
}
10. 两圆相交的面积
#include<cmath>
#include<iomanip>
#include<algorithm>
#include <iostream>
using namespace std;
int main()
{
double d,t,t1,area,x,y,xx,yy,r,rr;
while(cin>>x>>y>>r)
{
cin>>xx>>yy>>rr;
d=sqrt((x-xx)*(x-xx)+(y-yy)*(y-yy));
if(d>=r+rr)
area=0;
else if(d<=fabs(r-rr))
area=min(acos(-1.0)*r*r,acos(-1.0)*rr*rr);
else
{
t=(r*r+d*d-rr*rr)/2.0/d;
t1=sqrt(r*r-t*t);
area=-d*t1+r*r*acos(t/r)+rr*rr*acos((d-t)/rr);
}
cout<<fixed<<setprecision(3)<<area<<endl;
}
return 0;
}
