| 2 | #define L(p0, p1) P(p0), P(p1) |
| 3 | #define PL(p0, p1, p2) P(p0), P(p1), P(p2) |
| 4 | struct point3d { |
| 5 | double x, y, z; |
| 6 | point3d() : x(0), y(0), z(0) {} |
| 7 | point3d(double _x, double _y, double _z) |
| 8 | : x(_x), y(_y), z(_z) {} |
| 9 | point3d operator+(P(p)) const { |
| 10 | return point3d(x + p.x, y + p.y, z + p.z); } |
| 11 | point3d operator-(P(p)) const { |
| 12 | return point3d(x - p.x, y - p.y, z - p.z); } |
| 13 | point3d operator-() const { |
| 14 | return point3d(-x, -y, -z); } |
| 15 | point3d operator*(double k) const { |
| 16 | return point3d(x * k, y * k, z * k); } |
| 17 | point3d operator/(double k) const { |
| 18 | return point3d(x / k, y / k, z / k); } |
| 19 | double operator%(P(p)) const { |
| 20 | return x * p.x + y * p.y + z * p.z; } |
| 21 | point3d operator*(P(p)) const { |
| 22 | return point3d(y*p.z - z*p.y, |
| 23 | z*p.x - x*p.z, x*p.y - y*p.x); } |
| 24 | double length() const { |
| 25 | return sqrt(*this % *this); } |
| 26 | double distTo(P(p)) const { |
| 27 | return (*this - p).length(); } |
| 28 | double distTo(P(A), P(B)) const { |
| 29 | // A and B must be two different points |
| 30 | return ((*this - A) * (*this - B)).length() / A.distTo(B);} |
| 31 | double signedDistTo(PL(A,B,C)) const { |
| 32 | // A, B and C must not be collinear |
| 33 | point3d N = (B-A)*(C-A); double D = A%N; |
| 34 | return ((*this)%N - D)/N.length(); } |
| 35 | point3d normalize(double k = 1) const { |
| 36 | // length() must not return 0 |
| 37 | return (*this) * (k / length()); } |
| 38 | point3d getProjection(P(A), P(B)) const { |
| 39 | point3d v = B - A; |
| 40 | return A + v.normalize((v % (*this - A)) / v.length()); } |
| 41 | point3d rotate(P(normal)) const { |
| 42 | //normal must have length 1 and be orthogonal to the vector |
| 43 | return (*this) * normal; } |
| 44 | point3d rotate(double alpha, P(normal)) const { |
| 45 | return (*this) * cos(alpha) + rotate(normal) * sin(alpha);} |
| 46 | point3d rotatePoint(P(O), P(axe), double alpha) const{ |
| 47 | point3d Z = axe.normalize(axe % (*this - O)); |
| 48 | return O + Z + (*this - O - Z).rotate(alpha, O); } |
| 49 | bool isZero() const { |
| 50 | return abs(x) < EPS && abs(y) < EPS && abs(z) < EPS; } |
| 51 | bool isOnLine(L(A, B)) const { |
| 52 | return ((A - *this) * (B - *this)).isZero(); } |
| 53 | bool isInSegment(L(A, B)) const { |
| 54 | return isOnLine(A, B) && ((A - *this) % (B - *this))<EPS;} |
| 55 | bool isInSegmentStrictly(L(A, B)) const { |
| 56 | return isOnLine(A, B) && ((A - *this) % (B - *this))<-EPS;} |
| 57 | double getAngle() const { |
| 58 | return atan2(y, x); } |
| 59 | double getAngle(P(u)) const { |
| 60 | return atan2((*this * u).length(), *this % u); } |
| 61 | bool isOnPlane(PL(A, B, C)) const { |