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 [email protected]: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!

[neurolabusc]

@jcohenadad and @po09i I think this is resolved with the latest commit to the development branch. Thanks @jeffreyluci for XA data. You can also see dcm_qa for VB17 and dcm_qa_sag for E11 data.

Canonically, EPI data where slice times are important are saved as mosaics for systems VA through VE and as enhanced DICOM for XA. Mosaics are stored in the same spatial order regardless of order of temporal acquisition and reversing image numbering swaps the order in the mosaic but does not influence acquisitions. In contrast, for XA all slices from a single 3D volume are saved in a single file. From dcm2niix v1.0.20230411, the timing for V data stored as non-mosaics shows issues when the image has a positive determinant.

At the moment, I consider this fix experimental, as while we now have V* mosaic, V* 2D, and XA enhanced, I think we need some XA slices saved as classic 2D images to ensure we handle all the use cases.

[neurolabusc]

@po09i and @jcohenadad I have reverted this. The original solution is unintuitive but correct. I have described a procedure to validate this including your sample dataset.

The basic issues is that some but not all FSL tools mirror the first dimension if the volume has a posititve determinant. Since the FSL bvec format was adopted by BIDS, this confusion infuses our formats. In other words, to interpret an FSL bvec, one must flip the polarity of the X-component (or flip the order of columns) if the volume has a positive determinant. To minimize confusion, dcm2niix will attempt to create an image with a negative determinant by flipping the slice order, the slice timing and the spatial transforms. By ensuring that data has a negative determinant, one does not have to worry about transposing data. So in your DICOMs, the instance number, acquisition time and left-ness increase together, so the data are Right-to-left in both instance number and time. However, dcm2niix reorders the slices to be Left-to-Right and therefore the slice timing (which is relative to order stored on disk) is descending. The validation dataset makes this clear.

So the two challenges are:

• DICOM +Xd is Left, whereas for NIfTI +x is Right
• The BIDS and FSL formats are confusing (and indeed different FSL tools work inconsistently) when a volume has a positive determinant. Therefore, dcm2niix attempts to help by reversing the order slices are saved to disk to enforce negative determinants.

For full disclosure, when the BIDS format was being created I was a loud critic for using the FSL bvec format for precisely these reasons (and the fact that encoding vectors in world space rather than image space makes more sense to me). Any format is a compromise, and BIDS has proved a great unifying force in our field. However, it does have nuances that we must bear in mind.

As an aside, for your validation image I did burn in a bright white line on the first instance number (e.g. the first temporal scan). You will note that tools like FSLeyes show this bright line is on the final (not first) slice on the right side.

[po09i]

Thanks @neurolabusc for the investigation. I have tested the development branch and the slice timing is what I am expecting it to be. I also tested the validation procedure which also matches my expectations.

As an aside, for your validation image I did burn in a bright white line on the first instance number (e.g. the first temporal scan). You will note that tools like FSLeyes show this bright line is on the final (not first) slice on the right side.

This makes it much easier to see the solved issue.

Thanks again for the detailed explanation.