1.先来点对于精度要求的辅助函数
double PI = acos(-1);
double INF = 1e20;
double EPS = 1e-6;
bool IsZero(double x) {
return - EPS < x && x < EPS;
} // 是否为0
int dcmp(double x) {
if(fabs(x) < EPS) return 0;
else return x < 0 ? -1 : 1;
} // 判断是正是负
2.矢量、点、直线的定义以及四则运算
#define Point Vector // 点和向量的表示法都是一样的
struct Vector{
double x,y;
Vector(double x=0,double y=0) : x(x),y(y) {}
};
Vector operator + (Vector p,Vector q) {
return Vector(p.x + q.x, p.y + q.y);
} // 加
Vector operator -(Vector p, Vector q) {
return Vector(p.x - q.x, p.y - q.y);
} // 减
Vector operator *(double k, Vector p) {
return Vector(k * p.x, k * p.y);
} // 数乘
double operator *(Vector p, Vector q) {
return p.x * q.x + p.y * q.y;
} // 点乘
struct Line{
Point a,b;
Line(Point a,Point b) : a(a),b(b) {}
}; // 直线的定义 由两个点构成
3.点积的功能
求同向还是异向;(>0同向,=0垂直,<0异向)
求投影;
求出投影后用勾股定理求点到直线距离;
double dot(Vector p, Vector q) {
return p.x*q.x+p.y*q.y;
} // 点积
double project(Vector p, Vector n) {
return dot(p, unit(n));
} // p在n上的投影长度
4.叉积的功能
求面积(三角形面积等于任意两边求叉积/2);
求顺时针方向还是逆时针方向;(即a->b是顺时针转还是逆时针转)
若a 逆时针旋转小于180 度可到b, 则结果为正,否则结果为负
判断是否在半平面上
注意使用使用叉积的时候外面一定要加上(),也就是要写成(p^q)
double operator ^(Vector p, Vector q) {
return p.x * q.y - q.x * p.y;
} // 叉乘
double Cross(Vector p, Vector q) {
return p.x * q.y - q.x * p.y;
} // 叉积
5.一些其他向量/点之间的运算
1.给两个向量,求围成的三角形面积
2.求两点之间的距离
3.求一个点旋转alpha角度得到的点
double area(Vector p, Vector q) {
return (p^q)/2.0;
}
double dist(Point p, Point q) {
return length(p - q);
} // 两点之间的距离
Point rotate(Point b, Point a,double alpha) {
Vector p = b - a;
return Point(a.x + (p.x * cos(alpha)- p.y * sin(alpha)),a.y + (p.x * sin(alpha) + p.y * cos(alpha)) );
} // b绕a旋转alpha得到c点,顺时针旋转,用旋转矩阵理解
6.点与直线之间的一些操作
1.求点到直线的距离
2.判断点在向量AB的那一边
3.判断点是否在线段AB上
4.求点P在直线L上的垂足(用于求对称点)
5.求过点P关于L的垂线
//1.求点到直线的距离
//2.判断点在向量AB的那一边
//3.判断点是否在线段AB上
//4.求点P在直线L上的垂足
//5.求过点P关于L的垂线
double dist(Point p, Line l) {
return fabs((p - l.a) ^ (l.b - l.a)) / length(l.b-l.a);
} // 点到直线距离
int sideOfLine(Point p, Point a, Point b)
{
double result = (b - a) ^ (p - a);
if (IsZero(result)) return 0; //p 在 AB上
if (dcmp(result) > 0) return 1; //p在AB左侧
return -1; //p在AB右侧
} //判断p在直线 AB的哪一侧
bool OnSegment (Point p, Point a1, Point a2) {
if (dcmp(Cross(a1-p,a2-p))) return 0;
else
if (dcmp(p.x-min(a1.x,a2.x))>=0 && dcmp(p.x-max(a1.x,a2.x))<=0 &&
dcmp(p.y-min(a1.y,a2.y))>=0 && dcmp(p.y-max(a1.y,a2.y))<=0 ) return 1;
else return 0;
} // 点是否在线段上
Point foot(Point p, Line l) {
return l.a + project(p - l.a, l.b - l.a)*unit(l.b - l.a);
} // p点在l上的垂足
Line Vertical(Point p, Line l) {
return Line(p,p + (rotate(l.b, l.a, PI / 2) - l.a));
} // 求过p点的l的垂线
7.直线与直线之间的操作
1.求两个直线的交点
2.求两个直线的夹角
Point intersect(Line l, Line m, string msg) {
double x = area(m.a - l.a, l.b - l.a);
double y = area(l.b - l.a, m.b - l.a);
if (IsZero(x + y)) {
if (IsZero(dist(l.a, m))) msg = "重合";
else msg = "平行";
return NULL;
}
return m.a + x / (x + y) * (m.b - m.a);
} // 两条直线的交点,然后按照比例相加减即可
double angle(Line l, Line m) {
return acos( fabs(project(l.b - l.a, m.b - m.a)) / length(l.b - l.a));
} // 两条线的夹角
8.求多边形面积
使用有向面积计算
按逆时针方向顺次为各边指定方向
对于每条边AB,累加(A^B)/2
9.判断点是否在多边形内部
从这个点引一条射线,统计多边形穿过这条线正反多少次,设其在
逆时针穿过加一 顺时针穿过减一
最后wn!=0的话说明在内部
#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#define Point Vector
using namespace std;
/*
为了简略篇幅
这边省略了一些基本函数
到时候看自己代码库就行
*/
Point poly[10000];
Point p;
int n,i;
/*-----------------------------------------------------------------------------*/
int is_Point_in_Poly(Point p, Point *poly, int n)
{
int wn = 0;
for(int i = 1; i <= n; i++)
{
if (OnSegment(p, poly[i], poly[(i+1)%n])) return 0; // 注意这边不可以使用SideOfLine
int k = dcmp(Cross(poly[(i+1)%n]-poly[i], p-poly[i]));
int d1 = dcmp(poly[i].y-p.y);
int d2 = dcmp(poly[(i+1)%n].y-p.y);
if (k>0 && d1<=0 && d2>0) wn++;
if (k<0 && d2<=0 && d1>0) wn--;
}
if (wn!=0) return 1;
else return -1;
}
/*-----------------------------------------------------------------------------*/
int main()
{
scanf("%d",&n);
for (i=1; i<=n; i++) cin>>poly[i].x>>poly[i].y;
cin>>p.x>>p.y;
int ans=is_Point_in_Poly(p,poly,n);
if (ans>0) cout<<"内部"<<endl;
if (ans==0) cout<<"边上"<<endl;
if (ans<0) cout<<"外部"<<endl;
}
10.凸包
直观来看,把点集 S 中所有点钉在一个木板上,
用一个橡皮筋框起所有点来,一撒手,橡皮筋就表示出了 S 的凸包
任意两点的连线都在凸包内
凸包算是一个点集了
凸包的求法 基于水平序的Andrew算法(具体详解见刘汝佳p272)
1、先把所有点按x排序,x相同按y排序
2、然后依次加点,当新点在凸包左边就继续,否则就继续删除最新加入的点直到新点在左边
3、重复这一过程直到碰到最右边的点,这样求出了下凸包,然后反过来从p[n]做一次就求出了上凸包
#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#define Point Vector
using namespace std;
/*
为了简略篇幅
这边省略了一些基本函数
到时候看自己代码库就行
*/
Point p[10000];
Point ch[10000];
int n,m,i;
bool cmp(Point p,Point q) {
if (p.x<q.x) return 1;
if (p.x==q.x && p.y<q.y) return 1;
return 0;
}
/*-------------------------------------------------------------------------------*/
int ConvexHull(Point *p, int n, Point*ch ) //Andrew算法
{
sort(p+1,p+n+1,cmp);
int m = 0;
for (int i = 1; i <= n; i++)
{
while (m>1 && sideOfLine(p[i],ch[m-2],ch[m-1])<=0 ) m--;
// 或者写成 dcmp( Cross(ch[m-1]-ch[m-2], p[i]-ch[m-2]) )<=0
ch[m++]=p[i];
}
int k = m;
for (int i = n-1; i >= 1; i--)
{
while (m>k && sideOfLine(p[i],ch[m-2],ch[m-1])<=0 ) m--;
ch[m++]=p[i];
}
if (n>1) m--;
return m;
}
int main()
{
scanf("%d",&n);
for (i=1; i<=n; i++) cin>>p[i].x>>p[i].y;
m=ConvexHull(p,n,ch);
cout<<"凸包"<<endl;
for (i=0; i<m; i++) cout<<ch[i].x<<' '<<ch[i].y<<endl;
}