From d62e4715068ce6965a806bea0c7d87e8fdc4729c Mon Sep 17 00:00:00 2001 From: Chris Rorden Date: Fri, 9 Mar 2018 13:35:04 -0500 Subject: [PATCH 01/14] Fix for https://github.com/rordenlab/dcm2niix/issues/168 --- console/nii_dicom.cpp | 8 ++++++-- console/nii_dicom_batch.cpp | 30 +++++++++++++++++++++++++++++- 2 files changed, 35 insertions(+), 3 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 933782ad..9d8e901b 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -2932,6 +2932,7 @@ void _update_tvd(struct TVolumeDiffusion* ptvd) { struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, struct TDTI4D *dti4D) { struct TDICOMdata d = clear_dicom_data(); + d.imageNum = 0; //not set strcpy(d.protocolName, ""); //erase dummy with empty strcpy(d.protocolName, ""); //erase dummy with empty strcpy(d.seriesDescription, ""); //erase dummy with empty @@ -3574,8 +3575,7 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD d.acquNum = dcmStrInt(lLength, &buffer[lPos]); break; case kImageNum: - //int dx = 3; - if (d.imageNum <= 1) d.imageNum = dcmStrInt(lLength, &buffer[lPos]); //Philips renames each image as image1 in 2001,9000 + 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 kDimensionIndexValues: // kImageNum is not enough for 4D series from Philips 5.*. { // { necessary for initializing ndim. @@ -4169,6 +4169,10 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD printError("Unable to convert DTI [recompile with increased kMaxDTI4D] detected=%d, max = %d\n", d.CSA.numDti, kMaxDTI4D); d.CSA.numDti = 0; } + if (d.imageNum == 0) { + printError("Image number (0020,0013) not found.\n"); + d.imageNum = 1; //not set + } if (multiBandFactor > d.CSA.multiBandFactor) d.CSA.multiBandFactor = multiBandFactor; //SMS reported in 0051,1011 but not CSA header //d.isValid = false; //debug only - will not create output! diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 9e12cc0d..486ce914 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -2605,12 +2605,40 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis void fillTDCMsort(struct TDCMsort& tdcmref, const uint64_t indx, const struct TDICOMdata& dcmdata){ // Copy the relevant parts of dcmdata to tdcmref. tdcmref.indx = indx; + //printf("series/image %d %d\n", dcmdata.seriesNum, dcmdata.imageNum); tdcmref.img = ((uint64_t)dcmdata.seriesNum << 32) + dcmdata.imageNum; for(int i = 0; i < MAX_NUMBER_OF_DIMENSIONS; ++i) tdcmref.dimensionIndexValues[i] = dcmdata.dimensionIndexValues[i]; } // fillTDCMsort() int compareTDCMsort(void const *item1, void const *item2) { + //for quicksort http://blog.ablepear.com/2011/11/objective-c-tuesdays-sorting-arrays.html + struct TDCMsort const *dcm1 = (const struct TDCMsort *)item1; + struct TDCMsort const *dcm2 = (const struct TDCMsort *)item2; + + int retval = 0; // tie + + if (dcm1->img < dcm2->img) + retval = -1; + else if (dcm1->img > dcm2->img) + retval = 1; + if(retval != 0) return retval; //sorted images + // Check the dimensionIndexValues (useful for enhanced DICOM 4D series). + // ->img is basically behaving as a (seriesNum, imageNum) sort key + // concatenated into a (large) integer for qsort. That is unwieldy when + // dimensionIndexValues need to be compared, because the existence of + // uint128_t, uint256_t, etc. is not guaranteed. This sorts by + // (seriesNum, ImageNum, div[0], div[1], ...), or if you think of it as a + // number, the dimensionIndexValues come after the decimal point. + for(int i=0; i < MAX_NUMBER_OF_DIMENSIONS; ++i){ + if(dcm1->dimensionIndexValues[i] < dcm2->dimensionIndexValues[i]) + return -1; + else if(dcm1->dimensionIndexValues[i] > dcm2->dimensionIndexValues[i]) + return 1; + } + return retval; +} //compareTDCMsort() +/*int compareTDCMsort(void const *item1, void const *item2) { //for quicksort http://blog.ablepear.com/2011/11/objective-c-tuesdays-sorting-arrays.html struct TDCMsort const *dcm1 = (const struct TDCMsort *)item1; struct TDCMsort const *dcm2 = (const struct TDCMsort *)item2; @@ -2642,7 +2670,7 @@ int compareTDCMsort(void const *item1, void const *item2) { } } return retval; -} //compareTDCMsort() +} //compareTDCMsort()*/ int isSameFloatGE (float a, float b) { //Kludge for bug in 0002,0016="DIGITAL_JACKET", 0008,0070="GE MEDICAL SYSTEMS" DICOM data: Orient field (0020:0037) can vary 0.00604261 == 0.00604273 !!! From 32623a6236e8cd9a367edc64f032aeff42cd5d0f Mon Sep 17 00:00:00 2001 From: Chris Rorden Date: Fri, 9 Mar 2018 16:42:44 -0500 Subject: [PATCH 02/14] Experimental SQ skipping: Philips private use of public tags https://github.com/rordenlab/dcm2niix/issues/165 --- console/main_console.cpp | 6 +++++- console/nii_dicom.cpp | 36 ++++++++++++++++++++++++------------ 2 files changed, 29 insertions(+), 13 deletions(-) diff --git a/console/main_console.cpp b/console/main_console.cpp index 78cf75ba..19c2a2c5 100644 --- a/console/main_console.cpp +++ b/console/main_console.cpp @@ -324,7 +324,11 @@ int main(int argc, const char * argv[]) if (invalidParam(i, argv)) return 0; if ((argv[i][0] == 'n') || (argv[i][0] == 'N') || (argv[i][0] == '0')) //0: verbose OFF opts.isVerbose = 0; - else if ((argv[i][0] == 'h') || (argv[i][0] == 'H') || (argv[i][0] == '2')) //2: verbose HYPER + else if ((argv[i][0] == 'o') || (argv[i][0] == 'O')) {//-1: experimental SQ skipping: 'O'ptimized!? faster, does not get confused with Philips tags + opts.isVerbose = -1; + printf(">>Experimental SQ skipping enabled\n"); + + } else if ((argv[i][0] == 'h') || (argv[i][0] == 'H') || (argv[i][0] == '2')) //2: verbose HYPER opts.isVerbose = 2; else opts.isVerbose = 1; //1: verbose ON diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 9d8e901b..cf1db775 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -3183,8 +3183,7 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD float patientPosition[4] = {NAN, NAN, NAN, NAN}; //used to compute slice direction for Philips 4D float patientPositionEndPhilips[4] = {NAN, NAN, NAN, NAN}; float patientPositionStartPhilips[4] = {NAN, NAN, NAN, NAN}; - - while ((d.imageStart == 0) && ((lPos+8+lFileOffset) < fileLen)) { + while ((d.imageStart == 0) && ((lPos+8+lFileOffset) < fileLen)) { #ifndef myLoadWholeFileToReadHeader //read one segment at a time if ((size_t)(lPos + 128) > MaxBufferSz) { //avoid overreading the file lFileOffset = lFileOffset + lPos; @@ -3220,13 +3219,13 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD vr[0] = 'N'; vr[1] = 'A'; if (groupElement == kUnnest2) sqDepth--; + //if (sqDepth < 0) sqDepth = 0; //if (groupElement == kUnnest) geiisBug = false; //don't exit if there is a proprietary thumbnail lLength = 4; //Next: Skip nested items if explicit length provided - //int nestLength = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); - //if (nestLength < 0xFFFFFFFF) { - // printMessage(" >>>>> %d\n", nestLength); - //} + // int nestLength = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); + //if (isVerbose > 1) + // printMessage(" ---SQ\t%04x,%04x %d\n", groupElement & 65535,groupElement>>16, sqDepth); } else if (d.isExplicitVR) { vr[0] = buffer[lPos]; vr[1] = buffer[lPos+1]; if (buffer[lPos+1] < 'A') {//implicit vr with 32-bit length @@ -3275,17 +3274,27 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD } //if explicit else implicit VR if (lLength == 0xFFFFFFFF) { lLength = 8; //SQ (Sequences) use 0xFFFFFFFF [4294967295] to denote unknown length - vr[0] = 'S'; - vr[1] = 'Q'; + //09032018 - do not count these as SQs: Horos does not count even groups + uint32_t special = dcmInt(4,&buffer[lPos],d.isLittleEndian); + //http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_7.5.html + //uint32_t kItem = 0xFFFE + (0xE000 << 16 ); //#define kNest 0xFFFE +(0xE000 << 16 ) //Item follows SQ + //uint32_t kItemDelim = 0xFFFE + (0xE00D << 16); + uint32_t ksqDelim = 0xFFFE + (0xE0DD << 16); + if (special != ksqDelim) { + vr[0] = 'S'; + vr[1] = 'Q'; + } } if ((vr[0] == 'S') && (vr[1] == 'Q')) { sqDepth++; - //printMessage("SQstart %d\n", sqDepth); + //if (isVerbose > 1) + // printMessage(" +++SQ\t%04x,%04x %d\n", groupElement & 65535,groupElement>>16, sqDepth); } if ((groupElement == kNest) || ((vr[0] == 'S') && (vr[1] == 'Q'))) nest++; //if ((groupElement == kUnnest) || (groupElement == kUnnest2)) { // <- ? if (groupElement == kUnnest) { nest--; + if (nest < 0) nest = 0; } //next: look for required tags if ((groupElement == kNest) && (isEncapsulatedData)) { @@ -3298,7 +3307,10 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD } } if ((isIconImageSequence) && ((groupElement & 0x0028) == 0x0028 )) groupElement = kUnused; //ignore icon dimensions - if (true) { //(nest <= 0) { //some Philips images have different 0020,0013 + + if ((isVerbose >= 0) || (nest <= 0)) { //(nest <= 0) { //verbose=-1: turbo mode skips nested tags Philips images have different 0020,0013 + + //xx if (true) { //(nest <= 0) { //some Philips images have different 0020,0013 //verbose reporting : // printMessage("Pos %ld GroupElement %#08x,%#08x Length %d isLittle %d\n", lPos, (groupElement & 0xFFFF), (groupElement >> 16), lLength, d.isLittleEndian); @@ -4048,9 +4060,9 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD tagStr[pos] = 'x'; } printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\t%s\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos, tagStr); - //printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tnest=%d\t%s\n", groupElement & 65535,groupElement>>16, lLength, lPos, nest, tagStr); + //printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tnest=%d\tsq=%d\t%s\n", groupElement & 65535,groupElement>>16, lLength, lPos, nest, tagStr); } else - printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tnest=%d\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos, nest); + printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tnest=%d\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos, nest, sqDepth); //if (d.isExplicitVR) printMessage(" VR=%c%c\n", vr[0], vr[1]); } //printMessage(" tag=%04x,%04x length=%u pos=%ld %c%c nest=%d\n", groupElement & 65535,groupElement>>16, lLength, lPos,vr[0], vr[1], nest); lPos = lPos + (lLength); From 7aa3d18d48c0edcd6119638128c29f7760c9fcef Mon Sep 17 00:00:00 2001 From: Chris Rorden Date: Sun, 11 Mar 2018 22:42:53 -0400 Subject: [PATCH 03/14] Partial solution for https://github.com/rordenlab/dcm2niix/issues/165 Philips saves 4D datasets where the slice order of the first volume does not match subsequent volumes. Slices need to be re-ordered spatially. --- console/main_console.cpp | 2 +- console/nii_dicom.cpp | 140 +++++++++++++++++++++------------------ console/nii_dicom.h | 1 + 3 files changed, 79 insertions(+), 64 deletions(-) diff --git a/console/main_console.cpp b/console/main_console.cpp index 19c2a2c5..dc976cad 100644 --- a/console/main_console.cpp +++ b/console/main_console.cpp @@ -326,7 +326,7 @@ int main(int argc, const char * argv[]) opts.isVerbose = 0; else if ((argv[i][0] == 'o') || (argv[i][0] == 'O')) {//-1: experimental SQ skipping: 'O'ptimized!? faster, does not get confused with Philips tags opts.isVerbose = -1; - printf(">>Experimental SQ skipping enabled\n"); + printf(">> Experimental SQ skipping enabled\n"); } else if ((argv[i][0] == 'h') || (argv[i][0] == 'H') || (argv[i][0] == '2')) //2: verbose HYPER opts.isVerbose = 2; diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index cf1db775..c819baa9 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -1447,6 +1447,8 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D bool isIndexSequential = true; for (int i = 0; i < kMaxDTI4D; i++) dti4D->S[i].V[0] = -1.0; + dti4D->S[0].sliceNumberMrPhilips = -1.0; + dti4D->S[0].sliceNumberMrPhilipsVol2 = -1.0; //d.dti4D = (TDTI *)malloc(kMaxDTI4D * sizeof(TDTI)); while (p) { if (strlen(buff) < 1) @@ -1924,11 +1926,16 @@ unsigned char * nii_reorderSlices(unsigned char* bImg, struct nifti_1_header *h, if (sliceBytes < 0) return bImg; size_t volBytes = sliceBytes * h->dim[3]; unsigned char *srcImg = (unsigned char *)malloc(volBytes); + //printMessage("Reordering %d volumes\n", dim4to7); for (int v = 0; v < dim4to7; v++) { + //for (int v = 0; v < 1; v++) { + size_t volStart = v * volBytes; memcpy(&srcImg[0], &bImg[volStart], volBytes); //dest, src, size for (int z = 0; z < h->dim[3]; z++) { //for each slice int src = dti4D->S[z].sliceNumberMrPhilips - 1; //-1 as Philips indexes slices from 1 not 0 + if ((v > 0) && (dti4D->S[0].sliceNumberMrPhilipsVol2 >= 0)) + src = dti4D->S[z].sliceNumberMrPhilipsVol2 - 1; //printMessage("Reordering volume %d slice %d\n", v, dti4D->S[z].sliceNumberMrPhilips); if ((src < 0) || (src >= h->dim[3])) continue; memcpy(&bImg[volStart+(src*sliceBytes)], &srcImg[z*sliceBytes], sliceBytes); //dest, src, size @@ -2938,6 +2945,8 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru strcpy(d.seriesDescription, ""); //erase dummy with empty strcpy(d.sequenceName, ""); //erase dummy with empty //do not read folders - code specific to GCC (LLVM/Clang seems to recognize a small file size) + dti4D->S[0].sliceNumberMrPhilips = -1; + dti4D->S[0].sliceNumberMrPhilipsVol2 = -1; struct TVolumeDiffusion volDiffusion = initTVolumeDiffusion(&d, dti4D); @@ -3144,14 +3153,13 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kImageStart 0x7FE0+(0x0010 << 16 ) #define kImageStartFloat 0x7FE0+(0x0008 << 16 ) #define kImageStartDouble 0x7FE0+(0x0009 << 16 ) -uint32_t kNest = 0xFFFE +(0xE000 << 16 ); //#define kNest 0xFFFE +(0xE000 << 16 ) //Item follows SQ -uint32_t kUnnest = 0xFFFE +(0xE00D << 16 ); //#define kUnnest 0xFFFE +(0xE00D << 16 ) //ItemDelimitationItem [length defined] http://www.dabsoft.ch/dicom/5/7.5/ -uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD << 16 )//SequenceDelimitationItem [length undefined] - int nest = 0; - //double zSpacing = -1.0l; //includes slice thickness plus gap +uint32_t kItemTag = 0xFFFE +(0xE000 << 16 ); +uint32_t kItemDelimitationTag = 0xFFFE +(0xE00D << 16 ); +uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); int locationsInAcquisitionGE = 0; int locationsInAcquisitionPhilips = 0; int imagesInAcquisition = 0; + int sliceNumberMrPhilips = 0; uint32_t lLength; uint32_t groupElement; long lPos = 0; @@ -3180,6 +3188,8 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD bool isMosaic = false; int patientPositionNum = 0; int sqDepth = 0; + int sqDepthPrivate = 0; + int sqEndPrivate = -1; //used to skip private SQs that provide explicit length float patientPosition[4] = {NAN, NAN, NAN, NAN}; //used to compute slice direction for Philips 4D float patientPositionEndPhilips[4] = {NAN, NAN, NAN, NAN}; float patientPositionStartPhilips[4] = {NAN, NAN, NAN, NAN}; @@ -3213,19 +3223,11 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD } //transfer syntax requests switching VR after group 0001 //uint32_t group = (groupElement & 0xFFFF); lPos += 4; - if ((groupElement == kUnnest) || (groupElement == kUnnest2)) isIconImageSequence = false; - if (((groupElement == kNest) || (groupElement == kUnnest) || (groupElement == kUnnest2)) && (!isEncapsulatedData)) { - //if (((groupElement == kNest) || (groupElement == kUnnest) || (groupElement == kUnnest2)) ) { + if ((groupElement == kItemDelimitationTag) || (groupElement == kSequenceDelimitationItemTag)) isIconImageSequence = false; + if (((groupElement == kItemTag) || (groupElement == kItemDelimitationTag) || (groupElement == kSequenceDelimitationItemTag)) && (!isEncapsulatedData)) { vr[0] = 'N'; vr[1] = 'A'; - if (groupElement == kUnnest2) sqDepth--; - //if (sqDepth < 0) sqDepth = 0; - //if (groupElement == kUnnest) geiisBug = false; //don't exit if there is a proprietary thumbnail lLength = 4; - //Next: Skip nested items if explicit length provided - // int nestLength = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); - //if (isVerbose > 1) - // printMessage(" ---SQ\t%04x,%04x %d\n", groupElement & 65535,groupElement>>16, sqDepth); } else if (d.isExplicitVR) { vr[0] = buffer[lPos]; vr[1] = buffer[lPos+1]; if (buffer[lPos+1] < 'A') {//implicit vr with 32-bit length @@ -3275,29 +3277,34 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD if (lLength == 0xFFFFFFFF) { lLength = 8; //SQ (Sequences) use 0xFFFFFFFF [4294967295] to denote unknown length //09032018 - do not count these as SQs: Horos does not count even groups - uint32_t special = dcmInt(4,&buffer[lPos],d.isLittleEndian); + //uint32_t special = dcmInt(4,&buffer[lPos],d.isLittleEndian); + //http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_7.5.html - //uint32_t kItem = 0xFFFE + (0xE000 << 16 ); //#define kNest 0xFFFE +(0xE000 << 16 ) //Item follows SQ - //uint32_t kItemDelim = 0xFFFE + (0xE00D << 16); - uint32_t ksqDelim = 0xFFFE + (0xE0DD << 16); - if (special != ksqDelim) { + + //if (special != ksqDelim) { vr[0] = 'S'; vr[1] = 'Q'; - } + //} } if ((vr[0] == 'S') && (vr[1] == 'Q')) { - sqDepth++; - //if (isVerbose > 1) - // printMessage(" +++SQ\t%04x,%04x %d\n", groupElement & 65535,groupElement>>16, sqDepth); - } - if ((groupElement == kNest) || ((vr[0] == 'S') && (vr[1] == 'Q'))) nest++; - //if ((groupElement == kUnnest) || (groupElement == kUnnest2)) { // <- ? - if (groupElement == kUnnest) { - nest--; - if (nest < 0) nest = 0; - } + //http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_7.5.html + uint32_t special = dcmInt(4,&buffer[lPos],d.isLittleEndian); + uint32_t slen = dcmInt(4,&buffer[lPos+4],d.isLittleEndian); + uint32_t kUndefinedLen = 0xFFFFFFFF; + if (special == kSequenceDelimitationItemTag) { + //sqDepth--; + //printMessage(" SPECIAL >>>>t%04x,%04x %08x %d\n", groupElement & 65535,groupElement>>16, special, slen); + } else if (slen == kUndefinedLen) { + sqDepth++; + if ((sqDepthPrivate == 0) && ((groupElement & 65535) % 2)) + sqDepthPrivate = sqDepth; //in a private SQ: ignore contents + } else if ((groupElement & 65535) % 2) //private SQ of known length - lets jump over this! + slen = lFileOffset + slen; + if (slen > sqEndPrivate) + sqEndPrivate = slen; //if nested private SQs, remember the end address of the top parent SQ + } //next: look for required tags - if ((groupElement == kNest) && (isEncapsulatedData)) { + if ((groupElement == kItemTag) && (isEncapsulatedData)) { d.imageBytes = dcmInt(4,&buffer[lPos-4],d.isLittleEndian); //printMessage("compressed data %d-> %ld\n",d.imageBytes, lPos); if (d.imageBytes > 128) { @@ -3307,23 +3314,18 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD } } if ((isIconImageSequence) && ((groupElement & 0x0028) == 0x0028 )) groupElement = kUnused; //ignore icon dimensions - - if ((isVerbose >= 0) || (nest <= 0)) { //(nest <= 0) { //verbose=-1: turbo mode skips nested tags Philips images have different 0020,0013 - - //xx if (true) { //(nest <= 0) { //some Philips images have different 0020,0013 - //verbose reporting : - // printMessage("Pos %ld GroupElement %#08x,%#08x Length %d isLittle %d\n", lPos, (groupElement & 0xFFFF), (groupElement >> 16), lLength, d.isLittleEndian); - - // // Debugging - // int groupItem = groupElement >> 16; - // int groupGroup = groupElement - (groupItem << 16); - // printMessage("groupElement: (%04X, %04X)\n", groupElement & 65535,groupElement>>16); - + if (groupElement == kSequenceDelimitationItemTag) { + sqDepth--; + if (sqDepth < sqDepthPrivate) { + sqDepthPrivate = 0; //no longer in a private SQ + is2005140FSQ = false; + } + } + if (sqDepth < 0) sqDepth = 0; switch ( groupElement ) { case kTransferSyntax: { char transferSyntax[kDICOMStr]; dcmStr (lLength, &buffer[lPos], transferSyntax); - //printMessage("%d transfer syntax>>> '%s'\n", compressFlag, transferSyntax); if (strcmp(transferSyntax, "1.2.840.10008.1.2.1") == 0) ; //default isExplicitVR=true; //d.isLittleEndian=true else if (strcmp(transferSyntax, "1.2.840.10008.1.2.4.50") == 0) { @@ -3539,8 +3541,14 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD case kSeriesInstanceUID : // 0020, 000E dcmStr (lLength, &buffer[lPos], d.seriesInstanceUID); break; - case kPatientPosition : // 0020, 0032, ImagePositionPatient - if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (is2005140FSQ)) { + case kPatientPosition : { // 0020,0032, ImagePositionPatient + bool ignorePatientPos = false; + if (d.manufacturer == kMANUFACTURER_PHILIPS) { + if (is2005140FSQ) ignorePatientPos = true; + if ((lFileOffset + lPos) < sqEndPrivate) ignorePatientPos = true; //inside private SQ, SQ has defined length + if (sqDepthPrivate > 0) ignorePatientPos = true; //inside private SQ, SQ has undefined length + } + if (ignorePatientPos) { if (!is2005140FSQwarned) printWarning("Philips R3.2.2 can report different positions for the same slice. Attempting patch.\n"); is2005140FSQwarned = true; @@ -3567,9 +3575,8 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD set_isAtFirstPatientPosition_tvd(&volDiffusion, isAtFirstPatientPosition); if (isVerbose > 1) printMessage(" Patient Position 0020,0032 (#,@,X,Y,Z)\t%d\t%ld\t%g\t%g\t%g\n", patientPositionNum, lPos, patientPosition[1], patientPosition[2], patientPosition[3]); - } //not after 2005,140F - break; + break; } case kInPlanePhaseEncodingDirection: d.phaseEncodingRC = toupper(buffer[lPos]); //first character is either 'R'ow or 'C'ol break; @@ -3887,17 +3894,23 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD // // part of buffer[lPos]. // } // break; - case kSliceNumberMrPhilips : + case kSliceNumberMrPhilips : { if (d.manufacturer != kMANUFACTURER_PHILIPS) break; - /*if ((d.patientPositionNumPhilips >= kMaxDTI4D) || (d.patientPositionNumPhilips >= locationsInAcquisitionPhilips)) { - d.patientPositionNumPhilips++; - break; - }*/ - if ((locationsInAcquisitionPhilips > 0) && (d.patientPositionNumPhilips == locationsInAcquisitionPhilips)) - break; //we have acquired all slices in volume (e.g. all volumes after 1st for XYZT storage int sliceNumber; sliceNumber = dcmStrInt(lLength, &buffer[lPos]); + //printf("====> %d %d %d\n", locationsInAcquisitionPhilips, d.patientPositionNumPhilips, sliceNumber); + if ((sliceNumberMrPhilips >= locationsInAcquisitionPhilips) && (sliceNumberMrPhilips < (2*locationsInAcquisitionPhilips)) ) { + int p = sliceNumberMrPhilips - locationsInAcquisitionPhilips; + if (p < kMaxDTI4D) + dti4D->S[p].sliceNumberMrPhilipsVol2 = sliceNumber; + if ((p > 0) && (abs(dti4D->S[p].sliceNumberMrPhilipsVol2 - dti4D->S[p -1].sliceNumberMrPhilipsVol2) > 1) ) + d.isSlicesSpatiallySequentialPhilips = false; + } + sliceNumberMrPhilips ++; + if ((locationsInAcquisitionPhilips > 0) && (d.patientPositionNumPhilips == locationsInAcquisitionPhilips)) + break; //we have acquired all slices in volume (e.g. all volumes after 1st for XYZT storage + if ((d.patientPositionNumPhilips > 0) && (sliceNumber == dti4D->S[d.patientPositionNumPhilips-1].sliceNumberMrPhilips) ) break; //repeated spatial position (e.g. data saved XYTZ so several time points if (d.patientPositionNumPhilips >= kMaxDTI4D) { @@ -3919,7 +3932,7 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD } if (isVerbose > 1) printMessage("slice %d is spatial position %d\n", d.patientPositionNumPhilips, sliceNumber); - break; + break; } case kNumberOfSlicesMrPhilips : if (d.manufacturer != kMANUFACTURER_PHILIPS) break; @@ -4045,7 +4058,7 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD isIconImageSequence = false; break; } //switch/case for groupElement - } //if nest + if (isVerbose > 1) { if ((lLength > 12) && (lLength < 128)) { //if length is greater than 8 bytes (+4 hdr) the data must be a string [or image data] char tagStr[kDICOMStr]; @@ -4060,9 +4073,8 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD tagStr[pos] = 'x'; } printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\t%s\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos, tagStr); - //printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tnest=%d\tsq=%d\t%s\n", groupElement & 65535,groupElement>>16, lLength, lPos, nest, tagStr); } else - printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tnest=%d\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos, nest, sqDepth); + printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tsq=%d %d\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos, sqDepth, sqDepthPrivate); //if (d.isExplicitVR) printMessage(" VR=%c%c\n", vr[0], vr[1]); } //printMessage(" tag=%04x,%04x length=%u pos=%ld %c%c nest=%d\n", groupElement & 65535,groupElement>>16, lLength, lPos,vr[0], vr[1], nest); lPos = lPos + (lLength); @@ -4158,13 +4170,15 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD // else // printWarning("Spatial orientation ambiguous (tag 0020,0037 not found): %s\n", fname); // } - if ((d.numberOfDynamicScans < 2) && (!d.isSlicesSpatiallySequentialPhilips) && (!isnan(patientPositionStartPhilips[1])) && (!isnan(patientPositionEndPhilips[1]))) { - //to do: check for d.numberOfDynamicScans > 1 + if ( (!d.isSlicesSpatiallySequentialPhilips) && (!isnan(patientPositionStartPhilips[1])) && (!isnan(patientPositionEndPhilips[1]))) { for (int k = 0; k < 4; k++) { d.patientPosition[k] = patientPositionStartPhilips[k]; d.patientPositionLast[k] = patientPositionEndPhilips[k]; } - printMessage("Slices not spatially contiguous: please check output [new feature]\n"); + if (d.numberOfDynamicScans > 1) + printError("4D data with non-contiguous slices: please check output [newest feature]\n"); + else + printMessage("Slices not spatially contiguous: please check output [new feature]\n"); } if (isVerbose) { printMessage("%s\n patient position (0020,0032)\t%g\t%g\t%g\n",fname, d.patientPosition[1],d.patientPosition[2],d.patientPosition[3]); diff --git a/console/nii_dicom.h b/console/nii_dicom.h index fb513e0a..0d0c22d4 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -77,6 +77,7 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; struct TDTI { float V[4]; int sliceNumberMrPhilips; + int sliceNumberMrPhilipsVol2; }; struct TDTI4D { struct TDTI S[kMaxDTI4D]; From a1d2f67c3e13104dc4327f6e6f0ef550f5bd5650 Mon Sep 17 00:00:00 2001 From: Chris Rorden Date: Mon, 12 Mar 2018 13:08:39 -0400 Subject: [PATCH 04/14] More Philips enhanced support --- console/nii_dicom.cpp | 92 ++++++++++++++++++++++++++++++++----- console/nii_dicom.h | 7 +-- console/nii_dicom_batch.cpp | 26 ++++++----- 3 files changed, 99 insertions(+), 26 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index c819baa9..49635e92 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -693,6 +693,7 @@ struct TDICOMdata clear_dicom_data() { strcpy(d.bodyPartExamined,""); d.phaseEncodingLines = 0; d.patientPositionSequentialRepeats = 0; + d.patientPositionRepeats = 0; d.isHasPhase = false; d.isHasMagnitude = false; d.sliceOrient = kSliceOrientUnknown; @@ -704,6 +705,7 @@ struct TDICOMdata clear_dicom_data() { d.lastScanLoc = NAN; d.TR = 0.0; d.TE = 0.0; + d.TE2 = 0.0; d.TI = 0.0; d.flipAngle = 0.0; d.bandwidthPerPixelPhaseEncode = 0.0; @@ -750,6 +752,7 @@ struct TDICOMdata clear_dicom_data() { d.isSigned = false; //default is unsigned! d.isFloat = false; //default is for integers, not single or double precision d.isResampled = false; //assume data not resliced to remove gantry tilt problems + d.isLocalizer = false; d.compressionScheme = 0; //none d.isExplicitVR = true; d.isLittleEndian = true; //DICOM initially always little endian @@ -2710,7 +2713,13 @@ unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct hdr->dim[0] = 4; } if ((hdr->dim[0] > 3) && (dcm.patientPositionSequentialRepeats > 1)) //swizzle 3rd and 4th dimension (Philips stores time as 3rd dimension) - img = nii_XYTZ_XYZT(img, hdr,dcm.patientPositionSequentialRepeats ); + img = nii_XYTZ_XYZT(img, hdr,dcm.patientPositionSequentialRepeats); + if (((dcm.patientPositionSequentialRepeats * 2) == dcm.patientPositionRepeats) && (dcm.isHasPhase) && (dcm.isHasMagnitude)) { + hdr->dim[3] = hdr->dim[3] / 2; + hdr->dim[4] = hdr->dim[4] * 2; + hdr->dim[0] = 4; + printMessage("Splitting Phase+Magnitude into two volumes for %d slices (Z) and %d volumes (T).\n",hdr->dim[3], hdr->dim[4]); + } headerDcm2NiiSForm(dcm,dcm, hdr, false); return img; } //nii_loadImgXL() @@ -2726,6 +2735,12 @@ int isSQ(uint32_t groupElement) { //Detect sequence VR ("SQ") for implicit tags return 0; } //isSQ() +int isSameFloatGE (float a, float b) { +//Kludge for bug in 0002,0016="DIGITAL_JACKET", 0008,0070="GE MEDICAL SYSTEMS" DICOM data: Orient field (0020:0037) can vary 0.00604261 == 0.00604273 !!! + //return (a == b); //niave approach does not have any tolerance for rounding errors + return (fabs (a - b) <= 0.0001); +} + int isDICOMfile(const char * fname) { //0=NotDICOM, 1=DICOM, 2=Maybe(not Part 10 compliant) FILE *fp = fopen(fname, "rb"); if (!fp) return 0; @@ -3053,6 +3068,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kZThick 0x0018+(0x0050 << 16 ) #define kTR 0x0018+(0x0080 << 16 ) #define kTE 0x0018+(0x0081 << 16 ) +#define kEffectiveTE 0x0018+(0x9082 << 16 ) #define kTI 0x0018+(0x0082 << 16) // Inversion time #define kEchoNum 0x0018+(0x0086 << 16 ) //IS #define kMagneticFieldStrength 0x0018+(0x0087 << 16 ) //DS @@ -3290,7 +3306,11 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); //http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_7.5.html uint32_t special = dcmInt(4,&buffer[lPos],d.isLittleEndian); uint32_t slen = dcmInt(4,&buffer[lPos+4],d.isLittleEndian); + //if (d.isExplicitVR) + // slen = dcmInt(4,&buffer[lPos+8],d.isLittleEndian); uint32_t kUndefinedLen = 0xFFFFFFFF; + //printError(" SPECIAL >>>>t%04x,%04x %08x %08x\n", groupElement & 65535,groupElement>>16, special, slen); + //return d; if (special == kSequenceDelimitationItemTag) { //sqDepth--; //printMessage(" SPECIAL >>>>t%04x,%04x %08x %d\n", groupElement & 65535,groupElement>>16, special, slen); @@ -3435,8 +3455,11 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); break; case kComplexImageComponent: if (lLength < 2) break; - d.isHasPhase = (buffer[lPos]=='P') && (toupper(buffer[lPos+1]) == 'H'); - d.isHasMagnitude = (buffer[lPos]=='M') && (toupper(buffer[lPos+1]) == 'A'); + //not: it is possible for Philips enhanced DICOM to store BOTH magnitude and phase in the same image + if ((buffer[lPos]=='P') && (toupper(buffer[lPos+1]) == 'H')) + d.isHasPhase = true; + if ((buffer[lPos]=='M') && (toupper(buffer[lPos+1]) == 'A')) + d.isHasMagnitude = true; break; case kAcquisitionTime : char acquisitionTimeTxt[kDICOMStr]; @@ -3572,9 +3595,10 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); d.patientPositionSequentialRepeats = patientPositionNum-1; } //if different position from 1st slice in file } //if not first slice in file + if (isAtFirstPatientPosition) d.patientPositionRepeats ++; set_isAtFirstPatientPosition_tvd(&volDiffusion, isAtFirstPatientPosition); if (isVerbose > 1) - printMessage(" Patient Position 0020,0032 (#,@,X,Y,Z)\t%d\t%ld\t%g\t%g\t%g\n", patientPositionNum, lPos, patientPosition[1], patientPosition[2], patientPosition[3]); + printMessage(" Patient Position 0020,0032 (#,#repeat,@,X,Y,Z)\t%d\t%ld\t%g\t%g\t%g\n", patientPositionNum, d.patientPositionSequentialRepeats, lPos, patientPosition[1], patientPosition[2], patientPosition[3]); } //not after 2005,140F break; } case kInPlanePhaseEncodingDirection: @@ -3656,9 +3680,20 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); case kTR : d.TR = dcmStrFloat(lLength, &buffer[lPos]); break; - case kTE : - d.TE = dcmStrFloat(lLength, &buffer[lPos]); - break; + case kTE : { + double te = dcmStrFloat(lLength, &buffer[lPos]); + if (d.TE <= 0.0) + d.TE = te; + else if (d.TE != te) + d.TE2 = te; + break; } + case kEffectiveTE : { + double te = dcmFloatDouble(lLength, &buffer[lPos], d.isLittleEndian); + if (d.TE <= 0.0) + d.TE = te; + else if (d.TE != te) + d.TE2 = te; + break; } case kTI : d.TI = dcmStrFloat(lLength, &buffer[lPos]); break; @@ -3681,7 +3716,15 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); d.phaseFieldofView = dcmStrFloat(lLength, &buffer[lPos]); break; case kPixelBandwidth : + /*if (d.manufacturer == kMANUFACTURER_PHILIPS) { + //Private SQs can report different (more precise?) pixel bandwidth values than in the main header! + // https://github.com/rordenlab/dcm2niix/issues/170 + if (is2005140FSQ) break; + if ((lFileOffset + lPos) < sqEndPrivate) break; //inside private SQ, SQ has defined length + if (sqDepthPrivate > 0) break; //inside private SQ, SQ has undefined length + }*/ d.pixelBandwidth = dcmStrFloat(lLength, &buffer[lPos]); + //printWarning(" PixelBandwidth (0018,0095)====> %g @%d\n", d.pixelBandwidth, lPos); break; case kAcquisitionMatrix : if (lLength == 8) { @@ -3899,7 +3942,7 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); break; int sliceNumber; sliceNumber = dcmStrInt(lLength, &buffer[lPos]); - //printf("====> %d %d %d\n", locationsInAcquisitionPhilips, d.patientPositionNumPhilips, sliceNumber); + printError("====> %d %d %d\n", locationsInAcquisitionPhilips, d.patientPositionNumPhilips, sliceNumber); if ((sliceNumberMrPhilips >= locationsInAcquisitionPhilips) && (sliceNumberMrPhilips < (2*locationsInAcquisitionPhilips)) ) { int p = sliceNumberMrPhilips - locationsInAcquisitionPhilips; if (p < kMaxDTI4D) @@ -4013,10 +4056,24 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); case kOrientationACR : //use in emergency if kOrientation is not present! if (!isOrient) dcmMultiFloat(lLength, (char*)&buffer[lPos], 6, d.orient); break; - case kOrientation : + case kOrientation : { + if (isOrient) { //already read orient - read for this slice to see if it varies (localizer) + float orient[7]; + dcmMultiFloat(lLength, (char*)&buffer[lPos], 6, orient); + if ((!isSameFloatGE(d.orient[1], orient[1]) || !isSameFloatGE(d.orient[2], orient[2]) || !isSameFloatGE(d.orient[3], orient[3]) || + !isSameFloatGE(d.orient[4], orient[4]) || !isSameFloatGE(d.orient[5], orient[5]) || !isSameFloatGE(d.orient[6], orient[6]) ) ) { + d.isLocalizer = true; + if ((isVerbose > 1)) + printMessage("slices not stacked: orientation varies (localizer?) [%g %g %g %g %g %g] != [%g %g %g %g %g %g]\n", + d.orient[1], d.orient[2], d.orient[3],d.orient[4], d.orient[5], d.orient[6], + orient[1], orient[2], orient[3],orient[4], orient[5], orient[6]); + } + + } dcmMultiFloat(lLength, (char*)&buffer[lPos], 6, d.orient); + isOrient = true; - break; + break; } case kImagesInAcquisition : imagesInAcquisition = dcmStrInt(lLength, &buffer[lPos]); break; @@ -4074,7 +4131,7 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); } printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\t%s\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos, tagStr); } else - printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tsq=%d %d\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos, sqDepth, sqDepthPrivate); + printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tsq=%d %d %d %c%c\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos, sqDepth, sqDepthPrivate, sqEndPrivate, vr[0], vr[1]); //if (d.isExplicitVR) printMessage(" VR=%c%c\n", vr[0], vr[1]); } //printMessage(" tag=%04x,%04x length=%u pos=%ld %c%c nest=%d\n", groupElement & 65535,groupElement>>16, lLength, lPos,vr[0], vr[1], nest); lPos = lPos + (lLength); @@ -4170,6 +4227,19 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); // else // printWarning("Spatial orientation ambiguous (tag 0020,0037 not found): %s\n", fname); // } +/*if (d.isHasMagnitude) + printError("=====> mag %d %d\n", d.patientPositionRepeats, d.patientPositionSequentialRepeats ); +if (d.isHasPhase) + printError("=====> phase %d %d\n", d.patientPositionRepeats, d.patientPositionSequentialRepeats ); + + printError("=====> reps %d %d\n", d.patientPositionRepeats, d.patientPositionSequentialRepeats ); +*/ + if ((d.patientPositionSequentialRepeats > 1) && ( (d.xyzDim[3] % d.patientPositionSequentialRepeats) == 0 )) { + //will require Converting XYTZ to XYZT + d.numberOfDynamicScans = d.xyzDim[3] / d.patientPositionSequentialRepeats; + d.xyzDim[4] = d.xyzDim[3] / d.numberOfDynamicScans; + d.xyzDim[3] = d.numberOfDynamicScans; + } if ( (!d.isSlicesSpatiallySequentialPhilips) && (!isnan(patientPositionStartPhilips[1])) && (!isnan(patientPositionEndPhilips[1]))) { for (int k = 0; k < 4; k++) { d.patientPosition[k] = patientPositionStartPhilips[k]; diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 0d0c22d4..08ef2a81 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -118,8 +118,8 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; struct TDICOMdata { long seriesNum; int xyzDim[5]; - int modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, echoTrainLength, patientPositionNumPhilips, coilNum, echoNum, sliceOrient,numberOfDynamicScans, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,patientPositionSequentialRepeats,locationsInAcquisition, compressionScheme; - float patientWeight, zSpacing, zThick, pixelBandwidth, SAR, phaseFieldofView, accelFactPE, flipAngle, fieldStrength, TE, TI, TR, intenScale, intenIntercept, intenScalePhilips, gantryTilt, lastScanLoc, angulation[4]; + int modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, echoTrainLength, patientPositionNumPhilips, coilNum, echoNum, sliceOrient,numberOfDynamicScans, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,patientPositionSequentialRepeats,patientPositionRepeats,locationsInAcquisition, compressionScheme; + float patientWeight, zSpacing, zThick, pixelBandwidth, SAR, phaseFieldofView, accelFactPE, flipAngle, fieldStrength, TE, TE2, TI, TR, intenScale, intenIntercept, intenScalePhilips, gantryTilt, lastScanLoc, angulation[4]; float orient[7], patientPosition[4], patientPositionLast[4], xyzMM[4], stackOffcentre[4]; float radionuclidePositronFraction, radionuclideTotalDose, radionuclideHalfLife, doseCalibrationFactor; //PET ISOTOPE MODULE ATTRIBUTES (C.8-57) float ecat_isotope_halflife, ecat_dosage; @@ -129,11 +129,12 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; char imageComments[kDICOMStrLarge]; uint32_t dimensionIndexValues[MAX_NUMBER_OF_DIMENSIONS]; struct TCSAdata CSA; - bool isSegamiOasis, isDerived, isXRay, isMultiEcho, isSlicesSpatiallySequentialPhilips, isValid, is3DAcq, is2DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase,isHasMagnitude,isHasMixed, isFloat, isResampled; + bool isSegamiOasis, isDerived, isXRay, isMultiEcho, isSlicesSpatiallySequentialPhilips, isValid, is3DAcq, is2DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase,isHasMagnitude,isHasMixed, isFloat, isResampled, isLocalizer; char phaseEncodingRC, patientSex; }; size_t nii_ImgBytes(struct nifti_1_header hdr); + int isSameFloatGE (float a, float b); struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, struct TDTI4D *dti4D); struct TDICOMdata readDICOM(char * fname); struct TDICOMdata clear_dicom_data(); diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 486ce914..ff22c9ce 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -109,7 +109,6 @@ void dropTrailingFileSep(char *path) { // path[len] = '\0'; } - void getFileName( char *pathParent, const char *path) //if path is c:\d1\d2 then filename is 'd2' { const char *filename = strrchr(path, '/'); //UNIX @@ -784,7 +783,8 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, json_Float(fp, "\t\"SAR\": %g,\n", d.SAR ); if (d.echoNum > 1) fprintf(fp, "\t\"EchoNumber\": %d,\n", d.echoNum); if ((d.TE > 0.0) && (!d.isXRay)) fprintf(fp, "\t\"EchoTime\": %g,\n", d.TE / 1000.0 ); - json_Float(fp, "\t\"RepetitionTime\": %g,\n", d.TR / 1000.0 ); + if ((d.TE2 > 0.0) && (!d.isXRay)) fprintf(fp, "\t\"EchoTime2\": %g,\n", d.TE2 / 1000.0 ); + json_Float(fp, "\t\"RepetitionTime\": %g,\n", d.TR / 1000.0 ); json_Float(fp, "\t\"InversionTime\": %g,\n", d.TI / 1000.0 ); json_Float(fp, "\t\"FlipAngle\": %g,\n", d.flipAngle ); float pf = 1.0f; //partial fourier @@ -874,7 +874,6 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, bandwidthPerPixelPhaseEncode = d.CSA.bandwidthPerPixelPhaseEncode; json_Float(fp, "\t\"BandwidthPerPixelPhaseEncode\": %g,\n", bandwidthPerPixelPhaseEncode ); if (d.accelFactPE > 1.0) fprintf(fp, "\t\"ParallelReductionFactorInPlane\": %g,\n", d.accelFactPE); - //EffectiveEchoSpacing // Siemens bandwidthPerPixelPhaseEncode already accounts for the effects of parallel imaging, // interpolation, phaseOversampling, and phaseResolution, in the context of the size of the @@ -885,7 +884,6 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, if (d.effectiveEchoSpacingGE > 0.0) effectiveEchoSpacing = d.effectiveEchoSpacingGE / 1000000.0; json_Float(fp, "\t\"EffectiveEchoSpacing\": %g,\n", effectiveEchoSpacing); - // Calculate true echo spacing (should match what Siemens reports on the console) // i.e., should match "echoSpacing" extracted from the ASCII CSA header, when that exists double trueESfactor = 1.0; @@ -909,7 +907,6 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, // https://github.com/rordenlab/dcm2niix/issues/130 if ((reconMatrixPE > 0) && (effectiveEchoSpacing > 0.0)) fprintf(fp, "\t\"TotalReadoutTime\": %g,\n", effectiveEchoSpacing * (reconMatrixPE - 1.0)); - json_Float(fp, "\t\"PixelBandwidth\": %g,\n", d.pixelBandwidth ); if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.dwellTime > 0)) fprintf(fp, "\t\"DwellTime\": %g,\n", d.dwellTime * 1E-9); @@ -1445,8 +1442,11 @@ int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopt } if (f == 'N') strcat (outname,dcm.patientName); - if (f == 'P') + if (f == 'P') { strcat (outname,dcm.protocolName); + if (strlen(dcm.protocolName) < 1) + printWarning("Unable to append protocol name (0018,1030) to filename (it is empty).\n"); + } if (f == 'Q') strcat (outname,dcm.scanningSequence); if (f == 'S') { @@ -1514,9 +1514,11 @@ int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopt sprintf(newstr, "_e%d", dcm.echoNum); strcat (outname,newstr); } - if (dcm.isHasPhase) - strcat (outname,"_ph"); //manufacturer name not available - + if (dcm.isHasPhase) { + strcat (outname,"_ph"); //has phase map + if (dcm.isHasMagnitude) + strcat (outname,"Mag"); //Philips enhanced with BOTH phase and Magnitude in single file + } if (strlen(outname) < 1) strcpy(outname, "dcm2nii_invalidName"); if (outname[0] == '.') outname[0] = '_'; //make sure not a hidden file //eliminate illegal characters http://msdn.microsoft.com/en-us/library/windows/desktop/aa365247(v=vs.85).aspx @@ -2341,7 +2343,7 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis printMessage("Ignoring derived image(s) of series %ld %s\n", dcmList[indx].seriesNum, nameList->str[indx]); return EXIT_SUCCESS; } - if ((opts.isIgnoreDerivedAnd2D) && ((strcmp(dcmList[indx].sequenceName, "_fl3d1_ns")== 0) || (strcmp(dcmList[indx].sequenceName, "_fl2d1")== 0)) ) { + if ((opts.isIgnoreDerivedAnd2D) && ((dcmList[indx].isLocalizer) || (strcmp(dcmList[indx].sequenceName, "_fl3d1_ns")== 0) || (strcmp(dcmList[indx].sequenceName, "_fl2d1")== 0)) ) { printMessage("Ignoring localizer (sequence %s) of series %ld %s\n", dcmList[indx].sequenceName, dcmList[indx].seriesNum, nameList->str[indx]); return EXIT_SUCCESS; } @@ -2672,11 +2674,11 @@ int compareTDCMsort(void const *item1, void const *item2) { return retval; } //compareTDCMsort()*/ -int isSameFloatGE (float a, float b) { +/*int isSameFloatGE (float a, float b) { //Kludge for bug in 0002,0016="DIGITAL_JACKET", 0008,0070="GE MEDICAL SYSTEMS" DICOM data: Orient field (0020:0037) can vary 0.00604261 == 0.00604273 !!! //return (a == b); //niave approach does not have any tolerance for rounding errors return (fabs (a - b) <= 0.0001); -} +}*/ int isSameFloatDouble (double a, double b) { //Kludge for bug in 0002,0016="DIGITAL_JACKET", 0008,0070="GE MEDICAL SYSTEMS" DICOM data: Orient field (0020:0037) can vary 0.00604261 == 0.00604273 !!! From 3509af9de63e95355fac1739dbe8e0b238835f18 Mon Sep 17 00:00:00 2001 From: Yaroslav Halchenko Date: Fri, 16 Mar 2018 00:28:57 -0400 Subject: [PATCH 05/14] ENH: minor - consistent formatting (%?=description) of template keywords for -f --- console/main_console.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/console/main_console.cpp b/console/main_console.cpp index dc976cad..13ec9209 100644 --- a/console/main_console.cpp +++ b/console/main_console.cpp @@ -85,7 +85,7 @@ void showHelp(const char * argv[], struct TDCMopts opts) { #else #define kQstr "" #endif - printf(" -f : filename (%%a=antenna (coil) number, %%c=comments, %%d=description, %%e echo number, %%f=folder name, %%i ID of patient, %%j seriesInstanceUID, %%k studyInstanceUID, %%m=manufacturer, %%n=name of patient, %%p=protocol,%s %%s=series number, %%t=time, %%u=acquisition number, %%v=vendor, %%x=study ID; %%z sequence name; default '%s')\n", kQstr, opts.filename); + printf(" -f : filename (%%a=antenna (coil) number, %%c=comments, %%d=description, %%e=echo number, %%f=folder name, %%i=ID of patient, %%j=seriesInstanceUID, %%k=studyInstanceUID, %%m=manufacturer, %%n=name of patient, %%p=protocol,%s %%s=series number, %%t=time, %%u=acquisition number, %%v=vendor, %%x=study ID; %%z=sequence name; default '%s')\n", kQstr, opts.filename); printf(" -g : generate defaults file (y/n/o [o=only: reset and write defaults], default n)\n"); printf(" -h : show help\n"); printf(" -i : ignore derived, localizer and 2D images (y/n, default n)\n"); From 1ef50d09fabe65ba5bbb0d153468a087bfcbe912 Mon Sep 17 00:00:00 2001 From: Chris Rorden Date: Sat, 17 Mar 2018 15:34:02 -0400 Subject: [PATCH 06/14] Fixes for Philips https://github.com/rordenlab/dcm2niix/issues/165 https://github.com/rordenlab/dcm2niix/issues/170 --- Philips/README.md | 38 ++ README.md | 11 + console/nii_dicom.cpp | 748 ++++++++++++++++++++++++++++-------- console/nii_dicom.h | 27 +- console/nii_dicom_batch.cpp | 140 +++++-- 5 files changed, 757 insertions(+), 207 deletions(-) create mode 100644 Philips/README.md diff --git a/Philips/README.md b/Philips/README.md new file mode 100644 index 00000000..dcbd5dfe --- /dev/null +++ b/Philips/README.md @@ -0,0 +1,38 @@ +## About + +dcm2niix attempts to convert all DICOM images to NIfTI. The Philips enhanced DICOM images are elegantly able to save all images from a series as a single file. However, this format is necessarily complex. The usage of this format has evolved over time, and can become further complicated when DICOM are handled by DICOM tools (for example, anonymization, transfer which converts explicit VRs to implicit VRs, etc.). + +This web page describes some of the strategies handle these images. However, users should be vigilant when handling these datasets. If you encounter problems using dcm2niix you can explore [https://github.com/rordenlab/dcm2niix competing DICOM to NIfTI converters] or [ [https://github.com/rordenlab/dcm2niix report an issue]. + +## Image Patient Position + +The Image Patient Position (0020,0032) tag is required to determine the position of the slices with respect to each other (and with respect to the scanner bore). Philips scans often report two conflicting IPPs for each single slice: with one stored in the private sequence (SQ) 2005,140F while the other is in the public sequence. This is unusual, but is [ftp://medical.nema.org/medical/dicom/final/cp758_ft.pdf legal]. + +In practice, this complication has several major implications. First, not all software translate private SQs well. One potential problem is when the implicit VRs are saved as implicit VRs. This can obscure the fact that 2005,140F is an SQ. Indeed, some tools will convert the private SQ type as a "UN" unknown type and add another public sequence. This can confuse reading software. + +Furthermore, in the real world there are many Philips DICOM images that ONLY contain IPP inside SQ 2005,140F. It is not clear if this reflect the source Philips data or handling by other tools. + +Therefore, dcm2niix will ignore the IPP enclosed in 2005,140F unless no alternative exists. + +## Derived Apparent Diffusion Coefficient maps stored with raw diffusion data + +The DICOM standard requires derived data to be saved as a separate series number than the raw data. However, many Philips diffusion images append the derived ADC maps with the original data. Strangely, some of these ADC maps seem to be stored with specific diffusion directions, and with the same volume number as the B=0 image. In these cases, 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 ADC map (3rd dimension is "1" while the other images set this to "2"). + +``` +(0018,9075) CS [DIRECTIONAL] +(0020,9157) UL 1\1\2\32 +... +(0018,9075) CS [ISOTROPIC] +(0020,9157) UL 1\2\1\33 +... +(0018,9075) CS [NONE] +(0020,9157) UL 1\1\2\33 +``` + +## Diffusion Direction + +Some Philips enhanced scans use tag 0018,9089 to report the 3 gradient directions. Other files use the tags 2005,10b0, 2005,10b1, 2005,10b2. In general, dcm2niix will use the values that most closely precede the Dimension Index Values (0020,9157). + +## General variations + +The tags SliceNumberMrPhilips (2001,100A), NumberOfSlicesMrPhilips (2001,1018), and InStackPositionNumber (0020,9057) MRImageGradientOrientationNumber (2005,1413) are all potentially very useful for simplifying handling of enhanced DICOM images. However, none of these are reliably stored in all DICOM images. Therefore, dcm2niix does not depend on any of these, though it attempts to use several of these if they are available. \ No newline at end of file diff --git a/README.md b/README.md index 29b5ecc1..c3f115f8 100644 --- a/README.md +++ b/README.md @@ -59,6 +59,17 @@ The batch processing binary `dcm2niibatch` is optional. To build `dcm2niibatch` If you have any problems with the cmake build script described above or want to customize the software see the [COMPILE.md file for details on manual compilation](./COMPILE.md). +## Alternatives + + - [Valerio Luccio's dinifti](http://cbi.nyu.edu/software/dinifti.php) is focused on conversion of Siemens data + - [dcm2nii](http://www.mccauslandcenter.sc.edu/mricro/mricron/dcm2nii.htm) is the predecessor of dcm2niix. It is deprecated for modern images, but does handle image formats that predate DICOM (proprietary Elscint, GE and Siemens formats). + - [Xiangrui Li's dicm2nii](http://www.mathworks.com/matlabcentral/fileexchange/42997-dicom-to-nifti-converter) is written in Matlab. The Matlab language makes this very scriptable. + - [dicom2nifti](https://github.com/icometrix/dicom2nifti) uses the scriptable Python wrapper utilizes the [high performance GDCMCONV](http://gdcm.sourceforge.net/wiki/index.php/Gdcmconv) executables. + - [MRtrix mrconvert](http://mrtrix.readthedocs.io/en/latest/reference/commands/mrconvert.html) is a useful general purpose image converter and handles DTI data well. It is an outstanding tool for modern Philips enhanced images. + - [Jolinda Smith's mcverter](http://lcni.uoregon.edu/%7Ejolinda/MRIConvert/) has great support for various vendors. + - [mri_convert](https://surfer.nmr.mgh.harvard.edu/pub/docs/html/mri_convert.help.xml.html) is part of the popular FreeSurfer package. In my limited experience this tool works well for GE and Siemens data, but fails with Philips 4D datasets. + - [SPM12](http://www.fil.ion.ucl.ac.uk/spm/software/spm12/) is one of the most popular tools in the field. It includes DICOM to NIfTI conversion. Being based on Matlab it is easy to script. + ## Links - [Dcm2Bids](https://github.com/cbedetti/Dcm2Bids) uses dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets. diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 49635e92..5abf9371 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -692,8 +692,8 @@ struct TDICOMdata clear_dicom_data() { strcpy(d.studyInstanceUID, ""); strcpy(d.bodyPartExamined,""); d.phaseEncodingLines = 0; - d.patientPositionSequentialRepeats = 0; - d.patientPositionRepeats = 0; + //~ d.patientPositionSequentialRepeats = 0; + //~ d.patientPositionRepeats = 0; d.isHasPhase = false; d.isHasMagnitude = false; d.sliceOrient = kSliceOrientUnknown; @@ -705,7 +705,6 @@ struct TDICOMdata clear_dicom_data() { d.lastScanLoc = NAN; d.TR = 0.0; d.TE = 0.0; - d.TE2 = 0.0; d.TI = 0.0; d.flipAngle = 0.0; d.bandwidthPerPixelPhaseEncode = 0.0; @@ -714,7 +713,7 @@ struct TDICOMdata clear_dicom_data() { d.pixelBandwidth = 0.0; d.zSpacing = 0.0; d.zThick = 0.0; - d.numberOfDynamicScans = 0; + //~ d.numberOfDynamicScans = 0; d.echoNum = 1; d.echoTrainLength = 0; d.phaseFieldofView = 0.0; @@ -722,7 +721,7 @@ struct TDICOMdata clear_dicom_data() { d.phaseEncodingSteps = 0; d.coilNum = 0; d.accelFactPE = 0.0; - d.patientPositionNumPhilips = 0; + //d.patientPositionNumPhilips = 0; d.imageBytes = 0; d.intenScale = 1; d.intenScalePhilips = 0; @@ -740,9 +739,10 @@ struct TDICOMdata clear_dicom_data() { d.imageStart = 0; d.is3DAcq = false; //e.g. MP-RAGE, SPACE, TFE d.is2DAcq = false; // - d.isSlicesSpatiallySequentialPhilips = true; //Philips can save slices in random order, e.g. 4,5,6,1,2,3 + //~ d.isSlicesSpatiallySequentialPhilips = true; //Philips can save slices in random order, e.g. 4,5,6,1,2,3 d.isDerived = false; //0008,0008 = DERIVED,CSAPARALLEL,POSDISP d.isSegamiOasis = false; //these images do not store spatial coordinates + d.isScaleOrTEVaries = false; d.bitsAllocated = 16;//bits d.bitsStored = 0; d.samplesPerPixel = 1; @@ -1424,6 +1424,7 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D #define kv2 45 #define kv3 46 #define kASL 48 + printError("dcm2niix PAR is not actively supported. Use with extreme caution (hint: use dicm2nii)\n"); char buff[LINESZ]; int sliceNumberMrPhilipsB2[kMaxDTI4D], sliceNumberMrPhilipsVol2[kMaxDTI4D]; int patientPositionNumPhilipsB2 = 0; @@ -1442,16 +1443,20 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D int maxCardiac = 1; int nCols = 26; int slice = 0; + int patientPositionNumPhilips = 0; + bool isSlicesSpatiallySequentialPhilips = true; //int prevSliceIndex = 0; //index of prior slice: detect if images are not in order const int kMaxCols = 49; float *cols = (float *)malloc(sizeof(float) * kMaxCols); char *p = fgets (buff, LINESZ, fp); bool isIntenScaleVaries = false; bool isIndexSequential = true; + int prevSlice = 0; for (int i = 0; i < kMaxDTI4D; i++) dti4D->S[i].V[0] = -1.0; - dti4D->S[0].sliceNumberMrPhilips = -1.0; - dti4D->S[0].sliceNumberMrPhilipsVol2 = -1.0; + dti4D->sliceOrder[0] = -1; + //dti4D->S[0].sliceNumberMrPhilips = -1.0; + //dti4D->S[0].sliceNumberMrPhilipsVol2 = -1.0; //d.dti4D = (TDTI *)malloc(kMaxDTI4D * sizeof(TDTI)); while (p) { if (strlen(buff) < 1) @@ -1569,11 +1574,12 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D printError("Slice dimensions or bit depth varies %s\n", parname); return d; } - if ((d.patientPositionSequentialRepeats == 0) && ((!isSameFloat(d.patientPosition[1],cols[kPositionRL])) || + //~ + /*if ((d.patientPositionSequentialRepeats == 0) && ((!isSameFloat(d.patientPosition[1],cols[kPositionRL])) || (!isSameFloat(d.patientPosition[2],cols[kPositionAP])) || (!isSameFloat(d.patientPosition[3],cols[kPositionFH])) ) )//this is the first slice with different position d.patientPositionSequentialRepeats = slice-1; - + */ if ((d.intenScale != cols[kRS]) || (d.intenIntercept != cols[kRI])) isIntenScaleVaries = true; } @@ -1601,14 +1607,26 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D d.patientPosition[2] = cols[kPositionAP]; d.patientPosition[3] = cols[kPositionFH]; } + //~ + /* if (d.patientPositionNumPhilips < kMaxDTI4D) { dti4D->S[d.patientPositionNumPhilips].sliceNumberMrPhilips = round(cols[kSlice]); if ((d.patientPositionNumPhilips > 0) && (dti4D->S[d.patientPositionNumPhilips].sliceNumberMrPhilips < dti4D->S[d.patientPositionNumPhilips-1].sliceNumberMrPhilips)) { - d.isSlicesSpatiallySequentialPhilips = false; - //printMessage("slices are not contiguous\n"); + d.isSlicesSpatiallySequentialPhilips = false; } } - d.patientPositionNumPhilips++; + d.patientPositionNumPhilips++;*/ + /* + if (patientPositionNumPhilips < kMaxDTI4D) { + dti4D->S[patientPositionNumPhilips].sliceNumberMrPhilips = round(cols[kSlice]); + if ((patientPositionNumPhilips > 0) && (dti4D->S[patientPositionNumPhilips].sliceNumberMrPhilips < dti4D->S[patientPositionNumPhilips-1].sliceNumberMrPhilips)) { + isSlicesSpatiallySequentialPhilips = false; + } + }*/ + int slice = round(cols[kSlice]); + if (slice != (prevSlice + 1)) isSlicesSpatiallySequentialPhilips = false; + prevSlice = slice; + patientPositionNumPhilips++; } if (cols[kGradientNumber] > 0) { /*int dir = (int) cols[kGradientNumber]; @@ -1618,7 +1636,6 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D //it seems like the ADC is always saved as the final volume, so this solution SHOULD be foolproof. ADCwarning = true; }*/ - //666: test (cols[kDyn] == 1) //(cols[kImageType] == 0) means magnitude scan if ((cols[kImageType] == 0) && (cols[kDyn] == 1) && (cols[kEcho] == 1) && (cols[kCardiac] == 1) && (cols[kSlice] == 1)) { //only first slice d.CSA.numDti++; @@ -1650,14 +1667,17 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D for (int s = 0; s < patientPositionNumPhilipsVol2; s++) sliceNumberMrPhilipsVol2[s] = sliceNumberMrPhilipsB2[s]; } - if ((patientPositionNumPhilipsVol2 > 1) && (d.patientPositionNumPhilips == patientPositionNumPhilipsVol2)) { + if (!isSlicesSpatiallySequentialPhilips) + printError("PAR file order of slices is not sequential (hint: use dicm2nii)\n"); + //~ + /*if ((patientPositionNumPhilipsVol2 > 1) && (d.patientPositionNumPhilips == patientPositionNumPhilipsVol2)) { bool isSliceOrderConsistent = true; for (int s = 0; s < patientPositionNumPhilipsVol2; s++) if (sliceNumberMrPhilipsVol2[s] != dti4D->S[s].sliceNumberMrPhilips) isSliceOrderConsistent = false; if (!isSliceOrderConsistent) printError("PAR file order of slices varies between volumes (hint: use dicm2nii)\n"); d.isValid = false; - } + }*/ if (dynNotAscending) { printError("PAR file volumes not saved in ascending order (hint: use dicm2nii)\n"); d.isValid = false; @@ -1913,15 +1933,16 @@ unsigned char * nii_flipImgZ(unsigned char* bImg, struct nifti_1_header *hdr){ return bImg; } // nii_flipImgZ() -unsigned char * nii_reorderSlices(unsigned char* bImg, struct nifti_1_header *h, struct TDTI4D *dti4D){ +/*unsigned char * nii_reorderSlices(unsigned char* bImg, struct nifti_1_header *h, struct TDTI4D *dti4D){ //flip slice order - Philips scanners can save data in non-contiguous order //if ((h->dim[3] < 2) || (h->dim[4] > 1)) return bImg; + return bImg; if (h->dim[3] < 2) return bImg; if (h->dim[3] >= kMaxDTI4D) { printWarning("Unable to reorder slices (%d > %d)\n", h->dim[3], kMaxDTI4D); return bImg; } - //printMessage("<<< Slices not spatially contiguous: please check output [new feature]\n"); + printError("OBSOLETE<<< Slices not spatially contiguous: please check output [new feature]\n"); return bImg; int dim4to7 = 1; for (int i = 4; i < 8; i++) if (h->dim[i] > 1) dim4to7 = dim4to7 * h->dim[i]; @@ -1947,7 +1968,7 @@ unsigned char * nii_reorderSlices(unsigned char* bImg, struct nifti_1_header *h, free(srcImg); return bImg; }// nii_reorderSlices() - +*/ unsigned char * nii_flipZ(unsigned char* bImg, struct nifti_1_header *h){ //flip slice order if (h->dim[3] < 2) return bImg; @@ -2169,7 +2190,89 @@ unsigned char * nii_iVaries(unsigned char *img, struct nifti_1_header *hdr){ return (unsigned char*) img32; } //nii_iVaries() -unsigned char * nii_XYTZ_XYZT(unsigned char* bImg, struct nifti_1_header *hdr, int seqRepeats) { + +/*unsigned char * nii_reorderSlicesX(unsigned char* bImg, struct nifti_1_header *hdr, struct TDTI4D *dti4D) { +//Philips can save slices in any random order... rearrange all of them + int dim3to7 = 1; + for (int i = 3; i < 8; i++) + if (hdr->dim[i] > 1) dim3to7 = dim3to7 * hdr->dim[i]; + if (dim3to7 < 2) return bImg; + uint64_t sliceBytes = hdr->dim[1]*hdr->dim[2]*hdr->bitpix/8; + unsigned char *sliceImg = (unsigned char *)malloc( sliceBytes); + uint32_t *idx = (uint32_t *)malloc( dim3to7 * sizeof(uint32_t)); + for (int i = 0; i < dim3to7; i++) //for each volume + idx[i] = i; + for (int i = 0; i < dim3to7; i++) { //for each volume + int fromSlice = idx[dti4D->sliceOrder[i]]; + //if (i < 10) printMessage(" ===> Moving slice from/to positions\t%d\t%d\n", i, toSlice); + if (i != fromSlice) { + uint64_t inPos = fromSlice * sliceBytes; + uint64_t outPos = i * sliceBytes; + memcpy( &sliceImg[0], &bImg[outPos], sliceBytes); //dest, src -> copy slice about to be overwritten + memcpy( &bImg[outPos], &bImg[inPos], sliceBytes); //dest, src + memcpy( &bImg[inPos], &sliceImg[0], sliceBytes); // + idx[i] = fromSlice; + } + } + free(idx); + free(sliceImg); + return bImg; +}*/ + +unsigned char * nii_reorderSlicesX(unsigned char* bImg, struct nifti_1_header *hdr, struct TDTI4D *dti4D) { +//Philips can save slices in any random order... rearrange all of them + int dim3to7 = 1; + for (int i = 3; i < 8; i++) + if (hdr->dim[i] > 1) dim3to7 = dim3to7 * hdr->dim[i]; + if (dim3to7 < 2) return bImg; + uint64_t imgSz = nii_ImgBytes(*hdr); + uint64_t sliceBytes = hdr->dim[1]*hdr->dim[2]*hdr->bitpix/8; + unsigned char *outImg = (unsigned char *)malloc( imgSz); + memcpy( &outImg[0],&bImg[0], imgSz); + for (int i = 0; i < dim3to7; i++) { //for each volume + int fromSlice = dti4D->sliceOrder[i]; + //if (i < 10) printMessage(" ===> Moving slice from/to positions\t%d\t%d\n", i, toSlice); + if (i != fromSlice) { + uint64_t inPos = fromSlice * sliceBytes; + uint64_t outPos = i * sliceBytes; + memcpy( &bImg[outPos], &outImg[inPos], sliceBytes); + } + } + free(outImg); + return bImg; +} + +/*unsigned char * nii_reorderSlicesX(unsigned char* bImg, struct nifti_1_header *hdr, struct TDTI4D *dti4D) { +//Philips can save slices in any random order... rearrange all of them + //if (sliceOrder == NULL) return bImg; + int dim3to7 = 1; + for (int i = 3; i < 8; i++) + if (hdr->dim[i] > 1) dim3to7 = dim3to7 * hdr->dim[i]; + if (dim3to7 < 2) return bImg; + //printMessage("xxx NOT reordering %d Philips slices.\n", dim3to7); return bImg; + uint64_t sliceBytes = hdr->dim[1]*hdr->dim[2]*hdr->bitpix/8; + unsigned char *sliceImg = (unsigned char *)malloc( sliceBytes); + //for (int i = 0; i < dim3to7; i++) { //for each volume + uint64_t imgSz = nii_ImgBytes(*hdr); + //this uses a lot of RAM, someday this could be done in place... + unsigned char *outImg = (unsigned char *)malloc( imgSz); + memcpy( &outImg[0],&bImg[0], imgSz); + + for (int i = 0; i < dim3to7; i++) { //for each volume + int toSlice = dti4D->sliceOrder[i]; + //if (i < 10) printMessage(" ===> Moving slice from/to positions\t%d\t%d\n", i, toSlice); + if (i != toSlice) { + uint64_t inPos = i * sliceBytes; + uint64_t outPos = toSlice * sliceBytes; + memcpy( &bImg[outPos], &outImg[inPos], sliceBytes); + } + } + free(sliceImg); + free(outImg); + return bImg; +}*/ + +/*unsigned char * nii_XYTZ_XYZT(unsigned char* bImg, struct nifti_1_header *hdr, int seqRepeats) { //Philips can save time as 3rd dimensions, NIFTI requires time is 4th dimension int dim4to7 = 1; for (int i = 4; i < 8; i++) @@ -2203,6 +2306,7 @@ unsigned char * nii_XYTZ_XYZT(unsigned char* bImg, struct nifti_1_header *hdr, i free(bImg); return outImg; } //nii_XYTZ_XYZT() +*/ unsigned char * nii_byteswap(unsigned char *img, struct nifti_1_header *hdr){ if (hdr->bitpix < 9) return img; @@ -2655,7 +2759,7 @@ unsigned char * nii_loadImgRLE(char* imgname, struct nifti_1_header hdr, struct #endif #endif -unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct TDICOMdata dcm, bool iVaries, int compressFlag, int isVerbose) { +unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct TDICOMdata dcm, bool iVaries, int compressFlag, int isVerbose, struct TDTI4D *dti4D) { //provided with a filename (imgname) and DICOM header (dcm), creates NIfTI header (hdr) and img if (headerDcm2Nii(dcm, hdr, true) == EXIT_FAILURE) return NULL; //TOFU unsigned char * img; @@ -2712,14 +2816,17 @@ unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct hdr->dim[3] = nAcq; hdr->dim[0] = 4; } - if ((hdr->dim[0] > 3) && (dcm.patientPositionSequentialRepeats > 1)) //swizzle 3rd and 4th dimension (Philips stores time as 3rd dimension) - img = nii_XYTZ_XYZT(img, hdr,dcm.patientPositionSequentialRepeats); - if (((dcm.patientPositionSequentialRepeats * 2) == dcm.patientPositionRepeats) && (dcm.isHasPhase) && (dcm.isHasMagnitude)) { + //~ if ((hdr->dim[0] > 3) && (dcm.patientPositionSequentialRepeats > 1) && (dcm.sliceOrder == NULL)) //swizzle 3rd and 4th dimension (Philips stores time as 3rd dimension) + //~ img = nii_XYTZ_XYZT(img, hdr,dcm.patientPositionSequentialRepeats); + if (dti4D->sliceOrder[0] >= 0) + img = nii_reorderSlicesX(img, hdr, dti4D); + //~ + /*if (((dcm.patientPositionSequentialRepeats * 2) == dcm.patientPositionRepeats) && (dcm.isHasPhase) && (dcm.isHasMagnitude)) { hdr->dim[3] = hdr->dim[3] / 2; hdr->dim[4] = hdr->dim[4] * 2; hdr->dim[0] = 4; printMessage("Splitting Phase+Magnitude into two volumes for %d slices (Z) and %d volumes (T).\n",hdr->dim[3], hdr->dim[4]); - } + }*/ headerDcm2NiiSForm(dcm,dcm, hdr, false); return img; } //nii_loadImgXL() @@ -2952,6 +3059,27 @@ void _update_tvd(struct TVolumeDiffusion* ptvd) { }//_update_tvd() //END RIR +struct TDCMdim { //DimensionIndexValues + uint32_t dimIdx[MAX_NUMBER_OF_DIMENSIONS]; + uint32_t diskPos; + float TE, intenScale, intenIntercept, intenScalePhilips; + bool isPhase; +}; + +int compareTDCMdim(void const *item1, void const *item2) { + struct TDCMdim const *dcm1 = (const struct TDCMdim *)item1; + struct TDCMdim const *dcm2 = (const struct TDCMdim *)item2; + //for(int i=0; i < MAX_NUMBER_OF_DIMENSIONS; ++i){ + for(int i=MAX_NUMBER_OF_DIMENSIONS-1; i >=0; i--){ + + if(dcm1->dimIdx[i] < dcm2->dimIdx[i]) + return -1; + else if(dcm1->dimIdx[i] > dcm2->dimIdx[i]) + return 1; + } + return 0; +} //compareTDCMdim() + struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, struct TDTI4D *dti4D) { struct TDICOMdata d = clear_dicom_data(); d.imageNum = 0; //not set @@ -2960,9 +3088,9 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru strcpy(d.seriesDescription, ""); //erase dummy with empty strcpy(d.sequenceName, ""); //erase dummy with empty //do not read folders - code specific to GCC (LLVM/Clang seems to recognize a small file size) - dti4D->S[0].sliceNumberMrPhilips = -1; - dti4D->S[0].sliceNumberMrPhilipsVol2 = -1; - + //dti4D->S[0].sliceNumberMrPhilips = -1; + //dti4D->S[0].sliceNumberMrPhilipsVol2 = -1; + dti4D->sliceOrder[0] = -1; struct TVolumeDiffusion volDiffusion = initTVolumeDiffusion(&d, dti4D); struct stat s; @@ -3096,6 +3224,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kDiffusionOrientation 0x0018+uint32_t(0x9089<< 16 ) // FD, seen in enhanced // DICOM from Philips 5.* // and Siemens XA10. +//#define kMREchoSequence 0x0018+uint32_t(0x9114<< 16 ) //SQ #define kDwellTime 0x0019+(0x1018<< 16 ) //IS in NSec, see https://github.com/rordenlab/dcm2niix/issues/127 #define kLastScanLoc 0x0019+(0x101B<< 16 ) #define kDiffusionDirectionGEX 0x0019+(0x10BB<< 16 ) //DS @@ -3110,12 +3239,14 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kImageNum 0x0020+(0x0013 << 16 ) #define kStudyInstanceUID 0x0020+(0x000D << 16 ) #define kSeriesInstanceUID 0x0020+(0x000E << 16 ) -#define kPatientPosition 0x0020+(0x0032 << 16 ) // Actually ImagePositionPatient! +#define kImagePositionPatient 0x0020+(0x0032 << 16 ) // Actually ! #define kOrientationACR 0x0020+(0x0035 << 16 ) +#define kTemporalPositionIdentifier 0x0020+(0x0100 << 16 ) //IS #define kOrientation 0x0020+(0x0037 << 16 ) #define kImagesInAcquisition 0x0020+(0x1002 << 16 ) //IS #define kImageComments 0x0020+(0x4000<< 16 )// '0020' '4000' 'LT' 'ImageComments' #define kDimensionIndexValues 0x0020+uint32_t(0x9157<< 16 ) // UL n-dimensional index of frame. +#define kInStackPositionNumber 0x0020+uint32_t(0x9057<< 16 ) // UL can help determine slices in volume #define kLocationsInAcquisitionGE 0x0021+(0x104F<< 16 )// 'SS' 'LocationsInAcquisitionGE' #define kSamplesPerPixel 0x0028+(0x0002 << 16 ) #define kPhotometricInterpretation 0x0028+(0x0004 << 16 ) @@ -3148,9 +3279,9 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kIconImageSequence 0x0088+(0x0200 << 16 ) #define kDiffusionBFactor 0x2001+(0x1003 << 16 )// FL #define kSliceNumberMrPhilips 0x2001+(0x100A << 16 ) //IS Slice_Number_MR -#define kNumberOfSlicesMrPhilips 0x2001+(0x1018 << 16 )//SL 0x2001, 0x1018 ), "Number_of_Slices_MR" #define kSliceOrient 0x2001+(0x100B << 16 )//2001,100B Philips slice orientation (TRANSVERSAL, AXIAL, SAGITTAL) -//#define kLocationsInAcquisitionPhilips 0x2001+(0x1018 << 16 ) // +#define kNumberOfSlicesMrPhilips 0x2001+(0x1018 << 16 )//SL 0x2001, 0x1018 ), "Number_of_Slices_MR" +//#define kNumberOfLocationsPhilips 0x2001+(0x1015 << 16 ) //SS //#define kStackSliceNumber 0x2001+(0x1035 << 16 )//? Potential way to determine slice order for Philips? #define kNumberOfDynamicScans 0x2001+(0x1081 << 16 )//'2001' '1081' 'IS' 'NumberOfDynamicScans' #define kMRAcquisitionTypePhilips 0x2005+(0x106F << 16) @@ -3165,6 +3296,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kDiffusionDirectionAP 0x2005+(0x10B1 << 16) #define kDiffusionDirectionFH 0x2005+(0x10B2 << 16) #define kPrivatePerFrameSq 0x2005+(0x140F << 16) +#define kMRImageGradientOrientationNumber 0x2005+(0x1413 << 16) //IS #define kWaveformSq 0x5400+(0x0100 << 16) #define kImageStart 0x7FE0+(0x0010 << 16 ) #define kImageStartFloat 0x7FE0+(0x0008 << 16 ) @@ -3172,13 +3304,25 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru uint32_t kItemTag = 0xFFFE +(0xE000 << 16 ); uint32_t kItemDelimitationTag = 0xFFFE +(0xE00D << 16 ); uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); +double TE = 0.0; //most recent echo time recorded + bool is2005140FSQ = false; int locationsInAcquisitionGE = 0; + int inStackPositionNumber = 0; + int maxInStackPositionNumber = 0; + int temporalPositionIdentifier = 0; int locationsInAcquisitionPhilips = 0; int imagesInAcquisition = 0; + //int sumSliceNumberMrPhilips = 0; int sliceNumberMrPhilips = 0; + int numberOfFrames = 0; + int MRImageGradientOrientationNumber = 0; + //int maxMRImageGradientOrientationNumber = -1; + int maxGradNum = -1; + int numberOfDynamicScans = 0; uint32_t lLength; uint32_t groupElement; long lPos = 0; + bool isPhilipsADC = false; if (isPart10prefix) { //for part 10 files, skip preamble and prefix lPos = 128+4; //4-byte signature starts at 128 groupElement = buffer[lPos] | (buffer[lPos+1] << 8) | (buffer[lPos+2] << 16) | (buffer[lPos+3] << 24); @@ -3197,18 +3341,40 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); bool isIconImageSequence = false; bool isSwitchToImplicitVR = false; bool isSwitchToBigEndian = false; - //bool geiisBug = false; //for buggy GEIIS http://forum.dcmtk.org/viewtopic.php?p=7162&sid=3b516cc751aae51fbb5e73184abe37c2 - bool is2005140FSQ = false; //for buggy Philips - bool is2005140FSQwarned = false; //for buggy Philips bool isAtFirstPatientPosition = false; //for 3d and 4d files: flag is true for slices at same position as first slice bool isMosaic = false; int patientPositionNum = 0; - int sqDepth = 0; - int sqDepthPrivate = 0; - int sqEndPrivate = -1; //used to skip private SQs that provide explicit length + float B0Philips = -1.0; + float vRLPhilips = 0.0; + float vAPPhilips = 0.0; + float vFHPhilips = 0.0; + bool isPhase = false; + bool isMagnitude = false; + float patientPositionPrivate[4] = {NAN, NAN, NAN, NAN}; float patientPosition[4] = {NAN, NAN, NAN, NAN}; //used to compute slice direction for Philips 4D + float patientPositionPublic[4] = {NAN, NAN, NAN, NAN}; //used to compute slice direction for Philips 4D float patientPositionEndPhilips[4] = {NAN, NAN, NAN, NAN}; float patientPositionStartPhilips[4] = {NAN, NAN, NAN, NAN}; +/* //not required? + int sqDepth = 0; + int sqDepthPrivate = 0; + int sqEndPrivate = -1; //used to skip private SQs that provide explicit length + bool is2005140FSQ = false; //Philips stores 2D slice data here + bool is2005140FSQwarned = false; //for buggy Philips + +*/ + //array for storing Philips DTI when tags kMRImageGradientOrientationNumber and kDimensionIndexValues are available + struct TDTI philDTI[kMaxDTI4D]; + for (int i = 0; i < kMaxDTI4D; i++) + philDTI[i].V[0] = -1; + //array for storing DimensionIndexValues + int numDimensionIndexValues = 0; + TDCMdim dcmDim[kMaxSlice2D]; + for (int i = 0; i < kMaxSlice2D; i++) { + dcmDim[i].diskPos = i; + for (int j = 0; j < MAX_NUMBER_OF_DIMENSIONS; j++) + dcmDim[i].dimIdx[j] = 0; + } while ((d.imageStart == 0) && ((lPos+8+lFileOffset) < fileLen)) { #ifndef myLoadWholeFileToReadHeader //read one segment at a time if ((size_t)(lPos + 128) > MaxBufferSz) { //avoid overreading the file @@ -3266,7 +3432,6 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); } else if ((buffer[lPos] == 'S') && (buffer[lPos+1] == 'Q')) { lLength = 8; //Sequence Tag //printMessage(" !!!SQ\t%04x,%04x\n", groupElement & 65535,groupElement>>16); - is2005140FSQ = (groupElement == kPrivatePerFrameSq); } else { //explicit VR with 16-bit length if ((d.isLittleEndian) ) lLength = buffer[lPos+2] | (buffer[lPos+3] << 8); @@ -3275,8 +3440,8 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); lPos += 4; //skip 2 byte VR string and 2 length bytes = 4 bytes } } else { //implicit VR - vr[0] = 'N'; - vr[1] = 'A'; + vr[0] = 'U'; + vr[1] = 'N'; if (d.isLittleEndian) lLength = buffer[lPos] | (buffer[lPos+1] << 8) | (buffer[lPos+2] << 16) | (buffer[lPos+3] << 24); else @@ -3286,8 +3451,6 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); vr[0] = 'S'; vr[1] = 'Q'; lLength = 8; //Sequence Tag - is2005140FSQ = (groupElement == kPrivatePerFrameSq); - //if (is2005140FSQ) printMessage(" !!!SQ\t%04x,%04x\n", groupElement & 65535,groupElement>>16); } } //if explicit else implicit VR if (lLength == 0xFFFFFFFF) { @@ -3302,6 +3465,7 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); vr[1] = 'Q'; //} } + /* //Handle SQs: for explicit these have VR=SQ if ((vr[0] == 'S') && (vr[1] == 'Q')) { //http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_7.5.html uint32_t special = dcmInt(4,&buffer[lPos],d.isLittleEndian); @@ -3311,20 +3475,23 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); uint32_t kUndefinedLen = 0xFFFFFFFF; //printError(" SPECIAL >>>>t%04x,%04x %08x %08x\n", groupElement & 65535,groupElement>>16, special, slen); //return d; + is2005140FSQ = (groupElement == kPrivatePerFrameSq); + //if (isNextSQis2005140FSQ) is2005140FSQ = true; + //isNextSQis2005140FSQ = false; if (special == kSequenceDelimitationItemTag) { - //sqDepth--; - //printMessage(" SPECIAL >>>>t%04x,%04x %08x %d\n", groupElement & 65535,groupElement>>16, special, slen); + //unknown } else if (slen == kUndefinedLen) { sqDepth++; if ((sqDepthPrivate == 0) && ((groupElement & 65535) % 2)) sqDepthPrivate = sqDepth; //in a private SQ: ignore contents - } else if ((groupElement & 65535) % 2) //private SQ of known length - lets jump over this! - slen = lFileOffset + slen; - if (slen > sqEndPrivate) + } else if ((is2005140FSQ) || ((groupElement & 65535) % 2)) {//private SQ of known length - lets jump over this! + slen = lFileOffset + lPos + slen; + if ((sqEndPrivate < 0) || (slen > sqEndPrivate)) //sqEndPrivate is signed sqEndPrivate = slen; //if nested private SQs, remember the end address of the top parent SQ + } } //next: look for required tags - if ((groupElement == kItemTag) && (isEncapsulatedData)) { + if ((groupElement == kItemTag) && (isEncapsulatedData)) { d.imageBytes = dcmInt(4,&buffer[lPos-4],d.isLittleEndian); //printMessage("compressed data %d-> %ld\n",d.imageBytes, lPos); if (d.imageBytes > 128) { @@ -3334,14 +3501,17 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); } } if ((isIconImageSequence) && ((groupElement & 0x0028) == 0x0028 )) groupElement = kUnused; //ignore icon dimensions - if (groupElement == kSequenceDelimitationItemTag) { + if ((sqEndPrivate > 0) && ((lFileOffset + lPos) > sqEndPrivate)) + sqEndPrivate = -1; //end of private SQ with defined length + if (groupElement == kSequenceDelimitationItemTag) { //end of private SQ with undefined length sqDepth--; if (sqDepth < sqDepthPrivate) { sqDepthPrivate = 0; //no longer in a private SQ - is2005140FSQ = false; } } - if (sqDepth < 0) sqDepth = 0; + + if (sqDepth < 0) sqDepth = 0;*/ + if (groupElement == kSequenceDelimitationItemTag) is2005140FSQ = false; switch ( groupElement ) { case kTransferSyntax: { char transferSyntax[kDICOMStr]; @@ -3455,11 +3625,15 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); break; case kComplexImageComponent: if (lLength < 2) break; - //not: it is possible for Philips enhanced DICOM to store BOTH magnitude and phase in the same image + isPhase = false; + isMagnitude = false; if ((buffer[lPos]=='P') && (toupper(buffer[lPos+1]) == 'H')) - d.isHasPhase = true; + isPhase = true; if ((buffer[lPos]=='M') && (toupper(buffer[lPos+1]) == 'A')) - d.isHasMagnitude = true; + isMagnitude = true; + //not mutually exclusive: possible for Philips enhanced DICOM to store BOTH magnitude and phase in the same image + if (isPhase) d.isHasPhase = true; + if (isMagnitude) d.isHasMagnitude = true; break; case kAcquisitionTime : char acquisitionTimeTxt[kDICOMStr]; @@ -3529,9 +3703,14 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); case kPatientOrient : dcmStr (lLength, &buffer[lPos], d.patientOrient); break; - case kDiffusionDirectionality : // 0018, 9075 + case kDiffusionDirectionality : {// 0018, 9075 set_directionality0018_9075(&volDiffusion, (&buffer[lPos])); - break; + if ((d.manufacturer != kMANUFACTURER_PHILIPS) || (lLength < 10)) break; + char dir[kDICOMStr]; + dcmStr (lLength, &buffer[lPos], dir); + if (strcmp(dir, "ISOTROPIC") == 0) + isPhilipsADC = true; + break; } case kDwellTime : d.dwellTime = dcmStrInt(lLength, &buffer[lPos]); break; @@ -3564,43 +3743,41 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); case kSeriesInstanceUID : // 0020, 000E dcmStr (lLength, &buffer[lPos], d.seriesInstanceUID); break; - case kPatientPosition : { // 0020,0032, ImagePositionPatient - bool ignorePatientPos = false; - if (d.manufacturer == kMANUFACTURER_PHILIPS) { - if (is2005140FSQ) ignorePatientPos = true; - if ((lFileOffset + lPos) < sqEndPrivate) ignorePatientPos = true; //inside private SQ, SQ has defined length - if (sqDepthPrivate > 0) ignorePatientPos = true; //inside private SQ, SQ has undefined length - } - if (ignorePatientPos) { - if (!is2005140FSQwarned) - printWarning("Philips R3.2.2 can report different positions for the same slice. Attempting patch.\n"); - is2005140FSQwarned = true; - } else { - patientPositionNum++; - isAtFirstPatientPosition = true; - dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, &patientPosition[0]); //slice position - if (isnan(d.patientPosition[1])) { - //dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, &d.patientPosition[0]); //slice position - for (int k = 0; k < 4; k++) - d.patientPosition[k] = patientPosition[k]; - } else { - //dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, &d.patientPositionLast[0]); //slice direction for 4D - for (int k = 0; k < 4; k++) - d.patientPositionLast[k] = patientPosition[k]; - if ((isFloatDiff(d.patientPositionLast[1],d.patientPosition[1])) || - (isFloatDiff(d.patientPositionLast[2],d.patientPosition[2])) || - (isFloatDiff(d.patientPositionLast[3],d.patientPosition[3])) ) { - isAtFirstPatientPosition = false; //this slice is not at position of 1st slice - if (d.patientPositionSequentialRepeats == 0) //this is the first slice with different position - d.patientPositionSequentialRepeats = patientPositionNum-1; - } //if different position from 1st slice in file - } //if not first slice in file - if (isAtFirstPatientPosition) d.patientPositionRepeats ++; - set_isAtFirstPatientPosition_tvd(&volDiffusion, isAtFirstPatientPosition); - if (isVerbose > 1) - printMessage(" Patient Position 0020,0032 (#,#repeat,@,X,Y,Z)\t%d\t%ld\t%g\t%g\t%g\n", patientPositionNum, d.patientPositionSequentialRepeats, lPos, patientPosition[1], patientPosition[2], patientPosition[3]); - } //not after 2005,140F - break; } + case kImagePositionPatient : { + if (is2005140FSQ) { + dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, &patientPositionPrivate[0]); + break; + } + patientPositionNum++; + isAtFirstPatientPosition = true; + dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, &patientPosition[0]); //slice position + + /*if (!is2005140FSQ) { + printMessage("cxcPublic Patient Position 0020,0032 (#,@,X,Y,Z)\t%d\t%ld\t%g\t%g\t%g\n", patientPositionNum, lPos, patientPosition[1], patientPosition[2], patientPosition[3]); + for (int k = 0; k < 4; k++) + patientPositionPublic[k] = patientPosition[k]; + }*/ + + if (isnan(d.patientPosition[1])) { + //dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, &d.patientPosition[0]); //slice position + for (int k = 0; k < 4; k++) + d.patientPosition[k] = patientPosition[k]; + } else { + //dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, &d.patientPositionLast[0]); //slice direction for 4D + for (int k = 0; k < 4; k++) + d.patientPositionLast[k] = patientPosition[k]; + if ((isFloatDiff(d.patientPositionLast[1],d.patientPosition[1])) || + (isFloatDiff(d.patientPositionLast[2],d.patientPosition[2])) || + (isFloatDiff(d.patientPositionLast[3],d.patientPosition[3])) ) { + isAtFirstPatientPosition = false; //this slice is not at position of 1st slice + //if (d.patientPositionSequentialRepeats == 0) //this is the first slice with different position + // d.patientPositionSequentialRepeats = patientPositionNum-1; + } //if different position from 1st slice in file + } //if not first slice in file + set_isAtFirstPatientPosition_tvd(&volDiffusion, isAtFirstPatientPosition); + if (isVerbose > 1) + printMessage(" Patient Position 0020,0032 (#,@,X,Y,Z)\t%d\t%ld\t%g\t%g\t%g\n", patientPositionNum, lPos, patientPosition[1], patientPosition[2], patientPosition[3]); + break; } case kInPlanePhaseEncodingDirection: d.phaseEncodingRC = toupper(buffer[lPos]); //first character is either 'R'ow or 'C'ol break; @@ -3618,12 +3795,18 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); d.acquNum = dcmStrInt(lLength, &buffer[lPos]); break; case kImageNum: + //Enhanced Philips also uses this in once per file SQ 0008,1111 + //Enhanced Philips also uses this once per slice in SQ 2005,140f if (d.imageNum < 1) d.imageNum = dcmStrInt(lLength, &buffer[lPos]); //Philips renames each image again in 2001,9000, which can lead to duplicates - break; - case kDimensionIndexValues: // kImageNum is not enough for 4D series from Philips 5.*. - { // { necessary for initializing ndim. + break; + case kInStackPositionNumber: + if (d.manufacturer != kMANUFACTURER_PHILIPS) break; + inStackPositionNumber = dcmInt(4,&buffer[lPos],d.isLittleEndian); + if (inStackPositionNumber > maxInStackPositionNumber) maxInStackPositionNumber = inStackPositionNumber; + break; + case kDimensionIndexValues: { // kImageNum is not enough for 4D series from Philips 5.*. + if (lLength < 4) break; uint8_t ndim = lLength / 4; - if(ndim > MAX_NUMBER_OF_DIMENSIONS){ // Ideally degenerate axes would be cleverly handled. // They are commonly seen in NIfTIs, but perhaps not in DICOM. @@ -3632,20 +3815,106 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); ndim = MAX_NUMBER_OF_DIMENSIONS; // Truncate } dcmMultiLongs(4 * ndim, &buffer[lPos], ndim, d.dimensionIndexValues, d.isLittleEndian); + if (d.manufacturer != kMANUFACTURER_PHILIPS) break; + // + if (inStackPositionNumber > 0) { + //for images without SliceNumberMrPhilips (2001,100A) + int sliceNumber = inStackPositionNumber; + if ((sliceNumber == 1) && (!isnan(patientPosition[1])) ) { + for (int k = 0; k < 4; k++) + patientPositionStartPhilips[k] = patientPosition[k]; + } else if ((sliceNumber == 1) && (!isnan(patientPositionPrivate[1]))) { + for (int k = 0; k < 4; k++) + patientPositionStartPhilips[k] = patientPositionPrivate[k]; + } + if ((sliceNumber == maxInStackPositionNumber) && (!isnan(patientPosition[1])) ) { + for (int k = 0; k < 4; k++) + patientPositionEndPhilips[k] = patientPosition[k]; + } else if ((sliceNumber == maxInStackPositionNumber) && (!isnan(patientPositionPrivate[1])) ) { + for (int k = 0; k < 4; k++) + patientPositionEndPhilips[k] = patientPositionPrivate[k]; + } + patientPosition[1] = NAN; + patientPositionPrivate[1] = NAN; + } + inStackPositionNumber = 0; + if (numDimensionIndexValues >= kMaxSlice2D) { + printError("Too many slices to track dimensions. Only up to %d are supported\n", kMaxSlice2D); + break; + } + for (int i = 0; i < ndim; i++) + dcmDim[numDimensionIndexValues].dimIdx[i] = d.dimensionIndexValues[i]; + dcmDim[numDimensionIndexValues].TE = TE; + dcmDim[numDimensionIndexValues].intenScale = d.intenScale; + dcmDim[numDimensionIndexValues].intenIntercept = d.intenIntercept; + dcmDim[numDimensionIndexValues].isPhase = isPhase; + dcmDim[numDimensionIndexValues].intenScalePhilips = d.intenScalePhilips; + numDimensionIndexValues ++; + //next: add diffusion if reported + if (B0Philips < 0.0) break; //diffusion parameters not reported + // Philips does not always provide 2005,1413 (MRImageGradientOrientationNumber) and sometimes after dimensionIndexValues + int gradNum = 0; + for (int i = 0; i < ndim; i++) + if (d.dimensionIndexValues[i] > 0) gradNum = d.dimensionIndexValues[i]; + if (gradNum <= 0) break; + /* + With Philips 51.0 both ADC and B=0 are saved as same direction, though they vary in another dimension + (0018,9075) CS [ISOTROPIC] + (0020,9157) UL 1\2\1\33 << ADC MAP + (0018,9075) CS [NONE] + (0020,9157) UL 1\1\2\33 + next two lines attempt to skip ADC maps + we could also increment gradNum for ADC if we wanted... + */ + if (isPhilipsADC) { + gradNum ++; + B0Philips = 2000.0; + vRLPhilips = 0.0; + vAPPhilips = 0.0; + vFHPhilips = 0.0; + + } + //if ((MRImageGradientOrientationNumber > 0) && ((gradNum != MRImageGradientOrientationNumber)) break; + + if (gradNum >= maxGradNum) maxGradNum = gradNum; + if (gradNum >= kMaxDTI4D) { + printError("Number of DTI gradients exceeds 'kMaxDTI4D (%d).\n", kMaxDTI4D); + break; + } + gradNum = gradNum - 1;//index from 0 + philDTI[gradNum].V[0] = B0Philips; + philDTI[gradNum].V[1] = vRLPhilips; + philDTI[gradNum].V[2] = vAPPhilips; + philDTI[gradNum].V[3] = vFHPhilips; + isPhilipsADC = false; + //printMessage(" DimensionIndexValues grad %d b=%g vec=%gx%gx%g\n", gradNum, B0Philips, vRLPhilips, vAPPhilips, vFHPhilips); + + + /*if ((sliceNumberMrPhilips == locationsInAcquisitionPhilips) || ((locationsInAcquisitionPhilips == 0) &&(sliceNumberMrPhilips == maxInStackPositionNumber)) ) { + for (int k = 0; k < 4; k++) + patientPositionEndPhilips[k] = patientPosition[k]; + }*/ + //!!! 16032018 : next line as well as definition of B0Philips may need to be set to zero if Philips omits DiffusionBValue tag for B=0 + B0Philips = -1.0; //Philips may skip reporting B-values for B=0 volumes, so zero these + vRLPhilips = 0.0; + vAPPhilips = 0.0; + vFHPhilips = 0.0; + MRImageGradientOrientationNumber = 0; + break; } break; - case kPhotometricInterpretation: + case kPhotometricInterpretation: { char interp[kDICOMStr]; dcmStr (lLength, &buffer[lPos], interp); if (strcmp(interp, "PALETTE_COLOR") == 0) printError("Photometric Interpretation 'PALETTE COLOR' not supported\n"); - - break; + break; } case kPlanarRGB: d.isPlanarRGB = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); break; case kDim3: d.xyzDim[3] = dcmStrInt(lLength, &buffer[lPos]); + numberOfFrames = d.xyzDim[3]; break; case kSamplesPerPixel: d.samplesPerPixel = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); @@ -3680,20 +3949,16 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); case kTR : d.TR = dcmStrFloat(lLength, &buffer[lPos]); break; - case kTE : { - double te = dcmStrFloat(lLength, &buffer[lPos]); + case kTE : + TE = dcmStrFloat(lLength, &buffer[lPos]); if (d.TE <= 0.0) - d.TE = te; - else if (d.TE != te) - d.TE2 = te; - break; } + d.TE = TE; + break; case kEffectiveTE : { - double te = dcmFloatDouble(lLength, &buffer[lPos], d.isLittleEndian); + TE = dcmFloatDouble(lLength, &buffer[lPos], d.isLittleEndian); if (d.TE <= 0.0) - d.TE = te; - else if (d.TE != te) - d.TE2 = te; - break; } + d.TE = TE; + break; } case kTI : d.TI = dcmStrFloat(lLength, &buffer[lPos]); break; @@ -3765,7 +4030,6 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); if ((lLength == 4) && (d.manufacturer == kMANUFACTURER_PHILIPS)) d.intenScalePhilips = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); break; - case kIntercept : d.intenIntercept = dcmStrFloat(lLength, &buffer[lPos]); break; @@ -3814,7 +4078,9 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); break; }*/ case kNumberOfDynamicScans: - d.numberOfDynamicScans = dcmStrInt(lLength, &buffer[lPos]); + //~d.numberOfDynamicScans = dcmStrInt(lLength, &buffer[lPos]); + numberOfDynamicScans = dcmStrInt(lLength, &buffer[lPos]); + break; case kMRAcquisitionType: //detect 3D acquisition: we can reorient these without worrying about slice time correct or BVEC/BVAL orientation if (lLength > 1) d.is2DAcq = (buffer[lPos]=='2') && (toupper(buffer[lPos+1]) == 'D'); @@ -3873,6 +4139,10 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); else d.sliceOrient = kSliceOrientTra; //transverse (axial) break; } + case kDiffusionBFactor : + if (d.manufacturer != kMANUFACTURER_PHILIPS) break; + B0Philips = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); + break; // case kDiffusionBFactor: // 2001,1003 // if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) { // d.CSA.numDti++; //increment with BFactor: on Philips slices with B=0 have B-factor but no diffusion directions @@ -3900,9 +4170,9 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); // value *unset*! // GE started using this tag in 27, and annoyingly, NOT including // the b value if it is 0 for the slice. - if((d.manufacturer != kMANUFACTURER_PHILIPS) || !is2005140FSQ){ - // d.CSA.numDti++; + //if((d.manufacturer != kMANUFACTURER_PHILIPS) || !is2005140FSQ){ + // d.CSA.numDti++; // if (d.CSA.numDti == 2) { //First time we know that this is a 4D DTI dataset // //d.dti4D = (TDTI *)malloc(kMaxDTI4D * sizeof(TDTI)); // dti4D->S[0].V[0] = d.CSA.dtiV[0]; @@ -3910,12 +4180,14 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); // dti4D->S[0].V[2] = d.CSA.dtiV[2]; // dti4D->S[0].V[3] = d.CSA.dtiV[3]; // } + B0Philips = dcmFloatDouble(lLength, &buffer[lPos], d.isLittleEndian); + //d.CSA.dtiV[0] = dcmFloatDouble(lLength, &buffer[lPos], d.isLittleEndian); set_bVal(&volDiffusion, dcmFloatDouble(lLength, &buffer[lPos], d.isLittleEndian)); // if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D)) // dti4D->S[d.CSA.numDti-1].V[0] = d.CSA.dtiV[0]; - } + //} break; case kDiffusionOrientation: // 0018, 9089 // Note that this is ahead of kPatientPosition (0020,0032), so @@ -3927,9 +4199,19 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); // if (((d.manufacturer == kMANUFACTURER_SIEMENS) || // ((d.manufacturer == kMANUFACTURER_PHILIPS) && !is2005140FSQ)) && // (isAtFirstPatientPosition || isnan(d.patientPosition[1]))) - if((d.manufacturer == kMANUFACTURER_SIEMENS) || - ((d.manufacturer == kMANUFACTURER_PHILIPS) && !is2005140FSQ)) + + //if((d.manufacturer == kMANUFACTURER_SIEMENS) || ((d.manufacturer == kMANUFACTURER_PHILIPS) && !is2005140FSQ)) + if((d.manufacturer == kMANUFACTURER_SIEMENS) || (d.manufacturer == kMANUFACTURER_PHILIPS)) { + float v[4]; + //dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, v); + //dcmMultiFloatDouble(lLength, &buffer[lPos], 3, v, d.isLittleEndian); + dcmMultiFloatDouble(lLength, &buffer[lPos], 3, v, d.isLittleEndian); + vRLPhilips = v[0]; + vAPPhilips = v[1]; + vFHPhilips = v[2]; + set_orientation0018_9089(&volDiffusion, lLength, &buffer[lPos], d.isLittleEndian); + } break; // case kSharedFunctionalGroupsSequence: // if ((d.manufacturer == kMANUFACTURER_SIEMENS) && isAtFirstPatientPosition) { @@ -3937,20 +4219,41 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); // // part of buffer[lPos]. // } // break; + + //case kSliceNumberMrPhilips : + // sliceNumberMrPhilips3D = dcmStrInt(lLength, &buffer[lPos]); + // break; case kSliceNumberMrPhilips : { if (d.manufacturer != kMANUFACTURER_PHILIPS) break; - int sliceNumber; - sliceNumber = dcmStrInt(lLength, &buffer[lPos]); - printError("====> %d %d %d\n", locationsInAcquisitionPhilips, d.patientPositionNumPhilips, sliceNumber); + + sliceNumberMrPhilips = dcmStrInt(lLength, &buffer[lPos]); + int sliceNumber = sliceNumberMrPhilips; + /*//if (sliceNumberMrPhilips < 10) printMessage("%d ====> %d\n", sliceNumberMrPhilips, sliceNumber); + //printMessage("====> sliceNumberMrPhilips\t%d\t%d\t%d\n", locationsInAcquisitionPhilips, d.patientPositionNumPhilips, sliceNumber); if ((sliceNumberMrPhilips >= locationsInAcquisitionPhilips) && (sliceNumberMrPhilips < (2*locationsInAcquisitionPhilips)) ) { int p = sliceNumberMrPhilips - locationsInAcquisitionPhilips; if (p < kMaxDTI4D) dti4D->S[p].sliceNumberMrPhilipsVol2 = sliceNumber; - if ((p > 0) && (abs(dti4D->S[p].sliceNumberMrPhilipsVol2 - dti4D->S[p -1].sliceNumberMrPhilipsVol2) > 1) ) - d.isSlicesSpatiallySequentialPhilips = false; + //~ if ((p > 0) && (abs(dti4D->S[p].sliceNumberMrPhilipsVol2 - dti4D->S[p -1].sliceNumberMrPhilipsVol2) > 1) ) + //~ d.isSlicesSpatiallySequentialPhilips = false; + }*/ + //use public patientPosition if it exists - fall back to private patient position + if ((sliceNumber == 1) && (!isnan(patientPosition[1])) ) { + for (int k = 0; k < 4; k++) + patientPositionStartPhilips[k] = patientPosition[k]; + } else if ((sliceNumber == 1) && (!isnan(patientPositionPrivate[1]))) { + for (int k = 0; k < 4; k++) + patientPositionStartPhilips[k] = patientPositionPrivate[k]; + } + if ((sliceNumber == locationsInAcquisitionPhilips) && (!isnan(patientPosition[1])) ) { + for (int k = 0; k < 4; k++) + patientPositionEndPhilips[k] = patientPosition[k]; + } else if ((sliceNumber == locationsInAcquisitionPhilips) && (!isnan(patientPositionPrivate[1])) ) { + for (int k = 0; k < 4; k++) + patientPositionEndPhilips[k] = patientPositionPrivate[k]; } - sliceNumberMrPhilips ++; + /*sliceNumberMrPhilips ++; if ((locationsInAcquisitionPhilips > 0) && (d.patientPositionNumPhilips == locationsInAcquisitionPhilips)) break; //we have acquired all slices in volume (e.g. all volumes after 1st for XYZT storage @@ -3961,25 +4264,32 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); break; } dti4D->S[d.patientPositionNumPhilips].sliceNumberMrPhilips = sliceNumber; - if ((d.patientPositionNumPhilips > 0) && (abs(dti4D->S[d.patientPositionNumPhilips].sliceNumberMrPhilips - dti4D->S[d.patientPositionNumPhilips -1].sliceNumberMrPhilips) > 1) ) - d.isSlicesSpatiallySequentialPhilips = false; //slices not sequential (1,2,3,4 or 4,3,2,1) but 4,3,1,2 + //~ if ((d.patientPositionNumPhilips > 0) && (abs(dti4D->S[d.patientPositionNumPhilips].sliceNumberMrPhilips - dti4D->S[d.patientPositionNumPhilips -1].sliceNumberMrPhilips) > 1) ) + //~ d.isSlicesSpatiallySequentialPhilips = false; //slices not sequential (1,2,3,4 or 4,3,2,1) but 4,3,1,2 d.patientPositionNumPhilips++; //Philips can save 3D acquisitions in a single file with slices stored in non-sequential order. We need to know the first and final spatial position - if (sliceNumber == 1) { - for (int k = 0; k < 4; k++) - patientPositionStartPhilips[k] = patientPosition[k]; - } - if (sliceNumber == locationsInAcquisitionPhilips) { - for (int k = 0; k < 4; k++) - patientPositionEndPhilips[k] = patientPosition[k]; - } + //printMessage("xxxx====> - %d\n", sliceNumber); if (isVerbose > 1) printMessage("slice %d is spatial position %d\n", d.patientPositionNumPhilips, sliceNumber); + */ break; } case kNumberOfSlicesMrPhilips : if (d.manufacturer != kMANUFACTURER_PHILIPS) break; locationsInAcquisitionPhilips = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); + //printMessage("====> locationsInAcquisitionPhilips\t%d\n", locationsInAcquisitionPhilips); + break; + case kDiffusionDirectionRL: + if (d.manufacturer != kMANUFACTURER_PHILIPS) break; + vRLPhilips = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); + break; + case kDiffusionDirectionAP: + if (d.manufacturer != kMANUFACTURER_PHILIPS) break; + vAPPhilips = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); + break; + case kDiffusionDirectionFH: + if (d.manufacturer != kMANUFACTURER_PHILIPS) break; + vFHPhilips = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); break; // case kDiffusionDirectionRL: // if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) { @@ -4010,6 +4320,18 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); // d.CSA.dtiV[d.CSA.numDti-1][3] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);*/ // //http://www.na-mic.org/Wiki/index.php/NAMIC_Wiki:DTI:DICOM_for_DWI_and_DTI // break; + //~~ + case kPrivatePerFrameSq : + if (d.manufacturer != kMANUFACTURER_PHILIPS) break; + //if ((vr[0] == 'S') && (vr[1] == 'Q')) break; + //if (!is2005140FSQwarned) + // printWarning("expected VR of 2005,140F to be 'SQ' (prior DICOM->DICOM conversion error?)\n"); + is2005140FSQ = true; + //is2005140FSQwarned = true; + case kMRImageGradientOrientationNumber : + if (d.manufacturer == kMANUFACTURER_PHILIPS) + MRImageGradientOrientationNumber = dcmStrInt(lLength, &buffer[lPos]); + break; case kWaveformSq: d.imageStart = 1; //abort!!! printMessage("Skipping DICOM (audio not image) '%s'\n", fname); @@ -4056,22 +4378,23 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); case kOrientationACR : //use in emergency if kOrientation is not present! if (!isOrient) dcmMultiFloat(lLength, (char*)&buffer[lPos], 6, d.orient); break; + case kTemporalPositionIdentifier : + temporalPositionIdentifier = dcmStrInt(lLength, &buffer[lPos]); + break; case kOrientation : { if (isOrient) { //already read orient - read for this slice to see if it varies (localizer) float orient[7]; dcmMultiFloat(lLength, (char*)&buffer[lPos], 6, orient); if ((!isSameFloatGE(d.orient[1], orient[1]) || !isSameFloatGE(d.orient[2], orient[2]) || !isSameFloatGE(d.orient[3], orient[3]) || !isSameFloatGE(d.orient[4], orient[4]) || !isSameFloatGE(d.orient[5], orient[5]) || !isSameFloatGE(d.orient[6], orient[6]) ) ) { - d.isLocalizer = true; - if ((isVerbose > 1)) - printMessage("slices not stacked: orientation varies (localizer?) [%g %g %g %g %g %g] != [%g %g %g %g %g %g]\n", + if (!d.isLocalizer) + printMessage("slice orientation varies (localizer?) [%g %g %g %g %g %g] != [%g %g %g %g %g %g]\n", d.orient[1], d.orient[2], d.orient[3],d.orient[4], d.orient[5], d.orient[6], orient[1], orient[2], orient[3],orient[4], orient[5], orient[6]); + d.isLocalizer = true; } - } dcmMultiFloat(lLength, (char*)&buffer[lPos], 6, d.orient); - isOrient = true; break; } case kImagesInAcquisition : @@ -4131,7 +4454,8 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 ); } printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\t%s\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos, tagStr); } else - printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tsq=%d %d %d %c%c\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos, sqDepth, sqDepthPrivate, sqEndPrivate, vr[0], vr[1]); + printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos); + //printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tsq=%d %d %d %c%c\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos, sqDepth, sqDepthPrivate, sqEndPrivate, vr[0], vr[1]); //if (d.isExplicitVR) printMessage(" VR=%c%c\n", vr[0], vr[1]); } //printMessage(" tag=%04x,%04x length=%u pos=%ld %c%c nest=%d\n", groupElement & 65535,groupElement>>16, lLength, lPos,vr[0], vr[1], nest); lPos = lPos + (lLength); @@ -4234,33 +4558,143 @@ if (d.isHasPhase) printError("=====> reps %d %d\n", d.patientPositionRepeats, d.patientPositionSequentialRepeats ); */ - if ((d.patientPositionSequentialRepeats > 1) && ( (d.xyzDim[3] % d.patientPositionSequentialRepeats) == 0 )) { + /*if ((patientPositionSequentialRepeats > 1) && ( (d.xyzDim[3] % patientPositionSequentialRepeats) == 0 )) { //will require Converting XYTZ to XYZT - d.numberOfDynamicScans = d.xyzDim[3] / d.patientPositionSequentialRepeats; - d.xyzDim[4] = d.xyzDim[3] / d.numberOfDynamicScans; + //~ d.numberOfDynamicScans = d.xyzDim[3] / d.patientPositionSequentialRepeats; + //~ d.xyzDim[4] = d.xyzDim[3] / d.numberOfDynamicScans; + numberOfDynamicScans = d.xyzDim[3] / patientPositionSequentialRepeats; + d.xyzDim[4] = d.xyzDim[3] / numberOfDynamicScans; + d.xyzDim[3] = d.numberOfDynamicScans; + }*/ + if (numberOfFrames == 0) numberOfFrames = d.xyzDim[3]; + + if ((numberOfDynamicScans > 1) && (d.xyzDim[4] < 2) && (d.xyzDim[3] > 1) && ((d.xyzDim[3] % numberOfDynamicScans) == 0)) { + d.xyzDim[3] = d.xyzDim[3] / numberOfDynamicScans; + d.xyzDim[4] = numberOfDynamicScans; } - if ( (!d.isSlicesSpatiallySequentialPhilips) && (!isnan(patientPositionStartPhilips[1])) && (!isnan(patientPositionEndPhilips[1]))) { - for (int k = 0; k < 4; k++) { - d.patientPosition[k] = patientPositionStartPhilips[k]; - d.patientPositionLast[k] = patientPositionEndPhilips[k]; - } - if (d.numberOfDynamicScans > 1) + if ((maxInStackPositionNumber > 1) && (d.xyzDim[4] < 2) && (d.xyzDim[3] > 1) && ((d.xyzDim[3] % maxInStackPositionNumber) == 0)) { + d.xyzDim[4] = d.xyzDim[3] / maxInStackPositionNumber; + d.xyzDim[3] = maxInStackPositionNumber; + } + + if ((!isnan(patientPositionStartPhilips[1])) && (!isnan(patientPositionEndPhilips[1]))) { + for (int k = 0; k < 4; k++) { + d.patientPosition[k] = patientPositionStartPhilips[k]; + d.patientPositionLast[k] = patientPositionEndPhilips[k]; + } + /*if (d.numberOfDynamicScans > 1) printError("4D data with non-contiguous slices: please check output [newest feature]\n"); else - printMessage("Slices not spatially contiguous: please check output [new feature]\n"); + printMessage("Slices not spatially contiguous: please check output [new feature]\n");*/ } + if (!isnan(patientPositionStartPhilips[1])) //for Philips data without + for (int k = 0; k < 4; k++) + d.patientPosition[k] = patientPositionStartPhilips[k]; + if (isVerbose) { - printMessage("%s\n patient position (0020,0032)\t%g\t%g\t%g\n",fname, d.patientPosition[1],d.patientPosition[2],d.patientPosition[3]); - printMessage("%s\n orient (0020,0037)\t%g\t%g\t%g\t%g\t%g\t%g\n",fname, d.orient[1],d.orient[2],d.orient[3], d.orient[4],d.orient[5],d.orient[6]); - printMessage(" acq %d img %d ser %ld dim %dx%dx%d mm %gx%gx%g offset %d dyn %d loc %d valid %d ph %d mag %d posReps %d nDTI %d 3d %d bits %d littleEndian %d echo %d coil %d TE %g TR %g\n",d.acquNum,d.imageNum,d.seriesNum,d.xyzDim[1],d.xyzDim[2],d.xyzDim[3],d.xyzMM[1],d.xyzMM[2],d.xyzMM[3],d.imageStart, d.numberOfDynamicScans, d.locationsInAcquisition, d.isValid, d.isHasPhase, d.isHasMagnitude,d.patientPositionSequentialRepeats, d.CSA.numDti, d.is3DAcq, d.bitsAllocated, d.isLittleEndian, d.echoNum, d.coilNum, d.TE, d.TR); - if (d.CSA.dtiV[0] > 0) - printMessage(" DWI bxyz %g %g %g %g\n", d.CSA.dtiV[0], d.CSA.dtiV[1], d.CSA.dtiV[2], d.CSA.dtiV[3]); + printMessage("DICOM file %s:\n", fname); + printMessage(" patient position (0020,0032)\t%g\t%g\t%g\n", d.patientPosition[1],d.patientPosition[2],d.patientPosition[3]); + if (!isnan(patientPositionEndPhilips[1])) + printMessage(" patient position end (0020,0032)\t%g\t%g\t%g\n", patientPositionEndPhilips[1],patientPositionEndPhilips[2],patientPositionEndPhilips[3]); + printMessage(" orient (0020,0037)\t%g\t%g\t%g\t%g\t%g\t%g\n", d.orient[1],d.orient[2],d.orient[3], d.orient[4],d.orient[5],d.orient[6]); + printMessage(" acq %d img %d ser %ld dim %dx%dx%dx%d mm %gx%gx%g offset %d loc %d valid %d ph %d mag %d nDTI %d 3d %d bits %d littleEndian %d echo %d coil %d TE %g TR %g\n",d.acquNum,d.imageNum,d.seriesNum,d.xyzDim[1],d.xyzDim[2],d.xyzDim[3], d.xyzDim[4],d.xyzMM[1],d.xyzMM[2],d.xyzMM[3],d.imageStart, d.locationsInAcquisition, d.isValid, d.isHasPhase, d.isHasMagnitude, d.CSA.numDti, d.is3DAcq, d.bitsAllocated, d.isLittleEndian, d.echoNum, d.coilNum, d.TE, d.TR); + //if (d.CSA.dtiV[0] > 0) + // printMessage(" DWI bxyz %g %g %g %g\n", d.CSA.dtiV[0], d.CSA.dtiV[1], d.CSA.dtiV[2], d.CSA.dtiV[3]); } - if (d.patientPositionNumPhilips >= kMaxDTI4D) { + + if ((numDimensionIndexValues > 1) && (numDimensionIndexValues == numberOfFrames)) { + //Philips enhanced datasets can have custom slice orders and pack images with different TE, Phase/Magnitude/Etc. + if (isVerbose > 1) { // + int mn[MAX_NUMBER_OF_DIMENSIONS]; + int mx[MAX_NUMBER_OF_DIMENSIONS]; + for (int j = 0; j < MAX_NUMBER_OF_DIMENSIONS; j++) { + mx[j] = dcmDim[0].dimIdx[j]; + mn[j] = mx[j]; + for (int i = 0; i < numDimensionIndexValues; i++) { + if (mx[j] < dcmDim[i].dimIdx[j]) mx[j] = dcmDim[i].dimIdx[j]; + if (mn[j] > dcmDim[i].dimIdx[j]) mn[j] = dcmDim[i].dimIdx[j]; + } + } + printMessage(" DimensionIndexValues (0020,9157), dimensions with variability:\n"); + for (int i = 0; i < MAX_NUMBER_OF_DIMENSIONS; i++) + if (mn[i] != mx[i]) + printMessage(" Dimension %d Range: %d..%d\n", i, mn[i], mx[i]); + } //verbose > 1 + qsort(dcmDim, numberOfFrames, sizeof(struct TDCMdim), compareTDCMdim); + //for (int i = 0; i < numberOfFrames; i++) + // printf("%d -> %d %d %d %d\n", i, dcmDim[i].diskPos, dcmDim[i].dimIdx[1], dcmDim[i].dimIdx[2], dcmDim[i].dimIdx[3]); + + for (int i = 0; i < numberOfFrames; i++) + dti4D->sliceOrder[i] = dcmDim[i].diskPos; + if ((d.xyzDim[4] > 1) && (d.xyzDim[4] < kMaxDTI4D)) { //record variations in TE + d.isScaleOrTEVaries = false; + bool isTEvaries = false; + bool isScaleVaries = false; + //setting j = 1 in next few lines is a hack, just in case TE/scale/intercept listed AFTER dimensionIndexValues + int j = 0; + if (d.xyzDim[3] > 1) j = 1; + for (int i = 0; i < d.xyzDim[4]; i++) { + dti4D->TE[i] = dcmDim[j+(i * d.xyzDim[3])].TE; + dti4D->intenScale[i] = dcmDim[j+(i * d.xyzDim[3])].intenScale; + dti4D->intenIntercept[i] = dcmDim[j+(i * d.xyzDim[3])].intenIntercept; + dti4D->isPhase[i] = dcmDim[j+(i * d.xyzDim[3])].isPhase; + dti4D->intenScalePhilips[i] = dcmDim[j+(i * d.xyzDim[3])].intenScalePhilips; + if (dti4D->TE[i] != d.TE) isTEvaries = true; + if (dti4D->TE[i] != d.TE) ; + if (dti4D->intenScale[i] != d.intenScale) isScaleVaries = true; + if (dti4D->intenIntercept[i] != d.intenIntercept) isScaleVaries = true; + if (dti4D->isPhase[i] != isPhase) d.isScaleOrTEVaries = true; + } + if((isScaleVaries) || (isTEvaries)) d.isScaleOrTEVaries = true; + //if echoVaries,count number of echoes + /*int echoNum = 1; + for (int i = 1; i < d.xyzDim[4]; i++) { + if (dti4D->TE[i-1] != dti4D->TE[i]) + }*/ + if ((isVerbose) && (d.isScaleOrTEVaries)) { + printMessage("Parameters vary across 3D volumes packed in single DICOM file:\n"); + for (int i = 0; i < d.xyzDim[4]; i++) + printMessage(" %d TE=%g Slope=%g Inter=%g PhilipsScale=%g Phase=%d\n", i, dti4D->TE[i], dti4D->intenScale[i], dti4D->intenIntercept[i], dti4D->intenScalePhilips[i], dti4D->isPhase[i] ); + } + } + } + /* //Attempt to append ADC + printMessage("CXC grad %g %d %d\n", philDTI[0].V[0], maxGradNum, d.xyzDim[4]); + if ((maxGradNum > 1) && ((maxGradNum+1) == d.xyzDim[4]) ) { + //ADC map (non-zero b-value with zero vector length) + if (isVerbose) + printMessage("Final volume does not have an associated 0020,9157. Assuming final volume is an ADC map\n", philDTI[0].V[0], maxGradNum, d.xyzDim[4]); + philDTI[maxGradNum].V[0] = 1000.0; + philDTI[maxGradNum].V[1] = 0.0; + philDTI[maxGradNum].V[2] = 0.0; + philDTI[maxGradNum].V[3] = 0.0; + maxGradNum++; + }*/ + + + //printMessage("CXC grad %g %d %d\n", philDTI[0].V[0], maxGradNum, d.xyzDim[4]); + + if ((philDTI[0].V[0] >= 0) && (maxGradNum == d.xyzDim[4])) { + if (isVerbose) + printMessage("Using diffusion data coded by DimensionIndexValues\n"); + for (int i = 0; i < d.xyzDim[4]; i++) { + dti4D->S[i].V[0] = philDTI[i].V[0]; + dti4D->S[i].V[1] = philDTI[i].V[1]; + dti4D->S[i].V[2] = philDTI[i].V[2]; + dti4D->S[i].V[3] = philDTI[i].V[3]; + if (isVerbose > 1) + printMessage(" grad %d b=%g vec=%gx%gx%g\n", i, dti4D->S[i].V[0], dti4D->S[i].V[1], dti4D->S[i].V[2], dti4D->S[i].V[3]); + + } + d.CSA.numDti = maxGradNum; + } + //~~ + /*if (d.patientPositionNumPhilips >= kMaxDTI4D) { + printError("Too many 2D slices in a single file [recompile with increased kMaxDTI4D] detected=%d, max = %d\n", d.patientPositionNumPhilips, kMaxDTI4D); d.CSA.numDti = 0; - } + }*/ if (d.CSA.numDti >= kMaxDTI4D) { printError("Unable to convert DTI [recompile with increased kMaxDTI4D] detected=%d, max = %d\n", d.CSA.numDti, kMaxDTI4D); d.CSA.numDti = 0; diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 08ef2a81..4ee7e6fe 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -42,6 +42,8 @@ extern "C" { static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic static const int kMaxDTI4D = 4096; //maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images +static const int kMaxSlice2D = 16786; //maximum number of 2D slices in 4D (Philips) images + #define kDICOMStr 64 #define kDICOMStrLarge 256 #define kMANUFACTURER_UNKNOWN 0 @@ -76,11 +78,13 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; struct TDTI { float V[4]; - int sliceNumberMrPhilips; - int sliceNumberMrPhilipsVol2; + //int totalSlicesIn4DOrder; }; struct TDTI4D { struct TDTI S[kMaxDTI4D]; + int sliceOrder[kMaxSlice2D]; // [7,3,2] means the first slice on disk should be moved to 7th position + float TE[kMaxDTI4D], intenScale[kMaxDTI4D], intenIntercept[kMaxDTI4D], intenScalePhilips[kMaxDTI4D]; + bool isPhase[kMaxDTI4D]; }; #ifdef _MSC_VER //Microsoft nomenclature for packed structures is different... @@ -111,15 +115,18 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; } TCSAitem; //Siemens csa item structure #endif struct TCSAdata { - bool isPhaseMap; - float sliceTiming[kMaxEPI3D], dtiV[4], sliceNormV[4], bandwidthPerPixelPhaseEncode, sliceMeasurementDuration; + float sliceTiming[kMaxEPI3D], dtiV[4], sliceNormV[4], bandwidthPerPixelPhaseEncode, sliceMeasurementDuration; int numDti, SeriesHeader_offset, SeriesHeader_length, multiBandFactor, sliceOrder, slice_start, slice_end, mosaicSlices,protocolSliceNumber1,phaseEncodingDirectionPositive; + bool isPhaseMap; + }; struct TDICOMdata { long seriesNum; int xyzDim[5]; - int modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, echoTrainLength, patientPositionNumPhilips, coilNum, echoNum, sliceOrient,numberOfDynamicScans, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,patientPositionSequentialRepeats,patientPositionRepeats,locationsInAcquisition, compressionScheme; - float patientWeight, zSpacing, zThick, pixelBandwidth, SAR, phaseFieldofView, accelFactPE, flipAngle, fieldStrength, TE, TE2, TI, TR, intenScale, intenIntercept, intenScalePhilips, gantryTilt, lastScanLoc, angulation[4]; + //numberOfDynamicScans, patientPositionNumPhilips + //patientPositionSequentialRepeats,patientPositionRepeats, + int modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, echoTrainLength, coilNum, echoNum, sliceOrient, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,locationsInAcquisition, compressionScheme; + float patientWeight, zSpacing, zThick, pixelBandwidth, SAR, phaseFieldofView, accelFactPE, flipAngle, fieldStrength, TE, TI, TR, intenScale, intenIntercept, intenScalePhilips, gantryTilt, lastScanLoc, angulation[4]; float orient[7], patientPosition[4], patientPositionLast[4], xyzMM[4], stackOffcentre[4]; float radionuclidePositronFraction, radionuclideTotalDose, radionuclideHalfLife, doseCalibrationFactor; //PET ISOTOPE MODULE ATTRIBUTES (C.8-57) float ecat_isotope_halflife, ecat_dosage; @@ -129,8 +136,10 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; char imageComments[kDICOMStrLarge]; uint32_t dimensionIndexValues[MAX_NUMBER_OF_DIMENSIONS]; struct TCSAdata CSA; - bool isSegamiOasis, isDerived, isXRay, isMultiEcho, isSlicesSpatiallySequentialPhilips, isValid, is3DAcq, is2DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase,isHasMagnitude,isHasMixed, isFloat, isResampled, isLocalizer; + //isSlicesSpatiallySequentialPhilips + bool isSegamiOasis, isScaleOrTEVaries, isDerived, isXRay, isMultiEcho, isValid, is3DAcq, is2DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase,isHasMagnitude,isHasMixed, isFloat, isResampled, isLocalizer; char phaseEncodingRC, patientSex; + //uint32_t *totalSlicesIn4DOrder; //Reordering array for Philips slices }; size_t nii_ImgBytes(struct nifti_1_header hdr); @@ -141,14 +150,14 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D *dti4D); unsigned char * nii_flipY(unsigned char* bImg, struct nifti_1_header *h); unsigned char * nii_flipZ(unsigned char* bImg, struct nifti_1_header *h); - unsigned char * nii_reorderSlices(unsigned char* bImg, struct nifti_1_header *h, struct TDTI4D *dti4D); + //*unsigned char * nii_reorderSlices(unsigned char* bImg, struct nifti_1_header *h, struct TDTI4D *dti4D); void changeExt (char *file_name, const char* ext); unsigned char * nii_planar2rgb(unsigned char* bImg, struct nifti_1_header *hdr, int isPlanar); int isDICOMfile(const char * fname); //0=not DICOM, 1=DICOM, 2=NOTSURE(not part 10 compliant) void setQSForm(struct nifti_1_header *h, mat44 Q44i, bool isVerbose); int headerDcm2Nii2(struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1_header *h, int isVerbose); int headerDcm2Nii(struct TDICOMdata d, struct nifti_1_header *h, bool isComputeSForm) ; - unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct TDICOMdata dcm, bool iVaries, int compressFlag, int isVerbose); + unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct TDICOMdata dcm, bool iVaries, int compressFlag, int isVerbose, struct TDTI4D *dti4D); #ifdef __cplusplus } #endif diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index ff22c9ce..fd269c12 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -311,6 +311,8 @@ void siemensPhilipsCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI printWarning("Saving %d DTI gradients. Validate vectors (matrix had a negative determinant).\n", d->CSA.numDti); //perhaps Siemens sagittal } else if ( d->sliceOrient == kSliceOrientTra) { printMessage("Saving %d DTI gradients. Validate vectors.\n", d->CSA.numDti); + } else if ( d->sliceOrient == kSliceOrientUnknown) { + printWarning("Saving %d DTI gradients. Validate vectors (image slice orientation not reported, e.g. 2001,100B).\n", d->CSA.numDti); } else { printWarning("Saving %d DTI gradients. Validate vectors (images are not axial slices).\n", d->CSA.numDti); } @@ -783,7 +785,7 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, json_Float(fp, "\t\"SAR\": %g,\n", d.SAR ); if (d.echoNum > 1) fprintf(fp, "\t\"EchoNumber\": %d,\n", d.echoNum); if ((d.TE > 0.0) && (!d.isXRay)) fprintf(fp, "\t\"EchoTime\": %g,\n", d.TE / 1000.0 ); - if ((d.TE2 > 0.0) && (!d.isXRay)) fprintf(fp, "\t\"EchoTime2\": %g,\n", d.TE2 / 1000.0 ); + //if ((d.TE2 > 0.0) && (!d.isXRay)) fprintf(fp, "\t\"EchoTime2\": %g,\n", d.TE2 / 1000.0 ); json_Float(fp, "\t\"RepetitionTime\": %g,\n", d.TR / 1000.0 ); json_Float(fp, "\t\"InversionTime\": %g,\n", d.TI / 1000.0 ); json_Float(fp, "\t\"FlipAngle\": %g,\n", d.flipAngle ); @@ -1131,35 +1133,35 @@ int * nii_SaveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],str bvals[i] = bvals[i] + (0.5 * i/numDti); //add a small bias so ties are kept in sequential order } if (*numADC > 0) { - // DWIs (i.e. short diffusion scans with too few directions to - // calculate tensors...they typically acquire b=0 + 3 b > 0 so - // the isotropic trace or MD can be calculated) often come as - // b=0 and trace pairs, with the b=0 and trace in either order, - // and often as "ORIGINAL", even though the trace is not. - // The bval file is needed for downstream processing to know - // * which is the b=0 and which is the trace, and - // * what b is for the trace, - // so dcm2niix should *always* write the bval and bvec files, - // AND include the b for the trace for DWIs. - // One hackish way to accomplish that is to set *numADC = 0 - // when *numADC == 1 && numDti == 2. - // - Rob Reid, 2017-11-29. - if ((*numADC == 1) && ((numDti - *numADC) < 2)){ - *numADC = 0; - printMessage("Note: this appears to be a b=0+trace DWI; ADC/trace removal has been disabled.\n"); - } - else{ - if ((numDti - *numADC) < 2) { - if (!dcmList[indx0].isDerived) //no need to warn if images are derived Trace/ND pair - printWarning("No bvec/bval files created: only single value after ADC excluded\n"); - *numADC = 0; - free(bvals); - free(vx); - return NULL; - } - printMessage("Note: %d volumes appear to be ADC or trace images that will be removed to allow processing\n", - *numADC); - } + // DWIs (i.e. short diffusion scans with too few directions to + // calculate tensors...they typically acquire b=0 + 3 b > 0 so + // the isotropic trace or MD can be calculated) often come as + // b=0 and trace pairs, with the b=0 and trace in either order, + // and often as "ORIGINAL", even though the trace is not. + // The bval file is needed for downstream processing to know + // * which is the b=0 and which is the trace, and + // * what b is for the trace, + // so dcm2niix should *always* write the bval and bvec files, + // AND include the b for the trace for DWIs. + // One hackish way to accomplish that is to set *numADC = 0 + // when *numADC == 1 && numDti == 2. + // - Rob Reid, 2017-11-29. + if ((*numADC == 1) && ((numDti - *numADC) < 2)){ + *numADC = 0; + printMessage("Note: this appears to be a b=0+trace DWI; ADC/trace removal has been disabled.\n"); + } + else{ + if ((numDti - *numADC) < 2) { + if (!dcmList[indx0].isDerived) //no need to warn if images are derived Trace/ND pair + printWarning("No bvec/bval files created: only single value after ADC excluded\n"); + *numADC = 0; + free(bvals); + free(vx); + return NULL; + } + printMessage("Note: %d volumes appear to be ADC or trace images that will be removed to allow processing\n", + *numADC); + } } //sort ALL including ADC int * volOrderIndex = (int *) malloc(numDti * sizeof(int)); @@ -1200,7 +1202,6 @@ int * nii_SaveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],str if (!isSequential) printMessage("DTI volumes re-ordered by ascending b-value\n"); dcmList[indx0].CSA.numDti = numDti; //warning structure not changed outside scope! - geCorrectBvecs(&dcmList[indx0],sliceDir, vx); siemensPhilipsCorrectBvecs(&dcmList[indx0],sliceDir, vx); if (!opts.isFlipY ) { //!FLIP_Y&& (dcmList[indx0].CSA.mosaicSlices < 2) mosaics are always flipped in the Y direction @@ -2332,7 +2333,7 @@ void checkSliceTiming(struct TDICOMdata * d, struct TDICOMdata * d1) { printMessage("CSA slice timing based on 2nd volume, 1st volume corrupted (CMRR bug, range %g..%g, TR=%gms)\n", minT, maxT, d->TR); }//checkSliceTiming -int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmList[], struct TSearchList *nameList, struct TDCMopts opts, struct TDTI4D *dti4D) { +int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmList[], struct TSearchList *nameList, struct TDCMopts opts, struct TDTI4D *dti4D, int vol) { bool iVaries = intensityScaleVaries(nConvert,dcmSort,dcmList); float *sliceMMarray = NULL; //only used if slices are not equidistant uint64_t indx = dcmSort[0].indx; @@ -2353,7 +2354,7 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis } bool saveAs3D = dcmList[indx].isHasPhase; struct nifti_1_header hdr0; - unsigned char * img = nii_loadImgXL(nameList->str[indx], &hdr0,dcmList[indx], iVaries, opts.compressFlag, opts.isVerbose); + unsigned char * img = nii_loadImgXL(nameList->str[indx], &hdr0,dcmList[indx], iVaries, opts.compressFlag, opts.isVerbose, dti4D); if (strlen(opts.imageComments) > 0) { for (int i = 0; i < 24; i++) hdr0.aux_file[i] = 0; //remove dcm.imageComments snprintf(hdr0.aux_file,24,"%s",opts.imageComments); @@ -2443,6 +2444,7 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis hdr0.dim[5] = nConvert; hdr0.dim[0] = 5; } + /*if (nConvert > 1) { //next determine if TR is true time between volumes double startTime = dcmList[indx0].acquisitionTime; double endTime = startTime; @@ -2461,7 +2463,7 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis for (int i = 1; i < nConvert; i++) { //stack additional images indx = dcmSort[i].indx; //if (headerDcm2Nii(dcmList[indx], &hdrI) == EXIT_FAILURE) return EXIT_FAILURE; - img = nii_loadImgXL(nameList->str[indx], &hdrI, dcmList[indx],iVaries, opts.compressFlag, opts.isVerbose); + img = nii_loadImgXL(nameList->str[indx], &hdrI, dcmList[indx],iVaries, opts.compressFlag, opts.isVerbose, dti4D); if (img == NULL) return EXIT_FAILURE; if ((hdr0.dim[1] != hdrI.dim[1]) || (hdr0.dim[2] != hdrI.dim[2]) || (hdr0.bitpix != hdrI.bitpix)) { printError("Image dimensions differ %s %s",nameList->str[dcmSort[0].indx], nameList->str[indx]); @@ -2473,6 +2475,18 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis free(img); } } + if ((vol >= 0) && (hdr0.dim[4] > 1) && (vol < hdr0.dim[4])) { + size_t imgsz4D = imgsz; + hdr0.dim[0] = 3; //3D + hdr0.dim[4] = 1; + imgsz = nii_ImgBytes(hdr0); + unsigned char *img4D = (unsigned char *)malloc(imgsz4D); + memcpy(&img4D[0], &imgM[0], imgsz4D); + free(imgM); + imgM = (unsigned char *)malloc(imgsz); + memcpy(&imgM[0], &img4D[vol * imgsz], imgsz); + free(img4D); + } char pathoutname[2048] = {""}; if (nii_createFilename(dcmList[dcmSort[0].indx], pathoutname, opts) == EXIT_FAILURE) { free(imgM); @@ -2530,8 +2544,8 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis if ((hdr0.datatype == DT_UINT16) && (!dcmList[dcmSort[0].indx].isSigned)) nii_check16bitUnsigned(imgM, &hdr0); printMessage( "Convert %d DICOM as %s (%dx%dx%dx%d)\n", nConvert, pathoutname, hdr0.dim[1],hdr0.dim[2],hdr0.dim[3],hdr0.dim[4]); PhilipsPrecise(&dcmList[dcmSort[0].indx], opts.isPhilipsFloatNotDisplayScaling, &hdr0); - if (!dcmList[dcmSort[0].indx].isSlicesSpatiallySequentialPhilips) - nii_reorderSlices(imgM, &hdr0, dti4D); + //~ if (!dcmList[dcmSort[0].indx].isSlicesSpatiallySequentialPhilips) + //~ nii_reorderSlices(imgM, &hdr0, dti4D); if (hdr0.dim[3] < 2) printWarning("Check that 2D images are not mirrored.\n"); #ifndef HAVE_R @@ -2604,6 +2618,49 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis return returnCode;//EXIT_SUCCESS; }// saveDcm2Nii() +int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmList[], struct TSearchList *nameList, struct TDCMopts opts, struct TDTI4D *dti4D) { + //this wrapper does nothing if all the images share the same echo time and scale + // however, it segments images when these properties vary + uint64_t indx = dcmSort[0].indx; + if ((!dcmList[indx].isScaleOrTEVaries) || (dcmList[indx].xyzDim[4] < 2)) + return saveDcm2NiiCore(nConvert, dcmSort, dcmList, nameList, opts, dti4D, -1); + if ((dcmList[indx].xyzDim[4]) && (dti4D->sliceOrder[0] < 0)) { + printError("Unexpected error for image with varying echo time or intensity scaling\n"); + return EXIT_FAILURE; + } + + int echoNum[kMaxDTI4D]; + int echo = 1; + for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) + echoNum[i] = 0; + echoNum[0] = 1; + for (int i = 1; i < dcmList[indx].xyzDim[4]; i++) { + for (int j = 0; j < i; j++) + if (dti4D->TE[i] == dti4D->TE[j]) echoNum[i] = echoNum[j]; + if (echoNum[i] == 0) { + echo++; + echoNum[i] = echo; + } + + } + if (echo > 1) dcmList[indx].isMultiEcho = true; + for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) { + //printMessage("volume %d TE=%g Slope=%g Inter=%g Phase=%d\n", i, dti4D->TE[i], dti4D->intenScale[i], dti4D->intenIntercept[i], dti4D->isPhase[i] ); + dcmList[indx].TE = dti4D->TE[i]; + dcmList[indx].intenScale = dti4D->intenScale[i]; + dcmList[indx].intenIntercept = dti4D->intenIntercept[i]; + dcmList[indx].isHasPhase = dti4D->isPhase[i]; + dcmList[indx].intenScalePhilips = dti4D->intenScalePhilips[i]; + dcmList[indx].isHasMagnitude = false; + dcmList[indx].echoNum = echoNum[i]; + /*int echoNum = 1; + for (int i = 1; i < d.xyzDim[4]; i++) { + if (dti4D->TE[i-1] != dti4D->TE[i]) + }*/ + int ret = saveDcm2NiiCore(nConvert, dcmSort, dcmList, nameList, opts, dti4D, i); + } +} + void fillTDCMsort(struct TDCMsort& tdcmref, const uint64_t indx, const struct TDICOMdata& dcmdata){ // Copy the relevant parts of dcmdata to tdcmref. tdcmref.indx = indx; @@ -2794,7 +2851,8 @@ int singleDICOM(struct TDCMopts* opts, char *fname) { dcmList[0].converted2NII = 1; dcmList[0] = readDICOMv(nameList.str[0], opts->isVerbose, opts->compressFlag, &dti4D); //ignore compile warning - memory only freed on first of 2 passes fillTDCMsort(dcmSort[0], 0, dcmList[0]); - return saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, &dti4D); + int ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, &dti4D); + return ret; }// singleDICOM() size_t fileBytes(const char * fname) { @@ -3019,10 +3077,10 @@ int nii_loadDir(struct TDCMopts* opts) { bool convertError = false; for (int i = 0; i < (int)nDcm; i++ ) { dcmList[i] = readDICOMv(nameList.str[i], opts->isVerbose, opts->compressFlag, &dti4D); //ignore compile warning - memory only freed on first of 2 passes - if ((dcmList[i].isValid) &&((dcmList[i].patientPositionNumPhilips > 1) || (dcmList[i].CSA.numDti > 1))) { //4D dataset: dti4D arrays require huge amounts of RAM - write this immediately - struct TDCMsort dcmSort[1]; - - fillTDCMsort(dcmSort[0], i, dcmList[i]); + //~ if ((dcmList[i].isValid) &&((dcmList[i].totalSlicesIn4DOrder != NULL) ||(dcmList[i].patientPositionNumPhilips > 1) || (dcmList[i].CSA.numDti > 1))) { //4D dataset: dti4D arrays require huge amounts of RAM - write this immediately + if ((dcmList[i].isValid) &&((dti4D.sliceOrder[0] >= 0) || (dcmList[i].CSA.numDti > 1))) { //4D dataset: dti4D arrays require huge amounts of RAM - write this immediately + struct TDCMsort dcmSort[1]; + fillTDCMsort(dcmSort[0], i, dcmList[i]); dcmList[i].converted2NII = 1; int ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, &dti4D); if (ret == EXIT_SUCCESS) From eb567b39c530b6774039a58a1c8622119756ee06 Mon Sep 17 00:00:00 2001 From: Chris Rorden Date: Mon, 19 Mar 2018 12:34:49 -0400 Subject: [PATCH 07/14] Restore ability to handle DICOMs files with Icon Image --- console/main_console.cpp | 2 +- console/nii_dicom.cpp | 3 +-- console/nii_dicom.h | 2 +- console/nii_dicom_batch.cpp | 11 +++++++++-- console/nii_dicom_batch.h | 3 ++- 5 files changed, 14 insertions(+), 7 deletions(-) diff --git a/console/main_console.cpp b/console/main_console.cpp index 13ec9209..074bdea0 100644 --- a/console/main_console.cpp +++ b/console/main_console.cpp @@ -360,7 +360,7 @@ int main(int argc, const char * argv[]) strcpy(opts.outdir,argv[i]); } else if ((argv[i][1] == 'n') && ((i+1) < argc)) { i++; - int seriesNumber = atoi(argv[i]); + float seriesNumber = atof(argv[i]); if (seriesNumber < 0) opts.numSeries = -1; //report series: convert none else if ((opts.numSeries >= 0) && (opts.numSeries < MAX_NUM_SERIES)) { diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 5abf9371..370229c0 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -3500,7 +3500,6 @@ double TE = 0.0; //most recent echo time recorded encapsulatedDataFragmentStart = (int)lPos + (int)lFileOffset; } } - if ((isIconImageSequence) && ((groupElement & 0x0028) == 0x0028 )) groupElement = kUnused; //ignore icon dimensions if ((sqEndPrivate > 0) && ((lFileOffset + lPos) > sqEndPrivate)) sqEndPrivate = -1; //end of private SQ with defined length if (groupElement == kSequenceDelimitationItemTag) { //end of private SQ with undefined length @@ -3509,8 +3508,8 @@ double TE = 0.0; //most recent echo time recorded sqDepthPrivate = 0; //no longer in a private SQ } } - if (sqDepth < 0) sqDepth = 0;*/ + if ((isIconImageSequence) && ((groupElement & 0x0028) == 0x0028 )) groupElement = kUnused; //ignore icon dimensions if (groupElement == kSequenceDelimitationItemTag) is2005140FSQ = false; switch ( groupElement ) { case kTransferSyntax: { diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 4ee7e6fe..040e0127 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -38,7 +38,7 @@ extern "C" { #define kCCsuf " CompilerNA" //unknown compiler! #endif -#define kDCMvers "v1.0.20171215" kDCMsuf kCCsuf +#define kDCMvers "v1.0.20180303" kDCMsuf kCCsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic static const int kMaxDTI4D = 4096; //maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index fd269c12..73c50553 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -2500,7 +2500,10 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc for(int i = 0; i < nConvert; ++i) dcmList[dcmSort[i].indx].converted2NII = 1; if (opts.numSeries < 0) { //report series number but do not convert - printMessage("\t%ld\t%s\n", dcmList[dcmSort[0].indx].seriesNum, pathoutname); + if (vol >= 0) + printMessage("\t%ld.%d\t%s\n", dcmList[dcmSort[0].indx].seriesNum, vol, pathoutname); + else + printMessage("\t%ld\t%s\n", dcmList[dcmSort[0].indx].seriesNum, pathoutname); printMessage(" %s\n",nameList->str[dcmSort[0].indx]); return EXIT_SUCCESS; } @@ -2521,8 +2524,12 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc // skip converting if user has specified one or more series, but has not specified this one if (opts.numSeries > 0) { int i = 0; + float seriesNum = (float) dcmList[dcmSort[0].indx].seriesNum; + if (vol > 0) + seriesNum = seriesNum + (float) vol / 10.0; //n.b. we will have problems if vol > 9. However, 9 distinct TEs/scalings/PhaseMag seems unlikely for (; i < opts.numSeries; i++) { - if (opts.seriesNumber[i] == dcmList[dcmSort[0].indx].seriesNum) { + if (isSameFloatGE(opts.seriesNumber[i], seriesNum)) { + //if (opts.seriesNumber[i] == dcmList[dcmSort[0].indx].seriesNum) { break; } } diff --git a/console/nii_dicom_batch.h b/console/nii_dicom_batch.h index ced863e7..9cf0a7dc 100644 --- a/console/nii_dicom_batch.h +++ b/console/nii_dicom_batch.h @@ -28,7 +28,8 @@ extern "C" { bool isSave3D,isGz, isFlipY, isCreateBIDS, isSortDTIbyBVal, isAnonymizeBIDS, isOnlyBIDS, isCreateText, isIgnoreDerivedAnd2D, isPhilipsFloatNotDisplayScaling, isTiltCorrect, isRGBplanar, isOnlySingleFile, isForceStackSameSeries, isCrop; int isVerbose, compressFlag, gzLevel; //support for compressed data 0=none, char filename[512], outdir[512], indir[512], pigzname[512], optsname[512], indirParent[512], imageComments[24]; - long seriesNumber[MAX_NUM_SERIES], numSeries; + float seriesNumber[MAX_NUM_SERIES]; + long numSeries; #ifdef HAVE_R bool isScanOnly; void *imageList; From a72d18bb0fa34a86029dd148c0eef9b0ee4d5ca1 Mon Sep 17 00:00:00 2001 From: Chris Rorden Date: Tue, 20 Mar 2018 08:57:35 -0400 Subject: [PATCH 08/14] Restore support for compressed DICOM (Jon Clayden) --- README.md | 19 ++++++++++--------- console/nii_dicom.cpp | 9 +++++++++ console/nii_dicom.h | 2 +- 3 files changed, 20 insertions(+), 10 deletions(-) diff --git a/README.md b/README.md index c3f115f8..25c7c2b9 100644 --- a/README.md +++ b/README.md @@ -72,18 +72,19 @@ If you have any problems with the cmake build script described above or want to ## Links - - [Dcm2Bids](https://github.com/cbedetti/Dcm2Bids) uses dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets. + - [bidsify](https://github.com/spinoza-rec/bidsify) is a Python project that uses dcm2niix to convert DICOM and Philips PAR/REC images to the BIDS standard. - [bidskit](https://github.com/jmtyszka/bidskit) uses dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets. + - [boutiques-dcm2niix](https://github.com/lalet/boutiques-dcm2niix) is a dockerfile for installing and validating dcm2niix. - [DAC2BIDS](https://github.com/dangom/dac2bids) uses dcm2niibatch to create [BIDS](http://bids.neuroimaging.io/) datasets. - - [heudiconv](https://github.com/nipy/heudiconv) can use dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets. - - [nipype](https://github.com/nipy/nipype) can use dcm2niix to convert images. - - [pydcm2niix is a Python module for working with dcm2niix](https://github.com/jstutters/pydcm2niix). + - [Dcm2Bids](https://github.com/cbedetti/Dcm2Bids) uses dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets. - [dcm2niir](https://github.com/muschellij2/dcm2niir) R wrapper for dcm2niix/dcm2nii. - - [divest](https://github.com/jonclayden/divest) R interface to dcm2niix. - - [sci-tran dcm2niix](https://github.com/scitran-apps/dcm2niix) Flywheel Gear (docker). - - [neuro_docker](https://github.com/Neurita/neuro_docker) includes dcm2niix as part of a single, static Dockerfile. - - [boutiques-dcm2niix](https://github.com/lalet/boutiques-dcm2niix) is a dockerfile for installing and validating dcm2niix. - - [neurodocker](https://github.com/kaczmarj/neurodocker) generates [custom](https://github.com/rordenlab/dcm2niix/issues/138) Dockerfiles given specific versions of neuroimaging software. - [dcm2niix_afni](https://afni.nimh.nih.gov/pub/dist/doc/program_help/dcm2niix_afni.html) is a version of dcm2niix included with the [AFNI](https://afni.nimh.nih.gov/) distribution. + - [divest](https://github.com/jonclayden/divest) R interface to dcm2niix. + - [heudiconv](https://github.com/nipy/heudiconv) can use dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets. - [MRIcroGL](https://github.com/neurolabusc/MRIcroGL) is available for MacOS, Linux and Windows and provides a graphical interface for dcm2niix. You can get compiled copies from the [MRIcroGL NITRC web site](https://www.nitrc.org/projects/mricrogl/). + - [neuro_docker](https://github.com/Neurita/neuro_docker) includes dcm2niix as part of a single, static Dockerfile. - [NeuroDebian](http://neuro.debian.net/pkgs/dcm2niix.html) provides up-to-date version of dcm2niix for Debian-based systems. + - [neurodocker](https://github.com/kaczmarj/neurodocker) generates [custom](https://github.com/rordenlab/dcm2niix/issues/138) Dockerfiles given specific versions of neuroimaging software. + - [nipype](https://github.com/nipy/nipype) can use dcm2niix to convert images. + - [pydcm2niix is a Python module for working with dcm2niix](https://github.com/jstutters/pydcm2niix). + - [sci-tran dcm2niix](https://github.com/scitran-apps/dcm2niix) Flywheel Gear (docker). \ No newline at end of file diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 370229c0..8a9bfe4e 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -3509,6 +3509,15 @@ double TE = 0.0; //most recent echo time recorded } } if (sqDepth < 0) sqDepth = 0;*/ + if ((groupElement == kItemTag) && (isEncapsulatedData)) { + d.imageBytes = dcmInt(4,&buffer[lPos-4],d.isLittleEndian); + //printMessage("compressed data %d-> %ld\n",d.imageBytes, lPos); + if (d.imageBytes > 128) { + encapsulatedDataFragments++; + if (encapsulatedDataFragmentStart == 0) + encapsulatedDataFragmentStart = (int)lPos + (int)lFileOffset; + } + } if ((isIconImageSequence) && ((groupElement & 0x0028) == 0x0028 )) groupElement = kUnused; //ignore icon dimensions if (groupElement == kSequenceDelimitationItemTag) is2005140FSQ = false; switch ( groupElement ) { diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 040e0127..50ee8b2c 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -38,7 +38,7 @@ extern "C" { #define kCCsuf " CompilerNA" //unknown compiler! #endif -#define kDCMvers "v1.0.20180303" kDCMsuf kCCsuf +#define kDCMvers "v1.0.20180304" kDCMsuf kCCsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic static const int kMaxDTI4D = 4096; //maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images From 84a140e16db918f26176c5d5320e42199cf1e338 Mon Sep 17 00:00:00 2001 From: Chris Rorden Date: Tue, 20 Mar 2018 17:32:13 -0400 Subject: [PATCH 09/14] Preliminary attempt to support PAR/REC with arbitrary slice orders (e.g. order reconstructed, not acquired) --- console/nii_dicom.cpp | 90 +++++++++++++++++++++++++++++-------------- console/nii_dicom.h | 2 +- 2 files changed, 63 insertions(+), 29 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 8a9bfe4e..6425dcfc 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -1424,8 +1424,14 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D #define kv2 45 #define kv3 46 #define kASL 48 - printError("dcm2niix PAR is not actively supported. Use with extreme caution (hint: use dicm2nii)\n"); + printWarning("dcm2niix PAR is not actively supported. Use with extreme caution (hint: use dicm2nii)\n"); char buff[LINESZ]; + int maxNumberOfDiffusionValues = 1; + int maxNumberOfGradientOrients = 1; + int maxNumberOfCardiacPhases = 1; + int maxNumberOfEchoes = 1; + int maxNumberOfDynamics = 1; + int maxNumberOfMixes = 1; int sliceNumberMrPhilipsB2[kMaxDTI4D], sliceNumberMrPhilipsVol2[kMaxDTI4D]; int patientPositionNumPhilipsB2 = 0; int patientPositionNumPhilipsVol2 = 0; @@ -1436,6 +1442,7 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D int minDyn = 32767; int maxDyn = 0; bool ADCwarning = false; + int numSlice2D = 0; int prevDyn = -1; bool dynNotAscending = false; int parVers = 0; @@ -1529,6 +1536,26 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D if ((strcmp(Comment[0], "Max.") == 0) && (strcmp(Comment[3], "slices/locations") == 0)) { d.xyzDim[3] = atoi(Comment[5]); } + if ((strcmp(Comment[0], "Max.") == 0) && (strcmp(Comment[3], "diffusion") == 0)) { + maxNumberOfDiffusionValues = atoi(Comment[6]); + } + if ((strcmp(Comment[0], "Max.") == 0) && (strcmp(Comment[3], "gradient") == 0)) { + maxNumberOfGradientOrients = atoi(Comment[6]); + } + if ((strcmp(Comment[0], "Max.") == 0) && (strcmp(Comment[3], "cardiac") == 0)) { + maxNumberOfCardiacPhases = atoi(Comment[6]); + } + if ((strcmp(Comment[0], "Max.") == 0) && (strcmp(Comment[3], "echoes") == 0)) { + maxNumberOfEchoes = atoi(Comment[5]); + } + if ((strcmp(Comment[0], "Max.") == 0) && (strcmp(Comment[3], "dynamics") == 0)) { + maxNumberOfDynamics = atoi(Comment[5]); + } + if ((strcmp(Comment[0], "Max.") == 0) && (strcmp(Comment[3], "mixes") == 0)) { + maxNumberOfMixes = atoi(Comment[5]); + if (maxNumberOfMixes > 1) + printError("maxNumberOfMixes > 1. Please update this software to support these images\n"); + } p = fgets (buff, LINESZ, fp);//get next line continue; } //process '.' tag @@ -1607,26 +1634,26 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D d.patientPosition[2] = cols[kPositionAP]; d.patientPosition[3] = cols[kPositionFH]; } - //~ - /* - if (d.patientPositionNumPhilips < kMaxDTI4D) { - dti4D->S[d.patientPositionNumPhilips].sliceNumberMrPhilips = round(cols[kSlice]); - if ((d.patientPositionNumPhilips > 0) && (dti4D->S[d.patientPositionNumPhilips].sliceNumberMrPhilips < dti4D->S[d.patientPositionNumPhilips-1].sliceNumberMrPhilips)) { - d.isSlicesSpatiallySequentialPhilips = false; - } - } - d.patientPositionNumPhilips++;*/ - /* - if (patientPositionNumPhilips < kMaxDTI4D) { - dti4D->S[patientPositionNumPhilips].sliceNumberMrPhilips = round(cols[kSlice]); - if ((patientPositionNumPhilips > 0) && (dti4D->S[patientPositionNumPhilips].sliceNumberMrPhilips < dti4D->S[patientPositionNumPhilips-1].sliceNumberMrPhilips)) { - isSlicesSpatiallySequentialPhilips = false; - } - }*/ + patientPositionNumPhilips++; + } + if (true) { int slice = round(cols[kSlice]); if (slice != (prevSlice + 1)) isSlicesSpatiallySequentialPhilips = false; prevSlice = slice; - patientPositionNumPhilips++; + //(cols[kEcho] == 1) && (cols[kDyn] == 2) && (patientPositionNumPhilipsVol2 < kMaxDTI4D) && (cols[kCardiac] == 1) && (cols[kGradientNumber] == 1) + int volStep = maxNumberOfDynamics; + int vol = ((int)cols[kDyn] - 1); + vol = vol + (volStep * ((int)cols[kGradientNumber] - 1)); + volStep = volStep * maxNumberOfGradientOrients; + vol = vol + (volStep * ((int)cols[kEcho] - 1)); + volStep = volStep * maxNumberOfEchoes; + vol = vol + (volStep * ((int)cols[kCardiac] - 1)); + volStep = volStep * maxNumberOfCardiacPhases; + slice = slice + (vol * d.xyzDim[3]); + if (numSlice2D < kMaxSlice2D) { + dti4D->sliceOrder[slice -1] = numSlice2D; + } + numSlice2D++; } if (cols[kGradientNumber] > 0) { /*int dir = (int) cols[kGradientNumber]; @@ -1639,12 +1666,13 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D //(cols[kImageType] == 0) means magnitude scan if ((cols[kImageType] == 0) && (cols[kDyn] == 1) && (cols[kEcho] == 1) && (cols[kCardiac] == 1) && (cols[kSlice] == 1)) { //only first slice d.CSA.numDti++; - int dir = d.CSA.numDti; + //int dir = d.CSA.numDti; + int dir =(int)cols[kGradientNumber]; if (dir <= kMaxDTI4D) { - if (isVerbose ) { + //if (isVerbose ) { if (d.CSA.numDti == 1) printMessage("n\tdir\tbValue\tV1\tV2\tV3\n"); printMessage("%d\t%g\t%g\t%g\t%g\t%g\n", dir-1, cols[kGradientNumber], cols[kbval], cols[kv1], cols[kv2], cols[kv3]); - } + //} dti4D->S[dir-1].V[0] = cols[kbval]; dti4D->S[dir-1].V[1] = cols[kv1]; dti4D->S[dir-1].V[2] = cols[kv2]; @@ -1667,8 +1695,12 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D for (int s = 0; s < patientPositionNumPhilipsVol2; s++) sliceNumberMrPhilipsVol2[s] = sliceNumberMrPhilipsB2[s]; } - if (!isSlicesSpatiallySequentialPhilips) - printError("PAR file order of slices is not sequential (hint: use dicm2nii)\n"); + if (!isSlicesSpatiallySequentialPhilips) { + printWarning("PAR file order of slices is not sequential: please ensure correct slice re-ordering\n"); + if (numSlice2D > kMaxSlice2D) + printError("Overloaded slice re-ordering. Number of slices (%d) exceeds kMaxSlice2D (%d)\n", patientPositionNumPhilips, kMaxSlice2D); + } else + dti4D->sliceOrder[0] = -1; //sequential slices do not need to be re-ordered //~ /*if ((patientPositionNumPhilipsVol2 > 1) && (d.patientPositionNumPhilips == patientPositionNumPhilipsVol2)) { bool isSliceOrderConsistent = true; @@ -1679,8 +1711,8 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D d.isValid = false; }*/ if (dynNotAscending) { - printError("PAR file volumes not saved in ascending order (hint: use dicm2nii)\n"); - d.isValid = false; + printWarning("PAR file volumes not saved in ascending order (please check re-ordering)\n"); + //d.isValid = false; } if ((slice % d.xyzDim[3]) != 0) { printError("Total number of slices (%d) not divisible by slices per 3D volume (%d) [acquisition aborted]. Try dicm2nii or nii_rescue_par to fix this: %s\n", slice, d.xyzDim[3], parname); @@ -2232,7 +2264,9 @@ unsigned char * nii_reorderSlicesX(unsigned char* bImg, struct nifti_1_header *h for (int i = 0; i < dim3to7; i++) { //for each volume int fromSlice = dti4D->sliceOrder[i]; //if (i < 10) printMessage(" ===> Moving slice from/to positions\t%d\t%d\n", i, toSlice); - if (i != fromSlice) { + if ((i < 0) || (fromSlice >= dim3to7)) + printError("Re-ordered slice out-of-volume %d\n", fromSlice); + else if (i != fromSlice) { uint64_t inPos = fromSlice * sliceBytes; uint64_t outPos = i * sliceBytes; memcpy( &bImg[outPos], &outImg[inPos], sliceBytes); @@ -3509,7 +3543,7 @@ double TE = 0.0; //most recent echo time recorded } } if (sqDepth < 0) sqDepth = 0;*/ - if ((groupElement == kItemTag) && (isEncapsulatedData)) { + if ((groupElement == kItemTag) && (isEncapsulatedData)) { //use this to find image fragment for compressed datasets, e.g. JPEG transfer syntax d.imageBytes = dcmInt(4,&buffer[lPos-4],d.isLittleEndian); //printMessage("compressed data %d-> %ld\n",d.imageBytes, lPos); if (d.imageBytes > 128) { @@ -4276,7 +4310,7 @@ double TE = 0.0; //most recent echo time recorded //~ d.isSlicesSpatiallySequentialPhilips = false; //slices not sequential (1,2,3,4 or 4,3,2,1) but 4,3,1,2 d.patientPositionNumPhilips++; //Philips can save 3D acquisitions in a single file with slices stored in non-sequential order. We need to know the first and final spatial position - //printMessage("xxxx====> - %d\n", sliceNumber); + //printMessage("x====> - %d\n", sliceNumber); if (isVerbose > 1) printMessage("slice %d is spatial position %d\n", d.patientPositionNumPhilips, sliceNumber); */ diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 50ee8b2c..234ba481 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -38,7 +38,7 @@ extern "C" { #define kCCsuf " CompilerNA" //unknown compiler! #endif -#define kDCMvers "v1.0.20180304" kDCMsuf kCCsuf +#define kDCMvers "v1.0.20180305" kDCMsuf kCCsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic static const int kMaxDTI4D = 4096; //maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images From 34f44d1c69e0814e5c5c05ddb1e93a4bfbcdbab1 Mon Sep 17 00:00:00 2001 From: Chris Rorden Date: Wed, 21 Mar 2018 09:38:47 -0400 Subject: [PATCH 10/14] 0020,9157 does not necessarily index gradient directions from 1..n (e.g. 2...37 for ADNI 018_S_4868) --- Philips/README.md | 18 +++++++++---- console/nii_dicom.cpp | 63 ++++++++++++++++++++++++------------------- console/nii_dicom.h | 2 +- 3 files changed, 49 insertions(+), 34 deletions(-) diff --git a/Philips/README.md b/Philips/README.md index dcbd5dfe..b961f8d0 100644 --- a/Philips/README.md +++ b/Philips/README.md @@ -10,29 +10,37 @@ The Image Patient Position (0020,0032) tag is required to determine the position In practice, this complication has several major implications. First, not all software translate private SQs well. One potential problem is when the implicit VRs are saved as implicit VRs. This can obscure the fact that 2005,140F is an SQ. Indeed, some tools will convert the private SQ type as a "UN" unknown type and add another public sequence. This can confuse reading software. -Furthermore, in the real world there are many Philips DICOM images that ONLY contain IPP inside SQ 2005,140F. It is not clear if this reflect the source Philips data or handling by other tools. +Furthermore, in the real world there are many Philips DICOM images that ONLY contain IPP inside SQ 2005,140F. These situations appear to reflect modifications applied by a PACS to the DICOM data or attempts to anonymize the DICOM images (e.g. using simple Matlab scripts). Note that the private IPP differs from the public one by half a voxel. Therefore, in theory if one only has the private IPP, one can infer the public IPP location. Current versions of dcm2niix do not do this: the error is assumed small enough that it will not impact image processing steps such as coregistration. + +In general it is recommended that you archive and convert DICOM images as they are created from the scanner. If one does use an export tool such as the excellent dcmtk, it is recommended that you preserve the explicit VR, as implicit VR has the potential of obscuring private sequence (SQ) tags. Be aware that subsequent processing of DICOM data can disrupt data conversion. Therefore, dcm2niix will ignore the IPP enclosed in 2005,140F unless no alternative exists. -## Derived Apparent Diffusion Coefficient maps stored with raw diffusion data +## Derived parametric maps stored with raw diffusion data -The DICOM standard requires derived data to be saved as a separate series number than the raw data. However, many Philips diffusion images append the derived ADC maps with the original data. Strangely, some of these ADC maps seem to be stored with specific diffusion directions, and with the same volume number as the B=0 image. In these cases, 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 ADC 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. For example, ADC, Trace and Isotropic images can all be derived from the raw scans. As scientists, we want to discard these derived images, as they will disrupt data processing and we can generate better parametric maps after we have applied undistortion methods such as [Eddy and Topup](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/eddy/UsersGuide). The current version of dcm2niix uses the Diffusion Directionality (0018,9075) tag to detect B=0 unweighted ("NONE"), B-weighted ("DIRECTIONAL"), and derived ("ISOTROPIC") images. Note that the Dimension Index Values (0020,9157) tag provides an alternative approach to discriminate these images. Here are sample tags from a Philips enhanced image that includes and derived map (3rd dimension is "1" while the other images set this to "2"). ``` (0018,9075) CS [DIRECTIONAL] +(0018,9089) FD 1\0\0 +(0018,9087) FD 1000 (0020,9157) UL 1\1\2\32 ... (0018,9075) CS [ISOTROPIC] +(0018,9087) FD 1000 (0020,9157) UL 1\2\1\33 ... (0018,9075) CS [NONE] +(0018,9087) FD 0 (0020,9157) UL 1\1\2\33 ``` ## Diffusion Direction -Some Philips enhanced scans use tag 0018,9089 to report the 3 gradient directions. Other files use the tags 2005,10b0, 2005,10b1, 2005,10b2. In general, dcm2niix will use the values that most closely precede the Dimension Index Values (0020,9157). +Proper Philips enhanced scans use tag 0018,9089 to report the 3 gradient directions. However, in the wild, other files from Philips (presumably using older versions of Philips software) use the tags 2005,10b0, 2005,10b1, 2005,10b2. In general, dcm2niix will use the values that most closely precede the Dimension Index Values (0020,9157). + +The current version of dcm2niix uses Dimension Index Values (0020,9157) to determine gradient number, which can also be found in (2005,1413). However, while 2005,1413 is always indexed from one, this is not necessarily the case for 0020,9157. For example, the ADNI DWI dataset for participant 018_S_4868 has values of 2005,1413 that range from 1..36 for the 36 directions, while 0020,9157 ranges from 2..37. The current version of dcm2niix compensates for this by re-indexing the values of 0020,9157 after all the volumes have been read. ## General variations -The tags SliceNumberMrPhilips (2001,100A), NumberOfSlicesMrPhilips (2001,1018), and InStackPositionNumber (0020,9057) MRImageGradientOrientationNumber (2005,1413) are all potentially very useful for simplifying handling of enhanced DICOM images. However, none of these are reliably stored in all DICOM images. Therefore, dcm2niix does not depend on any of these, though it attempts to use several of these if they are available. \ No newline at end of file +Prior versions of dcm2niix used different methods to sort images. However, these have proved unreliable The undocumented tags SliceNumberMrPhilips (2001,100A), NumberOfSlicesMrPhilips (2001,1018) are not always present. In theory, InStackPositionNumber (0020,9057) should be present in all enhanced files, but has not proved reliable (perhaps not in older Philips images or DICOM images that were modified after leaving the scanner). MRImageGradientOrientationNumber (2005,1413) is complicated by the inclusion of derived images.Therefore, current versions of dcm2niix do not depend on any of these. diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 6425dcfc..ffb39347 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -1669,10 +1669,10 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D //int dir = d.CSA.numDti; int dir =(int)cols[kGradientNumber]; if (dir <= kMaxDTI4D) { - //if (isVerbose ) { + if (isVerbose ) { if (d.CSA.numDti == 1) printMessage("n\tdir\tbValue\tV1\tV2\tV3\n"); printMessage("%d\t%g\t%g\t%g\t%g\t%g\n", dir-1, cols[kGradientNumber], cols[kbval], cols[kv1], cols[kv2], cols[kv3]); - //} + } dti4D->S[dir-1].V[0] = cols[kbval]; dti4D->S[dir-1].V[1] = cols[kv1]; dti4D->S[dir-1].V[2] = cols[kv2]; @@ -3351,12 +3351,14 @@ double TE = 0.0; //most recent echo time recorded int numberOfFrames = 0; int MRImageGradientOrientationNumber = 0; //int maxMRImageGradientOrientationNumber = -1; + int minGradNum = kMaxDTI4D + 1; int maxGradNum = -1; int numberOfDynamicScans = 0; uint32_t lLength; uint32_t groupElement; long lPos = 0; - bool isPhilipsADC = false; + bool isPhilipsDerived = false; + bool isPhilipsDiffusion = false; if (isPart10prefix) { //for part 10 files, skip preamble and prefix lPos = 128+4; //4-byte signature starts at 128 groupElement = buffer[lPos] | (buffer[lPos+1] << 8) | (buffer[lPos+2] << 16) | (buffer[lPos+3] << 24); @@ -3751,7 +3753,7 @@ double TE = 0.0; //most recent echo time recorded char dir[kDICOMStr]; dcmStr (lLength, &buffer[lPos], dir); if (strcmp(dir, "ISOTROPIC") == 0) - isPhilipsADC = true; + isPhilipsDerived = true; break; } case kDwellTime : d.dwellTime = dcmStrInt(lLength, &buffer[lPos]); @@ -3908,7 +3910,7 @@ double TE = 0.0; //most recent echo time recorded next two lines attempt to skip ADC maps we could also increment gradNum for ADC if we wanted... */ - if (isPhilipsADC) { + if (isPhilipsDerived) { gradNum ++; B0Philips = 2000.0; vRLPhilips = 0.0; @@ -3916,8 +3918,14 @@ double TE = 0.0; //most recent echo time recorded vFHPhilips = 0.0; } + if (B0Philips == 0.0) { + //printMessage(" DimensionIndexValues grad %d b=%g vec=%gx%gx%g\n", gradNum, B0Philips, vRLPhilips, vAPPhilips, vFHPhilips); + vRLPhilips = 0.0; + vAPPhilips = 0.0; + vFHPhilips = 0.0; + } //if ((MRImageGradientOrientationNumber > 0) && ((gradNum != MRImageGradientOrientationNumber)) break; - + if (gradNum < minGradNum) minGradNum = gradNum; if (gradNum >= maxGradNum) maxGradNum = gradNum; if (gradNum >= kMaxDTI4D) { printError("Number of DTI gradients exceeds 'kMaxDTI4D (%d).\n", kMaxDTI4D); @@ -3928,7 +3936,7 @@ double TE = 0.0; //most recent echo time recorded philDTI[gradNum].V[1] = vRLPhilips; philDTI[gradNum].V[2] = vAPPhilips; philDTI[gradNum].V[3] = vFHPhilips; - isPhilipsADC = false; + isPhilipsDerived = false; //printMessage(" DimensionIndexValues grad %d b=%g vec=%gx%gx%g\n", gradNum, B0Philips, vRLPhilips, vAPPhilips, vFHPhilips); @@ -4609,6 +4617,7 @@ if (d.isHasPhase) d.xyzDim[3] = d.numberOfDynamicScans; }*/ + if (numberOfFrames == 0) numberOfFrames = d.xyzDim[3]; if ((numberOfDynamicScans > 1) && (d.xyzDim[4] < 2) && (d.xyzDim[3] > 1) && ((d.xyzDim[3] % numberOfDynamicScans) == 0)) { @@ -4713,30 +4722,28 @@ if (d.isHasPhase) philDTI[maxGradNum].V[3] = 0.0; maxGradNum++; }*/ + if ((minGradNum >= 1) && ((maxGradNum-minGradNum+1) == d.xyzDim[4])) { + //see ADNI DWI data for 018_S_4868 - the gradient numbers are in the range 2..37 for 36 volumes - no gradient number 1! + if (philDTI[minGradNum -1].V[0] >= 0) { + if (isVerbose) + printMessage("Using %d diffusion data directions coded by DimensionIndexValues\n", maxGradNum); + int off = 0; + if (minGradNum > 1) { + off = minGradNum - 1; + printWarning("DimensionIndexValues (0020,9157) is not indexed from 1 (range %d..%d). Please validate results\n", minGradNum, maxGradNum); + } + for (int i = 0; i < d.xyzDim[4]; i++) { + dti4D->S[i].V[0] = philDTI[i+off].V[0]; + dti4D->S[i].V[1] = philDTI[i+off].V[1]; + dti4D->S[i].V[2] = philDTI[i+off].V[2]; + dti4D->S[i].V[3] = philDTI[i+off].V[3]; + if (isVerbose > 1) + printMessage(" grad %d b=%g vec=%gx%gx%g\n", i, dti4D->S[i].V[0], dti4D->S[i].V[1], dti4D->S[i].V[2], dti4D->S[i].V[3]); - - //printMessage("CXC grad %g %d %d\n", philDTI[0].V[0], maxGradNum, d.xyzDim[4]); - - if ((philDTI[0].V[0] >= 0) && (maxGradNum == d.xyzDim[4])) { - if (isVerbose) - printMessage("Using diffusion data coded by DimensionIndexValues\n"); - for (int i = 0; i < d.xyzDim[4]; i++) { - dti4D->S[i].V[0] = philDTI[i].V[0]; - dti4D->S[i].V[1] = philDTI[i].V[1]; - dti4D->S[i].V[2] = philDTI[i].V[2]; - dti4D->S[i].V[3] = philDTI[i].V[3]; - if (isVerbose > 1) - printMessage(" grad %d b=%g vec=%gx%gx%g\n", i, dti4D->S[i].V[0], dti4D->S[i].V[1], dti4D->S[i].V[2], dti4D->S[i].V[3]); - + } + d.CSA.numDti = maxGradNum - off; } - d.CSA.numDti = maxGradNum; } - //~~ - /*if (d.patientPositionNumPhilips >= kMaxDTI4D) { - - printError("Too many 2D slices in a single file [recompile with increased kMaxDTI4D] detected=%d, max = %d\n", d.patientPositionNumPhilips, kMaxDTI4D); - d.CSA.numDti = 0; - }*/ if (d.CSA.numDti >= kMaxDTI4D) { printError("Unable to convert DTI [recompile with increased kMaxDTI4D] detected=%d, max = %d\n", d.CSA.numDti, kMaxDTI4D); d.CSA.numDti = 0; diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 234ba481..e4e7a3cc 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -38,7 +38,7 @@ extern "C" { #define kCCsuf " CompilerNA" //unknown compiler! #endif -#define kDCMvers "v1.0.20180305" kDCMsuf kCCsuf +#define kDCMvers "v1.0.20180306" kDCMsuf kCCsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic static const int kMaxDTI4D = 4096; //maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images From c18c6ef8fa90619dca3adf9592eccd967b25ba96 Mon Sep 17 00:00:00 2001 From: Chris Rorden Date: Wed, 21 Mar 2018 18:02:47 -0400 Subject: [PATCH 11/14] Better PAR/REC error reporting --- console/nii_dicom.cpp | 18 +++++++++++++++--- console/nii_dicom.h | 2 +- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index ffb39347..9776b933 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -1649,8 +1649,10 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D volStep = volStep * maxNumberOfEchoes; vol = vol + (volStep * ((int)cols[kCardiac] - 1)); volStep = volStep * maxNumberOfCardiacPhases; + //if (slice == 1) printWarning("%d\n", (int)cols[kEcho]); slice = slice + (vol * d.xyzDim[3]); - if (numSlice2D < kMaxSlice2D) { + + if ((slice >= 0) && (slice < kMaxSlice2D) && (numSlice2D < kMaxSlice2D) && (numSlice2D >= 0)) { dti4D->sliceOrder[slice -1] = numSlice2D; } numSlice2D++; @@ -1697,10 +1699,19 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D } if (!isSlicesSpatiallySequentialPhilips) { printWarning("PAR file order of slices is not sequential: please ensure correct slice re-ordering\n"); - if (numSlice2D > kMaxSlice2D) - printError("Overloaded slice re-ordering. Number of slices (%d) exceeds kMaxSlice2D (%d)\n", patientPositionNumPhilips, kMaxSlice2D); + if (numSlice2D > kMaxSlice2D) { + printError("Overloaded slice re-ordering. Number of slices (%d) exceeds kMaxSlice2D (%d)\n", numSlice2D, kMaxSlice2D); + dti4D->sliceOrder[0] = -1; + } } else dti4D->sliceOrder[0] = -1; //sequential slices do not need to be re-ordered + int numExpected = (d.xyzDim[3] * maxNumberOfDiffusionValues * maxNumberOfGradientOrients + * maxNumberOfCardiacPhases * maxNumberOfEchoes * maxNumberOfDynamics * maxNumberOfMixes); + if (numSlice2D != numExpected) { + printMessage(" found %d slices, but expected %d: slices*diff*grad*cardiac*echo*dynamic*mix = %d*%d*%d*%d*%d*%d*%d\n", numSlice2D, numExpected, + d.xyzDim[3], maxNumberOfDiffusionValues, maxNumberOfGradientOrients, + maxNumberOfCardiacPhases, maxNumberOfEchoes, maxNumberOfDynamics, maxNumberOfMixes); + } //~ /*if ((patientPositionNumPhilipsVol2 > 1) && (d.patientPositionNumPhilips == patientPositionNumPhilipsVol2)) { bool isSliceOrderConsistent = true; @@ -2257,6 +2268,7 @@ unsigned char * nii_reorderSlicesX(unsigned char* bImg, struct nifti_1_header *h for (int i = 3; i < 8; i++) if (hdr->dim[i] > 1) dim3to7 = dim3to7 * hdr->dim[i]; if (dim3to7 < 2) return bImg; + if (dim3to7 > kMaxSlice2D) return bImg; uint64_t imgSz = nii_ImgBytes(*hdr); uint64_t sliceBytes = hdr->dim[1]*hdr->dim[2]*hdr->bitpix/8; unsigned char *outImg = (unsigned char *)malloc( imgSz); diff --git a/console/nii_dicom.h b/console/nii_dicom.h index e4e7a3cc..8d974c3f 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -42,7 +42,7 @@ extern "C" { static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic static const int kMaxDTI4D = 4096; //maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images -static const int kMaxSlice2D = 16786; //maximum number of 2D slices in 4D (Philips) images +static const int kMaxSlice2D = 64000; //maximum number of 2D slices in 4D (Philips) images #define kDICOMStr 64 #define kDICOMStrLarge 256 From 806b1c1979c914d6ddcac4757a4ee840e7158e76 Mon Sep 17 00:00:00 2001 From: Chris Rorden Date: Sat, 24 Mar 2018 09:35:27 -0400 Subject: [PATCH 12/14] Improved Philips PARREC support --- console/nii_dicom.cpp | 118 ++++++++++++++++++++++++++++++------ console/nii_dicom.h | 6 +- console/nii_dicom_batch.cpp | 93 +++++++++++++++++++++++----- 3 files changed, 180 insertions(+), 37 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 9776b933..cbee2c49 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -696,6 +696,7 @@ struct TDICOMdata clear_dicom_data() { //~ d.patientPositionRepeats = 0; d.isHasPhase = false; d.isHasMagnitude = false; + //d.maxGradDynVol = -1; //PAR/REC only d.sliceOrient = kSliceOrientUnknown; d.dateTime = (double)19770703150928.0; d.acquisitionTime = 0.0f; @@ -1193,6 +1194,13 @@ int readCSAImageHeader(unsigned char *buff, int lLength, struct TCSAdata *CSA, i return EXIT_SUCCESS; } // readCSAImageHeader() +//int readCSAImageHeader(unsigned char *buff, int lLength, struct TCSAdata *CSA, int isVerbose) { +int readProtocolDataBlockGE(unsigned char *buff, int lLength, struct TCSAdata *CSA, int isVerbose) { + + return EXIT_SUCCESS; +} + + void dcmMultiShorts (int lByteLength, unsigned char lBuffer[], int lnShorts, uint16_t *lShorts, bool littleEndian) { //read array of unsigned shorts US http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_6.2.html if ((lnShorts < 1) || (lByteLength != (lnShorts * 2))) return; @@ -1382,7 +1390,7 @@ void changeExt (char *file_name, const char* ext) { } } //changeExt() -struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D *dti4D) { +struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D *dti4D, bool isReadPhase) { struct TDICOMdata d = clear_dicom_data(); strcpy(d.protocolName, ""); //erase dummy with empty strcpy(d.seriesDescription, ""); //erase dummy with empty @@ -1425,6 +1433,7 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D #define kv3 46 #define kASL 48 printWarning("dcm2niix PAR is not actively supported. Use with extreme caution (hint: use dicm2nii)\n"); + if (isReadPhase) printWarning(" Reading phase images from PAR/REC\n"); char buff[LINESZ]; int maxNumberOfDiffusionValues = 1; int maxNumberOfGradientOrients = 1; @@ -1450,6 +1459,9 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D int maxCardiac = 1; int nCols = 26; int slice = 0; + int numExpected = 0; + int maxVol = -1; + int patientPositionNumPhilips = 0; bool isSlicesSpatiallySequentialPhilips = true; //int prevSliceIndex = 0; //index of prior slice: detect if images are not in order @@ -1465,15 +1477,24 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D //dti4D->S[0].sliceNumberMrPhilips = -1.0; //dti4D->S[0].sliceNumberMrPhilipsVol2 = -1.0; //d.dti4D = (TDTI *)malloc(kMaxDTI4D * sizeof(TDTI)); + while (p) { if (strlen(buff) < 1) continue; if (buff[0] == '#') { //comment char Comment[7][50]; - sscanf(buff, "# %s %s\n", Comment[0], Comment[1]); + + //sscanf(buff, "# %s %s\n", Comment[0], Comment[1]); + + sscanf(buff, "# %s %s %s %s %s %s V%s\n", Comment[0], Comment[1], Comment[2], Comment[3],Comment[4], Comment[5],Comment[6]); + + if ((strcmp(Comment[0], "sl") == 0) && (strcmp(Comment[1], "ec") == 0) ) { + //printError("BINGO\n"); + numExpected = (d.xyzDim[3] * maxNumberOfDiffusionValues * maxNumberOfGradientOrients + * maxNumberOfCardiacPhases * maxNumberOfEchoes * maxNumberOfDynamics * maxNumberOfMixes); + } if (strcmp(Comment[1], "TRYOUT") == 0) { - sscanf(buff, "# %s %s %s %s %s %s V%s\n", Comment[0], Comment[1], Comment[2], Comment[3] - ,Comment[4], Comment[5],Comment[6]); + //sscanf(buff, "# %s %s %s %s %s %s V%s\n", Comment[0], Comment[1], Comment[2], Comment[3],Comment[4], Comment[5],Comment[6]); parVers = (int)round(atof(Comment[6])*10); //4.2 = 42 etc if (parVers < 40) { printMessage("This software is unable to convert ancient PAR files: please use legacy dcm2nii\n"); @@ -1547,6 +1568,7 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D } if ((strcmp(Comment[0], "Max.") == 0) && (strcmp(Comment[3], "echoes") == 0)) { maxNumberOfEchoes = atoi(Comment[5]); + if (maxNumberOfEchoes > 1) d.isMultiEcho = true; } if ((strcmp(Comment[0], "Max.") == 0) && (strcmp(Comment[3], "dynamics") == 0)) { maxNumberOfDynamics = atoi(Comment[5]); @@ -1570,6 +1592,21 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D } for (int i = 0; i <= nCols; i++) cols[i] = strtof(p, &p); // p+1 skip comma, read a float + //the initial header does not report if phase or magnitude are include, so we read each separately in two passes + /* + if ((!isReadPhase) && (cols[kImageType] != 0)) { //skip phase images + p = fgets (buff, LINESZ, fp);//get next line + numSlice2D++; + d.isHasPhase = true; + continue; + } + if ((isReadPhase) && (cols[kImageType] == 0)) { //skip magnitude images + p = fgets (buff, LINESZ, fp);//get next line + numSlice2D++; + d.isHasMagnitude = true; + continue; + } + */ if ((cols[kIndex]) != slice) isIndexSequential = false; //slices 0,1,2.. should have indices 0,1,2,3... slice ++; if (slice == 1) { @@ -1636,22 +1673,38 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D } patientPositionNumPhilips++; } - if (true) { + if (true) { //for every slice int slice = round(cols[kSlice]); - if (slice != (prevSlice + 1)) isSlicesSpatiallySequentialPhilips = false; - prevSlice = slice; //(cols[kEcho] == 1) && (cols[kDyn] == 2) && (patientPositionNumPhilipsVol2 < kMaxDTI4D) && (cols[kCardiac] == 1) && (cols[kGradientNumber] == 1) int volStep = maxNumberOfDynamics; int vol = ((int)cols[kDyn] - 1); + int gradDynVol = (int)cols[kGradientNumber] - 1; + if (vol > gradDynVol) gradDynVol = vol; //if fMRI series, else DWI vol = vol + (volStep * ((int)cols[kGradientNumber] - 1)); volStep = volStep * maxNumberOfGradientOrients; vol = vol + (volStep * ((int)cols[kEcho] - 1)); volStep = volStep * maxNumberOfEchoes; vol = vol + (volStep * ((int)cols[kCardiac] - 1)); volStep = volStep * maxNumberOfCardiacPhases; + if (cols[kImageType] != 0) //phase maps are not reported in the initial header! + vol = vol + volStep; + if (vol > maxVol) maxVol = vol; + if (vol <= kMaxDTI4D) { + // dti4D->S[vol].V[0] = cols[kbval]; + //dti4D->gradDynVol[vol] = gradDynVol; + dti4D->TE[vol] = cols[kTEcho]; + dti4D->intenIntercept[vol] = cols[kRI]; + dti4D->intenScale[vol] = cols[kRS]; + dti4D->intenScalePhilips[vol] = cols[kSS]; + dti4D->isPhase[vol] = (cols[kImageType] != 0); + } //if (slice == 1) printWarning("%d\n", (int)cols[kEcho]); slice = slice + (vol * d.xyzDim[3]); - + if (slice != (prevSlice + 1)) isSlicesSpatiallySequentialPhilips = false; + prevSlice = slice; + //if (cols[kImageType] != 0) //yikes - phase maps! + // slice = slice + numExpected; + //printWarning("%d\t%d\n", slice -1, numSlice2D); if ((slice >= 0) && (slice < kMaxSlice2D) && (numSlice2D < kMaxSlice2D) && (numSlice2D >= 0)) { dti4D->sliceOrder[slice -1] = numSlice2D; } @@ -1697,21 +1750,25 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D for (int s = 0; s < patientPositionNumPhilipsVol2; s++) sliceNumberMrPhilipsVol2[s] = sliceNumberMrPhilipsB2[s]; } - if (!isSlicesSpatiallySequentialPhilips) { + //even is sequential, we need to decode if multiple echoes/phasemaps,etc + //if (!isSlicesSpatiallySequentialPhilips) { printWarning("PAR file order of slices is not sequential: please ensure correct slice re-ordering\n"); if (numSlice2D > kMaxSlice2D) { printError("Overloaded slice re-ordering. Number of slices (%d) exceeds kMaxSlice2D (%d)\n", numSlice2D, kMaxSlice2D); dti4D->sliceOrder[0] = -1; } - } else - dti4D->sliceOrder[0] = -1; //sequential slices do not need to be re-ordered - int numExpected = (d.xyzDim[3] * maxNumberOfDiffusionValues * maxNumberOfGradientOrients - * maxNumberOfCardiacPhases * maxNumberOfEchoes * maxNumberOfDynamics * maxNumberOfMixes); - if (numSlice2D != numExpected) { + //} else + // dti4D->sliceOrder[0] = -1; //sequential slices do not need to be re-ordered + if ((d.isHasPhase) && (d.isHasMagnitude)) { + numExpected *= 2; //the initial header does not report that both types are stored + } + if (numSlice2D != numExpected) { printMessage(" found %d slices, but expected %d: slices*diff*grad*cardiac*echo*dynamic*mix = %d*%d*%d*%d*%d*%d*%d\n", numSlice2D, numExpected, d.xyzDim[3], maxNumberOfDiffusionValues, maxNumberOfGradientOrients, maxNumberOfCardiacPhases, maxNumberOfEchoes, maxNumberOfDynamics, maxNumberOfMixes); + d.isValid = false; } + //~ /*if ((patientPositionNumPhilipsVol2 > 1) && (d.patientPositionNumPhilips == patientPositionNumPhilipsVol2)) { bool isSliceOrderConsistent = true; @@ -1726,8 +1783,9 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D //d.isValid = false; } if ((slice % d.xyzDim[3]) != 0) { - printError("Total number of slices (%d) not divisible by slices per 3D volume (%d) [acquisition aborted]. Try dicm2nii or nii_rescue_par to fix this: %s\n", slice, d.xyzDim[3], parname); + printError("Total number of slices (%d) not divisible by slices per 3D volume (%d) [acquisition aborted]. Try dicm2nii: %s\n", slice, d.xyzDim[3], parname); d.isValid = false; + return d; } d.xyzDim[4] = slice/d.xyzDim[3]; d.locationsInAcquisition = d.xyzDim[3]; @@ -1823,6 +1881,20 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D if (fabs(TRms - d.TR) > 0.005f) printWarning("Reported TR=%gms, measured TR=%gms (prospect. motion corr.?)\n", d.TR, TRms); d.TR = TRms; + } + //check if dimensions vary + if (maxVol > 0) { //maxVol indexed from 0 + for (int i = 1; i <= maxVol; i++) { + //if (dti4D->gradDynVol[i] > d.maxGradDynVol) d.maxGradDynVol = dti4D->gradDynVol[i]; + if (dti4D->intenIntercept[i] != dti4D->intenIntercept[0]) d.isScaleOrTEVaries = true; + if (dti4D->intenScale[i] != dti4D->intenScale[0]) d.isScaleOrTEVaries = true; + if (dti4D->intenScalePhilips[i] != dti4D->intenScalePhilips[0]) d.isScaleOrTEVaries = true; + if (dti4D->isPhase[i] != dti4D->isPhase[0]) d.isScaleOrTEVaries = true; + } + if (d.isScaleOrTEVaries) + printWarning("Varying dimensions (echoes, phase maps, intensity scaling) will require volumes to be saved separately (hint: you may prefer dicm2nii output)\n"); + + } //check DTI makes sense if (d.CSA.numDti > 1) { @@ -2864,7 +2936,7 @@ unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct } //~ if ((hdr->dim[0] > 3) && (dcm.patientPositionSequentialRepeats > 1) && (dcm.sliceOrder == NULL)) //swizzle 3rd and 4th dimension (Philips stores time as 3rd dimension) //~ img = nii_XYTZ_XYZT(img, hdr,dcm.patientPositionSequentialRepeats); - if (dti4D->sliceOrder[0] >= 0) + if ((dti4D != NULL) && (dti4D->sliceOrder[0] >= 0)) img = nii_reorderSlicesX(img, hdr, dti4D); //~ /*if (((dcm.patientPositionSequentialRepeats * 2) == dcm.patientPositionRepeats) && (dcm.isHasPhase) && (dcm.isHasMagnitude)) { @@ -3242,7 +3314,8 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kZThick 0x0018+(0x0050 << 16 ) #define kTR 0x0018+(0x0080 << 16 ) #define kTE 0x0018+(0x0081 << 16 ) -#define kEffectiveTE 0x0018+(0x9082 << 16 ) +//#define kEffectiveTE 0x0018+(0x9082 << 16 ) +const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16); #define kTI 0x0018+(0x0082 << 16) // Inversion time #define kEchoNum 0x0018+(0x0086 << 16 ) //IS #define kMagneticFieldStrength 0x0018+(0x0087 << 16 ) //DS @@ -3294,6 +3367,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kDimensionIndexValues 0x0020+uint32_t(0x9157<< 16 ) // UL n-dimensional index of frame. #define kInStackPositionNumber 0x0020+uint32_t(0x9057<< 16 ) // UL can help determine slices in volume #define kLocationsInAcquisitionGE 0x0021+(0x104F<< 16 )// 'SS' 'LocationsInAcquisitionGE' +#define kProtocolDataBlockGE 0x0025+(0x101B<< 16 )// 'OB' #define kSamplesPerPixel 0x0028+(0x0002 << 16 ) #define kPhotometricInterpretation 0x0028+(0x0004 << 16 ) #define kPlanarRGB 0x0028+(0x0006 << 16 ) @@ -3400,7 +3474,7 @@ double TE = 0.0; //most recent echo time recorded bool isMagnitude = false; float patientPositionPrivate[4] = {NAN, NAN, NAN, NAN}; float patientPosition[4] = {NAN, NAN, NAN, NAN}; //used to compute slice direction for Philips 4D - float patientPositionPublic[4] = {NAN, NAN, NAN, NAN}; //used to compute slice direction for Philips 4D + //float patientPositionPublic[4] = {NAN, NAN, NAN, NAN}; //used to compute slice direction for Philips 4D float patientPositionEndPhilips[4] = {NAN, NAN, NAN, NAN}; float patientPositionStartPhilips[4] = {NAN, NAN, NAN, NAN}; /* //not required? @@ -3996,6 +4070,11 @@ double TE = 0.0; //most recent echo time recorded case kLocationsInAcquisitionGE: locationsInAcquisitionGE = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); break; + case kProtocolDataBlockGE : + if (d.manufacturer != kMANUFACTURER_GE) break; + //printError("ProtocolDataBlockGE to do %d %d %2x%2x\n", lPos, lLength, (uint8_t)buffer[lPos], (uint8_t)buffer[lPos+1]); + //readProtocolDataBlockGE(&buffer[lPos], lLength, &d.CSA, isVerbose); //, dti4D); + break; case kDoseCalibrationFactor : d.doseCalibrationFactor = dcmStrFloat(lLength, &buffer[lPos]); break; @@ -4698,18 +4777,19 @@ if (d.isHasPhase) int j = 0; if (d.xyzDim[3] > 1) j = 1; for (int i = 0; i < d.xyzDim[4]; i++) { + //dti4D->gradDynVol[i] = 0; //only PAR/REC dti4D->TE[i] = dcmDim[j+(i * d.xyzDim[3])].TE; dti4D->intenScale[i] = dcmDim[j+(i * d.xyzDim[3])].intenScale; dti4D->intenIntercept[i] = dcmDim[j+(i * d.xyzDim[3])].intenIntercept; dti4D->isPhase[i] = dcmDim[j+(i * d.xyzDim[3])].isPhase; dti4D->intenScalePhilips[i] = dcmDim[j+(i * d.xyzDim[3])].intenScalePhilips; if (dti4D->TE[i] != d.TE) isTEvaries = true; - if (dti4D->TE[i] != d.TE) ; if (dti4D->intenScale[i] != d.intenScale) isScaleVaries = true; if (dti4D->intenIntercept[i] != d.intenIntercept) isScaleVaries = true; if (dti4D->isPhase[i] != isPhase) d.isScaleOrTEVaries = true; } if((isScaleVaries) || (isTEvaries)) d.isScaleOrTEVaries = true; + if (isTEvaries) d.isMultiEcho = true; //if echoVaries,count number of echoes /*int echoNum = 1; for (int i = 1; i < d.xyzDim[4]; i++) { diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 8d974c3f..064c9d01 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -38,7 +38,7 @@ extern "C" { #define kCCsuf " CompilerNA" //unknown compiler! #endif -#define kDCMvers "v1.0.20180306" kDCMsuf kCCsuf +#define kDCMvers "v1.0.20180307" kDCMsuf kCCsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic static const int kMaxDTI4D = 4096; //maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images @@ -83,6 +83,7 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; struct TDTI4D { struct TDTI S[kMaxDTI4D]; int sliceOrder[kMaxSlice2D]; // [7,3,2] means the first slice on disk should be moved to 7th position + int gradDynVol[kMaxDTI4D]; //used to parse dimensions of Philips data, e.g. file with multiple dynamics, echoes, phase+magnitude float TE[kMaxDTI4D], intenScale[kMaxDTI4D], intenIntercept[kMaxDTI4D], intenScalePhilips[kMaxDTI4D]; bool isPhase[kMaxDTI4D]; }; @@ -125,6 +126,7 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; int xyzDim[5]; //numberOfDynamicScans, patientPositionNumPhilips //patientPositionSequentialRepeats,patientPositionRepeats, + //maxGradDynVol, gradDynVol, int modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, echoTrainLength, coilNum, echoNum, sliceOrient, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,locationsInAcquisition, compressionScheme; float patientWeight, zSpacing, zThick, pixelBandwidth, SAR, phaseFieldofView, accelFactPE, flipAngle, fieldStrength, TE, TI, TR, intenScale, intenIntercept, intenScalePhilips, gantryTilt, lastScanLoc, angulation[4]; float orient[7], patientPosition[4], patientPositionLast[4], xyzMM[4], stackOffcentre[4]; @@ -147,7 +149,7 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, struct TDTI4D *dti4D); struct TDICOMdata readDICOM(char * fname); struct TDICOMdata clear_dicom_data(); - struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D *dti4D); + struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D *dti4D, bool isReadPhase); unsigned char * nii_flipY(unsigned char* bImg, struct nifti_1_header *h); unsigned char * nii_flipZ(unsigned char* bImg, struct nifti_1_header *h); //*unsigned char * nii_reorderSlices(unsigned char* bImg, struct nifti_1_header *h, struct TDTI4D *dti4D); diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 73c50553..153f8ba2 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -1515,6 +1515,10 @@ int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopt sprintf(newstr, "_e%d", dcm.echoNum); strcat (outname,newstr); } + /*if (dcm.maxGradDynVol > 0) { //Philips segmented + sprintf(newstr, "_v%04d", dcm.gradDynVol+1); //+1 as indexed from zero + strcat (outname,newstr); + }*/ if (dcm.isHasPhase) { strcat (outname,"_ph"); //has phase map if (dcm.isHasMagnitude) @@ -2477,15 +2481,28 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc } if ((vol >= 0) && (hdr0.dim[4] > 1) && (vol < hdr0.dim[4])) { size_t imgsz4D = imgsz; - hdr0.dim[0] = 3; //3D + if (vol < 2) + hdr0.dim[0] = 3; //3D + int inVol = hdr0.dim[4]; hdr0.dim[4] = 1; - imgsz = nii_ImgBytes(hdr0); + size_t imgsz3D = nii_ImgBytes(hdr0); unsigned char *img4D = (unsigned char *)malloc(imgsz4D); memcpy(&img4D[0], &imgM[0], imgsz4D); free(imgM); - imgM = (unsigned char *)malloc(imgsz); - memcpy(&imgM[0], &img4D[vol * imgsz], imgsz); + imgM = (unsigned char *)malloc(imgsz3D * vol); + int outVol = 0; + for (int v = 0; v < inVol; v++) { + //dti4D cxc + if ((dti4D->gradDynVol[v] != 0) && (outVol < vol)) { + //printError("cxc -> %d %d %d\n", v, v * imgsz3D, imgsz4D); + memcpy(&imgM[outVol * imgsz3D], &img4D[v * imgsz3D], imgsz3D); + outVol ++; + } + } + hdr0.dim[4] = vol; + imgsz = nii_ImgBytes(hdr0); free(img4D); + saveAs3D = false; } char pathoutname[2048] = {""}; if (nii_createFilename(dcmList[dcmSort[0].indx], pathoutname, opts) == EXIT_FAILURE) { @@ -2623,8 +2640,11 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc #endif free(imgM); return returnCode;//EXIT_SUCCESS; -}// saveDcm2Nii() +}// saveDcm2NiiCore() +//int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmList[], struct TSearchList *nameList, struct TDCMopts opts, struct TDTI4D *dti4D) { +// return saveDcm2NiiCore(nConvert, dcmSort, dcmList, nameList, opts, dti4D, -1); +//} int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmList[], struct TSearchList *nameList, struct TDCMopts opts, struct TDTI4D *dti4D) { //this wrapper does nothing if all the images share the same echo time and scale // however, it segments images when these properties vary @@ -2635,7 +2655,10 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis printError("Unexpected error for image with varying echo time or intensity scaling\n"); return EXIT_FAILURE; } + int ret = EXIT_SUCCESS; + //check for repeated echoes - count unique number of echoes + /* //code below checks for multi-echoes - not required if maxNumberOfEchoes reported in PARREC int echoNum[kMaxDTI4D]; int echo = 1; for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) @@ -2648,11 +2671,50 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis echo++; echoNum[i] = echo; } + } + if (echo > 1) dcmList[indx].isMultiEcho = true;*/ + //check for repeated volumes + int seriesNum[kMaxDTI4D]; + int series = 1; + for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) + seriesNum[i] = 0; + seriesNum[0] = 1; + for (int i = 1; i < dcmList[indx].xyzDim[4]; i++) { + for (int j = 0; j < i; j++) + if ((dti4D->intenIntercept[i] == dti4D->intenIntercept[j]) && (dti4D->intenScale[i] == dti4D->intenScale[j]) && (dti4D->isPhase[i] == dti4D->isPhase[j]) && (dti4D->TE[i] == dti4D->TE[j])) seriesNum[i] = seriesNum[j]; + if (seriesNum[i] == 0) { + series++; + seriesNum[i] = series; + } + } + for (int s = 1; s <= series; s++) { + for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) + dti4D->gradDynVol[i] = 0; + int nVol = 0; + for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) + if (seriesNum[i] == s) { + dti4D->gradDynVol[i] = 1; + nVol ++; + dcmList[indx].TE = dti4D->TE[i]; + dcmList[indx].intenScale = dti4D->intenScale[i]; + dcmList[indx].intenIntercept = dti4D->intenIntercept[i]; + dcmList[indx].isHasPhase = dti4D->isPhase[i]; + dcmList[indx].intenScalePhilips = dti4D->intenScalePhilips[i]; + dcmList[indx].isHasMagnitude = false; + dcmList[indx].echoNum = echoNum[i]; + //dcmList[indx].gradDynVol = dti4D->gradDynVol[i]; + } + //printError("cxc files %d %d %g\n", nVol,echo, dcmList[indx].TE); + int ret2 = saveDcm2NiiCore(nConvert, dcmSort, dcmList, nameList, opts, dti4D, nVol); + if (ret2 != EXIT_SUCCESS) ret = ret2; //return EXIT_SUCCESS only if ALL are successful } - if (echo > 1) dcmList[indx].isMultiEcho = true; - for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) { - //printMessage("volume %d TE=%g Slope=%g Inter=%g Phase=%d\n", i, dti4D->TE[i], dti4D->intenScale[i], dti4D->intenIntercept[i], dti4D->isPhase[i] ); + // + //printError("cxc files %d %d %g\n", series,echo, dcmList[indx].TE); + //return ret; + + /*for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) { + printMessage("volume %d TE=%g Slope=%g Inter=%g Phase=%d\n", i, dti4D->TE[i], dti4D->intenScale[i], dti4D->intenIntercept[i], dti4D->isPhase[i] ); dcmList[indx].TE = dti4D->TE[i]; dcmList[indx].intenScale = dti4D->intenScale[i]; dcmList[indx].intenIntercept = dti4D->intenIntercept[i]; @@ -2660,13 +2722,12 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis dcmList[indx].intenScalePhilips = dti4D->intenScalePhilips[i]; dcmList[indx].isHasMagnitude = false; dcmList[indx].echoNum = echoNum[i]; - /*int echoNum = 1; - for (int i = 1; i < d.xyzDim[4]; i++) { - if (dti4D->TE[i-1] != dti4D->TE[i]) - }*/ - int ret = saveDcm2NiiCore(nConvert, dcmSort, dcmList, nameList, opts, dti4D, i); - } -} + dcmList[indx].gradDynVol = dti4D->gradDynVol[i]; + int ret2 = saveDcm2NiiCore(nConvert, dcmSort, dcmList, nameList, opts, dti4D, i); + if (ret2 != EXIT_SUCCESS) ret = ret2; //return EXIT_SUCCESS only if ALL are successful + }*/ + return ret; +}// saveDcm2Nii() void fillTDCMsort(struct TDCMsort& tdcmref, const uint64_t indx, const struct TDICOMdata& dcmdata){ // Copy the relevant parts of dcmdata to tdcmref. @@ -2978,7 +3039,7 @@ int convert_parRec(struct TDCMopts opts) { nameList.str[0] = (char *)malloc(strlen(opts.indir)+1); strcpy(nameList.str[0],opts.indir); TDTI4D dti4D; - dcmList[0] = nii_readParRec(nameList.str[0], opts.isVerbose, &dti4D); + dcmList[0] = nii_readParRec(nameList.str[0], opts.isVerbose, &dti4D, false); struct TDCMsort dcmSort[1]; dcmSort[0].indx = 0; if (dcmList[0].isValid) From 66ae53c5c294f69a639b271bc6889f408fe21c31 Mon Sep 17 00:00:00 2001 From: Chris Rorden Date: Sat, 24 Mar 2018 09:50:22 -0400 Subject: [PATCH 13/14] Further PARREC tuning --- console/nii_dicom.cpp | 8 ++++---- console/nii_dicom_batch.cpp | 6 ++---- 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index cbee2c49..7d619903 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -1752,7 +1752,8 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D } //even is sequential, we need to decode if multiple echoes/phasemaps,etc //if (!isSlicesSpatiallySequentialPhilips) { - printWarning("PAR file order of slices is not sequential: please ensure correct slice re-ordering\n"); + //printWarning("PAR file order of slices is not sequential: please ensure correct slice re-ordering\n"); + d.isScaleOrTEVaries = true; if (numSlice2D > kMaxSlice2D) { printError("Overloaded slice re-ordering. Number of slices (%d) exceeds kMaxSlice2D (%d)\n", numSlice2D, kMaxSlice2D); dti4D->sliceOrder[0] = -1; @@ -1893,8 +1894,6 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D } if (d.isScaleOrTEVaries) printWarning("Varying dimensions (echoes, phase maps, intensity scaling) will require volumes to be saved separately (hint: you may prefer dicm2nii output)\n"); - - } //check DTI makes sense if (d.CSA.numDti > 1) { @@ -1909,7 +1908,8 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D if ((!v1varies) || (!v2varies) || (!v3varies)) printError("Bizarre b-vectors %s\n", parname); } - if ((maxEcho > 1) || (maxCardiac > 1)) printWarning("Multiple Echo (%d) or Cardiac (%d). Segment output, e.g. nii_segment4d('img.nii', %d)\n", maxEcho, maxCardiac, maxEcho*maxCardiac); + if ((maxEcho > 1) || (maxCardiac > 1)) printWarning("Multiple Echo (%d) or Cardiac (%d). Carefully inspect output\n", maxEcho, maxCardiac); + if ((maxEcho > 1) || (maxCardiac > 1)) d.isScaleOrTEVaries = true; return d; } //nii_readParRec() diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 153f8ba2..bd782b5d 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -2656,9 +2656,8 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis return EXIT_FAILURE; } int ret = EXIT_SUCCESS; - //check for repeated echoes - count unique number of echoes - /* //code below checks for multi-echoes - not required if maxNumberOfEchoes reported in PARREC + //code below checks for multi-echoes - not required if maxNumberOfEchoes reported in PARREC int echoNum[kMaxDTI4D]; int echo = 1; for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) @@ -2672,7 +2671,7 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis echoNum[i] = echo; } } - if (echo > 1) dcmList[indx].isMultiEcho = true;*/ + if (echo > 1) dcmList[indx].isMultiEcho = true; //check for repeated volumes int seriesNum[kMaxDTI4D]; int series = 1; @@ -2687,7 +2686,6 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis seriesNum[i] = series; } } - for (int s = 1; s <= series; s++) { for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) dti4D->gradDynVol[i] = 0; From dd527d986c50bcc7fdd5980e6601382137c4d641 Mon Sep 17 00:00:00 2001 From: Chris Rorden Date: Sun, 25 Mar 2018 23:24:46 -0400 Subject: [PATCH 14/14] Improved GE support --- console/nii_dicom.cpp | 9 +- console/nii_dicom.h | 4 +- console/nii_dicom_batch.cpp | 225 +++++++++++++++++++++++++++--------- 3 files changed, 178 insertions(+), 60 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 7d619903..93c5a8c1 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -719,6 +719,8 @@ struct TDICOMdata clear_dicom_data() { d.echoTrainLength = 0; d.phaseFieldofView = 0.0; d.dwellTime = 0; + d.protocolBlockStartGE = 0; + d.protocolBlockLengthGE = 0; d.phaseEncodingSteps = 0; d.coilNum = 0; d.accelFactPE = 0.0; @@ -3444,7 +3446,7 @@ double TE = 0.0; //most recent echo time recorded uint32_t groupElement; long lPos = 0; bool isPhilipsDerived = false; - bool isPhilipsDiffusion = false; + //bool isPhilipsDiffusion = false; if (isPart10prefix) { //for part 10 files, skip preamble and prefix lPos = 128+4; //4-byte signature starts at 128 groupElement = buffer[lPos] | (buffer[lPos+1] << 8) | (buffer[lPos+2] << 16) | (buffer[lPos+3] << 24); @@ -4072,8 +4074,9 @@ double TE = 0.0; //most recent echo time recorded break; case kProtocolDataBlockGE : if (d.manufacturer != kMANUFACTURER_GE) break; - //printError("ProtocolDataBlockGE to do %d %d %2x%2x\n", lPos, lLength, (uint8_t)buffer[lPos], (uint8_t)buffer[lPos+1]); - //readProtocolDataBlockGE(&buffer[lPos], lLength, &d.CSA, isVerbose); //, dti4D); + d.protocolBlockLengthGE = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); + d.protocolBlockStartGE = lPos+lFileOffset+4; + //printError("ProtocolDataBlockGE %d @ %d\n", d.protocolBlockLengthGE, d.protocolBlockStartGE); break; case kDoseCalibrationFactor : d.doseCalibrationFactor = dcmStrFloat(lLength, &buffer[lPos]); diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 064c9d01..6fea7b85 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -38,7 +38,7 @@ extern "C" { #define kCCsuf " CompilerNA" //unknown compiler! #endif -#define kDCMvers "v1.0.20180307" kDCMsuf kCCsuf +#define kDCMvers "v1.0.20180325" kDCMsuf kCCsuf static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic static const int kMaxDTI4D = 4096; //maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images @@ -127,7 +127,7 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; //numberOfDynamicScans, patientPositionNumPhilips //patientPositionSequentialRepeats,patientPositionRepeats, //maxGradDynVol, gradDynVol, - int modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, echoTrainLength, coilNum, echoNum, sliceOrient, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,locationsInAcquisition, compressionScheme; + int protocolBlockStartGE, protocolBlockLengthGE, modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, echoTrainLength, coilNum, echoNum, sliceOrient, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,locationsInAcquisition, compressionScheme; float patientWeight, zSpacing, zThick, pixelBandwidth, SAR, phaseFieldofView, accelFactPE, flipAngle, fieldStrength, TE, TI, TR, intenScale, intenIntercept, intenScalePhilips, gantryTilt, lastScanLoc, angulation[4]; float orient[7], patientPosition[4], patientPositionLast[4], xyzMM[4], stackOffcentre[4]; float radionuclidePositronFraction, radionuclideTotalDose, radionuclideHalfLife, doseCalibrationFactor; //PET ISOTOPE MODULE ATTRIBUTES (C.8-57) diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index bd782b5d..15144a0e 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -407,6 +407,28 @@ const void * memmem(const char *l, size_t l_len, const char *s, size_t s_len) { //n.b. memchr returns "const void *" not "void *" for Windows C++ https://msdn.microsoft.com/en-us/library/d7zdhf37.aspx #endif //for systems without memmem +/*int readKeyX(const char * key, char * buffer, int remLength) { //look for text key in binary data stream, return subsequent integer value + int ret = 0; + char *keyPos = (char *)memmem(buffer, remLength, key, strlen(key)); + printWarning("<><><>\n"); + if (!keyPos) return 0; + printWarning("<><> %d\n", strlen(keyPos)); + int i = (int)strlen(key); + int numDigits = 0; + while( ( i< remLength) && (numDigits >= 0) ) { + printMessage("%c", keyPos[i]); + if( keyPos[i] >= '0' && keyPos[i] <= '9' ) { + ret = (10 * ret) + keyPos[i] - '0'; + numDigits ++; + } else if (numDigits > 0) + numDigits = -1; + i++; + } + printWarning("---> %d\n", ret); + return ret; +} //readKey() +*/ + int readKey(const char * key, char * buffer, int remLength) { //look for text key in binary data stream, return subsequent integer value int ret = 0; char *keyPos = (char *)memmem(buffer, remLength, key, strlen(key)); @@ -609,6 +631,81 @@ void siemensCsaAscii(const char * filename, int csaOffset, int csaLength, float } // siemensCsaAscii() #endif //myReadAsciiCsa() +#ifndef myDisableMiniZ + #define myReadGeProtocolBlock +#endif +#ifdef myReadGeProtocolBlock +int geProtocolBlock(const char * filename, int geOffset, int geLength, int isVerbose, int* sliceOrder, int* viewOrder) { + *sliceOrder = 0; + *viewOrder = 0; + int ret = EXIT_FAILURE; + if ((geOffset < 0) || (geLength < 20)) return ret; + FILE * pFile = fopen ( filename, "rb" ); + if(pFile==NULL) return ret; + fseek (pFile , 0 , SEEK_END); + long lSize = ftell (pFile); + if (lSize < (geOffset+geLength)) { + fclose (pFile); + return ret; + } + fseek(pFile, geOffset, SEEK_SET); + mz_uint8 * pCmp = (mz_uint8*) malloc (geLength); + if(pCmp == NULL) return ret; + size_t result = fread (pCmp,1,geLength,pFile); + if ((int)result != geLength) return ret; + int cmpSz = geLength; + //http://www.forensicswiki.org/wiki/Gzip + // always little endia! http://www.onicos.com/staff/iz/formats/gzip.html + if (cmpSz < 20) return ret; + if ((pCmp[0] != 31) || (pCmp[1] != 139) || (pCmp[2] != 8)) return ret; //check signature and deflate algorithm + uint8_t flags = pCmp[3]; + bool isFNAME = ((flags & 0x08) == 0x08); + bool isFCOMMENT = ((flags & 0x10) == 0x10); + size_t hdrSz = 10; + if (isFNAME) {//skip null-terminated string FNAME + for (hdrSz = hdrSz; hdrSz < cmpSz; hdrSz++) + if (pCmp[hdrSz] == 0) break; + hdrSz++; + } + if (isFCOMMENT) {//skip null-terminated string COMMENT + for (hdrSz = hdrSz; hdrSz < cmpSz; hdrSz++) + if (pCmp[hdrSz] == 0) break; + hdrSz++; + } + size_t unCmpSz = ((size_t)pCmp[cmpSz-4])+((size_t)pCmp[cmpSz-3] << 8)+((size_t)pCmp[cmpSz-2] << 16)+((size_t)pCmp[cmpSz-1] << 24); + //printf(">> %d %d %zu %zu %zu\n", isFNAME, isFCOMMENT, cmpSz, unCmpSz, hdrSz); + + z_stream s; + memset (&s, 0, sizeof (z_stream)); + inflateInit2(&s, -MZ_DEFAULT_WINDOW_BITS); + mz_uint8 *pUnCmp = (mz_uint8 *)malloc((size_t)unCmpSz); + s.avail_out = unCmpSz; + s.next_in = pCmp+ hdrSz; + s.avail_in = cmpSz-hdrSz-8; + s.next_out = (uint8_t *) pUnCmp; + ret = mz_inflate(&s, MZ_SYNC_FLUSH); + if (ret != MZ_STREAM_END) { + free(pUnCmp); + return EXIT_FAILURE; + } + //https://groups.google.com/forum/#!msg/comp.protocols.dicom/mxnCkv8A-i4/W_uc6SxLwHQJ + // DISCOVERY MR750 / 24\MX\MR Software release:DV24.0_R01_1344.a) are now storing an XML file + // + if ((pUnCmp[0] == '<') && (pUnCmp[1] == '?')) + printWarning("New XML-based GE Protocol Block is not yet supported: please report issue on dcm2niix Github page\n"); + char keyStrSO[] = "SLICEORDER"; + *sliceOrder = readKey(keyStrSO, (char *) pUnCmp, unCmpSz); + char keyStrVO[] = "VIEWORDER"; //"MATRIXX"; + *viewOrder = readKey(keyStrVO, (char *) pUnCmp, unCmpSz); + if (isVerbose > 1) { + printMessage("GE Protocol Block %s bytes %d compressed, %d uncompressed @ %d\n", filename, geLength, unCmpSz, geOffset); + printMessage(" ViewOrder %d SliceOrder %d\n", *viewOrder, *sliceOrder); + printMessage("%s\n", pUnCmp); + } + free(pUnCmp); + return EXIT_SUCCESS; +} +#endif //myReadGeProtocolBlock() void json_Str(FILE *fp, const char *sLabel, char *sVal) { if (strlen(sVal) < 1) return; //fprintf(fp, sLabel, sVal ); @@ -792,6 +889,17 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, float pf = 1.0f; //partial fourier bool interp = false; //2D interpolation float phaseOversampling = 0.0; + int viewOrderGE = -1; + #ifdef myReadGeProtocolBlock + int sliceOrderGE = -1; + if ((d.manufacturer == kMANUFACTURER_GE) && (d.protocolBlockStartGE> 0) && (d.protocolBlockLengthGE > 19)) { + printWarning("Using GE Protocol Data Block for BIDS data (beware: new feature)\n"); + int ok = geProtocolBlock(filename, d.protocolBlockStartGE, d.protocolBlockLengthGE, opts.isVerbose, &sliceOrderGE, &viewOrderGE); + if (ok != EXIT_SUCCESS) + printWarning("Unable to decode GE protocol block\n"); + printMessage(" ViewOrder %d SliceOrder %d\n", viewOrderGE, sliceOrderGE); + } //read protocolBlockGE + #endif #ifdef myReadAsciiCsa if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.CSA.SeriesHeader_offset > 0) && (d.CSA.SeriesHeader_length > 0)) { int baseResolution, interpInt, partialFourier, echoSpacing, parallelReductionFactorInPlane; @@ -913,7 +1021,10 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.dwellTime > 0)) fprintf(fp, "\t\"DwellTime\": %g,\n", d.dwellTime * 1E-9); // Phase encoding polarity - if (((d.phaseEncodingRC == 'R') || (d.phaseEncodingRC == 'C')) && (!d.is3DAcq) && ((d.CSA.phaseEncodingDirectionPositive == 1) || (d.CSA.phaseEncodingDirectionPositive == 0))) { + int phPos = d.CSA.phaseEncodingDirectionPositive; + if (viewOrderGE > -1) + phPos = viewOrderGE; + if (((d.phaseEncodingRC == 'R') || (d.phaseEncodingRC == 'C')) && (!d.is3DAcq) && (phPos >= 0)) { if (d.phaseEncodingRC == 'C') //Values should be "R"ow, "C"olumn or "?"Unknown fprintf(fp, "\t\"PhaseEncodingDirection\": \"j"); else if (d.phaseEncodingRC == 'R') @@ -924,17 +1035,40 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, //However, DICOM and NIfTI are reversed in the j (ROW) direction //Equivalent to dicm2nii's "if flp(iPhase), phPos = ~phPos; end" //for samples see https://github.com/rordenlab/dcm2niix/issues/125 - if (d.CSA.phaseEncodingDirectionPositive == -1) + if (phPos < 0) fprintf(fp, "?"); //unknown - else if ((d.CSA.phaseEncodingDirectionPositive == 0) && (d.phaseEncodingRC != 'C')) + else if ((phPos == 0) && (d.phaseEncodingRC != 'C')) fprintf(fp, "-"); - else if ((d.phaseEncodingRC == 'C') && (d.CSA.phaseEncodingDirectionPositive == 1) && (opts.isFlipY)) + else if ((d.phaseEncodingRC == 'C') && (phPos == 1) && (opts.isFlipY)) fprintf(fp, "-"); - else if ((d.phaseEncodingRC == 'C') && (d.CSA.phaseEncodingDirectionPositive == 0) && (!opts.isFlipY)) + else if ((d.phaseEncodingRC == 'C') && (phPos == 0) && (!opts.isFlipY)) fprintf(fp, "-"); fprintf(fp, "\",\n"); } //only save PhaseEncodingDirection if BOTH direction and POLARITY are known - // Slice Timing + // Slice Timing GE + if ((viewOrderGE > -1) && (h->dim[3] > 1) && (h->dim[4] > 1) && (d.TR > 0)) { // + //Warning: not correct for multiband sequences... not sure how these are stored + //Warning: will not create correct times for sparse acquisitions where DelayTimeInTR > 0 + float t = d.TR/ (float)h->dim[3] ; + fprintf(fp, "\t\"SliceTiming\": [\n"); + if (viewOrderGE == 1) {//interleaved ascending + for (int i = 0; i < h->dim[3]; i++) { + if (i != 0) + fprintf(fp, ",\n"); + int s = (i / 2); + if ((i % 2) != 0) s += (h->dim[3]+1)/2; + fprintf(fp, "\t\t%g", (float) s * t / 1000.0 ); + } + } else { //sequential ascending + for (int i = 0; i < h->dim[3]; i++) { + if (i != 0) + fprintf(fp, ",\n"); + fprintf(fp, "\t\t%g", (float) i * t / 1000.0 ); + } + } + fprintf(fp, "\t],\n"); + } + //Slice Timing Siemens if (d.CSA.sliceTiming[0] >= 0.0) { fprintf(fp, "\t\"SliceTiming\": [\n"); if (d.CSA.protocolSliceNumber1 > 1) { @@ -2337,7 +2471,7 @@ void checkSliceTiming(struct TDICOMdata * d, struct TDICOMdata * d1) { printMessage("CSA slice timing based on 2nd volume, 1st volume corrupted (CMRR bug, range %g..%g, TR=%gms)\n", minT, maxT, d->TR); }//checkSliceTiming -int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmList[], struct TSearchList *nameList, struct TDCMopts opts, struct TDTI4D *dti4D, int vol) { +int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmList[], struct TSearchList *nameList, struct TDCMopts opts, struct TDTI4D *dti4D, int segVol) { bool iVaries = intensityScaleVaries(nConvert,dcmSort,dcmList); float *sliceMMarray = NULL; //only used if slices are not equidistant uint64_t indx = dcmSort[0].indx; @@ -2448,7 +2582,6 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc hdr0.dim[5] = nConvert; hdr0.dim[0] = 5; } - /*if (nConvert > 1) { //next determine if TR is true time between volumes double startTime = dcmList[indx0].acquisitionTime; double endTime = startTime; @@ -2479,27 +2612,33 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc free(img); } } - if ((vol >= 0) && (hdr0.dim[4] > 1) && (vol < hdr0.dim[4])) { + if ((segVol >= 0) && (hdr0.dim[4] > 1)) { + int inVol = hdr0.dim[4]; + int nVol = 0; + for (int v = 0; v < inVol; v++) + if (dti4D->gradDynVol[v] == segVol) + nVol ++; + if (nVol < 1) { + printError("Series %d does not exist\n", segVol); + return EXIT_FAILURE; + } size_t imgsz4D = imgsz; - if (vol < 2) + if (nVol < 2) hdr0.dim[0] = 3; //3D - int inVol = hdr0.dim[4]; hdr0.dim[4] = 1; size_t imgsz3D = nii_ImgBytes(hdr0); unsigned char *img4D = (unsigned char *)malloc(imgsz4D); memcpy(&img4D[0], &imgM[0], imgsz4D); free(imgM); - imgM = (unsigned char *)malloc(imgsz3D * vol); + imgM = (unsigned char *)malloc(imgsz3D * nVol); int outVol = 0; for (int v = 0; v < inVol; v++) { - //dti4D cxc - if ((dti4D->gradDynVol[v] != 0) && (outVol < vol)) { - //printError("cxc -> %d %d %d\n", v, v * imgsz3D, imgsz4D); + if ((dti4D->gradDynVol[v] == segVol) && (outVol < nVol)) { memcpy(&imgM[outVol * imgsz3D], &img4D[v * imgsz3D], imgsz3D); outVol ++; } } - hdr0.dim[4] = vol; + hdr0.dim[4] = nVol; imgsz = nii_ImgBytes(hdr0); free(img4D); saveAs3D = false; @@ -2517,8 +2656,8 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc for(int i = 0; i < nConvert; ++i) dcmList[dcmSort[i].indx].converted2NII = 1; if (opts.numSeries < 0) { //report series number but do not convert - if (vol >= 0) - printMessage("\t%ld.%d\t%s\n", dcmList[dcmSort[0].indx].seriesNum, vol, pathoutname); + if (segVol >= 0) + printMessage("\t%ld.%d\t%s\n", dcmList[dcmSort[0].indx].seriesNum, segVol-1, pathoutname); else printMessage("\t%ld\t%s\n", dcmList[dcmSort[0].indx].seriesNum, pathoutname); printMessage(" %s\n",nameList->str[dcmSort[0].indx]); @@ -2542,8 +2681,8 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc if (opts.numSeries > 0) { int i = 0; float seriesNum = (float) dcmList[dcmSort[0].indx].seriesNum; - if (vol > 0) - seriesNum = seriesNum + (float) vol / 10.0; //n.b. we will have problems if vol > 9. However, 9 distinct TEs/scalings/PhaseMag seems unlikely + if (segVol > 0) + seriesNum = seriesNum + ((float) segVol - 1.0) / 10.0; //n.b. we will have problems if segVol > 9. However, 9 distinct TEs/scalings/PhaseMag seems unlikely for (; i < opts.numSeries; i++) { if (isSameFloatGE(opts.seriesNumber[i], seriesNum)) { //if (opts.seriesNumber[i] == dcmList[dcmSort[0].indx].seriesNum) { @@ -2642,9 +2781,6 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc return returnCode;//EXIT_SUCCESS; }// saveDcm2NiiCore() -//int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmList[], struct TSearchList *nameList, struct TDCMopts opts, struct TDTI4D *dti4D) { -// return saveDcm2NiiCore(nConvert, dcmSort, dcmList, nameList, opts, dti4D, -1); -//} int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmList[], struct TSearchList *nameList, struct TDCMopts opts, struct TDTI4D *dti4D) { //this wrapper does nothing if all the images share the same echo time and scale // however, it segments images when these properties vary @@ -2673,27 +2809,24 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis } if (echo > 1) dcmList[indx].isMultiEcho = true; //check for repeated volumes - int seriesNum[kMaxDTI4D]; int series = 1; for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) - seriesNum[i] = 0; - seriesNum[0] = 1; + dti4D->gradDynVol[i] = 0; + dti4D->gradDynVol[0] = 1; for (int i = 1; i < dcmList[indx].xyzDim[4]; i++) { for (int j = 0; j < i; j++) - if ((dti4D->intenIntercept[i] == dti4D->intenIntercept[j]) && (dti4D->intenScale[i] == dti4D->intenScale[j]) && (dti4D->isPhase[i] == dti4D->isPhase[j]) && (dti4D->TE[i] == dti4D->TE[j])) seriesNum[i] = seriesNum[j]; - if (seriesNum[i] == 0) { + if ((dti4D->intenIntercept[i] == dti4D->intenIntercept[j]) && (dti4D->intenScale[i] == dti4D->intenScale[j]) && (dti4D->isPhase[i] == dti4D->isPhase[j]) && (dti4D->TE[i] == dti4D->TE[j])) + dti4D->gradDynVol[i] = dti4D->gradDynVol[j]; + if (dti4D->gradDynVol[i] == 0) { series++; - seriesNum[i] = series; + dti4D->gradDynVol[i] = series; } } for (int s = 1; s <= series; s++) { for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) - dti4D->gradDynVol[i] = 0; - int nVol = 0; - for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) - if (seriesNum[i] == s) { - dti4D->gradDynVol[i] = 1; - nVol ++; + if (dti4D->gradDynVol[i] == s) { + //dti4D->gradDynVol[i] = s; + //nVol ++; dcmList[indx].TE = dti4D->TE[i]; dcmList[indx].intenScale = dti4D->intenScale[i]; dcmList[indx].intenIntercept = dti4D->intenIntercept[i]; @@ -2701,29 +2834,11 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis dcmList[indx].intenScalePhilips = dti4D->intenScalePhilips[i]; dcmList[indx].isHasMagnitude = false; dcmList[indx].echoNum = echoNum[i]; - //dcmList[indx].gradDynVol = dti4D->gradDynVol[i]; + break; } - //printError("cxc files %d %d %g\n", nVol,echo, dcmList[indx].TE); - int ret2 = saveDcm2NiiCore(nConvert, dcmSort, dcmList, nameList, opts, dti4D, nVol); + int ret2 = saveDcm2NiiCore(nConvert, dcmSort, dcmList, nameList, opts, dti4D, s); if (ret2 != EXIT_SUCCESS) ret = ret2; //return EXIT_SUCCESS only if ALL are successful } - // - //printError("cxc files %d %d %g\n", series,echo, dcmList[indx].TE); - //return ret; - - /*for (int i = 0; i < dcmList[indx].xyzDim[4]; i++) { - printMessage("volume %d TE=%g Slope=%g Inter=%g Phase=%d\n", i, dti4D->TE[i], dti4D->intenScale[i], dti4D->intenIntercept[i], dti4D->isPhase[i] ); - dcmList[indx].TE = dti4D->TE[i]; - dcmList[indx].intenScale = dti4D->intenScale[i]; - dcmList[indx].intenIntercept = dti4D->intenIntercept[i]; - dcmList[indx].isHasPhase = dti4D->isPhase[i]; - dcmList[indx].intenScalePhilips = dti4D->intenScalePhilips[i]; - dcmList[indx].isHasMagnitude = false; - dcmList[indx].echoNum = echoNum[i]; - dcmList[indx].gradDynVol = dti4D->gradDynVol[i]; - int ret2 = saveDcm2NiiCore(nConvert, dcmSort, dcmList, nameList, opts, dti4D, i); - if (ret2 != EXIT_SUCCESS) ret = ret2; //return EXIT_SUCCESS only if ALL are successful - }*/ return ret; }// saveDcm2Nii()