Advertisement
Mlxa

ALGO RayFloatGeometry

Mar 23rd, 2018
149
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 3.99 KB | None | 0 0
  1. #include <cmath>
  2. #include <cassert>
  3. #include <cstdio>
  4. #include <iostream>
  5. #include <iomanip>
  6. #include <vector>
  7. using namespace std;
  8.  
  9. #define all(x) begin(x), end(x)
  10. template<class A> void addlog(A a) { cerr << a << endl; }
  11. template<class A, class... B> void addlog(A a, B... b)
  12. { cerr << a << ' '; addlog(b...); }
  13.  
  14.  
  15. namespace floatGeometry {
  16.     // Types
  17.     #define pcout cout << fixed << setprecision(7)
  18.     typedef long long ll;
  19.     typedef long double ld;
  20.     const ld eps = 1e-12;
  21.     const ld inf = 1e+12;
  22.     bool eq (ld a, ld b) { return abs(a - b) < eps; }
  23.     ld abs (ld x) { return (x > 0? x : -x); }
  24.     struct Point { ld x, y; };
  25.     typedef Point Vector;
  26.    
  27.     // Constants
  28.     const Point     NO_POINT    {19, -inf},
  29.                     POINTS      {19, +inf};
  30.     // Point
  31.     ld abs (Point p) { return hypot(p.x, p.y); }
  32.     ld ang (Point p) { return atan2(p.y, p.x); }
  33.     bool operator== (Point a, Point b)
  34.     { return eq(a.x, b.x) && eq(a.y, b.y); }
  35.     bool operator!= (Point a, Point b)
  36.     { return !(a == b); }
  37.     // Vector
  38.     Vector operator- (Vector v)
  39.     { return Vector{ -v.x, -v.y }; }
  40.     Vector operator+ (Vector a, Vector b)
  41.     { return Vector{ a.x+b.x, a.y+b.y }; }
  42.     Vector operator- (Vector a, Vector b)
  43.     { return Vector{ a.x-b.x, a.y-b.y }; }
  44.     Vector operator* (Vector v, ld alpha)
  45.     { return Vector{ alpha * v.x, alpha * v.y }; }
  46.     ld operator* (Vector a, Vector b)
  47.     { return a.x*b.y - b.x*a.y; }
  48.     ld operator% (Vector a, Vector b)
  49.     { return a.x*b.x + a.y*b.y; }
  50.     Vector turn90 (Vector v)
  51.     { return Vector{v.y, -v.x}; }
  52.     // I/O.
  53.     istream& operator>> (istream& s, Point &p)
  54.     { return s >> p.x >> p.y; }
  55.     ostream& operator<< (ostream &s, Point p)
  56.     { return s << p.x << ' ' << p.y; }
  57.    
  58.    
  59.     struct Segment { Point p, q; };
  60.     ld abs (Segment s) { return abs(s.p - s.q); }
  61.    
  62.     istream& operator>> (istream &s, Segment &x)
  63.     { return s >> x.p >> x.q; }
  64.     ostream& operator<< (ostream &s, Segment x)
  65.     { return s << x.p << "  " << x.q; }
  66.    
  67.     struct Line { ld a, b, c; Line () {}
  68.         void normalize () {
  69.             ld z = hypot(a, b);
  70.             if (eq(a, 0)) { if (b < 0) z = -z; }
  71.             else if (a < 0) z = -z;
  72.             a /= z; b /= z; c /= z;
  73.                                                                         assert(eq(hypot(a, b), 1));
  74.         } Line (Point p, Point q) {
  75.             a = p.y - q.y; b = q.x - p.x;
  76.             c = - (a * p.x + b * p.y );
  77.             normalize();
  78.         } Line (Segment s) : Line(s.p, s.q) {}
  79.     };
  80.    
  81.     // Line
  82.     ld vMul (ld ax, ld ay, ld bx, ld by) { return ax*by - bx*ay; }
  83.     Point intersect (Line I, Line II) {
  84.         ld D0 = vMul(I.a, I.b, II.a, II.b);
  85.         ld Dx = vMul(I.b, I.c, II.b, II.c);
  86.         ld Dy = vMul(I.c, I.a, II.c, II.a);
  87.         if (eq(D0, 0)) {
  88.             assert(eq(Dx, Dy));
  89.             if (eq(Dx, 0))  return POINTS;
  90.             else            return NO_POINT;
  91.         } return Point{Dx / D0, Dy / D0};
  92.     }
  93.     // I/O
  94.     istream& operator>> (istream& s, Line &l)
  95.     { return s >> l.a >> l.b >> l.c; }
  96.     ostream& operator<< (ostream& s, Line l)
  97.     { return s << l.a << ' ' << l.b << ' ' << l.c; }
  98.    
  99.     struct Ray {
  100.         Point p; Vector v; Line l;
  101.         Ray (Point P, Vector V) : p(P), v(V), l(p, p + v) {}
  102.     };
  103.    
  104.     ld signDist (Line l, Point p) {
  105.         l.normalize();
  106.         return l.a * p.x + l.b * p.y + l.c;
  107.     } ld dist (Line l, Point p) { return abs(signDist(l, p)); }
  108.    
  109.     ld dist (Ray r, Point p) {
  110.         if ( (p - r.p) % r.v > -eps )
  111.             return dist(r.l, p);
  112.         else return abs(p - r.p);
  113.     }
  114.    
  115.     ld dist (Segment s, Point p) {
  116.         Ray r1 (s.p, s.q - s.p),
  117.             r2 (s.q, s.p - s.q);
  118.         return max(dist(r1, p), dist(r2, p));
  119.     }
  120.    
  121.     template <class G> bool
  122.     has (G g, Point p)
  123.     { return eq(dist(g, p), 0); }
  124.    
  125.     Point intersect (Segment a, Segment b) {
  126.         Line la(a), lb(b);
  127.         Point p = intersect(la, lb);
  128.         if (p == POINTS) {
  129.             // Return one of that points.
  130.             if (has(a, b.p)) return b.p;
  131.             if (has(a, b.q)) return b.q;
  132.             if (has(b, a.p)) return a.p;
  133.             if (has(b, a.q)) return a.q;
  134.                                                                         assert(false);
  135.         } else if (p != NO_POINT && has(a, p) && has(b, p)) {
  136.             return p;
  137.         } else
  138.             return NO_POINT;
  139.     }
  140. }
  141.  
  142.  
  143.  
  144.  
  145.  
  146.  
  147.  
  148. using namespace floatGeometry;
  149. int main () {
  150.    
  151.     for (Segment a, b; cin >> a >> b; )
  152.         pcout << intersect(a, b) << endl;
  153.    
  154.     return 0;
  155. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement