**
计算几何总结（使用教材算法设计与分析（第二版）李春葆 清华大学出版社）
**
（代码大部分为书中原有代码，如有雷同，实属正常。）
#include<bits/stdc++.h>
using namespace std;//万能头
#include<fundament.h>//向量计算的函数库
本次学习主要主要有以下几个方面：1、向量计算 2、求解凸包问题 3、求解最近点对问题
（所有向量点已做成头文件fundament.h，文章中头文件fundament.h都是代表引用向量的库函数，附带库函数需要自取用法直接搞进库然后include调用，向量计算中的所有函数都已包含！！）
链接：https://pan.baidu.com/s/1O61l1KOswRa6_KR87l7wEg //这个是头文件fundament.h
提取码：l7j6
一、向量计算
(总括) 全部向量函数：
class Point //点类
{
public:
double x; //行坐标
double y; //列坐标
Point() {} //默认构造函数
Point(double x1,double y1) //重载构造函数
{ x=x1;
y=y1;
}
void disp()
{ printf("(%g,%g) ",x,y); }
friend bool operator ==(Point &p1,Point &p2); //重载==运算符
friend Point operator +(Point &p1,Point &p2); //重载+运算符
friend Point operator -(Point &p1,Point &p2); //重载-运算符
friend double Dot(Point p1,Point p2); //两个向量的点积
friend double Length(Point &p); //求向量长度
friend int Angle(Point p0,Point p1,Point p2); //求两线段p0p1和p0p2的夹角
friend double Det(Point p1,Point p2); //两个向量的叉积
friend int Direction(Point p0,Point p1,Point p2); //判断两线段p0p1和p0p2的方向
friend double Distance(Point p1,Point p2); //求两个点的距离
friend double DistPtoSegment(Point p0,Point p1,Point p2); //求p0到p1p2线段的距离
friend bool InRectAngle(Point p0,Point p1,Point p2); //判断点p0是否在p1和p2表示的矩形内
friend bool OnSegment(Point p0,Point p1,Point p2); //判断点p0是否在p1p2线段上
friend bool Parallel(Point p1,Point p2,Point p3,Point p4); //判断p1p2和p3p4线段是否平行
friend bool SegIntersect(Point p1,Point p2,Point p3,Point p4);//判断p1p2和p3p4两线段是否相交
friend bool PointInPolygon(Point p0,vector<Point> a); //判断点p0是否在点集a所形成的多边形内
friend double triangleArea(Point p1,Point p2,Point p3); //求三角形面积
friend double polyArea(vector<Point> p); //求多边形面积
};
1、基本运算符的重载
bool operator ==(Point &p1,Point &p2) //重载==运算符
{
return p1.x==p2.x && p1.y==p2.y;
}
Point operator +(Point &p1,Point &p2) //重载+运算符
{
return Point(p1.x+p2.x,p1.y+p2.y);
}
Point operator -(Point &p1,Point &p2) //重载-运算符
{
return Point(p1.x-p2.x,p1.y-p2.y);
}
2、向量的点积运算
double Dot(Point p1,Point p2) //两个向量的点积
{
return p1.x*p2.x+p1.y*p2.y;
}
double Length(Point &p) //求向量长度
{
return sqrt(Dot(p,p));
}
int Angle(Point p0,Point p1,Point p2) //求两线段p0p1和p0p2的夹角
{
double d=Dot((p1-p0),(p2-p0));
if (d==0)
return 0; //两线段p1p0和p2p0夹角为直角
else if (d>0)
return 1; //两线段p1p0和p2p0的夹角为钝角
else
return -1; //两线段p1p0和p2p0的夹角为锐角
}
3、向量叉积运算
double Det(Point p1,Point p2) //两个向量的叉积
{
return p1.x*p2.y-p1.y*p2.x;
}
int Direction(Point p0,Point p1,Point p2) //判断两线段p0p1和p0p2的方向
{
double d=Det((p1-p0),(p2-p0));
if (d==0)
return 0; //3点共线
else if (d>0)
return 1; //p0p2在p0p1的顺时针方向上
else
return -1; //p0p2在p0p1的逆时针方向上
}
4、点间距离
double Distance(Point p1,Point p2) //两个点之间的距离
{
return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
}
5、点到线段的距离
double DistPtoSegment(Point p0,Point p1,Point p2) //求p0到p1p2线段的距离
{
Point v1=p2-p1,v2=p1-p2,v3=p0-p1,v4=p0-p2;
if (p1==p2)
{
Point poo=p0-p1;
return Length(poo);
}
if (Dot(v1,v3)<0)
return Length(v3);
else if (Dot(v2,v4)<0)
return Length(v4);
else
return fabs(Det(v1,v3))/Length(v1);
}
6、判断点是否在矩形里
bool InRectAngle(Point p0,Point p1,Point p2) //判断点p0是否在p1和p2表示的矩形内
{
return Dot(p1-p0,p2-p0)<=0;
}
7、判断一个点是否在一条线段上
bool OnSegment(Point p0,Point p1,Point p2) //判断点p0是否在p1p2线段上
{
return Det(p1-p0,p2-p0)==0 && Dot(p1-p0,p2-p0)<=0;
}
8、判断两条线段是否平行
bool Parallel(Point p1,Point p2,Point p3,Point p4) //判断p1p2和p3p4线段是否平行
{
return Det(p2-p1,p4-p3)==0;
}
9、判断两条线段是否相交
bool SegIntersect(Point p1,Point p2,Point p3,Point p4) //判断p1p2和p3p4两线段是否相交
{
int d1,d2,d3,d4;
d1=Direction(p3,p1,p4); //求p3p1在p3p4的哪个方向上
d2=Direction(p3,p2,p4); //求p3p2在p3p4的哪个方向上
d3=Direction(p1,p3,p2); //求p1p3在p1p2的哪个方向上
d4=Direction(p1,p4,p2); //求p1p4在p1p2的哪个方向上
if (d1*d2<0 && d3*d4<0)
return true;
if (d1==0 && OnSegment(p1,p3,p4)) //若d1为0且p1在p3p4线段上
return true;
else if (d2==0 && OnSegment(p2,p3,p4)) //若d2为0且p2在p3p4线段上
return true;
else if (d3==0 && OnSegment(p3,p1,p2)) //若d3为0且p3在p1p2线段上
return true;
else if (d4==0 && OnSegment(p4,p1,p2)) //若d4为0且p4在p1p2线段上
return true;
else
return false;
}
10、判断一个点是否在多边形内
bool PointInPolygon(Point p0,vector<Point> a) //判断点p0是否在点集a所形成的多边形内
{ int i,cnt=0; //cnt累加交点个数
double x;
Point p1,p2;
for (i=0;i<a.size();i++)
{ p1=a[i]; p2=a[i+1];
if (OnSegment(p0,p1,p2))
return true; //如果点p0在多边形边p1p2上,返回true
//以下求解y=p0.y与p1p2的交点
if (p1.y==p2.y) continue; //如果p1p2是水平线,直接跳过
//以下两种情况是交点在p1p2延长线上
if (p0.y<p1.y && p0.y<p2.y) continue; //p0在p1p2线段下方,直接跳过
if (p0.y>=p1.y && p0.y>=p2.y) continue; //p0在p1p2线段上方,直接跳过
x=(p0.y-p1.y)*(p2.x-p1.x)/(p2.y-p1.y)+p1.x; //求交点坐标的x值
if (x>p0.x) cnt++; //只统计射线的一边
}
return (cnt%2==1);
}
11、求三角形面积
double triangleArea(Point p1,Point p2,Point p3) //求三角形面积
{
return fabs(Det(p1-p2,p3-p2))/2;
}
12、求多边形面积
double polyArea(vector<Point> p) //求多边形面积
{
double ans=0.0;
for (int i=1;i<p.size()-1;i++)
ans+=Det(p[i]-p[0],p[i+1]-p[0]);
return fabs(ans)/2;
}
二、求解凸包问题
1、礼品包裹法
思路： 包裹法本质在于查最大顺时针，先找到最左下点（注意：左更优先）并记录然后换下一个点找下一个点的最大顺时针点，记录最大顺时针点的最大顺时针点。这些最大顺时针点就是凸点。
#include<bits/stdc++.h>
#include<fundament.h>
using namespace std;
void search(vector<Point> a,Point min,int xx);
bool cmp(Point aj,Point ai,Point ak);
int main()
{
Point mins;
int i,n,minx,miny,xx;
cin>>n;
vector<Point> a;
for(i=0;i<n;i++)
{
Point xy;
scanf("%lf%lf",&xy.x,&xy.y);
a.push_back(xy);
}
minx=a[0].x;
miny=a[0].y;
mins=a.front();
for(i=0;i<n;i++)//优先检索最下元素
{
if(a[i].x<=minx)//先最左，若都在最左则找最下
{
if(a[i].x==minx)
{
if(a[i].y<=mins.y)
{
mins=a[i];
xx=i;
minx=a[i].x;
}
}
if(a[i].x<minx)
{
mins=a[i];
xx=i;
minx=a[i].x;
}
}
//cout<<it->x<<" "<<it->y<<endl;
}
cout<<"第一个凸包点是：("<<mins.x<<","<<mins.y<<")"<<endl;
search(a,mins,xx);
}
void search(vector<Point> a,Point mins,int xx)
{
int tmp=xx,k,i,j=xx;
vector<int> ch;
while(true)// 每次循环换个点，找到的点是顺时针最大的！！！
{
ch.push_back(j);
k=-1;
for(i=0;i<a.size();i++)//土包巧妙循环QAQ
if(i!=j&&(k==-1||cmp(a[j],a[i],a[k])))//条件有主要两个
//（1、不取本身2、不是首个或者在顺时针）
k=i;
if(k==tmp)break;
j=k;
}
for(int i=0;i<ch.size();i++)
cout<<ch[i]<<" ";
}
bool cmp(Point aj,Point ai,Point ak)
{//一个比较函数通过方向判断顺逆时针
int d=Direction(aj,ai,ak);
if(d==0)
return Distance(aj,ak)<Distance(aj,ai);
else if(d>0)
return true;
else
return false;
}
2、葛立恒扫描法（Grahem）
思路：Greham扫描主要重点在于，先找最下左点（注意：最下优先）和右拐原则的判断（右拐原则说白了就是栈顶和栈顶前一个连成的线与栈顶前一个的点与第i个点的方向问题，及判断是否向量a[top-1]a[top]与向量a[top-1][i]夹角大于0大于0直接入栈，小于零弹出栈顶，再判断下一个新的栈顶直到角度大于0，将a[i]进栈）
#include<bits/stdc++.h>
#include<fundament.h>
using namespace std;
void Greheng(vector<Point> &a,Point min1,int pm,int n);
bool cmp(Point &a,Point &b);
Point p0;//一定要有全局记录初始点位的方便
//cmp中基点作为基线与后者的顺逆时针判断
//为了用sort排出由角度小到大的顺序
int main()
{
vector<Point> a;
int n;
cin>>n;
for(int i=0;i<n;i++)
{
Point x1;
cin>>x1.x>>x1.y;
a.push_back(x1);
}
Point min1=a[0];
int pm=0;
for(int i=0;i<n;i++)//找基础点
{//先最下，若都在最下则找最左
if(a[i].y<=min1.y)
{
if(a[i].y==min1.y&&a[i].x<=min1.x)
{
min1=a[i];
pm=i;
}
if(min1.y>a[i].y)
{
min1=a[i];
pm=i;
}
}
}
cout<<"这个凸包型的基础凸包点是("<<min1.x<<","<<min1.y<<") 是点"<<pm<<endl;
p0=min1;
Greheng(a,min1,pm,n);
}
bool cmp(Point &a,Point &b)
{//比较方向角大小 ，等于也不用管
if(Direction(p0,a,b)>0)
return true;
else
return false;
}
void Greheng(vector<Point> &a,Point min1,int pm,int n)
{
Point ch[n];//这里ch用数组加top代替栈更有效率
//，若用否则vector则要先用reserve分配预先资源
//否则及其之慢
Point change;
change=a[0];
a[0]=min1;
a[pm]=change;
sort(a.begin()+1,a.end(),cmp);
int top=0;
ch[0]=a[0];
ch[1]=a[1];top++;
ch[2]=a[2];top++;
for(int i=3;i<n;i++)
{//条件较多但是容易理解的判断循环条件
//首先栈不为空，并列的满足满足右拐则弹出栈顶
//一直判断并弹到不满足右拐，当相等时，长的和右拐一样
while(top>=0&&(Direction(ch[top-1],a[i],ch[top])>0||(Direction(ch[top-1],a[i],ch[top])==0&&Distance(ch[top-1],a[i])>Distance(ch[top-1],a[top]))))
{
top--;
}
top++;
ch[top]=a[i];
}
cout<<"除了基点外还有："<<endl;
for(int i=1;i<top+1;i++)//打印注意是打印到top+1之前
cout<<"("<<ch[i].x<<","<<ch[i].y<<")"<<endl;
}
三、求解最近点对问题
旋转卡壳法以及求最大三角形面积
思路： 最大面积求发葛立恒法或包裹法之后，暴力三循环遍历，最大点间距离求法也可以暴力法双循环，卡壳什么的好麻烦。。。QAQ但效率高，思路是由基点和i+1连成的便作为底边开始遍历其他点与这两个点的叉积（即求面积），记下最大记下最大面积时对应点分别到边两点的距离，然后换到i+1和i+2边再次遍历找叉积最大面积，记下最大面积时对应点分别到边两点的距离，以此类推，最后比较这些 最大点距得出最大点间距。
#include<bits/stdc++.h>
#include<fundament.h>
//这个代码是一个葛立恒代码基础上用卡壳法求凸包点
//maxS函数实在卡壳法上附加的求最大面积的函数
using namespace std;
void Greheng(vector<Point> &a,Point min1,int pm,int n,Point ch[],int &tt);
// 这里是葛立恒查询凸包点代码块
double maxJammed(Point ch[],int tt);
// 这里是卡壳查找最长距离的代码
double maxS(Point ch[],int tt);
//这个是找最大面积的函数
bool cmp(Point a,Point b);
Point p0;//这个全局变量用于标记凸包点中的基点，主要用于cmp函数的比较
bool cmp(Point a,Point b)
{//比较规则用于sort函数
if(Direction(p0,a,b)>0)
return true;
else
return false;
}
double maxS(Point ch[],int tt)
{//葛立恒扫描完后，传入ch寄存的凸包点，tt是栈顶top+1
double smax=0;
for(int i=0;i<tt;i++)
{//双循环暴力卡出最大叉积并记录作为最大面积
for(int l=i+1;l<tt;l++)
{
for(int j=l+1;j<tt;j++)
{
if(fabs(Det(ch[i]-ch[l],ch[j]-ch[l]))/2.00>smax)
{
smax=fabs(Det(ch[i]-ch[l],ch[j]-ch[l]))/2.00;
}
//cout<<fabs(Det(ch[i]-ch[j],ch[l]-ch[j]))<<endl;
}
}
}
return smax;
}
int main()
{
vector<Point> a;
int n;
cin>>n;
for(int i=0;i<n;i++)
{
Point x1;
cin>>x1.x>>x1.y;
a.push_back(x1);
}
Point min1=a[0];
int pm=0;
for(int i=0;i<n;i++)//找基础点
{
if(a[i].y<=min1.y)
{
if(a[i].y==min1.y&&a[i].x<=min1.x)
{
min1=a[i];
pm=i;
}
if(min1.y>a[i].y)
{
min1=a[i];
pm=i;
}
}
}
Point ch[10000];
cout<<"这个凸包型的基础凸包点是("<<min1.x<<","<<min1.y<<") 是点"<<pm<<endl;
p0=min1;
int tt=0;
Greheng(a,min1,pm,n,ch,tt);
//传入基础点和基础点所在数组位置，传入tt带出栈顶
//，传入ch带出凸包点坐标
//cout<<ch.size()<<"chang"<<endl;
//for(int i=0;i<ch.size();i++)
//cout<<"("<<ch[i].x<<","<<ch[i].y<<")"<<endl;
//cout<<tt<<" "<<endl;
cout<<"最长为："<<maxJammed(ch,tt)<<endl;
cout<<"最大组成面积为："<<maxS(ch,tt)<<endl;
//maxJammed(ch,min1);
}
void Greheng(vector<Point> &a,Point min1,int pm,int n,Point ch[],int &tt)
{//葛立恒找凸包点
p0=min1;
Point change;
change=a[0];
a[0]=a[pm];
a[pm]=change;
sort(a.begin()+1,a.end(),cmp);
int top=0;
ch[0]=a[0];
ch[1]=a[1];top++;
ch[2]=a[2];top++;
for(int i=3;i<n;i++)
{
while(top>=0&&(Direction(ch[top-1],a[i],ch[top])>0||((Direction(ch[top-1],a[i],ch[top])==0)&&Distance(ch[top-1],a[i])>Distance(ch[top-1],ch[top]))))
{
top--;
}
top++;ch[top]=a[i];
}
for(int i=0;i<top+1;i++)
cout<<"("<<ch[i].x<<","<<ch[i].y<<")"<<endl;
tt=top+1;
return;
}
double maxJammed(Point ch[],int tt)//旋转卡壳找最大点间距离
{
int j=0;
double max=0;
vector<double> lon;
lon.reserve(300);
ch[tt-1]=ch[0];
for(int i=0;i<tt;i++)
{
int l=(i+1)%tt;
while(fabs(Det(ch[j+1]-ch[l],ch[i]-ch[l])>Det(ch[j]-ch[l],ch[i]-ch[l])))
{
j=(j+1)%tt;
//cout<<max<<" "<<endl;
}
double d1=Distance(ch[i],ch[j]);
double d2=Distance(ch[l],ch[j]);
lon.push_back(d1);
lon.push_back(d2);
j=0;
}
max=0;
for(int i=0;i<lon.size();i++)
{
if(max<lon[i])
max=lon[i];
cout<<lon[i]<<endl;
}
return max;
}
扫一扫
442
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
