Merge pull request #7 from mjt320/develop

magnetisation transfer and inversion recovery T1

README.md

@@ -18,10 +18,11 @@ Created 28 September 2020
 - AIFs: including patient-specific (measured), Parker, bi-exponential Parker
 - Fitting free AIF time delay parameter
 - Relaxivity models: linear
-- Signal models: spoiled gradient echo
+- Signal models: spoiled gradient echo, inversion-recovery spin-echo
 - Water exchange models: FXL, NXL, NXL_be
-- T1 fitting using variable flip angle method, IR-SPGR and DESPOT1-HIFI
+- T1 fitting using variable flip angle method, IR-SPGR, DESPOT1-HIFI and inversion recovery 
 - T2(*) fitting for multi-TE acquisitions
+- MTR and MTSat calculation
 
 ### Not yet implemented/limitations:
 - Additional pharmacokinetic models (add by inheriting from PkModel class)
@@ -30,7 +31,4 @@ Created 28 September 2020
 - Additional signal models (add by inheriting from SignalModel class)
 - R2/R2* effects not included in fitting of enhancement curves (but is included for enhancement-to-concentration conversion)
 - Compartment-specific relaxivity parameters/models
-- Fitting free water exchange parameters
-
-### TODO:
-- inversion recovery T1 measurment
+- Fitting water exchange parameters

src/mt.py

@@ -0,0 +1,130 @@
+"""Functions to fit MRI SPGR signal to obtain MT parameters.
+
+Created 13 November 2023
+@authors: Michael Thrippleton
+@email: [email protected]
+@institution: University of Edinburgh, UK
+
+Classes:
+    MTR
+    MTSat
+
+"""
+
+import numpy as np
+from scipy.optimize import least_squares
+from fitting import Fitter
+
+
+class MTR(Fitter):
+    """Magnetisation transfer ratio calculation.
+
+    Calculates MTR based on MT ON and MT OFF acquisitions.
+    Subclass of Fitter.
+    """
+
+    def __init__(self):
+        """
+
+        Args:
+        """
+        pass
+
+    def output_info(self):
+        """Get output info. Overrides superclass method.
+        """
+        return ('mtr', False),
+
+    def proc(self, s):
+        """Calculate MTR. Overrides superclass method.
+
+        Args:
+            s (ndarray): 1D np array containing the MT ON, MT OFF signals
+
+        Returns:
+            tuple: (mtr)
+                mtr (float): magnetisation transfer ratio (percent)
+        """
+        if any(np.isnan(s)):
+            raise ArithmeticError(
+                f'Unable to calculate MTR: nan signal received.')
+
+        if s.shape != (2,):
+            raise ValueError(
+                f'Argument s should have shape (2,)'
+            )
+
+        return 100*(s[1]-s[0])/s[1]
+
+
+class MTSat(Fitter):
+    """Magnetisation transfer saturation calculation.
+
+    Calculates MTSat, MTR and T1 based on three acquisitions with MT, PD and
+    T1 weightings.
+    Subclass of Fitter.
+    """
+
+    def __init__(self, fa_mt, fa_pd, fa_t1, tr_pd, tr_t1):
+        """
+
+        Args:
+            fa_mt (float): excitation flip angle for MTw acquisition (deg)
+            fa_pd (float): excitation flip angle for PDw acquisition (deg)
+            fa_t1 (float): excitation flip angle for T1w acquisition (deg)
+            tr_pd (float): TR for PD and MT weighted acquisitions (s)
+            tr_t1 (float): TR for T1 weighted acquisition (s)
+        """
+        self.fa_mt = fa_mt * np.pi / 180
+        self.fa_pd = fa_pd * np.pi / 180
+        self.fa_t1 = fa_t1 * np.pi / 180
+        self.tr_pd = tr_pd
+        self.tr_t1 = tr_t1
+
+    def output_info(self):
+        """Get output info. Overrides superclass method.
+        """
+        return ('mtsat', False), ('mtr', False), ('t1', False), ('r1', False), \
+               ('a', False)
+
+    def proc(self, s):
+        """Calculate MTSat, MTR and T1. Overrides superclass method.
+        Uses equations from:
+            Helms et al., Magnetic Resonance in Medicine 60:1396–1407 (2008)
+            https://doi.org/10.1002/mrm.21732
+            and erratum https://doi.org/10.1002/mrm.22607
+
+        Args:
+            s (ndarray): 1D np array containing the MT, PD and T1
+            weighted signals in that sequence
+
+        Returns:
+            tuple: (mtsat, mtr and t1)
+                mtsat (float): magnetisation transfer saturation
+                mtr (float): magnetisation transfer ratio (percent)
+                t1 (float): T1 estimation (s)
+                r1 (float): R1 estimation (s^-1)
+                a (float): amplitude estimation
+        """
+        if any(np.isnan(s)):
+            raise ArithmeticError(
+                f'Unable to calculate MTSat: nan signal received.')
+
+        if s.shape != (3,):
+            raise ValueError(
+                f'Argument s should have shape (3,)'
+            )
+        s_mt, s_pd, s_t1 = s
+        mtr = 100*(s_pd-s_mt)/s_pd
+        a = (s_t1*(s_pd*self.tr_t1*self.fa_pd**2 -
+                   s_pd*self.tr_pd*self.fa_t1**2)) / (
+            self.fa_pd*self.fa_t1*(s_pd*self.tr_t1*self.fa_pd -
+                                 s_t1*self.tr_pd*self.fa_t1))
+        r1 = -(self.fa_pd*self.fa_t1*(s_pd*self.tr_t1*self.fa_pd -
+                                    s_t1*self.tr_pd*self.fa_t1)) / (
+                2*self.tr_t1*self.tr_pd*(s_pd*self.fa_t1 - s_t1*self.fa_pd))
+        mtsat = -self.fa_pd**2/2 + (a*r1*self.tr_pd*self.fa_pd)/s_mt - \
+                r1*self.tr_pd
+        t1 = 1/r1
+
+        return mtsat, mtr, t1, r1, a

src/signal_models.py

@@ -7,6 +7,7 @@
 
 Classes: SignalModel and derived subclasses:
     SPGR
+    IRSE
 """
 
 from abc import ABC, abstractmethod
@@ -74,3 +75,38 @@ def R_to_s(self, s0, R1, R2=None, R2s=0, k_fa=1):
                   (1.0-np.exp(-self.tr*R1)*np.cos(fa))
                   ) * np.exp(-self.te*R2s)
         return s
+
+
+class IRSE(SignalModel):
+    """Signal model subclass for inversion-recovery spin-echo pulse sequence.
+    Ignores impact of refocusing pulse on Mz.
+
+    """
+
+    def __init__(self, tr, ti, fa_ex, fa_inv, te, signed_signal=False):
+        """
+
+        Args:
+            tr (float): repetition time (s)
+            ti (float): inversion time (s)
+            fa_ex (float): excitation flip angle (deg)
+            fa_inv (float): inversion flip angle (deg)
+            te (float): echo time (s)
+            signed_signal (bool): whether signal can be + or -
+        """
+        self.tr = tr
+        self.ti = ti
+        self.fa_ex = fa_ex * np.pi/180
+        self.fa_inv = fa_inv * np.pi/180
+        self.te = te
+        self.signed_signal = signed_signal
+
+    def R_to_s(self, s0, R1, R2=0, R2s=None, k_fa=1):
+        """Get signal for this model. Overrides superclass method."""
+        fa_ex = k_fa * self.fa_ex
+        fa_inv = k_fa * self.fa_inv
+        s = s0 * np.exp(-self.te*R2) * np.sin(fa_ex) * (
+                1-np.cos(fa_inv)*np.exp(-self.tr*R1)-(1-np.cos(fa_inv))
+                * np.exp(-self.ti*R1)) / (
+                1-np.cos(fa_inv)*np.cos(fa_ex)*np.exp(-self.tr*R1))
+        return s if self.signed_signal else np.abs(s)

src/t1_fit.py

@@ -10,6 +10,7 @@
     VFALinear
     VFANonLinear
     HIFI
+    IRSE
 
 Functions:
     spgr_signal: get SPGR signal
@@ -19,6 +20,7 @@
 import numpy as np
 from scipy.optimize import least_squares
 from fitting import Fitter
+from utils.utilities import least_squares_global
 
 
 class VFA2Points(Fitter):
@@ -313,6 +315,102 @@ def __signal(self, x):
         return s
 
 
+class IRSE(Fitter):
+    """T1 estimation from inversion-recovery spin-echo.
+
+    Subclass of Fitter.
+    For fitting a series of acquisitions with the same TR/TE and different TI.
+    """
+
+    MIN_T1 = 0.01  # minimum allowable T1 (needed to avoid exceptions)
+
+    def __init__(self, tr, ti, pars_0=None, signed_signal=False, n_tries=1):
+        """
+
+        Args:
+            tr (float): TR for all acquisition (s)
+            ti (ndarray): TI for each acquisition (s)
+            pars_0 (list, optional): 2-element list of initial parameters s0 and
+                T1. If None (default), initial s0 is estimated as the maximum
+                absolute signal and initial T1 is estimated from the minimum
+                absolute signal.
+            signed_signal (bool, optional): image includes positive and negative
+                intensities, otherwise it is assumed to be a magnitude image.
+            n_tries (int, optional): if >1, repeats optimisation with
+                randomised starting values (uniform random values between 5%
+                and 2x initial estimates). Defaults to 1.
+        """
+        self.tr = tr
+        self.ti = ti
+        self.pars_0 = pars_0
+        self.signed_signal = signed_signal
+        self.n_tries = n_tries
+
+    def output_info(self):
+        """Get output info. Overrides superclass method.
+        """
+        return ('a', False), ('b', False), ('t1', False), ('s_opt', True)
+
+    def proc(self, s):
+        """Estimate T1 and s0. Overrides superclass method.
+        Uses the 3-parameter model: s = a + b exp(-TI/T1)
+
+        Args:
+            s (ndarray): 1D array containing the signals
+
+        Returns:
+            tuple: s0, t1, s_opt
+                a (float): estimated "a" parameter
+                b (float): estimated "b" parameter
+                t1 (float): estimated T1 (s)
+                s_opt (ndarray): fitted signal intensities
+        """
+        # if s_opts not given, generate starting estimates for s0, T1, a and b
+        if self.pars_0 is None:
+            s0_init = np.max(np.abs(s))  # order-of-magnitude estimate for s0
+            t1_init = self.ti[np.argmin(np.abs(s))]/np.log(2)  # ~null signal
+        else:
+            s0_init, t1_init = self.pars_0
+        t1_init = max(type(self).MIN_T1, t1_init)
+        a_init = s0_init * (1 + np.exp(-self.tr/t1_init))
+        b_init = -2 * s0_init
+        x_0 = [np.array((a_init, b_init, t1_init))]
+        # if >1 tries, generate additional randomised starting parameters
+        for n in range(1, self.n_tries):
+            s0_init_rnd = np.random.uniform(0.05, 2) * s0_init
+            t1_init_rnd = np.random.uniform(0.05, 2) * t1_init
+            t1_init_rnd = max(type(self).MIN_T1, t1_init_rnd)
+            a_init_rnd = s0_init_rnd * (1 + np.exp(-self.tr/t1_init_rnd))
+            b_init_rnd = -2 * s0_init_rnd
+            x_0_init_rnd = [np.array((a_init_rnd, b_init_rnd, t1_init_rnd))]
+            x_0 = x_0 + x_0_init_rnd
+
+
+        # optimise using the 3-parameter model
+        x_scale = (np.abs(a_init), np.abs(b_init), 1.0)
+        bounds = ((-np.inf, -np.inf, type(self).MIN_T1), (np.inf, np.inf, np.inf))
+        result = least_squares_global(self.__residuals, x_0,
+                                      args=(s,), method='trf',
+                                      bounds=bounds,
+                                      x_scale=x_scale)
+        if result.success is False:
+            raise ArithmeticError(
+                f'Unable to calculate T1'
+                f': {result.message}')
+
+        a_opt, b_opt, t1_opt = result.x
+        s_opt = self.__signal(result.x)
+
+        return a_opt, b_opt, t1_opt, s_opt
+
+    def __residuals(self, x, s):
+        return s - self.__signal(x)
+
+    def __signal(self, x):
+        signal = x[0] + x[1] * np.exp(-self.ti / x[2])
+        return signal if self.signed_signal else np.abs(signal)
+
+
 def spgr_signal(s0, t1, tr, fa):
     """Return signal for SPGR sequence.