Data on Regular Grids#
Imports#
[1]:
import sys
sys.path.insert(0, "../..")
sys.path.append("../")
import base64
from IPython.display import display, HTML
import pyvista
from pint import Quantity
import numpy as np
import imageio
import local_functions
import cmrsim
2022-12-09 13:35:32.046639: E tensorflow/stream_executor/cuda/cuda_driver.cc:271] failed call to cuInit: UNKNOWN ERROR (34)
2022-12-09 13:35:32.046695: I tensorflow/stream_executor/cuda/cuda_diagnostics.cc:156] kernel driver does not appear to be running on this host (89d9d482af0a): /proc/driver/nvidia/version does not exist
2022-12-09 13:35:32.046918: 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.
Offresonance in static Spherical Phantom#
In the following section, a spherical water phantom with air-filled holes is created. Subsequently, the RegularGridDataset is used to voxelize the unstructured grid and calculate the offresonsance in a static magentic field due to interfaces with different susceptibility. Finally, the phantom is returned as set of points as used for simulation.
Create and save phantom#
[ ]:
raster = pyvista.UniformGrid(dimensions=(120, 120, 120), spacing=(0.002, 0.002, 0.002), origin=(-0.11, -0.11, -0.11))
sphere_big = pyvista.Sphere(radius=0.1).extract_surface()
box = pyvista.Box(bounds=(-0.03, 0.03, -0.02, 0.02, -0.01, 0.01))
o = 0.07
spheres_small = [pyvista.Sphere(radius=0.005, center=c).extract_surface() for c in
((-o, 0., 0.), (o, 0., 0.), (0., -o, 0.), (0., o, 0.), (0., 0., -o), (0., 0., o))]
del o
filled_sphere = raster.clip_surface(sphere_big, progress_bar=True)
filled_sphere = filled_sphere.clip_surface(box, progress_bar=True, invert=False)
for sms in spheres_small:
filled_sphere = filled_sphere.clip_surface(sms, progress_bar=True, invert=False)
filled_sphere.save("spherical_phantom.vtk")
Assign physical properties#
[2]:
# load phantom
res = pyvista.read("spherical_phantom.vtk")
# Create field to signal if regular mesh-points are within original unstructured grid
res["in_mesh"] = np.ones(res.points.shape[0])
dataset = cmrsim.datasets.RegularGridDataset.from_unstructured_grid(res, pixel_spacing=Quantity([0.5, 0.5, 0.5], "mm")*3)
# Assign properties (everything not 'in-mesh' is assumed to be air)
dataset.mesh["chi"] = np.where(dataset.mesh["in_mesh"], np.ones_like(dataset.mesh["in_mesh"]) * (-9.05), np.ones_like(dataset.mesh["in_mesh"]) * 0.36) * 1e-6 # susceptibilty in ppm
dataset.mesh["M0"] = np.where(dataset.mesh["in_mesh"], np.ones_like(dataset.mesh["in_mesh"]), np.zeros_like(dataset.mesh["in_mesh"])) # density in percent
dataset.mesh["T1"] = np.where(dataset.mesh["in_mesh"], np.ones_like(dataset.mesh["in_mesh"]) * 1000, np.zeros_like(dataset.mesh["in_mesh"])) # time in ms
dataset.mesh["T2"] = np.where(dataset.mesh["in_mesh"], np.ones_like(dataset.mesh["in_mesh"]) * 300, np.zeros_like(dataset.mesh["in_mesh"])) # time in ms
dataset.mesh
[2]:
| Header | Data Arrays | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
Compute off-resonance map#
(The result is stored as field in the dataset.mesh with name ‘offres’)
[3]:
dataset.compute_offresonance(b0=Quantity(3, "T"), susceptibility_key="chi");
/code/cmrsim/notebooks/datasets/../../cmrsim/datasets/_regular_grid.py:204: RuntimeWarning: invalid value encountered in true_divide
kernel = 1/3 - kz2 / (kx2 + ky2 + kz2)
/usr/local/lib/python3.8/dist-packages/pyvista/core/dataset.py:1969: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
scalars = np.array(scalars)
Plot result#
[20]:
pyvista.close_all()
pyvista.start_xvfb()
plotter = pyvista.Plotter(notebook=False, off_screen=True, window_size=(1200, 600), shape=(1, 2))
plotter.add_mesh(res, opacity=0.6)
local_functions.add_custom_axes(plotter)
plotter.subplot(0, 1)
plotter.add_mesh(dataset.mesh.slice_orthogonal(), scalars="offres", opacity=0.8, cmap="twilight", scalar_bar_args=dict(title="Offresonsance (T)"))
local_functions.add_custom_axes(plotter)
plotter.screenshot("Offres.png")
b64 = base64.b64encode(open("Offres.png",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))
Resample to higher resolution#
[9]:
dataset.resample_spacing(Quantity([0.5, 0.5, 0.5], "mm"))
dataset.mesh
[9]:
| Header | Data Arrays | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
Select Slice#
MPS-coordinates
[7]:
slice_normal = np.array([0., 1., 1.])
readout_direction = np.array([1., 1., 0.])
slice_position = Quantity([0., 0., 1], "cm")
field_of_view = Quantity([30, 30, 1], "cm")
spacing = Quantity([1, 1, 1], "mm")
global_slice = dataset.select_slice(slice_normal, spacing, slice_position, field_of_view, readout_direction, in_mps=False)
mps_slice = dataset.select_slice(slice_normal, spacing, slice_position, field_of_view, readout_direction, in_mps=True)
/code/cmrsim/notebooks/datasets/../../cmrsim/datasets/_regular_grid.py:114: UserWarning: Optimal rotation is not uniquely or poorly defined for the given sets of vectors.
rot, _ = Rotation.align_vectors(original_basis, new_basis)
Illustrate slice selection as plots
[8]:
pyvista.close_all()
pyvista.start_xvfb()
plotter = pyvista.Plotter(notebook=False, off_screen=True, window_size=(300, 300))
plotter.add_mesh(dataset.mesh.slice_orthogonal(), scalars="offres",
opacity=0.8, cmap="twilight", scalar_bar_args=dict(title="Offresonsance (T)"))
local_functions.add_custom_axes(plotter)
plotter.add_title("Example input dataset", font_size=6)
plotter.screenshot("select_slice0.png");
plotter = pyvista.Plotter(notebook=False, off_screen=True, window_size=(300, 300))
plotter.add_mesh(dataset.mesh.slice_orthogonal(), scalars="offres", opacity=0.8,
cmap="twilight", scalar_bar_args=dict(title="Offresonsance (T)"))
plotter.add_mesh(global_slice, cmap="twilight", show_scalar_bar=False, opacity=0.8)
local_functions.add_custom_axes(plotter)
plotter.add_title("Probing slice in gobal coordinates", font_size=6)
plotter.screenshot("select_slice1.png");
plotter = pyvista.Plotter(notebook=False, off_screen=True, window_size=(300, 300))
plotter.add_mesh(mps_slice, scalars="offres", cmap="twilight", scalar_bar_args=dict(title="Offresonsance (T)"))
local_functions.add_custom_axes(plotter)
plotter.add_title("Slice back-transformed to MPS-coordinates", font_size=6)
plotter.screenshot("select_slice2.png");
plotter = pyvista.Plotter(notebook=False, off_screen=True, window_size=(300, 300))
plotter.add_mesh(mps_slice.threshold(0.01, scalars="M0"), scalars="offres",
cmap="twilight", scalar_bar_args=dict(title="Offresonsance (T)"))
local_functions.add_custom_axes(plotter)
plotter.add_title("Thresholded slice in MPS-coordinates", font_size=6)
plotter.screenshot("select_slice3.png");
pyvista.close_all()
image = np.concatenate([imageio.v3.imread(f"select_slice{i}.png") for i in range(4)], axis=1)
imageio.v3.imwrite("select_slice_example.png", image)
b64 = base64.b64encode(open("select_slice_example.png",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))
3D Shepp Logan-Phantom#
if the package phantominator is installed the RegularGridDataset can directly be used to create a 3D Shepp-Logan phantom via classmethod:
[2]:
dataset = cmrsim.datasets.RegularGridDataset.from_shepp_logan(map_size=(200, 200, 100), extent=Quantity([20, 20, 10], "cm"))
[3]:
import ipywidgets as widgets
hbox = widgets.HBox([widgets.Output(), widgets.Output()])
display(hbox)
with hbox.children[0]:
display(dataset.mesh)
with hbox.children[1]:
pyvista.close_all()
pyvista.start_xvfb()
plotter = pyvista.Plotter(notebook=False, off_screen=True, window_size=(600, 400))
plotter.add_mesh(dataset.mesh.slice_orthogonal(), scalars="rho", opacity=0.65, cmap="gray", show_scalar_bar=False)
local_functions.add_custom_axes(plotter)
plotter.screenshot("Offres.png")
b64 = base64.b64encode(open("Offres.png",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))
Get point-set dictionary#
[9]:
dataset.mesh["is_non_trivial"] = dataset.mesh["rho"] > 0.
phantom_dict = dataset.get_phantom_def(filter_by="is_non_trivial", keys=["rho", "T1", "T2"],
perturb_positions_std=Quantity(0.1, "mm").m_as("m"))
[10]:
polydata = pyvista.PolyData(phantom_dict["r_vectors"])
polydata["T1"] = phantom_dict["T1"]
pyvista.close_all()
pyvista.start_xvfb()
plotter = pyvista.Plotter(notebook=False, off_screen=True, window_size=(900, 400))
plotter.add_mesh(polydata, scalars="T1", opacity=0.65, cmap="gray", show_scalar_bar=False)
local_functions.add_custom_axes(plotter)
plotter.screenshot("polydata_temp.png")
b64 = base64.b64encode(open("polydata_temp.png",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))
[ ]: