Merge pull request #204 from GSTT-CSC/t2-rician-noise

Updates T2 relaxometry fitting with Rician noise model

hazenlib/relaxometry.py

@@ -93,14 +93,19 @@
     measurements.
 6. Determine relaxation time (T1 or T2) by fitting the decay equation to
     the ROI data for each sphere. The published values of the relaxation
-    times are used to seed the optimisation algorithm. For T2 fitting the
-    input data are truncated for TE > 5*T2 to avoid fitting Rician noise in
-    magnitude images with low signal intensity. Optionally plot and save the
-    decay curves.
+    times are used to seed the optimisation algorithm. A Rician nose model is
+    used for T2 fitting [1]_. Optionally plot and save the decay curves.
 7. Return plate number, relaxation type (T1 or T2), measured relaxation
     times, published relaxation times, and fractional differences in a
     dictionary.
 
+References
+==========
+.. [1] Raya, J.G., Dietrich, O., Horng, A., Weber, J., Reiser, M.F.
+and Glaser, C., 2010. T2 measurement in articular cartilage: impact of the
+fitting method on accuracy and precision at low SNR. Magnetic Resonance in
+Medicine: An Official Journal of the International Society for Magnetic
+Resonance in Medicine, 63(1), pp.181-193.
 
 Feature enhancements
 ====================
@@ -117,19 +122,24 @@
 Get r-squared measure of fit.
 
 """
-import pydicom
+import os.path
+
 import cv2 as cv
-import numpy as np
 import matplotlib.pyplot as plt
-import os
-import os.path
-import skimage.morphology
+import numpy as np
+import pydicom
 import scipy.ndimage
 import scipy.optimize
+import skimage.morphology
 from scipy.interpolate import UnivariateSpline
+from scipy.special import i0e, ive
 
 import hazenlib.exceptions
 
+# Parameters for Rician noise model
+MAX_RICIAN_NOISE = 20.0
+SEED_RICIAN_NOISE = 5.0
+
 # Use dict to store template and reference information
 # Coordinates are in array format (row,col), rather than plt.patches 
 # format (col,row)
@@ -322,7 +332,7 @@ def pixel_rescale(dcmfile):
 
     For Philips scanners the private DICOM fields 2005,100d (=SI) and 2005,100e
     (=SS) are used as inverse scaling factors to perform the inverse
-    transformation [1]_
+    transformation [1]_.
 
     Parameters
     ----------
@@ -454,12 +464,18 @@ def est_t1_s0(ti, tr, t1, pv):
     return -pv / (1 - 2 * np.exp(-ti / t1) + np.exp(-tr / t1))
 
 
-def t2_function(te, t2, s0):
-    """
-    Calculated pixel value given TE, T2, S0 and C.
+def t2_function(te, t2, s0, c):
+    r"""
+    Calculated pixel value with Rician noise model.
     
-    Calculates pixel value from::
-        S = S0 * np.exp(-te / t2)
+    Calculates pixel value from [1]_::
+        .. math::
+            S=\sqrt{\frac{\pi \alpha^2}{2}} \exp(- \alpha) \left( (1+ 2 \alpha)
+            \  \text{I_0}(\alpha) + 2 \alpha \ \text{I_1}(\alpha) \right)
+
+            \alpha() = \left( \frac{S_0}{2 \sigma} \ \exp{\left(-\frac{\text{TE}}{\text{T}_2}\right)} \right)^2
+
+            \text{I}_n() = n^\text{th} \ \text{order modified Bessel function of the first kind}
 
     Parameters
     ----------
@@ -469,46 +485,39 @@ def t2_function(te, t2, s0):
         T2 decay constant.
     S0 : float
         Initial pixel magnitude.
+    C : float
+        Noise parameter for Rician model (equivalent to st dev).
 
     Returns
     -------
     pv : array_like
         Theoretical pixel values (signal) at each TE.
 
-    """
-    pv = s0 * np.exp(-te / t2)
-    return pv
-
-
-def t2_jacobian(te, t2, s0):
-    """
-    Jacobian of ``t2_function`` used for curve fitting.
-
-    Parameters
+    References
     ----------
-    te : array_like
-        Echo times.
-    t2 : float
-        T2 decay constant.
-    s0 : float
-        Initial signal magnitude.
+    .. [1] Raya, J.G., Dietrich, O., Horng, A., Weber, J., Reiser, M.F. and
+    Glaser, C., 2010. T2 measurement in articular cartilage: impact of the
+    fitting method on accuracy and precision at low SNR. Magnetic Resonance in
+    Medicine: An Official Journal of the International Society for Magnetic
+    Resonance in Medicine, 63(1), pp.181-193.
+    """
 
-    Returns
-    -------
-    array
-        [t2_der, s0_der].T, where x_der is a 1-D array of the partial
-        derivatives at each ``te``.
+    s0 = s0
+    alpha = (s0 / (2 * c) * np.exp(-te / t2)) **2
+    # NB need to use `i0e` and `ive` below to avoid numeric inaccuracy from
+    # multiplying by huge exponentials then dividing by the same exponential
+    pv = np.sqrt(np.pi/2 * c ** 2) *  \
+         ((1 + 2 * alpha) * i0e(alpha) + 2 * alpha * ive(1, alpha))
 
-    """
-    t2_der = s0 * te / t2 ** 2 * np.exp(-te / t2)
-    s0_der = np.exp(-te / t2)
-    jacobian = np.array([t2_der, s0_der])
-    return jacobian.T
+    return pv
 
 
 def est_t2_s0(te, t2, pv, c=0.0):
     """
-    Initial guess for s0 to seed curve fitting.
+    Initial guess for s0 to seed curve fitting::
+        .. math::
+            S_0=\\frac{pv-c}{exp(-TE/T_2)}
+
 
     Parameters
     ----------
@@ -524,7 +533,7 @@ def est_t2_s0(te, t2, pv, c=0.0):
     Returns
     -------
     array_like
-        Initial s0 estimate ``s0 = (pv - c) / np.exp(-te/t2)``.
+        Initial s0 estimate.
 
     """
     return (pv - c) / np.exp(-te / t2)
@@ -999,8 +1008,8 @@ def __init__(self, image_slices, template_dcm=None, plate_number=None):
                          dicom_order_key='EchoTime')
 
         self.fit_function = t2_function
-        self.fit_jacobian = t2_jacobian
-        self.fit_eqn_str = 's0 * np.exp(-te / t2)'
+        self.fit_jacobian = None
+        self.fit_eqn_str = 'T2 with Rician noise (Raya et al 2010)'
 
     def initialise_fit_parameters(self, t2_estimates):
         """
@@ -1010,7 +1019,7 @@ def initialise_fit_parameters(self, t2_estimates):
         
         s0 is estimated using est_t2_s0(te, t2_est, mean_pv, c).
             
-        C is estimated as 0.0, the theoretical value assuming Gaussian noise.
+        C is estimated as 5.0.
    
         Parameters
         ----------
@@ -1026,25 +1035,23 @@ def initialise_fit_parameters(self, t2_estimates):
         self.t2_est = t2_estimates
         rois = self.ROI_time_series
         rois_second_mean = np.array([roi.means[1] for roi in rois])
-        self.c_est = np.full_like(self.t2_est, 0.0)
+        self.c_est = np.full_like(self.t2_est, SEED_RICIAN_NOISE)
         # estimate s0 from second image--first image is too low.
         self.s0_est = est_t2_s0(rois[0].times[1], t2_estimates,
                                 rois_second_mean, self.c_est)
-        # Get maximum time to use on fitting algorithm (5*t2_est)
-        # Truncating data after this avoids fitting Rician noise
-        self.max_fit_times = 5 * t2_estimates
 
     def find_relax_times(self):
         """
         Calculate T2 values. Access as ``image_stack.t2s``.
         
-        Uses the 'skip first echo' fit method [1]_. At times >> T2, the signal
-        is dwarfed by Rician noise (for magnitude images). This can lead to
-        inaccuracies in determining T2 as the measured signal does not tend to
-        zero. To counter this, the signal is truncated after
-        ``self.max_fit_times[i]``. At least three signals are used in the fit
-        even if this exceeds the above criteria.
-    
+        Uses the 'skip first echo' fit method [1]_ with a Rician noise model
+        [2]_. Ideally the Rician noise parameter should be determined from the
+        images rather than fitted. However, this is not possible as the noise
+        profile varies across the image due to spatial coil sensitivities and
+        whether the image is normalised or unfiltered. Fitting the noise
+        parameter makes this easier. It has an upper limit of MAX_RICIAN_NOISE,
+        currently set to 20.0.
+
         Returns
         -------
         None.
@@ -1054,25 +1061,30 @@ def find_relax_times(self):
         .. [1] McPhee, K. C., & Wilman, A. H. (2018). Limitations of skipping 
         echoes for exponential T2 fitting. Journal of Magnetic Resonance 
         Imaging, 48(5), 1432-1440. https://doi.org/10.1002/jmri.26052
+
+        .. [2] Raya, J.G., Dietrich, O., Horng, A., Weber, J., Reiser, M.F. and
+        Glaser, C., 2010. T2 measurement in articular cartilage: impact of the
+        fitting method on accuracy and precision at low SNR. Magnetic Resonance
+        in Medicine: An Official Journal of the International Society for
+         Magnetic Resonance in Medicine, 63(1), pp.181-193.
         """
-        # Require at least 4 samples (nb first sample is omitted)
-        min_number_times = 4
+
         rois = self.ROI_time_series
         #  Omit the first image data from the curve fit. This is achieved by
         #  slicing rois[i].times[1:] and rois[i].means[1:]. Skipping odd echoes
         #  can be implemented with rois[i].times[1::2] and .means[1::2]
-        bounds = ([0, 0], [np.inf, np.inf])
+
+        bounds = ([0, 0, 1], [np.inf, np.inf, MAX_RICIAN_NOISE])
+
         self.relax_fit = []
         for i in range(len(rois)):
-            times = [t for t in rois[i].times if t < self.max_fit_times[i]]
-            if len(times) < min_number_times:
-                times = rois[i].times[:min_number_times]
             self.relax_fit.append(
                 scipy.optimize.curve_fit(self.fit_function,
-                                         times[1:],
-                                         rois[i].means[1:len(times)],
+                                         rois[i].times[1:],
+                                         rois[i].means[1:],
                                          p0=[self.t2_est[i],
-                                             self.s0_est[i]],
+                                             self.s0_est[i],
+                                             self.c_est[i]],
                                          jac=self.fit_jacobian,
                                          bounds=bounds,
                                          method='trf'
@@ -1250,7 +1262,6 @@ def main(dcm_target_list, *, plate_number=None,
                           f'pub={relax_published[i]:.4g} '
                           f'({frac_time_diff[i] * 100:+.2f}%)',
                           fontsize=8)
-        # plt.tight_layout(rect=(0,0,0,0.95)) # Leave suptitle space at top
         if report_path:
             # Improve saved image quality
             old_dims = fig.get_size_inches()

tests/test_relaxometry.py

@@ -91,9 +91,9 @@ class TestRelaxometry(unittest.TestCase):
                  131.7, 93.3, 66.6, 45.1, 32.5, 22.4, 2632]
 
     # Values from testing (to check for variations)
-    PLATE4_T2 = [816.633328, 590.521423, 430.902617, 310.985158, 217.258533,
-                 156.093083, 109.839424, 79.269721, 55.998743, 39.044842,
-                 27.143183, 18.448322, 12.514836, 9.351455, 2376]
+    PLATE4_T2 = [ 816.33093, 590.26915, 430.69191, 310.801929, 216.998622,
+                  155.68107, 109.346141, 78.967168, 55.986545, 39.080988,
+                  27.59972,   18.389453, 12.336855, 8.035046, 2374.533555]
 
     TEMPLATE_PATH_T2 = os.path.join(TEST_DATA_DIR, 'relaxometry', 'T2',
                                     'Template_plate4_T2')