diff --git a/Philips/README.md b/Philips/README.md new file mode 100644 index 00000000..b961f8d0 --- /dev/null +++ b/Philips/README.md @@ -0,0 +1,46 @@ +## About + +dcm2niix attempts to convert all DICOM images to NIfTI. The Philips enhanced DICOM images are elegantly able to save all images from a series as a single file. However, this format is necessarily complex. The usage of this format has evolved over time, and can become further complicated when DICOM are handled by DICOM tools (for example, anonymization, transfer which converts explicit VRs to implicit VRs, etc.). + +This web page describes some of the strategies handle these images. However, users should be vigilant when handling these datasets. If you encounter problems using dcm2niix you can explore [https://github.com/rordenlab/dcm2niix competing DICOM to NIfTI converters] or [ [https://github.com/rordenlab/dcm2niix report an issue]. + +## Image Patient Position + +The Image Patient Position (0020,0032) tag is required to determine the position of the slices with respect to each other (and with respect to the scanner bore). Philips scans often report two conflicting IPPs for each single slice: with one stored in the private sequence (SQ) 2005,140F while the other is in the public sequence. This is unusual, but is [ftp://medical.nema.org/medical/dicom/final/cp758_ft.pdf legal]. + +In practice, this complication has several major implications. First, not all software translate private SQs well. One potential problem is when the implicit VRs are saved as implicit VRs. This can obscure the fact that 2005,140F is an SQ. Indeed, some tools will convert the private SQ type as a "UN" unknown type and add another public sequence. This can confuse reading software. + +Furthermore, in the real world there are many Philips DICOM images that ONLY contain IPP inside SQ 2005,140F. These situations appear to reflect modifications applied by a PACS to the DICOM data or attempts to anonymize the DICOM images (e.g. using simple Matlab scripts). Note that the private IPP differs from the public one by half a voxel. Therefore, in theory if one only has the private IPP, one can infer the public IPP location. Current versions of dcm2niix do not do this: the error is assumed small enough that it will not impact image processing steps such as coregistration. + +In general it is recommended that you archive and convert DICOM images as they are created from the scanner. If one does use an export tool such as the excellent dcmtk, it is recommended that you preserve the explicit VR, as implicit VR has the potential of obscuring private sequence (SQ) tags. Be aware that subsequent processing of DICOM data can disrupt data conversion. + +Therefore, dcm2niix will ignore the IPP enclosed in 2005,140F unless no alternative exists. + +## 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. For example, ADC, Trace and Isotropic images can all be derived from the raw scans. 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"). + +``` +(0018,9075) CS [DIRECTIONAL] +(0018,9089) FD 1\0\0 +(0018,9087) FD 1000 +(0020,9157) UL 1\1\2\32 +... +(0018,9075) CS [ISOTROPIC] +(0018,9087) FD 1000 +(0020,9157) UL 1\2\1\33 +... +(0018,9075) CS [NONE] +(0018,9087) FD 0 +(0020,9157) UL 1\1\2\33 +``` + +## Diffusion Direction + +Proper Philips enhanced scans use tag 0018,9089 to report the 3 gradient directions. However, in the wild, other files from Philips (presumably using older versions of Philips software) use the tags 2005,10b0, 2005,10b1, 2005,10b2. In general, dcm2niix will use the values that most closely precede the Dimension Index Values (0020,9157). + +The current version of dcm2niix uses Dimension Index Values (0020,9157) to determine gradient number, which can also be found in (2005,1413). However, while 2005,1413 is always indexed from one, this is not necessarily the case for 0020,9157. For example, the ADNI DWI dataset for participant 018_S_4868 has values of 2005,1413 that range from 1..36 for the 36 directions, while 0020,9157 ranges from 2..37. The current version of dcm2niix compensates for this by re-indexing the values of 0020,9157 after all the volumes have been read. + +## General variations + +Prior versions of dcm2niix used different methods to sort images. However, these have proved unreliable The undocumented tags SliceNumberMrPhilips (2001,100A), NumberOfSlicesMrPhilips (2001,1018) are not always present. 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 depend on any of these. diff --git a/README.md b/README.md index 29b5ecc1..25c7c2b9 100644 --- a/README.md +++ b/README.md @@ -59,20 +59,32 @@ The batch processing binary `dcm2niibatch` is optional. To build `dcm2niibatch` If you have any problems with the cmake build script described above or want to customize the software see the [COMPILE.md file for details on manual compilation](./COMPILE.md). +## Alternatives + + - [Valerio Luccio's dinifti](http://cbi.nyu.edu/software/dinifti.php) is focused on conversion of Siemens data + - [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). + - [Xiangrui Li's 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. + - [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. + - [Jolinda Smith's 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. + - [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. + ## Links - - [Dcm2Bids](https://github.com/cbedetti/Dcm2Bids) uses dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets. + - [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. - [bidskit](https://github.com/jmtyszka/bidskit) uses dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets. + - [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. - - [heudiconv](https://github.com/nipy/heudiconv) can use dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets. - - [nipype](https://github.com/nipy/nipype) can use dcm2niix to convert images. - - [pydcm2niix is a Python module for working with dcm2niix](https://github.com/jstutters/pydcm2niix). + - [Dcm2Bids](https://github.com/cbedetti/Dcm2Bids) uses dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets. - [dcm2niir](https://github.com/muschellij2/dcm2niir) R wrapper for dcm2niix/dcm2nii. - - [divest](https://github.com/jonclayden/divest) R interface to dcm2niix. - - [sci-tran dcm2niix](https://github.com/scitran-apps/dcm2niix) Flywheel Gear (docker). - - [neuro_docker](https://github.com/Neurita/neuro_docker) includes dcm2niix as part of a single, static Dockerfile. - - [boutiques-dcm2niix](https://github.com/lalet/boutiques-dcm2niix) is a dockerfile for installing and validating dcm2niix. - - [neurodocker](https://github.com/kaczmarj/neurodocker) generates [custom](https://github.com/rordenlab/dcm2niix/issues/138) Dockerfiles given specific versions of neuroimaging software. - [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. + - [divest](https://github.com/jonclayden/divest) R interface to dcm2niix. + - [heudiconv](https://github.com/nipy/heudiconv) can use dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets. - [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/). + - [neuro_docker](https://github.com/Neurita/neuro_docker) includes dcm2niix as part of a single, static Dockerfile. - [NeuroDebian](http://neuro.debian.net/pkgs/dcm2niix.html) provides up-to-date version of dcm2niix for Debian-based systems. + - [neurodocker](https://github.com/kaczmarj/neurodocker) generates [custom](https://github.com/rordenlab/dcm2niix/issues/138) Dockerfiles given specific versions of neuroimaging software. + - [nipype](https://github.com/nipy/nipype) can use dcm2niix to convert images. + - [pydcm2niix is a Python module for working with dcm2niix](https://github.com/jstutters/pydcm2niix). + - [sci-tran dcm2niix](https://github.com/scitran-apps/dcm2niix) Flywheel Gear (docker). \ No newline at end of file diff --git a/console/main_console.cpp b/console/main_console.cpp index 78cf75ba..074bdea0 100644 --- a/console/main_console.cpp +++ b/console/main_console.cpp @@ -85,7 +85,7 @@ void showHelp(const char * argv[], struct TDCMopts opts) { #else #define kQstr "" #endif - printf(" -f : filename (%%a=antenna (coil) number, %%c=comments, %%d=description, %%e echo number, %%f=folder name, %%i ID of patient, %%j seriesInstanceUID, %%k studyInstanceUID, %%m=manufacturer, %%n=name of patient, %%p=protocol,%s %%s=series number, %%t=time, %%u=acquisition number, %%v=vendor, %%x=study ID; %%z sequence name; default '%s')\n", kQstr, opts.filename); + printf(" -f : filename (%%a=antenna (coil) number, %%c=comments, %%d=description, %%e=echo number, %%f=folder name, %%i=ID of patient, %%j=seriesInstanceUID, %%k=studyInstanceUID, %%m=manufacturer, %%n=name of patient, %%p=protocol,%s %%s=series number, %%t=time, %%u=acquisition number, %%v=vendor, %%x=study ID; %%z=sequence name; default '%s')\n", kQstr, opts.filename); printf(" -g : generate defaults file (y/n/o [o=only: reset and write defaults], default n)\n"); printf(" -h : show help\n"); printf(" -i : ignore derived, localizer and 2D images (y/n, default n)\n"); @@ -324,7 +324,11 @@ int main(int argc, const char * argv[]) if (invalidParam(i, argv)) return 0; if ((argv[i][0] == 'n') || (argv[i][0] == 'N') || (argv[i][0] == '0')) //0: verbose OFF opts.isVerbose = 0; - else if ((argv[i][0] == 'h') || (argv[i][0] == 'H') || (argv[i][0] == '2')) //2: verbose HYPER + else if ((argv[i][0] == 'o') || (argv[i][0] == 'O')) {//-1: experimental SQ skipping: 'O'ptimized!? faster, does not get confused with Philips tags + opts.isVerbose = -1; + printf(">> Experimental SQ skipping enabled\n"); + + } else if ((argv[i][0] == 'h') || (argv[i][0] == 'H') || (argv[i][0] == '2')) //2: verbose HYPER opts.isVerbose = 2; else opts.isVerbose = 1; //1: verbose ON @@ -356,7 +360,7 @@ int main(int argc, const char * argv[]) strcpy(opts.outdir,argv[i]); } else if ((argv[i][1] == 'n') && ((i+1) < argc)) { i++; - int seriesNumber = atoi(argv[i]); + float seriesNumber = atof(argv[i]); if (seriesNumber < 0) opts.numSeries = -1; //report series: convert none else if ((opts.numSeries >= 0) && (opts.numSeries < MAX_NUM_SERIES)) { diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 933782ad..93c5a8c1 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -692,9 +692,11 @@ struct TDICOMdata clear_dicom_data() { strcpy(d.studyInstanceUID, ""); strcpy(d.bodyPartExamined,""); d.phaseEncodingLines = 0; - d.patientPositionSequentialRepeats = 0; + //~ d.patientPositionSequentialRepeats = 0; + //~ d.patientPositionRepeats = 0; d.isHasPhase = false; d.isHasMagnitude = false; + //d.maxGradDynVol = -1; //PAR/REC only d.sliceOrient = kSliceOrientUnknown; d.dateTime = (double)19770703150928.0; d.acquisitionTime = 0.0f; @@ -712,15 +714,17 @@ struct TDICOMdata clear_dicom_data() { d.pixelBandwidth = 0.0; d.zSpacing = 0.0; d.zThick = 0.0; - d.numberOfDynamicScans = 0; + //~ d.numberOfDynamicScans = 0; d.echoNum = 1; d.echoTrainLength = 0; d.phaseFieldofView = 0.0; d.dwellTime = 0; + d.protocolBlockStartGE = 0; + d.protocolBlockLengthGE = 0; d.phaseEncodingSteps = 0; d.coilNum = 0; d.accelFactPE = 0.0; - d.patientPositionNumPhilips = 0; + //d.patientPositionNumPhilips = 0; d.imageBytes = 0; d.intenScale = 1; d.intenScalePhilips = 0; @@ -738,9 +742,10 @@ struct TDICOMdata clear_dicom_data() { d.imageStart = 0; d.is3DAcq = false; //e.g. MP-RAGE, SPACE, TFE d.is2DAcq = false; // - d.isSlicesSpatiallySequentialPhilips = true; //Philips can save slices in random order, e.g. 4,5,6,1,2,3 + //~ d.isSlicesSpatiallySequentialPhilips = true; //Philips can save slices in random order, e.g. 4,5,6,1,2,3 d.isDerived = false; //0008,0008 = DERIVED,CSAPARALLEL,POSDISP d.isSegamiOasis = false; //these images do not store spatial coordinates + d.isScaleOrTEVaries = false; d.bitsAllocated = 16;//bits d.bitsStored = 0; d.samplesPerPixel = 1; @@ -750,6 +755,7 @@ struct TDICOMdata clear_dicom_data() { d.isSigned = false; //default is unsigned! d.isFloat = false; //default is for integers, not single or double precision d.isResampled = false; //assume data not resliced to remove gantry tilt problems + d.isLocalizer = false; d.compressionScheme = 0; //none d.isExplicitVR = true; d.isLittleEndian = true; //DICOM initially always little endian @@ -1190,6 +1196,13 @@ int readCSAImageHeader(unsigned char *buff, int lLength, struct TCSAdata *CSA, i return EXIT_SUCCESS; } // readCSAImageHeader() +//int readCSAImageHeader(unsigned char *buff, int lLength, struct TCSAdata *CSA, int isVerbose) { +int readProtocolDataBlockGE(unsigned char *buff, int lLength, struct TCSAdata *CSA, int isVerbose) { + + return EXIT_SUCCESS; +} + + void dcmMultiShorts (int lByteLength, unsigned char lBuffer[], int lnShorts, uint16_t *lShorts, bool littleEndian) { //read array of unsigned shorts US http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_6.2.html if ((lnShorts < 1) || (lByteLength != (lnShorts * 2))) return; @@ -1379,7 +1392,7 @@ void changeExt (char *file_name, const char* ext) { } } //changeExt() -struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D *dti4D) { +struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D *dti4D, bool isReadPhase) { struct TDICOMdata d = clear_dicom_data(); strcpy(d.protocolName, ""); //erase dummy with empty strcpy(d.seriesDescription, ""); //erase dummy with empty @@ -1421,7 +1434,15 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D #define kv2 45 #define kv3 46 #define kASL 48 + printWarning("dcm2niix PAR is not actively supported. Use with extreme caution (hint: use dicm2nii)\n"); + if (isReadPhase) printWarning(" Reading phase images from PAR/REC\n"); char buff[LINESZ]; + int maxNumberOfDiffusionValues = 1; + int maxNumberOfGradientOrients = 1; + int maxNumberOfCardiacPhases = 1; + int maxNumberOfEchoes = 1; + int maxNumberOfDynamics = 1; + int maxNumberOfMixes = 1; int sliceNumberMrPhilipsB2[kMaxDTI4D], sliceNumberMrPhilipsVol2[kMaxDTI4D]; int patientPositionNumPhilipsB2 = 0; int patientPositionNumPhilipsVol2 = 0; @@ -1432,6 +1453,7 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D int minDyn = 32767; int maxDyn = 0; bool ADCwarning = false; + int numSlice2D = 0; int prevDyn = -1; bool dynNotAscending = false; int parVers = 0; @@ -1439,24 +1461,42 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D int maxCardiac = 1; int nCols = 26; int slice = 0; + int numExpected = 0; + int maxVol = -1; + + int patientPositionNumPhilips = 0; + bool isSlicesSpatiallySequentialPhilips = true; //int prevSliceIndex = 0; //index of prior slice: detect if images are not in order const int kMaxCols = 49; float *cols = (float *)malloc(sizeof(float) * kMaxCols); char *p = fgets (buff, LINESZ, fp); bool isIntenScaleVaries = false; bool isIndexSequential = true; + int prevSlice = 0; for (int i = 0; i < kMaxDTI4D; i++) dti4D->S[i].V[0] = -1.0; + dti4D->sliceOrder[0] = -1; + //dti4D->S[0].sliceNumberMrPhilips = -1.0; + //dti4D->S[0].sliceNumberMrPhilipsVol2 = -1.0; //d.dti4D = (TDTI *)malloc(kMaxDTI4D * sizeof(TDTI)); + while (p) { if (strlen(buff) < 1) continue; if (buff[0] == '#') { //comment char Comment[7][50]; - sscanf(buff, "# %s %s\n", Comment[0], Comment[1]); + + //sscanf(buff, "# %s %s\n", Comment[0], Comment[1]); + + sscanf(buff, "# %s %s %s %s %s %s V%s\n", Comment[0], Comment[1], Comment[2], Comment[3],Comment[4], Comment[5],Comment[6]); + + if ((strcmp(Comment[0], "sl") == 0) && (strcmp(Comment[1], "ec") == 0) ) { + //printError("BINGO\n"); + numExpected = (d.xyzDim[3] * maxNumberOfDiffusionValues * maxNumberOfGradientOrients + * maxNumberOfCardiacPhases * maxNumberOfEchoes * maxNumberOfDynamics * maxNumberOfMixes); + } if (strcmp(Comment[1], "TRYOUT") == 0) { - sscanf(buff, "# %s %s %s %s %s %s V%s\n", Comment[0], Comment[1], Comment[2], Comment[3] - ,Comment[4], Comment[5],Comment[6]); + //sscanf(buff, "# %s %s %s %s %s %s V%s\n", Comment[0], Comment[1], Comment[2], Comment[3],Comment[4], Comment[5],Comment[6]); parVers = (int)round(atof(Comment[6])*10); //4.2 = 42 etc if (parVers < 40) { printMessage("This software is unable to convert ancient PAR files: please use legacy dcm2nii\n"); @@ -1519,6 +1559,27 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D if ((strcmp(Comment[0], "Max.") == 0) && (strcmp(Comment[3], "slices/locations") == 0)) { d.xyzDim[3] = atoi(Comment[5]); } + if ((strcmp(Comment[0], "Max.") == 0) && (strcmp(Comment[3], "diffusion") == 0)) { + maxNumberOfDiffusionValues = atoi(Comment[6]); + } + if ((strcmp(Comment[0], "Max.") == 0) && (strcmp(Comment[3], "gradient") == 0)) { + maxNumberOfGradientOrients = atoi(Comment[6]); + } + if ((strcmp(Comment[0], "Max.") == 0) && (strcmp(Comment[3], "cardiac") == 0)) { + maxNumberOfCardiacPhases = atoi(Comment[6]); + } + if ((strcmp(Comment[0], "Max.") == 0) && (strcmp(Comment[3], "echoes") == 0)) { + maxNumberOfEchoes = atoi(Comment[5]); + if (maxNumberOfEchoes > 1) d.isMultiEcho = true; + } + if ((strcmp(Comment[0], "Max.") == 0) && (strcmp(Comment[3], "dynamics") == 0)) { + maxNumberOfDynamics = atoi(Comment[5]); + } + if ((strcmp(Comment[0], "Max.") == 0) && (strcmp(Comment[3], "mixes") == 0)) { + maxNumberOfMixes = atoi(Comment[5]); + if (maxNumberOfMixes > 1) + printError("maxNumberOfMixes > 1. Please update this software to support these images\n"); + } p = fgets (buff, LINESZ, fp);//get next line continue; } //process '.' tag @@ -1533,6 +1594,21 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D } for (int i = 0; i <= nCols; i++) cols[i] = strtof(p, &p); // p+1 skip comma, read a float + //the initial header does not report if phase or magnitude are include, so we read each separately in two passes + /* + if ((!isReadPhase) && (cols[kImageType] != 0)) { //skip phase images + p = fgets (buff, LINESZ, fp);//get next line + numSlice2D++; + d.isHasPhase = true; + continue; + } + if ((isReadPhase) && (cols[kImageType] == 0)) { //skip magnitude images + p = fgets (buff, LINESZ, fp);//get next line + numSlice2D++; + d.isHasMagnitude = true; + continue; + } + */ if ((cols[kIndex]) != slice) isIndexSequential = false; //slices 0,1,2.. should have indices 0,1,2,3... slice ++; if (slice == 1) { @@ -1564,11 +1640,12 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D printError("Slice dimensions or bit depth varies %s\n", parname); return d; } - if ((d.patientPositionSequentialRepeats == 0) && ((!isSameFloat(d.patientPosition[1],cols[kPositionRL])) || + //~ + /*if ((d.patientPositionSequentialRepeats == 0) && ((!isSameFloat(d.patientPosition[1],cols[kPositionRL])) || (!isSameFloat(d.patientPosition[2],cols[kPositionAP])) || (!isSameFloat(d.patientPosition[3],cols[kPositionFH])) ) )//this is the first slice with different position d.patientPositionSequentialRepeats = slice-1; - + */ if ((d.intenScale != cols[kRS]) || (d.intenIntercept != cols[kRI])) isIntenScaleVaries = true; } @@ -1596,14 +1673,44 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D d.patientPosition[2] = cols[kPositionAP]; d.patientPosition[3] = cols[kPositionFH]; } - if (d.patientPositionNumPhilips < kMaxDTI4D) { - dti4D->S[d.patientPositionNumPhilips].sliceNumberMrPhilips = round(cols[kSlice]); - if ((d.patientPositionNumPhilips > 0) && (dti4D->S[d.patientPositionNumPhilips].sliceNumberMrPhilips < dti4D->S[d.patientPositionNumPhilips-1].sliceNumberMrPhilips)) { - d.isSlicesSpatiallySequentialPhilips = false; - //printMessage("slices are not contiguous\n"); - } + patientPositionNumPhilips++; + } + if (true) { //for every slice + int slice = round(cols[kSlice]); + //(cols[kEcho] == 1) && (cols[kDyn] == 2) && (patientPositionNumPhilipsVol2 < kMaxDTI4D) && (cols[kCardiac] == 1) && (cols[kGradientNumber] == 1) + int volStep = maxNumberOfDynamics; + int vol = ((int)cols[kDyn] - 1); + int gradDynVol = (int)cols[kGradientNumber] - 1; + if (vol > gradDynVol) gradDynVol = vol; //if fMRI series, else DWI + vol = vol + (volStep * ((int)cols[kGradientNumber] - 1)); + volStep = volStep * maxNumberOfGradientOrients; + vol = vol + (volStep * ((int)cols[kEcho] - 1)); + volStep = volStep * maxNumberOfEchoes; + vol = vol + (volStep * ((int)cols[kCardiac] - 1)); + volStep = volStep * maxNumberOfCardiacPhases; + if (cols[kImageType] != 0) //phase maps are not reported in the initial header! + vol = vol + volStep; + if (vol > maxVol) maxVol = vol; + if (vol <= kMaxDTI4D) { + // dti4D->S[vol].V[0] = cols[kbval]; + //dti4D->gradDynVol[vol] = gradDynVol; + dti4D->TE[vol] = cols[kTEcho]; + dti4D->intenIntercept[vol] = cols[kRI]; + dti4D->intenScale[vol] = cols[kRS]; + dti4D->intenScalePhilips[vol] = cols[kSS]; + dti4D->isPhase[vol] = (cols[kImageType] != 0); + } + //if (slice == 1) printWarning("%d\n", (int)cols[kEcho]); + slice = slice + (vol * d.xyzDim[3]); + if (slice != (prevSlice + 1)) isSlicesSpatiallySequentialPhilips = false; + prevSlice = slice; + //if (cols[kImageType] != 0) //yikes - phase maps! + // slice = slice + numExpected; + //printWarning("%d\t%d\n", slice -1, numSlice2D); + if ((slice >= 0) && (slice < kMaxSlice2D) && (numSlice2D < kMaxSlice2D) && (numSlice2D >= 0)) { + dti4D->sliceOrder[slice -1] = numSlice2D; } - d.patientPositionNumPhilips++; + numSlice2D++; } if (cols[kGradientNumber] > 0) { /*int dir = (int) cols[kGradientNumber]; @@ -1613,11 +1720,11 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D //it seems like the ADC is always saved as the final volume, so this solution SHOULD be foolproof. ADCwarning = true; }*/ - //666: test (cols[kDyn] == 1) //(cols[kImageType] == 0) means magnitude scan if ((cols[kImageType] == 0) && (cols[kDyn] == 1) && (cols[kEcho] == 1) && (cols[kCardiac] == 1) && (cols[kSlice] == 1)) { //only first slice d.CSA.numDti++; - int dir = d.CSA.numDti; + //int dir = d.CSA.numDti; + int dir =(int)cols[kGradientNumber]; if (dir <= kMaxDTI4D) { if (isVerbose ) { if (d.CSA.numDti == 1) printMessage("n\tdir\tbValue\tV1\tV2\tV3\n"); @@ -1645,21 +1752,43 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D for (int s = 0; s < patientPositionNumPhilipsVol2; s++) sliceNumberMrPhilipsVol2[s] = sliceNumberMrPhilipsB2[s]; } - if ((patientPositionNumPhilipsVol2 > 1) && (d.patientPositionNumPhilips == patientPositionNumPhilipsVol2)) { + //even is sequential, we need to decode if multiple echoes/phasemaps,etc + //if (!isSlicesSpatiallySequentialPhilips) { + //printWarning("PAR file order of slices is not sequential: please ensure correct slice re-ordering\n"); + d.isScaleOrTEVaries = true; + if (numSlice2D > kMaxSlice2D) { + printError("Overloaded slice re-ordering. Number of slices (%d) exceeds kMaxSlice2D (%d)\n", numSlice2D, kMaxSlice2D); + dti4D->sliceOrder[0] = -1; + } + //} else + // dti4D->sliceOrder[0] = -1; //sequential slices do not need to be re-ordered + if ((d.isHasPhase) && (d.isHasMagnitude)) { + numExpected *= 2; //the initial header does not report that both types are stored + } + if (numSlice2D != numExpected) { + printMessage(" found %d slices, but expected %d: slices*diff*grad*cardiac*echo*dynamic*mix = %d*%d*%d*%d*%d*%d*%d\n", numSlice2D, numExpected, + d.xyzDim[3], maxNumberOfDiffusionValues, maxNumberOfGradientOrients, + maxNumberOfCardiacPhases, maxNumberOfEchoes, maxNumberOfDynamics, maxNumberOfMixes); + d.isValid = false; + } + + //~ + /*if ((patientPositionNumPhilipsVol2 > 1) && (d.patientPositionNumPhilips == patientPositionNumPhilipsVol2)) { bool isSliceOrderConsistent = true; for (int s = 0; s < patientPositionNumPhilipsVol2; s++) if (sliceNumberMrPhilipsVol2[s] != dti4D->S[s].sliceNumberMrPhilips) isSliceOrderConsistent = false; if (!isSliceOrderConsistent) printError("PAR file order of slices varies between volumes (hint: use dicm2nii)\n"); d.isValid = false; - } + }*/ if (dynNotAscending) { - printError("PAR file volumes not saved in ascending order (hint: use dicm2nii)\n"); - d.isValid = false; + printWarning("PAR file volumes not saved in ascending order (please check re-ordering)\n"); + //d.isValid = false; } if ((slice % d.xyzDim[3]) != 0) { - printError("Total number of slices (%d) not divisible by slices per 3D volume (%d) [acquisition aborted]. Try dicm2nii or nii_rescue_par to fix this: %s\n", slice, d.xyzDim[3], parname); + printError("Total number of slices (%d) not divisible by slices per 3D volume (%d) [acquisition aborted]. Try dicm2nii: %s\n", slice, d.xyzDim[3], parname); d.isValid = false; + return d; } d.xyzDim[4] = slice/d.xyzDim[3]; d.locationsInAcquisition = d.xyzDim[3]; @@ -1756,6 +1885,18 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D printWarning("Reported TR=%gms, measured TR=%gms (prospect. motion corr.?)\n", d.TR, TRms); d.TR = TRms; } + //check if dimensions vary + if (maxVol > 0) { //maxVol indexed from 0 + for (int i = 1; i <= maxVol; i++) { + //if (dti4D->gradDynVol[i] > d.maxGradDynVol) d.maxGradDynVol = dti4D->gradDynVol[i]; + if (dti4D->intenIntercept[i] != dti4D->intenIntercept[0]) d.isScaleOrTEVaries = true; + if (dti4D->intenScale[i] != dti4D->intenScale[0]) d.isScaleOrTEVaries = true; + if (dti4D->intenScalePhilips[i] != dti4D->intenScalePhilips[0]) d.isScaleOrTEVaries = true; + if (dti4D->isPhase[i] != dti4D->isPhase[0]) d.isScaleOrTEVaries = true; + } + if (d.isScaleOrTEVaries) + printWarning("Varying dimensions (echoes, phase maps, intensity scaling) will require volumes to be saved separately (hint: you may prefer dicm2nii output)\n"); + } //check DTI makes sense if (d.CSA.numDti > 1) { bool v1varies = false; @@ -1769,7 +1910,8 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D if ((!v1varies) || (!v2varies) || (!v3varies)) printError("Bizarre b-vectors %s\n", parname); } - if ((maxEcho > 1) || (maxCardiac > 1)) printWarning("Multiple Echo (%d) or Cardiac (%d). Segment output, e.g. nii_segment4d('img.nii', %d)\n", maxEcho, maxCardiac, maxEcho*maxCardiac); + if ((maxEcho > 1) || (maxCardiac > 1)) printWarning("Multiple Echo (%d) or Cardiac (%d). Carefully inspect output\n", maxEcho, maxCardiac); + if ((maxEcho > 1) || (maxCardiac > 1)) d.isScaleOrTEVaries = true; return d; } //nii_readParRec() @@ -1908,15 +2050,16 @@ unsigned char * nii_flipImgZ(unsigned char* bImg, struct nifti_1_header *hdr){ return bImg; } // nii_flipImgZ() -unsigned char * nii_reorderSlices(unsigned char* bImg, struct nifti_1_header *h, struct TDTI4D *dti4D){ +/*unsigned char * nii_reorderSlices(unsigned char* bImg, struct nifti_1_header *h, struct TDTI4D *dti4D){ //flip slice order - Philips scanners can save data in non-contiguous order //if ((h->dim[3] < 2) || (h->dim[4] > 1)) return bImg; + return bImg; if (h->dim[3] < 2) return bImg; if (h->dim[3] >= kMaxDTI4D) { printWarning("Unable to reorder slices (%d > %d)\n", h->dim[3], kMaxDTI4D); return bImg; } - //printMessage("<<< Slices not spatially contiguous: please check output [new feature]\n"); + printError("OBSOLETE<<< Slices not spatially contiguous: please check output [new feature]\n"); return bImg; int dim4to7 = 1; for (int i = 4; i < 8; i++) if (h->dim[i] > 1) dim4to7 = dim4to7 * h->dim[i]; @@ -1924,11 +2067,16 @@ unsigned char * nii_reorderSlices(unsigned char* bImg, struct nifti_1_header *h, if (sliceBytes < 0) return bImg; size_t volBytes = sliceBytes * h->dim[3]; unsigned char *srcImg = (unsigned char *)malloc(volBytes); + //printMessage("Reordering %d volumes\n", dim4to7); for (int v = 0; v < dim4to7; v++) { + //for (int v = 0; v < 1; v++) { + size_t volStart = v * volBytes; memcpy(&srcImg[0], &bImg[volStart], volBytes); //dest, src, size for (int z = 0; z < h->dim[3]; z++) { //for each slice int src = dti4D->S[z].sliceNumberMrPhilips - 1; //-1 as Philips indexes slices from 1 not 0 + if ((v > 0) && (dti4D->S[0].sliceNumberMrPhilipsVol2 >= 0)) + src = dti4D->S[z].sliceNumberMrPhilipsVol2 - 1; //printMessage("Reordering volume %d slice %d\n", v, dti4D->S[z].sliceNumberMrPhilips); if ((src < 0) || (src >= h->dim[3])) continue; memcpy(&bImg[volStart+(src*sliceBytes)], &srcImg[z*sliceBytes], sliceBytes); //dest, src, size @@ -1937,7 +2085,7 @@ unsigned char * nii_reorderSlices(unsigned char* bImg, struct nifti_1_header *h, free(srcImg); return bImg; }// nii_reorderSlices() - +*/ unsigned char * nii_flipZ(unsigned char* bImg, struct nifti_1_header *h){ //flip slice order if (h->dim[3] < 2) return bImg; @@ -2159,7 +2307,92 @@ unsigned char * nii_iVaries(unsigned char *img, struct nifti_1_header *hdr){ return (unsigned char*) img32; } //nii_iVaries() -unsigned char * nii_XYTZ_XYZT(unsigned char* bImg, struct nifti_1_header *hdr, int seqRepeats) { + +/*unsigned char * nii_reorderSlicesX(unsigned char* bImg, struct nifti_1_header *hdr, struct TDTI4D *dti4D) { +//Philips can save slices in any random order... rearrange all of them + int dim3to7 = 1; + for (int i = 3; i < 8; i++) + if (hdr->dim[i] > 1) dim3to7 = dim3to7 * hdr->dim[i]; + if (dim3to7 < 2) return bImg; + uint64_t sliceBytes = hdr->dim[1]*hdr->dim[2]*hdr->bitpix/8; + unsigned char *sliceImg = (unsigned char *)malloc( sliceBytes); + uint32_t *idx = (uint32_t *)malloc( dim3to7 * sizeof(uint32_t)); + for (int i = 0; i < dim3to7; i++) //for each volume + idx[i] = i; + for (int i = 0; i < dim3to7; i++) { //for each volume + int fromSlice = idx[dti4D->sliceOrder[i]]; + //if (i < 10) printMessage(" ===> Moving slice from/to positions\t%d\t%d\n", i, toSlice); + if (i != fromSlice) { + uint64_t inPos = fromSlice * sliceBytes; + uint64_t outPos = i * sliceBytes; + memcpy( &sliceImg[0], &bImg[outPos], sliceBytes); //dest, src -> copy slice about to be overwritten + memcpy( &bImg[outPos], &bImg[inPos], sliceBytes); //dest, src + memcpy( &bImg[inPos], &sliceImg[0], sliceBytes); // + idx[i] = fromSlice; + } + } + free(idx); + free(sliceImg); + return bImg; +}*/ + +unsigned char * nii_reorderSlicesX(unsigned char* bImg, struct nifti_1_header *hdr, struct TDTI4D *dti4D) { +//Philips can save slices in any random order... rearrange all of them + int dim3to7 = 1; + for (int i = 3; i < 8; i++) + if (hdr->dim[i] > 1) dim3to7 = dim3to7 * hdr->dim[i]; + if (dim3to7 < 2) return bImg; + if (dim3to7 > kMaxSlice2D) return bImg; + uint64_t imgSz = nii_ImgBytes(*hdr); + uint64_t sliceBytes = hdr->dim[1]*hdr->dim[2]*hdr->bitpix/8; + unsigned char *outImg = (unsigned char *)malloc( imgSz); + memcpy( &outImg[0],&bImg[0], imgSz); + for (int i = 0; i < dim3to7; i++) { //for each volume + int fromSlice = dti4D->sliceOrder[i]; + //if (i < 10) printMessage(" ===> Moving slice from/to positions\t%d\t%d\n", i, toSlice); + if ((i < 0) || (fromSlice >= dim3to7)) + printError("Re-ordered slice out-of-volume %d\n", fromSlice); + else if (i != fromSlice) { + uint64_t inPos = fromSlice * sliceBytes; + uint64_t outPos = i * sliceBytes; + memcpy( &bImg[outPos], &outImg[inPos], sliceBytes); + } + } + free(outImg); + return bImg; +} + +/*unsigned char * nii_reorderSlicesX(unsigned char* bImg, struct nifti_1_header *hdr, struct TDTI4D *dti4D) { +//Philips can save slices in any random order... rearrange all of them + //if (sliceOrder == NULL) return bImg; + int dim3to7 = 1; + for (int i = 3; i < 8; i++) + if (hdr->dim[i] > 1) dim3to7 = dim3to7 * hdr->dim[i]; + if (dim3to7 < 2) return bImg; + //printMessage("xxx NOT reordering %d Philips slices.\n", dim3to7); return bImg; + uint64_t sliceBytes = hdr->dim[1]*hdr->dim[2]*hdr->bitpix/8; + unsigned char *sliceImg = (unsigned char *)malloc( sliceBytes); + //for (int i = 0; i < dim3to7; i++) { //for each volume + uint64_t imgSz = nii_ImgBytes(*hdr); + //this uses a lot of RAM, someday this could be done in place... + unsigned char *outImg = (unsigned char *)malloc( imgSz); + memcpy( &outImg[0],&bImg[0], imgSz); + + for (int i = 0; i < dim3to7; i++) { //for each volume + int toSlice = dti4D->sliceOrder[i]; + //if (i < 10) printMessage(" ===> Moving slice from/to positions\t%d\t%d\n", i, toSlice); + if (i != toSlice) { + uint64_t inPos = i * sliceBytes; + uint64_t outPos = toSlice * sliceBytes; + memcpy( &bImg[outPos], &outImg[inPos], sliceBytes); + } + } + free(sliceImg); + free(outImg); + return bImg; +}*/ + +/*unsigned char * nii_XYTZ_XYZT(unsigned char* bImg, struct nifti_1_header *hdr, int seqRepeats) { //Philips can save time as 3rd dimensions, NIFTI requires time is 4th dimension int dim4to7 = 1; for (int i = 4; i < 8; i++) @@ -2193,6 +2426,7 @@ unsigned char * nii_XYTZ_XYZT(unsigned char* bImg, struct nifti_1_header *hdr, i free(bImg); return outImg; } //nii_XYTZ_XYZT() +*/ unsigned char * nii_byteswap(unsigned char *img, struct nifti_1_header *hdr){ if (hdr->bitpix < 9) return img; @@ -2645,7 +2879,7 @@ unsigned char * nii_loadImgRLE(char* imgname, struct nifti_1_header hdr, struct #endif #endif -unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct TDICOMdata dcm, bool iVaries, int compressFlag, int isVerbose) { +unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct TDICOMdata dcm, bool iVaries, int compressFlag, int isVerbose, struct TDTI4D *dti4D) { //provided with a filename (imgname) and DICOM header (dcm), creates NIfTI header (hdr) and img if (headerDcm2Nii(dcm, hdr, true) == EXIT_FAILURE) return NULL; //TOFU unsigned char * img; @@ -2702,8 +2936,17 @@ unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct hdr->dim[3] = nAcq; hdr->dim[0] = 4; } - if ((hdr->dim[0] > 3) && (dcm.patientPositionSequentialRepeats > 1)) //swizzle 3rd and 4th dimension (Philips stores time as 3rd dimension) - img = nii_XYTZ_XYZT(img, hdr,dcm.patientPositionSequentialRepeats ); + //~ if ((hdr->dim[0] > 3) && (dcm.patientPositionSequentialRepeats > 1) && (dcm.sliceOrder == NULL)) //swizzle 3rd and 4th dimension (Philips stores time as 3rd dimension) + //~ img = nii_XYTZ_XYZT(img, hdr,dcm.patientPositionSequentialRepeats); + if ((dti4D != NULL) && (dti4D->sliceOrder[0] >= 0)) + img = nii_reorderSlicesX(img, hdr, dti4D); + //~ + /*if (((dcm.patientPositionSequentialRepeats * 2) == dcm.patientPositionRepeats) && (dcm.isHasPhase) && (dcm.isHasMagnitude)) { + hdr->dim[3] = hdr->dim[3] / 2; + hdr->dim[4] = hdr->dim[4] * 2; + hdr->dim[0] = 4; + printMessage("Splitting Phase+Magnitude into two volumes for %d slices (Z) and %d volumes (T).\n",hdr->dim[3], hdr->dim[4]); + }*/ headerDcm2NiiSForm(dcm,dcm, hdr, false); return img; } //nii_loadImgXL() @@ -2719,6 +2962,12 @@ int isSQ(uint32_t groupElement) { //Detect sequence VR ("SQ") for implicit tags return 0; } //isSQ() +int isSameFloatGE (float a, float b) { +//Kludge for bug in 0002,0016="DIGITAL_JACKET", 0008,0070="GE MEDICAL SYSTEMS" DICOM data: Orient field (0020:0037) can vary 0.00604261 == 0.00604273 !!! + //return (a == b); //niave approach does not have any tolerance for rounding errors + return (fabs (a - b) <= 0.0001); +} + int isDICOMfile(const char * fname) { //0=NotDICOM, 1=DICOM, 2=Maybe(not Part 10 compliant) FILE *fp = fopen(fname, "rb"); if (!fp) return 0; @@ -2930,14 +3179,38 @@ void _update_tvd(struct TVolumeDiffusion* ptvd) { }//_update_tvd() //END RIR +struct TDCMdim { //DimensionIndexValues + uint32_t dimIdx[MAX_NUMBER_OF_DIMENSIONS]; + uint32_t diskPos; + float TE, intenScale, intenIntercept, intenScalePhilips; + bool isPhase; +}; + +int compareTDCMdim(void const *item1, void const *item2) { + struct TDCMdim const *dcm1 = (const struct TDCMdim *)item1; + struct TDCMdim const *dcm2 = (const struct TDCMdim *)item2; + //for(int i=0; i < MAX_NUMBER_OF_DIMENSIONS; ++i){ + for(int i=MAX_NUMBER_OF_DIMENSIONS-1; i >=0; i--){ + + if(dcm1->dimIdx[i] < dcm2->dimIdx[i]) + return -1; + else if(dcm1->dimIdx[i] > dcm2->dimIdx[i]) + return 1; + } + return 0; +} //compareTDCMdim() + struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, struct TDTI4D *dti4D) { struct TDICOMdata d = clear_dicom_data(); + d.imageNum = 0; //not set strcpy(d.protocolName, ""); //erase dummy with empty strcpy(d.protocolName, ""); //erase dummy with empty strcpy(d.seriesDescription, ""); //erase dummy with empty 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->S[0].sliceNumberMrPhilips = -1; + //dti4D->S[0].sliceNumberMrPhilipsVol2 = -1; + dti4D->sliceOrder[0] = -1; struct TVolumeDiffusion volDiffusion = initTVolumeDiffusion(&d, dti4D); struct stat s; @@ -3043,6 +3316,8 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kZThick 0x0018+(0x0050 << 16 ) #define kTR 0x0018+(0x0080 << 16 ) #define kTE 0x0018+(0x0081 << 16 ) +//#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 @@ -3070,6 +3345,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kDiffusionOrientation 0x0018+uint32_t(0x9089<< 16 ) // FD, seen in enhanced // DICOM from Philips 5.* // and Siemens XA10. +//#define kMREchoSequence 0x0018+uint32_t(0x9114<< 16 ) //SQ #define kDwellTime 0x0019+(0x1018<< 16 ) //IS in NSec, see https://github.com/rordenlab/dcm2niix/issues/127 #define kLastScanLoc 0x0019+(0x101B<< 16 ) #define kDiffusionDirectionGEX 0x0019+(0x10BB<< 16 ) //DS @@ -3084,13 +3360,16 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kImageNum 0x0020+(0x0013 << 16 ) #define kStudyInstanceUID 0x0020+(0x000D << 16 ) #define kSeriesInstanceUID 0x0020+(0x000E << 16 ) -#define kPatientPosition 0x0020+(0x0032 << 16 ) // Actually ImagePositionPatient! +#define kImagePositionPatient 0x0020+(0x0032 << 16 ) // Actually ! #define kOrientationACR 0x0020+(0x0035 << 16 ) +#define kTemporalPositionIdentifier 0x0020+(0x0100 << 16 ) //IS #define kOrientation 0x0020+(0x0037 << 16 ) #define kImagesInAcquisition 0x0020+(0x1002 << 16 ) //IS #define kImageComments 0x0020+(0x4000<< 16 )// '0020' '4000' 'LT' 'ImageComments' #define kDimensionIndexValues 0x0020+uint32_t(0x9157<< 16 ) // UL n-dimensional index of frame. +#define kInStackPositionNumber 0x0020+uint32_t(0x9057<< 16 ) // UL can help determine slices in volume #define kLocationsInAcquisitionGE 0x0021+(0x104F<< 16 )// 'SS' 'LocationsInAcquisitionGE' +#define kProtocolDataBlockGE 0x0025+(0x101B<< 16 )// 'OB' #define kSamplesPerPixel 0x0028+(0x0002 << 16 ) #define kPhotometricInterpretation 0x0028+(0x0004 << 16 ) #define kPlanarRGB 0x0028+(0x0006 << 16 ) @@ -3122,9 +3401,9 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kIconImageSequence 0x0088+(0x0200 << 16 ) #define kDiffusionBFactor 0x2001+(0x1003 << 16 )// FL #define kSliceNumberMrPhilips 0x2001+(0x100A << 16 ) //IS Slice_Number_MR -#define kNumberOfSlicesMrPhilips 0x2001+(0x1018 << 16 )//SL 0x2001, 0x1018 ), "Number_of_Slices_MR" #define kSliceOrient 0x2001+(0x100B << 16 )//2001,100B Philips slice orientation (TRANSVERSAL, AXIAL, SAGITTAL) -//#define kLocationsInAcquisitionPhilips 0x2001+(0x1018 << 16 ) // +#define kNumberOfSlicesMrPhilips 0x2001+(0x1018 << 16 )//SL 0x2001, 0x1018 ), "Number_of_Slices_MR" +//#define kNumberOfLocationsPhilips 0x2001+(0x1015 << 16 ) //SS //#define kStackSliceNumber 0x2001+(0x1035 << 16 )//? Potential way to determine slice order for Philips? #define kNumberOfDynamicScans 0x2001+(0x1081 << 16 )//'2001' '1081' 'IS' 'NumberOfDynamicScans' #define kMRAcquisitionTypePhilips 0x2005+(0x106F << 16) @@ -3139,21 +3418,35 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kDiffusionDirectionAP 0x2005+(0x10B1 << 16) #define kDiffusionDirectionFH 0x2005+(0x10B2 << 16) #define kPrivatePerFrameSq 0x2005+(0x140F << 16) +#define kMRImageGradientOrientationNumber 0x2005+(0x1413 << 16) //IS #define kWaveformSq 0x5400+(0x0100 << 16) #define kImageStart 0x7FE0+(0x0010 << 16 ) #define kImageStartFloat 0x7FE0+(0x0008 << 16 ) #define kImageStartDouble 0x7FE0+(0x0009 << 16 ) -uint32_t kNest = 0xFFFE +(0xE000 << 16 ); //#define kNest 0xFFFE +(0xE000 << 16 ) //Item follows SQ -uint32_t kUnnest = 0xFFFE +(0xE00D << 16 ); //#define kUnnest 0xFFFE +(0xE00D << 16 ) //ItemDelimitationItem [length defined] http://www.dabsoft.ch/dicom/5/7.5/ -uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD << 16 )//SequenceDelimitationItem [length undefined] - int nest = 0; - //double zSpacing = -1.0l; //includes slice thickness plus gap +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 + bool is2005140FSQ = false; int locationsInAcquisitionGE = 0; + int inStackPositionNumber = 0; + int maxInStackPositionNumber = 0; + int temporalPositionIdentifier = 0; int locationsInAcquisitionPhilips = 0; int imagesInAcquisition = 0; + //int sumSliceNumberMrPhilips = 0; + int sliceNumberMrPhilips = 0; + int numberOfFrames = 0; + int MRImageGradientOrientationNumber = 0; + //int maxMRImageGradientOrientationNumber = -1; + int minGradNum = kMaxDTI4D + 1; + int maxGradNum = -1; + int numberOfDynamicScans = 0; uint32_t lLength; uint32_t groupElement; long lPos = 0; + bool isPhilipsDerived = false; + //bool isPhilipsDiffusion = false; if (isPart10prefix) { //for part 10 files, skip preamble and prefix lPos = 128+4; //4-byte signature starts at 128 groupElement = buffer[lPos] | (buffer[lPos+1] << 8) | (buffer[lPos+2] << 16) | (buffer[lPos+3] << 24); @@ -3172,18 +3465,41 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD bool isIconImageSequence = false; bool isSwitchToImplicitVR = false; bool isSwitchToBigEndian = false; - //bool geiisBug = false; //for buggy GEIIS http://forum.dcmtk.org/viewtopic.php?p=7162&sid=3b516cc751aae51fbb5e73184abe37c2 - bool is2005140FSQ = false; //for buggy Philips - bool is2005140FSQwarned = false; //for buggy Philips bool isAtFirstPatientPosition = false; //for 3d and 4d files: flag is true for slices at same position as first slice bool isMosaic = false; int patientPositionNum = 0; - int sqDepth = 0; + float B0Philips = -1.0; + float vRLPhilips = 0.0; + float vAPPhilips = 0.0; + float vFHPhilips = 0.0; + bool isPhase = false; + bool isMagnitude = false; + float patientPositionPrivate[4] = {NAN, NAN, NAN, NAN}; float patientPosition[4] = {NAN, NAN, NAN, NAN}; //used to compute slice direction for Philips 4D + //float patientPositionPublic[4] = {NAN, NAN, NAN, NAN}; //used to compute slice direction for Philips 4D float patientPositionEndPhilips[4] = {NAN, NAN, NAN, NAN}; float patientPositionStartPhilips[4] = {NAN, NAN, NAN, NAN}; +/* //not required? + int sqDepth = 0; + int sqDepthPrivate = 0; + int sqEndPrivate = -1; //used to skip private SQs that provide explicit length + bool is2005140FSQ = false; //Philips stores 2D slice data here + bool is2005140FSQwarned = false; //for buggy Philips - while ((d.imageStart == 0) && ((lPos+8+lFileOffset) < fileLen)) { +*/ + //array for storing Philips DTI when tags kMRImageGradientOrientationNumber and kDimensionIndexValues are available + struct TDTI philDTI[kMaxDTI4D]; + for (int i = 0; i < kMaxDTI4D; i++) + philDTI[i].V[0] = -1; + //array for storing DimensionIndexValues + int numDimensionIndexValues = 0; + TDCMdim dcmDim[kMaxSlice2D]; + for (int i = 0; i < kMaxSlice2D; i++) { + dcmDim[i].diskPos = i; + for (int j = 0; j < MAX_NUMBER_OF_DIMENSIONS; j++) + dcmDim[i].dimIdx[j] = 0; + } + while ((d.imageStart == 0) && ((lPos+8+lFileOffset) < fileLen)) { #ifndef myLoadWholeFileToReadHeader //read one segment at a time if ((size_t)(lPos + 128) > MaxBufferSz) { //avoid overreading the file lFileOffset = lFileOffset + lPos; @@ -3213,19 +3529,11 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD } //transfer syntax requests switching VR after group 0001 //uint32_t group = (groupElement & 0xFFFF); lPos += 4; - if ((groupElement == kUnnest) || (groupElement == kUnnest2)) isIconImageSequence = false; - if (((groupElement == kNest) || (groupElement == kUnnest) || (groupElement == kUnnest2)) && (!isEncapsulatedData)) { - //if (((groupElement == kNest) || (groupElement == kUnnest) || (groupElement == kUnnest2)) ) { + if ((groupElement == kItemDelimitationTag) || (groupElement == kSequenceDelimitationItemTag)) isIconImageSequence = false; + if (((groupElement == kItemTag) || (groupElement == kItemDelimitationTag) || (groupElement == kSequenceDelimitationItemTag)) && (!isEncapsulatedData)) { vr[0] = 'N'; vr[1] = 'A'; - if (groupElement == kUnnest2) sqDepth--; - //if (groupElement == kUnnest) geiisBug = false; //don't exit if there is a proprietary thumbnail lLength = 4; - //Next: Skip nested items if explicit length provided - //int nestLength = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); - //if (nestLength < 0xFFFFFFFF) { - // printMessage(" >>>>> %d\n", nestLength); - //} } else if (d.isExplicitVR) { vr[0] = buffer[lPos]; vr[1] = buffer[lPos+1]; if (buffer[lPos+1] < 'A') {//implicit vr with 32-bit length @@ -3248,7 +3556,6 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD } else if ((buffer[lPos] == 'S') && (buffer[lPos+1] == 'Q')) { lLength = 8; //Sequence Tag //printMessage(" !!!SQ\t%04x,%04x\n", groupElement & 65535,groupElement>>16); - is2005140FSQ = (groupElement == kPrivatePerFrameSq); } else { //explicit VR with 16-bit length if ((d.isLittleEndian) ) lLength = buffer[lPos+2] | (buffer[lPos+3] << 8); @@ -3257,8 +3564,8 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD lPos += 4; //skip 2 byte VR string and 2 length bytes = 4 bytes } } else { //implicit VR - vr[0] = 'N'; - vr[1] = 'A'; + vr[0] = 'U'; + vr[1] = 'N'; if (d.isLittleEndian) lLength = buffer[lPos] | (buffer[lPos+1] << 8) | (buffer[lPos+2] << 16) | (buffer[lPos+3] << 24); else @@ -3268,26 +3575,65 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD vr[0] = 'S'; vr[1] = 'Q'; lLength = 8; //Sequence Tag - is2005140FSQ = (groupElement == kPrivatePerFrameSq); - //if (is2005140FSQ) printMessage(" !!!SQ\t%04x,%04x\n", groupElement & 65535,groupElement>>16); } } //if explicit else implicit VR if (lLength == 0xFFFFFFFF) { lLength = 8; //SQ (Sequences) use 0xFFFFFFFF [4294967295] to denote unknown length - vr[0] = 'S'; - vr[1] = 'Q'; + //09032018 - do not count these as SQs: Horos does not count even groups + //uint32_t special = dcmInt(4,&buffer[lPos],d.isLittleEndian); + + //http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_7.5.html + + //if (special != ksqDelim) { + vr[0] = 'S'; + vr[1] = 'Q'; + //} } + /* //Handle SQs: for explicit these have VR=SQ if ((vr[0] == 'S') && (vr[1] == 'Q')) { - sqDepth++; - //printMessage("SQstart %d\n", sqDepth); + //http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_7.5.html + uint32_t special = dcmInt(4,&buffer[lPos],d.isLittleEndian); + uint32_t slen = dcmInt(4,&buffer[lPos+4],d.isLittleEndian); + //if (d.isExplicitVR) + // slen = dcmInt(4,&buffer[lPos+8],d.isLittleEndian); + uint32_t kUndefinedLen = 0xFFFFFFFF; + //printError(" SPECIAL >>>>t%04x,%04x %08x %08x\n", groupElement & 65535,groupElement>>16, special, slen); + //return d; + is2005140FSQ = (groupElement == kPrivatePerFrameSq); + //if (isNextSQis2005140FSQ) is2005140FSQ = true; + //isNextSQis2005140FSQ = false; + if (special == kSequenceDelimitationItemTag) { + //unknown + } else if (slen == kUndefinedLen) { + sqDepth++; + if ((sqDepthPrivate == 0) && ((groupElement & 65535) % 2)) + sqDepthPrivate = sqDepth; //in a private SQ: ignore contents + } else if ((is2005140FSQ) || ((groupElement & 65535) % 2)) {//private SQ of known length - lets jump over this! + slen = lFileOffset + lPos + slen; + if ((sqEndPrivate < 0) || (slen > sqEndPrivate)) //sqEndPrivate is signed + sqEndPrivate = slen; //if nested private SQs, remember the end address of the top parent SQ + } + } + //next: look for required tags + if ((groupElement == kItemTag) && (isEncapsulatedData)) { + d.imageBytes = dcmInt(4,&buffer[lPos-4],d.isLittleEndian); + //printMessage("compressed data %d-> %ld\n",d.imageBytes, lPos); + if (d.imageBytes > 128) { + encapsulatedDataFragments++; + if (encapsulatedDataFragmentStart == 0) + encapsulatedDataFragmentStart = (int)lPos + (int)lFileOffset; + } } - if ((groupElement == kNest) || ((vr[0] == 'S') && (vr[1] == 'Q'))) nest++; - //if ((groupElement == kUnnest) || (groupElement == kUnnest2)) { // <- ? - if (groupElement == kUnnest) { - nest--; + if ((sqEndPrivate > 0) && ((lFileOffset + lPos) > sqEndPrivate)) + sqEndPrivate = -1; //end of private SQ with defined length + if (groupElement == kSequenceDelimitationItemTag) { //end of private SQ with undefined length + sqDepth--; + if (sqDepth < sqDepthPrivate) { + sqDepthPrivate = 0; //no longer in a private SQ + } } - //next: look for required tags - if ((groupElement == kNest) && (isEncapsulatedData)) { + if (sqDepth < 0) sqDepth = 0;*/ + if ((groupElement == kItemTag) && (isEncapsulatedData)) { //use this to find image fragment for compressed datasets, e.g. JPEG transfer syntax d.imageBytes = dcmInt(4,&buffer[lPos-4],d.isLittleEndian); //printMessage("compressed data %d-> %ld\n",d.imageBytes, lPos); if (d.imageBytes > 128) { @@ -3297,20 +3643,11 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD } } if ((isIconImageSequence) && ((groupElement & 0x0028) == 0x0028 )) groupElement = kUnused; //ignore icon dimensions - if (true) { //(nest <= 0) { //some Philips images have different 0020,0013 - //verbose reporting : - // printMessage("Pos %ld GroupElement %#08x,%#08x Length %d isLittle %d\n", lPos, (groupElement & 0xFFFF), (groupElement >> 16), lLength, d.isLittleEndian); - - // // Debugging - // int groupItem = groupElement >> 16; - // int groupGroup = groupElement - (groupItem << 16); - // printMessage("groupElement: (%04X, %04X)\n", groupElement & 65535,groupElement>>16); - + if (groupElement == kSequenceDelimitationItemTag) is2005140FSQ = false; switch ( groupElement ) { case kTransferSyntax: { char transferSyntax[kDICOMStr]; dcmStr (lLength, &buffer[lPos], transferSyntax); - //printMessage("%d transfer syntax>>> '%s'\n", compressFlag, 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) { @@ -3420,8 +3757,15 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD break; case kComplexImageComponent: if (lLength < 2) break; - d.isHasPhase = (buffer[lPos]=='P') && (toupper(buffer[lPos+1]) == 'H'); - d.isHasMagnitude = (buffer[lPos]=='M') && (toupper(buffer[lPos+1]) == 'A'); + isPhase = false; + isMagnitude = false; + if ((buffer[lPos]=='P') && (toupper(buffer[lPos+1]) == 'H')) + isPhase = true; + if ((buffer[lPos]=='M') && (toupper(buffer[lPos+1]) == 'A')) + isMagnitude = true; + //not mutually exclusive: possible for Philips enhanced DICOM to store BOTH magnitude and phase in the same image + if (isPhase) d.isHasPhase = true; + if (isMagnitude) d.isHasMagnitude = true; break; case kAcquisitionTime : char acquisitionTimeTxt[kDICOMStr]; @@ -3491,9 +3835,14 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD case kPatientOrient : dcmStr (lLength, &buffer[lPos], d.patientOrient); break; - case kDiffusionDirectionality : // 0018, 9075 + case kDiffusionDirectionality : {// 0018, 9075 set_directionality0018_9075(&volDiffusion, (&buffer[lPos])); - break; + if ((d.manufacturer != kMANUFACTURER_PHILIPS) || (lLength < 10)) break; + char dir[kDICOMStr]; + dcmStr (lLength, &buffer[lPos], dir); + if (strcmp(dir, "ISOTROPIC") == 0) + isPhilipsDerived = true; + break; } case kDwellTime : d.dwellTime = dcmStrInt(lLength, &buffer[lPos]); break; @@ -3526,37 +3875,41 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD case kSeriesInstanceUID : // 0020, 000E dcmStr (lLength, &buffer[lPos], d.seriesInstanceUID); break; - case kPatientPosition : // 0020, 0032, ImagePositionPatient - if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (is2005140FSQ)) { - if (!is2005140FSQwarned) - printWarning("Philips R3.2.2 can report different positions for the same slice. Attempting patch.\n"); - is2005140FSQwarned = true; - } else { - patientPositionNum++; - isAtFirstPatientPosition = true; - 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 - for (int k = 0; k < 4; k++) - d.patientPosition[k] = patientPosition[k]; - } else { - //dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, &d.patientPositionLast[0]); //slice direction for 4D - for (int k = 0; k < 4; k++) - d.patientPositionLast[k] = patientPosition[k]; - if ((isFloatDiff(d.patientPositionLast[1],d.patientPosition[1])) || - (isFloatDiff(d.patientPositionLast[2],d.patientPosition[2])) || - (isFloatDiff(d.patientPositionLast[3],d.patientPosition[3])) ) { - isAtFirstPatientPosition = false; //this slice is not at position of 1st slice - if (d.patientPositionSequentialRepeats == 0) //this is the first slice with different position - d.patientPositionSequentialRepeats = patientPositionNum-1; - } //if different position from 1st slice in file - } //if not first slice in file - set_isAtFirstPatientPosition_tvd(&volDiffusion, isAtFirstPatientPosition); - if (isVerbose > 1) - 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]); - - } //not after 2005,140F - break; + case kImagePositionPatient : { + if (is2005140FSQ) { + dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, &patientPositionPrivate[0]); + break; + } + patientPositionNum++; + isAtFirstPatientPosition = true; + dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, &patientPosition[0]); //slice position + + /*if (!is2005140FSQ) { + printMessage("cxcPublic 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]); + for (int k = 0; k < 4; k++) + patientPositionPublic[k] = patientPosition[k]; + }*/ + + if (isnan(d.patientPosition[1])) { + //dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, &d.patientPosition[0]); //slice position + for (int k = 0; k < 4; k++) + d.patientPosition[k] = patientPosition[k]; + } else { + //dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, &d.patientPositionLast[0]); //slice direction for 4D + for (int k = 0; k < 4; k++) + d.patientPositionLast[k] = patientPosition[k]; + if ((isFloatDiff(d.patientPositionLast[1],d.patientPosition[1])) || + (isFloatDiff(d.patientPositionLast[2],d.patientPosition[2])) || + (isFloatDiff(d.patientPositionLast[3],d.patientPosition[3])) ) { + isAtFirstPatientPosition = false; //this slice is not at position of 1st slice + //if (d.patientPositionSequentialRepeats == 0) //this is the first slice with different position + // d.patientPositionSequentialRepeats = patientPositionNum-1; + } //if different position from 1st slice in file + } //if not first slice in file + set_isAtFirstPatientPosition_tvd(&volDiffusion, isAtFirstPatientPosition); + if (isVerbose > 1) + 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: d.phaseEncodingRC = toupper(buffer[lPos]); //first character is either 'R'ow or 'C'ol break; @@ -3574,13 +3927,18 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD d.acquNum = dcmStrInt(lLength, &buffer[lPos]); break; case kImageNum: - //int dx = 3; - if (d.imageNum <= 1) d.imageNum = dcmStrInt(lLength, &buffer[lPos]); //Philips renames each image as image1 in 2001,9000 - break; - case kDimensionIndexValues: // kImageNum is not enough for 4D series from Philips 5.*. - { // { necessary for initializing ndim. + //Enhanced Philips also uses this in once per file SQ 0008,1111 + //Enhanced Philips also uses this once per slice in SQ 2005,140f + if (d.imageNum < 1) d.imageNum = dcmStrInt(lLength, &buffer[lPos]); //Philips renames each image again in 2001,9000, which can lead to duplicates + break; + case kInStackPositionNumber: + if (d.manufacturer != kMANUFACTURER_PHILIPS) break; + inStackPositionNumber = dcmInt(4,&buffer[lPos],d.isLittleEndian); + if (inStackPositionNumber > maxInStackPositionNumber) maxInStackPositionNumber = inStackPositionNumber; + break; + case kDimensionIndexValues: { // kImageNum is not enough for 4D series from Philips 5.*. + if (lLength < 4) break; uint8_t ndim = lLength / 4; - if(ndim > MAX_NUMBER_OF_DIMENSIONS){ // Ideally degenerate axes would be cleverly handled. // They are commonly seen in NIfTIs, but perhaps not in DICOM. @@ -3589,20 +3947,112 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD ndim = MAX_NUMBER_OF_DIMENSIONS; // Truncate } dcmMultiLongs(4 * ndim, &buffer[lPos], ndim, d.dimensionIndexValues, d.isLittleEndian); + if (d.manufacturer != kMANUFACTURER_PHILIPS) break; + // + if (inStackPositionNumber > 0) { + //for images without SliceNumberMrPhilips (2001,100A) + int sliceNumber = inStackPositionNumber; + if ((sliceNumber == 1) && (!isnan(patientPosition[1])) ) { + for (int k = 0; k < 4; k++) + patientPositionStartPhilips[k] = patientPosition[k]; + } else if ((sliceNumber == 1) && (!isnan(patientPositionPrivate[1]))) { + for (int k = 0; k < 4; k++) + patientPositionStartPhilips[k] = patientPositionPrivate[k]; + } + if ((sliceNumber == maxInStackPositionNumber) && (!isnan(patientPosition[1])) ) { + for (int k = 0; k < 4; k++) + patientPositionEndPhilips[k] = patientPosition[k]; + } else if ((sliceNumber == maxInStackPositionNumber) && (!isnan(patientPositionPrivate[1])) ) { + for (int k = 0; k < 4; k++) + patientPositionEndPhilips[k] = patientPositionPrivate[k]; + } + patientPosition[1] = NAN; + patientPositionPrivate[1] = NAN; + } + inStackPositionNumber = 0; + if (numDimensionIndexValues >= kMaxSlice2D) { + printError("Too many slices to track dimensions. Only up to %d are supported\n", kMaxSlice2D); + break; + } + for (int i = 0; i < ndim; i++) + dcmDim[numDimensionIndexValues].dimIdx[i] = d.dimensionIndexValues[i]; + dcmDim[numDimensionIndexValues].TE = TE; + dcmDim[numDimensionIndexValues].intenScale = d.intenScale; + dcmDim[numDimensionIndexValues].intenIntercept = d.intenIntercept; + dcmDim[numDimensionIndexValues].isPhase = isPhase; + dcmDim[numDimensionIndexValues].intenScalePhilips = d.intenScalePhilips; + numDimensionIndexValues ++; + //next: add diffusion if reported + if (B0Philips < 0.0) break; //diffusion parameters not reported + // Philips does not always provide 2005,1413 (MRImageGradientOrientationNumber) and sometimes after dimensionIndexValues + int gradNum = 0; + for (int i = 0; i < ndim; i++) + if (d.dimensionIndexValues[i] > 0) gradNum = d.dimensionIndexValues[i]; + if (gradNum <= 0) break; + /* + With Philips 51.0 both ADC and B=0 are saved as same direction, though they vary in another dimension + (0018,9075) CS [ISOTROPIC] + (0020,9157) UL 1\2\1\33 << ADC MAP + (0018,9075) CS [NONE] + (0020,9157) UL 1\1\2\33 + next two lines attempt to skip ADC maps + we could also increment gradNum for ADC if we wanted... + */ + if (isPhilipsDerived) { + gradNum ++; + B0Philips = 2000.0; + vRLPhilips = 0.0; + vAPPhilips = 0.0; + vFHPhilips = 0.0; + + } + if (B0Philips == 0.0) { + //printMessage(" DimensionIndexValues grad %d b=%g vec=%gx%gx%g\n", gradNum, B0Philips, vRLPhilips, vAPPhilips, vFHPhilips); + vRLPhilips = 0.0; + vAPPhilips = 0.0; + vFHPhilips = 0.0; + } + //if ((MRImageGradientOrientationNumber > 0) && ((gradNum != MRImageGradientOrientationNumber)) break; + if (gradNum < minGradNum) minGradNum = gradNum; + if (gradNum >= maxGradNum) maxGradNum = gradNum; + if (gradNum >= kMaxDTI4D) { + printError("Number of DTI gradients exceeds 'kMaxDTI4D (%d).\n", kMaxDTI4D); + break; + } + gradNum = gradNum - 1;//index from 0 + philDTI[gradNum].V[0] = B0Philips; + philDTI[gradNum].V[1] = vRLPhilips; + philDTI[gradNum].V[2] = vAPPhilips; + philDTI[gradNum].V[3] = vFHPhilips; + isPhilipsDerived = false; + //printMessage(" DimensionIndexValues grad %d b=%g vec=%gx%gx%g\n", gradNum, B0Philips, vRLPhilips, vAPPhilips, vFHPhilips); + + + /*if ((sliceNumberMrPhilips == locationsInAcquisitionPhilips) || ((locationsInAcquisitionPhilips == 0) &&(sliceNumberMrPhilips == maxInStackPositionNumber)) ) { + for (int k = 0; k < 4; k++) + patientPositionEndPhilips[k] = patientPosition[k]; + }*/ + //!!! 16032018 : next line as well as definition of B0Philips may need to be set to zero if Philips omits DiffusionBValue tag for B=0 + B0Philips = -1.0; //Philips may skip reporting B-values for B=0 volumes, so zero these + vRLPhilips = 0.0; + vAPPhilips = 0.0; + vFHPhilips = 0.0; + MRImageGradientOrientationNumber = 0; + break; } break; - case kPhotometricInterpretation: + case kPhotometricInterpretation: { char interp[kDICOMStr]; dcmStr (lLength, &buffer[lPos], interp); if (strcmp(interp, "PALETTE_COLOR") == 0) printError("Photometric Interpretation 'PALETTE COLOR' not supported\n"); - - break; + break; } case kPlanarRGB: d.isPlanarRGB = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); break; case kDim3: d.xyzDim[3] = dcmStrInt(lLength, &buffer[lPos]); + numberOfFrames = d.xyzDim[3]; break; case kSamplesPerPixel: d.samplesPerPixel = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); @@ -3622,6 +4072,12 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD case kLocationsInAcquisitionGE: locationsInAcquisitionGE = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); break; + case kProtocolDataBlockGE : + if (d.manufacturer != kMANUFACTURER_GE) break; + d.protocolBlockLengthGE = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); + d.protocolBlockStartGE = lPos+lFileOffset+4; + //printError("ProtocolDataBlockGE %d @ %d\n", d.protocolBlockLengthGE, d.protocolBlockStartGE); + break; case kDoseCalibrationFactor : d.doseCalibrationFactor = dcmStrFloat(lLength, &buffer[lPos]); break; @@ -3638,8 +4094,15 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD d.TR = dcmStrFloat(lLength, &buffer[lPos]); break; case kTE : - d.TE = dcmStrFloat(lLength, &buffer[lPos]); - break; + TE = dcmStrFloat(lLength, &buffer[lPos]); + if (d.TE <= 0.0) + d.TE = TE; + break; + case kEffectiveTE : { + TE = dcmFloatDouble(lLength, &buffer[lPos], d.isLittleEndian); + if (d.TE <= 0.0) + d.TE = TE; + break; } case kTI : d.TI = dcmStrFloat(lLength, &buffer[lPos]); break; @@ -3662,7 +4125,15 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD d.phaseFieldofView = dcmStrFloat(lLength, &buffer[lPos]); break; case kPixelBandwidth : + /*if (d.manufacturer == kMANUFACTURER_PHILIPS) { + //Private SQs can report different (more precise?) pixel bandwidth values than in the main header! + // https://github.com/rordenlab/dcm2niix/issues/170 + if (is2005140FSQ) break; + if ((lFileOffset + lPos) < sqEndPrivate) break; //inside private SQ, SQ has defined length + if (sqDepthPrivate > 0) break; //inside private SQ, SQ has undefined length + }*/ d.pixelBandwidth = dcmStrFloat(lLength, &buffer[lPos]); + //printWarning(" PixelBandwidth (0018,0095)====> %g @%d\n", d.pixelBandwidth, lPos); break; case kAcquisitionMatrix : if (lLength == 8) { @@ -3703,7 +4174,6 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD if ((lLength == 4) && (d.manufacturer == kMANUFACTURER_PHILIPS)) d.intenScalePhilips = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); break; - case kIntercept : d.intenIntercept = dcmStrFloat(lLength, &buffer[lPos]); break; @@ -3752,7 +4222,9 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD break; }*/ case kNumberOfDynamicScans: - d.numberOfDynamicScans = dcmStrInt(lLength, &buffer[lPos]); + //~d.numberOfDynamicScans = dcmStrInt(lLength, &buffer[lPos]); + numberOfDynamicScans = dcmStrInt(lLength, &buffer[lPos]); + break; 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'); @@ -3811,6 +4283,10 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD else d.sliceOrient = kSliceOrientTra; //transverse (axial) break; } + case kDiffusionBFactor : + if (d.manufacturer != kMANUFACTURER_PHILIPS) break; + B0Philips = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); + break; // case kDiffusionBFactor: // 2001,1003 // if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) { // d.CSA.numDti++; //increment with BFactor: on Philips slices with B=0 have B-factor but no diffusion directions @@ -3838,9 +4314,9 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD // value *unset*! // GE started using this tag in 27, and annoyingly, NOT including // the b value if it is 0 for the slice. - if((d.manufacturer != kMANUFACTURER_PHILIPS) || !is2005140FSQ){ - // d.CSA.numDti++; + //if((d.manufacturer != kMANUFACTURER_PHILIPS) || !is2005140FSQ){ + // d.CSA.numDti++; // if (d.CSA.numDti == 2) { //First time we know that this is a 4D DTI dataset // //d.dti4D = (TDTI *)malloc(kMaxDTI4D * sizeof(TDTI)); // dti4D->S[0].V[0] = d.CSA.dtiV[0]; @@ -3848,12 +4324,14 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD // dti4D->S[0].V[2] = d.CSA.dtiV[2]; // dti4D->S[0].V[3] = d.CSA.dtiV[3]; // } + B0Philips = dcmFloatDouble(lLength, &buffer[lPos], d.isLittleEndian); + //d.CSA.dtiV[0] = dcmFloatDouble(lLength, &buffer[lPos], d.isLittleEndian); set_bVal(&volDiffusion, dcmFloatDouble(lLength, &buffer[lPos], d.isLittleEndian)); // if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D)) // dti4D->S[d.CSA.numDti-1].V[0] = d.CSA.dtiV[0]; - } + //} break; case kDiffusionOrientation: // 0018, 9089 // Note that this is ahead of kPatientPosition (0020,0032), so @@ -3865,9 +4343,19 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD // if (((d.manufacturer == kMANUFACTURER_SIEMENS) || // ((d.manufacturer == kMANUFACTURER_PHILIPS) && !is2005140FSQ)) && // (isAtFirstPatientPosition || isnan(d.patientPosition[1]))) - if((d.manufacturer == kMANUFACTURER_SIEMENS) || - ((d.manufacturer == kMANUFACTURER_PHILIPS) && !is2005140FSQ)) + + //if((d.manufacturer == kMANUFACTURER_SIEMENS) || ((d.manufacturer == kMANUFACTURER_PHILIPS) && !is2005140FSQ)) + if((d.manufacturer == kMANUFACTURER_SIEMENS) || (d.manufacturer == kMANUFACTURER_PHILIPS)) { + float v[4]; + //dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, v); + //dcmMultiFloatDouble(lLength, &buffer[lPos], 3, v, d.isLittleEndian); + dcmMultiFloatDouble(lLength, &buffer[lPos], 3, v, d.isLittleEndian); + vRLPhilips = v[0]; + vAPPhilips = v[1]; + vFHPhilips = v[2]; + set_orientation0018_9089(&volDiffusion, lLength, &buffer[lPos], d.isLittleEndian); + } break; // case kSharedFunctionalGroupsSequence: // if ((d.manufacturer == kMANUFACTURER_SIEMENS) && isAtFirstPatientPosition) { @@ -3875,17 +4363,44 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD // // part of buffer[lPos]. // } // break; - case kSliceNumberMrPhilips : + + //case kSliceNumberMrPhilips : + // sliceNumberMrPhilips3D = dcmStrInt(lLength, &buffer[lPos]); + // break; + case kSliceNumberMrPhilips : { if (d.manufacturer != kMANUFACTURER_PHILIPS) break; - /*if ((d.patientPositionNumPhilips >= kMaxDTI4D) || (d.patientPositionNumPhilips >= locationsInAcquisitionPhilips)) { - d.patientPositionNumPhilips++; - break; + + sliceNumberMrPhilips = dcmStrInt(lLength, &buffer[lPos]); + int sliceNumber = sliceNumberMrPhilips; + /*//if (sliceNumberMrPhilips < 10) printMessage("%d ====> %d\n", sliceNumberMrPhilips, sliceNumber); + //printMessage("====> sliceNumberMrPhilips\t%d\t%d\t%d\n", locationsInAcquisitionPhilips, d.patientPositionNumPhilips, sliceNumber); + if ((sliceNumberMrPhilips >= locationsInAcquisitionPhilips) && (sliceNumberMrPhilips < (2*locationsInAcquisitionPhilips)) ) { + int p = sliceNumberMrPhilips - locationsInAcquisitionPhilips; + if (p < kMaxDTI4D) + dti4D->S[p].sliceNumberMrPhilipsVol2 = sliceNumber; + //~ if ((p > 0) && (abs(dti4D->S[p].sliceNumberMrPhilipsVol2 - dti4D->S[p -1].sliceNumberMrPhilipsVol2) > 1) ) + //~ d.isSlicesSpatiallySequentialPhilips = false; }*/ + //use public patientPosition if it exists - fall back to private patient position + if ((sliceNumber == 1) && (!isnan(patientPosition[1])) ) { + for (int k = 0; k < 4; k++) + patientPositionStartPhilips[k] = patientPosition[k]; + } else if ((sliceNumber == 1) && (!isnan(patientPositionPrivate[1]))) { + for (int k = 0; k < 4; k++) + patientPositionStartPhilips[k] = patientPositionPrivate[k]; + } + if ((sliceNumber == locationsInAcquisitionPhilips) && (!isnan(patientPosition[1])) ) { + for (int k = 0; k < 4; k++) + patientPositionEndPhilips[k] = patientPosition[k]; + } else if ((sliceNumber == locationsInAcquisitionPhilips) && (!isnan(patientPositionPrivate[1])) ) { + for (int k = 0; k < 4; k++) + patientPositionEndPhilips[k] = patientPositionPrivate[k]; + } + /*sliceNumberMrPhilips ++; if ((locationsInAcquisitionPhilips > 0) && (d.patientPositionNumPhilips == locationsInAcquisitionPhilips)) break; //we have acquired all slices in volume (e.g. all volumes after 1st for XYZT storage - int sliceNumber; - sliceNumber = dcmStrInt(lLength, &buffer[lPos]); + if ((d.patientPositionNumPhilips > 0) && (sliceNumber == dti4D->S[d.patientPositionNumPhilips-1].sliceNumberMrPhilips) ) break; //repeated spatial position (e.g. data saved XYTZ so several time points if (d.patientPositionNumPhilips >= kMaxDTI4D) { @@ -3893,25 +4408,32 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD break; } dti4D->S[d.patientPositionNumPhilips].sliceNumberMrPhilips = sliceNumber; - if ((d.patientPositionNumPhilips > 0) && (abs(dti4D->S[d.patientPositionNumPhilips].sliceNumberMrPhilips - dti4D->S[d.patientPositionNumPhilips -1].sliceNumberMrPhilips) > 1) ) - d.isSlicesSpatiallySequentialPhilips = false; //slices not sequential (1,2,3,4 or 4,3,2,1) but 4,3,1,2 + //~ if ((d.patientPositionNumPhilips > 0) && (abs(dti4D->S[d.patientPositionNumPhilips].sliceNumberMrPhilips - dti4D->S[d.patientPositionNumPhilips -1].sliceNumberMrPhilips) > 1) ) + //~ d.isSlicesSpatiallySequentialPhilips = false; //slices not sequential (1,2,3,4 or 4,3,2,1) but 4,3,1,2 d.patientPositionNumPhilips++; //Philips can save 3D acquisitions in a single file with slices stored in non-sequential order. We need to know the first and final spatial position - if (sliceNumber == 1) { - for (int k = 0; k < 4; k++) - patientPositionStartPhilips[k] = patientPosition[k]; - } - if (sliceNumber == locationsInAcquisitionPhilips) { - for (int k = 0; k < 4; k++) - patientPositionEndPhilips[k] = patientPosition[k]; - } + //printMessage("x====> - %d\n", sliceNumber); if (isVerbose > 1) printMessage("slice %d is spatial position %d\n", d.patientPositionNumPhilips, sliceNumber); - break; + */ + break; } case kNumberOfSlicesMrPhilips : if (d.manufacturer != kMANUFACTURER_PHILIPS) break; locationsInAcquisitionPhilips = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); + //printMessage("====> locationsInAcquisitionPhilips\t%d\n", locationsInAcquisitionPhilips); + break; + case kDiffusionDirectionRL: + if (d.manufacturer != kMANUFACTURER_PHILIPS) break; + vRLPhilips = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); + break; + case kDiffusionDirectionAP: + if (d.manufacturer != kMANUFACTURER_PHILIPS) break; + vAPPhilips = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); + break; + case kDiffusionDirectionFH: + if (d.manufacturer != kMANUFACTURER_PHILIPS) break; + vFHPhilips = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); break; // case kDiffusionDirectionRL: // if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) { @@ -3942,6 +4464,18 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD // d.CSA.dtiV[d.CSA.numDti-1][3] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);*/ // //http://www.na-mic.org/Wiki/index.php/NAMIC_Wiki:DTI:DICOM_for_DWI_and_DTI // break; + //~~ + case kPrivatePerFrameSq : + if (d.manufacturer != kMANUFACTURER_PHILIPS) break; + //if ((vr[0] == 'S') && (vr[1] == 'Q')) break; + //if (!is2005140FSQwarned) + // printWarning("expected VR of 2005,140F to be 'SQ' (prior DICOM->DICOM conversion error?)\n"); + is2005140FSQ = true; + //is2005140FSQwarned = true; + case kMRImageGradientOrientationNumber : + if (d.manufacturer == kMANUFACTURER_PHILIPS) + MRImageGradientOrientationNumber = dcmStrInt(lLength, &buffer[lPos]); + break; case kWaveformSq: d.imageStart = 1; //abort!!! printMessage("Skipping DICOM (audio not image) '%s'\n", fname); @@ -3988,10 +4522,25 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD case kOrientationACR : //use in emergency if kOrientation is not present! if (!isOrient) dcmMultiFloat(lLength, (char*)&buffer[lPos], 6, d.orient); break; - case kOrientation : + case kTemporalPositionIdentifier : + temporalPositionIdentifier = dcmStrInt(lLength, &buffer[lPos]); + break; + case kOrientation : { + if (isOrient) { //already read orient - read for this slice to see if it varies (localizer) + float orient[7]; + dcmMultiFloat(lLength, (char*)&buffer[lPos], 6, orient); + if ((!isSameFloatGE(d.orient[1], orient[1]) || !isSameFloatGE(d.orient[2], orient[2]) || !isSameFloatGE(d.orient[3], orient[3]) || + !isSameFloatGE(d.orient[4], orient[4]) || !isSameFloatGE(d.orient[5], orient[5]) || !isSameFloatGE(d.orient[6], orient[6]) ) ) { + if (!d.isLocalizer) + printMessage("slice orientation varies (localizer?) [%g %g %g %g %g %g] != [%g %g %g %g %g %g]\n", + d.orient[1], d.orient[2], d.orient[3],d.orient[4], d.orient[5], d.orient[6], + orient[1], orient[2], orient[3],orient[4], orient[5], orient[6]); + d.isLocalizer = true; + } + } dcmMultiFloat(lLength, (char*)&buffer[lPos], 6, d.orient); isOrient = true; - break; + break; } case kImagesInAcquisition : imagesInAcquisition = dcmStrInt(lLength, &buffer[lPos]); break; @@ -4033,7 +4582,7 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD isIconImageSequence = false; break; } //switch/case for groupElement - } //if nest + if (isVerbose > 1) { if ((lLength > 12) && (lLength < 128)) { //if length is greater than 8 bytes (+4 hdr) the data must be a string [or image data] char tagStr[kDICOMStr]; @@ -4048,9 +4597,9 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD tagStr[pos] = 'x'; } printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\t%s\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos, tagStr); - //printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tnest=%d\t%s\n", groupElement & 65535,groupElement>>16, lLength, lPos, nest, tagStr); } else - printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tnest=%d\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos, nest); + printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos); + //printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tsq=%d %d %d %c%c\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos, sqDepth, sqDepthPrivate, sqEndPrivate, 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); lPos = lPos + (lLength); @@ -4146,29 +4695,158 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD // else // printWarning("Spatial orientation ambiguous (tag 0020,0037 not found): %s\n", fname); // } - if ((d.numberOfDynamicScans < 2) && (!d.isSlicesSpatiallySequentialPhilips) && (!isnan(patientPositionStartPhilips[1])) && (!isnan(patientPositionEndPhilips[1]))) { - //to do: check for d.numberOfDynamicScans > 1 - for (int k = 0; k < 4; k++) { - d.patientPosition[k] = patientPositionStartPhilips[k]; - d.patientPositionLast[k] = patientPositionEndPhilips[k]; - } - printMessage("Slices not spatially contiguous: please check output [new feature]\n"); +/*if (d.isHasMagnitude) + printError("=====> mag %d %d\n", d.patientPositionRepeats, d.patientPositionSequentialRepeats ); +if (d.isHasPhase) + printError("=====> phase %d %d\n", d.patientPositionRepeats, d.patientPositionSequentialRepeats ); + + printError("=====> reps %d %d\n", d.patientPositionRepeats, d.patientPositionSequentialRepeats ); +*/ + /*if ((patientPositionSequentialRepeats > 1) && ( (d.xyzDim[3] % patientPositionSequentialRepeats) == 0 )) { + //will require Converting XYTZ to XYZT + //~ d.numberOfDynamicScans = d.xyzDim[3] / d.patientPositionSequentialRepeats; + //~ d.xyzDim[4] = d.xyzDim[3] / d.numberOfDynamicScans; + numberOfDynamicScans = d.xyzDim[3] / patientPositionSequentialRepeats; + d.xyzDim[4] = d.xyzDim[3] / numberOfDynamicScans; + + d.xyzDim[3] = d.numberOfDynamicScans; + }*/ + + if (numberOfFrames == 0) numberOfFrames = d.xyzDim[3]; + + if ((numberOfDynamicScans > 1) && (d.xyzDim[4] < 2) && (d.xyzDim[3] > 1) && ((d.xyzDim[3] % numberOfDynamicScans) == 0)) { + d.xyzDim[3] = d.xyzDim[3] / numberOfDynamicScans; + d.xyzDim[4] = numberOfDynamicScans; + } + if ((maxInStackPositionNumber > 1) && (d.xyzDim[4] < 2) && (d.xyzDim[3] > 1) && ((d.xyzDim[3] % maxInStackPositionNumber) == 0)) { + d.xyzDim[4] = d.xyzDim[3] / maxInStackPositionNumber; + d.xyzDim[3] = maxInStackPositionNumber; + } + + if ((!isnan(patientPositionStartPhilips[1])) && (!isnan(patientPositionEndPhilips[1]))) { + for (int k = 0; k < 4; k++) { + d.patientPosition[k] = patientPositionStartPhilips[k]; + d.patientPositionLast[k] = patientPositionEndPhilips[k]; + } + /*if (d.numberOfDynamicScans > 1) + printError("4D data with non-contiguous slices: please check output [newest feature]\n"); + else + printMessage("Slices not spatially contiguous: please check output [new feature]\n");*/ } + if (!isnan(patientPositionStartPhilips[1])) //for Philips data without + for (int k = 0; k < 4; k++) + d.patientPosition[k] = patientPositionStartPhilips[k]; + if (isVerbose) { - printMessage("%s\n patient position (0020,0032)\t%g\t%g\t%g\n",fname, d.patientPosition[1],d.patientPosition[2],d.patientPosition[3]); - printMessage("%s\n orient (0020,0037)\t%g\t%g\t%g\t%g\t%g\t%g\n",fname, d.orient[1],d.orient[2],d.orient[3], d.orient[4],d.orient[5],d.orient[6]); - printMessage(" acq %d img %d ser %ld dim %dx%dx%d mm %gx%gx%g offset %d dyn %d loc %d valid %d ph %d mag %d posReps %d nDTI %d 3d %d bits %d littleEndian %d echo %d coil %d TE %g TR %g\n",d.acquNum,d.imageNum,d.seriesNum,d.xyzDim[1],d.xyzDim[2],d.xyzDim[3],d.xyzMM[1],d.xyzMM[2],d.xyzMM[3],d.imageStart, d.numberOfDynamicScans, d.locationsInAcquisition, d.isValid, d.isHasPhase, d.isHasMagnitude,d.patientPositionSequentialRepeats, d.CSA.numDti, d.is3DAcq, d.bitsAllocated, d.isLittleEndian, d.echoNum, d.coilNum, d.TE, d.TR); - 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]); + printMessage("DICOM file %s:\n", fname); + printMessage(" patient position (0020,0032)\t%g\t%g\t%g\n", d.patientPosition[1],d.patientPosition[2],d.patientPosition[3]); + if (!isnan(patientPositionEndPhilips[1])) + printMessage(" patient position end (0020,0032)\t%g\t%g\t%g\n", patientPositionEndPhilips[1],patientPositionEndPhilips[2],patientPositionEndPhilips[3]); + printMessage(" orient (0020,0037)\t%g\t%g\t%g\t%g\t%g\t%g\n", d.orient[1],d.orient[2],d.orient[3], d.orient[4],d.orient[5],d.orient[6]); + printMessage(" acq %d img %d ser %ld dim %dx%dx%dx%d mm %gx%gx%g offset %d loc %d valid %d ph %d mag %d nDTI %d 3d %d bits %d littleEndian %d echo %d coil %d TE %g TR %g\n",d.acquNum,d.imageNum,d.seriesNum,d.xyzDim[1],d.xyzDim[2],d.xyzDim[3], d.xyzDim[4],d.xyzMM[1],d.xyzMM[2],d.xyzMM[3],d.imageStart, d.locationsInAcquisition, d.isValid, d.isHasPhase, d.isHasMagnitude, d.CSA.numDti, d.is3DAcq, d.bitsAllocated, d.isLittleEndian, d.echoNum, d.coilNum, d.TE, d.TR); + //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.patientPositionNumPhilips >= kMaxDTI4D) { - printError("Too many 2D slices in a single file [recompile with increased kMaxDTI4D] detected=%d, max = %d\n", d.patientPositionNumPhilips, kMaxDTI4D); - d.CSA.numDti = 0; + + if ((numDimensionIndexValues > 1) && (numDimensionIndexValues == numberOfFrames)) { + //Philips enhanced datasets can have custom slice orders and pack images with different TE, Phase/Magnitude/Etc. + if (isVerbose > 1) { // + int mn[MAX_NUMBER_OF_DIMENSIONS]; + int mx[MAX_NUMBER_OF_DIMENSIONS]; + for (int j = 0; j < MAX_NUMBER_OF_DIMENSIONS; j++) { + mx[j] = dcmDim[0].dimIdx[j]; + mn[j] = mx[j]; + for (int i = 0; i < numDimensionIndexValues; i++) { + if (mx[j] < dcmDim[i].dimIdx[j]) mx[j] = dcmDim[i].dimIdx[j]; + if (mn[j] > dcmDim[i].dimIdx[j]) mn[j] = dcmDim[i].dimIdx[j]; + } + } + printMessage(" DimensionIndexValues (0020,9157), dimensions with variability:\n"); + for (int i = 0; i < MAX_NUMBER_OF_DIMENSIONS; i++) + if (mn[i] != mx[i]) + printMessage(" Dimension %d Range: %d..%d\n", i, mn[i], mx[i]); + } //verbose > 1 + qsort(dcmDim, numberOfFrames, sizeof(struct TDCMdim), compareTDCMdim); + //for (int i = 0; i < numberOfFrames; i++) + // printf("%d -> %d %d %d %d\n", i, dcmDim[i].diskPos, dcmDim[i].dimIdx[1], dcmDim[i].dimIdx[2], dcmDim[i].dimIdx[3]); + + for (int i = 0; i < numberOfFrames; i++) + dti4D->sliceOrder[i] = dcmDim[i].diskPos; + if ((d.xyzDim[4] > 1) && (d.xyzDim[4] < kMaxDTI4D)) { //record variations in TE + d.isScaleOrTEVaries = false; + bool isTEvaries = false; + bool isScaleVaries = false; + //setting j = 1 in next few lines is a hack, just in case TE/scale/intercept listed AFTER dimensionIndexValues + int j = 0; + if (d.xyzDim[3] > 1) j = 1; + for (int i = 0; i < d.xyzDim[4]; i++) { + //dti4D->gradDynVol[i] = 0; //only PAR/REC + dti4D->TE[i] = dcmDim[j+(i * d.xyzDim[3])].TE; + dti4D->intenScale[i] = dcmDim[j+(i * d.xyzDim[3])].intenScale; + dti4D->intenIntercept[i] = dcmDim[j+(i * d.xyzDim[3])].intenIntercept; + dti4D->isPhase[i] = dcmDim[j+(i * d.xyzDim[3])].isPhase; + dti4D->intenScalePhilips[i] = dcmDim[j+(i * d.xyzDim[3])].intenScalePhilips; + if (dti4D->TE[i] != d.TE) isTEvaries = true; + if (dti4D->intenScale[i] != d.intenScale) isScaleVaries = true; + if (dti4D->intenIntercept[i] != d.intenIntercept) isScaleVaries = true; + if (dti4D->isPhase[i] != isPhase) d.isScaleOrTEVaries = true; + } + if((isScaleVaries) || (isTEvaries)) d.isScaleOrTEVaries = true; + if (isTEvaries) d.isMultiEcho = true; + //if echoVaries,count number of echoes + /*int echoNum = 1; + for (int i = 1; i < d.xyzDim[4]; i++) { + if (dti4D->TE[i-1] != dti4D->TE[i]) + }*/ + if ((isVerbose) && (d.isScaleOrTEVaries)) { + printMessage("Parameters vary across 3D volumes packed in single DICOM file:\n"); + for (int i = 0; i < d.xyzDim[4]; i++) + printMessage(" %d TE=%g Slope=%g Inter=%g PhilipsScale=%g Phase=%d\n", i, dti4D->TE[i], dti4D->intenScale[i], dti4D->intenIntercept[i], dti4D->intenScalePhilips[i], dti4D->isPhase[i] ); + } + } } + /* //Attempt to append ADC + printMessage("CXC grad %g %d %d\n", philDTI[0].V[0], maxGradNum, d.xyzDim[4]); + if ((maxGradNum > 1) && ((maxGradNum+1) == d.xyzDim[4]) ) { + //ADC map (non-zero b-value with zero vector length) + if (isVerbose) + printMessage("Final volume does not have an associated 0020,9157. Assuming final volume is an ADC map\n", philDTI[0].V[0], maxGradNum, d.xyzDim[4]); + philDTI[maxGradNum].V[0] = 1000.0; + philDTI[maxGradNum].V[1] = 0.0; + philDTI[maxGradNum].V[2] = 0.0; + philDTI[maxGradNum].V[3] = 0.0; + maxGradNum++; + }*/ + if ((minGradNum >= 1) && ((maxGradNum-minGradNum+1) == d.xyzDim[4])) { + //see ADNI DWI data for 018_S_4868 - the gradient numbers are in the range 2..37 for 36 volumes - no gradient number 1! + if (philDTI[minGradNum -1].V[0] >= 0) { + if (isVerbose) + printMessage("Using %d diffusion data directions coded by DimensionIndexValues\n", maxGradNum); + int off = 0; + if (minGradNum > 1) { + off = minGradNum - 1; + printWarning("DimensionIndexValues (0020,9157) is not indexed from 1 (range %d..%d). Please validate results\n", minGradNum, maxGradNum); + } + for (int i = 0; i < d.xyzDim[4]; i++) { + dti4D->S[i].V[0] = philDTI[i+off].V[0]; + dti4D->S[i].V[1] = philDTI[i+off].V[1]; + dti4D->S[i].V[2] = philDTI[i+off].V[2]; + dti4D->S[i].V[3] = philDTI[i+off].V[3]; + if (isVerbose > 1) + printMessage(" grad %d b=%g vec=%gx%gx%g\n", i, dti4D->S[i].V[0], dti4D->S[i].V[1], dti4D->S[i].V[2], dti4D->S[i].V[3]); + + } + d.CSA.numDti = maxGradNum - off; + } + } if (d.CSA.numDti >= kMaxDTI4D) { printError("Unable to convert DTI [recompile with increased kMaxDTI4D] detected=%d, max = %d\n", d.CSA.numDti, kMaxDTI4D); d.CSA.numDti = 0; } + if (d.imageNum == 0) { + printError("Image number (0020,0013) not found.\n"); + d.imageNum = 1; //not set + } if (multiBandFactor > d.CSA.multiBandFactor) d.CSA.multiBandFactor = multiBandFactor; //SMS reported in 0051,1011 but not CSA header //d.isValid = false; //debug only - will not create output! diff --git a/console/nii_dicom.h b/console/nii_dicom.h index fb513e0a..6fea7b85 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -38,10 +38,12 @@ extern "C" { #define kCCsuf " CompilerNA" //unknown compiler! #endif -#define kDCMvers "v1.0.20171215" kDCMsuf kCCsuf +#define kDCMvers "v1.0.20180325" kDCMsuf kCCsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic static const int kMaxDTI4D = 4096; //maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images +static const int kMaxSlice2D = 64000; //maximum number of 2D slices in 4D (Philips) images + #define kDICOMStr 64 #define kDICOMStrLarge 256 #define kMANUFACTURER_UNKNOWN 0 @@ -76,10 +78,14 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; struct TDTI { float V[4]; - int sliceNumberMrPhilips; + //int totalSlicesIn4DOrder; }; struct TDTI4D { 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 TE[kMaxDTI4D], intenScale[kMaxDTI4D], intenIntercept[kMaxDTI4D], intenScalePhilips[kMaxDTI4D]; + bool isPhase[kMaxDTI4D]; }; #ifdef _MSC_VER //Microsoft nomenclature for packed structures is different... @@ -110,14 +116,18 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; } TCSAitem; //Siemens csa item structure #endif struct TCSAdata { - bool isPhaseMap; - float sliceTiming[kMaxEPI3D], dtiV[4], sliceNormV[4], bandwidthPerPixelPhaseEncode, sliceMeasurementDuration; + float sliceTiming[kMaxEPI3D], dtiV[4], sliceNormV[4], bandwidthPerPixelPhaseEncode, sliceMeasurementDuration; int numDti, SeriesHeader_offset, SeriesHeader_length, multiBandFactor, sliceOrder, slice_start, slice_end, mosaicSlices,protocolSliceNumber1,phaseEncodingDirectionPositive; + bool isPhaseMap; + }; struct TDICOMdata { long seriesNum; int xyzDim[5]; - int modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, echoTrainLength, patientPositionNumPhilips, coilNum, echoNum, sliceOrient,numberOfDynamicScans, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,patientPositionSequentialRepeats,locationsInAcquisition, compressionScheme; + //numberOfDynamicScans, patientPositionNumPhilips + //patientPositionSequentialRepeats,patientPositionRepeats, + //maxGradDynVol, gradDynVol, + int protocolBlockStartGE, protocolBlockLengthGE, modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, echoTrainLength, coilNum, echoNum, sliceOrient, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,locationsInAcquisition, compressionScheme; float patientWeight, zSpacing, zThick, pixelBandwidth, SAR, phaseFieldofView, accelFactPE, flipAngle, fieldStrength, TE, TI, TR, intenScale, intenIntercept, intenScalePhilips, gantryTilt, lastScanLoc, angulation[4]; float orient[7], patientPosition[4], patientPositionLast[4], xyzMM[4], stackOffcentre[4]; float radionuclidePositronFraction, radionuclideTotalDose, radionuclideHalfLife, doseCalibrationFactor; //PET ISOTOPE MODULE ATTRIBUTES (C.8-57) @@ -128,25 +138,28 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; char imageComments[kDICOMStrLarge]; uint32_t dimensionIndexValues[MAX_NUMBER_OF_DIMENSIONS]; struct TCSAdata CSA; - bool isSegamiOasis, isDerived, isXRay, isMultiEcho, isSlicesSpatiallySequentialPhilips, isValid, is3DAcq, is2DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase,isHasMagnitude,isHasMixed, isFloat, isResampled; + //isSlicesSpatiallySequentialPhilips + bool isSegamiOasis, isScaleOrTEVaries, isDerived, isXRay, isMultiEcho, isValid, is3DAcq, is2DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase,isHasMagnitude,isHasMixed, isFloat, isResampled, isLocalizer; char phaseEncodingRC, patientSex; + //uint32_t *totalSlicesIn4DOrder; //Reordering array for Philips slices }; size_t nii_ImgBytes(struct nifti_1_header hdr); + int isSameFloatGE (float a, float b); struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, struct TDTI4D *dti4D); struct TDICOMdata readDICOM(char * fname); struct TDICOMdata clear_dicom_data(); - struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D *dti4D); + struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D *dti4D, bool isReadPhase); unsigned char * nii_flipY(unsigned char* bImg, struct nifti_1_header *h); unsigned char * nii_flipZ(unsigned char* bImg, struct nifti_1_header *h); - unsigned char * nii_reorderSlices(unsigned char* bImg, struct nifti_1_header *h, struct TDTI4D *dti4D); + //*unsigned char * nii_reorderSlices(unsigned char* bImg, struct nifti_1_header *h, struct TDTI4D *dti4D); void changeExt (char *file_name, const char* ext); unsigned char * nii_planar2rgb(unsigned char* bImg, struct nifti_1_header *hdr, int isPlanar); int isDICOMfile(const char * fname); //0=not DICOM, 1=DICOM, 2=NOTSURE(not part 10 compliant) void setQSForm(struct nifti_1_header *h, mat44 Q44i, bool isVerbose); int headerDcm2Nii2(struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1_header *h, int isVerbose); int headerDcm2Nii(struct TDICOMdata d, struct nifti_1_header *h, bool isComputeSForm) ; - unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct TDICOMdata dcm, bool iVaries, int compressFlag, int isVerbose); + unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct TDICOMdata dcm, bool iVaries, int compressFlag, int isVerbose, struct TDTI4D *dti4D); #ifdef __cplusplus } #endif diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 9e12cc0d..15144a0e 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -109,7 +109,6 @@ void dropTrailingFileSep(char *path) { // path[len] = '\0'; } - void getFileName( char *pathParent, const char *path) //if path is c:\d1\d2 then filename is 'd2' { const char *filename = strrchr(path, '/'); //UNIX @@ -312,6 +311,8 @@ void siemensPhilipsCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI printWarning("Saving %d DTI gradients. Validate vectors (matrix had a negative determinant).\n", d->CSA.numDti); //perhaps Siemens sagittal } else if ( d->sliceOrient == kSliceOrientTra) { printMessage("Saving %d DTI gradients. Validate vectors.\n", d->CSA.numDti); + } else if ( d->sliceOrient == kSliceOrientUnknown) { + printWarning("Saving %d DTI gradients. Validate vectors (image slice orientation not reported, e.g. 2001,100B).\n", d->CSA.numDti); } else { printWarning("Saving %d DTI gradients. Validate vectors (images are not axial slices).\n", d->CSA.numDti); } @@ -406,6 +407,28 @@ const void * memmem(const char *l, size_t l_len, const char *s, size_t s_len) { //n.b. memchr returns "const void *" not "void *" for Windows C++ https://msdn.microsoft.com/en-us/library/d7zdhf37.aspx #endif //for systems without memmem +/*int readKeyX(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)); + printWarning("<><><>\n"); + if (!keyPos) return 0; + printWarning("<><> %d\n", strlen(keyPos)); + int i = (int)strlen(key); + int numDigits = 0; + while( ( i< remLength) && (numDigits >= 0) ) { + printMessage("%c", keyPos[i]); + if( keyPos[i] >= '0' && keyPos[i] <= '9' ) { + ret = (10 * ret) + keyPos[i] - '0'; + numDigits ++; + } else if (numDigits > 0) + numDigits = -1; + i++; + } + printWarning("---> %d\n", ret); + return ret; +} //readKey() +*/ + 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)); @@ -608,6 +631,81 @@ void siemensCsaAscii(const char * filename, int csaOffset, int csaLength, float } // siemensCsaAscii() #endif //myReadAsciiCsa() +#ifndef myDisableMiniZ + #define myReadGeProtocolBlock +#endif +#ifdef myReadGeProtocolBlock +int geProtocolBlock(const char * filename, int geOffset, int geLength, int isVerbose, int* sliceOrder, int* viewOrder) { + *sliceOrder = 0; + *viewOrder = 0; + int ret = EXIT_FAILURE; + if ((geOffset < 0) || (geLength < 20)) return ret; + FILE * pFile = fopen ( filename, "rb" ); + if(pFile==NULL) return ret; + fseek (pFile , 0 , SEEK_END); + long lSize = ftell (pFile); + if (lSize < (geOffset+geLength)) { + fclose (pFile); + return ret; + } + fseek(pFile, geOffset, SEEK_SET); + mz_uint8 * pCmp = (mz_uint8*) malloc (geLength); + if(pCmp == NULL) return ret; + size_t result = fread (pCmp,1,geLength,pFile); + if ((int)result != geLength) return ret; + int cmpSz = geLength; + //http://www.forensicswiki.org/wiki/Gzip + // always little endia! http://www.onicos.com/staff/iz/formats/gzip.html + if (cmpSz < 20) return ret; + if ((pCmp[0] != 31) || (pCmp[1] != 139) || (pCmp[2] != 8)) return ret; //check signature and deflate algorithm + uint8_t flags = pCmp[3]; + bool isFNAME = ((flags & 0x08) == 0x08); + bool isFCOMMENT = ((flags & 0x10) == 0x10); + size_t hdrSz = 10; + if (isFNAME) {//skip null-terminated string FNAME + for (hdrSz = hdrSz; hdrSz < cmpSz; hdrSz++) + if (pCmp[hdrSz] == 0) break; + hdrSz++; + } + if (isFCOMMENT) {//skip null-terminated string COMMENT + for (hdrSz = hdrSz; hdrSz < cmpSz; hdrSz++) + if (pCmp[hdrSz] == 0) break; + hdrSz++; + } + size_t unCmpSz = ((size_t)pCmp[cmpSz-4])+((size_t)pCmp[cmpSz-3] << 8)+((size_t)pCmp[cmpSz-2] << 16)+((size_t)pCmp[cmpSz-1] << 24); + //printf(">> %d %d %zu %zu %zu\n", isFNAME, isFCOMMENT, cmpSz, unCmpSz, hdrSz); + + z_stream s; + memset (&s, 0, sizeof (z_stream)); + inflateInit2(&s, -MZ_DEFAULT_WINDOW_BITS); + mz_uint8 *pUnCmp = (mz_uint8 *)malloc((size_t)unCmpSz); + s.avail_out = unCmpSz; + s.next_in = pCmp+ hdrSz; + s.avail_in = cmpSz-hdrSz-8; + s.next_out = (uint8_t *) pUnCmp; + ret = mz_inflate(&s, MZ_SYNC_FLUSH); + if (ret != MZ_STREAM_END) { + free(pUnCmp); + return EXIT_FAILURE; + } + //https://groups.google.com/forum/#!msg/comp.protocols.dicom/mxnCkv8A-i4/W_uc6SxLwHQJ + // DISCOVERY MR750 / 24\MX\MR Software release:DV24.0_R01_1344.a) are now storing an XML file + // + if ((pUnCmp[0] == '<') && (pUnCmp[1] == '?')) + printWarning("New XML-based GE Protocol Block is not yet supported: please report issue on dcm2niix Github page\n"); + char keyStrSO[] = "SLICEORDER"; + *sliceOrder = readKey(keyStrSO, (char *) pUnCmp, unCmpSz); + char keyStrVO[] = "VIEWORDER"; //"MATRIXX"; + *viewOrder = readKey(keyStrVO, (char *) pUnCmp, unCmpSz); + 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); + printMessage("%s\n", pUnCmp); + } + free(pUnCmp); + return EXIT_SUCCESS; +} +#endif //myReadGeProtocolBlock() void json_Str(FILE *fp, const char *sLabel, char *sVal) { if (strlen(sVal) < 1) return; //fprintf(fp, sLabel, sVal ); @@ -784,12 +882,24 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, json_Float(fp, "\t\"SAR\": %g,\n", d.SAR ); if (d.echoNum > 1) fprintf(fp, "\t\"EchoNumber\": %d,\n", d.echoNum); if ((d.TE > 0.0) && (!d.isXRay)) fprintf(fp, "\t\"EchoTime\": %g,\n", d.TE / 1000.0 ); - json_Float(fp, "\t\"RepetitionTime\": %g,\n", d.TR / 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\"InversionTime\": %g,\n", d.TI / 1000.0 ); json_Float(fp, "\t\"FlipAngle\": %g,\n", d.flipAngle ); float pf = 1.0f; //partial fourier bool interp = false; //2D interpolation float phaseOversampling = 0.0; + int viewOrderGE = -1; + #ifdef myReadGeProtocolBlock + int sliceOrderGE = -1; + if ((d.manufacturer == kMANUFACTURER_GE) && (d.protocolBlockStartGE> 0) && (d.protocolBlockLengthGE > 19)) { + printWarning("Using GE Protocol Data Block for BIDS data (beware: new feature)\n"); + int ok = geProtocolBlock(filename, d.protocolBlockStartGE, d.protocolBlockLengthGE, opts.isVerbose, &sliceOrderGE, &viewOrderGE); + if (ok != EXIT_SUCCESS) + printWarning("Unable to decode GE protocol block\n"); + printMessage(" ViewOrder %d SliceOrder %d\n", viewOrderGE, sliceOrderGE); + } //read protocolBlockGE + #endif #ifdef myReadAsciiCsa if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.CSA.SeriesHeader_offset > 0) && (d.CSA.SeriesHeader_length > 0)) { int baseResolution, interpInt, partialFourier, echoSpacing, parallelReductionFactorInPlane; @@ -874,7 +984,6 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, bandwidthPerPixelPhaseEncode = d.CSA.bandwidthPerPixelPhaseEncode; json_Float(fp, "\t\"BandwidthPerPixelPhaseEncode\": %g,\n", bandwidthPerPixelPhaseEncode ); if (d.accelFactPE > 1.0) fprintf(fp, "\t\"ParallelReductionFactorInPlane\": %g,\n", d.accelFactPE); - //EffectiveEchoSpacing // Siemens bandwidthPerPixelPhaseEncode already accounts for the effects of parallel imaging, // interpolation, phaseOversampling, and phaseResolution, in the context of the size of the @@ -885,7 +994,6 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, if (d.effectiveEchoSpacingGE > 0.0) effectiveEchoSpacing = d.effectiveEchoSpacingGE / 1000000.0; json_Float(fp, "\t\"EffectiveEchoSpacing\": %g,\n", effectiveEchoSpacing); - // Calculate true echo spacing (should match what Siemens reports on the console) // i.e., should match "echoSpacing" extracted from the ASCII CSA header, when that exists double trueESfactor = 1.0; @@ -909,12 +1017,14 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, // https://github.com/rordenlab/dcm2niix/issues/130 if ((reconMatrixPE > 0) && (effectiveEchoSpacing > 0.0)) fprintf(fp, "\t\"TotalReadoutTime\": %g,\n", effectiveEchoSpacing * (reconMatrixPE - 1.0)); - json_Float(fp, "\t\"PixelBandwidth\": %g,\n", d.pixelBandwidth ); if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.dwellTime > 0)) fprintf(fp, "\t\"DwellTime\": %g,\n", d.dwellTime * 1E-9); // Phase encoding polarity - if (((d.phaseEncodingRC == 'R') || (d.phaseEncodingRC == 'C')) && (!d.is3DAcq) && ((d.CSA.phaseEncodingDirectionPositive == 1) || (d.CSA.phaseEncodingDirectionPositive == 0))) { + int phPos = d.CSA.phaseEncodingDirectionPositive; + if (viewOrderGE > -1) + phPos = viewOrderGE; + if (((d.phaseEncodingRC == 'R') || (d.phaseEncodingRC == 'C')) && (!d.is3DAcq) && (phPos >= 0)) { if (d.phaseEncodingRC == 'C') //Values should be "R"ow, "C"olumn or "?"Unknown fprintf(fp, "\t\"PhaseEncodingDirection\": \"j"); else if (d.phaseEncodingRC == 'R') @@ -925,17 +1035,40 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, //However, DICOM and NIfTI are reversed in the j (ROW) direction //Equivalent to dicm2nii's "if flp(iPhase), phPos = ~phPos; end" //for samples see https://github.com/rordenlab/dcm2niix/issues/125 - if (d.CSA.phaseEncodingDirectionPositive == -1) + if (phPos < 0) fprintf(fp, "?"); //unknown - else if ((d.CSA.phaseEncodingDirectionPositive == 0) && (d.phaseEncodingRC != 'C')) + else if ((phPos == 0) && (d.phaseEncodingRC != 'C')) fprintf(fp, "-"); - else if ((d.phaseEncodingRC == 'C') && (d.CSA.phaseEncodingDirectionPositive == 1) && (opts.isFlipY)) + else if ((d.phaseEncodingRC == 'C') && (phPos == 1) && (opts.isFlipY)) fprintf(fp, "-"); - else if ((d.phaseEncodingRC == 'C') && (d.CSA.phaseEncodingDirectionPositive == 0) && (!opts.isFlipY)) + else if ((d.phaseEncodingRC == 'C') && (phPos == 0) && (!opts.isFlipY)) fprintf(fp, "-"); fprintf(fp, "\",\n"); } //only save PhaseEncodingDirection if BOTH direction and POLARITY are known - // Slice Timing + // Slice Timing GE + if ((viewOrderGE > -1) && (h->dim[3] > 1) && (h->dim[4] > 1) && (d.TR > 0)) { // + //Warning: not correct for multiband sequences... not sure how these are stored + //Warning: will not create correct times for sparse acquisitions where DelayTimeInTR > 0 + float t = d.TR/ (float)h->dim[3] ; + fprintf(fp, "\t\"SliceTiming\": [\n"); + if (viewOrderGE == 1) {//interleaved ascending + for (int i = 0; i < h->dim[3]; i++) { + if (i != 0) + fprintf(fp, ",\n"); + int s = (i / 2); + if ((i % 2) != 0) s += (h->dim[3]+1)/2; + fprintf(fp, "\t\t%g", (float) s * t / 1000.0 ); + } + } else { //sequential ascending + for (int i = 0; i < h->dim[3]; i++) { + if (i != 0) + fprintf(fp, ",\n"); + fprintf(fp, "\t\t%g", (float) i * t / 1000.0 ); + } + } + fprintf(fp, "\t],\n"); + } + //Slice Timing Siemens if (d.CSA.sliceTiming[0] >= 0.0) { fprintf(fp, "\t\"SliceTiming\": [\n"); if (d.CSA.protocolSliceNumber1 > 1) { @@ -1134,35 +1267,35 @@ int * nii_SaveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],str bvals[i] = bvals[i] + (0.5 * i/numDti); //add a small bias so ties are kept in sequential order } 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 - // the isotropic trace or MD can be calculated) often come as - // b=0 and trace pairs, with the b=0 and trace in either order, - // and often as "ORIGINAL", even though the trace is not. - // The bval file is needed for downstream processing to know - // * which is the b=0 and which is the trace, and - // * what b is for the trace, - // so dcm2niix should *always* write the bval and bvec files, - // AND include the b for the trace for DWIs. - // One hackish way to accomplish that is to set *numADC = 0 - // when *numADC == 1 && numDti == 2. - // - Rob Reid, 2017-11-29. - if ((*numADC == 1) && ((numDti - *numADC) < 2)){ - *numADC = 0; - printMessage("Note: this appears to be a b=0+trace DWI; ADC/trace removal has been disabled.\n"); - } - else{ - if ((numDti - *numADC) < 2) { - 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", - *numADC); - } + // DWIs (i.e. short diffusion scans with too few directions to + // calculate tensors...they typically acquire b=0 + 3 b > 0 so + // the isotropic trace or MD can be calculated) often come as + // b=0 and trace pairs, with the b=0 and trace in either order, + // and often as "ORIGINAL", even though the trace is not. + // The bval file is needed for downstream processing to know + // * which is the b=0 and which is the trace, and + // * what b is for the trace, + // so dcm2niix should *always* write the bval and bvec files, + // AND include the b for the trace for DWIs. + // One hackish way to accomplish that is to set *numADC = 0 + // when *numADC == 1 && numDti == 2. + // - Rob Reid, 2017-11-29. + if ((*numADC == 1) && ((numDti - *numADC) < 2)){ + *numADC = 0; + printMessage("Note: this appears to be a b=0+trace DWI; ADC/trace removal has been disabled.\n"); + } + else{ + if ((numDti - *numADC) < 2) { + 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", + *numADC); + } } //sort ALL including ADC int * volOrderIndex = (int *) malloc(numDti * sizeof(int)); @@ -1203,7 +1336,6 @@ int * nii_SaveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],str if (!isSequential) printMessage("DTI volumes re-ordered by ascending b-value\n"); dcmList[indx0].CSA.numDti = numDti; //warning structure not changed outside scope! - geCorrectBvecs(&dcmList[indx0],sliceDir, vx); siemensPhilipsCorrectBvecs(&dcmList[indx0],sliceDir, vx); if (!opts.isFlipY ) { //!FLIP_Y&& (dcmList[indx0].CSA.mosaicSlices < 2) mosaics are always flipped in the Y direction @@ -1445,8 +1577,11 @@ int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopt } if (f == 'N') strcat (outname,dcm.patientName); - if (f == 'P') + if (f == 'P') { strcat (outname,dcm.protocolName); + if (strlen(dcm.protocolName) < 1) + printWarning("Unable to append protocol name (0018,1030) to filename (it is empty).\n"); + } if (f == 'Q') strcat (outname,dcm.scanningSequence); if (f == 'S') { @@ -1514,9 +1649,15 @@ int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopt sprintf(newstr, "_e%d", dcm.echoNum); strcat (outname,newstr); } - if (dcm.isHasPhase) - strcat (outname,"_ph"); //manufacturer name not available - + /*if (dcm.maxGradDynVol > 0) { //Philips segmented + sprintf(newstr, "_v%04d", dcm.gradDynVol+1); //+1 as indexed from zero + strcat (outname,newstr); + }*/ + if (dcm.isHasPhase) { + strcat (outname,"_ph"); //has phase map + if (dcm.isHasMagnitude) + strcat (outname,"Mag"); //Philips enhanced with BOTH phase and Magnitude in single file + } if (strlen(outname) < 1) strcpy(outname, "dcm2nii_invalidName"); if (outname[0] == '.') outname[0] = '_'; //make sure not a hidden file //eliminate illegal characters http://msdn.microsoft.com/en-us/library/windows/desktop/aa365247(v=vs.85).aspx @@ -2330,7 +2471,7 @@ void checkSliceTiming(struct TDICOMdata * d, struct TDICOMdata * d1) { printMessage("CSA slice timing based on 2nd volume, 1st volume corrupted (CMRR bug, range %g..%g, TR=%gms)\n", minT, maxT, d->TR); }//checkSliceTiming -int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmList[], struct TSearchList *nameList, struct TDCMopts opts, struct TDTI4D *dti4D) { +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 uint64_t indx = dcmSort[0].indx; @@ -2341,7 +2482,7 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis printMessage("Ignoring derived image(s) of series %ld %s\n", dcmList[indx].seriesNum, nameList->str[indx]); return EXIT_SUCCESS; } - if ((opts.isIgnoreDerivedAnd2D) && ((strcmp(dcmList[indx].sequenceName, "_fl3d1_ns")== 0) || (strcmp(dcmList[indx].sequenceName, "_fl2d1")== 0)) ) { + if ((opts.isIgnoreDerivedAnd2D) && ((dcmList[indx].isLocalizer) || (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]); return EXIT_SUCCESS; } @@ -2351,7 +2492,7 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis } bool saveAs3D = dcmList[indx].isHasPhase; struct nifti_1_header hdr0; - unsigned char * img = nii_loadImgXL(nameList->str[indx], &hdr0,dcmList[indx], iVaries, opts.compressFlag, opts.isVerbose); + unsigned char * img = nii_loadImgXL(nameList->str[indx], &hdr0,dcmList[indx], iVaries, opts.compressFlag, opts.isVerbose, dti4D); if (strlen(opts.imageComments) > 0) { for (int i = 0; i < 24; i++) hdr0.aux_file[i] = 0; //remove dcm.imageComments snprintf(hdr0.aux_file,24,"%s",opts.imageComments); @@ -2459,7 +2600,7 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis for (int i = 1; i < nConvert; i++) { //stack additional images indx = dcmSort[i].indx; //if (headerDcm2Nii(dcmList[indx], &hdrI) == EXIT_FAILURE) return EXIT_FAILURE; - img = nii_loadImgXL(nameList->str[indx], &hdrI, dcmList[indx],iVaries, opts.compressFlag, opts.isVerbose); + img = nii_loadImgXL(nameList->str[indx], &hdrI, dcmList[indx],iVaries, opts.compressFlag, opts.isVerbose, dti4D); if (img == NULL) return EXIT_FAILURE; if ((hdr0.dim[1] != hdrI.dim[1]) || (hdr0.dim[2] != hdrI.dim[2]) || (hdr0.bitpix != hdrI.bitpix)) { printError("Image dimensions differ %s %s",nameList->str[dcmSort[0].indx], nameList->str[indx]); @@ -2471,6 +2612,37 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis free(img); } } + if ((segVol >= 0) && (hdr0.dim[4] > 1)) { + int inVol = hdr0.dim[4]; + int nVol = 0; + for (int v = 0; v < inVol; v++) + if (dti4D->gradDynVol[v] == segVol) + nVol ++; + if (nVol < 1) { + printError("Series %d does not exist\n", segVol); + return EXIT_FAILURE; + } + size_t imgsz4D = imgsz; + if (nVol < 2) + hdr0.dim[0] = 3; //3D + hdr0.dim[4] = 1; + size_t imgsz3D = nii_ImgBytes(hdr0); + unsigned char *img4D = (unsigned char *)malloc(imgsz4D); + memcpy(&img4D[0], &imgM[0], imgsz4D); + free(imgM); + imgM = (unsigned char *)malloc(imgsz3D * nVol); + int outVol = 0; + for (int v = 0; v < inVol; v++) { + if ((dti4D->gradDynVol[v] == segVol) && (outVol < nVol)) { + memcpy(&imgM[outVol * imgsz3D], &img4D[v * imgsz3D], imgsz3D); + outVol ++; + } + } + hdr0.dim[4] = nVol; + imgsz = nii_ImgBytes(hdr0); + free(img4D); + saveAs3D = false; + } char pathoutname[2048] = {""}; if (nii_createFilename(dcmList[dcmSort[0].indx], pathoutname, opts) == EXIT_FAILURE) { free(imgM); @@ -2484,7 +2656,10 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis for(int i = 0; i < nConvert; ++i) dcmList[dcmSort[i].indx].converted2NII = 1; if (opts.numSeries < 0) { //report series number but do not convert - printMessage("\t%ld\t%s\n", dcmList[dcmSort[0].indx].seriesNum, pathoutname); + if (segVol >= 0) + 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(" %s\n",nameList->str[dcmSort[0].indx]); return EXIT_SUCCESS; } @@ -2505,8 +2680,12 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis // 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 for (; i < opts.numSeries; i++) { - if (opts.seriesNumber[i] == dcmList[dcmSort[0].indx].seriesNum) { + if (isSameFloatGE(opts.seriesNumber[i], seriesNum)) { + //if (opts.seriesNumber[i] == dcmList[dcmSort[0].indx].seriesNum) { break; } } @@ -2528,8 +2707,8 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis if ((hdr0.datatype == DT_UINT16) && (!dcmList[dcmSort[0].indx].isSigned)) nii_check16bitUnsigned(imgM, &hdr0); 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]); PhilipsPrecise(&dcmList[dcmSort[0].indx], opts.isPhilipsFloatNotDisplayScaling, &hdr0); - if (!dcmList[dcmSort[0].indx].isSlicesSpatiallySequentialPhilips) - nii_reorderSlices(imgM, &hdr0, dti4D); + //~ if (!dcmList[dcmSort[0].indx].isSlicesSpatiallySequentialPhilips) + //~ nii_reorderSlices(imgM, &hdr0, dti4D); if (hdr0.dim[3] < 2) printWarning("Check that 2D images are not mirrored.\n"); #ifndef HAVE_R @@ -2600,17 +2779,106 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis #endif free(imgM); return returnCode;//EXIT_SUCCESS; +}// saveDcm2NiiCore() + +int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmList[], struct TSearchList *nameList, struct TDCMopts opts, struct TDTI4D *dti4D) { + //this wrapper does nothing if all the images share the same echo time and scale + // however, it segments images when these properties vary + uint64_t indx = dcmSort[0].indx; + if ((!dcmList[indx].isScaleOrTEVaries) || (dcmList[indx].xyzDim[4] < 2)) + return saveDcm2NiiCore(nConvert, dcmSort, dcmList, nameList, opts, dti4D, -1); + if ((dcmList[indx].xyzDim[4]) && (dti4D->sliceOrder[0] < 0)) { + printError("Unexpected error for image with varying echo time or intensity scaling\n"); + return EXIT_FAILURE; + } + int ret = EXIT_SUCCESS; + //check for repeated echoes - count unique number of echoes + //code below checks for multi-echoes - not required if maxNumberOfEchoes reported in PARREC + int echoNum[kMaxDTI4D]; + int echo = 1; + for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) + echoNum[i] = 0; + echoNum[0] = 1; + for (int i = 1; i < dcmList[indx].xyzDim[4]; i++) { + for (int j = 0; j < i; j++) + if (dti4D->TE[i] == dti4D->TE[j]) echoNum[i] = echoNum[j]; + if (echoNum[i] == 0) { + echo++; + echoNum[i] = echo; + } + } + if (echo > 1) dcmList[indx].isMultiEcho = true; + //check for repeated volumes + int series = 1; + for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) + dti4D->gradDynVol[i] = 0; + dti4D->gradDynVol[0] = 1; + for (int i = 1; i < dcmList[indx].xyzDim[4]; i++) { + for (int j = 0; j < i; j++) + if ((dti4D->intenIntercept[i] == dti4D->intenIntercept[j]) && (dti4D->intenScale[i] == dti4D->intenScale[j]) && (dti4D->isPhase[i] == dti4D->isPhase[j]) && (dti4D->TE[i] == dti4D->TE[j])) + dti4D->gradDynVol[i] = dti4D->gradDynVol[j]; + if (dti4D->gradDynVol[i] == 0) { + series++; + dti4D->gradDynVol[i] = series; + } + } + for (int s = 1; s <= series; s++) { + for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) + if (dti4D->gradDynVol[i] == s) { + //dti4D->gradDynVol[i] = s; + //nVol ++; + dcmList[indx].TE = dti4D->TE[i]; + dcmList[indx].intenScale = dti4D->intenScale[i]; + dcmList[indx].intenIntercept = dti4D->intenIntercept[i]; + dcmList[indx].isHasPhase = dti4D->isPhase[i]; + dcmList[indx].intenScalePhilips = dti4D->intenScalePhilips[i]; + dcmList[indx].isHasMagnitude = false; + dcmList[indx].echoNum = echoNum[i]; + break; + } + int ret2 = saveDcm2NiiCore(nConvert, dcmSort, dcmList, nameList, opts, dti4D, s); + if (ret2 != EXIT_SUCCESS) ret = ret2; //return EXIT_SUCCESS only if ALL are successful + } + return ret; }// saveDcm2Nii() void fillTDCMsort(struct TDCMsort& tdcmref, const uint64_t indx, const struct TDICOMdata& dcmdata){ // Copy the relevant parts of dcmdata to tdcmref. tdcmref.indx = indx; + //printf("series/image %d %d\n", dcmdata.seriesNum, dcmdata.imageNum); tdcmref.img = ((uint64_t)dcmdata.seriesNum << 32) + dcmdata.imageNum; for(int i = 0; i < MAX_NUMBER_OF_DIMENSIONS; ++i) tdcmref.dimensionIndexValues[i] = dcmdata.dimensionIndexValues[i]; } // fillTDCMsort() int compareTDCMsort(void const *item1, void const *item2) { + //for quicksort http://blog.ablepear.com/2011/11/objective-c-tuesdays-sorting-arrays.html + struct TDCMsort const *dcm1 = (const struct TDCMsort *)item1; + struct TDCMsort const *dcm2 = (const struct TDCMsort *)item2; + + int retval = 0; // tie + + if (dcm1->img < dcm2->img) + retval = -1; + else if (dcm1->img > dcm2->img) + retval = 1; + if(retval != 0) return retval; //sorted images + // Check the dimensionIndexValues (useful for enhanced DICOM 4D series). + // ->img is basically behaving as a (seriesNum, imageNum) sort key + // concatenated into a (large) integer for qsort. That is unwieldy when + // dimensionIndexValues need to be compared, because the existence of + // uint128_t, uint256_t, etc. is not guaranteed. This sorts by + // (seriesNum, ImageNum, div[0], div[1], ...), or if you think of it as a + // number, the dimensionIndexValues come after the decimal point. + for(int i=0; i < MAX_NUMBER_OF_DIMENSIONS; ++i){ + if(dcm1->dimensionIndexValues[i] < dcm2->dimensionIndexValues[i]) + return -1; + else if(dcm1->dimensionIndexValues[i] > dcm2->dimensionIndexValues[i]) + return 1; + } + return retval; +} //compareTDCMsort() +/*int compareTDCMsort(void const *item1, void const *item2) { //for quicksort http://blog.ablepear.com/2011/11/objective-c-tuesdays-sorting-arrays.html struct TDCMsort const *dcm1 = (const struct TDCMsort *)item1; struct TDCMsort const *dcm2 = (const struct TDCMsort *)item2; @@ -2642,13 +2910,13 @@ int compareTDCMsort(void const *item1, void const *item2) { } } return retval; -} //compareTDCMsort() +} //compareTDCMsort()*/ -int isSameFloatGE (float a, float b) { +/*int isSameFloatGE (float a, float b) { //Kludge for bug in 0002,0016="DIGITAL_JACKET", 0008,0070="GE MEDICAL SYSTEMS" DICOM data: Orient field (0020:0037) can vary 0.00604261 == 0.00604273 !!! //return (a == b); //niave approach does not have any tolerance for rounding errors return (fabs (a - b) <= 0.0001); -} +}*/ int isSameFloatDouble (double a, double b) { //Kludge for bug in 0002,0016="DIGITAL_JACKET", 0008,0070="GE MEDICAL SYSTEMS" DICOM data: Orient field (0020:0037) can vary 0.00604261 == 0.00604273 !!! @@ -2764,7 +3032,8 @@ int singleDICOM(struct TDCMopts* opts, char *fname) { dcmList[0].converted2NII = 1; dcmList[0] = readDICOMv(nameList.str[0], opts->isVerbose, opts->compressFlag, &dti4D); //ignore compile warning - memory only freed on first of 2 passes fillTDCMsort(dcmSort[0], 0, dcmList[0]); - return saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, &dti4D); + int ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, &dti4D); + return ret; }// singleDICOM() size_t fileBytes(const char * fname) { @@ -2883,7 +3152,7 @@ int convert_parRec(struct TDCMopts opts) { nameList.str[0] = (char *)malloc(strlen(opts.indir)+1); strcpy(nameList.str[0],opts.indir); TDTI4D dti4D; - dcmList[0] = nii_readParRec(nameList.str[0], opts.isVerbose, &dti4D); + dcmList[0] = nii_readParRec(nameList.str[0], opts.isVerbose, &dti4D, false); struct TDCMsort dcmSort[1]; dcmSort[0].indx = 0; if (dcmList[0].isValid) @@ -2989,10 +3258,10 @@ int nii_loadDir(struct TDCMopts* opts) { bool convertError = false; for (int i = 0; i < (int)nDcm; i++ ) { 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) &&((dcmList[i].patientPositionNumPhilips > 1) || (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]); + //~ if ((dcmList[i].isValid) &&((dcmList[i].totalSlicesIn4DOrder != NULL) ||(dcmList[i].patientPositionNumPhilips > 1) || (dcmList[i].CSA.numDti > 1))) { //4D dataset: dti4D arrays require huge amounts of RAM - write this immediately + 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]); dcmList[i].converted2NII = 1; int ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, &dti4D); if (ret == EXIT_SUCCESS) diff --git a/console/nii_dicom_batch.h b/console/nii_dicom_batch.h index ced863e7..9cf0a7dc 100644 --- a/console/nii_dicom_batch.h +++ b/console/nii_dicom_batch.h @@ -28,7 +28,8 @@ extern "C" { bool isSave3D,isGz, isFlipY, isCreateBIDS, isSortDTIbyBVal, isAnonymizeBIDS, isOnlyBIDS, isCreateText, isIgnoreDerivedAnd2D, isPhilipsFloatNotDisplayScaling, isTiltCorrect, isRGBplanar, isOnlySingleFile, isForceStackSameSeries, isCrop; int isVerbose, compressFlag, gzLevel; //support for compressed data 0=none, char filename[512], outdir[512], indir[512], pigzname[512], optsname[512], indirParent[512], imageComments[24]; - long seriesNumber[MAX_NUM_SERIES], numSeries; + float seriesNumber[MAX_NUM_SERIES]; + long numSeries; #ifdef HAVE_R bool isScanOnly; void *imageList;