Skip to content

Commit

Permalink
Respect pixel padding value during interpolation
Browse files Browse the repository at this point in the history 
* Fixes rordenlab#262
* Use a single float pixelPaddingValue in TDICOMData rather than
  pixelPaddingValue, floatPixelPaddingValue, isHasPixelPaddingValue and
  isHasFloatPixelPaddingValue
* Read kFloatPixelPaddingValue with dcmFloat() rather than dcmStrFloat()
* Use pixel padding value when available for padding added during tilt
  correction.
  • Loading branch information
@alund
alund committed Feb 1, 2019
1 parent c3f450d commit 7f378ce
Show file tree
Hide file tree
Showing 3 changed files with 46 additions and 27 deletions.
13 changes: 5 additions & 8 deletions console/nii_dicom.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -779,17 +779,14 @@ 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;
d.isScaleOrTEVaries = false;
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;
Expand Down Expand Up @@ -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]);
Expand Down
6 changes: 3 additions & 3 deletions console/nii_dicom.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
};

Expand Down
54 changes: 38 additions & 16 deletions console/nii_dicom_batch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand All @@ -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;
Expand Down Expand Up @@ -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];
Expand All @@ -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;
Expand Down

0 comments on commit 7f378ce

Please sign in to comment.