Merge pull request #149 from KVSlab/feature/geodesic-refinement

Add refinement based on Geodesic distance

src/vampy/automatedPreprocessing/automated_preprocessing.py

@@ -14,7 +14,8 @@
 from vampy.automatedPreprocessing.preprocessing_common import get_centers_for_meshing, \
     dist_sphere_diam, dist_sphere_curvature, dist_sphere_constant, get_regions_to_refine, add_flow_extension, \
     write_mesh, mesh_alternative, generate_mesh, find_boundaries, compute_flow_rate, setup_model_network, \
-    radiusArrayName, scale_surface, get_furtest_surface_point, check_if_closed_surface, remesh_surface
+    radiusArrayName, scale_surface, get_furtest_surface_point, check_if_closed_surface, remesh_surface, \
+    dist_sphere_geodesic
 from vampy.automatedPreprocessing.repair_tools import find_and_delete_nan_triangles, clean_surface, print_surface_info
 from vampy.automatedPreprocessing.simulate import run_simulation
 from vampy.automatedPreprocessing.visualize import visualize_model
@@ -24,14 +25,15 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f
                        meshing_method, refine_region, is_atrium, add_flow_extensions, visualize, config_path,
                        coarsening_factor, inlet_flow_extension_length, outlet_flow_extension_length, edge_length,
                        region_points, compress_mesh, add_boundary_layer, scale_factor, resampling_step,
-                       flow_rate_factor, moving_mesh, clamp_boundaries):
+                       flow_rate_factor, moving_mesh, clamp_boundaries, max_geodesic_distance):
     """
     Automatically generate mesh of surface model in .vtu and .xml format, including prescribed
     flow rates at inlet and outlet based on flow network model.
 
     Runs simulation of meshed case on a remote ssh server if server configuration is provided.
 
     Args:
+        max_geodesic_distance (float): Max distance when performing geodesic refinement (in mm)
         input_model (str): Name of case
         verbose_print (bool): Toggles verbose mode
         smoothing_method (str): Method for surface smoothing
@@ -66,9 +68,7 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f
     file_name_centerlines = base_path + "_centerlines.vtp"
     file_name_refine_region_centerlines = base_path + "_refine_region_centerline.vtp"
     file_name_region_centerlines = base_path + "_sac_centerline_{}.vtp"
-    file_name_distance_to_sphere_diam = base_path + "_distance_to_sphere_diam.vtp"
-    file_name_distance_to_sphere_const = base_path + "_distance_to_sphere_const.vtp"
-    file_name_distance_to_sphere_curv = base_path + "_distance_to_sphere_curv.vtp"
+    file_name_distance_to_sphere = base_path + f"_distance_to_sphere_{meshing_method}.vtp"
     file_name_distance_to_sphere_initial = base_path + "_distance_to_sphere_initial.vtp"
     file_name_probe_points = base_path + "_probe_point.json"
     file_name_voronoi = base_path + "_voronoi.vtp"
@@ -169,10 +169,10 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f
         # Extract the region centerline
         refine_region_centerline = []
         info = get_parameters(base_path)
-        num_anu = info["number_of_regions"]
+        number_of_refinement_regions = info["number_of_regions"]
 
         # Compute mean distance between points
-        for i in range(num_anu):
+        for i in range(number_of_refinement_regions):
             if not path.isfile(file_name_region_centerlines.format(i)):
                 line = extract_single_line(centerline_region, i)
                 locator = get_vtk_point_locator(centerlines)
@@ -193,7 +193,7 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f
             else:
                 refine_region_centerline.append(read_polydata(file_name_region_centerlines.format(i)))
 
-        # Merge the sac centerline
+        # Merge the refined region centerline
         region_centerlines = vtk_merge_polydata(refine_region_centerline)
 
         for region in refine_region_centerline:
@@ -354,25 +354,27 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f
 
     # Choose input for the mesh
     print("--- Computing distance to sphere\n")
-    if meshing_method == "constant":
-        if not path.isfile(file_name_distance_to_sphere_const):
+    if path.isfile(file_name_distance_to_sphere):
+        distance_to_sphere = read_polydata(file_name_distance_to_sphere)
+    else:
+        if meshing_method == "constant":
             distance_to_sphere = dist_sphere_constant(surface_extended, centerlines, region_center, misr_max,
-                                                      file_name_distance_to_sphere_const, edge_length)
-        else:
-            distance_to_sphere = read_polydata(file_name_distance_to_sphere_const)
-
-    elif meshing_method == "curvature":
-        if not path.isfile(file_name_distance_to_sphere_curv):
+                                                      file_name_distance_to_sphere, edge_length)
+        elif meshing_method == "curvature":
             distance_to_sphere = dist_sphere_curvature(surface_extended, centerlines, region_center, misr_max,
-                                                       file_name_distance_to_sphere_curv, coarsening_factor)
-        else:
-            distance_to_sphere = read_polydata(file_name_distance_to_sphere_curv)
-    elif meshing_method == "diameter":
-        if not path.isfile(file_name_distance_to_sphere_diam):
+                                                       file_name_distance_to_sphere, coarsening_factor)
+        elif meshing_method == "diameter":
             distance_to_sphere = dist_sphere_diam(surface_extended, centerlines, region_center, misr_max,
-                                                  file_name_distance_to_sphere_diam, coarsening_factor)
+                                                  file_name_distance_to_sphere, coarsening_factor)
+        elif meshing_method == "geodesic":
+            if edge_length is None:
+                print("Edge length needs to supplied when using Geodesic meshing")
+                sys.exit(0)
+            distance_to_sphere = dist_sphere_geodesic(surface_extended, region_center, max_geodesic_distance,
+                                                      file_name_distance_to_sphere, edge_length)
         else:
-            distance_to_sphere = read_polydata(file_name_distance_to_sphere_diam)
+            print(f"Method '{meshing_method}' is not a valid meshing method")
+            sys.exit(0)
 
     # Compute mesh
     if not path.isfile(file_name_vtu_mesh):
@@ -421,10 +423,9 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f
     print("--- Removing unused pre-processing files")
     files_to_remove = [
         file_name_centerlines, file_name_refine_region_centerlines, file_name_region_centerlines,
-        file_name_distance_to_sphere_diam, file_name_distance_to_sphere_const, file_name_distance_to_sphere_curv,
+        file_name_distance_to_sphere, file_name_remeshed, file_name_distance_to_sphere_initial,
         file_name_voronoi, file_name_voronoi_smooth, file_name_voronoi_surface, file_name_surface_smooth,
         file_name_model_flow_ext, file_name_clipped_model, file_name_flow_centerlines, file_name_surface_name,
-        file_name_remeshed, file_name_distance_to_sphere_initial
     ]
     for file in files_to_remove:
         if path.exists(file):
@@ -494,7 +495,7 @@ def read_command_line(input_path=None):
 
     parser.add_argument('-m', '--meshing-method',
                         type=str,
-                        choices=["diameter", "curvature", "constant"],
+                        choices=["diameter", "curvature", "constant", "geodesic"],
                         default="diameter",
                         help="Determines method of meshing. The method 'constant' is supplied with a constant edge " +
                              "length controlled by the -el argument, resulting in a constant density mesh. " +
@@ -587,6 +588,11 @@ def read_command_line(input_path=None):
                         default=False,
                         help="Clamps boundaries at inlet(s) and outlet(s) if true. Only used for moving mesh.")
 
+    parser.add_argument('-gd', '--max-geodesic-distance',
+                        default=10,
+                        type=float,
+                        help="Maximum distance when performing geodesic distance. In [mm].")
+
     # Parse path to get default values
     if required:
         args = parser.parse_args()
@@ -623,7 +629,7 @@ def verbose_print(*args):
                 outlet_flow_extension_length=args.outlet_flowextension, add_boundary_layer=args.add_boundary_layer,
                 scale_factor=args.scale_factor, resampling_step=args.resampling_step,
                 flow_rate_factor=args.flow_rate_factor, moving_mesh=args.moving_mesh,
-                clamp_boundaries=args.clamp_boundaries)
+                clamp_boundaries=args.clamp_boundaries, max_geodesic_distance=args.max_geodesic_distance)
 
 
 def main_meshing():

src/vampy/automatedPreprocessing/preprocessing_common.py

@@ -1,13 +1,15 @@
 import gzip
 import json
+from os import path
+
 import numpy as np
 import vtkmodules.numpy_interface.dataset_adapter as dsa
 from morphman import vtk_clean_polydata, vtk_triangulate_surface, get_parameters, write_parameters, read_polydata, \
     vmtkscripts, vtk, write_polydata, vtkvmtk, get_curvilinear_coordinate, vtk_compute_threshold, get_vtk_array, \
     get_distance, get_number_of_arrays, vmtk_smooth_surface, get_point_data_array, create_vtk_array, \
     get_vtk_point_locator, vtk_extract_feature_edges, get_uncapped_surface, vtk_compute_connectivity, \
     vtk_compute_mass_properties, extract_single_line, get_centerline_tolerance
-from os import path
+
 from vampy.automatedPreprocessing import ImportData
 from vampy.automatedPreprocessing.NetworkBoundaryConditions import FlowSplitting
 from vampy.automatedPreprocessing.vmtk_pointselector import vmtkPickPointSeedSelector
@@ -16,6 +18,7 @@
 distanceToSpheresArrayName = "DistanceToSpheres"
 radiusArrayName = 'MaximumInscribedSphereRadius'
 cellEntityArrayName = "CellEntityIds"
+dijkstraArrayName = "DijkstraDistanceToPoints"
 
 
 def get_regions_to_refine(surface, provided_points, dir_path):
@@ -426,7 +429,7 @@ def dist_sphere_diam(surface, centerlines, region_center, misr_max, save_path, f
     distance_to_sphere.GetPointData().AddArray(elements_vtk)
     element_size = diameter_array / element_size
 
-    # Reduce element size in aneurysm
+    # Reduce element size in refinement region
     for i in range(len(region_center)):
         distance_to_sphere = compute_distance_to_sphere(distance_to_sphere,
                                                         region_center[i],
@@ -922,3 +925,55 @@ def remesh_surface(surface, edge_length, element_size_mode="edgelength", exclude
     remeshed_surface = remeshing.Surface
 
     return remeshed_surface
+
+
+def dist_sphere_geodesic(surface, region_center, max_distance, save_path, edge_length):
+    """
+        Determines the target edge length for each cell on the surface, including
+        potential refinement or coarsening of certain user specified areas.
+        Level of refinement/coarseness is determined based on user selected region and the geodesic distance from it.
+        Args:
+            surface (vtkPolyData): Input surface model
+            region_center (list): Point representing region to refine
+            save_path (str): Location to store processed surface
+            edge_length (float): Target edge length
+            max_distance (float): Max distance of geodesic measure
+        Returns:
+            surface (vtkPolyData): Processed surface model with info on cell specific target edge length
+        """
+
+    # Define refine region point ID
+    locator = get_vtk_point_locator(surface)
+    point = region_center[0]
+    region_id = locator.FindClosestPoint(point)
+    region_ids = vtk.vtkIdList()
+    region_ids.InsertNextId(region_id)
+
+    # Compute geodesic distance to point
+    dijkstra = vtkvmtk.vtkvmtkPolyDataDijkstraDistanceToPoints()
+    dijkstra.SetInputData(surface)
+    dijkstra.SetSeedIds(region_ids)
+    dijkstra.SetDistanceOffset(0)
+    dijkstra.SetDistanceScale(1)
+    dijkstra.SetMinDistance(0)
+    dijkstra.SetMaxDistance(max_distance)
+    dijkstra.SetDijkstraDistanceToPointsArrayName(dijkstraArrayName)
+    dijkstra.Update()
+    geodesic_distance = dijkstra.GetOutput()
+
+    # Create smooth transition between LAA and LA lumen
+    N = surface.GetNumberOfPoints()
+    dist_array = geodesic_distance.GetPointData().GetArray(dijkstraArrayName)
+    for i in range(N):
+        dist = dist_array.GetComponent(i, 0) / max_distance
+        newDist = 1 / 3 * edge_length * (1 + 2 * dist ** 3)
+        dist_array.SetComponent(i, 0, newDist)
+
+    # Set element size based on geodesic distance
+    element_size = get_point_data_array(dijkstraArrayName, geodesic_distance)
+    vtk_array = create_vtk_array(element_size, "Size")
+    geodesic_distance.GetPointData().AddArray(vtk_array)
+
+    write_polydata(geodesic_distance, save_path)
+
+    return geodesic_distance

tests/test_pre_processing.py

@@ -1,5 +1,8 @@
-from dolfin import Mesh
+import os
 from os import path
+
+from dolfin import Mesh
+
 from vampy.automatedPreprocessing.automated_preprocessing import read_command_line, \
     run_pre_processing
 from vampy.automatedPreprocessing.preprocessing_common import read_polydata
@@ -39,6 +42,8 @@ def test_mesh_model_with_one_inlet():
     assert mesh_vtu.GetNumberOfPoints() == num_points
     assert mesh_xml.num_cells() == num_cells
 
+    tear_down(model_path)
+
 
 def test_mesh_model_with_one_inlet_and_two_outlets():
     model_path = "tests/test_data/artery/artery.stl"
@@ -53,7 +58,6 @@ def test_mesh_model_with_one_inlet_and_two_outlets():
                              outlet_flow_extension_length=1,
                              inlet_flow_extension_length=1
                              ))
-
     # Run pre processing
     run_pre_processing(**common_input)
 
@@ -74,6 +78,8 @@ def test_mesh_model_with_one_inlet_and_two_outlets():
     assert mesh_vtu.GetNumberOfPoints() == num_points
     assert mesh_xml.num_cells() == num_cells
 
+    tear_down(model_path)
+
 
 def test_mesh_model_with_one_inlet_and_one_outlet():
     model_path = "tests/test_data/vein/vein.stl"
@@ -108,8 +114,61 @@ def test_mesh_model_with_one_inlet_and_one_outlet():
     assert mesh_vtu.GetNumberOfPoints() == num_points
     assert mesh_xml.num_cells() == num_cells
 
+    tear_down(model_path)
+
+
+def test_mesh_model_with_geodesic_meshing():
+    model_path = "tests/test_data/artery/artery.stl"
+    # Get default input parameters
+    common_input = read_command_line(model_path)
+    region_point = [33.6612, 32.8443, 40.9184]
+    common_input.update(dict(meshing_method="geodesic",
+                             max_geodesic_distance=2.5,
+                             edge_length=0.5,
+                             refine_region=True,
+                             region_points=region_point,
+                             visualize=False,
+                             compress_mesh=False,
+                             outlet_flow_extension_length=1,
+                             inlet_flow_extension_length=1
+                             ))
+
+    # Run pre processing
+    run_pre_processing(**common_input)
+
+    # Check that mesh is created
+    mesh_path_vtu = model_path.replace("stl", "vtu")
+    mesh_path_xml = model_path.replace("stl", "xml")
+
+    assert path.isfile(mesh_path_vtu)
+    assert path.isfile(mesh_path_xml)
+
+    # Check that mesh is not empty with VTK/morphMan and FEniCS and contains correct amount of points and cells
+    mesh_vtu = read_polydata(mesh_path_vtu)
+    mesh_xml = Mesh(mesh_path_xml)
+
+    num_points = 6261
+    num_cells = 34849
+
+    assert mesh_vtu.GetNumberOfPoints() == num_points
+    assert mesh_xml.num_cells() == num_cells
+
+    tear_down(model_path)
+
+    # Remove additional file when performing refinement
+    os.remove(model_path.replace(".stl", "_sac_centerline_0.vtp"))
+
+
+def tear_down(model_path):
+    file_extensions = [".vtu", ".xml", "_info.json", "_probe_point.json"]
+    model_name = model_path.split(".")[0]
+    for file_extension in file_extensions:
+        file_to_remove = model_name + file_extension
+        os.remove(file_to_remove)
+
 
 if __name__ == "__main__":
     test_mesh_model_with_one_inlet()
     test_mesh_model_with_one_inlet_and_one_outlet()
     test_mesh_model_with_one_inlet_and_two_outlets()
+    test_mesh_model_with_geodesic_meshing()