diff --git a/.travis.yml b/.travis.yml index 6ecc20d5..80b8b2bb 100644 --- a/.travis.yml +++ b/.travis.yml @@ -34,8 +34,7 @@ script: - cd dcm_qa_uih && ./batch.sh && cd - before_deploy: - - export DATE=`date +%-d-%b-%Y` - - zip -j dcm2niix_${DATE}_${TARGET}.zip build/bin/* + - zip -j dcm2niix_${TARGET}.zip build/bin/* - sleep 300 # make sure appveyor deployment is done, thus proper release name is set deploy: diff --git a/PARREC/README.md b/PARREC/README.md index 78304b0d..fba3e5ee 100644 --- a/PARREC/README.md +++ b/PARREC/README.md @@ -5,7 +5,7 @@ dcm2niix attempts to convert Philips PAR/REC format images to NIfTI. While this According to [Matthew Clemence](https://www.nitrc.org/forum/forum.php?thread_id=9319&forum_id=4703) DICOM (classic and enhanced) and XML/REC are supported in the base product, NIFTI forms part of a Neuroscience commercial option from release 5 onwards. PAR/REC requires a research agreement to obtain. For the two formats XML/REC and PAR/REC, the "REC" part is identical but instead of a plain text file of the "par" format, the same information is now available as an XML file. This descision has been taken to allow the information to be more easily extended as the PAR file was getting increasingly limited. -## Detecting, Reporting and Fixing the V4 Image Offcente Bug +## Detecting, Reporting and Fixing the V4 Image offcentre Bug The PAR header contains a field 'image offcentre (ap,fh,rl in mm )' that we use to detect the spatial position of slices (e.g. for an axial scan is the first slice inferior or superior to the final slice). However, it appears that in some V4 images the values in these columns are actually stored in the order "rl,ap,fh". This has never been reported in V3, V4.1 and V4.2 images. A nice example of this is the ['philips_1_5T_intera' dataset provided with Rosetta Bit](https://www.nitrc.org/projects/rosetta/)(actually from a 3T MRI). This sample includes both DICOM and V4 PAR/REC data. Note the 'Off Centre midslice(ap,fh,rl) [mm]' field gives the volume center in the correct order. However, the subsequent 'image offcentre' fields are swizzled. The latest versions of dcm2niix will detect, report and correct this error. If you do see an error like the one below, please report it on Github as an issue, so we can have a better understanding of its prevalence. diff --git a/Philips/README.md b/Philips/README.md index 0616acbf..7faf6902 100644 --- a/Philips/README.md +++ b/Philips/README.md @@ -16,6 +16,23 @@ In general it is recommended that you archive and convert DICOM images as they a Therefore, dcm2niix will ignore the IPP enclosed in 2005,140F unless no alternative exists. +## Image Scaling + +dcm2niix losslessy copies the raw data from DICOM to NIfTI format. These values are typically stored as 16-bit integers in the range -32768..32767. Both the DICOM and NIfTI formats describe how scaling intercept and slope values can be used to convert these raw values into calibrated values. For example, with an intercept of 0 and slope of 0.01 the raw value of 50 would be converted to 0.5. + +Unlike other vendors, Philips can store different scaling factors in their DICOM header. For most MRI modalities where the intensity brightness is relative, this has no impact. However, for modalities like ASL it can have an impact. The NIfTI format requires a single intensity intercept and slope is chosen. Therefore, dcm2niix will choose the "Real World" values if provided. If these are not available, dcm2niix will choose either the "precise" (default) or "display" (if the user choose "-p n") value. dcm2niix will also populate the folllowing tags in the BIDS header that allow the user to select between different intensity scaling formats: "PhilipsRescaleSlope", "PhilipsRescaleIntercept", "PhilipsScaleSlope", "UsePhilipsFloatNotDisplayScaling" (where "1" indicates NIfTI uses precise value, and "0" indicates display values)., "PhilipsRWVSlope" and "PhilipsRWVIntercept". + +The relevant DICOM tags are +RS = rescale slope ([0028,1053](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0028,1053))) +RI = rescale intercept ([0028,1052](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0028,1052))) +SS = scale slope (2005,100E) +RealWorldIntercept = (0040,9224) +Real World Slope = (0040,9225) +The transformation formulas are: +R = raw value, P = precise value, D = displayed value +D = R * RS + RI +P = D/(RS * SS) + ## Derived parametric maps stored with raw diffusion data Some Philips diffusion DICOM images include derived image(s) along with the images. Other manufacturers save these derived images as a separate series number, and the DICOM standard seems ambiguous on whether it is allowable to mix raw and derived data in the same series (see PS 3.3-2008, C.7.6.1.1.2-3). In practice, many Philips diffusion images append [derived parametric maps](http://www.revisemri.com/blog/2008/diffusion-tensor-imaging/) with the original data. With Philips, appending the derived isotropic image is optional - it is only created for the 'clinical' DTI schemes for radiography analysis and is triggered if the first three vectors in the gradient table are the unit X,Y and Z vectors. For conventional DWI, the result is the conventional mean of the ADC X,Y,Z for DTI it the conventional mean of the 3 principle Eigen vectors. As scientists, we want to discard these derived images, as they will disrupt data processing and we can generate better parametric maps after we have applied undistortion methods such as [Eddy and Topup](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/eddy/UsersGuide). The current version of dcm2niix uses the Diffusion Directionality (0018,9075) tag to detect B=0 unweighted ("NONE"), B-weighted ("DIRECTIONAL"), and derived ("ISOTROPIC") images. Note that the Dimension Index Values (0020,9157) tag provides an alternative approach to discriminate these images. Here are sample tags from a Philips enhanced image that includes and derived map (3rd dimension is "1" while the other images set this to "2"). @@ -65,6 +82,20 @@ Likewise, the BIDS tag "PhaseEncodingDirection" allows tools like [eddy](https:/ Another value desirable for TOPUP is the "TotalReadoutTime". Again, one can not calculate this from Philips DICOMs. If you do decide to calculate this using values from the MRI console, be aware that the [FSL definition](https://github.com/rordenlab/dcm2niix/issues/130) is not intuitive for scans with interpolation, partial Fourier, parallel imaging, etc. However, it should be pointed out that the "TotalReadoutTime" only inlfuences TOPUP's calibrated validation images that are typically ignored. The data used in subsequent steps will not be influenced by this value. +## Non-Image DICOMs + +NIfTI is an image format, while DICOM is a multi-purpose format that can store videos (MPEG) or other data. Specifically, some Philips systems save Exam Cards and other non-image data as DICOM files. In these case, dcm2niix should skip these files, as they can not be represented in NIfTI. You can discriminate these files by reading the [MediaStorageSOPClassUID (0002,0002)](https://github.com/rordenlab/dcm2niix/issues/328). + +- MR Image Storage = 1.2.840.10008.5.1.4.1.1.4 +- Enhanced MR Image Storage = 1.2.840.10008.5.1.4.1.1.4.1 +- MR Spectroscopy Storage = 1.2.840.10008.5.1.4.1.1.4.2 +- Secondary Capture Image Storage = 1.2.840.10008.5.1.4.1.1.7 +- Grayscale Softcopy Presentation State = 1.2.840.10008.5.1.4.1.1.11.1 +- Raw Data Storage = 1.2.840.10008.5.1.4.1.1.66 +- (Old) Private MR Spectrum Storage = 1.3.46.670589.11.0.0.12.1 +- (Old) Private MR Series Data Storage = 1.3.46.670589.11.0.0.12.2 +- (Old) Private MR Examcard Storage = 1.3.46.670589.11.0.0.12.4 + ## General variations Prior versions of dcm2niix used different methods to sort images. However, these have proved unreliable The undocumented tags SliceNumberMrPhilips (2001,100A). In theory, InStackPositionNumber (0020,9057) should be present in all enhanced files, but has not proved reliable (perhaps not in older Philips images or DICOM images that were modified after leaving the scanner). MRImageGradientOrientationNumber (2005,1413) is complicated by the inclusion of derived images. Therefore, current versions of dcm2niix do not generally depend on any of these. diff --git a/README.md b/README.md index b24c9bb0..9d81925b 100644 --- a/README.md +++ b/README.md @@ -37,7 +37,8 @@ Command line usage is described in the [NITRC wiki](https://www.nitrc.org/plugin 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. - - [MRIcroGL](https://github.com/neurolabusc/MRIcroGL/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. + - Assuming the program `curl` is installed on your computer (e.g. Windows 10 build 1803 or later), you can download the latest Github release with the command. This allows you to download the latest stable release with a single line of code. For Linux (note prior comment regarding older distributions): `curl -fLO https://github.com/rordenlab/dcm2niix/releases/latest/download/dcm2niix_lnx.zip`. For MacOS: `curl -fLO https://github.com/rordenlab/dcm2niix/releases/latest/download/dcm2niix_mac.zip`. For Windows: `curl -fLO https://github.com/rordenlab/dcm2niix/releases/latest/download/dcm2niix_win.zip`. + - [MRIcroGL (stable)](https://www.nitrc.org/projects/mricrogl) or [MRIcroGL (development)](https://github.com/neurolabusc/MRIcroGL/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. - 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`. @@ -88,11 +89,11 @@ If you have any problems with the cmake build script described above or want to ## Alternatives - - [dcm2nii](http://www.mccauslandcenter.sc.edu/mricro/mricron/dcm2nii.htm) 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). + - [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). - [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/). - - [dinifti](http://cbi.nyu.edu/software/dinifti.php) is focused on conversion of Siemens data. + - [dinifti](http://as.nyu.edu/cbi/resources/Software/DINIfTI.html) is focused on conversion of Siemens data. - [DWIConvert](https://github.com/BRAINSia/BRAINSTools/tree/master/DWIConvert) converts DICOM images to NRRD and NIfTI formats. - [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. diff --git a/Siemens/README.md b/Siemens/README.md index b857f505..8f5918d9 100644 --- a/Siemens/README.md +++ b/Siemens/README.md @@ -30,7 +30,14 @@ Many crucial Siemens parameters are stored in the [proprietary CSA header](http: ## Slice Timing -For Siemens images created with software versions B15 through E11, the proprietary [CSA Image Header (0029,1010)](https://nipy.org/nibabel/dicom/siemens_csa.html) contains the array MosaicRefAcqTimes that provides [slice timing](https://www.mccauslandcenter.sc.edu/crnl/tools/stc). Earlier Siemens 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. Therefore dcm2niix typically is able to provide accurate slice timing information for non-Vida datasets (the prior section describes Vida slice timing issues seen with the XA software series). Some Siemens DICOMs stroe slice timings in the private tag [0019,1029](https://github.com/rordenlab/dcm2niix/issues/296). In theory, this could be used when the CSA header is missing. For archival studies, be aware that some sequences [incorrectly reported slice timing](https://github.com/rordenlab/dcm2niix/issues/126). The [SPM slice timing wiki](https://en.wikibooks.org/w/index.php?title=SPM/Slice_Timing&stable=0#Siemens_scanners) provides further information on Siemens 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. + +For software versions B15 through E11 where all slices of a volume are stored as a single mosaic file, the proprietary [CSA Image Header (0029,1010)](https://nipy.org/nibabel/dicom/siemens_csa.html) contains the array MosaicRefAcqTimes that provides [slice timing](https://www.mccauslandcenter.sc.edu/crnl/tools/stc). For volumes where each 2D slice is saved as a separate DICOM file, one can infer slice order from the DICOM tag Acquisition Time (0008,0032). + + The prior section describes Vida slice timing issues seen with the XA software series. In brief, dcm2niix will use Frame Acquisition Time (0018,9074) to determine slice times. Some Siemens DICOMs store slice timings in the private tag [0019,1029](https://github.com/rordenlab/dcm2niix/issues/296). In theory, this could be used when the CSA header is missing. For archival studies, be aware that some sequences [incorrectly reported slice timing](https://github.com/rordenlab/dcm2niix/issues/126). The [SPM slice timing wiki](https://en.wikibooks.org/w/index.php?title=SPM/Slice_Timing&stable=0#Siemens_scanners) provides further information on Siemens slice timing. ## Total Readout Time @@ -40,6 +47,58 @@ One often wants to determine [echo spacing, bandwidth, ](https://support.brainvo 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. +## 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. + +ep2d_pcasl, ep2d_pcasl_UI_PHC //pCASL 2D [Danny J.J. Wang](http://www.loft-lab.org) + - LabelOffset + - PostLabelDelay + - NumRFBlocks + - RFGap + - MeanGzx10 + - PhiAdjust + +tgse_pcasl //pCASL 3D [Danny J.J. Wang](http://www.loft-lab.org) + - RFGap + - MeanGzx10 + - T1 + +ep2d_pasl //PASL 2D Siemens Product + - InversionTime + - SaturationStopTime + +tgse_pasl //PASL 3D [Siemens Product](http://adni.loni.usc.edu/wp-content/uploads/2010/05/ADNI3_Basic_Siemens_Skyra_E11.pdf) + - BolusDuration + - InversionTime + +ep2d_fairest //PASL 2D http://www.pubmed.com/11746944 http://www.pubmed.com/21606572 + - PostInversionDelay + - PostLabelDelay + +to_ep2d_VEPCASL //pCASL 2D specific tags - Oxford (Thomas OKell) + - InversionTime + - BolusDuration + - TagRFFlipAngle + - TagRFDuration + - TagRFSeparation + - MeanTagGradient + - TagGradientAmplitude + - TagDuration + - MaximumT1Opt + - InitialPostLabelDelay [Array] + +jw_tgse_VEPCASL //pCASL 3D Oxford + - TagRFFlipAngle + - TagRFDuration + - TagRFSeparation + - MaximumT1Opt + - Tag0 + - Tag1 + - Tag2 + - Tag3 + - InitialPostLabelDelay [Array] + ## Sample Datasets - [Slice timing dataset](httphttps://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage#Slice_timing_corrections://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage). diff --git a/appveyor.yml b/appveyor.yml index e3c1863e..9f4a5d0f 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -45,7 +45,7 @@ after_build: - ps: >- cd c:\projects\dcm2niix - 7z a dcm2niix_$($env:DATE)_win.zip c:\projects\dcm2niix\build\bin\* >$null + 7z a dcm2niix_win.zip c:\projects\dcm2niix\build\bin\* >$null artifacts: - path: dcm2niix*.zip diff --git a/console/jpg_0XC3.cpp b/console/jpg_0XC3.cpp old mode 100755 new mode 100644 diff --git a/console/main_console.cpp b/console/main_console.cpp index 9dd02cab..9d7ef5dc 100644 --- a/console/main_console.cpp +++ b/console/main_console.cpp @@ -40,6 +40,7 @@ #include #include "nii_dicom_batch.h" #include "nii_dicom.h" +#include "nifti1_io_core.h" #include #if !defined(_WIN64) && !defined(_WIN32) @@ -95,7 +96,7 @@ void showHelp(const char * argv[], struct TDCMopts opts) { if (opts.isMaximize16BitRange) max16Ch = 'y'; printf(" -l : losslessly scale 16-bit integers to use dynamic range (y/n, default %c)\n", max16Ch); printf(" -m : merge 2D slices from same series regardless of echo, exposure, etc. (n/y or 0/1/2, default 2) [no, yes, auto]\n"); - printf(" -n : only convert this series number - can be used up to %i times (default convert all)\n", MAX_NUM_SERIES); + printf(" -n : only convert this series CRC number - can be used up to %i times (default convert all)\n", MAX_NUM_SERIES); printf(" -o : output directory (omit to save to input folder)\n"); printf(" -p : Philips precise float (not display) scaling (y/n, default y)\n"); printf(" -r : rename instead of convert DICOMs (y/n, default n)\n"); @@ -111,7 +112,7 @@ void showHelp(const char * argv[], struct TDCMopts opts) { //#define kNAME_CONFLICT_ADD_SUFFIX 2 //default 2 = write with new suffix as a new file printf(" -w : write behavior for name conflicts (0,1,2, default 2: 0=skip duplicates, 1=overwrite, 2=add suffix)\n"); printf(" -x : crop 3D acquisitions (y/n/i, default n, use 'i'gnore to neither crop nor rotate 3D acquistions)\n"); - char gzCh = 'n'; + char gzCh = 'n'; if (opts.isGz) gzCh = 'y'; #if defined(_WIN64) || defined(_WIN32) //n.b. the optimal use of pigz requires pipes that are not provided for Windows @@ -127,6 +128,10 @@ void showHelp(const char * argv[], struct TDCMopts opts) { printf(" -z : gz compress images (y/i/n/3, default %c) [y=pigz, i=internal:miniz, n=no, 3=no,3D]\n", gzCh); #endif #endif + printf(" --big-endian : byte order (y/n/o, default o) [y=big-end, n=little-end, o=optimal/native]\n"); + printf(" --progress : Slicer format progress information\n"); + printf(" --version : report version\n"); + printf(" --xml : Slicer format features\n"); printf(" Defaults stored in Windows registry\n"); printf(" Examples :\n"); printf(" %s c:\\DICOM\\dir\n", cstr); @@ -147,6 +152,10 @@ void showHelp(const char * argv[], struct TDCMopts opts) { printf(" -z : gz compress images (y/o/i/n/3, default %c) [y=pigz, o=optimal pigz, i=internal:miniz, n=no, 3=no,3D]\n", gzCh); #endif #endif + printf(" --big-endian : byte order (y/n/o, default o) [y=big-end, n=little-end, o=optimal/native]\n"); + printf(" --progress : Slicer format progress information\n"); + printf(" --version : report version\n"); + printf(" --xml : Slicer format features\n"); printf(" Defaults file : %s\n", opts.optsname); printf(" Examples :\n"); printf(" %s /Users/chris/dir\n", cstr); @@ -264,14 +273,27 @@ int main(int argc, const char * argv[]) if ((strlen(argv[i]) > 1) && (argv[i][0] == '-')) { //command if (argv[i][1] == 'h') { showHelp(argv, opts); - } else if (( ! strcmp(argv[1], "--progress")) && ((i+1) < argc)) { + } else if (( ! strcmp(argv[i], "--big-endian")) && ((i+1) < argc)) { + i++; + if ((littleEndianPlatform()) && ((argv[i][0] == 'y') || (argv[i][0] == 'Y'))) { + opts.isSaveNativeEndian = false; + printf("NIfTI data will be big-endian (byte-swapped)\n"); + } + if ((!littleEndianPlatform()) && ((argv[i][0] == 'n') || (argv[i][0] == 'N'))) { + opts.isSaveNativeEndian = false; + printf("NIfTI data will be little-endian\n"); + } + } else if ( ! strcmp(argv[i], "--version")) { + printf("%s\n", kDCMdate); + return kEXIT_REPORT_VERSION; + } else if (( ! strcmp(argv[i], "--progress")) && ((i+1) < argc)) { i++; if ((argv[i][0] == 'n') || (argv[i][0] == 'N') || (argv[i][0] == '0')) opts.isProgress = 0; else opts.isProgress = 1; if (argv[i][0] == '2') opts.isProgress = 2; //logorrheic - } else if ( ! strcmp(argv[1], "--xml")) { + } else if ( ! strcmp(argv[i], "--xml")) { showXML(); return EXIT_SUCCESS; } else if ((argv[i][1] == 'a') && ((i+1) < argc)) { //adjacent DICOMs @@ -281,7 +303,6 @@ int main(int argc, const char * argv[]) opts.isOneDirAtATime = false; else opts.isOneDirAtATime = true; - } else if ((argv[i][1] >= '1') && (argv[i][1] <= '9')) { opts.gzLevel = abs((int)strtol(argv[i], NULL, 10)); if (opts.gzLevel > 11) @@ -468,9 +489,9 @@ int main(int argc, const char * argv[]) strcpy(opts.outdir,argv[i]); } else if ((argv[i][1] == 'n') && ((i+1) < argc)) { i++; - float seriesNumber = atof(argv[i]); + double seriesNumber = atof(argv[i]); if (seriesNumber < 0) - opts.numSeries = -1; //report series: convert none + opts.numSeries = -1.0; //report series: convert none else if ((opts.numSeries >= 0) && (opts.numSeries < MAX_NUM_SERIES)) { opts.seriesNumber[opts.numSeries] = seriesNumber; opts.numSeries += 1; diff --git a/console/nifti1_io_core.cpp b/console/nifti1_io_core.cpp index f100c683..4320dd20 100644 --- a/console/nifti1_io_core.cpp +++ b/console/nifti1_io_core.cpp @@ -17,6 +17,7 @@ #include //requires VS 2015 or later #include "nifti1_io_core.h" +#include "nifti1.h" #include #include #include @@ -31,6 +32,7 @@ #include "print.h" + #ifndef USING_R void nifti_swap_8bytes( size_t n , void *ar ) // 4 bytes at a time { @@ -80,6 +82,68 @@ void nifti_swap_2bytes( size_t n , void *ar ) // 2 bytes at a time } #endif +/*-------------------------------------------------------------------------*/ +/*! Byte swap NIFTI-1 file header in various places and ways. + + If is_nifti, swap all (even UNUSED) fields of NIfTI header. + Else, swap as a nifti_analyze75 struct. +*//*---------------------------------------------------------------------- */ +void swap_nifti_header( struct nifti_1_header *h ) +{ + + /* otherwise, swap all NIFTI fields */ + + nifti_swap_4bytes(1, &h->sizeof_hdr); + nifti_swap_4bytes(1, &h->extents); + nifti_swap_2bytes(1, &h->session_error); + + nifti_swap_2bytes(8, h->dim); + nifti_swap_4bytes(1, &h->intent_p1); + nifti_swap_4bytes(1, &h->intent_p2); + nifti_swap_4bytes(1, &h->intent_p3); + + nifti_swap_2bytes(1, &h->intent_code); + nifti_swap_2bytes(1, &h->datatype); + nifti_swap_2bytes(1, &h->bitpix); + nifti_swap_2bytes(1, &h->slice_start); + + nifti_swap_4bytes(8, h->pixdim); + + nifti_swap_4bytes(1, &h->vox_offset); + nifti_swap_4bytes(1, &h->scl_slope); + nifti_swap_4bytes(1, &h->scl_inter); + nifti_swap_2bytes(1, &h->slice_end); + + nifti_swap_4bytes(1, &h->cal_max); + nifti_swap_4bytes(1, &h->cal_min); + nifti_swap_4bytes(1, &h->slice_duration); + nifti_swap_4bytes(1, &h->toffset); + nifti_swap_4bytes(1, &h->glmax); + nifti_swap_4bytes(1, &h->glmin); + + nifti_swap_2bytes(1, &h->qform_code); + nifti_swap_2bytes(1, &h->sform_code); + + nifti_swap_4bytes(1, &h->quatern_b); + nifti_swap_4bytes(1, &h->quatern_c); + nifti_swap_4bytes(1, &h->quatern_d); + nifti_swap_4bytes(1, &h->qoffset_x); + nifti_swap_4bytes(1, &h->qoffset_y); + nifti_swap_4bytes(1, &h->qoffset_z); + + nifti_swap_4bytes(4, h->srow_x); + nifti_swap_4bytes(4, h->srow_y); + nifti_swap_4bytes(4, h->srow_z); + + return ; +} + +bool littleEndianPlatform () +{ + uint32_t value = 1; + return (*((char *) &value) == 1); +} + int isSameFloat (float a, float b) { return (fabs (a - b) <= FLT_EPSILON); } diff --git a/console/nifti1_io_core.h b/console/nifti1_io_core.h index 1e7afd62..ff371ad5 100644 --- a/console/nifti1_io_core.h +++ b/console/nifti1_io_core.h @@ -8,6 +8,8 @@ #define STRICT_R_HEADERS #include "RNifti.h" #endif +#include "nifti1.h" +#include #ifdef __cplusplus extern "C" { @@ -52,6 +54,8 @@ float dotProduct(vec3 u, vec3 v); float nifti_mat33_determ( mat33 R ) ; int isSameFloat (float a, float b) ; int isSameDouble (double a, double b) ; +bool littleEndianPlatform (); + vec3 nifti_mat33_eig3(double bxx, double bxy, double bxz, double byy, double byz, double bzz); mat33 nifti_mat33_inverse( mat33 R ); @@ -66,6 +70,7 @@ 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); +void swap_nifti_header ( struct nifti_1_header *h) ; 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 diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 49777fca..456824ef 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -964,12 +964,7 @@ void dcmStr(int lLength, unsigned char lBuffer[], char* lOut, bool isStrLarge = //#endif } //dcmStr() -inline bool littleEndianPlatform () -{ - uint32_t value = 1; - return (*((char *) &value) == 1); -} - +#ifdef MY_OLD //this code works on Intel but not some older systems https://github.com/rordenlab/dcm2niix/issues/327 float dcmFloat(int lByteLength, unsigned char lBuffer[], bool littleEndian) {//read binary 32-bit float //http://stackoverflow.com/questions/2782725/converting-float-values-from-big-endian-to-little-endian bool swap = (littleEndian != littleEndianPlatform()); @@ -1010,6 +1005,51 @@ double dcmFloatDouble(const size_t lByteLength, const unsigned char lBuffer[], //printMessage("swapped val = %f\n",retVal); return retVal; } //dcmFloatDouble() +#else + +float dcmFloat(int lByteLength, unsigned char lBuffer[], bool littleEndian) {//read binary 32-bit float + //http://stackoverflow.com/questions/2782725/converting-float-values-from-big-endian-to-little-endian + if (lByteLength < 4) return 0.0; + bool swap = (littleEndian != littleEndianPlatform()); + union { + uint32_t i; + float f; + uint8_t c[4]; + } i,o; + memcpy(&i.i, (char*)&lBuffer[0], 4); + //printf("%02x%02x%02x%02x\n",i.c[0], i.c[1], i.c[2], i.c[3]); + if (!swap) return i.f; + o.c[0] = i.c[3]; + o.c[1] = i.c[2]; + o.c[2] = i.c[1]; + o.c[3] = i.c[0]; + //printf("swp %02x%02x%02x%02x\n",o.c[0], o.c[1], o.c[2], o.c[3]); + return o.f; +} //dcmFloat() + +double dcmFloatDouble(const size_t lByteLength, const unsigned char lBuffer[], + const bool littleEndian) {//read binary 64-bit float + //http://stackoverflow.com/questions/2782725/converting-float-values-from-big-endian-to-little-endian + if (lByteLength < 8) return 0.0; + bool swap = (littleEndian != littleEndianPlatform()); + union { + uint32_t i; + double d; + uint8_t c[8]; + } i,o; + memcpy(&i.i, (char*)&lBuffer[0], 8); + if (!swap) return i.d; + o.c[0] = i.c[7]; + o.c[1] = i.c[6]; + o.c[2] = i.c[5]; + o.c[3] = i.c[4]; + o.c[4] = i.c[3]; + o.c[5] = i.c[2]; + o.c[6] = i.c[1]; + o.c[7] = i.c[0]; + return o.d; +} //dcmFloatDouble() +#endif int dcmInt (int lByteLength, unsigned char lBuffer[], bool littleEndian) { //read binary 16 or 32 bit integer if (littleEndian) { @@ -2651,8 +2691,8 @@ unsigned char * nii_loadImgCore(char* imgname, struct nifti_1_header hdr, int bi // FileSize < (ImageSize+HeaderSize): 42399788 < ( 42399792.00) //note hdr.vox_offset is a float, and without a type-cast it can lead to unusual values //https://www.nitrc.org/forum/message.php?msg_id=27155 - printMessage("FileSize < (ImageSize+HeaderSize): %lu < (%lu+%lu) \n", fileLen, imgszRead, (long)hdr.vox_offset); - //printMessage("FileSize < (ImageSize+HeaderSize): %lu < (%lu) \n", fileLen, imgszRead+(long)hdr.vox_offset); + printMessage("FileSize < (ImageSize+HeaderSize): %ld < (%zu+%ld) \n", fileLen, imgszRead, (long)hdr.vox_offset); + //printMessage("FileSize < (ImageSize+HeaderSize): %ld < (%zu) \n", fileLen, imgszRead+(long)hdr.vox_offset); printWarning("File not large enough to store image data: %s\n", imgname); return NULL; } @@ -4005,6 +4045,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kTransferSyntax 0x0002+(0x0010 << 16) #define kImplementationVersionName 0x0002+(0x0013 << 16) #define kSourceApplicationEntityTitle 0x0002+(0x0016 << 16 ) +#define kDirectoryRecordSequence 0x0004+(0x1220 << 16 ) //#define kSpecificCharacterSet 0x0008+(0x0005 << 16 ) //someday we should handle foreign characters... #define kImageTypeTag 0x0008+(0x0008 << 16 ) #define kStudyDate 0x0008+(0x0020 << 16 ) @@ -4623,9 +4664,14 @@ double TE = 0.0; //most recent echo time recorded char mediaUID[kDICOMStr]; 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; + if (strstr(mediaUID, "1.3.46.670589.11.0.0.12.1") != NULL) d.isRawDataStorage = true; //Private MR Spectrum Storage + if (strstr(mediaUID, "1.3.46.670589.11.0.0.12.2") != NULL) d.isRawDataStorage = true; //Private MR Series Data Storage + if (strstr(mediaUID, "1.3.46.670589.11.0.0.12.4") != NULL) d.isRawDataStorage = true; //Private MR Examcard Storage if (d.isRawDataStorage) d.isDerived = true; - //Philips "PS_" files + if (d.isRawDataStorage) printMessage("Skipping non-image DICOM: %s\n", fname); + //Philips "PS_" files if (strstr(mediaUID, "1.2.840.10008.5.1.4.1.1.11.1") != NULL) d.isGrayscaleSoftcopyPresentationState = true; if (d.isGrayscaleSoftcopyPresentationState) d.isDerived = true; break; @@ -4716,6 +4762,10 @@ double TE = 0.0; //most recent echo time recorded if((slen < 5) || (strstr(saeTxt, "oasis") == NULL) ) break; d.isSegamiOasis = true; break; } + case kDirectoryRecordSequence: { + d.isRawDataStorage = true; + break; + } case kImageTypeTag: { dcmStr (lLength, &buffer[lPos], d.imageType); int slen; @@ -5058,6 +5108,12 @@ double TE = 0.0; //most recent echo time recorded } patientPositionNum++; isAtFirstPatientPosition = true; + + + //char dx[kDICOMStr]; + //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 @@ -5077,7 +5133,7 @@ double TE = 0.0; //most recent echo time recorded } //if not first slice in file set_isAtFirstPatientPosition_tvd(&volDiffusion, isAtFirstPatientPosition); //if (isAtFirstPatientPosition) numFirstPatientPosition++; - if (isVerbose == 1) //verbose > 1 will report full DICOM tag + if (isVerbose > 0) //verbose > 1 will report full DICOM tag printMessage(" Patient Position 0020,0032 (#,@,X,Y,Z)\t%d\t%ld\t%g\t%g\t%g\n", patientPositionNum, lPos, patientPosition[1], patientPosition[2], patientPosition[3]); break; } case kInPlanePhaseEncodingDirection: @@ -5535,6 +5591,9 @@ double TE = 0.0; //most recent echo time recorded isIconImageSequence = true; break; case kPMSCT_RLE1 : + //https://groups.google.com/forum/#!topic/comp.protocols.dicom/8HuP_aNy9Pc + //https://discourse.slicer.org/t/fail-to-load-pet-ct-gemini/8158/3 + // d.compressionScheme = kCompressPMSCT_RLE1; //force RLE if (d.compressionScheme != kCompressPMSCT_RLE1) break; d.imageStart = (int)lPos + (int)lFileOffset; d.imageBytes = lLength; @@ -6049,6 +6108,8 @@ double TE = 0.0; //most recent echo time recorded } //while d.imageStart == 0 free (buffer); + if (d.bitsStored < 0) d.isValid = false; + if (d.bitsStored == 1) printWarning("1-bit binary DICOMs not supported\n"); //maybe not valid - no examples to test //printf("%d bval=%g bvec=%g %g %g<<<\n", d.CSA.numDti, d.CSA.dtiV[0], d.CSA.dtiV[1], d.CSA.dtiV[2], d.CSA.dtiV[3]); //printMessage("><>< DWI bxyz %g %g %g %g\n", d.CSA.dtiV[0], d.CSA.dtiV[1], d.CSA.dtiV[2], d.CSA.dtiV[3]); if (encapsulatedDataFragmentStart > 0) { @@ -6214,9 +6275,10 @@ if (d.isHasPhase) if (d.CSA.dtiV[0] > 0) printMessage(" DWI bxyz %g %g %g %g\n", d.CSA.dtiV[0], d.CSA.dtiV[1], d.CSA.dtiV[2], d.CSA.dtiV[3]); } - if ((d.xyzDim[1] > 1) && (d.xyzDim[2] > 1) && (d.imageStart < 132) && (!d.isRawDataStorage)) { + if ((d.isValid) && (d.xyzDim[1] > 1) && (d.xyzDim[2] > 1) && (d.imageStart < 132) && (!d.isRawDataStorage)) { //20190524: Philips MR 55.1 creates non-image files that report kDim1/kDim2 - we can detect them since 0008,0016 reports "RawDataStorage" - printError("Conversion aborted due to corrupt file: %s %d %d\n", fname, d.xyzDim[1], d.xyzDim[2]); + //see https://neurostars.org/t/dcm2niix-error-from-philips-dicom-qsm-data-can-this-be-skipped/4883 + printError("Conversion aborted due to corrupt file: %s %dx%d %d\n", fname, d.xyzDim[1], d.xyzDim[2], d.imageStart); #ifdef USING_R Rf_error("Irrecoverable error during conversion"); #else @@ -6379,6 +6441,12 @@ if (d.isHasPhase) d.dimensionIndexValues[MAX_NUMBER_OF_DIMENSIONS-2] = d.echoNum; if (numDimensionIndexValues < MAX_NUMBER_OF_DIMENSIONS) //https://github.com/rordenlab/dcm2niix/issues/221 d.dimensionIndexValues[MAX_NUMBER_OF_DIMENSIONS-1] = mz_crc32X((unsigned char*) &d.seriesInstanceUID, strlen(d.seriesInstanceUID)); + if ((d.isValid) && (d.seriesUidCrc == 0)) { + if (d.seriesNum < 1) + d.seriesUidCrc = 1; //no series information + else + d.seriesUidCrc = d.seriesNum; //file does not have Series UID, use series number instead + } if (d.seriesNum < 1) //https://github.com/rordenlab/dcm2niix/issues/218 d.seriesNum = mz_crc32X((unsigned char*) &d.seriesInstanceUID, strlen(d.seriesInstanceUID)); getFileName(d.imageBaseName, fname); diff --git a/console/nii_dicom.h b/console/nii_dicom.h index d6657087..b9b80b21 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -43,7 +43,7 @@ extern "C" { #define kCCsuf " CompilerNA" //unknown compiler! #endif -#define kDCMdate "v1.0.20190720" +#define kDCMdate "v1.0.20190902" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic @@ -97,7 +97,7 @@ static const int kCompressC3 = 2; //obsolete JPEG lossless 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 = 5; //LoCo JPEG-LS +static const int kCompressJPEGLS = 6; //LoCo JPEG-LS #ifdef myEnableJasper static const int kCompressSupport = kCompressYes; //JASPER for JPEG2000 #else diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 2da90327..ccc3f07a 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -375,7 +375,7 @@ void nii_saveText(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, // this is not recommended: poorly documented // it is better to stick to the binary portion of the Siemens CSA image header -#if defined(_WIN64) || defined(_WIN32) || defined(__sun) +#if defined(_WIN64) || defined(_WIN32) || defined(__sun) || (defined(__APPLE__) && defined(__POWERPC__)) //https://opensource.apple.com/source/Libc/Libc-1044.1.2/string/FreeBSD/memmem.c /*- * Copyright (c) 2005 Pascal Gloor @@ -527,8 +527,8 @@ int phoenixOffsetCSASeriesHeader(unsigned char *buff, int lLength) { TCSAitem itemCSA; for (int lT = 1; lT <= lnTag; lT++) { memcpy(&tagCSA, &buff[lPos], sizeof(tagCSA)); //read tag - //if (!littleEndianPlatform()) - // nifti_swap_4bytes(1, &tagCSA.nitems); + if (!littleEndianPlatform()) + nifti_swap_4bytes(1, &tagCSA.nitems); //printf("%d CSA of %s %d\n",lPos, tagCSA.name, tagCSA.nitems); lPos +=sizeof(tagCSA); if (strcmp(tagCSA.name, "MrPhoenixProtocol") == 0) @@ -536,8 +536,8 @@ int phoenixOffsetCSASeriesHeader(unsigned char *buff, int lLength) { for (int lI = 1; lI <= tagCSA.nitems; lI++) { memcpy(&itemCSA, &buff[lPos], sizeof(itemCSA)); lPos +=sizeof(itemCSA); - //if (!littleEndianPlatform()) - // nifti_swap_4bytes(1, &itemCSA.xx2_Len); + if (!littleEndianPlatform()) + nifti_swap_4bytes(1, &itemCSA.xx2_Len); lPos += ((itemCSA.xx2_Len +3)/4)*4; } } @@ -558,7 +558,7 @@ 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->delayTimeInTR = -0.001; csaAscii->phaseOversampling = 0.0; csaAscii->phaseResolution = 0.0; csaAscii->txRefAmp = 0.0; @@ -618,7 +618,6 @@ void siemensCsaAscii(const char * filename, TCsaAscii* csaAscii, int csaOffset, if ((keyPosEnd) && ((keyPosEnd - keyPos) < csaLengthTrim)) //ignore binary data at end csaLengthTrim = (int)(keyPosEnd - keyPos); #endif - char keyStrLns[] = "sKSpace.lPhaseEncodingLines"; csaAscii->phaseEncodingLines = readKey(keyStrLns, keyPos, csaLengthTrim); char keyStrUcImg[] = "sSliceArray.ucImageNumb"; @@ -699,7 +698,7 @@ void siemensCsaAscii(const char * filename, TCsaAscii* csaAscii, int csaOffset, //check if ANY csaAscii.adFree tags exist keyPosFree = (char *)memmem(keyPos, csaLengthTrim, keyStrAdFree, strlen(keyStrAdFree)); if (!keyPosFree) { //"Wip" -> "WiP", modern -> old Siemens - strcpy(keyStrAdFree, "sWipMemBlock.adFree["); + strcpy(keyStrAdFree, "sWiPMemBlock.adFree["); keyPosFree = (char *)memmem(keyPos, csaLengthTrim, keyStrAdFree, strlen(keyStrAdFree)); } if (keyPosFree) { @@ -1167,7 +1166,7 @@ tse3d: T2*/ } //verbose */ //ASL specific tags - 2D pCASL Danny J.J. Wang http://www.loft-lab.org - if (strstr(pulseSequenceDetails,"ep2d_pcasl")) { + if ((strstr(pulseSequenceDetails,"ep2d_pcasl")) || (strstr(pulseSequenceDetails,"ep2d_pcasl_UI_PHC"))) { json_FloatNotNan(fp, "\t\"LabelOffset\": %g,\n", csaAscii.adFree[1]); //mm json_FloatNotNan(fp, "\t\"PostLabelDelay\": %g,\n", csaAscii.adFree[2] * (1.0/1000000.0)); //usec -> sec json_FloatNotNan(fp, "\t\"NumRFBlocks\": %g,\n", csaAscii.adFree[3]); @@ -1490,9 +1489,9 @@ tse3d: T2*/ if (d.phaseEncodingRC == 'R') fprintf(fp, "\t\"InPlanePhaseEncodingDirectionDICOM\": \"ROW\",\n" ); // Finish up with info on the conversion tool fprintf(fp, "\t\"ConversionSoftware\": \"dcm2niix\",\n"); - fprintf(fp, "\t\"ConversionSoftwareVersion\": \"%s\"\n", kDCMvers ); - //fprintf(fp, "\t\"DicomConversion\": [\"dcm2niix\", \"%s\"]\n", kDCMvers ); - fprintf(fp, "}\n"); + fprintf(fp, "\t\"ConversionSoftwareVersion\": \"%s\"\n", kDCMdate ); + //fprintf(fp, "\t\"ConversionSoftwareVersion\": \"%s\"\n", kDCMvers );kDCMdate + fprintf(fp, "}\n"); fclose(fp); }// nii_SaveBIDS() @@ -2337,7 +2336,7 @@ int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopt #ifdef myMultiEchoFilenameSkipEcho1 if ((!isEchoReported) && (dcm.isMultiEcho) && (dcm.echoNum >= 1)) { //multiple echoes saved as same series #else - if ((!isEchoReported) && (dcm.isMultiEcho)) { //multiple echoes saved as same series + if ((!isEchoReported) && ((dcm.isMultiEcho) || (dcm.echoNum > 1))) { //multiple echoes saved as same series #endif sprintf(newstr, "_e%d", dcm.echoNum); strcat (outname,newstr); @@ -2759,14 +2758,14 @@ void nii_saveAttributes (struct TDICOMdata &data, struct nifti_1_header &header, images->addAttribute("phaseEncodingSign", data.CSA.phaseEncodingDirectionPositive == 0 ? -1 : 1); } } - // Slice timing + + // Slice timing (stored in seconds) if (data.CSA.sliceTiming[0] >= 0.0 && (data.manufacturer == kMANUFACTURER_UIH || data.manufacturer == kMANUFACTURER_GE || (data.manufacturer == kMANUFACTURER_SIEMENS && !data.isXA10A))) { std::vector sliceTimes; - #pragma message ("Please test R specific code: at this stage all slice times should be in msec due to changes in checkSliceTiming() 20190704") - for (int i=header.dim[3]-1; i>=0; i--) { - if (data.CSA.sliceTiming[i] < 0.0) - break; - sliceTimes.push_back(data.CSA.sliceTiming[i]); //slice time in msec + for (int i=0; iaddAttribute("sliceTiming", sliceTimes); } @@ -2791,12 +2790,6 @@ void nii_saveAttributes (struct TDICOMdata &data, struct nifti_1_header &header, #else -inline bool littleEndianPlatformNRRD () -{ - uint32_t value = 1; - return (*((char *) &value) == 1); -} - int pigz_File(char * fname, struct TDCMopts opts, size_t imgsz) { //given "/dir/file.nii" creates "/dir/file.nii.gz" char blockSize[768]; @@ -2859,7 +2852,7 @@ int nii_saveNRRD(char * niiFilename, struct nifti_1_header hdr, unsigned char* i fprintf(fp,"NRRD0005\n"); fprintf(fp,"# Complete NRRD file format specification at:\n"); fprintf(fp,"# http://teem.sourceforge.net/nrrd/format.html\n"); - fprintf(fp,"# dcm2niix NRRD export transforms by Tashrif Billah\n"); + fprintf(fp,"# dcm2niix %s NRRD export transforms by Tashrif Billah\n", kDCMdate); char rgbNoneStr[10] = {""}; //type tag switch (hdr.datatype) { @@ -2914,7 +2907,7 @@ int nii_saveNRRD(char * niiFilename, struct nifti_1_header hdr, unsigned char* i fprintf(fp,"\n"); } //byteskip only for .nhdr, not .nrrd - if (littleEndianPlatformNRRD()) //raw data in native format + if (littleEndianPlatform()) //raw data in native format fprintf(fp,"endian: little\n"); else fprintf(fp,"endian: big\n"); @@ -3070,6 +3063,24 @@ int nii_saveNRRD(char * niiFilename, struct nifti_1_header hdr, unsigned char* i return pigz_File(fname, opts, imgsz); } // nii_saveNRRD() +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 + // one could also auto-detect: hdr->sizeof_hdr==348 + if (!isNative) swap_nifti_header(hdr); + int nVox = 1; + for (int i = 1; i < 8; i++) + if (hdr->dim[i] > 1) nVox = nVox * hdr->dim[i]; + int bitpix = hdr->bitpix; + int datatype = hdr->datatype; + if (isNative) swap_nifti_header(hdr); + if (datatype == DT_RGBA32) return; + //n.b. do not swap 8-bit, 24-bit RGB, and 32-bit RGBA + if (bitpix == 16) nifti_swap_2bytes(nVox, im); + if (bitpix == 32) nifti_swap_4bytes(nVox, im); + if (bitpix == 64) nifti_swap_8bytes(nVox, im); +} + int nii_saveNII(char * niiFilename, struct nifti_1_header hdr, unsigned char* im, struct TDCMopts opts, struct TDICOMdata d) { if (opts.isOnlyBIDS) return EXIT_SUCCESS; if (opts.isSaveNRRD) { @@ -3101,7 +3112,9 @@ int nii_saveNII(char * niiFilename, struct nifti_1_header hdr, unsigned char* im if ((imgsz+hdr.vox_offset) < kMaxPigz) printWarning(" Hint: using external compressor (pigz) should help.\n"); } else if ((opts.isGz) && (strlen(opts.pigzname) < 1) && ((imgsz+hdr.vox_offset) < kMaxGz) ) { //use internal compressor + if (!opts.isSaveNativeEndian) swapEndian(&hdr, im, true); //byte-swap endian (e.g. little->big) writeNiiGz (niiFilename, hdr, im, imgsz, opts.gzLevel, false); + if (!opts.isSaveNativeEndian) swapEndian(&hdr, im, false); //unbyte-swap endian (e.g. big->little) return EXIT_SUCCESS; } #endif @@ -3135,21 +3148,25 @@ int nii_saveNII(char * niiFilename, struct nifti_1_header hdr, unsigned char* im printError("Unable to open pigz pipe\n"); return EXIT_FAILURE; } + if (!opts.isSaveNativeEndian) swapEndian(&hdr, im, true); //byte-swap endian (e.g. little->big) fwrite(&hdr, sizeof(hdr), 1, pigzPipe); uint32_t pad = 0; fwrite(&pad, sizeof( pad), 1, pigzPipe); fwrite(&im[0], imgsz, 1, pigzPipe); pclose(pigzPipe); + if (!opts.isSaveNativeEndian) swapEndian(&hdr, im, false); //unbyte-swap endian (e.g. big->little) return EXIT_SUCCESS; } #endif FILE *fp = fopen(fname, "wb"); if (!fp) return EXIT_FAILURE; + if (!opts.isSaveNativeEndian) swapEndian(&hdr, im, true); //byte-swap endian (e.g. little->big) fwrite(&hdr, sizeof(hdr), 1, fp); uint32_t pad = 0; fwrite(&pad, sizeof( pad), 1, fp); fwrite(&im[0], imgsz, 1, fp); fclose(fp); + if (!opts.isSaveNativeEndian) swapEndian(&hdr, im, false); //unbyte-swap endian (e.g. big->little) if ((opts.isGz) && (strlen(opts.pigzname) > 0) ) { #ifndef myDisableGzSizeLimits if ((imgsz+hdr.vox_offset) > kMaxPigz) { @@ -4101,7 +4118,7 @@ void checkDateTimeOrder(struct TDICOMdata * d, struct TDICOMdata * d1) { printWarning("Images sorted by instance number [0020,0013](%d..%d), but AcquisitionTime [0008,0032] suggests a different order (%g..%g) \n", d->imageNum,d1->imageNum, d->acquisitionTime,d1->acquisitionTime); } -void checkSliceTiming(struct TDICOMdata * d, struct TDICOMdata * d1, int verbose) { +void checkSliceTiming(struct TDICOMdata * d, struct TDICOMdata * d1, int verbose, int isForceSliceTimeHHMMSS) { //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 @@ -4110,6 +4127,7 @@ void checkSliceTiming(struct TDICOMdata * d, struct TDICOMdata * d1, int verbose nSlices++; if (nSlices < 1) return; bool isSliceTimeHHMMSS = (d->manufacturer == kMANUFACTURER_UIH); + if (isForceSliceTimeHHMMSS) isSliceTimeHHMMSS = true; //if (d->isXA10A) isSliceTimeHHMMSS = true; //for XA10 use TimeAfterStart 0x0021,0x1104 -> Siemens de-identification can corrupt acquisition ties https://github.com/rordenlab/dcm2niix/issues/236 if (isSliceTimeHHMMSS) {//handle midnight crossing for (int i = 0; i < nSlices; i++) @@ -4125,7 +4143,7 @@ void checkSliceTiming(struct TDICOMdata * d, struct TDICOMdata * d1, int verbose float kNoonSec = 43200; if ((maxT - minT) > kNoonSec) { //volume started before midnight but ended next day! //identify and fix 'Cinderella error' where clock resets at midnight: untested - printWarning("UIH acquisition crossed midnight: check slice timing\n"); + printWarning("Acquisition crossed midnight: check slice timing\n"); for (int i = 0; i < nSlices; i++) if (d->CSA.sliceTiming[i] > kNoonSec) d->CSA.sliceTiming[i] = d->CSA.sliceTiming[i] - kMidnightSec; minT = d->CSA.sliceTiming[0]; @@ -4271,11 +4289,33 @@ void reverseSliceTiming(struct TDICOMdata * d, int verbose, int nSL) { d->CSA.sliceTiming[i] = sliceTiming[(nSL-1)-i]; } +int sliceTimingSiemens2D(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct nifti_1_header * hdr, int verbose, const char * filename, int nConvert) { + //only for Siemens 2D images, use acquisitionTime + uint64_t indx0 = dcmSort[0].indx; //first volume + if (!(dcmList[indx0].manufacturer == kMANUFACTURER_SIEMENS)) return 0; + if (dcmList[indx0].is3DAcq) return 0; //no need for slice times + if (dcmList[indx0].CSA.sliceTiming[0] >= 0.0) return 0; //slice times calculated + if (dcmList[indx0].CSA.mosaicSlices > 1) return 0; + if (nConvert != (hdr->dim[3]*hdr->dim[4])) return 0; + if (hdr->dim[3] > (kMaxEPI3D-1)) return 0; + int nZero = 0; //infer multiband: E11C may not populate kPATModeText + for (int v = 0; v < hdr->dim[3]; v++) { + 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() + if (dcmList[indx0].CSA.sliceTiming[v] == dcmList[indx0].CSA.sliceTiming[0]) nZero++; + } + if ((dcmList[indx0].CSA.multiBandFactor < 2) && (nZero > 1)) + dcmList[indx0].CSA.multiBandFactor = nZero; + return 1; +} + void rescueSliceTimingSiemens(struct TDICOMdata * d, int verbose, int nSL, const char * filename) { if (d->is3DAcq) return; //no need for slice times + if (d->CSA.multiBandFactor > 1) return; //pattern of multiband slice order unknown if (d->CSA.sliceTiming[0] >= 0.0) return; //slice times calculated + if (d->CSA.mosaicSlices < 2) return; //20190807 E11C 2D (not mosaic) files do not report mosaicAcqTimes or multi-band factor. if (nSL < 2) return; if ((d->manufacturer != kMANUFACTURER_SIEMENS) || (d->CSA.SeriesHeader_offset < 1) || (d->CSA.SeriesHeader_length < 1)) return; +#ifdef myReadAsciiCsa float shimSetting[8]; char protocolName[kDICOMStrLarge], fmriExternalInfo[kDICOMStrLarge], coilID[kDICOMStrLarge], consistencyInfo[kDICOMStrLarge], coilElements[kDICOMStrLarge], pulseSequenceDetails[kDICOMStrLarge], wipMemBlock[kDICOMStrLarge]; TCsaAscii csaAscii; @@ -4315,6 +4355,7 @@ void rescueSliceTimingSiemens(struct TDICOMdata * d, int verbose, int nSL, const } //if ucMode == 3 int //dicm2nii provides sSliceArray.ucImageNumb - similar to protocolSliceNumber1 //if asc_header(s, 'sSliceArray.ucImageNumb'), t = t(nSL:-1:1); end % rev-num +#endif } void sliceTimingUIH(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct nifti_1_header * hdr, int verbose, const char * filename, int nConvert) { @@ -4330,8 +4371,7 @@ void sliceTimingUIH(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct void sliceTimingGE(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 - //printf("XXX %d\n", indx0); - if (!(dcmList[indx0].manufacturer == kMANUFACTURER_GE)) return; + if (!(dcmList[indx0].manufacturer == kMANUFACTURER_GE)) return; bool GEsliceTiming_x0018x1060 = false; if ((hdr->dim[3] < (kMaxEPI3D-1)) && (hdr->dim[3] > 1) && (hdr->dim[4] > 1)) { //GE: 1st method for "epi" PSD @@ -4430,8 +4470,9 @@ int sliceTimingCore(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct struct TDICOMdata * d1 = &dcmList[indx1]; sliceTimingGE(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); + 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); @@ -4793,26 +4834,35 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc for(int i = 0; i < nConvert; ++i) dcmList[dcmSort[i].indx].converted2NII = 1; if (opts.numSeries < 0) { //report series number but do not convert - if (segVol >= 0) { - printMessage("\t%ld.%d\t%s\n", dcmList[dcmSort[0].indx].seriesNum, segVol-1, pathoutname); + int segVolEcho = segVol; + if ((dcmList[dcmSort[0].indx].echoNum > 1) && (segVolEcho <= 0)) + segVolEcho = dcmList[dcmSort[0].indx].echoNum+1; + if (segVolEcho >= 0) { + printMessage("\t%u.%d\t%s\n", dcmList[dcmSort[0].indx].seriesUidCrc, segVolEcho-1, pathoutname); + //printMessage("\t%ld.%d\t%s\n", dcmList[dcmSort[0].indx].seriesNum, segVol-1, pathoutname); } else { - printMessage("\t%ld\t%s\n", dcmList[dcmSort[0].indx].seriesNum, pathoutname); + printMessage("\t%u\t%s\n", dcmList[dcmSort[0].indx].seriesUidCrc, pathoutname); + //printMessage("\t%ld\t%s\n", dcmList[dcmSort[0].indx].seriesNum, pathoutname); } printMessage(" %s\n",nameList->str[dcmSort[0].indx]); return EXIT_SUCCESS; } - if (sliceDir < 0) { + if (sliceDir < 0) { 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! } // skip converting if user has specified one or more series, but has not specified this one if (opts.numSeries > 0) { int i = 0; - float seriesNum = (float) dcmList[dcmSort[0].indx].seriesNum; - if (segVol > 0) - seriesNum = seriesNum + ((float) segVol - 1.0) / 10.0; //n.b. we will have problems if segVol > 9. However, 9 distinct TEs/scalings/PhaseMag seems unlikely + //double seriesNum = (double) dcmList[dcmSort[0].indx].seriesNum; + double seriesNum = (double) dcmList[dcmSort[0].indx].seriesUidCrc; + int segVolEcho = segVol; + if ((dcmList[dcmSort[0].indx].echoNum > 1) && (segVolEcho <= 0)) + segVolEcho = dcmList[dcmSort[0].indx].echoNum+1; + if (segVolEcho > 0) + seriesNum = seriesNum + ((double) segVolEcho - 1.0) / 10.0; for (; i < opts.numSeries; i++) { - if (isSameFloatGE(opts.seriesNumber[i], seriesNum)) { + if (isSameDouble(opts.seriesNumber[i], seriesNum)) { //if (opts.seriesNumber[i] == dcmList[dcmSort[0].indx].seriesNum) { break; } @@ -4821,7 +4871,7 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc return EXIT_SUCCESS; } } - nii_saveText(pathoutname, dcmList[dcmSort[0].indx], opts, &hdr0, nameList->str[indx]); + nii_saveText(pathoutname, dcmList[dcmSort[0].indx], opts, &hdr0, nameList->str[indx]); int numADC = 0; int * volOrderIndex = nii_saveDTI(pathoutname,nConvert, dcmSort, dcmList, opts, sliceDir, dti4D, &numADC); PhilipsPrecise(&dcmList[dcmSort[0].indx], opts.isPhilipsFloatNotDisplayScaling, &hdr0, opts.isVerbose); @@ -4837,7 +4887,7 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc #ifndef USING_R fflush(stdout); //show immediately if run from MRIcroGL GUI #endif - //~ if (!dcmList[dcmSort[0].indx].isSlicesSpatiallySequentialPhilips) + //~ if (!dcmList[dcmSort[0].indx].isSlicesSpatiallySequentialPhilips) //~ nii_reorderSlices(imgM, &hdr0, dti4D); //hdr0.pixdim[3] = dxNoTilt; if (hdr0.dim[3] < 2) @@ -4879,10 +4929,9 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc } //avoid div/0: cosine not zero } //if gantry tilt //end: gantry tilt we need to save the shear in the transform - + int returnCode = EXIT_FAILURE; #ifndef myNoSave // Indicates success or failure of the (last) save - int returnCode = EXIT_FAILURE; //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 @@ -4939,7 +4988,7 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc #ifdef USING_R // Note that for R, only one image should be created per series // Hence this extra test - if (returnCode != EXIT_SUCCESS) + if (returnCode != EXIT_SUCCESS) returnCode = nii_saveNII(pathoutname, hdr0, imgM, opts, dcmList[dcmSort[0].indx]); if (returnCode == EXIT_SUCCESS) nii_saveAttributes(dcmList[dcmSort[0].indx], hdr0, opts, nameList->str[dcmSort[0].indx]); @@ -6276,6 +6325,7 @@ void setDefaultOpts (struct TDCMopts *opts, const char * argv[]) { //either "set opts->isCrop = false; opts->isRotate3DAcq = true; opts->isGz = false; + opts->isSaveNativeEndian = true; 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 c9880e02..9683e5dc 100644 --- a/console/nii_dicom_batch.h +++ b/console/nii_dicom_batch.h @@ -30,10 +30,10 @@ extern "C" { #define MAX_NUM_SERIES 16 struct TDCMopts { - bool isSaveNRRD, isOneDirAtATime, isRenameNotConvert, isMaximize16BitRange, isSave3D, isGz, isPipedGz, isFlipY, isCreateBIDS, isSortDTIbyBVal, isAnonymizeBIDS, isOnlyBIDS, isCreateText, isIgnoreDerivedAnd2D, isPhilipsFloatNotDisplayScaling, isTiltCorrect, isRGBplanar, isOnlySingleFile, isForceStackDCE, isRotate3DAcq, isCrop; + bool isSaveNativeEndian, isSaveNRRD, isOneDirAtATime, isRenameNotConvert, isMaximize16BitRange, isSave3D, isGz, isPipedGz, isFlipY, isCreateBIDS, isSortDTIbyBVal, isAnonymizeBIDS, isOnlyBIDS, isCreateText, isIgnoreDerivedAnd2D, isPhilipsFloatNotDisplayScaling, isTiltCorrect, isRGBplanar, isOnlySingleFile, isForceStackDCE, isRotate3DAcq, isCrop; int 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]; - float seriesNumber[MAX_NUM_SERIES]; + double seriesNumber[MAX_NUM_SERIES]; //requires double must store -1 (report but do not convert) as well as seriesUidCrc (uint32) long numSeries; #ifdef USING_R bool isScanOnly;