added geometry._Line._intersect_line()

freesurfer_surface/geometry.py

@@ -5,8 +5,7 @@ import numpy
 
 def _collinear(vector_a: numpy.ndarray, vector_b: numpy.ndarray) -> bool:
     # null vector: https://math.stackexchange.com/a/1772580
-    return numpy.allclose(numpy.cross(vector_a, vector_b),
-                          numpy.zeros(len(vector_a)))
+    return numpy.allclose(numpy.cross(vector_a, vector_b), 0)
 
 
 class _Line:
@@ -25,6 +24,17 @@ class _Line:
     def __repr__(self) -> str:
         return 'line(t) = {} + {} t'.format(self.point, self.vector)
 
+    def _intersect_line(self, other: '_Line') \
+            -> typing.Union[numpy.ndarray, bool]:
+        # https://en.wikipedia.org/wiki/Skew_lines#Distance
+        lines_normal_vector = numpy.cross(self.vector, other.vector)
+        if numpy.allclose(lines_normal_vector, 0):
+            return _collinear(self.vector, self.point - other.point)
+        plane_normal_vector = numpy.cross(other.vector, lines_normal_vector)
+        return self.point + self.vector \
+            * (numpy.inner(other.point - self.point, plane_normal_vector)
+               / numpy.inner(self.vector, plane_normal_vector))
+
 
 def _intersect_planes(normal_vector_a: numpy.ndarray,
                       constant_a: float,

tests/geometry/test_line.py

@@ -62,3 +62,54 @@ def test__equal(line_a, line_b, equal):
 def test_repr():
     line = _Line(point=[1, 2, 3], vector=[4, 5, 6])
     assert repr(line) == 'line(t) = [1. 2. 3.] + [4. 5. 6.] t'
+
+
+@pytest.mark.parametrize(('line_a', 'line_b', 'expected_point'), [
+    (_Line(point=(1, 2, 3), vector=(0, 0, 4)),
+     _Line(point=(1, 2, 3), vector=(0, 5, 0)),
+     [1, 2, 3]),
+    (_Line(point=(1, 2, 7), vector=(0, 0, 4)),
+     _Line(point=(1, -8, 3), vector=(0, 5, 0)),
+     [1, 2, 3]),
+    (_Line(point=(1, 2, 3), vector=(3, 2, 4)),
+     _Line(point=(1, 2, 3), vector=(4, -5, 9)),
+     [1, 2, 3]),
+    (_Line(point=(-2, 0, -1), vector=(3, 2, 4)),
+     _Line(point=(1, 2, 3), vector=(4, -5, 9)),
+     [1, 2, 3]),
+    (_Line(point=(-2, 0, -1), vector=(3, 2, 4)),
+     _Line(point=(9, -8, 21), vector=(4, -5, 9)),
+     [1, 2, 3]),
+    (_Line(point=(-7, 4, -2), vector=(2, 6, 3)),
+     _Line(point=(-7, 4, -2), vector=(-4, 8, -3)),
+     [-7, 4, -2]),
+    (_Line(point=(-5, 10, 1), vector=(2, 6, 3)),
+     _Line(point=(-15, 20, -8), vector=(-4, 8, -3)),
+     [-7, 4, -2]),
+    (_Line(point=(1, 2, 3), vector=(4, 8, 7)),
+     _Line(point=(1, 2, 3), vector=(4, 8, 7)),
+     True),
+    (_Line(point=(1, 2, 3), vector=(4, 8, 7)),
+     _Line(point=(1, 2, 3), vector=(8, 16, 14)),
+     True),
+    (_Line(point=(1, 2, 3), vector=(-4, -8, -7)),
+     _Line(point=(1, 2, 3), vector=(8, 16, 14)),
+     True),
+    (_Line(point=(-3, -6, -4), vector=(-4, -8, -7)),
+     _Line(point=(1, 2, 3), vector=(8, 16, 14)),
+     True),
+    (_Line(point=(-3, -6, -4), vector=(-4, -8, -7)),
+     _Line(point=(5, 10, 10), vector=(8, 16, 14)),
+     True),
+    (_Line(point=(-3, -6, -3), vector=(-4, -8, -7)),
+     _Line(point=(5, 10, 10), vector=(8, 16, 14)),
+     False),
+])
+def test__intersect_line(line_a, line_b, expected_point):
+    # pylint: disable=protected-access
+    point = line_a._intersect_line(line_b)
+    if isinstance(expected_point, bool):
+        assert isinstance(point, bool)
+        assert point == expected_point
+    else:
+        assert numpy.allclose(point, expected_point)