Merge pull request #224 from UKRIN-MAPS/rel/v0.7.2

rel/v0.7.2

.github/workflows/python_CI.yml

@@ -49,3 +49,5 @@ jobs:
         env_vars: OS,PYTHON
         fail_ci_if_error: true
         verbose: false
+      env:
+        CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}

CHANGELOG.md

@@ -1,3 +1,17 @@
+## [0.7.2] - 2024-07-05
+
+### Added
+* Export the fit signals from mapping functions e.g. the expected T1 recovery at times TI given the fit values of T1 and M0. #221
+* Warnings if there aren't many negative values in T1 mapping data. #222 #223
+* `get_fit_signal` methods have been added to most mapping classes to allow the user to get the fit signal at a given time point.
+
+### Fixed
+* Issue where resources sub-module was not included on PyPI
+
+### Changed
+* T1 mapping data is now assumed to have been magnitude corrected if the first percentile of the first inversion time is 
+  negative rather than the minimum value. This should make it more robust to noise/preprocessing artefacts. #222 #223
+
 ## [0.7.1] - 2024-02-23
 
 ### Added

CITATION.cff

@@ -17,7 +17,7 @@ authors:
   given-names: "Susan T"
   orcid: "https://orcid.org/0000-0003-0903-7507"
 title: "UKRIN Kidney Analysis Toolbox"
-version: 0.7.1
+version: 0.7.2
 doi: 10.5281/zenodo.4742470
-date-released: 2024-02-23
+date-released: 2024-07-05
 url: "https://github.com/UKRIN-MAPS/ukat"

setup.py

@@ -10,7 +10,7 @@
 
 setup(
     name="ukat",
-    version="0.7.1",
+    version="0.7.2",
     description="UKRIN Kidney Analysis Toolbox",
     long_description=long_description,
     long_description_content_type="text/markdown",

ukat/mapping/diffusion.py

@@ -302,6 +302,27 @@ def to_nifti(self, output_directory=os.getcwd(), base_file_name='Output',
                              'should be "all" or a list of maps from '
                              '"["adc", "adc_err", "mask"]".')
 
+    def get_fit_signal(self):
+        """
+        Get the fit signal from the model used to fit the data i.e. the
+        simulated signal at each b-value given the estimated ADC and S0.
+
+        Returns
+        -------
+        fit_signal : np.ndarray
+            An array containing the fit signal generated by the model
+        """
+        fit_signal = np.zeros((self.n_vox, self.n_bvals))
+        params = np.array([self.adc.reshape(-1),
+                           self.s0.reshape(-1)])
+
+        for n in range(self.n_vox):
+            fit_signal[n, :] = adc_eq(self.u_bvals,
+                                      params[0, n],
+                                      params[1, n])
+        fit_signal = fit_signal.reshape((*self.shape, self.n_bvals))
+        return fit_signal
+
 
 def adc_eq(bvals, adc, s0):
     """

ukat/mapping/t1.py

@@ -52,26 +52,56 @@ def __init__(self, pixel_array, ti, parameters=2, mask=None, tss=0,
         self.tss = tss
         self.tss_axis = tss_axis
 
-        if np.min(pixel_array) < 0:
+        # Assume the data has been magnitude corrected if the first
+        # percentile of the first inversion time is negative.
+        if np.percentile(pixel_array[..., 0], 1) < 0:
             self.mag_corr = True
+            neg_percent = (np.sum(pixel_array[..., 0] < 0)
+                           / pixel_array[..., 0].size)
+            if neg_percent < 0.05:
+                warnings.warn('Fitting data to a magnitude corrected '
+                              'inversion recovery curve however, less than 5% '
+                              'of the data from the first inversion is '
+                              'negative. If you have performed magnitude '
+                              'correction ignore this warning, otherwise the '
+                              'negative values could be due to noise or '
+                              'preprocessing steps  such as EPI distortion '
+                              'correction and  registration.\n'
+                              f'Percentage of first inversion data that is '
+                              f'negative = {neg_percent:.2%}')
         else:
             self.mag_corr = False
+            if np.nanmin(pixel_array) < 0:
+                warnings.warn('Negative values found in data from the first '
+                              'inversion but as the first percentile is not '
+                              'negative, it is assumed these are negative '
+                              'due to noise or preprocessing steps such as '
+                              'EPI distortion correction and registration. '
+                              'As such the data will be fit to the modulus of '
+                              'the recovery curve.\n'
+                              f'Min value = {np.nanmin(pixel_array[..., 0])}\n'
+                              '1st percentile = '
+                              f'{np.percentile(pixel_array[..., 0], 1)}')
 
         if self.parameters == 2:
             if self.mag_corr:
-                super().__init__(pixel_array, ti, two_param_eq, mask,
+                self.t1_eq = two_param_eq
+                super().__init__(pixel_array, ti, self.t1_eq, mask,
                                  multithread)
             else:
-                super().__init__(pixel_array, ti, two_param_abs_eq, mask,
+                self.t1_eq = two_param_abs_eq
+                super().__init__(pixel_array, ti, self.t1_eq, mask,
                                  multithread)
             self.bounds = ([0, 0], [5000, 1000000000])
             self.initial_guess = [1000, 30000]
         elif self.parameters == 3:
             if self.mag_corr:
-                super().__init__(pixel_array, ti, three_param_eq, mask,
+                self.t1_eq = three_param_eq
+                super().__init__(pixel_array, ti, self.t1_eq, mask,
                                  multithread)
             else:
-                super().__init__(pixel_array, ti, three_param_abs_eq, mask,
+                self.t1_eq = three_param_abs_eq
+                super().__init__(pixel_array, ti, self.t1_eq, mask,
                                  multithread)
             self.bounds = ([0, 0, 1], [5000, 1000000000, 2])
             self.initial_guess = [1000, 30000, 2]
@@ -224,10 +254,10 @@ def __init__(self, pixel_array, inversion_list, affine, tss=0, tss_axis=-2,
                               'using parameters=3.')
 
         # Fit Data
-        fitting_model = T1Model(self.pixel_array, self.inversion_list,
-                                self.parameters, self.mask, self.tss,
-                                self.tss_axis, self.multithread)
-        popt, error, r2 = fitting.fit_image(fitting_model)
+        self.fitting_model = T1Model(self.pixel_array, self.inversion_list,
+                                     self.parameters, self.mask, self.tss,
+                                     self.tss_axis, self.multithread)
+        popt, error, r2 = fitting.fit_image(self.fitting_model)
         self.t1_map = popt[0]
         self.m0_map = popt[1]
         self.t1_err = error[0]
@@ -242,8 +272,10 @@ def __init__(self, pixel_array, inversion_list, affine, tss=0, tss_axis=-2,
         # M0 out. Not filtering based on eff as this should ideally be at
         # the upper bound!
         threshold = 0.999  # 99.9% of the upper bound
-        bounds_mask = ((self.t1_map > fitting_model.bounds[1][0] * threshold) |
-                       (self.m0_map > fitting_model.bounds[1][1] * threshold))
+        bounds_mask = ((self.t1_map > self.fitting_model.bounds[1][0] *
+                        threshold) |
+                       (self.m0_map > self.fitting_model.bounds[1][1] *
+                        threshold))
         self.t1_map[bounds_mask] = 0
         self.m0_map[bounds_mask] = 0
         self.t1_err[bounds_mask] = 0
@@ -346,6 +378,44 @@ def to_nifti(self, output_directory=os.getcwd(), base_file_name='Output',
 
         return
 
+    def get_fit_signal(self):
+        """
+        Get the fit signal from the model used to fit the data i.e. the
+        simulated signal at each inversion time given the estimated T1, M0
+        (and inversion efficiency if applicable).
+
+        Returns
+        -------
+        fit_signal : np.ndarray
+            An array containing the fit signal generated by the model
+        """
+        if self.molli:
+            t1 = self.t1_map / ((self.m0_map * self.eff_map) / self.m0_map - 1)
+        else:
+            t1 = self.t1_map
+        t1_lin = t1.reshape(-1)
+        m0_lin = self.m0_map.reshape(-1)
+        if self.parameters == 3:
+            eff_lin = self.eff_map.reshape(-1)
+
+        fit_signal = np.zeros((self.n_vox, self.n_ti))
+
+        if self.parameters == 2:
+            for n, (ti, t1, m0) in enumerate(zip(self.fitting_model.x_list,
+                                                 t1_lin,
+                                                 m0_lin)):
+                fit_signal[n, :] = self.fitting_model.t1_eq(ti, t1, m0)
+        else:
+            for n, (ti, t1, m0, eff) in (
+                enumerate(zip(self.fitting_model.x_list,
+                              t1_lin,
+                              m0_lin,
+                              eff_lin))):
+                fit_signal[n, :] = self.fitting_model.t1_eq(ti, t1, m0, eff)
+
+        fit_signal = fit_signal.reshape((*self.shape, self.n_ti))
+        return fit_signal
+
 
 def two_param_abs_eq(t, t1, m0):
     """

ukat/mapping/t2.py

@@ -39,11 +39,13 @@ def __init__(self, pixel_array, te, method='2p_exp', mask=None,
         self.method = method
 
         if self.method == '2p_exp':
-            super().__init__(pixel_array, te, two_param_eq, mask, multithread)
+            self.t2_eq = two_param_eq
+            super().__init__(pixel_array, te, self.t2_eq, mask, multithread)
             self.bounds = ([0, 0], [1000, 100000000])
             self.initial_guess = [20, 10000]
         elif self.method == '3p_exp':
-            super().__init__(pixel_array, te, three_param_eq, mask,
+            self.t2_eq = three_param_eq
+            super().__init__(pixel_array, te, self.t2_eq, mask,
                              multithread)
             self.bounds = ([0, 0, 0], [1000, 100000000, 1000000])
             self.initial_guess = [20, 10000, 500]
@@ -176,12 +178,12 @@ def __init__(self, pixel_array, echo_list, affine, mask=None,
         self.multithread = multithread
 
         # Fit data
-        fitting_model = T2Model(self.pixel_array, self.echo_list,
+        self.fitting_model = T2Model(self.pixel_array, self.echo_list,
                                 self.method, self.mask, self.multithread)
 
         if self.noise_threshold > 0:
-            fitting_model.threshold_noise(self.noise_threshold)
-        popt, error, r2 = fitting.fit_image(fitting_model)
+            self.fitting_model.threshold_noise(self.noise_threshold)
+        popt, error, r2 = fitting.fit_image(self.fitting_model)
         self.t2_map = popt[0]
         self.m0_map = popt[1]
         self.t2_err = error[0]
@@ -195,8 +197,10 @@ def __init__(self, pixel_array, echo_list, affine, mask=None,
         # Filter values that are very close to models upper bounds of T2 or
         # M0 out.
         threshold = 0.999  # 99.9% of the upper bound
-        bounds_mask = ((self.t2_map > fitting_model.bounds[1][0] * threshold) |
-                       (self.m0_map > fitting_model.bounds[1][1] * threshold))
+        bounds_mask = ((self.t2_map > self.fitting_model.bounds[1][0] *
+                        threshold) |
+                       (self.m0_map > self.fitting_model.bounds[1][1] *
+                        threshold))
         self.t2_map[bounds_mask] = 0
         self.m0_map[bounds_mask] = 0
         self.t2_err[bounds_mask] = 0
@@ -268,6 +272,34 @@ def to_nifti(self, output_directory=os.getcwd(), base_file_name='Output',
                              '"mask"]".')
         return
 
+    def get_fit_signal(self):
+        """
+        Get the fit signal from the model used to fit the data i.e. the
+        simulated signal at each echo time given the estimated T2, M0
+        (and baseline noise floor (b) if applicable).
+
+        Returns
+        -------
+        fit_signal : np.ndarray
+            An array containing the fit signal generated by the model
+        """
+        fit_signal = np.zeros((self.n_vox, self.n_te))
+        if self.method == '2p_exp':
+            params = np.array([self.t2_map.reshape(-1),
+                               self.m0_map.reshape(-1)])
+        elif self.method == '3p_exp':
+            params = np.array([self.t2_map.reshape(-1),
+                               self.m0_map.reshape(-1),
+                               self.b_map.reshape(-1)])
+
+        for n in range(self.n_vox):
+            fit_signal[n] = self.fitting_model.t2_eq(self.echo_list,
+                                                     *params[:, n])
+
+        fit_signal = fit_signal.reshape((*self.shape, self.n_te))
+        return fit_signal
+
+
 
 def two_param_eq(t, t2, m0):
     """

ukat/mapping/t2_stimfit.py

@@ -442,6 +442,28 @@ def to_nifti(self, output_directory=os.getcwd(), base_file_name='Output',
                              '"["t2", "m0", "b1", "r2", '
                              '"mask"]".')
 
+    def get_fit_signal(self):
+        """
+        Get the fit signal from the model used to fit the data i.e. the
+        simulated signal at each echo time given the estimated T2, M0. Note
+        if norm=True, the fit signal will also be normalised.
+
+        Returns
+        -------
+        fit_signal : np.ndarray
+            An array containing the fit signal generated by the model
+        """
+        fit_signal = np.zeros((self.n_vox, self.model.opt['etl']))
+        params = np.array([self.t2_map.reshape(-1),
+                           self.m0_map.reshape(-1)])
+
+        for n in range(self.n_vox):
+            fit_signal[n, :] = two_param_eq(self.model.opt['te'] * 1000,
+                                            params[0, n],
+                                            params[1, n])
+        fit_signal = fit_signal.reshape((*self.shape, self.model.opt['etl']))
+        return fit_signal
+
     def _fit(self):
         mask = self.mask.flatten()
         signal = self.pixel_array.reshape(self.n_vox, self.model.opt['etl'])

ukat/mapping/t2star.py

@@ -358,6 +358,27 @@ def to_nifti(self, output_directory=os.getcwd(), base_file_name='Output',
 
         return
 
+    def get_fit_signal(self):
+        """
+        Get the fit signal from the model used to fit the data i.e. the
+        simulated signal at each echo time given the estimated T2* and M0.
+
+        Returns
+        -------
+        fit_signal : np.ndarray
+            An array containing the fit signal generated by the model
+        """
+        fit_signal = np.zeros((self.n_vox, self.n_te))
+        params = np.array([self.t2star_map.reshape(-1),
+                           self.m0_map.reshape(-1)])
+
+        for n in range(self.n_vox):
+            fit_signal[n, :] = two_param_eq(self.echo_list,
+                                            params[0, n],
+                                            params[1, n])
+        fit_signal = fit_signal.reshape((*self.shape, self.n_te))
+        return fit_signal
+
 
 def two_param_eq(t, t2star, m0):
     """

ukat/mapping/tests/test_b0.py

@@ -71,9 +71,6 @@ def test_mask(self):
         masked_pixels = B0(self.correct_array, self.correct_echo_list,
                            self.affine, mask=mask)
 
-        assert (all_pixels.phase_difference !=
-                masked_pixels.phase_difference).any()
-        assert (all_pixels.b0_map != masked_pixels.b0_map).any()
         assert (arraystats.ArrayStats(all_pixels.b0_map).calculate() !=
                 arraystats.ArrayStats(masked_pixels.b0_map).calculate())
 

ukat/mapping/tests/test_diffusion.py

@@ -197,6 +197,16 @@ def test_to_nifti(self):
         # Delete 'test_output' folder
         shutil.rmtree('test_output')
 
+    def test_get_fit_signal(self):
+        mapper = ADC(self.pixel_array, self.affine, self.bvals, self.mask)
+        fit_signal = mapper.get_fit_signal()
+        stats = arraystats.ArrayStats(fit_signal).calculate()
+        npt.assert_allclose([stats["mean"]["4D"], stats["std"]["4D"],
+                             stats["min"]["4D"], stats["max"]["4D"]],
+                            [33971.334576954156, 31856.26958366113,
+                             0.0, 241017.96908169912],
+                            rtol=1e-6, atol=1e-4)
+
 
 class TestDTI:
     pixel_array, affine, bvals, bvecs = fetch.dwi_philips()