Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Image re-ordering Philips with dynamic scans #529

Closed
drmclem opened this issue Jul 12, 2021 · 22 comments
Closed

Image re-ordering Philips with dynamic scans #529

drmclem opened this issue Jul 12, 2021 · 22 comments

Comments

@drmclem
Copy link

drmclem commented Jul 12, 2021

We have a classic DICOM data set which I think has both the images is indirect order and is also a dynamic scan (something I would expect to be typical for an fmri for example), However the converted files have the slices scrambled.

Looking at the logs it would appear that re-ordering is only being applied to the conventional scans not the dynamics.

A working one ... (single scan)

Convert 144 DICOM as non-anon_dcm2niix/_WIP_T2_brain_cor_20210709112826_1101 (336x336x144x1)
Converting ./DICOM/IM_0405
Warning: Interslice distance varies in this volume (incompatible with NIfTI format).
Dimensions 336 336 144 1 nAcq 1 nConvert 144
 Distance from first slice:
dx=[0 30 60 90 120 150 5 35 65 95 125 155 10 40 70 100 130 160 15 45 75 105 135 165 20 50 80 110 140 170 25 55 85 115 145 175 1.25 31.25 61.25 91.25 121.25 151.25 6.25 36.25 66.25 96.25 126.25 156.25 11.25 41.25 71.25 101.25 131.25 161.25 16.25 46.25 76.25 106.25 136.25 166.25 21.25 51.25 81.25 111.25 141.25 171.25 26.25 56.25 86.25 116.25 146.25 176.25 2.5 32.5 62.5 92.5 122.5 152.5 7.5 37.5 67.5 97.5 127.5 157.5 12.5 42.5 72.5 102.5 132.5 162.5 17.5 47.5 77.5 107.5 137.5 167.5 22.5 52.5 82.5 112.5 142.5 172.5 27.5 57.5 87.5 117.5 147.5 177.5 3.75001 33.75 63.75 93.75 123.75 153.75 8.75 38.75 68.75 98.75 128.75 158.75 13.75 43.75 73.75 103.75 133.75 163.75 18.75 48.75 78.75 108.75 138.75 168.75 23.75 53.75 83.75 113.75 143.75 173.75 28.75 58.75 88.75 118.75 148.75 178.75]
Warning: Order specified by DICOM instance number is not spatial (reordering).
Slice re-ordering resolved inter-slice distance variability.
Philips Scaling Values RS:RI:SS = 6.30281:0:0.0528026 (see PMC3998685)
 R = raw value, P = precise value, D = displayed value
 RS = rescale slope, RI = rescale intercept,  SS = scale slope
 D = R * RS + RI    , P = D/(RS * SS)
 D scl_slope:scl_inter = 6.30281:0
 P scl_slope:scl_inter = 18.9384:0
 Using P values ('-p n ' for D values)

Not working - no test done and no re-ordering ?

Convert 144 DICOM as non-anon_dcm2niix/_WIP_T2_brain_sag_20210709112826_1201 (336x336x144x1)
Converting ./DICOM/IM_0101
Philips Scaling Values RS:RI:SS = 6.88767:0:0.0524648 (see PMC3998685)
 R = raw value, P = precise value, D = displayed value
 RS = rescale slope, RI = rescale intercept,  SS = scale slope
 D = R * RS + RI    , P = D/(RS * SS)
 D scl_slope:scl_inter = 6.88767:0
 P scl_slope:scl_inter = 19.0604:0
 Using P values ('-p n ' for D values)

Would you expect this to work ?

Does this code imply that it only is applied to 4d scans with 1 in the 4th dimension ? (latest git hub version)

            if (hdr0.dim[4] < 2) {
                if (dxVaries) {
                    sliceMMarray = (float *) malloc(sizeof(float)*nConvert);
                    sliceMMarray[0] = 0.0f;
                    for (int i = 1; i < nConvert; i++)
						sliceMMarray[i] = intersliceDistance(dcmList[dcmSort[0].indx],dcmList[dcmSort[i].indx]);
                    printWarning("Interslice distance varies in this volume (incompatible with NIfTI format).\n");

Best wishes

Matthew

@neurolabusc
Copy link
Collaborator

Closing this issue as I do not have any access to a sample dataset that exhibits this issue. Either someone needs to validate a solution and provide a pull request, or I need to be provided with a sample dataset that exhibits this issue.

As an aside, while dcm2niix should handle these situations, Philips equipment should generate images in an order that is meaningful (e.g. sequential instances are sequential in space or time). I suspect the current jumbled order reflects the order of data returned from a parallel reconstruction system. While we can fix this one tool, such jumbled data is likely to cause grief for others. There are lots of legal ways to make a legal DICOM that is inscrutable, manufacturers should seek parsimonious solutions for this complicated format. Philips should be admonished for this behavior. We need to cure the root problem, not treating the symptoms.

@baxpr
Copy link

baxpr commented Jul 29, 2021

@neurolabusc I have emailed you an example of what is possibly the same problem - an fMRI where the slices are getting mis-sorted by v1.0.20210317. There is no relationship between InstanceNumber and slice position or TemporalPositionIdentifier in the DICOM files, which I gather might be the underlying issue.

These classic frames were exported on the scanner using a python interface to the DICOM leacher.

@neurolabusc
Copy link
Collaborator

@baxpr your issue seems different, but is indicative of how Philips assigns instance numbers with no rationale. My kludge generates an instance number based on three private Philips tags.

A sensible formula would be:

private volumeNumber (2005,1063)
private sliceNumber (2001,100A)
private slicesPerVolume (2001,1018)
instanceNumber = sliceNumber + ((volumeNumber-1) * slicesPerVolume)

However, your does not follow a logical pattern, and indeed uses a different pattern for each value. You can see this by running the development branch of dcm2niix in verbose mode (-v 1):

Warning: Philips instance number (10) does not make sense: slice 19 of 57, volume 1
Warning: Philips instance number (4) does not make sense: slice 13 of 57, volume 1
Warning: Philips instance number (38) does not make sense: slice 47 of 57, volume 1
Warning: Philips instance number (39) does not make sense: slice 48 of 57, volume 1
Warning: Philips instance number (5) does not make sense: slice 14 of 57, volume 1
Warning: Philips instance number (11) does not make sense: slice 1 of 57, volume 1
Warning: Philips instance number (7) does not make sense: slice 16 of 57, volume 1
Warning: Philips instance number (13) does not make sense: slice 3 of 57, volume 1
Warning: Philips instance number (12) does not make sense: slice 2 of 57, volume 1
Warning: Philips instance number (6) does not make sense: slice 15 of 57, volume 1
Warning: Philips instance number (2) does not make sense: slice 11 of 57, volume 1
Warning: Philips instance number (16) does not make sense: slice 6 of 57, volume 1
Warning: Philips instance number (17) does not make sense: slice 7 of 57, volume 1
Warning: Philips instance number (3) does not make sense: slice 12 of 57, volume 1
Warning: Philips instance number (29) does not make sense: slice 38 of 57, volume 1
Warning: Philips instance number (15) does not make sense: slice 5 of 57, volume 1
Warning: Philips instance number (1) does not make sense: slice 10 of 57, volume 1
Warning: Philips instance number (14) does not make sense: slice 4 of 57, volume 1
Warning: Philips instance number (28) does not make sense: slice 37 of 57, volume 1
Warning: Philips instance number (49) does not make sense: slice 20 of 57, volume 1
Warning: Philips instance number (48) does not make sense: slice 57 of 57, volume 1
Warning: Philips instance number (52) does not make sense: slice 23 of 57, volume 1
Warning: Philips instance number (46) does not make sense: slice 55 of 57, volume 1
Warning: Philips instance number (47) does not make sense: slice 56 of 57, volume 1
Warning: Philips instance number (53) does not make sense: slice 24 of 57, volume 1
Warning: Philips instance number (45) does not make sense: slice 54 of 57, volume 1
Warning: Philips instance number (51) does not make sense: slice 22 of 57, volume 1
Warning: Philips instance number (50) does not make sense: slice 21 of 57, volume 1
Warning: Philips instance number (44) does not make sense: slice 53 of 57, volume 1
Warning: Philips instance number (40) does not make sense: slice 49 of 57, volume 1
Warning: Philips instance number (54) does not make sense: slice 25 of 57, volume 1
Warning: Philips instance number (55) does not make sense: slice 26 of 57, volume 1
Warning: Philips instance number (41) does not make sense: slice 50 of 57, volume 1
Warning: Philips instance number (57) does not make sense: slice 28 of 57, volume 1
Warning: Philips instance number (43) does not make sense: slice 52 of 57, volume 1
Warning: Philips instance number (42) does not make sense: slice 51 of 57, volume 1
Warning: Philips instance number (56) does not make sense: slice 27 of 57, volume 1
Warning: Philips instance number (31) does not make sense: slice 40 of 57, volume 1
Warning: Philips instance number (25) does not make sense: slice 34 of 57, volume 1
Warning: Philips instance number (19) does not make sense: slice 9 of 57, volume 1
Warning: Philips instance number (18) does not make sense: slice 8 of 57, volume 1
Warning: Philips instance number (24) does not make sense: slice 33 of 57, volume 1
Warning: Philips instance number (30) does not make sense: slice 39 of 57, volume 1
Warning: Philips instance number (26) does not make sense: slice 35 of 57, volume 1
Warning: Philips instance number (32) does not make sense: slice 41 of 57, volume 1
Warning: Philips instance number (33) does not make sense: slice 42 of 57, volume 1
Warning: Philips instance number (27) does not make sense: slice 36 of 57, volume 1
Warning: Philips instance number (23) does not make sense: slice 32 of 57, volume 1
Warning: Philips instance number (37) does not make sense: slice 46 of 57, volume 1
Warning: Philips instance number (36) does not make sense: slice 45 of 57, volume 1
Warning: Philips instance number (22) does not make sense: slice 31 of 57, volume 1
Warning: Philips instance number (8) does not make sense: slice 17 of 57, volume 1
Warning: Philips instance number (34) does not make sense: slice 43 of 57, volume 1
Warning: Philips instance number (20) does not make sense: slice 29 of 57, volume 1
Warning: Philips instance number (21) does not make sense: slice 30 of 57, volume 1
Warning: Philips instance number (35) does not make sense: slice 44 of 57, volume 1
Warning: Philips instance number (9) does not make sense: slice 18 of 57, volume 1
Warning: Philips instance number (25847) does not make sense: slice 35 of 57, volume 2
Warning: Philips instance number (25853) does not make sense: slice 41 of 57, volume 2
Warning: Philips instance number (25852) does not make sense: slice 40 of 57, volume 2
Warning: Philips instance number (25846) does not make sense: slice 34 of 57, volume 2
Warning: Philips instance number (25850) does not make sense: slice 38 of 57, volume 2
Warning: Philips instance number (25844) does not make sense: slice 32 of 57, volume 2
Warning: Philips instance number (25878) does not make sense: slice 28 of 57, volume 2
Warning: Philips instance number (25845) does not make sense: slice 33 of 57, volume 2
Warning: Philips instance number (25851) does not make sense: slice 39 of 57, volume 2
Warning: Philips instance number (25869) does not make sense: slice 57 of 57, volume 2
Warning: Philips instance number (25855) does not make sense: slice 43 of 57, volume 2
Warning: Philips instance number (25841) does not make sense: slice 29 of 57, volume 2
Warning: Philips instance number (25840) does not make sense: slice 18 of 57, volume 2
Warning: Philips instance number (25854) does not make sense: slice 42 of 57, volume 2
Warning: Philips instance number (25868) does not make sense: slice 56 of 57, volume 2
Warning: Philips instance number (25842) does not make sense: slice 30 of 57, volume 2
Warning: Philips instance number (25856) does not make sense: slice 44 of 57, volume 2
Warning: Philips instance number (25857) does not make sense: slice 45 of 57, volume 2
Warning: Philips instance number (25843) does not make sense: slice 31 of 57, volume 2
Warning: Philips instance number (25824) does not make sense: slice 5 of 57, volume 2
Warning: Philips instance number (25830) does not make sense: slice 17 of 57, volume 2
Warning: Philips instance number (25831) does not make sense: slice 19 of 57, volume 2
Warning: Philips instance number (25825) does not make sense: slice 7 of 57, volume 2
Warning: Philips instance number (25833) does not make sense: slice 4 of 57, volume 2
Warning: Philips instance number (25827) does not make sense: slice 11 of 57, volume 2
Warning: Philips instance number (25826) does not make sense: slice 9 of 57, volume 2
Warning: Philips instance number (25832) does not make sense: slice 2 of 57, volume 2
Warning: Philips instance number (25836) does not make sense: slice 10 of 57, volume 2
Warning: Philips instance number (25822) does not make sense: slice 1 of 57, volume 2
Warning: Philips instance number (25823) does not make sense: slice 3 of 57, volume 2
Warning: Philips instance number (25837) does not make sense: slice 12 of 57, volume 2
Warning: Philips instance number (25835) does not make sense: slice 8 of 57, volume 2
Warning: Philips instance number (25834) does not make sense: slice 6 of 57, volume 2
Warning: Philips instance number (25839) does not make sense: slice 16 of 57, volume 2
Warning: Philips instance number (25838) does not make sense: slice 14 of 57, volume 2
Warning: Philips instance number (25828) does not make sense: slice 13 of 57, volume 2
Warning: Philips instance number (25829) does not make sense: slice 15 of 57, volume 2
Warning: Philips instance number (25866) does not make sense: slice 54 of 57, volume 2
Warning: Philips instance number (25872) does not make sense: slice 22 of 57, volume 2
Warning: Philips instance number (25873) does not make sense: slice 23 of 57, volume 2
Warning: Philips instance number (25867) does not make sense: slice 55 of 57, volume 2
Warning: Philips instance number (25871) does not make sense: slice 21 of 57, volume 2
Warning: Philips instance number (25865) does not make sense: slice 53 of 57, volume 2
Warning: Philips instance number (25859) does not make sense: slice 47 of 57, volume 2
Warning: Philips instance number (25858) does not make sense: slice 46 of 57, volume 2
Warning: Philips instance number (25864) does not make sense: slice 52 of 57, volume 2
Warning: Philips instance number (25870) does not make sense: slice 20 of 57, volume 2
Warning: Philips instance number (25848) does not make sense: slice 36 of 57, volume 2
Warning: Philips instance number (25874) does not make sense: slice 24 of 57, volume 2
Warning: Philips instance number (25860) does not make sense: slice 48 of 57, volume 2
Warning: Philips instance number (25861) does not make sense: slice 49 of 57, volume 2
Warning: Philips instance number (25875) does not make sense: slice 25 of 57, volume 2
Warning: Philips instance number (25849) does not make sense: slice 37 of 57, volume 2
Warning: Philips instance number (25863) does not make sense: slice 51 of 57, volume 2
Warning: Philips instance number (25877) does not make sense: slice 27 of 57, volume 2
Warning: Philips instance number (25876) does not make sense: slice 26 of 57, volume 2
Warning: Philips instance number (25862) does not make sense: slice 50 of 57, volume 2

The state of Philips DICOM images has been in disarray for several years now. They were the first to market with enhanced DICOM, but this was implemented in a very inefficient way both in terms of disk space and parsing speed. Since then, Canon and Siemens have released much cleaner enhanced DICOMs. While the initial XA10 Siemens had a lot of rough edges, Siemens quickly cleaned things up. In contrast, Philips seems to have not addressed the issues with either their enhanced or classic scans. I notice that your protocol is listed as a WIP. I am not sure I want to retain the kludge I wrote today in the final form. My sense is use should use the manufacturer provided NIfTI export functions until Philips is able to resolve their DICOMs.

To be clear, this is a limitation of Philips images, not my software. Their implementation impacts other tools like Horos and dicm2nii.

@drmclem
Copy link
Author

drmclem commented Jul 30, 2021

Hi

The instances of data that have been produced which show have been obtained from a research interface into the patient database and so is not widespread but is being relaunched as a research tool for our users (PRIDE 2.0 is our reference to it) and allows software to directly communicate to the scanner interface and obtain images for processing. The standard DICOM export does present them differently.

However both methods are unlikely to consider this a big issue and will probably not be resolved in the short term as there is no requirement to present images in any defined order in the DICOM and while its annoying (and I wish it were different) we will see this as a defect in the recieveing software making unwarrented assumptions. WIth increasingly multidimensional data sets, there is also a question as even if a defined order was chosen, the assumptions are different between different softwares (some use Z,T,Type, others Type,Z,T etc. the only solution is to correctly read the dimension information and act on it. - it is an unfortunate consequence that DICOM is principally a communications protocol bewteen two databases and the dicom file is an attempt to freeze this communcation stream to disk but without any negotiation or handshaking. All the information is correct and available to allow the images to be correctly identified - this does not necessarily become memory inefficient aside from storing the dimension vectors somewhere prior to writing image data to disk.

Instance number is just that - a unique identifier for the image - indeed it was renamed this in the earlier DICOM standard from "image number" to imply the unlinking of any physical meaning for this value.

I'm not sure I would agree with the terms "in disarray" (unless you mean this very literally in that the images are not provided in an array format :-) )

Also - unfortunatley our provided in built nifti export is limited to certain image types and sequences due to the requirements of only producing a single nifti file, so splitting and mutli file types are currently not supported.

M

@neurolabusc
Copy link
Collaborator

neurolabusc commented Jul 30, 2021

@drmclem here are four objective criteria that suggest disarray to me:

  1. Bloat of redundant data in Philips enhanced DICOMs that impact disk usage and read performance (compare to Siemens and Canon implementations).
  2. Missing information to replicate sequences or analyze data
  3. Jumbled order of instance number used by most tools to sort data. One could argue the merits of XYZT or XYTZ, but at least choosing one consistent behavior is better than random assignment.
  4. Failure to use available public tags (e.g. Temporal Position Index 0020,9128 and In-Stack Position Number Attribute 0020,9057) that would help in determining the desired slice ordering (e.g. XYZT vs XYTZ).
  5. Failure to evolve data over a span of years to meet needs of users. Contrast this with the responsiveness of engineers at Canon, GE, Siemens, UIH to address limitations to meet the needs of their customers.

Please do not get my tone wrong, you have been outstanding at providing example data and describing properties of Philips DICOM data. However, it is discouraging that none of the issues are making their way up stream. Looking at the dcm2niix code shows a huge number of kludges have been inserted to handle Philips data. I worry that since these are being done in isolation from the manufacturer, they are inherently fragile. Consider the kludge discussed above for @baxpr - how will this behave for multi-echo data? Does Philips only use the private volumeNumber tag (2005,1063) for fMRI, or will this also work for diffusion data?

I strongly encourage Philips to continue and extend the development of their NIfTI export to handle the broad range of data the instruments can acquire.

  1. This can help end users directly by providing a robust, rapid and compact method for getting data.
  2. This will help Philips engineers recognize the meta-data that should be retained to allow accurate conversion.
  3. This will provide developers like me with a reference solution.

With regards to the comment we will see this as a defect in the receiveing software. I agree that the images are technically legal DICOM, but it seems irresponsible to generate data where 100% of the tools I know of are defective in interpreting it. The generating software knows the number of slices in the volume and the temporal order, so providing a meaningful value is simple for Philips, but can be challenging for receiving tools, in particular when values that help decode the information is stored in private tags that are not defined. It seems really inconsiderate (and pragmatically unreasonable) to assume all other tools should develop algorithms to pick up the pieces. It reminds me of the joke of the British driver who takes his car through the Chunnel to France, and while driving he hears a radio report that there is an individual driving on the motorway against traffic. He calls the radio station and says Its not one driver, it is all of them. If all other manufacturers, and all other receiving tools are doing things differently to you, perhaps you should change your behavior.

@drmclem
Copy link
Author

drmclem commented Jul 30, 2021

Hi

I think part of the problem is that there as aspects of DICOM which are not fully defined, and certainly public tags we are incorporating as we move forward in software releases and any errors in those we will fix as we move forward (4) - we are looking how best to populate the defined timing blocks in the light of multiband sequences and the like.

I would also argue that using the instance number as anything other than a unique identifier is error prone.

A comprehensive solution based on public tags which works for philips data should also work for the other vendors, but I also your relcutance to move such solutions from Philips specific to general given that "if not broken don't change it".

I think we also need to be careful of confusing public tag requirements and private tag requirements - there is much additional information in there to support recall and rescanning of the patient and our internal useage, which i suspect is behind 1 and 2. While these can be stripped out on PACS communication (requesting "Pure" dicom for example) the disk file export cannot state its requirements in advance through handshake, but equally some of what you may consider bloat, it vital to other or our own users.

Specifically for 3 - the background is with the move to parallel moving that cost onto the scanner is detrimental to image reconstruction performance during the scan and is not particularly an issue for PACS - and as we have seen, the standard DICOM File export has resolved much of this - (but relying on instance number is inherently fragile) - this new research interface does not (currently, we have reported the issue and we will see what happens).

For 5 - the prinicple user is the PACS community, we understand the the needs of the research community often needs internal information which is not currently defined in the DICOM standard and therfore requires private tags - but in general private tags (as they are intended only to be used by the vendor who created them) are difficult to request and inevitably these lag behind the research community needs but we do try be responsive to them but I would argue we have introduced a lot of tooling to support this for Philips users, just not in the DICOM-nii axis. Even for (3) "There are no official DICOM tags to report this information - on Siemens scanners we have reverse engineered the proprietary CSA header 16."

For your final point on Nii export from the tools - it is there but we need good reference and fixed standards to work with (currently based on https://nifti.nimh.nih.gov/nifti-1 and support fMRI and structural formats). It is a testement to your software that it has almost become by definition the reference implementaion (without even considering the BIDS component) and I think support multiple evolving standards is always difficult - I think our efforts are better placed in ensuring the DICOM implementation meets the needs.

For 3 and 4 - the problem is that a modern data set potentially is (x,y,z,temporal, image type (R,I,M,P), state (label,control for example) with some elements not unique (for example many z same t as for current dyanmics scans, or unique z - t for fmri with slice timing etc. and different users have different needs - if you are going to sort them anyway ...

On the specifics for the solution to the above ill put in a separate comment if thats ok.

M

@drmclem
Copy link
Author

drmclem commented Jul 30, 2021

Hi @baxpr - any chance you could share the dataset with me so I can have a look at it and compare the dataset which started this convo ?

@neurolabusc
Copy link
Collaborator

Your comment regarding Siemens actually supports my point. Note the move from the proprietary CSA header in the V* systems (Trio, Prisma) to public tags in the XA systems (Vida, Sola). Between XA10 and XA11 there was a constructive dialog with users regarding the meta data that would allow replication and comprehensive analysis.

@baxpr
Copy link

baxpr commented Jul 30, 2021

Thanks y'all - really appreciating this discussion which I'm learning a lot from. @drmclem I have sent you the images as well. The spatial and temporal ordering can be reconstructed from ImagePositionPatient and TemporalPositionIdentifier, but I gather the assumptions that dcm2niix usually makes for Philips images are not met.

@baxpr
Copy link

baxpr commented Jul 30, 2021

fwiw I banged out some python code to rewrite instance numbers so that dcm2niix can handle the DICOMs we have:

#!/usr/bin/env python

import argparse
import pydicom
import os
import pandas
import numpy

# Input arguments
parser = argparse.ArgumentParser()
parser.add_argument('--dicom_dir', required=True) 
parser.add_argument('--out_dir', required=True) 
args = parser.parse_args()
dicom_dir = args.dicom_dir
out_dir = args.out_dir
 
# Search DICOM directory and get all filenames
allf = None
for r,d,f in os.walk(dicom_dir):
    thisf = [os.path.join(r,x) for x in f]
    if not allf:
        allf = thisf
    else:
        allf = allf + thisf

# Get relevant field values
DicomFile = list()
TempPos = list()
ImageOrientationPatient = list()
PixelSpacing = list()
ImagePositionPatient = list()
InstNum = list()
for f in allf:
    try:
        ds = pydicom.dcmread(f,stop_before_pixels=True)
    except pydicom.errors.InvalidDicomError:
        print(f'Invalid DICOM {f} - skipping')
        continue
    DicomFile.append(f)
    TempPos.append(int(ds.TemporalPositionIdentifier))
    InstNum.append(int(ds.InstanceNumber))
    ImageOrientationPatient.append(numpy.array(ds.ImageOrientationPatient))
    PixelSpacing.append(numpy.array(ds.PixelSpacing))
    ImagePositionPatient.append(numpy.array(ds.ImagePositionPatient))

# Verify we only have one ImageOrientationPatient and PixelSpacing for all frames
if not all( [all(x==ImageOrientationPatient[0]) for x in ImageOrientationPatient] ):
    raise Exception('Differences in ImageOrientationPatient')
if not all( [all(x==PixelSpacing[0]) for x in PixelSpacing] ):
    raise Exception('Differences in PixelSpacing')

# Compute the signed slice position on the through-slice axis
# https://nipy.org/nibabel/dicom/dicom_orientation.html
# https://itk.org/pipermail/insight-users/2003-September/004762.html
SliceNormal = numpy.cross(ImageOrientationPatient[0][0:3],ImageOrientationPatient[0][3:6])
SlicePos = [numpy.dot(SliceNormal,x) for x in ImagePositionPatient]

# Combine to data frame
info = pandas.DataFrame({
    'DicomFile': DicomFile,
    'InstNum': InstNum,
    'TempPos': TempPos,
    'ImageOrientationPatient': ImageOrientationPatient,
    'PixelSpacing': PixelSpacing,
    'ImagePositionPatient': ImagePositionPatient,
    'SlicePos': SlicePos,
    })

# Sort according to time, then slice
info = info.sort_values(['TempPos','SlicePos'])

# Update instance numbers
info['NewInstNum'] = range(1,info.shape[0]+1)

print('')
print(info[['DicomFile','InstNum','NewInstNum','TempPos','SlicePos']])

# Check some things
uTempPos = numpy.unique(info['TempPos'])
numT = len(uTempPos)
uSlicePos = numpy.unique(info['SlicePos'])
numS = len(uSlicePos)
numF = info.shape[0]
print(f'Found {numT} time points, {numS} slices, {numF} frames')
if len(uTempPos)<2:
    raise Exception('Only 1 temporal position - not an fmri series?')
if numF != numT*numS:
    raise Exception(f'Time and slice not completely sampled - expected {numT*numS} frames')

# Edit dicoms with new instance number
os.makedirs(out_dir,exist_ok=True)
for row in info.itertuples():
    ds = pydicom.dcmread(row.DicomFile)
    ds.InstanceNumber = row.NewInstNum
    ds.save_as(os.path.join(out_dir,f'{int(ds.InstanceNumber):06}.dcm'))

@captainnova
Copy link
Collaborator

Hi,
I've been following this as someone involved in several multicenter research studies, and hate to pile on Philips with a me too, but I agree with the points Chris enumerated.

  1. I think this refers to the repetition of settings that could theoretically change from frame to frame, but in practice never do. As far as I am aware there is nothing in the scanner console interface that would let users make them vary within a scan without programming a custom sequence, in which case it should (ideally) be up to the sequence programmer to write things as necessary. I might be wrong about that, but in any case it's ridiculous for about half of a 300MB dMRI to be metadata. I've asked one of our local Philips representatives to request a fix for this, and he eventually had to come back with the bad news that the DICOM team was opposed to making this change. It was really disappointing - Philips had been considered innovators in the DICOM field for rolling out enhanced DICOM, in part because 1 header for all the slices should use less disk space than classic DICOM.
  2. Can we pretty please have PhaseEncodingDirection? Philips makes it so easy to control the fat shift direction on the console, and its omission in the DICOM header is very odd. Is it perhaps extractable from the examcard (XX.dcm)?
  3. and 4. I can understand InstanceNumber being fairly arbitrary, and the loop order depending on the implementation details of the sequence and reconstruction, but that is exactly why public tags clearly giving Z (or the Z equivalent) and T need to be stored.
  4. Siemens switching from mosaic to enhanced DICOM in XA should be considered a feather in Philips' cap. Migrating a LOT of info from private to public tags was also good, but I think all the manufacturers agree on that goal in theory if not always in practice. The opportunity here for Philips is that Siemens XA is making the public tags for some of the newer features more established, i.e. cutting through the uncertainty and hesitation over what to call these things. The only other important difference is the practical one already discussed in 1. Whether to store the scan in 1 file or 1 file per volume seems to be slightly controversial, which I think means people greatly prefer either to 1 file per slice. I certainly don't like it when slices go missing during transmission.
  5. The research community is increasingly open to, or looking forward to, scanners directly saving .niis and BIDS, at least in neuroimaging. The transition will be difficult for those with existing pipelines and especially databases, but those people also know very well how difficult it is to deal with DICOM on a daily basis. I do think there will be some complications in practice, since as @drmclem mentioned DICOM is a more general format than NIFTI, but multple NIFTIs would be OK. dcm2niix does that anyway, which is annoying when you are only expecting 1 .nii to pop out, but that's life. How does BIDS handle the multiple .niis per series case? 1 .json per .nii? Or one .json with a manifest of the .niis keyed by whatever prevented the .niis from being merged? I have a feeling that has not been worked out yet, and it won't get much traction until all the scans in a session can be written in NIFTI format without resorting to DICOM.

Best regards,

Rob

neurolabusc added a commit that referenced this issue Jul 31, 2021
@neurolabusc
Copy link
Collaborator

neurolabusc commented Jul 31, 2021

@baxpr and @captainnova do you want to try out the latest commit? It attempts to generalize @baxpr Python code (which was not designed for diffusion data). Here are details regarding the implementation and the implications:

  1. Volume number is derived from largest value of the following DICOM tags: TemporalPosition (0020,0100), GradientOrientationNumber (2005,1413), or VolumeNumber (2005,1063).
  2. As @baxpr, slice position derived signed slice position on the through-slice axis.
  3. Philips often stores derived TRACE image with raw diffusion series. These are given the same GradientOrientationNumber as the B=0 volume. With this commit, these are identified an identified as volume N, where N is the total number of volumes in the series. Therefore, these should be placed at the end of the series, identical to the behavior of dcm2niix.
  4. In my experience, Philips often saves the B=0 volume first to disk, but gives it the highest GradientOrientationNumber. Therefore, the development version of dcm2niix may reorder DTI volumes differently than prior versions.
  5. I have no examples of Philips DTI data where a single series contains multiple averages (either repeated instances of the B=0, or repeating a gradient direction). I am not sure if the GradientOrientationNumber is reused or unique for each average. If the numbers are reused, we will need to develop logic to disambiguate the different instances (e.g. if a b=0 volume is acquired at the start and end of a series, we do not want to mix and match slices from these two volumes, as there may be head motion between them).
  6. The current solution only works for classic DICOM data. Once the logic is clear, we can extend to enhanced datasets.
  7. This solution should work with single volume datasets, such as the one that elicited this issue (WIP_T2_brain_cor_20210709112826_1101). @drmclem I would be grateful if you could test this.
  8. There is a small performance penalty for testing slice ordering for data where the the original slices are ordered correctly. I think this is negligible, but the testing is restricted to Philips data so the penalty will not impact other manufacturers.
  9. The sorting is done after a series is segmented based on type (phase, magnitude, real, imaginary), echo number, and trigger delay time. This should prevent contamination across these factors.
  10. Since Philips does not provide slice timing information, and this method orders slices as reported by the image plane, users will want to ensure that the order they save slice timing values matches the order that this algorithm saves the data to disk.

In brief, while this method seems to provide a robust solution for the problem (as long as point 5 does not pose an issue), the ordering for DTI data may be different than previous versions. This makes regression testing tricky. Also, since Philips often stores the b=0 volume as the first to disk, but gives it the highest GradientOrientationNumber, there may be complications with scripts that assume that the first volume is the reference volume (e.g. if one does not specify the reference_no to eddy_correct, it assumes that the initial volume is the reference volume.

@baxpr
Copy link

baxpr commented Jul 31, 2021

Thanks! Yes, that works for my test fmri.

@neurolabusc I've sent you a diffusion scan that might help. We quit making assumptions about ordering of diffusion volumes in the niftis a couple years ago due to these issues. Philips does not permit multiple b=0 to be specified, as far as I know - the data set I just sent you is how we have worked around that using a negligible b value and some arbitrary specified direction.

... That is an enhanced DICOM though, generated in a different way. I'll see if I can get a classic set for you.

@captainnova
Copy link
Collaborator

Unfortunately it broke diffusion. A scan with 1 b=0 followed by 32 b=999, and an apparent ADC probably at the end came out with the b=0 volume distributed among the b=999 volumes 2 slices at a time, i.e. 2 of the slices were replaced by the b=0 slices of the correct z. The bumped slices were placed at the highest z.

The "ADC" file is mostly the average of the b=999 volumes (which is normal - it wasn't really the ADC to begin with), except that two of its slices have also been replaced with b=0 and moved to the top of the image.

The array shape, .bval, and .bvec files were not changed.

Although the ImageType says ["ORIGINAL", "PRIMARY", "M", "SE", "M", "SE"], I think this series should really have been called DERIVED, since the SeriesDescription and ProtocolName say "Reg - AX_DTI SENSE", i.e. it was the MoCo version of the scan.
I don't know if that's relevant, though - it's the first scan I checked out of a large number of discrepancies. The (good) decision to stop repunctuating strings in the JSON made before and after look different for almost everything in my test suite.

Rob

@captainnova
Copy link
Collaborator

captainnova commented Aug 1, 2021

I should have given the dim field from the nii header:

 regular               38      1    r
  dim_info              39      1    57
  dim                   40      8    4 176 176 60 33 1 1 1
  srow_x               280      4    -1.363636 0.0 -0.0 120.572365
  srow_y               296      4    -0.0 1.363636 -0.0 -96.28479
  srow_z               312      4    0.0 0.0 2.699997 -76.36232

But again, it's probably best not to focus on this particular one yet. All I know is it used to work.

@captainnova
Copy link
Collaborator

Good news - enhanced DICOM (but a different scan, that really is ORIGINAL data) seems to be unharmed. The other series had been unenhanced with dcuncat.

@neurolabusc neurolabusc reopened this Aug 2, 2021
@neurolabusc
Copy link
Collaborator

neurolabusc commented Aug 3, 2021

@captainnova

  • the development branch only changes the sorting of classic DICOM. For enhanced DICOMs, dcm2niix will use the public tag Dimension Index Values (0020,9157). Unfortunately, Philips does not provide this public tag for their classic DICOM images, so we need to develop an alternative heuristic.
  • From your email, the problematic image is data that was data originally saved as Philips Enhanced DICOM, and converted to classic format using dicom3tools dcuncat.
    • I was unable to replicate this problem. I used the most recent version of dcuncat (20210306100017) to convert enhanced DTI data from both old (Achieva 3.2.2\3.2.2.0) and relatively modern (Achieva dStream 5.3.0\5.3.0.3) Philips scanners to classic format. The data was successfully converted using the development branch of dcm2niix. The original datasets included GradientOrientationNumber (2005,1413) which was used for sorting the volumes.
    • I realize you are converting data from different clinical sites, and therefore you do not have much control over the handling of the data, nor is it easy with you to share them. However, it is worth emphasizing that the DICOM format is very complex, and has been interpreted differently by different manufacturers. Further, each manufacturer has evolved their implementation over time. Therefore, automated DICOM conversion requires a complex set of heuristics. Modifying the data can have unintended consequences. Teams need to be very cautious regarding the tools used to manipulate DICOM images, this includes data anonymization tools, enhanced to classic conversion or even having the image touched by some PACS systems (e.g. the invalid thumbnail of a GEIIS PACS, or the obfuscation of an AGFA or dcm4che PACS).
    • I urge users with Philips hardware at their sites to lobby the Philips employees they work with to remove the redundancy from their Enhanced DICOMs (using Siemens XA and Canon as a reference for clean design and support of public tags). Neuroscientists like me are not so impacted by the Philips format: the only cost is the increased disk usage and the one-off time penalty of converting these DICOMs to NIfTI. The impact is far worse for those on the clinical side where the excessive meta data makes inspection of the headers unwieldy. To understand the issue, (after stopping all critical operations on your computer) open a Philips Enhanced scan in a popular DICOM viewer like Horos or Osirix and choose to 'View meta-data of this image`. The sheer size of the enhanced DICOM header will overwhelm the viewer. Your software may crash if it exceeds the RAM on your computer (notice RAM usage massively exceeds that of simply viewing the image data), and if the meta data window is shown, scrolling through the immense amount of data is sluggish and the sheer amount of redundant data makes the task of identifying key tags taxing. I suspect this is the reason that many clinical sites like the one you are working with are resorting to using dcuncat on their images. If Philips can streamline their enhanced implementation, users will appreciate the benefits of this enhanced format instead of begrudging the bloated implementation. I do think that many users worry about raising concerns with their manufacturers, as they feel that complaining may harm their relationship. Since I do not have Philips systems at my center, I do not think my voice carries a lot of weight. However, in my experience with GE and Siemens, the engineers who develop these tools appreciate constructive criticism.

@tomaroberts
Copy link

Hi all,

I'm the researcher who has been working with @drmclem, whose DICOM data set triggered this issue. I've caught up with the thread this morning.

I've just cloned and run the development branch of dcm2niix following the adjustments by @neurolabusc.

Can confirm that it has fixed the issue for our DICOMs. So, thank you very much @neurolabusc for the fix and @baxpr for your initial Python implementation too.

neurolabusc added a commit that referenced this issue Aug 5, 2021
@neurolabusc
Copy link
Collaborator

Paul Morgan at Nottingham University has provided images for a new validation repository dcm_qa_philips_asl. Paul provided both the enhanced and classic images, but the combined set seems to overwhelm Github, so this repository is only the classic images. I have updated dcm2niix to handle these images. This should now deal with diffusion, ASL and fMRI data. The ASL data has a few quirks (as described with the validation), so the dcm2niix commit includes a lot of heuristics to handle that modality. Hopefully others can test this on a broader set of images.

The ordering of ASL volumes with respect to repeat, phase and presence of label is arbitrary, so I just copied the order observed for the enhanced images.

It would be great if Philips could provide a bit more ASL sequence details for these DICOMs. BEP005 provides a community wishlist that will allow automated analyses for this modality.

@captainnova
Copy link
Collaborator

I no longer think dcuncat is the likely reason why the development version of dcm2niix missorts the slices in some of our test dMRI series. The affected series are from Philips 2, 3, and 4 software, and have been through at least one deidentification step (dcanon by me, and possibly others before they were sent to our lab), and I have lost track of whether the scanner stored them in classic or enhanced format, if I ever knew. However, in the old days it was standard practice for people here to run dcuncat on enhanced DICOM as soon as it arrived, so unfortunately for the data I have in hand the scanner software version is strongly conflated with whether it has been dcuncatted.

A Philips 5 dMRI that we scanned here works much better, although development dcm2niix prints some spurious error messages when working with the original enhanced version. It looks like it has well behaved volume indices, while the "volume indices" for the older scans go from 1 to 65537 and do not work with the new scheme.

@neurolabusc
Copy link
Collaborator

I am closing this issue. The current development version comprehensively addresses this issue. In cases where all the useful tags are stripped the algorithm resorts to assuming that the instance number is truthful, resolving the regression for the (over aggressively) anonymized data from @captainnova .

However, this fixes raises deeper concerns regarding how Philips ASL volumes should be ordered. This problem is introduced in issue 533.

yarikoptic added a commit to neurodebian/dcm2niix that referenced this issue Sep 9, 2021
* commit 'v1.0.20200427-207-g66a4fd0': (76 commits)
  Philips ASL in temporal order (rordenlab#533)
  Optional compilation to disable kludge for issue 532, e.g. "make CFLAGS=-DmyDisableGEPEPolarFlip" (rordenlab#532)
  0020,9157 if dcuncat, use 0020,0013 if 0019,10A2; 0020,0100; 2005,1063; 2005,1413 do not vary (rordenlab#529)
  Add BIDS ArterialSpinLabelingType field
  Diagnostics for issue 529 (rordenlab#529)
  Flip row order for kGE_EPI_PEPOLAR_REV
  Detect and correct GE epi_pepolar sequence (rordenlab#532)
  More CT meta data, detect manufacturer MEDISO
  Philips ASL volume order (rordenlab#529)
  ensureSequentialSlicePositions verbosity
  Handle Philips images where 2005,1063 is set to zero for all volumes, see Magdeburg_2014 from dcm_qa_philips Philips MR 51.0
  @baxpr kludge (rordenlab#529)
  Remove redundancy
  Kludge for Philips random instance numbers (rordenlab#529)
  BUG: Fix preprocessor conflict.
  Support ancient Linux, tested on holy-build-box (rordenlab#531)
  undefine DT_NONE (https://neurostars.org/t/adjusting-for-negative-mosaicrefacqtimes-issue-271-warning-but-with-normal-slice-times/19866/3)
  UIH new TotalReadoutTIme formula (rordenlab#531)
  UIH uodates
  Allow acquisition number (0020,0012) in filename (rordenlab#526)
  ...
@neurolabusc
Copy link
Collaborator

@drmclem and @sandeepganji please take my tone as constructive. However, I do think that the Philips enhanced DICOM places a tremendous burden on users. I was originally lured to South Carolina by the opportunity to use the brilliant Philips MRI systems. Years later, our entire hospital system is moving exclusively to Siemens. I do think Philips continues to make outstanding hardware, but data storage, transfer and display is a major factor for clinical sales. As a consumer, I have a vested interest in seeing competitive alternatives.

I do urge Philips to examine their competitors modern enhanced DICOM implementations. Here is data from Canon saved as both enhanced and classic DICOM. The uncompressed classic data requires 458mb, the enhanced is just 76mb. Likewise, dcm2niix is able to convert the clean enhanced DICOM in about 1/10th the time of classic data. The meta data is human readable, instead of overwhelming humans and machines with redundant data.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

5 participants