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

[FIX] Replace FSL mcflirt with AFNI 3dVolReg #1354

Closed
wants to merge 4 commits into from

Conversation

oesteban
Copy link
Member

@oesteban oesteban commented Oct 27, 2018

Closes #1535.

A second attempt to get #1283 finally merged in.

@oesteban
Copy link
Member Author

Testing this here https://github.com/oesteban/fmriprep/tree/tests/fix-1354 to avoid excessive Circle usage.

However, I'm not positive the removal of rotation and shear does not have any additional side effects we haven't envisaged.

if self.inputs.remove_rotation_and_shear:
out_fname = fname_presuffix(self.inputs.in_file, suffix='_valid',
newpath=runtime.cwd)
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.

Could it be that we're not correcting to closest canonical before removing rotation and shear in this interface?

Copy link
Member Author

Choose a reason for hiding this comment

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

Copy link
Member Author

Choose a reason for hiding this comment

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

You were right on that we were not correcting to closest canonical, but that does not solve the problem (https://circleci.com/gh/oesteban/fmriprep/2636#artifacts/containers/0).

Question is: it seems to be working with the pipeline. Should we just disable the remove_rotation_and_shear for the tests only? Or we want to update the masks to include that?

Copy link
Member

Choose a reason for hiding this comment

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

So I think if we're removing rotation and shear, we need to be updating to RAS. Perhaps it makes sense to actually put that step directly into remove_rotation_and_shear, to prevent the case where we destroy orientation information.

If the output masks are valid, just in a different space, then we should update the targets.

@effigies effigies mentioned this pull request Oct 27, 2018
3 tasks

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
@oesteban
Copy link
Member Author

oesteban commented Oct 29, 2018

Status report:

  • Smoke tests on ds005, ds054 and ds210 look good
  • Regression tests of masks don't work unless the remove_rotation_and_shear is disabled. I see two options for these tests to pass: keep remove_rotation_and_shear off, or update both test datasets and masks with the orientation fixed.

There seem to be no side effects to the remove_rotation_and_shear in the workflow. I think we can go ahead and merge this (currently, with remove_rotation_and_shear disabled in the regression tests).

Opinions? @chrisfilo @effigies

import numpy as np
import nibabel as nb

img = nb.as_closest_canonical(img)
Copy link
Member

Choose a reason for hiding this comment

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

Oh, I see you did put this in here. Ignore the above.


DEFAULT_MEMORY_MIN_GB = 0.01


# pylint: disable=R0914
def init_bold_hmc_wf(mem_gb, omp_nthreads, name='bold_hmc_wf'):
def init_bold_hmc_wf(use_mcflirt, mem_gb, omp_nthreads, name='bold_hmc_wf'):
Copy link
Member

Choose a reason for hiding this comment

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

If we eventually want to be able to swap in any MC algorithm, use_mcflirt is not a great name for the parameter. What about motion_correction that takes a string such as '3dvolreg' or 'mcflirt' (and I would accept 'afni' and 'fsl' as synonyms, at least at the command line).

@@ -58,7 +58,8 @@ def symmetric_overlap(img1, img2):
)
])
def test_masking(input_fname, expected_fname):
bold_reference_wf = init_bold_reference_wf(omp_nthreads=1)
bold_reference_wf = init_bold_reference_wf(
omp_nthreads=1, remove_rotation_and_shear=False)
Copy link
Member

Choose a reason for hiding this comment

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

Since you're not setting this to False in the actual pipeline, I would not set it False here. It's not obvious that a change to the algorithm that causes a 2% shift without this step will cause the same 2% shift with the step.

Copy link
Member Author

Choose a reason for hiding this comment

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

On the other hand, this test was placed because we wanted to have changes on the masking algorithm under control. Adding the reorientation just embeds complexity into the masking process but should not affect masking performance. I see that as an integration problem, and I'm not sure we want to upload test data and masks after reorientation to replace the current test data.

What if we decide to fall back into mcflirt? Then this test would require rolling data back again.

Copy link
Member Author

Choose a reason for hiding this comment

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

@chrisfilo where does the need from removing rotation and shear come from?

Copy link
Contributor

Choose a reason for hiding this comment

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

AFNI ignores those when calculating transformations. If you keep them in the reference and use ANTs to apply the transformations estimated using AFNI the results will be misaligned.

Copy link
Member Author

@oesteban oesteban Oct 30, 2018

Choose a reason for hiding this comment

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

That makes sense. Do you think it would be possible to include rotation and shear here in the Volreg2ITK conversion instead?

That way we would avoid changing the image's original header.

Copy link
Contributor

Choose a reason for hiding this comment

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

Potentially, but that would be harder to implement.

Copy link
Member

Choose a reason for hiding this comment

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

What if we decide to fall back into mcflirt? Then this test would require rolling data back again.

As long as we always apply this cleaning to the entire BOLD series and , MCFLIRT will also be working from the cardinal RAS+ orientation. This would actually end up finishing up #260.

@effigies
Copy link
Member

Should re-raise the following point from #260 (comment):

If we change the orientation of input files we need to adjust the phase encoding direction as well.

I would also add slice-encoding direction (#1345).

I think we need a quick interface (here as a function) to translate incoming directions to reoriented directions:

def update_encoding_axes(img, pe_dir, slice_dir):
    import nibabel as nb

    ornt = nb.orientations.io_orientation(img)

    axes = ornt[:, 0]
    flips = ornt[:, 1] # -1 if flip, 1 if no flip

    if pe_dir is not None:
        pe_idx = 'ijk'.index(pe_dir[0])
        pe_sign = -1 if pe_dir[1:] else 1
        pe_sign *= flips[pe_idx]
        pe_idx = axes[pe_idx]
        pe_dir = 'ijk'[pe_idx] + '-'[pe_sign:]
        
    if slice_dir is not None:
        ...

    return pe_dir, slice_dir

I can flesh out the interface and write tests, if you agree this is a good idea.

@oesteban
Copy link
Member Author

Isn't safer/easier to calculate the analytic inverse of remove_rotation_and_shear and apply it when building the ITK transforms?

@effigies
Copy link
Member

It's all passing around the correct metadata, so either way is fine with me. I've found applying transforms in the right order has been a bit of a pain, in the past. And transforms that involve adjusting the affine but not reslicing the data matrix are very finicky (see, e.g., #842), so I'm not sure how simple re-rotating-and-shearing will end up being.

@oesteban
Copy link
Member Author

oesteban commented Oct 30, 2018

I agree that adjusting the affine without reslicing is really finicky.

Another concern to take into account. Merging #1334 (or any derivative of it) will make it very easy to implement an orig keyword to generate BOLD series in the original space. However, these outputs will not match the original orientation somewhat defeating the purpose of the reference. Does it make sense?

This is less of a problem if we set a boldref as reference. Since that reference would be produced by fMRIPrep, it is not a problem that the shear and rotation have been removed.

@effigies
Copy link
Member

I guess the question is whether we ever really care about the "original space" so much as "absolutely minimal resampling". I'm inclined to say it's actually the latter, in which case boldref is the better choice, and we should just be explicit about it, but perhaps there is a use-case for orig?

@effigies
Copy link
Member

Is this okay to push off to the next release?

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.

This is going to require merging master. And I think we have some open discussion points. Maybe after OHBM abstracts are in, we should sit down and make a clear plan forward for swappable motion correction.

@oesteban
Copy link
Member Author

👍

@dangom
Copy link
Contributor

dangom commented Jan 9, 2020

Hi guys,

related to #1928, I would love to see this PR merged in order to start experimenting with fMRIPREP for the robust preprocessing of high-res/partial FOV data.

I noticed, and correct me if I'm wrong, that one of the obstacles towards getting this PR merged is thatAFNI messes up the nifti headers and therefore may cause potential issues in downstream resampling/registration steps - in particular when interacting with ANTs.

Interestingly enough, a solution to this problem seems to have already been implemented by Freesurfer. If we look at the code of freesurfer's 3dVolReg wrapper, mc-afni2, we see that they use mri_convert to fix the 3dVolReg output. Would it not simplify the implementation of this PR simply swapping mcflirt for mc-afni2 or porting freesurfer's solution by piping the output of 3dVolReg through a call to mri_convert as well?

# Always convert the output. AFNI mucks with the goemetry,
# so we have to run mri_convert with --like to keep it.
# 3dvolreg will remove off-diagonal components (but keeps
# the basic geometry correct).
echo "#@# --------------------------------" |& tee -a $LF
set cmd = (mri_convert $newoutvol $outvol --in_like $tempvol)
pwd |& tee -a $LF
echo $cmd |& tee -a $LF
$cmd |& tee -a $LF
if($status) then
  echo "ERROR: converting output " |& tee -a $LF
  if($cleanup) then
    echo "... cleaning up ...."; |& tee -a $LF
    rm -r $tmpdir
  endif
  exit 1;
endif

@oesteban
Copy link
Member Author

oesteban commented Jan 9, 2020

We are addressing the issue in nipy/nitransforms#9 without tricky hacks - just interpreting the "deobliquing" transforms as 3dVolreg does. An extra pair of hands will be more than very much appreciated there :)

@dangom
Copy link
Contributor

dangom commented Jan 12, 2020

Thanks for your answer.

I can't be of much help at the moment in interpreting AFNI transforms, since I don't have experience with AFNI. That said, I've been comparing 3dVolReg (via mc-afni2) and antsMotionCorr for high-res /partial FOV motion correction. Both work very well, and up to minor interpolation differences, results are almost identical.

3dVolReg has the disadvantage that it doesn't accept reference images with a different matrix size / resolution. antsMotionCorr does the right thing out of the box, so replacing an sbref for a whole-brain ref would work without further tweaks. Motion-correcting directly towards a whole brain (or large FOV) reference makes intuitively more sense for cases of extreme partial FOV ( < 1 cm coverage), where concepts such as a "mean volume" start to break down ( under motion each time-point may cover different, slightly non-overlapping anatomy).
The downside is that whole-brain references are generally separate acquisitions, so motion parameter estimates would potentially have to take that into account.

antsMotionCorr has the disadvantage that it's not currently supported by Nipype. An orphaned old-PR to add support exists, though.

@oesteban
Copy link
Member Author

I'd be 👍 with bringing back the idea of antsMotionCorr - that stub in nipype was created when we were first playing with it.

@oesteban
Copy link
Member Author

oesteban commented Jun 5, 2020

This PR has gone stale and almost impossible to merge right now. Also, we will hopefully go around the oblique x-forms issue in NiTransforms proper, and this migration (to AFNI's tooling) will become pretty trivial.

It is important to say that #1535 stands as relevant as when Chris initially posted it if not more (given the unreliability of FD Power traces we get from mcflirt).

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.

Replace mcflirt with 3dVolreg
4 participants