diff --git a/console/makefile b/console/makefile index 3fe5ea46..4b0c688b 100755 --- a/console/makefile +++ b/console/makefile @@ -1,9 +1,14 @@ -CFLAGS=-s +# Regular use +CFLAGS=-s -O3 + +# Debugging +#CFLAGS=-g + ifneq ($(OS),Windows_NT) OS = $(shell uname) ifeq "$(OS)" "Darwin" - CFLAGS=-dead_strip + CFLAGS=-dead_strip -O3 endif endif all: - g++ $(CFLAGS) -O3 -I. main_console.cpp nii_foreign.cpp nii_dicom.cpp jpg_0XC3.cpp ujpeg.cpp nifti1_io_core.cpp nii_ortho.cpp nii_dicom_batch.cpp -o dcm2niix -DmyDisableOpenJPEG + g++ $(CFLAGS) -I. main_console.cpp nii_foreign.cpp nii_dicom.cpp jpg_0XC3.cpp ujpeg.cpp nifti1_io_core.cpp nii_ortho.cpp nii_dicom_batch.cpp -o dcm2niix -DmyDisableOpenJPEG diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index ce470f84..03fc18cb 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -657,6 +657,8 @@ struct TDICOMdata clear_dicom_data() { d.angulation[i] = 0.0f; d.xyzMM[i] = 1; } + for (int i=0; i < MAX_NUMBER_OF_DIMENSIONS; ++i) + d.dimensionIndexValues[i] = 0; d.CSA.sliceTiming[0] = -1.0f; //impossible value denotes not known d.CSA.numDti = 0; for (int i=0; i < 5; i++) @@ -888,7 +890,8 @@ float dcmFloat(int lByteLength, unsigned char lBuffer[], bool littleEndian) {//r return swapVal; } //dcmFloat() -double dcmFloatDouble(int lByteLength, unsigned char lBuffer[], bool littleEndian) {//read binary 32-bit float +double dcmFloatDouble(const size_t lByteLength, const unsigned char lBuffer[], + const bool littleEndian) {//read binary 32-bit float //http://stackoverflow.com/questions/2782725/converting-float-values-from-big-endian-to-little-endian bool swap = (littleEndian != littleEndianPlatform()); double retVal = 0.0f; @@ -921,14 +924,14 @@ int dcmInt (int lByteLength, unsigned char lBuffer[], bool littleEndian) { //rea return lBuffer[3]+(lBuffer[2]<<8)+(lBuffer[1]<<16)+(lBuffer[0]<<24); //shortint vs word? } //dcmInt() -int dcmStrInt (int lByteLength, unsigned char lBuffer[]) {//read float stored as a string +int dcmStrInt (const int lByteLength, const unsigned char lBuffer[]) {//read float stored as a string //#ifdef _MSC_VER char * cString = (char *)malloc(sizeof(char) * (lByteLength + 1)); //#else // char cString[lByteLength + 1]; //#endif cString[lByteLength] =0; - memcpy(cString, (char*)&lBuffer[0], lByteLength); + memcpy(cString, (const unsigned char*)(&lBuffer[0]), lByteLength); //printMessage(" --> *%s* %s%s\n",cString, &lBuffer[0],&lBuffer[1]); int ret = atoi(cString); //#ifdef _MSC_VER @@ -937,7 +940,7 @@ int dcmStrInt (int lByteLength, unsigned char lBuffer[]) {//read float stored as return ret; } //dcmStrInt() -int dcmStrManufacturer (int lByteLength, unsigned char lBuffer[]) {//read float stored as a string +int dcmStrManufacturer (const int lByteLength, unsigned char lBuffer[]) {//read float stored as a string if (lByteLength < 2) return kMANUFACTURER_UNKNOWN; //#ifdef _MSC_VER char * cString = (char *)malloc(sizeof(char) * (lByteLength + 1)); @@ -1171,6 +1174,16 @@ void dcmMultiShorts (int lByteLength, unsigned char lBuffer[], int lnShorts, uin nifti_swap_2bytes(lnShorts, &lShorts[0]); } //dcmMultiShorts() +void dcmMultiLongs (int lByteLength, unsigned char lBuffer[], int lnLongs, uint32_t *lLongs, bool littleEndian) { + //read array of unsigned longs UL http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_6.2.html + if((lnLongs < 1) || (lByteLength != (lnLongs * 4))) + return; + memcpy(&lLongs[0], (uint32_t *)&lBuffer[0], lByteLength); + bool swap = (littleEndian != littleEndianPlatform()); + if (swap) + nifti_swap_4bytes(lnLongs, &lLongs[0]); +} //dcmMultiLongs() + void dcmMultiFloat (int lByteLength, char lBuffer[], int lnFloats, float *lFloats) { //warning: lFloats indexed from 1! will fill lFloats[1]..[nFloats] if ((lnFloats < 1) || (lByteLength < 1)) return; @@ -1205,7 +1218,7 @@ void dcmMultiFloat (int lByteLength, char lBuffer[], int lnFloats, float *lFloat //#endif } //dcmMultiFloat() -float dcmStrFloat (int lByteLength, unsigned char lBuffer[]) { //read float stored as a string +float dcmStrFloat (const int lByteLength, const unsigned char lBuffer[]) { //read float stored as a string //#ifdef _MSC_VER char * cString = (char *)malloc(sizeof(char) * (lByteLength + 1)); //#else @@ -2669,6 +2682,196 @@ int isDICOMfile(const char * fname) { //0=NotDICOM, 1=DICOM, 2=Maybe(not Part 10 return 0; } //isDICOMfile() +//START RIR 12/2017 Robert I. Reid + +// Gathering spot for all the info needed to get the b value and direction +// for a volume. +struct TVolumeDiffusion { + struct TDICOMdata* pdd; // The multivolume + struct TDTI4D* pdti4D; // permanent records. + + uint8_t manufacturer; // kMANUFACTURER_UNKNOWN, kMANUFACTURER_SIEMENS, etc. + + //void set_manufacturer(const uint8_t m) {manufacturer = m; update();} // unnecessary + + // Everything after this in the structure would be private if it were a C++ + // class, but it has been rewritten as a struct for C compatibility. I am + // using _ as a hint of that, although _ for privacy is not really a + // universal convention in C. Privacy is desired because immediately + // any of these are updated _update_tvd() should be called. + + bool _isAtFirstPatientPosition; // Limit b vals and vecs to 1 per volume. + + //float bVal0018_9087; // kDiffusion_b_value, always present in Philips/Siemens. + //float bVal2001_1003; // kDiffusionBFactor + // float dirRL2005_10b0; // kDiffusionDirectionRL + // float dirAP2005_10b1; // kDiffusionDirectionAP + // float dirFH2005_10b2; // kDiffusionDirectionFH + + // Philips diffusion scans tend to have a "trace" (average of the diffusion + // weighted volumes) volume tacked on, usually but not always at the end, + // so b is > 0, but the direction is meaningless. Most software versions + // explicitly set the direction to 0, but version 3 does not, making (0x18, + // 0x9075) necessary. + bool _isPhilipsNonDirectional; + + //char _directionality0018_9075[16]; // DiffusionDirectionality, not in Philips 2.6. + // float _orientation0018_9089[3]; // kDiffusionOrientation, always + // // present in Philips/Siemens for + // // volumes with a direction. + //char _seq0018_9117[64]; // MRDiffusionSequence, not in Philips 2.6. + + float _dtiV[4]; + //uint16_t numDti; +}; +struct TVolumeDiffusion initTVolumeDiffusion(struct TDICOMdata* ptdd, struct TDTI4D* dti4D); +void clear_volume(struct TVolumeDiffusion* ptvd); // Blank the volume-specific members or set them to impossible values. +void set_directionality0018_9075(struct TVolumeDiffusion* ptvd, unsigned char* inbuf); +void set_orientation0018_9089(struct TVolumeDiffusion* ptvd, int lLength, unsigned char* inbuf, + bool isLittleEndian); +void set_isAtFirstPatientPosition_tvd(struct TVolumeDiffusion* ptvd, bool iafpp); +void set_bValGE(struct TVolumeDiffusion* ptvd, int lLength, unsigned char* inbuf); +void set_diffusion_directionGE(struct TVolumeDiffusion* ptvd, int lLength, unsigned char* inbuf, int axis); +void set_bVal(struct TVolumeDiffusion* ptvd, float b); +void _update_tvd(struct TVolumeDiffusion* ptvd); + +struct TVolumeDiffusion initTVolumeDiffusion(struct TDICOMdata* ptdd, struct TDTI4D* dti4D) { + struct TVolumeDiffusion tvd; + tvd.pdd = ptdd; + tvd.pdti4D = dti4D; + clear_volume(&tvd); + return tvd; +} //initTVolumeDiffusion() + +void clear_volume(struct TVolumeDiffusion* ptvd) { + ptvd->_isAtFirstPatientPosition = false; + ptvd->manufacturer = kMANUFACTURER_UNKNOWN; + //bVal0018_9087 = -1; + //ptvd->_directionality0018_9075[0] = 0; + //ptvd->seq0018_9117[0] = 0; + //bVal2001_1003 = -1; + // dirRL2005_10b0 = 2; + // dirAP2005_10b1 = 2; + // dirFH2005_10b2 = 2; + ptvd->_isPhilipsNonDirectional = false; + ptvd->_dtiV[0] = -1; + for(int i = 1; i < 4; ++i) + ptvd->_dtiV[i] = 2; + //numDti = 0; +}//clear_volume() + +void set_directionality0018_9075(struct TVolumeDiffusion* ptvd, unsigned char* inbuf) { + if(strncmp(( char*)(inbuf), "DIRECTIONAL", 11) && // strncmp = 0 for ==. + strncmp(( char*)(inbuf), "BMATRIX", 7)){ // Siemens XA10 + ptvd->_isPhilipsNonDirectional = true; + // Explicitly set the direction to 0 now, because there may + // not be a 0018,9089 for this frame. + for(int i = 1; i < 4; ++i) // 1-3 is intentional. + ptvd->_dtiV[i] = 0.0; + } + else{ + ptvd->_isPhilipsNonDirectional = false; + // Wait for 0018,9089 to get the direction. + } + + _update_tvd(ptvd); +} //set_directionality0018_9075() + +void set_bValGE(struct TVolumeDiffusion* ptvd, int lLength, unsigned char* inbuf) { + //int dcmStrInt (int lByteLength, unsigned char lBuffer[]) {//read float stored as a string +dcmStrInt(lLength, inbuf); + + ptvd->_dtiV[0] = dcmStrInt(lLength, inbuf); + //dd.CSA.numDti = 1; // Always true for GE. + + _update_tvd(ptvd); +} //set_bValGE() + +// axis: 0 -> x, 1 -> y , 2 -> z +void set_diffusion_directionGE(struct TVolumeDiffusion* ptvd, int lLength, unsigned char* inbuf, const int axis){ + ptvd->_dtiV[axis + 1] = dcmStrFloat(lLength, inbuf); + + _update_tvd(ptvd); +}//set_diffusion_directionGE() + +void dcmMultiFloatDouble (size_t lByteLength, unsigned char lBuffer[], size_t lnFloats, float *lFloats, bool isLittleEndian) { + size_t floatlen = lByteLength / lnFloats; + + for(size_t i = 0; i < lnFloats; ++i) + lFloats[i] = dcmFloatDouble((int)floatlen, lBuffer + i * floatlen, isLittleEndian); +} //dcmMultiFloatDouble() + +void set_orientation0018_9089(struct TVolumeDiffusion* ptvd, int lLength, unsigned char* inbuf, bool isLittleEndian) { + if(ptvd->_isPhilipsNonDirectional){ + for(int i = 1; i < 4; ++i) // Deliberately ignore inbuf; it might be nonsense. + ptvd->_dtiV[i] = 0.0; + } + else + dcmMultiFloatDouble(lLength, inbuf, 3, ptvd->_dtiV + 1, isLittleEndian); + _update_tvd(ptvd); +}//set_orientation0018_9089() + +void set_bVal(struct TVolumeDiffusion* ptvd, const float b) { + ptvd->_dtiV[0] = b; + _update_tvd(ptvd); +}//set_bVal() + +void set_isAtFirstPatientPosition_tvd(struct TVolumeDiffusion* ptvd, const bool iafpp) { + ptvd->_isAtFirstPatientPosition = iafpp; + _update_tvd(ptvd); +}//set_isAtFirstPatientPosition_tvd() + +// Update the diffusion info in dd and *pdti4D for a volume once all the +// diffusion info for that volume has been read into pvd. +// +// Note that depending on the scanner software the diffusion info can arrive in +// different tags, in different orders (because of enclosing sequence tags), +// and the values in some tags may be invalid, or may be essential, depending +// on the presence of other tags. Thus it is best to gather all the diffusion +// info for a volume (frame) before taking action on it. +// +// On the other hand, dd and *pdti4D need to be updated as soon as the +// diffusion info is ready, before diffusion info for the next volume is read +// in. +void _update_tvd(struct TVolumeDiffusion* ptvd) { + // Figure out if we have both the b value and direction (if any) for this + // volume, and if isFirstPosition. + + // // GE (software version 27) is liable to NOT include kDiffusion_b_value for the + // // slice if it is 0, but should still have kDiffusionBFactor, which comes + // // after PatientPosition. + // if(isAtFirstPatientPosition && manufacturer == kMANUFACTURER_GE && dtiV[0] < 0) + // dtiV[0] = 0; // Implied 0. + + bool isReady = (ptvd->_isAtFirstPatientPosition && (ptvd->_dtiV[0] >= 0)); + if(isReady){ + for(int i = 1; i < 4; ++i){ + if(ptvd->_dtiV[i] > 1){ + isReady = false; + break; + } + } + } + if(!isReady) + return; + + // If still here, update dd and *pdti4D. + ptvd->pdd->CSA.numDti++; + if (ptvd->pdd->CSA.numDti == 2) { // First time we know that this is a 4D DTI dataset; + for(int i = 0; i < 4; ++i) // Start *pdti4D before ptvd->pdd->CSA.dtiV + ptvd->pdti4D->S[0].V[i] = ptvd->pdd->CSA.dtiV[i]; // is updated. + } + for(int i = 0; i < 4; ++i) // Update pdd + ptvd->pdd->CSA.dtiV[i] = ptvd->_dtiV[i]; + if((ptvd->pdd->CSA.numDti > 1) && (ptvd->pdd->CSA.numDti < kMaxDTI4D)){ // Update *pdti4D + //d.dti4D = (TDTI *)malloc(kMaxDTI4D * sizeof(TDTI)); + for(int i = 0; i < 4; ++i) + ptvd->pdti4D->S[ptvd->pdd->CSA.numDti - 1].V[i] = ptvd->_dtiV[i]; + } + clear_volume(ptvd); // clear the slate for the next volume. +}//_update_tvd() +//END RIR + struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, struct TDTI4D *dti4D) { struct TDICOMdata d = clear_dicom_data(); strcpy(d.protocolName, ""); //erase dummy with empty @@ -2676,6 +2879,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) + + struct TVolumeDiffusion volDiffusion = initTVolumeDiffusion(&d, dti4D); + struct stat s; if( stat(fname,&s) == 0 ) { if( !(s.st_mode & S_IFREG) ){ @@ -2800,12 +3006,19 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kInPlanePhaseEncodingDirection 0x0018+(0x1312<< 16 ) //CS #define kSAR 0x0018+(0x1316 << 16 ) //'DS' 'SAR' #define kPatientOrient 0x0018+(0x5100<< 16 ) //0018,5100. patient orientation - 'HFS' +#define kDiffusionDirectionality 0x0018+uint32_t(0x9075<< 16 ) // NONE, ISOTROPIC, or DIRECTIONAL //#define kDiffusionBFactorSiemens 0x0019+(0x100C<< 16 ) // 0019;000C;SIEMENS MR HEADER ;B_value +#define kDiffusion_bValue 0x0018+uint32_t(0x9087<< 16 ) // FD +#define kDiffusionOrientation 0x0018+uint32_t(0x9089<< 16 ) // FD, seen in enhanced + // DICOM from Philips 5.* + // and Siemens XA10. #define 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 #define kDiffusionDirectionGEY 0x0019+(0x10BC<< 16 ) //DS #define kDiffusionDirectionGEZ 0x0019+(0x10BD<< 16 ) //DS +#define kSharedFunctionalGroupsSequence 0x5200+uint32_t(0x9229<< 16 ) // SQ +#define kPerFrameFunctionalGroupsSequence 0x5200+uint32_t(0x9230<< 16 ) // SQ #define kBandwidthPerPixelPhaseEncode 0x0019+(0x1028<< 16 ) //FD #define kStudyID 0x0020+(0x0010 << 16 ) #define kSeriesNum 0x0020+(0x0011 << 16 ) @@ -2813,11 +3026,12 @@ 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 ) +#define kPatientPosition 0x0020+(0x0032 << 16 ) // Actually ImagePositionPatient! #define kOrientationACR 0x0020+(0x0035 << 16 ) #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 kLocationsInAcquisitionGE 0x0021+(0x104F<< 16 )// 'SS' 'LocationsInAcquisitionGE' #define kSamplesPerPixel 0x0028+(0x0002 << 16 ) #define kPhotometricInterpretation 0x0028+(0x0004 << 16 ) @@ -2910,6 +3124,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)) { #ifndef myLoadWholeFileToReadHeader //read one segment at a time if ((size_t)(lPos + 128) > MaxBufferSz) { //avoid overreading the file @@ -3021,6 +3236,12 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD 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", groupGroup, groupItem); + switch ( groupElement ) { case kTransferSyntax: { char transferSyntax[kDICOMStr]; @@ -3122,6 +3343,7 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD break; case kManufacturer: d.manufacturer = dcmStrManufacturer (lLength, &buffer[lPos]); + volDiffusion.manufacturer = d.manufacturer; break; case kInstitutionName: dcmStr(lLength, &buffer[lPos], d.institutionName); @@ -3205,6 +3427,9 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD case kPatientOrient : dcmStr (lLength, &buffer[lPos], d.patientOrient); break; + case kDiffusionDirectionality : // 0018, 9075 + set_directionality0018_9075(&volDiffusion, (&buffer[lPos])); + break; case kDwellTime : d.dwellTime = dcmStrInt(lLength, &buffer[lPos]); break; @@ -3217,27 +3442,27 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD break;*/ case kDiffusionDirectionGEX : - if (d.manufacturer == kMANUFACTURER_GE) d.CSA.dtiV[1] = dcmStrFloat(lLength, &buffer[lPos]); + if (d.manufacturer == kMANUFACTURER_GE) + set_diffusion_directionGE(&volDiffusion, lLength, (&buffer[lPos]), 0); break; case kDiffusionDirectionGEY : - if (d.manufacturer == kMANUFACTURER_GE) d.CSA.dtiV[2] = dcmStrFloat(lLength, &buffer[lPos]); + if (d.manufacturer == kMANUFACTURER_GE) + set_diffusion_directionGE(&volDiffusion, lLength, (&buffer[lPos]), 1); break; case kDiffusionDirectionGEZ : - if (d.manufacturer == kMANUFACTURER_GE) { - d.CSA.dtiV[3] = dcmStrFloat(lLength, &buffer[lPos]); - d.CSA.numDti = 1; - } + if (d.manufacturer == kMANUFACTURER_GE) + set_diffusion_directionGE(&volDiffusion, lLength, (&buffer[lPos]), 2); break; - case kBandwidthPerPixelPhaseEncode: - d.bandwidthPerPixelPhaseEncode = dcmFloatDouble(lLength, &buffer[lPos],d.isLittleEndian); - break; - case kStudyInstanceUID : + case kBandwidthPerPixelPhaseEncode: + d.bandwidthPerPixelPhaseEncode = dcmFloatDouble(lLength, &buffer[lPos],d.isLittleEndian); + break; + case kStudyInstanceUID : // 0020, 000D dcmStr (lLength, &buffer[lPos], d.studyInstanceUID); break; - case kSeriesInstanceUID : + case kSeriesInstanceUID : // 0020, 000E dcmStr (lLength, &buffer[lPos], d.seriesInstanceUID); break; - case kPatientPosition : + case kPatientPosition : // 0020, 0032, ImagePositionPatient if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (is2005140FSQ)) { if (!is2005140FSQwarned) printWarning("Philips R3.2.2 can report different positions for the same slice. Attempting patch.\n"); @@ -3262,6 +3487,9 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD 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]); @@ -3286,6 +3514,20 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD //int dx = 3; if (d.imageNum <= 1) d.imageNum = dcmStrInt(lLength, &buffer[lPos]); //Philips renames each image as image1 in 2001,9000 break; + case kDimensionIndexValues: // kImageNum is not enough for 4D series from Philips 5.*. + { // { necessary for initializing ndim. + 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. + printError("%d is too many dimensions. Only up to %d are supported\n", ndim, + MAX_NUMBER_OF_DIMENSIONS); + ndim = MAX_NUMBER_OF_DIMENSIONS; // Truncate + } + dcmMultiLongs(4 * ndim, &buffer[lPos], ndim, d.dimensionIndexValues, d.isLittleEndian); + } + break; case kPhotometricInterpretation: char interp[kDICOMStr]; dcmStr (lLength, &buffer[lPos], interp); @@ -3506,23 +3748,70 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD else d.sliceOrient = kSliceOrientTra; //transverse (axial) break; } - case kDiffusionBFactor: - 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 - 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]; - dti4D->S[0].V[1] = d.CSA.dtiV[1]; - dti4D->S[0].V[2] = d.CSA.dtiV[2]; - dti4D->S[0].V[3] = d.CSA.dtiV[3]; - } - d.CSA.dtiV[0] = dcmFloat(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]; - /*if ((d.CSA.numDti > 0) && (d.CSA.numDti <= kMaxDTIv)) - d.CSA.dtiV[d.CSA.numDti-1][0] = 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 + // 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]; + // dti4D->S[0].V[1] = d.CSA.dtiV[1]; + // dti4D->S[0].V[2] = d.CSA.dtiV[2]; + // dti4D->S[0].V[3] = d.CSA.dtiV[3]; + // } + // d.CSA.dtiV[0] = dcmFloat(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]; + // /*if ((d.CSA.numDti > 0) && (d.CSA.numDti <= kMaxDTIv)) + // d.CSA.dtiV[d.CSA.numDti-1][0] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);*/ + // } + // break; + case kDiffusion_bValue: // 0018, 9087 + // Note that this is ahead of kPatientPosition (0020,0032), so + // isAtFirstPatientPosition is not necessarily set yet. + // Philips uses this tag too, at least as of 5.1, but they also + // use kDiffusionBFactor (see above), and we do not want to + // double count. More importantly, with Philips this tag + // (sometimes?) gets repeated in a nested sequence with the + // 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.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]; + // dti4D->S[0].V[1] = d.CSA.dtiV[1]; + // dti4D->S[0].V[2] = d.CSA.dtiV[2]; + // dti4D->S[0].V[3] = d.CSA.dtiV[3]; + // } + //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 + // isAtFirstPatientPosition is not necessarily set yet. + // Philips uses this tag too, at least as of 5.1, but they also + // use kDiffusionDirectionRL, etc., and we do not want to double + // count. More importantly, with Philips this tag (sometimes?) + // gets repeated in a nested sequence with the value *unset*! + // 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)) + set_orientation0018_9089(&volDiffusion, lLength, &buffer[lPos], d.isLittleEndian); + break; + // case kSharedFunctionalGroupsSequence: + // if ((d.manufacturer == kMANUFACTURER_SIEMENS) && isAtFirstPatientPosition) { + // break; // For now - need to figure out how to get the nested + // // part of buffer[lPos]. + // } + // break; case kSliceNumberMrPhilips : if (d.manufacturer != kMANUFACTURER_PHILIPS) break; @@ -3561,35 +3850,35 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD break; locationsInAcquisitionPhilips = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); break; - case kDiffusionDirectionRL: - if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) { - d.CSA.dtiV[1] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); - if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D)) - dti4D->S[d.CSA.numDti-1].V[1] = d.CSA.dtiV[1]; - } - /*if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition) && (d.CSA.numDti > 0) && (d.CSA.numDti <= kMaxDTIv)) - d.CSA.dtiV[d.CSA.numDti-1][1] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);*/ - break; - case kDiffusionDirectionAP: - if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) { - d.CSA.dtiV[2] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); - if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D)) - dti4D->S[d.CSA.numDti-1].V[2] = d.CSA.dtiV[2]; - } - /*if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition) && (d.CSA.numDti > 0) && (d.CSA.numDti <= kMaxDTIv)) - d.CSA.dtiV[d.CSA.numDti-1][2] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);*/ - break; - case kDiffusionDirectionFH: - if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) { - d.CSA.dtiV[3] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); - if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D)) - dti4D->S[d.CSA.numDti-1].V[3] = d.CSA.dtiV[3]; - //printMessage("dti XYZ %g %g %g\n",d.CSA.dtiV[1],d.CSA.dtiV[2],d.CSA.dtiV[3]); - } - /*if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition) && (d.CSA.numDti > 0) && (d.CSA.numDti <= kMaxDTIv)) - 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 kDiffusionDirectionRL: + // if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) { + // d.CSA.dtiV[1] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); + // if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D)) + // dti4D->S[d.CSA.numDti-1].V[1] = d.CSA.dtiV[1]; + // } + // /*if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition) && (d.CSA.numDti > 0) && (d.CSA.numDti <= kMaxDTIv)) + // d.CSA.dtiV[d.CSA.numDti-1][1] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);*/ + // break; + // case kDiffusionDirectionAP: + // if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) { + // d.CSA.dtiV[2] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); + // if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D)) + // dti4D->S[d.CSA.numDti-1].V[2] = d.CSA.dtiV[2]; + // } + // /*if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition) && (d.CSA.numDti > 0) && (d.CSA.numDti <= kMaxDTIv)) + // d.CSA.dtiV[d.CSA.numDti-1][2] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);*/ + // break; + // case kDiffusionDirectionFH: + // if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) { + // d.CSA.dtiV[3] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian); + // if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D)) + // dti4D->S[d.CSA.numDti-1].V[3] = d.CSA.dtiV[3]; + // //printMessage("dti XYZ %g %g %g\n",d.CSA.dtiV[1],d.CSA.dtiV[2],d.CSA.dtiV[3]); + // } + // /*if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition) && (d.CSA.numDti > 0) && (d.CSA.numDti <= kMaxDTIv)) + // 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 kWaveformSq: d.imageStart = 1; //abort!!! printMessage("Skipping DICOM (audio not image) '%s'\n", fname); @@ -3619,7 +3908,8 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD if (d.manufacturer == kMANUFACTURER_GE) d.effectiveEchoSpacingGE = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); break; case kDiffusionBFactorGE : - if (d.manufacturer == kMANUFACTURER_GE) d.CSA.dtiV[0] = dcmStrInt(lLength, &buffer[lPos]); + if (d.manufacturer == kMANUFACTURER_GE) + set_bValGE(&volDiffusion, lLength, &buffer[lPos]); break; case kGeiisFlag: if ((lLength > 4) && (buffer[lPos]=='G') && (buffer[lPos+1]=='E') && (buffer[lPos+2]=='I') && (buffer[lPos+3]=='I')) { diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 67b7905f..5355c15d 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.20171204" kDCMsuf kCCsuf +#define kDCMvers "v1.0.20171208" 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 @@ -70,6 +70,10 @@ static const int kCompressC3 = 2; //obsolete JPEG lossless static const int kCompress50 = 3; //obsolete JPEG lossy static const int kCompressRLE = 4; //run length encoding +// Maximum number of dimensions for .dimensionIndexValues, i.e. possibly the +// number of axes in the output .nii. +static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; + struct TDTI { float V[4]; int sliceNumberMrPhilips; @@ -77,6 +81,7 @@ static const int kCompressRLE = 4; //run length encoding struct TDTI4D { struct TDTI S[kMaxDTI4D]; }; + #ifdef _MSC_VER //Microsoft nomenclature for packed structures is different... #pragma pack(2) typedef struct { @@ -121,6 +126,7 @@ static const int kCompressRLE = 4; //run length encoding //char mrAcquisitionType[kDICOMStr] char scanOptions[kDICOMStr], stationName[kDICOMStr], softwareVersions[kDICOMStr], deviceSerialNumber[kDICOMStr], institutionAddress[kDICOMStr], institutionName[kDICOMStr], referringPhysicianName[kDICOMStr], seriesInstanceUID[kDICOMStr], studyInstanceUID[kDICOMStr], bodyPartExamined[kDICOMStr], procedureStepDescription[kDICOMStr], imageType[kDICOMStr], institutionalDepartmentName[kDICOMStr], manufacturersModelName[kDICOMStr], patientID[kDICOMStr], patientOrient[kDICOMStr], patientName[kDICOMStr],seriesDescription[kDICOMStr], studyID[kDICOMStr], sequenceName[kDICOMStr], protocolName[kDICOMStr],sequenceVariant[kDICOMStr],scanningSequence[kDICOMStr], patientBirthDate[kDICOMStr], patientAge[kDICOMStr], studyDate[kDICOMStr],studyTime[kDICOMStr]; 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; char phaseEncodingRC, patientSex; diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 9b2a7d38..b012134d 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -70,7 +70,8 @@ #endif struct TDCMsort { - uint64_t indx, img; + uint64_t indx, img; + uint32_t dimensionIndexValues[MAX_NUMBER_OF_DIMENSIONS]; }; struct TSearchList { @@ -237,9 +238,9 @@ void geCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI *vx){ + (vx[i].V[2]*vx[i].V[2]) + (vx[i].V[3]*vx[i].V[3])); if ((vx[i].V[0] <= FLT_EPSILON)|| (vLen <= FLT_EPSILON) ) { //bvalue=0 - for (int v= 0; v < 4; v++) - vx[i].V[v] =0.0f; - continue; //do not normalize or reorient b0 vectors + for (int v= 1; v < 4; v++) + vx[i].V[v] = 0.0f; + continue; //do not normalize or reorient 0 vectors } if (!col) { //rows need to be swizzled //see Stanford dataset Ax_DWI_Tetrahedral_7 unable to resolve between possible solutions @@ -1130,15 +1131,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) { - 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 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)); @@ -2558,18 +2579,54 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis #endif free(imgM); + + // Prevent these DICOM files from being reused. + for(int i = 0; i < nConvert; ++i) + dcmList[dcmSort[i].indx].converted2NII = 1; + return returnCode;//EXIT_SUCCESS; }// saveDcm2Nii() +void fillTDCMsort(struct TDCMsort& tdcmref, const uint64_t indx, const struct TDICOMdata& dcmdata){ + // Copy the relevant parts of dcmdata to tdcmref. + tdcmref.indx = indx; + 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) - return -1; + retval = -1; else if (dcm1->img > dcm2->img) - return 1; - return 0; //tie + retval = 1; + + if(retval == 0){ + // 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]){ + retval = -1; + break; + } + else if(dcm1->dimensionIndexValues[i] > dcm2->dimensionIndexValues[i]){ + retval = 1; + break; + } + } + } + return retval; } //compareTDCMsort() int isSameFloatGE (float a, float b) { @@ -2689,10 +2746,9 @@ int singleDICOM(struct TDCMopts* opts, char *fname) { strcpy(nameList.str[nameList.numItems],filename); nameList.numItems++; struct TDCMsort dcmSort[1]; - dcmSort[0].indx = 0; - dcmSort[0].img = ((uint64_t)dcmList[0].seriesNum << 32) + dcmList[0].imageNum; 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); }// singleDICOM() @@ -2759,11 +2815,13 @@ int removeDuplicates(int nConvert, struct TDCMsort dcmSort[]){ if (nConvert < 2) return nConvert; int nDuplicates = 0; for (int i = 1; i < nConvert; i++) { - if (dcmSort[i].img == dcmSort[i-1].img) { + if (compareTDCMsort(&dcmSort[i], &dcmSort[i-1]) == 0) { nDuplicates ++; } else { dcmSort[i-nDuplicates].img = dcmSort[i].img; dcmSort[i-nDuplicates].indx = dcmSort[i].indx; + for(int j = 0; j < MAX_NUMBER_OF_DIMENSIONS; ++j) + dcmSort[i - nDuplicates].dimensionIndexValues[j] = dcmSort[i].dimensionIndexValues[j]; } } if (nDuplicates > 0) @@ -2776,18 +2834,20 @@ int removeDuplicatesVerbose(int nConvert, struct TDCMsort dcmSort[], struct TSea if (nConvert < 2) return nConvert; int nDuplicates = 0; for (int i = 1; i < nConvert; i++) { - if (dcmSort[i].img == dcmSort[i-1].img) { + if (compareTDCMsort(&dcmSort[i], &dcmSort[i-1]) == 0) { printMessage("\t%s\t=\t%s\n",nameList->str[dcmSort[i-1].indx],nameList->str[dcmSort[i].indx]); nDuplicates ++; } else { dcmSort[i-nDuplicates].img = dcmSort[i].img; dcmSort[i-nDuplicates].indx = dcmSort[i].indx; + for(int j = 0; j < MAX_NUMBER_OF_DIMENSIONS; ++j) + dcmSort[i - nDuplicates].dimensionIndexValues[j] = dcmSort[i].dimensionIndexValues[j]; } } if (nDuplicates > 0) printMessage("Some images have identical time, series, acquisition and image values. Duplicates removed.\n"); return nConvert - nDuplicates; -}// removeDuplicates() +}// removeDuplicatesVerbose() bool isExt (char *file_name, const char* ext) { char *p_extension; @@ -2914,8 +2974,8 @@ int nii_loadDir(struct TDCMopts* opts) { 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]; - dcmSort[0].indx = i; - dcmSort[0].img = ((uint64_t)dcmList[i].seriesNum << 32) + dcmList[i].imageNum; + + fillTDCMsort(dcmSort[0], i, dcmList[i]); dcmList[i].converted2NII = 1; int ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, &dti4D); if (ret == EXIT_SUCCESS) @@ -2975,9 +3035,7 @@ int nii_loadDir(struct TDCMopts* opts) { isMultiEcho = false; for (int j = i; j < (int)nDcm; j++) if (isSameSet(dcmList[i], dcmList[j], opts, &warnings, &isMultiEcho)) { - dcmSort[nConvert].indx = j; - dcmSort[nConvert].img = ((uint64_t)dcmList[j].seriesNum << 32) + dcmList[j].imageNum; - dcmList[j].converted2NII = 1; + fillTDCMsort(dcmSort[nConvert], j, dcmList[j]); nConvert++; } else { dcmList[i].isMultiEcho = isMultiEcho;