Surface normals, and pull in shapes from lace (#7)

Author: paulmelnikow

blmath/geometry/shapes.py

@@ -0,0 +1,96 @@
+import numpy as np
+
+
+def _maybe_flatten(vertices, faces, ret_unique_vertices_and_faces):
+    if ret_unique_vertices_and_faces:
+        return vertices, faces
+    else:
+        return vertices[faces]
+
+
+def create_rectangular_prism(origin, size, ret_unique_vertices_and_faces=False):
+    '''
+    Return vertices (or unique verties and faces) of an axis-aligned
+    rectangular prism. One vertex is `origin`; the diametrically opposite
+    vertex is `origin + size`.
+
+    size: 3x1 array.
+
+    '''
+    from lace.topology import quads_to_tris
+
+    lower_base_plane = np.array([
+        # Lower base plane
+        origin,
+        origin + np.array([size[0], 0, 0]),
+        origin + np.array([size[0], 0, size[2]]),
+        origin + np.array([0, 0, size[2]]),
+    ])
+    upper_base_plane = lower_base_plane + np.array([0, size[1], 0])
+
+    vertices = np.vstack([lower_base_plane, upper_base_plane])
+
+    faces = quads_to_tris(np.array([
+        [0, 1, 2, 3],  # lower base (-y)
+        [7, 6, 5, 4],  # upper base (+y)
+        [4, 5, 1, 0],  # +z face
+        [5, 6, 2, 1],  # +x face
+        [6, 7, 3, 2],  # -z face
+        [3, 7, 4, 0],  # -x face
+    ]))
+
+    return _maybe_flatten(vertices, faces, ret_unique_vertices_and_faces)
+
+
+def create_cube(origin, size, ret_unique_vertices_and_faces=False):
+    '''
+    Return vertices (or unique verties and faces) with an axis-aligned cube.
+    One vertex is `origin`; the diametrically opposite vertex is `size` units
+    along +x, +y, and +z.
+
+    size: int or float.
+
+    '''
+    return create_rectangular_prism(
+        origin,
+        np.repeat(size, 3),
+        ret_unique_vertices_and_faces=ret_unique_vertices_and_faces)
+
+
+def create_triangular_prism(p1, p2, p3, height, ret_unique_vertices_and_faces=False):
+    """
+    Return vertices (or unique verties and faces) of a triangular prism whose
+    base is the triangle p1, p2, p3. If the vertices are oriented in a
+    counterclockwise direction, the prism extends from behind them.
+
+    Imported from lace.
+    """
+    from . import Plane
+
+    base_plane = Plane.from_points(p1, p2, p3)
+    lower_base_to_upper_base = height * -base_plane.normal # pylint: disable=invalid-unary-operand-type
+    vertices = np.vstack(([p1, p2, p3], [p1, p2, p3] + lower_base_to_upper_base))
+
+    faces = np.array([
+        [0, 1, 2],  # base
+        [0, 3, 4], [0, 4, 1],  # side 0, 3, 4, 1
+        [1, 4, 5], [1, 5, 2],  # side 1, 4, 5, 2
+        [2, 5, 3], [2, 3, 0],  # side 2, 5, 3, 0
+        [5, 4, 3],  # base
+    ])
+
+    return _maybe_flatten(vertices, faces, ret_unique_vertices_and_faces)
+
+
+def create_horizontal_plane(ret_unique_vertices_and_faces=False):
+    '''
+    Creates a horizontal plane.
+    '''
+    vertices = np.array([
+        [1., 0., 0.],
+        [-1., 0., 0.],
+        [0., 0., 1.],
+        [0., 0., -1.]
+    ])
+    faces = [[0, 1, 2], [3, 1, 0]]
+    return _maybe_flatten(vertices, faces, ret_unique_vertices_and_faces)

blmath/geometry/surface_normals.py

@@ -0,0 +1,13 @@
+def surface_normal(points, normalize=True):
+    """
+    Compute the surface normal of a triangle.
+
+    Can provide three points or an array of points.
+    """
+    import numpy as np
+    from blmath.numerics import vx
+    p1 = points[..., 0, :]
+    p2 = points[..., 1, :]
+    p3 = points[..., 2, :]
+    normal = np.cross(p2 - p1, p3 - p1)
+    return vx.normalize(normal) if normalize else normal

blmath/geometry/test_surface_normals.py

@@ -0,0 +1,44 @@
+import unittest
+import numpy as np
+from ..numerics import vx
+from .surface_normals import surface_normal
+
+
+class TestSurfaceNormals(unittest.TestCase):
+    def test_surface_normal_single(self):
+        points = np.array([
+            [3.0, 0.0, 0.0],
+            [0.0, 3.0, 0.0],
+            [0.0, 0.0, 3.0],
+        ])
+
+        np.testing.assert_allclose(
+            surface_normal(points),
+            np.array([3 ** -0.5, 3 ** -0.5, 3 ** -0.5]))
+
+        np.testing.assert_allclose(
+            surface_normal(points, normalize=False),
+            np.array([9.0, 9.0, 9.0]))
+
+    def test_surface_normal_vectorized(self):
+        from .shapes import create_triangular_prism
+
+        p1 = np.array([3.0, 0.0, 0.0])
+        p2 = np.array([0.0, 3.0, 0.0])
+        p3 = np.array([0.0, 0.0, 3.0])
+        vertices = create_triangular_prism(p1, p2, p3, 1.0)
+
+        expected_normals = np.array([
+            vx.normalize(np.array([1.0, 1.0, 1.0])),
+            vx.normalize(np.array([1.0, 1.0, -2.0])),
+            vx.normalize(np.array([1.0, 1.0, -2.0])),
+            vx.normalize(np.array([-2.0, 1.0, 1.0])),
+            vx.normalize(np.array([-2.0, 1.0, 1.0])),
+            vx.normalize(np.array([1.0, -2.0, 1.0])),
+            vx.normalize(np.array([1.0, -2.0, 1.0])),
+            vx.normalize(np.array([-1.0, -1.0, -1.0])),
+        ])
+
+        np.testing.assert_allclose(
+            surface_normal(vertices),
+            expected_normals)