Susceptibility Distortion Correction - effective echo-spacing and total-readout time for Philips data

[treanus]

Dear all,

First, thank you for providing fmriprep - very useful!

I have questions about ESS and TRT for SDC of bold fMRI with a B0 fieldmap using Philips data.

In the documentation of fmriprep SDC and source code of fmriprep.interfaces.fmap there are the functions get_ees and get_trt. For Philips data it seems to be enough to define the WaterFatShift, MagneticFieldStrength, PhaseEncodingDirection and ParallelReductionFactorInPlane to calculate ess and trt.

First question:
a) does this signify that the EffectiveEchoSpacing does not need to be specified in the BIDS json file of the bold fmri, when WaterFatShift, MagneticFieldStrength, PhaseEncodingDirection and ParallelReductionFactorInPlane are given in the json file?
b) this then gives a warning in the bids-validation step ("You should define 'EffectiveEchoSpacing' for this file. If you don't provide this information field map correction will not be possible.") that can be safely ignored?

Second Question:
a) looking at the source code:

-
if wfs is not None:
        fstrength = in_meta['MagneticFieldStrength']
        wfd_ppm = 3.4  # water-fat diff in ppm
        g_ratio_mhz_t = 42.57  # gyromagnetic ratio for proton (1H) in MHz/T
        wfs_hz = fstrength * wfd_ppm * g_ratio_mhz_t
        return wfs / (wfs_hz * etl)
-

it seems that ParallelReductionFactorInPlane (the sense acceleration) is not used?

b) this ParallelReductionFactorInPlane is not used, since the echo-train-length (etl) in Philips images is already taking this into account? Correct?

c) I'm confused by other information on the web (among others: http://www.spinozacentre.nl/wiki/index.php/NeuroWiki:Current_developments#B0_correction) that also provide a calculation for the ESS as:
effective echo spacing = (((1000 * wfs)/(434.215 * (EPI factor+1))/acceleration)

Thus here 'acceleration' is used. Acceleration means the sense factor, which is also the same as the ParallelReductionFactorInPlane? Correct?
I guess it all depends on how 'etl' is defined. In the philips definition it already incorporates the ParallelReductionFactorInPlane (etl = NumberOfPhaseEncodingSteps/ParallelReductionFactorInPlane)

Last Question:

Starting from Philips dicom, I need to extract the correct information to put into the BIDS json file.

Using matlab dicominfo to extract relevant information out of Philips Dicom, one finds the following:

• MagneticFieldStrength in 'MagneticFieldStrength'
• etl in 'Private_2001_1013'
• WaterFatShift in 'Private_2001_1022'
• ParallelReductionFactorInPlane in 'Private_2005_140f.Item_1.ParallelReductionFactorInPlane';

EffectiveEchoSpacing is not in the philips dicom header, and is thus calculated.

Hope someone can help with (part of) my questions...

Best regards,

Stefan

Ps. Slightly of topic:

The PhaseEncodingDirection is not fully specified in the philips dicomheader? One cannot find e.g. 'j' or 'j-' as both are defined as 'COL' in the philips dicom header?
Also the slice timing info is not in the philips dicom header?
However ParallelReductionFactorOutOfPlane is specified in the philips dicom header, and this corresponds to the multiband factor?

[oesteban]

First question:

a) does this signify that the EffectiveEchoSpacing does not need to be specified in the BIDS json file of the bold fmri, when WaterFatShift, MagneticFieldStrength, PhaseEncodingDirection and ParallelReductionFactorInPlane are given in the json file?
b) this then gives a warning in the bids-validation step ("You should define 'EffectiveEchoSpacing' for this file. If you don't provide this information field map correction will not be possible.") that can be safely ignored?

PhaseEncodingDirection is mandatory for any type of field map correction. Only the experimental "fieldmap-less" (--use-syn command line argument) would work, but assuming that the PE is AP or PA.

The WaterFatShift and MagneticFieldStrength are vendor fields of DICOMS generated by Philips. This means that they are completely discouraged. If you can calculate the EffectiveEchoSpacing yourself, I'd say that would be the right path to choose. They are implemented in fmriprep as a back-up plan. This also replies to b), as we haven't been able to test the private fields thoroughly - not very safe.

That said, if you have Philips data and you want to help us make FMRIPREP robust to these details, we'll be more than happy to try your data ourselves and make the necessary improvements.

Second Question:

a) (...) it seems that ParallelReductionFactorInPlane (the sense acceleration) is not used?
b) this ParallelReductionFactorInPlane is not used, since the echo-train-length (etl) in Philips images is already taking this into account? Correct?

As you can see:

https://github.com/poldracklab/fmriprep/blob/260872273a1f4ef02de2cae20dd7d6948b531c4b/fmriprep/interfaces/fmap.py#L326

the ParallelReductionFactorInPlane is indeed taken into account.

c) I'm confused by other information on the web (among others: http://www.spinozacentre.nl/wiki/index.php/NeuroWiki:Current_developments#B0_correction) that also provide a calculation for the ESS as:
effective echo spacing = (((1000 * wfs)/(434.215 * (EPI factor+1))/acceleration)
Thus here 'acceleration' is used. Acceleration means the sense factor, which is also the same as the ParallelReductionFactorInPlane? Correct?

Correct.

I guess it all depends on how 'etl' is defined. In the philips definition it already incorporates the ParallelReductionFactorInPlane (etl = NumberOfPhaseEncodingSteps/ParallelReductionFactorInPlane)

Indeed:
https://github.com/poldracklab/fmriprep/blob/260872273a1f4ef02de2cae20dd7d6948b531c4b/fmriprep/interfaces/fmap.py#L328

Last Question:

Starting from Philips dicom, I need to extract the correct information to put into the BIDS json file.
Using matlab dicominfo to extract relevant information out of Philips Dicom, one finds the following:
MagneticFieldStrength in 'MagneticFieldStrength'
etl in 'Private_2001_1013'
WaterFatShift in 'Private_2001_1022'
ParallelReductionFactorInPlane in 'Private_2005_140f.Item_1.ParallelReductionFactorInPlane';
EffectiveEchoSpacing is not in the philips dicom header, and is thus calculated.

I added some pertinent comments to the current BIDS draft about this in this section. Meaning, there are no clear guidelines about how these private dicom fields should be made BIDS-compatible. And, as I said before, we are not even sure about our own implementation. But we would be happy to debug these fieldmaps if you have the time for that.

[treanus]

Thank you @oesteban for your reply.

If you can calculate the EffectiveEchoSpacing yourself, I'd say that would be the right path to choose.

This is exactly what I do. I use some matlab code, based on what is provided on http://support.brainvoyager.com/functional-analysis-preparation/27-pre-processing/459-epi-distortion-correction-echo-spacing.html (and is identical to the code in fmriprep/fmriprep/interfaces/fmap.py). Then I use the jooh fork of dcm2bids to put the ess and trt (and other data) in the json file. I was just confused about the correct formula.

the ParallelReductionFactorInPlane is indeed taken into account.

Indeed - sorry for my oversight. And npe is defined as npe = nb.load(in_file).shape[_get_pe_index(in_meta)]. Sorry, for asking once more, but I want to be sure: npe is the "PhaseEncodingSteps" in the json header?

But we would be happy to debug these fieldmaps if you have the time for that.

How can we get in touch?

[oesteban]

npe is the "PhaseEncodingSteps" in the json header?

Correct. Which is the number of voxels along the PE direction on the image matrix (without taking into account acceleration schemes)

How can we get in touch?

If you can share some of your data (e.g. one subject), you could upload it to OpenNeuro.org. Upload will require you valid BIDS data. Once uploaded, you'll see a large catalogue of applications available to you. One of them is the latest release of FMRIPREP. You can run it at no cost.

That way you can anticipate how the performance of FMRIPREP will be like on your data.

Otherwise, we are always open here or in neurostars.org to investigate the problems you may hit. Sometimes, we will ask you for some data to replicate the problem.

Finally, we always welcome new contributors. So if you find something that needs to be fixed, please do not hesitate to send us a pull request.

[treanus]

Looking again at the source code and at http://dbic.dartmouth.edu/pipermail/mrusers/attachments/20141112/eb1d20e6/attachment.pdf (from http://fmriprep.readthedocs.io/en/latest/sdc.html).

if wfs is not None:
        fstrength = in_meta['MagneticFieldStrength']
        wfd_ppm = 3.4  # water-fat diff in ppm
        g_ratio_mhz_t = 42.57  # gyromagnetic ratio for proton (1H) in MHz/T
        wfs_hz = fstrength * wfd_ppm * g_ratio_mhz_t
        return wfs / (wfs_hz * etl)

Shouldn't the last line be?

return wfs / (wfs_hz * (etl + 1) )

[oesteban]

One more interesting resource (recommended after a revision of the FSL FUGUE user guide) when we check on this - https://osf.io/hks7x/?show=view

[effigies]

This seems like the relevant issue, but I can open a new one if it's sufficiently different. I've recently converted some Philips phase1/phase2 fieldmaps with dcm2niix and gotten the following sidecars:

{
        "Manufacturer": "Philips",
        "PatientPosition": "HFS",
        "SeriesDescription": "ImageMRSERIES",
        "ProtocolName": "FieldMapping",
        "SeriesNumber": 16,
        "AcquisitionNumber": 16,
        "ImageComments": "Audio",
        "PhilipsRescaleSlope": 1.53455,
        "PhilipsRescaleIntercept": -3142,
        "PhilipsScaleSlope": 651.74,
        "UsePhilipsFloatNotDisplayScaling": 1,
        "EchoNumber": 1,
        "EchoTime": 0.0023,
        "RepetitionTime": 0.02,
        "ReconMatrixPE": 256,
        "ImageOrientationPatientDICOM": [
                1,
                0,
                0,
                0,
                0.996551,
                -0.0829821      ],
        "ConversionSoftware": "dcm2niix",
        "ConversionSoftwareVersion": "v1.0.20190902"
}

{
        "Manufacturer": "Philips",
        "PatientPosition": "HFS",
        "SeriesDescription": "ImageMRSERIES",
        "ProtocolName": "FieldMapping",
        "SeriesNumber": 16,
        "AcquisitionNumber": 16,
        "ImageComments": "Audio",
        "PhilipsRescaleSlope": 1.53455,
        "PhilipsRescaleIntercept": -3142,
        "PhilipsScaleSlope": 651.74,
        "UsePhilipsFloatNotDisplayScaling": 1,
        "EchoNumber": 2,
        "EchoTime": 0.0046,
        "RepetitionTime": 0.02,
        "ReconMatrixPE": 256,
        "ImageOrientationPatientDICOM": [
                1,
                0,
                0,
                0,
                0.996551,
                -0.0829821      ],
        "ConversionSoftware": "dcm2niix",
        "ConversionSoftwareVersion": "v1.0.20190902"
}

None of the methods for recovering EES seem to work with this metadata.

I'll also ping @neurolabusc in case there's something I need to do to finagle dcm2niix to retrieve the needed metadata.

For reference, the current get_ees code is:

https://github.com/poldracklab/sdcflows/blob/03ce021a792ee5a9ba23f99c7f9ddc3afcd494f2/sdcflows/interfaces/fmap.py#L382-L462

Happy to share data privately, but waiting for IRB approval to make public.

[neurolabusc]

Chris,

1. Your minimal BIDS sidecar for the phase map reflects the fact that Philips PAR/REC files are devoid of virtually all meta information. This is a limitation of the input data not dcm2niix.

I would urge all Philips users to save DICOM data direct from the console, the PAR/REC format is too lean to be considered archival quality. As @drmclem noted: the PAR-REC format was built for simpler times and its always been a tension between simplicity and comprehensiveness of the format.

2. With regards to the @treanus comments, it is worth noting that Philips DICOMs have evolved a lot over the years. I have tried to get as much information as possible, and some details (like phase encoding polarity and slice timing) are simply not in the Philips DICOM. I do not have access to Philips hardware, but if someone wants to develop a validation dataset with modern images that validate robust effective readout time estimates, I am happy to extend dcm2niix. As has been noted before, the "correct" value for TotalReadoutTime is only important in the context of FSL's topup if you need/want calibrated field maps, OR if for some reason your two acquisitions with differing polarities have differing TotalReadoutTimes. In most situations you can supply a bogus value, and the images used for subsequent processing are identical to if you had the precise values.

[effigies]

Thanks for weighing in, Chris. I need to integrate all of this a bit. Hopefully these fieldmaps are salvageable, but the ship has sailed in terms of getting DICOMs with the metadata we're currently expecting.

[neurolabusc]

@treanus I have updated dcm2niix to use the DICOM tags to populate the following BIDS tags:

	"ParallelReductionFactorInPlane": 2.5,
	"WaterFatShift": 7.95575,
	"EffectiveEchoSpacing": 0.000591038,

I used these samples to develop this. Since I do not have Philips data, I would appreciate if others would test this. Also, tell me if (wfs_hz * etl) or (wfs_hz * etl+1) is preferred. As noted, EffectiveEchoSpacing typically has no impact, but it is nice to include for completeness if we can.

Note that until Philips specifies the phase encoding polarity, dcm2niix will only be able to populate PhaseEncodingAxis (e.g. i or j)rather than PhaseEncodingDirection (disambiguating j from j-).

[oesteban]

Note that until Philips specifies the phase encoding polarity, dcm2niix will only be able to populate PhaseEncodingAxis (e.g. i or j)rather than PhaseEncodingDirection (disambiguating j from j-).

wow, that's really unfortunate!

[neurolabusc]

We do not have Philips hardware at my center, but it might be nice for someone with access to Philips equipment (@baxpr?) to acquire a validation dataset similar to what @mharms acquired on for Siemens, see his Excel spreadsheet for more details.

Glancing at your formula, I wonder if partial Fourier is handled correctly. As I recall, when calculating effective echo spacing using echo train length you ignore in plane acceleration (as that covers all of K-space, e.g. with SENSE=2 a 80 row scan with ETL=40 has the same distortion as an unaccelerated 40 row scan with the same FOV), but you must adjust based on partial Fourier (e.g. for a pF=5/8 with a 80-row scan and ETL=50 one has the same distortion as an unaccelerated 80 row scan).

Therefore, I think your formula needs to use the ParallelReductionFactorInPlane to determine if the scan has partial Fourier (or have some other measure to detect pF) and adjust your calculation.

Again, since in typical cases effective echo spacing has no impact on subsequent processing, one may ignore this. However, I hate adding new features to my software when I know that they make assumptions that are not always correct.

[baxpr]

@neurolabusc this question has come up here, too. I'll add this to my list - hopefully we can get some useful info or test scans.

[neurolabusc]

@treanus, @baxpr, @effigies and @oesteban I would be grateful if you could look at the formula and Excel spreadsheet provided here.

@baxpr, no need for test scans, as I have some from two sites, with permission to share from at least one site.

[oesteban]

Addressed in #132 - accepting patches on the dev/1.4.0 branch.