From 62d5e74212d6588b84c74b0c64e9200e4bd54034 Mon Sep 17 00:00:00 2001 From: plaresmedima Date: Tue, 1 Oct 2024 22:44:36 +0100 Subject: [PATCH 1/5] Added .venv to gitignore --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.gitignore b/.gitignore index 64659682..3dbcc034 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,6 @@ +# Virtual environments +.venv/ + # Visual Studio Code editor .vscode/ From 8ea0c184812991e252b621530c3821ae71c46523 Mon Sep 17 00:00:00 2001 From: plaresmedima Date: Tue, 1 Oct 2024 22:45:06 +0100 Subject: [PATCH 2/5] Add mdreg to requirements --- requirements.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/requirements.txt b/requirements.txt index a924e15d..16557d3d 100644 --- a/requirements.txt +++ b/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 From e8a27ed0d184dc3b8cbaff752fa8cf9351856144 Mon Sep 17 00:00:00 2001 From: plaresmedima Date: Tue, 1 Oct 2024 22:45:46 +0100 Subject: [PATCH 3/5] Add mdr keyword to T1() class --- ukat/mapping/t1.py | 79 ++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 73 insertions(+), 6 deletions(-) diff --git a/ukat/mapping/t1.py b/ukat/mapping/t1.py index 0fb63f47..0bf018c7 100644 --- a/ukat/mapping/t1.py +++ b/ukat/mapping/t1.py @@ -3,6 +3,8 @@ import os import warnings +import mdreg + from . import fitting @@ -148,8 +150,11 @@ class T1: 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,71 @@ def __init__(self, pixel_array, inversion_list, affine, tss=0, tss_axis=-2, 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}' + + # @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 +294,6 @@ def __init__(self, pixel_array, inversion_list, affine, tss=0, tss_axis=-2, 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 @@ -556,3 +615,11 @@ def magnitude_correct(pixel_array): 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 From dbfd7bd4fef40556974cc000d5fc24d8c1a30d21 Mon Sep 17 00:00:00 2001 From: plaresmedima Date: Tue, 1 Oct 2024 22:46:44 +0100 Subject: [PATCH 4/5] Added script to run mdreg examples Note this is not a proper test - this file probably belongs elsewhere - maybe can become the basis for a tutorial on motion correction. --- ukat/mapping/tests/test_mdreg.py | 50 ++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 ukat/mapping/tests/test_mdreg.py diff --git a/ukat/mapping/tests/test_mdreg.py b/ukat/mapping/tests/test_mdreg.py new file mode 100644 index 00000000..db672a1e --- /dev/null +++ b/ukat/mapping/tests/test_mdreg.py @@ -0,0 +1,50 @@ +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() + + + + From 1b392c026673984467867942b1a90bf6f99367d6 Mon Sep 17 00:00:00 2001 From: plaresmedima Date: Tue, 1 Oct 2024 23:31:04 +0100 Subject: [PATCH 5/5] autocorrect PEP8 violations --- ukat/mapping/t1.py | 57 ++++++++++++++++---------------- ukat/mapping/tests/test_mdreg.py | 8 ++--- 2 files changed, 30 insertions(+), 35 deletions(-) diff --git a/ukat/mapping/t1.py b/ukat/mapping/t1.py index 0bf018c7..5b2f01d7 100644 --- a/ukat/mapping/t1.py +++ b/ukat/mapping/t1.py @@ -150,7 +150,7 @@ class T1: apart from TI """ - #@Alex: suggestion: make affine a keyword parameter with default np.eye(4). + # @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, @@ -213,15 +213,15 @@ def __init__(self, pixel_array, inversion_list, affine, tss=0, tss_axis=-2, """ assert multithread is True \ - 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 + 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 + # 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]) @@ -229,37 +229,37 @@ def __init__(self, pixel_array, inversion_list, affine, tss=0, tss_axis=-2, multithread = True else: multithread = False - + if mdr: pixel_array, deform, _, _ = mdreg.fit( - pixel_array, - fit_image = { + pixel_array, + fit_image={ 'func': _T1_fit, 'inversion_list': inversion_list, 'affine': affine, - 'tss': tss, + 'tss': tss, 'tss_axis': tss_axis, - 'mask': mask, - 'parameters': parameters, - 'molli': molli, - 'multithread': multithread, + 'mask': mask, + 'parameters': parameters, + 'molli': molli, + 'multithread': multithread, }, - # @Alex: These coreg settings are default so technically speaking do not have + # @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 = { + # 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 + '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) + # 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. @@ -271,7 +271,6 @@ def __init__(self, pixel_array, inversion_list, affine, tss=0, tss_axis=-2, # @Alex: is this the expected default? self.deformation_field = None - self.pixel_array = pixel_array self.shape = pixel_array.shape[:-1] self.dimensions = len(pixel_array.shape) @@ -608,7 +607,7 @@ def magnitude_correct(pixel_array): 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) diff --git a/ukat/mapping/tests/test_mdreg.py b/ukat/mapping/tests/test_mdreg.py index db672a1e..39e226bd 100644 --- a/ukat/mapping/tests/test_mdreg.py +++ b/ukat/mapping/tests/test_mdreg.py @@ -10,7 +10,7 @@ def t1_maps(): # fetch data data = mdreg.fetch('MOLLI') - array = data['array'][:,:,0,:] + array = data['array'][:, :, 0, :] TI = np.array(data['TI']) # Calculate corrected uncorrected T1-map @@ -34,7 +34,7 @@ def plot_t1_maps(): 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) + 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() @@ -44,7 +44,3 @@ def plot_t1_maps(): # t1_maps() plot_t1_maps() - - - -