import pyvista as pv
import os
import numpy as np
import pandas as pd
import subprocess
import sys
matlabPath = "/Applications/MATLAB_R2023b.app/bin/matlab"
comsolPath = "/Applications/COMSOL61/Multiphysics/bin/comsol"
[docs]
def find_common_dir(start_dir, target_dir_name):
"""
Find the target directory by searching upwards from the start directory.
This function searches recursively upwards from a specified starting directory
to find a directory with the specified name.
Parameters:
start_dir (str): The directory path from which to start the search.
target_dir_name (str): The name of the directory to find.
Returns:
str: The absolute path to the target directory if found.
Raises:
FileNotFoundError: If the target directory is not found.
"""
current_dir = start_dir
while current_dir != os.path.dirname(current_dir): # Stop when reaching the root
if target_dir_name in os.listdir(current_dir):
return os.path.join(current_dir, target_dir_name)
current_dir = os.path.dirname(current_dir)
raise FileNotFoundError(
f"Could not find the directory '{target_dir_name}' from the start directory '{start_dir}'."
)
# dir_path = os.path.dirname(os.path.realpath(__file__))
# common_dir =find_common_dir(dir_path, 'T00_Common')
# sys.path.insert(0, os.path.join(common_dir, 'py/mesh-importer_v2'))
# import meshConversion as mConv
[docs]
class myMesh:
def __init__(self, paths, **kwargs):
"""
Initialize the myMesh object.
This constructor initializes the `myMesh` object by loading the specified meshes
from the given paths. It also prepares the object for further processing, such as
translation and scaling of the meshes.
Parameters:
paths (dict): A dictionary containing paths to the mesh files. Keys should include
'wall', 'calc', 'fibr', 'lrnc', and 'out' (for the output directory).
**kwargs: Additional keyword arguments.
Attributes:
loaded (dict): A dictionary holding the loaded meshes.
translated (dict): A dictionary holding the translated meshes.
scaled (dict): A dictionary holding the scaled meshes.
outPath (str): The output path where files will be saved.
script_dir (str): The directory path where the script is located.
common_dir (str): The path to the 'T00_Common' directory.
wls_script_path (str): The path to the Wolfram Language script.
"""
# Initialize the loaded attribute as a dictionary
self.loaded = {}
# Check if each path exists and is valid before reading the mesh
if "wall" in paths and os.path.exists(paths["wall"]):
self.loaded["wall"] = pv.read(paths["wall"])
if "calc" in paths and os.path.exists(paths["calc"]):
self.loaded["calc"] = pv.read(paths["calc"])
if "fibr" in paths and os.path.exists(paths["fibr"]):
self.loaded["fibr"] = pv.read(paths["fibr"])
if "lrnc" in paths and os.path.exists(paths["lrnc"]):
self.loaded["lrnc"] = pv.read(paths["lrnc"])
self.set_out_path(paths)
self.translated = {}
self.scaled = {}
self.translate_meshes_to_origin()
self.script_dir = dir_path
# Find the path to the T00_Common directory
self.common_dir = common_dir
# Construct the path to the .wls script in the T00_Common directory
self.wls_script_path = os.path.join(
self.common_dir, "generate-interpolation-materials.wls"
)
[docs]
def set_out_path(self, paths):
"""
Set the output path for saving the material map files.
This function checks if the output directory exists and creates it if necessary.
It then sets the `outPath` attribute to the specified output path.
Parameters:
paths (dict): A dictionary containing the 'out' key with the output directory path.
"""
mesh_folder_path = os.path.join(paths["out"], "mesh")
# Check if the folder exists, and create it if it doesn't
if not os.path.exists(mesh_folder_path):
os.makedirs(mesh_folder_path)
self.outPath = mesh_folder_path
[docs]
def translate_meshes_to_origin(self):
"""
Translate a list of meshes so that the smallest coordinate for each axis starts from zero.
This function translates the loaded meshes so that the smallest coordinates along
the x, y, and z axes are shifted to the origin (0, 0, 0). The translated meshes are
stored in the `translated` attribute.
Returns:
list: A list of translated meshes.
"""
global_min_coords = np.array([np.inf, np.inf, np.inf])
# Get the list of loaded meshes
meshes = [mesh for mesh in self.loaded.values()]
# Find the global minimum coordinates
for mesh in meshes:
bounds = mesh.bounds
min_coords = np.array([bounds[0], bounds[2], bounds[4]])
global_min_coords = np.minimum(global_min_coords, min_coords)
# Translate all meshes by the negative of the global minimum coordinates
translated_meshes = []
for mesh in meshes:
translated_mesh = mesh.translate(-global_min_coords)
translated_meshes.append(translated_mesh)
# Store the translated meshes back into self.translated dictionary
for key, translated_mesh in zip(self.loaded.keys(), translated_meshes):
self.translated[key] = translated_mesh
return translated_meshes
[docs]
def scale_meshes(self, scale_factor):
"""
Scale all loaded meshes by the specified scale factor and store them in self.scaled.
This function applies a uniform scaling transformation to all loaded meshes based
on the provided scale factor. The scaled meshes are stored in the `scaled` attribute.
Parameters:
scale_factor (float): The factor by which to scale the meshes.
"""
for key, mesh in self.translated.items():
scaled_mesh = mesh.copy()
scaled_mesh.scale([scale_factor, scale_factor, scale_factor], inplace=True)
self.scaled[key] = scaled_mesh
# check that they exist or assign them as empty meshes
if "wall" not in self.scaled:
self.scaled["wall"] = pv.PolyData()
if "calc" not in self.scaled:
self.scaled["calc"] = pv.PolyData()
if "fibr" not in self.scaled:
self.scaled["fibr"] = pv.PolyData()
if "lrnc" not in self.scaled:
self.scaled["lrnc"] = pv.PolyData()
[docs]
def visualize_mesh(self, type="scaled"):
"""
Visualize the meshes in the self.scaled dictionary.
This function creates a 3D plot of the meshes. The meshes can be visualized in their
original, translated, or scaled state based on the `type` parameter.
Parameters:
type (str): The type of meshes to visualize ('scaled', 'translated', or 'loaded').
"""
if type == "scaled":
meshes = self.scaled
elif type == "translated":
meshes = self.translated
else:
meshes = self.loaded
plotter = pv.Plotter()
# Check and add 'wall' mesh if it exists
if "wall" in meshes:
wall = meshes["wall"]
plotter.add_mesh(wall, color="red", opacity=0.3, label="Wall")
# Check and add 'calc' mesh if it exists
if "calc" in meshes:
calc = meshes["calc"]
plotter.add_mesh(calc, color="blue", opacity=0.3, label="Calcio")
# Check and add 'fibr' mesh if it exists
if "fibr" in meshes:
fibr = meshes["fibr"]
plotter.add_mesh(fibr, color="green", opacity=0.3, label="Fibrosa")
# Check and add 'lrnc' mesh if it exists
if "lrnc" in meshes:
lrnc = meshes["lrnc"]
plotter.add_mesh(lrnc, color="yellow", opacity=0.3, label="Lipid")
# Add legend and show the plot
plotter.add_legend()
plotter.show()
[docs]
def assign_material_map(self, spacing=1.0):
"""
Assign material IDs to each point in the spatial grid based on the mesh it is inside.
This function creates a structured grid of points within the bounds of the wall mesh
and checks which points are inside each of the provided meshes. The result is stored
in the `material_points` attribute.
Parameters:
spacing (float): The spacing between points in the grid. Default is 1.0.
"""
if self.scaled == {}:
self.scale_meshes(1.0)
bounds = self.scaled["wall"].bounds
grid = create_spatial_grid(bounds, spacing)
self.points_grid = grid
# check and define mesh or assign them as an empty mesh
wall = self.scaled["wall"] if "wall" in self.scaled else pv.PolyData()
calc = self.scaled["calc"] if "calc" in self.scaled else pv.PolyData()
fibr = self.scaled["fibr"] if "fibr" in self.scaled else pv.PolyData()
lrnc = self.scaled["lrnc"] if "lrnc" in self.scaled else pv.PolyData()
inside_points_wall = check_points_inside_mesh(wall, grid)
inside_points_lrnc = check_points_inside_mesh(lrnc, grid)
inside_points_calc = check_points_inside_mesh(calc, grid)
inside_points_fibr = check_points_inside_mesh(fibr, grid)
# Set to zero all points in inside_points_wall which are also inside any of the other three meshes
overlapping_points_wall = (
(inside_points_lrnc == 1)
| (inside_points_calc == 1)
| (inside_points_fibr == 1)
)
inside_points_wall[overlapping_points_wall] = 0
self.material_points = {}
self.material_points["wall"] = inside_points_wall
self.material_points["lrnc"] = inside_points_lrnc
self.material_points["calc"] = inside_points_calc
self.material_points["fibr"] = inside_points_fibr
[docs]
def visualize_material_map(self, plotter=None):
"""
Visualize the material map for the specified material.
This function creates a 3D plot that visualizes the distribution of material IDs
across the spatial grid. Each material type (e.g., wall, calcification, fibrous, lipid)
is displayed in a separate subplot.
Parameters:
plotter (pyvista.Plotter, optional): A PyVista plotter object for customized plotting.
If None, a new plotter is created.
"""
if self.material_points == {}:
raise ValueError("Invalid material specified")
if plotter is None:
plotter = pv.Plotter()
grid = self.points_grid
plotter = pv.Plotter(shape=(1, 4))
plotter.subplot(0, 0)
plotter = add_meshes_to_plotter(plotter, self.scaled)
plotter.add_mesh(
extract_plot(grid, self.material_points["wall"]),
color="black",
point_size=5,
)
plotter.add_text("WALL", font_size=14, position="upper_edge", color="black")
plotter.subplot(0, 1)
plotter = add_meshes_to_plotter(plotter, self.scaled)
plotter.add_mesh(
extract_plot(grid, self.material_points["calc"]),
color="black",
point_size=5,
)
plotter.add_text("CALC", font_size=14, position="upper_edge", color="black")
plotter.subplot(0, 2)
plotter = add_meshes_to_plotter(plotter, self.scaled)
plotter.add_mesh(
extract_plot(grid, self.material_points["fibr"]),
color="black",
point_size=5,
)
plotter.add_text("FIBR", font_size=14, position="upper_edge", color="black")
plotter.subplot(0, 3)
plotter = add_meshes_to_plotter(plotter, self.scaled)
plotter.add_mesh(
extract_plot(grid, self.material_points["lrnc"]),
color="black",
point_size=5,
)
plotter.add_text("LRNC", font_size=14, position="upper_edge", color="black")
plotter.show()
[docs]
def save_to_hdf(self, material, name, outPath=None):
"""
Save material points data to HDF5 format.
This function saves the material point data for a specific material to an HDF5 file.
Parameters:
material (str): The material for which to save the data (e.g., 'wall', 'calc', 'fibr', 'lrnc').
name (str): The name of the output HDF5 file (without extension).
outPath (str, optional): The directory where the HDF5 file will be saved.
If None, the default output path is used.
"""
if outPath is None:
outPath = self.outPath
points = self.points_grid.points
values = self.material_points[material].flatten(order="F")
data = np.column_stack((points, values))
# Convert the numpy array to a pandas DataFrame
df = pd.DataFrame(data, columns=["x", "y", "z", "value"])
# Save to HDF5 file
h5file_path = os.path.join(outPath, f"{name}.h5")
df.to_hdf(h5file_path, key="data", mode="w", format="table")
[docs]
def write_material_map(self, outPath=None):
"""
Save all material points to HDF5 files and then run the bash script.
This function saves the material point data for all materials ('wall', 'lrnc', 'calc', 'fibr')
to individual HDF5 files and then runs a predefined Wolfram Language script for further processing.
Parameters:
outPath (str, optional): The directory where the HDF5 files will be saved.
If None, the default output path is used.
"""
materials = ["wall", "lrnc", "calc", "fibr"]
names = {"wall": "matx", "lrnc": "lrnc", "calc": "calc", "fibr": "fibr"}
for material in materials:
if material in self.material_points:
self.save_to_hdf(material, names[material], outPath)
# Run the bash script
self.run_bash_script()
[docs]
def run_bash_script(self):
"""
Call the bash script.
This function executes a predefined Wolfram Language script using the `wolframscript`
command-line tool. The script processes the mesh data in the specified output directory.
"""
# Construct the command to run the Wolfram script
mesh_folder = self.outPath
command = ["wolframscript", "-file", self.wls_script_path, mesh_folder]
subprocess.run(command, check=True)
[docs]
def prepare_matlab_script(self, block_height_multiplier, hmax, hmin):
"""
Create a MATLAB script to generate the mesh.
This function generates a MATLAB script that automates the process of creating a mesh
using COMSOL Multiphysics. The script is saved in the specified output directory.
Parameters:
block_height_multiplier (float): Multiplier for the block height used in the mesh.
hmax (float): Maximum mesh element size.
hmin (float): Minimum mesh element size.
"""
# Example usage
create_matlab_script(
script_filename=os.path.join(self.outPath, "script.m"),
common_path=self.common_dir,
model_path=self.outPath,
stl_path=os.path.join(self.outPath, "vessel_wall.stl"),
block_height_multiplier=block_height_multiplier,
hmax=hmax,
hmin=hmin,
)
[docs]
def save_wall_mesh(self):
"""
Save the wall mesh to an STL file.
This function saves the wall mesh to an STL file named 'vessel_wall.stl' in the
specified output directory.
Raises:
ValueError: If no wall mesh is found to save.
"""
if "wall" in self.scaled:
wall = self.scaled["wall"]
wall.save(os.path.join(self.outPath, "vessel_wall.stl"))
self.vessel_wall_path = os.path.join(self.outPath, "vessel_wall.stl")
else:
raise ValueError("No wall mesh found to save")
[docs]
def create_mph_mesh(self, block_height_multiplier, hmax, hmin):
"""
Create the volumetric wall mesh using COMSOL and MATLAB.
This function orchestrates the process of creating a volumetric wall mesh by first
saving the wall mesh, then preparing and executing the required MATLAB script, and
finally processing the output mesh file using COMSOL and custom tools.
"""
self.save_wall_mesh()
self.prepare_matlab_script(block_height_multiplier, hmax, hmin)
comsol_command = [comsolPath, "mphserver", "-port", "2036"]
subprocess.Popen(comsol_command)
script_path = os.path.join(self.outPath, "script.m")
# Construct the MATLAB command
matlab_command = [
matlabPath,
"-batch",
f"run('{script_path}')", # Run the script
]
subprocess.run(matlab_command, check=True)
[docs]
def create_volumetric_wall_mesh(self, block_height_multiplier, hmax, hmin):
"""
Create the volumetric wall mesh using COMSOL and MATLAB.
This function orchestrates the process of creating a volumetric wall mesh by first
saving the wall mesh, then preparing and executing the required MATLAB script, and
finally processing the output mesh file using COMSOL and custom tools.
"""
self.save_wall_mesh()
self.create_mph_mesh(block_height_multiplier, hmax, hmin)
meshPath = os.path.join(self.outPath, "mesh.mphtxt")
mConv.main([meshPath])
h5_file_path = os.path.join(self.outPath, "mesh.h5")
mConv.classify_surfaces(h5_file_path)
return
[docs]
def add_meshes_to_plotter(plotter, meshes):
"""
Add multiple meshes to a PyVista plotter.
This function adds the provided meshes to a given PyVista plotter object for visualization.
Parameters:
plotter (pyvista.Plotter): The PyVista plotter object to which the meshes will be added.
meshes (dict): A dictionary containing the meshes to add, keyed by mesh type.
Returns:
pyvista.Plotter: The updated PyVista plotter object.
"""
plotter.add_mesh(meshes["wall"], color="red", opacity=0.3, label="Wall")
plotter.add_mesh(meshes["calc"], color="white", opacity=0.3, label="Calcio")
plotter.add_mesh(meshes["fibr"], color="white", opacity=0.3, label="Fibrosa")
plotter.add_mesh(meshes["lrnc"], color="white", opacity=0.3, label="Lipid")
return plotter
[docs]
def create_spatial_grid(bounds, spacing=1.0):
"""
Create a spatial grid of points that covers the size of the given bounds.
This function generates a 3D grid of points that spans the specified bounding box,
with points spaced according to the given spacing parameter.
Parameters:
bounds (tuple): The bounds of the region to cover (xmin, xmax, ymin, ymax, zmin, zmax).
spacing (float): The spacing between points in the grid. Default is 1.0.
Returns:
pyvista.StructuredGrid: The spatial grid of points.
"""
x = np.arange(bounds[0], bounds[1] + spacing, spacing)
y = np.arange(bounds[2], bounds[3] + spacing, spacing)
z = np.arange(bounds[4], bounds[5] + spacing, spacing)
xx, yy, zz = np.meshgrid(x, y, z, indexing="ij")
points = np.c_[xx.ravel(), yy.ravel(), zz.ravel()]
grid = pv.StructuredGrid()
grid.points = points
grid.dimensions = (len(x), len(y), len(z))
return grid
[docs]
def check_points_inside_mesh(mesh, grid):
"""
Check if each point in the grid is inside or outside the mesh.
Parameters:
mesh (pyvista.PolyData): The mesh to test against.
grid (pyvista.StructuredGrid): The grid of points to test.
Returns:
inside_points (np.ndarray): An array where 1 indicates the point is inside the mesh, and 0 indicates outside.
"""
point_cloud = pv.PolyData(grid.points)
selected = point_cloud.select_enclosed_points(
mesh, tolerance=0.001, inside_out=False, check_surface=True
)
inside_points = selected.point_data["SelectedPoints"].astype(int)
# Debug prints
# print(f"Number of points in the grid: {grid.points.shape[0]}")
# print(f"Number of points in inside_points array: {inside_points.size}")
# Ensure the length matches the number of grid points
if inside_points.size != grid.points.shape[0]:
raise ValueError(
"Mismatch in the number of points when reshaping inside_points array"
)
return inside_points.reshape(grid.dimensions, order="F")
# def check_points_inside_mesh(mesh, grid, method='select', tolerance=1e-6):
# """
# Check if each point in the grid is inside or outside the mesh.
# Parameters:
# mesh (pyvista.PolyData): The mesh to test against.
# grid (pyvista.StructuredGrid): The grid of points to test.
# method (str): The method to use for checking points ('select' or 'contains').
# tolerance (float): The tolerance for the 'select' method.
# Returns:
# inside_points (np.ndarray): An array where 1 indicates the point is inside the mesh, and 0 indicates outside.
# """
# if method == 'select':
# point_cloud = pv.PolyData(grid.points)
# selected = point_cloud.select_enclosed_points(mesh, tolerance=tolerance, inside_out=False, check_surface=True)
# inside_points = selected.point_data['SelectedPoints'].astype(int)
# elif method == 'contains':
# inside_points = np.array([mesh.contains(point) for point in grid.points], dtype=int)
# else:
# raise ValueError(f"Unknown method: {method}")
# # Ensure the length matches the number of grid points
# if inside_points.size != grid.points.shape[0]:
# raise ValueError("Mismatch in the number of points when reshaping inside_points array")
# return inside_points.reshape(grid.dimensions, order='F')
[docs]
def join_relative_paths(file_dir, relative_paths):
"""
Join the script directory with the relative paths.
This function joins the provided script directory with the relative paths of input files
and ensures that all paths are correctly formed.
Parameters:
file_dir (str): The directory where the script is located.
relative_paths (dict): A dictionary of relative paths.
Returns:
dict: A dictionary containing the absolute paths of the input files.
"""
# Join the script directory with the relative paths
paths = {
key: os.path.join(file_dir, value) if value else ""
for key, value in relative_paths.items()
}
if not paths["out"]:
paths["out"] = file_dir
for key, value in paths.items():
print(f"'{key}': '{value}'")
return paths
[docs]
def create_matlab_script(
script_filename,
common_path,
model_path,
stl_path,
block_height_multiplier,
hmax,
hmin,
):
"""
Create a MATLAB script to generate the mesh.
This function generates a MATLAB script that automates the creation of a mesh using
COMSOL Multiphysics, including setting the necessary parameters for the mesh generation.
Parameters:
script_filename (str): The path where the MATLAB script will be saved.
common_path (str): The path to the common directory containing necessary scripts.
model_path (str): The path to the model directory.
stl_path (str): The path to the STL file of the mesh.
block_height_multiplier (float): The multiplier for the block height.
hmax (str): The maximum element size for meshing.
hmin (str): The minimum element size for meshing.
"""
script_content = f"""\
addpath('{common_path}');
modelPath = '{model_path}';
blockHeightMultiplier = {block_height_multiplier};
hmax = '{hmax}';
hmin = '{hmin}';
% Set the STL file path
stl_path = '{stl_path}';
% Call the function defined in Untitled2.m
createMPHmesh(modelPath, stl_path, blockHeightMultiplier, hmax, hmin);
% Optionally, remove the path afterward to keep the environment clean
rmpath('{common_path}');
"""
# Write the content to the script file
with open(script_filename, "w") as file:
file.write(script_content)
print(f"MATLAB script saved as {script_filename}")