Skip to content

Add TetraMeshData class #441

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 5 commits into
base: master
Choose a base branch
from

Conversation

munechika-koyo
Copy link
Contributor

@munechika-koyo munechika-koyo commented Feb 15, 2025

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)
  1. Calculating the volume of a tetrahedron
>>> tetra.volume(0)
0.166666  # == 1 x 1 x 0.5 / 3.0
>>> tetra.volume_total()
0.166666
  1. Calculating the barycenter of a tetrahedron
>>> tetra.barycenter(0)
Point3D(0.25, 0.25, 0.25)
  1. Saving as .rsm
>>> tetra.save("tetra.rsm")
  1. 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
Loading

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.

@munechika-koyo munechika-koyo marked this pull request as ready for review February 15, 2025 12:02
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

1 participant