How-to: Use flowing particles as simulation input#
This notebook shows how velocity field used for particle simualtion is handled.
Import modules and set TF-GPU configuration#
[2]:
# Load TF and check for GPUs
import os
os.environ["TF_CPP_MIN_LOG_LEVEL"] = "2"
import tensorflow as tf
physical_devices = tf.config.list_physical_devices('GPU')
if physical_devices:
print("Available GPUS: \n\t", "\n\t".join([str(_) for _ in physical_devices]))
tf.config.experimental.set_visible_devices(physical_devices[0], 'GPU')
tf.config.experimental.set_memory_growth(physical_devices[0], True)
import base64
from IPython.display import display, HTML
import cmrseq
from pint import Quantity
import pyvista
import numpy as np
from tqdm.notebook import tqdm
import matplotlib.pyplot as plt
%matplotlib widget
# Project library cmrsim
import sys
sys.path.insert(0, "../../")
import cmrsim
import cmrsim.utils.particle_properties as part_factory
sys.path.append("..")
import local_functions
2024-02-26 16:17:41.932979: E external/local_xla/xla/stream_executor/cuda/cuda_dnn.cc:9261] Unable to register cuDNN factory: Attempting to register factory for plugin cuDNN when one has already been registered
2024-02-26 16:17:41.933916: E external/local_xla/xla/stream_executor/cuda/cuda_fft.cc:607] Unable to register cuFFT factory: Attempting to register factory for plugin cuFFT when one has already been registered
2024-02-26 16:17:42.079787: E external/local_xla/xla/stream_executor/cuda/cuda_blas.cc:1515] Unable to register cuBLAS factory: Attempting to register factory for plugin cuBLAS when one has already been registered
Available GPUS:
PhysicalDevice(name='/physical_device:GPU:0', device_type='GPU')
PhysicalDevice(name='/physical_device:GPU:1', device_type='GPU')
PhysicalDevice(name='/physical_device:GPU:2', device_type='GPU')
PhysicalDevice(name='/physical_device:GPU:3', device_type='GPU')
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[2], line 14
12 import base64
13 from IPython.display import display, HTML
---> 14 import cmrseq
15 from pint import Quantity
16 import pyvista
ModuleNotFoundError: No module named 'cmrseq'
Load Mesh#
Load provided aorta model from provided vtk file
Rotate & Translate the mesh to plausible scanner coordinates
Visualize the flow field
[ ]:
# Load Aorta Example.
mesh = pyvista.read("../example_resources/pv_computed_averaged_P14.vtk")
# The velocity [m/s] is stored in the mesh as U1. In cmrsim all definitions are done in [ms] therefore
# We add the velocity field as scaled data to the mesh
mesh["velocity"] = Quantity(mesh["U1"], "m/s").m_as("m/ms")
# To match the usual geometry (Head <-> Foot along z-axis) flip the axes of the mesh
mesh.points = mesh.points[:,[2,1,0]]
mesh.points[:, 1:] *= -1
mesh["velocity"] = mesh["velocity"][:,[2,1,0]]
mesh["velocity"][:, 1:] *= -1
global_mesh_offset = [0.02, 0.01, -0.02]
mesh.translate(global_mesh_offset, inplace=True)
mesh.rotate_z(60, (0., 0., 0.), inplace=True)
mesh = mesh.cell_data_to_point_data()
# clip to restricted length to reduce mesh-size
mesh.clip(normal='z', origin=[0., 0., 0.09], inplace=True)
mesh.clip(normal='-z', origin=[0., 0., -0.09], inplace=True)
Plot 1 - Input mesh#
(To enable interactive plots change to commented plotting instructions)
[ ]:
pyvista.close_all()
# pyvista.set_jupyter_backend('panel')
pyvista.start_xvfb()
# plotter = pyvista.Plotter(notebook=True, window_size=(500, 750))
plotter = pyvista.Plotter(off_screen=True, window_size=(750, 500), notebook=False)
local_functions.add_custom_axes(plotter)
plotter.add_mesh(mesh, scalars="velocity", opacity=0.5, scalar_bar_args={"title":"Velocity [m/ms]"}, cmap="twilight")
plotter.add_mesh(mesh.slice_orthogonal(), show_scalar_bar=False, cmap="twilight")
_ = plotter.show(screenshot="input_mesh.png");
b64 = base64.b64encode(open("input_mesh.png",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))
The re-filling Dataset#
When simulating MR-signals in moving particles, it is necessary to incorporate in-/out-flow of mangetization relative to the imaging volume. Given a seeding volume, that might be defined by e.g. a slice-selective excitation pulse, the refilling dataset uniformly fills this volume with particles. When passing a set of particles (e.g. for consequtive repetitions), the dataset object drops particles outside a defined distance to the seeding volume and refills the non-populated regions to match the initial particle density.
First, a re-seeding volume is defined as slice:
[3]:
# Define Reseeding slice parameters
seeding_slice_position = Quantity([0.03, 0.005, -0.025], "m")
seeding_slice_normal = np.array([0.3, 0., 1.])
seeding_slice_thickness = Quantity(25, "mm")
lookup_pixel_spacing = np.array((0.5e-3, 0.5e-3, 0.5e-3)) # resolution of seeding mesh, mm
seeding_pixel_spacing = np.array((0.5e-3, 0.5e-3, 0.5e-3))
seeding_slice_bb = Quantity([5, 5, 5], "cm")
particle_density = 60
Plot 2 - Seeding volume#
[4]:
pyvista.close_all()
pyvista.start_xvfb()
# pyvista.set_jupyter_backend('panel')
# plotter = pyvista.Plotter(notebook=True, window_size=(500, 750))
plotter = pyvista.Plotter(off_screen=True, window_size=(750, 500), notebook=False)
# Add input mesh and axes
local_functions.add_custom_axes(plotter)
plotter.add_mesh(mesh, scalars="velocity", opacity=0.75, scalar_bar_args={"title":"Velocity [m/ms]"}, cmap="twilight")
# Plot the Re-seeding slice as orange box
slice_mesh = pyvista.Box(bounds=(-seeding_slice_bb[0].m_as("m")/2, seeding_slice_bb[0].m_as("m")/2,
-seeding_slice_bb[1].m_as("m")/2, seeding_slice_bb[1].m_as("m") /2,
-seeding_slice_thickness.m_as("m")/2, seeding_slice_thickness.m_as("m")/2))
slice_mesh.translate(seeding_slice_position.m_as("m"), inplace=True)
plotter.add_mesh(slice_mesh, color="orange", opacity=0.8, show_edges=True)
_ = plotter.show(screenshot="seeding_volume.png");
b64 = base64.b64encode(open("seeding_volume.png",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))
Instantiate Dataset and seed particles#
To be flexible with including various signal effects that can be formulated as particle properties, the refilling dataset offers a generic way to add properties on particle creation, while also being compatible with the dropping/reseeding procedure. (This is implemented by passing a dictionary of callables to on initialization as shown below). All functions in the dictionary must take the positional argument n_new and return a numpy array of length n_new.
Furthermore, the dataset offers an initial randomly uniform filling of the entire mesh:
[5]:
# Setup flow dataset with seeding mesh information
particle_creation_fncs = {"magnetization": part_factory.norm_magnetization(),
"T1": part_factory.uniform(default_value=1300.),
"T2": part_factory.uniform(default_value=300.),
"M0": part_factory.uniform(default_value=1.)}
dummy_particles = particle_creation_fncs["magnetization"](n_new=1000)
print(f"Demonstration of particle creators: Shape[{dummy_particles.shape}], DType[{dummy_particles.dtype}]")
dataset = cmrsim.datasets.RefillingFlowDataset(mesh,
particle_creation_callables=particle_creation_fncs,
slice_position=seeding_slice_position.m_as("m"),
slice_normal=seeding_slice_normal,
slice_thickness=seeding_slice_thickness.m_as("m"),
slice_bounding_box=seeding_slice_bb.m_as("m"),
lookup_map_spacing=lookup_pixel_spacing,
seeding_vol_spacing=seeding_pixel_spacing,
field_list=[("velocity", 3), ])
# with a mesh spacing of 0.5x0.5x0.5 mm^3 the lowest achievable density per mm^3 is 8
positions, properties = dataset.initial_filling(particle_density,)
print(f"Initial particles seeded: {positions.shape}, {[(k, v.shape) for k, v in properties.items()]}")
# Call reseeding to get a density estimation
%time p, pr, _ = dataset(particle_density, residual_particle_pos=positions, particle_properties=properties, distance_tolerance=0.15, reseed_threshold=0.75)
print(f"Newly seeded particles: {dataset.n_new}, Dropped particles: {dataset.n_drop}, {p.shape}")
Demonstration of particle creators: Shape[(1000, 3)], DType[complex64]
Updated Mesh
/opt/conda/lib/python3.10/site-packages/pyvista/core/filters/data_set.py:3113: PyVistaDeprecationWarning: probe filter is deprecated and will be removed in a future version.
Use sample filter instead.
If using `mesh1.probe(mesh2)`, use `mesh2.sample(mesh1)`.
warnings.warn(
Updated Slice Position
Initial particles seeded: (7538515, 3), [('magnetization', (7538515, 3)), ('T1', (7538515,)), ('T2', (7538515,)), ('M0', (7538515,))]
CPU times: user 2.49 s, sys: 557 ms, total: 3.04 s
Wall time: 3.04 s
Newly seeded particles: 14689, Dropped particles: 0, (7553204, 3)
Plot 3 - Density estimation#
This plot shows the density estimation inside the re-seeding slice directly after initially seeding particles inside the entire mesh
[16]:
pyvista.close_all()
pyvista.start_xvfb()
# pyvista.set_jupyter_backend('panel')
# plotter = pyvista.Plotter(notebook=True, window_size=(500, 750))
plotter = pyvista.Plotter(off_screen=True, window_size=(750, 500), notebook=False)
local_functions.add_custom_axes(plotter)
# Add density estimation in reseeding slice:
# (the cell_data 'density' is evaluated on calling the refilling dataset with residual particles
density_projection_mesh = dataset.project_density_to_mesh(positions)
plotter.add_mesh(density_projection_mesh.slice_orthogonal(), scalars="density",
opacity=0.9, cmap="viridis", clim=[0, particle_density],
scalar_bar_args={"title":"Density [1/mm^3]"})
plotter.add_mesh(dataset.original_mesh, opacity=0.5, color="gray")
_ = plotter.show(screenshot="density_estimation.png");
b64 = base64.b64encode(open("density_estimation.png",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))
Trajectory Module#
To simulate the motion of particles inside the velocity field provides as mesh, incrementing the positions according to the euler-forward update-rule \(r^{t+1} = r^{t} + v(r^{t}) * dt\) is used in this notebook. The triliear interpolation lookup of the velocities for each step & and particle is implemented as GPU-compatible module cmrsim.dataset.flow.FlowTrajectory:
[6]:
# Setup lookup table for field information
velocity_field_3d, map_dimensions = dataset.get_lookup_table()
# Setup trajectory module
trajectory_module = cmrsim.trajectory.FlowTrajectory(velocity_field_3d.astype(np.float32),
map_dimensions=map_dimensions.astype(np.float32),
additional_dimension_mapping=dataset.field_list[1:],
device="GPU:0"
)
# As initial positions only use new particles from re seeding. Case using all initially seeded particles see next section
demonstration_positions, _, _ = dataset(particle_density=particle_density)
# Call the module to obtain the trajectory -> (#particles, n_steps, 3)
timing = np.arange(0, 20, 0.2)
trajectories, _ = trajectory_module(demonstration_positions.astype(np.float32), timing.astype(np.float32),
dt_max=np.array(0.05, np.float32), verbose=True)
print(f"Trajectories: Shape - {trajectories.shape}")
# For debugging & visualization purposes the momentary velocity can be returned along the position
_, additional_fields = trajectory_module(demonstration_positions.astype(np.float32), timing.astype(np.float32),
dt_max=np.array(0.05, np.float32), verbose=True, return_velocities=True)
print(additional_fields[0].keys())
velocities = np.stack([_["velocity"] for _ in additional_fields], axis=0)
velocities.shape
Trajectories: Shape - (1004304, 100, 3)
dict_keys(['velocity'])
[6]:
(100, 1004304, 3)
Plot 4 - Particle trajectories#
[14]:
pyvista.close_all()
pyvista.start_xvfb()
# pyvista.set_jupyter_backend('panel')
# plotter = pyvista.Plotter(notebook=True, window_size=(500, 750), shape=(1, 2))
# plotter = pyvista.Plotter(off_screen=True, window_size=(1000, 500), notebook=False, shape=(1, 2))
plotter = pyvista.Plotter(off_screen=True, window_size=(500, 500), notebook=False)
# plotter.subplot(0, 0)
# plotter.add_text("Initial Positions", font_size=16)
# local_functions.add_custom_axes(plotter)
# plotter.add_mesh(mesh, opacity=0.5, cmap="gray", show_scalar_bar=False)
# plotter.add_mesh(pyvista.PolyData(trajectories[::100, 0].numpy()),
# scalars=np.linalg.norm(velocities[0, ::100], axis=-1),
# scalar_bar_args={"title":"Velocity [m/ms]"}, cmap="rainbow")
# plotter.subplot(0, 1)
# plotter.add_text("Position after 20ms", font_size=16)
# local_functions.add_custom_axes(plotter)
plotter.add_mesh(mesh, opacity=0.5, cmap="gray", show_scalar_bar=False)
plotter.add_mesh(pyvista.PolyData(trajectories[::50, -1].numpy()),
scalars=np.linalg.norm(velocities[-1, ::50], axis=-1),
show_scalar_bar=False,
scalar_bar_args={"title":"Velocity [m/ms]"}, cmap="rainbow",
render_points_as_spheres=True)
plotter.camera.zoom(1.2)
_ = plotter.screenshot("trajectories.png", scale=3)
# _ = plotter.show(screenshot="trajectories.png");
b64 = base64.b64encode(open("trajectories.png",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))
Refill Seeding volumes#
This section demonstrates how to use the refilling functionality after advancing particle positions using the trajectory module.
[28]:
particle_density=50
positions, properties = dataset.initial_filling(particle_density=particle_density)
positions, properties, _ = dataset(particle_density=particle_density, residual_particle_pos=positions, particle_properties=properties, reseed_threshold=1., distance_tolerance=1.)
# positions, properties, _ = dataset(particle_density=particle_density, reseed_threshold=1, distance_tolerance=1.)
print(f"Density in seeding volume before re-seeding: {dataset.mean_density: 1.5f} $1/mm^3$")
print(f"Initial particles seeded: {positions.shape}, {[(k, v.shape) for k, v in properties.items()]}")
r_new = positions
for step in tqdm(range(1000)):
r_new, _ = trajectory_module.increment_particles(r_new, dt=0.01)
positions_new, properties_new, _ = dataset(particle_density=particle_density, residual_particle_pos=r_new.numpy(),
particle_properties=properties, reseed_threshold=0.7, distance_tolerance=0.5)
density_before = np.array(dataset.gridded_seeding_volume["density"])
print(f"New particles created: {dataset.n_new}\nParticles destroyed: {dataset.n_drop} "
f"\nNew array shape: {positions_new.shape}\nDensity in seeding volume before re-seeding: {dataset.mean_density: 1.5f} $1/mm^3$")
_ = dataset(particle_density=particle_density, residual_particle_pos=positions_new,
particle_properties=properties_new, reseed_threshold=0.7, distance_tolerance=0.5)
print(f"Density in seeding volume after re-seeding: {dataset.mean_density: 1.5f} $1/mm^3$")
Density in seeding volume before re-seeding: 0.10675 $1/mm^3$
Initial particles seeded: (6311345, 3), [('magnetization', (6311345, 3)), ('T1', (6311345,)), ('T2', (6311345,)), ('M0', (6311345,))]
New particles created: 247572
Particles destroyed: 187602
New array shape: (6371315, 3)
Density in seeding volume before re-seeding: 0.06731 $1/mm^3$
Density in seeding volume after re-seeding: 0.13860 $1/mm^3$
Plot 5 - Refilled particles#
Illustration of particle density before and after re-filling inside the specified re-seeding box. Note how the particle depleted regions are homogenously filled without increasing the density in non-depleted regions.
[29]:
from copy import deepcopy
d2 = deepcopy(dataset.gridded_seeding_volume)
d2["density"] = density_before
pyvista.close_all()
pyvista.start_xvfb()
# pyvista.set_jupyter_backend('panel')
# plotter = pyvista.Plotter(notebook=True, window_size=(500, 750), shape=(1, 2))
plotter = pyvista.Plotter(off_screen=True, window_size=(1000, 500), notebook=False, shape=(1, 2))
plotter.subplot(0, 0)
plotter.add_text("Density before re-seeding", font_size=16)
local_functions.add_custom_axes(plotter)
plotter.add_mesh(d2.slice_orthogonal(), scalars="density", cmap="jet")
plotter.subplot(0, 1)
plotter.add_text("Density after re-seeding", font_size=16)
local_functions.add_custom_axes(plotter)
plotter.add_mesh(dataset.gridded_seeding_volume.slice_orthogonal(), scalars="density", cmap="jet")
_ = plotter.show(screenshot="refilled_particles.png");
b64 = base64.b64encode(open("refilled_particles.png",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))
Additional Fields Lookup#
Spatial varying fields like off-resonance or turbulence parameters that influence the signal formation for moving particles can be looked up for every particle location. To increase eficiency, the lookup is implemented as additional channels of the velocity fields. The following cells demonstrates how the additional fields can be incorporated
[37]:
# This field is not realistic but it shows how the mechanism works
mesh["off_res"] = np.ones([mesh["velocity"].shape[0], 1])
# Setup flow dataset with seeding mesh information
particle_creation_fncs = {"magnetization": part_factory.norm_magnetization(),
"T1": part_factory.uniform(default_value=1300.),
"T2": part_factory.uniform(default_value=300.),
"M0": part_factory.uniform(default_value=1.)}
# Instantiate the dataset and specify the fields to be used in the field_list argument
dataset = cmrsim.datasets.RefillingFlowDataset(mesh,
particle_creation_callables=particle_creation_fncs,
slice_position=seeding_slice_position.m_as("m"),
slice_normal=seeding_slice_normal,
slice_thickness=seeding_slice_thickness.m_as("m"),
lookup_map_spacing=seeding_pixel_spacing,
seeding_vol_spacing=seeding_pixel_spacing,
field_list=[("velocity", 3), ("off_res", 1)])
# # Setup lookup table for field information
velocity_field_3d, map_dimensions = dataset.get_lookup_table()
# Setup trajectory module
trajectory_module = cmrsim.trajectory.FlowTrajectory(velocity_field_3d,
map_dimensions=map_dimensions.astype(np.float32),
additional_dimension_mapping=dataset.field_list[1:],
device="GPU:0")
positions, properties, _ = dataset(particle_density=1)
timing = Quantity(np.arange(0, 10, 0.1), "ms")
trajectories, additional_fields = trajectory_module(positions.astype(np.float32), timing=timing.m_as("ms").astype(np.float32),
dt_max=Quantity(0.05, "ms").m, return_velocities=True)
print(f"\nShape of trajectories: {trajectories.shape} \n"
f"Names and shapes of additional fields per step: {[(k, v.shape) for k, v in additional_fields[-1].items()]}\n")
Updated Mesh
Updated Slice Position
Shape of trajectories: (26529, 100, 3)
Names and shapes of additional fields per step: [('off_res', TensorShape([26529, 1])), ('velocity', TensorShape([26529, 3]))]
Sliced initial filling#
Special use-case: Only initially fill a slab withing the mesh
[ ]:
slice_position = Quantity([0.03, 0.00, -0.025], "m")
slice_normal = np.array([0., 1., -0.1])
slice_thickness = Quantity(10, "mm")
seeding_slice_bb = Quantity([6, 3, 3], "cm")
lookup_pixel_spacing = np.array((0.5e-3, 0.5e-3, 0.5e-3)) # resolution of seeding mesh, mm
seeding_pixel_spacing = np.array((0.5e-3, 0.5e-3, 0.5e-3))
particle_density = 128
particle_creation_fncs = {"magnetization": part_factory.norm_magnetization(),
"T1": part_factory.uniform(default_value=1300.),
"T2": part_factory.uniform(default_value=300.),
"M0": part_factory.uniform(default_value=1.)}
dataset = cmrsim.datasets.RefillingFlowDataset(mesh,
particle_creation_callables=particle_creation_fncs,
slice_position=slice_position.m_as("m"),
slice_normal=slice_normal,
slice_thickness=slice_thickness.m_as("m"),
slice_bounding_box=seeding_slice_bb.m_as("m"),
lookup_map_spacing=lookup_pixel_spacing,
seeding_vol_spacing=seeding_pixel_spacing,
field_list=[("velocity", 3), ])
positions, properties = dataset.initial_filling(particle_density,
slice_dictionary=dict(slice_normal=slice_normal,
slice_position=slice_position.m_as("m"),
slice_thickness=slice_thickness.m_as("m")))
# Setup lookup table for field information
velocity_field_3d, map_dimensions = dataset.get_lookup_table()
# Setup trajectory module
trajectory_module = cmrsim.trajectory.FlowTrajectory(velocity_field_3d,
map_dimensions=map_dimensions.astype(np.float32),
additional_dimension_mapping=dataset.field_list[1:],
device="GPU:0")
timing = Quantity(np.arange(0, 10, 0.1), "ms")
trajectories, additional_fields = trajectory_module(positions.astype(np.float32), timing=timing.m_as("ms").astype(np.float32),
dt_max=Quantity(0.05, "ms").m, return_velocities=True)
# positions_new, properties_new, _ = dataset(particle_density=particle_density, residual_particle_pos=r_new.numpy(),
# particle_properties=properties, reseed_threshold=0.9, distance_tolerance=0.5)
# print(dataset.n_new, dataset.mean_density)
[46]:
pyvista.close_all()
pyvista.start_xvfb()
# pyvista.set_jupyter_backend('panel')
# plotter = pyvista.Plotter(notebook=True, window_size=(500, 750), shape=(1, 2))
plotter = pyvista.Plotter(off_screen=True, window_size=(1000, 500), notebook=False, shape=(1, 2))
plotter.subplot(0, 0)
local_functions.add_custom_axes(plotter)
inital_seeding_slice = local_functions.transformed_box(reference_basis=np.eye(3, 3), slice_normal=slice_normal,
readout_direction=np.array([1., 0., 0.]), slice_position=slice_position,
extend=Quantity([10, 20, slice_thickness.m_as("cm")], "cm"))
plotter.add_mesh(inital_seeding_slice, opacity=0.5, show_edges=True, color="red")
plotter.add_mesh(mesh, opacity=0.75, color="gray")
plotter.add_mesh(dataset.gridded_seeding_volume, color="orange", show_edges=False, opacity=0.8)
plotter.subplot(0, 1)
local_functions.add_custom_axes(plotter)
plotter.add_mesh(mesh, opacity=0.75, color="gray")
plotter.add_mesh(pyvista.PolyData(positions), color="white")
_ = plotter.show(screenshot="sliced_inital_filling.png");
b64 = base64.b64encode(open("sliced_inital_filling.png",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))
/scratch/jweine/cmrsim/notebooks/datasets/../local_functions.py:43: UserWarning: Optimal rotation is not uniquely or poorly defined for the given sets of vectors.
rot, _ = Rotation.align_vectors(reference_basis, new_basis)
