Correction to Philips fieldmap preprocessing

[axkralj990]

Background

In this pull request, I corrected the issues from my previous pull request, #194, regarding Philips fieldmap image preprocessing. Upon further investigation, I discovered that I had made an incorrect assumption about the units of the input fieldmap images.

I had incorrectly assumed that the fieldmap image produced by the Philips Achieva 3.0T TX was a phase difference image in degrees. This assumption was based on the voxel values in the image, which ranged from about -180 to 180. However, upon closer inspection, the values were between about -200 and 200 (our actual dTE obtained from PAR/REC was 2.45 ms). To clarify the units, I contacted a Philips physicist and received confirmation that the fieldmap image was indeed in Hz.

The preprocessing steps are shown in the figure below, along with the old steps and current steps, in addition to the established Siemens processing, for reference.

Main Changes

In order to properly process the fieldmap, I have now corrected the previous steps. Since the input image is not a phase difference image, as previously assumed, it now needs to be first converted to a phase difference map, then unwrapped using FSL prelude, and finally converted back to a fieldmap.

I have also included a conditional check to determine if the echo difference parameter (dTE) has been specified. If it is not, the phase unwrapping is skipped. However, based on my research, I believe including the option for phase unwrapping is still useful. Some researchers skip phase unwrapping, especially when using the default echo difference of 3 ms during acquisition. In such cases, wrapping tends to be minimal and limited to the inferior-frontal region where BOLD signal loss usually occurs. When a larger dTE is set, the unwrapping option is necessary, as the phase is significantly wrapped (tested with dTE up to 10 ms).

Another factor contributing to the misconception that the earlier code was correctly preprocessing the fieldmaps was that the combination of all operations applied to the input image produced a similar scaling factor. Specifically, the input data was previously multiplied by pi/180 * 1000/(3 ms) = 5.81, whereas in the corrected code, the input data is multiplied by 2*pi = 6.28. The old cold produces similar values, but definitely lower than they should be, leading to under-correction.

The figure also shows an unwrapped fieldmap and three sagittal slices of the distortion-corrected BOLD image with the white matter surface outline. The effective echo spacing parameter for the correction was calculated as follows (Neurostars forum discussion):

effective_echo_spacing = (1000 * water_fat_shift)/(434.215 * (epi_factor+1))

I have made the necessary changes in the pull request to fix these issues and ensure proper distortion correction. Please let me know if there are any concerns or additional requests. Thank you.

Edit: Figure updated on August 8, 2023

[axkralj990]

I'm glad I could help. The correction works well on my end, and I also tested it yesterday on data collected with the Philips Achieva 3.0T after a dStream upgrade. I used the effective echo spacing estimate returned by the dcm2niix toolbox (EstimatedEffectiveEchoSpacing), and the correction results are practically the same as the correction with the topup method.

[glasserm]

I later saw you scene and I agree it looks good.

[mharms]

Couple things:

1. Is this particular approach for correction dependent on getting "echo spacing" correct? Because the issue of the ES for Philips is very complicated. See https://github.com/rordenlab/dcm2niix/tree/master/Philips and https://github.com/rordenlab/dcm2niix/issues/377
2. "real fieldmap" is probably not the best terminology to use for the input in the Siemens branch of your figure. It is actually a "phase diff(erence)" map, in Siemens arbitrary -4096:4096 convention
3. If the input for the Philips branch is a "real fieldmap (Hz)", then the range of that is not going to be [-500:500]/dTE, but rather whatever the actual field discrepancy (from B0) is in Hz, right?

[axkralj990]

@mharms, thank you for your feedback and important observations. Here are my thoughts; please correct me if my reasoning is incorrect in any way.

1. This is a valid point that I have thought about a lot. The success of this approach depends directly on using a correct value for the effective echo spacing. I am indeed familiar with the discussions to which you referred. After trying different formulas for calculating effective echo spacing (ees) on different EPI sequences (with different SENSE in-plane acceleration settings), I found that distortion correction using the ees estimates from the equations implemented in the dicm2niix toolbox gave very similar results to topup correction. On this basis, I assumed that the estimates were reasonable by simply comparing them to topup-corrected data across different sequences rather than trying to decipher the vendor-specific truth in the dicom files. To test this, I had access to a dataset collected on a Philips Achieva 3.0T before a dStream upgrade and to a few datasets collected on the same scanner after the dStream upgrade, where we collected data with different multiband-SENSE combinations. My goal was to analyze if different SENSE settings still lead to reasonable ees parameter estimates and consequently correct distortion correction. I think it is important to emphasize that future users perform a similar comparison to ensure that the data is properly corrected instead of simply assuming that the reported estimates for ees and total readout time are correct. If you have suggestions to further test the validity of the correction and estimates, I would be happy to explore them with our Philips data. We have collected some data sets (various short EPI sequences with both B0 maps and SE pairs) specifically for this purpose.

2. This is indeed a discrepancy that I have not changed since I first created the figure in 2020. I have updated the figure in my comment above. The input for the prelude must be a phase difference map anyway, so the figure might definitely confuse someone.

3. I assumed that at first as well, but the presence of wrapping indicates that the data must be bounded. Then I learned from the Philips physicist that the scanner converts a phase difference image into a fieldmap by dividing the data by the echo difference. The range of values shown in the figure spans one reciprocal of the echo difference. After the conversion, this leads to a phase difference image with values ranging from -pi to pi, where the wraps should occur. Philips apparently decided to omit phase unwrapping due to low wrapping when the default echo difference is used. I decided to take a more conservative approach and added an option to perform the unwrapping, which is especially critical in cases where the echo difference is increased from the default 3 ms, resulting in heavy wrapping. For example, I acquired a B0 map scan with an echo difference set to 10 ms, and the image had multiple wraps making it unusable for distortion correction without unwrapping.

[mharms]

Re (1): Did the EES (EffectiveEchoSpacing) value reported by dcm2niix change by the SENSE (in-plane) factor as you changed the SENSE factor? Can you document here which software version of Philips was used for the scanning and which dcm2niix version you used?

Re (2-3): In that case, might also be better to call the first two stages of the Philips branch of the figure a "(possibly) wrapped fieldmap (Hz)" rather than a "real fieldmap (Hz)" since a "real fieldmap" should not have any wrapping (at least for how I think about that phrasing).

[axkralj990]

@mharms, I have updated the figure so that it alerts users that we are not 100% sure of the specific image formats.

The software release is 5.7.1.2, and the scanner is a 2011 model. I have also included this information in the figure.

I am currently out of the office and will get back to you on Monday (July 31, 2023) on how SENSE factor changes affect dicm2niix effective echo spacing estimates, as I now only have access to images that vary in other parameters as well. I will take a few short EPI scans on a phantom with different SENSE values.

[axkralj990]

@mharms. Today, we did two pairs of BOLD scans on a phantom to see how the SENSE acceleration affects the estimated effective echo spacing.

If I use the equation (1000 * water_fat_shift)/(434.215 * (etl+1)), it appears that I need to divide the number by the in-plane acceleration factor (SENSE) to get values comparable to those calculated by dicm2niix. So in this case SENSE has an effect on the effective echo spacing.

The distortion correction also proved to be appropriate when evaluated on real data, one of which was acquired with SENSE = 1.9 and the other with SENSE = 1, both using the estimate from dicm2niix.

If you have any further questions, please let me know.

Here are the summary results of today's data.

1. MultiBand 3, SENSE 1.9

"EchoTime": 0.030001, "RepetitionTime": 1.3, "EchoTrainLength": 51, "WaterFatShift": 13.7557, "EstimatedEffectiveEchoSpacing": 0.000326937, "EstimatedTotalReadoutTime": 0.031059

ees = (1000 * water_fat_shift)/(434.215 * (etl+1))/SENSE
    = (1000 * 13.7557)/(434.215 * (51+1))/1.9
    = 0.32064 s
    = 0.00032064 ms

2. MultiBand 3, SENSE 1

"EchoTime": 0.031024, "RepetitionTime": 1.54233, "EchoTrainLength": 95, "WaterFatShift": 25.4869, "EstimatedEffectiveEchoSpacing": 0.0006112, "EstimatedTotalReadoutTime": 0.058064

ees = (1000 * water_fat_shift)/(434.215 * (etl+1))/SENSE
    = (1000 * 25.4869)/(434.215 * (95+1))/1
    = 0.6114 s
    = 0.0006114 ms

3. MultiBand 4, SENSE 1.9

"EchoTime": 0.03, "RepetitionTime": 1, "EchoTrainLength": 51, "WaterFatShift": 13.7557, "EstimatedEffectiveEchoSpacing": 0.000326937, "EstimatedTotalReadoutTime": 0.031059

ees = (1000 * water_fat_shift)/(434.215 * (etl+1))/SENSE
    = (1000 * 13.7557)/(434.215 * (51+1))/1.9
    = 0.32064 s
    = 0.00032064 ms

4. MultiBand 4, SENSE 1

"EchoTime": 0.03111, "RepetitionTime": 1.18339, "EchoTrainLength": 95, "WaterFatShift": 25.4869, "EstimatedEffectiveEchoSpacing": 0.0006112, "EstimatedTotalReadoutTime": 0.058064

ees = (1000 * water_fat_shift)/(434.215 * (etl+1))/SENSE
    = (1000 * 25.4869)/(434.215 * (95+1))/1
    = 0.6114 s
    = 0.0006114 ms

[mharms]

In the interest of being very specific, I think that the relevant equations used currently by dcm2niix are at the following link. Note that they use "EPI_Factor" (2001,1013) rather than "ETL", and the two are not always the same, which is an important consideration. https://github.com/rordenlab/dcm2niix/issues/377#issuecomment-1328530596

If you were to use the formulas at the above link, you should be able to duplicate the values returned by dcm2niix exactly (rather than just very close).

[axkralj990]

@mharms, thank you for pointing me to the actual equations implemented in the dcm2niix. Here are the reproduced effective echo spacing estimates.

Equations for estimating the Philips effective echo spacing parameter from dcm2niix:

ActualEchoSpacing = WaterFatShift / (ImagingFrequency * 3.4 * (EPI_Factor + 1)) 
TotalReadoutTIme = ActualEchoSpacing * EPI_Factor
EffectiveEchoSpacing = TotalReadoutTime / (ReconMatrixPE - 1)

1. MultiBand 3, SENSE 1.9

(2001, 1013) [EPI Factor]                        SL: 51
(2001, 1022) [Water Fat Shift]                   FL: 13.755722045898438
(2001, 1083) [Imaging Frequency]                 DS: '127.756744'
"ReconMatrixPE": 96                              # from dcm2niix .json

WaterFatShift = 13.755722045898438;
ImagingFrequency = 127.756744;
EPI_Factor = 51;
ReconMatrixPE = 96;

EffectiveEchoSpacing = 0.000326937
# dcm2niix "EstimatedEffectiveEchoSpacing": 0.000326937

2. MultiBand 3, SENSE 1

(2001, 1013) [EPI Factor]                        SL: 95
(2001, 1022) [Water Fat Shift]                   FL: 25.48691749572754
(2001, 1083) [Imaging Frequency]                 DS: '127.756744'
"ReconMatrixPE": 96                              # from dcm2niix .json

WaterFatShift = 25.48691749572754;
ImagingFrequency = 127.756744;
EPI_Factor = 95;
ReconMatrixPE = 96;

EffectiveEchoSpacing = 0.0006112
# dcm2niix "EstimatedEffectiveEchoSpacing": 0.0006112

3. MultiBand 4, SENSE 1.9

(2001, 1013) [EPI Factor]                        SL: 51
(2001, 1022) [Water Fat Shift]                   FL: 13.755722045898438
(2001, 1083) [Imaging Frequency]                 DS: '127.756744'
"ReconMatrixPE": 96                              # from dcm2niix .json

WaterFatShift = 13.755722045898438;
ImagingFrequency = 127.756744;
EPI_Factor = 51;
ReconMatrixPE = 96;

EffectiveEchoSpacing = 0.000326937
# dcm2niix "EstimatedEffectiveEchoSpacing": 0.000326937

4. MultiBand 4, SENSE 1

(2001, 1013) [EPI Factor]                        SL: 95
(2001, 1022) [Water Fat Shift]                   FL: 25.486919403076172
(2001, 1083) [Imaging Frequency]                 DS: '127.756744'
"ReconMatrixPE": 96                              # from dcm2niix .json

WaterFatShift = 25.486919403076172;
ImagingFrequency = 127.756744;
EPI_Factor = 95;
ReconMatrixPE = 96;

EffectiveEchoSpacing = 0.0006112
# dcm2niix "EstimatedEffectiveEchoSpacing": 0.0006112

---

The in-plane acceleration is accounted for in the EPI factor. For example, in both cases with SENSE = 1.9), I had to subtract 1 from the EPI factor reported in DICOM to get the exact acceleration factor for rescaling the actual echo spacing:

(ReconMatrixPE - 1)/(EPI_Factor - 1) = 95/50 = 1.9

[mharms]

Good to see that you can indeed (exactly) match the current output of dcm2niix.