"""Spoiled GRE simulation sub-routines."""
__all__ = ["SPGRModel"]
from ..base import AbstractModel
from ..base import autocast
import numpy.typing as npt
import torch
[docs]
class SPGRModel(AbstractModel):
"""
SPGR transverse signal at time TE after excitation.
This class models the transverse magnetization signal generated by the
spoiled gradient echo (SPGR) sequence, calculated at echo time (TE)
following RF excitation.
Methods
-------
set_properties(T1, T2star, M0=1.0, B0=0.0, chemshift=0.0):
Set tissue and system-specific properties for the SPGR model.
set_sequence(flip, TR, TE):
Set sequence parameters including flip angle, repetition time (TR),
and echo time (TE).
_engine(T1, T2star, TR, TE, flip, M0=1.0, field_map=0.0, delta_cs=0.0):
Compute the SPGR signal for given tissue, sequence, and field parameters.
Examples
--------
.. exec::
from torchsim.models import SPGRModel
model = SPGRModel()
model.set_properties(T1=1000, T2star=30)
model.set_sequence(flip=13.0, TR=10.0, TE=5.0)
signal = model()
"""
[docs]
@autocast
def set_properties(
self,
T1: float | npt.ArrayLike,
T2star: float | npt.ArrayLike,
M0: float | npt.ArrayLike = 1.0,
B0: float | npt.ArrayLike = 0.0,
chemshift: float | npt.ArrayLike = 0.0,
):
"""
Set tissue and system-specific properties for the SPGR model.
Parameters
----------
T1 : float | npt.ArrayLike
Longitudinal relaxation time in milliseconds.
T2star : float | npt.ArrayLike
Effective transverse relaxation time in milliseconds.
M0 : float | npt.ArrayLike, optional
Proton density scaling factor, default is ``1.0``.
B0 : float | npt.ArrayLike, optional
Frequency offset map in Hz, default is ``0.0.``
chemshift : float | npt.ArrayLik, optional
Chemical shift in Hz, default is ``0.0``.
"""
self.properties.T1 = T1
self.properties.T2star = T2star
self.properties.M0 = M0
self.properties.B0 = B0
self.properties.chemshift = chemshift
[docs]
@autocast
def set_sequence(
self,
flip: float | npt.ArrayLike,
TR: float | npt.ArrayLike,
TE: float | npt.ArrayLike,
):
"""
Set sequence parameters for the SPGR model.
Parameters
----------
flip : float | npt.ArrayLike
Flip angle in degrees.
TR : float | npt.ArrayLike
Repetition time in milliseconds.
TE : float | npt.ArrayLike
Echo time in milliseconds.
"""
self.sequence.flip = torch.pi * flip / 180.0
self.sequence.TR = TR * 1e-3 # ms -> s
self.sequence.TE = TE * 1e-3 # ms -> s
[docs]
@staticmethod
def _engine(
T1: float | npt.ArrayLike,
T2star: float | npt.ArrayLike,
TR: float | npt.ArrayLike,
TE: float | npt.ArrayLike,
flip: float | npt.ArrayLike,
M0: float | npt.ArrayLike = 1.0,
B0: float | npt.ArrayLike = 0.0,
chemshift: float | npt.ArrayLike = 0.0,
):
# Prepare relaxation parameters
R1, R2star = 1e3 / T1, 1e3 / T2star
# We are assuming Freeman-Hill convention for off-resonance map,
# so we need to negate to make use with this Ernst-Anderson-based implementation from Hoff
B0 = -B0
# Prepare off resonance
df = 2 * torch.pi * (B0 + chemshift)
# Divide-by-zero risk with PyTorch's nan_to_num
E1 = torch.exp(-R1 * TR)
E2 = torch.exp(-R2star * TE)
Phi = torch.exp(1j * df * TE)
# Precompute cos, sin
ca = torch.cos(flip)
sa = torch.sin(flip)
# Main calculation
den = 1 - E1 * ca
Mxy = M0 * ((1 - E1) * sa) / den
# Add decay
signal = Mxy * E2
# Add additional phase factor for readout at TE.
signal = signal * Phi
return signal