1. 程式人生 > >計算幾何 模板

計算幾何 模板

push_back 代碼 struct 大量 data 定義 getc names lan

大量功能未(lan)經(de)測試等待後續測試

好吧是不知道怎麽測試....

先把模板備份在這裏.

大佬們心情好的話可以幫忙查查錯呀

#include <cstdio>
#include <cstdlib>
#include <cctype>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <cassert>
#include <vector>
#define ri rd<int>
#define rep(i, a, b) for (int i = (a), _ = (b); i <= _; ++ i)
#define per(i, a, b) for (int i = (a), _ = (b); i >= _; -- i) #define For(i, a, b) for (int i = (a), _ = (b); i < _; ++ i) #define forea(it, v) for (__typeof(v.begin()) it = v.begin(); it != v.end(); ++ it) using namespace std; typedef long double db; const db pi = acos(-1.); const db eps = 1e-10
; template<class T> T rd() { bool f = 1; char c = getchar(); for (; !isdigit(c); c = getchar()) if (c == ‘-‘) f = 0; T x = 0; for (; isdigit(c); c = getchar()) x = x * 10 + c - 48; return f ? x : -x; } int sgn(db x) { return x < -eps ? -1 : x > eps; } bool le(db x, db y) { return sgn(y-x) != -1
; } bool lt(db x, db y) { return sgn(y-x) == 1; } bool ge(db x, db y) { return sgn(x-y) != -1; } bool gt(db x, db y) { return sgn(x-y) == 1; } bool eq(db x, db y) { return sgn(x-y) == 0; } struct Vector; typedef Vector Point; struct Vector { db x, y; Vector(db _x = 0, db _y = 0):x(_x), y(_y) { } friend Vector operator + (const Vector &a, const Vector &b) { return Vector(a.x + b.x, a.y + b.y); } friend Vector operator - (const Vector &a, const Vector &b) { return Vector(a.x - b.x, a.y - b.y); } friend Vector operator * (const Vector &a, const db k) { return Vector(a.x * k, a.y * k); } friend Vector operator * (const db k, const Vector &a) { return Vector(a.x * k, a.y * k); } friend Vector operator / (const Vector &a, const db k) { return Vector(a.x / k, a.y / k); } Vector operator - () { // 反向 return Vector(-x, -y); } friend bool operator < (const Point &a, const Point &b) { // 凸包排序用 return eq(a.x, b.x) ? lt(a.y, b.y) : lt(a.x, b.x); } db len() const { // 模長 return sqrt(x * x + y * y); } Vector unit() const { // 單位向量 return *this / len(); } Vector normal() const { // 法向量 return Vector(-y, x); } }; db dot(const Vector &a, const Vector &b) { return a.x * b.x + a.y * b.y; } db det(const Vector &a, const Vector &b) { return a.x * b.y - a.y * b.x; } db p_dis_p(const Point&, const Point&); Point p_middle_p(const Point&, const Point&); struct Line { // 直線 Point a, b; // 兩點確定一線, 方向向量定義為 a->b Line(Point _a, Point _b):a(_a), b(_b) { } bool chk_p(const Point &p) const { return 1; } }; struct Seg:Line { // 線段 Seg(Point _a, Point _b):Line(_a, _b) {} bool chk_p(const Point &p) const { return le(dot(p - a, p - b), 0); } Point middle() const { // 線段中點 return (a + b) / 2; } Line mperp() const { // 中垂線 Point u = middle(); return Line(u, u + (b - a).normal()); } db length() const { return p_dis_p(a, b); } }; struct Arrow:Line { // 射線 Arrow(Point _a, Point _b):Line(_a, _b) {} bool chk_p(const Point &p) const { return ge(dot(p - a, b - a), 0); } }; struct Circle { Point o; db r; Circle() {} Circle(Point _o, db _r):o(_o), r(_r) { } // 點徑型構造函數 Circle(Point a, Point b) { // 直徑型構造函數 o = p_middle_p(a, b); r = p_dis_p(a, o); } bool chk_p(const Point &p) const { return eq(p_dis_p(p, o), r); } }; struct Arc:Circle { Point a, b; // a 逆時針轉至 b Arc(Point _o, db _r, Point _a, Point _b):Circle(_o, _r), a(_a), b(_b) { } Arc(Circle _c, Point _a, Point _b):a(_a), b(_b) { o = _c.o, r = _c.r; } bool chk_p(const Point &p) const { return eq(p_dis_p(p, o), r) && le(dot(p - a, p - b), 0); } Arc bu() const { // 補弧 return Arc(o, r, b, a); } db bad() const { // 劣弧 | 半圓 ? return ge(det(a - o, b - o), 0); } db length() const { // 弧長 if (bad()) { Point u = p_middle_p(a, b); db x = p_dis_p(o, u); db y = p_dis_p(a, u); return atan2(y, x) * r; } else return 2 * pi * r - bu().length(); } db shan() const { // 扇形 return length() * r / 2; } db gong() const { // 弓形 Point u = p_middle_p(a, b); db x = p_dis_p(o, u); db y = p_dis_p(a, u); return atan2(y, x) * r * r / 2 - x * y; } }; template<class T1, class T2> bool l_parallel_l(const T1 a, const T2 b) { return eq(det(a.b - a.a, b.b - b.a), 0); } Point p_middle_p(const Point &a, const Point &b) { return (a + b) / 2; } template<class T> Point p_projection_l(const Point &p, const T &l) { Vector v = (l.b - l.a).unit(); return l.a + dot(p - l.a, v) * v; } template<class T> Seg seg_projection_l(const Seg &s, const T &l) { return Seg(p_projection_l(s.a, l), p_projection_l(s.b, l)); } Point p_over_p(const Point &a, const Point &b) { return 2 * b - a; } template<class T> Point p_over_l(const Point &p, const T &l) { return p_over_p(p_projection_l(p, l)); } db p_dis_p(const Point &a, const Point &b) { return (b - a).len(); } template<class T> db p_dis_l(const Point &p, const T &l) { Point q = p_projection_l(p, l); if (l.chk_p(q)) return p_dis_p(p, q); else return min(p_dis_p(p, l.a), p_dis_p(p, l.b)); } int c_relative_c(Circle a, Circle b) { // 代碼中統一不考慮兩圓重合的情況 if (lt(a.r, b.r)) swap(a, b); db d = p_dis_p(a.o, b.o); if (gt(d, a.r + b.r)) return 2; // 相離 if (eq(d, a.r + b.r)) return 1; // 外切 if (eq(d, a.r - b.r)) return -1; // 內切 if (lt(d, a.r - b.r)) return -2; // 內含 return 0; // 相交 } #define PUT(x) if (1) { Point _ = x; if (a.chk_p(_) && b.chk_p(_)) S.push_back(_); } template<class T1, class T2> vector<Point> l_intersection_l(const T1 &a, const T2 &b) { // 代碼中統一不考慮兩線在同一直線上的情況 vector<Point> S; if (l_parallel_l(a, b)) return S; db s1 = det(b.a - a.a, b.b - a.a); db s2 = det(b.b - a.b, b.a - a.b); PUT(a.a + (a.b - a.a) * s1 / (s1 + s2)); return S; } template<class T1, class T2> vector<Point> l_intersection_c(const T1 &a, const T2 &b) { vector<Point> S; Point u = p_projection_l(b.o, a); db x = p_dis_p(b.o, u); if (lt(b.r, x)) return S; if (eq(b.r, x)) { PUT(u); return S; } Vector v = (a.b - a.a).unit(); db y = sqrt(b.r * b.r - x * x); PUT(u + v * y); PUT(u - v * y); return S; } template<class T1, class T2> vector<Point> c_intersection_c(T1 a, T2 b) { if (lt(a.r, b.r)) swap(a, b); vector<Point> S; db d = p_dis_p(a.o, b.o); int kd = c_relative_c(a, b); if (abs(kd) == 2) return S; else if (abs(kd) == 1) { PUT(a.o + (b.o - a.o).unit() * a.r); return S; } db x = fabs((d + (a.r * a.r - b.r * b.r) / d) / 2); // 分類討論兩圓相交的各種情況, 式子相同, 符號不一 db z = sqrt(a.r * a.r - x * x); Vector u = (b.o - a.o).unit(), v = u.normal(); Point c = a.o + u * x; PUT(c + v * z); PUT(c - v * z); return S; } #undef PUT vector<Point> p_tangentp_c(const Point &p, const Circle &c) { Circle b(p, c.o); return c_intersection_c(b, c); } vector<Line> p_tangentl_c(const Point &p, const Circle &c) { vector<Point> S = p_tangentp_c(p, c); vector<Line> T; forea (it, S) T.push_back(Line(*it, p)); return T; } // 求兩圓的共切點的話, 四個點關系也不知道, 意義不大, 直接求線 vector<Line> c_itangentl_c(Circle a, Circle b) { // 內共切線 vector<Line> T; if (lt(a.r, b.r)) swap(a, b); int kd = c_relative_c(a, b); if (kd <= 0) return T; if (kd == 1) { Point u = a.o + (b.o - a.o).unit() * a.r; T.push_back(Line(u, u + (b.o - a.o).normal())); return T; } Circle c(a.o, a.r + b.r); Circle d(a.o, b.o); vector<Point> S = c_intersection_c(c, d); forea (it, S) { Point u = a.o + (*it - a.o).unit() * a.r; T.push_back(Line(u, u + (b.o - *it))); } return T; } vector<Line> c_extangentl_c(Circle a, Circle b) { // 外公切線 vector<Line> T; if (eq(a.r, b.r)) { Point u = a.o + (b.o - a.o).unit().normal() * a.r; T.push_back(Line(u, u + (b.o - a.o))); u = a.o - (b.o - a.o).unit().normal() * a.r; T.push_back(Line(u, u + (b.o - a.o))); return T; } if (lt(a.r, b.r)) swap(a, b); int kd = c_relative_c(a, b); if (kd == -2) return T; if (kd == -1) { Point u = a.o + (b.o - a.o).unit() * a.r; T.push_back(Line(u, u + (b.o - a.o).normal())); return T; } Circle c(a.o, a.r - b.r); Circle d(a.o, b.o); vector<Point> S = c_intersection_c(c, d); forea (it, S) { Point u = a.o + (*it - a.o).unit() * a.r; T.push_back(Line(u, u + (b.o - *it))); } return T; } int main() { // test data Circle a(Point(10, 10), 10); Circle b(Point(16, 14), 11.313708498984761); int kd = c_relative_c(a, b); kd = c_relative_c(b, a); vector<Point> S = c_intersection_c(a, b); S = c_intersection_c(b, a); vector<Line> T1 = c_extangentl_c(a, b); T1 = c_extangentl_c(b, a); vector<Line> T2 = c_itangentl_c(a, b); T2 = c_itangentl_c(b, a); return 0; }

計算幾何 模板