Susceptibility distortion correction without fieldmaps

chrisgorgo commented 7 years ago

Traditionally to perform susceptibility distortion correction (SDC) you need to have some sort of fieldmap data. However, there was some work showing that SDC can be performed without the need to acquire extra data.

Galiphour et. al 2008 proposed using a field distortion average template http://ieeexplore.ieee.org/document/4649099/ Another template can be found here: http://neurovault.org/collections/2426/
@juhuntenburg investigated non-linear coregistration between EPI and T1w of the same subject http://pubman.mpdl.mpg.de/pubman/item/escidoc:2327525:5/component/escidoc:2327523/master_thesis_huntenburg_4686947.pdf
Wang, @danjonpeterson et al. 2017 further improved and evaluated the non-linear coregistration technique http://journal.frontiersin.org/article/10.3389/fninf.2017.00017/full

It would be great if could add SDC to FMRIPREP that does not require fieldmaps. I would suggest starting, by implementing the approach from Wang et al. on a dataset with fieldmaps (for example Beast or UH2) and see how that compares. I suspect (judging from @juhuntenburg work) that some amount of regularization and use of priors will be required. Here are some ideas to try:

use 12 DOF (affine) boundary based coregistration as the initial transformation
use a field distortion template as the initial transformation
limit the non-linear coregistration only to the dimension of the phase encoding direction
create a mask of usually distorted areas from the distortion template an limit the coregistration using it

juhuntenburg commented 7 years ago

Hi! Something to consider is that, to my knowledge, this has only been shown successfully with DWI data, using the B0 image as a source for the registration step (e.g. in Wang et al. 2017). Based on my own work, I wouldn't recommend this as a standard preprocessing step for 3T BOLD data, which has less anatomical contrast, the results were extremely variable. It has worked really well for me though with 7T BOLD.

effigies commented 7 years ago

@chrisfilo So it seems like the place to do this would be with the reference EPI image, which should have the best contrast we can manage. Should we be re-using any prior registration steps, to keep it as apples-to-apples as possible? I'm not quite sure where we should put this in just yet.

effigies commented 7 years ago

Fieldmap unwarped on bias-corrected anatomical

fmap_unwarped

SyN-unwarped on bias-corrected anatomical

syn_unwarped

SyN-unwarped on fieldmap unwarped

syn_on_fmap

There are of course differences, but I don't know how to assess, on visual inspection, their severity.

This was run with:

/opt/ants/antsRegistrationSyN.sh -d 3 -f T1_inv.nii.gz -m /data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/epi_hmc_wf/gen_ref/ref_image.nii.gz -o EPItoT1 -t s

oesteban commented 7 years ago

How do you ensure that the deformation is only in the PE direction?

chrisgorgo commented 7 years ago

These figures are really hard to spot differences. Could you create set of animated reportlets via "SimpleBeforeAndAfter" interface in niworkflows?

On Wed, May 10, 2017 at 1:36 PM, Oscar Esteban notifications@github.com wrote:

How do you ensure that the deformation is only in the PE direction?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/poldracklab/fmriprep/issues/496#issuecomment-300604777, or mute the thread https://github.com/notifications/unsubscribe-auth/AAOkp9APhIRkaPf22gCcAVG_LP6QG2-1ks5r4h-2gaJpZM4NMxO0 .

effigies commented 7 years ago

@oesteban Haven't the slightest idea. ANTs documentation is almost entirely opaque to me, and I haven't seen any indication in what I have read that suggests ANTs can do this. I was starting from the command given in the Wang et al. supplementary materials.

@chrisfilo Sure.

effigies commented 7 years ago

@chrisfilo Huh. Well, they're at different resolutions, so those figures aren't much help.

 nib-ls /data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/sdc_unwarp_wf/apply_fov_mask/ref_image_corrected_fixhdr_trans_masked.nii.gz 
.../ref_image_corrected_fixhdr_trans_masked.nii.gz float32 [100, 100,  64] 2.20x2.20x2.20   sform

% nib-ls EPItoT1Warped.nii.gz 
EPItoT1Warped.nii.gz float64 [186, 256, 256] 0.90x0.90x0.90

I guess I'll need to try to figure out what's going on there.

danjonpeterson commented 7 years ago

The approach in Wang et al 2017 didn't restrict the deformation to the phase-encode direction, but this is something we're planning on looking into eventually, once we get around to it.

The -g (or --restrict-deformation) option on antsRegistration should do this (with -g 0x1x0 or 0.1x1x0.1)

The way to use this with the ANTS shell scripts (like antsRegistrationSyN.sh) is to ctrl-C the shell script call and look for the full antsRegistration call near the beginning of the output. These shell scrips are just wrappers for antsRegistration. Then you can add the more granular options like -g, etc. to the full call and run that. I've found the -v ( --write-interval-volumes ) options especially useful for troubleshooting registrations.

On Wed, May 10, 2017 at 1:54 PM, Chris Markiewicz notifications@github.com wrote:

@oesteban https://github.com/oesteban Haven't the slightest idea. ANTs documentation is almost entirely opaque to me, and I haven't seen any indication in what I have read that suggests ANTs can do this. I was starting from the command given in the Wang et al. supplementary materials.

@chrisfilo https://github.com/chrisfilo Sure.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/poldracklab/fmriprep/issues/496#issuecomment-300609480, or mute the thread https://github.com/notifications/unsubscribe-auth/AOj85BRSRfu9Csm_EIzTvj-6pnu6ptBmks5r4iQVgaJpZM4NMxO0 .

chrisgorgo commented 7 years ago

I think that's expected. You might need to reslice one of the images.

On May 10, 2017 2:06 PM, "Chris Markiewicz" notifications@github.com wrote:

@chrisfilo https://github.com/chrisfilo Huh. Well, they're at different resolutions, so those figures aren't much help.

nib-ls /data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/sdc_unwarp_wf/apply_fov_mask/ref_image_corrected_fixhdr_trans_masked.nii.gz .../ref_image_corrected_fixhdr_trans_masked.nii.gz float32 [100, 100, 64] 2.20x2.20x2.20 sform

% nib-ls EPItoT1Warped.nii.gz EPItoT1Warped.nii.gz float64 [186, 256, 256] 0.90x0.90x0.90

I guess I'll need to try to figure out what's going on there.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/poldracklab/fmriprep/issues/496#issuecomment-300612553, or mute the thread https://github.com/notifications/unsubscribe-auth/AAOkp3qD6e8jlUZsZ4qWX-tEgGuNyM1Yks5r4ibIgaJpZM4NMxO0 .

oesteban commented 7 years ago

Thanks a lot @danjonpeterson, I was right now looking at the --restrict-deformation option. https://github.com/stnava/ANTs/issues/216

chrisgorgo commented 7 years ago

Thanks @danjonpeterson - good to know! BTW have you tried this approach on non-diffusion data?

danjonpeterson commented 7 years ago

No, I haven't yet.

effigies commented 7 years ago

Resampled SyN-warped reference image to fieldmap-warped reference image (using -g 0x1x0):

https://imgh.us/report_1.svg

(Before = fieldmap-warped)

chrisgorgo commented 7 years ago

Thanks. It would be also good to see:

before/after the SyN-warp (analogous to before after applying fieldmap based distortion correction in current reports)
bbregister T1 coregistration reportlet (so we will see how well the corrected file overlaps with the anatomy)

In general there are some deformations that should not be there. Additionally the signal dropout from the ear canals got miraculously "recovered".

oesteban commented 7 years ago

I would have said that the miraculous recover of signal besides the ear canals is probably for the unrestricted deformation, but now I see that @effigies used the -g option. That's very surprising.

Also, if you look at the CC it is not corrected at all. What was the orientation of the input images to ants? It seems like maybe the j axis is not the second.

effigies commented 7 years ago

@oesteban Orientation LPS (inverted T1 is RAS). So j direction should be A/P in either case, which is what the json describes as the phase-encoding direction.

That said, perhaps ANTs cares about where exactly you place the -g flag? I put it at the end of the command, ending up with:

/opt/ants/antsRegistration --dimensionality 3 --float 0 --output '[EPItoT1,EPItoT1Warped.nii.gz]' --interpolation Linear --winsorize-image-intensities '[0.005,0.995]' --use-histogram-matching 0 --initial-moving-transform '[T1_inv.nii.gz,/data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/epi_hmc_wf/gen_ref/ref_image.nii.gz,1]' --transform 'Rigid[0.1]' --metric 'MI[T1_inv.nii.gz,/data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/epi_hmc_wf/gen_ref/ref_image.nii.gz,1,32,Regular,0.25]' --convergence '[1000x500x250x100,1e-6,10]' --shrink-factors 8x4x2x1 --smoothing-sigmas 3x2x1x0vox --transform 'Affine[0.1]' --metric 'MI[T1_inv.nii.gz,/data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/epi_hmc_wf/gen_ref/ref_image.nii.gz,1,32,Regular,0.25]' --convergence '[1000x500x250x100,1e-6,10]' --shrink-factors 8x4x2x1 --smoothing-sigmas 3x2x1x0vox --transform 'SyN[0.1,3,0]' --metric 'CC[T1_inv.nii.gz,/data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/epi_hmc_wf/gen_ref/ref_image.nii.gz,1,4]' --convergence '[100x70x50x20,1e-6,10]' --shrink-factors 8x4x2x1 --smoothing-sigmas 3x2x1x0vox -g 0x1x0

@chrisfilo

Before: Reference After: SyN-warped reference

I'll get the surface-overlaid images. I'm assuming you meant just overlay the surfaces, not actually run bbregister...

effigies commented 7 years ago

Surface-overlaid:

chrisgorgo commented 7 years ago

I would expect you have to specify the -g parameter per each stage: https://github.com/stnava/ANTs/wiki/Anatomy-of-an-antsRegistration-call - documentation is unclear.

This is actually not horrible, but the CC seems to be pulled into the signal dropout area in orbitofrontal cortex. I've seen similar errors when correcting images using blip up/down without restricting the deformations to the phase encoding direction (which is 3dQwarps default).

effigies commented 7 years ago

But don't we only want -g to apply to the SyN stage? If I'm reading the docs correctly, -g 0x1x0 on the rigid section means we're only allowing translation along the AP axis (or possibly rotation about it?). And I'm not sure what it would mean to restrict the affine transformation to 0x1x0.

Or should we not be performing the Rigid/Affine parts?

oesteban commented 7 years ago

You are right, -g should be defined only for the nonlinear stage.

That said, this could be the issue: we want to restrict the warp to the Y direction of the moving image. A way to do this:

Run first the rigid and affine registration levels using the T1 (or inverted T1) as fixed image. Then resample the T1 to the EPI space (possibly we want to keep resolution here, so you'll need to generate a reference image centered at the EPI volume).

Once you have your reference T1 in EPI space, then the Y axis will be aligned with the EPI, and you can run the nonlinear registration.

Sounds about right?

effigies commented 7 years ago

I'll give it a try. What could possibly go wrong?

effigies commented 7 years ago

Following @oesteban's advice, I mapped the reference EPI image to T1 using the rigid and affine steps:

antsRegistration --dimensionality 3 --float 0 --output '[Partial,Partial.nii.gz]' --interpolation Linear --winsorize-image-intensities '[0.005,0.995]' --use-histogram-matching 0 --initial-moving-transform '[T1_inv.nii.gz,/data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/epi_hmc_wf/gen_ref/ref_image.nii.gz,1]' --transform 'Rigid[0.1]' --metric 'MI[T1_inv.nii.gz,/data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/epi_hmc_wf/gen_ref/ref_image.nii.gz,1,32,Regular,0.25]' --convergence '[1000x500x250x100,1e-6,10]' --shrink-factors 8x4x2x1 --smoothing-sigmas 3x2x1x0vox --transform 'Affine[0.1]' --metric 'MI[T1_inv.nii.gz,/data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/epi_hmc_wf/gen_ref/ref_image.nii.gz,1,32,Regular,0.25]' --convergence '[1000x500x250x100,1e-6,10]' --shrink-factors 8x4x2x1 --smoothing-sigmas 3x2x1x0vox

I then upsampled the EPI into the T1 resolution:

from nilearn import image as nli
invt1 = nli.load_img('T1_inv.nii.gz')
nli.resample_img(ref_epi_image, target_affine=np.diag(invt1.header.get_zooms())).to_filename('epi_in_t1.nii.gz')

And performed the inverse affine transform to the inverted T1 to get a warp target:

antsApplyTransforms -i T1_inv.nii.gz -o T1_in_EPI.nii.gz --float -r epi_in_t1.nii.gz -t '[Partial0GenericAffine.mat,1]'

Finally, I applied the SyN step:

antsRegistration --dimensionality 3 --float 0 --output '[SyNOnly,SyNOnlyWarped.nii.gz]' --interpolation Linear --winsorize-image-intensities '[0.005,0.995]' --use-histogram-matching 0 --initial-moving-transform '[T1_in_EPI.nii.gz,/data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/epi_hmc_wf/gen_ref/ref_image.nii.gz,1]' --transform 'SyN[0.1,3,0]' --metric 'CC[T1_in_EPI.nii.gz,/data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/epi_hmc_wf/gen_ref/ref_image.nii.gz,1,4]' --convergence '[100x70x50x20,1e-6,10]' --shrink-factors 8x4x2x1 --smoothing-sigmas 3x2x1x0vox -g 0x1x0

Anatomical comparison

Before-after warp

Fieldmap vs SyN warp

oesteban commented 7 years ago

I would recommend doing the first resampling in the opposite way. Keeping your first command:

antsRegistration --dimensionality 3 --float 0 --output '[Partial,Partial.nii.gz]' --interpolation Linear --winsorize-image-intensities '[0.005,0.995]' --use-histogram-matching 0 --initial-moving-transform '[T1_inv.nii.gz,/data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/epi_hmc_wf/gen_ref/ref_image.nii.gz,1]' --transform 'Rigid[0.1]' --metric 'MI[T1_inv.nii.gz,/data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/epi_hmc_wf/gen_ref/ref_image.nii.gz,1,32,Regular,0.25]' --convergence '[1000x500x250x100,1e-6,10]' --shrink-factors 8x4x2x1 --smoothing-sigmas 3x2x1x0vox --transform 'Affine[0.1]' --metric 'MI[T1_inv.nii.gz,/data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/epi_hmc_wf/gen_ref/ref_image.nii.gz,1,32,Regular,0.25]' --convergence '[1000x500x250x100,1e-6,10]' --shrink-factors 8x4x2x1 --smoothing-sigmas 3x2x1x0vox

Then use the generated inverse transform and antsApplyTransform to move the T1w into EPI space. Here, I'd recommend to use as reference image an upsampled EPI which resolution is close to the resolution of the T1w image.

Once you have the T1w in the EPI space, then run the SyN registration.

The reason for this convoluted process is that, when you resample the EPI into the T1w space, you are implicitly rotating the axes. Thus the phase encoding direction is not aligned with an axis anymore and the -g flag does not have the expected impact.

chrisgorgo commented 7 years ago

What kind of skull stripping was applied to both images?

effigies commented 7 years ago

@chrisfilo The inverted T1 was masked with the ANTs skull stripping from fmriprep. The reference EPI wasn't skull-stripped... I suppose I ought to have done that.

@oesteban I think I did what you're saying I should have.

oesteban commented 7 years ago

In the first step you are resampling the EPI into T1w, if I understood correctly, right?

effigies commented 7 years ago

Yes, but only to get the affine transform to move the T1 into alignment with the EPI.

oesteban commented 7 years ago

Oh sorry, I can see that: T1_in_EPI.nii.gz and /data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/epi_hmc_wf/gen_ref/ref_image.nii.gz.

You are right: you are doing what I suggested. I'll double check the results, but it seems like the SyN step ~~is not doing any major warping right~~?

EDIT: actually yes, and in the opposite direction we would've expected.

oesteban commented 7 years ago

How does the T1_in_EPI.nii.gz look like? It is the T1 with the inverted contrast right?

Here it could help a lot to have not only the CC metric, but also the Mattes metric. To see if this idea may work, you can just replace the CC metric on your SyN command by the Mattes. It'll go a lot faster as well.

effigies commented 7 years ago

What would you like to see T1_in_EPI overlaid on? T1_inv? The reference EPI?

It looks fine. Aligns well with the reference EPI. And yes, it's the inverted contrast.

I'll try the final step with the Mattes metric.

oesteban commented 7 years ago

Yes I'd say that for the last before/after plot it'd be clearer to see the two images involved in the registration. Meaning, the inv T1 and the warped EPI

On May 12, 2017 09:48, "Chris Markiewicz" notifications@github.com wrote:

What would you like to see T1_in_EPI overlaid on? T1_inv? The reference EPI?

It looks fine. Aligns well with the reference EPI. And yes, it's the inverted contrast.

I'll try the final step with the Mattes metric.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/poldracklab/fmriprep/issues/496#issuecomment-301128125, or mute the thread https://github.com/notifications/unsubscribe-auth/AAkhxmzzIuGUzIQhwRJdICQNvHVKnWLlks5r5I1GgaJpZM4NMxO0 .

effigies commented 7 years ago

So at the end, do we need to do one more rigid transform from the warped EPI to the original T1, or just apply the Partial0GenericAffine.mat from earlier like so:

/opt/ants/antsApplyTransforms -i EPI_ref.nii.gz -o EPI_final.nii.gz -r T1_inv.nii.gz -t Partial0GenericAffine.mat -t Mattes1Warp.nii.gz

oesteban commented 7 years ago

If you want to use this exclusively for fieldmap correction, the answer would be no. You want the EPI in its EPI space but without correction. Maybe I don't understand the question.

If you were referring to my last comment: the before/after plot would have the T1w-inv_in_EPI and the EPI-warped. So the T1w-inv_in_EPI already has the Partial0GenericAffine.mat implicitly applied when you resampled, right?

effigies commented 7 years ago

Mattes metric produced blank images.

Here's the results with skull-stripping. This includes the T1w-inv_in_EPI to EPI-warped comparison.

EPI skull stripping, EPI resampling, and T1 skull-stripping and inversion:

#!/usr/bin/env python
import numpy as np
from nilearn import image as nli

T1_orig = '/data/out/uh2/derivatives/fmriprep/sub-s497/anat/sub-s497_T1w_preproc.nii.gz'
T1_mask = '/data/out/uh2/derivatives/fmriprep/sub-s497/anat/sub-s497_T1w_brainmask.nii.gz'
epi_ref = '/data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/epi_hmc_wf/gen_ref/ref_image.nii.gz'
epi_mask = '/data/scratch/uh2/fmriprep_wf/single_subject_s497_wf/func_preproc_ses_1_task_short_run_1_wf/epi_hmc_wf/skullstrip_epi/ref_image_corrected_fixhdr_brain.nii.gz'

T1_img = nli.load_img(T1_orig)
mask = nli.load_img(T1_mask).get_data().astype(bool)

# Mask EPI image
masked_epi = nli.math_img('epi * mask', epi=epi_ref, mask=epi_mask).get_data()
epi_ref = nli.new_img_like(epi_ref, masked_epi, copy_header=True)
epi_ref.to_filename('EPI_ref.nii.gz')

# Upsample EPI image
epi_rescaled = nli.resample_img(epi_ref, target_affine=np.diag(T1_img.header.get_zooms()))
epi_rescaled.to_filename('EPI_upsampled.nii.gz')

masked_t1 = T1_img.get_data() * mask

# Calculate scale factors
t1_min, t1_max = np.unique(masked_t1)[[1, -1]]
epi_min, epi_max = np.unique(masked_epi)[[1, -1]]
scale_factor = (epi_max - epi_min) / (t1_max - t1_min)

# Invert, rescale and mask T1 image
t1_inv = nli.new_img_like(T1_img, mask * ((t1_max -T1_img.get_data()) * scale_factor + epi_min), copy_header=True)
t1_inv.to_filename('T1_inv.nii.gz')

Creation of target (Inverted T1 registered to EPI):

#!/bin/bash

/opt/ants/antsRegistration --dimensionality 3 --float 0 --output '[Partial,Partial.nii.gz]' \
    --interpolation Linear --winsorize-image-intensities '[0.005,0.995]' --use-histogram-matching 0 \
    --initial-moving-transform '[T1_inv.nii.gz,EPI_ref.nii.gz,1]' \
    --transform 'Rigid[0.1]' --metric 'MI[T1_inv.nii.gz,EPI_ref.nii.gz,1,32,Regular,0.25]' --convergence '[1000x500x250x100,1e-6,10]' \
        --shrink-factors 8x4x2x1 --smoothing-sigmas 3x2x1x0vox \
    --transform 'Affine[0.1]' --metric 'MI[T1_inv.nii.gz,EPI_ref.nii.gz,1,32,Regular,0.25]' --convergence '[1000x500x250x100,1e-6,10]' \
        --shrink-factors 8x4x2x1 --smoothing-sigmas 3x2x1x0vox

/opt/ants/antsApplyTransforms -i T1_inv.nii.gz -o T1_in_EPI.nii.gz --float -r EPI_upsampled.nii.gz -t '[Partial0GenericAffine.mat,1]'

SyN registration (CC metric):

#!/bin/bash

/opt/ants/antsRegistration --dimensionality 3 --float 0 --output '[CC,CCWarped.nii.gz]' \
    --interpolation Linear --winsorize-image-intensities '[0.005,0.995]' --use-histogram-matching 0 \
    --initial-moving-transform '[T1_in_EPI.nii.gz,EPI_ref.nii.gz,1]' \
    --transform 'SyN[0.1,3,0]' --metric 'CC[T1_in_EPI.nii.gz,EPI_ref.nii.gz,1,4]' --convergence '[100x70x50x20,1e-6,10]' \
        --shrink-factors 8x4x2x1 --smoothing-sigmas 3x2x1x0vox -g 0x1x0

This produces the raw transform:

To compare to the original (inverted) anatomical, we additionally apply the affine transform:

Shown against the original anatomical:

For consistency, here it is shown against the fieldmap-corrected EPI:

chrisgorgo commented 7 years ago

It seems that we are seeing overfitting. I highlighted the issues below:

We need to regularize it better. Limit how much deformation can be applied or limit which regions should be affected by the deformation (using average fieldmap template)

nipreps / fmriprep