From ce038803a495d7fa18eac0d90d664e061aad3b97 Mon Sep 17 00:00:00 2001 From: Henrik Kjedsberg Date: Mon, 6 May 2024 01:26:08 +0200 Subject: [PATCH] Add geodesic refinement as meshing method and test --- .../automated_preprocessing.py | 60 ++++++++++-------- .../preprocessing_common.py | 59 ++++++++++++++++- tests/test_pre_processing.py | 63 ++++++++++++++++++- 3 files changed, 151 insertions(+), 31 deletions(-) diff --git a/src/vampy/automatedPreprocessing/automated_preprocessing.py b/src/vampy/automatedPreprocessing/automated_preprocessing.py index 87f04a66..6ed0a3a6 100644 --- a/src/vampy/automatedPreprocessing/automated_preprocessing.py +++ b/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,7 +25,7 @@ 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. @@ -32,6 +33,7 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f 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(): diff --git a/src/vampy/automatedPreprocessing/preprocessing_common.py b/src/vampy/automatedPreprocessing/preprocessing_common.py index a84eb0b4..a2a0d14a 100644 --- a/src/vampy/automatedPreprocessing/preprocessing_common.py +++ b/src/vampy/automatedPreprocessing/preprocessing_common.py @@ -1,5 +1,7 @@ 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, \ @@ -7,7 +9,7 @@ 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 diff --git a/tests/test_pre_processing.py b/tests/test_pre_processing.py index 3b9a5b3a..24045c11 100644 --- a/tests/test_pre_processing.py +++ b/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()