From b0d344e9b0fa1b1c61e34e975107dcac82e84d5a Mon Sep 17 00:00:00 2001 From: Igor Solovey Date: Wed, 15 May 2019 15:32:14 -0400 Subject: [PATCH] Bruker 4D datasets: correct slice ordering Use Dimension Index Sequence (0020,9222) to identify elements of Dimension Index Values (0020,9157) and re-order them such that In-Stack Position Number (0020,9057) is always first, and hence slices are correctly grouped together. (https://github.com/rordenlab/dcm2niix/issues/282) --- console/nii_dicom.cpp | 33 ++++++++++++++++++++++++++++++--- 1 file changed, 30 insertions(+), 3 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 13737d77..e927ef27 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -1015,6 +1015,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() + +uint32_t dcmAttributeTag (unsigned char lBuffer[], bool littleEndian) { + // read Attribute Tag (AT) value + // return in Group + (Element << 16) format + if (littleEndian) + return lBuffer[0]+(lBuffer[1]<<8)+(lBuffer[2]<<16)+(lBuffer[3]<<24); + return lBuffer[1]+(lBuffer[0]<<8)+(lBuffer[3]<<16)+(lBuffer[2]<<24); +} //dcmInt() /* //the code below trims strings after integer // does not appear required not http://en.cppreference.com/w/cpp/string/byte/atoi @@ -4094,6 +4102,7 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16); #define kTriggerDelayTime 0x0020+uint32_t(0x9153<< 16 ) //FD #define kDimensionIndexValues 0x0020+uint32_t(0x9157<< 16 ) // UL n-dimensional index of frame. #define kInStackPositionNumber 0x0020+uint32_t(0x9057<< 16 ) // UL can help determine slices in volume +#define kDimensionIndexPointer 0x0020+uint32_t(0x9165<< 16 ) //Private Group 21 as Used by Siemens: #define kSequenceVariant21 0x0021+(0x105B<< 16 )//CS #define kPATModeText 0x0021+(0x1009<< 16 )//LO, see kImaPATModeText @@ -4195,6 +4204,8 @@ double TE = 0.0; //most recent echo time recorded int locationsInAcquisitionGE = 0; int PETImageIndex = 0; int inStackPositionNumber = 0; + uint32_t dimensionIndexPointer[MAX_NUMBER_OF_DIMENSIONS]; + size_t dimensionIndexPointerCounter = 0; int maxInStackPositionNumber = 0; //int temporalPositionIdentifier = 0; int locationsInAcquisitionPhilips = 0; @@ -4348,9 +4359,24 @@ double TE = 0.0; //most recent echo time recorded printError("Too many slices to track dimensions. Only up to %d are supported\n", kMaxSlice2D); break; } + size_t dimensionIndexOrder[nDimIndxVal]; + for(size_t i = 0; i < nDimIndxVal; i++) + dimensionIndexOrder[i] = i; + + // Bruker Enhanced MR IOD: reorder dimensions to ensure InStackPositionNumber corresponds to the first one + // This will ensure correct ordering of slices in 4D datasets + if (d.manufacturer == kMANUFACTURER_BRUKER) { + for(size_t i = 1; i < dimensionIndexPointerCounter; i++){ + if (dimensionIndexPointer[i] == kInStackPositionNumber){ + //swap with first + dimensionIndexOrder[i] = 0; + dimensionIndexOrder[0] = i; + } + } + } int ndim = nDimIndxVal; for (int i = 0; i < ndim; i++) - dcmDim[numDimensionIndexValues].dimIdx[i] = d.dimensionIndexValues[i]; + dcmDim[numDimensionIndexValues].dimIdx[i] = d.dimensionIndexValues[dimensionIndexOrder[i]]; dcmDim[numDimensionIndexValues].TE = TE; dcmDim[numDimensionIndexValues].intenScale = d.intenScale; dcmDim[numDimensionIndexValues].intenIntercept = d.intenIntercept; @@ -5026,6 +5052,9 @@ double TE = 0.0; //most recent echo time recorded //printf("<%d>\n",inStackPositionNumber); if (inStackPositionNumber > maxInStackPositionNumber) maxInStackPositionNumber = inStackPositionNumber; break; + case kDimensionIndexPointer: + dimensionIndexPointer[dimensionIndexPointerCounter++] = dcmAttributeTag(&buffer[lPos],d.isLittleEndian); + break; case kFrameContentSequence : //if (!(d.manufacturer == kMANUFACTURER_BRUKER)) break; //see https://github.com/rordenlab/dcm2niix/issues/241 if (sqDepth == 0) sqDepth = 1; //should not happen, in case faulty anonymization @@ -6137,14 +6166,12 @@ if (d.isHasPhase) if (mn[i] != mx[i]) printMessage(" Dimension %d Range: %d..%d\n", i, mn[i], mx[i]); } //verbose > 1 - if (d.manufacturer != kMANUFACTURER_BRUKER) { //only single sample Bruker - perhaps use 0020,9057 to identify if space or time is 3rd dimension //sort dimensions #ifdef USING_R std::sort(dcmDim.begin(), dcmDim.begin() + numberOfFrames, compareTDCMdim); #else qsort(dcmDim, numberOfFrames, sizeof(struct TDCMdim), compareTDCMdim); #endif - } //for (int i = 0; i < numberOfFrames; i++) // printf("diskPos= %d dimIdx= %d %d %d %d TE= %g\n", i, dcmDim[i].diskPos, dcmDim[i].dimIdx[1], dcmDim[i].dimIdx[2], dcmDim[i].dimIdx[3], dti4D->TE[i]); for (int i = 0; i < numberOfFrames; i++)