Problem with new version of fmriprep (See update below): error while running processing of fieldmaps

Hello everyone!

Summary of what happened:

I want to run fmriprep on imaging data that has been acquired sagittally.
When I ignore the fieldmaps, everything runs nicely. However, with fieldmaps I get an error:
“ValueError: Shape of coefficients file [33 46 46] does not meet the expected shape [46. 46. 33.] (topup factors are [2. 2. 2.])”

For the fieldmaps I use two spin echo EPI sequences. For one of the two EPI sequences the phase encoding direction is inverted compared to the other one.
The error I get occurs during preprocessing of the fieldmaps.

Command used :

I use the following command stored in a bash script to run fmriprep:

 S_FolderOut=$HOME/data_nas04/VMatten/01_NameOfExperiment/ProcessedData/Experiment/fMRIprep
  S_FolderBIDS=$HOME/data_nas04/VMatten/01_NameOfExperiment/ProcessedData/Experiment/niftiUnprocessed
  S_FolderLicense=$HOME/data_nas04/VMatten/01_NameOfExperiment/Code/fMRIprep/license.txt
  S_FolderWorkingDir=/tmp/VMatten02/
  nSubj=1
  subj01=AVR024
  numCPUs=48
  nthreads=$((numCPUs / nSubj))
  mem_mb=45000


    docker run --rm -e DOCKER_VERSION_8395080871=20.10.21 -it \
    -v ${S_FolderLicense}:/opt/freesurfer/license.txt:ro \
    -v ${S_FolderBIDS}:/data:ro \
    -v ${S_FolderOut}:/out \
    -v ${S_FolderWorkingDir}:/scratch \
    -u $UID \
    nipreps/fmriprep:22.1.1 /data /out participant --participant-label  ${subj01} \
    --nthreads ${nthreads} \
    --stop-on-first-crash \
    --mem ${mem_mb} \
    --use-aroma \
    -w /scratch

Version and Environment:

I am running fmriprep via docker: DOCKER_VERSION_8395080871=20.10.21
The version of fmriprep is fmriprep:22.1.1

Data formatted according to a validatable standard?

I use dcm2niix and a custom matlab script to convert the DICOM data to BIDS standards. The correct formatting of the data is tested with this website: https://bids-standard.github.io/bids-validator/

Important note: for the fielmaps I need to add the tag PhaseEncodingDirection manually. Our Scanner does not give proper names to the fieldmaps such that dcm2niix automatically detects the phase encoding direction.
However, this is a known problem and I checked the correctness of this tag as far as I could.

Relevant log outputs (up to 20 lines):

File: /out/sub-AVR024/log/20230123-211135_78dd9716-9128-4c13-b20c-307a55d607aa/crash-20230124-012147-UID1136-fix_coeff-c5bedbd9-ba15-4dfc-ac55-080388ec52f5.txt
Working Directory: /scratch/fmriprep_22_1_wf/single_subject_AVR024_wf/fmap_preproc_wf/wf_auto_00000/fix_coeff
Inputs:
fmap_ref: /scratch/fmriprep_22_1_wf/single_subject_AVR024_wf/fmap_preproc_wf/wf_auto_00000/setwise_avg/sub-AVR024_ses-01_dir-PA_epi_average_merged_average.nii
in_coeff: ['/scratch/fmriprep_22_1_wf/single_subject_AVR024_wf/fmap_preproc_wf/wf_auto_00000/topup/reoriented_base_fieldcoef.nii.gz']
pe_dir: j
Traceback (most recent call last):
  File "/opt/conda/lib/python3.9/site-packages/nipype/pipeline/plugins/multiproc.py", line 344, in _send_procs_to_workers
    self.procs[jobid].run(updatehash=updatehash)
  File "/opt/conda/lib/python3.9/site-packages/nipype/pipeline/engine/nodes.py", line 527, in run
    result = self._run_interface(execute=True)
  File "/opt/conda/lib/python3.9/site-packages/nipype/pipeline/engine/nodes.py", line 645, in _run_interface
    return self._run_command(execute)
  File "/opt/conda/lib/python3.9/site-packages/nipype/pipeline/engine/nodes.py", line 771, in _run_command
    raise NodeExecutionError(msg)
nipype.pipeline.engine.nodes.NodeExecutionError: Exception raised while executing Node fix_coeff.

Traceback:
	Traceback (most recent call last):
	  File "/opt/conda/lib/python3.9/site-packages/nipype/interfaces/base/core.py", line 398, in run
	    runtime = self._run_interface(runtime)
	  File "/opt/conda/lib/python3.9/site-packages/sdcflows/interfaces/bspline.py", line 472, in _run_interface
	    self._results["out_coeff"] = [
	  File "/opt/conda/lib/python3.9/site-packages/sdcflows/interfaces/bspline.py", line 474, in 
	    _fix_topup_fieldcoeff(
	  File "/opt/conda/lib/python3.9/site-packages/sdcflows/interfaces/bspline.py", line 555, in _fix_topup_fieldcoeff
	    raise ValueError(
	ValueError: Shape of coefficients file [33 46 46] does not meet the expected shape [46. 46. 33.] (toupup factors are [2. 2. 2.]).

I am pretty new to fmriprep and at the moment I have no idea what to try next to solve the problem.
I am grateful for any suggestions or hints!
Does anyone know what the exact source of this error is and/or how to solve it? Could it be related to the sagittal image acquisition?

Best,
Gwen


Update:

After running the exact same command with an older version of fmriprep namely 22.0.2 instead of the current one 22.1.1 the error disappeared.

What would you advice me to do now? Simply run my data with the older version or is there an error in my data?

If the dataset works on 22.0.2, I would probably go with that for now. Would you be willing to share a subject of your dataset privately so we can identify the regression?

Hi @Gwen, so just looking at your dataset, I’m seeing:

Data shapes:

$ nib-ls **/*.nii
sub-AVR025/ses-01/anat/sub-AVR025_ses-01_T1w.nii                    int16 [256, 320, 320]      0.75x0.75x0.75        sform
sub-AVR025/ses-01/fmap/sub-AVR025_ses-01_dir-AP_epi.nii             int16 [ 86,  86,  60,   5] 2.50x2.50x2.50x0.57   sform
sub-AVR025/ses-01/fmap/sub-AVR025_ses-01_dir-PA_epi.nii             int16 [ 86,  86,  60,   5] 2.50x2.50x2.50x0.57   sform
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023a_bold.nii int16 [ 90,  90,  60, 540] 2.50x2.50x2.50x0.57   sform
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023b_bold.nii int16 [ 90,  90,  60, 540] 2.50x2.50x2.50x0.57   sform
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023_bold.nii  int16 [ 90,  90,  60, 540] 2.50x2.50x2.50x0.57   sform
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023c_bold.nii int16 [ 90,  90,  60, 540] 2.50x2.50x2.50x0.57   sform
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023d_bold.nii int16 [ 90,  90,  60, 540] 2.50x2.50x2.50x0.57   sform
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023e_bold.nii int16 [ 90,  90,  60, 540] 2.50x2.50x2.50x0.57   sform
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023f_bold.nii int16 [ 90,  90,  60, 540] 2.50x2.50x2.50x0.57   sform
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023g_bold.nii int16 [ 90,  90,  60, 540] 2.50x2.50x2.50x0.57   sform

PhaseEncodingDirections (missing for BOLD data):

$ grep -rI PhaseEncodingDirection **/*.json
sub-AVR025/ses-01/anat/sub-AVR025_ses-01_T1w.json:	"InPlanePhaseEncodingDirectionDICOM": "ROW",
sub-AVR025/ses-01/fmap/sub-AVR025_ses-01_dir-AP_epi.json:  "InPlanePhaseEncodingDirectionDICOM": "ROW",
sub-AVR025/ses-01/fmap/sub-AVR025_ses-01_dir-AP_epi.json:  "PhaseEncodingDirection": "j-"
sub-AVR025/ses-01/fmap/sub-AVR025_ses-01_dir-PA_epi.json:  "InPlanePhaseEncodingDirectionDICOM": "ROW",
sub-AVR025/ses-01/fmap/sub-AVR025_ses-01_dir-PA_epi.json:  "PhaseEncodingDirection": "j"
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023a_bold.json:	"InPlanePhaseEncodingDirectionDICOM": "ROW",
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023b_bold.json:	"InPlanePhaseEncodingDirectionDICOM": "ROW",
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023_bold.json:	"InPlanePhaseEncodingDirectionDICOM": "ROW",
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023c_bold.json:	"InPlanePhaseEncodingDirectionDICOM": "ROW",
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023d_bold.json:	"InPlanePhaseEncodingDirectionDICOM": "ROW",
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023e_bold.json:	"InPlanePhaseEncodingDirectionDICOM": "ROW",
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023f_bold.json:	"InPlanePhaseEncodingDirectionDICOM": "ROW",
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023g_bold.json:	"InPlanePhaseEncodingDirectionDICOM": "ROW",

Orientations:

In [1]: import nibabel as nb

In [2]: from glob import glob

In [3]: for fname in glob("**/*.nii", recursive=True):
    ...:     img = nb.load(fname)
    ...:     axcodes = ''.join(nb.aff2axcodes(img.affine))
    ...:     freq, phase, slc = img.header.get_dim_info()
    ...:     print(f"{fname:70s}\tOrientation:{axcodes}\t\tPE-Axis:{axcodes[phase] if phase is not None else None}")
    ...: 
sub-AVR025/ses-01/anat/sub-AVR025_ses-01_T1w.nii                      	Orientation:RAS		PE-Axis:R
sub-AVR025/ses-01/fmap/sub-AVR025_ses-01_dir-AP_epi.nii               	Orientation:PSL		PE-Axis:P
sub-AVR025/ses-01/fmap/sub-AVR025_ses-01_dir-PA_epi.nii               	Orientation:PSL		PE-Axis:P
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023a_bold.nii   	Orientation:PSL		PE-Axis:None
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023b_bold.nii   	Orientation:PSL		PE-Axis:None
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023c_bold.nii   	Orientation:PSL		PE-Axis:None
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023d_bold.nii   	Orientation:PSL		PE-Axis:None
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023e_bold.nii   	Orientation:PSL		PE-Axis:None
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023f_bold.nii   	Orientation:PSL		PE-Axis:None
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023g_bold.nii   	Orientation:PSL		PE-Axis:None
sub-AVR025/ses-01/func/sub-AVR025_ses-01_task-sssiNew2023_bold.nii    	Orientation:PSL		PE-Axis:None

So the first thing to note is that your fieldmaps have PhaseEncodingDirection set to j/j- in the BIDS metadata, which would correspond to inferior-superior and superior-inferior. The NIfTI headers agree with your labeling that these should be along the anterior-posterior axis, so you probably want to label them as i and i-. It seems likely that we will still have an orientation issue once they’re properly labeled, but you should update your metadata.

Also, your BOLD files do not have phase-encoding information either in the JSON or the NIfTI. Is this known? If not, you’ll be producing a preprocessed fieldmap, but there will be no way to apply them to your BOLD files.

While this is producing a workflow failure that we will debug as we have time, for now I will consider this a dataset issue as resolving it will not fix your preprocessing.

1 Like

Hi @effigies ,
Thank you so much for this helpful reply!
Thank you for being patient with such a mistake and in particular for the short python script. As a Newbie to fmri this really helps me to continue working with the data.
I will change things accordingly, test it and come back to you afterwards.

Hi @effigies
I corrected the PhaseEncodingDirection for the fieldmaps and fixed the problem with the missing PhaseEncodingDirection for the functional runs. Thank you for your help to correct these problem!
However, with the current version of fmriprep, I still receive the exact same error message.
Could there be further issues with my dataset? And if so in which direction should I search? Or should I continue with the older version of fmriprep?
Best,
Gwen

Hi @Gwen. If you wouldn’t mind sending me the corrected dataset, I can look at it next week.

Hi @effigies
Thank you so much for your offer to look at the data!
For your information, in case it might be helpful:
I successfully ran toptup from fsl to generate a fieldmap.
Best,
Gwen

Hi Gwen, just had another quick look. The PhaseEncodingDirection is still absent from the BOLD files. I should be able to debug the TOPUP issue with what I have, but the fieldmap won’t be able to be applied without knowing i or i- for the BOLD files.

Hi effigies,
Thank you for looking at the data!

I am not sure, I understand what you mean.
It is true, that for the fieldmaps dicom2niix still does not find the tag PhaseEncodingDirection. So I added it manually.

However, the problem regarding the orientation in the functional images should be solved.

When I run

In [1]:  import nibabel as nb

In [2]: from glob import glob

In [3]: for fname in glob("**/*.nii", recursive=True):
    ...:     img = nb.load(fname)
    ...:     axcodes = ''.join(nb.aff2axcodes(img.affine))
    ...:     freq, phase, slc = img.header.get_dim_info()
    ...:     print(f"{fname:70s}\tOrientation:{axcodes}\t\tPE-Axis:{axcodes[phase] if phase is not None else None}")

I get for the following:

sub-AVR025_ses-01_task-isss_bold.nii                                 Orientation:PSL   PE-Axis:P
sub-AVR025_ses-01_task-isssa_bold.nii                                Orientation:PSL  PE-Axis:P
sub-AVR025_ses-01_task-isssb_bold.nii                                Orientation:PSL  PE-Axis:P
sub-AVR025_ses-01_task-isssc_bold.nii                                Orientation:PSL  PE-Axis:P
sub-AVR025_ses-01_task-isssd_bold.nii                                Orientation:PSL  PE-Axis:P
sub-AVR025_ses-01_task-issse_bold.nii                                Orientation:PSL  PE-Axis:P
sub-AVR025_ses-01_task-isssf_bold.nii                                Orientation:PSL  PE-Axis:P
sub-AVR025_ses-01_task-isssg_bold.nii                                Orientation:PSL  PE-Axis:P

Best,
Gwen

Yes, the NIfTI can encode the phase-encoding axis, but it cannot tell us which direction (AP or PA) the image is in. I was using the header check to establish consistency. fMRIPrep will not use that, as it is unreliable. PhaseEncodingDirection in the JSON sidecar is mandatory for fieldmap correction.

ok, that was a misunderstanding on my side. I am sorry.
So I not only need to add the tag to the fieldmaps but also the functional images?
I can try that and start fmriprep again.

Yes. Note, however that it will not resolve the TOPUP issue, which I have successfully reproduced:

230223-14:32:10,504 nipype.workflow INFO:                                                                                                                                         
         [Node] Setting-up "fmriprep_23_0_wf.single_subject_AVR025_wf.fmap_preproc_wf.wf_auto_00000.fix_coeff" in "/scratch/fmriprep_23_0_wf/single_subject_AVR025_wf/fmap_preproc
_wf/wf_auto_00000/fix_coeff".                                                                                                                                                     
230223-14:32:10,512 nipype.workflow INFO:                                                                                                                                         
         [Node] Executing "fix_coeff" <sdcflows.interfaces.bspline.TOPUPCoeffReorient>                                                                                            230223-14:32:10,535 nipype.workflow INFO:                                                                                                                                                  [Node] Finished "fix_coeff", elapsed time 0.021339s.                                                                                                                     
230223-14:32:10,535 nipype.workflow WARNING:                                                                                                                                      
         Storing result file without outputs                                                                                                                                      
230223-14:32:10,536 nipype.workflow WARNING:                                                                                                                                      
         [Node] Error on "fmriprep_23_0_wf.single_subject_AVR025_wf.fmap_preproc_wf.wf_auto_00000.fix_coeff" (/scratch/fmriprep_23_0_wf/single_subject_AVR025_wf/fmap_preproc_wf/w
f_auto_00000/fix_coeff)                                                                                                                                                           
230223-14:32:10,536 nipype.workflow ERROR:                                                                                                                                                 Node fix_coeff failed to run on host 17a2da4ab666.                                                                                                                       
230223-14:32:10,537 nipype.workflow ERROR:                                                                                                                                        
         Saving crash info to /out/sub-AVR025/log/20230223-142713_804f3698-696d-4c0a-9bf1-1b6e1177bbc7/crash-20230223-143210-root-fix_coeff-17e3d614-83a1-41aa-bfeb-631fa6775953.t
xt                                                                                       
Traceback (most recent call last):                                                                                                                                                  File "/opt/conda/lib/python3.9/site-packages/nipype/pipeline/plugins/multiproc.py", line 344, in _send_procs_to_workers                                                         
    self.procs[jobid].run(updatehash=updatehash)                                                                                                                                  
  File "/opt/conda/lib/python3.9/site-packages/nipype/pipeline/engine/nodes.py", line 527, in run                                                                                 
    result = self._run_interface(execute=True)                                                                                                                                    
  File "/opt/conda/lib/python3.9/site-packages/nipype/pipeline/engine/nodes.py", line 645, in _run_interface                                                                      
    return self._run_command(execute)                                                                                                                                             
  File "/opt/conda/lib/python3.9/site-packages/nipype/pipeline/engine/nodes.py", line 771, in _run_command                                                                            raise NodeExecutionError(msg)                                                                                                                                                 
nipype.pipeline.engine.nodes.NodeExecutionError: Exception raised while executing Node fix_coeff.                                                                                 
                                                                                                                                                                                  
Traceback:                                                                                                                                                                        
        Traceback (most recent call last):                                                                                                                                                  File "/opt/conda/lib/python3.9/site-packages/nipype/interfaces/base/core.py", line 398, in run                                                                          
            runtime = self._run_interface(runtime)                                                                                                                                
          File "/opt/conda/lib/python3.9/site-packages/sdcflows/interfaces/bspline.py", line 472, in _run_interface                                                               
            self._results["out_coeff"] = [                                                                                                                                        
          File "/opt/conda/lib/python3.9/site-packages/sdcflows/interfaces/bspline.py", line 474, in <listcomp>                                                                   
            _fix_topup_fieldcoeff(                                                                                                                                                
          File "/opt/conda/lib/python3.9/site-packages/sdcflows/interfaces/bspline.py", line 555, in _fix_topup_fieldcoeff                                                                    raise ValueError(                                                                                                                                                     
        ValueError: Shape of coefficients file [33 46 46] does not meet the expected shape [46. 46. 33.] (toupup factors are [2. 2. 2.]). 

I will try to track that down as I can. I hope to get to it early next week.

@Gwen Just a heads up that I’m working on this here: FIX: Handle fieldmap orientations apart from R/LAS by effigies · Pull Request #339 · nipreps/sdcflows · GitHub

Looking at your images, they appear to be collected A>>P, since the smearing extends out the back of the head instead of the front. You probably want to set your PhaseEncodingDirection to "i".

ok, thank you so much for this information!
Then I will continue for now with manually adding this tag to the json files.
Thank you for working on the topic!

Dear effigies,
Thank you so much for your amazing support!
I downloaded the new version 23.0.0 and ran it for my current subjects. All passed successfully!
Again thank you so much for your fast and contructive help! A lot better than for many paid programs out there.
I am very grateful.
Best,
Gwen