We will wrap up the exercises for this session with a type of pulse design that visibly leverages some of the nicer abstractions within SigPy: parallel transmit pulse design.
Specifically, in this exercise we will focus on a spatial-domain spiral parallel transmit pulse design.
We start by writing our import statements one last time:
%matplotlib notebook
# typical sigpy and numpy imports
import numpy as np
import sigpy.mri as mr
import sigpy.mri.rf as rf # importing our rf tools separately
import sigpy.plot as pl
# to assist with importing data, we will also import scipy
import scipy.io as sio
import scipy.ndimage.filters as filt
import matplotlib.pyplot as mplib
For this problem, use a small problem geometry: a 2D 32x32 design grid with 8 channels, to keep computation simple for Binder! However, if you are running these exercises locally and have more RAM available, larger designs are no problem.
dim = 32
Nc = 8
sens_shape = [Nc, dim, dim]
img_shape = [dim, dim]
To begin, we will want some two-dimensional spatial magnetization pattern to design for. We provide one pattern in a matlab file for use in the data folder of this notebook, which we will load here. This is just to get you started - feel free to load other patterns or create your own!
mat_struct = sio.loadmat('data/smallv.mat')
d = mat_struct['d'].astype(np.single)
pl.ImagePlot(d, title='Target Excitation Pattern')
This is an interesting pattern, but the edges are extremely sharp. To avoid introducing Gibbs ringing into our design from these sharp edges, we will perform a very slight blur of the pattern using scipy. Our scipy filter requires real-valued numbers but we will want to use complex values in the design later, so we will cast to complex after the blurring is performed.
d = filt.gaussian_filter(d, 0.5)
d = d.astype(np.complex)
pl.ImagePlot(d, title='Blurred Target')
Here, we will break out the more general sigpy.mri submodule functions for the first time. Use sigpy.mri.birdcage_maps() to simulate our B1+ transmit sensitivities. Note that these simulated profiles would be equally valid to use as receive sensitivities for a reconstruction problem, such as a SENSE recon.
sens = mr.birdcage_maps(sens_shape)
pl.ImagePlot(sens)
If you haven't already, take a moment to play around with some of the hotkey features of sigpy.plot using the (complex-valued) sensitivities plot produced above. After clicking on the plot, try the following hotkeys:
You should be able to look at each of the coils' B1 sensitivities sequentially, and view their magnitude, phase, or real/imaginary components!
Before performing our design, we will also need to design time-varying gradient waveforms to produce a trajectory for our excitation.
Design a spiral-in trajectory using sigpy.mri.rf.spiral_arch(). Note that by default this trajectory designer produces a spiral-out trajectory - you will need to flip the waveforms to spiral inward! Numpy has functions that can assist with this, or you can use python built-ins (such as vector = vector[::-1])
Given:
Hint: to under- or oversample with the spiral use fov = fov/R
fov = 0.55 # FOV in m
gts = 6.4e-6 # hardware dwell time, s
gslew = 150 # gradient slew rate in mT/m/ms
gamp = 30 # maximum gradient amplitude in mT/m
R = 1/2 # degree of undersampling of outer region of trajectory- let's oversample by a factor of 2
dx = 0.025 # in m
# construct a trajectory
g, k, t, s = rf.spiral_arch(fov/R,dx,gts,gslew,gamp)
#Note that this trajectory is a spiral-out trajectory.
#Simply time-reverse it to create a spiral-in.
k = k[::-1]
g = g[::-1]
mplib.figure()
mplib.plot(k[:,0],k[:,1], color='orange')
mplib.title('Constant density spiral-in trajectory kx and ky')
mplib.figure()
mplib.plot(g)
mplib.title('Gradient waveforms')
Do this using the sigpy.mri.rf.stspa() small tip spatial-domain pulse designer. To keep the design fast, set explicit=False and perform ~ 10 iterations (explicit=True has not performed well in Binder). This should be sufficient to produce a good pattern.
Insert code for the pulse design here:
pulses = rf.stspa(d, sens, k, gts, explicit=False, max_iter=10)
pl.LinePlot(pulses)
One way of checking that our designed pulse is reasonable is by constructing the system matrix A used within the SmallTipSpatialPulseDesign app, and looking at $$ m = A b$$
Applying and transposing Linops is documented in detail on the Linop documentation page.
The linear operator is converted to its adjoint with A = A.H
You can apply the operation y = A(x) with y = A * x
A = mr.linop.Sense(sens, coord=k, ishape=d.shape).H
m = A * pulses
pl.ImagePlot(m, title = 'Mxy')
A better way of verifying the accuracy of the pulses is by performing a bloch simulation. Do so using the sigpy.mri.rf.sim.abrm_ptx() function.
Plot the Mxy and Mz profile once you're done.
#make spatial variables, in units of m
x, y = np.meshgrid(np.linspace(fov/2, -fov/2, dim), np.linspace(fov/2, -fov/2, dim))
spatial = np.fliplr(np.concatenate((np.reshape(x, (dim*dim, 1)), np.reshape(y, (dim*dim, 1))), axis=1))
gam = 42580 # Hz/mT
a, b, m, mz = rf.abrm_ptx(pulses/33, spatial, g * gam * gts * 2 * np.pi, gts, fmap=None, sens=sens)
pl.ImagePlot(np.reshape(2 * a * np.conj(b),img_shape), mode='m', title=('Mxy'))
pl.ImagePlot(np.reshape(1-2*np.absolute(b)**2,img_shape), mode='m', title=('Mz'))