Mdr t1 molli

.gitignore

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

requirements.txt

@@ -12,3 +12,4 @@ 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

ukat/mapping/t1.py

@@ -3,6 +3,8 @@
 import os
 import warnings
 
+import mdreg
+
 from . import fitting
 
 
@@ -148,8 +150,11 @@
         apart from TI
     """
 
+    # @Alex: suggestion: make affine a keyword parameter with default np.eye(4).
+    # As it is the last argument in the list this will not break existing code
+    # And it means the user is not forced to provide a dummy affine when it plays no role.
     def __init__(self, pixel_array, inversion_list, affine, tss=0, tss_axis=-2,
-                 mask=None, parameters=2, molli=False, multithread=True):
+                 mask=None, parameters=2, molli=False, multithread=True, mdr=False):
         """Initialise a T1 class instance.
 
         Parameters
@@ -201,12 +206,70 @@
             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 is True \
-               or multithread is False \
-               or multithread == 'auto', f'multithreaded must be True,' \
-                                         f'False or auto. You entered ' \
-                                         f'{multithread}'
+            or multithread is False \
+            or multithread == 'auto', f'multithreaded must be True,' \
+            f'False or auto. You entered ' \
+            f'{multithread}'
+
+        # @Alex: I have moved this up so multithreading
+        # settings can be reused in mdreg, which requires a True or False
+        # value. In this case (elastix) it is actually unnecessary as
+        # parallelization is not a real option anyway. But it would be relevant if
+        # wanted to run with skimage, for instamce.
+        if multithread == 'auto':
+            npixels = np.prod(pixel_array.shape[:-1])
+            if npixels > 20:
+                multithread = True
+            else:
+                multithread = False
+
+        if mdr:
+            pixel_array, deform, _, _ = mdreg.fit(
+                pixel_array,
+                fit_image={
+                    'func': _T1_fit,
+                    'inversion_list': inversion_list,
+                    'affine': affine,
+                    'tss': tss,
+                    'tss_axis': tss_axis,
+                    'mask': mask,
+                    'parameters': parameters,
+                    'molli': molli,
+                    'multithread': multithread,
+                },
+                # @Alex: These coreg settings are default so technically speaking do not have
+                # to be specified here. I am leaving it in nevertheless as a template in case we want
+                # to explore alternative coregistration settings later. The fit_coreg
+                # dictionary could also be exposed as a keyword argument in T1.__init__()
+                # in case we want to give the user the option to modify the detail.
+                # As it stands the user has no way of modifying the way mdr runs, eg.
+                # change verbosity level, stopping criteria or coreg options.
+                # I didn't change it as that may exactly be the intention of ukat?
+                fit_coreg={
+                    'package': 'elastix',
+                    'parallel': False,  # elastix is not parallelizable
+                }
+            )
+            # @Alex: At the moment the order of the dimensions in the deformation field returned by mdreg is awkward.
+            # Current dimensions are (x,y,2,t) for 2-dimensional pixel_arrays,
+            # and (x,y,z,3,t) for 3 dimensional. Better would be (x,y,t,2) and (x,y,z,t,3)
+            # - i.e add a new dimension at the end.
+            # We will probable change this in mdreg at some point so I put in an ad-hoc
+            # reordering at this stage, We can just take it out later if mdreg is updated.
+            self.deformation_field = np.swapaxes(deform, -2, -1)
+            # @Alex: Hack to avoid magnitude corrected model being selected
+            # This needs a better solution as this prevents a magnitude corrected mode
+            pixel_array = np.abs(pixel_array)
+        else:
+            # @Alex: is this the expected default?
+            self.deformation_field = None
 
         self.pixel_array = pixel_array
         self.shape = pixel_array.shape[:-1]
@@ -230,11 +293,6 @@
             self.tss = 0
         self.parameters = parameters
         self.molli = molli
-        if multithread == 'auto':
-            if self.n_vox > 20:
-                multithread = True
-            else:
-                multithread = False
         self.multithread = multithread
 
         # Some sanity checks
@@ -549,10 +607,18 @@
     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)
     return corrected_array
+
+
+# Private wrapper for use by mdreg
+def _T1_fit(pixel_array, inversion_list=None, affine=None, **kwargs):
+    # Alex: added the abs() here to avoid the bug with model selection
+    # This needs a better solution as this prevents a magnitude corrected mode
+    map = T1(np.abs(pixel_array), inversion_list, affine, **kwargs)
+    return map.get_fit_signal(), None

ukat/mapping/tests/test_mdreg.py

@@ -0,0 +1,46 @@
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import mdreg
+
+from ukat.mapping.t1 import T1
+
+
+def t1_maps():
+
+    # fetch data
+    data = mdreg.fetch('MOLLI')
+    array = data['array'][:, :, 0, :]
+    TI = np.array(data['TI'])
+
+    # Calculate corrected uncorrected T1-map
+    map = T1(array, TI, np.eye(4))
+    np.save(os.path.join(os.getcwd(), 't1_uncorr'), map.t1_map)
+
+    # Calculate corrected corrected T1-map
+    map = T1(array, TI, np.eye(4), mdr=True)
+    np.save(os.path.join(os.getcwd(), 't1_corr'), map.t1_map)
+
+
+def plot_t1_maps():
+
+    # Plot corrected and uncorrected side by side
+    t1_uncorr = np.load(os.path.join(os.getcwd(), 't1_uncorr.npy'))
+    t1_corr = np.load(os.path.join(os.getcwd(), 't1_corr.npy'))
+    fig, ax = plt.subplots(figsize=(10, 6), ncols=2, nrows=1)
+    for i in range(2):
+        ax[i].set_yticklabels([])
+        ax[i].set_xticklabels([])
+    ax[0].set_title('T1-map without MDR')
+    ax[1].set_title('T1-map with MDR')
+    im = ax[0].imshow(t1_uncorr.T, cmap='gray', vmin=0, vmax=2000)
+    im = ax[1].imshow(t1_corr.T, cmap='gray', vmin=0, vmax=2000)
+    fig.colorbar(im, ax=ax.ravel().tolist())
+    plt.savefig(os.path.join(os.getcwd(), 't1_mdr'))
+    plt.show()
+
+
+if __name__ == '__main__':
+
+    # t1_maps()
+    plot_t1_maps()