Add get_mesh_from_probability_map

More accurate sub-voxel mesh surface interpolation and reduces vtk-dep.

oai_analysis/mesh_processing.py

@@ -339,6 +339,15 @@ def get_mesh(itk_image, num_iterations=150):
     mesh = smooth_mesh(mesh, num_iterations=num_iterations)
     return mesh
 
+# To obtain mesh from itk image probability map
+def get_mesh_from_probability_map(image: itk.Image) -> itk.Mesh:
+    mesh = itk.cuberille_image_to_mesh_filter(image,
+                                              generate_triangle_faces=True,
+                                              iso_surface_value=0.5,
+                                              project_vertices_to_iso_surface=True,
+                                              project_vertex_surface_distance_threshold=0.05)
+    return mesh
+
 
 # To obtain inner and outer mesh splits given the mesh type
 def split_mesh(mesh, mesh_type="FC"):

pyproject.toml

@@ -29,6 +29,7 @@ classifiers = [
 dependencies = [
   "torch",
   "itk==5.3",
+  "itk-cuberille",
   "trimesh",
   "vtk",
   "icon_registration==0.3.4",

test/test_mesh_processing.py

@@ -4,9 +4,23 @@
 import trimesh
 from pathlib import Path
 import numpy as np
+import itk
 
 vtk = pytest.importorskip("vtk")
 
+def test_get_mesh():
+    test_prob_filepath = Path(__file__).parent / "test_files" / "colab_case" / "TC_probmap.nii.gz"
+    image = itk.imread(test_prob_filepath)
+    mesh = mp.get_mesh_from_probability_map(image)
+
+    baseline_mesh_filepath = Path(__file__).parent / "test_files" / "colab_case" / "TC_mesh.vtk"
+    baseline = itk.meshread(baseline_mesh_filepath)
+
+    trimesh = mp.get_trimesh(mesh)
+    baseline_trimesh = mp.get_trimesh(baseline)
+    np.testing.assert_allclose(trimesh.vertices, baseline_trimesh.vertices, atol=0.02)
+
+
 def test_get_cell_normals():
     test_filepath = Path(__file__).parent / "test_files" / "colab_case" / "avsm" / "TC_mesh_world.ply"