From e0e15a9ce549c20ec15530f1b60d74de01b0576e Mon Sep 17 00:00:00 2001 From: Tashrif Billah Date: Fri, 6 Nov 2020 23:13:59 -0500 Subject: [PATCH 01/41] correct name and source for PNL conversion library --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index b9d2bb1e..1ebe4d1b 100644 --- a/README.md +++ b/README.md @@ -161,7 +161,7 @@ The following tools exploit dcm2niix - [Neuroinformatics Database (NiDB)](https://github.com/gbook/nidb) is designed to store, retrieve, analyze, and share neuroimaging data. It uses dcm2niix for image QA and handling some formats. - [NiftyPET](https://niftypet.readthedocs.io/en/latest/install.html) provides PET image reconstruction and analysis, and uses dcm2niix to handle DICOM images. - [nipype](https://github.com/nipy/nipype) can use dcm2niix to convert images. - - [PNL-nipype](https://github.com/pnlbwh/Dummy-PNL-nipype) is a Python script that can convert dcm2niix created NIfTI files to the popular NRRD format (including DWI gradient tables). Note, recent versions of dcm2niix can directly convert DICOM images to NRRD. + - [conversion](https://github.com/pnlbwh/conversion) is a Python library that can convert dcm2niix created NIfTI files to the popular NRRD format (including DWI gradient tables). Note, recent versions of dcm2niix can directly convert DICOM images to NRRD. - [pydcm2niix is a Python module for working with dcm2niix](https://github.com/jstutters/pydcm2niix). - [pyBIDSconv provides a graphical format for converting DICOM images to the BIDS format](https://github.com/DrMichaelLindner/pyBIDSconv). It includes clever default heuristics for identifying Siemens scans. - [qsm](https://github.com/CAIsr/qsm) Quantitative Susceptibility Mapping software. From cb9948702936b0982adf6a5cebaa98e6a5d42064 Mon Sep 17 00:00:00 2001 From: Jon Clayden Date: Tue, 17 Nov 2020 16:17:34 +0000 Subject: [PATCH 02/41] Remove redundant underscores before creating subdirectories --- console/nii_dicom_batch.cpp | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 0b0a6a0f..6e227938 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -2714,6 +2714,18 @@ int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopt appendChar[0] = kPathSeparator; if ((strlen(pth) > 0) && (pth[strlen(pth)-1] != kPathSeparator) && (outname[0] != kPathSeparator)) strcat (baseoutname,appendChar); + + //remove redundant underscores + int len = strlen(outname); + int outpos = 0; + for (int inpos = 0; inpos < len; inpos ++) { + if ((outpos > 0) && (outname[inpos] == '_') && (outname[outpos-1] == '_')) + continue; + outname[outpos] = outname[inpos]; + outpos++; + } + outname[outpos] = 0; + //Allow user to specify new folders, e.g. "-f dir/%p" or "-f %s/%p/%m" // These folders are created if they do not exist char *sep = strchr(outname, kPathSeparator); @@ -2741,16 +2753,6 @@ int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopt strcat (newdir,ch); } } - //remove redundant underscores - int len = strlen(outname); - int outpos = 0; - for (int inpos = 0; inpos < len; inpos ++) { - if ((outpos > 0) && (outname[inpos] == '_') && (outname[outpos-1] == '_')) - continue; - outname[outpos] = outname[inpos]; - outpos++; - } - outname[outpos] = 0; //printMessage("path='%s' name='%s'\n", pathoutname, outname); //make sure outname is unique strcat (baseoutname,outname); From c57c37109f9e315d337038aa07a6c468438a1106 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Tue, 17 Nov 2020 13:48:10 -0500 Subject: [PATCH 03/41] Only report b-value for GE data with 0043,1039 if EPI (0018,9018) (https://github.com/rordenlab/dcm2niix/issues/449) --- console/nii_dicom.cpp | 3 --- console/nii_dicom_batch.cpp | 12 +++++++++++- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 9a6d77a4..160285aa 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -6708,9 +6708,6 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); } //printMessage(" tag=%04x,%04x length=%u pos=%ld %c%c nest=%d\n", groupElement & 65535,groupElement>>16, lLength, lPos,vr[0], vr[1], nest); #endif lPos = lPos + (lLength); - //printMessage("%d\n",d.imageStart); - //printMessage(" DWI bxyz %g %g %g %g %d\n", d.CSA.dtiV[0], d.CSA.dtiV[1], d.CSA.dtiV[2], d.CSA.dtiV[3], d.CSA.numDti); - } //while d.imageStart == 0 free (buffer); if (d.bitsStored < 0) d.isValid = false; diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 0b0a6a0f..eb1d3ea9 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -190,6 +190,7 @@ void geCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI *vx, int isV //COL then if swap the x and y value and reverse the sign on the z value. //If the phase encoding is not COL, then just reverse the sign on the x value. if ((d->manufacturer != kMANUFACTURER_GE) && (d->manufacturer != kMANUFACTURER_CANON)) return; + if ((!d->isEPI) && (d->CSA.numDti == 1)) d->CSA.numDti = 0; //issue449 if (d->CSA.numDti < 1) return; if ((toupper(d->patientOrient[0])== 'H') && (toupper(d->patientOrient[1])== 'F') && (toupper(d->patientOrient[2])== 'S')) ; //participant was head first supine @@ -1635,7 +1636,7 @@ tse3d: T2*/ } //only save PhaseEncodingDirection if BOTH direction and POLARITY are known //Slice Timing UIH or GE >>>> //in theory, we should also report XA10 slice times here, but see series 24 of https://github.com/rordenlab/dcm2niix/issues/236 - if ((d.modality != kMODALITY_PT) && (!d.is3DAcq) && (d.CSA.sliceTiming[0] >= 0.0)) { + if ((d.modality != kMODALITY_CT) && (d.modality != kMODALITY_PT) && (!d.is3DAcq) && (d.CSA.sliceTiming[0] >= 0.0)) { fprintf(fp, "\t\"SliceTiming\": [\n"); for (int i = 0; i < h->dim[3]; i++) { if (i != 0) @@ -2011,6 +2012,10 @@ int * nii_saveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],str dcmList[indx0].CSA.numDti = numDti; //warning structure not changed outside scope! geCorrectBvecs(&dcmList[indx0],sliceDir, vx, opts.isVerbose); siemensPhilipsCorrectBvecs(&dcmList[indx0],sliceDir, vx, opts.isVerbose); + if (dcmList[indx0].CSA.numDti < 1) { //issue449 + free(vx); + return NULL; + } if (!opts.isFlipY ) { //!FLIP_Y&& (dcmList[indx0].CSA.mosaicSlices < 2) mosaics are always flipped in the Y direction for (int i = 0; i < (numDti); i++) { if (fabs(vx[i].V[2]) > FLT_EPSILON) @@ -5011,7 +5016,12 @@ int sliceTimingCore(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct dcmList[dcmSort[0].indx].CSA.protocolSliceNumber1 = -1; } sliceTimingGE(d0, filename, opts, hdr, dcmSort, dcmList); + //ensure slice times have variability reverseSliceTiming(d0, verbose, hdr->dim[3]); + bool allSame = true; + for (int i = 0; i < hdr->dim[3]; i++) + if (!isSameFloatGE(d0->CSA.sliceTiming[i], d0->CSA.sliceTiming[0])) allSame = false; + if (allSame) d0->CSA.sliceTiming[0] = - 1.0; return sliceDir; } //sliceTiming() From ede3f49d4d78584f6694d6f2441480b7fc0d0b50 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Wed, 18 Nov 2020 08:17:08 -0500 Subject: [PATCH 04/41] Fix CMake for Apple Silicon (note similar change required for CloudFlare zlib) --- console/CMakeLists.txt | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/console/CMakeLists.txt b/console/CMakeLists.txt index e28566c9..0445c6fa 100644 --- a/console/CMakeLists.txt +++ b/console/CMakeLists.txt @@ -53,10 +53,15 @@ endif() # Compiler dependent flags include (CheckCXXCompilerFlag) if(UNIX) - check_cxx_compiler_flag(-msse2 HAS_SSE2) - if(HAS_SSE2) - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -msse2 -mfpmath=sse") - endif() + check_c_compiler_flag(-march=armv8-a+crc ARM_CRC) + if(ARM_CRC) + # wrong answer for Apple Silicon: check_cxx_compiler_flag(-msse2 HAS_SSE2) + else() + check_cxx_compiler_flag(-msse2 HAS_SSE2) + if(HAS_SSE2) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -msse2 -mfpmath=sse") + endif() + endif() endif() set(PROGRAMS dcm2niix) From ffa8435b4c28367d3e0166d5cd6c0dd5233fd109 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Wed, 18 Nov 2020 13:49:30 -0500 Subject: [PATCH 05/41] Apple. M1. use C++ not. C --- console/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/console/CMakeLists.txt b/console/CMakeLists.txt index 0445c6fa..05a16fbc 100644 --- a/console/CMakeLists.txt +++ b/console/CMakeLists.txt @@ -53,7 +53,7 @@ endif() # Compiler dependent flags include (CheckCXXCompilerFlag) if(UNIX) - check_c_compiler_flag(-march=armv8-a+crc ARM_CRC) + check_cxx_compiler_flag(-march=armv8-a+crc ARM_CRC) if(ARM_CRC) # wrong answer for Apple Silicon: check_cxx_compiler_flag(-msse2 HAS_SSE2) else() From 70f5fdf83a76413211cb0e6f5cc2a65afb6a2e86 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Tue, 24 Nov 2020 15:43:10 -0500 Subject: [PATCH 06/41] Use 1st study time if multiple times provided. --- console/nii_dicom.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 160285aa..902812f8 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -5216,7 +5216,8 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); dcmStr(lLength, &buffer[lPos], seriesTimeTxt); break; case kStudyTime : - dcmStr(lLength, &buffer[lPos], d.studyTime); + if (strlen(d.studyTime) < 2) + dcmStr(lLength, &buffer[lPos], d.studyTime); break; case kPatientName : dcmStr(lLength, &buffer[lPos], d.patientName); From d6acc85fad9bae989dde373b4bd3c32fc4374f2d Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Fri, 27 Nov 2020 06:29:45 -0500 Subject: [PATCH 07/41] When -n is specified, only save BIDS for requested series (https://github.com/rordenlab/dcm2niix/issues/453) --- console/nii_dicom.h | 2 +- console/nii_dicom_batch.cpp | 41 +++++++++++++++++-------------------- 2 files changed, 20 insertions(+), 23 deletions(-) diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 3fa01c4a..1574821d 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,7 +50,7 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20201102" +#define kDCMdate "v1.0.20201127" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index dd78fd8f..bd93b0e1 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -5509,8 +5509,25 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc free(imgM); return EXIT_FAILURE; } - //issue377(dcmList[indx0], &hdr0); return EXIT_SUCCESS; - nii_SaveBIDSX(pathoutname, dcmList[dcmSort[0].indx], opts, &hdr0, nameList->str[dcmSort[0].indx], dti4D); + + // skip converting if user has specified one or more series, but has not specified this one + if (opts.numSeries > 0) { //issue453: moved to before saveBIDS + int i = 0; + double seriesNum = (double) dcmList[dcmSort[0].indx].seriesUidCrc; + int segVolEcho = segVol; + if ((dcmList[dcmSort[0].indx].echoNum > 1) && (segVolEcho <= 0)) + segVolEcho = dcmList[dcmSort[0].indx].echoNum+1; + if (segVolEcho > 0) + seriesNum = seriesNum + ((double) segVolEcho - 1.0) / 10.0; + for (; i < opts.numSeries; i++) { + if (isSameDouble(opts.seriesNumber[i], seriesNum)) + break; + } + if (i == opts.numSeries) + return EXIT_SUCCESS; + } + if (opts.numSeries >= 0) //issue453 + nii_SaveBIDSX(pathoutname, dcmList[dcmSort[0].indx], opts, &hdr0, nameList->str[dcmSort[0].indx], dti4D); if (opts.isOnlyBIDS) { //note we waste time loading every image, however this ensures hdr0 matches actual output free(imgM); @@ -5572,26 +5589,6 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc isFlipZ = true; imgM = nii_flipZ(imgM, &hdr0); sliceDir = abs(sliceDir); //change this, we have flipped the image so GE DTI bvecs no longer need to be flipped! - } - // skip converting if user has specified one or more series, but has not specified this one - if (opts.numSeries > 0) { - int i = 0; - //double seriesNum = (double) dcmList[dcmSort[0].indx].seriesNum; - double seriesNum = (double) dcmList[dcmSort[0].indx].seriesUidCrc; - int segVolEcho = segVol; - if ((dcmList[dcmSort[0].indx].echoNum > 1) && (segVolEcho <= 0)) - segVolEcho = dcmList[dcmSort[0].indx].echoNum+1; - if (segVolEcho > 0) - seriesNum = seriesNum + ((double) segVolEcho - 1.0) / 10.0; - for (; i < opts.numSeries; i++) { - if (isSameDouble(opts.seriesNumber[i], seriesNum)) { - //if (opts.seriesNumber[i] == dcmList[dcmSort[0].indx].seriesNum) { - break; - } - } - if (i == opts.numSeries) { - return EXIT_SUCCESS; - } } nii_saveText(pathoutname, dcmList[dcmSort[0].indx], opts, &hdr0, nameList->str[indx]); int numADC = 0; From 25e314a68b1b0ffd956360465075d2fc19887805 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Wed, 2 Dec 2020 07:26:05 -0500 Subject: [PATCH 08/41] Single file mode memory allocation (https://github.com/rordenlab/dcm2niix/issues/454) --- console/nii_dicom_batch.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index bd93b0e1..d1b6ef3d 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -4996,8 +4996,9 @@ int sliceTimingCore(struct TDCMsort *dcmSort,struct TDICOMdata *dcmList, struct //uint64_t indx0 = dcmSort[0].indx; //uint64_t indx1 = dcmSort[1].indx; struct TDICOMdata * d0 = &dcmList[dcmSort[0].indx]; - uint64_t indx1 = dcmSort[1].indx; - if (nConvert < 2) indx1 = dcmSort[0].indx; + uint64_t indx1 = dcmSort[0].indx; + if (nConvert > 1) //use 2nd volume as CMRR bug can create bogus slice timing in first volume + indx1 = dcmSort[1].indx; struct TDICOMdata * d1 = &dcmList[indx1]; //oldSliceTimingGE(dcmSort, dcmList, hdr, verbose, filename, nConvert); sliceTimingUIH(dcmSort, dcmList, hdr, verbose, filename, nConvert); @@ -6166,12 +6167,13 @@ int singleDICOM(struct TDCMopts* opts, char *fname) { nameList.str[nameList.numItems] = (char *)malloc(strlen(fname)+1); strcpy(nameList.str[nameList.numItems],fname); nameList.numItems++; - struct TDCMsort dcmSort[1]; + TDCMsort * dcmSort = (TDCMsort *)malloc(sizeof(TDCMsort)); 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]); int ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, &dti4D); freeNameList(nameList); + free(dcmSort); return ret; }// singleDICOM() From 8725ec5b51486282e2aef37f66fc28558352f033 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Wed, 2 Dec 2020 13:21:45 -0500 Subject: [PATCH 09/41] leak (https://github.com/rordenlab/dcm2niix/issues/454) --- console/nii_dicom_batch.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index d1b6ef3d..b4928481 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -6174,6 +6174,7 @@ int singleDICOM(struct TDCMopts* opts, char *fname) { int ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, &dti4D); freeNameList(nameList); free(dcmSort); + free(dcmList); return ret; }// singleDICOM() From 14a66ae069c27e8588d95872e169acad2db68183 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Fri, 4 Dec 2020 06:20:58 -0500 Subject: [PATCH 10/41] Clear RepetitionTimeExcitation (https://github.com/rordenlab/dcm2niix/issues/457) --- console/nii_dicom.cpp | 4 ---- console/nii_dicom_batch.cpp | 3 +++ 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 902812f8..e9c0743c 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -1653,8 +1653,6 @@ dti4D->volumeOnsetTime[0] = -1; dti4D->decayFactor[0] = -1; dti4D->frameDuration[0] = -1; dti4D->intenScale[0] = 0.0; -dti4D->repetitionTimeExcitation = 0.0; -dti4D->repetitionTimeInversion = 0.0; strcpy(d.protocolName, ""); //erase dummy with empty strcpy(d.seriesDescription, ""); //erase dummy with empty strcpy(d.sequenceName, ""); //erase dummy with empty @@ -4064,8 +4062,6 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru dti4D->decayFactor[0] = -1; dti4D->frameDuration[0] = -1; dti4D->intenScale[0] = 0.0; - dti4D->repetitionTimeExcitation = 0.0; - dti4D->repetitionTimeInversion = 0.0; struct TVolumeDiffusion volDiffusion = initTVolumeDiffusion(&d, dti4D); struct stat s; if( stat(fname,&s) == 0 ) { diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index b4928481..113c4873 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -1670,6 +1670,7 @@ dti4D->volumeOnsetTime[0] = -1; dti4D->decayFactor[0] = -1; dti4D->intenScale[0] = 0.0; dti4D->repetitionTimeExcitation = 0.0; +dti4D->repetitionTimeInversion = 0.0; nii_SaveBIDSX(pathoutname, d, opts, h, filename, dti4D); }// nii_SaveBIDSX() @@ -5099,6 +5100,8 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc uint64_t indx1 = indx0; if (nConvert > 1) indx1 = dcmSort[1].indx; uint64_t indxEnd = dcmSort[nConvert-1].indx; + dti4D->repetitionTimeInversion = 0.0; //only set for Siemens and GE 3D T1 "TR" + dti4D->repetitionTimeExcitation = 0.0; //only set for Philips 3D T1 "TR" #ifdef newTilt //see issue 254 if (( (nConvert > 1) || (dcmList[indx0].xyzDim[3] > 1)) && ((dcmList[indx0].modality == kMODALITY_CT) || (dcmList[indx0].isXRay) || (dcmList[indx0].gantryTilt > 0.0))) { //issue372: enhanced DICOMs can also have gantry tilt dcmList[indx0].gantryTilt = computeGantryTiltPrecise(dcmList[indx0], dcmList[indxEnd], opts.isVerbose); From 8acdc7bcac48272d878162e04a877b939b78ed45 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Mon, 7 Dec 2020 12:51:29 -0500 Subject: [PATCH 11/41] Convert DICOM series where bit allocated (0028,0100) varies across slices (https://github.com/rordenlab/dcm2niix/issues/458) --- console/nii_dicom.h | 2 +- console/nii_dicom_batch.cpp | 26 ++++++++++++++------------ 2 files changed, 15 insertions(+), 13 deletions(-) diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 1574821d..810877ba 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,7 +50,7 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20201127" +#define kDCMdate "v1.0.20201207" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 113c4873..ddeb2864 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -2308,15 +2308,16 @@ bool intensityScaleVaries(int nConvert, struct TDCMsort dcmSort[],struct TDICOMd //some Siemens PET scanners generate 16-bit images where slice has its own scaling factor. // since NIfTI provides a single scaling factor for each file, these images require special consideration if (nConvert < 2) return false; - bool iVaries = false; + int dt = dcmList[dcmSort[0].indx].bitsAllocated; float iScale = dcmList[dcmSort[0].indx].intenScale; float iInter = dcmList[dcmSort[0].indx].intenIntercept; for (int i = 1; i < nConvert; i++) { //stack additional images uint64_t indx = dcmSort[i].indx; - if (fabs (dcmList[indx].intenScale - iScale) > FLT_EPSILON) iVaries = true; - if (fabs (dcmList[indx].intenIntercept- iInter) > FLT_EPSILON) iVaries = true; + if (dcmList[indx].bitsAllocated != dt) return true; + if (fabs (dcmList[indx].intenScale - iScale) > FLT_EPSILON) return true; + if (fabs (dcmList[indx].intenIntercept- iInter) > FLT_EPSILON) return true; } - return iVaries; + return false; } //intensityScaleVaries() /*unsigned char * nii_bgr2rgb(unsigned char* bImg, struct nifti_1_header *hdr) { @@ -5094,7 +5095,7 @@ void loadOverlay(char* imgname, unsigned char * img, int offset, int x, int y, i 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 + float *sliceMMarray = NULL; //only used if slices are not equidistant uint64_t indx = dcmSort[0].indx; uint64_t indx0 = dcmSort[0].indx; uint64_t indx1 = indx0; @@ -5513,7 +5514,6 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc free(imgM); return EXIT_FAILURE; } - // skip converting if user has specified one or more series, but has not specified this one if (opts.numSeries > 0) { //issue453: moved to before saveBIDS int i = 0; @@ -5735,6 +5735,8 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc #endif imgM = removeADC(&hdr0, imgM, numADC); #ifndef USING_R + if (iVaries) + printMessage("Saving as 32-bit float (slope, intercept or bits allocated varies).\n"); if (opts.isSaveNRRD) returnCode = nii_saveNRRD(pathoutname, hdr0, imgM, opts, dcmList[dcmSort[0].indx], dti4D, dcmList[indx0].CSA.numDti); else if (opts.isSave3D) @@ -6000,13 +6002,13 @@ int isSameFloatDouble (double a, double b) { } struct TWarnings { //generate a warning only once per set - bool acqNumVaries, bitDepthVaries, dateTimeVaries, echoVaries, triggerVaries, phaseVaries, coilVaries, forceStackSeries, seriesUidVaries, nameVaries, nameEmpty, orientVaries; + bool acqNumVaries, dimensionVaries, dateTimeVaries, echoVaries, triggerVaries, phaseVaries, coilVaries, forceStackSeries, seriesUidVaries, nameVaries, nameEmpty, orientVaries; }; TWarnings setWarnings() { TWarnings r; r.acqNumVaries = false; - r.bitDepthVaries = false; + r.dimensionVaries = false; r.dateTimeVaries = false; r.phaseVaries = false; r.echoVaries = false; @@ -6062,10 +6064,10 @@ bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, struct TDCMopts* opt if ((isSameStudyInstanceUID) && (d1.isXA10A) && (d2.isXA10A)) isSameTime = true; //kludge XA10A 0008,0030 incorrect https://github.com/rordenlab/dcm2niix/issues/236 if ((!isSameStudyInstanceUID) && (!isSameTime)) return false; - if ((d1.bitsAllocated != d2.bitsAllocated) || (d1.xyzDim[1] != d2.xyzDim[1]) || (d1.xyzDim[2] != d2.xyzDim[2]) || (d1.xyzDim[3] != d2.xyzDim[3]) ) { - if (!warnings->bitDepthVaries) - printMessage("Slices not stacked: dimensions or bit-depth varies\n"); - warnings->bitDepthVaries = true; + if ( (d1.xyzDim[1] != d2.xyzDim[1]) || (d1.xyzDim[2] != d2.xyzDim[2]) || (d1.xyzDim[3] != d2.xyzDim[3]) ) { + if (!warnings->dimensionVaries) + printMessage("Slices not stacked: dimensions vary across slices\n"); + warnings->dimensionVaries = true; return false; } #ifndef myIgnoreStudyTime From f827a591336a2ee2ae5c4861b97095a60ed48d2a Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Mon, 7 Dec 2020 15:37:48 -0500 Subject: [PATCH 12/41] Ignore LocationsInAcquisition (0020,1002) if number of slices is not evenly divisible with this value (https://github.com/rordenlab/dcm2niix/issues/459) --- console/nii_dicom.cpp | 1 - console/nii_dicom_batch.cpp | 9 +++++++-- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index e9c0743c..51550369 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -6895,7 +6895,6 @@ if (d.isHasPhase) if (!isnan(patientPositionStartPhilips[1])) //for Philips data without for (int k = 0; k < 4; k++) d.patientPosition[k] = patientPositionStartPhilips[k]; - //printMessage("%d %g\n", d.imageNum, sliceLocation); //shame this tag is optional if (isVerbose) { 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]); diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index ddeb2864..31d502eb 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -5190,8 +5190,13 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc dcmList[indx0].locationsInAcquisition = dcmList[indx0].locationsInAcquisitionConflict; } } - if ((nAcq == 1 ) && (dcmList[indx0].locationsInAcquisition > 0)) nAcq = nConvert/dcmList[indx0].locationsInAcquisition; - if (nAcq < 2 ) { + if ((nConvert > 1) && (nAcq == 1 ) && (dcmList[indx0].locationsInAcquisition > 0) ){ + if ((nConvert % dcmList[indx0].locationsInAcquisition) == 0) + nAcq = nConvert / dcmList[indx0].locationsInAcquisition; + else + printMessage("DICOM images may be missing, expected %d spatial locations per volume, but found %d slices.\n", dcmList[indx0].locationsInAcquisition, nConvert); + } + if (nAcq < 2 ) { nAcq = 0; for (int i = 0; i < nConvert; i++) if (isSamePosition(dcmList[dcmSort[0].indx],dcmList[dcmSort[i].indx])) nAcq++; From 7f9b266f416b7575e9f6728dacba07f7a77cfd5a Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Wed, 9 Dec 2020 08:57:51 -0500 Subject: [PATCH 13/41] PAR/REC field map calibrated in Hz given name _fieldmaphz (https://github.com/rordenlab/dcm2niix/issues/455) --- console/nii_dicom.cpp | 23 +++++++++++++++++------ console/nii_dicom.h | 4 ++-- console/nii_dicom_batch.cpp | 11 ++++++++--- 3 files changed, 27 insertions(+), 11 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 51550369..ea295aaa 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -848,6 +848,7 @@ struct TDICOMdata clear_dicom_data() { d.overlayStart[i] = 0; d.isHasOverlay = false; d.isPrivateCreatorRemap = false; + d.isRealIsPhaseMapHz = false; d.numberOfImagesInGridUIH = 0; d.phaseEncodingRC = '?'; d.patientSex = '?'; @@ -2017,9 +2018,10 @@ int kbval = 33; //V3: 27 } if (cols[kImageType] == 0) d.isHasMagnitude = true; if (cols[kImageType] != 0) d.isHasPhase = true; - if ((isSameFloat(cols[kImageType],18)) && (!isTypeWarning)) { - printWarning("Field map in Hz will be saved as the 'real' image.\n"); - isTypeWarning = true; + if (isSameFloat(cols[kImageType],18)) { + //printWarning("Field map in Hz will be saved as the 'real' image.\n"); + //isTypeWarning = true; + d.isRealIsPhaseMapHz = true; } else if (((cols[kImageType] < 0.0) || (cols[kImageType] > 4.0)) && (!isTypeWarning)) { printError("Unknown type %g: not magnitude[0], real[1], imaginary[2] or phase[3].\n", cols[kImageType]); isTypeWarning = true; @@ -2107,6 +2109,10 @@ int kbval = 33; //V3: 27 bool isReal = (cols[kImageType] == 1); bool isImaginary = (cols[kImageType] == 2); bool isPhase = (cols[kImageType] == 3); + if (cols[kImageType] == 18) { + isReal = true; + d.isRealIsPhaseMapHz = true; + } if (cols[kImageType] == 4) { if (!isType4Warning) { printWarning("Unknown image type (4). Be aware the 'phase' image is of an unknown type.\n"); @@ -2114,8 +2120,13 @@ int kbval = 33; //V3: 27 } isPhase = true; //2019 } - if ((cols[kImageType] < 0.0) || (cols[kImageType] > 3.0)) + if ((cols[kImageType] != 18) && ((cols[kImageType] < 0.0) || (cols[kImageType] > 3.0))) { + if (!isType4Warning) { + printWarning("Unknown image type (%g). Be aware the 'phase' image is of an unknown type.\n", round(cols[kImageType])); + isType4Warning = true; + } isReal = true; //<- this is not correct, kludge for bug in ROGERS_20180526_WIP_B0_NS_8_1.PAR + } if (isReal) vol += num3DExpected; if (isImaginary) vol += (2*num3DExpected); if (isPhase) vol += (3*num3DExpected); @@ -5171,7 +5182,7 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); //issue 256: Philips files report real ComplexImageComponent but Magnitude ImageType https://github.com/rordenlab/dcm2niix/issues/256 isPhase = false; isReal = false; - isImaginary = false; + isImaginary = false; isMagnitude = false; //see Table C.8-85 http://dicom.nema.org/medical/Dicom/2017c/output/chtml/part03/sect_C.8.13.3.html if ((buffer[lPos]=='R') && (toupper(buffer[lPos+1]) == 'E')) @@ -5185,7 +5196,7 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); //not mutually exclusive: possible for Philips enhanced DICOM to store BOTH magnitude and phase in the same image if (isPhase) d.isHasPhase = true; if (isReal) d.isHasReal = true; - if (isImaginary) d.isHasImaginary = true; + if (isImaginary) d.isHasImaginary = true; if (isMagnitude) d.isHasMagnitude = true; break; case kAcquisitionContrast: diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 810877ba..8acd9a7a 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -139,7 +139,7 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; bool isReal[kMaxDTI4D]; bool isImaginary[kMaxDTI4D]; bool isPhase[kMaxDTI4D]; - float repetitionTimeExcitation, repetitionTimeInversion; + float repetitionTimeExcitation, repetitionTimeInversion; }; #ifdef _MSC_VER //Microsoft nomenclature for packed structures is different... @@ -192,7 +192,7 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; char institutionAddress[kDICOMStrLarge], imageComments[kDICOMStrLarge]; uint32_t dimensionIndexValues[MAX_NUMBER_OF_DIMENSIONS]; struct TCSAdata CSA; - bool isPrivateCreatorRemap, isHasOverlay, isEPI, isIR, isPartialFourier, isDiffusion, isVectorFromBMatrix, isRawDataStorage, isGrayscaleSoftcopyPresentationState, isStackableSeries, isCoilVaries, isNonParallelSlices, isSegamiOasis, isXA10A, isScaleOrTEVaries, isScaleVariesEnh, isDerived, isXRay, isMultiEcho, isValid, is3DAcq, is2DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase, isHasImaginary, isHasReal, isHasMagnitude,isHasMixed, isFloat, isResampled, isLocalizer; + bool isRealIsPhaseMapHz, isPrivateCreatorRemap, isHasOverlay, isEPI, isIR, isPartialFourier, isDiffusion, isVectorFromBMatrix, isRawDataStorage, isGrayscaleSoftcopyPresentationState, isStackableSeries, isCoilVaries, isNonParallelSlices, isSegamiOasis, isXA10A, isScaleOrTEVaries, isScaleVariesEnh, isDerived, isXRay, isMultiEcho, isValid, is3DAcq, is2DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase, isHasImaginary, isHasReal, isHasMagnitude,isHasMixed, isFloat, isResampled, isLocalizer; char phaseEncodingRC, patientSex; }; diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 31d502eb..4c593e78 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -1140,6 +1140,8 @@ tse3d: T2*/ fprintf(fp,"\", \"REAL"); if ((d.isHasImaginary) && ((d.manufacturer == kMANUFACTURER_GE) || (strstr(d.imageType, "_IMAGINARY_") == NULL)) ) fprintf(fp,"\", \"IMAGINARY"); + if ((d.isRealIsPhaseMapHz))// && ((d.manufacturer == kMANUFACTURER_GE) || (strstr(d.imageType, "_IMAGINARY_") == NULL)) ) + fprintf(fp,"\", \"FIELDMAPHZ"); fprintf(fp, "\"],\n"); } if (d.isDerived) //DICOM is derived image or non-spatial file (sounds, etc) @@ -2667,7 +2669,10 @@ int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopt if ((isAddNamePostFixes) && (dcm.isHasImaginary)) { strcat (outname,"_imaginary"); //has phase map } - if ((isAddNamePostFixes) && (dcm.isHasReal)) { + if ((isAddNamePostFixes) && (dcm.isHasReal) && (dcm.isRealIsPhaseMapHz)) { + strcat (outname,"_fieldmaphz"); //has field map + } + if ((isAddNamePostFixes) && (dcm.isHasReal) && (!dcm.isRealIsPhaseMapHz)) { strcat (outname,"_real"); //has phase map } if ((isAddNamePostFixes) && (dcm.isHasPhase)) { @@ -6090,9 +6095,9 @@ bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, struct TDCMopts* opt //} return true; //we will stack these images, even if they differ in the following attributes } - if ((d1.isHasImaginary != d2.isHasImaginary) || (d1.isHasPhase != d2.isHasPhase) || ((d1.isHasReal != d2.isHasReal))) { + if ((d1.isHasImaginary != d2.isHasImaginary) || (d1.isHasPhase != d2.isHasPhase) || (d1.isHasReal != d2.isHasReal)) { if (!warnings->phaseVaries) - printMessage("Slices not stacked: some are phase/real/imaginary maps, others are not. Instances %d %d\n", d1.imageNum, d2.imageNum); + printMessage("Slices not stacked: some are phase/real/imaginary/phase maps, others are not. Instances %d %d\n", d1.imageNum, d2.imageNum); warnings->phaseVaries = true; return false; } From 2b96ee9fa26a4045434675b421d41725024937c9 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Thu, 10 Dec 2020 11:11:32 -0500 Subject: [PATCH 14/41] DICOM field map calibrated in Hz given name _fieldmaphz (https://github.com/rordenlab/dcm2niix/issues/455) --- console/nii_dicom.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index ea295aaa..0b642fbb 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -5074,6 +5074,8 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); //d.isDerived = true; //this would have 'i- y' skip MoCo images isMoCo = true; } + if ((slen > 5) && strstr(d.imageType, "B0") && strstr(d.imageType, "MAP")) + d.isRealIsPhaseMapHz = true; if((slen > 5) && strstr(d.imageType, "_ADC_") ) d.isDerived = true; if((slen > 5) && strstr(d.imageType, "_TRACEW_") ) From aa5dfdab76326f51db963f562d091b3d2f9653fb Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Thu, 10 Dec 2020 13:21:18 -0500 Subject: [PATCH 15/41] Fix PAR/REC ADC detection (https://github.com/rordenlab/dcm2niix/issues/462) --- FILENAMING.md | 1 + console/nii_dicom.cpp | 3 ++- console/nii_dicom_batch.cpp | 2 +- 3 files changed, 4 insertions(+), 2 deletions(-) diff --git a/FILENAMING.md b/FILENAMING.md index 4b79224c..620f9030 100644 --- a/FILENAMING.md +++ b/FILENAMING.md @@ -49,6 +49,7 @@ In general dcm2niix creates images with 3D dimensions, or 4 dimensions when the Some post-fixes are specific to Philips DICOMs - _ADC Philips specific case. A DWI image where derived isotropic, ADC or trace volume was appended to the series. Since this image will disrupt subsequent processing, and because subsequent processing (dwidenoise, topup, eddy) will yield better derived images, dcm2niix will also create an additional image without this volume. Therefore, the _ADC file should typically be discarded. If you want dcm2niix to discard these useless derived images, use the ignore feature ('-i y'). + - _fieldmaphz unwrapped B0 field map (in Hz) generated by a Philips scanner. Suggests Image Type (0008,0008) includes the terms 'B0' and 'MAP'. - _Raw Philips XX_* DICOMs (Raw Data Storage). - _PS Philips PS_* DICOMs (Grayscale Softcopy Presentation State). diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 0b642fbb..d85bcf9f 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -1969,7 +1969,7 @@ int kbval = 33; //V3: 27 } //diskSlice ++; bool isADC = false; - if ((maxNumberOfGradientOrients >= 2) && (cols[kbval] > 50) && isSameFloat(0.0, cols[kv1]) && isSameFloat(0.0, cols[kv2]) && isSameFloat(0.0, cols[kv2]) ) { + if ((maxNumberOfGradientOrients >= 2) && (cols[kbval] > 50) && isSameFloat(0.0, cols[kv1]) && isSameFloat(0.0, cols[kv2]) && isSameFloat(0.0, cols[kv3]) ) { isADC = true; ADCwarning = true; } @@ -2447,6 +2447,7 @@ int kbval = 33; //V3: 27 d.imageStart = 0; if (d.CSA.numDti >= kMaxDTI4D) { printError("Unable to convert DTI [increase kMaxDTI4D] found %d directions\n", d.CSA.numDti); + printMessage(" slices*grad*bval*cardiac*echo*dynamic*mix*label = %d*%d*%d*%d*%d*%d*%d*%d\n", d.xyzDim[3], maxNumberOfGradientOrients,maxNumberOfDiffusionValues, maxNumberOfCardiacPhases, maxNumberOfEchoes, maxNumberOfDynamics, maxNumberOfMixes, maxNumberOfLabels); d.CSA.numDti = 0; }; //check if dimensions vary diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 4c593e78..3e577421 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -6119,7 +6119,7 @@ bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, struct TDCMopts* opt } if (d1.coilCrc != d2.coilCrc) { if (!warnings->coilVaries) - printMessage("Slices not stacked: coil varies\n"); + printMessage("Slices not stacked: coil varies '%s' vs '%s'\n", d1.coilName, d2.coilName); warnings->coilVaries = true; *isCoilVaries = true; return false; From 28f63229fc3eb061548c0974b794c3e502fced52 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Fri, 11 Dec 2020 12:26:39 -0500 Subject: [PATCH 16/41] Prevent dti4D overflow --- console/nii_dicom.cpp | 23 +++++++++++++---------- console/nii_dicom.h | 6 +++++- 2 files changed, 18 insertions(+), 11 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index d85bcf9f..fa313140 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -1973,7 +1973,6 @@ int kbval = 33; //V3: 27 isADC = true; ADCwarning = true; } - //printMessage(">>%d %d\n", (int)cols[kSlice], diskSlice); if (numSlice2D < 1) { d.xyzMM[1] = cols[kXmm]; d.xyzMM[2] = cols[kYmm]; @@ -2138,7 +2137,7 @@ int kbval = 33; //V3: 27 free (cols); return d; } - // dti4D->S[vol].V[0] = cols[kbval]; + // dti4D->S[vol].V[0] = cols[kbval]; //dti4D->gradDynVol[vol] = gradDynVol; dti4D->TE[vol] = cols[kTEcho]; if (isSameFloatGE(cols[kTEcho], 0)) @@ -2161,7 +2160,7 @@ int kbval = 33; //V3: 27 if ((vol+1) > d.CSA.numDti) d.CSA.numDti = vol+1; } - if (numSlice2D < kMaxSlice2D) {//issue 363: intensity can vary with each 2D slice of 4D volume + if (numSlice2D < kMaxDTI4D) {//issue 363: intensity can vary with each 2D slice of 4D volume //printf("%d %g %g\n", numSlice2D, cols[kRI], cols[kRS]); dti4D->intenIntercept[numSlice2D] = cols[kRI]; dti4D->intenScale[numSlice2D] = cols[kRS]; @@ -2188,6 +2187,14 @@ int kbval = 33; //V3: 27 printError("Invalid PAR format header (unable to detect version or slices) %s\n", parname); return d; } + if (numSlice2D > kMaxSlice2D) { //check again after reading, as top portion of header does not report image types or isotropics + printError("Increase kMaxSlice2D from %d to at least %d (or use dicm2nii).\n", kMaxSlice2D, numSlice2D); + return d; + } + if (numSlice2D > kMaxDTI4D) { + printError("Increase kMaxDTI4D from %d to at least %d (or use dicm2nii).\n", kMaxDTI4D, numSlice2D); + return d; + } d.manufacturer = kMANUFACTURER_PHILIPS; d.isValid = true; d.isSigned = true; @@ -2212,10 +2219,6 @@ int kbval = 33; //V3: 27 } if (d.CSA.numDti > 0) d.CSA.numDti = maxVol; //e.g. gradient 2 can skip B=0 but include isotropic //remove unused slices - this will happen if unless we have all 4 image types: real, imag, mag, phase - if (numSlice2D > kMaxSlice2D) { //check again after reading, as top portion of header does not report image types or isotropics - printError("Increase kMaxSlice2D from %d to at least %d (or use dicm2nii).\n", kMaxSlice2D, numSlice2D); - d.isValid = false; - } int slice = 0; for (int i = 0; i < kMaxSlice2D; i++) { if (dti4D->sliceOrder[i] > -1) { //this slice was populated @@ -2251,7 +2254,7 @@ int kbval = 33; //V3: 27 dti4D->intenScalePhilips[i] = tmp.intenScalePhilips[j]; } } - d.isScaleOrTEVaries = true; + 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; @@ -2278,14 +2281,14 @@ int kbval = 33; //V3: 27 printWarning("Reported TR=%gms, measured TR=%gms (prospect. motion corr.?)\n", d.TR, TRms); d.TR = TRms; } - if ((isTypeWarning) && ((numSlice2D % num2DExpected) != 0) && ((numSlice2D % d.xyzDim[3]) == 0) ) { + if ((isTypeWarning) && ((numSlice2D % num2DExpected) != 0) && ((numSlice2D % d.xyzDim[3]) == 0) ) { num2DExpected = numSlice2D; } if ( ((numSlice2D % num2DExpected) != 0) && ((numSlice2D % d.xyzDim[3]) == 0) ) { num2DExpected = d.xyzDim[3] * (int)(numSlice2D / d.xyzDim[3]); if (!ADCwarning) printWarning("More volumes than described in header (ADC or isotropic?)\n"); } - if ((numSlice2D % num2DExpected) != 0) { + if ((numSlice2D % num2DExpected) != 0) { printMessage("Found %d slices, but expected divisible by %d: slices*grad*bval*cardiac*echo*dynamic*mix*labels = %d*%d*%d*%d*%d*%d*%d*%d %s\n", numSlice2D, num2DExpected, d.xyzDim[3], maxNumberOfGradientOrients, maxNumberOfDiffusionValues, maxNumberOfCardiacPhases, maxNumberOfEchoes, maxNumberOfDynamics, maxNumberOfMixes,maxNumberOfLabels, parname); diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 8acd9a7a..2bc99dfc 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -54,8 +54,12 @@ extern "C" { #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic -static const int kMaxDTI4D = 18000; //maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images static const int kMaxSlice2D = 64000; //maximum number of 2D slices in 4D (Philips) images +#if defined(_WIN64) || defined(_WIN32) //do not exceed windows stack +static const int kMaxDTI4D = 18000; //maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images +#else +static const int kMaxDTI4D = 18000; //issue460: maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images +#endif #define kDICOMStr 66 //64 characters plus NULL https://github.com/rordenlab/dcm2niix/issues/268 #define kDICOMStrLarge 256 From 28b72c239f08bacd6bcf80c2e1cb9e48bc7d91d2 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Fri, 11 Dec 2020 13:52:18 -0500 Subject: [PATCH 17/41] Support huge PAR/REC files (https://github.com/rordenlab/dcm2niix/issues/460) --- console/nii_dicom.cpp | 6 ++--- console/nii_dicom.h | 6 +---- console/nii_dicom_batch.cpp | 49 +++++++++++++++++++++---------------- 3 files changed, 32 insertions(+), 29 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index fa313140..d7c9ce03 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -1647,7 +1647,7 @@ int isSameFloatGE (float a, float b) { return (fabs (a - b) <= 0.0001); } -struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D *dti4D, bool isReadPhase) { +struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D *dti4D, bool isReadPhase) { struct TDICOMdata d = clear_dicom_data(); dti4D->sliceOrder[0] = -1; dti4D->volumeOnsetTime[0] = -1; @@ -1768,7 +1768,7 @@ int kbval = 33; //V3: 27 * maxNumberOfCardiacPhases * maxNumberOfEchoes * maxNumberOfDynamics * maxNumberOfMixes; num2DExpected = d.xyzDim[3] * num3DExpected; if ((num2DExpected ) >= kMaxSlice2D) { - printError("Use dicm2nii or increase kMaxDTI4D to be more than %d\n", num2DExpected); + printError("Use dicm2nii or increase kMaxSlice2D to be more than %d\n", num2DExpected); printMessage(" slices*grad*bval*cardiac*echo*dynamic*mix*label = %d*%d*%d*%d*%d*%d*%d*%d\n", d.xyzDim[3], maxNumberOfGradientOrients,maxNumberOfDiffusionValues, maxNumberOfCardiacPhases, maxNumberOfEchoes, maxNumberOfDynamics, maxNumberOfMixes, maxNumberOfLabels); @@ -2191,7 +2191,7 @@ int kbval = 33; //V3: 27 printError("Increase kMaxSlice2D from %d to at least %d (or use dicm2nii).\n", kMaxSlice2D, numSlice2D); return d; } - if (numSlice2D > kMaxDTI4D) { + if (numSlice2D > kMaxDTI4D) { //since issue460, kMaxSlice2D == kMaxSlice4D, so we should never get here printError("Increase kMaxDTI4D from %d to at least %d (or use dicm2nii).\n", kMaxDTI4D, numSlice2D); return d; } diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 2bc99dfc..76ca0383 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -55,11 +55,7 @@ extern "C" { static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic static const int kMaxSlice2D = 64000; //maximum number of 2D slices in 4D (Philips) images -#if defined(_WIN64) || defined(_WIN32) //do not exceed windows stack -static const int kMaxDTI4D = 18000; //maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images -#else -static const int kMaxDTI4D = 18000; //issue460: maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images -#endif +static const int kMaxDTI4D = kMaxSlice2D; //issue460: maximum number of DTI directions for 4D (Philips) images, also maximum number of 2D slices for Enhanced DICOM and PAR/REC #define kDICOMStr 66 //64 characters plus NULL https://github.com/rordenlab/dcm2niix/issues/268 #define kDICOMStrLarge 256 diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 3e577421..2ab943ce 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -3467,8 +3467,10 @@ void swapEndian(struct nifti_1_header* hdr, unsigned char* im, bool isNative) { int nii_saveNII(char * niiFilename, struct nifti_1_header hdr, unsigned char* im, struct TDCMopts opts, struct TDICOMdata d) { if (opts.isOnlyBIDS) return EXIT_SUCCESS; if (opts.isSaveNRRD) { - struct TDTI4D dti4D; - return nii_saveNRRD(niiFilename, hdr, im, opts, d, &dti4D, 0); + struct TDTI4D *dti4D = (struct TDTI4D *)malloc(sizeof(struct TDTI4D)); + int ret = nii_saveNRRD(niiFilename, hdr, im, opts, d, dti4D, 0); + free(dti4D); + return ret; } hdr.vox_offset = 352; size_t imgsz = nii_ImgBytes(hdr); @@ -5832,8 +5834,7 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis } } //bvec/bval saved for each series (real, phase, magnitude, imaginary) https://github.com/rordenlab/dcm2niix/issues/219 - TDTI4D dti4Ds; - dti4Ds = *dti4D; + TDTI4D dti4Ds = *dti4D; bool isHasDti = (dcmList[indx].CSA.numDti > 0); if ((isHasDti) && (dcmList[indx].CSA.numDti == dcmList[indx].xyzDim[4])) { int nDti = 0; @@ -6174,7 +6175,7 @@ int singleDICOM(struct TDCMopts* opts, char *fname) { return 0; } struct TDICOMdata *dcmList = (struct TDICOMdata *)malloc( sizeof(struct TDICOMdata)); - struct TDTI4D dti4D; + struct TDTI4D *dti4D = (struct TDTI4D *)malloc(sizeof(struct TDTI4D)); struct TSearchList nameList; nameList.maxItems = 1; // larger requires more memory, smaller more passes nameList.str = (char **) malloc((nameList.maxItems+1) * sizeof(char *)); //reserve one pointer (32 or 64 bits) per potential file @@ -6184,11 +6185,12 @@ int singleDICOM(struct TDCMopts* opts, char *fname) { nameList.numItems++; TDCMsort * dcmSort = (TDCMsort *)malloc(sizeof(TDCMsort)); 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 + 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]); - int ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, &dti4D); + int ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, dti4D); freeNameList(nameList); - free(dcmSort); + free(dti4D); + free(dcmSort); free(dcmList); return ret; }// singleDICOM() @@ -6229,7 +6231,7 @@ int textDICOM(struct TDCMopts* opts, char *fname) { #endif TDCMsort * dcmSort = (TDCMsort *)malloc(nConvert * sizeof(TDCMsort)); struct TDICOMdata *dcmList = (struct TDICOMdata *)malloc(nConvert * sizeof(struct TDICOMdata)); - struct TDTI4D dti4D; + struct TDTI4D *dti4D = (struct TDTI4D *)malloc(sizeof(struct TDTI4D)); struct TSearchList nameList; nameList.maxItems = nConvert; // larger requires more memory, smaller more passes nameList.str = (char **) malloc((nameList.maxItems+1) * sizeof(char *)); //reserve one pointer (32 or 64 bits) per potential file @@ -6243,15 +6245,16 @@ int textDICOM(struct TDCMopts* opts, char *fname) { nameList.str[nameList.numItems] = (char *)malloc(strlen(dcmname)+1); strcpy(nameList.str[nameList.numItems],dcmname); nameList.numItems++; - dcmList[nConvert] = readDICOMv(nameList.str[nConvert], opts->isVerbose, opts->compressFlag, &dti4D); //ignore compile warning - memory only freed on first of 2 passes + dcmList[nConvert] = readDICOMv(nameList.str[nConvert], opts->isVerbose, opts->compressFlag, dti4D); //ignore compile warning - memory only freed on first of 2 passes fillTDCMsort(dcmSort[nConvert], nConvert, dcmList[nConvert]); nConvert ++; } fclose(fp); qsort(dcmSort, nConvert, sizeof(struct TDCMsort), compareTDCMsort); //sort based on series and image numbers.... - int ret = saveDcm2Nii(nConvert, dcmSort, dcmList, &nameList, *opts, &dti4D); + int ret = saveDcm2Nii(nConvert, dcmSort, dcmList, &nameList, *opts, dti4D); free(dcmSort); free(dcmList); + free(dti4D); freeNameList(nameList); return ret; }//textDICOM() @@ -6364,12 +6367,13 @@ int convert_parRec(char * fnm, struct TDCMopts opts) { strcpy(nameList.str[0], fnm); //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, false); + struct TDTI4D *dti4D = (struct TDTI4D *)malloc(sizeof(struct TDTI4D)); + dcmList[0] = nii_readParRec(nameList.str[0], opts.isVerbose, dti4D, false); struct TDCMsort dcmSort[1]; dcmSort[0].indx = 0; if (dcmList[0].isValid) - ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, opts, &dti4D); + ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, opts, dti4D); + free(dti4D); free(dcmList);//if (nConvertTotal == 0) if (nameList.numItems < 1) printMessage("No valid PAR/REC files were found\n"); @@ -6655,7 +6659,7 @@ int nii_loadDirCore(char *indir, struct TDCMopts* opts) { progressPct = reportProgress(progressPct, kStage1Frac); //proportion correct, 0..100 // struct TDICOMdata dcmList [nameList.numItems]; //<- this exhausts the stack for large arrays struct TDICOMdata *dcmList = (struct TDICOMdata *)malloc(nameList.numItems * sizeof(struct TDICOMdata)); - struct TDTI4D dti4D; + struct TDTI4D *dti4D = (struct TDTI4D *)malloc(sizeof(struct TDTI4D)); int nConvertTotal = 0; bool compressionWarning = false; bool convertError = false; @@ -6674,13 +6678,13 @@ int nii_loadDirCore(char *indir, struct TDCMopts* opts) { convertError = true; continue; } - dcmList[i] = readDICOMv(nameList.str[i], opts->isVerbose, opts->compressFlag, &dti4D); //ignore compile warning - memory only freed on first of 2 passes + dcmList[i] = readDICOMv(nameList.str[i], opts->isVerbose, opts->compressFlag, dti4D); //ignore compile warning - memory only freed on first of 2 passes //if (!dcmList[i].isValid) printf(">>>>Not a valid DICOM %s\n", nameList.str[i]); - if ((dcmList[i].isValid) && ((dti4D.sliceOrder[0] >= 0) || (dcmList[i].CSA.numDti > 1))) { //4D dataset: dti4D arrays require huge amounts of RAM - write this immediately + 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); + int ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, dti4D); if (ret == EXIT_SUCCESS) nConvertTotal++; else @@ -6698,7 +6702,9 @@ int nii_loadDirCore(char *indir, struct TDCMopts* opts) { start = clock(); #endif if (opts->isRenameNotConvert) { - return EXIT_SUCCESS; + free(dcmList); + free(dti4D); + return EXIT_SUCCESS; } #ifdef USING_R if (opts->isScanOnly) { @@ -6783,7 +6789,7 @@ int nii_loadDirCore(char *indir, struct TDCMopts* opts) { else //nConvert = removeDuplicatesVerbose(nConvert, dcmSort, &nameList); nConvert = removeDuplicates(nConvert, dcmSort); - int ret = saveDcm2Nii(nConvert, dcmSort, dcmList, &nameList, *opts, &dti4D); + int ret = saveDcm2Nii(nConvert, dcmSort, dcmList, &nameList, *opts, dti4D); if (ret == EXIT_SUCCESS) nConvertTotal += nConvert; else @@ -6843,7 +6849,7 @@ int nii_loadDirCore(char *indir, struct TDCMopts* opts) { nConvert = removeDuplicatesVerbose(nConvert, dcmSort, &nameList); else nConvert = removeDuplicates(nConvert, dcmSort); - int ret = saveDcm2Nii(nConvert, dcmSort, dcmList, &nameList, *opts, &dti4D); + int ret = saveDcm2Nii(nConvert, dcmSort, dcmList, &nameList, *opts, dti4D); if (ret == EXIT_SUCCESS) nConvertTotal += nConvert; else @@ -6864,6 +6870,7 @@ int nii_loadDirCore(char *indir, struct TDCMopts* opts) { #endif if (opts->isProgress) progressPct = reportProgress(progressPct, 1); //proportion correct, 0..100 free(dcmList); + free(dti4D); freeNameList(nameList); if (convertError) { if (nConvertTotal == 0) From d7976467d6e4f8ce59fbd5accf68df4aed98fbfa Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Fri, 11 Dec 2020 16:22:38 -0500 Subject: [PATCH 18/41] Retain Philips Scaling Values for images where scaling does not vary but volumes must be separated (https://github.com/rordenlab/dcm2niix/issues/461) --- console/nii_dicom_batch.cpp | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 2ab943ce..e7bec953 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -5851,7 +5851,18 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis } //save each series bool isScaleVariesEnh = dcmList[indx].isScaleVariesEnh; //issue363: any variation in any image + float intenScale = dcmList[indx].intenScale; + float intenIntercept = dcmList[indx].intenIntercept; + float intenScalePhilips = dcmList[indx].intenScalePhilips; + float RWVIntercept = dcmList[indx].RWVIntercept; + float RWVScale = dcmList[indx].RWVScale; for (int s = 1; s <= series; s++) { + //issue461: assert these values as saveDcm2NiiCore modifies them when it applies Philips scaling + dcmList[indx].intenScale = intenScale; + dcmList[indx].intenIntercept = intenIntercept; + dcmList[indx].intenScalePhilips = intenScalePhilips; + dcmList[indx].RWVIntercept = RWVIntercept; + dcmList[indx].RWVScale = RWVScale; for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) { //for each volume if (dti4D->gradDynVol[i] == s) { //dti4D->gradDynVol[i] = s; @@ -5892,18 +5903,12 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis if (dti4Ds.intenScale[i] != dti4Ds.intenScale[0]) dcmList[indx].isScaleVariesEnh = true; if (dti4Ds.intenScalePhilips[i] != dti4Ds.intenScalePhilips[0]) dcmList[indx].isScaleVariesEnh = true; } - //in case scale doe not vary dcmList[indx].intenScale = dti4Ds.intenScale[0]; dcmList[indx].intenIntercept = dti4Ds.intenIntercept[0]; dcmList[indx].intenScalePhilips = dti4Ds.intenScalePhilips[0]; dcmList[indx].RWVIntercept = dti4Ds.RWVIntercept[0]; dcmList[indx].RWVScale = dti4Ds.RWVScale[0]; - } //if isScaleVariesEnh - - //if (dcmList[indx].isScaleVariesEnh) printf("Varies\n"); - //if (!dcmList[indx].isScaleVariesEnh) printf("no Varies\n"); - //printf("%g %g %g\n", dcmList[indx].intenScale, dcmList[indx].intenIntercept, dcmList[indx].intenScalePhilips); - //continue; + } if (s > 1) dcmList[indx].CSA.numDti = 0; //only save bvec for first type (magnitude) int ret2 = saveDcm2NiiCore(nConvert, dcmSort, dcmList, nameList, opts, &dti4Ds, s); if (ret2 != EXIT_SUCCESS) ret = ret2; //return EXIT_SUCCESS only if ALL are successful From 95b18f979d642f523c91ada2604fcad10770f434 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Fri, 11 Dec 2020 16:25:05 -0500 Subject: [PATCH 19/41] Bump version date --- console/nii_dicom.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 76ca0383..20fb26e4 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,7 +50,7 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20201207" +#define kDCMdate "v1.0.20201211" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic From 8b99a984c9f796b412c6874f26d529d3aad5ded7 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Wed, 16 Dec 2020 05:32:52 -0500 Subject: [PATCH 20/41] Increase details for series stacking, enhance merge override mode (https://github.com/rordenlab/dcm2niix/issues/463) --- console/nii_dicom.cpp | 9 +++++ console/nii_dicom.h | 4 +-- console/nii_dicom_batch.cpp | 66 +++++++++++++++++++++++++++++-------- 3 files changed, 63 insertions(+), 16 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index d7c9ce03..6b158fc7 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -5017,6 +5017,15 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); } else if ((compressFlag != kCompressNone) && (strcmp(transferSyntax, "1.2.840.10008.1.2.4.90") == 0)) { d.compressionScheme = kCompressYes; //printMessage("JPEG2000 Lossless support is new: please validate conversion\n"); + } else if ((strcmp(transferSyntax, "1.2.840.10008.1.2.1.99") == 0)){ + //n.b. Deflate compression applied applies to the encoding of the **entire** DICOM Data Set, not just image data + // see https://www.medicalconnections.co.uk/kb/Transfer-Syntax/ + //#ifndef myDisableZLib + //d.compressionScheme = kCompressDeflate; + //#else + printWarning("Unsupported transfer syntax '%s' (inflate files with 'dcmconv +te gz.dcm raw.dcm' or 'gdcmconv -w gz.dcm raw.dcm)'\n",transferSyntax); + d.imageStart = 1;//abort as invalid (imageStart MUST be >128) + //#endif } else if ((compressFlag != kCompressNone) && (strcmp(transferSyntax, "1.2.840.10008.1.2.4.91") == 0)) { d.compressionScheme = kCompressYes; //printMessage("JPEG2000 support is new: please validate conversion\n"); diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 20fb26e4..5360cdd2 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,11 +50,11 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20201211" +#define kDCMdate "v1.0.20201216" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic -static const int kMaxSlice2D = 64000; //maximum number of 2D slices in 4D (Philips) images +static const int kMaxSlice2D = 65535; //issue460 maximum number of 2D slices in 4D (Philips) images static const int kMaxDTI4D = kMaxSlice2D; //issue460: maximum number of DTI directions for 4D (Philips) images, also maximum number of 2D slices for Enhanced DICOM and PAR/REC #define kDICOMStr 66 //64 characters plus NULL https://github.com/rordenlab/dcm2niix/issues/268 diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index e7bec953..c3f6160d 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -6018,14 +6018,18 @@ int isSameFloatDouble (double a, double b) { } struct TWarnings { //generate a warning only once per set - bool acqNumVaries, dimensionVaries, dateTimeVaries, echoVaries, triggerVaries, phaseVaries, coilVaries, forceStackSeries, seriesUidVaries, nameVaries, nameEmpty, orientVaries; + bool manufacturerVaries, modalityVaries, derivedVaries, acqNumVaries, dimensionVaries, dateTimeVaries, studyUidVaries, echoVaries, triggerVaries, phaseVaries, coilVaries, forceStackSeries, seriesUidVaries, nameVaries, nameEmpty, orientVaries; }; TWarnings setWarnings() { TWarnings r; + r.manufacturerVaries = false; + r.modalityVaries = false; + r.derivedVaries = false; r.acqNumVaries = false; r.dimensionVaries = false; r.dateTimeVaries = false; + r.studyUidVaries = false; r.phaseVaries = false; r.echoVaries = false; r.triggerVaries = false; @@ -6041,14 +6045,29 @@ TWarnings setWarnings() { bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, struct TDCMopts* opts, struct TWarnings* warnings, bool *isMultiEcho, bool *isNonParallelSlices, bool *isCoilVaries) { //returns true if d1 and d2 should be stacked together as a single output if (!d1.isValid) return false; - if (!d2.isValid) return false; - if (d1.modality != d2.modality) return false; //do not stack MR and CT data! - if (d1.isDerived != d2.isDerived) return false; //do not stack raw and derived image types + if (!d2.isValid) return false; + if ((opts->isVerbose) && (d1.seriesNum == d2.seriesNum)) { + //one would never want to combine in these situations: only raise warning for verbose modes to help troubleshooting + if ((d1.manufacturer != d2.manufacturer) && (!warnings->manufacturerVaries)) { + printMessage("Volumes not stacked: manufacturer varies.\n"); + warnings->manufacturerVaries = true; + } + if ((d1.modality != d2.modality) && (!warnings->modalityVaries)) { + printMessage("Volumes not stacked: modality varies.\n"); + warnings->modalityVaries = true; + } + if ((d1.isDerived != d2.isDerived) && (!warnings->derivedVaries)) { + printMessage("Volumes not stacked: derived varies.\n"); + warnings->derivedVaries = true; + } + } if (d1.manufacturer != d2.manufacturer) return false; //do not stack data from different vendors - bool isForceStackSeries = false; + if (d1.modality != d2.modality) return false; //do not stack MR and CT data! + if (d1.isDerived != d2.isDerived) return false; //do not stack raw and derived image types + bool isForceStackSeries = false; if ((opts->isForceStackDCE) && (d1.isStackableSeries) && (d2.isStackableSeries) && (d1.seriesNum != d2.seriesNum)) { if (!warnings->forceStackSeries) - printMessage("Siemens volumes stacked despite varying series number (use '-m o' to turn off merging).\n"); + printMessage("Volumes stacked despite varying series number (use '-m o' to turn off merging).\n"); warnings->forceStackSeries = true; isForceStackSeries = true; } @@ -6079,8 +6098,20 @@ bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, struct TDCMopts* opt bool isSameTime = isSameFloatDouble(d1.dateTime, d2.dateTime); if ((isSameStudyInstanceUID) && (d1.isXA10A) && (d2.isXA10A)) isSameTime = true; //kludge XA10A 0008,0030 incorrect https://github.com/rordenlab/dcm2niix/issues/236 - if ((!isSameStudyInstanceUID) && (!isSameTime)) return false; - if ( (d1.xyzDim[1] != d2.xyzDim[1]) || (d1.xyzDim[2] != d2.xyzDim[2]) || (d1.xyzDim[3] != d2.xyzDim[3]) ) { + bool isDimensionVaries = ( (d1.xyzDim[1] != d2.xyzDim[1]) || (d1.xyzDim[2] != d2.xyzDim[2]) || (d1.xyzDim[3] != d2.xyzDim[3]) ); + if ((!isSameStudyInstanceUID) && (!isSameTime)) { + if (opts->isForceStackDCE) { + if (!warnings->studyUidVaries) + printMessage("Slices stacked despite Study Date/Time (0008,0020;0008,0030) and Study UID (0020,000E) variation %12.12f ~= %12.12f\n", d1.dateTime, d2.dateTime); + warnings->studyUidVaries = true; + } else { + if (!warnings->studyUidVaries) + printMessage("Slices not stacked: Study Date/Time (0008,0020;0008,0030) and Study UID (0020,000E) varies %12.12f ~= %12.12f\n", d1.dateTime, d2.dateTime); + warnings->studyUidVaries = true; + return false; + } + } + if (isDimensionVaries) { if (!warnings->dimensionVaries) printMessage("Slices not stacked: dimensions vary across slices\n"); warnings->dimensionVaries = true; @@ -6089,7 +6120,7 @@ bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, struct TDCMopts* opt #ifndef myIgnoreStudyTime if (!isSameTime) { //beware, some vendors incorrectly store Image Time (0008,0033) as Study Time (0008,0030). if (!warnings->dateTimeVaries) - printMessage("Slices not stacked: Study Date/Time (0008,0020 / 0008,0030) varies %12.12f ~= %12.12f\n", d1.dateTime, d2.dateTime); + printMessage("Slices not stacked: Study Date/Time (0008,0020;0008,0030) varies %12.12f ~= %12.12f\n", d1.dateTime, d2.dateTime); warnings->dateTimeVaries = true; return false; } @@ -6124,11 +6155,18 @@ bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, struct TDCMopts* opt return false; } if (d1.coilCrc != d2.coilCrc) { - if (!warnings->coilVaries) - printMessage("Slices not stacked: coil varies '%s' vs '%s'\n", d1.coilName, d2.coilName); - warnings->coilVaries = true; - *isCoilVaries = true; - return false; + if (opts->isForceStackDCE) { + if (!warnings->coilVaries) + printMessage("Slices stacked despite coil variation '%s' vs '%s'\n", d1.coilName, d2.coilName); + warnings->coilVaries = true; + *isCoilVaries = true; + } else { + if (!warnings->coilVaries) + printMessage("Slices not stacked: coil varies '%s' vs '%s'\n", d1.coilName, d2.coilName); + warnings->coilVaries = true; + *isCoilVaries = true; + return false; + } } if ((strlen(d1.protocolName) < 1) && (strlen(d2.protocolName) < 1)) { if (!warnings->nameEmpty) From 829f541d686997f6c98d4bf309791daef3c28f68 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Wed, 23 Dec 2020 08:35:23 -0500 Subject: [PATCH 21/41] MacOS notarization, minor tweaks --- console/nii_dicom.cpp | 49 ++++++++++++++++------------- console/nii_dicom.h | 2 +- console/nii_dicom_batch.cpp | 8 +++-- console/notarize.sh | 63 +++++++++++++++++++++++++++++++++++++ 4 files changed, 97 insertions(+), 25 deletions(-) create mode 100755 console/notarize.sh diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 6b158fc7..e187cf82 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -1238,7 +1238,8 @@ void checkSliceTimes(struct TCSAdata *CSA, int itemsOK, int isVerbose, bool is3D // Nov 1, 2018 wrote: // If you have an interleaved dataset we can more definitively validate this formula (aka sliceTime(i) - min(sliceTimes())). if (minTimeValue < 0) { - printWarning("Adjusting for negative MosaicRefAcqTimes (issue 271).\n"); + //printWarning("Adjusting for negative MosaicRefAcqTimes (issue 271).\n"); //if uncommented, overwhelming number of warnings (one per DICOM input), better once per series + CSA->sliceTiming[kMaxEPI3D-1] = -2.0; //issue 271: flag for unified warning for (int z = 0; z < itemsOK; z++) CSA->sliceTiming[z] = CSA->sliceTiming[z] - minTimeValue; } @@ -7022,29 +7023,31 @@ if (d.isHasPhase) int j = 0; if (d.xyzDim[3] > 1) j = 1; for (int i = 0; i < d.xyzDim[4]; i++) { + int slice = j+(i * d.xyzDim[3]); //dti4D->gradDynVol[i] = 0; //only PAR/REC - dti4D->TE[i] = dcmDim[j+(i * d.xyzDim[3])].TE; - dti4D->isPhase[i] = dcmDim[j+(i * d.xyzDim[3])].isPhase; - dti4D->isReal[i] = dcmDim[j+(i * d.xyzDim[3])].isReal; - dti4D->isImaginary[i] = dcmDim[j+(i * d.xyzDim[3])].isImaginary; - dti4D->triggerDelayTime[i] = dcmDim[j+(i * d.xyzDim[3])].triggerDelayTime; - dti4D->S[i].V[0] = dcmDim[j+(i * d.xyzDim[3])].V[0]; - dti4D->S[i].V[1] = dcmDim[j+(i * d.xyzDim[3])].V[1]; - dti4D->S[i].V[2] = dcmDim[j+(i * d.xyzDim[3])].V[2]; - dti4D->S[i].V[3] = dcmDim[j+(i * d.xyzDim[3])].V[3]; + dti4D->TE[i] = dcmDim[slice].TE; + dti4D->isPhase[i] = dcmDim[slice].isPhase; + dti4D->isReal[i] = dcmDim[slice].isReal; + dti4D->isImaginary[i] = dcmDim[slice].isImaginary; + dti4D->triggerDelayTime[i] = dcmDim[slice].triggerDelayTime; + dti4D->S[i].V[0] = dcmDim[slice].V[0]; + dti4D->S[i].V[1] = dcmDim[slice].V[1]; + dti4D->S[i].V[2] = dcmDim[slice].V[2]; + dti4D->S[i].V[3] = dcmDim[slice].V[3]; //printf("te=\t%g\tscl=\t%g\tintercept=\t%g\n",dti4D->TE[i], dti4D->intenScale[i],dti4D->intenIntercept[i]); if ((!isSameFloatGE(dti4D->TE[i],0.0)) && (dti4D->TE[i] != d.TE)) isTEvaries = true; if (dti4D->isPhase[i] != isPhase) d.isScaleOrTEVaries = true; if (dti4D->triggerDelayTime[i] != d.triggerDelayTime) d.isScaleOrTEVaries = true; if (dti4D->isReal[i] != isReal) d.isScaleOrTEVaries = true; if (dti4D->isImaginary[i] != isImaginary) d.isScaleOrTEVaries = true; - //dti4D->intenScale[i] = dcmDim[j+(i * d.xyzDim[3])].intenScale; - //dti4D->intenIntercept[i] = dcmDim[j+(i * d.xyzDim[3])].intenIntercept; - //dti4D->intenScalePhilips[i] = dcmDim[j+(i * d.xyzDim[3])].intenScalePhilips; - //dti4D->RWVIntercept[i] = dcmDim[j+(i * d.xyzDim[3])].RWVIntercept; - //dti4D->RWVScale[i] = dcmDim[j+(i * d.xyzDim[3])].RWVScale; - //if (dti4D->intenScale[i] != d.intenScale) isScaleVaries = true; - //if (dti4D->intenIntercept[i] != d.intenIntercept) isScaleVaries = true; + /*Philips can vary intensity scalings for separate slices within a volume! + dti4D->intenScale[i] = dcmDim[slice].intenScale; + dti4D->intenIntercept[i] = dcmDim[slice].intenIntercept; + dti4D->intenScalePhilips[i] = dcmDim[slice].intenScalePhilips; + dti4D->RWVIntercept[i] = dcmDim[slice].RWVIntercept; + dti4D->RWVScale[i] = dcmDim[slice].RWVScale; + if (dti4D->intenScale[i] != d.intenScale) isScaleVaries = true; + if (dti4D->intenIntercept[i] != d.intenIntercept) isScaleVaries = true;*/ } if((isScaleVaries) || (isTEvaries)) d.isScaleOrTEVaries = true; if (isTEvaries) d.isMultiEcho = true; @@ -7055,8 +7058,10 @@ if (d.isHasPhase) }*/ 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] ); + for (int i = 0; i < d.xyzDim[4]; i++) { + int slice = (i * d.xyzDim[3]); + printMessage(" %d TE=%g Slope=%g Inter=%g PhilipsScale=%g Phase=%d\n", i, dti4D->TE[i], dti4D->intenScale[slice], dti4D->intenIntercept[slice], dti4D->intenScalePhilips[slice], dti4D->isPhase[i] ); + } } } if ((d.xyzDim[3] == maxInStackPositionNumber) && (maxInStackPositionNumber > 1) && (d.zSpacing <= 0.0)) { @@ -7249,7 +7254,9 @@ if (d.isHasPhase) } // readDICOM() struct TDICOMdata readDICOM(char * fname) { - TDTI4D unused; - return readDICOMv(fname, false, kCompressSupport, &unused); + struct TDTI4D *dti4D = (struct TDTI4D *)malloc(sizeof(struct TDTI4D)); //unused + TDICOMdata ret = readDICOMv(fname, false, kCompressSupport, dti4D); + free(dti4D); + return ret; } // readDICOM() diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 5360cdd2..c31a24a9 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,7 +50,7 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20201216" +#define kDCMdate "v1.0.20201223" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index c3f6160d..bf6695ab 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -4467,6 +4467,8 @@ void checkSliceTiming(struct TDICOMdata * d, struct TDICOMdata * d1, int verbose while ((nSlices < kMaxEPI3D) && (d->CSA.sliceTiming[nSlices] >= 0.0)) nSlices++; if (nSlices < 1) return; + if (d->CSA.sliceTiming[kMaxEPI3D-1] < 1.0) + printWarning("Adjusting for negative MosaicRefAcqTimes (issue 271).\n"); bool isSliceTimeHHMMSS = (d->manufacturer == kMANUFACTURER_UIH); if (isForceSliceTimeHHMMSS) isSliceTimeHHMMSS = true; //if (d->isXA10A) isSliceTimeHHMMSS = true; //for XA10 use TimeAfterStart 0x0021,0x1104 -> Siemens de-identification can corrupt acquisition ties https://github.com/rordenlab/dcm2niix/issues/236 @@ -5516,9 +5518,11 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc printMessage("|%d|%s\n", i, nameList->str[dcmSort[i].indx]); } #endif + if (strlen(dcmList[dcmSort[0].indx].protocolName) < 1) //beware: tProtocolName can vary within a series "t1+AF8-mpr+AF8-ns+AF8-sag+AF8-p2+AF8-iso" vs "T1_mprage_ns_sag_p2_iso 1.0mm_192" + rescueProtocolName(&dcmList[dcmSort[0].indx], nameList->str[dcmSort[0].indx]); //move before headerDcm2Nii2 checkSliceTiming(&dcmList[indx0], &dcmList[indx1]); char pathoutname[2048] = {""}; - if (nii_createFilename(dcmList[dcmSort[0].indx], pathoutname, opts) == EXIT_FAILURE) { + if (nii_createFilename(dcmList[dcmSort[0].indx], pathoutname, opts) == EXIT_FAILURE) { //make sure call is subsequent to rescueProtocolName() free(imgM); return EXIT_FAILURE; } @@ -5580,8 +5584,6 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc free(img4D); saveAs3D = false; } - if (strlen(dcmList[dcmSort[0].indx].protocolName) < 1) //beware: tProtocolName can vary within a series "t1+AF8-mpr+AF8-ns+AF8-sag+AF8-p2+AF8-iso" vs "T1_mprage_ns_sag_p2_iso 1.0mm_192" - rescueProtocolName(&dcmList[dcmSort[0].indx], nameList->str[dcmSort[0].indx]); // Prevent these DICOM files from being reused. for(int i = 0; i < nConvert; ++i) dcmList[dcmSort[i].indx].converted2NII = 1; diff --git a/console/notarize.sh b/console/notarize.sh new file mode 100755 index 00000000..9a61a0b7 --- /dev/null +++ b/console/notarize.sh @@ -0,0 +1,63 @@ +#!/bin/bash +set -e + +#com.${COMPANY_NAME}.${APP_NAME} e.g. com.mricro.niimath +COMPANY_NAME=mycompany +APP_NAME=dcm2niix +APP_SPECIFIC_PASSWORD=abcd-efgh-ijkl-mnop +APPLE_ID_USER=myname@gmail.com +APPLE_ID_INSTALL="Developer ID Installer: My Name" +APPLE_ID_APP="Developer ID Application: My Name" + +if [[ "$APPLE_ID_USER" == "myname@gmail.com" ]] +then + echo "You need to set your personal IDs and password" + exit 1 +fi + +g++ -O3 -sectcreate __TEXT __info_plist Info.plist -I. main_console.cpp nii_foreign.cpp nii_dicom.cpp jpg_0XC3.cpp ujpeg.cpp nifti1_io_core.cpp nii_ortho.cpp nii_dicom_batch.cpp -o dcm2niixX86 -DmyDisableOpenJPEG -target x86_64-apple-macos10.12 -mmacosx-version-min=10.12 +g++ -O3 -sectcreate __TEXT __info_plist Info.plist -I. main_console.cpp nii_foreign.cpp nii_dicom.cpp jpg_0XC3.cpp ujpeg.cpp nifti1_io_core.cpp nii_ortho.cpp nii_dicom_batch.cpp -o dcm2niixARM -DmyDisableOpenJPEG -target arm64-apple-macos11 -mmacosx-version-min=11.0 +# Create the universal binary. +strip ./dcm2niixARM; strip ./dcm2niixX86 +lipo -create -output ${APP_NAME} dcm2niixARM dcm2niixX86 +rm ./dcm2niixARM; rm ./dcm2niixX86 +# Create a staging area for the installer package. +mkdir -p usr/local/bin +# Move the binary into the staging area. +mv ${APP_NAME} usr/local/bin +# Sign the binary. +codesign --timestamp --options=runtime -s "${APPLE_ID_APP}" -v usr/local/bin/${APP_NAME} +# Build the package. +pkgbuild --identifier "com.${COMPANY_NAME}.${APP_NAME}.pkg" --sign "${APPLE_ID_INSTALL}" --timestamp --root usr/local --install-location /usr/local/ ${APP_NAME}.pkg +# Submit the package to the notarization service. + +xcrun altool --notarize-app --primary-bundle-id "com.${COMPANY_NAME}.${APP_NAME}.pkg" --username $APPLE_ID_USER --password $APP_SPECIFIC_PASSWORD --file ${APP_NAME}.pkg --output-format xml > upload_log_file.txt + +# now we need to query apple's server to the status of notarization +# when the "xcrun altool --notarize-app" command is finished the output plist +# will contain a notarization-upload->RequestUUID key which we can use to check status +echo "Checking status..." +sleep 50 +REQUEST_UUID=`/usr/libexec/PlistBuddy -c "Print :notarization-upload:RequestUUID" upload_log_file.txt` +while true; do + xcrun altool --notarization-info $REQUEST_UUID -u $APPLE_ID_USER -p $APP_SPECIFIC_PASSWORD --output-format xml > request_log_file.txt + # parse the request plist for the notarization-info->Status Code key which will + # be set to "success" if the package was notarized + STATUS=`/usr/libexec/PlistBuddy -c "Print :notarization-info:Status" request_log_file.txt` + if [ "$STATUS" != "in progress" ]; then + break + fi + # echo $STATUS + echo "$STATUS" + sleep 10 +done + +# download the log file to view any issues +/usr/bin/curl -o log_file.txt `/usr/libexec/PlistBuddy -c "Print :notarization-info:LogFileURL" request_log_file.txt` + +# staple +echo "Stapling..." +xcrun stapler staple ${APP_NAME}.pkg +xcrun stapler validate ${APP_NAME}.pkg + +open log_file.txt From 0587941bb939fa46843917740f2c652a4219911f Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Wed, 23 Dec 2020 09:43:16 -0500 Subject: [PATCH 22/41] Overflow for Siemens data [missing protocol name] (https://github.com/rordenlab/dcm2niix/issues/466) --- console/nii_dicom.h | 2 +- console/nii_dicom_batch.cpp | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/console/nii_dicom.h b/console/nii_dicom.h index c31a24a9..cdfbeb05 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,7 +50,7 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20201223" +#define kDCMdate "v1.0.20201224" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index bf6695ab..8473eb06 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -968,6 +968,8 @@ void rescueProtocolName(struct TDICOMdata *d, const char * filename) { char protocolName[kDICOMStrLarge], fmriExternalInfo[kDICOMStrLarge], coilID[kDICOMStrLarge], consistencyInfo[kDICOMStrLarge], coilElements[kDICOMStrLarge], pulseSequenceDetails[kDICOMStrLarge], wipMemBlock[kDICOMStrLarge]; TCsaAscii csaAscii; siemensCsaAscii(filename, &csaAscii, d->CSA.SeriesHeader_offset, d->CSA.SeriesHeader_length, shimSetting, coilID, consistencyInfo, coilElements, pulseSequenceDetails, fmriExternalInfo, protocolName, wipMemBlock); + if (strlen(protocolName) >= kDICOMStr) + protocolName[kDICOMStr-1] = 0; strcpy(d->protocolName, protocolName); #endif } From 05a5d23c67d04df2c144f90266a782190860ff9c Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Thu, 24 Dec 2020 08:04:40 -0500 Subject: [PATCH 23/41] Update --- console/nii_dicom.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/console/nii_dicom.h b/console/nii_dicom.h index cdfbeb05..0eb6d9ba 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -44,7 +44,7 @@ extern "C" { #endif #if defined(__arm__) || defined(__ARM_ARCH) #define kCPUsuf " ARM" -#elif defined(__x86_64) +#elif defined(__x86_64) #define kCPUsuf " x86-64" #else #define kCPUsuf " " //unknown CPU From a7e384392066030e99befc9510b8398749250cab Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Wed, 13 Jan 2021 07:56:01 -0500 Subject: [PATCH 24/41] Preserve Series Description and Protocol Name for GE MRI (https://github.com/rordenlab/dcm2niix/issues/473) --- console/nii_dicom.cpp | 13 ++++++++----- console/nii_dicom.h | 2 +- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index e187cf82..adf58136 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -5317,7 +5317,6 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); break; } case kProtocolName : { - //if ((strlen(d.protocolName) < 1) || (d.manufacturer != kMANUFACTURER_GE)) //GE uses a generic session name here: do not overwrite kProtocolNameGE dcmStr(lLength, &buffer[lPos], d.protocolName); //see also kSequenceName break; } case kPatientOrient : @@ -6040,7 +6039,6 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); dcmStr(lLength, &buffer[lPos], d.scanOptions); break; case kSequenceName : { - //if (strlen(d.protocolName) < 1) //precedence given to kProtocolName and kProtocolNameGE dcmStr(lLength, &buffer[lPos], d.sequenceName); break; } @@ -6862,8 +6860,13 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); printWarning("0008,0008=MOSAIC but number of slices not specified: %s\n", fname); if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.CSA.dtiV[1] < -1.0) && (d.CSA.dtiV[2] < -1.0) && (d.CSA.dtiV[3] < -1.0)) d.CSA.dtiV[0] = 0; //SiemensTrio-Syngo2004A reports B=0 images as having impossible b-vectors. - if ((d.manufacturer == kMANUFACTURER_GE) && (d.modality == kMODALITY_MR) && (strlen(d.seriesDescription) > 1)) //GE uses a generic session name here: do not overwrite kProtocolNameGE + if ((d.manufacturer == kMANUFACTURER_GE) && (d.modality == kMODALITY_MR) && (strlen(d.seriesDescription) > 1)) {//GE uses a generic session name here: do not overwrite kProtocolNameGE + //issue 473: swap protocolName and seriesDescription + char tempTxt[kDICOMStr] = ""; + strcpy(tempTxt, d.protocolName); strcpy(d.protocolName, d.seriesDescription); + strcpy(d.seriesDescription, tempTxt); + } if ((strlen(d.protocolName) < 1) && (strlen(d.seriesDescription) > 1)) strcpy(d.protocolName, d.seriesDescription); if ((strlen(d.protocolName) > 1) && (isMoCo)) @@ -7145,7 +7148,7 @@ if (d.isHasPhase) d.triggerDelayTime = 0.0; //printf("%d\t%g\t%g\t%g\n", d.imageNum, d.acquisitionTime, d.triggerDelayTime, MRImageDynamicScanBeginTime); if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (strlen(seriesTimeTxt) > 1) && (d.isXA10A) && (d.xyzDim[3] == 1) && (d.xyzDim[4] < 2)) { - //printWarning("Ignoring series number of XA data saved as classic DICOM (issue 394)\n"); + //printWarning(">>Ignoring series number of XA data saved as classic DICOM (issue 394)\n"); d.isStackableSeries = true; d.imageNum += (d.seriesNum * 1000); strcpy(d.seriesInstanceUID, seriesTimeTxt); // dest <- src @@ -7193,7 +7196,7 @@ if (d.isHasPhase) } d.seriesUidCrc = mz_crc32X((unsigned char*) &d.seriesInstanceUID, strlen(d.seriesInstanceUID)); } - if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (strstr(d.sequenceName, "fl2d1") != NULL)) { + if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (strstr(d.sequenceName, "fl2d1") != NULL)) { d.isLocalizer = true; } //UIH 3D T1 scans report echo train length, which is interpreted as 3D EPI diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 0eb6d9ba..443e0975 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,7 +50,7 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20201224" +#define kDCMdate "v1.0.20210112" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic From e32fb1b958194c9ea6633724e18b3dd62b18b8c2 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Sun, 24 Jan 2021 14:45:42 -0500 Subject: [PATCH 25/41] GE TotalReadouTime, BEP009 fixes (https://github.com/rordenlab/dcm2niix/issues/476) --- GE/README.md | 10 +++++++++- console/nii_dicom.cpp | 5 +++++ console/nii_dicom.h | 2 +- console/nii_dicom_batch.cpp | 11 +++++++++-- 4 files changed, 24 insertions(+), 4 deletions(-) diff --git a/GE/README.md b/GE/README.md index cc05df1e..0761295b 100644 --- a/GE/README.md +++ b/GE/README.md @@ -34,7 +34,11 @@ Some sequences allow the user to interpolate images in plane (e.g. saving a 2D 6 ## Total Readout Time -One often wants to determine [echo spacing, bandwidth](https://support.brainvoyager.com/brainvoyager/functional-analysis-preparation/29-pre-processing/78-epi-distortion-correction-echo-spacing-and-bandwidth) and total read-out time for EPI data so they can be undistorted. Total readout time is influenced by parallel acceleration factor, bandwidth, number of EPI lines, and partial Fourier. Not all of these parameters are available from the GE DICOM images, so a user needs to check the scanner console. +One often wants to determine [echo spacing, bandwidth](https://support.brainvoyager.com/brainvoyager/functional-analysis-preparation/29-pre-processing/78-epi-distortion-correction-echo-spacing-and-bandwidth) and total read-out time for EPI data so they can be undistorted. Speifically, we are interested in FSL's definition of total read-out time, which may differ from the actual read-out time. FSL expects “the time from the middle of the first each to the middle of the last echo, as it would have been had partial k-space not been used”. So total read-out time is influenced by parallel acceleration factor, bandwidth, number of EPI lines, but not partial Fourier. For GE data we can use the Acquisition Matrix (0018,1310) in the phase-encoding direction, the in-plane acceleration ASSET R factor (0043,1083) and the Effective Echo Spacing (0043,102C). + +``` +TotalReadoutTime = ( ( ceil (0.5 * PE_AcquisitionMatrix / Asset_R_factor ) * 2) - 1 ] * EchoSpacing +``` ## Image Acceleration @@ -42,6 +46,10 @@ The BIDS [ParallelReductionFactorInPlane](https://bids-specification.readthedocs The BIDS [MultibandAccelerationFactor](https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html#rf-contrast) can be determined from the private tag `Multiband Parameters` (0043,10B6). This is an array with at least three values, the first is the Multiband (aka HyperBand) factor, the second is the Slice FOV Shift Factor, and the final is the Calibration method. For example, a scan using a multiband factor of 2 could be reported as `(0043,10b6) LO [2\4\19\\\\]`. +## Missing Information + +While GE DICOMs will report if the image uses partial Fourier, it will not reveal the partial Fourier fraction. + ## Phase-Encoding Polarity All EPI scans have spatial distortion, particularly those with longer readout times. Tools like [FSL topup](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/topup/TopupUsersGuide) can leverage data where two spin-echo images are acquired that are identical except for using opposite phase-encoding polarity (e.g. one uses A>P, the other P>A). Each image is distorted with the same magnitude, but in the opposite direction. GE's Rx27 software version and later populate the [Rectilinear Phase Encode Reordering (0018,9034)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9034)) tag which for EPI is set to either LINEAR or REVERSE_LINEAR. diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index adf58136..a7e76112 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -739,6 +739,11 @@ struct TDICOMdata clear_dicom_data() { strcpy(d.coilName, ""); strcpy(d.coilElements, ""); strcpy(d.radiopharmaceutical, ""); + strcpy(d.convolutionKernel, ""); + strcpy(d.unitsPT, ""); + strcpy(d.decayCorrection, ""); + strcpy(d.attenuationCorrectionMethod, ""); + strcpy(d.reconstructionMethod, ""); d.phaseEncodingLines = 0; //~ d.patientPositionSequentialRepeats = 0; //~ d.patientPositionRepeats = 0; diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 443e0975..a45585ee 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,7 +50,7 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20210112" +#define kDCMdate "v1.0.20210124" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 8473eb06..c60eb353 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -1547,8 +1547,6 @@ tse3d: T2*/ // effectiveEchoSpacing = d.CSA.sliceMeasurementDuration / (reconMatrixPE * 1000.0); if ((reconMatrixPE > 0) && (bandwidthPerPixelPhaseEncode > 0.0)) effectiveEchoSpacing = 1.0 / (bandwidthPerPixelPhaseEncode * reconMatrixPE); - if (d.effectiveEchoSpacingGE > 0.0) - effectiveEchoSpacing = d.effectiveEchoSpacingGE / 1000000.0; if ((effectiveEchoSpacing == 0.0) && (d.fieldStrength > 0) && (d.waterFatShift != 0.0) && (d.echoTrainLength > 0) && (reconMatrixPE > 1)) { json_Float(fp, "\t\"WaterFatShift\": %g,\n", d.waterFatShift); //https://github.com/rordenlab/dcm2niix/issues/377 @@ -1571,6 +1569,15 @@ tse3d: T2*/ json_Float(fp, "\t\"EstimatedEffectiveEchoSpacing\": %g,\n", effectiveEchoSpacingPhil); fprintf(fp, "\t\"EstimatedTotalReadoutTime\": %g,\n", totalReadoutTime); } + if (d.effectiveEchoSpacingGE > 0.0) { + //TotalReadoutTime = [ ceil (PE_AcquisitionMatrix / Asset_R_factor) - 1] * ESP + float totalReadoutTime = ((ceil (0.5 * d.phaseEncodingLines / d.accelFactPE) * 2.0) - 1.0) * d.effectiveEchoSpacingGE * 0.000001; + printf("ASSET= %g PE_AcquisitionMatrix= %d ESP= %d TotalReadoutTime= %g\n", d.accelFactPE, d.phaseEncodingLines, d.effectiveEchoSpacingGE, totalReadoutTime); + json_Float(fp, "\t\"TotalReadoutTime\": %g,\n", totalReadoutTime); + float effectiveEchoSpacingGE = totalReadoutTime / (reconMatrixPE - 1); + json_Float(fp, "\t\"EstimatedEffectiveEchoSpacing\": %g,\n", effectiveEchoSpacingGE); + effectiveEchoSpacing = 0.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 From 275969c7c3459e2f3bd943a3439fa4a646fac9d3 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Fri, 29 Jan 2021 17:05:06 -0500 Subject: [PATCH 26/41] Update explicit naming of DICOMs (https://github.com/rordenlab/dcm2niix/issues/252), GE file naming changes (https://github.com/rordenlab/dcm2niix/issues/476) --- GE/README.md | 20 ++++- console/main_console.cpp | 3 +- console/nifti1_io_core.cpp | 12 +++ console/nifti1_io_core.h | 1 + console/nii_dicom.cpp | 13 ++-- console/nii_dicom.h | 2 +- console/nii_dicom_batch.cpp | 149 +++++++++++++----------------------- 7 files changed, 97 insertions(+), 103 deletions(-) diff --git a/GE/README.md b/GE/README.md index 0761295b..13a6ef40 100644 --- a/GE/README.md +++ b/GE/README.md @@ -34,10 +34,26 @@ Some sequences allow the user to interpolate images in plane (e.g. saving a 2D 6 ## Total Readout Time -One often wants to determine [echo spacing, bandwidth](https://support.brainvoyager.com/brainvoyager/functional-analysis-preparation/29-pre-processing/78-epi-distortion-correction-echo-spacing-and-bandwidth) and total read-out time for EPI data so they can be undistorted. Speifically, we are interested in FSL's definition of total read-out time, which may differ from the actual read-out time. FSL expects “the time from the middle of the first each to the middle of the last echo, as it would have been had partial k-space not been used”. So total read-out time is influenced by parallel acceleration factor, bandwidth, number of EPI lines, but not partial Fourier. For GE data we can use the Acquisition Matrix (0018,1310) in the phase-encoding direction, the in-plane acceleration ASSET R factor (0043,1083) and the Effective Echo Spacing (0043,102C). +One often wants to determine [echo spacing, bandwidth](https://support.brainvoyager.com/brainvoyager/functional-analysis-preparation/29-pre-processing/78-epi-distortion-correction-echo-spacing-and-bandwidth) and total read-out time for EPI data so they can be undistorted. Speifically, we are interested in FSL's definition of total read-out time, which may differ from the actual read-out time. FSL expects “the time from the middle of the first each to the middle of the last echo, as it would have been had partial k-space not been used”. So total read-out time is influenced by parallel acceleration factor, bandwidth, number of EPI lines, but not partial Fourier. For GE data we can use the Acquisition Matrix (0018,1310) in the phase-encoding direction, the in-plane acceleration ASSET R factor (the reciprocal of this is stored as the first element of 0043,1083) and the Effective Echo Spacing (0043,102C). While GE does not tell us the [partial Fourier fraction](https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html), is does reveal if it is present with the ScanOptions (0018,1022) reporting [PFF](http://dicomlookup.com/lookup.asp?sw=Ttable&q=C.8-4) (in my experience, GE does not populate [(0018,9081)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9081))). While partial Fourier does not impact FSL's totalReadoutTime directly, it can interact with the number of lines acquired when combined with parallel imaging (the `Round_factor`). + +The formula for FSL's definition of TotalReadoutTime (in seconds) is: + +``` +TotalReadoutTime = ( ( ceil ((1/Round_factor) * PE_AcquisitionMatrix / Asset_R_factor ) * Round_factor) - 1 ] * EchoSpacing * 0.000001 +``` + +Consider an example: +``` +(0018,1310) US 128\0\0\128 # 8, 4 AcquisitionMatrix +(0018,0022) CS [SAT_GEMS\MP_GEMS\EPI_GEMS\ACC_GEMS\PFF\FS] # 42, 6 ScanOptions +(0043,102c) SS 636 # 2, 1 EchoSpacing +(0043,1083) DS [0.666667\1] # 10, 2 Acceleration +``` + +From this we can derive: ``` -TotalReadoutTime = ( ( ceil (0.5 * PE_AcquisitionMatrix / Asset_R_factor ) * 2) - 1 ] * EchoSpacing +ASSET= 1.5 PE_AcquisitionMatrix= 128 EchoSpacing= 636 Round_Factor= 4 TotalReadoutTime= 0.055332 ``` ## Image Acceleration diff --git a/console/main_console.cpp b/console/main_console.cpp index 1c135f2d..a049830e 100644 --- a/console/main_console.cpp +++ b/console/main_console.cpp @@ -102,7 +102,8 @@ void showHelp(const char * argv[], struct TDCMopts opts) { printf(" -p : Philips precise float (not display) scaling (y/n, default y)\n"); printf(" -r : rename instead of convert DICOMs (y/n, default n)\n"); printf(" -s : single file mode, do not convert other images in folder (y/n, default n)\n"); - printf(" -t : text notes includes private patient details (y/n, default n)\n"); + //text notes replaced with BIDS: this function is deprecated + //printf(" -t : text notes includes private patient details (y/n, default n)\n"); #if !defined(_WIN64) && !defined(_WIN32) //shell script for Unix only printf(" -u : up-to-date check\n"); #endif diff --git a/console/nifti1_io_core.cpp b/console/nifti1_io_core.cpp index 758e2c01..9935af68 100644 --- a/console/nifti1_io_core.cpp +++ b/console/nifti1_io_core.cpp @@ -192,6 +192,18 @@ vec3 nifti_vect33_norm (vec3 v) { //normalize vector length return vO; } +vec4 nifti_vect44_norm (vec4 v) { //normalize vector length + vec4 vO = v; + float vLen = sqrt( (v.v[0]*v.v[0]) + + (v.v[1]*v.v[1]) + + (v.v[2]*v.v[2]) + + (v.v[3]*v.v[3])); + if (vLen <= FLT_EPSILON) return vO; //avoid divide by zero + for (int i = 0; i < 4; i++) + vO.v[i] = v.v[i]/vLen; + return vO; +} + vec3 nifti_vect33mat33_mul(vec3 v, mat33 m ) { //multiply vector * 3x3matrix vec3 vO; for (int i=0; i<3; i++) { //multiply Pcrs * m diff --git a/console/nifti1_io_core.h b/console/nifti1_io_core.h index 1e15cbe3..98607888 100644 --- a/console/nifti1_io_core.h +++ b/console/nifti1_io_core.h @@ -68,6 +68,7 @@ mat44 nifti_mat44_inverse( mat44 R ); mat44 nifti_mat44_mul( mat44 A , mat44 B ); vec3 crossProduct(vec3 u, vec3 v); vec3 nifti_vect33_norm (vec3 v); +vec4 nifti_vect44_norm (vec4 v); vec3 nifti_vect33mat33_mul(vec3 v, mat33 m ); ivec3 setiVec3(int x, int y, int z); vec3 setVec3(float x, float y, float z); diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index a7e76112..ef71ae0e 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -5322,7 +5322,7 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); break; } case kProtocolName : { - dcmStr(lLength, &buffer[lPos], d.protocolName); //see also kSequenceName + dcmStr(lLength, &buffer[lPos], d.protocolName); break; } case kPatientOrient : dcmStr(lLength, &buffer[lPos], d.patientOrient); @@ -6027,7 +6027,7 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); //Canon/Toshiba omits 0018,9018, but reports [SE\EP];[GR\EP] for 0018,0020 //GE omits 0018,9018, but reports [EP\GR];[EP\SE] for 0018,0020 //Philips reports 0018,9018, but reports [SE];[GR] for 0018,0020 - if ((lLength > 1) && (strstr(d.scanningSequence, "IR") != NULL)) + if ((lLength > 1) && (strstr(d.scanningSequence, "IR") != NULL)) d.isIR = true; if ((lLength > 1) && (strstr(d.scanningSequence, "EP") != NULL)) d.isEPI = true; @@ -6042,6 +6042,8 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); } case kScanOptions: dcmStr(lLength, &buffer[lPos], d.scanOptions); + if ((lLength > 1) && (strstr(d.scanOptions, "PFF") != NULL)) + d.isPartialFourier = true; //e.g. GE does not populate (0018,9081) break; case kSequenceName : { dcmStr(lLength, &buffer[lPos], d.sequenceName); @@ -6865,13 +6867,14 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); printWarning("0008,0008=MOSAIC but number of slices not specified: %s\n", fname); if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.CSA.dtiV[1] < -1.0) && (d.CSA.dtiV[2] < -1.0) && (d.CSA.dtiV[3] < -1.0)) d.CSA.dtiV[0] = 0; //SiemensTrio-Syngo2004A reports B=0 images as having impossible b-vectors. - if ((d.manufacturer == kMANUFACTURER_GE) && (d.modality == kMODALITY_MR) && (strlen(d.seriesDescription) > 1)) {//GE uses a generic session name here: do not overwrite kProtocolNameGE - //issue 473: swap protocolName and seriesDescription + /* + if ((d.manufacturer == kMANUFACTURER_GE) && (d.modality == kMODALITY_MR) && (strlen(d.seriesDescription) > 1)) { + //issue 473, 476: swap protocolName and seriesDescription char tempTxt[kDICOMStr] = ""; strcpy(tempTxt, d.protocolName); strcpy(d.protocolName, d.seriesDescription); strcpy(d.seriesDescription, tempTxt); - } + }*/ if ((strlen(d.protocolName) < 1) && (strlen(d.seriesDescription) > 1)) strcpy(d.protocolName, d.seriesDescription); if ((strlen(d.protocolName) > 1) && (isMoCo)) diff --git a/console/nii_dicom.h b/console/nii_dicom.h index a45585ee..cc6e51e8 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,7 +50,7 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20210124" +#define kDCMdate "v1.0.20210129" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index c60eb353..a33d32bc 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -57,11 +57,9 @@ #endif #if defined(_WIN64) || defined(_WIN32) - #define myTextFileInputLists //comment out to disable this feature: https://github.com/rordenlab/dcm2niix/issues/288 const char kPathSeparator ='\\'; const char kFileSep[2] = "\\"; #else - #define myTextFileInputLists const char kPathSeparator ='/'; const char kFileSep[2] = "/"; #endif @@ -6249,77 +6247,6 @@ int singleDICOM(struct TDCMopts* opts, char *fname) { return ret; }// singleDICOM() -#ifdef myTextFileInputLists //https://github.com/rordenlab/dcm2niix/issues/288 -int textDICOM(struct TDCMopts* opts, char *fname) { - //check input file - FILE *fp = fopen(fname, "r"); - if (fp == NULL) -#ifdef USING_R - return EXIT_FAILURE; -#else - exit(EXIT_FAILURE); -#endif - int nConvert = 0; - char dcmname[2048]; - while (fgets(dcmname, sizeof(dcmname), fp)) { - int sz = strlen(dcmname); - if (sz > 0 && dcmname[sz-1] == '\n') dcmname[sz-1] = 0; //Unix LF - if (sz > 1 && dcmname[sz-2] == '\r') dcmname[sz-2] = 0; //Windows CR/LF - //if (isDICOMfile(dcmname) == 0) { //<- this will reject DICOM metadata not wrapped with a header - if ((!is_fileexists(dcmname)) || (!is_fileNotDir(dcmname)) ) { //<-this will accept meta data - fclose(fp); - printError("Problem with file '%s'\n", dcmname); - return EXIT_FAILURE; - } - //printf("%s\n", dcmname); - nConvert ++; - } - fclose(fp); - if (nConvert < 1) { - printError("No DICOM files found '%s'\n", dcmname); - return EXIT_FAILURE; - } - printMessage("Found %d DICOM file(s)\n", nConvert); - #ifndef USING_R - fflush(stdout); //show immediately if run from MRIcroGL GUI - #endif - TDCMsort * dcmSort = (TDCMsort *)malloc(nConvert * sizeof(TDCMsort)); - struct TDICOMdata *dcmList = (struct TDICOMdata *)malloc(nConvert * sizeof(struct TDICOMdata)); - struct TDTI4D *dti4D = (struct TDTI4D *)malloc(sizeof(struct TDTI4D)); - struct TSearchList nameList; - nameList.maxItems = nConvert; // larger requires more memory, smaller more passes - nameList.str = (char **) malloc((nameList.maxItems+1) * sizeof(char *)); //reserve one pointer (32 or 64 bits) per potential file - nameList.numItems = 0; - nConvert = 0; - fp = fopen(fname, "r"); - while (fgets(dcmname, sizeof(dcmname), fp)) { - int sz = strlen(dcmname); - if (sz > 0 && dcmname[sz-1] == '\n') dcmname[sz-1] = 0; //Unix LF - if (sz > 1 && dcmname[sz-2] == '\r') dcmname[sz-2] = 0; //Windows CR/LF - nameList.str[nameList.numItems] = (char *)malloc(strlen(dcmname)+1); - strcpy(nameList.str[nameList.numItems],dcmname); - nameList.numItems++; - dcmList[nConvert] = readDICOMv(nameList.str[nConvert], opts->isVerbose, opts->compressFlag, dti4D); //ignore compile warning - memory only freed on first of 2 passes - fillTDCMsort(dcmSort[nConvert], nConvert, dcmList[nConvert]); - nConvert ++; - } - fclose(fp); - qsort(dcmSort, nConvert, sizeof(struct TDCMsort), compareTDCMsort); //sort based on series and image numbers.... - int ret = saveDcm2Nii(nConvert, dcmSort, dcmList, &nameList, *opts, dti4D); - free(dcmSort); - free(dcmList); - free(dti4D); - freeNameList(nameList); - return ret; -}//textDICOM() - -#else //ifdef myTextFileInputLists -int textDICOM(struct TDCMopts* opts, char *fname) { - printError("Unable to parse txt files: re-compile with 'myTextFileInputLists' (see issue 288)"); - return EXIT_FAILURE; -} -#endif - size_t fileBytes(const char * fname) { FILE *fp = fopen(fname, "rb"); if (!fp) return 0; @@ -6684,25 +6611,57 @@ int nii_loadDirCore(char *indir, struct TDCMopts* opts) { #ifdef myTimer clock_t start = clock(); #endif - //1: find filenames of dicom files: up to two passes if we found more files than we allocated memory - for (int i = 0; i < 2; i++ ) { - nameList.str = (char **) malloc((nameList.maxItems+1) * sizeof(char *)); //reserve one pointer (32 or 64 bits) per potential file - nameList.numItems = 0; - searchDirForDICOM(indir, &nameList, opts->dirSearchDepth, 0, opts); - if (nameList.numItems <= nameList.maxItems) - break; - freeNameList(nameList); - nameList.maxItems = nameList.numItems+1; - //printMessage("Second pass required, found %ld images\n", nameList.numItems); - } - if (nameList.numItems < 1) { - if (opts->dirSearchDepth > 0) - printError("Unable to find any DICOM images in %s (or subfolders %d deep)\n", indir, opts->dirSearchDepth); - else //keep silent for dirSearchDepth = 0 - presumably searching multiple folders - {}; - free(nameList.str); //ignore compile warning - memory only freed on first of 2 passes - return kEXIT_NO_VALID_FILES_FOUND; - } + if ((is_fileNotDir(opts->indir)) && isExt(opts->indir, ".txt") ){ + nameList.str = (char **) malloc((nameList.maxItems+1) * sizeof(char *)); //reserve one pointer (32 or 64 bits) per potential file + nameList.numItems = 0; + FILE *fp = fopen(opts->indir, "r"); //textDICOM + if (fp == NULL) + return EXIT_FAILURE; + char dcmname[2048]; + while (fgets(dcmname, sizeof(dcmname), fp)) { + int sz = (int)strlen(dcmname); + if (sz > 0 && dcmname[sz-1] == '\n') dcmname[sz-1] = 0; //Unix LF + if (sz > 1 && dcmname[sz-2] == '\r') dcmname[sz-2] = 0; //Windows CR/LF + if ((!is_fileexists(dcmname)) || (!is_fileNotDir(dcmname)) ) { //<-this will accept meta data + fclose(fp); + printError("Problem with file '%s'\n", dcmname); + return EXIT_FAILURE; + } + if (nameList.numItems < nameList.maxItems) { + nameList.str[nameList.numItems] = (char *)malloc(strlen(dcmname)+1); + strcpy(nameList.str[nameList.numItems],dcmname); + } + nameList.numItems ++; + } + fclose(fp); + if (nameList.numItems >= nameList.maxItems) { + printError("Too many file names in '%s'\n", opts->indir); + return EXIT_FAILURE; + } + if (nameList.numItems < 1) + return kEXIT_NO_VALID_FILES_FOUND; + printMessage("Found %lu files in '%s'\n", nameList.numItems, opts->indir); + } else { + //1: find filenames of dicom files: up to two passes if we found more files than we allocated memory + for (int i = 0; i < 2; i++ ) { + nameList.str = (char **) malloc((nameList.maxItems+1) * sizeof(char *)); //reserve one pointer (32 or 64 bits) per potential file + nameList.numItems = 0; + searchDirForDICOM(indir, &nameList, opts->dirSearchDepth, 0, opts); + if (nameList.numItems <= nameList.maxItems) + break; + freeNameList(nameList); + nameList.maxItems = nameList.numItems+1; + //printMessage("Second pass required, found %ld images\n", nameList.numItems); + } + if (nameList.numItems < 1) { + if (opts->dirSearchDepth > 0) + printError("Unable to find any DICOM images in %s (or subfolders %d deep)\n", indir, opts->dirSearchDepth); + else //keep silent for dirSearchDepth = 0 - presumably searching multiple folders + {}; + free(nameList.str); //ignore compile warning - memory only freed on first of 2 passes + return kEXIT_NO_VALID_FILES_FOUND; + } + } size_t nDcm = nameList.numItems; printMessage( "Found %lu DICOM file(s)\n", nameList.numItems); //includes images and other non-image DICOMs #ifdef myTimer @@ -7043,8 +7002,10 @@ int nii_loadDir(struct TDCMopts* opts) { return convert_parRec(pname, *opts); }; } - if (isFile && (opts->isOnlySingleFile) && isExt(indir, ".txt") ) - return textDICOM(opts, indir); + if (isFile && (opts->isOnlySingleFile) && isExt(indir, ".txt") ) { + strcpy(opts->indir,indir); + return nii_loadDirCore(opts->indir, opts); + } if (opts->isRenameNotConvert) { int nConvert = searchDirRenameDICOM(opts->indir, opts->dirSearchDepth, 0, opts); if (nConvert < 0) return kEXIT_RENAME_ERROR; From 6445c1c4999fc5df544e1dd061cc1a6a214fbeb9 Mon Sep 17 00:00:00 2001 From: Jaemin Shin Date: Wed, 3 Feb 2021 10:36:24 -0500 Subject: [PATCH 27/41] GE Round factor for Total Readout Time --- GE/README.md | 2 +- console/nii_dicom_batch.cpp | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/GE/README.md b/GE/README.md index 13a6ef40..03453129 100644 --- a/GE/README.md +++ b/GE/README.md @@ -34,7 +34,7 @@ Some sequences allow the user to interpolate images in plane (e.g. saving a 2D 6 ## Total Readout Time -One often wants to determine [echo spacing, bandwidth](https://support.brainvoyager.com/brainvoyager/functional-analysis-preparation/29-pre-processing/78-epi-distortion-correction-echo-spacing-and-bandwidth) and total read-out time for EPI data so they can be undistorted. Speifically, we are interested in FSL's definition of total read-out time, which may differ from the actual read-out time. FSL expects “the time from the middle of the first each to the middle of the last echo, as it would have been had partial k-space not been used”. So total read-out time is influenced by parallel acceleration factor, bandwidth, number of EPI lines, but not partial Fourier. For GE data we can use the Acquisition Matrix (0018,1310) in the phase-encoding direction, the in-plane acceleration ASSET R factor (the reciprocal of this is stored as the first element of 0043,1083) and the Effective Echo Spacing (0043,102C). While GE does not tell us the [partial Fourier fraction](https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html), is does reveal if it is present with the ScanOptions (0018,1022) reporting [PFF](http://dicomlookup.com/lookup.asp?sw=Ttable&q=C.8-4) (in my experience, GE does not populate [(0018,9081)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9081))). While partial Fourier does not impact FSL's totalReadoutTime directly, it can interact with the number of lines acquired when combined with parallel imaging (the `Round_factor`). +One often wants to determine [echo spacing, bandwidth](https://support.brainvoyager.com/brainvoyager/functional-analysis-preparation/29-pre-processing/78-epi-distortion-correction-echo-spacing-and-bandwidth) and total read-out time for EPI data so they can be undistorted. Speifically, we are interested in FSL's definition of total read-out time, which may differ from the actual read-out time. FSL expects “the time from the middle of the first each to the middle of the last echo, as it would have been had partial k-space not been used”. So total read-out time is influenced by parallel acceleration factor, bandwidth, number of EPI lines, but not partial Fourier. For GE data we can use the Acquisition Matrix (0018,1310) in the phase-encoding direction, the in-plane acceleration ASSET R factor (the reciprocal of this is stored as the first element of 0043,1083) and the Effective Echo Spacing (0043,102C). While GE does not tell us the [partial Fourier fraction](https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html), is does reveal if it is present with the ScanOptions (0018,1022) reporting [PFF](http://dicomlookup.com/lookup.asp?sw=Ttable&q=C.8-4) (in my experience, GE does not populate [(0018,9081)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9081))). While partial Fourier does not impact FSL's totalReadoutTime directly, it can interact with the number of lines acquired when combined with parallel imaging (the `Round_factor` 2 (Full Fourier) or 4 (Partial Fourier)). The formula for FSL's definition of TotalReadoutTime (in seconds) is: diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index a33d32bc..524e5ca1 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -1569,7 +1569,9 @@ tse3d: T2*/ } if (d.effectiveEchoSpacingGE > 0.0) { //TotalReadoutTime = [ ceil (PE_AcquisitionMatrix / Asset_R_factor) - 1] * ESP - float totalReadoutTime = ((ceil (0.5 * d.phaseEncodingLines / d.accelFactPE) * 2.0) - 1.0) * d.effectiveEchoSpacingGE * 0.000001; + float roundFactor = 2.0; + if (d.isPartialFourier) roundFactor = 4.0; + float totalReadoutTime = ((ceil (1/roundFactor * d.phaseEncodingLines / d.accelFactPE) * roundFactor) - 1.0) * d.effectiveEchoSpacingGE * 0.000001; printf("ASSET= %g PE_AcquisitionMatrix= %d ESP= %d TotalReadoutTime= %g\n", d.accelFactPE, d.phaseEncodingLines, d.effectiveEchoSpacingGE, totalReadoutTime); json_Float(fp, "\t\"TotalReadoutTime\": %g,\n", totalReadoutTime); float effectiveEchoSpacingGE = totalReadoutTime / (reconMatrixPE - 1); From d0a99861c5960b1fd3160160f40ec6ffe5df6074 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Wed, 3 Feb 2021 15:32:57 -0500 Subject: [PATCH 28/41] EstimatedEffectiveEchoSpacing -> EffectiveEchoSpacing (https://github.com/rordenlab/dcm2niix/pull/480) --- FILENAMING.md | 2 +- console/nii_dicom.h | 2 +- console/nii_dicom_batch.cpp | 6 ++---- 3 files changed, 4 insertions(+), 6 deletions(-) diff --git a/FILENAMING.md b/FILENAMING.md index 620f9030..e08b01c2 100644 --- a/FILENAMING.md +++ b/FILENAMING.md @@ -19,7 +19,7 @@ You request the output file name with the `-f` argument. For example, consider y - %m=manufacturer short name (from 0008,0070: GE, Ph, Si, To, UI, NA) - %n=name of patient (from 0010,0010) - %o=mediaObjectInstanceUID (0002,0003)* - - %p=protocol name (from 0018,1030). If 0018,1030 is empty, or if the Manufacturer (0008,0070) is GE with Modality (0008,0060) of MR, then the SequenceName (0018,0024) is used if it is not empty. + - %p=protocol name (from 0018,1030). If 0018,1030 is empty, the SequenceName (0018,0024) is used. - %r=instance number (from 0020,0013)* - %s=series number (from 0020,0011) - %t=time of study (from 0008,0020 and 0008,0030) diff --git a/console/nii_dicom.h b/console/nii_dicom.h index cc6e51e8..c272823e 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,7 +50,7 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20210129" +#define kDCMdate "v1.0.20210202" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 524e5ca1..8c1d9b67 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -1572,11 +1572,9 @@ tse3d: T2*/ float roundFactor = 2.0; if (d.isPartialFourier) roundFactor = 4.0; float totalReadoutTime = ((ceil (1/roundFactor * d.phaseEncodingLines / d.accelFactPE) * roundFactor) - 1.0) * d.effectiveEchoSpacingGE * 0.000001; - printf("ASSET= %g PE_AcquisitionMatrix= %d ESP= %d TotalReadoutTime= %g\n", d.accelFactPE, d.phaseEncodingLines, d.effectiveEchoSpacingGE, totalReadoutTime); + //printf("ASSET= %g PE_AcquisitionMatrix= %d ESP= %d TotalReadoutTime= %g\n", d.accelFactPE, d.phaseEncodingLines, d.effectiveEchoSpacingGE, totalReadoutTime); json_Float(fp, "\t\"TotalReadoutTime\": %g,\n", totalReadoutTime); - float effectiveEchoSpacingGE = totalReadoutTime / (reconMatrixPE - 1); - json_Float(fp, "\t\"EstimatedEffectiveEchoSpacing\": %g,\n", effectiveEchoSpacingGE); - effectiveEchoSpacing = 0.0; + effectiveEchoSpacing = totalReadoutTime / (reconMatrixPE - 1); } json_Float(fp, "\t\"EffectiveEchoSpacing\": %g,\n", effectiveEchoSpacing); // Calculate true echo spacing (should match what Siemens reports on the console) From 2ad99528ba6a10501dbd05250c2a6abab02b1e9f Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Thu, 4 Feb 2021 10:41:29 -0500 Subject: [PATCH 29/41] report totalReadoutTime once --- GE/README.md | 1 + console/nii_dicom_batch.cpp | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/GE/README.md b/GE/README.md index 03453129..1dad29e9 100644 --- a/GE/README.md +++ b/GE/README.md @@ -40,6 +40,7 @@ The formula for FSL's definition of TotalReadoutTime (in seconds) is: ``` TotalReadoutTime = ( ( ceil ((1/Round_factor) * PE_AcquisitionMatrix / Asset_R_factor ) * Round_factor) - 1 ] * EchoSpacing * 0.000001 +EffectiveEchoSpacing = TotalReadoutTime/ (reconMatrixPE - 1) ``` Consider an example: diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 8c1d9b67..78178f70 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -1573,7 +1573,7 @@ tse3d: T2*/ if (d.isPartialFourier) roundFactor = 4.0; float totalReadoutTime = ((ceil (1/roundFactor * d.phaseEncodingLines / d.accelFactPE) * roundFactor) - 1.0) * d.effectiveEchoSpacingGE * 0.000001; //printf("ASSET= %g PE_AcquisitionMatrix= %d ESP= %d TotalReadoutTime= %g\n", d.accelFactPE, d.phaseEncodingLines, d.effectiveEchoSpacingGE, totalReadoutTime); - json_Float(fp, "\t\"TotalReadoutTime\": %g,\n", totalReadoutTime); + //json_Float(fp, "\t\"TotalReadoutTime\": %g,\n", totalReadoutTime); effectiveEchoSpacing = totalReadoutTime / (reconMatrixPE - 1); } json_Float(fp, "\t\"EffectiveEchoSpacing\": %g,\n", effectiveEchoSpacing); From 4a4ec3ed7f09d099e948498dca337fd544672aa7 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Tue, 9 Feb 2021 17:48:28 -0500 Subject: [PATCH 30/41] Forbid semicolon in filenames as Linux uses this for command concatenation and Windows function GetOpenFileName will truncate (https://github.com/rordenlab/dcm2niix/issues/425) --- FILENAMING.md | 5 +++-- GE/README.md | 3 ++- console/nii_dicom.h | 2 +- console/nii_dicom_batch.cpp | 2 +- 4 files changed, 7 insertions(+), 5 deletions(-) diff --git a/FILENAMING.md b/FILENAMING.md index e08b01c2..8a3a4a31 100644 --- a/FILENAMING.md +++ b/FILENAMING.md @@ -68,9 +68,9 @@ dcm2niix will attempt to write your image using the naming scheme you specify wi ## Special Characters -[Some characters are not permitted](https://stackoverflow.com/questions/1976007/what-characters-are-forbidden-in-windows-and-linux-directory-names) in file names. The following characters will be replaced with underscorces (`_`). Note that the forbidden characters vary between operating systems (Linux only forbids the forward slash, MacOS forbids forward slash and colon, while Windows forbids any of the characters listed below). To ensure that files can be easily copied between file systems, [dcm2niix restricts file names to characters allowed by Windows](https://github.com/rordenlab/dcm2niix/issues/237). +[Some characters are not permitted](https://stackoverflow.com/questions/1976007/what-characters-are-forbidden-in-windows-and-linux-directory-names) in file names. The following characters will be replaced with underscorces (`_`). Note that the forbidden characters vary between operating systems (Linux only forbids the forward slash, MacOS forbids forward slash and colon, while Windows forbids any of the characters listed below). To ensure that files can be easily copied between file systems, [dcm2niix restricts file names to characters allowed by Windows](https://github.com/rordenlab/dcm2niix/issues/237). While technically legal in all filesystems, the semicolon can wreak havoc in [Windows](https://stackoverflow.com/questions/3869594/semi-colons-in-windows-filenames) and [Linux](https://forums.plex.tv/t/linux-hates-semicolons-in-file-names/49098/2). -### List of Forbidden Characters (based on Windows) +### List of Forbidden Characters ``` < (less than) > (greater than) @@ -81,6 +81,7 @@ dcm2niix will attempt to write your image using the naming scheme you specify wi | (vertical bar or pipe) ? (question mark) * (asterisk) +; (semicolon) ``` [Control characters](https://en.wikipedia.org/wiki/ASCII#Control_characters) like backspace and tab are also forbidden. diff --git a/GE/README.md b/GE/README.md index 1dad29e9..6f6b8226 100644 --- a/GE/README.md +++ b/GE/README.md @@ -34,7 +34,7 @@ Some sequences allow the user to interpolate images in plane (e.g. saving a 2D 6 ## Total Readout Time -One often wants to determine [echo spacing, bandwidth](https://support.brainvoyager.com/brainvoyager/functional-analysis-preparation/29-pre-processing/78-epi-distortion-correction-echo-spacing-and-bandwidth) and total read-out time for EPI data so they can be undistorted. Speifically, we are interested in FSL's definition of total read-out time, which may differ from the actual read-out time. FSL expects “the time from the middle of the first each to the middle of the last echo, as it would have been had partial k-space not been used”. So total read-out time is influenced by parallel acceleration factor, bandwidth, number of EPI lines, but not partial Fourier. For GE data we can use the Acquisition Matrix (0018,1310) in the phase-encoding direction, the in-plane acceleration ASSET R factor (the reciprocal of this is stored as the first element of 0043,1083) and the Effective Echo Spacing (0043,102C). While GE does not tell us the [partial Fourier fraction](https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html), is does reveal if it is present with the ScanOptions (0018,1022) reporting [PFF](http://dicomlookup.com/lookup.asp?sw=Ttable&q=C.8-4) (in my experience, GE does not populate [(0018,9081)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9081))). While partial Fourier does not impact FSL's totalReadoutTime directly, it can interact with the number of lines acquired when combined with parallel imaging (the `Round_factor` 2 (Full Fourier) or 4 (Partial Fourier)). +One often wants to determine [echo spacing, bandwidth](https://support.brainvoyager.com/brainvoyager/functional-analysis-preparation/29-pre-processing/78-epi-distortion-correction-echo-spacing-and-bandwidth) and total read-out time for EPI data so they can be undistorted. Specifically, we are interested in FSL's definition of total read-out time, which may differ from the actual read-out time. FSL expects “the time from the middle of the first echo to the middle of the last echo, as it would have been had partial k-space not been used”. So total read-out time is influenced by parallel acceleration factor, bandwidth, number of EPI lines, but not partial Fourier. For GE data we can use the Acquisition Matrix (0018,1310) in the phase-encoding direction, the in-plane acceleration ASSET R factor (the reciprocal of this is stored as the first element of 0043,1083) and the Effective Echo Spacing (0043,102C). While GE does not tell us the [partial Fourier fraction](https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html), is does reveal if it is present with the ScanOptions (0018,1022) reporting [PFF](http://dicomlookup.com/lookup.asp?sw=Ttable&q=C.8-4) (in my experience, GE does not populate [(0018,9081)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9081))). While partial Fourier does not impact FSL's totalReadoutTime directly, it can interact with the number of lines acquired when combined with parallel imaging (the `Round_factor` 2 (Full Fourier) or 4 (Partial Fourier)). The formula for FSL's definition of TotalReadoutTime (in seconds) is: @@ -44,6 +44,7 @@ EffectiveEchoSpacing = TotalReadoutTime/ (reconMatrixPE - 1) ``` Consider an example: + ``` (0018,1310) US 128\0\0\128 # 8, 4 AcquisitionMatrix (0018,0022) CS [SAT_GEMS\MP_GEMS\EPI_GEMS\ACC_GEMS\PFF\FS] # 42, 6 ScanOptions diff --git a/console/nii_dicom.h b/console/nii_dicom.h index c272823e..0dac0574 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,7 +50,7 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20210202" +#define kDCMdate "v1.0.20210209" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 78178f70..39bd6763 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -2709,7 +2709,7 @@ int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopt #endif #if defined(_WIN64) || defined(_WIN32) || defined(kMASK_WINDOWS_SPECIAL_CHARACTERS)//https://stackoverflow.com/questions/1976007/what-characters-are-forbidden-in-windows-and-linux-directory-names for (size_t pos = 0; pos') || (outname[pos] == ':') + if ((outname[pos] == '\\') || (outname[pos] == '/') || (outname[pos] == ' ') || (outname[pos] == '<') || (outname[pos] == '>') || (outname[pos] == ':') || (outname[pos] == ';') || (outname[pos] == '"') // || (outname[pos] == '/') || (outname[pos] == '\\') //|| (outname[pos] == '^') issue398 || (outname[pos] == '*') || (outname[pos] == '|') || (outname[pos] == '?')) From 4f27d508680c4e59b19b2544f981981177071314 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Wed, 10 Feb 2021 07:45:48 -0500 Subject: [PATCH 31/41] Tell user to override merging (-m o) to separate coils, specific with fmrib usage that disrupts Siemens DCE processing (https://github.com/rordenlab/dcm2niix/issues/187) --- console/nii_dicom_batch.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 39bd6763..436c2aa6 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -6166,7 +6166,7 @@ bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, struct TDCMopts* opt if (d1.coilCrc != d2.coilCrc) { if (opts->isForceStackDCE) { if (!warnings->coilVaries) - printMessage("Slices stacked despite coil variation '%s' vs '%s'\n", d1.coilName, d2.coilName); + printMessage("Slices stacked despite coil variation '%s' vs '%s' (use '-m o' to turn off merging)\n", d1.coilName, d2.coilName); warnings->coilVaries = true; *isCoilVaries = true; } else { From 050e7a3bf5a32e6833462e80626855392a4074ca Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Mon, 22 Feb 2021 17:26:58 -0500 Subject: [PATCH 32/41] Kludge to separate vNav files as 3D (https://github.com/rordenlab/dcm2niix/issues/481) --- console/nii_dicom.h | 2 +- console/nii_dicom_batch.cpp | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 0dac0574..0ea9c9a9 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,7 +50,7 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20210209" +#define kDCMdate "v1.0.20210222" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 436c2aa6..0ef46f16 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -6187,6 +6187,7 @@ bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, struct TDCMopts* opt warnings->nameVaries = true; return false; } + if (( *isNonParallelSlices) && (d1.CSA.mosaicSlices > 1 )) return false; if ((!isSameFloatGE(d1.orient[1], d2.orient[1]) || !isSameFloatGE(d1.orient[2], d2.orient[2]) || !isSameFloatGE(d1.orient[3], d2.orient[3]) || !isSameFloatGE(d1.orient[4], d2.orient[4]) || !isSameFloatGE(d1.orient[5], d2.orient[5]) || !isSameFloatGE(d1.orient[6], d2.orient[6]) ) ) { if ((!warnings->orientVaries) && (!d1.isNonParallelSlices) && (!d1.isLocalizer)) @@ -6772,9 +6773,9 @@ int nii_loadDirCore(char *indir, struct TDCMopts* opts) { if (nConvert < 1) nConvert = 1; //prevents compiler warning for next line: never executed since j=i always causes nConvert ++ TDCMsort * dcmSort = (TDCMsort *)malloc(nConvert * sizeof(TDCMsort)); nConvert = 0; + isNonParallelSlices = false; for (int j = i; j < (int)nDcm; j++) { isMultiEcho = false; - isNonParallelSlices = false; isCoilVaries = false; if (isSameSet(dcmList[i], dcmList[j], opts, &warnings, &isMultiEcho, &isNonParallelSlices, &isCoilVaries)) { dcmList[j].converted2NII = 1; //do not reprocess repeats From ae16a61c81e158c9207501fdc12c322709ffe40b Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Mon, 22 Feb 2021 17:35:34 -0500 Subject: [PATCH 33/41] Add notes --- console/nii_dicom_batch.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 0ef46f16..da2c4106 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -6187,7 +6187,7 @@ bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, struct TDCMopts* opt warnings->nameVaries = true; return false; } - if (( *isNonParallelSlices) && (d1.CSA.mosaicSlices > 1 )) return false; + if (( *isNonParallelSlices) && (d1.CSA.mosaicSlices > 1 )) return false; //issue481 if ((!isSameFloatGE(d1.orient[1], d2.orient[1]) || !isSameFloatGE(d1.orient[2], d2.orient[2]) || !isSameFloatGE(d1.orient[3], d2.orient[3]) || !isSameFloatGE(d1.orient[4], d2.orient[4]) || !isSameFloatGE(d1.orient[5], d2.orient[5]) || !isSameFloatGE(d1.orient[6], d2.orient[6]) ) ) { if ((!warnings->orientVaries) && (!d1.isNonParallelSlices) && (!d1.isLocalizer)) @@ -6773,7 +6773,7 @@ int nii_loadDirCore(char *indir, struct TDCMopts* opts) { if (nConvert < 1) nConvert = 1; //prevents compiler warning for next line: never executed since j=i always causes nConvert ++ TDCMsort * dcmSort = (TDCMsort *)malloc(nConvert * sizeof(TDCMsort)); nConvert = 0; - isNonParallelSlices = false; + //isNonParallelSlices = false; //issue481 for (int j = i; j < (int)nDcm; j++) { isMultiEcho = false; isCoilVaries = false; From b7708063e50a11462bff82e461239137ae97888c Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Sun, 28 Feb 2021 16:00:18 -0500 Subject: [PATCH 34/41] If mosaics where ANY volume is non-planar, save ALL volumes as 3D(https://github.com/rordenlab/dcm2niix/issues/481) --- console/nii_dicom.h | 2 +- console/nii_dicom_batch.cpp | 29 ++++++++++++++++++++++++----- 2 files changed, 25 insertions(+), 6 deletions(-) diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 0ea9c9a9..bcaac4d1 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,7 +50,7 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20210222" +#define kDCMdate "v1.0.20210228" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index da2c4106..61592500 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -1375,7 +1375,7 @@ tse3d: T2*/ } //for k */ for (int k = 3; k < 11; k++) { //vessel locations char newstr[256]; - sprintf(newstr, "\t\"sWipMemBlock.AdFree%d\": %%g,\n", k); + sprintf(newstr, "\t\"sWipMemBlockAdFree%d\": %%g,\n", k); //issue483: sWipMemBlock.AdFree -> sWipMemBlockAdFree json_FloatNotNan(fp, newstr, csaAscii.adFree[k]); } } @@ -2665,8 +2665,8 @@ int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopt strcat (outname,newstr); isEchoReported = true; } - if ((dcm.isNonParallelSlices) && (!isImageNumReported)) { - sprintf(newstr, "_i%05d", dcm.imageNum); + if ((dcm.isNonParallelSlices) && (!isImageNumReported)) { + sprintf(newstr, "_i%05d", dcm.imageNum); strcat (outname,newstr); } /*if (dcm.maxGradDynVol > 0) { //Philips segmented @@ -4474,7 +4474,7 @@ void checkSliceTiming(struct TDICOMdata * d, struct TDICOMdata * d1, int verbose while ((nSlices < kMaxEPI3D) && (d->CSA.sliceTiming[nSlices] >= 0.0)) nSlices++; if (nSlices < 1) return; - if (d->CSA.sliceTiming[kMaxEPI3D-1] < 1.0) + if (d->CSA.sliceTiming[kMaxEPI3D-1] < 1.0) printWarning("Adjusting for negative MosaicRefAcqTimes (issue 271).\n"); bool isSliceTimeHHMMSS = (d->manufacturer == kMANUFACTURER_UIH); if (isForceSliceTimeHHMMSS) isSliceTimeHHMMSS = true; @@ -6773,10 +6773,10 @@ int nii_loadDirCore(char *indir, struct TDCMopts* opts) { if (nConvert < 1) nConvert = 1; //prevents compiler warning for next line: never executed since j=i always causes nConvert ++ TDCMsort * dcmSort = (TDCMsort *)malloc(nConvert * sizeof(TDCMsort)); nConvert = 0; - //isNonParallelSlices = false; //issue481 for (int j = i; j < (int)nDcm; j++) { isMultiEcho = false; isCoilVaries = false; + isNonParallelSlices = false; if (isSameSet(dcmList[i], dcmList[j], opts, &warnings, &isMultiEcho, &isNonParallelSlices, &isCoilVaries)) { dcmList[j].converted2NII = 1; //do not reprocess repeats fillTDCMsort(dcmSort[nConvert], j, dcmList[j]); @@ -6839,6 +6839,25 @@ int nii_loadDirCore(char *indir, struct TDCMopts* opts) { nConvert++; } } //for all images with same seriesUID as first one + if ((isNonParallelSlices) && (dcmList[ii].CSA.mosaicSlices > 1) && (nConvert > 0)) { //issue481: if ANY volumes are non-parallel, save ALL as 3D + printWarning("Saving mosaics with non-parallel slices as 3D (issue 481)\n"); + for (int j = i; j < (int)nDcm; j++) { + int ji = crcSort[j].indx; + if (dcmList[ii].seriesUidCrc != dcmList[ji].seriesUidCrc) break; + dcmList[ji].converted2NII = 1; + dcmList[ji].isNonParallelSlices = true; + if (isMultiEcho) dcmList[ji].isMultiEcho = true; + if (isCoilVaries) dcmList[ji].isCoilVaries = true; + struct TDCMsort dcmSort[1]; + fillTDCMsort(dcmSort[0], ji, dcmList[ji]); + int ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, dti4D); + if (ret == EXIT_SUCCESS) + nConvertTotal++; + else + convertError = true; + } + continue; + } //issue481 //issue 381: ensure all images are informed if there are variations in echo, parallel slices, coil name: if (isMultiEcho) for (int j = i; j <= jMax; j++) { From 3e2e3429275335105f4585166031dc3182de3ef1 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Thu, 4 Mar 2021 06:59:50 -0800 Subject: [PATCH 35/41] Use first recognized manufacturer tag (https://github.com/rordenlab/dcm2niix/issues/487) --- console/nii_dicom.cpp | 8 +++++--- console/nii_dicom.h | 2 +- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index ef71ae0e..b72d75db 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -1139,9 +1139,8 @@ int dcmStrManufacturer (const int lByteLength, unsigned char lBuffer[]) {//read // char cString[lByteLength + 1]; //#endif int ret = kMANUFACTURER_UNKNOWN; - cString[lByteLength] =0; + cString[lByteLength] = 0; memcpy(cString, (char*)&lBuffer[0], lByteLength); - //printMessage("MANU %s\n",cString); if ((toupper(cString[0])== 'S') && (toupper(cString[1])== 'I')) ret = kMANUFACTURER_SIEMENS; if ((toupper(cString[0])== 'G') && (toupper(cString[1])== 'E')) @@ -1159,6 +1158,8 @@ int dcmStrManufacturer (const int lByteLength, unsigned char lBuffer[]) {//read ret = kMANUFACTURER_UIH; if ((toupper(cString[0])== 'B') && (toupper(cString[1])== 'R')) ret = kMANUFACTURER_BRUKER; + if (ret == kMANUFACTURER_UNKNOWN) + printWarning("Unknown manufacturer %s\n",cString); //#ifdef _MSC_VER free(cString); //#endif @@ -5185,7 +5186,8 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); d.modality = kMODALITY_US; break; case kManufacturer: - d.manufacturer = dcmStrManufacturer (lLength, &buffer[lPos]); + if (d.manufacturer == kMANUFACTURER_UNKNOWN) + d.manufacturer = dcmStrManufacturer (lLength, &buffer[lPos]); volDiffusion.manufacturer = d.manufacturer; break; case kInstitutionName: diff --git a/console/nii_dicom.h b/console/nii_dicom.h index bcaac4d1..10780a92 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,7 +50,7 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20210228" +#define kDCMdate "v1.0.20210304" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic From 61f970de737f74a11c7d4e1e132cd1a14098d479 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Tue, 9 Mar 2021 10:48:35 -0500 Subject: [PATCH 36/41] Support Canon Enhanced DICOM (https://github.com/rordenlab/dcm2niix/issues/491) --- Canon/README.md | 10 ++++- README.md | 2 +- console/nii_dicom.cpp | 82 +++++++++++++++++++++++++++++-------- console/nii_dicom.h | 4 +- console/nii_dicom_batch.cpp | 19 +++++---- 5 files changed, 88 insertions(+), 29 deletions(-) diff --git a/Canon/README.md b/Canon/README.md index fdf51e98..ee651505 100644 --- a/Canon/README.md +++ b/Canon/README.md @@ -4,15 +4,21 @@ dcm2niix can convert Canon (né Toshiba) DICOM format images to NIfTI. This page ## Diffusion Weighted Imaging Notes -In contrast to several other vendors, Toshiba used public tags to report diffusion properties. Specifically, [DiffusionBValue (0018,9087)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9087)) and [DiffusionGradientOrientation (0018,9089)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9089)). Be aware that these tags are only populated for images where a diffusion gradient is applied. Consider a typical diffusion series where some volumes are acquired with B=0 while others have B=1000. In this case, only the volumes with B>0 will report a DiffusionBValue. +In contrast to several other vendors, Toshiba used public tags to report diffusion properties. Specifically, [DiffusionBValue (0018,9087)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9087)) and [DiffusionGradientOrientation (0018,9089)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9089)). Be aware that these tags are only populated for images where a diffusion gradient is applied. Consider a typical diffusion series where some volumes are acquired with B=0 while others have B=1000. In this case, only the volumes with B>0 will report a DiffusionBValue. These coordinates are with respect to the scanner bore, not image space. -Since the acquisition by Canon, these public tags are no longer populated. The diffusion gradient directions are now stored in the ASCII Image Comments tag. Like GE (but unlike [Siemens, GE and Toshiba](https://www.na-mic.org/wiki/NAMIC_Wiki:DTI:DICOM_for_DWI_and_DTI)), these directions are with respect to the image space, not the scanner bore. For empirical data see the Sample Datasets section. +Since the acquisition by Canon, these public tags are no longer populated for images saved in classic 2D DICOM format. The diffusion gradient directions are now stored in the ASCII Image Comments tag. Like GE (but unlike [Siemens, GE and Toshiba](https://www.na-mic.org/wiki/NAMIC_Wiki:DTI:DICOM_for_DWI_and_DTI)), these directions are with respect to the image space, not the scanner bore. Further, gradient direction is not adjusted for phase encoding polarity, and it is impossible to determine phase encoding polarity. For detailed discussion and a validation dataset that exhibits these attributes please see [dcm_qa_canon](https://github.com/neurolabusc/dcm_qa_canon). A Canon classic DICOM DWI image may report: ``` (0018,9087) FD 1500 # 8, 1 DiffusionBValue (0020,4000) LT [b=1500(0.445,0.000,0.895)] # 26, 1 ImageComments ``` +In contrast, when exporting images as enhanced (4D) DICOM, information is stored in public tags and does appear to compensate for phase encode polarity. These coordinates are with respect to the scanner bore, not image space. A Canon classic DICOM DWI image may report: + +``` +(0018,9087) FD 1500 # 8, 1 DiffusionBValue +(0018,9089) FD 0.29387456178665161\-0.95365142822265625\-0.064700603485107422 # 24, 3 DiffusionGradientOrientation +``` ## Unknown Properties diff --git a/README.md b/README.md index 1ebe4d1b..62ab7874 100644 --- a/README.md +++ b/README.md @@ -121,7 +121,7 @@ The following tools exploit dcm2niix - [abcd-dicom2bids](https://github.com/DCAN-Labs/abcd-dicom2bids) selectively downloads high quality ABCD datasets. - [autobids](https://github.com/khanlab/autobids) automates dcm2bids which uses dcm2niix. - [BiDirect_BIDS_Converter](https://github.com/wulms/BiDirect_BIDS_Converter) for conversion from DICOM to the BIDS standard. - - [BIDScoin](https://github.com/Donders-Institute/bidscoin) is a DICOM to BIDS converter with thorough [documentation](https://bircibrain.github.io/computingguide/docs/bids/bidscoin). + - [BIDScoin](https://github.com/Donders-Institute/bidscoin) is a DICOM to BIDS converter with a GUI and thorough [documentation](https://bidscoin.readthedocs.io). - [BIDS Toolbox](https://github.com/cardiff-brain-research-imaging-centre/bids-toolbox) is a web service for the creation and manipulation of BIDS datasets, using dcm2niix for importing DICOM data. - [birc-bids](https://github.com/bircibrain/birc-bids) provides a Docker/Singularity container with various BIDS conversion utilities. - [BOLD5000_autoencoder](https://github.com/nmningmei/BOLD5000_autoencoder) uses dcm2niix to pipe imaging data into an unsupervised machine learning algorithm. diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index b72d75db..9194a637 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -810,6 +810,7 @@ struct TDICOMdata clear_dicom_data() { d.is2DAcq = false; // d.isDerived = false; //0008,0008 = DERIVED,CSAPARALLEL,POSDISP d.isSegamiOasis = false; //these images do not store spatial coordinates + d.isBVecWorldCoordinates = false; //bvecs can be in image space (GE) or world coordinates (Siemens) d.isGrayscaleSoftcopyPresentationState = false; d.isRawDataStorage = false; d.isPartialFourier = false; @@ -4053,12 +4054,22 @@ bool compareTDCMdim (const TDCMdim &dcm1, const TDCMdim &dcm2) { return false; } //compareTDCMdim() +bool compareTDCMdimRev (const TDCMdim &dcm1, const TDCMdim &dcm2) { + for (int i = 0; i < MAX_NUMBER_OF_DIMENSIONS; i++) { + if(dcm1.dimIdx[i] < dcm2.dimIdx[i]) + return true; + else if(dcm1.dimIdx[i] > dcm2.dimIdx[i]) + return false; + } + return false; +} //compareTDCMdimRev() + #else 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 = 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]) @@ -4069,6 +4080,18 @@ int compareTDCMdim(void const *item1, void const *item2) { return 0; } //compareTDCMdim() +int compareTDCMdimRev(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++) { + if(dcm1->dimIdx[i] < dcm2->dimIdx[i]) + return -1; + else if(dcm1->dimIdx[i] > dcm2->dimIdx[i]) + return 1; + } + return 0; +} //compareTDCMdimRev() + #endif // USING_R struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, struct TDTI4D *dti4D) { @@ -4294,6 +4317,8 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16); #define kTriggerDelayTime 0x0020+uint32_t(0x9153<< 16 ) //FD #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 kTemporalPositionIndex 0x0020+uint32_t(0x9128<< 16 ) // UL #define kDimensionIndexPointer 0x0020+uint32_t(0x9165<< 16 ) //Private Group 21 as Used by Siemens: #define kSequenceVariant21 0x0021+(0x105B<< 16 )//CS @@ -4442,6 +4467,8 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); uint32_t dimensionIndexPointer[MAX_NUMBER_OF_DIMENSIONS]; size_t dimensionIndexPointerCounter = 0; int maxInStackPositionNumber = 0; + int temporalPositionIndex = 0; + int maxTemporalPositionIndex = 0; //int temporalPositionIdentifier = 0; int locationsInAcquisitionPhilips = 0; int imagesInAcquisition = 0; @@ -4605,7 +4632,7 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); is2005140FSQ = false; if (sqDepth < 0) sqDepth = 0; //should not happen, but protect for faulty anonymization //if we leave the folder MREchoSequence 0018,9114 - if (( nDimIndxVal > 0) && ((d.manufacturer == kMANUFACTURER_BRUKER) || (d.manufacturer == kMANUFACTURER_PHILIPS)) && (sqDepth00189114 >= sqDepth)) { + if (( nDimIndxVal > 0) && ((d.manufacturer == kMANUFACTURER_CANON) || (d.manufacturer == kMANUFACTURER_BRUKER) || (d.manufacturer == kMANUFACTURER_PHILIPS)) && (sqDepth00189114 >= sqDepth)) { sqDepth00189114 = -1; //triggered //printf("slice %d---> 0020,9157 = %d %d %d\n", inStackPositionNumber, d.dimensionIndexValues[0], d.dimensionIndexValues[1], d.dimensionIndexValues[2]); if (inStackPositionNumber > 0) { @@ -4640,7 +4667,8 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); // Bruker Enhanced MR IOD: reorder dimensions to ensure InStackPositionNumber corresponds to the first one // This will ensure correct ordering of slices in 4D datasets - if (d.manufacturer == kMANUFACTURER_BRUKER) { + /* + if (d.manufacturer == kMANUFACTURER_BRUKER) { for(size_t i = 1; i < dimensionIndexPointerCounter; i++){ if (dimensionIndexPointer[i] == kInStackPositionNumber){ //swap with first @@ -4648,8 +4676,9 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); dimensionIndexOrder[0] = i; } } - } + }*/ //Canon and Bruker reverse dimensionIndexItem order relative to Philips: new versions introduce compareTDCMdimRev int ndim = nDimIndxVal; + //printf("%d: %d %d %d %d\n", ndim, d.dimensionIndexValues[0], d.dimensionIndexValues[1], d.dimensionIndexValues[2], d.dimensionIndexValues[3]); for (int i = 0; i < ndim; i++) dcmDim[numDimensionIndexValues].dimIdx[i] = d.dimensionIndexValues[dimensionIndexOrder[i]]; dcmDim[numDimensionIndexValues].TE = TE; @@ -5593,12 +5622,16 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); 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_HITACHI) && (d.manufacturer != kMANUFACTURER_UNKNOWN) && (d.manufacturer != kMANUFACTURER_PHILIPS) && (d.manufacturer != kMANUFACTURER_BRUKER)) break; + if ((d.manufacturer != kMANUFACTURER_CANON) && (d.manufacturer != kMANUFACTURER_HITACHI) && (d.manufacturer != kMANUFACTURER_UNKNOWN) && (d.manufacturer != kMANUFACTURER_PHILIPS) && (d.manufacturer != kMANUFACTURER_BRUKER)) break; inStackPositionNumber = dcmInt(4,&buffer[lPos],d.isLittleEndian); //if (inStackPositionNumber == 1) numInStackPositionNumber1 ++; //printf("<%d>\n",inStackPositionNumber); if (inStackPositionNumber > maxInStackPositionNumber) maxInStackPositionNumber = inStackPositionNumber; break; + case kTemporalPositionIndex: + temporalPositionIndex = dcmInt(4,&buffer[lPos],d.isLittleEndian); + if (temporalPositionIndex > maxTemporalPositionIndex) maxTemporalPositionIndex = temporalPositionIndex; + break; case kDimensionIndexPointer: dimensionIndexPointer[dimensionIndexPointerCounter++] = dcmAttributeTag(&buffer[lPos],d.isLittleEndian); break; @@ -5617,13 +5650,13 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); case kDimensionIndexValues: { // kImageNum is not enough for 4D series from Philips 5.*. if (lLength < 4) break; nDimIndxVal = lLength / 4; - if(nDimIndxVal > MAX_NUMBER_OF_DIMENSIONS){ + if(nDimIndxVal > MAX_NUMBER_OF_DIMENSIONS){ printError("%d is too many dimensions. Only up to %d are supported\n", nDimIndxVal, MAX_NUMBER_OF_DIMENSIONS); nDimIndxVal = MAX_NUMBER_OF_DIMENSIONS; // Truncate } dcmMultiLongs(4 * nDimIndxVal, &buffer[lPos], nDimIndxVal, d.dimensionIndexValues, d.isLittleEndian); - break; } + break; } case kPhotometricInterpretation: { char interp[kDICOMStr]; dcmStr(lLength, &buffer[lPos], interp); @@ -6201,6 +6234,7 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); //d.CSA.dtiV[3] = v[2]; //printMessage("><>< 0018,9089: DWI bxyz %g %g %g %g\n", d.CSA.dtiV[0], d.CSA.dtiV[1], d.CSA.dtiV[2], d.CSA.dtiV[3]); hasDwiDirectionality = true; + d.isBVecWorldCoordinates = true; //e.g. Canon saved image space coordinates in Comments, world space in 0018, 9089 set_orientation0018_9089(&volDiffusion, lLength, &buffer[lPos], d.isLittleEndian); } break; @@ -6991,9 +7025,10 @@ if (d.isHasPhase) // printWarning("3D EPI with FrameAcquisitionDuration = %gs volumes = %d (see issue 369)\n", frameAcquisitionDuration/1000.0, d.xyzDim[4]); if (numDimensionIndexValues > 1) strcpy(d.imageType, imageType1st); //for multi-frame datasets, return name of book, not name of last chapter - if ((numDimensionIndexValues > 1) && (numDimensionIndexValues == numberOfFrames)) { + 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 maxVariableItem = 0; + if (true) { // int mn[MAX_NUMBER_OF_DIMENSIONS]; int mx[MAX_NUMBER_OF_DIMENSIONS]; for (int j = 0; j < MAX_NUMBER_OF_DIMENSIONS; j++) { @@ -7003,20 +7038,35 @@ if (d.isHasPhase) 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]; } + if (mx[j] != mn[j]) maxVariableItem = 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]); + if (isVerbose > 1) { + 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 + //see http://dicom.nema.org/medical/Dicom/2018d/output/chtml/part03/sect_C.8.24.3.3.html + //Philips puts spatial position as lower item than temporal position, the reverse is true for Bruker and Canon + int stackPositionItem = 0; + if (dimensionIndexPointerCounter > 0) + for(size_t i = 0; i < dimensionIndexPointerCounter; i++) + if (dimensionIndexPointer[i] == kInStackPositionNumber) stackPositionItem = i; //sort dimensions #ifdef USING_R - std::sort(dcmDim.begin(), dcmDim.begin() + numberOfFrames, compareTDCMdim); + if (stackPositionItem < maxVariableItem) + std::sort(dcmDim.begin(), dcmDim.begin() + numberOfFrames, compareTDCMdim); + else + std::sort(dcmDim.begin(), dcmDim.begin() + numberOfFrames, compareTDCMdimRev); #else - qsort(dcmDim, numberOfFrames, sizeof(struct TDCMdim), compareTDCMdim); + if (stackPositionItem < maxVariableItem) + qsort(dcmDim, numberOfFrames, sizeof(struct TDCMdim), compareTDCMdim); + else + qsort(dcmDim, numberOfFrames, sizeof(struct TDCMdim), compareTDCMdimRev); #endif //for (int i = 0; i < numberOfFrames; i++) - // printf("diskPos= %d dimIdx= %d %d %d %d TE= %g\n", i, dcmDim[i].diskPos, dcmDim[i].dimIdx[1], dcmDim[i].dimIdx[2], dcmDim[i].dimIdx[3], dti4D->TE[i]); + // printf("i %d diskPos= %d dimIdx= %d %d %d %d TE= %g\n", i, dcmDim[i].diskPos, dcmDim[i].dimIdx[0], dcmDim[i].dimIdx[1], dcmDim[i].dimIdx[2], dcmDim[i].dimIdx[3], dti4D->TE[i]); for (int i = 0; i < numberOfFrames; i++) { dti4D->sliceOrder[i] = dcmDim[i].diskPos; dti4D->intenScale[i] = dcmDim[i].intenScale; diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 10780a92..6ce08bf9 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,7 +50,7 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20210304" +#define kDCMdate "v1.0.20210308" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic @@ -192,7 +192,7 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; char institutionAddress[kDICOMStrLarge], imageComments[kDICOMStrLarge]; uint32_t dimensionIndexValues[MAX_NUMBER_OF_DIMENSIONS]; struct TCSAdata CSA; - bool isRealIsPhaseMapHz, isPrivateCreatorRemap, isHasOverlay, isEPI, isIR, isPartialFourier, isDiffusion, isVectorFromBMatrix, isRawDataStorage, isGrayscaleSoftcopyPresentationState, isStackableSeries, isCoilVaries, isNonParallelSlices, isSegamiOasis, isXA10A, isScaleOrTEVaries, isScaleVariesEnh, isDerived, isXRay, isMultiEcho, isValid, is3DAcq, is2DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase, isHasImaginary, isHasReal, isHasMagnitude,isHasMixed, isFloat, isResampled, isLocalizer; + bool isRealIsPhaseMapHz, isPrivateCreatorRemap, isHasOverlay, isEPI, isIR, isPartialFourier, isDiffusion, isVectorFromBMatrix, isRawDataStorage, isGrayscaleSoftcopyPresentationState, isStackableSeries, isCoilVaries, isNonParallelSlices, isBVecWorldCoordinates, isSegamiOasis, isXA10A, isScaleOrTEVaries, isScaleVariesEnh, isDerived, isXRay, isMultiEcho, isValid, is3DAcq, is2DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase, isHasImaginary, isHasReal, isHasMagnitude,isHasMixed, isFloat, isResampled, isLocalizer; char phaseEncodingRC, patientSex; }; diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 61592500..9d99ce3f 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -188,6 +188,7 @@ void geCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI *vx, int isV //COL then if swap the x and y value and reverse the sign on the z value. //If the phase encoding is not COL, then just reverse the sign on the x value. if ((d->manufacturer != kMANUFACTURER_GE) && (d->manufacturer != kMANUFACTURER_CANON)) return; + if (d->isBVecWorldCoordinates) return; //Canon classic DICOMs use image space, enhanced use world space! if ((!d->isEPI) && (d->CSA.numDti == 1)) d->CSA.numDti = 0; //issue449 if (d->CSA.numDti < 1) return; if ((toupper(d->patientOrient[0])== 'H') && (toupper(d->patientOrient[1])== 'F') && (toupper(d->patientOrient[2])== 'S')) @@ -204,7 +205,7 @@ void geCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI *vx, int isV return; } if (abs(sliceDir) != 3) - printWarning("GE DTI only tested for axial acquisitions (solution: use Xiangrui Li's dicm2nii)\n"); + printWarning("Limited validation for non-Axial DTI: confirm gradient vector transformation.\n"); //GE vectors from Xiangrui Li' dicm2nii, validated with datasets from https://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage#Diffusion_Tensor_Imaging ivec3 flp; if (abs(sliceDir) == 1) @@ -291,9 +292,9 @@ void siemensPhilipsCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI //convert DTI vectors from scanner coordinates to image frame of reference //Uses 6 orient values from ImageOrientationPatient (0020,0037) // requires PatientPosition 0018,5100 is HFS (head first supine) - if ((d->manufacturer != kMANUFACTURER_BRUKER) && (d->manufacturer != kMANUFACTURER_TOSHIBA) && (d->manufacturer != kMANUFACTURER_HITACHI) && (d->manufacturer != kMANUFACTURER_UIH) && (d->manufacturer != kMANUFACTURER_SIEMENS) && (d->manufacturer != kMANUFACTURER_PHILIPS)) return; + if ((!d->isBVecWorldCoordinates) && (d->manufacturer != kMANUFACTURER_BRUKER) && (d->manufacturer != kMANUFACTURER_TOSHIBA) && (d->manufacturer != kMANUFACTURER_HITACHI) && (d->manufacturer != kMANUFACTURER_UIH) && (d->manufacturer != kMANUFACTURER_SIEMENS) && (d->manufacturer != kMANUFACTURER_PHILIPS)) return; if (d->CSA.numDti < 1) return; - if (d->manufacturer == kMANUFACTURER_UIH) { + if (d->manufacturer == kMANUFACTURER_UIH) { for (int i = 0; i < d->CSA.numDti; i++) { vx[i].V[2] = -vx[i].V[2]; for (int v= 0; v < 4; v++) @@ -1545,9 +1546,10 @@ tse3d: T2*/ // effectiveEchoSpacing = d.CSA.sliceMeasurementDuration / (reconMatrixPE * 1000.0); if ((reconMatrixPE > 0) && (bandwidthPerPixelPhaseEncode > 0.0)) effectiveEchoSpacing = 1.0 / (bandwidthPerPixelPhaseEncode * reconMatrixPE); - if ((effectiveEchoSpacing == 0.0) && (d.fieldStrength > 0) && (d.waterFatShift != 0.0) && (d.echoTrainLength > 0) && (reconMatrixPE > 1)) { - json_Float(fp, "\t\"WaterFatShift\": %g,\n", d.waterFatShift); - //https://github.com/rordenlab/dcm2niix/issues/377 + json_Float(fp, "\t\"WaterFatShift\": %g,\n", d.waterFatShift); + if ((effectiveEchoSpacing == 0.0) && (d.imagingFrequency > 0.0) && (d.waterFatShift != 0.0) && (d.echoTrainLength > 0) && (reconMatrixPE > 1)) { + //in theory we could use either fieldStrength or imagingFrequency, but the former is typically provided with low precision + //https://github.com/rordenlab/dcm2niix/issues/377 // EchoSpacing 1/BW/EPI_factor https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=ind1308&L=FSL&D=0&P=113520 // this formula from https://support.brainvoyager.com/brainvoyager/functional-analysis-preparation/29-pre-processing/78-epi-distortion-correction-echo-spacing-and-bandwidth // https://neurostars.org/t/consolidating-epi-echo-spacing-and-readout-time-for-philips-scanner/4406 @@ -1561,11 +1563,12 @@ tse3d: T2*/ ReconMatrixPE = 0028,0010 or 0028,0011 depending on 0018,1312 */ - float actualEchoSpacing = d.waterFatShift / (d.imagingFrequency * 3.4 * (d.echoTrainLength + 1)); + float actualEchoSpacing = d.waterFatShift / (d.imagingFrequency * 3.4 * (d.echoTrainLength + 1)); float totalReadoutTime = actualEchoSpacing * d.echoTrainLength; float effectiveEchoSpacingPhil = totalReadoutTime / (reconMatrixPE - 1); json_Float(fp, "\t\"EstimatedEffectiveEchoSpacing\": %g,\n", effectiveEchoSpacingPhil); - fprintf(fp, "\t\"EstimatedTotalReadoutTime\": %g,\n", totalReadoutTime); + + fprintf(fp, "\t\"EstimatedTotalReadoutTime\": %g,\n", totalReadoutTime); } if (d.effectiveEchoSpacingGE > 0.0) { //TotalReadoutTime = [ ceil (PE_AcquisitionMatrix / Asset_R_factor) - 1] * ESP From 3dcc308543f1e1e99962779e5d163a3250b423dd Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Wed, 10 Mar 2021 18:54:39 -0500 Subject: [PATCH 37/41] Update dcm_qa submodule. --- Philips/README.md | 45 +++++++++++++++++++++++++++++++------------ console/nii_dicom.cpp | 1 + 2 files changed, 34 insertions(+), 12 deletions(-) diff --git a/Philips/README.md b/Philips/README.md index 6813e9b6..a5a0b817 100644 --- a/Philips/README.md +++ b/Philips/README.md @@ -20,19 +20,40 @@ Therefore, dcm2niix will ignore the IPP enclosed in 2005,140F unless no alternat dcm2niix losslessy copies the raw data from DICOM to NIfTI format. These values are typically stored as 16-bit integers in the range -32768..32767. Both the DICOM and NIfTI formats describe how scaling intercept and slope values can be used to convert these raw values into calibrated values. For example, with an intercept of 0 and slope of 0.01 the raw value of 50 would be converted to 0.5. -Unlike other vendors, Philips can store different scaling factors in their DICOM header. For most MRI modalities where the intensity brightness is relative, this has no impact. However, for modalities like ASL it can have an impact. The NIfTI format requires a single intensity intercept and slope is chosen. Therefore, dcm2niix will choose the "Real World" values if provided. If these are not available, dcm2niix will choose either the "precise" (default) or "display" (if the user choose "-p n") value. dcm2niix will also populate the folllowing tags in the BIDS header that allow the user to select between different intensity scaling formats: "PhilipsRescaleSlope", "PhilipsRescaleIntercept", "PhilipsScaleSlope", "UsePhilipsFloatNotDisplayScaling" (where "1" indicates NIfTI uses precise value, and "0" indicates display values)., "PhilipsRWVSlope" and "PhilipsRWVIntercept". - -The relevant DICOM tags are -RS = rescale slope ([0028,1053](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0028,1053))) -RI = rescale intercept ([0028,1052](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0028,1052))) -SS = scale slope (2005,100E) -RealWorldIntercept = (0040,9224) -Real World Slope = (0040,9225) -The transformation formulas are: -R = raw value, P = precise value, D = displayed value -D = R * RS + RI -P = D/(RS * SS) +The [NIfTI](https://nifti.nimh.nih.gov/pub/dist/src/niftilib/nifti1.h) header provides the `scl_slope` and `scl_inter` fields so each voxel value in the dataset is scaled as: +``` +I = scl_slope * R + scl_inter +``` + +where `R` is the voxel value stored and `I` is the "true" transformed voxel intensity. + +For all manufacturers except Philips, the `scl_slope` and `scl_inter` correspond to the DICOM tags [0028,1053](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0028,1053)) and [0028,1052](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0028,1053)), respectively. + +Philips has three possible intensity transforms for their DICOM images (world (`W`), display (`D`),precise (`P`)). All of these transforms might be provided in a single DICOM image, while the NIfTI standard only designates a single `scl_slope` and `scl_inter` for each image. dcm2niix will retain the raw value (`R`) and sets the NIfTI `scl_inter` and `scl_slope` values for the desired intensity transform. dcm2niix will use `W` if possible. If this is not possible it will use `P` unless the user specifies `-p n` to disable the precise calculation resulting in `D`. + +The formulates are provided below. The DICOM tags are in brackets (e.g. `(0040,9225)`) and the BIDS tag is in double quotes (e.g. `"PhilipsRWVSlope"`). Since all the scaling values are stored in the BIDS sidecar, you can always use these to generate your preferred intensity transform. + + +``` +Inputs: + R = raw stored value of voxel in DICOM without scaling + WS = RealWorldValue slope (0040,9225) "PhilipsRWVSlope" + WI = RealWorldValue intercept (0040,9224) "PhilipsRWVIntercept" + RS = rescale slope (0028,1053) "PhilipsRescaleSlope" + RI = rescale intercept (0028,1052) "PhilipsRescaleIntercept" + SS = scale slope (2005,100E) "PhilipsScaleSlope" +Outputs: + W = real world value + P = precise value + D = displayed value +Formulas: + W = R * WS + WI + D = R * RS + RI + P = D / (RS * SS) +``` + + ## Derived parametric maps stored with raw diffusion data Some Philips diffusion DICOM images include derived image(s) along with the images. Other manufacturers save these derived images as a separate series number, and the DICOM standard seems ambiguous on whether it is allowable to mix raw and derived data in the same series (see PS 3.3-2008, C.7.6.1.1.2-3). In practice, many Philips diffusion images append [derived parametric maps](http://www.revisemri.com/blog/2008/diffusion-tensor-imaging/) with the original data. With Philips, appending the derived isotropic image is optional - it is only created for the 'clinical' DTI schemes for radiography analysis and is triggered if the first three vectors in the gradient table are the unit X,Y and Z vectors. For conventional DWI, the result is the conventional mean of the ADC X,Y,Z for DTI it the conventional mean of the 3 principle Eigen vectors. As scientists, we want to discard these derived images, as they will disrupt data processing and we can generate better parametric maps after we have applied undistortion methods such as [Eddy and Topup](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/eddy/UsersGuide). The current version of dcm2niix uses the Diffusion Directionality (0018,9075) tag to detect B=0 unweighted ("NONE"), B-weighted ("DIRECTIONAL"), and derived ("ISOTROPIC") images. Note that the Dimension Index Values (0020,9157) tag provides an alternative approach to discriminate these images. Here are sample tags from a Philips enhanced image that includes and derived map (3rd dimension is "1" while the other images set this to "2"). diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 9194a637..0bf88ec0 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -6013,6 +6013,7 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); break; case kUnitsPT: //CS dcmStr(lLength, &buffer[lPos], d.unitsPT); + printf(">>>>>>>>%s<<<%s\n", d.unitsPT, fname); break; case kAttenuationCorrectionMethod: //LO dcmStr(lLength, &buffer[lPos], d.attenuationCorrectionMethod); From e9eda83b6feba59e778bbdb36b77353e14d20cb6 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Thu, 11 Mar 2021 11:37:22 -0500 Subject: [PATCH 38/41] Update Philips notes (https://github.com/rordenlab/dcm2niix/issues/493) --- Philips/README.md | 51 +++++++++++++++++++++++++++++++++++------------ README.md | 1 + UIH/README.md | 14 ++++++------- 3 files changed, 46 insertions(+), 20 deletions(-) diff --git a/Philips/README.md b/Philips/README.md index a5a0b817..afc9f660 100644 --- a/Philips/README.md +++ b/Philips/README.md @@ -1,6 +1,6 @@ ## 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.). +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 modified by DICOM tools (for example, anonymization, mangled by a [dcm4che/AGFA PACS](https://github.com/neurolabusc/dcm_qa_agfa), conversion of 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 [alternative DICOM to NIfTI converters](https://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage#Alternatives) or [report an issue](https://github.com/rordenlab/dcm2niix). @@ -18,26 +18,48 @@ Therefore, dcm2niix will ignore the IPP enclosed in 2005,140F unless no alternat ## Image Scaling -dcm2niix losslessy copies the raw data from DICOM to NIfTI format. These values are typically stored as 16-bit integers in the range -32768..32767. Both the DICOM and NIfTI formats describe how scaling intercept and slope values can be used to convert these raw values into calibrated values. For example, with an intercept of 0 and slope of 0.01 the raw value of 50 would be converted to 0.5. +How data is represented in DICOM for MR has several challenges and the technology and standard has evolved over the years to accommodate new uses. Unlike CT, where the signal is naturally displayed in Hounsfield units, MR has no natural signal units and the magnitude is influenced by the electronics and the software processing required to bring this to the final image. Secondly most of the original DICOM implementations used small bit number integers to store the underlying images for economy of storage. As a result it is necessary to apply scaling from the internal DICOM storage to a form suitable for radiographic display or quantitative measurement. There remain several challenges with this process, ensuring that the mapping to the integer values makes best use of the available bit depth for images with large dynamic range, or large changes between images, without clipping the data while also preserving the appearance of the noise field which is demanded by the needs of radiographic visual review. + +At its simplest this requires a rescale slope and intercept defined by the DICOM standard tags [0028,1053](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0028,1053)) and [0028,1052](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0028,1053)). Whether these values are the same for all images, or image specific depends upon the implementation and potentially the location of these tags withing the DICOM tag structure. + +In addition, the DICOM standard introduced the concept of [`real world units`](http://dicom.nema.org/dicom/2013/output/chtml/part03/sect_A.46.html). This allows the storage of one or more mappings to allow selective viewing of the data mapped into different value ranges (which may also be non-linear mappings). + +Clearly, scaling is a crucial aspect to be aware of for quantitative methods and which representation is used depends upon your needs. + +Philips thinks in terms of three different representations (using the terminology of the documentation available to Philips collaborators): + +| Name | ID | Description| +| ---------------- | ------------- | ------------- | +| Stored Value | SV | Raw data stored in DICOM tag PIXEL DATA tag (7FE0,0010)| +| Displayed Value | DV | The value which is shown to the user when using scanner interface, ROIS, measurements etc. | +| Floating Point | FP | An internal value at a point earlier in the reconstruction chain before the conversion to DICOM/integer for image presentation. | + +In general SV should not be used for quantitative measurements as it is an integer format. In practice, if the Rescale values are the same for all images (the typical case, but not guaranteed) SV can be used to compare signal intensities between images from the same scan. + +DV can be used for quantitative comparison of signal intensities between images in the same scan as long as the relevant rescale values are taken into account. These rescale values may come from the tags standard tags 0028,1053 and 0028,1052 or from a relevant RealWorld block if present. + +If the DV is derived from a RealWorld block with defined units (tag (0008,0104) such as Hz or ms rather than “no units”) or a RescaleType (0028,1054) with a non-US type (not defined by the standard), then the DV is already quantitative and cross scan comparison may be done. + +However, in general DV is not sufficient to compare images from different scans, especially if the signal intensity varies a lot (eg multiple inversion recovery scans) in which case the FP value may be used as this may be compared (with some caveats) across scans and across timescales. This scaling requires an additional scale factor on top of the DV value, the Scale Slope (private tag (2005,100E)) + +dcm2niix copies the raw pixel data from the DICOM tag 7FE0,0010) to NIfTI image. These values are typically stored as 16-bit integers in the range -32768..32767. Both the DICOM and NIfTI formats describe how scaling intercept and slope values can be used to convert these raw values into calibrated values. For example, with an intercept of 0 and slope of 0.01 the raw value of 50 would be converted to 0.5. The [NIfTI](https://nifti.nimh.nih.gov/pub/dist/src/niftilib/nifti1.h) header provides the `scl_slope` and `scl_inter` fields so each voxel value in the dataset is scaled as: ``` -I = scl_slope * R + scl_inter +I = scl_slope * SV + scl_inter ``` -where `R` is the voxel value stored and `I` is the "true" transformed voxel intensity. - -For all manufacturers except Philips, the `scl_slope` and `scl_inter` correspond to the DICOM tags [0028,1053](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0028,1053)) and [0028,1052](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0028,1053)), respectively. +where `SV` is the raw stored value and `I` is the "true" transformed voxel intensity. -Philips has three possible intensity transforms for their DICOM images (world (`W`), display (`D`),precise (`P`)). All of these transforms might be provided in a single DICOM image, while the NIfTI standard only designates a single `scl_slope` and `scl_inter` for each image. dcm2niix will retain the raw value (`R`) and sets the NIfTI `scl_inter` and `scl_slope` values for the desired intensity transform. dcm2niix will use `W` if possible. If this is not possible it will use `P` unless the user specifies `-p n` to disable the precise calculation resulting in `D`. +Philips has three possible intensity transforms for their DICOM images (world (`W`), display (`D`), precise (`P`)). All of these transforms might be provided in a single DICOM image, while the [NIfTI](https://nifti.nimh.nih.gov/pub/dist/src/niftilib/nifti1.h) header only designates a single `scl_slope` and `scl_inter` for each image. dcm2niix will retain the stored values (`SV`) and sets the NIfTI `scl_inter` and `scl_slope values` for the desired intensity transform. dcm2niix will use `W` if possible. If this is not possible it will use `P` unless the user specifies `-p n` to disable the precise calculation resulting in `D`. -The formulates are provided below. The DICOM tags are in brackets (e.g. `(0040,9225)`) and the BIDS tag is in double quotes (e.g. `"PhilipsRWVSlope"`). Since all the scaling values are stored in the BIDS sidecar, you can always use these to generate your preferred intensity transform. +The formulas are provided below. The DICOM tags are in brackets (e.g. `(0040,9225)`) and the BIDS tag is in double quotes (e.g. `"PhilipsRWVSlope"`). Since all the scaling values are stored in the BIDS sidecar, you can always use these to generate your preferred intensity transform. ``` Inputs: - R = raw stored value of voxel in DICOM without scaling + SV = stored value of DICOM PIXEL DATA without scaling WS = RealWorldValue slope (0040,9225) "PhilipsRWVSlope" WI = RealWorldValue intercept (0040,9224) "PhilipsRWVIntercept" RS = rescale slope (0028,1053) "PhilipsRescaleSlope" @@ -48,8 +70,8 @@ Outputs: P = precise value D = displayed value Formulas: - W = R * WS + WI - D = R * RS + RI + W = SV * WS + WI + D = SV * RS + RI P = D / (RS * SS) ``` @@ -146,5 +168,8 @@ Prior versions of dcm2niix used different methods to sort images. However, these ## Sample Datasets - [National Alliance for Medical Image Computing (NAMIC) samples](http://www.insight-journal.org/midas/collection/view/194) - - [Unusual Philips Examples](https://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage#Unusual_MRI). - - [Diffusion Examples](https://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage#Diffusion_Tensor_Imaging). \ No newline at end of file + - [Unusual Philips Examples (e.g. multi-echo)](https://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage#Unusual_MRI) + - [Archival samples](https://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage#Archival_MRI) + - [Diffusion Examples](https://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage#Diffusion_Tensor_Imaging) + - [Additional Diffusion Examples](https://github.com/neurolabusc/dcm_qa_philips) + - [Enhanced DICOMs](https://github.com/neurolabusc/dcm_qa_enh) \ No newline at end of file diff --git a/README.md b/README.md index 62ab7874..45e63545 100644 --- a/README.md +++ b/README.md @@ -40,6 +40,7 @@ There are a couple ways to install dcm2niix - Run the following command to get the latest version for Linux, Macintosh or Windows: * `curl -fLO https://github.com/rordenlab/dcm2niix/releases/latest/download/dcm2niix_lnx.zip` * `curl -fLO https://github.com/rordenlab/dcm2niix/releases/latest/download/dcm2niix_mac.zip` + * `curl -fLO https://github.com/rordenlab/dcm2niix/releases/latest/download/dcm2niix_mac_arm.pkg` * `curl -fLO https://github.com/rordenlab/dcm2niix/releases/latest/download/dcm2niix_win.zip` - [MRIcroGL (NITRC)](https://www.nitrc.org/projects/mricrogl) or [MRIcroGL (GitHub)](https://github.com/rordenlab/MRIcroGL12/releases) includes dcm2niix that can be run from the command line or from the graphical user interface (select the Import menu item). The Linux version of dcm2niix is compiled on a [holy build box](https://github.com/phusion/holy-build-box), so it should run on any Linux distribution. - If you have a MacOS computer with Homebrew you can run `brew install dcm2niix`. diff --git a/UIH/README.md b/UIH/README.md index 1c2f37a1..fbbf4fb7 100644 --- a/UIH/README.md +++ b/UIH/README.md @@ -12,7 +12,7 @@ UIH supports two ways of archiving the DWI/DTI and fMRI data. One way is one DIC Tag ID | Tag Name | VR | VM | Description | Sample -- | -- | -- | -- | -- | -- 0061,1002 | Generate Private | US | 1 | Flag to generate private format file | 1 -**0061,4002** | **FOV** | SH | 1 | FOV(mm) | 224*224 +0061,4002 | FOV | SH | 1 | FOV(mm) | 224*224 0065,1000 | MeasurmentUID | UL | 1 | Measurement UID of Protocol | 12547865 0065,1002 | ImageOrientationDisplayed | SH | 1 | Image Orientation Displayed | Sag or Sag>Cor 0065,1003 | ReceiveCoil | LO | 1 | Receive Coil Information | H 8 @@ -21,9 +21,9 @@ Tag ID | Tag Name | VR | VM | Description | Sample 0065,1006 | Slice Group ID | IS | 1 | Slice Group ID | 1 0065,1007 | Uprotocol | OB | 1 | Uprotocol value |   0065,1009 | BActualValue | FD | 1 | Actual B-Value from sequence | 1000.0 -**0065,100A** | **BUserValue** | FD | 1 | User Choose B-Value from UI | 1000.0 -**0065,100B** | **Block Size** | DS | 1 | Size of the paradigm/block | 10 -**0065,100C** | **Experimental status** | SH | 1 | fMRI | rest/active +0065,100A | BUserValue | FD | 1 | User Choose B-Value from UI | 1000.0 +0065,100B | Block Size | DS | 1 | Size of the paradigm/block | 10 +0065,100C | Experimental status | SH | 1 | fMRI | rest/active 0065,100D | Parallel Information | SH | 1 | ratio of parallel acquisition and acceleration |   0065,100F | Slice Position | SH | 1 | Slice location displayed on the screen | H23.4 0065,1011 | Sections | SH | 1 |   |   @@ -38,17 +38,17 @@ Tag ID | Tag Name | VR | VM | Description | Sample 0065,1029 | AcquisitionDuration | SH | 1 | Acquisition Duration | 0.03 0065,102B | ApplicationCategory | LT | 1 | Application names available | DTI\Func 0065,102C | RepeatitionIndex | IS | 1 |   | 0 -**0065,102D** | **SequenceDisplayName** | ST | 1 | Sequence display name | Epi_dti_b0 +0065,102D | SequenceDisplayName | ST | 1 | Sequence display name | Epi_dti_b0 0065,102E | NoiseDecovarFlag | LO | 1 | Noise decorrelation flag | PreWhite 0065,102F | ScaleFactor | FL | 1 | scale factor | 2.125 0065,1031 | MRSequenceVariant | SH | 1 | SequenceVariant |   0065,1032 | MRKSpaceFilter | SH | 1 | K space filter |   0065,1033 | MRTableMode | SH | 1 | Table mode | Fix 0065,1036 | MRDiscoParameter | OB | 1 |   |   -**0065,1037** | **MRDiffusionGradOrientation** | FD | 3 | Diffusion gradient orientation | 0\0\0 +0065,1037 | MRDiffusionGradOrientation | FD | 3 | Diffusion gradient orientation | 0\0\0 0065,1038 | MRPerfusionNoiseLevel | FD | 1 | epi_dwi/perfusion noise level | 40 0065,1039 | MRGradRange | SH | 6 | linear range of gradient | 0.0\157\0.0\157\0.0\125 -**0065,1050** | **MR Number Of Slice In Volume** | DS | 1 | Number Of Frames In a Volume,Columns of each frame: cols =ceil(sqrt(total)) ; Rows of each frame: rows =ceil(total/cols) ; appeared when image type (00080008) has VFRAME | 27 +0065,1050 | MR Number Of Slice In Volume | DS | 1 | Number Of Frames In a Volume,Columns of each frame: cols =ceil(sqrt(total)) ; Rows of each frame: rows =ceil(total/cols) ; appeared when image type (00080008) has VFRAME | 27 0065,1051 | MR VFrame Sequence | SQ | 1 | 1 |   ->0008,0022 | Acquisition Date | DA | 1 |   |   ->0008,0032 | Acquisition Time | TM | 1 |   |   From 93e1373e6dec869c330a2dc5747ce7ad80bdb81b Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Wed, 17 Mar 2021 11:53:28 -0400 Subject: [PATCH 39/41] Describe and provide kludge for mangled Canon DICOMs (https://github.com/rordenlab/dcm2niix/issues/495) --- Canon/README.md | 4 ++++ Philips/README.md | 39 ++++++++++++++++--------------------- console/main_console.cpp | 9 +++++++-- console/nii_dicom.h | 2 +- console/nii_dicom_batch.cpp | 10 +++++++--- console/nii_dicom_batch.h | 2 +- 6 files changed, 37 insertions(+), 29 deletions(-) diff --git a/Canon/README.md b/Canon/README.md index ee651505..7ee79845 100644 --- a/Canon/README.md +++ b/Canon/README.md @@ -2,6 +2,10 @@ dcm2niix can convert Canon (né Toshiba) DICOM format images to NIfTI. This page notes vendor specific conversion details. +## Avoid Classic DICOM + +Users of Canon MRI equipment are strongly advised to export data from their scanners as enhanced DICOM (with all images from the series stored as a single file) rather than classic DICOM (each 2D slice stored as a separate file). Limitations of the Canon classic DICOMs are described [here](https://github.com/rordenlab/dcm2niix/issues/495) and [here](https://github.com/neurolabusc/dcm_qa_canon). + ## Diffusion Weighted Imaging Notes In contrast to several other vendors, Toshiba used public tags to report diffusion properties. Specifically, [DiffusionBValue (0018,9087)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9087)) and [DiffusionGradientOrientation (0018,9089)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9089)). Be aware that these tags are only populated for images where a diffusion gradient is applied. Consider a typical diffusion series where some volumes are acquired with B=0 while others have B=1000. In this case, only the volumes with B>0 will report a DiffusionBValue. These coordinates are with respect to the scanner bore, not image space. diff --git a/Philips/README.md b/Philips/README.md index afc9f660..97fa3138 100644 --- a/Philips/README.md +++ b/Philips/README.md @@ -18,13 +18,11 @@ Therefore, dcm2niix will ignore the IPP enclosed in 2005,140F unless no alternat ## Image Scaling -How data is represented in DICOM for MR has several challenges and the technology and standard has evolved over the years to accommodate new uses. Unlike CT, where the signal is naturally displayed in Hounsfield units, MR has no natural signal units and the magnitude is influenced by the electronics and the software processing required to bring this to the final image. Secondly most of the original DICOM implementations used small bit number integers to store the underlying images for economy of storage. As a result it is necessary to apply scaling from the internal DICOM storage to a form suitable for radiographic display or quantitative measurement. There remain several challenges with this process, ensuring that the mapping to the integer values makes best use of the available bit depth for images with large dynamic range, or large changes between images, without clipping the data while also preserving the appearance of the noise field which is demanded by the needs of radiographic visual review. +How data is represented in DICOM for MR has several challenges and the technology and standard has evolved over the years to accommodate new uses. Unlike CT, where the signal is naturally displayed in Hounsfield units, MR has no natural signal units and the magnitude is influenced by the electronics and the software processing required to bring this to the final image. Secondly most of the original DICOM implementations used small bit number integers to store the underlying images for economy of storage. As a result it is necessary to apply scaling from the internal DICOM storage to a form suitable for radiographic display or quantitative measurement. There remain several challenges with this process, ensuring that the mapping to the integer values makes best use of the available bit depth for images with large dynamic range, or large changes between images, without clipping the data while also preserving the appearance of the noise field which is demanded by the needs of radiographic visual review. Note that for most MRI modalities these concerns do not impact analyses: the intensity is assumed arbitrary, the statistics treat signal offset and scaling as nuisance regressors when fitting models, and cacluations are computed with high precision floating point numbers. However, there are some situations such as arterial spin labeling where image scaling is important. In these situations, scaling is a crucial aspect to be aware of for quantitative methods and which representation is used depends upon your needs. -At its simplest this requires a rescale slope and intercept defined by the DICOM standard tags [0028,1053](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0028,1053)) and [0028,1052](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0028,1053)). Whether these values are the same for all images, or image specific depends upon the implementation and potentially the location of these tags withing the DICOM tag structure. +At its simplest image scaling requires a rescale slope and intercept defined by the DICOM standard tags [0028,1053](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0028,1053)) and [0028,1052](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0028,1053)). Whether these values are the same for all images, or image specific depends upon the implementation and potentially the location of these tags withing the DICOM tag structure. For manufacturers other than Philips, these are the only intensity scaling values provided, so there is no concern regarding which scaling values should be used. -In addition, the DICOM standard introduced the concept of [`real world units`](http://dicom.nema.org/dicom/2013/output/chtml/part03/sect_A.46.html). This allows the storage of one or more mappings to allow selective viewing of the data mapped into different value ranges (which may also be non-linear mappings). - -Clearly, scaling is a crucial aspect to be aware of for quantitative methods and which representation is used depends upon your needs. +However, the DICOM standard introduced the concept of [`real world units`](http://dicom.nema.org/dicom/2013/output/chtml/part03/sect_A.46.html). This allows the storage of one or more mappings to allow selective viewing of the data mapped into different value ranges (which may also be non-linear mappings). Philips thinks in terms of three different representations (using the terminology of the documentation available to Philips collaborators): @@ -33,16 +31,15 @@ Philips thinks in terms of three different representations (using the terminolog | Stored Value | SV | Raw data stored in DICOM tag PIXEL DATA tag (7FE0,0010)| | Displayed Value | DV | The value which is shown to the user when using scanner interface, ROIS, measurements etc. | | Floating Point | FP | An internal value at a point earlier in the reconstruction chain before the conversion to DICOM/integer for image presentation. | +| Real World Value | WV | DICOM defined real world units| -In general SV should not be used for quantitative measurements as it is an integer format. In practice, if the Rescale values are the same for all images (the typical case, but not guaranteed) SV can be used to compare signal intensities between images from the same scan. - -DV can be used for quantitative comparison of signal intensities between images in the same scan as long as the relevant rescale values are taken into account. These rescale values may come from the tags standard tags 0028,1053 and 0028,1052 or from a relevant RealWorld block if present. +In general SV should not be used for quantitative measurements as it is an integer format. In practice, if the Rescale values are the same for all images (the typical case, but not guaranteed) SV can be used to compare signal intensities between images from the same scan. Note that the NIfTI format only provides a single `scl_slope` and `scl_inter` for the entire file, whereas in DICOM rescale values can in theory differ across 2D slices. Therefore, in situations where the rescale values do differ across slices, dcm2niix will apply the requested rescale to each slice and save the scaled data as the 32-bit float NIfTI dataset. This preserves the varibility reported by the rescale tags, at the cost of disk space. -If the DV is derived from a RealWorld block with defined units (tag (0008,0104) such as Hz or ms rather than “no units”) or a RescaleType (0028,1054) with a non-US type (not defined by the standard), then the DV is already quantitative and cross scan comparison may be done. +DV can be used for quantitative comparison of signal intensities between images in the same scan as long as the relevant rescale values are taken into account. These rescale values may come from the tags standard tags 0028,1053 and 0028,1052 or from a relevant RealWorld block if present. If the DV is derived from a RealWorld block with defined units (tag (0008,0104) such as Hz or ms rather than “no units”) or a RescaleType (0028,1054) with a non-US type (not defined by the standard), then the DV is already quantitative and cross scan comparison may be done. However, in general DV is not sufficient to compare images from different scans, especially if the signal intensity varies a lot (eg multiple inversion recovery scans) in which case the FP value may be used as this may be compared (with some caveats) across scans and across timescales. This scaling requires an additional scale factor on top of the DV value, the Scale Slope (private tag (2005,100E)) -dcm2niix copies the raw pixel data from the DICOM tag 7FE0,0010) to NIfTI image. These values are typically stored as 16-bit integers in the range -32768..32767. Both the DICOM and NIfTI formats describe how scaling intercept and slope values can be used to convert these raw values into calibrated values. For example, with an intercept of 0 and slope of 0.01 the raw value of 50 would be converted to 0.5. +As long as rescale values are identical across all DICOM slices, dcm2niix losslessly copies the raw pixel data from the DICOM tag (7FE0,0010) to NIfTI image. These values are typically stored as 16-bit integers in the range -32768..32767. Both the DICOM and NIfTI formats describe how scaling intercept and slope values can be used to convert these raw values into calibrated values. For example, with an intercept of 0 and slope of 0.01 the raw value of 50 would be converted to 0.5. The [NIfTI](https://nifti.nimh.nih.gov/pub/dist/src/niftilib/nifti1.h) header provides the `scl_slope` and `scl_inter` fields so each voxel value in the dataset is scaled as: @@ -52,10 +49,9 @@ I = scl_slope * SV + scl_inter where `SV` is the raw stored value and `I` is the "true" transformed voxel intensity. -Philips has three possible intensity transforms for their DICOM images (world (`W`), display (`D`), precise (`P`)). All of these transforms might be provided in a single DICOM image, while the [NIfTI](https://nifti.nimh.nih.gov/pub/dist/src/niftilib/nifti1.h) header only designates a single `scl_slope` and `scl_inter` for each image. dcm2niix will retain the stored values (`SV`) and sets the NIfTI `scl_inter` and `scl_slope values` for the desired intensity transform. dcm2niix will use `W` if possible. If this is not possible it will use `P` unless the user specifies `-p n` to disable the precise calculation resulting in `D`. +Philips has three possible intensity transforms for their DICOM images (world (`W`), display (`D`), precise (`P`)). All of these transforms might be provided in a single DICOM image, while the [NIfTI](https://nifti.nimh.nih.gov/pub/dist/src/niftilib/nifti1.h) header only designates a single `scl_slope` and `scl_inter` for each image. dcm2niix will attempt to retain the stored values (`SV`) and sets the NIfTI `scl_inter` and `scl_slope values` for the desired intensity transform. dcm2niix will use `FP` if possible. If this is not possibleor the user specifies `-p n` dcm2niix will use the transforms for `DV`. -The formulas are provided below. The DICOM tags are in brackets (e.g. `(0040,9225)`) and the BIDS tag is in double quotes (e.g. `"PhilipsRWVSlope"`). Since all the scaling values are stored in the BIDS sidecar, you can always use these to generate your preferred intensity transform. - +The formulas are provided below. The DICOM tags are in brackets (e.g. `(0040,9225)`) and the BIDS tag is in double quotes (e.g. `"PhilipsRWVSlope"`). Since all the scaling values are stored in the BIDS sidecar, you can always use these to later your preferred intensity transform (assume all slices used the same scaling values). ``` Inputs: @@ -66,19 +62,18 @@ Inputs: RI = rescale intercept (0028,1052) "PhilipsRescaleIntercept" SS = scale slope (2005,100E) "PhilipsScaleSlope" Outputs: - W = real world value - P = precise value - D = displayed value + WV = real world value + FP = precise value + DV = displayed value Formulas: - W = SV * WS + WI - D = SV * RS + RI - P = D / (RS * SS) + WV = SV * WS + WI + DV = SV * RS + RI + FP = DV / (RS * SS) ``` - ## Derived parametric maps stored with raw diffusion data -Some Philips diffusion DICOM images include derived image(s) along with the images. Other manufacturers save these derived images as a separate series number, and the DICOM standard seems ambiguous on whether it is allowable to mix raw and derived data in the same series (see PS 3.3-2008, C.7.6.1.1.2-3). In practice, many Philips diffusion images append [derived parametric maps](http://www.revisemri.com/blog/2008/diffusion-tensor-imaging/) with the original data. With Philips, appending the derived isotropic image is optional - it is only created for the 'clinical' DTI schemes for radiography analysis and is triggered if the first three vectors in the gradient table are the unit X,Y and Z vectors. For conventional DWI, the result is the conventional mean of the ADC X,Y,Z for DTI it the conventional mean of the 3 principle Eigen vectors. As scientists, we want to discard these derived images, as they will disrupt data processing and we can generate better parametric maps after we have applied undistortion methods such as [Eddy and Topup](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/eddy/UsersGuide). The current version of dcm2niix uses the Diffusion Directionality (0018,9075) tag to detect B=0 unweighted ("NONE"), B-weighted ("DIRECTIONAL"), and derived ("ISOTROPIC") images. Note that the Dimension Index Values (0020,9157) tag provides an alternative approach to discriminate these images. Here are sample tags from a Philips enhanced image that includes and derived map (3rd dimension is "1" while the other images set this to "2"). +Some Philips diffusion DICOM images include derived image(s) along with the images. Other manufacturers save these derived images as a separate series number, and the DICOM standard seems ambiguous on whether it is allowable to mix raw and derived data in the same series (see PS 3.3-2008, C.7.6.1.1.2-3). In practice, many Philips diffusion images append [derived parametric maps](http://www.revisemri.com/blog/2008/diffusion-tensor-imaging/) with the original data. With Philips, appending the derived isotropic image is optional - it is only created for the 'clinical' DTI schemes for radiography analysis and is triggered if the first three vectors in the gradient table are the unit X,Y and Z vectors. For conventional DWI, the result is the conventional mean of the ADC X,Y,Z for DTI it the conventional mean of the 3 principle Eigen vectors. As scientists, we want to discard these derived images, as they will disrupt data processing and we can generate better parametric maps after we have applied undistortion methods such as [Eddy and Topup](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/eddy/UsersGuide). The current version of dcm2niix uses the Diffusion Directionality (0018,9075) tag to detect B=0 unweighted ("NONE"), B-weighted ("DIRECTIONAL"), and derived ("ISOTROPIC") images. Note that the Dimension Index Values (0020,9157) tag provides an alternative approach to discriminate these images. Here are sample tags from a Philips enhanced image that includes and derived map (3rd dimension is "1" while the other images set this to "2"). ``` (0018,9075) CS [DIRECTIONAL] @@ -145,7 +140,7 @@ Another value desirable for TOPUP is the "TotalReadoutTime". Again, one can not ## Partial Volumes -NIfTI expects all 3D volumes of a 4D series to have the same number of series (e.g. a time series of 3D fMRI volumes, or a diffusion set with 3D volumes with different gradients applied). If a fMRI sequence is aborted part way through, it is possible that a Philips scanner will only save part of the final volume. An example would be where the total slices (9970) does not equal Dynamics (290) x slices (35) = 10150. Current versions of dcm2niix expect complete volumes. You can repair your data using the console or a Python script, as discussed in [issue 357](https://github.com/rordenlab/dcm2niix/issues/357). To resolve this situation by hand you could also [rename](RENAMING.md) your DICOM files with a call like `./dcm2niix -r y -f %t/%s_%p_%4y_%2r.dcm ~/out 0020,0100`. In this example, the [`%4y`](FILENAMING.md) parameter adds the volume (Temporal Position, 0020,0100) to the filename, allowing you to identify volumes with missing slices. +NIfTI expects all 3D volumes of a 4D series to have the same number of series (e.g. a time series of 3D fMRI volumes, or a diffusion set with 3D volumes with different gradients applied). If a fMRI sequence is aborted part way through, it is possible that a Philips scanner will only save part of the final volume. An example would be where the total slices (9970) does not equal Dynamics (290) x slices (35) = 10150. Current versions of dcm2niix expect complete volumes. You can repair your data using the console or a Python script, as discussed in [issue 357](https://github.com/rordenlab/dcm2niix/issues/357). To resolve this situation by hand you could also [rename](RENAMING.md) your DICOM files with a call like `./dcm2niix -r y -f %t/%s_%p_%4y_%2r.dcm ~/out 0020,0100`. In this example, the [`%4y`](FILENAMING.md) parameter adds the volume (Temporal Position, 0020,0100) to the filename, allowing you to identify volumes with missing slices. ## Non-Image DICOMs diff --git a/console/main_console.cpp b/console/main_console.cpp index a049830e..b3e05613 100644 --- a/console/main_console.cpp +++ b/console/main_console.cpp @@ -400,9 +400,14 @@ int main(int argc, const char * argv[]) opts.isForceStackSameSeries = 1; if ((argv[i][0] == '2')) opts.isForceStackSameSeries = 2; - if ((argv[i][0] == 'o') || (argv[i][0] == 'O')) + if ((argv[i][0] == 'o') || (argv[i][0] == 'O')) { opts.isForceStackDCE = false; - + //printf("Advanced feature: '-m o' merges images despite varying series number\n"); + } + if ((argv[i][0] == '2')) { + opts.isIgnoreSeriesInstanceUID = true; + printf("Advanced feature: '-m 2' ignores Series Instance UID.\n"); + } } else if ((argv[i][1] == 'p') && ((i+1) < argc)) { i++; if (invalidParam(i, argv)) return 0; diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 6ce08bf9..fa3605bf 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -50,7 +50,7 @@ extern "C" { #define kCPUsuf " " //unknown CPU #endif -#define kDCMdate "v1.0.20210308" +#define kDCMdate "v1.0.20210317" #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 9d99ce3f..32233b78 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -4308,7 +4308,8 @@ float PhilipsPreciseVal (float lPV, float lRS, float lRI, float lSS) { void PhilipsPrecise(struct TDICOMdata * d, bool isPhilipsFloatNotDisplayScaling, struct nifti_1_header *h, int verbose) { if (d->manufacturer != kMANUFACTURER_PHILIPS) return; //not Philips if (d->isScaleVariesEnh) return; //issue363 rescaled before slice reordering - if (!isSameFloatGE(0.0, d->RWVScale)) { + /* + if (!isSameFloatGE(0.0, d->RWVScale)) { //https://github.com/rordenlab/dcm2niix/issues/493 h->scl_slope = d->RWVScale; h->scl_inter = d->RWVIntercept; printMessage("Using RWVSlope:RWVIntercept = %g:%g\n",d->RWVScale,d->RWVIntercept); @@ -4319,7 +4320,7 @@ void PhilipsPrecise(struct TDICOMdata * d, bool isPhilipsFloatNotDisplayScaling, printMessage(" RS = rescale slope, RI = rescale intercept, SS = scale slope\n"); printMessage(" D = R * RS + RI , P = D/(RS * SS)\n"); return; - } + }*/ if (d->intenScalePhilips == 0) return; //no Philips Precise //we will report calibrated "FP" values http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3998685/ float l0 = PhilipsPreciseVal (0, d->intenScale, d->intenIntercept, d->intenScalePhilips); @@ -6696,7 +6697,9 @@ int nii_loadDirCore(char *indir, struct TDCMopts* opts) { continue; } dcmList[i] = readDICOMv(nameList.str[i], opts->isVerbose, opts->compressFlag, dti4D); //ignore compile warning - memory only freed on first of 2 passes - //if (!dcmList[i].isValid) printf(">>>>Not a valid DICOM %s\n", nameList.str[i]); + if (opts->isIgnoreSeriesInstanceUID) + dcmList[i].seriesUidCrc = dcmList[i].seriesNum; + //if (!dcmList[i].isValid) printf(">>>>Not a valid DICOM %s\n", nameList.str[i]); if ((dcmList[i].isValid) && ((dti4D->sliceOrder[0] >= 0) || (dcmList[i].CSA.numDti > 1))) { //4D dataset: dti4D arrays require huge amounts of RAM - write this immediately struct TDCMsort dcmSort[1]; fillTDCMsort(dcmSort[0], i, dcmList[i]); @@ -7235,6 +7238,7 @@ void setDefaultOpts (struct TDCMopts *opts, const char * argv[]) { //either "set opts->isRenameNotConvert = false; opts->isForceStackSameSeries = 2; //automatic: stack CTs, do not stack MRI opts->isForceStackDCE = true; + opts->isIgnoreSeriesInstanceUID = false; opts->isIgnoreDerivedAnd2D = false; opts->isForceOnsetTimes = true; opts->isPhilipsFloatNotDisplayScaling = true; diff --git a/console/nii_dicom_batch.h b/console/nii_dicom_batch.h index e8cad56e..39e66c89 100644 --- a/console/nii_dicom_batch.h +++ b/console/nii_dicom_batch.h @@ -35,7 +35,7 @@ extern "C" { #define MAX_NUM_SERIES 16 struct TDCMopts { - bool isTestx0021x105E, isAddNamePostFixes, isSaveNativeEndian, isSaveNRRD, isOneDirAtATime, isRenameNotConvert, isSave3D, isGz, isPipedGz, isFlipY, isCreateBIDS, isSortDTIbyBVal, isAnonymizeBIDS, isOnlyBIDS, isCreateText, isForceOnsetTimes,isIgnoreDerivedAnd2D, isPhilipsFloatNotDisplayScaling, isTiltCorrect, isRGBplanar, isOnlySingleFile, isForceStackDCE, isRotate3DAcq, isCrop; + bool isTestx0021x105E, isAddNamePostFixes, isSaveNativeEndian, isSaveNRRD, isOneDirAtATime, isRenameNotConvert, isSave3D, isGz, isPipedGz, isFlipY, isCreateBIDS, isSortDTIbyBVal, isAnonymizeBIDS, isOnlyBIDS, isCreateText, isForceOnsetTimes,isIgnoreDerivedAnd2D, isPhilipsFloatNotDisplayScaling, isTiltCorrect, isRGBplanar, isOnlySingleFile, isForceStackDCE, isIgnoreSeriesInstanceUID, isRotate3DAcq, isCrop; int isMaximize16BitRange, isForceStackSameSeries, nameConflictBehavior, isVerbose, isProgress, compressFlag, dirSearchDepth, gzLevel; //support for compressed data 0=none, char filename[512], outdir[512], indir[512], pigzname[512], optsname[512], indirParent[512], imageComments[24]; double seriesNumber[MAX_NUM_SERIES]; //requires double must store -1 (report but do not convert) as well as seriesUidCrc (uint32) From 9a5fa2627bc57180dd34a9261e2e942659a549e7 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Fri, 19 Mar 2021 17:53:34 -0400 Subject: [PATCH 40/41] help should describe accession number (https://github.com/rordenlab/dcm2niix/issues/496) --- FILENAMING.md | 1 + console/main_console.cpp | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/FILENAMING.md b/FILENAMING.md index 8a3a4a31..e1090ab8 100644 --- a/FILENAMING.md +++ b/FILENAMING.md @@ -12,6 +12,7 @@ You request the output file name with the `-f` argument. For example, consider y - %d=description (from 0008,103E) - %e=echo number (from 0018,0086) - %f=folder name (name of folder containing first DICOM) + - %g=accession number (0008,0050) - %i=ID of patient (from 0010,0020) - %j=series instance UID (from 0020,000E) - %k=study instance UID (from 0020,000D) diff --git a/console/main_console.cpp b/console/main_console.cpp index b3e05613..2adc4982 100644 --- a/console/main_console.cpp +++ b/console/main_console.cpp @@ -88,7 +88,7 @@ void showHelp(const char * argv[], struct TDCMopts opts) { #else #define kQstr "" #endif - printf(" -f : filename (%%a=antenna (coil) name, %%b=basename, %%c=comments, %%d=description, %%e=echo number, %%f=folder name, %%i=ID of patient, %%j=seriesInstanceUID, %%k=studyInstanceUID, %%m=manufacturer, %%n=name of patient, %%o=mediaObjectInstanceUID, %%p=protocol,%s %%r=instance number, %%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) name, %%b=basename, %%c=comments, %%d=description, %%e=echo number, %%f=folder name, %%g=accession number, %%i=ID of patient, %%j=seriesInstanceUID, %%k=studyInstanceUID, %%m=manufacturer, %%n=name of patient, %%o=mediaObjectInstanceUID, %%p=protocol,%s %%r=instance number, %%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/i [o=only: reset and write defaults; i=ignore: reset defaults], default n)\n"); printf(" -h : show help\n"); printf(" -i : ignore derived, localizer and 2D images (y/n, default n)\n"); From 7831666a4aa128c41f1fd8fba7ebccb3099e1273 Mon Sep 17 00:00:00 2001 From: neurolabusc Date: Tue, 23 Mar 2021 15:06:51 -0400 Subject: [PATCH 41/41] Remove diagnostic message --- console/nii_dicom.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 0bf88ec0..9194a637 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -6013,7 +6013,6 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); break; case kUnitsPT: //CS dcmStr(lLength, &buffer[lPos], d.unitsPT); - printf(">>>>>>>>%s<<<%s\n", d.unitsPT, fname); break; case kAttenuationCorrectionMethod: //LO dcmStr(lLength, &buffer[lPos], d.attenuationCorrectionMethod);