二维计算几何学习笔记

2023-10-26 13:16:22

三态函数dcmp

dcmp用来减少精度问题。

const double eps =1e-10;
int dcmp(double x) {
	if(fabs(x)<eps) return 0;
	else return x<0?-1:1;
}

点与向量

点与向量的关系

令点为\(D\)，向量为\(V\)，常数为\(C\)

则有：

\[D-D=V \]

\[V+V=V \]

\[D+V=D \]

\[D+D=\textbf{无意义} \]

\[V \times C=C \]

\[\frac{V}{C}=V \]

点与向量的实现

struct Point{
    double x,y;
    Point(double x=0,double y=0):x(x),y(y){ }
};

typedef Point Vector; // 偷个懒

Vector operator+(Vector A,Vector B){
    return Vector(A.x+B.x,A.y+B.y);
}

Vector operator-(Point A,Point B){
    return Vector(A.x-B.x,A.y-B.y);
}

Vector operator*(Vector A,double p){
    return Vector(A.x*p,A.y*p);
}

Vector operator/(Vector A,double p){
    return A*(1/p);
}
bool operator==(const Point& a,const Point &b){
	return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}
bool operator<(const Point &a,const Point &b){
	return a.x<b.x||(a.x==b.x&&a.y<b.y);
}

点积、长度与夹角

// 计算点积
double Dot(Vector A,Vector B){
	return A.x*B.x+A.y*B.y;
}

//长度
double Length(Vector A){
	return sqrt(Dot(A,A));
}

// 夹角
double Angle(Vector A,Vector B){
	return acos(Dot(A,B)/Length(A)/Length(B));
}

内积（点乘|点积）的几何意义包括：

表征或计算两个向量之间的夹角
b向量在a向量方向上的投影

from：Here

叉积与三角形面积

// 差积
double Cross(Vector A,Vector B){
	return A.x*B.y-A.y*B.x;
}
// 三角形面积
double Area2(Point A, Point B,Point C){
	return Cross(B-A,C-A);
}

向量旋转、单位法线

// 向量旋转，rad是弧度
Vector Rotate(Vector A,double rad){
	return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));
}

// 单位法线，A不能是零向量
Vector Normal(Vector A){
	double L = Length(A);
	return Vector(-A.y/L,A.x/L);
}

点与线的关系

直线交点

// 直线交点 两直线为f(t)=P+tv，g(t)=Q+tw。
Point GetLineIntersection(Point P,Vector v,Point Q,Vector w){
	Vector u=P-Q;
	double t=Cross(w,u)/Cross(v,w);
	return P+v*t;
}

点到直线的距离

// 点到直线的距离，P到直线AB的距离，directed确定是否为有向距离
double DistanceToLine(Point P,Point A,Point B,bool directed=false){
	Vector v1=B-A;
	Vector v2=P-A;
	return directed?(fabs(Cross(v1,v2))/Length(v1)):((Cross(v1,v2))/Length(v1));
}

点到线段的距离

double DistanceToSegment(Point P,Point A,Point B){
	if(A==B)return Length(P-A);
	Vector v1=B-A,v2=P-A,v3=P-B;
	if(dcmp(Dot(v1,v2))<0){
		return Length(v2);
	}
	if(dcmp(Dot(v1,v3))>0){
		return Length(v3);
	}
	return fabs(Cross(v1,v2))/Length(v1);
}

点在直线上的投影

// 点P在直线AB上的投影
Point GetLineProjection(Point P,Point A,Point B){
	Vector v=B-A;
	return A+v*(Dot(v,P-A)/Dot(v,v));
}

线段相交判定、点在直线上判定

// 线段相交判定，线段(a1,a2)与(b1,b2)
// 不含端点相交
bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2){
	double c1=Cross(a2-a1,b1-a1);
	double c2=Cross(a2-a1,b2-a1);
	double c3=Cross(b2-b1,a1-b1);
	double c4=Cross(b2-b1,a2-b1);
	return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}

// 点p是否在线段(a1,a2)上
bool OnSegment(Point p,Point a1,Point a2){
	return dcmp(Cross(a1-p,a2-p))==0&&dcmp(Dot(a1-p,a2-p))<0;
}

多边形

多边形有向面积

// 多边形有向面积：n边形面积
double ConvexPolygonArea(Point *p,int n){
	double area=0;
	for(int i=2;i<n;i++){
		area+=Cross(p[i]-p[1],p[i+1]-p[1]);
	}
	return area/2;
}

例题

UVA11178 Morley's Theorem

莫雷定理，关键在于计算\(D\)点。

首先需要计算\(\angle ABC\)的值\(\alpha\)，然后把射线\(BC\)逆时针旋转\(\frac{\alpha}{3}\)，得到直线\(BD\)，同理可得直线\(CD\)，然后求交点。

AC代码如下：

#include <bits/stdc++.h>
using namespace std;

struct Point{
	double x,y;
	Point(double x=0,double y=0):x(x),y(y){ }
};

typedef Point Vector; // 偷个懒

Vector operator+(Vector A,Vector B){
	return Vector(A.x+B.x,A.y+B.y);
}

Vector operator-(Point A,Point B){
	return Vector(A.x-B.x,A.y-B.y);
}

Vector operator*(Vector A,double p){
	return Vector(A.x*p,A.y*p);
}

Vector operator/(Vector A,double p){
	return A*(1/p);
}

// 计算点积
double Dot(Vector A,Vector B){
	return A.x*B.x+A.y*B.y;
}

//长度
double Length(Vector A){
	return sqrt(Dot(A,A));
}

// 夹角
double Angle(Vector A,Vector B){
	return acos(Dot(A,B)/Length(A)/Length(B));
}
// 差积
double Cross(Vector A,Vector B){
	return A.x*B.y-A.y*B.x;
}
// 向量旋转，rad是弧度
Vector Rotate(Vector A,double rad){
	return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));
}
// 直线交点 两直线为f(t)=P+tv，g(t)=Q+tw。
Point GetLineIntersection(Point P,Vector v,Point Q,Vector w){
	Vector u=P-Q;
	double t=Cross(w,u)/Cross(v,w);
	return P+v*t;
}

Point getD(Point A, Point B, Point C){
	Vector v1=C-B;
	double a1=Angle(A-B,v1);
	v1=Rotate(v1,a1/3);
	Vector v2=B-C;
	double a2=Angle(A-C,v2);
	v2=Rotate(v2,-a2/3);
	return GetLineIntersection(B,v1,C,v2);
}

void read_point(Point &poi){
	int x,y;
	cin>>x>>y;
	poi.x=x;
	poi.y=y;
}

int main(){
	int t;
	Point A,B,C,D,E,F;
	cin>>t;
	while(t--){
		read_point(A);
		read_point(B);
		read_point(C);
		D=getD(A,B,C);
		E=getD(B,C,A);
		F=getD(C,A,B);
		printf("%.6lf %.6lf %.6lf %.6lf %.6lf %.6lf\n",D.x,D.y,E.x,E.y,F.x,F.y);
	}
	return 0;
}

UVA11796 Dog Distance

主要代码：

const int MX=60;
int t,a,b;
Point p[MX],q[MX];
double minn,maxn;

void update(Point P,Point A,Point B){
	minn=min(minn,DistanceToSegment(P,A,B));
	maxn=max(maxn,Length(P-A));
	maxn=max(maxn,Length(P-B));
}
void read_point(Point &poi){
	int x,y;
	cin>>x>>y;
	poi.x=x;
	poi.y=y;
}

int main(){
	cin>>t;
	for(int kase=1;kase<=t;kase++){
		cin>>a>>b;
		for(int i=0;i<a;i++){
			read_point(p[i]);
		}
		for(int i=0;i<b;i++){
			read_point(q[i]);
		}
		double lena=0,lenb=0;
		for(int i=0;i<a-1;i++){
			lena+=Length(p[i+1]-p[i]);
		}
		for(int i=0;i<b-1;i++){
			lenb=Length(q[i+1]-q[i]);
		}
		int sa=0,sb=0;
		Point pa=p[0],pb=q[0];
		minn=1e9;maxn=-minn;
		while(sa<a-1&&sb<b-1){
			double la=Length(p[sa+1]-pa);
			double lb=Length(q[sb+1]-pb);
			double T = min(la/lena,lb/lenb);
			Vector va=(p[sa+1]-pa)/la*T*lena;
			Vector vb=(q[sb+1]-pb)/lb*T*lenb;
		
		update(pa,pb,pb+vb-va);
		pa=pa+va;
		pb=pb+vb;
		if(pa==p[sa+1]){
			sa++;
		}
		if(pb==q[sb+1]){
			sb++;
		}
		}
		printf("Case %d: %.0lf\n",kase,maxn-minn);
	}
	return 0;
}

结果\(\color{red}\textbf{ Unaccepted}\)了。但是样例过了，~~求放过~~。

码农公寓