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

UniformGrid	Information
N Cells	2282544
N Points	2334948
X Bounds	-9.985e-02, 9.815e-02
Y Bounds	-9.931e-02, 9.719e-02
Z Bounds	-1.000e-01, 9.800e-02
Dimensions	133, 132, 133
Spacing	1.500e-03, 1.500e-03, 1.500e-03
N Arrays	6

Name	Field	Type	N Comp	Min	Max
in_mesh	Points	float64	1	0.000e+00	1.000e+00
vtkValidPointMask	Points	int8	1	0.000e+00	1.000e+00
chi	Points	float64	1	-9.050e-06	3.600e-07
M0	Points	float64	1	0.000e+00	1.000e+00
T1	Points	float64	1	0.000e+00	1.000e+03
T2	Points	float64	1	0.000e+00	3.000e+02

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

UniformGrid	Information
N Cells	61161800
N Points	61628688
X Bounds	-9.985e-02, 9.765e-02
Y Bounds	-9.931e-02, 9.669e-02
Z Bounds	-1.000e-01, 9.750e-02
Dimensions	396, 393, 396
Spacing	5.000e-04, 5.000e-04, 5.000e-04
N Arrays	9

Name	Field	Type	N Comp	Min	Max
M0	Points	float64	1	0.000e+00	1.000e+00
T1	Points	float64	1	0.000e+00	1.000e+03
T2	Points	float64	1	0.000e+00	3.000e+02
chi	Points	float64	1	-9.050e-06	3.600e-07
in_mesh	Points	float64	1	0.000e+00	1.000e+00
offres	Points	float32	1	-3.505e-02	1.895e-02
vtkValidPointMask	Points	int8	1	1.000e+00	1.000e+00
vtkGhostType	Points	uint8	1	0.000e+00	0.000e+00
vtkGhostType	Cells	uint8	1	0.000e+00	0.000e+00

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}" />'))

[ ]: