diff --git a/COMPILE.md b/COMPILE.md index 53cb8265..fec61b75 100644 --- a/COMPILE.md +++ b/COMPILE.md @@ -137,34 +137,26 @@ lipo -create dcm2niix32 dcm2niix64 -o dcm2niix file ./dcm2niix ``` -##### OSX GRAPHICAL INTERFACE BUILD - -You can building the OSX graphical user interface using Xcode. First, Copy contents of "console" folder to /xcode/dcm2/core. Next, open and compile the project "dcm2.xcodeproj" with Xcode 4.6 or later - -##### THE QT AND wxWIDGETS GUIs ARE NOT YET SUPPORTED - FOLLOWING LINES FOR FUTURE VERSIONS - -Building QT graphical user interface: - Open "dcm2.pro" with QTCreator. This should work on OSX and Linux. On Windows the printf information is not redirected to the user interface - If compile gives you grief look at the .pro file which has notes for different operating systems. - -Building using wxWidgets -wxWdigets makefiles are pretty complex and specific for your operating system. For simplicity, we will build the "clipboard" example that comes with wxwidgets and then substitute our own code. The process goes something like this. - a.) Install wxwdigets - b.) successfully "make" the samples/clipboard program - c.) DELETE console/makefile. WE DO NOT WANT TO OVERWRITE the WX MAKEFILE - d.) with the exception of "makefile", copy the contents of console to /samples/clipboard - e.) overwrite the original /samples/clipboard.cpp with the dcm2niix file of the same name - f.) Older Xcodes have problems with .cpp files, whereas wxWidgets's makefiles do not compile with "-x cpp". So the core files are called '.c' but we will rename them to .cpp for wxWidgets: - rename 's/\.c$/\.cpp/' * - g.) edit the /samples/clipboard makefile: Add "nii_dicom.o nifti1_io_core.o nii_ortho.o nii_dicom_batch.o \" to CLIPBOARD_OBJECTS: -CLIPBOARD_OBJECTS = \ - nii_dicom.o nifti1_io_core.o nii_ortho.o nii_dicom_batch.o \ - $(__clipboard___win32rc) \ - $(__clipboard_os2_lib_res) \ - clipboard_clipboard.o - h.) edit the /samples/clipboard makefile: With wxWidgets we will capture std::cout comments, not printf, so we need to add "-DDmyUseCOut" to CXXFLAGS: -CXXFLAGS = -DmyUseCOut -DWX_PRECOMP .... - i.) For a full refresh -rm clipboard -rm *.o -make +##### CMAKE INSTALLATION + +While most of this page describes how to use `make` to compile dcm2niix, `cmake` can automatically aid complex builds. The [home page](https://github.com/rordenlab/dcm2niix) describes typical cmake options. The cmake command will attempt to pull additional code from git as needed for zlib, OpenJPEG etc. If you get the following error: + +``` +fatal: unable to connect to github.com: +github.com[0: 140.82.121.4]: errno=Connection timed out +``` + +This suggests git is unable to connect using ssh. You have two options, first you can disable the cmake option USE_GIT_PROTOCOL (which is on by default). Alternatively, to use https instead using the following lines prior to running cmake: + +``` +git config --global url."https://github.com/".insteadOf git@github.com: +git config --global url."https://".insteadOf git:// +``` + +Once the installation is completed, you can revert these changes: + +``` +git config --global --unset-all url.https://github.com/.insteadof +git config --global --unset-all url.https://.insteadof +``` + diff --git a/Toshiba/README.md b/Canon/README.md similarity index 59% rename from Toshiba/README.md rename to Canon/README.md index 8e36949d..fdf51e98 100644 --- a/Toshiba/README.md +++ b/Canon/README.md @@ -1,14 +1,22 @@ ## About -dcm2niix attempts to convert Toshiba DICOM format images to NIfTI. This page notes vendor specific conversion details. +dcm2niix can convert Canon (né Toshiba) DICOM format images to NIfTI. This page notes vendor specific conversion details. ## Diffusion Weighted Imaging Notes -In contrast to several other vendors, Toshiba uses public tags to report diffusion properties. Specifically, [DiffusionBValue (0018,9087)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9087)) and [DiffusionGradientOrientation (0018,9089)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9089)). Be aware that these tags are only populated for images where a diffusion gradient is applied. Consider a typical diffusion series where some volumes are acquired with B=0 while others have B=1000. In this case, only the volumes with B>0 will report a DiffusionBValue. +In contrast to several other vendors, Toshiba used public tags to report diffusion properties. Specifically, [DiffusionBValue (0018,9087)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9087)) and [DiffusionGradientOrientation (0018,9089)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9089)). Be aware that these tags are only populated for images where a diffusion gradient is applied. Consider a typical diffusion series where some volumes are acquired with B=0 while others have B=1000. In this case, only the volumes with B>0 will report a DiffusionBValue. + +Since the acquisition by Canon, these public tags are no longer populated. The diffusion gradient directions are now stored in the ASCII Image Comments tag. Like GE (but unlike [Siemens, GE and Toshiba](https://www.na-mic.org/wiki/NAMIC_Wiki:DTI:DICOM_for_DWI_and_DTI)), these directions are with respect to the image space, not the scanner bore. For empirical data see the Sample Datasets section. + +``` + (0018,9087) FD 1500 # 8, 1 DiffusionBValue + (0020,4000) LT [b=1500(0.445,0.000,0.895)] # 26, 1 ImageComments +``` + ## Unknown Properties -The (BIDS format)[https://bids.neuroimaging.io] can record several sequence properties that are useful for processing MRI data. The DICOM headers created by Toshiba scanners are very clean and minimalistic, and do not report several of these advanced properties. Therefore, dcm2niix is unable to populate these properties of the JSON file. This reflects a limitation of the DICOM images, not of dcm2niix. +The [BIDS format](https://bids.neuroimaging.io) can record several sequence properties that are useful for processing MRI data. The DICOM headers created by Toshiba scanners are very clean and minimalistic, and do not report several of these advanced properties. Therefore, dcm2niix is unable to populate these properties of the JSON file. This reflects a limitation of the DICOM images, not of dcm2niix. - SliceTiming is not recorded. This can be useful for slice time correction. - Phase encoding polarity is not record. This is useful for undistortion with [TOPUP](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/topup). @@ -18,3 +26,4 @@ The (BIDS format)[https://bids.neuroimaging.io] can record several sequence prop - [Toshiba Aquilion](https://www.aliza-dicom-viewer.com/download/datasets). - [Toshiba 3T Galan Diffusion Dataset](https://github.com/neurolabusc/dcm_qa_toshiba). + - [Canon 3T Galan Diffusion Dataset](https://github.com/neurolabusc/dcm_qa_canon). \ No newline at end of file diff --git a/ERRORS.md b/ERRORS.md new file mode 100644 index 00000000..89d270b3 --- /dev/null +++ b/ERRORS.md @@ -0,0 +1,29 @@ +## About + +dcm2niix will return an exit status to allow scripts to determine if a conversion was successful. Following Unix convention, the value value 0 is used for [EXIT_SUCCESS](https://www.gnu.org/software/libc/manual/html_node/Exit-Status.html). In contrast, any non-zero result suggests an error. In the Unix terminal and with shell scripts, the variable `$?` reports the most recent exit status. Here is an example of a successful conversion: + +``` +>dcm2niix ~/dcm +Chris Rorden's dcm2niiX version v1.0.20200331 Clang11.0.0 (64-bit MacOS) +Found 2 DICOM file(s) +Convert 2 DICOM as ~/dcm/dcm_ax_asc_6 (64x64x35x2) +Conversion required 0.015866 seconds (0.012676 for core code). +>echo $? +0 +``` + +Below is a list of possible return values from running dcm2niix. + +| Exit Status | Meaning | +| ----------- | ----------------------------------------------------------- | +| 0 | Success | +| 1 | Unspecified error (see console output for details) | +| 2 | No DICOM images found in input folder | +| 3 | Exit from report version (result of `dcm2niix -v`) | +| 4 | Corrupt DICOM file (Irrecoverable error during conversion) | +| 5 | Input folder invalid | +| 6 | Output folder invalid | +| 7 | Unable to write to output folder (check file permissions) | +| 8 | Converted some but not all of the input DICOMs | +| 9 | Unable to rename files (result of `dcm2niix -r y ~/in`) | + diff --git a/FILENAMING.md b/FILENAMING.md index 6ebb2370..4b79224c 100644 --- a/FILENAMING.md +++ b/FILENAMING.md @@ -19,7 +19,7 @@ You request the output file name with the `-f` argument. For example, consider y - %m=manufacturer short name (from 0008,0070: GE, Ph, Si, To, UI, NA) - %n=name of patient (from 0010,0010) - %o=mediaObjectInstanceUID (0002,0003)* - - %p=protocol name (from 0018,1030) + - %p=protocol name (from 0018,1030). If 0018,1030 is empty, or if the Manufacturer (0008,0070) is GE with Modality (0008,0060) of MR, then the SequenceName (0018,0024) is used if it is not empty. - %r=instance number (from 0020,0013)* - %s=series number (from 0020,0011) - %t=time of study (from 0008,0020 and 0008,0030) @@ -54,6 +54,10 @@ Some post-fixes are specific to Philips DICOMs If you do not want post-fixes, run dcm2niix in the terse mode (`--terse`). In this mode, most post-fixes will be omitted. Beware that this mode can have name clashes, and images from a series may over write each other. +## Overlays + +DICOM images can have up to [16](https://www.medicalconnections.co.uk/kb/Number-Of-Overlays-In-Image/) binary (black or white) overlays as described by the [Overlay Plane Module](http://dicom.nema.org/dicom/2013/output/chtml/part03/sect_C.9.html). dcm2niix will save these regions of interest with the post-fix "_ROIn" where N is the overlay number (1..16). + ## File Name Conflicts dcm2niix will attempt to write your image using the naming scheme you specify with the '-f' parameter. However, if an image already exists with the specified output name, dcm2niix will append a letter (e.g. 'a') to the end of a file name to avoid overwriting existing images. Consider a situation where dcm2niix is run with '-f %t'. This will name images based on the study time. If a single study has multiple series (for example, a T1 sequence and a fMRI scan, the reulting file names will conflict with each other. In order to avoid overwriting images, dcm2niix will resort to adding the post fix 'a', 'b', etc. There are a few solutions to avoiding these situations. You may want to consider using both of these: @@ -78,4 +82,6 @@ dcm2niix will attempt to write your image using the naming scheme you specify wi * (asterisk) ``` +[Control characters](https://en.wikipedia.org/wiki/ASCII#Control_characters) like backspace and tab are also forbidden. + Be warned that dcm2niix will copy all allowed characters verbatim, which can cause problems for some other tools. Consider this [sample dataset](https://github.com/neurolabusc/dcm_qa_nih/tree/master/In/20180918GE/mr_0004) where the DICOM Protocol Name (0018,1030) is 'Axial_EPI-FMRI_(Interleaved_I_to_S)'. The parentheses ("round brackets") may cause other tools issues. Consider converting this series with the command 'dcm2niix -f %s_%p ~/DICOM' to create the file '4_Axial_EPI-FMRI_(Interleaved_I_to_S).nii'.If you now run the command 'fslhd 4_Axial_EPI-FMRI_(Interleaved_I_to_S).nii' you will get the error '-bash: syntax error near unexpected token `(''. Therefore, it is often a good idea to use double quotes to specify the names of files. In this example 'fslhd "4_Axial_EPI-FMRI_(Interleaved_I_to_S).nii"' will work correctly. \ No newline at end of file diff --git a/GE/README.md b/GE/README.md index b8a7e1e5..cc05df1e 100644 --- a/GE/README.md +++ b/GE/README.md @@ -14,8 +14,12 @@ Knowing the relative timing of the acquisition for each 2D slice in a 3D volume In general, fMRI acquired using GE product sequence (PSD) “epi” with the multiphase option will store slice timing in the Trigger Time (DICOM 0018,1060) element. In contrast, the popular PSD “epiRT” (BrainWave RT, fMRI/DTI package provided by Medical Numerics) does not save this tag (though in some cases it saves the RTIA Timer). Examples are [available](https://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage#Slice_timing_correction) for both the “epiRT” and “epi” sequences. +The “epiRT” sequences also allow the user to specify the `Group Delay`, which is 0 msec by default. Increasing this value will create a pause at the end of each volume, and this value is recorded in the DICOM header(as 0043,107C, reported in seconds). This option can be used for sparse designs, where one wants a pause after each value. Be aware that the `RepetitionTime` (0018,0080) reported in this header omits the group delay. So a study with a TR of 2000ms and a Group Delay of 55ms will report the values (0018,0080 = 2000, 0043,107C = 0.055), while the actual sampling rate will be 2055ms. This is unintuitive, the TR with respect to tissue contrast is 2055ms, not the reported 2000. + If neither Trigger Time (DICOM 0018,1060) or RTIA Timer (0021,105E) store slice timing information, a final option is to decode the GE Protocol Data Block as described below. At best, this block only reports whether the acquisition was interleaved or sequential. As long as one assumes the acquisition was continuous (with no temporal gap between volumes, e.g. sparse images) on can use this value, the number of slices in the volume and the repetition time to infer slice times. +Due to these various methods, recent releases of dcm2niix read the Protocol Data Block to determine the multi-bandFactor, number of samples, sampling rate (TR), interleaved or sequential slices, group delay and software version. This information is used to estimate the slice timing directly, without requiring the previously described tags. The [dcm_qa_ge](https://github.com/neurolabusc/dcm_qa_ge) provides validation images and more details. Note that this slice timing approach does not support GE's diffusion weighted imaging sequences. + ## User Define Data GE (0043,102A) This private element of the DICOM header is used to determine the phase encoding polarity. Specifically, we need to know the "Ky traversal direction" (top-down, or bottom up) and the phase encoding polarity. Unfortunately, this data is stored in a complicated, proprietary structure, that has changed with different releases of GE software. [Click here to see the definition for this structure](https://github.com/ScottHaileRobertson/GE-MRI-Tools/blob/master/GePackage/%2BGE/%2BPfile/%2BHeader/%2BRDB15/rdbm.h). @@ -24,26 +28,44 @@ This private element of the DICOM header is used to determine the phase encoding Current GE software (DV26.0_R03_1831.b) running research multi-echo sequences create invalid DICOM images. The required public [EchoTime (0018,0081)](https://dicom.innolitics.com/ciods/mr-image/mr-image/00180081) attribute lists the shortest echo time for the series, rather than the actual echo time for the given DICOM image. The public tag [EchoNumber (0018,0086)](https://dicom.innolitics.com/ciods/mr-image/mr-image/00180086) reports `1` for all echoes. These limitations in GE's DICOM images disrupt dcm2niix's image conversion. Hopefully future product sequences will generate valid DICOM data. In the meantime, [issue 359](https://github.com/rordenlab/dcm2niix/issues/359) provides a kludge for image conversion. - -## Imkage Interpolation +## Image Interpolation Some sequences allow the user to interpolate images in plane (e.g. saving a 2D 64x64 EPI image as 128x128) or between slices (e.g. saving a 126 slice T1-weighted image as 252 images). The resulting files require much more disk space, add no new information, are slower to process and can [disrupt some tools](https://mrtrix.readthedocs.io/en/latest/reference/commands/mrdegibbs.html). Users are strongly discouraged from interpolating raw data. However, dcm2niix should correctly detect this interpolation, resolving apparent discrepancies between tags (0020,1002; 0021,104F; 0054,0081). [Issue 355](https://github.com/rordenlab/dcm2niix/issues/355) provides details. ## Total Readout Time -One often wants to determine [echo spacing, bandwidth](https://support.brainvoyager.com/brainvoyager/functional-analysis-preparation/29-pre-processing/78-epi-distortion-correction-echo-spacing-and-bandwidth) and total read-out time for EPI data so they can be undistorted. Total readout time is influence by parallel acceleration factor, bandwidth, number of EPI lines, and partial Fourier. Not all of these parameters are available from the GE DICOM images, so a user needs to check the scanner console. +One often wants to determine [echo spacing, bandwidth](https://support.brainvoyager.com/brainvoyager/functional-analysis-preparation/29-pre-processing/78-epi-distortion-correction-echo-spacing-and-bandwidth) and total read-out time for EPI data so they can be undistorted. Total readout time is influenced by parallel acceleration factor, bandwidth, number of EPI lines, and partial Fourier. Not all of these parameters are available from the GE DICOM images, so a user needs to check the scanner console. + +## Image Acceleration + +The BIDS [ParallelReductionFactorInPlane](https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html#in-plane-spatial-encoding) can be determined from the reciprocal of the first element of the private tag `Asset R Factors` (0043,1083). The first value is for acceleration in the phase direction equivalent to the public tag 0018,9069), the second for the slice direction (equivalent to the public tag 0018,9155). For 2D EPI scans, the second value will always be 1, whereas for 3D acquisitions [acceleration can take place in both the phase-encoding and the slice-encoding directions](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4459721/). For example, a scan using a acceleration factor of 1.5 would be reported as `(0043,1083) DS [0.666667\1]`. + +The BIDS [MultibandAccelerationFactor](https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html#rf-contrast) can be determined from the private tag `Multiband Parameters` (0043,10B6). This is an array with at least three values, the first is the Multiband (aka HyperBand) factor, the second is the Slice FOV Shift Factor, and the final is the Calibration method. For example, a scan using a multiband factor of 2 could be reported as `(0043,10b6) LO [2\4\19\\\\]`. + +## Phase-Encoding Polarity + +All EPI scans have spatial distortion, particularly those with longer readout times. Tools like [FSL topup](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/topup/TopupUsersGuide) can leverage data where two spin-echo images are acquired that are identical except for using opposite phase-encoding polarity (e.g. one uses A>P, the other P>A). Each image is distorted with the same magnitude, but in the opposite direction. GE's Rx27 software version and later populate the [Rectilinear Phase Encode Reordering (0018,9034)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9034)) tag which for EPI is set to either LINEAR or REVERSE_LINEAR. ## GE Protocol Data Block In addition to the public DICOM tags, previous versions of dcm2niix attempted to decode the proprietary GE Protocol Data Block (0025,101B). This is essentially a [GZip format](http://www.onicos.com/staff/iz/formats/gzip.html) file embedded inside the DICOM header. Unfortunately, this data seems to be [unreliable](https://github.com/rordenlab/dcm2niix/issues/163) and therefore this strategy is not used anymore. The notes below regarding the usage of this data block are provided for historical purposes. - - The VIEWORDER tag is used to set the polarity of the BIDS tag PhaseEncodingDirection, with VIEWORDER of 1 suggesting bottom up phase encoding. Unfortunately, users can separately reverse the phase encoding direction making this tag unreliable. + - The VIEWORDER tag is used to set the polarity of the BIDS tag PhaseEncodingDirection, with VIEWORDER of 1 suggesting bottom up phase encoding. Unfortunately, users can separately reverse the phase encoding direction making this tag unreliable. Thankfully, recent scanners provide 0018,9034. - The SLICEORDER tag could be used to set the SliceTiming for the BIDS tag PhaseEncodingDirection, with a SLICEORDER of 1 suggesting interleaved acquisition. - - There are reports that newer versions of GE equipement (e.g. DISCOVERY MR750 / 24\MX\MR Software release:DV24.0_R01_1344.a) are now storing an [XML](https://groups.google.com/forum/#!msg/comp.protocols.dicom/mxnCkv8A-i4/W_uc6SxLwHQJ) file within the Protocolo Data Block (compressed). In theory this might also provide useful information. + - There are reports that newer versions of GE equipement (e.g. DISCOVERY MR750 / 24\MX\MR Software release:DV24.0_R01_1344.a) are now storing an [XML](https://groups.google.com/forum/#!msg/comp.protocols.dicom/mxnCkv8A-i4/W_uc6SxLwHQJ) file within the Protocol Data Block (compressed). In theory this might also provide useful information. + +## Complex Image Component + +Most vendors store the complex image component (magnitude, phase, real or imaginary) in the [Complex Image Component (0008,9208)](http://dicom.nema.org/medical/dicom/2016e/output/chtml/part03/sect_C.8.13.3.html) tag or in the [Image Type (0008,0008)](http://incenter.medical.philips.com/doclib/enc/fetch/2000/4504/577242/577256/588723/5144873/5144488/5144982/DICOM_Conformance_Statement_Intera_R7%2c_R8_and_R9.pdf%3fnodeid%3d5147977%26vernum%3d-2) tag. GE stores this information as a signed short of the Private Image Type(0043,102F) tag. The values 0, 1, 2, 3 correspond to magnitude, phase, real, and imaginary (respectively). + +## Detecting Anatomical Localizers + +Anatomical localizers (e.g. scout images) are quick-and-dirty scans used to position subsequent slower but higher quality images. These scans are typically discarded in subsequent analyses. The dcm2niix argument `-i y` will ignore these scans. For GE, these sequences are detected based on the SeriesPlane (0019,1017) tag, which is of type [SS](http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_6.2.html) and can report the values 2 (Axial), 4 (Sagittal), 8 (Coronal), 16 (Oblique) or 256 (3plane). The 3plane value is consistent with an anatomical localizer. ## Sample Datasets - [A validation dataset for dcm2niix commits](https://github.com/neurolabusc/dcm_qa_nih). + - [Slice Timing and Phase Encoding examples](https://github.com/jannikadon/cc-dcm2bids-wrapper/tree/main/dicom-qa-examples) - [Slice timing validation](https://github.com/neurolabusc/dcm_qa_stc) for different varieties of GE EPI sequences. - - [Examples of phase encoding polarity, slice timing and diffusion gradients](https://github.com/nikadon/cc-dcm2bids-wrapper/tree/master/dicom-qa-examples/). + - [Examples of phase encoding polarity, slice timing and diffusion gradients](https://github.com/neurolabusc/dcm_qa_ge). - The dcm2niix [wiki](https://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage) includes examples of diffusion data, slice timing, and other variations. \ No newline at end of file diff --git a/Philips/README.md b/Philips/README.md index eaae15ff..6813e9b6 100644 --- a/Philips/README.md +++ b/Philips/README.md @@ -90,7 +90,7 @@ MyCustomDirections ``` -## Missing Information. +## Missing Information Philips DICOMs do not contain all the information desired by many neuroscientists. Due to this, the [BIDS](http://bids.neuroimaging.io/) files created by dcm2niix are impoverished relative to data from other vendors. This reflects a limitation in the Philips DICOMs, not dcm2niix. diff --git a/README.md b/README.md index ba006775..b9d2bb1e 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@ ## About -dcm2niix is a designed to convert neuroimaging data from the DICOM format to the NIfTI format. This web page hosts the developmental source code - a compiled version for Linux, MacOS, and Windows of the most recent stable release is included with [MRIcroGL](https://www.nitrc.org/projects/mricrogl/). A full manual for this software is available in the form of a [NITRC wiki](http://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage). +dcm2niix is designed to convert neuroimaging data from the DICOM format to the NIfTI format. This web page hosts the developmental source code - a compiled version for Linux, MacOS, and Windows of the most recent stable release is included with [MRIcroGL](https://www.nitrc.org/projects/mricrogl/). A full manual for this software is available in the form of a [NITRC wiki](http://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage). ## License @@ -25,7 +25,7 @@ DICOM provides many ways to store/compress image data, known as [transfer syntax ## Versions -[See the VERSIONS.md file for details on releases](./VERSIONS.md). +[See releases](https://github.com/rordenlab/dcm2niix/releases) for recent release notes. [See the VERSIONS.md file for details on earlier releases](./VERSIONS.md). ## Running @@ -36,12 +36,12 @@ Command line usage is described in the [NITRC wiki](https://www.nitrc.org/plugin ## Install There are a couple ways to install dcm2niix - - [Github Releases](https://github.com/rordenlab/dcm2niix/releases) provides the latest compiled executables. This is an excellent option for MacOS and Windows users. However, the provided Linux executable requires a recent version of Linux, so the provided Unix executable is not suitable for all distributions. + - [Github Releases](https://github.com/rordenlab/dcm2niix/releases) provides the latest compiled executables. This is an excellent option for MacOS and Windows users. However, the provided Linux executable requires a recent version of Linux (e.g. Ubuntu 14.04 or later), so the provided Unix executable is not suitable for very old distributions. Specifically, it requires Glibc 2.19 (from 2014) or later. Users of older systems can compile their own copy of dcm2niix or download the compiled version included with MRIcroGL Glibc 2.12 (from 2011, see below). - Run the following command to get the latest version for Linux, Macintosh or Windows: * `curl -fLO https://github.com/rordenlab/dcm2niix/releases/latest/download/dcm2niix_lnx.zip` * `curl -fLO https://github.com/rordenlab/dcm2niix/releases/latest/download/dcm2niix_mac.zip` * `curl -fLO https://github.com/rordenlab/dcm2niix/releases/latest/download/dcm2niix_win.zip` - - [MRIcroGL (NITRC)](https://www.nitrc.org/projects/mricrogl) or [MRIcroGL (GitHub)](https://github.com/rordenlab/MRIcroGL12/releases) includes dcm2niix that can be run from the command line or from the graphical user interface (select the Import menu item). The Linux version of dcm2niix is compiled on a holy build box, so it should run on any Linux distribution. + - [MRIcroGL (NITRC)](https://www.nitrc.org/projects/mricrogl) or [MRIcroGL (GitHub)](https://github.com/rordenlab/MRIcroGL12/releases) includes dcm2niix that can be run from the command line or from the graphical user interface (select the Import menu item). The Linux version of dcm2niix is compiled on a [holy build box](https://github.com/phusion/holy-build-box), so it should run on any Linux distribution. - If you have a MacOS computer with Homebrew you can run `brew install dcm2niix`. - If you have Conda, [`conda install -c conda-forge dcm2niix`](https://anaconda.org/conda-forge/dcm2niix) on Linux, MacOS or Windows. - On Debian Linux computers you can run `sudo apt-get install dcm2niix`. @@ -98,6 +98,7 @@ If you have any problems with the cmake build script described above or want to ## Alternatives - [dcm2nii](https://people.cas.sc.edu/rorden/mricron/dcm2nii.html) is the predecessor of dcm2niix. It is deprecated for modern images, but does handle image formats that predate DICOM (proprietary Elscint, GE and Siemens formats). + - Python [dcmstack](https://github.com/moloney/dcmstack) DICOM to Nifti conversion with meta data preservation. - [dicm2nii](http://www.mathworks.com/matlabcentral/fileexchange/42997-dicom-to-nifti-converter) is written in Matlab. The Matlab language makes this very scriptable. - [dicom2nifti](https://github.com/icometrix/dicom2nifti) uses the scriptable Python wrapper utilizes the [high performance GDCMCONV](http://gdcm.sourceforge.net/wiki/index.php/Gdcmconv) executables. - [dicomtonifti](https://github.com/dgobbi/vtk-dicom/wiki/dicomtonifti) leverages [VTK](https://www.vtk.org/). @@ -106,9 +107,12 @@ If you have any problems with the cmake build script described above or want to - [mcverter](http://lcni.uoregon.edu/%7Ejolinda/MRIConvert/) has great support for various vendors. - [mri_convert](https://surfer.nmr.mgh.harvard.edu/pub/docs/html/mri_convert.help.xml.html) is part of the popular FreeSurfer package. In my limited experience this tool works well for GE and Siemens data, but fails with Philips 4D datasets. - [MRtrix mrconvert](http://mrtrix.readthedocs.io/en/latest/reference/commands/mrconvert.html) is a useful general purpose image converter and handles DTI data well. It is an outstanding tool for modern Philips enhanced images. + - [PET CT viewer](http://petctviewer.org/index.php/feature/results-exports/nifti-export) for [Fiji](https://fiji.sc) can load DICOM images and export as NIfTI. - [Plastimatch](https://www.plastimatch.org/) is a Swiss Army knife - it computes registration, image processing, statistics and it has a basic image format converter that can convert some DICOM images to NIfTI or NRRD. + - [Simple Dicom Reader 2 (Sdr2)](http://ogles.sourceforge.net/sdr2-doc/index.html) uses [dcmtk](https://dicom.offis.de/dcmtk.php.en) to read DICOM images and convert them to the NIfTI format. + - [SlicerHeart extension](https://github.com/SlicerHeart/SlicerHeart) is specifically designed to help 3D Slicer support ultra sound (US) images stored as DICOM. + - [spec2nii](https://github.com/wexeee/spec2nii) converts MR spectroscopy to NIFTI. - [SPM12](http://www.fil.ion.ucl.ac.uk/spm/software/spm12/) is one of the most popular tools in the field. It includes DICOM to NIfTI conversion. Being based on Matlab it is easy to script. - - dcm2niix is largely used for converting MRI and CT DICOMs. While it should in theory support ultra sound (US) images stored as DICOMs, vendors tend to use private tags. The [SlicerHeart extension](https://github.com/SlicerHeart/SlicerHeart) is specifically designed to help 3D Slicer support this modality. ## Links @@ -117,10 +121,13 @@ The following tools exploit dcm2niix - [abcd-dicom2bids](https://github.com/DCAN-Labs/abcd-dicom2bids) selectively downloads high quality ABCD datasets. - [autobids](https://github.com/khanlab/autobids) automates dcm2bids which uses dcm2niix. - [BiDirect_BIDS_Converter](https://github.com/wulms/BiDirect_BIDS_Converter) for conversion from DICOM to the BIDS standard. + - [BIDScoin](https://github.com/Donders-Institute/bidscoin) is a DICOM to BIDS converter with thorough [documentation](https://bircibrain.github.io/computingguide/docs/bids/bidscoin). - [BIDS Toolbox](https://github.com/cardiff-brain-research-imaging-centre/bids-toolbox) is a web service for the creation and manipulation of BIDS datasets, using dcm2niix for importing DICOM data. + - [birc-bids](https://github.com/bircibrain/birc-bids) provides a Docker/Singularity container with various BIDS conversion utilities. + - [BOLD5000_autoencoder](https://github.com/nmningmei/BOLD5000_autoencoder) uses dcm2niix to pipe imaging data into an unsupervised machine learning algorithm. - [brainnetome DiffusionKit](http://diffusion.brainnetome.org/en/latest/) uses dcm2niix to convert images. - [Brain imAgiNg Analysis iN Arcana (Banana)](https://pypi.org/project/banana/) is a collection of brain imaging analysis workflows, it uses dcm2niix for format conversions. - - [BOLD5000_autoencoder](https://github.com/nmningmei/BOLD5000_autoencoder) uses dcm2niix to pipe imaging data into an unsupervised machine learning algorithm. + - [BraTS-Preprocessor](https://neuronflow.github.io/BraTS-Preprocessor/) uses dcm2niix to import files for [Brain Tumor Segmentation](https://www.frontiersin.org/articles/10.3389/fnins.2020.00125/full). - [clinica](https://github.com/aramis-lab/clinica) is a software platform for clinical neuroimaging studies that uses dcm2niix to convert DICOM images. - [dcm2niix can help convert data from the Adolescent Brain Cognitive Development (ABCD) DICOM to BIDS](https://github.com/ABCD-STUDY/abcd-dicom2bids) - [bidsify](https://github.com/spinoza-rec/bidsify) is a Python project that uses dcm2niix to convert DICOM and Philips PAR/REC images to the BIDS standard. @@ -128,20 +135,23 @@ The following tools exploit dcm2niix - [BioImage Suite Web Project](https://github.com/bioimagesuiteweb/bisweb) is a JavaScript project that uses dcm2niix for its DICOM conversion module. - [boutiques-dcm2niix](https://github.com/lalet/boutiques-dcm2niix) is a dockerfile for installing and validating dcm2niix. - [DAC2BIDS](https://github.com/dangom/dac2bids) uses dcm2niibatch to create [BIDS](http://bids.neuroimaging.io/) datasets. - - [Dcm2Bids](https://github.com/cbedetti/Dcm2Bids) uses dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets. + - [Dcm2Bids](https://github.com/cbedetti/Dcm2Bids) uses dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets. Here is a [tutorial](https://andysbrainbook.readthedocs.io/en/latest/OpenScience/OS/BIDS_Overview.html) describing usage. - [dcm2niir](https://github.com/muschellij2/dcm2niir) R wrapper for dcm2niix/dcm2nii. - [dcm2niix_afni](https://afni.nimh.nih.gov/pub/dist/doc/program_help/dcm2niix_afni.html) is a version of dcm2niix included with the [AFNI](https://afni.nimh.nih.gov/) distribution. - [dcm2niiXL](https://github.com/neurolabusc/dcm2niiXL) is a shell script and tuned compilation of dcm2niix designed for accelerated conversion of extra large datasets. - [DICOM2BIDS](https://github.com/klsea/DICOM2BIDS) is a Python 2 script for creating BIDS files. + - [dicom2bids](https://github.com/Jolinda/lcnimodules) includes python modules for converting dicom files to nifti in a bids-compatible file structure that use dcm2niix. - [dcmwrangle](https://github.com/jbteves/dcmwrangle) a Python interactive and static tool for organizing dicoms. - [dicom2nifti_batch](https://github.com/scanUCLA/dicom2nifti_batch) is a Matlab script for automating dcm2niix. - [divest](https://github.com/jonclayden/divest) R interface to dcm2niix. + - [ezBIDS](https://github.com/brainlife/ezbids) is a web service for converting directory full of DICOM images into BIDS without users having to learn python nor custom configuration file. - [fmrif tools](https://github.com/nih-fmrif/fmrif_tools) uses dcm2niix for its [oxy2bids](https://fmrif-tools.readthedocs.io/en/latest/#) tool. - [fsleyes](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FSLeyes) is a powerful Python-based image viewer. It uses dcm2niix to handle DICOM files through its fslpy libraries. - [Functional Real-Time Interactive Endogenous Neuromodulation and Decoding (FRIEND) Engine](https://github.com/InstitutoDOr/FriendENGINE) uses dcm2niix. - - [heudiconv](https://github.com/nipy/heudiconv) can use dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets. + - [heudiconv](https://github.com/nipy/heudiconv) can use dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets. Data acquired using the [reproin](https://github.com/ReproNim/reproin) convention can be easily converted to BIDS. - [kipettools](https://github.com/mathesong/kipettools) uses dcm2niix to load PET data. - [LEAD-DBS](http://www.lead-dbs.org/) uses dcm2niix for [DICOM import](https://github.com/leaddbs/leaddbs/blob/master/ea_dicom_import.m). + - [lin4neuro](http://www.lin4neuro.net/lin4neuro/18.04bionic/vm/) releases such as the English l4n-18.04.4-amd64-20200801-en.ova include MRIcroGL and dcm2niix pre-installed. This allows user with VirtualBox or VMWarePlayer to use these tools (and many other neuroimaging tools) in a graphical virtual machine. - [MRIcroGL](https://github.com/neurolabusc/MRIcroGL) is available for MacOS, Linux and Windows and provides a graphical interface for dcm2niix. You can get compiled copies from the [MRIcroGL NITRC web site](https://www.nitrc.org/projects/mricrogl/). - [neurodocker](https://github.com/kaczmarj/neurodocker) includes dcm2niix as a lean, minimal install Dockerfile. - [neuro_docker](https://github.com/Neurita/neuro_docker) includes dcm2niix as part of a single, static Dockerfile. diff --git a/Siemens/README.md b/Siemens/README.md index 8f5918d9..51975e7c 100644 --- a/Siemens/README.md +++ b/Siemens/README.md @@ -2,9 +2,15 @@ dcm2niix attempts to convert Siemens DICOM format images to NIfTI. This page describes some vendor-specific details. -## Vida XA10 and XA11 +## Siemens X-Series -The Siemens Vida introduced the new XA10 DICOM format. Users are strongly encouraged to export data using the "Enhanced" format and to not use any of the "Anonymize" features on the console. The consequences of these options is discussed in detail in [issue 236](https://github.com/rordenlab/dcm2niix/issues/236). In brief, the Vida can export to enhanced, mosaic or classic 2D. Note that the mosaics are considered secondary capture images intended for quality assurance only. The mosaic scans lack several "Type 1" DICOM properties, necessarily limiting conversion. The non-mosaic 2D enhanced DICOMs are compact and efficient, but appear to have limited details relative to the enhanced output. Finally, each of the formats (enhanced, mosaic, classic) can be exported as anonymized. The Siemens console anonymization of current XA10A (Fall 2018) strips many useful tags. Siemens suggests "the use an offline/in-house anonymization software instead." Another limitation of the current XA10A format is that it retains no versioning details for software and hardware stepping, despite the fact that the data format is rapidly evolving. If you use a Vida, you are strongly encouraged to log every hardware or software upgrade to allow future analyses to identify and regress out any effects of these modifications. Since the XA10A format does not have a CSA header, dcm2niix will attempt to use the new private DICOM tags to populate the BIDS file. These tags are described in [issue 240](https://github.com/rordenlab/dcm2niix/issues/240). +Siemens MR is named by Series, Generation, Major Version and Minor Version. Prior to the Siemens Vida, all contemporary Siemens MRI systems (Trio, Prisma, Skyra, etc) were part of the V series. So a Trio might be on VB17, and a Prisma on VE11 (series 'V', generation 'E', major version '1', minor version '1'). The 3T Vida and 1.5T Sola introduce the X-series (XA10, XA11, XA20). Since the V-series was dominant for so long, most users simply omit the series, e.g. referring to a system as `B19`. However, Siemens has recently introduced a new X-series. + +The DICOM images exported by the X-series is radically different than the V-series. The images lack the proprietary CSA header with its rich meta data. + +X-series users are strongly encouraged to export data using the "Enhanced" format and to not use any of the "Anonymize" features on the console. The consequences of these options is discussed in detail in [issue 236](https://github.com/rordenlab/dcm2niix/issues/236). Siemens notes `We highly recommend that the Enhanced DICOM format be used. This is because this format retains far more information in the header`. Failure to export data in this format has led to catastrophic data loss for numerous users (for publicly reported details see issues [203](https://github.com/rordenlab/dcm2niix/issues/203), [236](https://github.com/rordenlab/dcm2niix/issues/236), [240](https://github.com/rordenlab/dcm2niix/issues/240), [274](https://github.com/rordenlab/dcm2niix/issues/274), [303](https://github.com/rordenlab/dcm2niix/issues/303), [370](https://github.com/rordenlab/dcm2niix/issues/370), [394](https://github.com/rordenlab/dcm2niix/issues/394)). This reflects limitations of the DICOM data, not dcm2niix. + +While X-series consoles allow users to export data as enhanced, mosaic or classic 2D formats, choosing an option other than enhanced dramatically degrades the meta data. Note that the Siemens considers mosaic images `secondary capture` data intended for quality assurance only. The mosaic scans lack several "Type 1" DICOM properties, necessarily limiting conversion. The non-mosaic 2D enhanced DICOMs are compact and efficient, but appear to have limited details relative to the enhanced output. This is unfortunate, as for the V-series the mosaic format has major benefits, so users may be in the habit of prefering mosaic export. Finally, each of the formats (enhanced, mosaic, classic) can be exported as anonymized. The Siemens console anonymization of current XA10A (Fall 2018) strips many useful tags. Siemens suggests `the use an offline/in-house anonymization software instead`. Another limitation of the current X-series format is that it retains no versioning details beyond the minor version for software and hardware stepping (e.g. versions are merely XA10 or XA11 with no details for service packs). If you use a X-series, you are strongly encouraged to manually log every hardware or software upgrade to allow future analyses to identify and regress out any effects of these modifications. Since the X-series format does not have a CSA header, dcm2niix will attempt to use the new private DICOM tags to populate the BIDS file. These tags are described in [issue 240](https://github.com/rordenlab/dcm2niix/issues/240). When creating enhanced DICOMs diffusion information is provided in public tags. Based on a limited sample, it seems that classic DICOMs do not store diffusion data for XA10, and use private tags for [XA11](https://www.nitrc.org/forum/forum.php?thread_id=10013&forum_id=4703). @@ -26,11 +32,10 @@ In theory, the public DICOM tag 'Frame Acquisition Date Time' (0018,9074) and th ## CSA Header -Many crucial Siemens parameters are stored in the [proprietary CSA header](http://nipy.org/nibabel/dicom/siemens_csa.html). This has a binary section that allows quick reading for many useful parameters. It also includes an ASCII text portion that includes a lot of information but is slow to parse and poorly curated. Be aware that Siemens Vida scanners do not generate a CSA header. +Many crucial Siemens parameters are stored in the [proprietary CSA header](http://nipy.org/nibabel/dicom/siemens_csa.html), in particular the CSA Image Header Info (0029, 1010) and CSA Series Header Info (0029, 1020). These have binary sections that allows quick reading for many useful parameters. They also include an ASCII text portion that includes a lot of information but is slow to parse and poorly curated. Be aware that Siemens Vida scanners do not generate a CSA header. ## Slice Timing - See the [dcm_qa_stc](https://github.com/neurolabusc/dcm_qa_stc) repository with sample data that exhibits different methods used by Siemens to record slice timing. Older software (e.g. A25 through B13) sometimes populates the tag sSliceArray.ucMode in the [CSA Series Header (0029, 1020)](https://nipy.org/nibabel/dicom/siemens_csa.html) where the values [1, 2, and 4](https://github.com/xiangruili/dicm2nii/issues/18) correspond to Ascending, Descending and Interleaved acquisitions. @@ -45,11 +50,13 @@ One often wants to determine [echo spacing, bandwidth, ](https://support.brainvo ## Diffusion Tensor Notes -Diffusion specific parameters are described on the [NA-MIC](https://www.na-mic.org/wiki/NAMIC_Wiki:DTI:DICOM_for_DWI_and_DTI#Private_vendor:_Siemens) website. Gradient vectors are reported with respect to the scanner bore, and dcm2niix will attempt to re-orient these to [FSL format](http://justinblaber.org/brief-introduction-to-dwmri/) [bvec files](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FDT/FAQ#What_conventions_do_the_bvecs_use.3F). For the Vida, see the Vida section for specific diffusion details. +Diffusion specific parameters are described on the [NA-MIC](https://www.na-mic.org/wiki/NAMIC_Wiki:DTI:DICOM_for_DWI_and_DTI#Private_vendor:_Siemens) website. Gradient vectors are reported with respect to the scanner bore, and dcm2niix will attempt to re-orient these to [FSL format](http://justinblaber.org/brief-introduction-to-dwmri/) [bvec files](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FDT/FAQ#What_conventions_do_the_bvecs_use.3F). + +For Siemens V-series systems from the B-generation onward (around 2005), the most reliable way to read diffusion gradients is from the [CSA header](https://nipy.org/nibabel/dicom/siemens_csa.html). Specially, the CSA's 'DiffusionGradientDirection' and 'B_value' tags. For the X-series, the private DICOM tags B_value (0019,100c) and DiffusionGradientDirection (0019,100e) are used. ## Arterial Spin Labeling -Tools like [ExploreASL](https://sites.google.com/view/exploreasl) and [FSL BASIL](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/BASIL) can help process arterial spin labeling data. These tools require sequence details. These details differ between different sequences. If you create a BIDS JSON file with dcm2niix, the following tags will be created, using the same names used in the Siemens sequence PDFs. Note different sequences provide different values. +Tools like [ExploreASL](https://sites.google.com/view/exploreasl) and [FSL BASIL](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/BASIL) can help process arterial spin labeling data. These tools require sequence details. These details differ between different sequences. If you create a BIDS JSON file with dcm2niix, the following tags will be created, using the same names used in the Siemens sequence PDFs. Note different sequences provide different values. The [dcm_qa_asl](https://github.com/neurolabusc/dcm_qa_asl) repository provides example DICOM ASL datasets. ep2d_pcasl, ep2d_pcasl_UI_PHC //pCASL 2D [Danny J.J. Wang](http://www.loft-lab.org) - LabelOffset diff --git a/console/main_console.cpp b/console/main_console.cpp index 25452ffc..1c135f2d 100644 --- a/console/main_console.cpp +++ b/console/main_console.cpp @@ -80,7 +80,7 @@ void showHelp(const char * argv[], struct TDCMopts opts) { printf(" -a : adjacent DICOMs (images from same series always in same folder) for faster conversion (n/y, default n)\n"); printf(" -b : BIDS sidecar (y/n/o [o=only: no NIfTI], default %c)\n", bool2Char(opts.isCreateBIDS)); printf(" -ba : anonymize BIDS (y/n, default %c)\n", bool2Char(opts.isAnonymizeBIDS)); - printf(" -c : comment stored in NIfTI aux_file (provide up to 24 characters)\n"); + printf(" -c : comment stored in NIfTI aux_file (provide up to 24 characters e.g. '-c first_visit')\n"); printf(" -d : directory search depth. Convert DICOMs in sub-folders of in_folder? (0..9, default %d)\n", opts.dirSearchDepth); printf(" -e : export as NRRD instead of NIfTI (y/n, default n)\n"); #ifdef mySegmentByAcq @@ -374,6 +374,13 @@ int main(int argc, const char * argv[]) opts.isIgnoreDerivedAnd2D = false; else opts.isIgnoreDerivedAnd2D = true; + } else if ((argv[i][1] == 'j') && ((i+1) < argc)) { + i++; + if (invalidParam(i, argv)) return 0; + if ((argv[i][0] == 'y') || (argv[i][0] == 'Y') || (argv[i][0] == '1')) { + opts.isTestx0021x105E = true; + printf("undocumented '-j y' compares GE slice timing from 0021,105E\n"); + } } else if ((argv[i][1] == 'l') && ((i+1) < argc)) { i++; if (invalidParam(i, argv)) return 0; diff --git a/console/nifti1.h b/console/nifti1.h index 7176c968..680d1637 100644 --- a/console/nifti1.h +++ b/console/nifti1.h @@ -1498,4 +1498,4 @@ typedef struct { unsigned char r,g,b; } rgb_byte ; #endif /*=================*/ -#endif /* _NIFTI_HEADER_ */ \ No newline at end of file +#endif /* _NIFTI_HEADER_ */ diff --git a/console/nifti1_io_core.cpp b/console/nifti1_io_core.cpp index 4320dd20..758e2c01 100644 --- a/console/nifti1_io_core.cpp +++ b/console/nifti1_io_core.cpp @@ -17,7 +17,6 @@ #include //requires VS 2015 or later #include "nifti1_io_core.h" -#include "nifti1.h" #include #include #include @@ -80,7 +79,6 @@ void nifti_swap_2bytes( size_t n , void *ar ) // 2 bytes at a time } return ; } -#endif /*-------------------------------------------------------------------------*/ /*! Byte swap NIFTI-1 file header in various places and ways. @@ -137,6 +135,7 @@ void swap_nifti_header( struct nifti_1_header *h ) return ; } +#endif bool littleEndianPlatform () { diff --git a/console/nifti1_io_core.h b/console/nifti1_io_core.h index ff371ad5..1e15cbe3 100644 --- a/console/nifti1_io_core.h +++ b/console/nifti1_io_core.h @@ -6,10 +6,12 @@ #ifdef USING_R #define STRICT_R_HEADERS +#include #include "RNifti.h" -#endif +#else #include "nifti1.h" #include +#endif #ifdef __cplusplus extern "C" { @@ -70,7 +72,10 @@ vec3 nifti_vect33mat33_mul(vec3 v, mat33 m ); ivec3 setiVec3(int x, int y, int z); vec3 setVec3(float x, float y, float z); vec4 setVec4(float x, float y, float z); +#ifndef USING_R +// This declaration differs from the equivalent function in the current nifti1_io.h, so avoid the clash void swap_nifti_header ( struct nifti_1_header *h) ; +#endif vec4 nifti_vect44mat44_mul(vec4 v, mat44 m ); void nifti_swap_2bytes( size_t n , void *ar ); // 2 bytes at a time void nifti_swap_4bytes( size_t n , void *ar ); // 4 bytes at a time @@ -87,4 +92,4 @@ mat44 nifti_quatern_to_mat44( float qb, float qc, float qd, } #endif -#endif \ No newline at end of file +#endif diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index bed326e9..9a6d77a4 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -738,6 +738,7 @@ struct TDICOMdata clear_dicom_data() { strcpy(d.bodyPartExamined,""); strcpy(d.coilName, ""); strcpy(d.coilElements, ""); + strcpy(d.radiopharmaceutical, ""); d.phaseEncodingLines = 0; //~ d.patientPositionSequentialRepeats = 0; //~ d.patientPositionRepeats = 0; @@ -770,6 +771,8 @@ struct TDICOMdata clear_dicom_data() { d.echoNum = 1; d.echoTrainLength = 0; d.waterFatShift = 0.0; + d.groupDelay = 0.0; + d.decayFactor = 0.0; d.percentSampling = 0.0; d.phaseFieldofView = 0.0; d.dwellTime = 0; @@ -780,6 +783,7 @@ struct TDICOMdata clear_dicom_data() { d.seriesUidCrc = 0; d.instanceUidCrc = 0; d.accelFactPE = 0.0; + d.accelFactOOP = 0.0; //d.patientPositionNumPhilips = 0; d.imageBytes = 0; d.intenScale = 1; @@ -787,11 +791,12 @@ struct TDICOMdata clear_dicom_data() { d.intenIntercept = 0; d.gantryTilt = 0.0; d.radionuclidePositronFraction = 0.0; - d.radionuclideTotalDose = 0.0; d.radionuclideHalfLife = 0.0; d.doseCalibrationFactor = 0.0; d.ecat_isotope_halflife = 0.0; + d.frameDuration = -1.0; d.ecat_dosage = 0.0; + d.radionuclideTotalDose = 0.0; d.seriesNum = 1; d.acquNum = 0; d.imageNum = 1; @@ -803,6 +808,8 @@ struct TDICOMdata clear_dicom_data() { d.isGrayscaleSoftcopyPresentationState = false; d.isRawDataStorage = false; d.isPartialFourier = false; + d.isIR = false; + d.isEPI = false; d.isDiffusion = false; d.isVectorFromBMatrix = false; d.isStackableSeries = false; //combine DCE series https://github.com/rordenlab/dcm2niix/issues/252 @@ -834,6 +841,13 @@ struct TDICOMdata clear_dicom_data() { d.rtia_timerGE = -1.0; d.rawDataRunNumber = -1; d.maxEchoNumGE = -1; + d.epiVersionGE = -1; + d.internalepiVersionGE = -1; + d.interp3D = -1; + for (int i = 0; i < kMaxOverlay; i++) + d.overlayStart[i] = 0; + d.isHasOverlay = false; + d.isPrivateCreatorRemap = false; d.numberOfImagesInGridUIH = 0; d.phaseEncodingRC = '?'; d.patientSex = '?'; @@ -912,19 +926,14 @@ uint32_t mz_crc32X(unsigned char *ptr, size_t buf_len) return ~crcu32; } - - -void dcmStr(int lLength, unsigned char lBuffer[], char* lOut, bool isStrLarge = false, bool isSpaceToUnderscore = true) { +void dcmStr(int lLength, unsigned char lBuffer[], char* lOut, bool isStrLarge = false) { if (lLength < 1) return; -//#ifdef _MSC_VER char * cString = (char *)malloc(sizeof(char) * (lLength + 1)); -//#else -// char cString[lLength + 1]; -//#endif cString[lLength] =0; memcpy(cString, (char*)&lBuffer[0], lLength); //memcpy(cString, test, lLength); //printMessage("X%dX\n", (unsigned char)d.patientName[1]); + #ifdef ISO8859 for (int i = 0; i < lLength; i++) //assume specificCharacterSet (0008,0005) is ISO_IR 100 http://en.wikipedia.org/wiki/ISO/IEC_8859-1 if (cString[i]< 1) { @@ -952,22 +961,10 @@ void dcmStr(int lLength, unsigned char lBuffer[], char* lOut, bool isStrLarge = if (c == 253) cString[i] = 'y'; if (c == 255) cString[i] = 'y'; } - for (int i = 0; i < lLength; i++) - if ((cString[i]<1) || (cString[i]==',') || (cString[i]=='^') || (cString[i]=='/') || (cString[i]=='\\') || (cString[i]=='%') || (cString[i]=='*') || (cString[i] == 9) || (cString[i] == 10) || (cString[i] == 11) || (cString[i] == 13)) cString[i] = '_'; - if (!isSpaceToUnderscore) { - if (cString[lLength-1]==' ') cString[lLength-1] = '_'; //only remove trailing space, e.g. Philips Image Type with odd length - } else { - for (int i = 0; i < lLength; i++) - if (cString[i]==' ') cString[i] = '_'; - } - int len = 1; - for (int i = 1; i < lLength; i++) { //remove repeated "_" - if ((cString[i-1]!='_') || (cString[i]!='_')) { - cString[len] =cString[i]; - len++; - } - } //for each item - if (cString[len-1] == '_') len--; + #endif + //we no longer sanitize strings, see issue 425 + int len = lLength; + if (cString[len-1] == ' ') len--; //while ((len > 0) && (cString[len]=='_')) len--; //remove trailing '_' cString[len] = 0; //null-terminate, strlcpy does this anyway int maxLen = kDICOMStr; @@ -978,9 +975,7 @@ void dcmStr(int lLength, unsigned char lBuffer[], char* lOut, bool isStrLarge = } memcpy(lOut,cString,len-1); lOut[len-1] = 0; -//#ifdef _MSC_VER free(cString); -//#endif } //dcmStr() #ifdef MY_OLD //this code works on Intel but not some older systems https://github.com/rordenlab/dcm2niix/issues/327 @@ -1523,7 +1518,7 @@ int headerDcm2Nii(struct TDICOMdata d, struct nifti_1_header *h, bool isComputeS h->pixdim[1] = d.xyzMM[1]; h->pixdim[2] = d.xyzMM[2]; h->pixdim[3] = d.xyzMM[3]; - h->pixdim[4] = d.TR/1000; //TR reported in msec, time is in sec + h->pixdim[4] = d.TR / 1000.0; //TR reported in msec, time is in sec h->dim[1] = d.xyzDim[1]; h->dim[2] = d.xyzDim[2]; h->dim[3] = d.xyzDim[3]; @@ -1624,7 +1619,8 @@ void cleanStr(char* lOut) { if (c == 255) cString[i] = 'y'; } for (int i = 0; i < lLength; i++) - if ((cString[i]<1) || (cString[i]==' ') || (cString[i]==',') || (cString[i]=='^') || (cString[i]=='/') || (cString[i]=='\\') || (cString[i]=='%') || (cString[i]=='*') || (cString[i] == 9) || (cString[i] == 10) || (cString[i] == 11) || (cString[i] == 13)) cString[i] = '_'; + if ((cString[i]<1) || (cString[i]==' ') || (cString[i]==',') || (cString[i]=='/') || (cString[i]=='\\') || (cString[i]=='%') || (cString[i]=='*') || (cString[i] == 9) || (cString[i] == 10) || (cString[i] == 11) || (cString[i] == 13)) cString[i] = '_'; //issue398 + //if ((cString[i]<1) || (cString[i]==' ') || (cString[i]==',') || (cString[i]=='^') || (cString[i]=='/') || (cString[i]=='\\') || (cString[i]=='%') || (cString[i]=='*') || (cString[i] == 9) || (cString[i] == 10) || (cString[i] == 11) || (cString[i] == 13)) cString[i] = '_'; int len = 1; for (int i = 1; i < lLength; i++) { //remove repeated "_" if ((cString[i-1]!='_') || (cString[i]!='_')) { @@ -1653,7 +1649,12 @@ int isSameFloatGE (float a, float b) { struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D *dti4D, bool isReadPhase) { struct TDICOMdata d = clear_dicom_data(); dti4D->sliceOrder[0] = -1; +dti4D->volumeOnsetTime[0] = -1; +dti4D->decayFactor[0] = -1; +dti4D->frameDuration[0] = -1; dti4D->intenScale[0] = 0.0; +dti4D->repetitionTimeExcitation = 0.0; +dti4D->repetitionTimeInversion = 0.0; strcpy(d.protocolName, ""); //erase dummy with empty strcpy(d.seriesDescription, ""); //erase dummy with empty strcpy(d.sequenceName, ""); //erase dummy with empty @@ -2264,7 +2265,6 @@ int kbval = 33; //V3: 27 num2DExpected = d.xyzDim[3] * num3DExpected; } float TRms = 1000.0f * (maxDynTime - minDynTime) / (float)(numDyn-1); //-1 for fence post - //float TRms = 1000.0f * (maxDynTime - minDynTime) / (float)(d.CSA.numDti-1); if (fabs(TRms - d.TR) > 0.005f) printWarning("Reported TR=%gms, measured TR=%gms (prospect. motion corr.?)\n", d.TR, TRms); d.TR = TRms; @@ -4060,7 +4060,12 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru strcpy(d.sequenceName, ""); //erase dummy with empty //do not read folders - code specific to GCC (LLVM/Clang seems to recognize a small file size) dti4D->sliceOrder[0] = -1; + dti4D->volumeOnsetTime[0] = -1; + dti4D->decayFactor[0] = -1; + dti4D->frameDuration[0] = -1; dti4D->intenScale[0] = 0.0; + dti4D->repetitionTimeExcitation = 0.0; + dti4D->repetitionTimeInversion = 0.0; struct TVolumeDiffusion volDiffusion = initTVolumeDiffusion(&d, dti4D); struct stat s; if( stat(fname,&s) == 0 ) { @@ -4141,6 +4146,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kAcquisitionDate 0x0008+(0x0022 << 16 ) #define kAcquisitionDateTime 0x0008+(0x002A << 16 ) #define kStudyTime 0x0008+(0x0030 << 16 ) +#define kSeriesTime 0x0008+(0x0031 << 16 ) #define kAcquisitionTime 0x0008+(0x0032 << 16 ) //TM //#define kContentTime 0x0008+(0x0033 << 16 ) //TM #define kModality 0x0008+(0x0060 << 16 ) //CS @@ -4171,15 +4177,14 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kScanOptions 0x0018+(0x0022 << 16) #define kMRAcquisitionType 0x0018+(0x0023 << 16) #define kSequenceName 0x0018+(0x0024 << 16) +#define kRadiopharmaceutical 0x0018+(0x0031 << 16 ) //LO #define kZThick 0x0018+(0x0050 << 16 ) #define kTR 0x0018+(0x0080 << 16 ) #define kTE 0x0018+(0x0081 << 16 ) +#define kTI 0x0018+(0x0082 << 16) // Inversion time #define kNumberOfAverages 0x0018+(0x0083 << 16 ) //DS #define kImagingFrequency 0x0018+(0x0084 << 16 ) //DS -#define kTriggerTime 0x0018+(0x1060 << 16 ) //DS //#define kEffectiveTE 0x0018+(0x9082 << 16 ) -const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16); -#define kTI 0x0018+(0x0082 << 16) // Inversion time #define kEchoNum 0x0018+(0x0086 << 16 ) //IS #define kMagneticFieldStrength 0x0018+(0x0087 << 16 ) //DS #define kZSpacing 0x0018+(0x0088 << 16 ) //'DS' 'SpacingBetweenSlices' @@ -4191,29 +4196,37 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16); #define kDeviceSerialNumber 0x0018+(0x1000 << 16 ) //LO #define kSoftwareVersions 0x0018+(0x1020 << 16 ) //LO #define kProtocolName 0x0018+(0x1030<< 16 ) +#define kTriggerTime 0x0018+(0x1060 << 16 ) //DS #define kRadionuclideTotalDose 0x0018+(0x1074<< 16 ) #define kRadionuclideHalfLife 0x0018+(0x1075<< 16 ) #define kRadionuclidePositronFraction 0x0018+(0x1076<< 16 ) #define kGantryTilt 0x0018+(0x1120 << 16 ) #define kXRayExposure 0x0018+(0x1152 << 16 ) +#define kConvolutionKernel 0x0018+(0x1210 << 16 ) //SH +#define kFrameDuration 0x0018+(0x1242 << 16 ) //IS #define kReceiveCoilName 0x0018+(0x1250 << 16 ) // SH #define kAcquisitionMatrix 0x0018+(0x1310 << 16 ) //US #define kFlipAngle 0x0018+(0x1314 << 16 ) #define kInPlanePhaseEncodingDirection 0x0018+(0x1312<< 16 ) //CS #define kSAR 0x0018+(0x1316 << 16 ) //'DS' 'SAR' #define kPatientOrient 0x0018+(0x5100<< 16 ) //0018,5100. patient orientation - 'HFS' +#define kInversionRecovery 0x0018+uint32_t(0x9009 << 16 ) //'CS' 'YES'/'NO' +#define kEchoPlanarPulseSequence 0x0018+uint32_t(0x9018 << 16 ) //'CS' 'YES'/'NO' +#define kRectilinearPhaseEncodeReordering 0x0018+uint32_t(0x9034 << 16) //'CS' 'REVERSE_LINEAR'/'LINEAR' #define kParallelReductionFactorInPlane 0x0018+uint32_t(0x9069<< 16 ) //FD #define kAcquisitionDuration 0x0018+uint32_t(0x9073<< 16 ) //FD //#define kFrameAcquisitionDateTime 0x0018+uint32_t(0x9074<< 16 ) //DT "20181019212528.232500" #define kDiffusionDirectionality 0x0018+uint32_t(0x9075<< 16 ) // NONE, ISOTROPIC, or DIRECTIONAL #define kInversionTimes 0x0018+uint32_t(0x9079<< 16 ) //FD #define kPartialFourier 0x0018+uint32_t(0x9081<< 16 ) //CS +const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16); //#define kDiffusionBFactorSiemens 0x0019+(0x100C<< 16 ) // 0019;000C;SIEMENS MR HEADER ;B_value #define kDiffusion_bValue 0x0018+uint32_t(0x9087<< 16 ) // FD #define kDiffusionOrientation 0x0018+uint32_t(0x9089<< 16 ) // FD, seen in enhanced // DICOM from Philips 5.* // and Siemens XA10. #define kImagingFrequency2 0x0018+uint32_t(0x9098 << 16 ) //FD +#define kParallelReductionFactorOutOfPlane 0x0018+uint32_t(0x9155<< 16 ) //FD //#define kFrameAcquisitionDuration 0x0018+uint32_t(0x9220 << 16 ) //FD #define kDiffusionBValueXX 0x0018+uint32_t(0x9602 << 16 ) //FD #define kDiffusionBValueXY 0x0018+uint32_t(0x9603 << 16 ) //FD @@ -4224,10 +4237,11 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16); #define kMREchoSequence 0x0018+uint32_t(0x9114<< 16 ) //SQ #define kMRAcquisitionPhaseEncodingStepsInPlane 0x0018+uint32_t(0x9231<< 16 ) //US #define kNumberOfImagesInMosaic 0x0019+(0x100A<< 16 ) //US NumberOfImagesInMosaic +#define kSeriesPlaneGE 0x0019+(0x1017<< 16 ) //SS #define kDwellTime 0x0019+(0x1018<< 16 ) //IS in NSec, see https://github.com/rordenlab/dcm2niix/issues/127 //https://nmrimaging.wordpress.com/2011/12/20/when-we-process/ // https://nciphub.org/groups/qindicom/wiki/DiffusionrelatedDICOMtags:experienceacrosssites?action=pdf -#define kDiffusionBValueSiemens 0x0019+(0x100C<< 16 ) //IS +#define kDiffusion_bValueSiemens 0x0019+(0x100C<< 16 ) //IS #define kSliceTimeSiemens 0x0019+(0x1029<< 16) ///FD #define kDiffusionGradientDirectionSiemens 0x0019+(0x100E<< 16 ) //FD #define kNumberOfDiffusionDirectionGE 0x0019+(0x10E0<< 16) ///DS NumberOfDiffusionDirection:UserData24 @@ -4237,6 +4251,8 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16); #define kDiffusionDirectionGEX 0x0019+(0x10BB<< 16 ) //DS phase diffusion direction #define kDiffusionDirectionGEY 0x0019+(0x10BC<< 16 ) //DS frequency diffusion direction #define kDiffusionDirectionGEZ 0x0019+(0x10BD<< 16 ) //DS slice diffusion direction +#define kPulseSequenceNameGE 0x0019+(0x109C<< 16 ) //LO 'epiRT' or 'epi' +#define kInternalPulseSequenceNameGE 0x0019+(0x109E<< 16 ) //LO 'EPI' or 'EPI2' #define kSharedFunctionalGroupsSequence 0x5200+uint32_t(0x9229<< 16 ) // SQ #define kPerFrameFunctionalGroupsSequence 0x5200+uint32_t(0x9230<< 16 ) // SQ #define kBandwidthPerPixelPhaseEncode 0x0019+(0x1028<< 16 ) //FD @@ -4250,7 +4266,9 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16); #define kOrientationACR 0x0020+(0x0035 << 16 ) //#define kTemporalPositionIdentifier 0x0020+(0x0100 << 16 ) //IS #define kOrientation 0x0020+(0x0037 << 16 ) -#define kTemporalPosition 0x0020+(0x0100 << 16 ) //IS +//#define kTemporalPosition 0x0020+(0x0100 << 16 ) //IS +//#define kNumberOfTemporalPositions 0x0020+(0x0105 << 16 ) //IS public tag for NumberOfDynamicScans +#define kTemporalResolution 0x0020+(0x0110 << 16 ) //DS #define kImagesInAcquisition 0x0020+(0x1002 << 16 ) //IS //#define kSliceLocation 0x0020+(0x1041 << 16 ) //DS would be useful if not optional type 3 #define kImageComments 0x0020+(0x4000<< 16 )// '0020' '4000' 'LT' 'ImageComments' @@ -4263,7 +4281,7 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16); #define kSequenceVariant21 0x0021+(0x105B<< 16 )//CS #define kPATModeText 0x0021+(0x1009<< 16 )//LO, see kImaPATModeText #define kTimeAfterStart 0x0021+(0x1104<< 16 )//DS -#define kPhaseEncodingDirectionPositive 0x0021+(0x111C<< 16 )//IS +#define kPhaseEncodingDirectionPositiveSiemens 0x0021+(0x111C<< 16 )//IS //#define kRealDwellTime 0x0021+(0x1142<< 16 )//IS #define kBandwidthPerPixelPhaseEncode21 0x0021+(0x1153<< 16 )//FD #define kCoilElements 0x0021+(0x114F<< 16 )//LO @@ -4276,7 +4294,7 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16); #define kPhotometricInterpretation 0x0028+(0x0004 << 16 ) #define kPlanarRGB 0x0028+(0x0006 << 16 ) #define kDim3 0x0028+(0x0008 << 16 ) //number of frames - for Philips this is Dim3*Dim4 -#define kDim2 0x0028+(0x0010 << 16 ) +#define kDim2 0x0028+(0x0010 << 16 ) #define kDim1 0x0028+(0x0011 << 16 ) #define kXYSpacing 0x0028+(0x0030 << 16 ) //DS 'PixelSpacing' #define kBitsAllocated 0x0028+(0x0100 << 16 ) @@ -4297,11 +4315,20 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16); #define kRealWorldSlope 0x0040+uint32_t(0x9225 << 16 ) //IS dicm2nii's SlopInt_6_9 #define kUserDefineDataGE 0x0043+(0x102A << 16 ) //OB #define kEffectiveEchoSpacingGE 0x0043+(0x102C << 16 ) //SS -#define kDiffusionBFactorGE 0x0043+(0x1039 << 16 ) //IS dicm2nii's SlopInt_6_9 +#define kImageTypeGE 0x0043+(0x102F << 16 ) //SS 0/1/2/3 for magnitude/phase/real/imaginary +#define kDiffusion_bValueGE 0x0043+(0x1039 << 16 ) //IS dicm2nii's SlopInt_6_9 +#define kEpiRTGroupDelayGE 0x0043+(0x107C << 16 ) //FL +#define kAssetRFactorsGE 0x0043+(0x1083 << 16 ) //DS +#define kMultiBandGE 0x0043+(0x10B6 << 16 ) //LO #define kAcquisitionMatrixText 0x0051+(0x100B << 16 ) //LO #define kCoilSiemens 0x0051+(0x100F << 16 ) #define kImaPATModeText 0x0051+(0x1011 << 16 ) #define kLocationsInAcquisition 0x0054+(0x0081 << 16 ) +#define kUnitsPT 0x0054+(0x1001<< 16 ) //CS +#define kAttenuationCorrectionMethod 0x0054+(0x1101<< 16 ) //LO +#define kDecayCorrection 0x0054+(0x1102<< 16 ) //CS +#define kReconstructionMethod 0x0054+(0x1103<< 16 ) //LO +#define kDecayFactor 0x0054+(0x1321<< 16 ) //LO //ftp://dicom.nema.org/MEDICAL/dicom/2014c/output/chtml/part03/sect_C.8.9.4.html //If ImageType is REPROJECTION we slice direction is reversed - need example to test // #define kSeriesType 0x0054+(0x1000 << 16 ) @@ -4311,17 +4338,22 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16); #define kDiffusion_bValueUIH 0x0065+(0x1009<< 16 ) //FD #define kParallelInformationUIH 0x0065+(0x100D<< 16 ) //SH #define kNumberOfImagesInGridUIH 0x0065+(0x1050<< 16 ) //DS -#define kMRVFrameSequenceUIH 0x0065+(0x1050<< 16 ) //SQ #define kDiffusionGradientDirectionUIH 0x0065+(0x1037<< 16 ) //FD +//#define kMRVFrameSequenceUIH 0x0065+(0x1050<< 16 ) //SQ +#define kPhaseEncodingDirectionPositiveUIH 0x0065+(0x1058<< 16 )//IS issue410 #define kIconImageSequence 0x0088+(0x0200 << 16 ) #define kElscintIcon 0x07a3+(0x10ce << 16 ) //see kGeiisFlag and https://github.com/rordenlab/dcm2niix/issues/239 #define kPMSCT_RLE1 0x07a1+(0x100a << 16 ) //Elscint/Philips compression -#define kPrivateCreator 0x2001+(0x0010 << 16 )// LO -#define kDiffusionBFactor 0x2001+(0x1003 << 16 )// FL +#define kPrivateCreator 0x2001+(0x0010 << 16 )// LO (Private creator is any tag where group is odd and element is x0010-x00FF +#define kDiffusion_bValuePhilips 0x2001+(0x1003 << 16 )// FL +#define kCardiacSync 0x2001+(0x1010 << 16 ) //CS //#define kDiffusionDirectionPhilips 0x2001+(0x1004 << 16 )//CS Diffusion Direction #define kSliceNumberMrPhilips 0x2001+(0x100A << 16 ) //IS Slice_Number_MR #define kSliceOrient 0x2001+(0x100B << 16 )//2001,100B Philips slice orientation (TRANSVERSAL, AXIAL, SAGITTAL) #define kEPIFactorPhilips 0x2001+(0x1013 << 16 ) //SL +#define kPrepulseDelay 0x2001+(0x101B << 16 ) //FL +#define kPrepulseType 0x2001+(0x101C << 16 ) //CS +#define kRespirationSync 0x2001+(0x101F << 16 ) //CS #define kNumberOfSlicesMrPhilips 0x2001+(0x1018 << 16 )//SL 0x2001, 0x1018 ), "Number_of_Slices_MR" #define kPartialMatrixScannedPhilips 0x2001+(0x1019 << 16 )// CS #define kWaterFatShiftPhilips 0x2001+(0x1022 << 16 ) //FL @@ -4352,8 +4384,28 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16); uint32_t kItemTag = 0xFFFE +(0xE000 << 16 ); uint32_t kItemDelimitationTag = 0xFFFE +(0xE00D << 16 ); uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); -double TE = 0.0; //most recent echo time recorded +#define salvageAgfa +#ifdef salvageAgfa //issue435 + // handle PrivateCreator renaming e.g. 0021,10xx -> 0021,11xx + // https://github.com/dcm4che/dcm4che/blob/master/dcm4che-dict/src/main/dicom3tools/libsrc/standard/elmdict/siemens.tpl + // https://github.com/neurolabusc/dcm_qa_agfa + // http://dicom.nema.org/medical/dicom/current/output/chtml/part05/sect_7.8.html + #define kMaxRemaps 16 //no vendor uses more than 5 private creator groups + //we need to keep track of multiple remappings, e.g. issue 437 2005,0014->2005,0012; 2005,0015->2005,0011 + int nRemaps = 0; + uint32_t privateCreatorMasks[kMaxRemaps]; //0 -> none + uint32_t privateCreatorRemaps[kMaxRemaps]; //0 -> none + //uint32_t privateCreatorMask = 0; //0 -> none + //uint32_t privateCreatorRemap = 0; //0 -> none +#endif + double TE = 0.0; //most recent echo time recorded + float temporalResolutionMS = 0.0; + float MRImageDynamicScanBeginTime = 0.0; bool is2005140FSQ = false; + bool overlayOK = true; + int overlayRows = 0; + int overlayCols = 0; + bool isTriggerSynced = false; bool isDICOMANON = false; //issue383 bool isMATLAB = false; //issue383 //double contentTime = 0.0; @@ -4406,6 +4458,7 @@ double TE = 0.0; //most recent echo time recorded } char vr[2]; //float intenScalePhilips = 0.0; + char seriesTimeTxt[kDICOMStr] = ""; char acquisitionDateTimeTxt[kDICOMStr] = ""; char imageType1st[kDICOMStr] = ""; bool isEncapsulatedData = false; @@ -4421,6 +4474,7 @@ double TE = 0.0; //most recent echo time recorded bool isPaletteColor = false; bool isInterpolated = false; bool isIconImageSequence = false; + int sqDepthIcon = -1; bool isSwitchToImplicitVR = false; bool isSwitchToBigEndian = false; bool isAtFirstPatientPosition = false; //for 3d and 4d files: flag is true for slices at same position as first slice @@ -4506,18 +4560,29 @@ double TE = 0.0; //most recent echo time recorded } //transfer syntax requests switching VR after group 0001 //uint32_t group = (groupElement & 0xFFFF); lPos += 4; - if ((groupElement == kItemDelimitationTag) || (groupElement == kSequenceDelimitationItemTag)) isIconImageSequence = false; + //issue409 - icons can have their own sub-sections... keep reading until we get to the icon image? + //if ((groupElement == kItemDelimitationTag) || (groupElement == kSequenceDelimitationItemTag)) isIconImageSequence = false; //if (groupElement == kItemTag) sqDepth++; bool unNest = false; while ((nNestPos > 0) && (nestPos[nNestPos] <= (lFileOffset+lPos))) { - nNestPos--; - sqDepth--; - unNest = true; + nNestPos--; + sqDepth--; + unNest = true; + if ((sqDepthIcon >= 0) && (sqDepth <= sqDepthIcon)) { //issue415 + sqDepthIcon = -1; + isIconImageSequence = false; + } + } if (groupElement == kItemDelimitationTag) { //end of item with undefined length sqDepth--; unNest = true; + if ((sqDepthIcon >= 0) && (sqDepth <= sqDepthIcon)) { //issue415 + sqDepthIcon = -1; + isIconImageSequence = false; + } } + if (unNest) { is2005140FSQ = false; if (sqDepth < 0) sqDepth = 0; //should not happen, but protect for faulty anonymization @@ -4578,7 +4643,10 @@ double TE = 0.0; //most recent echo time recorded dcmDim[numDimensionIndexValues].intenScalePhilips = d.intenScalePhilips; dcmDim[numDimensionIndexValues].RWVScale = d.RWVScale; dcmDim[numDimensionIndexValues].RWVIntercept = d.RWVIntercept; - dcmDim[numDimensionIndexValues].triggerDelayTime = d.triggerDelayTime; + if (isSameFloat(MRImageDynamicScanBeginTime * 1000.0, d.triggerDelayTime)) + dcmDim[numDimensionIndexValues].triggerDelayTime = 0.0; //issue395 + else + dcmDim[numDimensionIndexValues].triggerDelayTime = d.triggerDelayTime; dcmDim[numDimensionIndexValues].V[0] = -1.0; #ifdef MY_DEBUG if (numDimensionIndexValues < 19) { @@ -4788,10 +4856,98 @@ double TE = 0.0; //most recent echo time recorded } } if ((isIconImageSequence) && ((groupElement & 0x0028) == 0x0028 )) groupElement = kUnused; //ignore icon dimensions - switch ( groupElement ) { + #ifdef salvageAgfa //issue435 + //Handle remapping using integers, and slower but simpler approach is with strings: + // https://github.com/pydicom/pydicom/blob/master/pydicom/_private_dict.py + if (((groupElement & 65535) % 2) == 0) goto skipRemap; //remap odd (private) groups + //printf("tag %04x,%04x\n", groupElement & 65535, groupElement >> 16); + + if (((groupElement>>16) >= 0x10) && ((groupElement>>16) <= 0xFF)) { //tags (gggg,0010-00FF) may define new remapping + //if remapping tag + //first: see if this remapping overwrites existing tag + uint32_t privateCreatorMask = 0; //0 -> none + uint32_t privateCreatorRemap = 0; //0 -> none + privateCreatorMask = (groupElement & 65535) + ((groupElement & 0xFFFF0000) << 8); + if (nRemaps > 0) { + int j = 0; + for (int i = 0; i < nRemaps; i++) //remove duplicate remapping + //copy all remaps except exact match + if (privateCreatorMasks[i] != privateCreatorMask) { + privateCreatorMasks[j] = privateCreatorMasks[i]; + privateCreatorRemaps[j] = privateCreatorRemaps[i]; + j++; + } + nRemaps = j; + } + //see if this is known private vendor tag + privateCreatorRemap = 0; + char privateCreator[kDICOMStr]; + dcmStr(lLength, &buffer[lPos], privateCreator); + //next lines determine remapping, append as needed + //Siemens https://github.com/dcm4che/dcm4che/blob/master/dcm4che-dict/src/main/dicom3tools/libsrc/standard/elmdict/siemens.tpl + if (strstr(privateCreator, "SIEMENS MR HEADER") != NULL) privateCreatorRemap = 0x0019 +(0x1000 << 16 ); + if (strstr(privateCreator, "SIEMENS MR SDS 01") != NULL) privateCreatorRemap = 0x0021 +(0x1000 << 16 ); + if (strstr(privateCreator, "SIEMENS MR SDI 02") != NULL) privateCreatorRemap = 0x0021 +(0x1100 << 16 ); + if (strstr(privateCreator, "SIEMENS CSA HEADER") != NULL) privateCreatorRemap = 0x0029 +(0x1000 << 16 ); + if (strstr(privateCreator, "SIEMENS MR HEADER") != NULL) privateCreatorRemap = 0x0051 +(0x1000 << 16 ); + //GE https://github.com/dcm4che/dcm4che/blob/master/dcm4che-dict/src/main/dicom3tools/libsrc/standard/elmdict/gems.tpl + if (strstr(privateCreator, "GEMS_ACQU_01") != NULL) privateCreatorRemap = 0x0019 +(0x1000 << 16 ); + if (strstr(privateCreator, "GEMS_RELA_01") != NULL) privateCreatorRemap = 0x0021 +(0x1000 << 16 ); + if (strstr(privateCreator, "GEMS_SERS_01") != NULL) privateCreatorRemap = 0x0025 +(0x1000 << 16 ); + if (strstr(privateCreator, "GEMS_PARM_01") != NULL) privateCreatorRemap = 0x0043 +(0x1000 << 16 ); + //ELSCINT https://github.com/dcm4che/dcm4che/blob/master/dcm4che-dict/src/main/dicom3tools/libsrc/standard/elmdict/elscint.tpl + int grp = (groupElement & 65535); + if ((grp == 0x07a1) && (strstr(privateCreator, "ELSCINT1") != NULL)) privateCreatorRemap = 0x07a1 +(0x1000 << 16 ); + if ((grp == 0x07a3) && (strstr(privateCreator, "ELSCINT1") != NULL)) privateCreatorRemap = 0x07a3 +(0x1000 << 16 ); + //Philips https://github.com/dcm4che/dcm4che/blob/master/dcm4che-dict/src/main/dicom3tools/libsrc/standard/elmdict/philips.tpl + if (strstr(privateCreator, "PHILIPS IMAGING DD 001") != NULL) privateCreatorRemap = 0x2001 +(0x1000 << 16 ); + if (strstr(privateCreator, "Philips Imaging DD 001") != NULL) privateCreatorRemap = 0x2001 +(0x1000 << 16 ); + if (strstr(privateCreator, "PHILIPS MR IMAGING DD 001") != NULL) privateCreatorRemap = 0x2005 +(0x1000 << 16 ); + if (strstr(privateCreator, "Philips MR Imaging DD 001") != NULL) privateCreatorRemap = 0x2005 +(0x1000 << 16 ); + if (strstr(privateCreator, "PHILIPS MR IMAGING DD 005") != NULL) privateCreatorRemap = 0x2005 +(0x1400 << 16 ); + if (strstr(privateCreator, "Philips MR Imaging DD 005") != NULL) privateCreatorRemap = 0x2005 +(0x1400 << 16 ); + //UIH https://github.com/neurolabusc/dcm_qa_uih + if (strstr(privateCreator, "Image Private Header") != NULL) privateCreatorRemap = 0x0065 +(0x1000 << 16 ); + //sanity check: group should match + if (grp != (privateCreatorRemap & 65535)) privateCreatorRemap = 0; + if (privateCreatorRemap == 0) goto skipRemap; //this is not a known private group + if (privateCreatorRemap == privateCreatorMask) goto skipRemap; //the remapping and mask are identical 2005,1000 -> 2005,1000 + if ((nRemaps + 1) >=kMaxRemaps) goto skipRemap; //all slots full (should never happen) + //add new remapping + privateCreatorMasks[nRemaps] = privateCreatorMask; + privateCreatorRemaps[nRemaps] = privateCreatorRemap; + //printf("new remapping %04x,%04x -> %04x,%04x\n", privateCreatorMask & 65535, privateCreatorMask >> 16, privateCreatorRemap & 65535, privateCreatorRemap >> 16); + if (isVerbose > 1) + printf("new remapping (%d) %04x,%02xxy -> %04x,%02xxy\n", nRemaps, privateCreatorMask & 65535, privateCreatorMask >> 24, privateCreatorRemap & 65535, privateCreatorRemap >> 24); + nRemaps += 1; + //for (int i = 0; i < nRemaps; i++) + // printf(" %d = %04x,%02xxy -> %04x,%02xxy\n", i, privateCreatorMasks[i] & 65535, privateCreatorMasks[i] >> 24, privateCreatorRemaps[i] & 65535, privateCreatorRemaps[i] >> 24); + goto skipRemap; + } + + if (nRemaps < 1) goto skipRemap; + { + uint32_t remappedGroupElement = 0; + for (int i = 0; i < nRemaps; i++) + if ((groupElement & 0xFF00FFFF) == (privateCreatorMasks[i] & 0xFF00FFFF)) + remappedGroupElement = privateCreatorRemaps[i] + (groupElement & 0x00FF0000); + if (remappedGroupElement == 0) goto skipRemap; + if (isVerbose > 1) + printf("remapping %04x,%04x -> %04x,%04x\n", groupElement & 65535, groupElement >> 16, remappedGroupElement & 65535, remappedGroupElement >> 16); + groupElement = remappedGroupElement; + } + skipRemap: + #endif // salvageAgfa + if ((lLength % 2) != 0) { //https://www.nitrc.org/forum/forum.php?thread_id=11827&forum_id=4703 + printf("Illegal DICOM tag %04x,%04x (odd element length %d): %s\n", groupElement & 65535,groupElement>>16, lLength, fname); + //proper to return here, but we can carry on as a hail mary + // d.isValid = false; + //return d; + } + switch ( groupElement ) { case kMediaStorageSOPClassUID: { char mediaUID[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], mediaUID); + dcmStr(lLength, &buffer[lPos], mediaUID); //Philips "XX_" files //see https://github.com/rordenlab/dcm2niix/issues/328 if (strstr(mediaUID, "1.2.840.10008.5.1.4.1.1.66") != NULL) d.isRawDataStorage = true; @@ -4807,7 +4963,7 @@ double TE = 0.0; //most recent echo time recorded } case kMediaStorageSOPInstanceUID : {// 0002, 0003 //char SOPInstanceUID[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], d.instanceUID); + dcmStr(lLength, &buffer[lPos], d.instanceUID); //printMessage(">>%s\n", d.seriesInstanceUID); d.instanceUidCrc = mz_crc32X((unsigned char*) &d.instanceUID, strlen(d.instanceUID)); break; @@ -4815,7 +4971,7 @@ double TE = 0.0; //most recent echo time recorded case kTransferSyntax: { char transferSyntax[kDICOMStr]; strcpy(transferSyntax, ""); - dcmStr (lLength, &buffer[lPos], transferSyntax); + dcmStr(lLength, &buffer[lPos], transferSyntax); if (strcmp(transferSyntax, "1.2.840.10008.1.2.1") == 0) ; //default isExplicitVR=true; //d.isLittleEndian=true else if (strcmp(transferSyntax, "1.2.840.10008.1.2.4.50") == 0) { @@ -4870,14 +5026,14 @@ double TE = 0.0; //most recent echo time recorded break;} //{} provide scope for variable 'transferSyntax /*case kImplementationVersionName: { char impTxt[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], impTxt); + dcmStr(lLength, &buffer[lPos], impTxt); int slen = (int) strlen(impTxt); if((slen < 6) || (strstr(impTxt, "OSIRIX") == NULL) ) break; printError("OSIRIX Detected\n"); break; }*/ case kImplementationVersionName: { char impTxt[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], impTxt); + dcmStr(lLength, &buffer[lPos], impTxt); int slen = (int) strlen(impTxt); if ((slen > 5) && (strstr(impTxt, "MATLAB") != NULL) ) isMATLAB = true; @@ -4887,7 +5043,7 @@ double TE = 0.0; //most recent echo time recorded break; } case kSourceApplicationEntityTitle: { char saeTxt[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], saeTxt); + dcmStr(lLength, &buffer[lPos], saeTxt); int slen = (int) strlen(saeTxt); if((slen < 5) || (strstr(saeTxt, "oasis") == NULL) ) break; d.isSegamiOasis = true; @@ -4898,11 +5054,15 @@ double TE = 0.0; //most recent echo time recorded } case kImageTypeTag: { bool is1st = strlen(d.imageType) == 0; - dcmStr (lLength, &buffer[lPos], d.imageType, false, false); //<-distinguish spaces from pathdelim: [ORIGINAL\PHASE MAP\FFE] should return "PHASE MAP" not "PHASE_MAP" - if (is1st) - strcpy(imageType1st, d.imageType); + dcmStr(lLength, &buffer[lPos], d.imageType, false); //<-distinguish spaces from pathdelim: [ORIGINAL\PHASE MAP\FFE] should return "PHASE MAP" not "PHASE_MAP" int slen; slen = (int) strlen(d.imageType); + if (slen > 1) { + for (int i = 0; i 5) && strstr(d.imageType, "_MOCO_") ) { //d.isDerived = true; //this would have 'i- y' skip MoCo images isMoCo = true; @@ -4972,16 +5132,16 @@ double TE = 0.0; //most recent echo time recorded break; } case kAcquisitionDate: char acquisitionDateTxt[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], acquisitionDateTxt); + dcmStr(lLength, &buffer[lPos], acquisitionDateTxt); d.acquisitionDate = atof(acquisitionDateTxt); break; case kAcquisitionDateTime: //char acquisitionDateTimeTxt[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], acquisitionDateTimeTxt); + dcmStr(lLength, &buffer[lPos], acquisitionDateTimeTxt); //printMessage("%s\n",acquisitionDateTimeTxt); break; case kStudyDate: - dcmStr (lLength, &buffer[lPos], d.studyDate); + dcmStr(lLength, &buffer[lPos], d.studyDate); break; case kModality: if (lLength < 2) break; @@ -5004,7 +5164,7 @@ double TE = 0.0; //most recent echo time recorded dcmStr(lLength, &buffer[lPos], d.institutionName); break; case kInstitutionAddress: //VR is "ST": 1024 chars maximum - dcmStr(lLength, &buffer[lPos], d.institutionAddress, true); + dcmStr(lLength, &buffer[lPos], d.institutionAddress); break; case kReferringPhysicianName: dcmStr(lLength, &buffer[lPos], d.referringPhysicianName); @@ -5038,97 +5198,125 @@ double TE = 0.0; //most recent echo time recorded if (((int) strlen(acqContrast) > 8) && (strstr(acqContrast, "DIFFUSION") != NULL)) d.isDiffusion = true; break; - case kAcquisitionTime : + case kAcquisitionTime : { char acquisitionTimeTxt[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], acquisitionTimeTxt); + dcmStr(lLength, &buffer[lPos], acquisitionTimeTxt); d.acquisitionTime = atof(acquisitionTimeTxt); if (d.manufacturer != kMANUFACTURER_UIH) break; //UIH slice timing- do not use for Siemens as Siemens de-identification can corrupt this field https://github.com/rordenlab/dcm2niix/issues/236 d.CSA.sliceTiming[acquisitionTimesGE_UIH] = d.acquisitionTime; acquisitionTimesGE_UIH ++; - break; + break; } //case kContentTime : // char contentTimeTxt[kDICOMStr]; - // dcmStr (lLength, &buffer[lPos], contentTimeTxt); + // dcmStr(lLength, &buffer[lPos], contentTimeTxt); // contentTime = atof(contentTimeTxt); // break; + case kSeriesTime : + dcmStr(lLength, &buffer[lPos], seriesTimeTxt); + break; case kStudyTime : - dcmStr (lLength, &buffer[lPos], d.studyTime); + dcmStr(lLength, &buffer[lPos], d.studyTime); break; case kPatientName : - dcmStr (lLength, &buffer[lPos], d.patientName); + dcmStr(lLength, &buffer[lPos], d.patientName); break; case kAnatomicalOrientationType: { char aotTxt[kDICOMStr]; //ftp://dicom.nema.org/MEDICAL/dicom/2015b/output/chtml/part03/sect_C.7.6.2.html#sect_C.7.6.2.1.1 - dcmStr (lLength, &buffer[lPos], aotTxt); + dcmStr(lLength, &buffer[lPos], aotTxt); int slen = (int) strlen(aotTxt); if((slen < 9) || (strstr(aotTxt, "QUADRUPED") == NULL) ) break; printError("Anatomical Orientation Type (0010,2210) is QUADRUPED: rotate coordinates accordingly\n"); break; } case kDeidentificationMethod: { //issue 383 char anonTxt[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], anonTxt); + dcmStr(lLength, &buffer[lPos], anonTxt); int slen = (int) strlen(anonTxt); if((slen < 10) || (strstr(anonTxt, "DICOMANON") == NULL) ) break; isDICOMANON = true; printWarning("Matlab DICOMANON can scramble SeriesInstanceUID (0020,000e) and remove crucial data (see issue 383). \n"); break; } case kPatientID : - dcmStr (lLength, &buffer[lPos], d.patientID); + if (strlen(d.patientID) > 1) break; + dcmStr(lLength, &buffer[lPos], d.patientID); break; case kAccessionNumber : - dcmStr (lLength, &buffer[lPos], d.accessionNumber); + dcmStr(lLength, &buffer[lPos], d.accessionNumber); break; case kPatientBirthDate : - dcmStr (lLength, &buffer[lPos], d.patientBirthDate); + dcmStr(lLength, &buffer[lPos], d.patientBirthDate); break; - case kPatientSex : - d.patientSex = toupper(buffer[lPos]); //first character is either 'R'ow or 'C'ol + case kPatientSex : { + //must be M,F,O: http://dicom.nema.org/dicom/2013/output/chtml/part03/sect_C.2.html + char patientSex = toupper(buffer[lPos]); + if ((patientSex == 'M') || (patientSex == 'F') || (patientSex == 'O')) + d.patientSex = patientSex; break; + } case kPatientAge : - dcmStr (lLength, &buffer[lPos], d.patientAge); + dcmStr(lLength, &buffer[lPos], d.patientAge); break; case kPatientWeight : d.patientWeight = dcmStrFloat(lLength, &buffer[lPos]); break; case kStationName : - dcmStr (lLength, &buffer[lPos], d.stationName); + dcmStr(lLength, &buffer[lPos], d.stationName); + break; + case kSeriesDescription: + dcmStr(lLength, &buffer[lPos], d.seriesDescription); break; - case kSeriesDescription: { - dcmStr (lLength, &buffer[lPos], d.seriesDescription); - break; } case kInstitutionalDepartmentName: - dcmStr (lLength, &buffer[lPos], d.institutionalDepartmentName); + dcmStr(lLength, &buffer[lPos], d.institutionalDepartmentName); break; case kManufacturersModelName : - dcmStr (lLength, &buffer[lPos], d.manufacturersModelName); + dcmStr(lLength, &buffer[lPos], d.manufacturersModelName); break; case kDerivationDescription : { //strcmp(transferSyntax, "1.2.840.10008.1.2") char derivationDescription[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], derivationDescription);//strcasecmp, strcmp + dcmStr(lLength, &buffer[lPos], derivationDescription);//strcasecmp, strcmp if (strcasecmp(derivationDescription, "MEDCOM_RESAMPLED") == 0) d.isResampled = true; break; } case kDeviceSerialNumber : { - dcmStr (lLength, &buffer[lPos], d.deviceSerialNumber); + dcmStr(lLength, &buffer[lPos], d.deviceSerialNumber); break; } case kSoftwareVersions : { - dcmStr (lLength, &buffer[lPos], d.softwareVersions); + dcmStr(lLength, &buffer[lPos], d.softwareVersions); int slen = (int) strlen(d.softwareVersions); if((slen > 4) && (strstr(d.softwareVersions, "XA11") != NULL) ) d.isXA10A = true; + if((slen > 4) && (strstr(d.softwareVersions, "XA20") != NULL) ) d.isXA10A = true; + if((slen > 4) && (strstr(d.softwareVersions, "XA30") != NULL) ) d.isXA10A = true; if((slen < 5) || (strstr(d.softwareVersions, "XA10") == NULL) ) break; - d.isXA10A = true; + d.isXA10A = true; break; } case kProtocolName : { //if ((strlen(d.protocolName) < 1) || (d.manufacturer != kMANUFACTURER_GE)) //GE uses a generic session name here: do not overwrite kProtocolNameGE - dcmStr (lLength, &buffer[lPos], d.protocolName); //see also kSequenceName + dcmStr(lLength, &buffer[lPos], d.protocolName); //see also kSequenceName break; } case kPatientOrient : - dcmStr (lLength, &buffer[lPos], d.patientOrient); + dcmStr(lLength, &buffer[lPos], d.patientOrient); + break; + case kInversionRecovery : // CS [YES],[NO] + if (lLength < 2) break; + if (toupper(buffer[lPos]) == 'Y') + d.isIR = true; + break; + case kEchoPlanarPulseSequence : // CS [YES],[NO] + if (lLength < 2) break; + if (toupper(buffer[lPos]) == 'Y') + d.isEPI = true; break; + case kRectilinearPhaseEncodeReordering : { //'CS' [REVERSE_LINEAR],[LINEAR],[CENTRIC],[REVERSE_CENTRIC] + if (d.manufacturer != kMANUFACTURER_GE) break; //only found in GE software beginning with RX27 + if (lLength < 2) break; + if (toupper(buffer[lPos]) == 'L') + d.phaseEncodingGE = kGE_PHASE_ENCODING_POLARITY_UNFLIPPED; + if (toupper(buffer[lPos]) == 'R') + d.phaseEncodingGE = kGE_PHASE_ENCODING_POLARITY_FLIPPED; + break; } case kParallelReductionFactorInPlane: if (d.manufacturer == kMANUFACTURER_SIEMENS) break; d.accelFactPE = dcmFloatDouble(lLength, &buffer[lPos],d.isLittleEndian); @@ -5142,14 +5330,14 @@ double TE = 0.0; //most recent echo time recorded // //(0018,9074) DT [20190621095516.140000] YYYYMMDDHHMMSS // //see https://github.com/rordenlab/dcm2niix/issues/303 // char dateTime[kDICOMStr]; - // dcmStr (lLength, &buffer[lPos], dateTime); + // dcmStr(lLength, &buffer[lPos], dateTime); // printf("%s\tkFrameAcquisitionDateTime\n", dateTime); //} case kDiffusionDirectionality : {// 0018, 9075 set_directionality0018_9075(&volDiffusion, (&buffer[lPos])); if ((d.manufacturer != kMANUFACTURER_PHILIPS) || (lLength < 10)) break; char dir[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], dir); + dcmStr(lLength, &buffer[lPos], dir); if (strcmp(dir, "ISOTROPIC") == 0) isPhilipsDerived = true; break; } @@ -5186,13 +5374,18 @@ double TE = 0.0; //most recent echo time recorded if (d.manufacturer == kMANUFACTURER_SIEMENS) numberOfImagesInMosaic = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); break; + case kSeriesPlaneGE : //SS 2=Axi, 4=Sag, 8=Cor, 16=Obl, 256=3plane + if (d.manufacturer != kMANUFACTURER_GE) break; + if (dcmInt(lLength,&buffer[lPos],d.isLittleEndian) == 256) d.isLocalizer = true; + break; case kDwellTime : d.dwellTime = dcmStrInt(lLength, &buffer[lPos]); break; - case kDiffusionBValueSiemens : + case kDiffusion_bValueSiemens : if (d.manufacturer != kMANUFACTURER_SIEMENS) break; - d.CSA.dtiV[0] = dcmStrInt(lLength, &buffer[lPos]); - d.CSA.numDti = 1; + //issue409 + B0Philips = dcmStrInt(lLength, &buffer[lPos]); + set_bVal(&volDiffusion, B0Philips); break; case kSliceTimeSiemens : {//Array of FD (64-bit double) if (d.manufacturer != kMANUFACTURER_SIEMENS) break; @@ -5261,14 +5454,44 @@ double TE = 0.0; //most recent echo time recorded if (d.manufacturer == kMANUFACTURER_GE) set_diffusion_directionGE(&volDiffusion, lLength, (&buffer[lPos]), 2); break; + case kPulseSequenceNameGE : { //LO 'epi'/'epiRT' + if (d.manufacturer != kMANUFACTURER_GE) break; + char epiStr[kDICOMStr]; + dcmStr(lLength, &buffer[lPos], epiStr); + if ((strstr(epiStr, "epi") != NULL) && (strstr(epiStr, "epi2") == NULL)){ + d.epiVersionGE = 0; //-1 = not epi, 0 = epi, 1 = epiRT + } + if (strstr(epiStr, "epi2") != NULL){ + d.epiVersionGE = 2; //-1 = not epi, 0 = epi, 1 = epiRT, 2 = epi2 + } + if (strstr(epiStr, "epiRT") != NULL){ + d.epiVersionGE = 1; //-1 = not epi, 0 = epi, 1 = epiRT + } + break; + } + case kInternalPulseSequenceNameGE : { //LO 'EPI'(gradient echo)/'EPI2'(spin echo): + if (d.manufacturer != kMANUFACTURER_GE) break; + char epiStr[kDICOMStr]; + dcmStr(lLength, &buffer[lPos], epiStr); + if (strcmp(epiStr, "EPI") == 0){ + d.internalepiVersionGE = 1; //-1 = not EPI, 1 = EPI, 2 = EPI2 + if (d.epiVersionGE != 1){ // 1 = epiRT by kEpiRTGroupDelayGE or kPulseSequenceNameGE + d.epiVersionGE = 0; // 0 = epi (multi-phase epi) + } + } + if (strcmp(epiStr, "EPI2") == 0){ + d.internalepiVersionGE = 2; //-1 = not epi, 1 = EPI, 2 = EPI2 + } + break; + } case kBandwidthPerPixelPhaseEncode: d.bandwidthPerPixelPhaseEncode = dcmFloatDouble(lLength, &buffer[lPos],d.isLittleEndian); break; - case kStudyInstanceUID : // 0020, 000D - dcmStr (lLength, &buffer[lPos], d.studyInstanceUID); + case kStudyInstanceUID : // 0020,000D + dcmStr(lLength, &buffer[lPos], d.studyInstanceUID); break; - case kSeriesInstanceUID : // 0020, 000E - dcmStr (lLength, &buffer[lPos], d.seriesInstanceUID); + case kSeriesInstanceUID : // 0020,000E + dcmStr(lLength, &buffer[lPos], d.seriesInstanceUID); //printMessage(">>%s\n", d.seriesInstanceUID); d.seriesUidCrc = mz_crc32X((unsigned char*) &d.seriesInstanceUID, strlen(d.seriesInstanceUID)); break; @@ -5280,9 +5503,8 @@ double TE = 0.0; //most recent echo time recorded patientPositionNum++; isAtFirstPatientPosition = true; //char dx[kDICOMStr]; - //dcmStr (lLength, &buffer[lPos], dx); + //dcmStr(lLength, &buffer[lPos], dx); //printMessage("*%s*", dx); - dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, &patientPosition[0]); //slice position if (isnan(d.patientPosition[1])) { //dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, &d.patientPosition[0]); //slice position @@ -5319,7 +5541,6 @@ double TE = 0.0; //most recent echo time recorded } nSliceMM++; } - break; } case kInPlanePhaseEncodingDirection: d.phaseEncodingRC = toupper(buffer[lPos]); //first character is either 'R'ow or 'C'ol @@ -5328,7 +5549,7 @@ double TE = 0.0; //most recent echo time recorded d.SAR = dcmStrFloat(lLength, &buffer[lPos]); break; case kStudyID: - dcmStr (lLength, &buffer[lPos], d.studyID); + dcmStr(lLength, &buffer[lPos], d.studyID); break; case kSeriesNum: d.seriesNum = dcmStrInt(lLength, &buffer[lPos]); @@ -5375,7 +5596,7 @@ double TE = 0.0; //most recent echo time recorded break; } case kPhotometricInterpretation: { char interp[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], interp); + dcmStr(lLength, &buffer[lPos], interp); if (strcmp(interp, "PALETTE_COLOR") == 0) isPaletteColor = true; //printError("Photometric Interpretation 'PALETTE COLOR' not supported\n"); @@ -5407,20 +5628,20 @@ double TE = 0.0; //most recent echo time recorded // dcmMultiFloat(lLength, (char*)&buffer[lPos], 2, d.xyzMM); // break; case kImageComments: - dcmStr (lLength, &buffer[lPos], d.imageComments, true); + dcmStr(lLength, &buffer[lPos], d.imageComments, true); break; //group 21: siemens //g21 case kPATModeText : { //e.g. Siemens iPAT x2 listed as "p2" char accelStr[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], accelStr); + dcmStr(lLength, &buffer[lPos], accelStr); char *ptr; dcmStrDigitsOnlyKey('p', accelStr); //e.g. if "p2s4" return "2", if "s4" return "" d.accelFactPE = (float)strtof(accelStr, &ptr); if (*ptr != '\0') d.accelFactPE = 0.0; //between slice accel - dcmStr (lLength, &buffer[lPos], accelStr); + dcmStr(lLength, &buffer[lPos], accelStr); dcmStrDigitsOnlyKey('s', accelStr); //e.g. if "p2s4" return "4", if "p2" return "" multiBandFactor = (int)strtol(accelStr, &ptr, 10); if (*ptr != '\0') @@ -5438,7 +5659,7 @@ double TE = 0.0; //most recent echo time recorded //printf("x\t%d\t%g\tkTimeAfterStart\n", acquisitionTimesGE_UIH, d.CSA.sliceTiming[acquisitionTimesGE_UIH]); acquisitionTimesGE_UIH ++; break; - case kPhaseEncodingDirectionPositive: { + case kPhaseEncodingDirectionPositiveSiemens: { if (d.manufacturer != kMANUFACTURER_SIEMENS) break; int ph = dcmStrInt(lLength, &buffer[lPos]); if (ph == 0) d.phaseEncodingGE = kGE_PHASE_ENCODING_POLARITY_FLIPPED; @@ -5454,7 +5675,7 @@ double TE = 0.0; //most recent echo time recorded break; case kCoilElements: if (d.manufacturer != kMANUFACTURER_SIEMENS) break; - dcmStr (lLength, &buffer[lPos], d.coilElements); + dcmStr(lLength, &buffer[lPos], d.coilElements); break; //group 21: GE case kLocationsInAcquisitionGE: @@ -5480,20 +5701,19 @@ double TE = 0.0; //most recent echo time recorded break; case kPEDirectionDisplayedUIH : if (d.manufacturer != kMANUFACTURER_UIH) break; - dcmStr (lLength, &buffer[lPos], d.phaseEncodingDirectionDisplayedUIH); + dcmStr(lLength, &buffer[lPos], d.phaseEncodingDirectionDisplayedUIH); break; case kDiffusion_bValueUIH : { if (d.manufacturer != kMANUFACTURER_UIH) break; float v[4]; dcmMultiFloatDouble(lLength, &buffer[lPos], 1, v, d.isLittleEndian); - d.CSA.dtiV[0] = v[0]; - d.CSA.numDti = 1; - //printf("%d>>>%g\n", lPos, v[0]); + B0Philips = v[0]; + set_bVal(&volDiffusion, v[0]); break; } case kParallelInformationUIH: {//SENSE factor (0065,100d) SH [F:2S] if (d.manufacturer != kMANUFACTURER_UIH) break; char accelStr[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], accelStr); + dcmStr(lLength, &buffer[lPos], accelStr); //char *ptr; dcmStrDigitsDotOnlyKey(':', accelStr); //e.g. if "p2s4" return "2", if "s4" return "" d.accelFactPE = atof(accelStr); @@ -5503,6 +5723,12 @@ double TE = 0.0; //most recent echo time recorded d.numberOfImagesInGridUIH = dcmStrFloat(lLength, &buffer[lPos]); d.CSA.mosaicSlices = d.numberOfImagesInGridUIH; break; + case kPhaseEncodingDirectionPositiveUIH: { + if (d.manufacturer != kMANUFACTURER_UIH) break; + int ph = dcmStrInt(lLength, &buffer[lPos]); + if (ph == 1) d.phaseEncodingGE = kGE_PHASE_ENCODING_POLARITY_FLIPPED; + if (ph == 0) d.phaseEncodingGE = kGE_PHASE_ENCODING_POLARITY_UNFLIPPED; + break; } case kDiffusionGradientDirectionUIH : { //0065,1037 //0.03712929804225321\-0.5522387869760447\-0.8328587749392602 if (d.manufacturer != kMANUFACTURER_UIH) break; @@ -5559,6 +5785,9 @@ double TE = 0.0; //most recent echo time recorded //printf("%g\n", d.CSA.sliceTiming[acquisitionTimesGE_UIH]); acquisitionTimesGE_UIH ++; break; } + case kRadionuclideTotalDose : + d.radionuclideTotalDose = dcmStrFloat(lLength, &buffer[lPos]); + break; case kEffectiveTE : { TE = dcmFloatDouble(lLength, &buffer[lPos], d.isLittleEndian); if (d.TE <= 0.0) @@ -5616,9 +5845,6 @@ double TE = 0.0; //most recent echo time recorded case kFlipAngle : d.flipAngle = dcmStrFloat(lLength, &buffer[lPos]); break; - case kRadionuclideTotalDose : - d.radionuclideTotalDose = dcmStrFloat(lLength, &buffer[lPos]); - break; case kRadionuclideHalfLife : d.radionuclideHalfLife = dcmStrFloat(lLength, &buffer[lPos]); break; @@ -5634,8 +5860,14 @@ double TE = 0.0; //most recent echo time recorded d.TE = dcmStrFloat(lLength, &buffer[lPos]); } break; + case kConvolutionKernel: //CS + dcmStr(lLength, &buffer[lPos], d.convolutionKernel); + break; + case kFrameDuration : + d.frameDuration = dcmStrInt(lLength, &buffer[lPos]); + break; case kReceiveCoilName : - dcmStr (lLength, &buffer[lPos], d.coilName); + dcmStr(lLength, &buffer[lPos], d.coilName); if (strlen(d.coilName) < 1) break; d.coilCrc = mz_crc32X((unsigned char*) &d.coilName, strlen(d.coilName)); break; @@ -5652,14 +5884,17 @@ double TE = 0.0; //most recent echo time recorded case kMRImageDynamicScanBeginTime: { //FL if (lLength != 4) break; - float dyn = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); - if (dyn < minDynamicScanBeginTime) minDynamicScanBeginTime = dyn; - if (dyn > maxDynamicScanBeginTime) maxDynamicScanBeginTime = dyn; + MRImageDynamicScanBeginTime = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); + if (MRImageDynamicScanBeginTime < minDynamicScanBeginTime) minDynamicScanBeginTime = MRImageDynamicScanBeginTime; + if (MRImageDynamicScanBeginTime > maxDynamicScanBeginTime) maxDynamicScanBeginTime = MRImageDynamicScanBeginTime; break; } case kIntercept : d.intenIntercept = dcmStrFloat(lLength, &buffer[lPos]); break; + case kRadiopharmaceutical : + dcmStr(lLength, &buffer[lPos], d.radiopharmaceutical); + break; case kZThick : d.xyzMM[3] = dcmStrFloat(lLength, &buffer[lPos]); d.zThick = d.xyzMM[3]; @@ -5669,7 +5904,7 @@ double TE = 0.0; //most recent echo time recorded case kAcquisitionMatrixText : { if (d.manufacturer == kMANUFACTURER_SIEMENS) { char matStr[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], matStr); + dcmStr(lLength, &buffer[lPos], matStr); char* pPosition = strchr(matStr, 'I'); if (pPosition != NULL) isInterpolated = true; @@ -5680,7 +5915,7 @@ double TE = 0.0; //most recent echo time recorded //see if image from single coil "H12" or an array "HEA;HEP" //char coilStr[kDICOMStr]; //int coilNum; - dcmStr (lLength, &buffer[lPos], d.coilName); + dcmStr(lLength, &buffer[lPos], d.coilName); if (strlen(d.coilName) < 1) break; //printf("-->%s\n", coilStr); //d.coilName = coilStr; @@ -5697,14 +5932,14 @@ double TE = 0.0; //most recent echo time recorded break; } case kImaPATModeText : { //e.g. Siemens iPAT x2 listed as "p2" char accelStr[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], accelStr); + dcmStr(lLength, &buffer[lPos], accelStr); char *ptr; dcmStrDigitsOnlyKey('p', accelStr); //e.g. if "p2s4" return "2", if "s4" return "" d.accelFactPE = (float)strtof(accelStr, &ptr); if (*ptr != '\0') d.accelFactPE = 0.0; //between slice accel - dcmStr (lLength, &buffer[lPos], accelStr); + dcmStr(lLength, &buffer[lPos], accelStr); dcmStrDigitsOnlyKey('s', accelStr); //e.g. if "p2s4" return "4", if "p2" return "" multiBandFactor = (int)strtol(accelStr, &ptr, 10); if (*ptr != '\0') @@ -5713,8 +5948,24 @@ double TE = 0.0; //most recent echo time recorded case kLocationsInAcquisition : d.locationsInAcquisition = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); break; + case kUnitsPT: //CS + dcmStr(lLength, &buffer[lPos], d.unitsPT); + break; + case kAttenuationCorrectionMethod: //LO + dcmStr(lLength, &buffer[lPos], d.attenuationCorrectionMethod); + break; + case kDecayCorrection: //CS + dcmStr(lLength, &buffer[lPos], d.decayCorrection); + break; + case kReconstructionMethod: //LO + dcmStr(lLength, &buffer[lPos], d.reconstructionMethod); + break; + case kDecayFactor : + d.decayFactor = dcmStrFloat(lLength, &buffer[lPos]); + break; case kIconImageSequence: isIconImageSequence = true; + if (sqDepthIcon < 0) sqDepthIcon = sqDepth; break; /*case kStackSliceNumber: { //https://github.com/Kevin-Mattheus-Moerman/GIBBON/blob/master/dicomDict/PMS-R32-dict.txt int stackSliceNumber = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); @@ -5731,29 +5982,42 @@ double TE = 0.0; //most recent echo time recorded case kMRAcquisitionType: //detect 3D acquisition: we can reorient these without worrying about slice time correct or BVEC/BVAL orientation if (lLength > 1) d.is2DAcq = (buffer[lPos]=='2') && (toupper(buffer[lPos+1]) == 'D'); if (lLength > 1) d.is3DAcq = (buffer[lPos]=='3') && (toupper(buffer[lPos+1]) == 'D'); - //dcmStr (lLength, &buffer[lPos], d.mrAcquisitionType); + //dcmStr(lLength, &buffer[lPos], d.mrAcquisitionType); break; case kBodyPartExamined : { - dcmStr (lLength, &buffer[lPos], d.bodyPartExamined); + dcmStr(lLength, &buffer[lPos], d.bodyPartExamined); break; } case kScanningSequence : { - dcmStr (lLength, &buffer[lPos], d.scanningSequence); - break; + dcmStr(lLength, &buffer[lPos], d.scanningSequence); + //According to the DICOM standard 0018,9018 is REQUIRED for EPI raw data + // http://dicom.nema.org/MEDICAL/Dicom/2015c/output/chtml/part03/sect_C.8.13.4.html + //In practice, this is not the case for all vendors + //Fortunately, the combination of 0018,0020 and 0018,9018 appears to reliably detect EPI data + //Siemens (pre-XA) omits 0018,9018, but reports [EP] for 0018,0020 (regardless of SE/GR) + //Siemens (XA) reports 0018,9018 but omits 0018,0020 + //Canon/Toshiba omits 0018,9018, but reports [SE\EP];[GR\EP] for 0018,0020 + //GE omits 0018,9018, but reports [EP\GR];[EP\SE] for 0018,0020 + //Philips reports 0018,9018, but reports [SE];[GR] for 0018,0020 + if ((lLength > 1) && (strstr(d.scanningSequence, "IR") != NULL)) + d.isIR = true; + if ((lLength > 1) && (strstr(d.scanningSequence, "EP") != NULL)) + d.isEPI = true; + break; //warp } case kSequenceVariant21 : if (d.manufacturer != kMANUFACTURER_SIEMENS) break; //see GE dataset in dcm_qa_nih //fall through... case kSequenceVariant : { - dcmStr (lLength, &buffer[lPos], d.sequenceVariant); + dcmStr(lLength, &buffer[lPos], d.sequenceVariant); break; } case kScanOptions: - dcmStr (lLength, &buffer[lPos], d.scanOptions); + dcmStr(lLength, &buffer[lPos], d.scanOptions); break; case kSequenceName : { //if (strlen(d.protocolName) < 1) //precedence given to kProtocolName and kProtocolNameGE - dcmStr (lLength, &buffer[lPos], d.sequenceName); + dcmStr(lLength, &buffer[lPos], d.sequenceName); break; } case kMRAcquisitionTypePhilips: //kMRAcquisitionType @@ -5780,7 +6044,7 @@ double TE = 0.0; //most recent echo time recorded case kSliceOrient: { char orientStr[kDICOMStr]; orientStr[0] = 'X'; //avoid compiler warning: orientStr filled by dcmStr - dcmStr (lLength, &buffer[lPos], orientStr); + dcmStr(lLength, &buffer[lPos], orientStr); if (toupper(orientStr[0])== 'S') d.sliceOrient = kSliceOrientSag; //sagittal else if (toupper(orientStr[0])== 'C') @@ -5791,6 +6055,7 @@ double TE = 0.0; //most recent echo time recorded case kElscintIcon : printWarning("Assuming icon SQ 07a3,10ce.\n"); isIconImageSequence = true; + if (sqDepthIcon < 0) sqDepthIcon = sqDepth; break; case kPMSCT_RLE1 : //https://groups.google.com/forum/#!topic/comp.protocols.dicom/8HuP_aNy9Pc @@ -5804,9 +6069,8 @@ double TE = 0.0; //most recent echo time recorded if (d.manufacturer != kMANUFACTURER_UNKNOWN) break; d.manufacturer = dcmStrManufacturer (lLength, &buffer[lPos]); volDiffusion.manufacturer = d.manufacturer; - //printf(">>>>%d\n", d.manufacturer); break; } - case kDiffusionBFactor : + case kDiffusion_bValuePhilips : if (d.manufacturer != kMANUFACTURER_PHILIPS) break; B0Philips = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); set_bVal(&volDiffusion, B0Philips); @@ -5836,18 +6100,23 @@ double TE = 0.0; //most recent echo time recorded //CS: Possible values: P (PreparationDirection), M (MeasurementDirection),S (Selection Direction),O(Oblique Direction),I (Isotropic),Only applicable for diffusion scans. if (d.manufacturer != kMANUFACTURER_PHILIPS) break; char diffDir[kDICOMStr]; - dcmStr (lLength, &buffer[lPos], diffDir); + dcmStr(lLength, &buffer[lPos], diffDir); printf(">>%s %s\n", diffDir, fname); break; } */ + case kCardiacSync : //CS [TRIGGERED],[NO] + if (lLength < 2) break; + if (toupper(buffer[lPos]) != 'N') + isTriggerSynced = true; + break; case kDiffusion_bValue: // 0018,9087 if (d.manufacturer == kMANUFACTURER_UNKNOWN ) { d.manufacturer = kMANUFACTURER_PHILIPS; printWarning("Found 0018,9087 but manufacturer (0008,0070) unknown: assuming Philips.\n"); } - // Note that this is ahead of kPatientPosition (0020,0032), so + // Note that this is ahead of kImagePositionPatient (0020,0032), so // isAtFirstPatientPosition is not necessarily set yet. // Philips uses this tag too, at least as of 5.1, but they also // use kDiffusionBFactor (see above), and we do not want to @@ -5874,7 +6143,7 @@ double TE = 0.0; //most recent echo time recorded //} break; case kDiffusionOrientation: // 0018, 9089 - // Note that this is ahead of kPatientPosition (0020,0032), so + // Note that this is ahead of kImagePositionPatient (0020,0032), so // isAtFirstPatientPosition is not necessarily set yet. // Philips uses this tag too, at least as of 5.1, but they also // use kDiffusionDirectionRL, etc., and we do not want to double @@ -5917,6 +6186,10 @@ double TE = 0.0; //most recent echo time recorded case kImagingFrequency2 : d.imagingFrequency = dcmFloatDouble(lLength, &buffer[lPos], d.isLittleEndian); break; + case kParallelReductionFactorOutOfPlane: + if (d.manufacturer == kMANUFACTURER_SIEMENS) break; + d.accelFactOOP = dcmFloatDouble(lLength, &buffer[lPos],d.isLittleEndian); + break; //case kFrameAcquisitionDuration : // frameAcquisitionDuration = dcmFloatDouble(lLength, &buffer[lPos], d.isLittleEndian); //issue369 // break; @@ -5977,6 +6250,23 @@ double TE = 0.0; //most recent echo time recorded break; echoTrainLengthPhil = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); break; + case kPrepulseDelay : //FL + if (d.manufacturer != kMANUFACTURER_PHILIPS) + break; + d.TI = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); + break; + case kPrepulseType : //CS [INV] + if (d.manufacturer != kMANUFACTURER_PHILIPS) + break; + if (lLength < 3) break; + if ((toupper(buffer[lPos]) != 'I') && (toupper(buffer[lPos+1]) != 'N') && (toupper(buffer[lPos+2]) != 'V')) + d.isIR = true; + break; + case kRespirationSync : //CS [TRIGGERED],[NO] + if (lLength < 2) break; + if (toupper(buffer[lPos]) != 'N') + isTriggerSynced = true; + break; case kNumberOfSlicesMrPhilips : if (d.manufacturer != kMANUFACTURER_PHILIPS) break; @@ -6086,8 +6376,9 @@ double TE = 0.0; //most recent echo time recorded case kRealWorldSlope: if (d.manufacturer != kMANUFACTURER_PHILIPS) break; d.RWVScale = dcmFloatDouble(lLength, &buffer[lPos],d.isLittleEndian); - //printMessage("RWVScale %g\n", d.RWVScale); - if (isSameFloat(1.0, d.intenScale)) //give precedence to standard value + if (d.RWVScale > 1.0E38) + d.RWVScale = 0.0; + else if (isSameFloat(1.0, d.intenScale)) //give precedence to standard value d.intenScale = d.RWVScale; break; case kUserDefineDataGE: { //0043,102A @@ -6178,28 +6469,62 @@ double TE = 0.0; //most recent echo time recorded case kEffectiveEchoSpacingGE: if (d.manufacturer == kMANUFACTURER_GE) d.effectiveEchoSpacingGE = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); break; - case kDiffusionBFactorGE : + case kImageTypeGE: { //0/1/2/3 for magnitude/phase/real/imaginary + if (d.manufacturer != kMANUFACTURER_GE) break; + int dt = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); + if (dt == 0) d.isHasMagnitude = true; + if (dt == 1) d.isHasPhase = true; + if (dt == 2) d.isHasReal = true; + if (dt == 3) d.isHasImaginary = true; + break; } + case kDiffusion_bValueGE : if (d.manufacturer == kMANUFACTURER_GE) { d.CSA.dtiV[0] = (float)set_bValGE(&volDiffusion, lLength, &buffer[lPos]); d.CSA.numDti = 1; } break; + case kEpiRTGroupDelayGE : //FL + if (d.manufacturer != kMANUFACTURER_GE) + break; + d.groupDelay = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); + d.groupDelay *= 1000.0; //sec -> ms + // If kEpiRTGroupDelayGE (0043,107C) exists, epiRT + d.epiVersionGE = 1; //-1 = not epi, 0 = epi, 1 = epiRT + break; + case kAssetRFactorsGE: { //DS issue427GE + if (d.manufacturer != kMANUFACTURER_GE) break; + float PhaseSlice[3]; + dcmMultiFloat(lLength, (char*)&buffer[lPos], 2, PhaseSlice); + if (PhaseSlice[1] > 0.0) + d.accelFactPE = 1.0f / PhaseSlice[1]; + if (PhaseSlice[2] > 0.0) + d.accelFactOOP = 1.0f / PhaseSlice[2]; + break; + } + case kMultiBandGE: { //LO issue427GE + if (d.manufacturer != kMANUFACTURER_GE) break; + //LO array: Value 1 = Multiband factor, Value 2 = Slice FOV Shift Factor, Value 3 = Calibration method + int mb = dcmStrInt(lLength, &buffer[lPos]); + if (mb > 1) d.CSA.multiBandFactor = mb; + break; + } case kGeiisFlag: if ((lLength > 4) && (buffer[lPos]=='G') && (buffer[lPos+1]=='E') && (buffer[lPos+2]=='I') && (buffer[lPos+3]=='I')) { //read a few digits, as bug is specific to GEIIS, while GEMS are fine printWarning("GEIIS violates the DICOM standard. Inspect results and admonish your vendor.\n"); isIconImageSequence = true; + if (sqDepthIcon < 0) sqDepthIcon = sqDepth; //geiisBug = true; //compressed thumbnails do not follow transfer syntax! GE should not re-use pulbic tags for these proprietary images http://sonca.kasshin.net/gdcm/Doc/GE_ImageThumbnails } break; case kStudyComments: { //char commentStr[kDICOMStr]; - //dcmStr (lLength, &buffer[lPos], commentStr); + //dcmStr(lLength, &buffer[lPos], commentStr); //printf(">> %s\n", commentStr); break; } case kProcedureStepDescription: - dcmStr (lLength, &buffer[lPos], d.procedureStepDescription); + dcmStr(lLength, &buffer[lPos], d.procedureStepDescription); break; case kOrientationACR : //use in emergency if kOrientation is not present! if (!isOrient) dcmMultiFloat(lLength, (char*)&buffer[lPos], 6, d.orient); @@ -6227,12 +6552,10 @@ double TE = 0.0; //most recent echo time recorded //printf("sliceV %g %g %g\n", sliceV.v[0], sliceV.v[1], sliceV.v[2]); isOrient = true; break; } - case kTemporalPosition : - if (d.manufacturer == kMANUFACTURER_GE) - break; //for GE use kRawDataRunNumber - d.rawDataRunNumber = dcmStrInt(lLength, &buffer[lPos]); + case kTemporalResolution : + temporalResolutionMS = dcmStrFloat(lLength, &buffer[lPos]); break; - case kImagesInAcquisition : + case kImagesInAcquisition : imagesInAcquisition = dcmStrInt(lLength, &buffer[lPos]); break; //case kSliceLocation : //optional so useless, infer from image position patient (0020,0032) and image orientation (0020,0037) @@ -6275,13 +6598,58 @@ double TE = 0.0; //most recent echo time recorded isIconImageSequence = false; break; } //switch/case for groupElement - + if ((((groupElement >>8) & 0xFF) == 0x60) && (groupElement % 2 == 0) && ((groupElement & 0xFF) < 0x1E)) { //Group 60xx: OverlayGroup http://dicom.nema.org/dicom/2013/output/chtml/part03/sect_C.9.html + //even group numbers 0x6000..0x601E + int overlayN = ((groupElement & 0xFF) >> 1); + //printf("%08x %d %d\n", groupElement, (groupElement & 0xFF), overlayN); + int element = groupElement>>16; + switch(element) { + case 0x0010: //US OverlayRows + overlayRows = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); + break; + case 0x0011: //US OverlayColumns + overlayCols = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); + break; + case 0x0050: {//SSx2! OverlayOrigin + if (lLength != 4) break; + int row = dcmInt(2,&buffer[lPos],d.isLittleEndian); + int col = dcmInt(2,&buffer[lPos+2],d.isLittleEndian); + if ((row == 1) && (col == 1)) break; + printMessage("Unsupported overlay origin %d/%d\n", row, col); + overlayOK = false; + break; + } + case 0x0100: {//US OverlayBitsAllocated + int bits = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); + if (bits == 1) break; + //old style Burned-In + printMessage("Illegal/Obsolete DICOM: Overlay Bits Allocated must be 1, not %d\n", bits); + overlayOK = false; + break; + } + case 0x0102: {//US OverlayBitPosition + int pos = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); + if (pos == 0) break; + //old style Burned-In + printMessage("Illegal/Obsolete DICOM: Overlay Bit Position shall be 0, not %d\n", pos); + overlayOK = false; + break; + } + case 0x3000: { + d.overlayStart[overlayN] = (int)lPos + (int)lFileOffset; + d.isHasOverlay = true; + break; + } + } + }//Group 60xx even values 0x6000..0x601E https://www.medicalconnections.co.uk/kb/Number-Of-Overlays-In-Image/ + #ifndef USING_R if (isVerbose > 1) { //dcm2niix i fast because it does not use a dictionary. // this is a very incomplete DICOM header report, and not a substitute for tools like dcmdump // the purpose is to see how dcm2niix has parsed the image for diagnostics // this section will report very little for implicit data + //if (d.isHasReal) printf("r");else printf("m"); char str[kDICOMStr]; sprintf(str, "%*c%04x,%04x %u@%ld ", sqDepth+1, ' ', groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos); bool isStr = false; @@ -6323,19 +6691,20 @@ double TE = 0.0; //most recent echo time recorded //tagStr[0] = 'X'; //avoid compiler warning: orientStr filled by dcmStr strcpy(tagStr,""); if (lLength > 0) - dcmStr (lLength, &buffer[lPos], tagStr); + dcmStr(lLength, &buffer[lPos], tagStr); if (strlen(tagStr) > 1) { for (size_t pos = 0; pos') || (tagStr[pos] == ':') || (tagStr[pos] == '"') || (tagStr[pos] == '\\') || (tagStr[pos] == '/') - || (tagStr[pos] == '^') || (tagStr[pos] < 33) + || (tagStr[pos] < 32) //issue398 + //|| (tagStr[pos] == '^') || (tagStr[pos] < 33) || (tagStr[pos] == '*') || (tagStr[pos] == '|') || (tagStr[pos] == '?')) - tagStr[pos] = 'x'; + tagStr[pos] = '_'; } printMessage("%s %s\n", str, tagStr); } else printMessage("%s\n", str); - //if (d.isExplicitVR) printMessage(" VR=%c%c\n", vr[0], vr[1]); + //if (d.isExplicitVR) printMessage(" VR=%c%c\n", vr[0], vr[1]); } //printMessage(" tag=%04x,%04x length=%u pos=%ld %c%c nest=%d\n", groupElement & 65535,groupElement>>16, lLength, lPos,vr[0], vr[1], nest); #endif lPos = lPos + (lLength); @@ -6360,10 +6729,18 @@ double TE = 0.0; //most recent echo time recorded printWarning("DICOM violation (contact vendor): compressed image without image fragments, assuming image offset defined by 0x7FE0,x0010: %s\n", fname); d.imageStart = encapsulatedDataImageStart; } + if ((d.manufacturer == kMANUFACTURER_GE) && (d.groupDelay > 0.0)) + d.TR += d.groupDelay; //Strangely, for GE the sample rate is (0018,0080) + ((0043,107c) * 1000.0) if ((d.modality == kMODALITY_PT) && (PETImageIndex > 0)) { d.imageNum = PETImageIndex; //https://github.com/rordenlab/dcm2niix/issues/184 //printWarning("PET scan using 0054,1330 for image number %d\n", PETImageIndex); } + if (d.isHasOverlay) { + if ((overlayCols > 0) && (d.xyzDim[1] != overlayCols)) overlayOK = false; + if ((overlayRows > 0) && (d.xyzDim[2] != overlayRows)) overlayOK = false; + if (!overlayOK) + d.isHasOverlay = false; + } //Recent Philips images include DateTime (0008,002A) but not separate date and time (0008,0022 and 0008,0032) #define kYYYYMMDDlen 8 //how many characters to encode year,month,day in "YYYYDDMM" format if ((strlen(acquisitionDateTimeTxt) > (kYYYYMMDDlen+5)) && (!isFloatDiff(d.acquisitionTime, 0.0f)) && (!isFloatDiff(d.acquisitionDate, 0.0f)) ) { @@ -6391,7 +6768,8 @@ double TE = 0.0; //most recent echo time recorded /* SAH.start: Fix for ZIP2 */ int zipFactor = (int) roundf(d.xyzMM[3] / d.zSpacing); if (zipFactor > 1) { - printMessage("Issue 373: Check for ZIP2 Factor: %d SliceThickness+SliceGap: %f, SpacingBetweenSlices: %f \n", zipFactor, d.xyzMM[3], d.zSpacing); + d.interp3D = zipFactor; + //printMessage("Issue 373: Check for ZIP2 Factor: %d SliceThickness+SliceGap: %f, SpacingBetweenSlices: %f \n", zipFactor, d.xyzMM[3], d.zSpacing); locationsInAcquisitionGE *= zipFactor; // Multiply number of slices by ZIP factor. Do this prior to checking for conflict below (?). } /* SAH.end */ @@ -6444,7 +6822,6 @@ double TE = 0.0; //most recent echo time recorded printError("Unable to decode %d-bit images with Transfer Syntax 1.2.840.10008.1.2.4.51, decompress with dcmdjpg or gdcmconv\n", d.bitsAllocated); d.isValid = false; } - if ((isMosaic) && (d.CSA.mosaicSlices < 1) && (numberOfImagesInMosaic < 1) && (!isInterpolated) && (d.phaseEncodingLines > 0) && (frequencyRows > 0) && ((d.xyzDim[1] % frequencyRows) == 0) && ((d.xyzDim[1] / frequencyRows) > 2) && ((d.xyzDim[2] % d.phaseEncodingLines) == 0) && ((d.xyzDim[2] / d.phaseEncodingLines) > 2) ) { //n.b. in future check if frequency is in row or column direction (and same with phase) // >2 avoids detecting interpolated as mosaic, in future perhaps check "isInterpolated" @@ -6464,7 +6841,7 @@ double TE = 0.0; //most recent echo time recorded printWarning("0008,0008=MOSAIC but number of slices not specified: %s\n", fname); if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.CSA.dtiV[1] < -1.0) && (d.CSA.dtiV[2] < -1.0) && (d.CSA.dtiV[3] < -1.0)) d.CSA.dtiV[0] = 0; //SiemensTrio-Syngo2004A reports B=0 images as having impossible b-vectors. - if ((d.manufacturer == kMANUFACTURER_GE) && (strlen(d.seriesDescription) > 1)) //GE uses a generic session name here: do not overwrite kProtocolNameGE + if ((d.manufacturer == kMANUFACTURER_GE) && (d.modality == kMODALITY_MR) && (strlen(d.seriesDescription) > 1)) //GE uses a generic session name here: do not overwrite kProtocolNameGE strcpy(d.protocolName, d.seriesDescription); if ((strlen(d.protocolName) < 1) && (strlen(d.seriesDescription) > 1)) strcpy(d.protocolName, d.seriesDescription); @@ -6517,6 +6894,10 @@ if (d.isHasPhase) d.patientPositionLast[k] = patientPositionEndPhilips[k]; } } + if ((B0Philips >= 0) && (d.CSA.numDti == 0)) { + d.CSA.dtiV[0] = B0Philips; + d.CSA.numDti = 1; + } //issue409 Siemens XA saved as classic 2D not enhanced if (!isnan(patientPositionStartPhilips[1])) //for Philips data without for (int k = 0; k < 4; k++) d.patientPosition[k] = patientPositionStartPhilips[k]; @@ -6633,7 +7014,7 @@ if (d.isHasPhase) dti4D->S[i].V[2] = dcmDim[j+(i * d.xyzDim[3])].V[2]; dti4D->S[i].V[3] = dcmDim[j+(i * d.xyzDim[3])].V[3]; //printf("te=\t%g\tscl=\t%g\tintercept=\t%g\n",dti4D->TE[i], dti4D->intenScale[i],dti4D->intenIntercept[i]); - if (dti4D->TE[i] != d.TE) isTEvaries = true; + if ((!isSameFloatGE(dti4D->TE[i],0.0)) && (dti4D->TE[i] != d.TE)) isTEvaries = true; if (dti4D->isPhase[i] != isPhase) d.isScaleOrTEVaries = true; if (dti4D->triggerDelayTime[i] != d.triggerDelayTime) d.isScaleOrTEVaries = true; if (dti4D->isReal[i] != isReal) d.isScaleOrTEVaries = true; @@ -6734,6 +7115,18 @@ if (d.isHasPhase) strcpy(d.seriesInstanceUID, d.studyInstanceUID); d.seriesUidCrc = mz_crc32X((unsigned char*) &d.protocolName, strlen(d.protocolName)); } + if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (!isTriggerSynced)) //issue408 + d.triggerDelayTime = 0.0; + if (isSameFloat(MRImageDynamicScanBeginTime * 1000.0, d.triggerDelayTime)) //issue395 + d.triggerDelayTime = 0.0; + //printf("%d\t%g\t%g\t%g\n", d.imageNum, d.acquisitionTime, d.triggerDelayTime, MRImageDynamicScanBeginTime); + if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (strlen(seriesTimeTxt) > 1) && (d.isXA10A) && (d.xyzDim[3] == 1) && (d.xyzDim[4] < 2)) { + //printWarning("Ignoring series number of XA data saved as classic DICOM (issue 394)\n"); + d.isStackableSeries = true; + d.imageNum += (d.seriesNum * 1000); + strcpy(d.seriesInstanceUID, seriesTimeTxt); // dest <- src + d.seriesUidCrc = mz_crc32X((unsigned char*) &seriesTimeTxt, strlen(seriesTimeTxt)); + } if (((d.manufacturer == kMANUFACTURER_TOSHIBA) || (d.manufacturer == kMANUFACTURER_CANON))&& (B0Philips > 0.0)) {//issue 388 char txt[1024] = {""}; sprintf(txt, "b=%d(", (int) round(B0Philips)); @@ -6749,23 +7142,34 @@ if (d.isHasPhase) } float v[4]; dcmMultiFloat(len,(char*)&txt[0], 3, &v[0]); - //printf(">>>%g = %g %g %g : %s %s\n", v[0], v[1], v[2], v[3], txt,fname); - d.CSA.dtiV[0] = B0Philips; - d.CSA.dtiV[1] = v[1]; - d.CSA.dtiV[2] = v[2]; - d.CSA.dtiV[3] = v[3]; + d.CSA.dtiV[0] = B0Philips; + #ifdef swizzleCanon //see issue422 and dcm_qa_canon + d.CSA.dtiV[1] = v[2]; + d.CSA.dtiV[2] = v[1]; + d.CSA.dtiV[3] = -v[3]; + #else + d.CSA.dtiV[1] = v[2]; + d.CSA.dtiV[2] = v[1]; + d.CSA.dtiV[3] = v[3]; + d.manufacturer = kMANUFACTURER_CANON; + #endif + //d.CSA.dtiV[1] = v[1]; + //d.CSA.dtiV[2] = v[2]; + //d.CSA.dtiV[3] = v[3]; d.CSA.numDti = 1; } } if ((isDICOMANON) && (isMATLAB)) { //issue 383 - d.seriesInstanceUID[0] = '\0'; - strcat(d.seriesInstanceUID, d.studyDate); - strcat(d.seriesInstanceUID, d.studyTime); + strcpy(d.seriesInstanceUID, d.studyDate); + // This check is unlikely to be important in practice, but it silences a warning from GCC with -Wrestrict + if (strlen(d.studyDate) < kDICOMStr) { + strncat(d.seriesInstanceUID, d.studyTime, kDICOMStr-strlen(d.studyDate)); + } d.seriesUidCrc = mz_crc32X((unsigned char*) &d.seriesInstanceUID, strlen(d.seriesInstanceUID)); } - if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (strstr(d.sequenceName, "_fl2d1") != NULL)) { + if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (strstr(d.sequenceName, "fl2d1") != NULL)) { d.isLocalizer = true; } //UIH 3D T1 scans report echo train length, which is interpreted as 3D EPI @@ -6798,6 +7202,11 @@ if (d.isHasPhase) #ifndef myLoadWholeFileToReadHeader fclose(file); #endif + if ((temporalResolutionMS > 0.0) && (isSameFloatGE(d.TR,temporalResolutionMS)) ) { + //do something profound + //in practice 0020,0110 not used + //https://github.com/bids-standard/bep001/blob/repetitiontime/Proposal_RepetitionTime.md + } if (hasDwiDirectionality) d.isVectorFromBMatrix = false; //issue 265: Philips/Siemens have both directionality and bmatrix, Bruker only has bmatrix /* fixed 2/2019 by modifying to kDiffusionBFactor, kDiffusionDirectionRL, kDiffusionDirectionAP, kDiffusionDirectionFH @@ -6813,6 +7222,8 @@ if (d.isHasPhase) */ //printf("%s\t%s\t%s\t%s\t%s_%s\n",d.patientBirthDate, d.procedureStepDescription,d.patientName, fname, d.studyDate, d.studyTime); //d.isValid = false; + //printMessage(" patient position (0020,0032)\t%g\t%g\t%g\n", d.patientPosition[1],d.patientPosition[2],d.patientPosition[3]); + //printf("%d\t%g\t%g\t%g\t%g\n", d.imageNum, d.rtia_timerGE, d.patientPosition[1],d.patientPosition[2],d.patientPosition[3]); //printf("%g\t\t%g\t%g\t%g\t%s\n", d.CSA.dtiV[0], d.CSA.dtiV[1], d.CSA.dtiV[2], d.CSA.dtiV[3], fname); //printMessage("buffer usage %d %d %d\n",d.imageStart, lPos+lFileOffset, MaxBufferSz); return d; diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 37d57167..3fa01c4a 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -42,9 +42,16 @@ extern "C" { #else #define kCCsuf " CompilerNA" //unknown compiler! #endif +#if defined(__arm__) || defined(__ARM_ARCH) + #define kCPUsuf " ARM" +#elif defined(__x86_64) + #define kCPUsuf " x86-64" +#else + #define kCPUsuf " " //unknown CPU +#endif -#define kDCMdate "v1.0.20200331" -#define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf +#define kDCMdate "v1.0.20201102" +#define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic static const int kMaxDTI4D = 18000; //maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images @@ -90,8 +97,8 @@ static const int kMaxSlice2D = 64000; //maximum number of 2D slices in 4D (Phili #define kEXIT_INPUT_FOLDER_INVALID 5 #define kEXIT_OUTPUT_FOLDER_INVALID 6 #define kEXIT_OUTPUT_FOLDER_READ_ONLY 7 -#define kEXIT_RENAME_ERROR 8 #define kEXIT_SOME_OK_SOME_BAD 8 +#define kEXIT_RENAME_ERROR 9 static const int kSliceOrientUnknown = 0; @@ -106,6 +113,7 @@ static const int kCompress50 = 3; //obsolete JPEG lossy static const int kCompressRLE = 4; //run length encoding static const int kCompressPMSCT_RLE1 = 5; //see rle2img: Philips/ELSCINT1 run-length compression 07a1,1011= PMSCT_RLE1 static const int kCompressJPEGLS = 6; //LoCo JPEG-LS +static const int kMaxOverlay = 16; //even group values 0x6000..0x601E #ifdef myEnableJasper static const int kCompressSupport = kCompressYes; //JASPER for JPEG2000 #else @@ -127,10 +135,11 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; struct TDTI S[kMaxDTI4D]; int sliceOrder[kMaxSlice2D]; // [7,3,2] means the first slice on disk should be moved to 7th position int gradDynVol[kMaxDTI4D]; //used to parse dimensions of Philips data, e.g. file with multiple dynamics, echoes, phase+magnitude - float triggerDelayTime[kMaxDTI4D], TE[kMaxDTI4D], RWVScale[kMaxDTI4D], RWVIntercept[kMaxDTI4D], intenScale[kMaxDTI4D], intenIntercept[kMaxDTI4D], intenScalePhilips[kMaxDTI4D]; + float frameDuration[kMaxDTI4D], decayFactor[kMaxDTI4D], volumeOnsetTime[kMaxDTI4D], triggerDelayTime[kMaxDTI4D], TE[kMaxDTI4D], RWVScale[kMaxDTI4D], RWVIntercept[kMaxDTI4D], intenScale[kMaxDTI4D], intenIntercept[kMaxDTI4D], intenScalePhilips[kMaxDTI4D]; bool isReal[kMaxDTI4D]; bool isImaginary[kMaxDTI4D]; bool isPhase[kMaxDTI4D]; + float repetitionTimeExcitation, repetitionTimeInversion; }; #ifdef _MSC_VER //Microsoft nomenclature for packed structures is different... @@ -170,18 +179,20 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; long seriesNum; int xyzDim[5]; uint32_t coilCrc, seriesUidCrc, instanceUidCrc; - int maxEchoNumGE, rawDataRunNumber, numberOfImagesInGridUIH, numberOfDiffusionDirectionGE, phaseEncodingGE, protocolBlockStartGE, protocolBlockLengthGE, modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, echoTrainLength, echoNum, sliceOrient, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,locationsInAcquisition, locationsInAcquisitionConflict, compressionScheme; - float percentSampling,waterFatShift, numberOfAverages, imagingFrequency, patientWeight, zSpacing, zThick, pixelBandwidth, SAR, phaseFieldofView, accelFactPE, flipAngle, fieldStrength, TE, TI, TR, intenScale, intenIntercept, intenScalePhilips, gantryTilt, lastScanLoc, angulation[4]; + int overlayStart[kMaxOverlay]; + int interp3D, epiVersionGE, internalepiVersionGE, maxEchoNumGE, rawDataRunNumber, numberOfImagesInGridUIH, numberOfDiffusionDirectionGE, phaseEncodingGE, protocolBlockStartGE, protocolBlockLengthGE, modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, echoTrainLength, echoNum, sliceOrient, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,locationsInAcquisition, locationsInAcquisitionConflict, compressionScheme; + float groupDelay, decayFactor, percentSampling,waterFatShift, numberOfAverages, imagingFrequency, patientWeight, zSpacing, zThick, pixelBandwidth, SAR, phaseFieldofView, accelFactPE, accelFactOOP, flipAngle, fieldStrength, TE, TI, TR, intenScale, intenIntercept, intenScalePhilips, gantryTilt, lastScanLoc, angulation[4]; float orient[7], patientPosition[4], patientPositionLast[4], xyzMM[4], stackOffcentre[4]; float rtia_timerGE, radionuclidePositronFraction, radionuclideTotalDose, radionuclideHalfLife, doseCalibrationFactor; //PET ISOTOPE MODULE ATTRIBUTES (C.8-57) - float ecat_isotope_halflife, ecat_dosage; + float frameDuration, ecat_isotope_halflife, ecat_dosage; float pixelPaddingValue; // used for both FloatPixelPaddingValue (0028, 0122) and PixelPaddingValue (0028, 0120); NaN if not present. double acquisitionDuration, triggerDelayTime, RWVScale, RWVIntercept, dateTime, acquisitionTime, acquisitionDate, bandwidthPerPixelPhaseEncode; + char radiopharmaceutical[kDICOMStr], convolutionKernel[kDICOMStr], unitsPT[kDICOMStr], decayCorrection[kDICOMStr], attenuationCorrectionMethod[kDICOMStr],reconstructionMethod[kDICOMStr]; char coilElements[kDICOMStr], coilName[kDICOMStr], phaseEncodingDirectionDisplayedUIH[kDICOMStr], imageBaseName[kDICOMStr], scanOptions[kDICOMStr], stationName[kDICOMStr], softwareVersions[kDICOMStr], deviceSerialNumber[kDICOMStr], institutionName[kDICOMStr], referringPhysicianName[kDICOMStr], instanceUID[kDICOMStr], seriesInstanceUID[kDICOMStr], studyInstanceUID[kDICOMStr], bodyPartExamined[kDICOMStr], procedureStepDescription[kDICOMStr], imageType[kDICOMStr], institutionalDepartmentName[kDICOMStr], manufacturersModelName[kDICOMStr], patientID[kDICOMStr], patientOrient[kDICOMStr], patientName[kDICOMStr], accessionNumber[kDICOMStr], seriesDescription[kDICOMStr], studyID[kDICOMStr], sequenceName[kDICOMStr], protocolName[kDICOMStr],sequenceVariant[kDICOMStr],scanningSequence[kDICOMStr], patientBirthDate[kDICOMStr], patientAge[kDICOMStr], studyDate[kDICOMStr],studyTime[kDICOMStr]; char institutionAddress[kDICOMStrLarge], imageComments[kDICOMStrLarge]; uint32_t dimensionIndexValues[MAX_NUMBER_OF_DIMENSIONS]; struct TCSAdata CSA; - bool isPartialFourier, isDiffusion, isVectorFromBMatrix, isRawDataStorage, isGrayscaleSoftcopyPresentationState, isStackableSeries, isCoilVaries, isNonParallelSlices, isSegamiOasis, isXA10A, isScaleOrTEVaries, isScaleVariesEnh, isDerived, isXRay, isMultiEcho, isValid, is3DAcq, is2DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase, isHasImaginary, isHasReal, isHasMagnitude,isHasMixed, isFloat, isResampled, isLocalizer; + bool isPrivateCreatorRemap, isHasOverlay, isEPI, isIR, isPartialFourier, isDiffusion, isVectorFromBMatrix, isRawDataStorage, isGrayscaleSoftcopyPresentationState, isStackableSeries, isCoilVaries, isNonParallelSlices, isSegamiOasis, isXA10A, isScaleOrTEVaries, isScaleVariesEnh, isDerived, isXRay, isMultiEcho, isValid, is3DAcq, is2DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase, isHasImaginary, isHasReal, isHasMagnitude,isHasMixed, isFloat, isResampled, isLocalizer; char phaseEncodingRC, patientSex; }; diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 0dee754e..0b0a6a0f 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -189,7 +189,7 @@ void geCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI *vx, int isV //phase encoding. This is in DICOM message 0018,1312. If this has value //COL then if swap the x and y value and reverse the sign on the z value. //If the phase encoding is not COL, then just reverse the sign on the x value. - if (d->manufacturer != kMANUFACTURER_GE) return; + if ((d->manufacturer != kMANUFACTURER_GE) && (d->manufacturer != kMANUFACTURER_CANON)) return; if (d->CSA.numDti < 1) return; if ((toupper(d->patientOrient[0])== 'H') && (toupper(d->patientOrient[1])== 'F') && (toupper(d->patientOrient[2])== 'S')) ; //participant was head first supine @@ -236,16 +236,30 @@ void geCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI *vx, int isV continue; //do not normalize or reorient 0 vectors } if ((vLen > 0.03) && (vLen < 0.97)) { - //bVal scaled by norm(g)^2 https://github.com/rordenlab/dcm2niix/issues/163 - float bVal = vx[i].V[0] * (vLen * vLen); + //bVal scaled by norm(g)^2 issue163,245 + float bValtemp = 0, bVal = 0, bVecScale=0; + // rounding by 5 with mimimum of 5 if b-value > 0 + bValtemp = vx[i].V[0] * (vLen * vLen); + if (bValtemp > 0 && bValtemp < 5) { + bVal = 5; + } + else { + bVal = (int)((bValtemp + 2.5f)/5)*5; + } + if(bVal == 0) + bVecScale = 0; + else + { + bVecScale = sqrt((float)vx[i].V[0]/bVal); + } if (!scaledBValWarning) { printMessage("GE BVal scaling (e.g. %g -> %g s/mm^2)\n", vx[i].V[0], bVal); scaledBValWarning = true; } vx[i].V[0] = bVal; - vx[i].V[1] = vx[i].V[1]/vLen; - vx[i].V[2] = vx[i].V[2]/vLen; - vx[i].V[3] = vx[i].V[3]/vLen; + vx[i].V[1] = vx[i].V[1]*bVecScale; + vx[i].V[2] = vx[i].V[2]*bVecScale; + vx[i].V[3] = vx[i].V[3]*bVecScale; } if (!col) { //rows need to be swizzled //see Stanford dataset Ax_DWI_Tetrahedral_7 unable to resolve between possible solutions @@ -278,7 +292,7 @@ void siemensPhilipsCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI //convert DTI vectors from scanner coordinates to image frame of reference //Uses 6 orient values from ImageOrientationPatient (0020,0037) // requires PatientPosition 0018,5100 is HFS (head first supine) - if ((d->manufacturer != kMANUFACTURER_BRUKER) && (d->manufacturer != kMANUFACTURER_TOSHIBA) && (d->manufacturer != kMANUFACTURER_CANON) && (d->manufacturer != kMANUFACTURER_HITACHI) && (d->manufacturer != kMANUFACTURER_UIH) && (d->manufacturer != kMANUFACTURER_SIEMENS) && (d->manufacturer != kMANUFACTURER_PHILIPS)) return; + if ((d->manufacturer != kMANUFACTURER_BRUKER) && (d->manufacturer != kMANUFACTURER_TOSHIBA) && (d->manufacturer != kMANUFACTURER_HITACHI) && (d->manufacturer != kMANUFACTURER_UIH) && (d->manufacturer != kMANUFACTURER_SIEMENS) && (d->manufacturer != kMANUFACTURER_PHILIPS)) return; if (d->CSA.numDti < 1) return; if (d->manufacturer == kMANUFACTURER_UIH) { for (int i = 0; i < d->CSA.numDti; i++) { @@ -286,8 +300,8 @@ void siemensPhilipsCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI for (int v= 0; v < 4; v++) if (vx[i].V[v] == -0.0f) vx[i].V[v] = 0.0f; //remove sign from values that are virtually zero } - //for (int i = 0; i < 3; i++) - // printf("%g %g %g\n", vx[i].V[1], vx[i].V[2], vx[i].V[3]); + for (int i = 0; i < 3; i++) + printf("%g = %g %g %g\n", vx[i].V[0], vx[i].V[1], vx[i].V[2], vx[i].V[3]); return; } //https://github.com/rordenlab/dcm2niix/issues/225 if ((toupper(d->patientOrient[0])== 'H') && (toupper(d->patientOrient[1])== 'F') && (toupper(d->patientOrient[2])== 'S')) @@ -426,9 +440,9 @@ const void * memmem(const char *l, size_t l_len, const char *s, size_t s_len) { #endif //for systems without memmem int readKeyN1(const char * key, char * buffer, int remLength) { //look for text key in binary data stream, return subsequent integer value - int ret = 0; char *keyPos = (char *)memmem(buffer, remLength, key, strlen(key)); if (!keyPos) return -1; + int ret = 0; int i = (int)strlen(key); while( ( i< remLength) && (keyPos[i] != 0x0A) ) { if( keyPos[i] >= '0' && keyPos[i] <= '9' ) @@ -441,7 +455,7 @@ int readKeyN1(const char * key, char * buffer, int remLength) { //look for text int readKey(const char * key, char * buffer, int remLength) { //look for text key in binary data stream, return subsequent integer value int ret = 0; char *keyPos = (char *)memmem(buffer, remLength, key, strlen(key)); - if (!keyPos) return 0; + if (!keyPos) return ret; int i = (int)strlen(key); while( ( i< remLength) && (keyPos[i] != 0x0A) ) { if( keyPos[i] >= '0' && keyPos[i] <= '9' ) @@ -449,7 +463,7 @@ int readKey(const char * key, char * buffer, int remLength) { //look for text k i++; } return ret; -} //readKey() +} //readKey() //return 0 if key not found float readKeyFloatNan(const char * key, char * buffer, int remLength) { //look for text key in binary data stream, return subsequent integer value char *keyPos = (char *)memmem(buffer, remLength, key, strlen(key)); @@ -460,7 +474,7 @@ float readKeyFloatNan(const char * key, char * buffer, int remLength) { //look tmpstr[1] = 0; int i = (int)strlen(key); while( ( i< remLength) && (keyPos[i] != 0x0A) ) { - if( (keyPos[i] >= '0' && keyPos[i] <= '9') || (keyPos[i] <= '.') || (keyPos[i] <= '-') ) { + if( (keyPos[i] >= '0' && keyPos[i] <= '9') || (keyPos[i] == '.') || (keyPos[i] == '-') ) { tmpstr[0] = keyPos[i]; strcat (str, tmpstr); } @@ -479,7 +493,7 @@ float readKeyFloat(const char * key, char * buffer, int remLength) { //look for tmpstr[1] = 0; int i = (int)strlen(key); while( ( i< remLength) && (keyPos[i] != 0x0A) ) { - if( (keyPos[i] >= '0' && keyPos[i] <= '9') || (keyPos[i] <= '.') || (keyPos[i] <= '-') ) { + if( (keyPos[i] >= '0' && keyPos[i] <= '9') || (keyPos[i] == '.') || (keyPos[i] == '-') ) { tmpstr[0] = keyPos[i]; strcat (str, tmpstr); } @@ -544,7 +558,7 @@ int phoenixOffsetCSASeriesHeader(unsigned char *buff, int lLength) { #define kMaxWipFree 64 typedef struct { - float delayTimeInTR, phaseOversampling, phaseResolution, txRefAmp; + float TE0, TE1, delayTimeInTR, phaseOversampling, phaseResolution, txRefAmp; int phaseEncodingLines, existUcImageNumb, ucMode, baseResolution, interp, partialFourier,echoSpacing, difBipolar, parallelReductionFactorInPlane, refLinesPE; float alFree[kMaxWipFree] ; @@ -556,7 +570,9 @@ typedef struct { void siemensCsaAscii(const char * filename, TCsaAscii* csaAscii, int csaOffset, int csaLength, float* shimSetting, char* coilID, char* consistencyInfo, char* coilElements, char* pulseSequenceDetails, char* fmriExternalInfo, char* protocolName, char* wipMemBlock) { //reads ASCII portion of CSASeriesHeaderInfo and returns lEchoTrainDuration or lEchoSpacing value // returns 0 if no value found - csaAscii->delayTimeInTR = -0.001; + csaAscii->TE0 = 0.0; + csaAscii->TE1 = 0.0; + csaAscii->delayTimeInTR = -0.001; csaAscii->phaseOversampling = 0.0; csaAscii->phaseResolution = 0.0; csaAscii->txRefAmp = 0.0; @@ -659,6 +675,10 @@ void siemensCsaAscii(const char * filename, TCsaAscii* csaAscii, int csaOffset, readKeyStr(keyStrWipMemBlock, keyPos, csaLengthTrim, wipMemBlock); char keyStrPn[] = "tProtocolName"; readKeyStr(keyStrPn, keyPos, csaLengthTrim, protocolName); + char keyStrTE0[] = "alTE[0]"; + csaAscii->TE0 = readKeyFloatNan(keyStrTE0, keyPos, csaLengthTrim); + char keyStrTE1[] = "alTE[1]"; + csaAscii->TE1 = readKeyFloatNan(keyStrTE1, keyPos, csaLengthTrim); //read ALL alTI[*] values for (int k = 0; k < kMaxWipFree; k++) csaAscii->alTI[k] = NAN; @@ -762,9 +782,12 @@ void siemensCsaAscii(const char * filename, TCsaAscii* csaAscii, int csaOffset, #define myReadGeProtocolBlock #endif #ifdef myReadGeProtocolBlock -int geProtocolBlock(const char * filename, int geOffset, int geLength, int isVerbose, int* sliceOrder, int* viewOrder, int* mbAccel) { +int geProtocolBlock(const char * filename, int geOffset, int geLength, int isVerbose, int* sliceOrder, int* viewOrder, int* mbAccel, int* nSlices, float* groupDelay, char ioptGE[]) { *sliceOrder = -1; *viewOrder = 0; + *mbAccel = 0; + *nSlices = 0; + *groupDelay = 0.0; int ret = EXIT_FAILURE; if ((geOffset < 0) || (geLength < 20)) return ret; FILE * pFile = fopen ( filename, "rb" ); @@ -836,6 +859,30 @@ int geProtocolBlock(const char * filename, int geOffset, int geLength, int isV *viewOrder = readKey(keyStrVO, (char *) pUnCmp, unCmpSz); char keyStrMB[] = "MBACCEL"; *mbAccel = readKey(keyStrMB, (char *) pUnCmp, unCmpSz); + char keyStrNS[] = "NOSLC"; + *nSlices = readKey(keyStrNS, (char *) pUnCmp, unCmpSz); + char keyStrDELACQ[] = "DELACQ"; + char DELACQ[100]; + readKeyStr(keyStrDELACQ, (char *) pUnCmp, unCmpSz, DELACQ); + char keyStrGD[] = "DELACQNOAV"; + *groupDelay = readKeyFloat(keyStrGD, (char *) pUnCmp, unCmpSz); + char keyStrIOPT[] = "IOPT"; + readKeyStr(keyStrIOPT, (char *) pUnCmp, unCmpSz, ioptGE); + char PHASEDELAYS1[10000]; + char keyStrPHASEDELAYS1[] = "PHASEDELAYS1"; + readKeyStr(keyStrPHASEDELAYS1, (char *) pUnCmp, unCmpSz, PHASEDELAYS1); + if (strstr(ioptGE,"MPh") != NULL) { + if (strcmp(DELACQ, "Minimum") == 0) { + *groupDelay = 0; + } + if (strstr(ioptGE,"MPhVar") != NULL) { + *groupDelay = -1; + // Multiphase EPI with Variable Delays + // TO-DO + // NEED TO rescue ALL_PHASES case (=Group delay) + // IF values in PHASEDELAYS1 are all same except 1st value (0), this case should be same as Group Delay + } + } if (isVerbose > 1) { printMessage("GE Protocol Block %s bytes %d compressed, %d uncompressed @ %d\n", filename, geLength, unCmpSz, geOffset); printMessage(" ViewOrder %d SliceOrder %d\n", *viewOrder, *sliceOrder); @@ -846,25 +893,49 @@ int geProtocolBlock(const char * filename, int geOffset, int geLength, int isV } #endif //myReadGeProtocolBlock() -void json_Str(FILE *fp, const char *sLabel, char *sVal) { +void json_Str(FILE *fp, const char *sLabel, char *sVal) { // issue131,425 if (strlen(sVal) < 1) return; - //fprintf(fp, sLabel, sVal ); - //convert \ ' " characters to _ see https://github.com/rordenlab/dcm2niix/issues/131 - for (size_t pos = 0; pos < strlen(sVal); pos ++) { - if ((sVal[pos] == '\'') || (sVal[pos] == '"') || (sVal[pos] == '\\')) - sVal[pos] = '_'; + unsigned char sValEsc[2048] = {""}; + unsigned char *iVal = (unsigned char *) sVal; + int o = 0; + for (int i = 0; i < strlen(sVal); i ++) { + //escape double quote (") and Backslash + if ((sVal[i] == '"') || (sVal[i] == '\\')) { //escape double quotes and back slash + sValEsc[o] = '\\'; + o++; + } + if ((sVal[i] >= 0x01) && (sVal[i] <= 0x07)) continue; //control characters like "bell" + //http://dicom.nema.org/medical/dicom/current/output/html/part05.html + //0x08 Backspace is replaced with \b + //0x09 Tab is replaced with \t + //0x0A Newline is replaced with \n + //0x0B Escape ?? + //0x0C Form feed is replaced with \f + //0x0D Carriage return is replaced with \r + if ((sVal[i] >= 0x08) && (sVal[i] <= 0x0D)) { + sValEsc[o] = '\\'; + o++; + if (sVal[i] == 0x08) sValEsc[o] = 'b'; + if (sVal[i] == 0x09) sValEsc[o] = '9'; + if (sVal[i] == 0x0A) sValEsc[o] = 'n'; + if (sVal[i] == 0x0B) sValEsc[o] = '\\'; + if (sVal[i] == 0x0C) sValEsc[o] = 'f'; + if (sVal[i] == 0x0D) sValEsc[o] = 'r'; + o++; + continue; + } + //https://stackoverflow.com/questions/4059775/convert-iso-8859-1-strings-to-utf-8-in-c-c + if (iVal[i] >= 128) { + sValEsc[o]=0xc2+(iVal[i]>0xbf); + o++; + sValEsc[o]=(iVal[i]&0x3f)+0x80; + } else { + sValEsc[o] = sVal[i]; + } + o++; } - fprintf(fp, sLabel, sVal ); - /*char outname[PATH_MAX] = {""}; - char appendChar[2] = {"\\"}; - char passChar[2] = {"\\"}; - for (int pos = 0; posdecayFactor[0] >= 0.0) { //see BEP009 PET https://docs.google.com/document/d/1mqMLnxVdLwZjDd4ZiWFqjEAmOmfcModA_R535v3eQs0 + fprintf(fp, "\t\"DecayFactor\": [\n"); + for (int i = 0; i < h->dim[4]; i++) { + if (i != 0) + fprintf(fp, ",\n"); + if (dti4D->decayFactor[i] < 0) break; + fprintf(fp, "\t\t%g", dti4D->decayFactor[i] ); + } + fprintf(fp, "\t],\n"); + } + if (dti4D->volumeOnsetTime[0] >= 0.0) { //see BEP009 PET https://docs.google.com/document/d/1mqMLnxVdLwZjDd4ZiWFqjEAmOmfcModA_R535v3eQs0 + fprintf(fp, "\t\"FrameTimesStart\": [\n"); + for (int i = 0; i < h->dim[4]; i++) { + if (i != 0) + fprintf(fp, ",\n"); + if (dti4D->volumeOnsetTime[i] < 0) break; + fprintf(fp, "\t\t%g", dti4D->volumeOnsetTime[i] ); + } + fprintf(fp, "\t],\n"); + } + if (dti4D->frameDuration[0] >= 0.0) { //see BEP009 PET https://docs.google.com/document/d/1mqMLnxVdLwZjDd4ZiWFqjEAmOmfcModA_R535v3eQs0 + fprintf(fp, "\t\"FrameDuration\": [\n"); + for (int i = 0; i < h->dim[4]; i++) { + if (i != 0) + fprintf(fp, ",\n"); + if (dti4D->frameDuration[i] < 0) break; + fprintf(fp, "\t\t%g", dti4D->frameDuration[i] / 1000.0 ); // from 0018,1242 ms -> sec + } + fprintf(fp, "\t],\n"); + } + + //CT parameters if ((d.TE > 0.0) && (d.isXRay)) fprintf(fp, "\t\"XRayExposure\": %g,\n", d.TE ); //MRI parameters @@ -1128,6 +1247,8 @@ tse3d: T2*/ if ((d.TE > 0.0) && (!d.isXRay)) fprintf(fp, "\t\"EchoTime\": %g,\n", d.TE / 1000.0 ); //if ((d.TE2 > 0.0) && (!d.isXRay)) fprintf(fp, "\t\"EchoTime2\": %g,\n", d.TE2 / 1000.0 ); json_Float(fp, "\t\"RepetitionTime\": %g,\n", d.TR / 1000.0 ); + json_Float(fp, "\t\"RepetitionTimeExcitation\": %g,\n", dti4D->repetitionTimeExcitation); + json_Float(fp, "\t\"RepetitionTimeInversion\": %g,\n", dti4D->repetitionTimeInversion); json_Float(fp, "\t\"InversionTime\": %g,\n", d.TI / 1000.0 ); json_Float(fp, "\t\"FlipAngle\": %g,\n", d.flipAngle ); bool interp = false; //2D interpolation @@ -1144,7 +1265,6 @@ tse3d: T2*/ if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.CSA.SeriesHeader_offset > 0) && (d.CSA.SeriesHeader_length > 0)) { float pf = 1.0f; //partial fourier float shimSetting[8]; - char protocolName[kDICOMStrLarge], fmriExternalInfo[kDICOMStrLarge], coilID[kDICOMStrLarge], consistencyInfo[kDICOMStrLarge], coilElements[kDICOMStrLarge], pulseSequenceDetails[kDICOMStrLarge], wipMemBlock[kDICOMStrLarge]; TCsaAscii csaAscii; siemensCsaAscii(filename, &csaAscii, d.CSA.SeriesHeader_offset, d.CSA.SeriesHeader_length, shimSetting, coilID, consistencyInfo, coilElements, pulseSequenceDetails, fmriExternalInfo, protocolName, wipMemBlock); @@ -1153,6 +1273,11 @@ tse3d: T2*/ //if (d.phaseEncodingLines != csaAscii.phaseEncodingLines) //e.g. phaseOversampling // printWarning("PhaseEncodingLines reported in DICOM (%d) header does not match value CSA-ASCII (%d) %s\n", d.phaseEncodingLines, csaAscii.phaseEncodingLines, pathoutname); delayTimeInTR = csaAscii.delayTimeInTR; + if ((d.isHasPhase) && (csaAscii.TE0 > 0.0) && (csaAscii.TE1 > 0.0)) { //issue400 + //https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html + json_Float(fp, "\t\"EchoTime1\": %g,\n", csaAscii.TE0 / 1000000.0 ); + json_Float(fp, "\t\"EchoTime2\": %g,\n", csaAscii.TE1 / 1000000.0 ); + } phaseOversampling = csaAscii.phaseOversampling; if (csaAscii.existUcImageNumb > 0) { if (d.CSA.protocolSliceNumber1 < 2) { @@ -1294,7 +1419,7 @@ tse3d: T2*/ fprintf(fp, "\t\"TagPlaneSNormalDTra\": %g,\n", csaAscii.sNormalDTra); } //general properties - if (csaAscii.partialFourier > 0) { + if ((csaAscii.partialFourier > 0) && ((d.modality == kMODALITY_MR))) { //check MR, e.g. do not report for Siemens PET //https://github.com/ismrmrd/siemens_to_ismrmrd/blob/master/parameter_maps/IsmrmrdParameterMap_Siemens_EPI_FLASHREF.xsl if (csaAscii.partialFourier == 1) pf = 0.5; // 4/8 if (csaAscii.partialFourier == 2) pf = 0.625; // 5/8 @@ -1302,10 +1427,12 @@ tse3d: T2*/ if (csaAscii.partialFourier == 8) pf = 0.875; fprintf(fp, "\t\"PartialFourier\": %g,\n", pf); } - if (csaAscii.interp > 0) { + if (csaAscii.interp > 0) { //in-plane interpolation interp = true; fprintf(fp, "\t\"Interpolation2D\": %d,\n", interp); } + if (d.interp3D > 1) //through-plane interpolation e.g. GE ZIP2 through-plane http://mriquestions.com/zip.html + fprintf(fp, "\t\"Interpolation3D\": %d,\n", d.interp3D); if (csaAscii.baseResolution > 0) fprintf(fp, "\t\"BaseResolution\": %d,\n", csaAscii.baseResolution ); if (shimSetting[0] != 0.0) { fprintf(fp, "\t\"ShimSetting\": [\n"); @@ -1324,15 +1451,15 @@ tse3d: T2*/ // https://groups.google.com/forum/#!topic/bids-discussion/nmg1BOVH1SU // https://groups.google.com/forum/#!topic/bids-discussion/seD7AtJfaFE json_Float(fp, "\t\"DelayTime\": %g,\n", delayTimeInTR/ 1000000.0); //DelayTimeInTR usec -> sec - json_Float(fp, "\t\"TxRefAmp\": %g,\n", csaAscii.txRefAmp); - json_Float(fp, "\t\"PhaseResolution\": %g,\n", csaAscii.phaseResolution); + if (d.modality == kMODALITY_MR) json_Float(fp, "\t\"TxRefAmp\": %g,\n", csaAscii.txRefAmp); + if (d.modality == kMODALITY_MR) json_Float(fp, "\t\"PhaseResolution\": %g,\n", csaAscii.phaseResolution); json_Float(fp, "\t\"PhaseOversampling\": %g,\n", phaseOversampling); json_Float(fp, "\t\"VendorReportedEchoSpacing\": %g,\n", csaAscii.echoSpacing / 1000000.0); //usec -> sec //ETD and epiFactor not useful/reliable https://github.com/rordenlab/dcm2niix/issues/127 //if (echoTrainDuration > 0) fprintf(fp, "\t\"EchoTrainDuration\": %g,\n", echoTrainDuration / 1000000.0); //usec -> sec //if (epiFactor > 0) fprintf(fp, "\t\"EPIFactor\": %d,\n", epiFactor); json_Str(fp, "\t\"ReceiveCoilName\": \"%s\",\n", coilID); - json_Str(fp, "\t\"ReceiveCoilActiveElements\": \"%s\",\n", coilElements); + if (d.modality == kMODALITY_MR) json_Str(fp, "\t\"ReceiveCoilActiveElements\": \"%s\",\n", coilElements); if (strcmp(coilElements,d.coilName) != 0) json_Str(fp, "\t\"CoilString\": \"%s\",\n", d.coilName); strcpy(d.coilName, ""); @@ -1361,6 +1488,7 @@ tse3d: T2*/ float pf = (float)d.phaseEncodingLines; if (d.accelFactPE > 1) pf = (float)pf / (float)d.accelFactPE; //estimate: not sure if we round up or down + pf = (float)d.echoTrainLength / (float)pf; if (pf < 1.0) //e.g. if difference between lines and echo length not all explained by iPAT (SENSE/GRAPPA) fprintf(fp, "\t\"PartialFourier\": %g,\n", pf); @@ -1368,7 +1496,7 @@ tse3d: T2*/ //printf("PhaseLines=%d EchoTrainLength=%d SENSE=%g\n", d.phaseEncodingLines, d.echoTrainLength, d.accelFactPE); //n.b. we can not distinguish pF from SENSE/GRAPPA for UIH } #endif - if (d.CSA.multiBandFactor > 1) //AccelFactorSlice + if ((d.CSA.multiBandFactor > 1) && (d.modality == kMODALITY_MR)) //AccelFactorSlice fprintf(fp, "\t\"MultibandAccelerationFactor\": %d,\n", d.CSA.multiBandFactor); json_Float(fp, "\t\"PercentPhaseFOV\": %g,\n", d.phaseFieldofView); json_Float(fp, "\t\"PercentSampling\": %g,\n", d.percentSampling); @@ -1380,7 +1508,7 @@ tse3d: T2*/ fprintf(fp, "\t\"PhaseEncodingStepsNoPartialFourier\": %d,\n", d.phaseEncodingSteps ); } else if (d.phaseEncodingSteps > 0) fprintf(fp, "\t\"PhaseEncodingSteps\": %d,\n", d.phaseEncodingSteps ); - if (d.phaseEncodingLines > 0) fprintf(fp, "\t\"AcquisitionMatrixPE\": %d,\n", d.phaseEncodingLines ); + if ((d.phaseEncodingLines > 0) && (d.modality == kMODALITY_MR)) fprintf(fp, "\t\"AcquisitionMatrixPE\": %d,\n", d.phaseEncodingLines ); //Compute ReconMatrixPE // Actual size of the *reconstructed* data in the PE dimension, which does NOT match @@ -1388,7 +1516,6 @@ tse3d: T2*/ // We'll need this for generating a value for effectiveEchoSpacing that is consistent // with the *reconstructed* data. int reconMatrixPE = d.phaseEncodingLines; - if ((h->dim[2] > 0) && (h->dim[1] > 0)) { if (h->dim[1] == h->dim[2]) //phase encoding does not matter reconMatrixPE = h->dim[2]; @@ -1397,13 +1524,14 @@ tse3d: T2*/ else if (d.phaseEncodingRC =='R') reconMatrixPE = h->dim[1]; } - if (reconMatrixPE > 0) fprintf(fp, "\t\"ReconMatrixPE\": %d,\n", reconMatrixPE ); + if ((d.modality == kMODALITY_MR) && (reconMatrixPE > 0)) fprintf(fp, "\t\"ReconMatrixPE\": %d,\n", reconMatrixPE ); double bandwidthPerPixelPhaseEncode = d.bandwidthPerPixelPhaseEncode; if (bandwidthPerPixelPhaseEncode == 0.0) bandwidthPerPixelPhaseEncode = d.CSA.bandwidthPerPixelPhaseEncode; json_Float(fp, "\t\"BandwidthPerPixelPhaseEncode\": %g,\n", bandwidthPerPixelPhaseEncode ); //if ((!d.is3DAcq) && (d.accelFactPE > 1.0)) fprintf(fp, "\t\"ParallelReductionFactorInPlane\": %g,\n", d.accelFactPE); if (d.accelFactPE > 1.0) fprintf(fp, "\t\"ParallelReductionFactorInPlane\": %g,\n", d.accelFactPE); //https://github.com/rordenlab/dcm2niix/issues/314 + if (d.accelFactOOP > 1.0) fprintf(fp, "\t\"ParallelReductionOutOfPlane\": %g,\n", d.accelFactOOP); //EffectiveEchoSpacing // Siemens bandwidthPerPixelPhaseEncode already accounts for the effects of parallel imaging, // interpolation, phaseOversampling, and phaseResolution, in the context of the size of the @@ -1471,6 +1599,7 @@ tse3d: T2*/ //next two conditionals updated: make GE match Siemens if (d.phaseEncodingGE == kGE_PHASE_ENCODING_POLARITY_UNFLIPPED) phPos = 1; if (d.phaseEncodingGE == kGE_PHASE_ENCODING_POLARITY_FLIPPED) phPos = 0; + //if ((phPos >= 0) && (d.phaseEncodingRC == 'R') && (d.manufacturer == kMANUFACTURER_UIH)) phPos = 1 - phPos; //issue410 bool isSkipPhaseEncodingAxis = d.is3DAcq; if (d.echoTrainLength > 1) isSkipPhaseEncodingAxis = false; //issue 371: ignore phaseEncoding for 3D MP-RAGE/SPACE, but report for 3D EPI if (((d.phaseEncodingRC == 'R') || (d.phaseEncodingRC == 'C')) && (!isSkipPhaseEncodingAxis) && (phPos < 0)) { @@ -1506,7 +1635,7 @@ tse3d: T2*/ } //only save PhaseEncodingDirection if BOTH direction and POLARITY are known //Slice Timing UIH or GE >>>> //in theory, we should also report XA10 slice times here, but see series 24 of https://github.com/rordenlab/dcm2niix/issues/236 - if ((!d.is3DAcq) && (d.CSA.sliceTiming[0] >= 0.0)) { + if ((d.modality != kMODALITY_PT) && (!d.is3DAcq) && (d.CSA.sliceTiming[0] >= 0.0)) { fprintf(fp, "\t\"SliceTiming\": [\n"); for (int i = 0; i < h->dim[3]; i++) { if (i != 0) @@ -1531,7 +1660,17 @@ tse3d: T2*/ //fprintf(fp, "\t\"ConversionSoftwareVersion\": \"%s\"\n", kDCMvers );kDCMdate fprintf(fp, "}\n"); fclose(fp); -}// nii_SaveBIDS() +}// nii_SaveBIDSX() + +void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, struct nifti_1_header *h, const char * filename) { +struct TDTI4D *dti4D; +dti4D->sliceOrder[0] = -1; +dti4D->volumeOnsetTime[0] = -1; +dti4D->decayFactor[0] = -1; +dti4D->intenScale[0] = 0.0; +dti4D->repetitionTimeExcitation = 0.0; +nii_SaveBIDSX(pathoutname, d, opts, h, filename, dti4D); +}// nii_SaveBIDSX() bool isADCnotDTI(TDTI bvec) { //returns true if bval!=0 but all bvecs == 0 (Philips code for derived ADC image) return ((!isSameFloat(bvec.V[0],0.0f)) && //not a B-0 image @@ -1620,17 +1759,23 @@ int * nii_saveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],str uint64_t indx0 = dcmSort[0].indx; //first volume int numDti = dcmList[indx0].CSA.numDti; #ifdef USING_R - //R might want to pad zeros -#else + ImageList *images = (ImageList *) opts.imageList; +#endif //https://github.com/rordenlab/dcm2niix/issues/352 bool allB0 = dcmList[indx0].isDiffusion; if (dcmList[indx0].isDerived) allB0 = false; //e.g. FA map if ((numDti == numVol) && (numDti > 1)) allB0 = false; if (numDti > 1) allB0 = false; + if (nConvert > 1) allB0 = false; if ((numDti == 1) && (dti4D->S[0].V[0] > 50.0)) allB0 = false; - if (allB0) { - if (opts.isVerbose) + if (allB0) { + if (opts.isVerbose) printMessage("Diffusion image without gradients: assuming %d volume B=0 series\n", numVol); +#ifdef USING_R + // The image hasn't been created yet, so the attributes must be deferred + images->addDeferredAttribute("bValues", std::vector(numVol,0.0)); + images->addDeferredAttribute("bVectors", std::vector(numVol*3,0.0), numVol, 3); +#else char sep = '\t'; if (opts.isCreateBIDS) sep = ' '; @@ -1653,8 +1798,8 @@ int * nii_saveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],str fprintf(fp, "\n"); } fclose(fp); +#endif } -#endif if (numDti < 1) return NULL; if ((numDti < 3) && (nConvert < 3)) return NULL; TDTI * vx = NULL; @@ -1684,7 +1829,10 @@ int * nii_saveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],str for (int i = 1; i < numDti; i++) //check if all bvalues match first volume if (vx[i].V[0] < minBval) minBval = vx[i].V[0]; if (minBval > 50.0) bValueVaries = true; - //do not save files without variability + //start issue394: experimental, single volume per series, Siemens XA + if (!isAllZeroFloat(vx[0].V[1], vx[0].V[2], vx[0].V[3])) + bValueVaries = true; + //end issue394 if (!bValueVaries) { bool bVecVaries = false; for (int i = 1; i < numDti; i++) {//check if all bvalues match first volume @@ -1763,13 +1911,16 @@ int * nii_saveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],str } if (numGEwarn > 0) printWarning("Some images had bval>0 but bvec=0 (either Trace or b=0, see issue 245)\n"); - if ((*numADC == numDti) || (numGEwarn == numDti)) { + /*if ((*numADC == numDti) || (numGEwarn == numDti)) { //issue 405: we now save bvals file for isotropic series //all isotropic/ADC images - no valid bvecs *numADC = 0; free(bvals); free(vx); return NULL; - } + }*/ + bool isIsotropic = false; + if ((*numADC == numDti) || (numGEwarn == numDti)) //issue 405: we now save bvals file for isotropic series + isIsotropic = true; if (*numADC > 0) { // DWIs (i.e. short diffusion scans with too few directions to // calculate tensors...they typically acquire b=0 + 3 b > 0 so @@ -1790,14 +1941,16 @@ int * nii_saveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],str } else{ if ((numDti - *numADC) < 2) { - if (!dcmList[indx0].isDerived) //no need to warn if images are derived Trace/ND pair + /*if (!dcmList[indx0].isDerived) //no need to warn if images are derived Trace/ND pair printWarning("No bvec/bval files created: only single value after ADC excluded\n"); *numADC = 0; free(bvals); free(vx); - return NULL; - } - printMessage("Note: %d volumes appear to be ADC or trace images that will be removed to allow processing\n", + return NULL;*/ + printMessage("Warning: Isotropic DWI series, all bvecs are zero (issue 405)\n"); + *numADC = 0; + } else + printMessage("Note: %d volumes appear to be ADC or trace images that will be removed to allow processing\n", *numADC); } } @@ -1839,6 +1992,22 @@ int * nii_saveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],str } if (!isSequential) printMessage("DTI volumes re-ordered by ascending b-value\n"); + #ifdef RAW_BVEC //save raw bvecs as reported by DICOM, not rotated to image space + char txtname[2048] = {""}; + strcpy(txtname,pathoutname); + strcat (txtname,".rvec"); + FILE *fp = fopen(txtname, "w"); + for (int i = 0; i < numDti; i++) + fprintf(fp, "%g\t",vx[i].V[1]); + fprintf(fp, "\n"); + for (int i = 0; i < numDti; i++) + fprintf(fp, "%g\t",vx[i].V[2]); + fprintf(fp, "\n"); + for (int i = 0; i < numDti; i++) + fprintf(fp, "%g\t",vx[i].V[3]); + fprintf(fp, "\n"); + fclose(fp); + #endif dcmList[indx0].CSA.numDti = numDti; //warning structure not changed outside scope! geCorrectBvecs(&dcmList[indx0],sliceDir, vx, opts.isVerbose); siemensPhilipsCorrectBvecs(&dcmList[indx0],sliceDir, vx, opts.isVerbose); @@ -1866,7 +2035,6 @@ int * nii_saveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],str bVectors[i+j*numDti] = vx[i].V[j+1]; } // The image hasn't been created yet, so the attributes must be deferred - ImageList *images = (ImageList *) opts.imageList; images->addDeferredAttribute("bValues", bValues); images->addDeferredAttribute("bVectors", bVectors, numDti, 3); #else @@ -1898,6 +2066,10 @@ int * nii_saveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],str } fprintf(fp, "%g\n", vx[numDti-1].V[0]); fclose(fp); + if (isIsotropic) { //issue 405: ISOTROPIC images have bval but not bvec + free(vx); + return volOrderIndex; + } strcpy(txtname,pathoutname); if (dcmList[indx0].isVectorFromBMatrix) strcat (txtname,".mvec"); @@ -2225,7 +2397,37 @@ bool isExt (char *file_name, const char* ext) { return false; }// isExt() - +void cleanISO8859(char * cString) { + int len = strlen(cString); + if (len < 1) return; + for (int i = 0; i < len; i++) + //assume specificCharacterSet (0008,0005) is ISO_IR 100 http://en.wikipedia.org/wiki/ISO/IEC_8859-1 + if (cString[i]< 1) { + unsigned char c = (unsigned char)cString[i]; + if ((c >= 192) && (c <= 198)) cString[i] = 'A'; + if (c == 199) cString[i] = 'C'; + if ((c >= 200) && (c <= 203)) cString[i] = 'E'; + if ((c >= 204) && (c <= 207)) cString[i] = 'I'; + if (c == 208) cString[i] = 'D'; + if (c == 209) cString[i] = 'N'; + if ((c >= 210) && (c <= 214)) cString[i] = 'O'; + if (c == 215) cString[i] = 'x'; + if (c == 216) cString[i] = 'O'; + if ((c >= 217) && (c <= 220)) cString[i] = 'O'; + if (c == 221) cString[i] = 'Y'; + if ((c >= 224) && (c <= 230)) cString[i] = 'a'; + if (c == 231) cString[i] = 'c'; + if ((c >= 232) && (c <= 235)) cString[i] = 'e'; + if ((c >= 236) && (c <= 239)) cString[i] = 'i'; + if (c == 240) cString[i] = 'o'; + if (c == 241) cString[i] = 'n'; + if ((c >= 242) && (c <= 246)) cString[i] = 'o'; + if (c == 248) cString[i] = 'o'; + if ((c >= 249) && (c <= 252)) cString[i] = 'u'; + if (c == 253) cString[i] = 'y'; + if (c == 255) cString[i] = 'y'; + } +} int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopts opts) { char pth[PATH_MAX] = {""}; if (strlen(opts.outdir) > 0) { @@ -2260,6 +2462,10 @@ int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopt if (strlen(inname) < 1) { strcpy(inname, "T%t_N%n_S%s"); } + const char kTempPathSeparator ='\a'; + for (size_t pos = 0; pos') || (outname[pos] == ':') + if ((outname[pos] == '\\') || (outname[pos] == '/') || (outname[pos] == ' ') || (outname[pos] == '<') || (outname[pos] == '>') || (outname[pos] == ':') || (outname[pos] == '"') // || (outname[pos] == '/') || (outname[pos] == '\\') - || (outname[pos] == '^') + //|| (outname[pos] == '^') issue398 || (outname[pos] == '*') || (outname[pos] == '|') || (outname[pos] == '?')) outname[pos] = '_'; - #if defined(_WIN64) || defined(_WIN32) - const char kForeignPathSeparator ='/'; - #else - const char kForeignPathSeparator ='\\'; - #endif - for (int pos = 0; pos 0) && (outname[inpos] == '_') && (outname[outpos-1] == '_')) + continue; + outname[outpos] = outname[inpos]; + outpos++; + } + outname[outpos] = 0; + //printMessage("path='%s' name='%s'\n", pathoutname, outname); //make sure outname is unique strcat (baseoutname,outname); char pathoutname[2048] = {""}; @@ -2737,38 +2953,21 @@ void nii_saveAttributes (struct TDICOMdata &data, struct nifti_1_header &header, { ImageList *images = (ImageList *) opts.imageList; switch (data.modality) { - case kMODALITY_CR: - images->addAttribute("modality", "CR"); - break; - case kMODALITY_CT: - images->addAttribute("modality", "CT"); - break; - case kMODALITY_MR: - images->addAttribute("modality", "MR"); - break; - case kMODALITY_PT: - images->addAttribute("modality", "PT"); - break; - case kMODALITY_US: - images->addAttribute("modality", "US"); - break; + case kMODALITY_CR: images->addAttribute("modality", "CR"); break; + case kMODALITY_CT: images->addAttribute("modality", "CT"); break; + case kMODALITY_MR: images->addAttribute("modality", "MR"); break; + case kMODALITY_PT: images->addAttribute("modality", "PT"); break; + case kMODALITY_US: images->addAttribute("modality", "US"); break; } switch (data.manufacturer) { - case kMANUFACTURER_SIEMENS: - images->addAttribute("manufacturer", "Siemens"); - break; - case kMANUFACTURER_GE: - images->addAttribute("manufacturer", "GE"); - break; - case kMANUFACTURER_PHILIPS: - images->addAttribute("manufacturer", "Philips"); - break; - case kMANUFACTURER_TOSHIBA: - images->addAttribute("manufacturer", "Toshiba"); - break; - case kMANUFACTURER_Canon: - images->addAttribute("manufacturer", "Canon"); - break; + case kMANUFACTURER_SIEMENS: images->addAttribute("manufacturer", "Siemens"); break; + case kMANUFACTURER_GE: images->addAttribute("manufacturer", "GE"); break; + case kMANUFACTURER_PHILIPS: images->addAttribute("manufacturer", "Philips"); break; + case kMANUFACTURER_TOSHIBA: images->addAttribute("manufacturer", "Toshiba"); break; + case kMANUFACTURER_UIH: images->addAttribute("manufacturer", "UIH"); break; + case kMANUFACTURER_BRUKER: images->addAttribute("manufacturer", "Bruker"); break; + case kMANUFACTURER_HITACHI: images->addAttribute("manufacturer", "Hitachi"); break; + case kMANUFACTURER_CANON: images->addAttribute("manufacturer", "Canon"); break; } if (strlen(data.manufacturersModelName) > 0) images->addAttribute("scannerModelName", data.manufacturersModelName); @@ -2814,11 +3013,11 @@ void nii_saveAttributes (struct TDICOMdata &data, struct nifti_1_header &header, // See the nii_SaveBIDS() function for details int reconMatrixPE = data.phaseEncodingLines; if ((header.dim[2] > 0) && (header.dim[1] > 0)) { - if (header.dim[2] == header.dim[2]) //phase encoding does not matter - reconMatrixPE = header.dim[2]; - else if (data.phaseEncodingRC =='R') + if (header.dim[1] == header.dim[2]) //phase encoding does not matter reconMatrixPE = header.dim[2]; else if (data.phaseEncodingRC =='C') + reconMatrixPE = header.dim[2]; + else if (data.phaseEncodingRC =='R') reconMatrixPE = header.dim[1]; } @@ -3142,7 +3341,6 @@ int nii_saveNRRD(char * niiFilename, struct nifti_1_header hdr, unsigned char* i } else fprintf(fp,"DWMRI_gradient_%04d:=%.17g %.17g %.17g\n", i, factor*dti4D->S[i].V[1], factor*dti4D->S[i].V[2], factor*dti4D->S[i].V[3]); //printf("%g = %g %g %g>>>>\n",dti4D->S[i].V[0], dti4D->S[i].V[1],dti4D->S[i].V[2],dti4D->S[i].V[3]); - } } fprintf(fp,"\n"); //blank line: end of NRRD header @@ -3170,6 +3368,70 @@ int nii_saveNRRD(char * niiFilename, struct nifti_1_header hdr, unsigned char* i return pigz_File(fname, opts, imgsz); } // nii_saveNRRD() + +#ifndef max + #define max(a,b) \ + ({ __typeof__ (a) _a = (a); \ + __typeof__ (b) _b = (b); \ + _a > _b ? _a : _b; }) +#endif + +#ifndef min + #define min(a,b) \ + ({ __typeof__ (a) _a = (a); \ + __typeof__ (b) _b = (b); \ + _a < _b ? _a : _b; }) +#endif + +void removeSclSlopeInter(struct nifti_1_header* hdr, unsigned char* img) { + //NRRD does not have scl_slope scl_inter. Adjust data if possible + // https://discourse.slicer.org/t/preserve-image-rescale-and-slope-when-saving-in-nrrd-file/13357 + if (isSameFloat(hdr->scl_inter,0.0) && isSameFloat(hdr->scl_slope,1.0)) return; + if ((!isSameFloat(fmod(hdr->scl_inter, 1.0),0.0)) || (!isSameFloat(fmod(hdr->scl_slope, 1.0),0.0))) return; + int nVox = 1; + for (int i = 1; i < 8; i++) + if (hdr->dim[i] > 1) nVox = nVox * hdr->dim[i]; + if (hdr->datatype == DT_INT16) { + int16_t * img16 = (int16_t*) img; + int16_t mn, mx; + mn = img16[0]; + mx = mn; + for (int i=0; i < nVox; i++) { + mn = min(mn, img16[i]); + mx = max(mx, img16[i]); + } + float v = (mn * hdr->scl_slope) + hdr->scl_inter; + if ((v < -32768) || (v > 32767)) return; + v = (mx * hdr->scl_slope) + hdr->scl_inter; + if ((v < -32768) || (v > 32767)) return; + for (int i=0; i < nVox; i++) + img16[i] = round((img16[i] * hdr->scl_slope) + hdr->scl_inter); + hdr->scl_slope = 1.0; + hdr->scl_inter = 0.0; + return; + } + if (hdr->datatype == DT_UINT16) { + uint16_t * img16 = (uint16_t*) img; + uint16_t mn, mx; + mn = img16[0]; + mx = mn; + for (int i=0; i < nVox; i++) { + mn = min(mn, img16[i]); + mx = max(mx, img16[i]); + } + float v = (mn * hdr->scl_slope) + hdr->scl_inter; + if ((v < 0) || (v > 65535)) return; + v = (mx * hdr->scl_slope) + hdr->scl_inter; + if ((v < 0) || (v > 65535)) return; + for (int i=0; i < nVox; i++) + img16[i] = round((img16[i] * hdr->scl_slope) + hdr->scl_inter); + hdr->scl_slope = 1.0; + hdr->scl_inter = 0.0; + return; + } + //printWarning("NRRD unable to record scl_slope/scl_inter %g/%g\n", hdr->scl_slope, hdr->scl_inter); +} + void swapEndian(struct nifti_1_header* hdr, unsigned char* im, bool isNative) { //swap endian from big->little or little->big // must be told which is native to detect datatype and number of voxels @@ -3510,13 +3772,11 @@ void adjustOriginForNegativeTilt(struct nifti_1_header * hdr, float shiftPxY) { hdr->srow_y[3] -= shiftPxY * hdr->srow_y[1]; hdr->srow_z[3] -= shiftPxY * hdr->srow_y[2]; } - - if (hdr->qform_code > 0) { + if (hdr->qform_code > 0) { // Adjust the quaternion offsets using quatern_* and pixdim mat44 mat = nifti_quatern_to_mat44(hdr->quatern_b, hdr->quatern_c, hdr->quatern_d, hdr->qoffset_x, hdr->qoffset_y, hdr->qoffset_z, hdr->pixdim[1], hdr->pixdim[2], hdr->pixdim[3], hdr->pixdim[0]); - hdr->qoffset_x -= shiftPxY * mat.m[1][0]; hdr->qoffset_y -= shiftPxY * mat.m[1][1]; hdr->qoffset_z -= shiftPxY * mat.m[1][2]; @@ -4151,7 +4411,7 @@ int nii_saveCrop(char * niiFilename, struct nifti_1_header hdr, unsigned char* i return returnCode; }// nii_saveCrop() -float dicomTimeToSec (float dicomTime) { +double dicomTimeToSec (double dicomTime) { //convert HHMMSS to seconds, 135300.024 -> 135259.731 are 0.293 sec apart char acqTimeBuf[64]; snprintf(acqTimeBuf, sizeof acqTimeBuf, "%+013.5f", (double)dicomTime); @@ -4163,14 +4423,15 @@ float dicomTimeToSec (float dicomTime) { return (ahour * 3600)+(amin * 60) + asec; } -float acquisitionTimeDifference(struct TDICOMdata * d1, struct TDICOMdata * d2) { +double acquisitionTimeDifference(struct TDICOMdata * d1, struct TDICOMdata * d2) { if (d1->acquisitionDate != d2->acquisitionDate) return -1; //to do: scans running across midnight - float sec1 = dicomTimeToSec(d1->acquisitionTime); - float sec2 = dicomTimeToSec(d2->acquisitionTime); + double sec1 = dicomTimeToSec(d1->acquisitionTime); + double sec2 = dicomTimeToSec(d2->acquisitionTime); //printMessage("%g\n",d2->acquisitionTime); if ((sec1 < 0) || (sec2 < 0)) return -1; return (sec2 - sec1); } + void checkDateTimeOrder(struct TDICOMdata * d, struct TDICOMdata * d1) { if (d->acquisitionDate < d1->acquisitionDate) return; //d1 occurred on later date if (d->acquisitionTime <= d1->acquisitionTime) return; //d1 occurred on later (or same) time @@ -4184,6 +4445,8 @@ void checkSliceTiming(struct TDICOMdata * d, struct TDICOMdata * d1, int verbose //detect images with slice timing errors. https://github.com/rordenlab/dcm2niix/issues/126 //modified 20190704: this function now ensures all slice times are in msec if ((d->TR < 0.0) || (d->CSA.sliceTiming[0] < 0.0)) return; //no slice timing + if (d->manufacturer == kMANUFACTURER_GE) return; //compute directly from Protocol Block + if (d->modality == kMODALITY_PT) return; //issue407 int nSlices = 0; while ((nSlices < kMaxEPI3D) && (d->CSA.sliceTiming[nSlices] >= 0.0)) nSlices++; @@ -4234,6 +4497,7 @@ void checkSliceTiming(struct TDICOMdata * d, struct TDICOMdata * d1, int verbose if ((minT == maxT) && (d->is3DAcq)) return; //fine: 3D EPI if ((minT == maxT) && (d->CSA.multiBandFactor == d->CSA.mosaicSlices)) return; //fine: all slices single excitation if ((strlen(d->seriesDescription) > 0) && (strstr(d->seriesDescription, "SBRef") != NULL)) return; //fine: single-band calibration data, the slice timing WILL exceed the TR + if (verbose > 1) printMessage("Slice timing range of first volume: range %g..%g, TR=%g ms)\n", minT, maxT, TRms); //check if 2nd image has valid slice timing float minT1 = d1->CSA.sliceTiming[0]; float maxT1 = minT1; @@ -4242,9 +4506,17 @@ void checkSliceTiming(struct TDICOMdata * d, struct TDICOMdata * d1, int verbose if (d1->CSA.sliceTiming[i] < minT1) minT1 = d1->CSA.sliceTiming[i]; if (d1->CSA.sliceTiming[i] > maxT1) maxT1 = d1->CSA.sliceTiming[i]; } + if (verbose > 1) printMessage("Slice timing range of 2nd volume: range %g..%g, TR=%g ms)\n", minT, maxT, TRms); + if ((minT1 < maxT1) && (minT1 > 0.0) && ((maxT1-minT1) <= TRms) ) { //issue 429: 2nd volume may not start from zero + for (int i = 0; i < nSlices; i++) + d1->CSA.sliceTiming[i] -= minT1; + maxT1 -= minT1; + minT1 -= minT1; + } if ((minT1 < 0.0) && (d->rtia_timerGE >= 0.0)) return; //use rtia timer if (minT1 < 0.0) { //https://github.com/neurolabusc/MRIcroGL/issues/31 - printWarning("Siemens MoCo? Bogus slice timing (range %g..%g, TR=%g seconds)\n", minT1, maxT1, TRms); + if (d->modality == kMODALITY_MR) + printWarning("Siemens MoCo? Bogus slice timing (range %g..%g, TR=%g seconds)\n", minT1, maxT1, TRms); return; } if ((minT1 == maxT1) || (maxT1 >= TRms)) { //both first and second image corrupted @@ -4265,7 +4537,7 @@ void sliceTimingXA(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct n // Ignore first volume: For an example of erroneous first volume timing, see series 10 (Functional_w_SMS=3) https://github.com/rordenlab/dcm2niix/issues/240 // an alternative would be to use 0018,9074 - this would need to be converted from DT to Secs, and is scrambled if de-identifies data see enhanced de-identified series 26 from issue 236 uint64_t indx0 = dcmSort[0].indx; //first volume - if (!dcmList[indx0].isXA10A) return; + if ((!dcmList[indx0].isXA10A) || (hdr->dim[3] < 1)) return; if ( (nConvert == (hdr->dim[3]*hdr->dim[4])) && (hdr->dim[3] < (kMaxEPI3D-1)) && (hdr->dim[3] > 1) && (hdr->dim[4] > 1)) { //XA11 2D classic for (int v = 0; v < hdr->dim[3]; v++) @@ -4286,21 +4558,172 @@ void sliceTimingXA(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct n } if ((dcmList[indx0].CSA.multiBandFactor < 2) && (mb > 1)) dcmList[indx0].CSA.multiBandFactor = mb; - //for (int v = 0; v < hdr->dim[3]; v++) - // printf("XA10sliceTiming\t%d\t%g\n", v, dcmList[indx0].CSA.sliceTiming[v]); + return; //we have subtracted min } + //issue429: subtract min + float mn = dcmList[indx0].CSA.sliceTiming[0]; + for (int v = 0; v < hdr->dim[3]; v++) + mn = min(mn, dcmList[indx0].CSA.sliceTiming[v]); + if (isSameFloatGE(mn, 0.0)) return; + for (int v = 0; v < hdr->dim[3]; v++) + dcmList[indx0].CSA.sliceTiming[v] -= mn; } //sliceTimingXA() -void rescueSliceTimingGE(struct TDICOMdata * d, int verbose, int nSL, const char * filename) { +void sliceTimeGE (struct TDICOMdata * d, int mb, int dim3, float TR, bool isInterleaved, bool is27r3, float groupDelaysec) { +//mb : multiband factor +//dim3 : number of slices in volume +//TRsec : repetition time in seconds +//isInterleaved : interleaved or sequential slice order +//is27r3 : software release 27.0 R03 or later + float sliceTiming[kMaxEPI3D]; + //multiband can be fractional! 'extra' slices discarded + int nExcitations = ceil(float(dim3) / float(mb)); + int nDiscardedSlices = (nExcitations * mb) - dim3; + float secPerSlice = (TR - groupDelaysec) / (nExcitations); + if (!isInterleaved) { + for (int i = 0; i < nExcitations; i++) + sliceTiming[i] = i * secPerSlice; + } else { + int nOdd = (nExcitations - 1) / 2; + for (int i = 0; i < nExcitations; i++) { + if (i % 2 == 0) //ODD slices since we index from 0! + sliceTiming[i] = (i/2) * secPerSlice; + else + sliceTiming[i] = (nOdd+((i+1)/2)) * secPerSlice; + } //for each slice + if ((mb > 1) && (is27r3) && (isInterleaved) && (nExcitations > 2) && ((nExcitations % 2) == 0)) { + float tmp = sliceTiming[nExcitations - 1]; + sliceTiming[nExcitations - 1] = sliceTiming[nExcitations - 3]; + sliceTiming[nExcitations - 3] = tmp; + //printf("SWAP!\n"); + } + } //if interleaved + for (int i = 0; i < dim3; i++) + sliceTiming[i] = sliceTiming[i % nExcitations]; + //#define testSliceTimesGE //note that early GE HyperBand sequences reported single-band values in x0021x105E + #ifdef testSliceTimesGE + float maxErr = 0.0; + for (int i = 0; i < dim3; i++) + maxErr = max(maxErr, fabs(sliceTiming[i] - d->CSA.sliceTiming[i])); + if ((d->CSA.sliceTiming[0] >= 0.0) && (maxErr > 1.0) ) { //allow a 1.0 msec tolerance for rounding + printMessage("GE estimated slice times differ from reported (max error: %g)\n", maxErr); + printMessage("Slice\tEstimated\tReported\n"); + for (int i = 0; i < dim3; i++) { + printMessage("%d %g %g\n", i, sliceTiming[i], d->CSA.sliceTiming[i]); + maxErr = max(maxErr, fabs(sliceTiming[i] - d->CSA.sliceTiming[i])); + } + } + #endif + for (int i = 0; i < dim3; i++) + d->CSA.sliceTiming[i] = sliceTiming[i]; +} // sliceTimeGE() + +void readSoftwareVersionsGE(char softwareVersionsGE[], int verbose,char geVersionPrefix[], float* geMajorVersion, int* geMajorVersionInt, int* geMinorVersionInt, int* geReleaseVersionInt, bool* is27r3) { + // softwareVersionsGE + // "27\LX\MR Software release:RX27.0_R02_1831.a" -> 27 + // "28\LX\MR29.1_EA_2039.g" -> 29 + // geVersionPrefix + // RX27.0_R02_1831.a -> RX + // MR29.1_EA_2039.g -> MR + // geMajorVersion + // RX27.0_R02_1831.a -> 27.0 + // MR29.1_EA_2039.g -> 29.1 + // geMajorVersionInt + // RX27.0_R02_1831.a -> 27 + // MR29.1_EA_2039.g -> 29 + // geMinorVersionInt + // RX27.0_R02_1831.a -> 0 + // MR29.1_EA_2039.g -> 1 + // geReleaseVersionInt + // EA->0, R01->1, R02->2, R03->4 + // RX27.0_R02_1831.a -> 2 + // MR29.1_EA_2039.g -> 0 + + int len = 0; + // If softwareVersionsGE is 27\LX\MR Software release:RX27.0_R02_1831.a + char *sepStart = strchr(softwareVersionsGE, ':'); + if (sepStart == NULL) { + // If softwareVersionsGE is 28\LX\MR29.1_EA_2039.g + sepStart = strrchr(softwareVersionsGE, '\\'); + if (sepStart == NULL) return; + } + sepStart += 1; + len = 11; + char * versionString = (char *)malloc(sizeof(char) * len); + versionString[len] =0; + memcpy(versionString, sepStart, len); + int ver1, ver2, ver3; + char c1, c2, c3, c4; + // RX27.0_R02_ or MR29.1_EA_2 + sscanf( versionString, "%c%c%d.%d_%c%c%d", &c1, &c2, geMajorVersionInt, geMinorVersionInt, &c3, &c4, geReleaseVersionInt); + memcpy(geVersionPrefix, &c1, 1); + memcpy(geVersionPrefix+1, &c2, 1); + if ( (c3 == 'E') && (c4 == 'A') ){ + *geReleaseVersionInt = 0; + } + free(versionString); + *geMajorVersion = (float) *geMajorVersionInt + (float) 0.1 * (float) *geMinorVersionInt; + *is27r3 = ((*geMajorVersion >= 27.1) || ((*geMajorVersionInt == 27) && (*geReleaseVersionInt >= 3))); + if (verbose > 1) { + printMessage("GE Software VersionPrefix: %s\n", geVersionPrefix); + printMessage("GE Software MajorVersion: %d\n", *geMajorVersionInt); + printMessage("GE Software MinorVersion: %d\n", *geMinorVersionInt); + printMessage("GE Software ReleaseVersion: %d\n", *geReleaseVersionInt); + printMessage("GE Software is27r3: %d\n", *is27r3); + } +} // readSoftwareVersionsGE() + +void sliceTimingGE_Testx0021x105E(struct TDICOMdata * d, struct TDCMopts opts, struct nifti_1_header * hdr, struct TDCMsort *dcmSort,struct TDICOMdata *dcmList) { + if ((!opts.isTestx0021x105E) || (hdr->dim[3] < 2) || (hdr->dim[4] < 1)) return; + if (d->rtia_timerGE <= 0.0) { + printMessage("DICOM images do not report RTIA timer(0021,105E)\n"); + return; + } + int j = hdr->dim[3]; + float sliceTiming[kMaxEPI3D]; + float mn = INFINITY; + for (int v = 0; v < hdr->dim[3]; v++) { + sliceTiming[v] = dcmList[dcmSort[v+j].indx].rtia_timerGE; + mn = min(mn, sliceTiming[v]); + } + if (mn < 0.0) return; + float mxErr = 0.0; + for (int v = 0; v < hdr->dim[3]; v++) { + sliceTiming[v] = (sliceTiming[v] - mn) * 1000.0; //subtract offset, convert sec -> ms + mxErr = max(mxErr, fabs(sliceTiming[v] - d->CSA.sliceTiming[v])); + } + printMessage("Slice Timing Error between calculated and RTIA timer(0021,105E): %gms\n", mxErr); + if ((mxErr < 1.0) && (opts.isVerbose < 1)) return; + for (int v = 0; v < hdr->dim[3]; v++) + printMessage("\t%g\t%g\n", d->CSA.sliceTiming[v], sliceTiming[v]); +} + +void sliceTimingGE(struct TDICOMdata * d, const char * filename, struct TDCMopts opts, struct nifti_1_header * hdr, struct TDCMsort *dcmSort,struct TDICOMdata *dcmList) { //we can often read GE slice timing from TriggerTime (0018,1060) or RTIA Timer (0021,105E) // if both of these methods fail, we can often guess based on slice order stored in the Private Protocol Data Block (0025,101B) // this is referred to as "rescue" as we only know the TR, not the TA. So assumes continuous scans with no gap - if (d->is3DAcq) return; //no need for slice times - if (nSL < 2) return; + if ((d->is3DAcq) || (d->isLocalizer) || (hdr->dim[4] < 2)) return; //no need for slice times + if (hdr->dim[3] < 2) return; if (d->manufacturer != kMANUFACTURER_GE) return; - if (d->maxEchoNumGE > 0) - printWarning("GE sequence with %d echoes. See https://github.com/rordenlab/dcm2niix/issues/359\n", d->maxEchoNumGE); - + //start version check: + float geMajorVersion = 0; + int geMajorVersionInt = 0, geMinorVersionInt = 0, geReleaseVersionInt = 0; + char geVersionPrefix[2] = " "; + bool is27r3 = false; + readSoftwareVersionsGE(d->softwareVersions, opts.isVerbose, geVersionPrefix, &geMajorVersion, &geMajorVersionInt, &geMinorVersionInt, &geReleaseVersionInt, &is27r3); + //readSoftwareVersionsGE(&geMajorVersion); + /* + //used for oldSliceTimingGE + if (!opts.isIgnorex0021x105E) { + if ((geMajorVersionInt >= 28) && (d->CSA.sliceTiming[0] >= 0.0)) { + //if (opts.isVerbose > 1) + printMessage("GEversion %.1f, slice timing from DICOM (0021,105E).\n", geMajorVersion); + return; //trust slice timings for versions > 27, see issue 336 + } + }*/ + //end: version check + if (d->maxEchoNumGE > 0) + printWarning("GE sequence with %d echoes. See issue 359\n", d->maxEchoNumGE); if ((d->protocolBlockStartGE < 1) || (d->protocolBlockLengthGE < 19)) return; #ifdef myReadGeProtocolBlock //GE final desperate attempt to determine slice order @@ -4310,54 +4733,78 @@ void rescueSliceTimingGE(struct TDICOMdata * d, int verbose, int nSL, const char int viewOrderGE = -1; int sliceOrderGE = -1; int mbAccel = -1; + int nSlices = -1; + float groupDelay = 0.0; + char ioptGE[3000] = ""; //printWarning("Using GE Protocol Data Block for BIDS data (beware: new feature)\n"); - int ok = geProtocolBlock(filename, d->protocolBlockStartGE, d->protocolBlockLengthGE, verbose, &sliceOrderGE, &viewOrderGE, &mbAccel); + int ok = geProtocolBlock(filename, d->protocolBlockStartGE, d->protocolBlockLengthGE, opts.isVerbose, &sliceOrderGE, &viewOrderGE, &mbAccel, &nSlices, &groupDelay, ioptGE); if (ok != EXIT_SUCCESS) { - printWarning("Unable to decode GE protocol block\n"); + printWarning("Unable to estimate slice times: issue decoding GE protocol block.\n"); return; } - if (mbAccel > 1) { - d->CSA.multiBandFactor = mbAccel; - printWarning("Unabled to compute slice times for GE multi-band. SliceOrder=%d (seq=0, int=1)\n", sliceOrderGE); + mbAccel = max(mbAccel, 1); + if (nSlices != hdr->dim[3]) //redundant with locationsInAcquisition check? + printWarning("Missing DICOMs, number of slices estimated (%d) differs from Protocol Block (0025,101B) report (%d).\n", hdr->dim[3], nSlices); + d->CSA.multiBandFactor = max(d->CSA.multiBandFactor, mbAccel); + bool isInterleaved = (sliceOrderGE != 0); + groupDelay *= 1000.0; //sec -> ms + // + // Estimate GE Slice Time only for EPI Multi-Phase (epi) or EPI fMRI/BrainWave (epiRT) + // + // BrainWave (epiRT) + if ((d->epiVersionGE == 1) || (strstr(ioptGE,"FMRI") != NULL)) { //-1 = not epi, 0 = epi, 1 = epiRT + d->epiVersionGE = 1; + d->internalepiVersionGE = 1; // 'EPI'(gradient echo)/'EPI2'(spin echo) + if (!isSameFloatGE(groupDelay, d->groupDelay)) + printWarning("With epiRT (i.e. FMRI option), Group delay reported in private tag (0043,107C = %g) and Protocol Block (0025,101B = %g) differ.\n", d->groupDelay, groupDelay); + } + // EPI Multi-Phase (epi) + else if ((d->epiVersionGE == 0) || (strstr(ioptGE,"MPh") != NULL)) { //-1 = not epi, 0 = epi, 1 = epiRT + d->epiVersionGE = 0; + d->internalepiVersionGE = 1; // 'EPI'(gradient echo)/'EPI2'(spin echo) + if (groupDelay > 0) { + d->TR += groupDelay; + d->groupDelay = groupDelay; + } + // EPI Multi-Phase (epi) with Variable Delays (Unsupported) + if (groupDelay < -0.5) { + printWarning("SliceTiming Unspported: GE Multi-Phase EPI with Variable Delays\n"); + d->CSA.sliceTiming[0] = -1; + return; + } + } + // Diffusion (Unsupported) + else if ( (d->epiVersionGE == 2) || (d->internalepiVersionGE == 2) || (strstr(ioptGE,"DIFF") != NULL) ) { + printWarning("Unable to compute slice times for GE Diffusion\n"); d->CSA.sliceTiming[0] = -1.0; return; - } - if (d->CSA.sliceTiming[0] >= 0.0) return; //slice times calculated - moved here to detect multiband, see issue 336 - - if ((sliceOrderGE < 0) || (sliceOrderGE > 1)) return; - // 0=sequential/1=interleaved - printWarning("Guessing slice times using ProtocolBlock SliceOrder=%d (seq=0, int=1)\n", sliceOrderGE); - int nOdd = nSL / 2; - float secPerSlice = d->TR/nSL; //should be TA not TR! We do not know TR - if (sliceOrderGE == 0) { - for (int i = 0; i < nSL; i++) - d->CSA.sliceTiming[i] = i * secPerSlice; - } else { - for (int i = 0; i < nSL; i++) { - if (i % 2 == 0) { //ODD slices since we index from 0! - d->CSA.sliceTiming[i] = (i/2) * secPerSlice; - //printf("%g\n", d->CSA.sliceTiming[i]); - } else { - d->CSA.sliceTiming[i] = (nOdd+((i+1)/2)) * secPerSlice; - //printf("%g\n", d->CSA.sliceTiming[i]); - } - } //for each slice - } //if interleaved + // Others (Unsupported) + else { + printWarning("Unable to compute slice times for this GE dataset\n"); + d->CSA.sliceTiming[0] = -1.0; + return; + } + if (opts.isVerbose > 1) { + printMessage("GEiopt: %s, groupDelay (%g), internalepiVersionGE (%d), epiVersionGE(%d)\n", ioptGE, groupDelay, d->internalepiVersionGE, d->epiVersionGE); + printMessage("GEversion %s%.1f_R0%d, TRms %g, interleaved %d, multiband %d, groupdelayms %g\n", geVersionPrefix, geMajorVersion, geReleaseVersionInt, d->TR, isInterleaved, d->CSA.multiBandFactor, d->groupDelay); + } + sliceTimeGE(d, d->CSA.multiBandFactor, hdr->dim[3], d->TR, isInterleaved, is27r3, d->groupDelay); + sliceTimingGE_Testx0021x105E(d, opts, hdr, dcmSort, dcmList); #endif -} +} //sliceTimingGE() -void reverseSliceTiming(struct TDICOMdata * d, int verbose, int nSL) { +void reverseSliceTiming(struct TDICOMdata * d, int verbose, int nSL) { if ((d->CSA.protocolSliceNumber1 == 0) || ((d->CSA.protocolSliceNumber1 == 1))) return; //slices not flipped if (d->is3DAcq) return; //no need for slice times if (d->CSA.sliceTiming[0] < 0.0) return; //no slice times if (nSL > kMaxEPI3D) return; if (nSL < 2) return; if (verbose) - printMessage("Slices were spatially flipped, so slice times are flipped\n"); - d->CSA.protocolSliceNumber1 = 0; - float sliceTiming[kMaxEPI3D]; - for (int i = 0; i < nSL; i++) + printMessage("Slices were spatially flipped, so slice times are flipped\n"); + d->CSA.protocolSliceNumber1 = 0; + float sliceTiming[kMaxEPI3D]; + for (int i = 0; i < nSL; i++) sliceTiming[i] = d->CSA.sliceTiming[i]; for (int i = 0; i < nSL; i++) d->CSA.sliceTiming[i] = sliceTiming[(nSL-1)-i]; @@ -4442,7 +4889,8 @@ void sliceTimingUIH(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct dcmList[indx0].CSA.sliceTiming[v] = dcmList[dcmSort[v].indx].acquisitionTime; //nb format is HHMMSS we need to handle midnight-crossing and convert to ms, see checkSliceTiming() } -void sliceTimingGE(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct nifti_1_header * hdr, int verbose, const char * filename, int nConvert) { +/* +void oldSliceTimingGE(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct nifti_1_header * hdr, int verbose, const char * filename, int nConvert) { //GE check slice timing >>> uint64_t indx0 = dcmSort[0].indx; //first volume if (!(dcmList[indx0].manufacturer == kMANUFACTURER_GE)) return; @@ -4532,9 +4980,10 @@ void sliceTimingGE(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct n } //verbose > 1 } //if maxTime != minTIme } //GE slice timing from 0021,105E -} //sliceTimingGE() +} //oldSliceTimingGE() +*/ -int sliceTimingCore(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct nifti_1_header * hdr, int verbose, const char * filename, int nConvert) { +int sliceTimingCore(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct nifti_1_header * hdr, int verbose, const char * filename, int nConvert, struct TDCMopts opts) { int sliceDir = 0; if (hdr->dim[3] < 2) return sliceDir; //uint64_t indx0 = dcmSort[0].indx; @@ -4543,13 +4992,12 @@ int sliceTimingCore(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct uint64_t indx1 = dcmSort[1].indx; if (nConvert < 2) indx1 = dcmSort[0].indx; struct TDICOMdata * d1 = &dcmList[indx1]; - sliceTimingGE(dcmSort, dcmList, hdr, verbose, filename, nConvert); + //oldSliceTimingGE(dcmSort, dcmList, hdr, verbose, filename, nConvert); sliceTimingUIH(dcmSort, dcmList, hdr, verbose, filename, nConvert); int isSliceTimeHHMMSS = sliceTimingSiemens2D(dcmSort, dcmList, hdr, verbose, filename, nConvert); sliceTimingXA(dcmSort, dcmList, hdr, verbose, filename, nConvert); checkSliceTiming(d0, d1, verbose, isSliceTimeHHMMSS); rescueSliceTimingSiemens(d0, verbose, hdr->dim[3], filename); //desperate attempts if conventional methods fail - rescueSliceTimingGE(d0, verbose, hdr->dim[3], filename); //desperate attempts if conventional methods fail if (hdr->dim[3] > 1)sliceDir = headerDcm2Nii2(dcmList[dcmSort[0].indx], dcmList[indx1] , hdr, true); //UNCOMMENT NEXT TWO LINES TO RE-ORDER MOSAIC WHERE CSA's protocolSliceNumber does not start with 1 if (dcmList[dcmSort[0].indx].CSA.protocolSliceNumber1 > 1) { @@ -4562,6 +5010,7 @@ int sliceTimingCore(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct if ((dcmList[dcmSort[0].indx].manufacturer == kMANUFACTURER_UIH) || (dcmList[dcmSort[0].indx].manufacturer == kMANUFACTURER_GE)) dcmList[dcmSort[0].indx].CSA.protocolSliceNumber1 = -1; } + sliceTimingGE(d0, filename, opts, hdr, dcmSort, dcmList); reverseSliceTiming(d0, verbose, hdr->dim[3]); return sliceDir; } //sliceTiming() @@ -4588,6 +5037,47 @@ int sliceTimingCore(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct return 0; }*/ +void loadOverlay(char* imgname, unsigned char * img, int offset, int x, int y, int z) { + int nvox = x * y * z; + size_t imgszRead = (nvox+7) >> 3; //overlay stored as 1 bit per voxel + FILE *file = fopen(imgname , "rb"); + if (!file) { + printError("Unable to open '%s'\n", imgname); + return; + } + fseek(file, 0, SEEK_END); + long fileLen=ftell(file); + if (fileLen < (imgszRead+offset)) { + printWarning("File not large enough to store overlay: %s\n", imgname); + return; + } + fseek(file, (long) offset, SEEK_SET); + unsigned char *bImg = (unsigned char *)malloc(imgszRead); + size_t sz = fread(bImg, 1, imgszRead, file); + //static unsigned char mask[] = {128, 64, 32, 16, 8, 4, 2, 1}; + static unsigned char mask[] = {1, 2, 4, 8, 16, 32, 64, 128}; + for (int i = 0; i < nvox; i++) { + int byt = (i >> 3); + int bit = (i % 8); + img[i] = ((bImg[byt] & mask[bit]) != 0); + } + /* + if (isFlipY) { + unsigned char *tImg = (unsigned char *)malloc(nvox); + memcpy(&tImg[0], &img[0], nvox); + int i = 0; + for (int yi = y-1; yi >= 0; yi--) + for (int xi = 0; xi < x; xi++) { + img[(yi*x)+xi] = tImg[i]; + i++; + } + free(tImg); + }*/ + free(bImg); + fclose(file); + return; +} //loadOverlay() + int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmList[], struct TSearchList *nameList, struct TDCMopts opts, struct TDTI4D *dti4D, int segVol) { bool iVaries = intensityScaleVaries(nConvert,dcmSort,dcmList); float *sliceMMarray = NULL; //only used if slices are not equidistant @@ -4602,6 +5092,8 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc if (isnan(dcmList[indx0].gantryTilt)) return EXIT_FAILURE; } #endif //newTilt see issue 254 + if (dcmList[indx0].isPrivateCreatorRemap) + printWarning("PrivateCreator remapping detected. DICOMs are not archival quality (issue 435).\n"); if (dcmList[indx0].isScaleVariesEnh) //issue363 iVaries = true; if ((dcmList[indx].isXA10A) && (dcmList[indx].CSA.mosaicSlices < 0)) { @@ -4621,7 +5113,7 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc return EXIT_SUCCESS; } if ((opts.isIgnoreDerivedAnd2D) && ((dcmList[indx].isLocalizer) || (strcmp(dcmList[indx].sequenceName, "_tfl2d1")== 0) || (strcmp(dcmList[indx].sequenceName, "_fl3d1_ns")== 0) || (strcmp(dcmList[indx].sequenceName, "_fl2d1")== 0)) ) { - printMessage("Ignoring localizer (sequence %s) of series %ld %s\n", dcmList[indx].sequenceName, dcmList[indx].seriesNum, nameList->str[indx]); + printMessage("Ignoring localizer (sequence '%s') of series %ld %s\n", dcmList[indx].sequenceName, dcmList[indx].seriesNum, nameList->str[indx]); return EXIT_SUCCESS; } if ((opts.isIgnoreDerivedAnd2D) && (nConvert < 2) && (dcmList[indx].CSA.mosaicSlices < 2) && (dcmList[indx].xyzDim[3] < 2)) { @@ -4649,6 +5141,7 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc memcpy(&imgM[0], &img[0], imgsz); free(img); //printMessage(" %d %d %d %d %lu\n", hdr0.dim[1], hdr0.dim[2], hdr0.dim[3], hdr0.dim[4], (unsigned long)[imgM length]); + bool isHasOverlay = dcmList[indx0].isHasOverlay; if (nConvert > 1) { //next: detect trigger time see example https://www.slicer.org/wiki/Documentation/4.4/Modules/MultiVolumeExplorer double triggerDx = dcmList[dcmSort[nConvert-1].indx].triggerDelayTime - dcmList[indx0].triggerDelayTime; @@ -4695,6 +5188,8 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc hdr0.dim[3] = nConvert/nAcq; hdr0.dim[4] = nAcq; hdr0.dim[0] = 4; + if ((dcmList[indx0].locationsInAcquisition > 0) && (dcmList[indx0].locationsInAcquisition != hdr0.dim[3])) + printMessage("DICOM images may be missing, expected %d spatial locations per volume, but found %d.\n", dcmList[indx0].locationsInAcquisition, hdr0.dim[3]); } else if ((dcmList[indx0].isXA10A) && (nConvert > nAcq) && (nAcq > 1) ) { nAcq -= 1; hdr0.dim[3] = nConvert/nAcq; @@ -4708,11 +5203,12 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc if ((nAcq > 1) && (nConvert != nAcq)) { printMessage("Slice positions repeated, but number of slices (%d) not divisible by number of repeats (%d): missing images?\n", nConvert, nAcq); if (dcmList[indx0].locationsInAcquisition > 0) - printMessage("Hint: expected %d locations", dcmList[indx0].locationsInAcquisition); + printMessage("Hint: expected %d locations\n", dcmList[indx0].locationsInAcquisition); } } //next options removed: features now thoroughly detected in nii_loadDir() - for (int i = 0; i < nConvert; i++) { //make sure 1st volume describes shared features + for (int i = 0; i < nConvert; i++) { //make sure 1st volume describes shared features + if (dcmList[dcmSort[i].indx].isHasOverlay) isHasOverlay = true; if (dcmList[dcmSort[i].indx].isCoilVaries) dcmList[indx0].isCoilVaries = true; if (dcmList[dcmSort[i].indx].isMultiEcho) dcmList[indx0].isMultiEcho = true; if (dcmList[dcmSort[i].indx].isNonParallelSlices) dcmList[indx0].isNonParallelSlices = true; @@ -4721,26 +5217,72 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc if (dcmList[dcmSort[i].indx].isHasImaginary) dcmList[indx0].isHasImaginary = true; } //next: detect variable inter-volume time https://github.com/rordenlab/dcm2niix/issues/184 - if (dcmList[indx0].modality == kMODALITY_PT) { + //if ((nConvert > 1) && ((dcmList[indx0].modality == kMODALITY_PT)|| (opts.isForceOnsetTimes))) { + if ((nConvert > 1) && ((dcmList[indx0].modality == kMODALITY_PT)|| ((opts.isForceOnsetTimes) && (dcmList[indx0].manufacturer != kMANUFACTURER_GE)))) { + //note: GE 0008,0032 unreliable, see mb=6 data from sw27.0 20201026 + //issue 407 + int nTR = 0; + for (int i = 0; i < nConvert; i++) + if (isSamePosition(dcmList[indx0],dcmList[dcmSort[i].indx])) { + dti4D->frameDuration[nTR] = dcmList[dcmSort[i].indx].frameDuration; + nTR += 1; + if (nTR >= kMaxDTI4D) break; + } + bool trVaries = false; bool dayVaries = false; - float tr = -1; + float tr = -1.0; + float mintr = 1000000; + float maxtr = -1.0; + float volumeTimeStartFirstStartLast = -1.0; + int nVol = 0; uint64_t prevVolIndx = indx0; for (int i = 0; i < nConvert; i++) if (isSamePosition(dcmList[indx0],dcmList[dcmSort[i].indx])) { float trDiff = acquisitionTimeDifference(&dcmList[prevVolIndx], &dcmList[dcmSort[i].indx]); prevVolIndx = dcmSort[i].indx; + nVol ++; if (trDiff <= 0) continue; + mintr = min(mintr, trDiff); + maxtr = max(maxtr, trDiff); if (tr < 0) tr = trDiff; if (trDiff < 0) dayVaries = true; - if (!isSameFloatGE(tr,trDiff)) - trVaries = true; + float trDiff0 = acquisitionTimeDifference(&dcmList[indx0], &dcmList[dcmSort[i].indx]); + volumeTimeStartFirstStartLast = max(volumeTimeStartFirstStartLast, trDiff0); } + float toleranceSec = 50.0/1000.0; //e.g. 50/1000 = 50ms + if ((nVol > 1) && (volumeTimeStartFirstStartLast > 0.0)) { + tr = volumeTimeStartFirstStartLast / (nVol - 1.0); + if (fabs(tr - hdr0.pixdim[4]) > toleranceSec) { + if ((dcmList[indx0].isIR) && (dcmList[indx0].manufacturer != kMANUFACTURER_PHILIPS)) + dti4D->repetitionTimeInversion = hdr0.pixdim[4]; + else + dti4D->repetitionTimeExcitation = hdr0.pixdim[4]; + hdr0.pixdim[4] = tr; + dcmList[indx0].TR = tr * 1000.0; //as msec + } + } + if ((maxtr - mintr) > toleranceSec) trVaries = true; if (trVaries) { if (dayVaries) printWarning("Seconds between volumes varies (perhaps run through midnight)\n"); else printWarning("Seconds between volumes varies\n"); + //issue 407 + int nTR = 0; + for (int i = 0; i < nConvert; i++) + if (isSamePosition(dcmList[indx0],dcmList[dcmSort[i].indx])) { + float trDiff = acquisitionTimeDifference(&dcmList[indx0], &dcmList[dcmSort[i].indx]); + dti4D->volumeOnsetTime[nTR] = trDiff; + dti4D->decayFactor[nTR] = dcmList[dcmSort[i].indx].decayFactor; + //printf("%d %g\n", i, dcmList[dcmSort[i].indx].decayFactor); + nTR += 1; + if (nTR >= kMaxDTI4D) break; + } + if (dcmList[indx0].modality != kMODALITY_PT) + dti4D->decayFactor[0] = -1.0; //only for PET + else + hdr0.pixdim[4] = 0.0; // saveAs3D = true; // printWarning("Creating independent volumes as time between volumes varies\n"); if (opts.isVerbose) { @@ -4754,7 +5296,6 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc } } //if trVaries } //if PET - //next: detect variable inter-slice distance float dx = intersliceDistance(dcmList[dcmSort[0].indx],dcmList[dcmSort[1].indx]); #ifdef myInstanceNumberOrderIsNotSpatial @@ -4848,7 +5389,10 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc #endif int imageNumRange = 1 + abs( dcmList[dcmSort[nConvert-1].indx].imageNum - dcmList[dcmSort[0].indx].imageNum); if ((imageNumRange > 1) && (imageNumRange != nConvert)) { - printWarning("Missing images? Expected %d images, but instance number (0020,0013) ranges from %d to %d\n", nConvert, dcmList[dcmSort[0].indx].imageNum, dcmList[dcmSort[nConvert-1].indx].imageNum); + if ((dcmList[dcmSort[0].indx].locationsInAcquisition > 0) && ((nConvert % dcmList[dcmSort[0].indx].locationsInAcquisition) != 0) ) + printError("Missing images. Found %d images, expected %d slices per volume and instance number (0020,0013) ranges from %d to %d\n", nConvert, dcmList[dcmSort[0].indx].locationsInAcquisition, dcmList[dcmSort[0].indx].imageNum, dcmList[dcmSort[nConvert-1].indx].imageNum); + else + printWarning("Missing images? Expected %d images, but instance number (0020,0013) ranges from %d to %d\n", nConvert, dcmList[dcmSort[0].indx].imageNum, dcmList[dcmSort[nConvert-1].indx].imageNum); if (opts.isVerbose) { printMessage("instance=["); for (int i = 0; i < nConvert; i++) { @@ -4888,7 +5432,6 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc if (dcmList[dcmSort[i].indx].CSA.numDti > 0) dcmList[indx0].CSA.numDti =1; } - /*if (nConvert > 1) { //next determine if TR is true time between volumes double startTime = dcmList[indx0].acquisitionTime; double endTime = startTime; @@ -4933,7 +5476,7 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc if (hdr0.dim[4] > 1) //for 4d datasets, last volume should be acquired before first checkDateTimeOrder(&dcmList[dcmSort[0].indx], &dcmList[dcmSort[nConvert-1].indx]); } - int sliceDir = sliceTimingCore(dcmSort, dcmList, &hdr0, opts.isVerbose, nameList->str[dcmSort[0].indx], nConvert); + int sliceDir = sliceTimingCore(dcmSort, dcmList, &hdr0, opts.isVerbose, nameList->str[dcmSort[0].indx], nConvert, opts); #ifdef myReportSliceFilenames if (sliceDir < 0) { for (int i = nConvert; i > 0; --i) @@ -4955,7 +5498,7 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc return EXIT_FAILURE; } //issue377(dcmList[indx0], &hdr0); return EXIT_SUCCESS; - nii_SaveBIDS(pathoutname, dcmList[dcmSort[0].indx], opts, &hdr0, nameList->str[dcmSort[0].indx]); + nii_SaveBIDSX(pathoutname, dcmList[dcmSort[0].indx], opts, &hdr0, nameList->str[dcmSort[0].indx], dti4D); if (opts.isOnlyBIDS) { //note we waste time loading every image, however this ensures hdr0 matches actual output free(imgM); @@ -5011,7 +5554,10 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc printMessage(" %s\n",nameList->str[dcmSort[0].indx]); return EXIT_SUCCESS; } - if (sliceDir < 0) { + struct nifti_1_header hdrrx = hdr0; + bool isFlipZ = false; + if (sliceDir < 0) { + isFlipZ = true; imgM = nii_flipZ(imgM, &hdr0); sliceDir = abs(sliceDir); //change this, we have flipped the image so GE DTI bvecs no longer need to be flipped! } @@ -5047,6 +5593,10 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc nii_scale16bitUnsigned(imgM, &hdr0, opts.isVerbose); //allow UINT16 to use full dynamic range } else if ((opts.isMaximize16BitRange == kMaximize16BitRange_False) && (hdr0.datatype == DT_UINT16) && (!dcmList[dcmSort[0].indx].isSigned)) nii_check16bitUnsigned(imgM, &hdr0, opts.isVerbose); //save UINT16 as INT16 if we can do this losslessly + if ((dcmList[dcmSort[0].indx].isXA10A) && (nConvert < 2)) + printWarning("Siemens XA DICOM inadequate for robust conversion (issue 236)\n"); + if ((dcmList[dcmSort[0].indx].isXA10A) && (nConvert > 1)) + printWarning("Siemens XA exported as classic not enhanced DICOM (issue 236)\n"); printMessage( "Convert %d DICOM as %s (%dx%dx%dx%d)\n", nConvert, pathoutname, hdr0.dim[1],hdr0.dim[2],hdr0.dim[3],hdr0.dim[4]); #ifndef USING_R fflush(stdout); //show immediately if run from MRIcroGL GUI @@ -5060,11 +5610,16 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc else fflush(stdout); //GUI buffers printf, display all results #endif - if ((opts.isRotate3DAcq) && (dcmList[dcmSort[0].indx].is3DAcq) && (dcmList[dcmSort[0].indx].echoTrainLength <= 1) && (hdr0.dim[3] > 1) && (hdr0.dim[0] < 4)) - imgM = nii_setOrtho(imgM, &hdr0); //printMessage("ortho %d\n", echoInt (33)); - else if (opts.isFlipY)//(FLIP_Y) //(dcmList[indx0].CSA.mosaicSlices < 2) && + //3D-EPI vs 3D SPACE/MPRAGE/ETC + bool isFlipY = false; + bool isSetOrtho = false; + if ((opts.isRotate3DAcq) && (dcmList[dcmSort[0].indx].is3DAcq) && (!dcmList[dcmSort[0].indx].isEPI) && (hdr0.dim[3] > 1) && (hdr0.dim[0] < 4)) { + imgM = nii_setOrtho(imgM, &hdr0); + isSetOrtho = true; + } else if (opts.isFlipY){//(FLIP_Y) //(dcmList[indx0].CSA.mosaicSlices < 2) && imgM = nii_flipY(imgM, &hdr0); - else + isFlipY = true; + } else printMessage("DICOM row order preserved: may appear upside down in tools that ignore spatial transforms\n"); //begin: gantry tilt we need to save the shear in the transform mat44 sForm; @@ -5096,6 +5651,8 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc int returnCode = EXIT_FAILURE; #ifndef myNoSave // Indicates success or failure of the (last) save + if (opts.isSaveNRRD) + removeSclSlopeInter(&hdr0, imgM); //printMessage(" x--> %d ----\n", nConvert); if (! opts.isRGBplanar) //save RGB as packed RGBRGBRGB... instead of planar RRR..RGGG..GBBB..B imgM = nii_planar2rgb(imgM, &hdr0, true); //NIfTI is packed while Analyze was planar @@ -5116,6 +5673,52 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc else nii_saveNII(pathoutnameADC, hdr0, imgM, opts, dcmList[dcmSort[0].indx]); } + if (isHasOverlay) { //each series can have up to 16 overlays, overlays may not be on all slices + for (int j = 0; j < kMaxOverlay; j++) { + bool isOverlay = false; + for (int i = 0; i < nConvert; i++) + if (dcmList[dcmSort[i].indx].overlayStart[j] > 0) isOverlay = true; + if (!isOverlay) continue; + char pathoutnameROI[2048] = {""}; + strcat(pathoutnameROI,pathoutname); + char append[128] = {""}; + sprintf(append,"_ROI%d",j+1); + strcat(pathoutnameROI,append); + struct nifti_1_header hdrr = hdrrx; + hdrr.dim[0] = 3; + if (hdrr.dim[1] < 1) hdrr.dim[1] = 1; + if (hdrr.dim[2] < 1) hdrr.dim[2] = 1; + if (hdrr.dim[3] < 1) hdrr.dim[3] = 1; + hdrr.dim[4] = 1; + hdrr.bitpix = 8; + hdrr.datatype = 2; + hdrr.scl_inter = 0.0; + hdrr.scl_slope = 1.0; + int nvox = hdrr.dim[1] * hdrr.dim[2] * hdrr.dim[3]; + unsigned char *imgR = (unsigned char *)malloc(nvox); + for (int v = 0; v < nvox; v++) + imgR[v] = 0; + if (nConvert == 1) { + int indx = dcmSort[0].indx; + loadOverlay(nameList->str[indx], imgR, dcmList[indx].overlayStart[j], hdrr.dim[1], hdrr.dim[2], hdrr.dim[3]); + } else if (nConvert == hdrr.dim[3]) { + for (int i = 0; i < nConvert; i++) { + int indx = dcmSort[i].indx; + if (dcmList[indx].overlayStart[j] > 0) { + unsigned char *imgRS = imgR + (i * hdrr.dim[1] * hdrr.dim[2]); + loadOverlay(nameList->str[indx], imgRS, dcmList[indx].overlayStart[j], hdrr.dim[1], hdrr.dim[2], 1); + } //if overlay on slice + } //for each volume + } // + if (isFlipZ) + imgR = nii_flipZ(imgR, &hdrr); + if (isSetOrtho) + imgR = nii_setOrtho(imgR, &hdrr); + if (isFlipY) + imgR = nii_flipY(imgR, &hdrr); + nii_saveNII(pathoutnameROI, hdrr, imgR, opts, dcmList[dcmSort[0].indx]); + } + } #endif imgM = removeADC(&hdr0, imgM, numADC); #ifndef USING_R @@ -5147,7 +5750,8 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc returnCode = nii_saveNII3Deq(pathoutname, hdr0, imgM,opts, dcmList[dcmSort[0].indx], sliceMMarray); free(sliceMMarray); } - if ((opts.isRotate3DAcq) && (opts.isCrop) && (dcmList[indx0].is3DAcq) && (hdr0.dim[3] > 1) && (hdr0.dim[0] < 4))//for T1 scan: && (dcmList[indx0].TE < 25) + //3D-EPI vs 3D SPACE/MPRAGE/ETC + if ((opts.isRotate3DAcq) && (opts.isCrop) && (dcmList[indx0].is3DAcq) && (!dcmList[indx0].isEPI) && (hdr0.dim[3] > 1) && (hdr0.dim[0] < 4))//for T1 scan: && (dcmList[indx0].TE < 25) returnCode = nii_saveCrop(pathoutname, hdr0, imgM, opts, dcmList[dcmSort[0].indx]); //n.b. must be run AFTER nii_setOrtho()! #ifdef USING_R // Note that for R, only one image should be created per series @@ -5270,11 +5874,10 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis dcmList[indx].RWVScale = dti4Ds.RWVScale[0]; } //if isScaleVariesEnh - if (dcmList[indx].isScaleVariesEnh) printf("Varies\n"); - if (!dcmList[indx].isScaleVariesEnh) printf("no Varies\n"); - printf("%g %g %g\n", dcmList[indx].intenScale, dcmList[indx].intenIntercept, dcmList[indx].intenScalePhilips); + //if (dcmList[indx].isScaleVariesEnh) printf("Varies\n"); + //if (!dcmList[indx].isScaleVariesEnh) printf("no Varies\n"); + //printf("%g %g %g\n", dcmList[indx].intenScale, dcmList[indx].intenIntercept, dcmList[indx].intenScalePhilips); //continue; - if (s > 1) dcmList[indx].CSA.numDti = 0; //only save bvec for first type (magnitude) int ret2 = saveDcm2NiiCore(nConvert, dcmSort, dcmList, nameList, opts, &dti4Ds, s); if (ret2 != EXIT_SUCCESS) ret = ret2; //return EXIT_SUCCESS only if ALL are successful @@ -5414,7 +6017,7 @@ bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, struct TDCMopts* opt bool isForceStackSeries = false; if ((opts->isForceStackDCE) && (d1.isStackableSeries) && (d2.isStackableSeries) && (d1.seriesNum != d2.seriesNum)) { if (!warnings->forceStackSeries) - printMessage("DCE volumes stacked despite varying series number (use '-m o' to turn off merging).\n"); + printMessage("Siemens volumes stacked despite varying series number (use '-m o' to turn off merging).\n"); warnings->forceStackSeries = true; isForceStackSeries = true; } @@ -5462,14 +6065,14 @@ bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, struct TDCMopts* opt #endif if ((opts->isForceStackSameSeries == 1) || ((opts->isForceStackSameSeries == 2) && (d1.isXRay) )) { // "isForceStackSameSeries == 2" will automatically stack CT scans but not MR - //if ((d1.TE != d2.TE) || (d1.echoNum != d2.echoNum)) - if ((!(isSameFloat(d1.TE, d2.TE))) || (d1.echoNum != d2.echoNum)) - *isMultiEcho = true; + //if (((!(isSameFloat(d1.TE, d2.TE))) || (d1.echoNum != d2.echoNum)) && (!d1.isXRay)) { + // *isMultiEcho = true; + //} return true; //we will stack these images, even if they differ in the following attributes } if ((d1.isHasImaginary != d2.isHasImaginary) || (d1.isHasPhase != d2.isHasPhase) || ((d1.isHasReal != d2.isHasReal))) { if (!warnings->phaseVaries) - printMessage("Slices not stacked: some are phase/real/imaginary maps, others are not. Use 'f 2D slices' option to force stacking\n"); + printMessage("Slices not stacked: some are phase/real/imaginary maps, others are not. Instances %d %d\n", d1.imageNum, d2.imageNum); warnings->phaseVaries = true; return false; } @@ -5508,7 +6111,7 @@ bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, struct TDCMopts* opt } if ((!isSameFloatGE(d1.orient[1], d2.orient[1]) || !isSameFloatGE(d1.orient[2], d2.orient[2]) || !isSameFloatGE(d1.orient[3], d2.orient[3]) || !isSameFloatGE(d1.orient[4], d2.orient[4]) || !isSameFloatGE(d1.orient[5], d2.orient[5]) || !isSameFloatGE(d1.orient[6], d2.orient[6]) ) ) { - if ((!warnings->orientVaries) && (!d1.isNonParallelSlices)) + if ((!warnings->orientVaries) && (!d1.isNonParallelSlices) && (!d1.isLocalizer)) printMessage("Slices not stacked: orientation varies (vNav or localizer?) [%g %g %g %g %g %g] != [%g %g %g %g %g %g]\n", d1.orient[1], d1.orient[2], d1.orient[3],d1.orient[4], d1.orient[5], d1.orient[6], d2.orient[1], d2.orient[2], d2.orient[3],d2.orient[4], d2.orient[5], d2.orient[6]); @@ -6045,6 +6648,7 @@ int nii_loadDirCore(char *indir, struct TDCMopts* opts) { continue; } dcmList[i] = readDICOMv(nameList.str[i], opts->isVerbose, opts->compressFlag, &dti4D); //ignore compile warning - memory only freed on first of 2 passes + //if (!dcmList[i].isValid) printf(">>>>Not a valid DICOM %s\n", nameList.str[i]); if ((dcmList[i].isValid) && ((dti4D.sliceOrder[0] >= 0) || (dcmList[i].CSA.numDti > 1))) { //4D dataset: dti4D arrays require huge amounts of RAM - write this immediately struct TDCMsort dcmSort[1]; fillTDCMsort(dcmSort[0], i, dcmList[i]); @@ -6290,8 +6894,7 @@ int nii_loadDir(struct TDCMopts* opts) { if (!is_dir(opts->indir,true)) { printError("Input folder invalid: %s\n",opts->indir); return kEXIT_INPUT_FOLDER_INVALID; - } - + } #ifdef USING_R // Full file paths are only used by R/divest when reorganising DICOM files if (opts->isRenameNotConvert) { @@ -6561,12 +7164,14 @@ void setDefaultOpts (struct TDCMopts *opts, const char * argv[]) { //either "set opts->isForceStackSameSeries = 2; //automatic: stack CTs, do not stack MRI opts->isForceStackDCE = true; opts->isIgnoreDerivedAnd2D = false; + opts->isForceOnsetTimes = true; opts->isPhilipsFloatNotDisplayScaling = true; opts->isCrop = false; opts->isRotate3DAcq = true; opts->isGz = false; opts->isSaveNativeEndian = true; opts->isAddNamePostFixes = true; //e.g. "_e2" added for second echo + opts->isTestx0021x105E = false; //GE test slice times stored in 0021,105E opts->isSaveNRRD = false; opts->isPipedGz = false; //e.g. pipe data directly to pigz instead of saving uncompressed to disk opts->isSave3D = false; diff --git a/console/nii_dicom_batch.h b/console/nii_dicom_batch.h index 59fecf72..e8cad56e 100644 --- a/console/nii_dicom_batch.h +++ b/console/nii_dicom_batch.h @@ -35,7 +35,7 @@ extern "C" { #define MAX_NUM_SERIES 16 struct TDCMopts { - bool isAddNamePostFixes, isSaveNativeEndian, isSaveNRRD, isOneDirAtATime, isRenameNotConvert, isSave3D, isGz, isPipedGz, isFlipY, isCreateBIDS, isSortDTIbyBVal, isAnonymizeBIDS, isOnlyBIDS, isCreateText, isIgnoreDerivedAnd2D, isPhilipsFloatNotDisplayScaling, isTiltCorrect, isRGBplanar, isOnlySingleFile, isForceStackDCE, isRotate3DAcq, isCrop; + bool isTestx0021x105E, isAddNamePostFixes, isSaveNativeEndian, isSaveNRRD, isOneDirAtATime, isRenameNotConvert, isSave3D, isGz, isPipedGz, isFlipY, isCreateBIDS, isSortDTIbyBVal, isAnonymizeBIDS, isOnlyBIDS, isCreateText, isForceOnsetTimes,isIgnoreDerivedAnd2D, isPhilipsFloatNotDisplayScaling, isTiltCorrect, isRGBplanar, isOnlySingleFile, isForceStackDCE, isRotate3DAcq, isCrop; int isMaximize16BitRange, isForceStackSameSeries, nameConflictBehavior, isVerbose, isProgress, compressFlag, dirSearchDepth, gzLevel; //support for compressed data 0=none, char filename[512], outdir[512], indir[512], pigzname[512], optsname[512], indirParent[512], imageComments[24]; double seriesNumber[MAX_NUM_SERIES]; //requires double must store -1 (report but do not convert) as well as seriesUidCrc (uint32) diff --git a/console/nii_foreign.h b/console/nii_foreign.h index cf4bc466..9249a4c9 100644 --- a/console/nii_foreign.h +++ b/console/nii_foreign.h @@ -16,4 +16,4 @@ int convert_foreign (const char *fn, struct TDCMopts opts); } #endif -#endif \ No newline at end of file +#endif diff --git a/console/tinydir.h b/console/tinydir.h index cc65b7f7..82d41837 100644 --- a/console/tinydir.h +++ b/console/tinydir.h @@ -347,7 +347,9 @@ int tinydir_readfile(const tinydir_dir *dir, tinydir_file *file) dir->_e->d_name #endif ); - strcat(file->path, file->name); + /* Limit the number of bytes copied to the maximum length of the name, + to avoid spurious compiler warnings about possible overlap */ + strncat(file->path, file->name, _TINYDIR_FILENAME_MAX); #ifndef _MSC_VER if (stat(file->path, &file->_s) == -1) { diff --git a/dcm_qa b/dcm_qa index cb8f40e9..2f59a065 160000 --- a/dcm_qa +++ b/dcm_qa @@ -1 +1 @@ -Subproject commit cb8f40e991bd75db8b0d89611fc76ec4805499c6 +Subproject commit 2f59a065686f1bdab52f8c4c529e1e2e7e4c133f diff --git a/dcm_qa_nih b/dcm_qa_nih index 9460104a..a6fffa39 160000 --- a/dcm_qa_nih +++ b/dcm_qa_nih @@ -1 +1 @@ -Subproject commit 9460104acfeed6acf7c8139f530e785f8e64e96c +Subproject commit a6fffa392766b7bd816fe0e822e6469f044af9ee diff --git a/dcm_qa_uih b/dcm_qa_uih index d1329e1b..78f21294 160000 --- a/dcm_qa_uih +++ b/dcm_qa_uih @@ -1 +1 @@ -Subproject commit d1329e1b4bc151d7c05c03197e768b937c456843 +Subproject commit 78f212941ca7adc13bd3810667084b48a0047484