`omitted axis "4"` Error in dwipreproc of Siemens Vida (XA11A) data

error
preprocessing

#1

Dear mrtrix community,

In preprocessing diffusion data from our new Siemens Vida I got stuck all with the following error in the dwipreproc

mrconvert: [ERROR] omitted axis "4" has dimension greater than 1

I’ll explain the background and how I got there.

Some while ago I used the mrtrix preprocessing on SMS data from an older Trio. It all worked out fine. We got a new Siemens Vida and had to deal with significant delays (associated with severe problems from the Numaris-syngo Version XA10). Updates (Numaris-syngo Version XA11A) fixed that now. I was starting back up with my Kurtosis pipeline.

The Data is acquired with a product SMS MDDW with b = 0 1k 2.5k s/mm² in 30 standard Siemens half sphere directions and additional 9 b0.
Everything of the above is acquired in both AP and PA.
(I know - not optimal, full sphere would be better, I promise change :wink: By the way: Siemens seems to have largely reworked the interface for using a custom b-table: If anyone out there already has experience with the XA11A please let me know!)
The b-table I got for now hence for either PE direction looks like:

  • 1xb0
  • 30xb1000
  • 30xb2500
  • 9xb0

I then do a mrconvert to put the DICOMs into an AP.mif and a PA.mif, followed by a mrcat AP.mif PA.mif AP_PA.mif-axis 3

followed by:

$ dwipreproc -pe_dir AP AP_PA.mif AP_PA_preproc.mif -rpe_all -topup_options "--verbose "

Although I am doing the work on an i7 then takes topup foreeeeeeeeeeeeeever (eddy is finished astonishingly quickly thanks to CUDA …Any way to parallelize topup?) and gives the following output …

dwipreproc: 
dwipreproc: Note that this script makes use of commands / algorithms that have relevant articles for citation; INCLUDING FROM EXTERNAL SOFTWARE PACKAGES. Please consult the help page (-help option) for more information.
dwipreproc: 
dwipreproc: Generated temporary directory: /home/ingomar/Desktop/from_stick/dwipreproc-tmp-6XKITW/
Command:  mrconvert "/home/ingomar/Desktop/from_stick/AP_PA.mif" "/home/ingomar/Desktop/from_stick/dwipreproc-tmp-6XKITW/dwi.mif"
dwipreproc: Changing to temporary directory (/home/ingomar/Desktop/from_stick/dwipreproc-tmp-6XKITW/)
Command:  dirstat dwi.mif -output asym
dwipreproc: Total readout time not provided at command-line; assuming sane default of 0.1
Command:  mrconvert dwi.mif -import_pe_table dwi_manual_pe_scheme.txt - | dwiextract - se_epi.mif -bzero
Command:  mrinfo dwi.mif -export_grad_mrtrix grad.b
Command:  mrconvert se_epi.mif topup_in.nii -strides -1,+2,+3,+4 -export_pe_table topup_datain.txt
Command:  topup --imain=topup_in.nii --datain=topup_datain.txt --out=field --fout=field_map.nii.gz --config=/usr/local/fsl/etc/flirtsch/b02b0.cnf --verbose
Command:  mrconvert dwi.mif -import_pe_table dwi_manual_pe_scheme.txt - | mrinfo - -export_pe_eddy applytopup_config.txt applytopup_indices.txt
Command:  dwiextract dwi.mif -import_pe_table dwi_manual_pe_scheme.txt -pe 0,-1,0,0.1 - | mrconvert - dwi_pe_1.nii -json_export dwi_pe_1.json
Command:  applytopup --imain=dwi_pe_1.nii --datain=applytopup_config.txt --inindex=1 --topup=field --out=dwi_pe_1_applytopup.nii --method=jac
Command:  mrconvert dwi_pe_1_applytopup.nii.gz dwi_pe_1_applytopup.mif -json_import dwi_pe_1.json
Command:  dwiextract dwi.mif -import_pe_table dwi_manual_pe_scheme.txt -pe 0,1,0,0.1 - | mrconvert - dwi_pe_2.nii -json_export dwi_pe_2.json
Command:  applytopup --imain=dwi_pe_2.nii --datain=applytopup_config.txt --inindex=2 --topup=field --out=dwi_pe_2_applytopup.nii --method=jac
Command:  mrconvert dwi_pe_2_applytopup.nii.gz dwi_pe_2_applytopup.mif -json_import dwi_pe_2.json
Command:  mrcat dwi_pe_1_applytopup.mif dwi_pe_2_applytopup.mif - -axis 3 | dwi2mask - - | maskfilter - dilate - | mrconvert - eddy_mask.nii -datatype float32 -strides -1,+2,+3
Command:  mrconvert dwi.mif -import_pe_table dwi_manual_pe_scheme.txt eddy_in.nii -strides -1,+2,+3,+4 -export_grad_fsl bvecs bvals -export_pe_eddy eddy_config.txt eddy_indices.txt
Command:  eddy_cuda --imain=eddy_in.nii --mask=eddy_mask.nii --acqp=eddy_config.txt --index=eddy_indices.txt --bvecs=bvecs --bvals=bvals --topup=field --out=dwi_post_eddy
dwipreproc: Detected matching DWI volumes with opposing phase encoding; performing explicit volume recombination

… before throwing an error:


Command:  mrcalc field_map.nii.gz 0.1 -mult -1.0 -mult - | mrfilter - gradient - | mrconvert - field_deriv_pe_1.mif -coord 3 1 -axes 0,1,2
dwipreproc:
dwipreproc: [ERROR] Command failed: mrcalc field_map.nii.gz 0.1 -mult -1.0 -mult - | mrfilter - gradient - | mrconvert - field_deriv_pe_1.mif -coord 3 1 -axes 0,1,2 (dwipreproc:996)
dwipreproc: Output of failed command:
mrconvert: [ERROR] omitted axis "4" has dimension greater than 1
dwipreproc: 
dwipreproc: Changing back to original directory (/home/ingomar/Desktop/from_stick/)
dwipreproc: Script failed while executing the command: mrcalc field_map.nii.gz 0.1 -mult -1.0 -mult - | mrfilter - gradient - | mrconvert - field_deriv_pe_1.mif -coord 3 1 -axes 0,1,2
dwipreproc: For debugging, inspect contents of temporary directory: /home/ingomar/Desktop/from_stick/dwipreproc-tmp-6XKITW/

As a first step in working towards handling the error I tried to speed up the thing by providing it with a -se_epi b0 image pair from closely adjoint b0 (last b0 from AP and first from PA) and align it with -align_seepi, that however ended in:

$ dwipreproc -pe_dir AP AP_PA.mif AP_PA_prepreproc.mif -rpe_all -se_epi se_epi_2vol.mif 
dwipreproc: 
dwipreproc: Note that this script makes use of commands / algorithms that have relevant articles for citation; INCLUDING FROM EXTERNAL SOFTWARE PACKAGES. Please consult the help page (-help option) for more information.
dwipreproc: 
dwipreproc: Generated temporary directory: /home/ingomar/Desktop/from_stick/dwipreproc-tmp-GVVHAU/
Command:  mrconvert "/home/ingomar/Desktop/from_stick/AP_PA.mif" "/home/ingomar/Desktop/from_stick/dwipreproc-tmp-GVVHAU/dwi.mif"
Command:  mrconvert "/home/ingomar/Desktop/from_stick/se_epi_2vol.mif" "/home/ingomar/Desktop/from_stick/dwipreproc-tmp-GVVHAU/se_epi.mif"
dwipreproc: Changing to temporary directory (/home/ingomar/Desktop/from_stick/dwipreproc-tmp-GVVHAU/)
Command:  dirstat dwi.mif -output asym
dwipreproc: Total readout time not provided at command-line; assuming sane default of 0.1

dwipreproc: [ERROR] If explicitly including SE EPI images when using -rpe_all option, they must come with their own associated phase-encoding information in the image header

dwipreproc: Changing back to original directory (/home/ingomar/Desktop/from_stick/)
dwipreproc: Contents of temporary directory kept, location: /home/ingomar/Desktop/from_stick/dwipreproc-tmp-GVVHAU/)

I then by using
mrconvert AP_PA.mif -import_pe_table dwi_manual_pe_scheme.txt - | dwiextract - - -bzero | mrconvert - se_epi_2vol.mif -coord 3 9,10 -force included the following PE table I got from one of the tmp directories:

  • 70x 0.0 -1.0 0.0 0.1
  • 70x 0.0 1.0 0.0 0.1

However then I got this Error

)$ dwipreproc -pe_dir AP AP_PA.mif AP_PA_prepreproc.mif -rpe_all -se_epi se_epi_2vol.mif -align_seepi
dwipreproc: 
dwipreproc: Note that this script makes use of commands / algorithms that have relevant articles for citation; INCLUDING FROM EXTERNAL SOFTWARE PACKAGES. Please consult the help page (-help option) for more information.
dwipreproc: 
dwipreproc: Generated temporary directory: /home/cercare/Desktop/from_stick/New Folder (14)/dwipreproc-tmp-NR98QH/
Command:  mrconvert "/home/cercare/Desktop/from_stick/New Folder (14)/AP_PA.mif" "/home/cercare/Desktop/from_stick/New Folder (14)/dwipreproc-tmp-NR98QH/dwi.mif"
Command:  mrconvert "/home/cercare/Desktop/from_stick/New Folder (14)/se_epi_2vol.mif" "/home/cercare/Desktop/from_stick/New Folder (14)/dwipreproc-tmp-NR98QH/se_epi.mif"
dwipreproc: Changing to temporary directory (/home/cercare/Desktop/from_stick/New Folder (14)/dwipreproc-tmp-NR98QH/)
Command:  dirstat dwi.mif -output asym
dwipreproc: Total readout time not provided at command-line; assuming sane default of 0.1
Command:  mrconvert dwi.mif dwi_first_bzero.mif -coord 3 0 -axes 0,1,2
dwipreproc: Balanced phase-encoding scheme detected in SE-EPI series; volume 0 will be removed and replaced with first b=0 from DWIs
Command:  mrconvert se_epi.mif - -coord 3 1 | mrcat dwi_first_bzero.mif - se_epi_firstdwibzero.mif -axis 3
Command:  mrinfo dwi.mif -export_grad_mrtrix grad.b
Command:  mrconvert se_epi_firstdwibzero.mif topup_in.nii -strides -1,+2,+3,+4 -export_pe_table topup_datain.txt
dwipreproc: 
dwipreproc: [ERROR] Command failed: mrconvert se_epi_firstdwibzero.mif topup_in.nii -strides -1,+2,+3,+4 -export_pe_table topup_datain.txt (dwipreproc:736)
dwipreproc: Output of failed command:
            mrconvert: [ERROR] no phase-encoding information found within image "topup_in.nii"
dwipreproc: 
dwipreproc: Changing back to original directory (/home/cercare/Desktop/from_stick/New Folder (14))
dwipreproc: Script failed while executing the command: mrconvert se_epi_firstdwibzero.mif topup_in.nii -strides -1,+2,+3,+4 -export_pe_table topup_datain.txt
dwipreproc: For debugging, inspect contents of temporary directory: /home/cercare/Desktop/from_stick/New Folder (14)/dwipreproc-tmp-NR98QH/

Desperately grasping that It is better to have slightly misaligned than no data I dropped the -align_seepi option.
(Any Ideas how to proceed without dropping the -align_seepi?!)

Topup with the -se_epi is much quicker (minutes instead of hours) but I still get the same error as in the beginning:

Command:  mrconvert dwi.mif -import_pe_table dwi_manual_pe_scheme.txt eddy_in.nii -strides -1,+2,+3,+4 -export_grad_fsl bvecs bvals -export_pe_eddy eddy_config.txt eddy_indices.txt
Command:  eddy_cuda --imain=eddy_in.nii --mask=eddy_mask.nii --acqp=eddy_config.txt --index=eddy_indices.txt --bvecs=bvecs --bvals=bvals --topup=field --out=dwi_post_eddy
dwipreproc: Detected matching DWI volumes with opposing phase encoding; performing explicit volume recombination
Command:  mrcalc field_map.nii.gz 0.1 -mult -1.0 -mult - | mrfilter - gradient - | mrconvert - field_deriv_pe_1.mif -coord 3 1 -axes 0,1,2
dwipreproc: 
dwipreproc: [ERROR] Command failed: mrcalc field_map.nii.gz 0.1 -mult -1.0 -mult - | mrfilter - gradient - | mrconvert - field_deriv_pe_1.mif -coord 3 1 -axes 0,1,2 (dwipreproc:996)
dwipreproc: Output of failed command:
mrconvert: [ERROR] omitted axis "4" has dimension greater than 1

So I am back to square one (however I little quicker :face_with_raised_eyebrow: ) - this brings me here.

I most grateful for any help!:grinning:

p.s.: One additional related question:

When trying to use eddy’s slice2volume correction via $ dwipreproc -pe_dir AP AP_PA.mif AP_PA_prepreproc.mif -rpe_all -se_epi se_epi_2vol.mif -eddy_options "--mporder=# "I get the error …

dwipreproc: No slice encoding direction information present; assuming third axis corresponds to slices
dwipreproc: [ERROR] Cannot perform slice-to-volume correction in eddy: No slspec file provided, and no slice timing information present in header

… although I am using a full .mif workflow within mrtrix. In the manual for dwipreproc it is indicated that one need not provide a --slspec file if doing so.
Anyone having a guess if this is due to an oversight on my side or an incompatiblity with the Vida-XA11A-DICOMs.

I very much appreciate your help!

Thank you all,
Ingomar


#2

Thanks for the very complete problem description – I’m going to use this as a model of how it’s done… :+1:

In terms of the underlying problem, there’s a quite a few things going on here, and you’re trying quite a few different strategies that @rsmith may be better placed to advise on (dwipreproc is primarily his baby – even if he sometimes wishes it wasn’t… :wink:). But the main error:

mrcalc field_map.nii.gz 0.1 -mult -1.0 -mult - | mrfilter - gradient - | mrconvert - field_deriv_pe_1.mif -coord 3 1 -axes 0,1,2
mrconvert: [ERROR] omitted axis "4" has dimension greater than 1

suggests that the field_map.nii.gz file is 4D, with more than one volume – rather than the expected single 3D volume. This cause mrfilter gradient to output a 5D image (dx,dy,dz for each volume), and when mrconvert tries to extract the second volume (dy), but selects only the x,y,z, axes to output a 3D image, mrconvert complains that the image has size > 1 along the 5th dimension (axis 4, indexed from zero). I can’t think of many other ways you’d get this error.

Looking further up, the field_map.nii.gz is generated directly by the topup call using the --fout option:

Command:  topup --imain=topup_in.nii --datain=topup_datain.txt --out=field --fout=field_map.nii.gz --config=/usr/local/fsl/etc/flirtsch/b02b0.cnf --verbose

Unfortunately, I can’t figure out from the topup help page why (or indeed whether) topup would ever produce a 4D output.

So I guess the first thing is to inspect the field_map.nii.gz in the temporary folder created by dwipreproc – assuming you still have access to it. If it isn’t 4D, then I’m not sure what’s going on…

Sorry, can’t help on that front…

Not that I know, unfortunately. We also find it can be very slow, particularly when there’s lots of b=0 images and all of them are used in the estimation (this also explains why your explicit -se_epi call was quicker: only 2 b=0 volumes instead of 10). It’s single-threaded, but I don’t know of any option to parallelise it – if anyone knows about this, let us know and we’ll switch it on!

This relies on the MRtrix3 DICOM handling code having imported the slice timings correctly to figure out the correct slice grouping. Check whether your AP_PA_prepreproc.mif file contains the SliceTiming entry. If not, then the DICOM data either don’t contain the information, or Siemens have changed the way it’s stored and we’d need to update our code to support that (we may need to get access to the data to work this out…).


#3

Thank you very much for your quick reply so close to the ISMRM deadline :wink:

Yes, indeed, I can confirm that it is a 4d file:thinking:

$ mrinfo field_map.nii.gz 
************************************************
Image:               "field_map.nii.gz"
************************************************
  Dimensions:        88 x 88 x 60 x 20
  Voxel size:        2.5 x 2.5 x 2.5 x 1
  Data strides:      [ -1 2 3 4 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:                0.993      0.0477      0.1081      -114.5
                         -0.01882       0.967     -0.2539      -107.7
                          -0.1166      0.2501      0.9612      -96.71
  comments:          6.0.0

It comprises 1 fieldmap and 20 (transformed?) T2 low res (==b0?) images.

So I thought it could have to do with me putting in concatenated AP PA data with the following structure

  • 1xb0
  • 30xb1000
  • 30xb2500
  • 9xb0
  • 1xb0
  • 30xb1000
  • 30xb2500
  • 9xb0

Trying to circumvent any potential related issues I did cut out the middle b0 volumes via

mrconvert -coord 3 0:60,70:130 AP_PA.mif AP_PA_concatenated.mif

and did rerun the preprocessing …

$ dwipreproc -pe_dir AP AP_PA_concatenated.mif AP_PA_concatenated_preproc.mif -rpe_all
dwipreproc: 
dwipreproc: Note that this script makes use of commands / algorithms that have relevant articles for citation; INCLUDING FROM EXTERNAL SOFTWARE PACKAGES. Please consult the help page (-help option) for more information.
dwipreproc: 
dwipreproc: Generated temporary directory: /home/ingomar/from_stick/dwipreproc-tmp-HLRCGM/
Command:  mrconvert " /home/ingomar/from_stick/AP_PA_concatenated.mif" " /home/ingomar/from_stick/dwipreproc-tmp-HLRCGM/dwi.mif"
dwipreproc: Changing to temporary directory ( /home/ingomar/from_stick/dwipreproc-tmp-HLRCGM/)
Command:  dirstat dwi.mif -output asym
dwipreproc: Total readout time not provided at command-line; assuming sane default of 0.1
Command:  mrconvert dwi.mif -import_pe_table dwi_manual_pe_scheme.txt - | dwiextract - se_epi.mif -bzero
Command:  mrinfo dwi.mif -export_grad_mrtrix grad.b
Command:  mrconvert se_epi.mif topup_in.nii -strides -1,+2,+3,+4 -export_pe_table topup_datain.txt
Command:  topup --imain=topup_in.nii --datain=topup_datain.txt --out=field --fout=field_map.nii.gz --config=/usr/local/fsl/etc/flirtsch/b02b0.cnf
Command:  mrconvert dwi.mif -import_pe_table dwi_manual_pe_scheme.txt - | mrinfo - -export_pe_eddy applytopup_config.txt applytopup_indices.txt
Command:  dwiextract dwi.mif -import_pe_table dwi_manual_pe_scheme.txt -pe 0,-1,0,0.1 - | mrconvert - dwi_pe_1.nii -json_export dwi_pe_1.json
Command:  applytopup --imain=dwi_pe_1.nii --datain=applytopup_config.txt --inindex=1 --topup=field --out=dwi_pe_1_applytopup.nii --method=jac
Command:  mrconvert dwi_pe_1_applytopup.nii.gz dwi_pe_1_applytopup.mif -json_import dwi_pe_1.json
Command:  dwiextract dwi.mif -import_pe_table dwi_manual_pe_scheme.txt -pe 0,1,0,0.1 - | mrconvert - dwi_pe_2.nii -json_export dwi_pe_2.json
Command:  applytopup --imain=dwi_pe_2.nii --datain=applytopup_config.txt --inindex=2 --topup=field --out=dwi_pe_2_applytopup.nii --method=jac
Command:  mrconvert dwi_pe_2_applytopup.nii.gz dwi_pe_2_applytopup.mif -json_import dwi_pe_2.json
Command:  mrcat dwi_pe_1_applytopup.mif dwi_pe_2_applytopup.mif - -axis 3 | dwi2mask - - | maskfilter - dilate - | mrconvert - eddy_mask.nii -datatype float32 -strides -1,+2,+3
Command:  mrconvert dwi.mif -import_pe_table dwi_manual_pe_scheme.txt eddy_in.nii -strides -1,+2,+3,+4 -export_grad_fsl bvecs bvals -export_pe_eddy eddy_config.txt eddy_indices.txt
Command:  eddy_cuda --imain=eddy_in.nii --mask=eddy_mask.nii --acqp=eddy_config.txt --index=eddy_indices.txt --bvecs=bvecs --bvals=bvals --topup=field --out=dwi_post_eddy

… but it ran into the same problem as before …

dwipreproc: Detected matching DWI volumes with opposing phase encoding; performing explicit volume recombination
Command:  mrcalc field_map.nii.gz 0.1 -mult -1.0 -mult - | mrfilter - gradient - | mrconvert - field_deriv_pe_1.mif -coord 3 1 -axes 0,1,2
dwipreproc: 
dwipreproc: [ERROR] Command failed: mrcalc field_map.nii.gz 0.1 -mult -1.0 -mult - | mrfilter - gradient - | mrconvert - field_deriv_pe_1.mif -coord 3 1 -axes 0,1,2 (dwipreproc:996)
dwipreproc: Output of failed command:
            mrconvert: [ERROR] omitted axis "4" has dimension greater than 1
dwipreproc: 
dwipreproc: Changing back to original directory (/home/ingomar/from_stick/)
dwipreproc: Script failed while executing the command: mrcalc field_map.nii.gz 0.1 -mult -1.0 -mult - | mrfilter - gradient - | mrconvert - field_deriv_pe_1.mif -coord 3 1 -axes 0,1,2
dwipreproc: For debugging, inspect contents of temporary directory: /home/ingomar/from_stick/dwipreproc-tmp-HLRCGM/

The structure of the resulting fielmapfile is simmilar:

$ mrinfo field_map.nii.gz 
************************************************
Image:               "field_map.nii.gz"
************************************************
  Dimensions:        88 x 88 x 60 x 2
  Voxel size:        2.5 x 2.5 x 2.5 x 1
  Data strides:      [ -1 2 3 4 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:                0.993      0.0477      0.1081      -114.5
                         -0.01882       0.967     -0.2539      -107.7
                          -0.1166      0.2501      0.9612      -96.71
  comments:          6.0.0

Again one fielmap Image and one low res T2 weighted image.

In the end I thought that an additional source of error could be the new FSL 6.0.0 I have installed- anyone else out there with feedback for that?

I am very grateful for any help or a quick fix (small alterations in the code to filter out the nasty second volume in the field_map.nii.gz), espeacially with the ISMRM deadline approaching.

Thank you so much! :grinning: @rsmith @jdtournier


#4

@XA10 / XA11A DICOMs / Slice Timing

The dcm2niix community has figured out how to optain proper slice timing info from XA10 :grinning:: https://github.com/rordenlab/dcm2niix/issues/240
XA11A data looks good too!
(Weĺl get XA11B in 2 weeks … Siemens promised not to change anything in the DICOM-subsystem … I`ll keept you updated as it could, maybe, at sometime, possibly not to far in the furture, get out for Prisma too.)
(btw: All challanges with DICOM set aside – data quality itself is excellent, especially in the neck with the 3rd order local Coil Shim)

@dwipreprocess

Any ideas for a quick fix? (ducktape & WD40 style and any tricks fully accepted :wink:)

Would safe my day :grinning:

Thank you all for making this great community and software possible!


#5

p.s.: using the slice timing from dcm2niix converted with the script on the eddy user page works great even for terrible head motion with parameters

-eddy_options "--niter=8 --fwhm=10,6,4,2,0,0,0,0 --repol --ol_type=both --mporder=8 --s2v_niter=8 --slspec=AP_slspec.txt"

If only the AP PA combination would work … :roll_eyes:


#6

Thanks for confirming. I guess a proper fix involve detecting when the fieldmap is 4D, and picking the first volume only. But to get you going, the simplest ‘fix’ might be to fix up the field map in the temporary folder:

# move to temp folder:
cd /home/ingomar/from_stick/dwipreproc-tmp-HLRCGM/
# extract first volume:
mrconvert field_map.nii.gz -coord 3 0 field_map2.nii.gz
# overwrite field map with the 3D version:
mv field_map2.nii.gz field_map.nii.gz
# return to original folder:
cd -

and then running your command again with dwipreproc’s -continue option. I think that might work as-is… :crossed_fingers: