Fix line intersects line, line contains point. Add more tests.

math/include/game/math/compare.hpp

@@ -4,6 +4,18 @@
 #include <cmath>
 
 namespace math {
+  struct lessinf {
+    template <typename T>
+    bool operator()(T const & lhs, T const & rhs) {
+      return !isinf(lhs) && (isinf(rhs) || lhs < rhs);
+    }
+  };
+
+  template <typename T>
+  T safe_div(T num, T denom) {
+    return (num == denom && denom == 0) ? 0 : num / denom;
+  }
+
   template <typename T>
   bool between(T val, T min, T max) {
     return val >= min && (min == max || val < max);

math/src/common.cpp

@@ -33,17 +33,24 @@ namespace math {
     };
   }
   
+  enum orientation {
+    colinear = 0, clockwise = 1, anticlockwise = 2
+  };
+  
+  orientation orient(dim2::point const & p, dim2::point const & q,
+                     dim2::point const & r) {
+    float val = (q[1] - p[1]) * (r[0] - q[0]) - (q[0] - p[0]) * (r[1] - q[1]);
+    
+    if (val == 0) return colinear;
+    return (val > 0) ? clockwise : anticlockwise;
+  }
+  
   bool contains(dim2::line const & ln, dim2::point const & pt) {
-    if (ln.first[0] == ln.second[0]) {
-      return pt[0] == ln.first[0] && between(pt[1], ln.first[1], ln.second[1]);
-    } else if (ln.first == pt || ln.second == pt) {
-      return true;
-    } else {
-      auto slope = dim2::line{ln.first, pt}.slope();
-      return (slope == ln.slope() || slope == -1/ln.slope()) &&
-        between(pt[0], ln.first[0], ln.second[0]) &&
-        between(pt[1], ln.first[1], ln.second[1]);
-    }
+    auto xs = std::minmax(ln.first[0], ln.second[0], lessinf());
+    auto ys = std::minmax(ln.first[1], ln.second[1], lessinf());
+    return orient(ln.first, ln.second, pt) == colinear &&
+        between(pt[0], xs.first, xs.second) &&
+        between(pt[1], ys.first, ys.second);
   }
   
   bool contains(dim2::circle const & shape, dim2::point const & pt) {
@@ -51,18 +58,17 @@ namespace math {
     return delta.dot(delta) <= std::pow(shape.radius, 2);
   }
   
-  static dim2::line ray_x(dim2::point const & pt) {
-    return {pt, {{std::numeric_limits<float>::infinity(), pt[1]}}};
+  static dim2::line ray_x(dim2::point const & pt, dim2::line const & l) {
+    auto x_inf = std::max({l.first[0], l.second[0], pt[0]});
+    return {pt, {{x_inf, pt[1]}}};
   }
   
   static bool contains(std::vector<dim2::line> const & segments,
                        dim2::point const & pt) {
-    dim2::line const ray = ray_x(pt);
     int hits = 0;
     for (auto l : segments) {
-      dim2::point const inter = lines::intersection(l, ray);
-      if (inter == pt) return true;
-      else if (contains(l, inter)) ++hits;
+      if (contains(l, pt)) return true;
+      else if (intersects(l, ray_x(pt, l))) ++hits;
     }
     return (hits & 1) == 1;
   }
@@ -72,8 +78,28 @@ namespace math {
   }
   
   bool intersects(dim2::line const & lhs, dim2::line const & rhs) {
-    dim2::point const inter = lines::intersection(lhs, rhs);
-    return contains(rhs, inter);
+    // Find the four orientations needed for general and
+    // special cases
+    orientation o1 = orient(lhs.first, lhs.second, rhs.first);
+    orientation o2 = orient(lhs.first, lhs.second, rhs.second);
+    orientation o3 = orient(rhs.first, rhs.second, lhs.first);
+    orientation o4 = orient(rhs.first, rhs.second, lhs.second);
+    
+    // General case
+    if (o1 != o2 && o3 != o4)
+      return true;
+    
+    // Special Cases
+    // p1, q1 and p2 are colinear and p2 lies on segment p1q1
+    if (o1 == colinear && contains(lhs, rhs.second)) { return true; }
+    // p1, q1 and q2 are colinear and q2 lies on segment p1q1
+    if (o2 == colinear && contains(lhs, rhs.second)) { return true; }
+    // p2, q2 and p1 are colinear and p1 lies on segment p2q2
+    if (o3 == colinear && contains(rhs, lhs.first)) { return true; }
+    // p2, q2 and q1 are colinear and q1 lies on segment p2q2
+    if (o4 == colinear && contains(rhs, lhs.second)) { return true; }
+    
+    return false; // Doesn't fall in any of the above cases
   }
   
   bool intersects(dim2::line const & lhs, dim2::circle const & rhs) {

math/src/shape.cpp

@@ -10,13 +10,15 @@
 #include <iostream>
 #include <vector>
 
+#include "game/math/compare.hpp"
+
 namespace math { namespace dim2 {
   float line::length() const {
     return (second - first).magnitude();
   }
-  
+    
   float line::slope() const {
-    return (second[1] - first[1]) / (second[0] - first[0]);
+    return safe_div(second[1] - first[1], second[0] - first[0]);
   }
   
   rectangle::operator quad() const {
@@ -57,13 +59,16 @@ namespace math { namespace lines {
   }
   
   inline dim2::point intersection(float b1, float s1, float b2, float s2) {
-    float const x = (b1 + b2) / (s1 - s2);
+    float const x = safe_div(b1 + b2, s1 - s2);
     return {{x, b1 + x * s1}};
   }
 
-  inline dim2::point intersection(dim2::point p1, float s1, dim2::point p2,
-                                  float s2) {
-    return intersection(p1[1] - s1 * p1[0], s1, p2[1] - s2 * p2[0], s2);
+  inline dim2::point intersection(dim2::point const & p1, float s1,
+                                  dim2::point const & p2, float s2) {
+    // Solve for Inf * NaN
+    float y1 = p1[0] == 0 ? 0 : s1 * p1[0];
+    float y2 = p2[0] == 0 ? 0 : s2 * p2[0];
+    return intersection(p1[1] - y1, s1, p2[1] - y2, s2);
   }
   
   dim2::point intersection(dim2::line const & lhs, dim2::line const & rhs) {
@@ -72,6 +77,7 @@ namespace math { namespace lines {
 
   dim2::line orthogonal(dim2::line const & from, dim2::point const & to) {
     float const slope = from.slope();
+    if (slope == 0 || isinf(slope)) { return {to, {{to[0], from.first[1]}}}; }
     return {to, intersection(from.first, slope, to, -1/slope)};
   }
 } }

math/test/shape_test.cxx

@@ -53,6 +53,15 @@ TEST_P(XUnitTest, OrthoOnIntersection) {
               Eq(expected));
 }
 
+TEST_P(XUnitTest, OrthoOnIntersectionY) {
+  line const ln{{{0, 0}}, {{0, 1}}};
+  point const pt = GetParam();
+  line const expected{pt, {{pt[0], 0}}};
+  math::lines::orthogonal(ln, pt);
+  EXPECT_THAT(math::lines::orthogonal(ln, pt),
+              Eq(expected));
+}
+
 std::vector<point> x_orthos{
   {{0, 1}}, {{1, 1}}, {{1, 0}}, {{-1, 0}}, {{0, -1}}, {{-1, -1}},
   {{0, 2}}, {{2, 2}}, {{2, 0}}, {{-2, 0}}, {{0, -2}}, {{-2, -2}},