Add TetraMeshData class

raysect/primitive/mesh/__init__.pxd

@@ -30,3 +30,4 @@
 # POSSIBILITY OF SUCH DAMAGE.
 
 from raysect.primitive.mesh.mesh cimport Mesh, MeshIntersection
+from raysect.primitive.mesh.tetra_mesh cimport TetraMeshData

raysect/primitive/mesh/__init__.py

@@ -30,6 +30,7 @@
 # POSSIBILITY OF SUCH DAMAGE.
 
 from .mesh import Mesh, MeshIntersection
+from .tetra_mesh import TetraMeshData
 from .stl import import_stl, export_stl, STL_AUTOMATIC, STL_ASCII, STL_BINARY
 from .obj import import_obj, export_obj
 from .ply import import_ply, export_ply, PLY_AUTOMATIC, PLY_ASCII, PLY_BINARY

raysect/primitive/mesh/tests/test_tetra_mesh.py

@@ -0,0 +1,111 @@
+import pickle
+import unittest
+
+from raysect.core.math import AffineMatrix3D, Point3D
+from raysect.primitive.mesh.tetra_mesh import TetraMeshData
+
+# Configure Test Framework
+
+
+class TestTetraMeshData(unittest.TestCase):
+    def setUp(self):
+        # Define sample points and a single tetrahedron.
+        # Points of a unit tetrahedron: volume should be 1/6.
+        self.points = [
+            (0.0, 0.0, 0.0),
+            (1.0, 0.0, 0.0),
+            (0.0, 1.0, 0.0),
+            (0.0, 0.0, 1.0),
+        ]
+        self.tetrahedra = [
+            (0, 1, 2, 3),
+        ]
+        self.mesh = TetraMeshData(self.points, self.tetrahedra)
+
+    def test_initialization(self):
+        # Test that the mesh is created with the expected number of points and tetrahedra.
+        self.assertEqual(len(self.mesh.vertices), 4)
+        self.assertEqual(len(self.mesh.tetrahedra), 1)
+
+    def test_invalid_tetrahedron_indices(self):
+        # Test that constructing a mesh with tetrahedron indices out of bounds raises an error.
+        invalid_tetrahedra = [(0, 1, 2, 5)]  # Index 5 is out-of-range.
+        with self.assertRaises(IndexError):
+            TetraMeshData(self.points, invalid_tetrahedra)
+
+    def test_vertex_method(self):
+        # Test that the vertex method returns the correct point.
+        vertex0 = self.mesh.vertex(0)
+        self.assertAlmostEqual(vertex0.x, 0.0, places=5)
+        self.assertAlmostEqual(vertex0.y, 0.0, places=5)
+        self.assertAlmostEqual(vertex0.z, 0.0, places=5)
+
+    def test_invalid_vertex_index(self):
+        # Test that accessing a vertex with an out-of-range index raises an IndexError.
+        with self.assertRaises(IndexError):
+            self.mesh.vertex(10)
+
+    def test_barycenter(self):
+        # Test that barycenter of the tetrahedron is correctly computed.
+        # Barycenter is the average of the four vertices.
+        expected_barycenter = (
+            (0.0 + 1.0 + 0.0 + 0.0) / 4,
+            (0.0 + 0.0 + 1.0 + 0.0) / 4,
+            (0.0 + 0.0 + 0.0 + 1.0) / 4,
+        )
+        barycenter = self.mesh.barycenter(0)
+        self.assertAlmostEqual(barycenter.x, expected_barycenter[0], places=5)
+        self.assertAlmostEqual(barycenter.y, expected_barycenter[1], places=5)
+        self.assertAlmostEqual(barycenter.z, expected_barycenter[2], places=5)
+
+    def test_compute_volume(self):
+        # Assume TetraMeshData has a volume method that returns the volume of specified tetrahedron
+        # and a volume_total method that returns the total volume of all tetrahedra.
+        # The volume of the tetrahedron with the provided vertices is 1/6.
+        expected_volume = 1 / 6
+        volume = self.mesh.volume(0)
+        self.assertAlmostEqual(volume, expected_volume, places=5)
+
+        total_volume = self.mesh.volume_total()
+        self.assertAlmostEqual(total_volume, expected_volume, places=5)
+
+    def test_is_contained(self):
+        # Test that is_contained returns True for a point inside the tetrahedron
+        # and False for a point outside.
+        inside_point = Point3D(0.1, 0.1, 0.1)
+        outside_point = Point3D(1.0, 1.0, 1.0)
+        self.assertTrue(self.mesh.is_contained(inside_point))
+        self.assertFalse(self.mesh.is_contained(outside_point))
+
+    def test_bounding_box(self):
+        # Test the bounding box computation using an identity transformation.
+        # Construct an identity affine matrix.
+        identity = AffineMatrix3D([1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1])
+        bbox = self.mesh.bounding_box(identity)
+        # The tetrahedron vertices span from (0,0,0) to (1,1,1)
+        # and the bounding box is padded by BOX_PADDING (1e-6) on each side.
+        padding = 1e-6
+        self.assertAlmostEqual(bbox.lower.x, 0.0 - padding, places=5)
+        self.assertAlmostEqual(bbox.lower.y, 0.0 - padding, places=5)
+        self.assertAlmostEqual(bbox.lower.z, 0.0 - padding, places=5)
+        self.assertAlmostEqual(bbox.upper.x, 1.0 + padding, places=5)
+        self.assertAlmostEqual(bbox.upper.y, 1.0 + padding, places=5)
+        self.assertAlmostEqual(bbox.upper.z, 1.0 + padding, places=5)
+
+    def test_pickle_state(self):
+        # Test that the mesh can be pickled and unpickled without loss of state.
+        state = pickle.dumps(self.mesh)
+        new_mesh = pickle.loads(state)
+        self.assertEqual(len(new_mesh.vertices), len(self.mesh.vertices))
+        self.assertEqual(len(new_mesh.tetrahedra), len(self.mesh.tetrahedra))
+        self.assertAlmostEqual(new_mesh.volume(0), self.mesh.volume(0), places=5)
+        # Verify barycenter consistency.
+        bary_orig = self.mesh.barycenter(0)
+        bary_new = new_mesh.barycenter(0)
+        self.assertAlmostEqual(bary_new.x, bary_orig.x, places=5)
+        self.assertAlmostEqual(bary_new.y, bary_orig.y, places=5)
+        self.assertAlmostEqual(bary_new.z, bary_orig.z, places=5)
+
+
+if __name__ == "__main__":
+    unittest.main()

raysect/primitive/mesh/tetra_mesh.pxd

@@ -0,0 +1,73 @@
+# cython: language_level=3
+
+# Copyright (c) 2014-2023, Dr Alex Meakins, Raysect Project
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+#
+#     1. Redistributions of source code must retain the above copyright notice,
+#        this list of conditions and the following disclaimer.
+#
+#     2. Redistributions in binary form must reproduce the above copyright
+#        notice, this list of conditions and the following disclaimer in the
+#        documentation and/or other materials provided with the distribution.
+#
+#     3. Neither the name of the Raysect Project nor the names of its
+#        contributors may be used to endorse or promote products derived from
+#        this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+# POSSIBILITY OF SUCH DAMAGE.
+
+from numpy cimport ndarray, int32_t, uint8_t
+from raysect.core cimport BoundingBox3D, Point3D, AffineMatrix3D
+from raysect.core.math.spatial.kdtree3d cimport KDTree3DCore
+
+
+cdef class TetraMeshData(KDTree3DCore):
+
+    cdef:
+        ndarray _vertices
+        ndarray _tetrahedra
+        double[:, ::1] vertices_mv
+        int32_t[:, ::1] tetrahedra_mv
+        int32_t tetrahedra_id
+        int32_t i1, i2, i3, i4
+        double alpha, beta, gamma, delta
+        bint _cache_available
+        double _cached_x
+        double _cached_y
+        double _cached_z
+        bint _cached_result
+
+    cpdef Point3D vertex(self, int index)
+
+    cpdef ndarray tetrahedron(self, int index)
+
+    cpdef Point3D barycenter(self, int index)
+
+    cpdef double volume(self, int index)
+
+    cpdef double volume_total(self)
+
+    cdef double _volume(self, int index)
+
+    cdef object _filter_tetrahedra(self)
+
+    cdef BoundingBox3D _generate_bounding_box(self, int32_t tetrahedra)
+
+    cpdef BoundingBox3D bounding_box(self, AffineMatrix3D to_world)
+
+    cdef uint8_t _read_uint8(self, object file)
+
+    cdef bint _read_bool(self, object file)