二维平面计算几何模板
1 //author Eterna 2 #include<iostream> 3 #include<algorithm> 4 #include<cstdio> 5 #include<vector> 6 #include<cstring> 7 #include<string> 8 #include<cmath> 9 #include<cstdlib> 10 #include<utility> 11 #include<deque> 12 #include<queue> 13 using namespace std; 14 //------------------------------点与直线部分----------------------------- 15 const double pi = acos(-1.0); 16 const double eps = 1e-10; 17 //double比较和0的大小关系 18 inline int Dcmp(double x) { 19 if (fabs(x) < eps)return 0; 20 else return x < 0 ? -1 : 1; 21 } 22 struct point { 23 friend istream& operator >>(istream& in, point& rhs) { 24 in >> rhs.x >> rhs.y; 25 return in; 26 } 27 friend ostream& operator <<(ostream& os, const point& rhs) { 28 os << rhs.x << ' ' << rhs.y; 29 return os; 30 } 31 bool operator == (const point& rhs)const { 32 return Dcmp(x - rhs.x) == 0 && Dcmp(y - rhs.y) == 0; 33 } 34 bool operator < (const point& rhs)const { 35 return x < rhs.x || (x == rhs.x && y < rhs.y); 36 } 37 bool operator >(const point& rhs)const { 38 return !((*this) < rhs) && !((*this) == rhs); 39 } 40 point operator + (const point& p)const { 41 return point(x + p.x, y + p.y); 42 } 43 point operator - (const point& p)const { 44 return point(x - p.x, y - p.y); 45 } 46 point operator * (double p)const { 47 return point(x * p, y * p); 48 } 49 point operator / (double p)const { 50 return point(x / p, y / p); 51 } 52 //点乘 53 double dot(const point& p) const { 54 return x * p.x + y * p.y; 55 } 56 //叉乘 57 double det(const point& p)const { 58 return x * p.y - y * p.x; 59 } 60 point() {} 61 point(double _x, double _y) :x(_x), y(_y) {} 62 double x, y; 63 }; 64 struct line { 65 friend istream& operator >>(istream& in, line& rhs) { 66 in >> rhs.p >> rhs.v; 67 return in; 68 } 69 friend ostream& operator <<(ostream& os, const line& rhs) { 70 os << rhs.p << ' ' << rhs.v; 71 return os; 72 } 73 bool operator < (const line& rhs)const { 74 if (Dcmp(ang - rhs.ang))return ang < rhs.ang; 75 else return Dcmp(v.det(rhs.v)) == -1; 76 } 77 point Get_Point(double t)const { 78 return p + v * t; 79 } 80 bool IsParallel(const line& rhs)const { return Dcmp(v.det(rhs.v)) == 0; } 81 point p, v; 82 double ang; 83 line() {} 84 line(point _p, point _v) :p(_p), v(_v) { ang = atan2(_v.y, _v.x); } 85 }; 86 //不想写scanf 87 inline point read_point() { 88 double x, y; 89 scanf("%lf %lf", &x, &y); 90 return point(x, y); 91 } 92 //向量长度 93 inline double length(const point& a) { return sqrt(a.dot(a)); } 94 //两点之间长度的平方 95 inline double Dist(const point& a, const point& b) { return (a - b).dot(a - b); } 96 //两向量的夹角 97 inline double Angle(const point& a, const point& b) { return acos(a.dot(b) / (length(a) * length(b))); } 98 //三点无向面积公式 99 inline double Area2(const point& a, const point& b, const point& c) { return fabs((b - a).det(c - a)) / 2; } 100 //判断p q两点是否在线段ab的同侧 101 inline bool SameSide(const point& a, const point& b, const point& p, const point& q) { 102 point ab = b - a, ac = p - a, ap = q - a; 103 double res = ab.det(ac) * ab.det(ap); 104 return Dcmp(res) >= 0; 105 } 106 //判断点是否在三角形内 107 inline bool InTriangle(const point& a, const point& b, const point& c, const point& p) { return SameSide(a, b, c, p) && SameSide(b, c, a, p) && SameSide(c, a, b, p); } 108 //逆时针旋转向量 109 inline point Rotate(const point& rhs, double rad) { return point(rhs.x * cos(rad) - rhs.y *sin(rad), rhs.x * sin(rad) + rhs.y * cos(rad)); } 110 //求向量单位法向量 111 inline point Normal(const point& rhs) { 112 if (Dcmp(rhs.x) == 0 && Dcmp(rhs.y) == 0)return point(0.0, 0.0); 113 double len = length(rhs); 114 return point(-rhs.y / len, rhs.x / len); 115 } 116 //求单位向量 117 inline point Unit_vector(const point& rhs) { 118 double len = length(rhs); 119 return point(rhs.x / len, rhs.y / len); 120 } 121 //判断q点是否在线段p1, p1上 122 inline bool On_seg(const point& p1, const point& p2, const point& q) { return Dcmp((p1 - q).det(p2 - q)) == 0 && Dcmp((p1 - q).dot(p2 - q)) < 0; } 123 //判断p是否在直线左边 124 inline bool On_Left(const line& L, const point& p) { return L.v.det(p - L.p) > 0; } 125 //通过4点坐标求两直线交点 126 inline point Intersection_point(const point& p1, const point& p2, const point& q1, const point& q2) { 127 return p1 + (p2 - p1) * ((q2 - q1).det(q1 - p1) / (q2 - q1).det(p2 - p1)); 128 } 129 //通过直线上一点及其方向向量求交点 130 inline point Intersection_line(const line& L1, const line& L2) { 131 point u = L1.p - L2.p; 132 double t = L2.v.det(u) / L1.v.det(L2.v); 133 return L1.p + L1.v * t; 134 } 135 //判断两线段是否规范相交 136 inline bool SegmentProperIntersection(const point& a1, const point& a2, const point& b1, const point& b2) { 137 double c1 = (a2 - a1).det(b1 - a1), c2 = (a2 - a1).det(b2 - a1), c3 = (b2 - b1).det(a1 - b1), c4 = (b2 - b1).det(a2 - b1); 138 return Dcmp(c1) * Dcmp(c2) < 0 && Dcmp(c3) * Dcmp(c4) < 0; 139 } 140 //判断两线段是否相交，判断所有情况 141 inline bool SegmentIntersection(const point& a1, const point& a2, const point& b1, const point& b2) { 142 if (Dcmp((a2 - a1).det(b2 - b1)) == 0)return On_seg(a1, a2, b1) || On_seg(a1, a2, b2) || On_seg(b1, b2, a1) || On_seg(b2, b2, a2); 143 double c1 = (a2 - a1).det(b1 - a1), c2 = (a2 - a1).det(b2 - a1), c3 = (b2 - b1).det(a1 - b1), c4 = (b2 - b1).det(a2 - b1); 144 return Dcmp(c1) * Dcmp(c2) < 0 && Dcmp(c3) * Dcmp(c4) < 0; 145 } 146 //求点p到直线ab的垂直距离 147 inline double DistanceToLine(const point& p, const point& a, const point& b) { return 2 * Area2(p, a, b) / length(a - b); } 148 //求点p到线段ab的距离 149 inline double DistanceToSegment(const point& p, const point& a, const point& b) { 150 if (a == b)return length(p - a); 151 point v1 = b - a, v2 = p - a, v3 = p - b; 152 if (Dcmp(v1.dot(v2) < 0))return length(v2); 153 else if (Dcmp(v1.dot(v3) > 0))return length(v3); 154 else return DistanceToLine(p, a, b); 155 } 156 //求点p在直线ab上的投影点 157 inline point GetLineProjection(const point& p, const point& a, const point& b) { 158 point v = b - a; 159 return a + v * (v.dot(p - a) / v.dot(v)); 160 } 161 //计算凸多边形面积 162 double ConvexPolygonArea(point* p, int n) { 163 double res = 0; 164 for (int i = 1; i < n - 1; ++i)res += Area2(p[0], p[i], p[i + 1]); 165 return res; 166 } 167 //计算任意多边形面积公式 168 double PolygonArea(point* p, int n) { 169 double res = 0; 170 for (int i = 1; i < n - 1; ++i)res += (p[i] - p[0]).det(p[i + 1] - p[0]); 171 return res / 2; 172 } 173 //判断点是否在多边形内 174 int isPointinPolygon(const point& p, vector<point>& v) { 175 int cnt = 0, n = v.size(); 176 for (int i = 0; i < n; ++i) { 177 if (On_seg(v[i], v[(i + 1) % n], p))return -1; 178 int k = Dcmp((v[(i + 1) % n] - v[i]).det(p - v[i])); 179 int d1 = Dcmp(v[i].y - p.y), d2 = Dcmp(v[(i + 1) % n].y - p.y); 180 if (k > 0 && d1 <= 0 && d2 > 0)++cnt; 181 if (k < 0 && d2 <= 0 && d1 > 0)--cnt; 182 } 183 if (cnt)return 1;//内部 184 else return 0;//外部 185 } 186 vector<point> v; 187 //构造凸包 188 inline void Get_convex_hall(point* arr, int n) { 189 v.clear(); 190 v.resize(n << 1); 191 sort(arr + 1, arr + 1 + n); 192 int cnt = 0; 193 for (int i = 1; i <= n; ++i) { 194 while (cnt > 1 && (v[cnt - 1] - v[cnt - 2]).det(arr[i] - v[cnt - 1]) <= 0)--cnt; 195 v[cnt++] = arr[i]; 196 } 197 for (int i = n - 1, t = cnt; i > 0; --i) { 198 while (cnt > t && (v[cnt - 1] - v[cnt - 2]).det(arr[i] - v[cnt - 1]) <= 0)--cnt; 199 v[cnt++] = arr[i]; 200 } 201 v.resize(cnt - 1); 202 } 203 //旋转卡壳 204 double Rotate_Calipers() { 205 double res = 0.0; 206 int n = v.size(); 207 v.push_back(v[0]); 208 int now = 1; 209 for (int i = 0; i < n; ++i) { 210 while ((v[now] - v[i + 1]).det(v[i] - v[i + 1]) < (v[now + 1] - v[i + 1]).det(v[i] - v[i + 1]))now = (now + 1) % n; 211 res = max(res, max(Dist(v[now], v[i]), Dist(v[now + 1], v[i + 1]))); 212 } 213 v.pop_back(); 214 return res; 215 } 216 //求半平面交 217 bool HalfPlaneIntersection(line* arr, int n, vector<point>& v) { 218 deque<line> qLine; 219 deque<point> qPoint; 220 sort(arr + 1, arr + 1 + n); 221 qLine.push_back(arr[1]); 222 for (int i = 2; i <= n; ++i) { 223 if (Dcmp(arr[i].ang - arr[i - 1].ang)) { 224 if (qLine.size() > 1 && (qLine[0].IsParallel(qLine[1]) || qLine[qLine.size() - 1].IsParallel(qLine[qLine.size() - 2])))return false; 225 while (qLine.size() > 1 && !On_Left(arr[i], qPoint[qPoint.size() - 1]))qPoint.pop_back(), qLine.pop_back(); 226 while (qLine.size() > 1 && !On_Left(arr[i], qPoint[0]))qPoint.pop_front(), qLine.pop_front(); 227 qLine.push_back(arr[i]); 228 if (qLine.size() > 1)qPoint.push_back(Intersection_line(qLine[qLine.size() - 1], qLine[qLine.size() - 2])); 229 } 230 } 231 while (qLine.size() > 1 && !On_Left(qLine[0], qPoint[qPoint.size() - 1]))qLine.pop_back(), qPoint.pop_back(); 232 while (qLine.size() > 1 && !On_Left(qLine[qLine.size() - 1], qPoint[0]))qLine.pop_front(), qPoint.pop_front(); 233 if (qPoint.size() < 2)return false; 234 qPoint.push_back(Intersection_line(qLine[0], qLine[qLine.size() - 1])); 235 for (int i = 0; i < qPoint.size(); ++i)v.push_back(qPoint[i]); 236 return true; 237 } 238 //------------------------------点与直线部分----------------------------- 239 240 //------------------------------圆与球部分------------------------------- 241 struct circle { 242 friend istream& operator >>(istream& in, circle& rhs) { 243 in >> rhs.c >> rhs.r; 244 return in; 245 } 246 friend ostream& operator <<(ostream& os, const circle& rhs) { 247 os << rhs.c << " " << rhs.r; 248 return os; 249 } 250 point c; 251 double r; 252 circle() {} 253 circle(point _c, double _r) :c(_c), r(_r) {} 254 //通过圆心角求圆上点 255 point Get_point(double a)const { return point(c.x + r * cos(a), c.y + r * sin(a)); } 256 }; 257 //求极角 258 inline double polar_angle(const point& rhs) { return atan2(rhs.y, rhs.x); } 259 //将角度转换为弧度 260 inline double torad(double deg) { return deg / 180. * pi; } 261 //解方程求直线与圆的交点 262 int GetLineCircleIntersection(const line& L, const circle& C, double& t1, double& t2, vector<point>& v) { 263 double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y - C.c.y; 264 double e = a * a + c * c, f = 2 * (a * b + c * d), g = b * b + d * d - C.r * C.r; 265 double delta = f * f - 4 * e * g; 266 if (Dcmp(delta) < 0)return 0; 267 if (Dcmp(delta) == 0) { 268 t1 = t2 = -f / (2 * e); 269 v.push_back(L.Get_Point(t1)); 270 return 1; 271 } 272 t1 = (-f - sqrt(delta)) / (2 * e); 273 t2 = (-f + sqrt(delta)) / (2 * e); 274 v.push_back(L.Get_Point(t1)); 275 v.push_back(L.Get_Point(t2)); 276 return 2; 277 } 278 //勾股定理求直线与圆的交点 279 int GetLineCircleIntersection(const point& A, const point& B, const circle&C, double& L, vector<point>& v) { 280 double d = DistanceToLine(C.c, A, B); 281 if (Dcmp(C.r - d) < 0)return 0; 282 point p = GetLineProjection(C.c, A, B); 283 double len = sqrt(C.r * C.r - d * d); 284 if (Dcmp(len) == 0) { 285 L = 0; 286 v.push_back(p); 287 return 1; 288 } 289 point unit = Unit_vector(B - A); 290 L = len; 291 v.push_back(p - unit * len); 292 v.push_back(p + unit * len); 293 return 2; 294 } 295 //求两圆交点 296 int GetCircleCircleIntersection(const circle& C1, const circle& C2, vector<point>& v) { 297 double d = length(C1.c - C2.c); 298 if (Dcmp(d) == 0) { 299 if (Dcmp(C1.r - C2.r) == 0)return -1;//两圆重合 300 return 0; 301 } 302 if (Dcmp(C1.r + C2.r - d) < 0 || Dcmp(fabs(C1.r - C1.r) - d) > 0)return 0; 303 double a = polar_angle(C2.c - C1.c); 304 double da = acos((C1.r * C1.r + d * d - C2.r * C2.r) / (2 * C1.r * d)); 305 point p1 = C1.Get_point(a - da), p2 = C1.Get_point(a + da); 306 v.push_back(p1); 307 if (p1 == p2)return 1; 308 v.push_back(p2); 309 return 2; 310 } 311 //求点到圆的切线 312 int GetTangents(const point& p, const circle& c, vector<point>& v) { 313 point u = c.c - p; 314 double dist = length(u); 315 if (dist < c.r)return 0; 316 else if (Dcmp(dist - c.r) == 0) { 317 v.push_back(Rotate(u, pi / 2)); 318 return 1; 319 } 320 else { 321 double ang = asin(c.r / dist); 322 v.push_back(Rotate(u, -ang)); 323 v.push_back(Rotate(u, +ang)); 324 return 2; 325 } 326 } 327 //求两圆的切线，返回切线个数 328 int GetTangents(circle A, circle B, vector<point>& a, vector<point>& b) { 329 if (Dcmp(A.r - B.r) < 0) { 330 swap(A, B); 331 swap(a, b); 332 } 333 double d2 = Dist(A.c, B.c), rdiff = A.r - B.r, rsum = A.r + B.r; 334 if (Dcmp(d2 - rdiff * rdiff) == -1)return 0;//内含 335 double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x); 336 if (Dcmp(d2) == 0 && Dcmp(A.r - B.r) == 0)return -1; 337 if (Dcmp(d2 - rdiff * rdiff) == 0) { 338 a.push_back(A.Get_point(base)); 339 b.push_back(B.Get_point(base)); 340 return 1; 341 }//内切 342 double ang = acos((A.r - B.r) / sqrt(d2)); 343 a.push_back(A.Get_point(base + ang)), a.push_back(A.Get_point(base - ang)); 344 b.push_back(B.Get_point(base + ang)), b.push_back(B.Get_point(base - ang)); 345 if (Dcmp(d2 - rsum * rsum) == 0) { 346 a.push_back(A.Get_point(base)); 347 b.push_back(B.Get_point(pi + base)); 348 } 349 else if (Dcmp(d2 - rsum * rsum) > 0) { 350 ang = acos((A.r + B.r) / sqrt(d2)); 351 a.push_back(A.Get_point(base + ang)), a.push_back(A.Get_point(base - ang)); 352 b.push_back(B.Get_point(pi + base + ang)), b.push_back(B.Get_point(pi + base - ang)); 353 } 354 return (int)a.size(); 355 } 356 //求三角形外接圆 357 circle CircumscribedCircle(const point& a, const point& b, const point& c) { 358 double Bx = b.x - a.x, By = b.y - a.y, Cx = c.x - a.x, Cy = c.y - a.y; 359 double D = 2 * (Bx * Cy - Cx * By); 360 double cx = (Cy * (Bx * Bx + By * By) - By *(Cx * Cx + Cy * Cy)) / D + a.x; 361 double cy = (Bx * (Cx * Cx + Cy * Cy) - Cx *(Bx * Bx + By * By)) / D + a.y; 362 point p = point(cx, cy); 363 return circle(p, length(a - p)); 364 } 365 //求三角形内接圆 366 circle InscribedCircle(const point& a, const point& b, const point& c) { 367 double lena = length(b - c), lenb = length(c - a), lenc = length(a - b); 368 point p = (a * lena + b * lenb + c * lenc) / (lena + lenb + lenc); 369 return circle(p, DistanceToLine(p, a, b)); 370 } 371 //------------------------------圆与球部分-------------------------------
posted on 2019-03-16 16:38 Eterna_King 阅读(222) 评论(0) 编辑 收藏 举报
