added geometry._intersect_planes

freesurfer_surface/geometry.py

@@ -1,3 +1,5 @@
+import typing
+
 import numpy
 
 
@@ -19,3 +21,20 @@ class _Line:
         if not _collinear(self.vector, other.vector):
             return False
         return _collinear(self.vector, self.point - other.point)
+
+    def __repr__(self) -> str:
+        return 'line(t) = {} + {} t'.format(self.point, self.vector)
+
+
+def _intersect_planes(normal_vector_a: numpy.ndarray,
+                      constant_a: float,
+                      normal_vector_b: numpy.ndarray,
+                      constant_b: float) -> typing.Union[_Line, bool]:
+    line_vector = numpy.cross(normal_vector_a, normal_vector_b)
+    if numpy.allclose(line_vector, 0):
+        return constant_a == constant_b
+    point = numpy.linalg.solve(
+        numpy.vstack((normal_vector_a, normal_vector_b, line_vector)),
+        numpy.vstack((constant_a, constant_b, 0)),
+    )
+    return _Line(point=point.reshape(3), vector=line_vector)

tests/geometry/test_intersect.py

@@ -0,0 +1,39 @@
+import numpy
+import pytest
+
+from freesurfer_surface.geometry import _intersect_planes, _Line
+
+
+@pytest.mark.parametrize(
+    ('normal_vector_a', 'constant_a',
+     'normal_vector_b', 'constant_b',
+     'expected_line'),
+    [([0, 0, 1], 0, [0, 1, 0], 0, _Line([0, 0, 0], [1, 0, 0])),
+     ([0, 0, 1], 0, [0, 3, 0], 0, _Line([0, 0, 0], [1, 0, 0])),
+     ([0, 0, 1], 0, [4, 0, 0], 0, _Line([0, 0, 0], [0, 1, 0])),
+     ([0, 2, 2], 0, [3, 0, 3], 0, _Line([0, 0, 0], [1, 1, -1])),
+     ([1, 2, 4], 0, [2, 3, 5], 0, _Line([0, 0, 0], [-2, 3, -1])),
+     ([1, 2, 4], 0, [2, 3, 5], 2, _Line([2, 1, -1], [-2, 3, -1])),
+     ([2, 3, 5], 2, [1, 2, 4], 0, _Line([2, 1, -1], [-2, 3, -1])),
+     ([2, 3, 5], 2, [1, 2, 4], 7, _Line([-7, -3, 5], [-2, 3, -1])),
+     ([1, 2, 4], 0, [2, 4, 8], 1, False),
+     ([1, 2, 4], 0, [2, 4, 8], 0, True)],
+)
+def test__intersect_planes(normal_vector_a, constant_a,
+                           normal_vector_b, constant_b,
+                           expected_line):
+    line = _intersect_planes(normal_vector_a, constant_a,
+                             normal_vector_b, constant_b)
+    assert line == expected_line
+    if not isinstance(expected_line, bool):
+        assert numpy.isclose(numpy.inner(normal_vector_a, line.vector), 0)
+        assert numpy.isclose(numpy.inner(normal_vector_b, line.vector), 0)
+        assert numpy.isclose(numpy.inner(normal_vector_a, line.point),
+                             constant_a)
+        assert numpy.isclose(numpy.inner(normal_vector_b, line.point),
+                             constant_b)
+        other_point = line.point + line.vector
+        assert numpy.isclose(numpy.inner(normal_vector_a, other_point),
+                             constant_a)
+        assert numpy.isclose(numpy.inner(normal_vector_b, other_point),
+                             constant_b)

tests/geometry/test_line.py

@@ -57,3 +57,8 @@ def test_init_numpy_array():
 ])
 def test__equal(line_a, line_b, equal):
     assert (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'