Source code for sigpy.mri.rf.ptx

# -*- coding: utf-8 -*-
"""MRI RF excitation pulse design functions,
    including SLR and small tip spatial design
"""

import sigpy as sp
import numpy as np
from sigpy.mri import rf as rf
from sigpy import backend
from scipy.interpolate import interp1d


__all__ = ['stspa', 'stspk']


[docs]def stspa(target, sens, coord, dt, roi=None, alpha=0, b0=None, tseg=None, st=None, phase_update_interval=float('inf'), explicit=False, max_iter=1000, tol=1E-6): """Small tip spatial domain method for multicoil parallel excitation. Allows for constrained or unconstrained designs. Args: target (array): desired magnetization profile. [dim dim] sens (array): sensitivity maps. [Nc dim dim] coord (array): coordinates for noncartesian trajectories. [Nt 2] dt (float): hardware sampling dwell time. roi (array): array for error weighting, specify spatial ROI. [dim dim] alpha (float): regularization term, if unconstrained. b0 (array): B0 inhomogeneity map [dim dim]. For explicit matrix building. tseg (None or Dictionary): parameters for time-segmented off-resonance correction. Parameters are 'b0' (array), 'dt' (float), 'lseg' (int), and 'n_bins' (int). Lseg is the number of time segments used, and n_bins is the number of histogram bins. st (None or Dictionary): 'subject to' constraint parameters. Parameters are avg power 'cNorm' (float), peak power 'cMax' (float), 'mu' (float), 'rhoNorm' (float), 'rhoMax' (float), 'cgiter' (int), 'max_iter' (int), 'L' (list of arrays), 'c' (float), 'rho' (float), and 'lam' (float). These parameters are explained in detail in the SDMM documentation. phase_update_interval (int): number of iters between exclusive phase updates. If 0, no phase updates performed. explicit (bool): Use explicit matrix. max_iter (int): max number of iterations. tol (float): allowable error. Returns: array: pulses out. References: Grissom, W., Yip, C., Zhang, Z., Stenger, V. A., Fessler, J. A. & Noll, D. C.(2006). Spatial Domain Method for the Design of RF Pulses in Multicoil Parallel Excitation. Magnetic resonance in medicine, 56, 620-629. """ Nc = sens.shape[0] Nt = coord.shape[0] device = backend.get_device(target) xp = device.xp with device: pulses = xp.zeros((Nc, Nt), xp.complex64) # set up the system matrix if explicit: A = rf.linop.PtxSpatialExplicit(sens, coord, dt, target.shape, b0) else: A = sp.mri.linop.Sense(sens, coord, weights=None, tseg=tseg, ishape=target.shape).H # handle the Ns * Ns error weighting ROI matrix W = sp.linop.Multiply(A.oshape, xp.ones(target.shape)) if roi is not None: W = sp.linop.Multiply(A.oshape, roi) # apply ROI A = W * A # Unconstrained, use conjugate gradient if st is None: I = sp.linop.Identity((Nc, coord.shape[0])) b = A.H * W * target alg_method = sp.alg.ConjugateGradient(A.H * A + alpha * I, b, pulses, P=None, max_iter=max_iter, tol=tol) # Constrained case, use SDMM else: # vectorize target for SDMM target = W * target d = xp.expand_dims(target.flatten(), axis=0) alg_method = sp.alg.SDMM(A, d, st['lam'], st['L'], st['c'], st['mu'], st['rho'], st['rhoMax'], st['rhoNorm'], 10**-5, 10**-2, st['cMax'], st['cNorm'], st['cgiter'], st['max_iter']) # perform the design: apply optimization method to find solution pulse while not alg_method.done(): # phase_update switch if (alg_method.iter > 0) and \ (alg_method.iter % phase_update_interval == 0): target = xp.abs(target) * xp.exp( 1j * xp.angle( xp.reshape(A * alg_method.x, target.shape))) b = A.H * target alg_method.b = b alg_method.update() if st is not None: pulses = xp.reshape(alg_method.x, [Nc, Nt]) return pulses
[docs]def stspk(mask, sens, n_spokes, fov, dx_max, gts, sl_thick, tbw, dgdtmax, gmax, alpha=1, iter_dif=0.01): """Small tip spokes parallel transmit pulse designer. Args: mask (ndarray): region in which to optimize flip angle uniformity in slice. [dim dim] sens (ndarray): sensitivity maps. [nc dim dim] n_spokes (int): number of spokes to be created in the design. fov (float): excitation FOV (cm). dx_max (float): max. resolution of the trajectory (cm). gts (float): hardware sampling dwell time (s). sl_thick (float): slice thickness (mm). tbw (int): time-bandwidth product. dgdtmax (float): max gradient slew (g/cm/s). gmax (float): max gradient amplitude (g/cm). alpha (float): regularization parameter. iter_dif (float): for each spoke, the difference in cost btwn. successive iterations at which to terminate MLS iterations. Returns: 2-element tuple containing - **pulses** (*array*): RF waveform out. - **g** (*array*): corresponding gradient, in g/cm. References: Grissom, W., Khalighi, M., Sacolick, L., Rutt, B. & Vogel, M (2012). Small-tip-angle spokes pulse design using interleaved greedy and local optimization methods. Magnetic Resonance in Medicine, 68(5), 1553-62. """ device = backend.get_device(sens) xp = device.xp with device: nc = sens.shape[0] kmax = 1 / dx_max # /cm, max spatial freq of trajectory # greedy kx, ky grid kxs, kys = xp.meshgrid(xp.linspace(-kmax / 2, kmax / 2 - 1 / fov, int(fov * kmax)), xp.linspace(-kmax / 2, kmax / 2 - 1 / fov, int(fov * kmax))) # vectorize the grid kxs = kxs.flatten() kys = kys.flatten() # remove DC dc = xp.intersect1d(xp.where((kxs == 0)), xp.where((kys == 0)))[0] kxs = xp.concatenate([kxs[:dc], kxs[dc+1:]]) kys = xp.concatenate([kys[:dc], kys[dc+1:]]) # step 2: design the weights # initial kx/ky location is DC k = xp.expand_dims(xp.array([0, 0]), 0) # initial target phase phs = xp.zeros((xp.count_nonzero(mask), 1), dtype=xp.complex64) for ii in range(n_spokes): # build Afull (and take only 0 locations into matrix) Anum = rf.PtxSpatialExplicit(sens, k, gts, mask.shape, ret_array=True) Anum = Anum[~(Anum == 0).all(1)] # design wfull using MLS: # initialize wfull sys_a = (Anum.conj().T @ Anum + alpha * xp.eye((ii+1)*nc)) sys_b = (Anum.conj().T @ xp.exp(1j*phs)) w_full = xp.linalg.solve(sys_a, sys_b) err = Anum @ w_full - xp.exp(1j * phs) cost = err.conj().T @ err + alpha * w_full.conj().T @ w_full cost = xp.real(cost) cost_old = 10 * cost # to get the loop going while xp.absolute(cost - cost_old) > iter_dif * cost_old: cost_old = cost phs = xp.angle(Anum @ w_full) w_full = xp.linalg.solve( (Anum.conj().T @ Anum + alpha * xp.eye((ii + 1) * nc)), (Anum.conj().T @ xp.exp(1j * phs))) err = Anum @ w_full - xp.exp(1j * phs) cost = xp.real(err.conj().T @ err + alpha * w_full.conj().T @ w_full) # add a spoke using greedy method if ii < n_spokes - 1: r = xp.exp(1j * phs) - Anum @ w_full rfnorm = xp.zeros(kxs.shape, dtype=xp.complex64) for jj in range(kxs.size): ks_test = xp.expand_dims(xp.array([kxs[jj], kys[jj]]), 0) Anum = rf.PtxSpatialExplicit(sens, ks_test, gts, mask.shape, ret_array=True) Anum = Anum[~(Anum == 0).all(1)] rfm = xp.linalg.solve((Anum.conj().T @ Anum), (Anum.conj().T @ r)) rfnorm[jj] = xp.linalg.norm(rfm) ind = xp.argmax(rfnorm) k_new = xp.expand_dims(xp.array([kxs[ind], kys[ind]]), 0) if ii % 2 != 0: # add to end of pulse k = xp.concatenate((k, k_new)) else: # add to beginning of pulse k = xp.concatenate((k_new, k)) # remove chosen point from candidates kxs = xp.concatenate([kxs[:ind], kxs[ind + 1:]]) kys = xp.concatenate([kys[:ind], kys[ind + 1:]]) # from our spoke selections, build the whole waveforms # first, design our gradient waveforms: g = rf.spokes_grad(k, tbw, sl_thick, gmax, dgdtmax, gts) # design our rf # calc. the size of the traps in our gz waveform- will use to calc rf area = tbw / (sl_thick / 10) / 4257 # thick*kwid=twb, kwid=gam*area [subgz, nramp] = rf.min_trap_grad(area, gmax, dgdtmax, gts) npts = 128 subrf = rf.dzrf(npts, tbw, 'st') n_plat = subgz.size - 2 * nramp # time points on trap plateau # interpolate to stretch out waveform to appropriate length f = interp1d(np.arange(0, npts, 1) / npts, subrf, fill_value='extrapolate') subrf = f(xp.arange(0, n_plat, 1) / n_plat) subrf = xp.concatenate((xp.zeros(nramp), subrf, xp.zeros(nramp))) pulses = xp.kron(xp.reshape(w_full, (nc, n_spokes)), subrf) # add zeros for gzref rf_ref = xp.zeros((nc, g.shape[1] - pulses.shape[1])) pulses = xp.concatenate((pulses, rf_ref), 1) return pulses, g