Skip to content

Commit

Permalink
More XA10A Vida header info (#240)
Browse files Browse the repository at this point in the history
  • Loading branch information
@neurolabusc
neurolabusc committed Oct 20, 2018
1 parent cd08175 commit 32bdca6
Show file tree
Hide file tree
Showing 3 changed files with 116 additions and 18 deletions.
88 changes: 78 additions & 10 deletions console/nii_dicom.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -613,7 +613,7 @@ int headerDcm2NiiSForm(struct TDICOMdata d, struct TDICOMdata d2, struct nifti_
if ((d.isDerived) || ((d.bitsAllocated == 8) && (d.samplesPerPixel == 3) && (d.manufacturer == kMANUFACTURER_SIEMENS))) {
printMessage("Unable to determine spatial orientation: 0020,0037 missing (probably not a problem: derived image)\n");
} else {
printMessage("Unable to determine spatial orientation: 0020,0037 missing!\n");
printMessage("Unable to determine spatial orientation: 0020,0037 missing (Type 1 attribute: not a valid DICOM) Series %d\n", d.seriesNum);
}
}
mat44 Q44 = set_nii_header_x(d, d2, h, &sliceDir, isVerbose);
Expand Down Expand Up @@ -715,6 +715,7 @@ struct TDICOMdata clear_dicom_data() {
strcpy(d.studyInstanceUID, "");
strcpy(d.bodyPartExamined,"");
strcpy(d.coilName, "");
strcpy(d.coilElements, "");
d.phaseEncodingLines = 0;
//~ d.patientPositionSequentialRepeats = 0;
//~ d.patientPositionRepeats = 0;
Expand Down Expand Up @@ -3875,6 +3876,7 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16);
#define kRadionuclidePositronFraction 0x0018+(0x1076<< 16 )
#define kGantryTilt 0x0018+(0x1120 << 16 )
#define kXRayExposure 0x0018+(0x1152 << 16 )
#define kReceiveCoilName 0x0018+(0x1250 << 16 ) // SH
#define kAcquisitionMatrix 0x0018+(0x1310 << 16 ) //US
#define kFlipAngle 0x0018+(0x1314 << 16 )
#define kInPlanePhaseEncodingDirection 0x0018+(0x1312<< 16 ) //CS
Expand All @@ -3888,6 +3890,7 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16);
// DICOM from Philips 5.*
// and Siemens XA10.
#define kMREchoSequence 0x0018+uint32_t(0x9114<< 16 ) //SQ
#define kMRAcquisitionPhaseEncodingStepsInPlane 0x0018+uint32_t(0x9231<< 16 ) //US
#define kNumberOfImagesInMosaic 0x0019+(0x100A<< 16 ) //US NumberOfImagesInMosaic
#define kDwellTime 0x0019+(0x1018<< 16 ) //IS in NSec, see https://github.com/rordenlab/dcm2niix/issues/127
#define kLastScanLoc 0x0019+(0x101B<< 16 )
Expand All @@ -3913,18 +3916,28 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16);
#define kTriggerDelayTime 0x0020+uint32_t(0x9153<< 16 ) //FD
#define kDimensionIndexValues 0x0020+uint32_t(0x9157<< 16 ) // UL n-dimensional index of frame.
#define kInStackPositionNumber 0x0020+uint32_t(0x9057<< 16 ) // UL can help determine slices in volume
#define kLocationsInAcquisitionGE 0x0021+(0x104F<< 16 )// 'SS' 'LocationsInAcquisitionGE'
#define kRTIA_timer 0x0021+(0x105E<< 16 )// 'DS'
#define kProtocolDataBlockGE 0x0025+(0x101B<< 16 )// 'OB'
//Private Group 21 as Used by Siemens:
#define kPATModeText 0x0021+(0x1009<< 16 )//LO, see kImaPATModeText
#define kTimeAfterStart 0x0021+(0x1104<< 16 )//DS
#define kPhaseEncodingDirectionPositive 0x0021+(0x111C<< 16 )//IS
#define kRealDwellTime 0x0021+(0x1142<< 16 )//IS
#define kBandwidthPerPixelPhaseEncode21 0x0021+(0x1153<< 16 )//FD
#define kCoilElements 0x0021+(0x114F<< 16 )//LO

//g21
//Private Group 21 as used by GE:
#define kLocationsInAcquisitionGE 0x0021+(0x104F<< 16 )//SS 'LocationsInAcquisitionGE'
#define kRTIA_timer 0x0021+(0x105E<< 16 )//DS
#define kProtocolDataBlockGE 0x0025+(0x101B<< 16 )//OB
#define kSamplesPerPixel 0x0028+(0x0002 << 16 )
#define kPhotometricInterpretation 0x0028+(0x0004 << 16 )
#define kPlanarRGB 0x0028+(0x0006 << 16 )
#define kDim3 0x0028+(0x0008 << 16 ) //number of frames - for Philips this is Dim3*Dim4
#define kDim2 0x0028+(0x0010 << 16 )
#define kDim1 0x0028+(0x0011 << 16 )
#define kXYSpacing 0x0028+(0x0030 << 16 ) //'0028' '0030' 'DS' 'PixelSpacing'
#define kXYSpacing 0x0028+(0x0030 << 16 ) //DS 'PixelSpacing'
#define kBitsAllocated 0x0028+(0x0100 << 16 )
#define kBitsStored 0x0028+(0x0101 << 16 )//'0028' '0101' 'US' 'BitsStored'
#define kBitsStored 0x0028+(0x0101 << 16 )//US 'BitsStored'
#define kIsSigned 0x0028+(0x0103 << 16 )
#define kIntercept 0x0028+(0x1052 << 16 )
#define kSlope 0x0028+(0x1053 << 16 )
Expand Down Expand Up @@ -3954,6 +3967,7 @@ const uint32_t kEffectiveTE = 0x0018+ (0x9082 << 16);
#define kMRVFrameSequenceUIH 0x0065+(0x1050<< 16 ) //SQ
#define kDiffusionGradientDirectionUIH 0x0065+(0x1037<< 16 ) //FD
#define kIconImageSequence 0x0088+(0x0200 << 16 )
#define kElscintIcon 0x07a3+(0x10ce << 16 ) //see kGeiisFlag and https://github.com/rordenlab/dcm2niix/issues/239
#define kPMSCT_RLE1 0x07a1+(0x100a << 16 ) //Elscint/Philips compression
#define kDiffusionBFactor 0x2001+(0x1003 << 16 )// FL
#define kSliceNumberMrPhilips 0x2001+(0x100A << 16 ) //IS Slice_Number_MR
Expand Down Expand Up @@ -4539,9 +4553,8 @@ double TE = 0.0; //most recent echo time recorded
char acquisitionTimeTxt[kDICOMStr];
dcmStr (lLength, &buffer[lPos], acquisitionTimeTxt);
d.acquisitionTime = atof(acquisitionTimeTxt);
//printf("%s\n", acquisitionTimeTxt);
if (d.manufacturer != kMANUFACTURER_UIH) break;
//UIH slice timing
//UIH slice timing- do not use for Siemens as Siemens de-identification can corrupt this field https://github.com/rordenlab/dcm2niix/issues/236
d.CSA.sliceTiming[acquisitionTimesGE_UIH] = d.acquisitionTime;
acquisitionTimesGE_UIH ++;
break;
Expand Down Expand Up @@ -4631,6 +4644,9 @@ double TE = 0.0; //most recent echo time recorded
if (sqDepth == 0) sqDepth = 1; //should not happen, in case faulty anonymization
sqDepth00189114 = sqDepth - 1;
break;
case kMRAcquisitionPhaseEncodingStepsInPlane :
d.phaseEncodingLines = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
break;
case kNumberOfImagesInMosaic :
if (d.manufacturer == kMANUFACTURER_SIEMENS)
numberOfImagesInMosaic = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
Expand Down Expand Up @@ -4776,6 +4792,50 @@ double TE = 0.0; //most recent echo time recorded
case kImageComments:
dcmStr (lLength, &buffer[lPos], d.imageComments, true);
break;
//group 21: siemens
//g21
case kPATModeText : { //e.g. Siemens iPAT x2 listed as "p2"
char accelStr[kDICOMStr];
dcmStr (lLength, &buffer[lPos], accelStr);
char *ptr;
dcmStrDigitsOnlyKey('p', accelStr); //e.g. if "p2s4" return "2", if "s4" return ""
d.accelFactPE = (float)strtof(accelStr, &ptr);
if (*ptr != '\0')
d.accelFactPE = 0.0;
//between slice accel
dcmStr (lLength, &buffer[lPos], accelStr);
dcmStrDigitsOnlyKey('s', accelStr); //e.g. if "p2s4" return "4", if "p2" return ""
multiBandFactor = (int)strtol(accelStr, &ptr, 10);
if (*ptr != '\0')
multiBandFactor = 0.0;
//printMessage("p%gs%d\n", d.accelFactPE, multiBandFactor);
break; }
case kTimeAfterStart:
if (d.manufacturer != kMANUFACTURER_SIEMENS) break;
if (acquisitionTimesGE_UIH >= kMaxEPI3D) break;
d.CSA.sliceTiming[acquisitionTimesGE_UIH] = dcmStrFloat(lLength, &buffer[lPos]);
//printf("%d %g\n", acquisitionTimesGE_UIH, d.CSA.sliceTiming[acquisitionTimesGE_UIH]);
acquisitionTimesGE_UIH ++;
break;
case kPhaseEncodingDirectionPositive: {
if (d.manufacturer != kMANUFACTURER_SIEMENS) break;
int ph = dcmStrInt(lLength, &buffer[lPos]);
if (ph == 1) d.phaseEncodingGE = kGE_PHASE_ENCODING_POLARITY_FLIPPED;
if (ph == 0) d.phaseEncodingGE = kGE_PHASE_ENCODING_POLARITY_UNFLIPPED;
break; }
case kRealDwellTime :
if (d.manufacturer != kMANUFACTURER_SIEMENS) break;
d.dwellTime = dcmStrInt(lLength, &buffer[lPos]);
break;
case kBandwidthPerPixelPhaseEncode21:
if (d.manufacturer != kMANUFACTURER_SIEMENS) break;
d.bandwidthPerPixelPhaseEncode = dcmFloatDouble(lLength, &buffer[lPos],d.isLittleEndian);
break;
case kCoilElements:
if (d.manufacturer != kMANUFACTURER_SIEMENS) break;
dcmStr (lLength, &buffer[lPos], d.coilElements);
break;
//group 21: GE
case kLocationsInAcquisitionGE:
locationsInAcquisitionGE = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
break;
Expand Down Expand Up @@ -4926,6 +4986,11 @@ double TE = 0.0; //most recent echo time recorded
d.TE = dcmStrFloat(lLength, &buffer[lPos]);
}
break;
case kReceiveCoilName :
dcmStr (lLength, &buffer[lPos], d.coilName);
if (strlen(d.coilName) < 1) break;
d.coilCrc =(long)abs( (long)mz_crc32((unsigned char*) &d.coilName, strlen(d.coilName)));
break;
case kSlope :
d.intenScale = dcmStrFloat(lLength, &buffer[lPos]);
break;
Expand Down Expand Up @@ -4986,8 +5051,7 @@ double TE = 0.0; //most recent echo time recorded
multiBandFactor = (int)strtol(accelStr, &ptr, 10);
if (*ptr != '\0')
multiBandFactor = 0.0;
//printMessage("p%gs%d\n", d.accelFactPE, multiBandFactor);
break; }
break; }
case kLocationsInAcquisition :
d.locationsInAcquisition = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
break;
Expand Down Expand Up @@ -5061,6 +5125,10 @@ double TE = 0.0; //most recent echo time recorded
else
d.sliceOrient = kSliceOrientTra; //transverse (axial)
break; }
case kElscintIcon :
printWarning("Assuming icon SQ 07a3,10ce.\n");
isIconImageSequence = true;
break;
case kPMSCT_RLE1 :
if (d.compressionScheme != kCompressPMSCT_RLE1) break;
d.imageStart = (int)lPos + (int)lFileOffset;
Expand Down
4 changes: 2 additions & 2 deletions console/nii_dicom.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ extern "C" {
#define kCCsuf " CompilerNA" //unknown compiler!
#endif

#define kDCMvers "v1.0.20181013 " kJP2suf kLSsuf kCCsuf
#define kDCMvers "v1.0.20181020 " kJP2suf kLSsuf kCCsuf

static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic
static const int kMaxDTI4D = 18000; //maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images
Expand Down 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 acquisitionDuration, triggerDelayTime, RWVScale, RWVIntercept, dateTime, acquisitionTime, acquisitionDate, bandwidthPerPixelPhaseEncode;
char 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 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;
Expand Down
42 changes: 36 additions & 6 deletions console/nii_dicom_batch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -976,6 +976,18 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts,
if (parallelReductionFactorInPlane != (int)(d.accelFactPE))
printWarning("ParallelReductionFactorInPlane reported in DICOM [0051,1011] (%d) does not match CSA series value %d\n", (int)(d.accelFactPE), parallelReductionFactorInPlane);
}
} else { //e.g. Siemens Vida does not have CSA header, but has many attributes
json_Str(fp, "\t\"ReceiveCoilActiveElements\": \"%s\",\n", d.coilElements);
if (strcmp(d.coilElements,d.coilName) != 0)
json_Str(fp, "\t\"CoilString\": \"%s\",\n", d.coilName);
if ((d.phaseEncodingLines > d.echoTrainLength) && (d.echoTrainLength > 0)) {
float pF = (float)d.phaseEncodingLines;
if (d.accelFactPE > 1)
pF = pF / (float)d.accelFactPE; //estimate: not sure if we round up or down
pF = (float)d.echoTrainLength / pF;
if (pF < 1.0) //e.g. if difference between lines and echo length not all explained by iPAT (SENSE/GRAPPA)
fprintf(fp, "\t\"PartialFourier\": %g,\n", pf);
} //compute partial Fourier: not reported in XA10, so infer
}
#endif
if (d.CSA.multiBandFactor > 1) //AccelFactorSlice
Expand Down Expand Up @@ -1083,7 +1095,7 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts,
} //only save PhaseEncodingDirection if BOTH direction and POLARITY are known
//Slice Timing UIH or GE >>>>
//in theory, we should also report XA10 slice times here, but see series 24 of https://github.com/rordenlab/dcm2niix/issues/236
if (((d.manufacturer == kMANUFACTURER_UIH) || (d.manufacturer == kMANUFACTURER_GE)) && (d.CSA.sliceTiming[0] >= 0.0)) {
if (((d.manufacturer == kMANUFACTURER_UIH) || (d.manufacturer == kMANUFACTURER_GE) || (d.isXA10A)) && (d.CSA.sliceTiming[0] >= 0.0)) {
fprintf(fp, "\t\"SliceTiming\": [\n");
for (int i = 0; i < h->dim[3]; i++) {
if (i != 0)
Expand Down Expand Up @@ -2827,7 +2839,7 @@ void checkSliceTiming(struct TDICOMdata * d, struct TDICOMdata * d1) {
nSlices++;
if (nSlices < 1) return;
bool isSliceTimeHHMMSS = (d->manufacturer == kMANUFACTURER_UIH);
if (d->isXA10A) isSliceTimeHHMMSS = true;
//if (d->isXA10A) isSliceTimeHHMMSS = true; //for XA10 use TimeAfterStart 0x0021,0x1104 -> Siemens de-identification can corrupt acquisition ties https://github.com/rordenlab/dcm2niix/issues/236
if (isSliceTimeHHMMSS) {//convert HHMMSS to Sec
for (int i = 0; i < nSlices; i++)
d->CSA.sliceTiming[i] = dicomTimeToSec(d->CSA.sliceTiming[i]);
Expand All @@ -2841,7 +2853,7 @@ void checkSliceTiming(struct TDICOMdata * d, struct TDICOMdata * d1) {
float kNoonSec = 43200;
if ((maxT - minT) > kNoonSec) { //volume started before midnight but ended next day!
//identify and fix 'Cinderella error' where clock resets at midnight: untested
printWarning("XA10A/UIH acquisition crossed midnight: check slice timing\n");
printWarning("UIH acquisition crossed midnight: check slice timing\n");
for (int i = 0; i < nSlices; i++)
if (d->CSA.sliceTiming[i] > kNoonSec) d->CSA.sliceTiming[i] = d->CSA.sliceTiming[i] - kMidnightSec;
minT = d->CSA.sliceTiming[0];
Expand Down Expand Up @@ -3112,9 +3124,27 @@ 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]);
}
//XA10/UIH 2D slice timing
bool isXA2D = ((dcmList[dcmSort[0].indx].CSA.mosaicSlices < 2) && (dcmList[dcmSort[0].indx].manufacturer == kMANUFACTURER_SIEMENS) && (dcmList[dcmSort[0].indx].isXA10A) );
if (((isXA2D) || (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)) {
//Siemens XA10 slice timing - first https://github.com/rordenlab/dcm2niix/issues/240
if (((dcmList[dcmSort[0].indx].isXA10A)) && (nConvert == (hdr0.dim[4])) && (hdr0.dim[3] < (kMaxEPI3D-1)) && (hdr0.dim[3] > 1) && (hdr0.dim[4] > 1)) {
float mn = dcmList[dcmSort[1].indx].CSA.sliceTiming[0];
for (int v = 0; v < hdr0.dim[3]; v++) {
dcmList[dcmSort[0].indx].CSA.sliceTiming[v] = dcmList[dcmSort[1].indx].CSA.sliceTiming[v];
if (dcmList[dcmSort[0].indx].CSA.sliceTiming[v] < mn) mn = dcmList[dcmSort[0].indx].CSA.sliceTiming[v];
}
if (mn < 0.0) mn = 0.0;
int mb = 0;
for (int v = 0; v < hdr0.dim[3]; v++) {
dcmList[dcmSort[0].indx].CSA.sliceTiming[v] -= mn;
if (isSameFloatGE(dcmList[dcmSort[0].indx].CSA.sliceTiming[v], 0.0)) mb ++;
}
if ((dcmList[dcmSort[0].indx].CSA.multiBandFactor < 2) && (mb > 1))
dcmList[dcmSort[0].indx].CSA.multiBandFactor = mb;
//for (int v = 0; v < hdr0.dim[3]; v++)
// printf("XA10sliceTiming\t%d\t%g\n", v, dcmList[dcmSort[0].indx].CSA.sliceTiming[v]);

}
//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;
}
Expand Down

0 comments on commit 32bdca6

Please sign in to comment.