计算几何.docx

1、计算几何/计算几何(二维) #include #include #include #include using namespace std;const double EPSILON = 1e-8;const double INFINITY = 1e10;const double PI = 3.14159265358979323846;inline double DegToRad(double deg)return (deg * PI / 180.0);inline double RadToDeg(double rad)return (rad * 180.0 / PI);inline doubl

2、e Sin(double deg)return sin(DegToRad(deg);inline double Cos(double deg)return cos(DegToRad(deg);inline double Tan(double deg)return tan(DegToRad(deg);inline double ArcSin(double val)return RadToDeg(asin(val);inline double ArcCos(double val)return RadToDeg(acos(val);inline double ArcTan(double val)re

3、turn RadToDeg(atan(val);struct POINT double x, y; POINT(); POINT(double _x, double _y) x = _x; y = _y; bool operator = (const POINT& _P) const return fabs(x - _P.x) = EPSILON & fabs(y - _P.y) = EPSILON; bool operator != (const POINT & _P) const return !(fabs(x - _P.x) = EPSILON & fabs(y - _P.y) = EP

4、SILON); bool operator (const POINT& _P) const if(y != _P.y) return y _P.y; else return x _P.x; ; POINT operator + (const POINT& _P) const return POINT(x + _P.x, y + _P.y); POINT operator - (const POINT& _P) const return POINT(x - _P.x, y - _P.y); double Angle(const POINT& O) /返回点关于O点的极坐标 double x1 =

5、 x - O.x; double y1 = y - O.y; if(fabs(x1) 0.0) return 90.0; else return 270.0; if(fabs(y1) 0.0) return 0.0; else return 180.0; double angle = ArcTan(fabs(y1 / x1); if(x1 0) return y1 0 ? angle : 360.0 - angle; else return y1 0 ? 180.0 - angle : 180.0 + angle; double Angle() /默认O点为(0，0) return Angle

6、(POINT(0,0); ;inline double Distance(const POINT & a, const POINT & b) return sqrt(a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y); struct SEGMENT POINT a, b; SEGMENT(); SEGMENT(POINT _a, POINT _b) a = _a; b = _b; ; bool operator = (const SEGMENT& _SEGMENT) const return a = _SEGMENT.a & b = _SE

7、GMENT.b; ;/Ax + By + C = 0struct LINE POINT a, b; double A, B, C; LINE(POINT _a, POINT _b) a = _a; b = _b; A = b.y - a.y; B = a.x - b.x; C = b.x * a.y - a.x * b.y; ; LINE(double _A, double _B, double _C) A = _A; B = _B; C = _C; ; bool operator=(const LINE& _L) const return fabs(A * _L.B - _L.A * B)

8、= EPSILON & fabs(A * _L.C - _L.A * C) = EPSILON & fabs(B * _L.C - _L.B * C) EPSILON) return 1; else if(t EPSILON) return 1; else if(t -EPSILON) return -1; else return 0;/点到直线的垂足POINT Pedal_Point(const POINT& P, const LINE& L) if(fabs(L.B) = EPSILON) /如果直线平行Y轴 return POINT(-L.C / L.A, P.y); else if(f

9、abs(L.A) = EPSILON) /如果直线平行X轴 return POINT(P.x, -L.C / L.B); else double k1 = - L.A / L.B; double k2 = - 1 / k1; double b1 = - L.C / L.B; double b2 = P.y - k2 * P.x; double x = (b2 - b1) / (k1 - k2); double y = k2 * x + b2; return POINT(x, y); /p点关于直线L的对称点POINT Symmetrical_Point_Of_Line(POINT P, LIN

10、E L) POINT _P; double d; d = L.A * L.A + L.B * L.B; _P.x = (L.B * L.B * P.x - L.A * L.A * P.x - 2.0 * L.A * L.B * P.y - 2.0 * L.A * L.C) / d; _P.y = (L.A * L.A * P.y - L.B * L.B * P.y - 2.0 * L.A * L.B * P.x - 2.0 * L.B * L.C) / d; return _P;bool Is_On_Segment(const POINT & P, const SEGMENT & S) if(

11、P = S.a | P = S.b) return true; else if(P.x min(S.a.x, S.b.x) - EPSILON & P.y min(S.a.y, S.b.y) - EPSILON & fabs(Cross(S.b, P, S.a) = EPSILON) return true; return false;/点到线段的最近点POINT Nearest_Point_To_Segment(const POINT& P, const SEGMENT& S) LINE L(S.a, S.b); POINT I = Pedal_Point(P, L); if(Is_On_S

12、egment(I, S) return I; else double t1 = Distance(P, S.a); double t2 = Distance(P, S.b); return t1 min(v.a.x, v.b.x) - EPSILON) & (max(v.a.x, v.b.x) min(u.a.x, u.b.x) - EPSILON) & (max(u.a.y, u.b.y) min(v.a.y, v.b.y) - EPSILON) & (max(v.a.y, v.b.y) min(u.a.y, u.b.y) - EPSILON) int t1 = Sign_Cross(v.a

13、, u.b, u.a); int t2 = Sign_Cross(u.b, v.b, u.a); int t3 = Sign_Cross(u.a, v.b, v.a); int t4 = Sign_Cross(v.b, u.b, v.a); return (t1 * t2 0) & (t3 * t4 0); return false;bool Is_Touch(const SEGMENT &u, const SEGMENT & v) if( (max(u.a.x, u.b.x) = min(v.a.x, v.b.x) - EPSILON) & (max(v.a.x, v.b.x) = min(

14、u.a.x, u.b.x) - EPSILON) & (max(u.a.y, u.b.y) = min(v.a.y, v.b.y) - EPSILON) & (max(v.a.y, v.b.y) = min(u.a.y, u.b.y) - EPSILON) int t1 = Sign_Cross(v.a, u.b, u.a); int t2 = Sign_Cross(u.b, v.b, u.a); int t3 = Sign_Cross(u.a, v.b, v.a); int t4 = Sign_Cross(v.b, u.b, v.a); if(t1 * t2 = 0) return (t3

15、* t4 = 0); else return (t1 * t2 0 & t3 * t4 = 0); return false;/判断线断P和直线Q是否相交，端点重合算相交bool Is_Intersect(const SEGMENT & p, const LINE & q) SEGMENT p1,p2,q0; p1.a=q.a;p1.b=p.a; p2.a=q.a;p2.b=p.b; q0.a=q.a;q0.b=q.b; int t1 = Sign_Cross(p1, q0); int t2 = Sign_Cross(q0, p2); return ( t1 * t2 = 0); /判断两条直

16、线是否平行bool Is_Parallel(const LINE & L1, const LINE & L2) /共线不算平行/ return (L1.A * L2.B = L2.A * L1.B) & / (L1.A * L2.C != L2.A * L1.C) | (L1.B * L2.C != L2.B * L1.C); /共线算平行 return (fabs(L1.A * L2.B - L2.A * L1.B) = EPSILON); /判断两条直线是否相交 bool Is_Intersect(const LINE & L1, const LINE & L2) return !(L1

17、= L2) & !(Is_Parallel(L1, L2);/角(AOB)的平分线LINE Angle_Bisector(POINT A, const POINT& O, POINT B) double angle1 = A.Angle(O); double angle2 = B.Angle(O); double angle3 = (angle2 - angle1) / 2.0 + angle1; if(fabs(angle3 - 90.0) = EPSILON | fabs(angle3 - 270.0) = EPSILON) POINT P1(O.x, O.y + 1.0); return

18、 LINE(P1, O); double k = Tan(angle3); POINT P2(O.x + 1.0, O.y + k); return LINE(P2, O);/线段的垂直平分线LINE SEGMENT_Bisector(const SEGMENT& S) POINT O(S.a.x + S.b.x) / 2.0, (S.a.y + S.b.y) / 2.0); return Angle_Bisector(S.a, O, S.b);/从X轴旋转到S的角度double Angle(const SEGMENT& S) double x1 = S.b.x - S.a.x; double

19、 y1 = S.b.y - S.a.y; if(fabs(x1) 0.0) return 0.5 * PI; else return 1.5 * PI; if(fabs(y1) 0.0) return 0.0; else return PI; double angle = atan(fabs(y1 / x1); if(x1 0) return y1 0 ? angle : 2.0 * PI - angle; else return y1 0 ? PI - angle : PI + angle; /线段s1与s2的夹角double Angle_Between(const SEGMENT& s1,

20、 const SEGMENT& s2) /注意误差，若太接近，就另为0 double angle = Angle(s2) - Angle(s1); if(fabs(angle) 2.0 * PI) angle = 0.0; else if(angle 0.0) angle += 2.0 * PI; return angle;/返回两线段的重合部分vector Coincidence_With_Segment(const SEGMENT & s1, const SEGMENT & s2) vector S; SEGMENT s; if(Is_Touch(s1, s2) if(Is_On_Segm

21、ent(s1.a, s2) & Is_On_Segment(s1.b, s2) S.push_back(s1); else if(Is_On_Segment(s2.a, s1) & Is_On_Segment(s2.b, s1) S.push_back(s2); else if(Is_On_Segment(s2.b, s1) & Is_On_Segment(s1.b, s2) & s2.b != s1.b) s.a = s2.b; s.b = s1.b; S.push_back(s); else if(Is_On_Segment(s2.a, s1) & Is_On_Segment(s1.a,

22、s2) & s2.a != s1.a) s.a = s2.a; s.b = s1.a; S.push_back(s); else if(Is_On_Segment(s1.a, s2) s.a = s1.a; s.b = s1.a; S.push_back(s); else if(Is_On_Segment(s2.a, s1) s.a = s2.a; s.b = s2.a; S.push_back(s); return S;/ 圆/struct CIRCLE POINT C; double R; CIRCLE() CIRCLE(POINT _P, double _R) C = _P; R = _

23、R; bool operator=(const CIRCLE& _C) const return fabs(R - _C.R) B.R) ? A : B; CIRCLE N = (A.R B.R) ? B : A; double D = Distance(M.C, N.C); if (D M.R - N.R) double cosM = (M.R * M.R + D * D - N.R * N.R) / (2.0 * M.R * D); double cosN = (N.R * N.R + D * D - M.R * M.R) / (2.0 * N.R * D); double alpha = 2.0 * ArcCos(cosM); double beta = 2.0 * ArcCos(cosN); double TM = 0.5 * M.R * M.R * Sin(alpha); double TN = 0.5 * N.R * N.R * Sin(beta); double FM = (alpha / 360.0) * M.Area(); double FN = (beta / 360.0) * N.Area(); area = FM + FN - TM - TN; else if (D = M.R - N.R) area = N.Ar