From f24bda7a2ae1d43e903cde922a6f2dfb371987e8 Mon Sep 17 00:00:00 2001 From: munechika-koyo Date: Sat, 15 Feb 2025 11:15:02 +0100 Subject: [PATCH 1/5] Add `TetraMeshData` class --- raysect/primitive/mesh/__init__.pxd | 1 + raysect/primitive/mesh/__init__.py | 1 + raysect/primitive/mesh/tetra_mesh.pxd | 73 +++ raysect/primitive/mesh/tetra_mesh.pyx | 620 ++++++++++++++++++++++++++ 4 files changed, 695 insertions(+) create mode 100644 raysect/primitive/mesh/tetra_mesh.pxd create mode 100644 raysect/primitive/mesh/tetra_mesh.pyx diff --git a/raysect/primitive/mesh/__init__.pxd b/raysect/primitive/mesh/__init__.pxd index 7fa0f6bc..3cd7938f 100644 --- a/raysect/primitive/mesh/__init__.pxd +++ b/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 diff --git a/raysect/primitive/mesh/__init__.py b/raysect/primitive/mesh/__init__.py index a0d71a23..d4944233 100644 --- a/raysect/primitive/mesh/__init__.py +++ b/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 diff --git a/raysect/primitive/mesh/tetra_mesh.pxd b/raysect/primitive/mesh/tetra_mesh.pxd new file mode 100644 index 00000000..d3e56d0c --- /dev/null +++ b/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) \ No newline at end of file diff --git a/raysect/primitive/mesh/tetra_mesh.pyx b/raysect/primitive/mesh/tetra_mesh.pyx new file mode 100644 index 00000000..c920593d --- /dev/null +++ b/raysect/primitive/mesh/tetra_mesh.pyx @@ -0,0 +1,620 @@ +# cython: language_level=3 + +# Copyright (c) 2014-2021, 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. + +import io +import struct + +import numpy as np +cimport numpy as np +from libc.math cimport fabs +from raysect.core.math.spatial.kdtree3d cimport Item3D +from raysect.core.boundingbox cimport BoundingBox3D, new_boundingbox3d +from raysect.core.math cimport AffineMatrix3D, Point3D, new_point3d, Vector3D +from raysect.core.math.cython cimport barycentric_inside_tetrahedra, barycentric_coords_tetra +from cpython.bytes cimport PyBytes_AsString +cimport cython + +# bounding box is padded by a small amount to avoid numerical accuracy issues +DEF BOX_PADDING = 1e-6 + +# convenience defines +DEF V1 = 0 +DEF V2 = 1 +DEF V3 = 2 +DEF V4 = 3 + +DEF X = 0 +DEF Y = 1 +DEF Z = 2 + +# raysect mesh format constants +DEF RSM_VERSION_MAJOR = 1 +DEF RSM_VERSION_MINOR = 1 + + +cdef class TetraMeshData(KDTree3DCore): + """ + Holds the 3D tetrahedral mesh data and acceleration structures. + + This arrangement simplifies tetrahedral mesh instancing and the load/dump methods. + + The mesh vertices are supplied as an Nx3 list/array of floating point + values. For each Vertex, x, y and z coordinates must be supplied. e.g. + + vertices = [[0.0, 0.0, 1.0], [1.0, 0.0, 0.0], ...] + + The tetrahedral array is an Mx4. Tetrahedra are defined by indexing i.e: + + tetrahedra = [[v1, v2, v3, v4], ...] + + where v1, v2, v3, v4 are the vertex array indices specifying the tetrahedral vertices. + + :param object vertices: A list/array or tetrahedral vertices with shape Nx3, + where N is the number of vertices. + :param object tetrahedra: A list/array of tetrahedra with shape Mx4, + where M is the number of tetrahedra in the mesh. For each tetrahedra there + must be four integers identifying the tetrahedral vertices in the vertices array. + :param bool tolerant: Toggles filtering out of degenerate tetrahedra + (default=True). + + :ivar ndarray vertices: tetrahedral vertices with shape Nx3, where N is the number of vertices. + :ivar ndarray tetrahedra: tetrahedra with shape Mx4, where M is the number + of tetrahedra in the mesh. + """ + + def __init__(self, object vertices not None, object tetrahedra not None, bint tolerant=True): + + vertices = np.array(vertices, dtype=np.double) + tetrahedra = np.array(tetrahedra, dtype=np.int32) + + # check dimensions are correct + if vertices.ndim != 2 or vertices.shape[1] != 3: + raise ValueError("The vertex array must have dimensions Nx3.") + + if tetrahedra.ndim != 2 or tetrahedra.shape[1] != 4: + raise ValueError("The tetrahedra array must have dimensions Mx4.") + + # check tetrahedra contains only valid indices + invalid = (tetrahedra[:, 0:4] < 0) | (tetrahedra[:, 0:4] >= vertices.shape[0]) + if invalid.any(): + raise ValueError("The tetrahedra array references non-existent vertices.") + + # assign to internal attributes + self._vertices = vertices + self._tetrahedra = tetrahedra + + # assign to memory views + self.vertices_mv = vertices + self.tetrahedra_mv = tetrahedra + + # initialise hit state attributes + self.tetrahedra_id = -1 + self.i1 = -1 + self.i2 = -1 + self.i3 = -1 + self.i4 = -1 + self.alpha = 0.0 + self.beta = 0.0 + self.gamma = 0.0 + self.delta = 0.0 + + # filter out degenerate tetrahedra if we are being tolerant + if tolerant: + self._filter_tetrahedra() + + # kd-Tree init + items = [] + for tetrahedra in range(self._tetrahedra.shape[0]): + items.append(Item3D(tetrahedra, self._generate_bounding_box(tetrahedra))) + super().__init__(items, max_depth=0, min_items=1, hit_cost=50.0, empty_bonus=0.2) + + # TODO: (possible enhancement) check if tetrahedra are overlapping? + # (any non-owned vertex lying inside another tetrahedra) + + # init cache + self._cache_available = False + self._cached_x = 0.0 + self._cached_y = 0.0 + self._cached_z = 0.0 + self._cached_result = False + + def __getstate__(self): + state = io.BytesIO() + self.save(state) + return state.getvalue() + + def __setstate__(self, state): + self.load(io.BytesIO(state)) + + def __reduce__(self): + return self.__new__, (self.__class__, ), self.__getstate__() + + @property + def vertices(self): + return self._vertices.copy() + + @property + def tetrahedra(self): + return self._tetrahedra.copy() + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cpdef Point3D vertex(self, int index): + """ + Returns the specified vertex. + + :param int index: The vertex index. + :return: A Point3D object. + :rtype: Point3D + """ + + if index < 0 or index >= self.vertices_mv.shape[0]: + raise ValueError('Vertex index is out of range: [0, {}].'.format(self.vertices_mv.shape[0])) + + return new_point3d( + self.vertices_mv[index, X], + self.vertices_mv[index, Y], + self.vertices_mv[index, Z] + ) + + cpdef ndarray tetrahedron(self, int index): + """ + Returns the specified tetrahedron. + + The returned data will be a 4 element numpy array which are the tetrahedral vertex indices. + + :param int index: The tetrahedral index. + :return: A numpy array. + :rtype: ndarray + """ + + if index < 0 or index >= self.vertices_mv.shape[0]: + raise ValueError('Tetrahedral index is out of range: [0, {}].'.format(self.tetrahedra_mv.shape[0])) + + return self._tetrahedra[index, :].copy() + + cpdef Point3D barycenter(self, int index): + """ + Returns the barycenter point of the specified tetrahedron. + + :param int index: The tetrahedral index. + :return: A barycenter point + :rtype: Point3D + """ + cdef: + np.int32_t i1, i2, i3, i4 + + if index < 0 or index >= self.tetrahedra_mv.shape[0]: + raise ValueError('Tetrahedral index is out of range: [0, {}].'.format(self.tetrahedra_mv.shape[0])) + + i1 = self.tetrahedra_mv[index, V1] + i2 = self.tetrahedra_mv[index, V2] + i3 = self.tetrahedra_mv[index, V3] + i4 = self.tetrahedra_mv[index, V4] + + return new_point3d( + 0.25 * (self.vertices_mv[i1, X] + self.vertices_mv[i2, X] + self.vertices_mv[i3, X] + self.vertices_mv[i4, X]), + 0.25 * (self.vertices_mv[i1, Y] + self.vertices_mv[i2, Y] + self.vertices_mv[i3, Y] + self.vertices_mv[i4, Y]), + 0.25 * (self.vertices_mv[i1, Z] + self.vertices_mv[i2, Z] + self.vertices_mv[i3, Z] + self.vertices_mv[i4, Z]), + ) + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cpdef double volume(self, int index): + """ + Calculate a volume of the specified tetrahedron + + :param int index: The tetrahedral index. + :return: A volume of specified tetrahedron + :rtype: double + """ + if index < 0 or index >= self.tetrahedra_mv.shape[0]: + raise ValueError('Tetrahedral index is out of range: [0, {}].'.format(self.tetrahedra_mv.shape[0])) + + return self._volume(index) + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cpdef double volume_total(self): + """ + Calculate a total volume of all tetrahedra + + :return: total volume of all tetrahedra + :rtype: double + """ + cdef: + np.int32_t i + double volume = 0.0 + + for i in range(self.tetrahedra_mv.shape[0]): + volume += self._volume(i) + + return volume + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + @cython.cdivision(True) + cdef double _volume(self, int index): + """ + Fast calculation of volume of a tetrahedron + """ + cdef: + np.int32_t i1, i2, i3, i4 + Point3D p1, p2, p3, p4 + Vector3D v1, v2, v3, V4 + double area, height + + i1 = self.tetrahedra_mv[index, V1] + i2 = self.tetrahedra_mv[index, V2] + i3 = self.tetrahedra_mv[index, V3] + i4 = self.tetrahedra_mv[index, V4] + + p1 = new_point3d(self.vertices_mv[i1, X], self.vertices_mv[i1, Y], self.vertices_mv[i1, Z]) + p2 = new_point3d(self.vertices_mv[i2, X], self.vertices_mv[i2, Y], self.vertices_mv[i2, Z]) + p3 = new_point3d(self.vertices_mv[i3, X], self.vertices_mv[i3, Y], self.vertices_mv[i3, Z]) + p4 = new_point3d(self.vertices_mv[i4, X], self.vertices_mv[i4, Y], self.vertices_mv[i4, Z]) + + v1 = p1.vector_to(p2) + v2 = p1.vector_to(p3) + v3 = p1.vector_to(p4) + + # area of the base + v4 = v1.cross(v2) + area = v4.get_length() * 0.5 + + # height from the base to the apex + height = fabs(v4.normalise().dot(v3)) + + return (area * height) / 3.0 + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cdef object _filter_tetrahedra(self): + + cdef: + np.int32_t i, valid + + # scan tetrahedra and make valid tetrahedra contiguous + valid = 0 + for i in range(self.tetrahedra_mv.shape[0]): + + if self._volume(i) == 0.0: + + # tetrahedron is degenerate, skip + continue + + # shift tetrahedra + self.tetrahedra_mv[valid, :] = self.tetrahedra_mv[i, :] + valid += 1 + + # reslice array to contain only valid tetrahedra + self.tetrahedra_mv = self.tetrahedra_mv[:valid, :] + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cdef BoundingBox3D _generate_bounding_box(self, np.int32_t tetrahedra): + """ + Generates a bounding box for the specified tetrahedron. + + A small degree of padding is added to the bounding box to provide the + conservative bounds required by the watertight mesh algorithm. + + :param tetrahedra: tetrahedral array index. + :return: A BoundingBox3D object. + :rtype: BoundingBox3D + """ + + cdef: + np.int32_t i1, i2, i3, i4 + BoundingBox3D bbox + + i1 = self.tetrahedra_mv[tetrahedra, V1] + i2 = self.tetrahedra_mv[tetrahedra, V2] + i3 = self.tetrahedra_mv[tetrahedra, V3] + i4 = self.tetrahedra_mv[tetrahedra, V4] + + bbox = new_boundingbox3d( + new_point3d( + min(self.vertices_mv[i1, X], self.vertices_mv[i2, X], self.vertices_mv[i3, X], self.vertices_mv[i4, X]), + min(self.vertices_mv[i1, Y], self.vertices_mv[i2, Y], self.vertices_mv[i3, Y], self.vertices_mv[i4, Y]), + min(self.vertices_mv[i1, Z], self.vertices_mv[i2, Z], self.vertices_mv[i3, Z], self.vertices_mv[i4, Z]), + ), + new_point3d( + max(self.vertices_mv[i1, X], self.vertices_mv[i2, X], self.vertices_mv[i3, X], self.vertices_mv[i4, X]), + max(self.vertices_mv[i1, Y], self.vertices_mv[i2, Y], self.vertices_mv[i3, Y], self.vertices_mv[i4, Y]), + max(self.vertices_mv[i1, Z], self.vertices_mv[i2, Z], self.vertices_mv[i3, Z], self.vertices_mv[i4, Z]), + ), + ) + + # The bounding box and tetrahedral vertices may not align following coordinate + # transforms in the water tight mesh algorithm, therefore a small bit of padding + # is added to avoid numerical representation issues. + bbox.pad(max(BOX_PADDING, bbox.largest_extent() * BOX_PADDING)) + + return bbox + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cdef bint _is_contained_leaf(self, np.int32_t id, Point3D point): + + cdef: + np.int32_t index, tetrahedra_id, i1, i2, i3, i4 + double alpha, beta, gamma, delta + + # identify the first tetrahedra that contains the point, if any + for index in range(self._nodes[id].count): + + # obtain vertex indices + tetrahedra_id = self._nodes[id].items[index] + i1 = self.tetrahedra_mv[tetrahedra_id, V1] + i2 = self.tetrahedra_mv[tetrahedra_id, V2] + i3 = self.tetrahedra_mv[tetrahedra_id, V3] + i4 = self.tetrahedra_mv[tetrahedra_id, V4] + + barycentric_coords_tetra(self.vertices_mv[i1, X], self.vertices_mv[i1, Y], self.vertices_mv[i1, Z], + self.vertices_mv[i2, X], self.vertices_mv[i2, Y], self.vertices_mv[i2, Z], + self.vertices_mv[i3, X], self.vertices_mv[i3, Y], self.vertices_mv[i3, Z], + self.vertices_mv[i4, X], self.vertices_mv[i4, Y], self.vertices_mv[i4, Z], + point.x, point.y, point.z, &alpha, &beta, &gamma, &delta) + + if barycentric_inside_tetrahedra(alpha, beta, gamma, delta): + + # store id of tetrahedral hit + self.tetrahedra_id = tetrahedra_id + + # store vertex indices and barycentric coords + self.i1 = i1 + self.i2 = i2 + self.i3 = i3 + self.i4 = i4 + self.alpha = alpha + self.beta = beta + self.gamma = gamma + self.delta = delta + + return True + + return False + + cpdef bint is_contained(self, Point3D point): + """ + Traverses the kd-Tree to identify if the point is contained by an item. + :param Point3D point: A Point3D object. + :return: True if the point lies inside an item, false otherwise. + :rtype: bool + """ + + cdef bint result + + if self._cache_available and point.x == self._cached_x and point.y == self._cached_y and point.z == self._cached_z: + return self._cached_result + + result = self._is_contained(point) + + # add cache + self._cache_available = True + self._cached_x = point.x + self._cached_y = point.y + self._cached_z = point.z + self._cached_result = result + + return result + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.initializedcheck(False) + cpdef BoundingBox3D bounding_box(self, AffineMatrix3D to_world): + """ + Returns a bounding box that encloses the mesh. + + The box is padded by a small margin to reduce the risk of numerical + accuracy problems between the mesh and box representations following + coordinate transforms. + + :param to_world: Local to world space transform matrix. + :return: A BoundingBox3D object. + """ + + cdef: + np.int32_t i + BoundingBox3D bbox + Point3D vertex + + # TODO: padding should really be a function of mesh extent + # convert vertices to world space and grow a bounding box around them + bbox = BoundingBox3D() + for i in range(self.vertices_mv.shape[0]): + vertex = new_point3d(self.vertices_mv[i, X], self.vertices_mv[i, Y], self.vertices_mv[i, Z]) + bbox.extend(vertex.transform(to_world), BOX_PADDING) + + return bbox + + @cython.boundscheck(False) + @cython.wraparound(False) + def save(self, object file): + """ + Save the mesh's kd-Tree representation to a binary Raysect mesh file (.rsm). + + :param object file: File stream or string file name to save state. + """ + + cdef: + np.int32_t i, j + double[:, ::1] vertices + np.int32_t[:, ::1] tetrahedra + + close = False + + # treat as a filename if a stream is not supplied + if not isinstance(file, io.IOBase): + file = open(file, mode="wb") + close = True + + # hold local references to avoid repeated memory view object checks + vertices = self.vertices_mv + tetrahedra = self.tetrahedra_mv + + # write header + file.write(b"RSM") + file.write(struct.pack(">> from raysect.primitive.mesh import TetraMeshData + >>> + >>> mesh = TetraMeshData.from_file("test.rsm") + """ + + m = cls.__new__(cls) + m.load(file) + return m + + cdef uint8_t _read_uint8(self, object file): + return ( PyBytes_AsString(file.read(sizeof(uint8_t))))[0] + + cdef bint _read_bool(self, object file): + return self._read_uint8(file) != 0 From 1234389053e259ce01b5c5c6586dd8aa1a5e1a10 Mon Sep 17 00:00:00 2001 From: munechika-koyo Date: Sat, 15 Feb 2025 12:30:51 +0100 Subject: [PATCH 2/5] Change ValueError to IndexError for invalid tetrahedra references --- raysect/primitive/mesh/tetra_mesh.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/raysect/primitive/mesh/tetra_mesh.pyx b/raysect/primitive/mesh/tetra_mesh.pyx index c920593d..82cee09d 100644 --- a/raysect/primitive/mesh/tetra_mesh.pyx +++ b/raysect/primitive/mesh/tetra_mesh.pyx @@ -105,7 +105,7 @@ cdef class TetraMeshData(KDTree3DCore): # check tetrahedra contains only valid indices invalid = (tetrahedra[:, 0:4] < 0) | (tetrahedra[:, 0:4] >= vertices.shape[0]) if invalid.any(): - raise ValueError("The tetrahedra array references non-existent vertices.") + raise IndexError("The tetrahedra array references non-existent vertices.") # assign to internal attributes self._vertices = vertices From 857fc33461ffbd0ada1103d89968949f285318ad Mon Sep 17 00:00:00 2001 From: munechika-koyo Date: Sat, 15 Feb 2025 12:40:56 +0100 Subject: [PATCH 3/5] Change ValueError to IndexError for out-of-range vertex and tetrahedral index checks --- raysect/primitive/mesh/tetra_mesh.pyx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/raysect/primitive/mesh/tetra_mesh.pyx b/raysect/primitive/mesh/tetra_mesh.pyx index 82cee09d..103d0223 100644 --- a/raysect/primitive/mesh/tetra_mesh.pyx +++ b/raysect/primitive/mesh/tetra_mesh.pyx @@ -178,7 +178,7 @@ cdef class TetraMeshData(KDTree3DCore): """ if index < 0 or index >= self.vertices_mv.shape[0]: - raise ValueError('Vertex index is out of range: [0, {}].'.format(self.vertices_mv.shape[0])) + raise IndexError('Vertex index is out of range: [0, {}].'.format(self.vertices_mv.shape[0])) return new_point3d( self.vertices_mv[index, X], @@ -198,7 +198,7 @@ cdef class TetraMeshData(KDTree3DCore): """ if index < 0 or index >= self.vertices_mv.shape[0]: - raise ValueError('Tetrahedral index is out of range: [0, {}].'.format(self.tetrahedra_mv.shape[0])) + raise IndexError('Tetrahedral index is out of range: [0, {}].'.format(self.tetrahedra_mv.shape[0])) return self._tetrahedra[index, :].copy() @@ -214,7 +214,7 @@ cdef class TetraMeshData(KDTree3DCore): np.int32_t i1, i2, i3, i4 if index < 0 or index >= self.tetrahedra_mv.shape[0]: - raise ValueError('Tetrahedral index is out of range: [0, {}].'.format(self.tetrahedra_mv.shape[0])) + raise IndexError('Tetrahedral index is out of range: [0, {}].'.format(self.tetrahedra_mv.shape[0])) i1 = self.tetrahedra_mv[index, V1] i2 = self.tetrahedra_mv[index, V2] @@ -239,7 +239,7 @@ cdef class TetraMeshData(KDTree3DCore): :rtype: double """ if index < 0 or index >= self.tetrahedra_mv.shape[0]: - raise ValueError('Tetrahedral index is out of range: [0, {}].'.format(self.tetrahedra_mv.shape[0])) + raise IndexError('Tetrahedral index is out of range: [0, {}].'.format(self.tetrahedra_mv.shape[0])) return self._volume(index) From 5ada94122d18dd487afb9d79866d57b72e105e0f Mon Sep 17 00:00:00 2001 From: munechika-koyo Date: Sat, 15 Feb 2025 13:01:40 +0100 Subject: [PATCH 4/5] Add unit tests for `TetraMeshData` class --- .../primitive/mesh/tests/test_tetra_mesh.py | 111 ++++++++++++++++++ 1 file changed, 111 insertions(+) create mode 100644 raysect/primitive/mesh/tests/test_tetra_mesh.py diff --git a/raysect/primitive/mesh/tests/test_tetra_mesh.py b/raysect/primitive/mesh/tests/test_tetra_mesh.py new file mode 100644 index 00000000..ca5931d4 --- /dev/null +++ b/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() From 78ae8a4bd741f2524bd5292e7815ec7a2d850c8d Mon Sep 17 00:00:00 2001 From: munechika-koyo Date: Sat, 15 Feb 2025 20:07:41 +0100 Subject: [PATCH 5/5] Refactor comments for consistency in `TetraMeshData` class --- raysect/primitive/mesh/tetra_mesh.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/raysect/primitive/mesh/tetra_mesh.pyx b/raysect/primitive/mesh/tetra_mesh.pyx index 103d0223..bf5852c1 100644 --- a/raysect/primitive/mesh/tetra_mesh.pyx +++ b/raysect/primitive/mesh/tetra_mesh.pyx @@ -495,7 +495,7 @@ cdef class TetraMeshData(KDTree3DCore): file.write(struct.pack("