rordenlab / dcm2niix

dcm2nii DICOM to NIfTI converter: compiled versions available from NITRC
https://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage
Other
899 stars 229 forks source link

dcm2niix does not support data from United Imaging Healthcare(UIH), China #225

Closed chunshan closed 6 years ago

chunshan commented 6 years ago

UIH is a noted MR manufacturer from China which is currently exploring the worldwide markets. By now many institutes over the world have been using UIH MRscanners to acquire DWI/DTI and fMRI data. The users also want to convert their dicom data to NIfTI format using dcm2niix which is very popular in the brain and neurology research community. But the dcm2niix currently does not recognize the UIH data format correctly. As an example, We tried to convert the fMRI data using dcm2niix and got the following result of improper format conversion.

UIH supports two ways of archiving the DWI/DTI and fMRI data. One way is one dicom file per slice and the other is one dicom file per volume (We call it GRID format). The private tags used in the images are shown in the following table.

Tag ID Tag Name VR VM Description Sample
0061,1002 Generate Private US 1 Flag to generate private format file 1
0061,4002 FOV SH 1 FOV(mm) 224*224
0065,1000 MeasurmentUID UL 1 Measurement UID of Protocol 12547865
0065,1002 ImageOrientationDisplayed SH 1 Image Orientation Displayed Sag or Sag>Cor
0065,1003 ReceiveCoil LO 1 Receive Coil Information H 8
0065,1004 Interpolation SH 1 Interpolation I
0065,1005 PE Direction Displayed SH 1 Phase encoding diretion displayed A->P or H->F
0065,1006 Slice Group ID IS 1 Slice Group ID 1
0065,1007 Uprotocol OB 1 Uprotocol value  
0065,1009 BActualValue FD 1 Actual B-Value from sequence 1000.0
0065,100A BUserValue FD 1 User Choose B-Value from UI 1000.0
0065,100B Block Size DS 1 Size of the paradigm/block 10
0065,100C Experimental status SH 1 fMRI rest/active
0065,100D Parallel Information SH 1 ratio of parallel acquisition and acceleration  
0065,100F Slice Position SH 1 Slice location displayed on the screen H23.4
0065,1011 Sections SH 1    
0065,1013 InPlaneRotAngle FD(°) 1 Rotation angle in the plane -0.5936
0065,1014 SliceNormalVector DS 3 Normal vector of the slice 0\0\1
0065,1015 SliceCenterPosition DS 3 Center position of the slice 0\0\0
0065,1016 PixelRotateModel UL 1 Pixel Rotation Model 4
0065,1017 SAR Model LO 1 Calculation model of SAR value Normal:WHBST
0065,1018 dB/dt Model LO 1 Calculation model of dB/dt Normal
0065,1023 TablePosition LO 1 Table Position 0
0065,1025 Slice Gap DS 1 Slice Gap 0.0
0065,1029 AcquisitionDuration SH 1 Acquisition Duration 0.03
0065,102B ApplicationCategory LT 1 Application names available DTI\Func
0065,102C RepeatitionIndex IS 1   0
0065,102D SequenceDisplayName ST 1 Sequence display name Epi_dti_b0
0065,102E NoiseDecovarFlag LO 1 Noise decorrelation flag PreWhite
0065,102F ScaleFactor FL 1 scale factor 2.125
0065,1031 MRSequenceVariant SH 1 SequenceVariant  
0065,1032 MRKSpaceFilter SH 1 K space filter  
0065,1033 MRTableMode SH 1 Table mode Fix
0065,1036 MRDiscoParameter OB 1    
0065,1037 MRDiffusionGradOrientation FD 3 Diffusion gradient orientation 0\0\0
0065,1038 MRPerfusionNoiseLevel FD 1 epi_dwi/perfusion noise level 40
0065,1039 MRGradRange SH 6 linear range of gradient 0.0\157\0.0\157\0.0\125
0065,1050 MR Number Of Slice In Volume DS 1 Number Of Frames In a Volume,Columns of each frame: cols =ceil(sqrt(total)) ; Rows of each frame: rows =ceil(total/cols) ; appeared when image type (00080008) has VFRAME 27
0065,1051 MR VFrame Sequence SQ 1 1  
->0008,0022 Acquisition Date DA 1    
->0008,0032 Acquisition Time TM 1    
->0008,002A Acquisition DateTime DT 1    
->0020,0032 ImagePosition(Patient) DS 3    
->00201041 Slice Location DS 1    
->0018,9073 Acquisition Duration FD 1    
->0065,100C MRExperimental Status SH 1   rest/active

The example data can be downloaded from https://1drv.ms/f/s!Avf7THyflzj1gnO37GL2I8Hk-0MV . The gradient table for the DTI data is also provided.

If any questions, please give me a reply or contact me via email chunshan.yang@united-imaging.com

neurolabusc commented 6 years ago

I have uploaded a new patch that adds experimental support for some UIH features. I would suggest the engineers at UIH make a fork of dcm2niix and fill out support for the non-DICOM aspects of their format. Once complete, a pull request can be made that can insert the UIH-specific code into the main branch.

At the moment, the format seems underspecified, and the engineers can resolve this with their patches. A few initial comments:

  1. The grid format is underspecified Consider when 0065,1050 is 20 - is this saved as a 4x5 grid, a 5x4 grid or a 5x5 grid. At the moment my software uses the default Siemens expansion, but it is unclear if UIH uses the same algorithm.
  2. The format is underspecified with regards to the image coordinates of the Grid images. The current code has the comment "UIH Grid to matrix transform unknown\n" - UIH engineers will need to insert the correct transform here (the samples demonstrate that this is not identical to the Siemens scheme, but the code for Siemens provides a clear template for a UIH specific solution).
  3. It is unclear how diffusion vectors are stored - with respect to the imaging plane or the scanner bore? Are the values never stored in the DICOM files, and if they are the format needs to be specified.
neurolabusc commented 6 years ago

I have ported Xiangrui Li's modifications for UIH support. This provides experimental support for de-interlacing grids, slice timing, diffusion gradient directions, and phase polarity. As I do not have access to UIH hardware, I encourage users to test this carefully.

chunshan commented 6 years ago

@neurolabusc ,For your comments,

• Is 0065,102F ScaleFactor identical to 0028,1053? Which gets priority if they differ. [UIH-MR] The two tags are different. (0065,102F) is compensation for coil or sequence, only for UIH engineers’ internal use. (0028,1053) is the scaling factor between display value and stored value, as written in DICOM standard definition.

• Is 0065,1025 identical to the difference between 0018,0050 and 0018,0088? Which gets priority if they differ. [UIH-MR] (0065,1025) indicates gap between slices, which doesn’t calculate slice thickness. Generally, spacingbetweenslices(0018,0088)=MRSliceGap(0065,1025) + thickness(0018,0050)

• How is phase encoding polarity stored (e.g. A->P vs P->A)? this is important for tools like TOPUP. [UIH-MR] The phase encoding polarity is stored in (0065,1005),such as A->P,H->F,R->L.

• How is the slice timing for 2D slices of a 3D EPI image stored - e.g. how does one determine ascending, descending, interleaved. How does one determine if a temporal gap is inserted between the last slice of one volume and the first slice of the next 3D volume (e.g. sparse imaging, dynamic correction, etc)? [UIH-MR] Acquisition timing of each 2D slice is written in tag Acquisition Time(0008,0032). With reference to tag Image Position Patient (0020,0032), we can determine if the acquisition is ascending, descending or interleaved. For calculating temporal gap inserted between the last slice of one volume and the first slice of the next 3D volume, tag (0018, 9073) Acquisition Duration can be used.

• The Grid layout is in theory similar to the deprecated Siemens Mosaic and has the same benefits and weaknesses. I strongly suggest UIH consider moving to the enhanced DICOM as illustrated by the Siemens XA images. Neither the grid or mosaic formats are truly DICOM conformant, and provide incorrect coordinates with valid DICOM viewing software. In addition, the Grid/Mosaic formats impose a performance penalty for de-interlacing [UIH-MR] Thanks for the suggestions. We have currently acquired large amount of data in GRID format. We hope it can be support at current stage. Will consider moving to enhanced MR DICOM in the future.

Another 3 questions:

  1. The grid format is underspecified Consider when 0065,1050 is 20 - is this saved as a 4x5 grid, a 5x4 grid or a 5x5 grid. At the moment my software uses the default Siemens expansion, but it is unclear if UIH uses the same algorithm. [UIH-MR] Our algorithm is slightly different from Siemens. Our number of columns for each frame: cols=ceil(sqrt(total)). Number of rows for each frame: rows= ceil(total/cols number); if 065,1050 is 20, the frame will be in 4*5 grid.

  2. The format is underspecified with regards to the image coordinates of the Grid images. The current code has the comment "UIH Grid to matrix transform unknown\n" - UIH engineers will need to insert the correct transform here (the samples demonstrate that this is not identical to the Siemens scheme, but the code for Siemens provides a clear template for a UIH specific solution). [UIH-MR] Please refer to the above answer.

  3. It is unclear how diffusion vectors are stored - with respect to the imaging plane or the scanner bore? Are the values never stored in the DICOM files, and if they are the format needs to be specified. [UIH-MR] Current private Tag 0065,1037 stores the gradient direction of each image in patient coordinates, x\y\z.

chunshan commented 6 years ago

@xiangruili, for your comments in the email, I also want to give our comments as below. Thanks very much for all your comments. @neurolabusc @xiangruili

  1. The phase encoding polarity is stored in (0065,1005), something like A->P. Very good.

  2. Slice timing is stored as timestamp in MRVFrameSequence for each mosaic frame, or in regular tag of regular dicom. This serves the purpose, but I would suggest to save this as a dicom tag like Siemens mosaic file does. It has nSlices entries in milliseconds since the start of the TR. I believe this will be easier to both the UIH engineer and converter developers. [UIH-MR] We have used this method to store time information and collected large amount of data. We hope this approach can be supported at current stage.

  3. The /VFRAME in ImageType is the flag for mosaic format, right? [UIH-MR] YES

  4. The /DIFFUSION is missing in ImageType. I am now using BActualValue as a flag for DTI dicom. This should be fixed. [UIH-MR] UIH uses 0065,102b Image Category to illustrate the type of images. For fMRI, this tag will contain BOLD. For DTI, this tag will contain DTI.

  5. B value and vector are stored in BActualValue and MRDiffusionGradOrientation. I suppose the vector is in dicom reference like Siemens and Philips, but please confirm. The picture below shows the result, and both lines and color display look right to me. I wonder why UIH stores these in private tags. There are standard dicom tags for both B value and vector. [UIH-MR] For compatibility, we will write BActualValue and MRDiffusionGradOrientation into both public and private tags.

  6. VFRAME or mosaic format. As Chris wrote, we expect this is similar to Siemens mosaic. There should be always n by n grids, where n is ceil(sqrt(nSlices)). Please confirm this. As converter developers, we also like the ImagePositionPatient to be consistent with Siemens, to simplify the special cases. But Siemens's treatment to put mosaic center as slice center was unfortunate. In my opinion, UIH method to keep original ImagePositionPatient makes more sense, although both treatment are not dicom compliant, as Chris already pointed out. [UIH-MR] Our algorithm is slightly different from Siemens. Our number of columns for each frame: cols=ceil(sqrt(total)). Number of rows for each frame: rows= ceil(total/cols number); if 065,1050 is 20, the frame will be in 4*5 grid. ImagePositionPatient is indicating the upper left position of the image, which is a more common approach as defined in DICOM standard.

  7. It is a very nice idea to have SliceNormalVector tag. This would avoid a common problem for converters to determine the slice stacking direction. Unfortunately, it seems to me the vector has wrong sign, indicating how confusing this might be. Please correct this if I am right. [UIH-MR] We suggest ignore this private tag when converting, Use public tags of Image Position Patient and Image Orientation Patient for ordering slices instead.

neurolabusc commented 6 years ago

@chunshan thanks for the clarifications.

Would it be possible for your team to acquire a reference dataset that demonstrates these features. Ideally, it would be great to have an open source dataset that I (and any other developer) could use like the Siemens dcm_qa where every commit is validated against this dataset using Travis/AppVeyor. An ideal dataset would have the following features:

Thanks for your efforts and clarifications. This should allow dcm2niix, dicm2nii and others to provide robust support for these images.

xiangruili commented 6 years ago

@chunshan When you have a chance to collect a dataset, it will be useful to include a short scan with sparse-like imaging, so we can check AcquisitionDuration tag. For the images you provided, I don't see this tag, supposing whole duration of RepetitionTime is used for acquisition?

chunshan commented 6 years ago

@xiangruili In the images (BOLD and DTI) we provide, AcquisitionDuration Tag (0018, 9073) (VR - FD) is provided. In the single slice file, this tag is the first level. In the GRID file, this tag is embedded in the tag (0065, 1051), which means MR VFrame Sequence.

@neurolabusc We are planning to acquire the reference dataset. Is the data added into the repository dcm_qa ?

neurolabusc commented 6 years ago

@chunshan if you agree, I would like to create a new repository named dcm_qa_uih that would provide a reference dataset. Like the current dcm_qa (Siemens) and dcm_qa_nih (GE), these would be automatically used to detect regressions in new dcm2niix commits, and will be openly available for other developers.

One more feature request for a sample dataset:

xiangruili commented 6 years ago

@chunshan RE: AcquisitionDuration I found the tag, but with several issues.

First, since both standard and private tags are used for this parameter, the private one should have a different name, otherwise the later tag will overwrite earlier one. I am now calling the private one AcquisitionDurationUIH.

Second, the meaning of the parameter is unclear to me. For the example of BOLD_resting scan, it is 34 or 33 (AcquisitionDurationUIH='0.03'). I guess it is in ms (DICOM definition is a mess for the unit of this). For BOLD_resting sequence, TR=2000, nSlice=33. My understanding is that AcquisitionDuration should be close to 2000/33=60.6ms for a slice. This can be verified by AcquisitionTime (~60ms difference). What does 34 or 33 for AcquisitionDuration mean here?

Last, AcquisitionDuration is supposed to be the same for all slices, so it seems to me it does not belong to MRVFrameSequence.

chunshan commented 6 years ago

@xiangruili , For the questions you mentioned, 1, We recommend you use the public tag(0018,9073). These two tags have the meaning but for different scenarios. As you know, the unit of the private tag is ms, and the unit of the private tag is s. The private tag is designed for the viewer to display the corner information directly.

2, Acquisition duration in UIH MR is the acquisition time window from the first line in k-space to the last line while the time of RF, gradient and others is not included. So you could notice that the acquisition duration is shorter than the time TR/nSlice. Which value is used in your correction?

3, Theoretically, acquisition duration is the same for all slices. But if we examine the value recorded in the dicom file, there is minor difference between slices (about 1~2ms).

chunshan commented 6 years ago

@neurolabusc Data acquired for qa has been uploaded. The gradient tables are also provided. Please help to check it and create a repository. https://1drv.ms/f/s!Avf7THyflzj16SYrIQWziLOvQ7LY

neurolabusc commented 6 years ago

Thanks for the example images. In general these are very clear, and you have provided a nice, compact dataset.

Unfortunately, these seem to report phase-encoding axis but not phase-encoding polarity. Specifically, consider series 6 (dti_tra_dir16_PA) and series 9 (dti_tra_dir16_AP). Both report (0065,1005) SH [A->P]. I had expected the PA file to report P->A, but this is not the case. I did not see any other DICOM elements that encode the fact that these two series have the same phase encoding axis, but opposite polarity from each other. Knowing the polarity is useful for image undistortion methods and both axis and polarity are required to determine the BIDS tag PhaseEncodingDirection. chunshan , Can you tell me if there is any way to distinguish A->P from P->A files based on the DICOM header?

chunshan commented 6 years ago

@neurolabusc

Thank you for the reply. Your understanding is correct about the phase- encoding direction tag of UIH. Currently, (0065,1005) only reports the phased-encoding axis without the polarity information. UIH uses two parameters on user interface to control the phase-encoding axis and polarity separately. Unfortunately, only axis info is recorded in dicom header at this point. We plan to save the polarity info in a new tag soon. Besides, if there is any other missing information in dicom header necessary for data processing, please let us know ASAP. Thanks again.

xiangruili commented 6 years ago

@chunshan
A new tag for polarity is fine, but something like A->P and P->A in the existing tag, as Chris mentioned, is preferred. This will be consistent with as least Philips, and it also corrects the possibly wrong polarity in A->P tag. Actually we need this tag only for polarity, since the axis is already stored in public tag InPlanePhaseEncodingDirection.

xiangruili commented 6 years ago

FYI @chunshan Siemens has no tag such as AP or PA, but uses a private tag PhaseEncodingDirectionPositive to indicate the polarity for ROW or COL in InPlanePhaseEncodingDirection.

xiangruili commented 6 years ago

@chunshan RE: AcquisitionDuration I was figuring out if it can be used for "Effective EPI echo spacing" called by FSL. It is the time between successive k-space lines, taking the parallel imaging acceleration into account if applicable.

You stated: "For calculating temporal gap inserted between the last slice of one volume and the first slice of the next 3D volume, tag (0018, 9073) Acquisition Duration can be used" and "Acquisition duration in UIH MR is the acquisition time window from the first line in k-space to the last line while the time of RF, gradient and others is not included"

I don't know how to intemperate the first statement. From the second statement, it seems AcquisitionDuration / dimAtPhaseDirection is "Effective EPI echo spacing" we wanted for distortion correction. Could @chunshan and @neurolabusc confirm or correct me?

xiangruili commented 6 years ago

@neurolabusc RE: bvec with larger rotation My quick test on series dti_tra_dir16_PA_rot_SaveBySlc indicates something is wrong, as shown in the image below. Have you tested it?

The bvec from dicm2nii.m is different from Grad_PA_rot.csv, and neither gives correct result.

image

neurolabusc commented 6 years ago

@xiangruili 1) I agree it would be nice to calculate "Effective EPI echo spacing" for FSL. I do have two small comments. First, for the typical study that uses identical sequences just with reversed phase encoding polarity the value for this only influences the calibrated scan used for visualization. If the values are wrong, the calibration numbers will be higher or lower but the warping will be the same and the end product is identical. So for most studies this value is nice but not critical. Second, to calculate this one typically needs to know a lot of values: bandwidth, EPI lines, partial Fourier and in-plane acceleration. @mharms did a great job with ensuring we can compute these with Siemens and created both good sample data and an excel spreadsheet.

2) Looking at the colors is a challenge when the colors are shown with respect to the image slice but the head is rotated in plane. For these images, I much prefer using FSLeyes or FSLview to show the vector lines. These look pretty good to my eyes. Please see attached view for this dataset:

dti

xiangruili commented 6 years ago

Thank @neurolabusc .

RE: Effective echo spacing I knew for topup, the readout time is not critical. Good to know the same for fieldmap correction. Then we don't have to be picky about its absolute value.

RE: bvec with rotation That was silly me :). I wrote my nii_viewer, but forgot the color is along display axis.

@chunshan It seems the bvec provided in csv files has some sign problem. This is not an issue since our converters will save bvec files based on dicom info. The sign depends on how nifti image is saved (e.g. left-handed or not), so in theory, a standalone bvec table for nifti is not meaningful.

mharms commented 6 years ago

While the value of EffectiveEchoSpacing may not be critical to getting the unwarping correct in the context of 'topup' applied using blip up/down scans, I would caution against including it if you aren't confident that the value is an accurate representation of the actual scan. My perspective is that users should be able to assume that actual field maps (e.g., in Hz) derived from blip up/down scans are accurate, and that is only true if EffectiveEchoSpacing is accurately specified.

neurolabusc commented 6 years ago

@chunshan

chunshan commented 6 years ago

@neurolabusc @xiangruili Thanks for giving us so much constructive advice and pointing out the errors of the tags.

RE:Acquisition Duration @xiangruili , I am sorry that the the first state is not correct.

RE: Phase Encoding Polarity I will discuss how to record this information in the DICOM tags.

I am pleased to colse this issue and to make the tags more reasonable.

xiangruili commented 6 years ago

@neurolabusc UIH group seemed unsure about the bvec reference (image or patient coordinate system, PCS). I thought it is in PCS like Siemens/Philips, but they suggest it is in image reference. I did dtifit again for the two in-plane rotated dataset, and the results look much better to my eyes if I treat the bvec as image reference. The .bvec file will have opposite sign at y-axis compared to bvec in dicom if we save left-handed nifti for these two dataset. Could you verify the image reference?

I just updated dicm2nii for this in case you want to compare.

In my opinion, PCS reference is deterministic, while the image reference may cause confusion for signs for the cases with different phase encoding direction, non-axial slices, and reversed slice order in mosaic etc.

neurolabusc commented 6 years ago
dcm2niix_uih
xiangruili commented 6 years ago

@neurolabusc RE: bvec sign relation to qform determinant Because of the special treatment by FSL, dicm2nii.m default chooses left-handed storage. I used positive determinant long time ago, and it worked with FSL 5.0.8 (the version may not be accurate). It seemed that early FSL did not care the determinant, and worked for both +/- determinant (sound a little strange though). Later FSL seems strict on this. I remember you told me I need to reverse x sign with right-handed storage. So for the current FSL, the answer is yes, we need to reverse bvec x sign if the determinant is positive (causing inconsistency for sign and voxel indexing, which may have problem with other packages). The negative determinant should work for all I think.

issakomi commented 6 years ago

@chunshan

FYI, applicable to GRID format DICOM files:

Private SQ items with private data elements are missing private creator inside http://dicom.nema.org/medical/dicom/current/output/chtml/part05/sect_7.8.html

0054,0081 (Number of Slices) is available, but MR Image Module, IMHO, can not contain "Number of Slices" http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.8.3.html#sect_C.8.3.1

Image Plane Module missing Image Position (Patient) 0020,0032 (there are multiple positions in private SQ) http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.7.6.2.html

Patient Age 0010,1010 - VR AS has invalid character (*) http://dicom.nema.org/medical/dicom/current/output/chtml/part05/sect_6.2.html

BTW, there is e.g. Enhanced MR Image IOD with MR Diffusion Functional Group Macro, might be it were an option?

CC @dclunie @malaterre

Here is the above link again for simplicity, s. folder epi_dti_tra_dir64_SaveBySlc https://1drv.ms/f/s!Avf7THyflzj1gnO37GL2I8Hk-0MV

Regards

neurolabusc commented 6 years ago

@chunshan Can you tell me how one can detect if in plane acceleration (e.g. SENSE/GRAPPA) and/or partial Fourier was applied? As an example, consider the series dti_tra_dir16_PA_SaveBySlc__135612you provided. This reports:

  (0018,0091) IS [32]       #   2, 1 EchoTrainLength
  (0018,1310) US 64\0\0\64  #   8, 4 AcquisitionMatrix

I assume this means that the image was acquired with SENSE x2 (as partial Fourier might be none, 7/8, 6/8, 5/8 but NOT 4/8). However, for the BIDS files it is important to distinguish all combinations of parallel imaging and partial Fourier.

(0065,100d) SH [F:2S]

chunshan commented 6 years ago

@neurolabusc I asked some other colleagues about the acceleration information. They told me that the in plane acceleration factor is recorded in (0065,100d) while the partial Fourier acceleration is not recorded in the DICOM header.