Skip to content

Commit

Permalink
Experimental UIH support (#225)
Browse files Browse the repository at this point in the history
  • Loading branch information
@neurolabusc
neurolabusc committed Sep 17, 2018
1 parent c28e331 commit cb05c3a
Show file tree
Hide file tree
Showing 3 changed files with 83 additions and 10 deletions.
49 changes: 45 additions & 4 deletions console/nii_dicom.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -556,9 +556,9 @@ mat44 set_nii_header_x(struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1
}
return Q44;
}
if ((d.manufacturer == kMANUFACTURER_UIH) && (d.CSA.mosaicSlices > 1) ) {
printWarning("UIH Grid to matrix transform unknown\n");
} else if (d.CSA.mosaicSlices > 1) {
//next line only for Siemens mosaic: ignore for UIH grid
// https://github.com/xiangruili/dicm2nii/commit/47ad9e6d9bc8a999344cbd487d602d420fb1509f
if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.CSA.mosaicSlices > 1)) {
double nRowCol = ceil(sqrt((double) d.CSA.mosaicSlices));
double lFactorX = (d.xyzDim[1] -(d.xyzDim[1]/nRowCol) )/2.0;
double lFactorY = (d.xyzDim[2] -(d.xyzDim[2]/nRowCol) )/2.0;
Expand Down Expand Up @@ -674,7 +674,9 @@ struct TDICOMdata clear_dicom_data() {
}
for (int i=0; i < MAX_NUMBER_OF_DIMENSIONS; ++i)
d.dimensionIndexValues[i] = 0;
d.CSA.sliceTiming[0] = -1.0f; //impossible value denotes not known
//d.CSA.sliceTiming[0] = -1.0f; //impossible value denotes not known
for (int z = 0; z < kMaxEPI3D; z++)
d.CSA.sliceTiming[z] = -1.0;
d.CSA.numDti = 0;
for (int i=0; i < 5; i++)
d.xyzDim[i] = 1;
Expand All @@ -685,6 +687,7 @@ struct TDICOMdata clear_dicom_data() {
strcpy(d.imageType,"");
strcpy(d.imageComments, "");
strcpy(d.imageBaseName, "");
strcpy(d.phaseEncodingDirectionDisplayedUIH, "");
strcpy(d.studyDate, "");
strcpy(d.studyTime, "");
strcpy(d.protocolName, "");
Expand Down Expand Up @@ -3926,7 +3929,11 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16);
// #define kSeriesType 0x0054+(0x1000 << 16 )
#define kDoseCalibrationFactor 0x0054+(0x1322<< 16 )
#define kPETImageIndex 0x0054+(0x1330<< 16 )
#define kPEDirectionDisplayedUIH 0x0065+(0x1005<< 16 )//SH
#define kDiffusion_bValueUIH 0x0065+(0x1009<< 16 ) //FD
#define kNumberOfImagesInGridUIH 0x0065+(0x1050<< 16 ) //DS
#define kMRVFrameSequenceUIH 0x0065+(0x1050<< 16 ) //SQ
#define kDiffusionGradientDirectionUIH 0x0065+(0x1037<< 16 ) //FD
#define kIconImageSequence 0x0088+(0x0200 << 16 )
#define kPMSCT_RLE1 0x07a1+(0x100a << 16 ) //Elscint/Philips compression
#define kDiffusionBFactor 0x2001+(0x1003 << 16 )// FL
Expand Down Expand Up @@ -3962,6 +3969,7 @@ double TE = 0.0; //most recent echo time recorded
double contentTime = 0.0;
int philMRImageDiffBValueNumber = 0;
int sqDepth = 0;
int acquisitionTimesUIH = 0;
int sqDepth00189114 = -1;
int nDimIndxVal = -1; //tracks Philips kDimensionIndexValues
int locationsInAcquisitionGE = 0;
Expand Down Expand Up @@ -4492,6 +4500,11 @@ double TE = 0.0; //most recent echo time recorded
char acquisitionTimeTxt[kDICOMStr];
dcmStr (lLength, &buffer[lPos], acquisitionTimeTxt);
d.acquisitionTime = atof(acquisitionTimeTxt);
if (d.manufacturer != kMANUFACTURER_UIH) break;
//UIH slice timing
d.CSA.sliceTiming[acquisitionTimesUIH] = d.acquisitionTime;
//printf("%g\n", d.CSA.sliceTiming[acquisitionTimesUIH]);
acquisitionTimesUIH ++;
break;
case kContentTime :
char contentTimeTxt[kDICOMStr];
Expand Down Expand Up @@ -4738,11 +4751,38 @@ double TE = 0.0; //most recent echo time recorded
case kPETImageIndex :
PETImageIndex = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
break;
case kPEDirectionDisplayedUIH :
if (d.manufacturer != kMANUFACTURER_UIH) break;
dcmStr (lLength, &buffer[lPos], d.phaseEncodingDirectionDisplayedUIH);
break;
case kDiffusion_bValueUIH : {
if (d.manufacturer != kMANUFACTURER_UIH) break;
float v[4];
dcmMultiFloatDouble(lLength, &buffer[lPos], 1, v, d.isLittleEndian);
d.CSA.dtiV[0] = v[0];
d.CSA.numDti = 1;
//printf("%d>>>%g\n", lPos, v[0]);
break; }
case kNumberOfImagesInGridUIH :
if (d.manufacturer != kMANUFACTURER_UIH) break;
d.numberOfImagesInGridUIH = dcmStrFloat(lLength, &buffer[lPos]);
d.CSA.mosaicSlices = d.numberOfImagesInGridUIH;
break;
case kDiffusionGradientDirectionUIH : { //0065,1037
//0.03712929804225321\-0.5522387869760447\-0.8328587749392602
if (d.manufacturer != kMANUFACTURER_UIH) break;
float v[4];
dcmMultiFloatDouble(lLength, &buffer[lPos], 3, v, d.isLittleEndian);
//dcmMultiFloat(lLength, (char*)&buffer[lPos], 3, v);
//printf(">>>%g %g %g\n", v[0], v[1], v[2]);
d.CSA.dtiV[1] = v[0];
d.CSA.dtiV[2] = v[1];
d.CSA.dtiV[3] = v[2];
//vRLPhilips = v[0];
//vAPPhilips = v[1];
//vFHPhilips = v[2];
break; }

case kBitsAllocated :
d.bitsAllocated = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
break;
Expand Down Expand Up @@ -5654,6 +5694,7 @@ if (d.isHasPhase)
#ifndef myLoadWholeFileToReadHeader
fclose(file);
#endif
//printf("%g\t\t%g\t%g\t%g\n", d.CSA.dtiV[0], d.CSA.dtiV[1], d.CSA.dtiV[2], d.CSA.dtiV[3]);
//printMessage("buffer usage %d %d %d\n",d.imageStart, lPos+lFileOffset, MaxBufferSz);
return d;
} // readDICOM()
Expand Down
2 changes: 1 addition & 1 deletion console/nii_dicom.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -160,7 +160,7 @@ static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8;
float rtia_timerGE, radionuclidePositronFraction, radionuclideTotalDose, radionuclideHalfLife, doseCalibrationFactor; //PET ISOTOPE MODULE ATTRIBUTES (C.8-57)
float ecat_isotope_halflife, ecat_dosage;
double triggerDelayTime, RWVScale, RWVIntercept, dateTime, acquisitionTime, acquisitionDate, bandwidthPerPixelPhaseEncode;
char 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 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;
Expand Down
42 changes: 37 additions & 5 deletions console/nii_dicom_batch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -264,7 +264,7 @@ void siemensPhilipsCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI
//convert DTI vectors from scanner coordinates to image frame of reference
//Uses 6 orient values from ImageOrientationPatient (0020,0037)
// requires PatientPosition 0018,5100 is HFS (head first supine)
if ((d->manufacturer != kMANUFACTURER_SIEMENS) && (d->manufacturer != kMANUFACTURER_PHILIPS)) return;
if ((d->manufacturer != kMANUFACTURER_UIH) && (d->manufacturer != kMANUFACTURER_SIEMENS) && (d->manufacturer != kMANUFACTURER_PHILIPS)) return;
if (d->CSA.numDti < 1) return;
if ((toupper(d->patientOrient[0])== 'H') && (toupper(d->patientOrient[1])== 'F') && (toupper(d->patientOrient[2])== 'S'))
; //participant was head first supine
Expand Down Expand Up @@ -912,6 +912,7 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts,
int viewOrderGE = -1;
int sliceOrderGE = -1;
//n.b. https://neurostars.org/t/getting-missing-ge-information-required-by-bids-for-common-preprocessing/1357/7
json_Str(fp, "\t\"PhaseEncodingDirectionDisplayed\": \"%s\",\n", d.phaseEncodingDirectionDisplayedUIH);
if (d.phaseEncodingGE != kGE_PHASE_ENCODING_POLARITY_UNKNOWN) { //only set for GE
if (d.phaseEncodingGE == kGE_PHASE_ENCODING_POLARITY_UNFLIPPED) fprintf(fp, "\t\"PhaseEncodingPolarityGE\": \"Unflipped\",\n" );
if (d.phaseEncodingGE == kGE_PHASE_ENCODING_POLARITY_FLIPPED) fprintf(fp, "\t\"PhaseEncodingPolarityGE\": \"Flipped\",\n" );
Expand Down Expand Up @@ -1086,6 +1087,19 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts,
fprintf(fp, "-");
fprintf(fp, "\",\n");
} //only save PhaseEncodingDirection if BOTH direction and POLARITY are known
//Slice Timing UIH >>>>
if ((d.manufacturer == kMANUFACTURER_UIH) && (d.CSA.sliceTiming[0] >= 0.0)) {
fprintf(fp, "\t\"SliceTiming\": [\n");
for (int i = 0; i < h->dim[3]; i++) {
if (i != 0)
fprintf(fp, ",\n");
if (d.CSA.protocolSliceNumber1 < 0)
fprintf(fp, "\t\t%g", d.CSA.sliceTiming[(h->dim[3]-1) - i]);
else
fprintf(fp, "\t\t%g", d.CSA.sliceTiming[i]);
}
fprintf(fp, "\t],\n");
}
//Slice Timing GE >>>>
if ((d.sliceOrderGE != kGE_SLICE_ORDER_UNKNOWN) &&(d.manufacturer == kMANUFACTURER_GE) && (d.CSA.sliceTiming[0] >= 0.0)) {
fprintf(fp, "\t\"SliceTiming\": [\n");
Expand Down Expand Up @@ -1319,7 +1333,8 @@ int * nii_SaveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],str
for (int i = 0; i < numDti; i++) {
bvals[i] = vx[i].V[0];
//printMessage("---bxyz %g %g %g %g\n",vx[i].V[0],vx[i].V[1],vx[i].V[2],vx[i].V[3]);
if (isADCnotDTI(vx[i])) {
//Philips includes derived isotropic images
if ((dcmList[indx0].manufacturer == kMANUFACTURER_PHILIPS) && (isADCnotDTI(vx[i]))) {
*numADC = *numADC + 1;
bvals[i] = kADCval;
//printMessage("+++bxyz %d\n",i);
Expand Down Expand Up @@ -2703,6 +2718,18 @@ void checkSliceTiming(struct TDICOMdata * d, struct TDICOMdata * d1) {
if (d->CSA.sliceTiming[i] < minT) minT = d->CSA.sliceTiming[i];
if (d->CSA.sliceTiming[i] > maxT) maxT = d->CSA.sliceTiming[i];
}
if (d->manufacturer == kMANUFACTURER_UIH) {
maxT = 0;
for (int i = 0; i < kMaxEPI3D; i++) {
if (d->CSA.sliceTiming[i] < 0.0) break;
d->CSA.sliceTiming[i] = d->CSA.sliceTiming[i] - minT;
d->CSA.sliceTiming[i] = dicomTimeToSec(d->CSA.sliceTiming[i]);
if (d->CSA.sliceTiming[i] > maxT)
maxT = d->CSA.sliceTiming[i];
//printf("::>>>%g\n", d->CSA.sliceTiming[i]);
}
//t = (t - min(t)) * 24 * 3600 * 1000; % day to ms
}
if ((minT != maxT) && (maxT <= d->TR)) return; //looks fine
if ((minT == maxT) && (d->is3DAcq)) return; //fine: 3D EPI
if ((minT == maxT) && (d->CSA.multiBandFactor == d->CSA.mosaicSlices)) return; //fine: all slices single excitation
Expand Down Expand Up @@ -2928,6 +2955,11 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc
if (hdr0.dim[4] > 1) //for 4d datasets, last volume should be acquired before first
checkDateTimeOrder(&dcmList[dcmSort[0].indx], &dcmList[dcmSort[nConvert-1].indx]);
}
//UIH 2D slice timing
if ((dcmList[dcmSort[0].indx].manufacturer == kMANUFACTURER_UIH) && (nConvert == (hdr0.dim[3]*hdr0.dim[4])) && (hdr0.dim[3] < (kMaxEPI3D-1)) && (hdr0.dim[3] > 1) && (hdr0.dim[4] > 1)) {
for (int v = 0; v < hdr0.dim[3]; v++)
dcmList[dcmSort[0].indx].CSA.sliceTiming[v] = dcmList[dcmSort[v].indx].acquisitionTime;
}
//GE check slice timing >>>
if ((dcmList[dcmSort[0].indx].manufacturer == kMANUFACTURER_GE) && (hdr0.dim[3] < (kMaxEPI3D-1)) && (hdr0.dim[3] > 1) && (hdr0.dim[4] > 1)) {
//ignore bogus values of first volume https://neurostars.org/t/getting-missing-ge-information-required-by-bids-for-common-preprocessing/1357/6
Expand Down Expand Up @@ -3015,10 +3047,10 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc
if (sliceDir < 0) {
imgM = nii_flipZ(imgM, &hdr0);
sliceDir = abs(sliceDir); //change this, we have flipped the image so GE DTI bvecs no longer need to be flipped!
if ((dcmList[dcmSort[0].indx].manufacturer == kMANUFACTURER_GE) && (dcmList[dcmSort[0].indx].CSA.sliceTiming[0] >= 0.0) ) {
if ((dcmList[dcmSort[0].indx].manufacturer == kMANUFACTURER_UIH) && (dcmList[dcmSort[0].indx].CSA.sliceTiming[0] >= 0.0) )
dcmList[dcmSort[0].indx].CSA.protocolSliceNumber1 = -1;
if ((dcmList[dcmSort[0].indx].manufacturer == kMANUFACTURER_GE) && (dcmList[dcmSort[0].indx].CSA.sliceTiming[0] >= 0.0) )
dcmList[dcmSort[0].indx].CSA.protocolSliceNumber1 = -1;
//printWarning("GE slices flipped: check slice timing\n"); //>>>
}
}
// skip converting if user has specified one or more series, but has not specified this one
if (opts.numSeries > 0) {
Expand Down

0 comments on commit cb05c3a

Please sign in to comment.