started to split 4d philips into individual files due to "Slices not stacked: trigger time varies"

[yarikoptic]

Data is at https://github.com/datalad/example-dicom-functional and we use it for testing of https://github.com/datalad/datalad-neuroimaging/ where we run heudiconv with dcm2niix to convert into bids dataset first. Our tests started to fail, boiled down to newer dcm2niix (1.0.20200331) stopped placing them into a single 4d .nii.gz

here is how looks with elderly 1.0.20181125

$> dcm2niix -9 -b y -o . -z y .                
Compression will be faster with 'pigz' installed
Chris Rorden's dcm2niiX version v1.0.20181125  (JP2:OpenJPEG) GCC8.2.0 (64-bit Linux)
Found 630 DICOM file(s)
swizzling 3rd and 4th dimensions (XYTZ -> XYZT), assuming interslice distance is 3.300000
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 630 DICOM as ./__func_task-oneback_run-1_20140425155335_401 (80x80x35x18)
Conversion required 1.171342 seconds (1.170979 for core code).

how looks with 1.0.20200331

$> dcm2niix -9 -b y -o . -z y .
Chris Rorden's dcm2niiX version v1.0.20200331  (JP2:OpenJPEG) GCC8.3.0 (64-bit Linux)
Found 630 DICOM file(s)
Slices not stacked: trigger time varies
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t2000 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t4000 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t6000 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t8000 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t10000 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t12000 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t14000 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t16000 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t18000 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t20000 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t22000 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t24000 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t26000 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t28000 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t30000 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t32000 (80x80x35x1)
Using RWVSlope:RWVIntercept = 4.00757:0
 Philips Scaling Values RS:RI:SS = 4.00757:0:0.0132383 (see PMC3998685)
Convert 35 DICOM as ./__func_task-oneback_run-1_20140425155335_401_t34000 (80x80x35x1)
Conversion required 1.183462 seconds (1.183202 for core code).

at least according to .json, seems to be consistent in TriggerDelayTime increment at 2000 msec interval:

$> grep TriggerDelayTime *json
__func_task-oneback_run-1_20140425155335_401_t10000.json:	"TriggerDelayTime": 10000,
__func_task-oneback_run-1_20140425155335_401_t12000.json:	"TriggerDelayTime": 12000,
__func_task-oneback_run-1_20140425155335_401_t14000.json:	"TriggerDelayTime": 14000,
__func_task-oneback_run-1_20140425155335_401_t16000.json:	"TriggerDelayTime": 16000,
__func_task-oneback_run-1_20140425155335_401_t18000.json:	"TriggerDelayTime": 18000,
__func_task-oneback_run-1_20140425155335_401_t20000.json:	"TriggerDelayTime": 20000,
__func_task-oneback_run-1_20140425155335_401_t2000.json:	"TriggerDelayTime": 2000,
__func_task-oneback_run-1_20140425155335_401_t22000.json:	"TriggerDelayTime": 22000,
__func_task-oneback_run-1_20140425155335_401_t24000.json:	"TriggerDelayTime": 24000,
__func_task-oneback_run-1_20140425155335_401_t26000.json:	"TriggerDelayTime": 26000,
__func_task-oneback_run-1_20140425155335_401_t28000.json:	"TriggerDelayTime": 28000,
__func_task-oneback_run-1_20140425155335_401_t30000.json:	"TriggerDelayTime": 30000,
__func_task-oneback_run-1_20140425155335_401_t32000.json:	"TriggerDelayTime": 32000,
__func_task-oneback_run-1_20140425155335_401_t34000.json:	"TriggerDelayTime": 34000,
__func_task-oneback_run-1_20140425155335_401_t4000.json:	"TriggerDelayTime": 4000,
__func_task-oneback_run-1_20140425155335_401_t6000.json:	"TriggerDelayTime": 6000,
__func_task-oneback_run-1_20140425155335_401_t8000.json:	"TriggerDelayTime": 8000,

[neurolabusc]

@yarikoptic I believe this is an unintended consequence of the changes made for issue 189. Specifically, if you look at the Multi-echo Philips 3T Achieva dStream 5.3.0.3 images from Baxter Rogers, you will see that the images must be segmented based on the DICOM Trigger Time (0018,1060) tags (e.g. series 701 and 1201 of sample dataset).

I really do not know how to proceed, as dcm2niix seems to be following the recipe provided by your DICOMs. I note your images have been handled by a few tools, so it is unclear if the trigger time tag was initially generated by Philips software or as a time-stamp that was inserted later by some other tool. While my aim is to resolve your issue automatically, I note that this data was touched by gdcmanon, and it would be easy enough to use that tool to simply remove the mysterious tag (e.g. --remove 18,1060).

Maybe @baxpr or @drmclem can provide input. Do we ignore 0018,1060 if the maximum value is far too long to be the time between two successive heart beats? In your case, the 155 volumes have the trigger time increment with TR, so should a tool ignore the trigger time if it precisely matches the TR? I really do not know what the intention of the trigger time is in your data, as the R-wave should be occurring many times in the 310 seconds of the series. In brief, your data seems to have a creative interpretation for the meaning of 0018,1060.

[baxpr]

Issue 189 is the only time I have encountered TriggerDelayTime. In that case it was extremely helpful for dcm2niix to sort and name nii files based on the trigger delays, and I believe the implemented solution was necessary to avoid ambiguity or errors in the nii files.

I haven't seen this use of TriggerDelayTime in any of our Philips fMRIs, for what that's worth. Anybody know if these values came from the scanner, or were inserted later? I agree this use doesn't seem to be in line with the DICOM spec.

[yarikoptic]

May be @mih could also shine some light on these data he shared if there was anything special about it (dicoms were anonimized but I guess that should have not changed any other field). I also won't be surprised if it ends up being scanner software version specific behavior which was changed later to accommodate new paradigms.

[yarikoptic]

BTW in my case trigger times are all at a regular interval. May be dcm2niix should split into multiple when difference varies, like it does for slices iirc?

[baxpr]

In the specific case in our DICOMs that prompted the change, the trigger delay times were evenly spaced at 1010, 3310, 5610 ms. Thinking about it, that's even in a plausible range for fMRI volume times too.

[baxpr]

Our pulse sequence is doing this by the way:
https://pubmed.ncbi.nlm.nih.gov/19819338/

[mih]

FTR: The gdcmanon call used on the files was

    gdcmanon --dumb --remove 400,500 --remove 12,62 --remove 12,63 \
      --replace 0008,103E,func_task-oneback_run-1 \
      --replace 0018,1030,func_task-oneback_run-1 \
      --replace 0008,1030,Hanke_Stadler^0083_transrep2 \
      --replace 0010,0010,Jane_Doe \
      --replace 0010,0020,02 \
      --replace 0010,0040,F \
      --replace 0010,1010,42 \
      --replace 0010,1030,75 \
      --replace 0010,0030,19660101 \
      -i ... -o ...

They should be unmodified otherwise. I could investigate further, if needed.

[drmclem]

Our pulse sequence is doing this by the way:
https://pubmed.ncbi.nlm.nih.gov/19819338/

Hi - are you doing dynamic MP2Rage ? MP2RAge is not defualt philips seqence and will be a custom implementation. One way of perfoming the necessary two different inversion times is to re-use cardiac inversion functionality which could result in the trigger time value being populated with the inversion delays or somethng like it. Or course if you are not doing this sequence - ignore this!

[baxpr]

Not dynamic. We are as you say re-using the cardiac inversion functionality. It's a local patch, but I think the only change from product is the shape of the pulses.

[drmclem]

So - on the list above the trigger time is the inversion time ?

[baxpr]

That's right.

[neurolabusc]

So it seems like neither of these examples are actually using the public tag TriggerTime as specified by the DICOM standard. I have committed a patch that should resolve this. It leverages the fact that for @yarikoptic's data the TriggerTime (in milliseconds) is identical to the Philips private tag MRImageDynamicScanBeginTime (2005,10a0) which is specified in seconds. If these two values match, dcm2niix ignores the TriggerTime. It is a kludge, but I think it is a robust one.

@yarikoptic data where TriggerTime actually stores Dynamic Scan Begin Time:

(0018,1020) LO [5.1.1\5.1.1.0]                          #  14, 2 SoftwareVersions
(0018,1060) DS [308000]                                 #   6, 1 TriggerTime
(0018,1081) IS [0]                                      #   2, 1 LowRRValue
(0018,1082) IS [0]                                      #   2, 1 HighRRValue
(0018,1083) IS [0]                                      #   2, 1 IntervalsAcquired
(0018,1084) IS [0]                                      #   2, 1 IntervalsRejected
(0018,1088) IS [0]                                      #   2, 1 HeartRate
...
(2005,10a0) FL 308                                      #   4, 1 Unknown Tag & Data

@baxpr example where TriggerTime actually stores Inversion Time:

(0018,1020) LO [5.3.0\5.3.0.3]                          #  14, 2 SoftwareVersions
(0018,1060) DS [511]                                    #   4, 1 TriggerTime
(0018,1081) IS [0]                                      #   2, 1 LowRRValue
(0018,1082) IS [0]                                      #   2, 1 HighRRValue
(0018,1083) IS [0]                                      #   2, 1 IntervalsAcquired
(0018,1084) IS [0]                                      #   2, 1 IntervalsRejected
(0018,1088) IS [60]                                     #   2, 1 HeartRate
...
(2005,10a0) FL 0                                        #   4, 1 Unknown Tag & Data

[drmclem]

So just to be clear - these are basically single customer solutions because of implementation details of a custom MP2RAGE sequence and don't form part of the Philips product.

[baxpr]

If there is a standard triggered cardiac scan from the product sequence, don't we still need to split on the TriggerTime field? I don't have test data for this though.

[neurolabusc]

I think this should be robust. It should segment based on TriggerTime regardless of whether this truly reflect TriggerTime or some other variable. The only situation where data will not be segmented based on TriggerTime is if it matches the MRImageDynamicScanBeginTime. So while it is based on a single customer's experiences, it should work across sequences. @yarikoptic if you are happy with this solution, can you close the issue.

[yarikoptic]

as long as it works for me ;) that field (was not) extracted by dcm2niix so can't check quickly if it is the same... will dig with dcmdump later on

[yarikoptic]

FWIW, in 'my' data -- I see no other field matching the Time, e.g.

$> dcmdump MR.1.3.46.670589.11.38317.5.0.4476.2014042516045946836 | grep 34000
(0018,1060) DS [34000]                                  #   6, 1 TriggerTime

[neurolabusc]

Remember, 0018,1060 is in milliseconds, but 2005,10a0 is in seconds. In your example you would want to look for "[34]"

[yarikoptic]

;-) thank you @neurolabusc - sorry that I had missed your original comment which had all the details. Is a new release coming? or would you bless me to include the patch into neurodebian pkg?

[neurolabusc]

Stable releases are on typically about a 6 month cadence. Unless this issue impacts many users, my sense is to not generate a new release, as it seems like an edge case. I try to test releases on a wide set of images.

On the other hand, since there are no major changes in the development branch. Therefore, I think the chances for regression are minor. Feel free to update neurodebian.

[yarikoptic]

I try to test releases on a wide set of images.

Isn't that testing automated already or there are some additional sweeps of cases?

Feel free to update neurodebian.

Ok, I will go directly to current state of development branch v1.0.20190410-94-gad0e1c8 (according to git describe --tags) but in effect it is v1.0.20200331-3-gad0e1c8 or so if only release tag was merged into development which I think would be a good strategy in general. But it is actually the v1.0.20200427 as version constant was updated if I see it right. So would you bless me to claim it to be the 1.0.20200427 version or should I make it 1.0.20200331+git3-gad0e1c8 to signify that there no "official" release?

[neurolabusc]

I would suggest calling it 1.0.20200427. While I test on a battery of images, I also send release candidates to power users who have large datasets of unusual data that can not be publicly shared. You are correct that each commit is tested on a small battery that helps with regression testing. However, each vendor is constantly evolving their interpretation of the DICOM format.

[yarikoptic]

I just uploaded 1.0.20200427 to neurodebian, and consider this issue resolved. Thank you @neurolabusc !