计算几何细节梳理&模板 - Flash_Hu - 博客园 (cnblogs.com) 必看,还有一些知识点我漏掉了(细节)
1、二维几何
点和向量 point Vector(直线、线段、射线)
#include<bits/stdc++.h>
using namespace std;
const double pi=acos(-1.0);
const double eps=1e-8;
int sgn(double x){ //判断x是不是等于0
if(fabs(x)<eps) return 0;
else return x<0? -1:1;
}
int dcmp(double x,double y){ //比较两个浮点数,0为相等,-1为小于,1为大于
if(fabs(x-y)<eps) return 0;
else return x<y? -1:1;
}
//点,点用坐标(x,y)表示
struct point{
double x,y;
point(){}
point(double x,double y):x(x),y(y){}
point operator + (point B){return point(x+B.x,y+B.y);}
point operator - (point B){return point(x-B.x,y-B.y);}
point operator * (double k){return point(x*k,y*k);}
point operator / (double k){return point(x/k,y/k);}
bool operator == (point B){return sgn(x-B.x)==0&&sgn(y-B.y)==0;}
bool operator < (point B){
return sgn(x-B.x)<0||(sgn(x-B.x)==0&&sgn(y-B.y)<0);}
};
//两点之间的距离
double Distance(point a,point b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
//或者用库函数hypot()
return hypot(a.x-b.x,a.y-b.y);
}
点积和叉积
点积dot(a,b) 叉积cross(a,b)
//向量:有大小有方向,如果把其中一个点移到源点,那么就可以用点来表示了
//向量的运算:运算符重载,加减乘除
typedef point Vector;
//点积和叉积
//点积是A*B=|A||B|cosC C是A,B之间的夹角
//A*B=A.x*B.x+A.y*B.y
double dot(Vector A,Vector B) {
return A.x*B.x+A.y*B.y;
}
//应用
//1.判断A,B的夹角是钝角还是锐角
//如果dot(A,B)>0 ,那么夹角为锐角 如果 dot(A,B)<0 ,那么夹角为钝角 如果dot(A,B)=0 ,那么夹角为直角
//2.求向量A的长度
double len(Vector A){
return sqrt(dot(A,A));
//或者求长度的平方
}
double len2(Vector A){
return dot(A,A);
}
//3.求向量A,B的夹角大小
double angle(Vector A,Vector B){
return acos(dot(A,B)/len(A)/len(B));
}
//叉积是A*B=|A||B|sinC C表示A旋转到B的夹角
//A*B=A.x*B.y-A.y*B.x 有正负号,几何意义是A和B形成的平行四边形的“有向”面积
//计算过程中A,B是有顺序的,A*B与B*A是相反的 ------小窍门,左边的边向逆时针转到第二条边,就是正的
double cross(Vector A,Vector B){
return A.x*B.y-A.y*B.x;
}
//应用:
//1.判断方向关系:如果A*B>0 那么B在A的逆时针方向
// 如果A*B<0 那么B在A的顺时针方向
// 如果A*B=0 那么B与A共线,可能同向,也可能反向
//2。计算两向量构成的平行四边形的有向面积
//3个点A,B,C以A为公共点,得到两个向量B-A,C-A,构成的平行四边形面积为:
double area2(Vector A,Vector B,Vector C){
return cross(B-A,C-A);
}
//3.三个点构成的三角形面积=平行四边形的一半
//4.向量旋转
//使向量(x,y)绕起点逆时针旋转,角度为C,那么旋转后的向量为
//x`=xcosC-ysinC
//y`=xsinC+ycosC
Vector rotate(Vector A,double rad){ //rad是旋转角度
return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));
}
//如果逆时针旋转90,rotate(a,pi/2),返回vector(-a.y,a.x)
//如果顺时针旋转90,rotate(a,-pi/2),返回vector(a.y,-a.x)
//如果求单位法向量:逆时针转90然后取单位值
Vector normal(Vector A){
return Vector(-A.y/len(A),A.x/len(A));
}
//5.用叉积检查两个向量是不是平行或者重合
bool parallel(Vector A,Vector B){
return sgn(cross(A,B))==0;
}
点和线:
1.直线的表示、线段的表示
2、点和直线(位置关系、距离、投影、对称点) 3、点和线段(位置关系、距离)
4、两条直线(位置关系、交点) 5、两条线段(是否相交、交点)
//直线的表示:
//1.用直线上的两个点来表示
//2.ax+by+c=0
//3.y=kx+b
//4.P=P0+vt 点向式:比较适合计算机处理,表示直线(t无限制)、线段(t=0~1)、射线(t>0)
struct Line{
point p1,p2;
Line(){}
//两个点
Line(point p1,point p2):p1(p1),p2(p2){}
//一个点,一个倾斜角angle在0,pi之间
Line(point p,double angle){
p1=p;
if(sgn(angle-pi/2)==0) {p2=(p1+point(0,1));}
else{p2=(p1+point(1,tan(angle)));}
}
Line(double a,double b,double c){
if(sgn(a)==0){
p1=point(0,-c/b);
p2=point(1,-c/b);
}
if(sgn(b)==0){
p1=point(-c/a,0);
p2=point(-c/a,1);
}
else{
p1=point(0,-c/b);
p2=point(1,((-c-a)/b));
}
}
};
//点和直线的关系
int point_line_relation(point p,Line v){
int c=sgn(cross(p-v.p1,v.p2-v.p1));
if(c<0) return 1; //1 p在v的左边
if(c>0) return 2; //2 p在v的右边
return 0; //p在v上
}
//点和直线的距离 p v(p1,p2)
//用叉积求p,p1,p2构成的平行四边形面积,然后除以平行四边形的底边长,就得到了高
double dis_point_line(point p,Line v){
return fabs(cross(p-v.p1,v.p2-v.p1))/Distance(v.p1,v.p2);
}
//点在直线上的投影
point point_line_proj(point p,Line v){
double k=dot(v.p2-v.p1,p-v.p1)/len2(v.p2-v.p1);
return v.p1+(v.p2-v.p1)*k;
}
//点关于直线的对称点
point point_line_symmetry(point p,Line v){
point q=point_line_proj(p,v);
return point(2*q.x-p.x,2*q.y-p.y);
}
//两条直线的关系
int line_relation(Line v1,Line v2){
if(sgn(cross(v1.p2-v1.p1,v2.p2-v2.p1))==0){
if(point_line_relation(v1.p1,v2)==0) return 1; //重合
else return 0; //平行
}
return 2; //相交
} x
//两条直线的交点
point cross_point(point a,point b,point c,point d){
double s1=cross(b-a,c-a);
double s2=cross(b-a,d-a);
//在调用这个函数之前要确保s1 s2不相等
return point(c.x*s2-d.x*s1,c.y*s2-d.y*s1)/(s2-s1);
}
//线段的表示:直接用直线的结构
typedef Line Segment;
//点和线段的关系
//先用叉积看是不是共线,然后点积看p与v的两个端点是不是形成的钝角
bool point_on_seg(point p,Line v){
return sgn(cross(p-v.p1,v.p2-v.p1))==0&&sgn(dot(p-v.p1,p-v.p2))<=0;
}
//点到线段的距离:从p出发对A,B做垂线,如果就在AB上,那就是最小值,否则就是到A或者到B的最小值
double dis_point_seg(point p,Segment v){
if(sgn(dot(p-v.p1,v.p2-v.p1))<0||sgn(dot(p-v.p2,v.p1-v.p2))<0) return min(Distance(p,v.p1),Distance(p,v.p2));
return dis_point_line(point(p,v));
}
//判断两个线段是否相交
//如果一条线段的两端在另一条线段的两侧,那么两个端点与另一线段产生的两个叉积正负相反
bool cross_segment(point a,point b,point c,point d){
double c1=cross(b-a,c-a),c2=cross(b-a,d-a);
double d1=cross(d-c,a-c),d2=cross(d-c,b-c);
return sgn(c1)*sgn(c2)<0&&sgn(d1)*sgn(d2)<0; //1:相交 0:不想交
}
//两个线段的交点:先判断两条线段是否相交,如果相交转化为两条直线求交点
多边形
1.点在三角形的内部还是外部 2.求多边形的面积 3.求多边形的重心
//多边形
//1.判断点在多边形的内部:转角法,以P为起点引起一条水平线,检查与多边形每条边的相交情况,统计P穿过这些边的次数
//c=cross(p-j,i-j) u=i.y-p.y v=j.y-p.y
int point_in_polygon(point pt,point *p,int n){
for(int i=0;i<n;i++){
if(p[i]==pt) return 3; //3.点在多边形的顶点上
}
for(int i=0;i<n;i++){
Line v=Line(p[i],p[(i+1)%n]);
if(point_on_seg(pt,v)) return 2; //2.点在多边形的边上
}
int num=0;
for(int i=0;i<n;i++){
int j=(i+1)%n;
int c=sgn(cross(pt-p[j],p[i]-p[j]));
int u=sgn(p[i].y-pt.y);
int v=sgn(p[j].y-pt.y);
if(c>0&&u<0&&v>=0) num++;
if(c<0&&u>=0&&v<0) num--;
}
return num!=0; //1:点在内部 0:点在外部
}
//求多边形的面积
//以原点为中心点划分三角形,然后求多边形的面积
double polygon_area(point *p,int n){
double area=0;
for(int i=0;i<n;i++) area+=cross(p[i],p[(i+1)%n]);
return area/2;
}
//求多边形的重心
//将多边形进行三角剖分,算出每个三角形的重心,三角形的重心是3点坐标的平均值,然后对每个三角形的邮箱面积加权平均???
point polygon_center(point *p,int n){
point ans(0,0);
if(polygon_area(p,n)==0) return ans;
for(int i=0;i<n;i++){
ans=ans+(p[i]+p[(i+1)%n])*cross(p[i],p[(i+1)%n]);
}
return ans/polygon_area(p,n)/6;
}
hdu 1115 求N(3<N<1000000)边形的重心
#include<iostream>
#include<cmath>
using namespace std;
const double pi=acos(-1.0);
const double eps=1e-8;
int sgn(double x){ //判断x是不是等于0
if(fabs(x)<eps) return 0;
else return x<0? -1:1;
}
struct point{
double x,y;
point(){}
point(double x,double y):x(x),y(y){}
point operator + (point B){return point(x+B.x,y+B.y);}
point operator - (point B){return point(x-B.x,y-B.y);}
point operator * (double k){return point(x*k,y*k);}
point operator / (double k){return point(x/k,y/k);}
bool operator == (point B){return sgn(x-B.x)==0&&sgn(y-B.y)==0;}
bool operator < (point B){
return sgn(x-B.x)<0||(sgn(x-B.x)==0&&sgn(y-B.y)<0);}
};
typedef point Vector;
//点积和叉积
//点积是A*B=|A||B|cosC C是A,B之间的夹角
//A*B=A.x*B.x+A.y*B.y
double dot(Vector A,Vector B) {
return A.x*B.x+A.y*B.y;
}
double cross(Vector A,Vector B){
return A.x*B.y-A.y*B.x;
}
//求多边形的面积
//以原点为中心点划分三角形,然后求多边形的面积
double polygon_area(point *p,int n){
double area=0;
for(int i=0;i<n;i++) area+=cross(p[i],p[(i+1)%n]);
return area/2;
}
point polygon_center(point *p,int n){
point ans(0,0);
if(polygon_area(p,n)==0) return ans;
for(int i=0;i<n;i++){
ans=ans+(p[i]+p[(i+1)%n])*cross(p[i],p[(i+1)%n]);
}
return ans/polygon_area(p,n)/6;
}
int main(){
int t,n,i;
point center;
point p[100000];
scanf("%d",&t);
while(t--){
scanf("%d",&n);
for(int i=0;i<n;i++) scanf("%lf %lf",&p[i].x,&p[i].y);
center=polygon_center(p,n);
printf("%.2f %.2f\n",center.x,center.y);
}
return 0;
}
凸包
搞懂概念:给定一些点,求所有能把这些点包含在内的面积最小的多边形
主要是用点求凸包的算法:Graham算法O(nlong2n)、jarvis步进法O(nh),h是凸包上的顶点数
andrew算法(算法复杂度:nlogn):从做左下扫描到右上得到下凸包,从右到左扫描得到上凸包
hdu 1392 求凸包的周长
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<stack>
#include<cstdio>
#include<queue>
#include<map>
#include<vector>
#include<set>
using namespace std;
//hdu 1392 求凸包的周长
const double pi=acos(-1.0);
const double eps=1e-8;
const int maxn=104;
int sgn(double x){ //判断x是不是等于0
if(fabs(x)<eps) return 0;
else return x<0? -1:1;
}
struct point{
double x,y;
point(){}
point(double x,double y):x(x),y(y){}
point operator + (point B){return point(x+B.x,y+B.y);}
point operator - (point B){return point(x-B.x,y-B.y);}
point operator * (double k){return point(x*k,y*k);}
point operator / (double k){return point(x/k,y/k);}
bool operator == (point B){return sgn(x-B.x)==0&&sgn(y-B.y)==0;}
bool operator < (point B){
return sgn(x-B.x)<0||(sgn(x-B.x)==0&&sgn(y-B.y)<0);}
};
typedef point Vector;
//点积和叉积
//点积是A*B=|A||B|cosC C是A,B之间的夹角
//A*B=A.x*B.x+A.y*B.y
double dot(Vector A,Vector B) {
return A.x*B.x+A.y*B.y;
}
double cross(Vector A,Vector B){
return A.x*B.y-A.y*B.x;
}
double Distance(point a,point b){
return hypot(a.x-b.x,a.y-b.y);
}
int convex_hull(point *p,int n,point *ch){
sort(p,p+n);
n=unique(p,p+n)-p;
int v=0;
//求下凸包,如果p[i]是右拐弯的,那么这个点不在凸包上, 往回退
for(int i=0;i<n;i++){
while(v>1&&sgn(cross(ch[v-1]-ch[v-2],p[i]-ch[v-2]))<=0) v--;
ch[v++]=p[i];
}
int j=v;
//求上凸包
for(int i=n-2;i>=0;i--){
while(v>j&&sgn(cross(ch[v-1]-ch[v-2],p[i]-ch[v-2]))<=0) v--;
ch[v++]=p[i];
}
if(n>1) v--;
return v;
}
int main(){
int n;
point p[maxn],ch[maxn];
while(scanf("%d",&n)&&n){
for(int i=0;i<n;i++) scanf("%lf %lf",&p[i].x,&p[i].y);
int v=convex_hull(p,n,ch);
double ans=0;
if(v==1) ans=0;
else if(v==2) ans=Distance(ch[0],ch[1]);
else{
for(int i=0;i<v;i++) ans+=Distance(ch[i],ch[(i+1)%v]);
}
printf("%.2f\n",ans);
}
return 0;
}
Graham扫描法:
凸包--Graham扫描法_小蒟蒻yyb的博客-CSDN博客
首先找到最靠近左下的那个点,这个点一定在凸包上(不难理解吧。。。画个图就知道了)
以这个点为极点,其他点按照极角排序
然后按照顺序依次访问所有点,判断可行性
struct Node
{
int x,y;
}p[MAX],S[MAX];//p储存节点的位置,S是凸包的栈
inline bool cmp(Node a,Node b)//比较函数,对点的极角进行排序
{
double A=atan2((a.y-p[1].y),(a.x-p[1].x));
double B=atan2((b.y-p[1].y),(b.x-p[1].x));
if(A!=B)return A<B;
else return a.x<b.x; //这里注意一下,如果极角相同,优先放x坐标更小的点
}
long long Cross(Node a,Node b,Node c)//计算叉积
{
return 1LL*(b.x-a.x)*(c.y-a.y)-1LL*(b.y-a.y)*(c.x-a.x);
}
void Get()//求出凸包
{
p[0]=(Node){INF,INF};int k;
for(int i=1;i<=n;++i)//找到最靠近左下的点
if(p[0].y>p[i].y||(p[0].y==p[i].y&&p[i].x<p[0].x))
{p[0]=p[i];k=i;}
swap(p[k],p[1]);
sort(&p[2],&p[n+1],cmp);//对于剩余点按照极角进行排序
S[0]=p[1],S[1]=p[2];top=1;//提前在栈中放入节点
for(int i=3;i<=n;)//枚举其他节点
{
if(top&&Cross(S[top-1],p[i],S[top])>=0)
top--;//如果当前栈顶不是凸包上的节点则弹出
else S[++top]=p[i++];//加入凸包的栈中
}
//底下这个玩意用来输出凸包上点的坐标
//for(int i=0;i<=top;++i)
// printf("(%d,%d)\n",S[i].x,S[i].y);
}
最近点对
分治法O(nlog2n)
合并的时候,如果两个点在两个集合里,就不能简单合并,如果集合s1中最小距离是d1,集合s2中最小距离是d2,那么在两个集合中间的点找距离他们
小于d1,d2的,记录在tmp_p[]中,但不能用暴力法直接列出点集tmp_p[]中的所有点对,否则会超时,按照y坐标对tmp_p[]中的点排序,再加上剪枝
hdu 1007 求最近点对的一半
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<stack>
#include<cstdio>
#include<queue>
#include<map>
#include<vector>
#include<set>
//最近点对
using namespace std;
const double pi=acos(-1.0);
const double eps=1e-8;
const int maxn=100010;
const double INF=1e20;
int sgn(double x){ //判断x是不是等于0
if(fabs(x)<eps) return 0;
else return x<0? -1:1;
}
struct point{
double x,y;
};
double Distance(point a,point b){
return hypot(a.x-b.x,a.y-b.y);
}
bool cmpxy(point a,point b){
return sgn(a.x-b.x)<0||sgn(a.x-b.x)==0&&sgn(a.y-b.y)<0;
}
bool cmpy(point a,point b){
return sgn(a.y-b.y)<0;
}
point p[maxn],tmp_p[maxn];
double closest_pair(int left,int right){
double dis=INF;
if(left==right) return dis;
if(left+1==right) return Distance(p[left],p[right]);
int mid=(left+right)/2;
double d1=closest_pair(left,mid);
double d2=closest_pair(mid+1,right);
dis=min(d1,d2);
int k=0;
for(int i=left;i<=right;i++){
if(fabs(p[mid].x-p[i].x)<=dis){
tmp_p[k++]=p[i];
}
}
sort(tmp_p,tmp_p+k,cmpy);
for(int i=0;i<k;i++){
for(int j=i+1;j<k;j++){
if(tmp_p[j].y-tmp_p[i].y>=dis) break; //剪枝
dis=min(dis,Distance(tmp_p[i],tmp_p[j]));
}
}
return dis;
}
int main(){
int n;
while(scanf("%d",&n)&&n){
for(int i=0;i<n;i++) scanf("%lf %lf",&p[i].x,&p[i].y);
sort(p,p+n,cmpxy);
printf("%.2f\n",closest_pair(0,n-1)/2);
}
return 0;
}
旋转卡壳
是一种思想,可以用来求出凸包上最远点对距离
凸包的直径——旋转卡壳 - 小蒟蒻yyb - 博客园 (cnblogs.com)
poj会有编译错误,为啥
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<stack>
#include<cstdio>
#include<queue>
#include<map>
#include<vector>
#include<set>
using namespace std;
const int maxn=50010;
const int INF=0x3f3f3f3f;
struct point{
int x,y;
}p[maxn],p0,s[maxn];
int n,top,T;
bool cmp(point a,point b){
double aa=atan2(a.y-p0.y,a.x-p0.x);
double bb=atan2(b.y-p0.y,b.x-p0.x);
if(aa!=bb) return aa<bb;
else return a.x<b.x;
}
long long cross(int x1,int y1,int x2,int y2){
return (1LL*x1*y2-1LL*x2*y1);
}
long long cross2(point a,point b,point c){ //以a为中心点
return cross((b.x-a.x),(b.y-a.y),(c.x-a.x),(c.y-a.y));
}
void convex_hull(){ //寻找凸包
p0=(point){INF,INF};
int k=0;
for(int i=0;i<n;i++){ //寻找最下方的点
if(p0.y>p[i].y||(p0.y==p[i].y)&&(p0.x>p[i].x))
p0=p[i],k=i;
}
swap(p[k],p[0]);
sort(&p[1],&p[n],cmp); //按照极角对下方的点进行排序
s[0]=p[0];
s[1]=p[1];
top=1; //栈顶
for(int i=2;i<n;){ //求出凸包
if(top&&cross2(s[top-1],p[i],s[top])>=0) top--;
else s[++top]=p[i++]; //s里面放的是凸包边上的点
}
}
long long distan(point a,point b){
return 1LL*(a.x-b.x)*(a.x-b.x)+1LL*(a.y-b.y)*(a.y-b.y);
}
long long getmax(){ //求出直径
long long re=0;
if(top==1) //仅有两个点
return distan(s[0],s[1]);
s[++top]=s[0]; //把第一个点放到最后
int j=2;
for(int i=0;i<top;i++){ //枚举边
while(cross2(s[i],s[i+1],s[j])<cross2(s[i],s[i+1],s[j+1])) j=(j+1)%top;
re=max(re,max(distan(s[i],s[j]),distan(s[i+1],s[j])));
}
return re;
}
int main(){
scanf("%d",&n);
for(int i=0;i<n;i++) scanf("%lld %lld",&p[i].x,&p[i].y);
long long ans=INF,ss;
convex_hull();
ans=getmax();
printf("%lld\n",ans);
return 0;
}
半平面交
给出一堆半平面,然后求出凸包先
算法:增量法O(n^2)
这个算法:O(nlong2n)
主要是增加新板平面的时候根据情况进行队列的增删改查
应用:hdu 2297 "RUN"(真看不出来是半平面,没学过的话)
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<stack>
#include<cstdio>
#include<queue>
#include<map>
#include<vector>
#include<set>
using namespace std;
const int maxn=50010;
const double INF=1e12;
const double pi=acos(-1.0);
const double eps=1e-8;
int sgn(double x){
if(fabs(x)<eps) return 0;
else return x<0? -1:1;
}
struct point{
double x,y;
point(){}
point(double x,double y):x(x),y(y){}
point operator + (point B){return point(x+B.x,y+B.y);}
point operator - (point B){return point(x-B.x,y-B.y);}
point operator * (double k){return point(x*k,y*k);}
};
typedef point Vector;
double cross(Vector A,Vector B){
return A.x*B.y-A.y*B.x;
}
struct Line{
point p;
Vector v;
double ang;
Line(){}
Line(point p,Vector v):p(p),v(v){ang=atan2(v.y,v.x);}
bool operator < (Line &L){ //用于极角排序
return ang<L.ang;
}
};
bool onleft(Line L,point p){
return sgn(cross(L.v,p-L.p))>0; //点p在L的左边:即点p在L的外面
}
point cross_point(Line a,Line b){ //两条直线的交点
Vector u=a.p-b.p;
double t=cross(b.v,u)/cross(a.v,b.v);
return a.p+a.v*t;
}
vector<point> HPI(vector<Line> L){ //求半平面交,返回凸多边形
int n=L.size();
sort(L.begin(),L.end()); //将所有半平面按照极角排序
int first,last;
vector<point> p(n); //两个相交半平面的交点
vector<Line> q(n); //双端队列
vector<point> ans; //半平面形成的凸包
q[first=last=0]=L[0];
for(int i=0;i<n;i++){
//情况1:删除尾部的半平面
while(first<last&&!onleft(L[i],p[last-1])) last--;
//情况2:删除首部的半平面
while(first<last&&!onleft(L[i],p[first])) first++;
q[++last]=L[i]; //加入到双端队列的尾部
//极角相同的两个半平面保留左边
if(fabs(cross(q[last].v,q[last-1].v))<eps){
last--;
if(onleft(q[last],L[i].p)) q[last]=L[i];
}
//计算队列尾部的半平面交点
if(first<last) p[last-1]=cross_point(q[last-1],q[last]);
}
//情况3:删除队列尾部的无用半平面:尾部的交点在第一条线段外面
while(first<last&&!onleft(q[first],p[last-1])) last--;
if(last-first<=1) return ans;
p[last]=cross_point(q[last],q[first]); //计算队列尾首部的交点
for(int i=first;i<=last;i++) ans.push_back(p[i]);
return ans; //返回凸多边形的交点
}
int main(){
int t,n;
cin>>t;
while(t--){
cin>>n;
vector<Line> L;
L.push_back(Line(point(0,0),Vector(0,-1)));
L.push_back(Line(point(0,INF),Vector(-1,0))); //添加反向y轴,y极大的向左的直线
while(n--){
double a,b;
scanf("%lf %lf",&a,&b);
L.push_back(Line(point(0,a),Vector(1,b)));
}
vector<point> ans=HPI(L);
printf("%d\n",ans.size()-2); //去掉认为加上的两个点
}
return 0;
}
圆
圆的定义、点和圆的关系、直线和圆的关系、线段和圆的关系、直线和圆的交点(会用到前面的点、线知识)
#include<bits/stdc++.h>
using namespace std;
const double pi=acos(-1.0);
const double eps=1e-8;
//圆:圆的定义、直线和圆的关系、线段和圆的关系、直线和圆的交点(会用到前面的点、线知识)
int sgn(double x){
if(fabs(x)<eps) return 0;
else return x<0? -1:1;
}
struct point{
double x,y;
point(){}
point(double x,double y):x(x),y(y){}
point operator + (point B){return point(x+B.x,y+B.y);}
point operator - (point B){return point(x-B.x,y-B.y);}
point operator * (double k){return point(x*k,y*k);}
};
double Distance(point a,point b){
//return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
//或者用库函数hypot()
return hypot(a.x-b.x,a.y-b.y);
}
struct Line{
point p1,p2;
Line(){}
//两个点
Line(point p1,point p2):p1(p1),p2(p2){}
//一个点,一个倾斜角angle在0,pi之间
Line(point p,double angle){
p1=p;
if(sgn(angle-pi/2)==0) {p2=(p1+point(0,1));}
else{p2=(p1+point(1,tan(angle)));}
}
Line(double a,double b,double c){
if(sgn(a)==0){
p1=point(0,-c/b);
p2=point(1,-c/b);
}
if(sgn(b)==0){
p1=point(-c/a,0);
p2=point(-c/a,1);
}
else{
p1=point(0,-c/b);
p2=point(1,((-c-a)/b));
}
}
};
double dis_point_line(point p,Line v){
return fabs(cross(p-v.p1,v.p2-v.p1))/Distance(v.p1,v.p2);
}
//点在直线上的投影
point point_line_proj(point p,Line v){
double k=dot(v.p2-v.p1,p-v.p1)/len2(v.p2-v.p1);
return v.p1+(v.p2-v.p1)*k;
}
//点关于直线的对称点
point point_line_symmetry(point p,Line v){
point q=point_line_proj(p,v);
return point(2*q.x-p.x,2*q.y-p.y);
}
//线段的表示:直接用直线的结构
typedef Line Segment;
//点和线段的关系
//先用叉积看是不是共线,然后点积看p与v的两个端点是不是形成的钝角
bool point_on_seg(point p,Line v){
return sgn(cross(p-v.p1,v.p2-v.p1))==0&&sgn(dot(p-v.p1,p-v.p2))<=0;
}
//点到线段的距离:从p出发对A,B做垂线,如果就在AB上,那就是最小值,否则就是到A或者到B的最小值
double dis_point_seg(point p,Segment v){
if(sgn(dot(p-v.p1,v.p2-v.p1))<0||sgn(dot(p-v.p2,v.p1-v.p2))<0) return min(Distance(p,v.p1),Distance(p,v.p2));
return dis_point_line(point(p,v));
}
//圆的定义
struct circle{
point c;
double r;
circle(){}
circle(point c,double r):c(r),r(r){}
circle(double a,double b,double _r){
c=point(a,b);r=_r;
}
};
//点和圆的关系:根据点到圆心的距离判断
int point_circle_relation(point p,circle c){
double dst=Distance(p,c.c);
if(sgn(dst-c.r)<0) return 0; //在园内
if(sgn(dst-c.r)==0) return 1; //在圆上
else return 2; //点在圆外
}
//直线和圆的关系:根据圆心到直线的距离判断
int line_circle_relation(Line v,circle c){
double dst=dis_point_line(c.c,v);
if(sgn(dst-c.r)<0) return 0; //直线和圆相交
if(sgn(dst-c.r)==0) return 1; //直线和圆相切
return 2; //在圆外
}
//线段和圆的关系
int seg_circle_relation(Segment v,circle c){
double dst=dis_point_seg(c.c,v);
if(sgn(dst-c.r)<0) return 0; //线段在圆内
if(sgn(dst-c.r)==0) return 1; //线段和圆相切
return 2; //线段在圆外
}
//直线和圆的交点 pa和pb是交点,返回的是交点的个数
int line_cross_circle(Line v,circle c,point &pa,point &pb){
if(line_circle_relation(v,c)==2) return 0 ; //无交点
point q=point_line_proj(c.c,v); //圆心在直线上的投影点
double d=dis_point_line(c.c,v); //圆心到直线的距离
double k=sqrt(c.r*c.r-d*d);
if(sgn(k)==0){
pa=q;pb=q;return 1; //相切的情况
}
point n=(v.p2-v.p1)/len(v.p2-v.p1); //单位向量
pa=q+n*k; pb=q-n*k;
return 2; //两个交点
}
圆的模板的应用:hdu 5572
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<stack>
#include<cstdio>
#include<queue>
#include<map>
#include<vector>
#include<set>
//hdu 5572 圆的模板应用
using namespace std;
const double pi=acos(-1.0);
const double eps=1e-8;
int sgn(double x){ //判断x是不是等于0
if(fabs(x)<eps) return 0;
else return x<0? -1:1;
}
struct point{
double x,y;
point(){}
point(double x,double y):x(x),y(y){}
point operator + (point B){return point(x+B.x,y+B.y);}
point operator - (point B){return point(x-B.x,y-B.y);}
point operator * (double k){return point(x*k,y*k);}
point operator / (double k){return point(x/k,y/k);}
};
typedef point Vector;
double dot(Vector A,Vector B) {
return A.x*B.x+A.y*B.y;
}
double len(Vector A){
return sqrt(dot(A,A));
//或者求长度的平方
}
double len2(Vector A){
return dot(A,A);
}
double cross(Vector A,Vector B){
return A.x*B.y-A.y*B.x;
}
//两点之间的距离
double Distance(point a,point b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
//或者用库函数hypot()
//return hypot(a.x-b.x,a.y-b.y);
}
struct Line{
point p1,p2;
Line(){}
//两个点
Line(point p1,point p2):p1(p1),p2(p2){}
};
typedef Line Segment;
//点和直线的关系
int point_line_relation(point p,Line v){
int c=sgn(cross(p-v.p1,v.p2-v.p1));
if(c<0) return 1; //1 p在v的左边
if(c>0) return 2; //2 p在v的右边
return 0; //p在v上
}
//点和直线的距离 p v(p1,p2)
//用叉积求p,p1,p2构成的平行四边形面积,然后除以平行四边形的底边长,就得到了高
double dis_point_line(point p,Line v){
return fabs(cross(p-v.p1,v.p2-v.p1))/Distance(v.p1,v.p2);
}
//点到线段的距离:从p出发对A,B做垂线,如果就在AB上,那就是最小值,否则就是到A或者到B的最小值
double dis_point_seg(point p,Segment v){
if(sgn(dot(p-v.p1,v.p2-v.p1))<0||sgn(dot(p-v.p2,v.p1-v.p2))<0) return min(Distance(p,v.p1),Distance(p,v.p2));
return dis_point_line(p,v);
}
//点在直线上的投影
point point_line_proj(point p,Line v){
double k=dot(v.p2-v.p1,p-v.p1)/len2(v.p2-v.p1);
return v.p1+(v.p2-v.p1)*k;
}
//点关于直线的对称点
point point_line_symmetry(point p,Line v){
point q=point_line_proj(p,v);
return point(2*q.x-p.x,2*q.y-p.y);
}
//圆的定义
struct circle{
point c;
double r;
circle(){}
circle(point c,double r):c(c),r(r){}
circle(double a,double b,double _r){
c=point(a,b);r=_r;
}
};
//直线和圆的关系:根据圆心到直线的距离判断
int line_circle_relation(Line v,circle c){
double dst=dis_point_line(c.c,v);
if(sgn(dst-c.r)<0) return 0; //直线和圆相交
if(sgn(dst-c.r)==0) return 1; //直线和圆相切
return 2; //在圆外
}
//线段和圆的关系
int seg_circle_relation(Segment v,circle c){
double dst=dis_point_seg(c.c,v);
if(sgn(dst-c.r)<0) return 0; //线段在圆内
if(sgn(dst-c.r)==0) return 1; //线段和圆相切
return 2; //线段在圆外
}
//直线和圆的交点 pa和pb是交点,返回的是交点的个数
int line_cross_circle(Line v,circle c,point &pa,point &pb){
if(line_circle_relation(v,c)==2) return 0 ; //无交点
point q=point_line_proj(c.c,v); //圆心在直线上的投影点
double d=dis_point_line(c.c,v); //圆心到直线的距离
double k=sqrt(c.r*c.r-d*d);
if(sgn(k)==0){
pa=q;pb=q;return 1; //相切的情况
}
point n=(v.p2-v.p1)/len(v.p2-v.p1); //单位向量
pa=q+n*k; pb=q-n*k;
return 2; //两个交点
}
int main(){
int t;
scanf("%d",&t);
for(int cas=1;cas<=t;cas++){
circle o;
point a,b,v;
scanf("%lf %lf %lf",&o.c.x,&o.c.y,&o.r);
scanf("%lf %lf %lf %lf",&a.x,&a.y,&v.x,&v.y);
scanf("%lf %lf",&b.x,&b.y);
Line l(a,a+v);
Line t(a,b);
//情况1:直线和圆不想交,而且直线经过点
if(point_line_relation(b,l)==0&&seg_circle_relation(t,o)>=1&&sgn(cross(b-a,v))==0) printf("Case #%d: Yes\n",cas);
else{
point pa,pb;
//情况2:直线和圆相切,不经过点
if(line_cross_circle(l,o,pa,pb)!=2) printf("Case #%d: No\n",cas);
else{
//情况3:直线和圆相交
point cut; //直线和圆的碰撞点
if(Distance(pa,a)>Distance(pb,a)) cut=pb;
else cut=pa;
Line mid(cut,o.c); //圆心到碰撞点的线
point en=point_line_symmetry(a,mid); //镜像点
Line light(cut,en); //反射线
if(Distance(light.p2,b)>Distance(light.p1,b)) swap(light.p1,light.p2);
if(sgn(cross(light.p2-light.p1,point(b.x-cut.x,b.y-cut.y)))==0) printf("Case #%d: Yes\n",cas);
else printf("Case #%d: No\n",cas) ;
}
}
}
return 0;
}
最小圆覆盖:就是找个点能够覆盖平面上所有的点
两种方法:几何算法(增量法)和模拟退火(很慢)
hdu 3007 裸题
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<stack>
#include<cstdio>
#include<queue>
#include<map>
#include<vector>
#include<set>
//hdu 5572 圆的模板应用
using namespace std;
const double pi=acos(-1.0);
const int maxn=505;
const double eps=1e-8;
int sgn(double x){ //判断x是不是等于0
if(fabs(x)<eps) return 0;
else return x<0? -1:1;
}
struct point{
double x,y;
point(){}
point(double x,double y):x(x),y(y){}
point operator + (point B){return point(x+B.x,y+B.y);}
point operator - (point B){return point(x-B.x,y-B.y);}
point operator * (double k){return point(x*k,y*k);}
point operator / (double k){return point(x/k,y/k);}
};
//两点之间的距离
double Distance(point a,point b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
//或者用库函数hypot()
//return hypot(a.x-b.x,a.y-b.y);
}
//求三角形abc的外接圆的圆心
point circle_center(const point a,const point b,const point c){
point center;
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;
center.x=a.x+(c1*b2-c2*b1)/d;
center.y=a.y+(a1*c2-a2*c1)/d;
return center;
}
//求最小圆覆盖,返回圆心c,半径r
void min_cover_circle(point *p,int n,point &c,double &r){
//随机函数,打乱所有点,很重要,降低算法复杂度
random_shuffle(p,p+n);
c=p[0];r=0; //初始值
for(int i=1;i<n;i++){
if(sgn(Distance(p[i],c)-r)>0){ //点i在圆的外部,需要拓展
c=p[i];r=0; //这里就直接更新圆(相当于直接进队首)
for(int j=0;j<i;j++){ //重新检查前面所有的点
if(sgn(Distance(p[j],c)-r)>0) { //两点定圆
c.x=(p[j].x+p[i].x)/2;
c.y=(p[j].y+p[i].y)/2;
r=Distance(p[j],c);
for(int k=0;k<j;k++){ //两点不能定圆那就三点定圆
if(sgn(Distance(p[k],c)-r)>0){
c=circle_center(p[i],p[j],p[k]);
r=Distance(p[i],c);
}
}
}
}
}
}
}
int main(){
int n;
point p[maxn];
point c;double r;
while(~scanf("%d",&n)&&n){
for(int i=0;i<n;i++) scanf("%lf %lf",&p[i].x,&p[i].y);
min_cover_circle(p,n,c,r);
printf("%.2f %.2f %.2f\n",c.x,c.y,r);
}
return 0;
}
模拟退火
//模拟退火算法 :算法复杂度很高
void min_cover_circle2(point *p,int n,point &c,double &r){
double T=100.0; //初识温度
double delta=0.98 //降温系数
c=p[0];
int pos;
while(T>eps){ //eps是终止温度
pos=0;r=0; //初识:p[0]是圆心,半径是0
for(int i=0;i<=n-1;i++){ //找距离圆心最远的点
if(Distance(c,p[i])>r){
r=Distance(c,p[i]); //距圆心最远的点肯定在圆周上
pos=i;
}
}
c.x+=(p[pos].x-c.x)/r*T; //逼近最后的解
c.y+=(p[pos].y-c.y)/r*T;
T*=delta;
}
}
三维几何
除了一些基础的点、向量、直线、线段的表示以外,还有三维点积、叉积、三角形的面积、平面(3个点)、平面法向量、直线和平面的关系
四面体的有向体积
最小球覆盖、三维凸包(懒得写了)
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<stack>
#include<cstdio>
#include<queue>
#include<map>
#include<vector>
#include<set>
//三维几何
int sgn(double x){ //判断x是不是等于0
if(fabs(x)<eps) return 0;
else return x<0? -1:1;
}
struct point3{
double x,y,z;
point3(){}
point3(double x,double y,double z) x(x),y(y),z(z){}
point3 operator + (point3 b){return point3(x+b.x,y+b.y,z+b.z);}
point3 operator - (point3 b){return point3(x-b.x,y-b.y,z-b.z);}
point3 operator * (double b){return point3(x*b,y*b,z*b);}
point3 operator / (point3 b){return point3(x/k,y/k,z/k);}
bool operator == (point3 b){
return sgn(x-b.x)==0&&sgn(y-b.y)==0&&sgn(z-b.z)==0;
}
};
typedef point3 Vector3; //三维向量
double Distance(Vector3 a,Vector3 b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+(a.z-b.z)*(a.z-b.z));
}
//线和线段
struct Line3{
point3 p1,p2;
Line3(){}
Line3(point3 p1,point3 p2):p1(p1),p2(p2){}
};
typedef Line3 Segment3;
//三维点积
double dot3(Vector3 a,Vector3 b){
return a.x*b.x+a.y*b.y+a.z*b.z;
}
//求向量的长度
double len3(Vector3 a){
return sqrt(dot3(a,a));
}
//或者是求长度的平方
double len32(Vector3 a){
return dot3(a,a);
}
//求向量A,B的夹角大小
double angle3(Vector3 a,Vector3 b){
return acos(dot3(a,b)/len3(a,a)/len3(b,b));
}
//三维叉积:是一个向量,可以把三维叉积A,B的叉积看成垂直于A和B的向量
Vector3 cross3(Vector3 a,Vector3 b){
return point3(a.y*b.z-a.z-b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y-b.x);
}
//三角形面积:有向面积,先求三维叉积,然后取叉积的长度值
double area32(point3 a,point3 b,point3 c){
return len3(cross3(b-a,c-a));
}
// 判断点p在不在三角形ABC内:在内的话进行三角剖分后的面积相加为原本ABC的面积
dcmp(area32(p,a,b)+area32(p,b,c)+area32(p,a,c),area32(a,b,c))==0
//平面:用三个点可以确定一个平面
struct plane{
point3 p1,p2,p3;
plane(){}
plane(point3 p1,point3 p2,point3 p3):p1(p1),p2(p2),p3(p3){}
};
//平面法向量
point3 pvec(point3 a,point3 b,point3 c){
return cross3(b-a,c-a);
}
//或者这样也行
point3 pvec(plane x){
return cross(x.p2-x.p1,x-p3-x.p1);
}
//直线和平面的关系
int line_corss_plane(Line3 u,plane f,point3 &p){
point3 v=pvec(f);
double x=dot3(v,u.p2-f.p1);
double y=dot3(v,u.p1-f.p1);
double d=x-y;
if(sgn(x)==0&&sgn(y)==0) return -1; //直线在平面上
if(sgn(d)==0) return 0; //平行
p=((u.p1*x)-(u.p2*y))/d; //相交的点
return 1;
}
//四面体的有向体积
double volume4(point3 a,point3 b,point3 c,point3 d){
return dot3(cross(b-a,c-a),d-a);
}
//最小球覆盖
//三维凸包:增量法、懒得写了
posted on
浙公网安备 33010602011771号