peltorator

!Стереома Рубинчика

Dec 8th, 2017
79
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <cmath>
  2. #include <iostream>
  3. #include <cstdio>
  4.  
  5.  
  6. using namespace std;
  7.  
  8.  
  9. double sqr(double a) // квадрат числа
  10. {
  11.         return a * a;
  12. }
  13.  
  14.  
  15. bool doubleEqual(double a, double b) // сравниваем на равенство с eps
  16. {
  17.         return fabs(a - b) < 1e-9;
  18. }
  19.  
  20.  
  21. bool doubleLessOrEqual(double a, double b) // <= с eps
  22. {
  23.         return a < b || doubleEqual(a, b);
  24. }
  25.  
  26.  
  27. bool doubleLess(double a, double b) // < с eps
  28. {
  29.         return a < b && !doubleEqual(a, b);
  30. }
  31.  
  32.  
  33. bool doubleGreaterOrEqual(double a, double b) // >= с eps
  34. {
  35.         return a > b || doubleEqual(a, b);
  36. }
  37.  
  38.  
  39. bool doubleGreater(double a, double b) // > с eps
  40. {
  41.         return a > b && !doubleEqual(a, b);
  42. }
  43.  
  44.  
  45. double mySqrt(double a) // mySqrt с проверкой, что корректен аргумент
  46. {
  47.         if(doubleLess(a, 0) )
  48.         {
  49.                 throw "sqrt(-1)";
  50.         }
  51.         if(a < 0) return 0;
  52.         return sqrt(a);
  53. }
  54.  
  55.  
  56. struct Point{ // класс точки или вектора, далее мы эти понятия разделять не будем
  57.                                 // но условимся (для удобного чтения) большими буквами обозначать точки (A, B, C, D, ...)
  58.                                 // маленькими - вектора(v, u, w,...)
  59. private: double x, y, z; // 3 поля, других не будет
  60. public:
  61.         Point(): x(0), y(0), z(0) {} // конструктор по умолчанию
  62.  
  63.         Point(double x, double y, double z): x(x), y(y), z(z) {} // намеренно сделаем 2 конструктора вместо Point(x = 0...)
  64.                                                                                         //Это изавит нас от ошибок вида Point A = 2;
  65.  
  66.         void scan() // читаем координаты точки
  67.         {
  68.                 scanf("%lf %lf %lf", &x, &y, &z);
  69.         }
  70.  
  71.         void print() const // выводим координаты точки
  72.         {
  73.                 printf("%.10lf %.10lf %.10lf\n", x, y, z);
  74.         }
  75.  
  76.         Point operator+(const Point & p) const // сложение 2х векторов
  77.         {
  78.                 return Point(x + p.x, y + p.y, z + p.z);
  79.         }
  80.  
  81.         Point operator-(const Point & p) const // вычитание 2х векторов
  82.         {
  83.                 return Point(x - p.x, y - p.y, z - p.z);
  84.         }
  85.  
  86.         Point operator-() const // вычитание 2х векторов
  87.         {
  88.                 return Point(-x, -y, -z);
  89.         }
  90.  
  91.         Point operator*(double k) const // умножение вектора на скаляр
  92.         {
  93.                 return Point(x * k, y * k, z * k);
  94.         }
  95.  
  96.         Point operator/(double k) const // деление вектора на скаляр
  97.         {
  98.                 return Point(x / k, y / k, z / k);
  99.         }
  100.  
  101.         double operator%(const Point & p) const // скалярное произведение
  102.         {
  103.                 return x * p.x + y * p.y + z * p.z;
  104.         }
  105.  
  106.         Point operator*(const Point & p) const // векторное произведение
  107.         {
  108.                 return
  109.                         Point
  110.                         (
  111.                                  y * p.z - z * p.y,
  112.                                  z * p.x - x * p.z,
  113.                                  x * p.y - y * p.x
  114.                         );
  115.         }
  116.  
  117.         double length() const // длина вектора по определению из ан.гема (корень из скалярного квадрата)
  118.         {
  119.                 return mySqrt(*this % *this);
  120.         }
  121.  
  122.         double distTo(const Point & p) const //расстояние между 2мя точками - модуль вектора между ними
  123.         {
  124.                 return (*this - p).length();
  125.         }
  126.  
  127.         double distTo(const Point & A, const Point & B) const // расстояние от точки до прямой (всегда >= 0)
  128.         {
  129.                 double d = A.distTo(B);
  130.                 if(doubleEqual(d, 0) ) // прямая должна задаваться 2мя точками
  131.                 {
  132.                         throw "A = B";
  133.                 }
  134.                 double s = ( (*this - A) * (*this - B) ).length(); // площадь треугольника
  135.                 return s / d; // метод площадей
  136.         }
  137.  
  138.         Point normalize(double k = 1) const // выставляет длину вектора в k
  139.         {
  140.                 double len = length(); // Текущая длина
  141.                 if(doubleEqual(len, 0) )
  142.                 {
  143.                         if(doubleEqual(k, 0) )
  144.                         {
  145.                                 return Point();
  146.                         }
  147.                         throw "zero-size vector";
  148.                 }
  149.                 return *this * (k / len);
  150.         }
  151.  
  152.         Point getH(const Point & A, const Point & B) const // опускаем высоту из точки на прямую (A, B)
  153.         {
  154.                 Point C = *this;
  155.                 Point v = B - A; // направляющий вектор прямой
  156.                 Point u = C - A; // вектор, проекция которого нам нужна
  157.                 double k = v % u / v.length(); // нашли длину проекции
  158.                 v = v.normalize(k); // нашли проекцию u на v
  159.                 Point H = A + v; // подвинули точку A в проекцию конца вектора u, отложенного из A
  160.                 return H;
  161.         }
  162.  
  163.         Point rotate(Point normal) const // поворот на 90 градусов против часовой стрелки (положительное направление)
  164.         {
  165.                 Point v = *this;
  166.                 if(!doubleEqual(v % normal, 0) || !doubleEqual( normal.length(), 1) )
  167.                 {
  168.                         throw "no axe";
  169.                 }
  170.                 return v * normal;
  171.         }
  172.  
  173.         Point rotate(double alpha, const Point & normal) const // поворот на угол alpha против часовой стрелки
  174.                                                                         // (по часовой стрелке, если alpha < 0)
  175.         {
  176.                 return rotate(cos(alpha), sin(alpha), normal ); // делегируем задачу другому экземпляру функции
  177.         }
  178.  
  179.         Point rotate(double cosa, double sina, const Point & normal) const // поворот с заданными косинусом и синусом
  180.         {
  181.                 Point v = *this;
  182.                 Point u = v.rotate(normal); // вектор, перпендикулярный v, теперь (v, u) - базис, в котором мы знаем ответ
  183.                 Point w = v * cosa + u * sina; // зная координаты в базисе (v, u), нашли w
  184.                 return w;
  185.         }
  186.  
  187.         bool isZero() const // проверка на то, что точка нулевая без сложных операций
  188.         {
  189.                 return doubleEqual(x, 0) && doubleEqual(y, 0) && doubleEqual(z, 0);
  190.         }
  191.  
  192.         bool isOnLine(const Point & A, const Point & B) const // точка на прямой?
  193.         {
  194.                 return ( (A - *this) * (B - *this) ).isZero();
  195.         }
  196.  
  197.         bool isInSegment(const Point & A, const Point & B) const // точка внутри отрезка?
  198.         {
  199.                 return isOnLine(A, B) && doubleLessOrEqual( (A - *this) % (B - *this), 0 );
  200.         }
  201.  
  202.         bool isInSegmentStrictly(const Point & A, const Point & B) const // точка внутри отрезка строго?
  203.         {
  204.                 return isOnLine(A, B) && doubleLess( (A - *this) % (B - *this), 0 );
  205.         }
  206.  
  207.         double getAngle() const
  208.         {
  209.                 return atan2(y, x); // угол между вектором и осью ОХ
  210.         }
  211.  
  212.         double getAngle(Point u) const
  213.         {
  214.                 Point v = *this;
  215.                 return atan2((v * u).length() , v % u); // угол между векторами (Не ориентированный)
  216.         }
  217.  
  218.         bool isOnPlane(const Point & A, const Point & B, const Point & C) const // проверка принадлежности точки к плоскости
  219.         {
  220.                 return doubleEqual( (A - *this) * (B - *this) % (C - *this), 0);
  221.         }
  222.  
  223. };
  224.  
  225.  
  226. int getIntersection // ищем пересечение прямой (A, B) и прямой (C, D)
  227.         (
  228.                 const Point & A,
  229.                 const Point & B,
  230.                 const Point & C,
  231.                 const Point & D,
  232.                 Point & O
  233.         )
  234. {
  235.         if( !doubleEqual( (B - A) * (C - A) % (D - A), 0)  )
  236.         {
  237.                 throw "It's not plane";
  238.         }
  239.         if( doubleEqual( ( (A - B) * (C - D) ).length(), 0) )
  240.         {
  241.                 if(A.isOnLine(C, D) ) return 2;
  242.                 return 0;
  243.         }
  244.         Point normal = ( (A - B) * (C - B) ).normalize();
  245.         Point v = B - A; // направляющий вектор прямой (A, B)
  246.         double s1 = (C - A) * (D - A) % normal; // площадь треугольника A, C, D
  247.         double s2 = (D - B) * (C - B) % normal; // площадь треугольника B, D, C
  248.         double s = s1 + s2; // площадь четурёхугольника
  249.         v = v / s;
  250.         v = v * s1; // вектора AO и AB пропорциональны площадям s1 и s
  251.         O = A + v; // нашли точку пересечения
  252.         return 1; // 1 точка в пересечении
  253. }
  254.  
  255.  
  256.  
  257. int getIntersection // ищем пересечение прямой (A, B) и плоскости (C, D, E)
  258.         (
  259.                 const Point & A,
  260.                 const Point & B,
  261.                 const Point & C,
  262.                 const Point & D,
  263.                 const Point & E,
  264.                 Point & O
  265.         )
  266. {
  267.         Point v = B - A; // направляющий вектор прямой (A, B)
  268.         double V1 = (C - A) * (D - A) % (E - A); // объём тетраэдра A, C, D, E
  269.         double V2 = (D - B) * (C - B) % (E - B); // объём тетраэдра B, D, C, E
  270.         double V = V1 + V2; // площадь четурёхугольника
  271.         v = v / V;
  272.         if(doubleEqual(V, 0) )
  273.         {
  274.                 if(A.isOnPlane(C, D, E) ) return 2;
  275.                 return 0;
  276.         }
  277.         v = v * V1; // вектора AO и AB пропорциональны площадям s1 и s
  278.         O = A + v; // нашли точку пересечения
  279.         return 1; // 1 точка в пересечении
  280. }
  281.  
  282.  
  283. bool getIntersection // ищем пересечение плоскости (A, nA) и плоскости (B, nB). Плоскость даём точкой и нормалью
  284.         (
  285.                 const Point & A,
  286.                 const Point & nA,
  287.                 const Point & B,
  288.                 const Point & nB,
  289.                 Point & P,
  290.                 Point & Q
  291.         )
  292. {
  293.         Point n = nA * nB;
  294.         if(n.isZero() ) return false;
  295.         Point v = n * nA;
  296.         double k = (B - A) % nB / (v % nB);
  297.        
  298.         v = v * k;
  299.         P = A + v;
  300.         Q = P + n;
  301.         return true;
  302. }
RAW Paste Data