Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

mrconvert: can not determine phase encoding direction from GE SIGNA Premier v29.0 DTI #2517

Open
blezek opened this issue Oct 17, 2022 · 22 comments
Labels

Comments

@blezek
Copy link
Contributor

blezek commented Oct 17, 2022

Describe the bug

mrconvert (and other mrtrix tools) do not identify the phase encoding direction on GE 3T SIGNA Premier imaging version v29.0. Results in incorrect dec images (and really bad tractography!).

To Reproduce

mrinfo -export_pe_table freq-encoding-ap.txt MR-INNOVATION-20221014/Ax-DTI-70-dir-default-Freq-A-P-2/
mrinfo: [done] scanning DICOM folder "MR-INNOVAT...-70-dir-default-Freq-A-P-2/"
mrinfo: [100%] reading DICOM series "Ax DTI (70 dir-default)Freq A/P"
mrinfo: [ERROR] no phase-encoding information found within image "<REDACTED> [MR] Ax DTI (70 dir-default)Freq A/P"

mrinfo -export_pe_table freq-encoding-lr.txt MR-INNOVATION-20221014/Ax-DTI-70-dir-default-Freq-R-L-5/
mrinfo: [done] scanning DICOM folder "MR-INNOVAT...-70-dir-default-Freq-R-L-5/"
mrinfo: [100%] reading DICOM series "Ax DTI (70 dir-default)Freq R/L"
mrinfo: [ERROR] no phase-encoding information found within image "<REDACTED> [MR] Ax DTI (70 dir-default)Freq R/L"

Output of dcminfo -all

these are randomly selected slices from the series, but all images in the respective series have the same phase / frequency encoding directions.

Platform/Environment/Version

OS:

LSB Version:	:core-4.1-amd64:core-4.1-noarch
Distributor ID:	CentOS
Description:	CentOS Linux release 7.9.2009 (Core)
Release:	7.9.2009
Codename:	Core

MRtrix3 version (example: mrinfo -version)

== mrinfo 3.0_RC2 ==
64 bit release version, built Mar 16 2022, using Eigen 3.3.7
Author(s): J-Donald Tournier (d.tournier@brain.org.au) and Robert E. Smith (robert.smith@florey.edu.au)
Copyright (c) 2008-2017 the MRtrix3 contributors.
@blezek blezek added the bug label Oct 17, 2022
@bjeurissen bjeurissen reopened this Oct 18, 2022
@bjeurissen
Copy link
Member

Thanks for reporting this issue. I don’t think it really qualifies as a bug. None of the issues with wrong dec or tractography are directly related to missing PE info, they are the result of displaying the dec within a viewer that uses different coordinate conventions. How those different conventions manifest in this particular scenario happens to coincide with successful extraction of the PE direction.

That said, I think your DICOMs might be of use to help us extract the PE directions in currently unsupported scenarios.

@blezek
Copy link
Contributor Author

blezek commented Oct 18, 2022

@bjeurissen based on your suggestions I did some further investigation visualizing the tensors acquired by the scanner. Agree that this is not a bug in mrtrix, rather something funny is going on, likely on our side. I've engaged our site MR application scientist to help, will see what he comes up with. In the meantime, I processed the raw DICOM using this code for both frequency encoding directions (L/R and A/P):

  mrconvert -force "$dicom" dti-$suffix.mif
  dwi2mask -force dti-$suffix.mif mask-$suffix.mif
  dwidenoise -force -mask mask-$suffix.mif dti-$suffix.mif dti_denoised-$suffix.mif
  dwiextract -force dti-$suffix.mif - -bzero | mrmath -quiet -axis 3 -force - mean b0-$suffix.mif

  dwi2tensor -force -mask mask-$suffix.mif dti-$suffix.mif dt-$suffix.mif
  tckgen \
    -force \
    -nthreads 16 \
    -algorithm Tensor_Prob \
    -select "$NSTREAMLINES" \
    -minlength 60 \
    -seed_image mask-$suffix.mif \
    -mask mask-$suffix.mif \
    dti_denoised-$suffix.mif \
    streamlines-$suffix.tck

I then displayed using mrview:

mrview b0-ap.mif -odf.load_tensor dt-ap.mif -tractography.load streamlines-ap.tck &
mrview b0-lr.mif -odf.load_tensor dt-lr.mif  -tractography.load streamlines-lr.tck &

The tensor directions appear flipped in plane for A/P frequency encoding, while the L/R frequency encoding are as expected.

2022-10-18_08-28-19

The tractography likewise looks correct for L/R frequency encoding.

2022-10-18_08-36-38

mrtrix reports the gradient directions are the same in both acquisitions (I can post if it would be helpful).

@bjeurissen
Copy link
Member

@bjeurissen based on your suggestions I did some further investigation visualizing the tensors acquired by the scanner. Agree that this is not a bug in mrtrix, rather something funny is going on, likely on our side.

Based on what you posted above I would say it is a bug. Your DICOMs might be of interest to @jdtournier to fix edge cases in the DICOM conversion.

@blezek
Copy link
Contributor Author

blezek commented Oct 25, 2022

The dcm2niix wiki has a bit more information (see https://github.com/rordenlab/dcm2niix/tree/master/GE#diffusion-tensor-notes). Because GE encodes directions in "physics" orientation, the DICOM tag Phase Encoding Direction (0018,1312) must be used to convert to scanner or physical orientation. I had a look through the DICOM conversion code, but couldn't figure out where the reading was happening.

I do have data that I could share with phase encode and L/R and A/P directions from the same subject on two different GE scanners.

For the moment, I'll convert using dcm2niix, then use mrconvert to get into mif format and all should be well.

@bjeurissen
Copy link
Member

The dcm2niix wiki has a bit more information (see https://github.com/rordenlab/dcm2niix/tree/master/GE#diffusion-tensor-notes). Because GE encodes directions in "physics" orientation, the DICOM tag Phase Encoding Direction (0018,1312) must be used to convert to scanner or physical orientation. I had a look through the DICOM conversion code, but couldn't figure out where the reading was happening.

@jdtournier : how does the code currently deal with this?

@jdtournier
Copy link
Member

This is where this happens for the DW scheme. I guess things could be handled similarly for the PE direction, though that might require some refactoring...

@Lestropie
Copy link
Member

Somewhat relates to #1038.

The issue is not with the reference axes relative to which "directions" are defined. The DICOM standard specifies for a given cartesian slice whether the phase encoding was applied along the row or column. The problem is that it does not specify the sign of phase encoding along such. In similar cases I've seen how dcm2niix will export field PhaseEncodingAxis instead of PhaseEncodingDirection, specifically because the sign information required for the latter cannot be determined from the header data. Unless @neurolabusc has found a trick I'm not aware of (or have since forgotten), what's required is pressure on the vendor to actually encode this information in their DICOMs.

@jdtournier
Copy link
Member

The issue is not with the reference axes relative to which "directions" are defined.

Yes, you're right. I got sidetracked by the original report which was about the DEC maps being wrong when the PE direction is not AP - which I think will relate to the fact that GE tend to rely on the PRS frame rather than the DICOM frame for directional information (at least it does in other contexts). But this shouldn't impact on extracting the correct PE direction, as you say.

@neurolabusc
Copy link

neurolabusc commented Oct 27, 2022

@Lestropie

For modern GE scanners you can determine phase encoding polarity from the public RectilinearPhaseEncodeReordering (0018,9034). So you can detect:

(0018,9034) CS [LINEAR] 

versus

(0018,9034) CS [REVERSE_LINEAR]

For older datasets, you can decode the private UserDefineDataGE (0043,102A). The dcm2niix code for this is here. Albeit, this code was written long ago by reverse engineering this tag which was not publicly documented at the time (I know some of the people who signed the GE NDA leaked documents onto Github, but I did not look at those when I wrote this code). However, my solution has proven robust. @mr-jaemin might be able to provide a clearer definition of this tag.

@mr-jaemin has provided sample DICOM datasets that illustrate both methods.

@Lestropie
Copy link
Member

Sorry I myself overlooked the fact that there was indeed an issue with interpreted diffusion directions; I commented based on the deviation of the discussion from the topic title. So we need to explicitly consider two separate issues:

  1. Phase encoding direction not being extracted. This is something that could potentially be implemented thanks to @neurolabusc's contribution. It is however unrelated to the way in which GE encodes its diffusion sensitisation directions.

  2. Diffusion sensitisation directions not being encoded / decoded correctly, leading to poor tractography / DEC. MRtrix3 is supposed to handle these correctly based on the PRS vs. XYZ distinction, as per the code linked by @jdtournier. Again here I would suggest running dcm2niix and see whether MRtrix3 commands using the bvec produced by that software looks correct: if yes, then MRtrix3's DICOM handling has some residual issue; if no, there might be a DICOM creation issue.

@blezek
Copy link
Contributor Author

blezek commented Nov 2, 2022

2. Again here I would suggest running dcm2niix and see whether MRtrix3 commands using the bvec produced by that software looks correct: if yes, then MRtrix3's DICOM handling has some residual issue; if no, there might be a DICOM creation issue.

@Lestropie, I can confirm that converting from DICOM using dcm2niix, then using MRTrix3 performs as expected. For the record, our standard protocol is frequency encoding in the l/r direction, but according to our lead GE tech, the scanner sometimes flips the frequency encoding (for reasons known only to GE engineers), and there is no way to lock it in. We discovered this issue during a protocol development scan when our best tech noticed that the DTI was acquired in the a/p frequency encoding direction and repeated the scan in the l/r frequency encoding direction.

processing steps

Processed both l/r and a/p frequency encoding directions using "fsl" (really dcm2niix) and directly from DICOM ("dicom").

~/Source/dcm2niix/build/bin/dcm2niix -o . -f "fsl-$suffix" "$dicom"
mrconvert -force -fslgrad fsl-$suffix.bvec fsl-$suffix.bval fsl-$suffix.nii dti-fsl-$suffix.mif  
mrconvert -force "$dicom" dti-dicom-$suffix.mif

dwi2mask -quiet -force dti-fsl-$suffix.mif mask-$suffix.mif	
dwi2mask -quiet -force dti-dicom-$suffix.mif mask-dicom-$suffix.mif

tckgen \
  -force \
  -nthreads 16 \
  -algorithm Tensor_Prob \
  -select 20000 \
  -minlength 60 \
  -seed_image mask-$suffix.mif \
  -mask mask-$suffix.mif \
  dti-fsl-$suffix.mif \
  streamlines-fsl-$suffix.tck

tckgen \
  -force \
  -nthreads 16 \
  -algorithm Tensor_Prob \
  -select 20000 \
  -minlength 60 \
  -seed_image mask-dicom-$suffix.mif \
  -mask mask-dicom-$suffix.mif \
  dti-dicom-$suffix.mif \
  streamlines-dicom-$suffix.tck

dwiextract -quiet -force dti-fsl-$suffix.mif - -bzero | mrmath -quiet -axis 3 -force - mean b0-fsl-$suffix.nii.gz
dwiextract -quiet -force dti-dicom-$suffix.mif - -bzero | mrmath -quiet -axis 3 -force - mean b0-dicom-$suffix.nii.gz

visualization

Left column is frequency encoding in the l/r direction, right column is frequency encoding in the a/p direction. Top row is loading using dcm2niix, then mrconvert, bottom row is mrconvert directly from the DICOM. Streamlines are overlaid on the b0 image. The mrconvert directly from DICOM does not produce correct streamlines (lower left image).

2022-11-02_07-41-23

Thanks for all the attention on this issue.

@jdtournier
Copy link
Member

@blezek, I might have some idea about what might be going wrong on the import of the DW directions. I can try to implement some fixes and you can check whether that fixes the issue, or if you can share one example dataset with me, I might be able to make faster progress...

@blezek
Copy link
Contributor Author

blezek commented Nov 2, 2022

or if you can share one example dataset with me

I'm working through the proper channels to get my images out of my institution, unfortunately, it may take a some time (2 weeks+). Will happily test any fixes in the mean time!

@jdtournier
Copy link
Member

OK, I've had a go at fixing what I think might be the issue. Can you try the fix_GE_grad_import branch I just pushed (also in PR #2523) and see whether it fixes that side of things?

@mr-jaemin
Copy link

I quickly tested two datasets: "R/L" frequency encoding and "A/P" frequency encoding using dcm2niix and mrconvert. The issue would be related to GE convention. As noted dcm2niix/GE or ge-mri for more details, the GE convention for reported diffusion gradient direction has always been in “MR physics” logical coordinate, i.e Freq (X), Phase (Y), Slice (Z). Note that this is neither “with reference to the scanner bore” (like Siemens or Philips) nor “with reference to the imaging plane” (as expected by FSL tools). This is the main source of confusion.

"R/L" frequency encoding (i.e. "A/P" phase encoding)
dcm2niix output:

0 0.654875 -0.271924 -0.958586 0.11881 0.326879 -0.955003
0 -0.355659 0.933965 -0.187305 0.110437 0.866547 -0.186604
0 -0.666817 0.231877 0.22031 0.986756 0.377155 0.219487

mrconvert output:

0 0.654875162745229 -0.271923961120244 -0.957386756364346 0.118810027790566 0.326879092309074 -0.957386903620765
-0 -0.355659083615388 0.933964896447314 -0.187071218460557 0.110437025997087 0.866547224504841 -0.187070016939661
0 0.66681716945632 -0.231876975035755 -0.220033992607682 -0.986756221408296 -0.377155096304502 -0.220034373385653

"A/P" frequency encoding (i.e. "R/L" phase encoding)
dcm2niix output:

0 0.355659 -0.933965 0.187305 -0.110437 -0.866547 0.186604
0 0.654875 -0.271924 -0.958586 0.11881 0.326879 -0.955003
0 -0.666817 0.231877 0.22031 0.986756 0.377155 0.219487

mrconvert output:

0 0.654875161727358 -0.271923945456056 -0.957386753335705 0.118810029903769 0.326879105692582 -0.95738690059212
-0 -0.355659074039292 0.933964900449171 -0.187071235934739 0.11043702827826 0.866547221357538 -0.187070034412982
0 0.66681717553057 -0.231876977344428 -0.220033991048607 -0.98675622099648 -0.377155092087822 -0.220034371823866

@jdtournier
Copy link
Member

Thanks for looking into it, @mr-jaemin. MRtrix has been able to deal with GE data for years though, so there's more to it than that. My gut feeling is it'll be due to there being a mix of GE private tags and DICOM standard tags both carrying information about the gradient table, and our code might be getting confused about whether it needs to re-orient the gradients or not. The way it has worked so far is to assume that the gradients need to be reoriented if we encounter one of the GE-specific private tags, but if we encounter the DICOM standard ones later in the parse and record those gradients, we'll be rotating them when we shouldn't. The pull request I just pushed (#2523) should fix that - assuming that is indeed the problem... 🤞

@neurolabusc
Copy link

@blezek with regards to your comment For the record, our standard protocol is frequency encoding in the l/r direction, but according to our lead GE tech, the scanner sometimes flips the frequency encoding (for reasons known only to GE engineers), and there is no way to lock it in. a friend who is an expert with GE told me:

When planning scans on a localiser, often the interface places them with a different orientation – e.g., you originally saved the protocol with slightly tilted axial slices, but when loading it back in again and clicking on a new localiser image, the default slices appear as a coronal stack. The scan operator can then click and rotate the stack back to axial – but depending on whether they rotate the stack clockwise or anticlockwise, will affect whether slice 1 is at the top or bottom of the now axial stack, i.e., ascending or descending. I know that you had a correction for this in dcm2niix, if I remember whether the z coordinate of the first slice was greater or less than the last slice, or something similar.

I suppose that as GE stores its diffusion directions relative to freq/phase/slice, then ascending or descending slice ordering could effectively be changing the sign of the z diffusion direction? And obviously a change in sign of the z diffusion direction could manifest as flips in x or y as well, once the tensor is calculated?

Also – although I’ve never seen this – if when clicking on a localiser image to plan slices, results in a coronal or sagittal block which is then manually rotated back to axial – then the scanner could “helpfully” swap the phase encoding direction to avoid wrap, and this phase encoding setting doesn’t revert back when rotating the stack back to axial? But this should be picked up from the DICOM tag 0018,1312 as ROW or COL as you’ve already mentioned.

Perhaps @mr-jaemin can confirm whether these user actions change the protocol property. If this behavior can be replicated, perhaps GE can update the software to allow these slice parameters to be locked.

@jdtournier
Copy link
Member

@blezek, any chance you could check out branch fix_GE_grad_import (#2523) and see if that fixes the RGB issue?

It won't fix the PE issue yet, that'll need more investigation...

@Lestropie
Copy link
Member

the scanner sometimes flips the frequency encoding

To be clear, what you're describing elsewhere is not so much a "flip" (as in 180-degree rotation) but a 90-degree rotation (is. exchanging frequency encoding and phase encoding axes)?

An alternative explanation for happening here is that the system is determining that a L>>R or R>>L phase encoding direction would result in peripheral nerve stimulation; but rather than refusing to progress with the scan it instead rotates the phase encoding direction to A>>P or P>>A. In my own protocol development I really wanted to use 4 phase encoding directions, but to use L>>R or R>>L it was necessary to drastically reduce the bandwidth, increasing EPI distortions and extending TE. It may only require a slight difference in FoV angulation to trigger this.

People may want to use L>>R and R>>L because that's what HCP did to reduce TE. But they were able to get away with it only because the Connectom gradient system doesn't have the same spatial extent as a typical whole-body system, so the gradients roll off below the shoulders and PNS is less of a problem. This is taken to the extreme with the latest head-only gradient coil systems, where they can push the gradient slew rates higher for the same reason.

Only a hypothesis though.

@jdtournier
Copy link
Member

My gut feeling is it'll be due to there being a mix of GE private tags and DICOM standard tags both carrying information about the gradient table, and our code might be getting confused about whether it needs to re-orient the gradients or not.

Turns out my gut was wrong, and we've simply never properly handled GE data when the phase-encode direction wasn't AP... Preliminary fix on #2523, though it'll require more extensive testing...

@jdtournier
Copy link
Member

jdtournier commented Nov 16, 2022

For modern GE scanners you can determine phase encoding polarity from the public RectilinearPhaseEncodeReordering (0018,9034).

Thanks @neurolabusc! I've just pushed a commit (31f819d) to parse this field, and that does seem to do the trick - though currently I only have data with a single PE direction, so I won't be able to fully validate till I get hold of some reversed PE data... 😁

[EDIT: I've just spotted the links posted by @neurolabusc with sample validation data - I'll have a look at those now]

@jdtournier
Copy link
Member

jdtournier commented Nov 17, 2022

@neurolabusc, @mr-jaemin: I've been looking into the sample dcm_qa_ge test data, which is very useful. I have however hit this weird issue, I wouldn't mind your thoughts:

The 04_Ax_DWI_TENSOR_R2MB2 contains both RectilinearPhaseEncodeReordering and UserDefineDataGE tags, but when I parse the information, the PE directions don't match. I guess I can prefer the standard RectilinearPhaseEncodeReordering tag in this case, but it does suggest there may be something wrong in the interpretation of the UserDefineDataGE tag in general?

I note all the other Ax_DTI_TENSOR_XX datasets also have both tags defined, but they all match. The odd one out is 04_Ax_DWI_TENSOR_R2MB2...

Could this be due to this bit of code, where the polarity is inverted if the slice order is bottom-up? Why would this have anything to do with the PE direction...?!?

Also, my DW gradients come out reversed compared to dcm2niix. I'm not sure that's necessarily an issue though. I can't think of any application where that would be a problem, but would this be worth getting to grips with...?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

6 participants