Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[ENH] Replace FSL mcflirt with AFNI 3dVolReg #1283

Merged
merged 30 commits into from
Oct 26, 2018

Conversation

chrisgorgo
Copy link
Contributor

Copy link
Member

@effigies effigies left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is a 2005 comparison paper still a valid criterion for making such a move?

import nibabel as nb
import numpy as np
affine = img.affine
affine[:3, :3] = np.diag(np.diag(affine[:3, :3]))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This will simply zero out the off-diagonal, right? That will end up producing an affine with smaller zooms than the original. Probably the correct way to do this is to do a rotation/zoom/shear/translation decomposition, and then to remove the rotations and shears.

That said, this is less a de-oblique than a removal of information that describes obliqueness. Is this what you want, and what's the reason for that?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes the intention was to remove off diagonal values (because otherwise AFNI doesn't generate transforms that can be applied by ANTs). Is there a python code snippet for rotation/zoom/shear/translation decomposition I could use?

import numpy as np
affine = img.affine
affine[:3, :3] = np.diag(np.diag(affine[:3, :3]))
return nb.Nifti1Image(np.asanyarray(img.dataobj), affine, img.header)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would prefer img.__class__ to nb.Nifti1Image.

mcflirt = pe.Node(fsl.MCFLIRT(save_mats=True, save_plots=True),
name='mcflirt', mem_gb=mem_gb * 3)
mc = pe.Node(afni.Volreg(args='-prefix NULL -twopass',
zpad=4, outputtype='NIFTI_GZ'), name="mc", mem_gb=mem_gb * 3)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we need to specify outputtype if -prefix NULL?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, but we never use the output so there is no need for saving it (even if Nipype will delete it anyway).

@chrisgorgo
Copy link
Contributor Author

chrisgorgo commented Sep 14, 2018

Is a 2005 comparison paper still a valid criterion for making such a move?

Probably not alone.

@chrisgorgo
Copy link
Contributor Author

Tests are passing, but before I finish the PR by updating the documentation lets decide if we want to do this switch. @oesteban @effigies

@effigies
Copy link
Member

I'm okay with moving to 3dvolreg by default, primarily on the grounds of decreasing dependence on FSL-licensed code.

However, given that we've gone through the trouble of implementing a MCFLIRT pathway, does it make sense to leave it as an option? Particularly from the perspective of people eventually using fMRIPrep to compare the quality of different algorithms/implementations, it seems worth making this a swappable component.

@oesteban oesteban changed the title [WIP] Replace MCFLIRT with 3dVolReg [ENH] Replace FSL mcflirt with AFNI 3dVolReg Oct 3, 2018
@oesteban
Copy link
Member

oesteban commented Oct 4, 2018

I'm putting together the last release before RC1. My idea was to upgrade the version minor after that last release (i.e. start 1.2.x).

I think this is a somewhat large change that would fit better within the next 1.2.0, rather than 1.1.8.

WDYT?

@chrisgorgo
Copy link
Contributor Author

I agree it's fairly substantial.

Copy link
Member

@oesteban oesteban left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good overall. Left few comments.

@@ -105,6 +105,7 @@
'scikit-image',
'versioneer',
'pyyaml',
'transforms3d'
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please add trailing comma

requirements.txt Outdated
@@ -1,4 +1,4 @@
niworkflows>=0.4.3
git+https://github.com/poldracklab/niworkflows.git@fdf20103b9dedd4632bb8e26a5eb3b896e349fd2#egg=niworkflows
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

0.4.3 has this commit already, please revert this change (I guess the niworkflows dependency has been updated since you sent the PR)


fsl2itk = pe.Node(MCFLIRT2ITK(), name='fsl2itk',
mem_gb=0.05, n_procs=omp_nthreads)
afni2itk = pe.Node(niu.Function(function=afni2itk_func,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If there's nothing against it, I'll convert this function node into a VOLREG2ITK fully-fledged nipype interface before merging in, similarly to MCFLIRT2ITK. I can remember that parallelizing that one was tricky but paid off with an important speed up on long datasets.

@oesteban
Copy link
Member

This is waiting on my offline tests comparing master vs. this PR on a real dataset.

@effigies
Copy link
Member

@oesteban Any results from those tests?

@oesteban
Copy link
Member

Not yet. I've been blocked with other things. I should have a decent run on ds005 with current master, so I just need to test out this branch.

@oesteban
Copy link
Member

I think this is a thumbs up, see below (cc @effigies):

Master:
mcflirt-sub-01_task-mixedgamblestask_run-01_carpetplot

3dVolReg:
volreg-sub-01_task-mixedgamblestask_run-01_carpetplot

Master:
mcflirt-sub-01_task-mixedgamblestask_run-02_carpetplot

3dVolReg:
volreg-sub-01_task-mixedgamblestask_run-02_carpetplot

Master:
mcflirt-sub-01_task-mixedgamblestask_run-03_carpetplot

3dVolReg:
volreg-sub-01_task-mixedgamblestask_run-03_carpetplot

@effigies effigies added this to the 1.2.0 milestone Oct 24, 2018
Copy link
Member

@effigies effigies left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good. That said, I think there's a fairly major bug here. See below comment.

Also, did we want to permit selection of MCFLIRT/3dVolReg with a CLI flag (I agree let's make 3dVolReg default)? If someone did want to run an apples-to-apples comparison of the tools (such as a replication of Oakes, et al., 2005, and potentially relevant to @miykael's Neurostars post), that would be useful infrastructure.

imgs = nb.four_to_three(nb.load(in_file))
for i, img in enumerate(imgs):
out_file = os.path.abspath('vol%04d.nii.gz' % i)
img = remove_rotation_and_shear(img)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we need to do img = nb.as_closest_canonical(img) before this. Consider an image with the affine:

>>> img.affine
array([[  0.        ,   3.        ,   0.        , -88.4322052 ],
       [ -1.66738188,   0.        ,   2.6871357 , -47.71268463],
       [  3.35891914,   0.        ,   1.33390486, -77.85404205],
       [  0.        ,   0.        ,   0.        ,   1.        ]])
>>> img.shape
(33, 64, 64)
>>> remove_rotation_and_shear(img).affine
array([[  3.75000002,   0.        ,   0.        , -88.4322052 ],
       [  0.        ,   3.        ,   0.        , -47.71268463],
       [  0.        ,   0.        ,   3.00000007, -77.85404205],
       [  0.        ,   0.        ,   0.        ,   1.        ]])

>>> remove_rotation_and_shear(nb.as_closest_canonical(img)).affine
array([[  3.        ,   0.        ,   0.        , -88.4322052 ],
       [  0.        ,   3.00000007,   0.        , -47.71268463],
       [  0.        ,   0.        ,   3.75000002, -77.85404205],
       [  0.        ,   0.        ,   0.        ,   1.        ]])

Notice that both are RAS matrices, but with different zooms. And if you look at the shapes:

>>> remove_rotation_and_shear(img).shape
(33, 64, 64)

>>> remove_rotation_and_shear(nb.as_closest_canonical(img)).shape
(64, 64, 33)

We clearly have a rotated data matrix, which cannot be true if both correctly have an RAS orientation. So we end up destroying all orientation information, not just the "off-diagonal" information, unless we rotate the data matrix first.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch

@effigies
Copy link
Member

@oesteban
Copy link
Member

Good job folks

@oesteban oesteban merged commit b9576a5 into nipreps:master Oct 26, 2018
@oesteban
Copy link
Member

Sorry I was too fast:

https://5616-53175327-gh.circle-artifacts.com/0/tmp/data/reports/ds001362/sub-01_task-taskname_run-01_bold_masks.svg

With the new bold references, it seems data are changed the orientation. I'm unsure whether we just need to update the masks or this regression reflects a deeper problem. But it seems clear that it was introduced here: https://circleci.com/gh/chrisfilo/fmriprep/616#artifacts/containers/0

Sorry about that :(

@effigies
Copy link
Member

Hmm. Looks like the tests aren't failing properly with 0 overlap. Or possibly they're close enough voxel-wise but the affines are far off.

@oesteban
Copy link
Member

oesteban commented Oct 27, 2018 via email

oesteban added a commit that referenced this pull request Oct 27, 2018
oesteban added a commit that referenced this pull request Oct 27, 2018
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants