ACM計算幾何總結
阿新 • • 發佈：2018-11-13
文章目錄
1 幾何基礎
#include <cstdio>
#include <cmath>
using namespace std;
const double pi = acos(-1.0);
const double inf = 1e100;
const double eps = 1e-6;
int sgn(double d){
if(fabs(d) < eps)
return 0;
if(d > 0)
return 2 點與向量
2.1 手動實現
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 Vector(A.x/p, A.y/p);
}
bool operator < (const Point& a, const Point& b){
if(a.x == b.x)
return a.y < b.y;
return a.x < b.x;
}
const double eps = 1e-6;
int sgn(double x){
if(fabs(x) < eps)
return 0;
if(x < 0)
return -1;
return 1;
}
bool operator == (const Point& a, const Point& b){
if(sgn(a.x-b.x) == 0 && sgn(a.y-b.y) == 0)
return true;
return false;
}
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;
}
double Area2(Point A, Point B, Point C){
return Cross(B-A, C-A);
}
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));
}
Vector Normal(Vector A){//向量A左轉90°的單位法向量
double L = Length(A);
return Vector(-A.y/L, A.x/L);
}
bool ToLeftTest(Point a, Point b, Point c){
return Cross(b - a, c - b) > 0;
}
2.2 複數黑科技
#include <complex>
using namespace std;
typedef complex<double> Point;
typedef Point Vector;//複數定義向量後，自動擁有建構函式、加減法和數量積
const double eps = 1e-9;
int sgn(double x){
if(fabs(x) < eps)
return 0;
if(x < 0)
return -1;
return 1;
}
double Length(Vector A){
return abs(A);
}
double Dot(Vector A, Vector B){//conj(a+bi)返回共軛複數a-bi
return real(conj(A)*B);
}
double Cross(Vector A, Vector B){
return imag(conj(A)*B);
}
Vector Rotate(Vector A, double rad){
return A*exp(Point(0, rad));//exp(p)返回以e為底複數的指數
}
3 點與線
3.1 直線定義
struct Line{//直線定義
Point v, p;
Line(Point v, Point p):v(v), p(p) {}
Point point(double t){//返回點P = v + (p - v)*t
return v + (p - v)*t;
}
};
3.2 求兩直線交點
//呼叫前需保證 Cross(v, w) != 0
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;
}
3.3 求點到直線距離
//點P到直線AB距離
double DistanceToLine(Point P, Point A, Point B){
Vector v1 = B-A, v2 = P-A;
return fabs(Cross(v1, v2)/Length(v1));
}//不去絕對值，得到的是有向距離
3.4 求點到線段距離
//點P到線段AB距離公式
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 DistanceToLine(P, A, B);
}
3.5 求點在直線上的投影點
//點P在直線AB上的投影點
Point GetLineProjection(Point P, Point A, Point B){
Vector v = B-A;
return A+v*(Dot(v, P-A)/Dot(v, v));
}
3.6 判斷點是否線上段上
//判斷p點是否線上段a1a2上
bool OnSegment(Point p, Point a1, Point a2){
return dcmp(Cross(a1-p, a2-p)) == 0 && dcmp(Dot(a1-p, a2-p)) < 0;
}
3.7 判斷兩線段是否相交
//判斷兩線段是否相交
bool SegmentProperIntersection(Point a1, Point a2, Point b1, Point b2){
double c1 = Cross(a2-a1, b1-a1), c2 = Cross(a2-a1, b2-a1);
double c3 = Cross(b2-b1, a1-b1), c4 = Cross(b2-b1, a2-b1);
//if判斷控制是否允許線段在端點處相交，根據需要新增
if(!sgn(c1) || !sgn(c2) || !sgn(c3) || !sgn(c4)){
bool f1 = OnSegment(b1, a1, a2);
bool f2 = OnSegment(b2, a1, a2);
bool f3 = OnSegment(a1, b1, b2);
bool f4 = OnSegment(a2, b1, b2);
bool f = (f1|f2|f3|f4);
return f;
}
return (sgn(c1)*sgn(c2) < 0 && sgn(c3)*sgn(c4) < 0);
}
4 三角形
4.1 求三角形外心
struct Point {
double x, y;
};
struct Line{
Point a, b;
};
Point Intersection(Line u, Line v){
Point ret = u.a;
double t1 = (u.a.x - v.a.x)*(v.a.y - v.b.y) - (u.a.y - v.a.y)*(v.a.x - v.b.x);
double t2 = (u.a.x - u.b.x)*(v.a.y - v.b.y) - (u.a.y - u.b.y)*(v.a.x - v.b.x);
double t = t1/t2;
ret.x += (u.b.x - u.a.x)*t;
ret.y += (u.b.y - u.a.y)*t;
return ret;
}
//外心
Point Circumcenter(Point a, Point b, Point c){
Line u, v;
u.a.x = (a.x + b.x)/2;
u.a.y = (a.y + b.y)/2;
u.b.x = u.a.x - a.y + b.y;
u.b.y = u.a.y + a.x - b.x;
v.a.x = (a.x + c.x)/2;
v.a.y = (a.y + c.y)/2;
v.b.x = v.a.x - a.y + c.y;
v.b.y = v.a.y + a.x - c.x;
return Intersection(u, v);
}
4.1 求三角形內心
struct Point {
double x, y;
};
struct Line{
Point a, b;
};
Point Intersection(Line u, Line v){
Point ret = u.a;
double t1 = (u.a.x - v.a.x)*(v.a.y - v.b.y) - (u.a.y - v.a.y)*(v.a.x - v.b.x);
double t2 = (u.a.x - u.b.x)*(v.a.y - v.b.y) - (u.a.y - u.b.y)*(v.a.x - v.b.x);
double t = t1/t2;
ret.x += (u.b.x - u.a.x)*t;
ret.y += (u.b.y - u.a.y)*t;
return ret;
}
//三角形內心
Point Incenter(Point a, Point b, Point c){
Line u, v;
double m, n;
u.a = a;
m = atan2(b.y - a.y, b.x - a.x);
n = atan2(c.y - a.y, c.x - a.x);
u.b.x = u.a.x + cos((m + n)/2);
u.b.y = u.a.y + sin((m + n)/2);
v.a = b;
m = atan2(a.y - b.y, a.x - b.x);
n = atan2(c.y - b.y, c.x - b.x);
v.b.x = v.a.x + cos((m + n)/2);
v.b.y = v.a.y + sin((m + n)/2);
return Intersection(u, v);
}
4.2 求三角形垂心
struct Point {