Siemens Diffusion Phase Encoding Directionality

[matthew-brett]

By @ baratzz on October 27, 2021, 3:11pm

Hi all,

I'm trying to determine the phase encoding directionality for some diffusion-weighted scans acquired by a Siemens Prisma scanner. dcm2niix (@ neurolabusc) seems to determine phase encoding smoothly, however, looking at the private tags, I couldn't find the private tag (0019, 100e) or any other tag that adequately differentiates the two.
Looking deeper into the Series' CSA header revealed some differences under the "Diffusion" tag, however, I could not determine with confidence which of these differences is a reliable source for deducing the directionality. E.g. for an AP-encoded series:

{'Mode': 1,
 'Scheme': 1,
 'QSpaceCoverage': 1,
 'QSpaceSampling': 1,
 'QSpaceSteps': 1,
 'BValue': {'__attribute__': {'size': 128}},
 'Averages': {'__attribute__': {'size': 128}},
 'FreeDiffusionData': {'CoordinateSystem': 1,
  'Normalization': 1,
  'DiffDirVector': {'__attribute__': {'size': 0}}}}

and for a PA-encoded series:

{'DiffWeightings': 2,
 'NoiseLevel': 40,
 'DiffDirections': 6,
 'Mode': 128,
 'Scheme': 2,
 'QSpaceCoverage': 1,
 'QSpaceSampling': 1,
 'QSpaceMaxBValue': 5,
 'QSpaceSteps': 1,
 'BValue': [4000],
 'Averages': [1],
 'Averages[0]': 1,
 'FreeDiffusionData': {'Comment': {'0': " = 0x62    # 'b'",
   '1': " = 0x30    # '0'",
   '2': " = 0x73    # 's'",
   '3': " = 0x5f    # '_'",
   '4': " = 0x64    # 'd'",
   '5': " = 0x69    # 'i'",
   '6': " = 0x72    # 'r'",
   '7': " = 0x65    # 'e'",
   '8': " = 0x63    # 'c'",
   '9': " = 0x74    # 't'",
   '10': " = 0x69    # 'i'",
   '11': " = 0x6f    # 'o'",
   '12': " = 0x6e    # 'n'",
   '13': " = 0x5f    # '_'",
   '14': " = 0x69    # 'i'",
   '15': " = 0x64    # 'd'",
   '16': " = 0x6f    # 'o'",
   '17': " = 0x2e    # '.'",
   '18': " = 0x64    # 'd'",
   '19': " = 0x76    # 'v'",
   '20': " = 0x73    # 's'",
   '21': ' = 0xa',
   '22': " = 0x53    # 'S'",
   '23': " = 0x49    # 'I'",
   '24': " = 0x45    # 'E'",
   '25': " = 0x4d    # 'M'",
   '26': " = 0x45    # 'E'",
   '27': " = 0x4e    # 'N'",
   '28': " = 0x53    # 'S'",
   '29': ' = 0xa',
   '30': " = 0x49    # 'I'",
   '31': " = 0x6e    # 'n'",
   '32': " = 0x74    # 't'",
   '33': " = 0x65    # 'e'",
   '34': " = 0x72    # 'r'",
   '35': " = 0x6e    # 'n'",
   '36': " = 0x61    # 'a'",
   '37': " = 0x6c    # 'l'",
   '38': " = 0x20    # ' '",
   '39': " = 0x36    # '6'",
   '40': " = 0x2d    # '-'",
   '41': " = 0x64    # 'd'",
   '42': " = 0x69    # 'i'",
   '43': " = 0x72    # 'r'",
   '44': " = 0x65    # 'e'",
   '45': " = 0x63    # 'c'",
   '46': " = 0x74    # 't'",
   '47': " = 0x69    # 'i'",
   '48': " = 0x6f    # 'o'",
   '49': " = 0x6e    # 'n'",
   '50': " = 0x73    # 's'",
   '51': " = 0x20    # ' '",
   '52': " = 0x73    # 's'",
   '53': " = 0x65    # 'e'",
   '54': " = 0x74    # 't'"},
  'DiffDirections': 6,
  'CoordinateSystem': 1,
  'Normalization': 1,
  'DiffDirVector': [{'Sag': 1, 'Tra': 1},
   {'Sag': -1, 'Tra': 1},
   {'Cor': 1, 'Tra': 1},
   {'Cor': 1, 'Tra': -1},
   {'Sag': 1, 'Cor': 1},
   {'Sag': -1, 'Cor': 1}]}}

Any help figuring this out will be greatly appreciated!
Best,
Zvi

[matthew-brett]

By @ neurolabusc on October 27, 2021, 4:23pm

Since you refer to the CSA header, I assume this Prisma is running V* (e.g. VE11C) not X* (e.g. XA30). You can get DICOMs for validating phase encoding polarity here. Specifically, the last characters of the protocol name reveal phase encoding direction (e.g. AP, PA, RL, LR).

dcm2niix uses the CSA tag PhaseEncodingDirectionPositive from nii_dicom.cpp.

You can use the logorrheic vebosity (-v 2) to see how dcm2niix sees the DICOM header. This is very extensive, so you may want to apply it to a single DICOM file rather than an entire sessions. The minimal call would be dcm2niix -v 2 /path/to/DICOM

For Siemens X* series (enhanced DICOM, without the CSA header), you will want to read the private DICOM tag PhaseEncodingDirectionPositiveSiemens 0021,111C. See dcm_qa_enh for sample datasets.

[matthew-brett]

By @ baratzz on November 16, 2021, 5:38pm

@ neurolabusc Thank you so much for following up!

Indeed, it turned out my implementation of the CsaHeader class was partial, and returned only the ASCII header. I resolved the issue and am now returning the full CSA header, including the required PhaseEncodingDirectionPositive value.

I've started work on a method to infer phase encoding as it is specified in the BIDS standard, but am unfortunately still struggling with the translation of the DICOM metadata to reliable i/j/k/-/+ and AP/PA/RL/LR(/IS/SI?) values.

For now, I'm running the experimental version of the code locally because my lab requires it for automated BIDS-compatible path generation (following the conversion to NIfTI using dcm2niix, but I would be very grateful for any help figuring this out and coming up with a more robust solution.

[matthew-brett]

By @ neurolabusc on November 16, 2021, 6:38pm

You can look at dcm_qa for concrete examples of AP/PA/RL/LR examples. Here is the algorithm

1. Use 0018,1312 to determine COL or ROW. Use this to infer BIDS axis (COL → j ROW → i)
2. Use phase encoding direction positive to determine if the direction is inverted (-) or not. Be aware, that this depends on how the image data is stored to disk. The first line of DICOM data is at the top of the screen, the way we read English. With NIfTI, the first line is at the bottom of the screen, the way we draw a cartesian graph with increasing Y moving up the screen. Most DICOM to NIfTI converters flip row order as saved to disk (though this lossless flip could also be stored in the S-Form/Q-Form).

DICOM data is stored as 2D slices, with row and columns always referring to either frequency or phase direction. The third dimension is used for slices. Therefore, raw DICOM will never yield a BIDS k. It could be done if the row, column and slice order was changed when converting DICOM to NIfTI, but this would disrupt the BIDS slicetiming field. Therefore, in practice DICOM and BIDS agree that the 3rd dimension saved to disk is always slices.

[matthew-brett]

By @ baratzz on November 16, 2021, 7:26pm

Thank you so much for the clear explanation!

I implemented the simple COL → j / ROW → i and 0 → - / 1 → + conversion (see here and here) based on similar advice you gave in some older post concerning this, but had trouble translating that to the AP/PA/LR/RL. It was my understanding the is determined by the plane of the image (axial/sagittal/coronal, which I determine based on this SO answer), but I still haven't managed to reliably infer it over a large dataset. I will definitely take a closer look at the examples in dcm_qa when I manage to find the time to dig a little deeper.

Edit: I realize the dictionary in the first link shows opposite i/j mapping, for some reason this is the configuration that seems to work best with most of my scans.

[matthew-brett]

By @ neurolabusc on November 17, 2021, 12:31pm

Note that the file name label is arbitrary. In other words, you could always use FWD and REV for forward and reverse phase encoding. Brain images are often acquired oblique to the scanner bore to match human anatomy, and oblique to anatomy to maximize coverage (e.g. an axial scan often includes a pitch rotation relative to the AC-PC plane to provide more coverage of the cerebellum). You can compute a cross-product of the canonical slice directions and the DICOM ImageOrientationPatient, but one worries that this will sometimes be misleading.

[matthew-brett]

By @ baratzz on November 17, 2021, 1:30pm

Brain images are often acquired oblique to the scanner bore to match
human anatomy, and oblique to anatomy to maximize coverage

I was hoping rouding the ImageOrientationPatient will return sufficiently reliable results, I guess it's still a little
risky.

you could always use FWD and REV for forward and reverse phase
encoding.

That does seem like a much better solution! I think I'll try that. If BIDS-compatible apps generally rely on the JSON sidecar's IntendedFor field, it definitely sounds like the reasonable way to go. Thank you!