diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index e35e2218..69c4b3e0 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -779,8 +779,6 @@ 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; @@ -788,8 +786,7 @@ struct TDICOMdata clear_dicom_data() { d.bitsAllocated = 16;//bits d.bitsStored = 0; d.samplesPerPixel = 1; - d.pixelPaddingValue = -1; - d.floatPixelPaddingValue = -1.0; + d.pixelPaddingValue = NAN; d.isValid = false; d.isXRay = false; d.isMultiEcho = false; @@ -5035,12 +5032,12 @@ double TE = 0.0; //most recent echo time recorded d.isSigned = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); break; case kPixelPaddingValue : - d.pixelPaddingValue = (short) dcmInt(lLength,&buffer[lPos],d.isLittleEndian); - d.isHasPixelPaddingValue = true; + // According to the DICOM standard, this can be either unsigned (US) or signed (SS). Currently this + // is used only in nii_saveNII3Dtilt() which only allows DT_INT16, so treat it as signed. + d.pixelPaddingValue = (float) (short) dcmInt(lLength,&buffer[lPos],d.isLittleEndian); break; case kFloatPixelPaddingValue : - d.floatPixelPaddingValue = dcmStrFloat(lLength, &buffer[lPos]); - d.isHasFloatPixelPaddingValue = true; + d.pixelPaddingValue = dcmFloat(lLength, &buffer[lPos], d.isLittleEndian); break; case kTR : d.TR = dcmStrFloat(lLength, &buffer[lPos]); diff --git a/console/nii_dicom.h b/console/nii_dicom.h index d2698e8f..24c574ae 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -158,16 +158,16 @@ 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], floatPixelPaddingValue; + float orient[7], patientPosition[4], patientPositionLast[4], xyzMM[4], stackOffcentre[4]; float rtia_timerGE, radionuclidePositronFraction, radionuclideTotalDose, radionuclideHalfLife, doseCalibrationFactor; //PET ISOTOPE MODULE ATTRIBUTES (C.8-57) float ecat_isotope_halflife, ecat_dosage; + float pixelPaddingValue; // used for both FloatPixelPaddingValue (0028, 0122) and PixelPaddingValue (0028, 0120); NaN if not present. double acquisitionDuration, triggerDelayTime, RWVScale, RWVIntercept, dateTime, acquisitionTime, acquisitionDate, bandwidthPerPixelPhaseEncode; char coilElements[kDICOMStr], coilName[kDICOMStr], phaseEncodingDirectionDisplayedUIH[kDICOMStr], imageBaseName[kDICOMStr], 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; - 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; + 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; char phaseEncodingRC, patientSex; }; diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index a783e864..3551a37e 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -2960,13 +2960,24 @@ unsigned char * nii_saveNII3DtiltFloat32(char * niiFilename, struct nifti_1_head int nVox2D = hdr->dim[1]*hdr->dim[2]; unsigned char * imOut = (unsigned char *)malloc(nVox2D * hdrIn.dim[3] * 4);// *4 as 32-bits per voxel, sizeof(float) ); float * imOut32 = ( float*) imOut; - //set surrounding voxels to darkest observed value - float minVoxVal = imIn32[0]; - for (int v = 0; v < (nVox2DIn * hdrIn.dim[3]); v++) - if (imIn32[v] < minVoxVal) - minVoxVal = imIn32[v]; + + //set surrounding voxels to padding (if present) or darkest observed value + bool hasPixelPaddingValue = !isnan(d.pixelPaddingValue); + float pixelPaddingValue; + if (hasPixelPaddingValue) { + pixelPaddingValue = d.pixelPaddingValue; + } + else { + // Find darkest pixel value. Note that `hasPixelPaddingValue` remains false so that the darkest value + // will not trigger nearest neighbor interpolation below when this value is found in the image. + pixelPaddingValue = imIn32[0]; + for (int v = 0; v < (nVox2DIn * hdrIn.dim[3]); v++) + if (imIn32[v] < pixelPaddingValue) + pixelPaddingValue = imIn32[v]; + } for (int v = 0; v < (nVox2D * hdrIn.dim[3]); v++) - imOut32[v] = minVoxVal; + imOut32[v] = pixelPaddingValue; + //copy skewed voxels for (int s = 0; s < hdrIn.dim[3]; s++) { //for each slice float sliceMM = s * hdrIn.pixdim[3]; @@ -2989,7 +3000,7 @@ unsigned char * nii_saveNII3DtiltFloat32(char * niiFilename, struct nifti_1_head 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)) { + if (hasPixelPaddingValue && (valLo == pixelPaddingValue || valHi == pixelPaddingValue)) { // 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; @@ -3057,13 +3068,24 @@ unsigned char * nii_saveNII3Dtilt(char * niiFilename, struct nifti_1_header * hd int nVox2D = hdr->dim[1]*hdr->dim[2]; unsigned char * imOut = (unsigned char *)malloc(nVox2D * hdrIn.dim[3] * 2);// *2 as 16-bits per voxel, sizeof( short) ); short * imOut16 = ( short*) imOut; - //set surrounding voxels to darkest observed value - int minVoxVal = imIn16[0]; - for (int v = 0; v < (nVox2DIn * hdrIn.dim[3]); v++) - if (imIn16[v] < minVoxVal) - minVoxVal = imIn16[v]; + + //set surrounding voxels to padding (if present) or darkest observed value + bool hasPixelPaddingValue = !isnan(d.pixelPaddingValue); + short pixelPaddingValue; + if (hasPixelPaddingValue) { + pixelPaddingValue = (short) round(d.pixelPaddingValue); + } + else { + // Find darkest pixel value. Note that `hasPixelPaddingValue` remains false so that the darkest value + // will not trigger nearest neighbor interpolation below when this value is found in the image. + pixelPaddingValue = imIn16[0]; + for (int v = 0; v < (nVox2DIn * hdrIn.dim[3]); v++) + if (imIn16[v] < pixelPaddingValue) + pixelPaddingValue = imIn16[v]; + } for (int v = 0; v < (nVox2D * hdrIn.dim[3]); v++) - imOut16[v] = minVoxVal; + imOut16[v] = pixelPaddingValue; + //copy skewed voxels for (int s = 0; s < hdrIn.dim[3]; s++) { //for each slice float sliceMM = s * hdrIn.pixdim[3]; @@ -3084,9 +3106,9 @@ 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 - int valLo = imIn16[rLo+c]; - int valHi = imIn16[rHi+c]; - if (d.isHasPixelPaddingValue && (valLo == d.pixelPaddingValue || valHi == d.pixelPaddingValue)) { + short valLo = imIn16[rLo+c]; + short valHi = imIn16[rHi+c]; + if (hasPixelPaddingValue && (valLo == pixelPaddingValue || valHi == 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;