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#
[1]:
# Load TF and check for GPUs
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)
# 3rd Party dependencies
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.append("../../")
import cmrsim
import cmrsim.utils.particle_properties as part_factory
2022-12-19 19:27:27.817151: I tensorflow/core/platform/cpu_feature_guard.cc:193] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations: AVX2 FMA
To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.
2022-12-19 19:27:27.937542: W tensorflow/compiler/xla/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory
2022-12-19 19:27:27.937569: I tensorflow/compiler/xla/stream_executor/cuda/cudart_stub.cc:29] Ignore above cudart dlerror if you do not have a GPU set up on your machine.
2022-12-19 19:27:28.869890: W tensorflow/compiler/xla/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libnvinfer.so.7'; dlerror: libnvinfer.so.7: cannot open shared object file: No such file or directory
2022-12-19 19:27:28.870018: W tensorflow/compiler/xla/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libnvinfer_plugin.so.7'; dlerror: libnvinfer_plugin.so.7: cannot open shared object file: No such file or directory
2022-12-19 19:27:28.870031: W tensorflow/compiler/tf2tensorrt/utils/py_utils.cc:38] TF-TRT Warning: Cannot dlopen some TensorRT libraries. If you would like to use Nvidia GPU with TensorRT, please make sure the missing libraries mentioned above are installed properly.
2022-12-19 19:27:29.796475: W tensorflow/compiler/xla/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcuda.so.1'; dlerror: libcuda.so.1: cannot open shared object file: No such file or directory
2022-12-19 19:27:29.796535: W tensorflow/compiler/xla/stream_executor/cuda/cuda_driver.cc:265] failed call to cuInit: UNKNOWN ERROR (303)
2022-12-19 19:27:29.796577: I tensorflow/compiler/xla/stream_executor/cuda/cuda_diagnostics.cc:156] kernel driver does not appear to be running on this host (runner-xwyxd7yv-project-23798-concurrent-1): /proc/driver/nvidia/version does not exist
2022-12-19 19:27:31.571818: I tensorflow/core/platform/cpu_feature_guard.cc:193] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations: AVX2 FMA
To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.
Load Mesh#
Load provided aorta model from provided vtk file
Rotate & Translate the mesh to plausible scanner coordinates
Visualize the flow field
[2]:
# 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)
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In [2], line 2
1 # Load Aorta Example.
----> 2 mesh = pyvista.read("../example_resources/pv_computed_averaged_P14.vtk")
4 # The velocity [m/s] is stored in the mesh as U1. In cmrsim all definitions are done in [ms] therefore
5 # We add the velocity field as scaled data to the mesh
6 mesh["velocity"] = Quantity(mesh["U1"], "m/s").m_as("m/ms")
File /opt/conda/lib/python3.9/site-packages/pyvista/utilities/fileio.py:166, in read(filename, attrs, force_ext, file_format, progress_bar)
164 filename = os.path.abspath(os.path.expanduser(str(filename)))
165 if not os.path.isfile(filename):
--> 166 raise FileNotFoundError(f'File ({filename}) not found')
168 # Read file using meshio.read if file_format is present
169 if file_format:
FileNotFoundError: File (/builds/jweine/cmrsim/docs/source/example_gallery/example_resources/pv_computed_averaged_P14.vtk) not found
Plot 1 - Input mesh#
[3]:
pyvista.close_all()
pyvista.set_jupyter_backend('panel')
pyvista.start_xvfb()
plotter = pyvista.Plotter(notebook=True, window_size=(500, 750))
# Plot axes at origin
origin_arrows = [pyvista.Arrow(start=(0., 0., 0.), direction=direction, tip_length=0.1, tip_radius=0.025,
tip_resolution=10, shaft_radius=0.005, shaft_resolution=10, scale=0.1)
for direction in [(1., 0., 0.), (0., 1., 0.), (0., 0., 1.)]]
plotter.add_mesh(origin_arrows[0], color="blue", label="x")
plotter.add_mesh(origin_arrows[1], color="green", label="y")
plotter.add_mesh(origin_arrows[2], color="red", label="z")
# Add mesh
plotter.add_mesh(mesh, scalars="velocity", opacity=0.5, scalar_bar_args={"title":"Velocity [m/ms]"})
plotter.add_mesh(mesh.slice_orthogonal())
plotter.show()
---------------------------------------------------------------------------
OSError Traceback (most recent call last)
Cell In [3], line 3
1 pyvista.close_all()
2 pyvista.set_jupyter_backend('panel')
----> 3 pyvista.start_xvfb()
4 plotter = pyvista.Plotter(notebook=True, window_size=(500, 750))
6 # Plot axes at origin
File /opt/conda/lib/python3.9/site-packages/pyvista/utilities/xvfb.py:47, in start_xvfb(wait, window_size)
44 raise OSError('`start_xvfb` is only supported on Linux')
46 if os.system('which Xvfb > /dev/null'):
---> 47 raise OSError(XVFB_INSTALL_NOTES)
49 # use current default window size
50 if window_size is None:
OSError: Please install Xvfb with:
Debian
$ sudo apt install libgl1-mesa-glx xvfb
CentOS / RHL
$ sudo yum install libgl1-mesa-glx xvfb
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:
[4]:
# 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#
[5]:
pyvista.close_all()
pyvista.set_jupyter_backend('panel')
pyvista.start_xvfb()
plotter = pyvista.Plotter(notebook=True, window_size=(500, 750))
# Plot axes at origin
origin_arrows = [pyvista.Arrow(start=(0., 0., 0.), direction=direction, tip_length=0.1, tip_radius=0.025,
tip_resolution=10, shaft_radius=0.005, shaft_resolution=10, scale=0.1)
for direction in [(1., 0., 0.), (0., 1., 0.), (0., 0., 1.)]]
plotter.add_mesh(origin_arrows[0], color="blue", label="x")
plotter.add_mesh(origin_arrows[1], color="green", label="y")
plotter.add_mesh(origin_arrows[2], color="red", label="z")
# Add mesh
plotter.add_mesh(mesh, scalars="velocity", opacity=0.75, scalar_bar_args={"title":"Velocity [m/ms]"})
# 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()
---------------------------------------------------------------------------
OSError Traceback (most recent call last)
Cell In [5], line 3
1 pyvista.close_all()
2 pyvista.set_jupyter_backend('panel')
----> 3 pyvista.start_xvfb()
4 plotter = pyvista.Plotter(notebook=True, window_size=(500, 750))
6 # Plot axes at origin
File /opt/conda/lib/python3.9/site-packages/pyvista/utilities/xvfb.py:47, in start_xvfb(wait, window_size)
44 raise OSError('`start_xvfb` is only supported on Linux')
46 if os.system('which Xvfb > /dev/null'):
---> 47 raise OSError(XVFB_INSTALL_NOTES)
49 # use current default window size
50 if window_size is None:
OSError: Please install Xvfb with:
Debian
$ sudo apt install libgl1-mesa-glx xvfb
CentOS / RHL
$ sudo yum install libgl1-mesa-glx xvfb
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:
[6]:
# 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]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [6], line 10
7 dummy_particles = particle_creation_fncs["magnetization"](n_new=1000)
8 print(f"Demonstration of particle creators: Shape[{dummy_particles.shape}], DType[{dummy_particles.dtype}]")
---> 10 dataset = cmrsim.datasets.RefillingFlowDataset(mesh,
11 particle_creation_callables=particle_creation_fncs,
12 slice_position=seeding_slice_position.m_as("m"),
13 slice_normal=seeding_slice_normal,
14 slice_thickness=seeding_slice_thickness.m_as("m"),
15 slice_bounding_box=seeding_slice_bb.m_as("m"),
16 lookup_map_spacing=lookup_pixel_spacing,
17 seeding_vol_spacing=seeding_pixel_spacing,
18 field_list=[("velocity", 3), ])
20 # with a mesh spacing of 0.5x0.5x0.5 mm^3 the lowest achievable density per mm^3 is 8
21 positions, properties = dataset.initial_filling(particle_density,)
NameError: name 'mesh' is not defined
Plot 3 - Density estimation#
This plot shows the density estimation inside the re-seeding slice directly after initially seeding particles inside the entire mesh
[7]:
pyvista.close_all()
pyvista.set_jupyter_backend('panel')
pyvista.start_xvfb()
plotter = pyvista.Plotter(notebook=True, window_size=(500, 750))
# Plot axes at origin
origin_arrows = [pyvista.Arrow(start=(0., 0., 0.), direction=direction, tip_length=0.1, tip_radius=0.025,
tip_resolution=10, shaft_radius=0.005, shaft_resolution=10, scale=0.1)
for direction in [(1., 0., 0.), (0., 1., 0.), (0., 0., 1.)]]
plotter.add_mesh(origin_arrows[0], color="blue", label="x")
plotter.add_mesh(origin_arrows[1], color="green", label="y")
plotter.add_mesh(origin_arrows[2], color="red", label="z")
# Add seeded particles
plotter.add_mesh(pyvista.PolyData(positions[::100]), color="white", opacity=0.1)
# 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, scalars="density",
opacity=0.5, cmap="viridis", clim=[0, particle_density],
scalar_bar_args={"title":"Density [1/mm^3]"})
plotter.show()
---------------------------------------------------------------------------
OSError Traceback (most recent call last)
Cell In [7], line 3
1 pyvista.close_all()
2 pyvista.set_jupyter_backend('panel')
----> 3 pyvista.start_xvfb()
4 plotter = pyvista.Plotter(notebook=True, window_size=(500, 750))
5 # Plot axes at origin
File /opt/conda/lib/python3.9/site-packages/pyvista/utilities/xvfb.py:47, in start_xvfb(wait, window_size)
44 raise OSError('`start_xvfb` is only supported on Linux')
46 if os.system('which Xvfb > /dev/null'):
---> 47 raise OSError(XVFB_INSTALL_NOTES)
49 # use current default window size
50 if window_size is None:
OSError: Please install Xvfb with:
Debian
$ sudo apt install libgl1-mesa-glx xvfb
CentOS / RHL
$ sudo yum install libgl1-mesa-glx xvfb
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:
[8]:
# Setup lookup table for field information
velocity_field_3d, map_dimensions = dataset.get_lookup_table()
# Setup trajectory module
trajectory_module = cmrsim.datasets.FlowTrajectory(velocity_field_3d,
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)
trajectories, _ = trajectory_module(initial_positions=demonstration_positions, dt=0.1, n_steps=200, 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(initial_positions=demonstration_positions, dt=0.1, n_steps=200, verbose=True, return_velocities=True)
print(additional_fields[0].keys())
velocities = np.stack([_["velocity"] for _ in additional_fields], axis=0)
velocities.shape
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [8], line 2
1 # Setup lookup table for field information
----> 2 velocity_field_3d, map_dimensions = dataset.get_lookup_table()
4 # Setup trajectory module
5 trajectory_module = cmrsim.datasets.FlowTrajectory(velocity_field_3d,
6 map_dimensions=map_dimensions.astype(np.float32),
7 additional_dimension_mapping=dataset.field_list[1:],
8 device="GPU:0")
NameError: name 'dataset' is not defined
Plot 4 - Particle trajectories#
[9]:
pyvista.close_all()
pyvista.set_jupyter_backend('panel')
pyvista.start_xvfb()
plotter = pyvista.Plotter(notebook=True, window_size=(500, 750), shape=(1, 2))
plotter.subplot(0, 0)
plotter.add_text("Initial Positions", font_size=16)
plotter.add_mesh(mesh, opacity=0.5, cmap="gray")
plotter.add_mesh(pyvista.PolyData(trajectories[::10, 0].numpy()),
scalars=np.linalg.norm(velocities[0, ::10], axis=-1),
scalar_bar_args={"title":"Velocity [m/ms]"}, cmap="jet")
plotter.subplot(0, 1)
plotter.add_text("Position after 20ms", font_size=16)
plotter.add_mesh(mesh, opacity=0.5, cmap="gray")
plotter.add_mesh(pyvista.PolyData(trajectories[::10, -1].numpy()),
scalars=np.linalg.norm(velocities[-1, ::10], axis=-1),
scalar_bar_args={"title":"Velocity [m/ms]"}, cmap="jet")
plotter.show()
---------------------------------------------------------------------------
OSError Traceback (most recent call last)
Cell In [9], line 3
1 pyvista.close_all()
2 pyvista.set_jupyter_backend('panel')
----> 3 pyvista.start_xvfb()
4 plotter = pyvista.Plotter(notebook=True, window_size=(500, 750), shape=(1, 2))
6 plotter.subplot(0, 0)
File /opt/conda/lib/python3.9/site-packages/pyvista/utilities/xvfb.py:47, in start_xvfb(wait, window_size)
44 raise OSError('`start_xvfb` is only supported on Linux')
46 if os.system('which Xvfb > /dev/null'):
---> 47 raise OSError(XVFB_INSTALL_NOTES)
49 # use current default window size
50 if window_size is None:
OSError: Please install Xvfb with:
Debian
$ sudo apt install libgl1-mesa-glx xvfb
CentOS / RHL
$ sudo yum install libgl1-mesa-glx xvfb
Refill Seeding volumes#
This section demonstrates how to use the refilling and trajectory module without saving the entire trajectory. In this demonstration the initial filling is used as input to the trajectory module.
[10]:
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$")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [10], line 2
1 particle_density=50
----> 2 positions, properties = dataset.initial_filling(particle_density=particle_density)
3 positions, properties, _ = dataset(particle_density=particle_density, residual_particle_pos=positions, particle_properties=properties, reseed_threshold=1., distance_tolerance=1.)
4 # positions, properties, _ = dataset(particle_density=particle_density, reseed_threshold=1, distance_tolerance=1.)
NameError: name 'dataset' is not defined
Plot 5 - Refilled particles#
[11]:
pyvista.close_all()
pyvista.set_jupyter_backend('panel')
pyvista.start_xvfb()
plotter = pyvista.Plotter(notebook=True, window_size=(500, 750), shape=(1, 2))
# Plot axes at origin
for i in range(2):
plotter.subplot(0, i)
origin_arrows = [pyvista.Arrow(start=(0., 0., 0.), direction=direction, tip_length=0.1, tip_radius=0.025,
tip_resolution=10, shaft_radius=0.005, shaft_resolution=10, scale=0.1)
for direction in [(1., 0., 0.), (0., 1., 0.), (0., 0., 1.)]]
plotter.add_mesh(origin_arrows[0], color="blue", label="x")
plotter.add_mesh(origin_arrows[1], color="green", label="y")
plotter.add_mesh(origin_arrows[2], color="red", label="z")
# Add the original mesh to the plot
plotter.subplot(0, 0)
from copy import deepcopy
d2 = deepcopy(dataset.gridded_seeding_volume)
d2["density"] = density_before
plotter.add_text("Density before re-seeding", font_size=16)
plotter.add_mesh(d2.slice_orthogonal(), scalars="density", cmap="jet")
plotter.subplot(0, 1)
plotter.add_text("Density after re-seeding", font_size=16)
plotter.add_mesh(dataset.gridded_seeding_volume.slice_orthogonal(), scalars="density", cmap="jet")
plotter.show()
---------------------------------------------------------------------------
OSError Traceback (most recent call last)
Cell In [11], line 3
1 pyvista.close_all()
2 pyvista.set_jupyter_backend('panel')
----> 3 pyvista.start_xvfb()
4 plotter = pyvista.Plotter(notebook=True, window_size=(500, 750), shape=(1, 2))
6 # Plot axes at origin
File /opt/conda/lib/python3.9/site-packages/pyvista/utilities/xvfb.py:47, in start_xvfb(wait, window_size)
44 raise OSError('`start_xvfb` is only supported on Linux')
46 if os.system('which Xvfb > /dev/null'):
---> 47 raise OSError(XVFB_INSTALL_NOTES)
49 # use current default window size
50 if window_size is None:
OSError: Please install Xvfb with:
Debian
$ sudo apt install libgl1-mesa-glx xvfb
CentOS / RHL
$ sudo yum install libgl1-mesa-glx xvfb
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
[12]:
# 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.datasets.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)
trajectories, additional_fields = trajectory_module(positions, dt=0.1, n_steps=100)
print(f"\nShape of trajectories: {trajectories.shape} \n"
f"Names and shapes of additional fields: {[(k, v.shape) for k, v in additional_fields[-1].items()]}\n")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [12], line 2
1 # This field is not realistic but it shows how the mechanism works
----> 2 mesh["off_res"] = np.ones([mesh["velocity"].shape[0], 1])
4 # Setup flow dataset with seeding mesh information
5 particle_creation_fncs = {"magnetization": part_factory.norm_magnetization(),
6 "T1": part_factory.uniform(default_value=1300.),
7 "T2": part_factory.uniform(default_value=300.),
8 "M0": part_factory.uniform(default_value=1.)}
NameError: name 'mesh' is not defined
Sliced initial filling#
[13]:
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.datasets.FlowTrajectory(velocity_field_3d,
map_dimensions=map_dimensions.astype(np.float32),
additional_dimension_mapping=dataset.field_list[1:],
device="GPU:0")
r_new = positions
for step in tqdm(range(500)):
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.9, distance_tolerance=0.5)
print(dataset.n_new, dataset.mean_density)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [13], line 15
7 particle_density = 128
9 particle_creation_fncs = {"magnetization": part_factory.norm_magnetization(),
10 "T1": part_factory.uniform(default_value=1300.),
11 "T2": part_factory.uniform(default_value=300.),
12 "M0": part_factory.uniform(default_value=1.)}
---> 15 dataset = cmrsim.datasets.RefillingFlowDataset(mesh,
16 particle_creation_callables=particle_creation_fncs,
17 slice_position=slice_position.m_as("m"),
18 slice_normal=slice_normal,
19 slice_thickness=slice_thickness.m_as("m"),
20 slice_bounding_box=seeding_slice_bb.m_as("m"),
21 lookup_map_spacing=lookup_pixel_spacing,
22 seeding_vol_spacing=seeding_pixel_spacing,
23 field_list=[("velocity", 3), ])
25 positions, properties = dataset.initial_filling(particle_density,
26 slice_dictionary=dict(slice_normal=slice_normal,
27 slice_position=slice_position.m_as("m"),
28 slice_thickness=slice_thickness.m_as("m")))
30 # Setup lookup table for field information
NameError: name 'mesh' is not defined
[14]:
pyvista.close_all()
pyvista.set_jupyter_backend('panel')
pyvista.start_xvfb()
plotter = pyvista.Plotter(notebook=True, window_size=(750, 1100), shape=(1, 3))
# Plot axes at origin
origin_arrows = [pyvista.Arrow(start=(0., 0., 0.), direction=direction, tip_length=0.1, tip_radius=0.025,
tip_resolution=10, shaft_radius=0.005, shaft_resolution=10, scale=0.1)
for direction in [(1., 0., 0.), (0., 1., 0.), (0., 0., 1.)]]
plotter.add_mesh(origin_arrows[0], color="blue", label="x")
plotter.add_mesh(origin_arrows[1], color="green", label="y")
plotter.add_mesh(origin_arrows[2], color="red", label="z")
_scale = 0.1
inital_seeding_slice = pyvista.Box(bounds=(-_scale, _scale, -_scale, _scale, -slice_thickness.m_as("m")/2, slice_thickness.m_as("m")/2))
inital_seeding_slice.points = cmrseq.utils.mps_to_xyz(inital_seeding_slice.points[np.newaxis], slice_normal=slice_normal, readout_direction=np.array([1., 0., 0.]))[0]
inital_seeding_slice.translate(slice_position, inplace=True)
plotter.add_mesh(inital_seeding_slice, opacity=0.1, show_edges=True, color="red")
plotter.add_mesh(mesh, opacity=0.75)
plotter.add_mesh(pyvista.PolyData(positions), color="white")
plotter.add_mesh(dataset.gridded_seeding_volume, color="orange", show_edges=False, opacity=0.8)
plotter.subplot(0, 1)
# Plot axes at origin
origin_arrows = [pyvista.Arrow(start=(0., 0., 0.), direction=direction, tip_length=0.1, tip_radius=0.025,
tip_resolution=10, shaft_radius=0.005, shaft_resolution=10, scale=0.1)
for direction in [(1., 0., 0.), (0., 1., 0.), (0., 0., 1.)]]
plotter.add_mesh(origin_arrows[0], color="blue", label="x")
plotter.add_mesh(origin_arrows[1], color="green", label="y")
plotter.add_mesh(origin_arrows[2], color="red", label="z")
# plotter.add_mesh(pyvista.PolyData(positions_new[:-dataset.n_new:10]), color="white")
plotter.add_mesh(pyvista.PolyData(positions_new[-dataset.n_new:]), color="green")
plotter.subplot(0, 2)
plotter.add_mesh(dataset.gridded_seeding_volume, scalars="density")
plotter.show()
---------------------------------------------------------------------------
OSError Traceback (most recent call last)
Cell In [14], line 3
1 pyvista.close_all()
2 pyvista.set_jupyter_backend('panel')
----> 3 pyvista.start_xvfb()
4 plotter = pyvista.Plotter(notebook=True, window_size=(750, 1100), shape=(1, 3))
6 # Plot axes at origin
File /opt/conda/lib/python3.9/site-packages/pyvista/utilities/xvfb.py:47, in start_xvfb(wait, window_size)
44 raise OSError('`start_xvfb` is only supported on Linux')
46 if os.system('which Xvfb > /dev/null'):
---> 47 raise OSError(XVFB_INSTALL_NOTES)
49 # use current default window size
50 if window_size is None:
OSError: Please install Xvfb with:
Debian
$ sudo apt install libgl1-mesa-glx xvfb
CentOS / RHL
$ sudo yum install libgl1-mesa-glx xvfb
[ ]: