From c3f450d845063b840f7d90e84a5ce67f421037ff Mon Sep 17 00:00:00 2001 From: Alan Lund Date: Thu, 31 Jan 2019 16:26:11 -0600 Subject: [PATCH] Respect pixel padding value during interpolation * Fixes https://github.com/rordenlab/dcm2niix/issues/262 --- console/nii_dicom.cpp | 14 ++++++++++++++ console/nii_dicom.h | 5 +++-- console/nii_dicom_batch.cpp | 24 +++++++++++++++++++++--- 3 files changed, 38 insertions(+), 5 deletions(-) diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 02b1bee5..e35e2218 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -779,6 +779,8 @@ struct TDICOMdata clear_dicom_data() { d.isRawDataStorage = false; d.isStackableSeries = false; //combine DCE series https://github.com/rordenlab/dcm2niix/issues/252 d.isXA10A = false; //https://github.com/rordenlab/dcm2niix/issues/236 + d.isHasPixelPaddingValue = false; + d.isHasFloatPixelPaddingValue = false; d.triggerDelayTime = 0.0; d.RWVScale = 0.0; d.RWVIntercept = 0.0; @@ -786,6 +788,8 @@ struct TDICOMdata clear_dicom_data() { d.bitsAllocated = 16;//bits d.bitsStored = 0; d.samplesPerPixel = 1; + d.pixelPaddingValue = -1; + d.floatPixelPaddingValue = -1.0; d.isValid = false; d.isXRay = false; d.isMultiEcho = false; @@ -3999,6 +4003,8 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16); #define kBitsAllocated 0x0028+(0x0100 << 16 ) #define kBitsStored 0x0028+(0x0101 << 16 )//US 'BitsStored' #define kIsSigned 0x0028+(0x0103 << 16 ) //PixelRepresentation +#define kPixelPaddingValue 0x0028+(0x0120 << 16 ) // https://github.com/rordenlab/dcm2niix/issues/262 +#define kFloatPixelPaddingValue 0x0028+(0x0122 << 16 ) // https://github.com/rordenlab/dcm2niix/issues/262 #define kIntercept 0x0028+(0x1052 << 16 ) #define kSlope 0x0028+(0x1053 << 16 ) //#define kSpectroscopyDataPointColumns 0x0028+(0x9002 << 16 ) //IS @@ -5028,6 +5034,14 @@ double TE = 0.0; //most recent echo time recorded case kIsSigned : //http://dicomiseasy.blogspot.com/2012/08/chapter-12-pixel-data.html d.isSigned = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); break; + case kPixelPaddingValue : + d.pixelPaddingValue = (short) dcmInt(lLength,&buffer[lPos],d.isLittleEndian); + d.isHasPixelPaddingValue = true; + break; + case kFloatPixelPaddingValue : + d.floatPixelPaddingValue = dcmStrFloat(lLength, &buffer[lPos]); + d.isHasFloatPixelPaddingValue = true; + break; case kTR : d.TR = dcmStrFloat(lLength, &buffer[lPos]); break; diff --git a/console/nii_dicom.h b/console/nii_dicom.h index 4ca932c1..d2698e8f 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -158,7 +158,7 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; uint32_t coilCrc, seriesUidCrc; int numberOfImagesInGridUIH, numberOfDiffusionDirectionGE, phaseEncodingGE, protocolBlockStartGE, protocolBlockLengthGE, modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, echoTrainLength, echoNum, sliceOrient, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,locationsInAcquisition, compressionScheme; float numberOfAverages, imagingFrequency, 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 orient[7], patientPosition[4], patientPositionLast[4], xyzMM[4], stackOffcentre[4], floatPixelPaddingValue; float rtia_timerGE, radionuclidePositronFraction, radionuclideTotalDose, radionuclideHalfLife, doseCalibrationFactor; //PET ISOTOPE MODULE ATTRIBUTES (C.8-57) float ecat_isotope_halflife, ecat_dosage; double acquisitionDuration, triggerDelayTime, RWVScale, RWVIntercept, dateTime, acquisitionTime, acquisitionDate, bandwidthPerPixelPhaseEncode; @@ -166,7 +166,8 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8; char imageComments[kDICOMStrLarge]; uint32_t dimensionIndexValues[MAX_NUMBER_OF_DIMENSIONS]; struct TCSAdata CSA; - bool isRawDataStorage, isGrayscaleSoftcopyPresentationState, isStackableSeries, isCoilVaries, isNonParallelSlices, isSegamiOasis, isXA10A, isScaleOrTEVaries, isDerived, isXRay, isMultiEcho, isValid, is3DAcq, is2DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase, isHasImaginary, isHasReal, isHasMagnitude,isHasMixed, isFloat, isResampled, isLocalizer; + short pixelPaddingValue; + bool isRawDataStorage, isGrayscaleSoftcopyPresentationState, isStackableSeries, isCoilVaries, isNonParallelSlices, isSegamiOasis, isXA10A, isScaleOrTEVaries, isDerived, isXRay, isMultiEcho, isValid, is3DAcq, is2DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase, isHasImaginary, isHasReal, isHasMagnitude,isHasMixed, isFloat, isResampled, isLocalizer, isHasPixelPaddingValue, isHasFloatPixelPaddingValue; char phaseEncodingRC, patientSex; }; diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 8bdbdd0a..a783e864 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -2986,8 +2986,17 @@ unsigned char * nii_saveNII3DtiltFloat32(char * niiFilename, struct nifti_1_head rLo = (rLo * hdrIn.dim[1]) + (s * nVox2DIn); //offset to start of row below rHi = (rHi * hdrIn.dim[1]) + (s * nVox2DIn); //offset to start of row above int rOut = (r * hdrIn.dim[1]) + (s * nVox2D); //offset to output row - for (int c = 0; c < hdrIn.dim[1]; c++) { //for each row - imOut32[rOut+c] = round( ( ((float)imIn32[rLo+c])*fracLo) + ((float)imIn32[rHi+c])*fracHi); + for (int c = 0; c < hdrIn.dim[1]; c++) { //for each column + float valLo = (float) imIn32[rLo+c]; + float valHi = (float) imIn32[rHi+c]; + if (d.isHasFloatPixelPaddingValue && (valLo == d.floatPixelPaddingValue || valHi == d.floatPixelPaddingValue)) { + // https://github.com/rordenlab/dcm2niix/issues/262 - Use nearest neighbor interpolation + // when at least one of the values is padding. + imOut32[rOut+c] = fracHi >= 0.5 ? valHi : valLo; + } + else { + imOut32[rOut+c] = round(valLo*fracLo + valHi*fracHi); + } } //for c (each column) } //rI (input row) in range } //for r (each row) @@ -3075,7 +3084,16 @@ unsigned char * nii_saveNII3Dtilt(char * niiFilename, struct nifti_1_header * hd rHi = (rHi * hdrIn.dim[1]) + (s * nVox2DIn); //offset to start of row above int rOut = (r * hdrIn.dim[1]) + (s * nVox2D); //offset to output row for (int c = 0; c < hdrIn.dim[1]; c++) { //for each row - imOut16[rOut+c] = round( ( ((float)imIn16[rLo+c])*fracLo) + ((float)imIn16[rHi+c])*fracHi); + int valLo = imIn16[rLo+c]; + int valHi = imIn16[rHi+c]; + if (d.isHasPixelPaddingValue && (valLo == d.pixelPaddingValue || valHi == d.pixelPaddingValue)) { + // https://github.com/rordenlab/dcm2niix/issues/262 - Use nearest neighbor interpolation + // when at least one of the values is padding. + imOut16[rOut+c] = fracHi >= 0.5 ? valHi : valLo; + } + else { + imOut16[rOut+c] = round((((float) valLo)*fracLo) + ((float) valHi)*fracHi); + } } //for c (each column) } //rI (input row) in range } //for r (each row)