Detect variability in first frame's image position patient (#888)

console/nii_dicom.cpp

@@ -874,6 +874,7 @@ struct TDICOMdata clear_dicom_data() {
 	d.radionuclideTotalDose = 0.0;
 	d.seriesNum = 1;
 	d.acquNum = 0;
+	d.frameNum = 0; //first shall be one
 	d.imageNum = 1;
 	d.imageStart = 0;
 	d.is3DAcq = false;				  // e.g. MP-RAGE, SPACE, TFE
@@ -4373,6 +4374,7 @@ struct TDICOMdata readDICOMx(char *fname, struct TDCMprefs *prefs, struct TDTI4D
 #define kStudyDescription 0x0008 + (0x1030 << 16)  // LO
 #define kSeriesDescription 0x0008 + (0x103E << 16) // '0008' '103E' 'LO' 'SeriesDescription'
 #define kInstitutionalDepartmentName 0x0008 + (0x1040 << 16)
+#define kFrameNum 0x0008 + (0x1160 << 16) // IS
 #define kManufacturersModelName 0x0008 + (0x1090 << 16)
 #define kDerivationDescription 0x0008 + (0x2111 << 16)
 #define kReferencedImageEvidenceSQ (uint32_t)0x0008 + (0x9092 << 16)
@@ -5801,6 +5803,9 @@ struct TDICOMdata readDICOMx(char *fname, struct TDCMprefs *prefs, struct TDTI4D
 		case kInstitutionalDepartmentName:
 			dcmStr(lLength, &buffer[lPos], d.institutionalDepartmentName);
 			break;
+		case kFrameNum:
+			d.frameNum = dcmStrInt(lLength, &buffer[lPos]);
+			break;
 		case kManufacturersModelName:
 			dcmStr(lLength, &buffer[lPos], d.manufacturersModelName);
 			break;

console/nii_dicom.h

@@ -251,7 +251,7 @@ struct TDICOMdata {
 	int xyzDim[5];
 	uint32_t coilCrc, seriesUidCrc, instanceUidCrc;
 	int overlayStart[kMaxOverlay];
-	int postLabelDelay, shimGradientX, shimGradientY, shimGradientZ, phaseNumber, spoiling, mtState, partialFourierDirection, interp3D, aslFlags, durationLabelPulseGE, epiVersionGE, internalepiVersionGE, maxEchoNumGE, rawDataRunNumber, numberOfTR, numberOfImagesInGridUIH, numberOfDiffusionT2GE, numberOfDiffusionDirectionGE, tensorFileGE, diffCyclingModeGE, phaseEncodingGE, protocolBlockStartGE, protocolBlockLengthGE, modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, frequencyEncodingSteps, phaseEncodingStepsOutOfPlane, echoTrainLength, echoNum, sliceOrient, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel, locationsInAcquisition, locationsInAcquisitionConflict, compressionScheme;
+	int postLabelDelay, shimGradientX, shimGradientY, shimGradientZ, phaseNumber, spoiling, mtState, partialFourierDirection, interp3D, aslFlags, durationLabelPulseGE, epiVersionGE, internalepiVersionGE, maxEchoNumGE, rawDataRunNumber, numberOfTR, numberOfImagesInGridUIH, numberOfDiffusionT2GE, numberOfDiffusionDirectionGE, tensorFileGE, diffCyclingModeGE, phaseEncodingGE, protocolBlockStartGE, protocolBlockLengthGE, modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, frequencyEncodingSteps, phaseEncodingStepsOutOfPlane, echoTrainLength, echoNum, sliceOrient, manufacturer, converted2NII, acquNum, frameNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel, locationsInAcquisition, locationsInAcquisitionConflict, compressionScheme;
 	float compressedSensingFactor, xRayTubeCurrent, exposureTimeMs, numberOfExcitations, numberOfArms, numberOfPointsPerArm, groupDelay, decayFactor, scatterFraction, percentSampling, waterFatShift, numberOfAverages, patientSize, patientWeight, zSpacing, zThick, pixelBandwidth, SAR, phaseFieldofView, accelFactPE, accelFactOOP, flipAngle, fieldStrength, TE, TI, TR, intenScale, intenIntercept, intenScalePhilips, gantryTilt, lastScanLoc, angulation[4], velocityEncodeScaleGE;
 	float orient[7], patientPosition[4], patientPositionLast[4], xyzMM[4], stackOffcentre[4];
 	float rtia_timerGE, radionuclidePositronFraction, radionuclideTotalDose, radionuclideHalfLife, doseCalibrationFactor, injectedVolume, reconFilterSize; // PET ISOTOPE MODULE ATTRIBUTES (C.8-57)

console/nii_dicom_batch.cpp

@@ -3264,8 +3264,14 @@ bool ensureSequentialSlicePositions(int d3, int d4, struct TDCMsort dcmSort[], s
 	int minInstance = dcmList[dcmSort[0].indx].imageNum;
 	int maxInstance = minInstance;
 	int maxPhase = 1; // Philips Multi-Phase
+	int idxFrame1 = 0;
+	int nFrameIs1 = 0;
 	for (int i = 0; i < nConvert; i++) {
 		int vol = dcmList[dcmSort[i].indx].rawDataRunNumber;
+		if (dcmList[dcmSort[i].indx].frameNum == 1) {
+			nFrameIs1 ++;
+			idxFrame1 = i;
+		}
 		minVol = min(minVol, vol);
 		maxVol = max(maxVol, vol);
 		if (vol < kMaxDTI4D)
@@ -3275,6 +3281,20 @@ bool ensureSequentialSlicePositions(int d3, int d4, struct TDCMsort dcmSort[], s
 		maxInstance = max(maxInstance, instance);
 		maxPhase = max(maxPhase, dcmList[dcmSort[i].indx].phaseNumber);
 	}
+	if (nFrameIs1 > 1) {
+		//all samples of ReferencedFrameNumber (0008,1160) should have identical ImagePositionPatient
+		int lastVol = idxFrame1;
+		float maxDx = 0.0;
+		for (int i = 0; i < idxFrame1; i++) {
+			if (dcmList[dcmSort[i].indx].frameNum != 1)
+				continue;
+			float dx = fabs(intersliceDistanceSigned(dcmList[dcmSort[i].indx], dcmList[dcmSort[idxFrame1].indx]));
+			maxDx = max(dx, maxDx);
+		}
+		if (!isSameFloatGE(0.0, maxDx)) {
+			printError("%d images report ReferencedFrameNumber (0008,1160) of 1, but ImagePositionPatient varies by %gmm (issue 888).\n", nFrameIs1, maxDx);
+		}
+	}
 	bool isUseFrameReferenceTimeForVolume = false;
 	if (dcmList[dcmSort[0].indx].frameReferenceTime >= 0.0)
 		isUseFrameReferenceTimeForVolume = true;