[ENH] Add phase1/phase2 SDC support

[mattcieslak]

Changes proposed in this pull request

Add support for phase1 and phase2 fieldmaps, related to #1169. This is just one idea for how to support them, another possibility would be to make a separate workflow.

Documentation that should be reviewed

Documentation probably doesn't need to be changed.

[welcome]

Thanks for opening this pull request! We have detected this is the first time for you to contribute to fMRIPrep. Please check out our contributing guidelines.
We invite you to list yourself as a fMRIPrep contributor, so if your name is not already mentioned, please modify the .zenodo.json file with your data right above Russ' entry. Example:

{
   "name": "Contributor, New FMRIPrep",
   "affiliation": "Department of fMRI prep'ing, Open Science Made-Up University",
   "orcid": "<your id>"
},
{
   "name": "Poldrack, Russell A.",
   "affiliation": "Department of Psychology, Stanford University",
   "orcid": "0000-0001-6755-0259"
},

Of course, if you want to opt-out this time there is no problem at all with adding your name later. You will be always welcome to add it in the future whenever you feel it should be listed.

[effigies]

Thanks for this. Some initial thoughts, though again @oesteban will be the more competent reviewer, here.

[oesteban]

Please remove the [WIP] tag to indicate this is almost there, and I'll gladly review it :)

Thanks for the contribution!!

[mattcieslak]

@oesteban have you ever worked with phase1, phase2 images? I tried running some directly through the siemens2rads function, which looks like it handled these at some point, but the output looks wrong.

There is a script used here internally that handles these images, but it has two different equations depending on whether the readout was bipolar or monopolar (which isn't specified in the json sidecars). I'm wondering if supporting this kind of fieldmap is worth the trouble.

[mattcieslak]

Got it to work! The subtraction done in the nipype workflow definitely doesn't work for calculating a phasediff. Instead the math was a little more complicated (https://github.com/mattcieslak/fmriprep/blob/phase1phase2/fmriprep/interfaces/fmap.py#L578). Currently getting a data set on openneuro that has phase1/phase2 fieldmaps.

[effigies]

Will phases always range from 0 to 4096?

[mattcieslak]

We're trying to recreate these scanning sequences to see if this is still true for the current Siemens product. This has been true for all the older versions of the sequence I've run into.

[effigies]

Thanks for the diligence on this. Do you know if any other scanners produce phase1/phase2 fieldmaps? It would be good to see if this convention extends beyond Siemens.

Perhaps @neurolabusc might have the experience or know someone who does.

[neurolabusc]

I have a web page that is relevant. For Siemens images I think the scaled range is -4096..4092 (scl_slope = 2.000000; scl_inter = -4096.000000). For Philips I think the scaled range is ~-500..496 for some systems and -pi..+pi for others. For GE I think the convention is to save real/imaginary images so you would use fslcomplex to get the phase map. As a good estimate, you can assume min..max should be scaled very close to -pi..+pi.

[effigies]

Thanks for the rundown. If the slope and intercept are 2, -4096, and the original range is 0..4096, then this should be:

rad0 = np.pi * phase0 / 4096
rad1 = np.pi * phase1 / 4096

This is because img.get_fdata() already applies the scaling, so we should work from the scaled values, not the raw values.

I think we can ignore the GE case for now; that seems like it should probably not be considered "phase1"/"phase2" fieldmaps, and it's a job for BIDS to establish what the correct naming scheme is to indicate those files.

I think in the short term, we could have some a heuristic where we check for values above 500 and pi, to establish the correct denominator.

In the long term, I'm leery of putting too many heuristics for interpreting data into fMRIPrep that are not in the BIDS spec. We may want to consider a metadata field that unambiguously indicates the correct interpretation. (Though the heuristic would probably need to stay, to keep processing datasets that conform with BIDS 1.0.)

[effigies]

I see. So I guess we need to watch for cases where the scale and intercept are missing.

[mattcieslak]

We scanned all the available GRE sequences on a prisma and uploaded them on openneuro in case you want to take a look. Here's the link: https://openneuro.org/datasets/ds001600. These were converted with a recent dcm2niix.

I'm open to suggestions on how to handle these robustly. Also, sorry about the formatting mess. I ran a misconfigured yapf and am still putting things back.

[neurolabusc]

@mattcieslak can you send me the raw DICOM data and permission to distribute the DICOMs on the dcm2niix wiki? This would provide a nice reference set for all DICOM-to-BIDS conversion tools to test, and it would provide nice validation for any future changes or features.

[neurolabusc]

Despite being acquired on a E11C Prisma that supports high dynamic range (16-bit DICOM), all the samples you acquired appear to use 12-bit precision (0..4095). My guess is 12-bit is more than sufficient for the SNR provided, but I would make sure that any code generates an error if image intensities are outside this range. This will detect any sequences that provide more than 12-bit precision.

• sub-1_acq-v1_magnitude1 0..1821
• sub-1_acq-v1_magnitude2 0..1705
• sub-1_acq-v1_phase1 0..4095
• sub-1_acq-v1_phase2 0..4095
• sub-1_acq-v2_magnitude1 0..1835
• sub-1_acq-v2_magnitude2 0..1738
• sub-1_acq-v2_phase1 0..4095
• sub-1_acq-v2_phase2 0..4095
• sub-1_acq-v4_magnitude1 0..1598
• sub-1_acq-v4_magnitude2 0..1566
• sub-1_acq-v4_phasediff 0..4095
• sub-1_dir-AP_epi 0..2341
• sub-1_dir-PA_epi 0..1796

[mattcieslak]

@neurolabusc sent a link to your mailbox.sc.edu from my pennmedicine email

[neurolabusc]

@mattcieslak please check the latest dcm2niix development branch commit. It will attempt to read the CSA tag sDiffusion.dsScheme (where 1=bipolar, 2=monopolar) and if that fails it will attempt to read ucReadOutMode (where 0x1=monopolar and 0x2=bipolar). If either is found, the BIDS file will report either "DiffusionScheme": "Monopolar" or "DiffusionScheme": "Bipolar".

I believe sDiffusion.dsScheme is used on E11 for both CMRR and Siemens product. I think the ucReadOutMode was used by CMRR prior to E11. I am pretty sure that prior to E11, Siemens product diffusion sequences were always bipolar (e.g. used a bipolar twice-refocused spin echo sequence).

[mattcieslak]

@neurolabusc I added a check to see if "DiffusionScheme" is specified in the json sidecar. If it is "Bipolar", the workflow does not do SDC. If it is specified as "Monopolar" it proceeds without a warning and if it's not included in the json sidecar it will warn that it's assuming Monopolar. Does that seem reasonable @effigies ?

[neurolabusc]

@mattcieslak I think your method will work with all the Siemens images I have seen. I do not know how to detect these properties for GE and Philips. My sense is that monopolar is the default for Philips. @drmclem might know the answer.

[neurolabusc]

@mattcieslak if it is Siemens and not detected by my software, it seems that the it is most likely bipolar. I think the Siemens sequences used to be only bipolar, and I think my new version detects bipolar/monopolar settings for CMRR sequences as well as the new Siemens sequences that allow one to select between bipolar and monopolar. For many DWI scans, one should be able to reliably infer bipolar or monopolar once knows the matrix size, partial Fourier and iPAT. For most DTI sequences users tend to select the minimum possible echo time to maximize SNR. While bipolar sequences show less spatial distortion and are less sensitive to background gradients, they pay a large penalty with regards to TE. The one excpetion might be scans where users want to measure kurtosis, as these groups tend to prefer longer echo times, on the other hand they are also concerned with background gradients so the kurtosis sequences I have seen follow the trend of long TE and bipolar sequence.

[effigies]

I'm okay in either direction. Whether we reject anything that isn't explicitly monopolar or only reject things that are explicitly bipolar both seem like arguable choices. If we run into the problem that a lot of people with unlabeled bipolar fieldmaps are having issues, we can revisit it or perhaps encourage users to better label their data.

[mattcieslak]

@neurolabusc that makes sense for diffusion-weighted image series, but in this case the difference between monopolar and bipolar phase images doesn't just change the SNR/TE but the biploar readout also changes the values in the phase images. Bipolar phase images require that you scale the non-phase-encoding direction by its coordinate. This was previously corrected using

NP=`fslval ${outfroot}_rad2raw dim1`
3dcalc -a ${outfroot}_rad2raw.nii -expr "(i-$NP/2-1)*3.14159/$NP/2" -prefix ${outfroot}_rad2err.nii

which was checked against the monopolar sequence using a phantom. I have no idea how general this solution would be, which is why I don't think it should go into fmriprep. It's also probably very rare to see a bipolar GRE scan.

The phasediff calculation in this patch works easily on the monopolar sequences, but the displacements that come out of the bipolar sequence greatly shift one hemisphere (you can see the shearing here: https://github.com/mattcieslak/example_fieldmaps/blob/master/figures/sub-1_task-rest_acq-v1_sdc_phase.svg. right click and "view image" for the animation).

[neurolabusc]

OK, thanks for the clarification.

[mattcieslak]

I'm hoping to finish this up. I merged in the latest fmriprep and am running into an issue with fmriprep-docker where --patch-fmriprep results in smriprep being missing. I'm using version 1.3.2. It this expected?

[effigies]

@mattcieslak We should probably add a --patch-smriprep option to fmriprep-docker. But you can use -i poldracklab/fmriprep:1.4.0a1 to use the latest pre-release, which should have a compatible version of smriprep.

[mattcieslak]

Ok, this is working now on all my test datasets. Anything else I should add?

[oesteban]

Thanks a lot. I'll need some time to go through this PR with enough depth. Since I'm at ISMRM right now, and we are about to release a new version with big changes, the review may take a little while.

[oesteban]

I'm concerned that phase unwrapping should be performed on each phase map separately. Since we are meeting tomorrow, and considering that all the fieldmap correction infrastructure in fMRIPrep is needing a deep overhaul, I would hold off merging this PR and try to pick up as much as I can from here during the refactor of fieldmaps (which should start after OHBM, at the latest).

Of course, I'd like to add Matt's name to the .zenodo.json file ASAP because regardless of the fate of this particular PR, it is clear that the main contribution (diagnosing and drafting a solution to the problem) has already been done. Very much appreciated 👍

[oesteban]

I don't think there exists such a suffix in BIDS. This works because fmap_key gets renamed to phase1 either way (not sure about how then it digests the phase2 suffix, I'll check below).

[oesteban]

Along the lines of the previous comment.

Suggested change

	if fmap['suffix'] in ['fieldmap', 'phasediff', 'phase']:
	if fmap['suffix'] in ['fieldmap', 'phasediff', 'phase1']:

[oesteban]

Okay, I have confirmed this point: as per the FUGUE guide, each phase map should be unwrapped with PRELUDE separately (step 3) of the guide. Then they can be subtracted as they are not wrapped within [-pi, pi) anymore.

After that point, the pipeline should be exactly the same as for the phasediff type of fieldmap.

[mattcieslak]

This might also solve some of the issues we've been running into with these fieldmaps when they come from our 7T scanner. Also, I'd be happy to help with a refactor of the fieldmap workflows. I'm maintaining a diffusion pipeline that shamelessly copies them and would like to keep them up to date. Should we close this PR then?

[oesteban]

This will be continued on https://github.com/poldracklab/sdcflows