Skip to content

Commit

Permalink
Handle enahnced DICOMs lacking 0020,9157 (#372)
Browse files Browse the repository at this point in the history
  • Loading branch information
@neurolabusc
neurolabusc committed Feb 7, 2020
1 parent 0be430f commit 5390c66
Show file tree
Hide file tree
Showing 2 changed files with 62 additions and 19 deletions.
67 changes: 52 additions & 15 deletions console/nii_dicom.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3993,6 +3993,24 @@ void getFileName( char *pathParent, const char *path) {//if path is c:\d1\d2 the
getFileNameX(pathParent, path, kDICOMStr);
}

struct fidx
{
float value;
int index;
};

int cmp(const void *a, const void *b)
{
struct fidx *a1 = (struct fidx *)a;
struct fidx *a2 = (struct fidx *)b;
if ((*a1).value > (*a2).value)
return -1;
else if ((*a1).value < (*a2).value)
return 1;
else
return 0;
}

#ifdef USING_R

// True iff dcm1 sorts *before* dcm2
Expand Down Expand Up @@ -4389,11 +4407,15 @@ double TE = 0.0; //most recent echo time recorded
bool isImaginary = false;
bool isMagnitude = false;
d.seriesNum = -1;
vec3 sliceV; //cross-product of kOrientation 0020,0037
sliceV.v[0] = NAN;
float patientPositionPrivate[4] = {NAN, NAN, NAN, NAN};
float patientPosition[4] = {NAN, NAN, NAN, NAN}; //used to compute slice direction for Philips 4D
//float patientPositionPublic[4] = {NAN, NAN, NAN, NAN}; //used to compute slice direction for Philips 4D
float patientPositionEndPhilips[4] = {NAN, NAN, NAN, NAN};
float patientPositionStartPhilips[4] = {NAN, NAN, NAN, NAN};
float sliceMM[kMaxSlice2D];
int nSliceMM = 0;
//struct TDTI philDTI[kMaxDTI4D];
//for (int i = 0; i < kMaxDTI4D; i++)
// philDTI[i].V[0] = -1;
Expand Down Expand Up @@ -4656,6 +4678,11 @@ double TE = 0.0; //most recent echo time recorded
vr[1] = 'Q';
lLength = 0; //Do not skip kItemTag - required to determine nesting of Philips Enhanced
}
if ((d.manufacturer != kMANUFACTURER_PHILIPS) && (isSQ(groupElement))) { //https://github.com/rordenlab/dcm2niix/issues/144
vr[0] = 'S';
vr[1] = 'Q';
lLength = 0; //Do not skip kItemTag - required to determine nesting of Philips Enhanced
}
} //if explicit else implicit VR
if (lLength == 0xFFFFFFFF) {
lLength = 8; //SQ (Sequences) use 0xFFFFFFFF [4294967295] to denote unknown length
Expand Down Expand Up @@ -5184,8 +5211,6 @@ double TE = 0.0; //most recent echo time recorded
}
patientPositionNum++;
isAtFirstPatientPosition = true;


//char dx[kDICOMStr];
//dcmStr (lLength, &buffer[lPos], dx);
//printMessage("*%s*", dx);
Expand All @@ -5211,6 +5236,12 @@ double TE = 0.0; //most recent echo time recorded
//if (isAtFirstPatientPosition) numFirstPatientPosition++;
if (isVerbose > 0) //verbose > 1 will report full DICOM tag
printMessage(" Patient Position 0020,0032 (#,@,X,Y,Z)\t%d\t%ld\t%g\t%g\t%g\n", patientPositionNum, lPos, patientPosition[1], patientPosition[2], patientPosition[3]);
if ((isOrient) && (nSliceMM < kMaxSlice2D)) {
vec3 pos = setVec3(patientPosition[1], patientPosition[2], patientPosition[3]);
sliceMM[nSliceMM] = dotProduct(pos, sliceV);
nSliceMM++;
}

break; }
case kInPlanePhaseEncodingDirection:
d.phaseEncodingRC = toupper(buffer[lPos]); //first character is either 'R'ow or 'C'ol
Expand Down Expand Up @@ -6079,6 +6110,10 @@ double TE = 0.0; //most recent echo time recorded
}
}
dcmMultiFloat(lLength, (char*)&buffer[lPos], 6, d.orient);
vec3 readV = setVec3(d.orient[1],d.orient[2],d.orient[3]);
vec3 phaseV = setVec3(d.orient[4],d.orient[5],d.orient[6]);
sliceV = crossProduct(readV ,phaseV);
//printf("sliceV %g %g %g\n", sliceV.v[0], sliceV.v[1], sliceV.v[2]);
isOrient = true;
break; }
case kTemporalPosition :
Expand Down Expand Up @@ -6384,6 +6419,20 @@ if (d.isHasPhase)
exit (kEXIT_CORRUPT_FILE_FOUND);
#endif
}
if ((numberOfFrames > 1) && (numDimensionIndexValues == 0) && (numberOfFrames == nSliceMM)) { //issue 372
fidx* objects = (fidx*)malloc(sizeof(struct fidx) * numberOfFrames);
for (int i = 0; i < numberOfFrames; i++) {
objects[i].value = sliceMM[i];
objects[i].index = i;
}
qsort(objects, numberOfFrames, sizeof(struct fidx), cmp);
numDimensionIndexValues = numberOfFrames;
for (int i = 0; i < numberOfFrames; i++) {
// printf("%d > %g\n", objects[i].index, objects[i].value);
dcmDim[objects[i].index].dimIdx[0] = i;
}
free(objects);
} //issue 372
if (numDimensionIndexValues > 1)
strcpy(d.imageType, imageType1st); //for multi-frame datasets, return name of book, not name of last chapter
if ((numDimensionIndexValues > 1) && (numDimensionIndexValues == numberOfFrames)) {
Expand All @@ -6404,19 +6453,7 @@ if (d.isHasPhase)
if (mn[i] != mx[i])
printMessage(" Dimension %d Range: %d..%d\n", i, mn[i], mx[i]);
} //verbose > 1
/* for (int i = 0; i < numberOfFrames; i++) {
//issue 363: Philips can provide unique slope/intercept for each slice of volume. Do this BEFORE sort
dti4D->intenScale[i] = dcmDim[i].intenScale;
dti4D->intenIntercept[i] = dcmDim[i].intenIntercept;
dti4D->intenScalePhilips[i] = dcmDim[i].intenScalePhilips;
dti4D->RWVIntercept[i] = dcmDim[i].RWVIntercept;
dti4D->RWVScale[i] = dcmDim[i].RWVScale;
if (dti4D->intenIntercept[i] != dti4D->intenIntercept[0]) d.isScaleVariesEnh = true;
if (dti4D->intenScale[i] != dti4D->intenScale[0]) d.isScaleVariesEnh = true;
if (dti4D->intenScalePhilips[i] != dti4D->intenScalePhilips[0]) d.isScaleVariesEnh = true;
//printf("%g >>> %g\n", dcmDim[i].intenScale, dcmDim[i].intenIntercept);
}*/
//sort dimensions
//sort dimensions
#ifdef USING_R
std::sort(dcmDim.begin(), dcmDim.begin() + numberOfFrames, compareTDCMdim);
#else
Expand Down
14 changes: 10 additions & 4 deletions console/nii_dicom_batch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1440,7 +1440,9 @@ tse3d: T2*/
//next two conditionals updated: make GE match Siemens
if (d.phaseEncodingGE == kGE_PHASE_ENCODING_POLARITY_UNFLIPPED) phPos = 1;
if (d.phaseEncodingGE == kGE_PHASE_ENCODING_POLARITY_FLIPPED) phPos = 0;
if (((d.phaseEncodingRC == 'R') || (d.phaseEncodingRC == 'C')) && (!d.is3DAcq) && (phPos < 0)) {
bool isSkipPhaseEncodingAxis = d.is3DAcq;
if (d.echoTrainLength > 1) isSkipPhaseEncodingAxis = false; //issue 371: ignore phaseEncoding for 3D MP-RAGE/SPACE, but report for 3D EPI
if (((d.phaseEncodingRC == 'R') || (d.phaseEncodingRC == 'C')) && (!isSkipPhaseEncodingAxis) && (phPos < 0)) {
//when phase encoding axis is known but we do not know phase encoding polarity
// https://github.com/rordenlab/dcm2niix/issues/163
// This will typically correspond with InPlanePhaseEncodingDirectionDICOM
Expand All @@ -1449,7 +1451,8 @@ tse3d: T2*/
else if (d.phaseEncodingRC == 'R')
fprintf(fp, "\t\"PhaseEncodingAxis\": \"i\",\n");
}
if (((d.phaseEncodingRC == 'R') || (d.phaseEncodingRC == 'C')) && (!d.is3DAcq) && (phPos >= 0)) {
if (((d.phaseEncodingRC == 'R') || (d.phaseEncodingRC == 'C')) && (!isSkipPhaseEncodingAxis) && (phPos >= 0)) {
//printf("%ld %d %d %c %d\n", d.seriesNum, d.echoTrainLength, isSkipPhaseEncodingAxis, d.phaseEncodingRC, phPos); //test issue 371
if (d.phaseEncodingRC == 'C') //Values should be "R"ow, "C"olumn or "?"Unknown
fprintf(fp, "\t\"PhaseEncodingDirection\": \"j");
else if (d.phaseEncodingRC == 'R')
Expand Down Expand Up @@ -2801,7 +2804,10 @@ void nii_saveAttributes (struct TDICOMdata &data, struct nifti_1_header &header,

// Phase encoding polarity
// We only save these attributes if both direction and polarity are known
if (((data.phaseEncodingRC == 'R') || (data.phaseEncodingRC == 'C')) && (!data.is3DAcq) && ((data.CSA.phaseEncodingDirectionPositive == 1) || (data.CSA.phaseEncodingDirectionPositive == 0))) {
bool isSkipPhaseEncodingAxis = data.is3DAcq;
if (data.echoTrainLength > 1) isSkipPhaseEncodingAxis = false; //issue 371: ignore phaseEncoding for 3D MP-RAGE/SPACE, but report for 3D EPI

if (((data.phaseEncodingRC == 'R') || (data.phaseEncodingRC == 'C')) && (!isSkipPhaseEncodingAxis) && ((data.CSA.phaseEncodingDirectionPositive == 1) || (data.CSA.phaseEncodingDirectionPositive == 0))) {
if (data.phaseEncodingRC == 'C') {
images->addAttribute("phaseEncodingDirection", "j");
// Notice the XOR (^): the sense of phaseEncodingDirectionPositive
Expand Down Expand Up @@ -4998,7 +5004,7 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc
else
fflush(stdout); //GUI buffers printf, display all results
#endif
if ((opts.isRotate3DAcq) && (dcmList[dcmSort[0].indx].is3DAcq) && (hdr0.dim[3] > 1) && (hdr0.dim[0] < 4))
if ((opts.isRotate3DAcq) && (dcmList[dcmSort[0].indx].is3DAcq) && (dcmList[dcmSort[0].indx].echoTrainLength <= 1) && (hdr0.dim[3] > 1) && (hdr0.dim[0] < 4))
imgM = nii_setOrtho(imgM, &hdr0); //printMessage("ortho %d\n", echoInt (33));
else if (opts.isFlipY)//(FLIP_Y) //(dcmList[indx0].CSA.mosaicSlices < 2) &&
imgM = nii_flipY(imgM, &hdr0);
Expand Down

0 comments on commit 5390c66

Please sign in to comment.