sct_fmri_moco: Generalize motion correction for sagittal acquisitions and other improvements #2022

jcohenadad · 2018-09-10T19:29:14Z

Currently, sct_fmri_moco is designed such that the input image is assumed to be in axial orientation (i.e. 3rd dimension is S-I), and the algorithm estimates a set of Tx,Ty for each slice, regularized using polynom functions along the S-I direction.

However, for images acquired sagitally (or coronally), this approach does not work anymore. Beyond the simple orientation issue described in #1983, the type of motion-related corruption is not addressed by the current algorithm. More specifically, if motion occurs between sagittal slices, then the cord is non-linearly distorted in the axial plane and a simple translation within the axial plane is not sufficient. This is shown in the animation below. Without motion correction (left), a few slices are translated in the AP direction, leaving ugly distortion of the cord in the axial plane. With motion correction (right), the algorithm tries to find an average translation, but this is of course not sufficient for this type of motion.

The implemented solution is to split the data along Z (R-L dimension) and use antsRegistration in 2D mode for each slice independently. For now we use affine transfo, which is fast and robust, but later we could investigate the use of other transfo. What helps in particular is to use a mask to focus on a region of interest for registration (esp. because motion is non-affine). See below a comparison without (left) and with the new moco algorithm (right):

This PR also provides a few other fixes, refactoring and improvements (runs faster).

Fixes #1984, Fixes #820, Fixes #2017, Fixes #1878

zougloub

Just a question about concurrent branches

zougloub · 2018-09-12T02:50:07Z

scripts/sct_create_mask.py

@@ -110,12 +101,11 @@ def create_mask(param):

    path_tmp = sct.tmp_create(basename="create_mask", verbose=param.verbose)

-    sct.printv('\nCheck orientation...', param.verbose)
+    sct.printv('\nOrientation:', param.verbose)


Hmmm, which PR do you want to merge first?

i knew there would be a problem there... i started #2022, then found a bug which is solved by #2021, but then i needed the changes from #2021 to be able to continue working on #2022, so i did a rebase 2021-->2022 and here we are.

so i guess we should merge #2021 first

This fixes issue #2017

Oherwise 4D data with singleton will become 3D data

Now, if the splitting dimension corresponds to the last dim of the image, then it reshapes to remove the last dim (e.g. 4d --> 3d). Otherwise, it keeps the singleton.

For faster computation

This is sometimes useful in case singletons need to be removed

via param squeeze_data (True by default to keep current behavior)

These modifications include the possibility to split along Z in case the input image is in sagittal orientation. The rest of the code is untouched, so that the motion correction as implemented previously is now encapsulated inside the z-splitted loop. UNFINISHED

need to be singleton

Fixes #1984

To avoid incessant warning messages.

After several unsuccessful testing using antsSliceRegularizedRegistration in 2D for moco applied to sagittal images, i realized that using antsRegistration in 2D mode, with affine transformation, gave best results. UNFINISHED: need to take care of the .mat output when using -g different from 1 (or maybe we could avoid supporting it)

Did some refactoring, notably made is_sagittal a field in param structure.

There are still a few checks to do, notably make sure the code runs when using mask with g>1 in 2D mode

Fixes #1878

...Which were introduced with the recent changes to msct_moco

Some cleaning, and miminized verbose (because it conflicts with the recent introduction of tqdm).

zougloub · 2018-10-04T13:26:26Z

scripts/sct_deepseg_sc.py

+
+    # compute z_centerline in image coordinates with continuous precision
+    voxel_coordinates = image_labels.transfo_phys2pix([[x_centerline_fit_rescorr[i], y_centerline_fit_rescorr[i], z_centerline_rescorr[i]] for i in range(len(z_centerline_rescorr))], real=False)
+    x_centerline_voxel_cont = [coord[0] for coord in voxel_coordinates]


Note that since transfo_phys2pix now returns a numpy array, you could have used something like x, y, z = voxel_coordinates.T

zougloub · 2018-10-04T13:29:52Z

scripts/sct_dmri_moco.py

@@ -369,36 +372,30 @@ def dmri_moco(param):

    # DWI groups
    file_dwi_mean = []
+    tqdm_bar = tqdm(total=nb_groups, unit='iter', unit_scale=False, desc="Merge within groups", ascii=True, ncols=80)
    for iGroup in range(nb_groups):


You could have used a tqdm iterator object instead of the range() here (avoids to update() and close())

i was not aware of that-- i'll look into it

fixed in c68600d

zougloub · 2018-10-04T13:31:58Z

scripts/sct_fmri_moco.py

 import time
 import math
-
+import shutil


There is a copy() function in sctutils

Fixed in 44e02e4

zougloub · 2018-10-04T13:35:17Z

scripts/msct_moco.py

+                im_maskz.save()
+                file_mask_splitZ.append(im_maskz.absolutepath)
+        # initialize file list for output matrices
+        file_mat = np.chararray([nz, nt],


chararray is deprecated and also not a good idea
https://docs.scipy.org/doc/numpy/reference/generated/numpy.chararray.html

after spending 10min on this i need help-- i've tried

file_mat = np.array([nz, nt], dtype="str")

but pb because string is immutable. solution?

zougloub · 2018-10-04T13:36:17Z

scripts/msct_moco.py

+        # Replace failed transformation with the closest good one
+        fT = [i for i, j in enumerate(failed_transfo) if j == 1]
+        gT = [i for i, j in enumerate(failed_transfo) if j == 0]
+        for it in range(len(fT)):


range(len(x)) -> enumerate(fT)

Fixed in 14f6794

zougloub · 2018-10-04T13:36:40Z

scripts/msct_moco.py

+        fT = [i for i, j in enumerate(failed_transfo) if j == 1]
+        gT = [i for i, j in enumerate(failed_transfo) if j == 0]
+        for it in range(len(fT)):
+            abs_dist = [abs(gT[i] - fT[it]) for i in range(len(gT))]


fixed in 7b3647c

zougloub · 2018-10-04T13:39:13Z

scripts/sct_fmri_moco.py

+        # Copy registration matrices
+        sct.printv('\nCopy transformations...', param.verbose)
+        for iGroup in range(nb_groups):
+            for data in range(len(group_indexes[iGroup])):


fixed in ef52dd3

actually we can't here-- fixed in c9aaa84

zougloub · 2018-10-04T13:41:02Z

scripts/sct_image.py

        im_out.absolutepath = sct.add_suffix(im_in.absolutepath, "_{}{}".format(dim_list[dim].upper(), str(idx_img).zfill(4)))
        im_out_list.append(im_out)

    return im_out_list


-def concat_data(fname_in_list, dim, pixdim=None):
+def concat_data(fname_in_list, dim, pixdim=None, squeeze_data=False):


Why do you absolutely want to have squeezing performed within this function and not elsewhere (caller or otherwise)?

squeeze=False, so I guess it is not performed-- or am i missing something?

just a question of design. why would a function named concat_data perform a squeeze.

We cannot use enumerate because group_indexes has 2 dim. So changing it back to the uglier range(len)

Otherwise I get the following error when running sct_testing: UnicodeEncodeError: 'ascii' codec can't encode characters in position 27-28: ordinal not in range(128)

sct_fmri_moco: Generalize motion correction for sagittal acquisitions and other improvements Former-commit-id: 55756cf

jcohenadad added sct_fmri_moco context: feature category: new functionality labels Sep 11, 2018

jcohenadad added this to the v3.2.4 milestone Sep 11, 2018

jcohenadad requested a review from zougloub September 11, 2018 20:19

jcohenadad changed the title ~~(WIP) sct_fmri_moco: Generalize motion correction for sagittal acquisitions~~ sct_fmri_moco: Generalize motion correction for sagittal acquisitions and other improvements. Sep 11, 2018

jcohenadad changed the title ~~sct_fmri_moco: Generalize motion correction for sagittal acquisitions and other improvements.~~ sct_fmri_moco: Generalize motion correction for sagittal acquisitions and other improvements Sep 11, 2018

jcohenadad modified the milestones: v3.2.4, v3.2.5 Sep 11, 2018

zougloub requested changes Sep 12, 2018

View reviewed changes

jcohenadad force-pushed the jca/1984-moco-sagittal branch from 4eefb4f to 604c68e Compare September 12, 2018 15:39

jcohenadad added enhancement category: improves performance/results of an existing feature and removed feature category: new functionality labels Oct 4, 2018

jcohenadad added 18 commits October 4, 2018 09:28

cleanup

df58e24

sct_image: Does not reshape on split

57ea9ee

This fixes issue #2017

sct_fmri_moco: no more data squeezing upon conversion

31740ec

Oherwise 4D data with singleton will become 3D data

sct_image: Fixed singleton removal based on dim

c195889

Now, if the splitting dimension corresponds to the last dim of the image, then it reshapes to remove the last dim (e.g. 4d --> 3d). Otherwise, it keeps the singleton.

msct_moco: removed old code

5262fb8

sct_fmri_moco: Does not average anymore if g=1

1d5d9fc

For faster computation

sct_image/concat: Introduced possibility to squeeze data

44076da

This is sometimes useful in case singletons need to be removed

sct_image: added possibility to not remove last singleton

747550e

via param squeeze_data (True by default to keep current behavior)

sct_maths: Fixed ineffective verbose option

65a5c0a

sct_image: added a condition on squeeze data

639d718

need to be singleton

sct_fmri_moco: (almost) finished implementing sagittal feature

edba5ca

Fixes #1984

sct_fmri_moco: Fixed wrong naming of copied matrix

2a3e09a

msct_moco: set verbose=0 when saving temp file

0db36c0

To avoid incessant warning messages.

sct_fmri_moco: Added more param flags accessible

ac60362

msct_moco: Now managing input mask with sagittal

bf4d64d

sct_utils: Added antsRegistration syntax to get_interpolation()

a05d6a0

jcohenadad added 8 commits October 4, 2018 09:28

sct_fmri_moco: Imposing group_size=1 for sagittal

bee0112

Did some refactoring, notably made is_sagittal a field in param structure.

msct_moco: Now accepting masks for 2D mode

276e4d3

There are still a few checks to do, notably make sure the code runs when using mask with g>1 in 2D mode

sct_image: Fixes wrong output file for copy-header

bc979dd

Fixes #1878

sct_fmri_moco: Updated default params

7134f85

sct_dmri_moco: Added missing field in param struct

7058551

sct_dmri_moco: Fixed compatibility issue

60b73bf

...Which were introduced with the recent changes to msct_moco

msct_moco: Cosmetic changes (use tqdm)

e5251a6

sct_dmri_moco/sct_fmri_moco/msct_moco: cleaning

272d3e2

Some cleaning, and miminized verbose (because it conflicts with the recent introduction of tqdm).

jcohenadad force-pushed the jca/1984-moco-sagittal branch from 7efa2bd to 272d3e2 Compare October 4, 2018 13:29

zougloub reviewed Oct 4, 2018

View reviewed changes

jcohenadad added 6 commits October 4, 2018 09:42

*_moco: Replaced tqdm by iterator

c68600d

sct_fmri_moco: replaced shutil.copy() by sct_utils

44e02e4

msct_moco: Use enumerate for loop

14f6794

sct_dmri/fmri_moco: use enumerate for loops

ef52dd3

msct_moco: Use numpy for abs

7b3647c

sct_dmri_moco,sct_fmri_moco: Change back from enumerate

d784521

We cannot use enumerate because group_indexes has 2 dim. So changing it back to the uglier range(len)

jcohenadad force-pushed the jca/1984-moco-sagittal branch from c9aaa84 to d784521 Compare October 4, 2018 17:28

*_moco: Had to change tqdm(ascii=True)

6fafce8

Otherwise I get the following error when running sct_testing: UnicodeEncodeError: 'ascii' codec can't encode characters in position 27-28: ordinal not in range(128)

zougloub approved these changes Oct 5, 2018

View reviewed changes

zougloub merged commit 55756cf into master Oct 5, 2018

jcohenadad deleted the jca/1984-moco-sagittal branch October 26, 2018 20:21

jcohenadad pushed a commit that referenced this pull request Dec 18, 2019

Merge pull request #2022 from neuropoly/jca/1984-moco-sagittal

32b973b

sct_fmri_moco: Generalize motion correction for sagittal acquisitions and other improvements Former-commit-id: 55756cf

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

sct_fmri_moco: Generalize motion correction for sagittal acquisitions and other improvements #2022

sct_fmri_moco: Generalize motion correction for sagittal acquisitions and other improvements #2022

jcohenadad commented Sep 10, 2018 •

edited

Loading

zougloub left a comment

zougloub Sep 12, 2018

jcohenadad Sep 12, 2018 •

edited

Loading

zougloub Oct 4, 2018

zougloub Oct 4, 2018

jcohenadad Oct 4, 2018

jcohenadad Oct 4, 2018

zougloub Oct 4, 2018

jcohenadad Oct 4, 2018

zougloub Oct 4, 2018

jcohenadad Oct 4, 2018

zougloub Oct 4, 2018

jcohenadad Oct 4, 2018

zougloub Oct 4, 2018

jcohenadad Oct 4, 2018

zougloub Oct 4, 2018

jcohenadad Oct 4, 2018

jcohenadad Oct 4, 2018

zougloub Oct 4, 2018

jcohenadad Oct 4, 2018 •

edited

Loading

zougloub Oct 5, 2018

sct_fmri_moco: Generalize motion correction for sagittal acquisitions and other improvements #2022

sct_fmri_moco: Generalize motion correction for sagittal acquisitions and other improvements #2022

Conversation

jcohenadad commented Sep 10, 2018 • edited Loading

zougloub left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jcohenadad Sep 12, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jcohenadad Oct 4, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jcohenadad commented Sep 10, 2018 •

edited

Loading

jcohenadad Sep 12, 2018 •

edited

Loading

jcohenadad Oct 4, 2018 •

edited

Loading