fmriprep icon indicating copy to clipboard operation
fmriprep copied to clipboard

Bad phasediff SDC with fmriprep-21.0.0rc1

Open Avijit-Chowdhury1 opened this issue 3 years ago • 18 comments

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:

image

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 avatar Oct 18 '21 15:10 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):

image

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

image

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

Avijit-Chowdhury1 avatar Oct 19 '21 00:10 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?

oesteban avatar Oct 22 '21 14:10 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?

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?

Avijit-Chowdhury1 avatar Oct 22 '21 14:10 Avijit-Chowdhury1

@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.

effigies avatar Oct 22 '21 15:10 effigies

@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.

Avijit-Chowdhury1 avatar Oct 24 '21 17:10 Avijit-Chowdhury1

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?

oesteban avatar Oct 31 '21 16:10 oesteban

@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 avatar Nov 01 '21 01:11 Avijit-Chowdhury1

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

Avijit-Chowdhury1 avatar Nov 05 '21 13:11 Avijit-Chowdhury1

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.

oesteban avatar Nov 05 '21 13:11 oesteban

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

Avijit-Chowdhury1 avatar Nov 15 '21 21:11 Avijit-Chowdhury1

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.

oesteban avatar Nov 16 '21 08:11 oesteban

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

Avijit-Chowdhury1 avatar Nov 16 '21 12:11 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?

Avijit-Chowdhury1 avatar Nov 17 '21 16:11 Avijit-Chowdhury1

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.

oesteban avatar Nov 17 '21 16:11 oesteban

@oesteban @effigies Hi, any update on this?

Thanks again

Avijit-Chowdhury1 avatar Nov 22 '21 17:11 Avijit-Chowdhury1

@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.

oesteban avatar Nov 29 '21 15:11 oesteban

Hi, any updates on this? @oesteban @effigies

Avijit-Chowdhury1 avatar Jan 20 '22 03:01 Avijit-Chowdhury1

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 avatar Jan 21 '22 14:01 effigies

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

effigies avatar Dec 07 '22 19:12 effigies

Please reopen if this isn't resolved by 22.1.0.

effigies avatar Dec 14 '22 03:12 effigies