Mdr t1 molli

.github/workflows/python_CI.yml

@@ -21,14 +21,15 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, windows-latest, macos-latest]
-        python-version: [3.8, 3.9, "3.10", "3.11"]
+        python-version: [3.9, "3.10", "3.11"]
 
     steps:
     - uses: actions/checkout@v2
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v2
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
+        cache: 'pip'
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip

.gitignore

@@ -1,3 +1,6 @@
+# Virtual environments
+.venv/
+
 # Visual Studio Code editor
 .vscode/
 

requirements.txt

@@ -12,3 +12,7 @@ scipy>=1.9.1, <1.12
 scikit-learn~=1.2.1  # Changing versions is known to cause changes in iSNR results
 tabulate>=0.9.0
 tqdm>=4.64.1
+mdreg>=0.4.1
+itk>=5.4.0; sys_platform == "linux"
+itk>=5.4.0; sys_platform == "darwin"
+itk<5.4.0; sys_platform == "win32"

setup.py

@@ -17,7 +17,7 @@
     url="https://github.com/UKRIN-MAPS/ukat",
     license="GPL-3.0",
 
-    python_requires='>=3.8, <4',
+    python_requires='>=3.9, <3.12',
     packages=find_packages(),
     install_requires=requirements,
     include_package_data=True,
@@ -29,7 +29,6 @@
         'Topic :: Scientific/Engineering',
         'Environment :: Console',
         'Operating System :: OS Independent',
-        'Programming Language :: Python :: 3.8',
         'Programming Language :: Python :: 3.9',
         'Programming Language :: Python :: 3.10',
         'Programming Language :: Python :: 3.11',

ukat/mapping/t1.py

@@ -3,6 +3,8 @@
 import os
 import warnings
 
+import mdreg
+
 from . import fitting
 
 
@@ -130,6 +132,9 @@
     r2 : np.ndarray
         The R-Squared value of the fit, values close to 1 indicate a good
         fit, lower values indicate a poorer fit
+    deformation_field : np.ndarray
+        The deformation field generated by the model-driven registration
+        process.
     shape : tuple
         The shape of the T1 map
     n_ti : int
@@ -141,7 +146,7 @@
 
     def __init__(self, pixel_array, inversion_list, affine, tss=0, tss_axis=-2,
                  mask=None, parameters=2, mag_corr=False, molli=False,
-                 multithread=True):
+                 multithread=True, mdr=False):
         """Initialise a T1 class instance.
 
         Parameters
@@ -200,15 +205,11 @@
             increase in speed distributing the calculation would generate.
             'auto' attempts to apply multithreading where appropriate based
             on the number of voxels being fit.
+        mdr : bool, optional
+            Default 'False`
+            If True, this performs a motion correction with model-driven
+            registration before performing the final fit to the model function.
         """
-        assert multithread in [True,
-                               False,
-                               'auto'], (f'multithreaded must '
-                                         f'be True, False or auto. You '
-                                         f'entered { multithread}.')
-        assert mag_corr in [True,
-                            False], (f'mag_corr must be True or False. '
-                                     f'You entered {mag_corr}.')
         # Normalise the data so its roughly in the same range across vendors
         self.scale = np.nanmax(pixel_array)
         self.pixel_array = pixel_array / self.scale
@@ -222,6 +223,9 @@
         if mask is None:
             self.mask = np.ones(self.shape, dtype=bool)
         else:
+            if mdr is True:
+                raise ValueError('Masking is not supported when using '
+                                 'model-driven registration.')
             self.mask = mask.astype(bool)
         # Don't process any nan values
         self.mask[np.isnan(np.sum(pixel_array, axis=-1))] = False
@@ -241,8 +245,20 @@
             else:
                 multithread = False
         self.multithread = multithread
+        self.mdr = mdr
 
         # Some sanity checks
+        assert multithread in [True,
+                               False,
+                               'auto'], (f'multithreaded must '
+                                         f'be True, False or auto. You '
+                                         f'entered {multithread}.')
+        assert mag_corr in [True,
+                            False], (f'mag_corr must be True or False. '
+                                     f'You entered {mag_corr}.')
+
+        assert mdr in [True, False], (f'mdr must be True or False. '
+                                      f'You entered {mdr}.')
         assert (pixel_array.shape[-1]
                 == len(inversion_list)), 'Number of inversions does not ' \
                                          'match the number of time frames ' \
@@ -252,12 +268,92 @@
                 'Temporal slice spacing can\'t be applied to the TI axis.'
             assert (tss_axis < self.dimensions), \
                 'tss_axis must be less than the number of spatial dimensions'
+            if (self.tss_axis != 2) & (self.mdr is True):
+                raise ValueError('Temporal slice spacing only supported '
+                                 'along the z direction when using '
+                                 'model-driven registration.')
+
         if self.molli:
             if self.parameters == 2:
                 self.parameters = 3
                 warnings.warn('MOLLI requires a three parameter fit, '
                               'using parameters=3.')
 
+        if mdr:
+            # The mdreg package can handle register each slice of an M2D
+            # image separately provided the inversion times of each slice
+            # are the same.
+            if self.tss == 0:
+                pixel_array, deform, _, _ = mdreg.fit(
+                    self.pixel_array,
+                    force_2d=True,
+                    verbose=1,
+                    fit_image={
+                        'func': _t1_fit,
+                        'inversion_list': self.inversion_list,
+                        'affine': self.affine,
+                        'tss': self.tss,
+                        'tss_axis': self.tss_axis,
+                        'mask': self.mask,
+                        'parameters': self.parameters,
+                        'mag_corr': self.mag_corr,
+                        # MOLLI-correction is not relevant for MDR
+                        'molli': False,
+                        'multithread': self.multithread,
+                    },
+                    # All default settings but kept here as a template for if
+                    # we decide to expose coreg options to ukat users in the
+                    # future.
+                    fit_coreg={
+                        'package': 'elastix',
+                        'parallel': False,  # elastix is not parallelizable
+                    }
+                )
+            else:
+                pixel_array = np.zeros(self.pixel_array.shape)
+                deform = np.zeros((*self.pixel_array.shape[:3], 2,
+                                   self.pixel_array.shape[3]))
+
+                # The following for loop is a workaround to allow a
+                # different inversion list for each slice of data. Most of
+                # the code comes from the mdreg.fit() function.
+                for slice in range(self.shape[-1]):
+                    print('-----------------')
+                    print('Fitting slice ' + str(slice).zfill(3))
+                    print('-----------------')
+                    inversion_list = (np.array(self.inversion_list)
+                                      + self.tss * slice)
+                    (pixel_array[..., slice, :], deform[..., slice, :, :], _,
+                     _) = mdreg.fit(
+                        self.pixel_array[..., slice, :],
+                        force_2d=True,
+                        verbose=1,
+                        fit_image={
+                            'func': _t1_fit,
+                            'inversion_list': inversion_list,
+                            'affine': self.affine,
+                            'tss': 0,
+                            'tss_axis': None,
+                            'mask': self.mask[..., slice],
+                            'parameters': self.parameters,
+                            'mag_corr': self.mag_corr,
+                            # MOLLI-correction is not relevant for MDR
+                            'molli': False,
+                            'multithread': self.multithread,
+                        },
+                        # All default settings but kept here as a template for
+                        # if we decide to expose coreg options to ukat users
+                        # in the future.
+                        fit_coreg={
+                            'package': 'elastix',
+                            'parallel': False,  # elastix is not parallelizable
+                        }
+                    )
+            # Changing the dimensions of the deformation field to a more
+            # intuitive order.
+            self.deformation_field = np.swapaxes(deform, -2, -1)
+            self.pixel_array = pixel_array
+
         # Fit Data
         self.fitting_model = T1Model(self.pixel_array, self.inversion_list,
                                      self.parameters, self.mask, self.tss,
@@ -338,13 +434,13 @@
         maps : list or 'all', optional
             List of maps to save to NIFTI. This should either the string "all"
             or a list of maps from ["t1", "t1_err", "m0", "m0_err", "eff",
-            "eff_err", "r1", "r2", "mask"]
+            "eff_err", "deformation_field", "r1", "r2", "mask"]
         """
         os.makedirs(output_directory, exist_ok=True)
         base_path = os.path.join(output_directory, base_file_name)
         if maps == 'all' or maps == ['all']:
-            maps = ['t1', 't1_err', 'm0', 'm0_err', 'eff', 'eff_err', 'r1_map',
-                    'r2', 'mask']
+            maps = ['t1', 't1_err', 'm0', 'm0_err', 'eff', 'eff_err',
+                    'deformation_field', 'r1_map', 'r2', 'mask']
         if isinstance(maps, list):
             for result in maps:
                 if result == 't1' or result == 't1_map':
@@ -370,6 +466,11 @@
                     eff_err_nifti = nib.Nifti1Image(self.eff_err,
                                                     affine=self.affine)
                     nib.save(eff_err_nifti, base_path + '_eff_err.nii.gz')
+                elif self.mdr is True and result == 'deformation_field':
+                    deformation_nifti = nib.Nifti1Image(self.deformation_field,
+                                                        affine=self.affine)
+                    nib.save(deformation_nifti,
+                             base_path + '_deformation_field.nii.gz')
                 elif result == 'r1' or result == 'r1_map':
                     r1_nifti = nib.Nifti1Image(T1.r1_map(self),
                                                affine=self.affine)
@@ -386,7 +487,7 @@
             raise ValueError('No NIFTI file saved. The variable "maps" '
                              'should be "all" or a list of maps from '
                              '"["t1", "t1_err", "m0", "m0_err", "eff", '
-                             '"eff_err", "r1", "mask"]".')
+                             '"eff_err", "deformation_field", "r1", "mask"]".')
 
         return
 
@@ -428,6 +529,19 @@
         fit_signal = fit_signal.reshape((*self.shape, self.n_ti))
         return fit_signal
 
+    def get_pixel_array(self):
+        """
+        Get the pixel array from the T1 class. This method should be used
+        rather than T1.pixel_array as it will return the data in the
+        original scale.
+
+        Returns
+        -------
+        pixel_array : np.ndarray
+            An array containing the pixel data in the original scale.
+        """
+        return self.pixel_array * self.scale
+
 
 def two_param_abs_eq(t, t1, m0):
     """
@@ -561,10 +675,24 @@
     for ti in range(pixel_array.shape[-1]):
         pixel_array_prime[..., ti] = (pixel_array[..., ti] *
                                       pixel_array[..., -1].conjugate()) \
-                                     / np.abs(pixel_array[..., -1])
+            / np.abs(pixel_array[..., -1])
 
     phase_factor = np.imag(np.log(pixel_array_prime / np.abs(pixel_array)))
     phase_offset = np.abs(phase_factor) - (np.pi / 2)
     sign = -(phase_offset / np.abs(phase_offset))
     corrected_array = sign * np.abs(pixel_array)
+    corrected_array = np.nan_to_num(corrected_array)
     return corrected_array
+
+
+# Private wrapper for use by mdreg
+def _t1_fit(pixel_array, inversion_list=None, affine=None, **kwargs):
+    """
+    Private wrapper for use by mdreg
+    """
+    map = T1(pixel_array, inversion_list, affine, **kwargs)
+    if map.parameters == 2:
+        pars = np.stack((map.t1_map, map.m0_map), axis=-1)
+    else:
+        pars = np.stack((map.t1_map, map.m0_map, map.eff_map), axis=-1)
+    return map.get_fit_signal(), pars

ukat/mapping/tests/test_t1.py

@@ -301,14 +301,28 @@ def test_mag_corr_options(self):
                         mag_corr='yes please',
                         multithread=False)
 
-
     def test_molli_2p_warning(self):
         signal_array = np.tile(self.correct_signal_three_param, (10, 10, 3, 1))
         with pytest.warns(UserWarning):
             mapper = T1(pixel_array=signal_array,
                         inversion_list=self.t,
                         affine=self.affine, parameters=2, molli=True)
 
+    def test_tss_mdr_error(self):
+        signal_array = np.tile(self.correct_signal_three_param, (10, 10, 3, 1))
+        with pytest.raises(ValueError):
+            mapper = T1(pixel_array=signal_array,
+                        inversion_list=self.t,
+                        affine=self.affine, tss=10, tss_axis=1, mdr=True)
+
+    def test_mask_mdr_error(self):
+        signal_array = np.tile(self.correct_signal_three_param, (10, 10, 3, 1))
+        mask = np.ones(signal_array.shape[:-1])
+        with pytest.raises(ValueError):
+            mapper = T1(pixel_array=signal_array,
+                        inversion_list=self.t,
+                        affine=self.affine, mdr=True, mask=mask)
+
     def test_real_data(self):
         # Get test data
         magnitude, phase, affine, ti, tss = fetch.t1_philips(2)
@@ -329,6 +343,10 @@ def test_real_data(self):
         gold_standard_3p_single = [1347.824169, 657.254769, 0.0, 3948.24018]
         gold_standard_molli = [1554.586501,  606.863022, -170.611303,
                                6025.763663]
+        gold_standard_molli_mdr = [1528.876958, 659.720555, -209.721654,
+                                   5707.662715]
+        gold_standard_2p_mdr = [1038.024629, 427.33669, 223.047457,
+                                2600.325215]
 
         # Two parameter method
         mapper = T1(magnitude, ti, affine, parameters=2, tss=tss)
@@ -360,6 +378,35 @@ def test_real_data(self):
                              t1_stats['min']['3D'], t1_stats['max']['3D']],
                             gold_standard_molli, rtol=1e-6, atol=5e-3)
 
+        # MDR TSS == 0
+        mapper = T1(image_molli[:, :, :2, :], ti_molli, affine_molli,
+                    parameters=3, molli=True, mdr=True)
+        t1_stats = arraystats.ArrayStats(mapper.t1_map).calculate()
+        # Large tolerance as ITK performs differently on MacOS, Linux and
+        # Windows
+        npt.assert_allclose([t1_stats['mean']['3D'], t1_stats['std']['3D'],
+                             t1_stats['min']['3D'], t1_stats['max']['3D']],
+                            gold_standard_molli_mdr, rtol=0.1, atol=50)
+
+        # MDR TSS != 0
+        mapper = T1(magnitude[:, :, :2, :], ti, affine,
+                    parameters=2, tss=tss, mdr=True)
+        t1_stats = arraystats.ArrayStats(mapper.t1_map).calculate()
+        # Large tolerance as ITK performs differently on MacOS, Linux and
+        # Windows
+        npt.assert_allclose([t1_stats['mean']['3D'], t1_stats['std']['3D'],
+                             t1_stats['min']['3D'], t1_stats['max']['3D']],
+                            gold_standard_2p_mdr, rtol=0.1, atol=50)
+
+    def test_get_pixel_array(self):
+        # Create a T1 map instance and test different export to NIFTI scenarios
+        signal_array = np.tile(self.correct_signal_two_param, (10, 10, 3, 1))
+        mapper = T1(signal_array, self.t, self.affine, parameters=2)
+
+        # Check that the pixel array is returned
+        pixel_array = mapper.get_pixel_array()
+        npt.assert_array_almost_equal(pixel_array, signal_array)
+
     def test_to_nifti(self):
         # Create a T1 map instance and test different export to NIFTI scenarios
         signal_array = np.tile(self.correct_signal_three_param, (10, 10, 3, 1))