计算几何学习笔记
叉乘（叉积）
定义（二维）：
\[\vec{p_1} \times \vec{p_2} = x_1y_2-x_2y_1=-\vec{p_2}\times\vec{p_1}
\]
性质：
若 \(\vec{p_2}\) 在 \(\vec{p_1}\) 的逆时针方向，则结果大于 0 。
若 \(\vec{p_2}\) 在 \(\vec{p_1}\) 的顺时针方向，则结果小于 0 。
若 \(\vec{p_2}\) 与 \(\vec{p_1}\) 共线，则结果为 0 。
应用
- 叉积等于两个向量对应平行四边形的有向面积，三角形面积等于叉积的一般
- 判断一条折线是向左拐还是向右拐
- 判断两个线段 ab 和 cd 是否相交：
- \((\vec{ab}\times\vec{ac})\cdot(\vec{ab}\times\vec{ad})<0\) (点 c、d 分别在直线 ab 两侧) 并且 \((\vec{cd}\times\vec{ca})\cdot(\vec{cd}\times\vec{cb})<0\) (点 a、b 分别在直线 cd 两侧)
- \(\vec{ab}\times\vec{ac}=0\) （某端点在另一条直线上） 且 \(min(x_a,x_b)<=x_c<=max(x_a,x_b) \&\& min(y_a,y_b)<=y_c<=max(y_a,y_b)\) （该端点在另一条线段上）
- 多边形的面积：取一个点（通常为原点）作为起点，以起点到各顶点的连线为边分割成的三角形的面积（带正负）之和。
求多边形的重心
关于重心的定义：物理意义上的重心
三角形的重心公式：
\[x_0=\frac{x_1+x_2+x_3}{3}\\
y_0=\frac{y_1+y_2+y_3}{3}
\]
此公式若想应用到多边形，当且仅当多边形的重量全部分布在各个顶点上（公式内的系数为顶点的质量）。
质量均匀分布的多边形重心的求法：任取一个点与各个顶点连线组成三角形，各个三角形的重心的位置（三角形重心公式）和大小（三角形面积）可求。于是各三角形可以等价成各个重心，最后用加权的公式求解。
凸包
最小的凸多边形，满足覆盖了所有的点。
性质：凸包的所有点一定属于原来的点集。
格雷汉姆扫描法
取y值最小的点 \(p_0\) 作为原点。按照 \(p_0p_i\) 的斜率从小到大排序（用向量叉乘写 cmp 函数，放止精度问题）。开一个栈存凸包上的点，依次遍历每个点，检查栈顶的前两个元素与这个点构成的折线段是否“拐”向右侧（叉积 ≤ 0 ），是的话则不断弹出栈顶，直到折线段向左拐，此时将此点入栈。
时间复杂度为 \(O(n\log n)\)
贾维斯步进法
观察到确定原点后，斜率最小的点一定在凸包内（下一步一定向左转），那么可以每次找到斜率最小的点作为新的原点，不断更新凸包。
时间复杂度为 \(O(nh)\) ，h为凸包上的点个数（最坏情况为 n）（有一种死去的spfa的既视感）
旋转卡壳
用途：求凸包直径、用正方形框住所有点的最小面积等等。
核心思路：逆时针地遍历凸包上的边，对于每条边都找到离这条边最远的点，那么这时随着边的转动，对应的最远点也在逆时针旋转，不会有反向的情况。这样就保证了复杂度是O(n)的。
旋转卡壳 - OI Wiki yyds
习题
gxx's Problem
恶心了我足足两个周的题
先判断两个线段是否有交点，有的话输出交点，如果非整数则需要输出最简分数。
找两个不共线的线段的交点采用的是等底三角形面积之比等于高之比，进而再用相似三角形求算。详见此文章：CSDN博客
但是这个题的恶心之处在于需要判断两条线段是否有交点（一个还是无穷多个）以及是否共线。
总之就是分各种情况讨论，然后发现还是有的点wa了。
对拍正解，上千组数据也没拍出来。
最后在电梯里一想，会不会有线段是一个点的情况，试了试果然。
不得不说，好漂亮的数据点（微笑
#include<iostream>
#include<algorithm>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<queue>
#include<set>
#include<map>
#include<vector>
#include<iomanip>
#include<ctime>
#include<stack>
#define int long long
using namespace std;
inline int read(){
int x=0,f=1;char c=getchar();
while(!(c>='0'&&c<='9')) {if(c=='-') f=-1;c=getchar();}
while(c>='0'&&c<='9') {x=x*10+c-'0';c=getchar();}
return x*f;
}
long long cross(int x0,int y0,int x1,int y1){
return 1ll*x0*y1-1ll*x1*y0;
}
long long direction(int x0,int y0,int x1,int y1,int x2,int y2){
return cross(x1-x0,y1-y0,x2-x0,y2-y0);
}
bool onseg(int xa,int ya,int xb,int yb,int xc,int yc){
return min(xa,xb)<=xc&&xc<=max(xa,xb)&&min(ya,yb)<=yc&&yc<=max(ya,yb);
}
bool inseg(int xa,int ya,int xb,int yb,int xc,int yc){
return min(xa,xb)<xc&&xc<max(xa,xb)||min(ya,yb)<yc&&yc<max(ya,yb);
}
int intersection(int xa,int ya,int xb,int yb,int xc,int yc,int xd,int yd){
long long d1=direction(xa,ya,xb,yb,xc,yc);
long long d2=direction(xa,ya,xb,yb,xd,yd);
long long d3=direction(xc,yc,xd,yd,xa,ya);
long long d4=direction(xc,yc,xd,yd,xb,yb);
if(d1*d2<0&&d3*d4<0) return 1;
if(d1==0&&d2==0){
if(onseg(xa,ya,xb,yb,xc,yc)&&onseg(xa,ya,xb,yb,xd,yd)) return -1;
if(onseg(xc,yc,xd,yd,xa,ya)&&onseg(xc,yc,xd,yd,xb,yb)) return -1;
if(inseg(xa,ya,xb,yb,xc,yc)||inseg(xa,ya,xb,yb,xd,yd)) return -1;
if(inseg(xc,yc,xd,yd,xa,ya)||inseg(xc,yc,xd,yd,xb,yb)) return -1;
if(onseg(xa,ya,xb,yb,xc,yc)||onseg(xa,ya,xb,yb,xd,yd)||onseg(xc,yc,xd,yd,xa,ya)||onseg(xc,yc,xd,yd,xb,yb)) return 2;
return 0;
}
if(d1==0&&onseg(xa,ya,xb,yb,xc,yc)) return 1;
if(d2==0&&onseg(xa,ya,xb,yb,xd,yd)) return 1;
if(d3==0&&onseg(xc,yc,xd,yd,xa,ya)) return 1;
if(d4==0&&onseg(xc,yc,xd,yd,xb,yb)) return 1;
return 0;
}
long long gcd(long long a,long long b){
if(b==0) return a;
return gcd(b,a%b);
}
void Prin(long long a,long long b){
int f=1;
if(a<0) f*=-1,a=-a;
if(b<0) f*=-1,b=-b;
long long c=gcd(a,b);
a/=c;b/=c;
if(b==1) printf("%lld",f*a);
else printf("%lld/%lld",f*a,b);
}
void putinter(int x1,int y1,int px1,int py1,int x2,int y2,int px2,int py2){
int px3=x1-x2,py3=y1-y2;
long long fr1=cross(px2,py2,px3,py3),fr2=cross(px1,py1,px2,py2);
Prin(1ll*x1*fr2+1ll*fr1*px1,fr2);
putchar(' ');
Prin(1ll*y1*fr2+1ll*py1*fr1,fr2);
puts("");
}
signed main()
{
int T=read();
while(T--){
int xa=read(),ya=read(),xb=read(),yb=read();
int xc=read(),yc=read(),xd=read(),yd=read();
int num=intersection(xa,ya,xb,yb,xc,yc,xd,yd);
if(xa==xb&&ya==yb){
if(onseg(xc,yc,xd,yd,xa,ya)) printf("1\n%lld %lld\n",xa,ya);
else puts("0");
continue;
}
if(xc==xd&&yc==yd){
if(onseg(xa,ya,xb,yb,xc,yc)) printf("1\n%lld %lld\n",xc,yc);
else puts("0");
continue;
}
if(num==0) puts("0");
else if(num==-1) puts("INF");
else{
puts("1");
if(num==1) putinter(xa,ya,xb-xa,yb-ya,xc,yc,xd-xc,yd-yc);
else{
if(onseg(xa,ya,xb,yb,xc,yc)) printf("%lld %lld\n",xc,yc);
else printf("%lld %lld\n",xd,yd);
}
}
}
return 0;
}
Surround the Trees
凸包板子题
#include<iostream>
#include<algorithm>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<queue>
#include<set>
#include<map>
#include<vector>
#include<iomanip>
#include<ctime>
#include<stack>
using namespace std;
inline int read(){
int x=0,f=1;char c=getchar();
while(!(c>='0'&&c<='9')) {if(c=='-') f=-1;c=getchar();}
while(c>='0'&&c<='9') {x=x*10+c-'0';c=getchar();}
return x*f;
}
const int maxn=105;
struct node{
int x,y;
node operator +(const node &b){
node c;
c.x=x+b.x;
c.y=y+b.y;
return c;
}
node operator -(const node &b){
node c;
c.x=x-b.x;
c.y=y-b.y;
return c;
}
}a[maxn];
node ans[maxn];
stack<node> s;
int cross(node a,node b){
return a.x*b.y-a.y*b.x;
}
bool cmp1(node a,node b){
return a.y<b.y;
}
bool cmp2(node b,node c){
int cr=cross(b-a[1],c-a[1]);
return cr!=0?cr>0:b.x<c.x;
}
double cal(node a,node b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
int main()
{
int n=read();
while(n){
while(!s.empty()) s.pop();
for(int i=1;i<=n;i++) a[i].x=read(),a[i].y=read();
sort(a+1,a+n+1,cmp1);
sort(a+2,a+n+1,cmp2);
for(int i=1;i<=n;i++){
while(s.size()>=2){
node pre=s.top();s.pop();
node ppre=s.top();
if(cross(pre-ppre,a[i]-pre)>0){
s.push(pre);
break;
}
}
s.push(a[i]);
}
int cnt=0;
double anss=0;
while(!s.empty()){
ans[++cnt]=s.top();
s.pop();
}
for(int i=1;i<cnt;i++) anss+=cal(ans[i],ans[i+1]);
if(cnt>2) anss+=cal(ans[cnt],ans[1]);
printf("%.2lf\n",anss);
n=read();
}
return 0;
}
最大三角形
凸包上旋转卡壳板子题。
易得最大的三角形的顶点一定在凸包上。所以先 \(O(n\log n)\) 求出凸包。
然后类比找直径的过程，易得当三角形底边跨的点数一定时，对面的顶点一定也满足凸的性质。
所以最外层枚举三角形底边跨的点数 len，然后模拟指针指向对面的最远点，不断更新答案。
对于每一个 len 时间复杂度为 \(O(n)\) 的，所以总复杂度为 \(O(n^2)\) 。
#include<iostream>
#include<algorithm>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<queue>
#include<set>
#include<map>
#include<vector>
#include<iomanip>
#include<ctime>
#include<stack>
using namespace std;
inline int read(){
int x=0,f=1;char c=getchar();
while(!(c>='0'&&c<='9')) {if(c=='-') f=-1;c=getchar();}
while(c>='0'&&c<='9') {x=x*10+c-'0';c=getchar();}
return x*f;
}
const int maxn=200005;
int n;
struct node{
int x,y;
node operator +(const node &b)const{
return {x+b.x,y+b.y};
}
node operator -(const node &b)const{
return {x-b.x,y-b.y};
}
}a[maxn];
int cross(node a,node b){
return a.x*b.y-a.y*b.x;
}
bool cmp1(node a,node b){
return a.y<b.y;
}
bool cmp2(node aa,node bb){
return cross(aa-a[1],bb-a[1])>0;
}
double S(node a,node b){
return 0.5*abs(cross(a,b));
}
int main()
{
while(scanf("%d",&n)!=EOF){
for(int i=1;i<=n;i++){
a[i].x=read(),a[i].y=read();
}
sort(a+1,a+n+1,cmp1);
sort(a+2,a+n+1,cmp2);
stack<node> q;
for(int i=1;i<=n;i++){
while(q.size()>=2){
node pre=q.top();q.pop();
node pree=q.top();
if(cross(pre-pree,a[i]-pre)>0){
q.push(pre);
break;
}
}
q.push(a[i]);
}
int cnt=0;
while(!q.empty()){
a[++cnt]=q.top();
q.pop();
}
for(int i=1;i<=cnt;i++){
a[i+cnt*2]=a[i+cnt]=a[i];
}
double ans=0;
for(int len=1;len<cnt;len++){
int now=len+2;
for(int i=1;i<=cnt;i++){
int j=i+len;
while(S(a[j]-a[i],a[now]-a[i])<S(a[j]-a[i],a[now+1]-a[i])) now++;
ans=max(ans,S(a[j]-a[i],a[now]-a[i]));
}
}
printf("%.2lf\n",ans);
}
return 0;
}
