Add `TetraMeshData` class

[munechika-koyo]

Proposal of new feature: TetraMeshData

TL;DR: Reading files in the created tree structure will greatly enhance the speed of the process with this PR.

Description

I would propose a new feature, the TetraMeshData class, modeled after the existing MeshData definition. (I have already suggested #407 but would like to split the function.)

It was implemented with the following objectives:

• Saving/Loading from the file formatted .rsm, which has constructed a K-D tree structure.
• Offering useful methods like calculating a tetrahedral volume and barycenter.

Additionally, this class will be used when loading mesh files formatted as .vtk or .obj, which contain tetrahedral (3-D unstructured) mesh data in the future.

Examples of the use:

1. Creating the instance of TetraMeshData

>>> from raysect.primitive.mesh.tetra_mesh import TetraMeshData
>>> verts = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]]
>>> tets = [[0, 1, 2, 3]]
>>> tetra = TetraMeshData(verts, tets)

2. Calculating the volume of a tetrahedron

>>> tetra.volume(0)
0.166666  # == 1 x 1 x 0.5 / 3.0
>>> tetra.volume_total()
0.166666

3. Calculating the barycenter of a tetrahedron

>>> tetra.barycenter(0)
Point3D(0.25, 0.25, 0.25)

4. Saving as .rsm

>>> tetra.save("tetra.rsm")

5. Loading from a .rsm file

>>> TetraMeshData.from_file("tetra.rsm")
<raysect.primitive.mesh.tetra_mesh.TetraMeshData at 0x1336e7440>

Speed test

I attempted to compare the instancing of Discrete3DMesh with TetraMeshData loaded from a file. As a mesh file, I used the stanford_bunny.mesh file located in demos/resources/.

	elapse time (sec)
Create Discrete3DMesh	15.441215
Create TetraMeshData	15.894989
Load TetraMeshData from file	01.766619

I attained a loading speed for tetra mesh that is approximately 9 times faster than generating TetraMeshData using vertices and tetrahedral index arrays. This difference is expected to widen as the vertex data size increases. The test script that I used is in the collapsed section below.

Test script

from datetime import timedelta
from pathlib import Path
from timeit import timeit

import numpy as np

from raysect.core.math.function.float.function3d.interpolate import Discrete3DMesh
from raysect.primitive.mesh.tetra_mesh import TetraMeshData


def load_mesh(mesh_filepath):
    vertices = []
    tetrahedra = []

    with open(mesh_filepath, "r") as f:
        lines = f.readlines()

    i = 0
    while i < len(lines):
        line = lines[i].strip()
        if not line or line.startswith("#"):
            i += 1
            continue

        # Check section headers (case-sensitive)
        if line == "Vertices":
            i += 1
            vertex_count = int(lines[i].strip())
            i += 1
            for _ in range(vertex_count):
                parts = lines[i].strip().split()
                # take only the first 3 coordinates
                vertices.append([float(v) for v in parts[:3]])
                i += 1
            continue
        elif line == "Tetrahedra":
            i += 1
            tet_count = int(lines[i].strip())
            i += 1
            for _ in range(tet_count):
                parts = lines[i].strip().split()
                # take only the first 4 indices, converting them to int.
                tetrahedra.append([int(idx) for idx in parts[:4]])
                i += 1
            continue

        i += 1

    vertices = np.array(vertices)
    tetrahedra = np.array(tetrahedra, dtype=np.int32) - 1  # 0-based indexing
    return vertices, tetrahedra


def create_discrete3dmesh(vertices, tetrahedra):
    return Discrete3DMesh(
        vertices, tetrahedra, np.ones((tetrahedra.shape[0])), False, 0
    )


def create_tetrameshdata(vertices, tetrahedra):
    return TetraMeshData(vertices, tetrahedra)


def load_tetrameshdata(mesh_filepath):
    return TetraMeshData.from_file(mesh_filepath)


if __name__ == "__main__":
    ROOT = Path(__file__).parent
    mesh_file = ROOT / "demos" / "resources" / "stanford_bunny.mesh"
    vertices, tetrahedra = load_mesh(mesh_file)

    # Save tetramesh data in advance
    tetra = TetraMeshData(vertices, tetrahedra)
    tetra_file = ROOT / "temp_tetra.rsm"
    tetra.save(tetra_file)

    # === Measure time ===
    loop = 5

    # create Discrete3DMesh
    result = timeit(
        "create_discrete3dmesh(vertices, tetrahedra)", globals=globals(), number=loop
    )
    elapsed_time = timedelta(seconds=result / loop)
    print(f"Create Discrete3DMesh: {elapsed_time}")

    # create TetraMeshData
    result = timeit(
        "create_tetrameshdata(vertices, tetrahedra)", globals=globals(), number=loop
    )
    elapsed_time = timedelta(seconds=result / loop)
    print(f"Create TetraMeshData: {elapsed_time}")

    # load TetraMeshData from file
    result = timeit("load_tetrameshdata(tetra_file)", globals=globals(), number=loop)
    elapsed_time = timedelta(seconds=result / loop)
    print(f"Load TetraMeshData from file: {elapsed_time}")

Class structure

The structure of the TetraMeshData I implemented is as follows.

classDiagram
    class KDTree3DCore {
        +bint is_contained(Point3D point)
        +void save(object file)
        +void load(object file)
    }

    class TetraMeshData {
        +__init__(object vertices, object tetrahedra, bint tolerant=True)
        +__getstate__()
        +__setstate__(state)
        +__reduce__()
        +vertices
        +tetrahedra
        +Point3D vertex(int index)
        +ndarray tetrahedron(int index)
        +Point3D barycenter(int index)
        +double volume(int index)
        +double volume_total()
        +BoundingBox3D bounding_box(AffineMatrix3D to_world)
        +bint is_contained(Point3D point)
        +void save(object file)
        +void load(object file)
        +classmethod from_file(cls, file)
    }

    KDTree3DCore <|-- TetraMeshData

Details of unit test

I also implemented a unit test for TestMeshData in raysect/primitive/mesh/tests/test_tetra_mesh.py. Below is a table showing the correspondence between the methods of the TestTetraMeshData class in unit tests and the methods tested in the TestMeshData class.

Test Method Name in TestTetraMeshData	Tested Method/Function in TetraMeshData
test_initialization()	__init__()
test_invalid_tetrahedron_indices()	__init__()
test_vertex_method()	vertex()
test_invalid_vertex_index()	vertex()
test_barycenter()	barycenter()
test_compute_volume()	volume(), volume_total()
test_is_contained()	is_contained()
test_bounding_box()	bounding_box()
test_pickle_state ()	__getstate__(), __setstate__(), save(), load()

---

I would appreciate it if you would review it and make any suggestions and comments.