Correction to Philips fieldmap preprocessing #271

axkralj990 · 2023-07-19T15:34:28Z

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

…rrection

coalsont · 2023-07-19T23:04:37Z

global/scripts/FieldMapPreprocessingAll.sh

-		    PhaseInputName=`getopt1 "--fmapphase" $@`
-		    DeltaTE=`getopt1 "--echodiff" $@`
+	    PhaseInputName=`getopt1 "--fmapphase" $@`
+		DeltaTE=`getopt1 "--echodiff" $@`


In terms of style, these whitespace changes show that the file has a mix of tabs and spaces (even within a single line), and the current changes also leave edited lines inconsistent with their neighbors (these two changed lines use a tab and 4 spaces, and then two tabs for the next line). I suggest trying to match the indentation style of at least consecutive lines within lines you were going to edit anyway.

I changed this in my most recent commit. There were several lines in the FieldmapPreprocessingAll.sh script with such inconsistencies. I changed all the indentations in the file to 4 spaces, as this seemed as the most common pattern. There were also a few lines with 2-space indentation that I increased to 4 spaces.

…apPreprocessingAll.sh

glasserm

If this is tested and working well in your hands, I am happy to merge. Thanks for your work on this!

axkralj990 · 2023-07-23T11:48:27Z

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 · 2023-07-23T13:02:02Z

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

mharms · 2023-07-24T22:23:06Z

Couple things:

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 Determine Effective Echo Spacing for Philips rordenlab/dcm2niix#377
"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
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 · 2023-07-25T15:42:51Z

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

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.
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.
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 · 2023-07-25T16:53:14Z

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 · 2023-07-28T12:00:12Z

@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 · 2023-07-31T18:09:39Z

@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 · 2023-08-07T19:10:53Z

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.
rordenlab/dcm2niix#377 (comment)

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 · 2023-08-08T09:02:03Z

@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 · 2023-08-08T18:37:50Z

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

Aleksij Kraljic and others added 3 commits May 4, 2023 03:49

implemented the proper Philips fieldmap preprocessing and distort. co…

7bea1e4

…rrection

Merge remote-tracking branch 'upstream/master' into philips_fieldmap_fix

c34e9e0

remove debug/support statements from fieldmap scripts

0afe56d

coalsont reviewed Jul 19, 2023

View reviewed changes

fix the indentation incosistencies with tabs-spaces throughout FieldM…

cc24fd2

…apPreprocessingAll.sh

glasserm reviewed Jul 22, 2023

View reviewed changes

glasserm self-requested a review July 23, 2023 13:02

glasserm approved these changes Jul 23, 2023

View reviewed changes

coalsont merged commit 17686f3 into Washington-University:master Jul 24, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Correction to Philips fieldmap preprocessing #271

Correction to Philips fieldmap preprocessing #271

axkralj990 commented Jul 19, 2023 •

edited

Loading

coalsont Jul 19, 2023

axkralj990 Jul 20, 2023

glasserm left a comment

axkralj990 commented Jul 23, 2023

glasserm commented Jul 23, 2023

mharms commented Jul 24, 2023

axkralj990 commented Jul 25, 2023 •

edited

Loading

mharms commented Jul 25, 2023

axkralj990 commented Jul 28, 2023

axkralj990 commented Jul 31, 2023

mharms commented Aug 7, 2023

axkralj990 commented Aug 8, 2023

mharms commented Aug 8, 2023

Correction to Philips fieldmap preprocessing #271

Correction to Philips fieldmap preprocessing #271

Conversation

axkralj990 commented Jul 19, 2023 • edited Loading

Background

Main Changes

coalsont Jul 19, 2023

Choose a reason for hiding this comment

axkralj990 Jul 20, 2023

Choose a reason for hiding this comment

glasserm left a comment

Choose a reason for hiding this comment

axkralj990 commented Jul 23, 2023

glasserm commented Jul 23, 2023

mharms commented Jul 24, 2023

axkralj990 commented Jul 25, 2023 • edited Loading

mharms commented Jul 25, 2023

axkralj990 commented Jul 28, 2023

axkralj990 commented Jul 31, 2023

1. MultiBand 3, SENSE 1.9

2. MultiBand 3, SENSE 1

3. MultiBand 4, SENSE 1.9

4. MultiBand 4, SENSE 1

mharms commented Aug 7, 2023

axkralj990 commented Aug 8, 2023

1. MultiBand 3, SENSE 1.9

2. MultiBand 3, SENSE 1

3. MultiBand 4, SENSE 1.9

4. MultiBand 4, SENSE 1

mharms commented Aug 8, 2023

axkralj990 commented Jul 19, 2023 •

edited

Loading

axkralj990 commented Jul 25, 2023 •

edited

Loading