Possible incorrect SliceTiming

[po09i]

Describe the bug

The acquisition was acquired on Siemens with 5 purely sagittal slices (slices are RL/LR).

The NIfTI scanner coordinates increase with slice indexes (index 0: -6.270688, index 4: 13.72931)
The SliceTiming reported is: "SliceTiming": [0, 0.51562, 1.03125, 1.53125, 2.04688]

The DICOMs for a single echo/volume reports:

AcquisitionTime	SliceLocation
16:01:01.210000	-13.729311943054
16:01:01.717500	-8.7293119430542
16:01:02.227500	-3.7293121814728
16:01:02.737500	1.2706878185272
16:01:03.245000	6.2706880569458

I believe the DICOMs are in LPS+ coordinates whereas NIfTIs are RAS+ so I would expect the Slice location of DICOMs and NIFTIs to be going with opposite polarity with respect to acquisition time. If DICOM slice coordinates increase with higher acquisition times, then slices in NIfTI coordinates should decrease with higher acquisition times. Another way of looking at it would be that NIfTI coordinate -6.270688 (dcm2niix reports a SliceTiming of 0s) refers to DICOM coordinate 6.270688 which is the last slice acquired in the DICOMs.

Am I missing something or is this not intended?

To reproduce
Link to dataset

Steps to reproduce the behavior:

1. Run the command './bin/dcm2niix -o niftis/ -z y 02-gre_field_mapping_PMUlog'
2. Compare SliceTiming from JSON file and from DICOMs

Expected behavior

I would expect a SliceTiming similar to: [2.04688, 1.53125, 1.03125, 0.51562, 0]

Output log

./bin/dcm2niix -o niftis/dcm2niix_test -z y niftis/sourcedata/02-gre_field_mapping_PMUlog  
Chris Rorden's dcm2niiX version v1.0.20240202  (JP2:OpenJPEG) (JP-LS:CharLS) Clang15.0.0 ARM (64-bit MacOS)
Found 300 DICOM file(s)
Slices not stacked: echo varies (TE 2.46, 4.92; echo 1, 2). Use 'merge 2D slices' option to force stacking
Warning: Discrepancy between reported (0.0067s) and estimated (2.54603s) repetition time (issue 560).
Convert 150 DICOM as niftis/dcm2niix_test/02-gre_field_mapping_PMUlog_gre_field_mapping_PMUlog_20231128152351_2_e1e (42x64x5x30)
Warning: Discrepancy between reported (0.0067s) and estimated (2.54603s) repetition time (issue 560).
Convert 150 DICOM as niftis/dcm2niix_test/02-gre_field_mapping_PMUlog_gre_field_mapping_PMUlog_20231128152351_2_e2e (42x64x5x30)
Conversion required 0.236481 seconds (0.138406 for core code).

Note: I removed full paths to show as relative paths in the log

Version

v1.0.20240202

[neurolabusc]

@po09i dcm2niix does not rotate EPI data out of plane. An image that is acquired as sagittal DICOM slices will be saved as NIfTI sagittal images: the columns [i] will be A<->P, the rows [j] will be I<->S and the slices [k] will be L<->R. You can use fslhd to show the s-form of the NIfTI header to confirm this. The NIfTI s-form transforms voxels as stored to disk to the NIfTI coordinate system (which uses RAS as its identity matrix).

As an aside, your volumes are gre_field_mapping which suggests field maps where slice timing is not important.

As a concrete example, here is a sagittal NIfTI image created by dcm2niix from dcm_qa

git clone git@github.com:neurolabusc/dcm_qa.git
cd ./dcm_qa/Ref
fslhd sag_asc_35sl_22.nii
...
sto_xyz:1	0.000000 0.000000 -3.600000 61.200001 
sto_xyz:2	-3.250000 0.000000 0.000000 140.319641 
sto_xyz:3	0.000000 3.250000 0.000000 -126.173706 
sto_xyz:4	0.000000 0.000000 0.000000 1.000000 
sform_xorient	Anterior-to-Posterior
sform_yorient	Inferior-to-Superior
sform_zorient	Right-to-Left

[po09i]

@neurolabusc Thanks for your answer!

Indeed the NIfTI files are still sagittal, as expected (i.e.: the slices are in the last dimension but with coordinates RAS+ rather than LPS+ when converted using the sform).

Here is the sform from the provided dataset:

sto_xyz:1	0.000000 0.000000 5.000000 -6.270688 
sto_xyz:2	-4.375000 -0.000000 0.000000 98.774040 
sto_xyz:3	0.000000 4.375000 0.000000 -78.311218 
sto_xyz:4	0.000000 0.000000 0.000000 1.000000 
sform_xorient	Anterior-to-Posterior
sform_yorient	Inferior-to-Superior
sform_zorient	Left-to-Right

As an aside, your volumes are gre_field_mapping which suggests field maps where slice timing is not important.

This is precisely what I am interested in. I don't think the sliceTiming for these slices are correct, which is what I am describing in this issue. For context, we are trying to extract SliceTiming information from field maps, for real-time shimming purposes. The timing is therefore really important to us.

[po09i]

@neurolabusc I also tested a transverse dataset and the slice timings from the JSON file were correct. SliceLocation coordinates in the DICOMs increased with acquisition times and slice coordinates also increased with acquisition time in the NIfTI/JSON. The coordinates of RAS+ and LPS+ in the slice axis for a transverse acquisition (I<-->S) are the same so this is expected.

[neurolabusc]

Great. So it sounds like dcm2niix is working correctly and as defined in the BIDS specification.

[po09i]

Great. So it sounds like dcm2niix is working correctly and as defined in the BIDS specification.

It seems right with transverse acquisitions, but sagittal acquisitions are not as explained above.

[neurolabusc]

Can I suggest you acquire slice timing validation datasets for your sequences. Here are the ones I used when developing dcm2niix, albeit they use mosaics rather than one-slice per volume. Feel free to share datasets with my institutional email.

[jcohenadad]

Hi Chris, I'm not sure I see the point in acquiring a validation dataset to address this issue. To clarify: what we observe is that the DICOM timestamp of each sagittal slice is incorrectly associated with the SliceTiming vector. In Alex's original comment #797 (comment), you can note that:

• the first slice acquired (at 16:01:01.210000) corresponds to index 4 (in physical coordinate: 13.72931), whereas
• the last slice acquired (at 16:01:03.245000) corresponds to index 0 (in physical coordinate: -6.270688)

Yet, the output SliceTiming gives:

• SliceTiming[0]=0 (should be 2.04688)
• SliceTiming[4]=2.04688 (should be 0)

[neurolabusc]

@jcohenadad I can not replicate this. My validation data come are available here but are saved as mosaics. Feel free to share a dataset that illustrates the issue with my institutional email or submit a PR with a patch. It is unclear how to resolve the issue without a concrete exemplar.

[jcohenadad]

OK, we will acquire a couple of scans (2D sagittal gre_field_mapping with interleave on/off, ascending/descending), and remove the phantom half-way and upload here.

[po09i]

A sagittal dataset is linked in the original post. Here it is for convenience:

To reproduce
Link to dataset

[po09i]

@neurolabusc I have acquired a dataset linked here similar to what you had suggested (sagittal, coronal and axial acquired using ascending/descending/interleaved) where we removed the phantom halfway through the acquisition.

Findings: What I am noticing from viewing in FSLeyes and looking at the reported slice timing in the JSON sidecar is that all sagittal SliceTimings are erroneously converted while all other acquisition SliceTimings (transverse and coronal) seem to be right.

Note that I changed to a product sequence (gre_field_mapping).

[po09i]

Hey @neurolabusc, are you able to replicate my findings using the validation dataset?

[neurolabusc]

I have an exceptionally heavy teaching, service, research and admin load this term. It's on my radar and I have booked a slot on our Prisma.

[po09i]

Of course no worries! I was making sure you had seen this :)

I have booked a slot on our Prisma.

I did not think you would need to acquire a new dataset since I acquired one, but I it makes sense to see if we can replicate on multiple scanners. Thanks!

[neurolabusc]

@po09i and @jcohenadad can you please test the latest commits to the development branch (v1.0.20240327):

git clone --branch development https://github.com/rordenlab/dcm2niix.git
cd dcm2niix/console
make

The crucial change is here. It adds the conditional:

if((dcmList[dcmSort[0].indx].manufacturer == kMANUFACTURER_SIEMENS) && (dcmList[dcmSort[0].indx].CSA.mosaicSlices < 2))
  dcmList[dcmSort[0].indx].CSA.protocolSliceNumber1 = -1;

The logic here is pretty complicated. The core issue is that many FSL tools will flip the order of image columns if an image has a positive determinant. The column flipping confuses many and is not even consistent within FSL (e.g. fslstats vs fslmaths). Unfortunately, this behavior is baked into the FSL/BIDS bvec format. Therefore, any BIDS compatiblte tools must either flip image columns or flip the x component of the bvec file when handling an image with a positive determinant.

To combat this, dcm2niix will always save EPI data with a negative determinant. This explains why fslstats and fslmaths always give the same results for images created by dcm2niix, and the x component of the bvec file is always in register with the image data.

However, unlike GE and UIH, for Siemens mosaics image numbering is always anatomical, regardless of the temporal order (determined by the image numbering parameter). Therefore, for Siemens we use the CSA headerProtocolSliceNumber to infer slice order. Your sample dataset suggests this behavior is specific to Siemens Mosaics, and not images stored in non-mosaic format. Your data suggests that dcm2niix reports the BIDS slice timing in reverse order when the data is acquired with sagittal acquisition (regardless of ascending or descending slice order) and saved as non-mosaic data. I do think we need a dataset from Siemens enhanced XA to see if the suggested change helps or hurts those images.

[po09i]

Awesome! Thanks for the detailed explanation.

I tested the development branch and it does indeed solve the issue. Thanks a lot!