Philips dicom with interspersed b0s gets converted incorrectly

[vuiiscci]

It appears the nifti image looks correct (b0s are interspersed correctly). But, the bvals don't look correct (they look ordered and all b0s are in the front).

Example dicom here: https://www.dropbox.com/s/g7op2sseml3h109/1.3.46.670589.11.17240.5.0.5636.2018061108033918000-901-1-1xjxrsf.dcm?dl=0

[baxpr]

I confirmed that this is a change in behavior for release v1.0.20180614 compared to release v1.0.20171215 (at least, for default options).

[baxpr]

As far as I know, Philips re-orders these at the scanner so they are sorted first by bvalue, then by time. So what we've seen in the past is the ordering that's in the bval file - but now, the nii doesn't match that.

[neurolabusc]

Historically, dcm2niix converted Philips volumes by relying on the instance number (0020,0013). However, modern Philips software may now store the 2D slices in these 4D formats in a completely jumbled order (it is truly jumbled, not the order images come off the reconstructor: you can find derived isotropic images stored prior to raw data). Volumes, slices and derived data can be completely intermingled. The instance number often has no bearing with acquisition order. Unfortunately, Philips does not provide a varying DICOM FrameAcquisitionDateTime (0018,9074) or PAR dyn_scan_begin_time to determine the actual order these images were acquired. Therefore, these files are sorted based on the DICOM Dimension Index Values (0020,9157) or PAR diffusion b value number/gradient orientation number. This is the only way to handle modern Philips data.

So this is a feature not a defect of recent releases of dcm2niix: older versions will fail with modern jumbled Philips files. The modern Philips files store no information regarding actual acquisition order. The old version preserved acquisition because it assumed instance number (0020,0013) could provide reliable information about order. However, this is not required by DICOM and this assumption will often have catastrophic consequences when applied to Philips data. I have documented the dcm2niix strategy for handling Philips DICOM and PAR/REC. For your DICOM images, the recent version faithfully converts these in the order specified by Dimension Index Values.

If you know the actual acquisition order, you can always write a simple script to reorder the NIfTI, bval and bvec files. Another option would be to make a fork of dcm2niix that sorts based on instance number - you would change this line of nii_dicom:
qsort(dcmDim, numberOfFrames, sizeof(struct TDCMdim), compareTDCMdim);
However, you should be aware that this will fail on many modern Philips datasets.

@drmclem might note that as a feature request it would be great if Philips could make Dimension Index Values correspond to acquisition time. It would also be great if instance number (0020,0013) could be trusted to have some relationship with acquisition order. The current Philips DICOM is legal, but it provides no information regarding the actual acquisition order.

The ordering of these values probably has little implication. At one time, dcm2niix included an option to sort by b-value, as this is the format required by DKE. However, @mharms raised some concern. Here were Jesper Andersson's comments on this topic from November 16, 2017 :

eddy is in principle fine with any re-ordering as long as you 
1. Make use that the first volume in your topup input is also the first volume in your eddy input
2. Re-order everything else accordingly (bvecs, bvals, index etc).
The continuous movement model in the slice-to-vol is only continuous within a volume, so it won’t affect that.
But yes, Mike has a good point that it will make the initial volume-to-volume movement correction a little trickier. For that reason you might want to add something like
  --niter=8 --fwhm=10,6,4,2,0,0,0,0
to your eddy command.

It should be noted that this is also why modern Philips datasets look jumbled when you view them with your favorite DICOM viewer like Horos. In your example the first few slices are very non-sequential

(0020,9157) UL 1\1\1\66
(0020,9157) UL 1\1\2\2  
(0020,9157) UL 1\1\2\9  
(0020,9157) UL 1\1\2\17  

[neurolabusc]

I did look at the sample data and I agree that dcm2niix is not parsing the bval in the same order as the other images. I will investigate this, as I have not seen this before. But all my comments above are accurate.

[neurolabusc]

@BaxterPRogers and @vuiiscci can you please test latest commit (v1.0.20180622). This ensures DTI get reordered along with all other jumbled Philips volume information. Since this issue impacts the recent 'stable' release (v1.0.20180614), I am eager to test this thorough for a new release. I have tried to make my changes and minimal as possible to address the issue, so I do not think there are any unintended consequences (indeed, it simplifies the code). Since I have a limited number of Philips DTI samples, the real trick is making sure we still identify the derived isotropic volumes that can be embedded in the raw diffusion series.

[baxpr]

@neurolabusc @vuiiscci I tested v1.0.20180622 -

For the scan linked above, the bvals/vecs and nifti now appear to be ordered in the same way, solving the issue we were having. There is no derived volume in this scan I think? I also did a cursory test of DTI reconstruction and fiber tracking with the dcm2niix nii/bvecs, and those results looked sane.

I also tested on a scan from last year that does have a derived volume. Seems to be working fine there - it finds the derived ADC volume and produces .nii both with and without it. The ordering of the images is different from past behavior (B0 image is at the end of the series instead of the beginning) but I gather that's related to the issues you described above? And a warning to that effect was generated.

[drmclem]

Hi We have been looking at what is the most sensible way of implementing this – especially with the commercial multiband (SMS) product now out - Siemens I believe have a private tag (0019,0010) but I’m not sure how it is defined. I think there is some interesting comments here https://wiki.humanconnectome.org/display/PublicData/HCP+fMRI+slice-timing+acquisition+parameters on applicability of slice timing. Below is the relevant DICOM standard timing references so 0018,9074 would seem the most applicable,  There are others but have historically been used for the whole set of data rather than on a per image basis (in classic) 0008,0022/0032 Acquisition Date/Time I have seen various definitions such as “The date/time the acquisition of data (set of slices) that resulted in this image started” – however since the enhanced changes this now states “The date the acquisition of data that resulted in this image started” but there may be some historical clashing. 0008,0023/0033 Content Date/Time The date/time the image pixel data creation started. Required if image is part of a series in which the images are temporally related. May be present otherwise. 0008,0023 could imply this should be the reconstruction time (which would be unhelpful) Matthew

[vuiiscci]

@neurolabusc Is there any way we could put the b0 in the beginning? For eddy/topup its desired to have the scans ordered in which they were acquired... Although I do understand maybe this is impossible to obtain from the dicom tags.

[neurolabusc]

@vuiiscci it is impossible to determine the acquisition order from the DICOM tags. If you want, you can edit nii_dicom_batch.cpp so that the line that is currently opts->isSortDTIbyBVal = false; reads opts->isSortDTIbyBVal = true; and recompile. This will sort by bvalues, regardless of the order specified by 0020,9157.

You do this at your own risk. I do not want to add this to the new release, as I have not tested code in a while and as a rule dcm2niix does not second guess the ordering data provided by the vendor. The goal is to be lossless. I understand that the order provided by Philips 0020,9157 may not always match the temporal acquisition order. But the job of dcm2niix is to extract the data as described by the vendor. The code is already complicated enough and hard enough to maintain without adding further demands. I really think you are better off writing a separate script to do this re-ordering for you. For example, if you look at my Eddy/TOPUP shell scripts you will see that I identify minBvalIdx from the bval file.

[neurolabusc]

@drmclem thanks for your feedback. Would be ideal if Philips could have new public tags defined, but using private tags is probably the most expedient. I think the typical wish list for data from Philips data is

• Slice timing: for 2D EPI, how long after the first slice was this slice acquired. This value will allow the user to determine slice order (ascending, descending, interleaved), multi-band factor (which slices were acquired simultaneously) and if successive volumes are separated by a pause (sparse imaging or acquisitions with prospective motion correction).
• Phase encoding polarity. For tools like TOPUP we want to known not only the dimension of the phase encoding (e.g. A<->P vs R<->L)but the polarity (A->P vs P->A).
• Effective echo spacing for calibrated maps from FUGUE. This is a want rather than a need: the spatial undistortion is not influenced by this value, but since it provides calibrated maps it allows the user to assess different sequences.
• PartialFourier : nice to check for mrdegibbs
• ReconMatrixPE : hopefully no scientist willingly acquires interpolated data, but if the data is interpolated it should be encoded. This has many implications including mrdegibbs
• ParallelReductionFactorInPlane: e.g. SENSE factor
• MultibandAccelerationFactor: this is encoded by Slice timing
• Volume order: critical for fMRI so we know that data is in correct temporal order, but also useful for DTI.

With your regards to slice timing and multiband (SMS), I think this issue is beyond the scope of this issue. I think it is important to provide the users with these parameters, even if in typical cases it may not be useful. In particular, it is possible that future tools may leverage this information in ways we do not at the moment. For example, all current tools compute motion correction and slice timing as two separate stages, but an optimal motion correction should have knowledge about slice timing. Therefore, including these data may help future studies.

For completeness, here are a few thoughts for those interested in slice timing and SMS. I also have a web page on this topic.

• The basic paradox: we need slice timing correction when the TR is long (where timing differences between slices in a volume is large), yet we can only accurately infer slice timing when the TR is short (to minimize interpolation errors).
• SMS dramatically reduces TR, so slice timing does less.
• If your TR is very fast (~1s), the temporal derivative should be able to account for slice time errors, so diminishing returns for slice time correction with short.
• As your TR gets very short, you have the SNR has diminishing returns: there is less time for T1 relaxation (less signal) and more autocorrelation between signals. While the HCP sequences use very agressive multi-band, there are tradeoffs to this, including aliasing. In many cases, it may be worthwhile to relax the multi-band factor a bit, getting a slightly longer TR but avoiding some of the penalties. In other words, one should consider what the motivation for extremely short TRs is (e.g. are you trying to get below the Nyquist of physiological artifacts).

[neurolabusc]

This issue is addressed in the new release v1.0.20180622.