Add function to voxelise a mesh

examples/voxelise-mesh.py

@@ -0,0 +1,26 @@
+# OpenEP
+# Copyright (c) 2021 OpenEP Collaborators
+#
+# This file is part of OpenEP.
+#
+# OpenEP is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# OpenEP is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program (LICENSE.txt).  If not, see <http://www.gnu.org/licenses/>
+
+import pyvista
+
+import openep
+from openep._datasets.meshes import MESH_2
+
+mesh = pyvista.read(MESH_2)
+voxels = openep.mesh.voxelise(mesh, extract_myocardium=True)
+voxels.plot()

openep/_datasets/meshes.py

@@ -0,0 +1,16 @@
+"""
+Meshes created for testing openep
+=================================
+
+MESH_2 was created from the OpenEP DATASET_2. The largest surface was extracted
+and the surface then smoothed.
+
+"""
+__all__ = [
+    "MESH_2",
+]
+
+from pkg_resources import resource_filename
+
+MESH_2 = resource_filename(__name__,
+                           "VTK/openep_mesh_2_smooth.vtk")

openep/mesh/__init__.py

@@ -8,4 +8,5 @@
     calculate_vertex_path,
     get_free_boundaries,
     repair_mesh,
+    voxelise
 )

openep/mesh/mesh_routines.py

@@ -59,8 +59,10 @@
 """
 
 from attr import attrs
+from typing import Union
 
 import numpy as np
+import scipy.stats
 
 import pyvista
 import pymeshfix
@@ -74,6 +76,7 @@
     "calculate_field_area",
     "calculate_vertex_distance",
     "calculate_vertex_path",
+    "voxelise",
 ]
 
 
@@ -463,3 +466,140 @@ def calculate_vertex_path(
         path = np.array([])
 
     return path
+
+
+def _get_unreferenced_points(mesh):
+    """Determine indices of points not referenced in the triangulation"""
+
+    indices = np.arange(mesh.n_points)
+    referenced_indices = np.unique(mesh.faces.reshape(mesh.n_faces, 4)[:, 1:].ravel())
+    unreferenced_indices = np.isin(indices, referenced_indices, assume_unique=True, invert=True)
+
+    return unreferenced_indices
+
+
+def _determine_voxel_bins(mesh, edge_length, border=10):
+    """Determine the bins to voxelise a mesh.
+
+    Args:
+        mesh (pyvista.PolyData): Mesh to be voxelised.
+        edge_length (float): Edge length of each voxel, in mm.
+        border (float): Minimum border around the mesh, in mm. 
+
+    Returns:
+        (np.ndarray): Bins that can be used to construct arrays of voxels.
+
+    """
+
+    def _round_up(value, nearest):
+        return np.ceil(value / nearest) * nearest
+
+    def _round_down(value, nearest):
+        return np.floor(value / nearest) * nearest
+
+    low_values = np.asarray(mesh.bounds[::2])
+    high_values = np.asarray(mesh.bounds[1::2])
+
+    low_values = _round_down(low_values, nearest=border) - border / 2
+    high_values = _round_up(high_values, nearest=border) + border / 2
+
+    x_low, y_low, z_low = low_values
+    x_high, y_high, z_high = high_values
+
+    x_bins = np.arange(x_low, x_high + edge_length, edge_length)
+    y_bins = np.arange(y_low, y_high + edge_length, edge_length)
+    z_bins = np.arange(z_low, z_high + edge_length, edge_length)
+
+    bin_edges = [x_bins, y_bins, z_bins]
+    bin_centres = [
+        x_bins[:-1] + edge_length / 2,
+        y_bins[:-1] + edge_length / 2,
+        z_bins[:-1] + edge_length / 2,
+    ]
+
+    return bin_edges, bin_centres
+
+
+def voxelise(
+    mesh: pyvista.PolyData,
+    thickness: Union[float, np.ndarray] = 2,
+    n_surfaces: int = 11,
+    edge_length: float = 1,
+    extract_myocardium: bool = False,
+) -> pyvista.PolyData:
+    """Voxelise a surface mesh.
+
+    Args:
+        mesh (PolyData): Surface mesh to be voxelised.
+        thickness (float or np.ndarray): If a float, this defines to thickness of the myocardium.
+            An array of thicknesses - one per point in the mesh - can be passed to create a voxelised
+            mesh with heterogenous thickness.
+        edge_length (float): Length of the voxel edges, in mm.
+        n_surfaces (int): A series of surface meshes are created by interpolating points between the
+            endocardium and epicardium. Points from this series of meshes are used to determine which voxels
+            should be filled. The smaller the voxel edge length, the larger :attr:`n_surfaces`
+            should be.
+        extract_myocardium (bool, optional): If True the voxelised myocardium will be extracted and
+            returned. If False, the voxels in a StructuredGrid will be labelled as filled (1) or empty (0),
+            and this data stored as point data in the returned mesh.
+
+    Returns:
+        StructuredGrid: The voxelised mesh.
+    """
+
+    # Don't make any changes to the mesh
+    mesh = mesh.copy(deep=True)
+
+    # Remove points not referenced by the triangulation
+    not_referenced = _get_unreferenced_points(mesh)
+    mesh.remove_points(not_referenced, inplace=True)
+
+    # Compute normals and set thicknesses
+    mesh.compute_normals(inplace=True, auto_orient_normals=True, cell_normals=False, point_normals=True)
+    thickness = np.full(mesh.n_points, fill_value=thickness, dtype=float) if isinstance(thickness, float) else thickness
+    mesh.point_data['Thickness'] = thickness
+
+    # Calculate voxel bins and create output mesh
+    bin_edges, bin_centres = _determine_voxel_bins(mesh, edge_length=edge_length)
+    bin_centres_x, bin_centres_y, bin_centres_z = bin_centres
+
+    XX, YY, ZZ = np.meshgrid(bin_centres_x, bin_centres_y, bin_centres_z, indexing='xy')  # must use bin centres, must use Cartesian indexing
+    voxels = pyvista.StructuredGrid(XX, YY, ZZ)
+
+    voxel_filled = np.zeros(voxels.n_points, dtype=int)  # keep track of which voxels are filled. 0: empty, 1: filled
+
+    # Determine how much to 
+    n_voxels_x = np.ptp(bin_centres_x)
+    n_voxels_y = np.ptp(bin_centres_y)
+    n_voxels_z = np.ptp(bin_centres_z)
+
+    # Create a series - of open meshes and voxelise each mesh
+    for shell_distance in np.linspace(0, 1, n_surfaces):
+
+        shell = mesh.copy(deep=True)
+        shell.points += mesh.point_data['Normals'] * mesh.point_data['Thickness'][:, np.newaxis] * shell_distance
+        shell.subdivide_adaptive(max_edge_len=edge_length, inplace=True)
+
+        mesh_binned = scipy.stats.binned_statistic_dd(
+            sample=np.asarray(shell.points),
+            values=np.zeros(shell.n_points),
+            statistic='count',
+            bins=bin_edges,
+            expand_binnumbers=True,
+        )
+
+        x_indices, y_indices, z_indices = mesh_binned.binnumber - 1
+        bin_indices = np.ravel_multi_index(
+            [y_indices, x_indices, z_indices],  # we're using Cartesian indexing (as does PyVista)
+            dims=np.asarray([n_voxels_y+1, n_voxels_x+1, n_voxels_z+1], dtype=int),
+            order='F',
+        )
+
+        voxel_filled[bin_indices] = 1
+
+    voxels.point_data['Filled'] = voxel_filled
+
+    if extract_myocardium:
+        voxels = voxels.extract_points(voxel_filled.astype(bool))
+
+    return voxels