Issue with ROMEO Phase Unwrapping in SEPIA Pipeline for GE Scanner QSM Data

[LuuuXG]

Hi,
Thank you for the handy QSM preprocessing GUI! I successfully ran the pipeline and obtained an acceptable QSM chimap. However, after reviewing the discussion section on the GitHub page, I noticed the recommendation of specific parameters for QSM consensus data processing (as described here). My QSM data was acquired using a GE scanner, and I found that both the MEDI and FANSI pipelines employ the ROMEO total field calculation for echo phase combination.
Unfortunately, the ROMEO algorithm did not seem to produce the desired output in my case, specifically for the fieldmap and the subsequent images.

I have read through issue #85 and tried performing phase unwrapping with ROMEO, with both bipolar eddy correction and MCPC-3D-S set to off. However, the result still appears unreliable.
I would appreciate any insight into what might be causing this issue, as I’ve obtained good results using other methods.
I have attached an example dataset and the SEPIA log file for reference. Looking forward to your suggestions!
https://drive.google.com/file/d/1Jb2TdxUY5hC88nqMQpluTXM6nNE1EFvK/view?usp=sharing

[kschan0214]

Hi Youjie,

The issue you encountered is different from Issue #85. The artefact you observed is caused by the inter-slice phase polarity differences in the GE data (see the discontinuity between even and odd number slices):

and causes the subsequent phase unwrapping method fail.

In the consensus paper's tutorial, I did not obatin the phase data from dcm2niix on the GE data. Instead, I got only the real and imaginary images. In Step 4 of the tutorial, the interslice polarity is corrected in the preparation script. If you have the real/imaginary NIFTI data, you may try to follow the tutorial. Alternatively, you can try the attached script here:
correct_GE_interslice_phase_offset.zip which should fix your issue. In my preliminary testing, you will need to turn the "MCPC-3D-S phase offset correction" option to "On" for your data when you use ROMEO for total field calculation.

Let me know how it goes!

Kwok

[LuuuXG]

Thank you for your help! I used the script and configuration file from the zip you provided, and I successfully obtained the susceptibility map.

Before closing this issue, I have a few additional questions to clarify:

1. You mentioned "inter-slice phase polarity differences" in your previous response. Does this imply that my QSM data was acquired using a "bipolar" readout? After reading this paper and two discussions on Neurostars (topic 1, topic 2), I understand that a bipolar acquisition is somewhat similar to an EPI readout, where the direction alternates between echoes. I suspect the "inter-slice phase polarity differences" could be related to an issue during the raw data saving process, as outlined in the consensus paper. Are there any resources you would recommend for learning about this concept? I'm still unfamiliar with some basic MRI principles.
2. For my dataset, I used Optimum weights and MEDI nonlinear fit for total field recovery and phase unwrapping. The output Chimap seemed to be fine. However, after realizing the issue with discrepancy of interslice polarity, I wonder whether this discrepancy was automatically corrected by these methods, or if I may have missed something critical.
3. As demonstrated in Figures 5 and 6 of the consensus paper (and in the SEPIA 101 exercise), applying stricter masks improves the quality of QSM reconstruction. However, this also results in the erosion of brain tissue, which poses a challenge when registering to standard MNI space for whole-brain analysis. Is this a common issue, or is there a recommended strategy for handling this paradox?

Apologies for the questions beyond the scope of the original issue, but I would greatly appreciate any further clarification on these points!

[kschan0214]

Hi! I converted this issue to the Discussion as it is a more appropriate place.

To answer your questions:

1. The phase issue you encountered is more related to the vendor than the way how the data was collected. This issue exists when using dcm2niix to convert GE scanner data into NIFTI format. Somehow the data in the odd-numbered slices are missing a negative sign vs the even-numbered slices (or vice versa). I haven't looked into whether this is truly the way the DICOM images are stored in GE or if there is any compatibility issue on the dcm2niix side. Regarding the bipolar readout, its effect is more on the time domain of the data. When you look into the plot of the phase images in time series (e.g. loading the 4D images into fsleyes), you could observe a zig-zag pattern between odd and even numbered echoes which are often absent in the monopolar data. It is because the eddy current induced by the bipolar readout gradient introduced an offset between successive echoes. You may check out this paper for more information. You can also download the example data from the consensus paper which contains both monopolar and bipolar images.

2. The short answer is yes. The principle of both methods is to first perform temporal phase unwrapping, followed by spatial unwrapping. In this case, the negative sign is handled by the temporal phase unwrapping step and does not cause any issue in the subsequent processing. For multi-echo data, ROMEO takes both spatial and temporal phase coherence into account. Therefore, the inconsistent phase error between slices (i.e. spatial phase) will produce incorrect output.

3. Unfortunately, I don't have a good solution for that (let me know if you find one!). For volumetric analysis, I usually use the magnitude of the GRE data for registration and apply the transformation matrix to bring the Chimap to the target space. While this ensures a good registration result, the eroded voxels close to the pial surface will still be missing. One way to potentially recover the missing cortical voxels would be using a two-pass pipeline. You can generate the Chimap using a tight mask in the first run and then use a loose mask in the second run. You can then combine the two Chimaps retrospectively to fill out the missing voxels. This would be similar to the concept of Star-QSM and QSMxT. In this way, at least you can be sure the QSM recon inside the brain is good and you also get the missing voxels back (but those voxels are less reliable).

[kschan0214]

Forgot to mention...you might also find the mask section (Section 5) of the consensus paper relevant for cortical analysis.

[LuuuXG]

Thank you so much for the detailed response!