Source code for sigpy.mri.rf.sim

# -*- coding: utf-8 -*-
"""RF Pulse Simulation Functions.

"""
from sigpy import backend

__all__ = ['abrm', 'abrm_nd', 'abrm_hp', 'abrm_ptx']


[docs]def abrm(rf, x, balanced=False): r"""1D RF pulse simulation, with simultaneous RF + gradient rotations. Args: rf (array): rf waveform input. x (array): spatial locations. balanced (bool): toggles application of rewinder. Returns: 2-element tuple containing - **a** (*array*): SLR alpha parameter. - **b** (*array*): SLR beta parameter. References: Pauly, J., Le Roux, Patrick., Nishimura, D., and Macovski, A.(1991). 'Parameter Relations for the Shinnar-LeRoux Selective Excitation Pulse Design Algorithm'. IEEE Transactions on Medical Imaging, Vol 10, No 1, 53-65. """ device = backend.get_device(rf) xp = device.xp with device: eps = 1e-16 g = xp.ones(xp.size(rf)) * 2 * xp.pi / xp.size(rf) a = xp.ones(xp.size(x), dtype=complex) b = xp.zeros(xp.size(x), dtype=complex) for mm in range(xp.size(rf)): om = x * g[mm] phi = xp.sqrt(xp.abs(rf[mm]) ** 2 + om ** 2) + eps n = xp.column_stack((xp.real(rf[mm]) / phi, xp.imag(rf[mm]) / phi, om / phi)) av = xp.cos(phi / 2) - 1j * n[:, 2] * xp.sin(phi / 2) bv = -1j * (n[:, 0] + 1j * n[:, 1]) * xp.sin(phi / 2) at = av * a - xp.conj(bv) * b bt = bv * a + xp.conj(av) * b a = at b = bt if balanced: # apply a rewinder g = -2 * xp.pi / 2 om = x * g phi = xp.abs(om) + eps nz = om / phi av = xp.cos(phi / 2) - 1j * nz * xp.sin(phi / 2) a = av * a b = xp.conj(av) * b return a, b
[docs]def abrm_nd(rf, x, g): r"""N-dim RF pulse simulation Assumes that x has inverse spatial units of g, and g has gamma*dt applied. Assumes dimensions x = [...,Ndim], g = [Ndim,Nt]. Args: rf (array): rf waveform input. x (array): spatial locations. g (array): gradient array. Returns: 2-element tuple containing - **a** (*array*): SLR alpha parameter. - **b** (*array*): SLR beta parameter. References: Pauly, J., Le Roux, Patrick., Nishimura, D., and Macovski, A.(1991). 'Parameter Relations for the Shinnar-LeRoux Selective Excitation Pulse Design Algorithm'. IEEE Transactions on Medical Imaging, Vol 10, No 1, 53-65. """ device = backend.get_device(rf) xp = device.xp with device: eps = 1e-16 a = xp.ones(xp.shape(x)[0], dtype=complex) b = xp.zeros(xp.shape(x)[0], dtype=complex) for mm in range(xp.size(rf)): om = x @ g[mm, :] phi = xp.sqrt(xp.abs(rf[mm]) ** 2 + om ** 2) n = xp.column_stack((xp.real(rf[mm]) / (phi + eps), xp.imag(rf[mm]) / (phi + eps), om / (phi + eps))) av = xp.cos(phi / 2) - 1j * n[:, 2] * xp.sin(phi / 2) bv = -1j * (n[:, 0] + 1j * n[:, 1]) * xp.sin(phi / 2) at = av * a - xp.conj(bv) * b bt = bv * a + xp.conj(av) * b a = at b = bt return a, b
[docs]def abrm_hp(rf, gamgdt, xx, dom0dt=0): r"""1D RF pulse simulation, with non-simultaneous RF + gradient rotations. Args: rf (array): rf pulse samples in radians. gamdt (array): gradient samples in radians/(units of xx). xx (array): spatial locations. dom0dt (float): off-resonance phase in radians. Returns: 2-element tuple containing - **a** (*array*): SLR alpha parameter. - **b** (*array*): SLR beta parameter. References: Pauly, J., Le Roux, Patrick., Nishimura, D., and Macovski, A.(1991). 'Parameter Relations for the Shinnar-LeRoux Selective Excitation Pulse Design Algorithm'. IEEE Transactions on Medical Imaging, Vol 10, No 1, 53-65. """ device = backend.get_device(rf) xp = device.xp with device: Ns = xp.shape(xx) Ns = Ns[0] # Ns: # of spatial locs Nt = xp.shape(gamgdt) Nt = Nt[0] # Nt: # time points a = xp.ones((Ns,)) b = xp.zeros((Ns,)) for ii in xp.arange(Nt): # apply phase accural z = xp.exp(-1j * (xx * gamgdt[ii, ] + dom0dt)) b = b * z # apply rf C = xp.cos(xp.abs(rf[ii]) / 2) S = 1j * xp.exp(1j * xp.angle(rf[ii])) * xp.sin(xp.abs(rf[ii]) / 2) at = a * C - b * xp.conj(S) bt = a * S + b * C a = at b = bt z = xp.exp(1j / 2 * (xx * xp.sum(gamgdt, axis=0) + Nt * dom0dt)) a = a * z b = b * z return a, b
[docs]def abrm_ptx(b1, x, g, dt, fmap=None, sens=None): r"""N-dim RF pulse simulation Assumes that x has inverse spatial units of g, and g has gamma*dt applied. Assumes dimensions rf = [Nc, Nt], x = [...,Ndim], g = [Ndim,Nt], and sens = [Nc, dim, dim]. Args: b1 (array): rf waveform input samples in radians. x (array): spatial locations (m). g (array): gradient array (mT/m with gamma*dt applied). dt (float): hardware dwell time (s). fmap (array): off-resonance map (Hz). sens (array or None): B1+ sensitivity matrix. If None, creates matrix of 1's. Input size [Nc dim dim] Returns: 4-element tuple containing - **a** (*array*): SLR alpha parameter. - **b** (*array*): SLR beta parameter. - **m** (*array*): transverse magnetization. - **mz** (*array*): longitudinal magnetization. References: Pauly, J., Le Roux, Patrick., Nishimura, D., and Macovski, A.(1991). 'Parameter Relations for the Shinnar-LeRoux Selective Excitation Pulse Design Algorithm'. IEEE Transactions on Medical Imaging, Vol 10, No 1, 53-65. Grissom, W., Xu, D., Kerr, A., Fessler, J. and Noll, D. (2009). 'Fast large-tip-angle multidimensional and parallel RF pulse design in MRI' IEEE Trans Med Imaging, Vol 28, No 10, 1548-59. """ device = backend.get_device(b1) xp = device.xp with device: gam = 267.522 * 1e6 / 1000 # rad/s/mT dim = int(xp.sqrt(x.shape[0])) Ns = dim * dim Nc = b1.shape[0] Nt = b1.shape[1] dim = int(xp.sqrt(x.shape[0])) if sens is None: sens = xp.ones((dim*dim, Nc)) else: sens = xp.transpose(sens) sens = xp.reshape(sens, (dim*dim, Nc)) bxy = sens @ b1 bz = x @ xp.transpose(g) if fmap is not None and xp.sum(xp.abs(fmap)) != 0: rep_b0 = xp.repeat(xp.expand_dims(fmap.flatten(), 0), Nt, axis=0) bz += xp.transpose(rep_b0/gam*2*xp.pi) statea = xp.ones((Ns, 1)) stateb = xp.zeros((Ns, 1)) a = xp.ones(xp.shape(x)[0], dtype=complex) b = xp.zeros(xp.shape(x)[0], dtype=complex) for mm in range(Nt): phi = dt*gam*xp.sqrt(xp.abs(bxy[:, mm]) ** 2 + bz[:, mm] ** 2) with xp.errstate(divide='ignore'): normfact = dt*gam*(phi ** -1) normfact[xp.isinf(normfact)] = 0 nxy = normfact * bxy[:, mm] nxy[xp.isinf(nxy)] = 0 nz = normfact * bz[:, mm] nz[xp.isinf(nz)] = 0 cp = xp.cos(phi/2) sp = xp.sin(phi/2) alpha = xp.expand_dims(cp + 1j * nz * sp, 1) beta = xp.expand_dims(1j * xp.conj(nxy) * sp, 1) tmpa = xp.multiply(alpha, statea) + xp.multiply(beta, stateb) tmpb = -xp.conj(beta) * statea + xp.conj(alpha) * stateb statea, stateb = tmpa, tmpb # NOT returning all states: a = statea b = -xp.conj(stateb) mxy0 = 0 + 1j * 0 mz0 = 1 m = mz0 * xp.conj(statea) * stateb m += mxy0*xp.conj(statea) ** 2 m -= xp.conj(mxy0)*(stateb ** 2) mz = mz0 * (statea * xp.conj(statea) - stateb * xp.conj(stateb)) mz += 2 * xp.real(mxy0 * xp.conj(statea)*xp.negative(xp.conj(stateb))) return a, b, m, mz