Bad phasediff SDC with fmriprep-21.0.0rc1

[Avijit-Chowdhury1]

What version of fMRIPrep are you using?

fmriprep-21.0.0rc1

What kind of installation are you using? Containers (Singularity, Docker), or "bare-metal"?

Singularity

What is the exact command-line you used?

#!/bin/bash

export SINGULARITYENV_TEMPLATEFLOW_HOME=/home/fmriprep/.cache/templateflow
/usr/bin/singularity run \
    --contain \
    --cleanenv \
    -B /tmp -B /data -B /data1 -B /data2 -B /data3 -B /cm/shared \
    /cm/shared/singularity/images/fmriprep-21.0.0rc1.simg \
    /data/achowdhury/CaseStudy/ \
    /data/achowdhury/CaseStudy/derivatives \
    participant \
    --fs-license-file /cm/shared/freesurfer-6.0.1/license.txt \
    --participant_label ID01 \
    --output-spaces MNI152NLin2009cAsym:res-2 anat func fsaverage \
    --n_cpus 8 \
    --mem-mb 65536 \
    --notrack \
    --skip_bids_validation \
    --fs-no-reconall \
    --t nirod \
    -vv \
    -w  /data/fmriprep-workdir
<place your command line here>

Have you checked that your inputs are BIDS valid?

Skip BIDS validation

Did fMRIPrep generate the visual report for this particular subject? If yes, could you share it?

https://www.dropbox.com/sh/poznbevp7dbmgpw/AACfaPHZhdbCehSnLGFp1BFda?dl=0

Can you find some traces of the error reported in the visual report (at the bottom) or in crashfiles?

No error
<please copy&paste all the information you can gather about errors>

Are you reusing previously computed results (e.g., FreeSurfer, Anatomical derivatives, work directory of previous run)?

Yes

fMRIPrep log

If you have access to the output logged by fMRIPrep, please make sure to attach it as a text file to this issue.
attached

I am preprocessing 7T fMRI data with BOLD resolution :

pixdim1 1.099138
pixdim2 1.099138
pixdim3 1.100000

Using phase-difference-based fieldmap correction.

T1 to MNI transform looks perfect. But BOLD to MNI normalization looks way off:

Here is normalized bold MASK overlayed on top of standard MNI:

I tried using --output-spaces MNI152NLin2009cAsym:res-native and the output is even worse.

Any idea why this is happening and how to fix it?

Thanks all

[Avijit-Chowdhury1]

update:

We also collected PE-polar fieldmap and using that mode of fieldmap correction results in slightly better normalization (although not perfect - cerebellum is out of normalized brain mask):

Further update: Using fMRIPREP version 20.2.5, using PE-polar fieldmap, normalization is better than 21.0.0.rc1 (not perfect):

Using fMRIPREP version 20.2.5, using phase-difference fieldmap, normalization is as bad as 21.0.0.rc1:

[oesteban]

The reports you provided show that GRE/phasediff fieldmaps were provided. There seems to be a problem in the direction of correction (opposite to what it should do, so enlarging the distortion instead of correcting for it) that has been fixed in the development branch, but not yet released as rc2 (#2576).

Could you share the visual report corresponding to the PEPOLAR attempt?

[Avijit-Chowdhury1]

The reports you provided show that GRE/phasediff fieldmaps were provided. There seems to be a problem in the direction of correction (opposite to what it should do, so enlarging the distortion instead of correcting for it) that has been fixed in the development branch, but not yet released as rc2 (#2576).

Could you share the visual report corresponding to the PEPOLAR attempt?

I resorted to using version 20.2.5 with PEPOLAR fieldmap, which seems to be giving the best (not accurate) output so far: https://www.dropbox.com/sh/evw50h7os8voceu/AAAl0JX9nlWPD0W8vdKPrAy2a?dl=0.

Also, the direction of the phase-encoding is reported wrongly as LEFT-RIGHT (when it should be A-P), as noted in question number #2553.

When is the rc2 release going to be available?

[effigies]

@Avijit-Chowdhury1 You can pull the nipreps/fmriprep:unstable docker image for the latest development build.

We're trying not to snow people under with RCs so are doing some extensive testing on our end before RC2.

[Avijit-Chowdhury1]

@effigies @oesteban

Tried using the rc2 unstable release, still getting the same results as rc1: https://www.dropbox.com/sh/t1njjxbh3nw11wa/AADMNM8NDh3uphphl5V-92ama?dl=0
.
I think overall phase-difference fieldmap correction is not working. Using PEpolar fieldmap the OFC area is being corrected for distortion while using phase-difference fieldmap it gets worse.

[oesteban]

I think overall phase-difference fieldmap correction is not working.

I disagree. Estimation of the fieldmap seems to be working, and the fact that the distortion is accentuated (i.e., correction is happening in the opposite direction it should) means that the field we estimate is physically plausible.

Now, since the estimation side of things is okay, what could be wrong? There seems to be occurring an undesired flip in the direction of correction. That can be due to the following factors:

• The PhaseEncodingDirection in the BIDS metadata is incorrect.
• fMRIPrep is incorrectly introducing the flip.

For us to debug this, we'll need access to some data. Would you be allowed to share with us some? To avoid issues, let's hold the T1w image off, this subject's fieldmap (and corresponding metadata) plus some BOLD you want to correct (it would be fine if you crop the BOLD after 50 timepoints so we genuinely don't have any other use of the data than resolving this bug. Would that be possible?

[Avijit-Chowdhury1]

@oesteban @effigies

Hi, sure, I have shared data of one subject with phasediff fieldmap here: https://www.dropbox.com/sh/2jzefzhpf0ez0nr/AAB0kp6_tFeJ8fJ4G3_nOTsna?dl=0

Please note that I raised this issue before in #2553 for the same dataset where PE direction was reported wrongly in the output HTML for fMRIPrep version: 20.2.3. Also, this is highres 7T data.

Your help is much appreciated.

[Avijit-Chowdhury1]

@oesteban @effigies
Hi, was wondering if you had a chance to look at the data?

[oesteban]

I haven't yet, I've been swamped this week. But I intend to do it as soon as next Monday. Thanks for your patience.

[Avijit-Chowdhury1]

@oesteban @effigies Would there be a way to debug this at my end while you are looking at the data?

[oesteban]

If you are using containers, the docker image nipreps/fmriprep:unstable should contain the latest code. Because we have seen this flip on some other data, I suspect it will remain with that docker image - but still definitely worth the try.

[Avijit-Chowdhury1]

@oesteban Thanks for your reply. I have tried using the latest code - the results are the same (bad SDC).

[Avijit-Chowdhury1]

@oesteban @effigies

Is there a document which specifies the intermediate files in the fmriprep working directory? For example, if I want to run the normalization by myself using FSL, which file should I selected from the working directory that has most preprocessing steps done (slice-timing, motion-correction, fieldmap correction), and I can just proceed to do normalization using FSL/SPM?

[oesteban]

Although possible, that would not be recommended and even less encouraged by us. Unless we are aware of a limitation that forces users to poke through the temporary directory, we would discourage you from doing so. The reason being that, when the time comes to write your paper, the methodological section will be extremely difficult to write so what you did is reproducible in any reasonable way.

Further reasons (which are related to the previous) include that it might seem appealing to run some stuff with fMRIPrep and then patch the stuff I don't like with other tools. The next reasonable step is to apply the patch only on those subjects where I like better my patch over fMRIPrep's result, and so on so forth. fMRIPrep intendedly limits the researcher degrees of freedom, surely at the cost of maybe not always getting the best result (which is again contingent on the viewer's eyes, rather than an objective data point).

So, long story short, unless it is very well justified, that practice is definitely discouraged.

[Avijit-Chowdhury1]

@oesteban @effigies Hi, any update on this?

Thanks again

[oesteban]

@Avijit-Chowdhury1

I finally managed to run your data through the following script (cc @effigies and @mgxd to check whether they can see some important mistake in it):

from pathlib import Path
import json
from nipype.pipeline import engine as pe
from niworkflows.interfaces.nibabel import SplitSeries
from niworkflows.workflows.epi.refmap import init_epi_reference_wf
from niworkflows.func.util import init_enhance_and_skullstrip_bold_wf
from sdcflows.workflows.outputs import init_fmap_derivatives_wf, init_fmap_reports_wf
from sdcflows.workflows.apply.registration import init_coeff2epi_wf
from sdcflows.workflows.apply.correction import init_unwarp_wf
from sdcflows.fieldmaps import FieldmapEstimation, FieldmapFile

datadir = Path.home() / "Downloads"
outdir = (Path() / "out").absolute()
in_phasediff = datadir / "ses-01/fmap/sub-ID01_ses-01_run-01_phasediff.nii.gz"
in_meta = datadir / "ses-01/fmap/sub-ID01_ses-01_run-01_phasediff.json"
in_magnitude1 = datadir / "ses-01/fmap/sub-ID01_ses-01_run-01_magnitude1.nii.gz"
in_magnitude2 = datadir / "ses-01/fmap/sub-ID01_ses-01_run-01_magnitude2.nii.gz"
in_epi = datadir / "ses-01/func/sub-ID01_ses-01_task-jhana_run-01_bold.nii.gz"
in_epi_meta = datadir / "ses-01/func/sub-ID01_ses-01_task-jhana_run-01_bold.json"

ref_wf = init_epi_reference_wf(omp_nthreads=32, auto_bold_nss=True)
ref_wf.inputs.inputnode.in_files = [str(in_epi)]
enhance_and_skullstrip_bold_wf = init_enhance_and_skullstrip_bold_wf()

phdiff_wf = FieldmapEstimation(
    [in_phasediff, in_magnitude1, in_magnitude2],
    bids_id="phasediff01",
).get_workflow()

fmap_derivatives_wf = init_fmap_derivatives_wf(
    output_dir=str(outdir),
    write_coeff=True,
    bids_fmap_id="phasediff01",
)
fmap_derivatives_wf.inputs.inputnode.source_files = [str(in_phasediff)]
fmap_derivatives_wf.inputs.inputnode.fmap_meta = [json.loads(in_meta.read_text())]

fmap_reports_wf = init_fmap_reports_wf(
    output_dir=str(outdir),
    fmap_type="phasediff",
)
fmap_reports_wf.inputs.inputnode.source_files = [str(in_phasediff)]

coreg_wf = init_coeff2epi_wf(omp_nthreads=32, write_coeff=True)

unwarp_wf = init_unwarp_wf(omp_nthreads=32)
unwarp_wf.inputs.inputnode.metadata = json.loads(in_epi_meta.read_text())
unwarp_wf.inputs.inputnode.hmc_xforms = [str(Path("identity.txt").absolute())] * 3

split = pe.Node(SplitSeries(in_file=str(in_epi)), name="split")

wf = pe.Workflow(name="my_wf")
wf.connect([
    (phdiff_wf, coreg_wf, [
        ("outputnode.fmap_ref", "inputnode.fmap_ref"),
        ("outputnode.fmap_mask", "inputnode.fmap_mask"),
        ("outputnode.fmap_coeff", "inputnode.fmap_coeff")]),
    (ref_wf, coreg_wf, [
        ("outputnode.epi_ref_file", "inputnode.target_ref")]),
    (ref_wf, enhance_and_skullstrip_bold_wf, [
        ("outputnode.epi_ref_file", "inputnode.in_file")
    ]),
    (enhance_and_skullstrip_bold_wf, coreg_wf, [
        ("outputnode.mask_file", "inputnode.target_mask")]),
    (phdiff_wf, fmap_reports_wf, [
        ("outputnode.fmap", "inputnode.fieldmap"),
        ("outputnode.fmap_ref", "inputnode.fmap_ref"),
        ("outputnode.fmap_mask", "inputnode.fmap_mask")]),
    (phdiff_wf, fmap_derivatives_wf, [
        ("outputnode.fmap", "inputnode.fieldmap"),
        ("outputnode.fmap_ref", "inputnode.fmap_ref"),
        ("outputnode.fmap_coeff", "inputnode.fmap_coeff"),
    ]),
    (coreg_wf, unwarp_wf, [
        ("outputnode.fmap_coeff", "inputnode.fmap_coeff")]),
    (split, unwarp_wf, [
        ("out_files", "inputnode.distorted")]),
])
wf.base_dir = str(Path() / "work")
wf.run()

I can confirm that the preprocessed fieldmap looks good, but then things fall apart when projecting the fieldmap into EPI space and applying the unwarping.

First thing I noticed is that coregistration between the magnitude and the EPI reference does not work with your parameters. We were initializing the registration with the affines found in the NIfTI headers and ITK/ANTs definitely doesn't like them. nipreps/sdcflows#253 seems to do the trick and I managed to get the magnitude and the boldref aligned.

Once the magnitude-boldref alignment worked out, the problem with the affines in the headers pops up again in correction. I will need more time to debug this.

[Avijit-Chowdhury1]

Hi, any updates on this? @oesteban @effigies

[effigies]

Unfortunately not, our focus has been on fixing more straightforward bugs with less chance of modifying the behavior on existing datasets.

One thing you can try doing to see if it works for you is to patch the modified SDCflows into your container:

git clone https://github.com/nipreps/sdcflows.git
git -C sdcflows checkout fix/magnitude-epi-coreg-params
singularity run -B $PWD/sdcflows/sdcflows:/opt/conda/lib/python3.8/site-packages/sdcflows ...

[effigies]

@Avijit-Chowdhury1 Have you by chance tried out the 22.0.x series?

[effigies]

Please reopen if this isn't resolved by 22.1.0.