From 49c3ecd7c6f052ff94aaad45ea608178d1712585 Mon Sep 17 00:00:00 2001
From: Alex Weiss <algrs@cacography.net>
Date: Fri, 30 Jun 2017 10:08:25 -0400
Subject: [PATCH] vectorize line intersection (#7)

* vectorize line intersection

* fix lint

* bump versiom
---
 blmath/__init__.py                       |  2 +-
 blmath/geometry/primitives/plane.py      | 26 ++++++++++++++++++++----
 blmath/geometry/primitives/test_plane.py | 26 ++++++++++++++++++++++++
 3 files changed, 49 insertions(+), 5 deletions(-)

diff --git a/blmath/__init__.py b/blmath/__init__.py
index 3f262a6..923b987 100644
--- a/blmath/__init__.py
+++ b/blmath/__init__.py
@@ -1 +1 @@
-__version__ = '1.2.1'
+__version__ = '1.2.2'
diff --git a/blmath/geometry/primitives/plane.py b/blmath/geometry/primitives/plane.py
index c19a060..e1952b7 100644
--- a/blmath/geometry/primitives/plane.py
+++ b/blmath/geometry/primitives/plane.py
@@ -229,6 +229,22 @@ def _line_segment_xsection(self, a, b):
                 return None
         return pt
 
+    def line_xsections(self, pts, rays):
+        denoms = np.dot(rays, self.normal)
+        denom_is_zero = denoms == 0
+        denoms[denom_is_zero] = np.nan
+        p = np.dot(self.reference_point - pts, self.normal) / denoms
+        return np.vstack([p, p, p]).T * rays + pts, ~denom_is_zero
+
+    def line_segment_xsections(self, a, b):
+        pts, pt_is_valid = self.line_xsections(a, b-a)
+        pt_is_out_of_bounds = np.logical_or(np.any(np.logical_and(pts[pt_is_valid] > a[pt_is_valid], pts[pt_is_valid] > b[pt_is_valid]), axis=1),
+                                            np.any(np.logical_and(pts[pt_is_valid] < a[pt_is_valid], pts[pt_is_valid] < b[pt_is_valid]), axis=1))
+        pt_is_valid[pt_is_valid] = ~pt_is_out_of_bounds
+        pts[~pt_is_valid] = np.nan
+        return pts, pt_is_valid
+
+
     def mesh_xsection(self, m, neighborhood=None):
         '''
         Backwards compatible.
@@ -309,11 +325,13 @@ def get(self, u, v):
                     return None
 
         intersection_map = EdgeMap()
-        for e in es:
+
+        pts, pt_is_valid = self.line_segment_xsections(m.v[es[:, 0]], m.v[es[:, 1]])
+        valid_pts = pts[pt_is_valid]
+        valid_es = es[pt_is_valid]
+        for val, e in zip(valid_pts, valid_es):
             if not intersection_map.contains(e[0], e[1]):
-                val = self._line_segment_xsection(m.v[e[0]], m.v[e[1]])
-                if val is not None:
-                    intersection_map.add(e[0], e[1], val)
+                intersection_map.add(e[0], e[1], val)
         verts = np.array(intersection_map.values)
 
         class Graph(object):
diff --git a/blmath/geometry/primitives/test_plane.py b/blmath/geometry/primitives/test_plane.py
index 3bdda4f..dfeda0c 100644
--- a/blmath/geometry/primitives/test_plane.py
+++ b/blmath/geometry/primitives/test_plane.py
@@ -245,6 +245,19 @@ def test_line_plane_intersection(self):
         np.testing.assert_array_equal(plane.line_xsection(pt=[0., -1., 0.], ray=[0., 1., 0.]), [0., 0., 0.])
         np.testing.assert_array_equal(plane.line_xsection(pt=[0., -1., 0.], ray=[1., 1., 0.]), [1., 0., 0.])
 
+    def test_line_plane_intersections(self):
+        # x-z plane
+        normal = np.array([0., 1., 0.])
+        sample = np.array([0., 0., 0.])
+
+        plane = Plane(sample, normal)
+        pts = np.array([[0., -1., 0.], [0., 0., 0.], [0., -1., 0.], [0., -1., 0.]])
+        rays = np.array([[1., 0., 0.], [1., 0., 0.], [0., 1., 0.], [1., 1., 0.]])
+        expected = np.array([[np.nan, np.nan, np.nan], [np.nan, np.nan, np.nan], [0., 0., 0.], [1., 0., 0.]])
+        intersections, is_intserseting = plane.line_xsections(pts, rays)
+        np.testing.assert_array_equal(intersections, expected)
+        np.testing.assert_array_equal(is_intserseting, [False, False, True, True])
+
     def test_line_segment_plane_intersection(self):
         # x-z plane
         normal = np.array([0., 1., 0.])
@@ -257,6 +270,19 @@ def test_line_segment_plane_intersection(self):
         np.testing.assert_array_equal(plane.line_segment_xsection([0., -1., 0.], [2., 1., 0.]), [1., 0., 0.])
         self.assertIsNone(plane.line_segment_xsection([0., 1., 0.], [0., 2., 0.])) # line intersecting, but not in segment
 
+    def test_line_segment_plane_intersections(self):
+        # x-z plane
+        normal = np.array([0., 1., 0.])
+        sample = np.array([0., 0., 0.])
+
+        plane = Plane(sample, normal)
+        a = np.array([[0., -1., 0.], [0., 0., 0.], [0., -1., 0.], [0., -1., 0.], [0., 1., 0.]])
+        b = np.array([[1., -1., 0.], [1., 0., 0.], [0., 1., 0.], [2., 1., 0.], [0., 2., 0.]])
+        expected = np.array([[np.nan, np.nan, np.nan], [np.nan, np.nan, np.nan], [0., 0., 0.], [1., 0., 0.], [np.nan, np.nan, np.nan]])
+        intersections, is_intserseting = plane.line_segment_xsections(a, b)
+        np.testing.assert_array_equal(intersections, expected)
+        np.testing.assert_array_equal(is_intserseting, [False, False, True, True, False])
+
     def test_mesh_plane_intersection(self):
         # x-z plane
         normal = np.array([0., 1., 0.])