From 518bc77928fa11c5c92cfd1b7d86b6c98ee1c22b Mon Sep 17 00:00:00 2001 From: Chris Rorden Date: Fri, 9 Jun 2017 07:56:14 -0400 Subject: [PATCH] Computer TotalReadoutDuration for Siemens --- COMPILE.md | 4 +-- console/nii_dicom.cpp | 56 +++++++++++++++++++++++++++++-------- console/nii_dicom.h | 6 ++-- console/nii_dicom_batch.cpp | 52 ++++++++++++++++++++++++++-------- 4 files changed, 89 insertions(+), 29 deletions(-) diff --git a/COMPILE.md b/COMPILE.md index ed69b487..cdeb516a 100644 --- a/COMPILE.md +++ b/COMPILE.md @@ -59,13 +59,13 @@ g++ -dead_strip -O3 -I. main_console.cpp nii_dicom.cpp jpg_0XC3.cpp ujpeg.cpp ni You should then be able to run then run: ``` -g++ -O3 -dead_strip -I. main_console.cpp nii_dicom.cpp nifti1_io_core.cpp nii_ortho.cpp nii_dicom_batch.cpp jpg_0XC3.cpp ujpeg.cpp -o dcm2niix -lopenjp2 +g++ -O3 -dead_strip -I. main_console.cpp nii_dicom.cpp nifti1_io_core.cpp nii_ortho.cpp nii_dicom_batch.cpp jpg_0XC3.cpp ujpeg.cpp nii_foreign.cpp -o dcm2niix -lopenjp2 ``` But in my experience this works best if you explicitly tell the software how to find the libraries, so your compile will probably look like one of these two options: ``` -g++ -O3 -dead_strip -I. main_console.cpp nii_dicom.cpp nifti1_io_core.cpp nii_ortho.cpp nii_dicom_batch.cpp jpg_0XC3.cpp ujpeg.cpp -o dcm2niix -I/usr/local/include /usr/local/lib/libopenjp2.a +g++ -O3 -dead_strip -I. main_console.cpp nii_dicom.cpp nifti1_io_core.cpp nii_ortho.cpp nii_dicom_batch.cpp jpg_0XC3.cpp ujpeg.cpp nii_foreign.cpp -o dcm2niix -I/usr/local/include/openjpeg-2.1 /usr/local/lib/libopenjp2.a ``` ``` diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp index 375d5c41..95b27a11 100644 --- a/console/nii_dicom.cpp +++ b/console/nii_dicom.cpp @@ -637,7 +637,7 @@ struct TDICOMdata clear_dicom_data() { strcpy(d.sequenceName, ""); strcpy(d.scanningSequence, ""); strcpy(d.sequenceVariant, ""); - strcpy(d.manufacturersModelName, "N/A"); + strcpy(d.manufacturersModelName, ""); strcpy(d.procedureStepDescription, ""); strcpy(d.institutionName, ""); strcpy(d.referringPhysicianName, ""); @@ -648,6 +648,7 @@ struct TDICOMdata clear_dicom_data() { strcpy(d.studyID, ""); strcpy(d.studyInstanceUID, ""); strcpy(d.bodyPartExamined,""); + d.phaseEncodingLines = 0; d.patientPositionSequentialRepeats = 0; d.isHasPhase = false; d.isHasMagnitude = false; @@ -662,11 +663,14 @@ struct TDICOMdata clear_dicom_data() { d.TE = 0.0; d.TI = 0.0; d.flipAngle = 0.0; + d.bandwidthPerPixelPhaseEncode = 0.0; d.fieldStrength = 0.0; d.numberOfDynamicScans = 0; d.echoNum = 1; d.echoTrainLength = 0; + d.phaseEncodingSteps = 0; d.coilNum = 1; + d.accelFactPE = 0.0; d.patientPositionNumPhilips = 0; d.imageBytes = 0; d.intenScale = 1; @@ -724,7 +728,6 @@ void dcmStrDigitsOnly(char* lStr) { for (int i = 0; i < (int) len; i++) if (!isdigit(lStr[i]) ) lStr[i] = ' '; - } void dcmStr(int lLength, unsigned char lBuffer[], char* lOut) { @@ -966,11 +969,9 @@ int readCSAImageHeader(unsigned char *buff, int lLength, struct TCSAdata *CSA, i for (int lT = 1; lT <= lnTag; lT++) { memcpy(&tagCSA, &buff[lPos], sizeof(tagCSA)); //read tag lPos +=sizeof(tagCSA); - // Storage order is always little-endian, so byte-swap required values if necessary if (!littleEndianPlatform()) nifti_swap_4bytes(1, &tagCSA.nitems); - if (isVerbose > 1) //extreme verbosity: show every CSA tag printMessage("%d CSA of %s %d\n",lPos, tagCSA.name, tagCSA.nitems); if (tagCSA.nitems > 0) { @@ -1092,6 +1093,15 @@ int readCSAImageHeader(unsigned char *buff, int lLength, struct TCSAdata *CSA, i return EXIT_SUCCESS; } // readCSAImageHeader() +void dcmMultiShorts (int lByteLength, unsigned char lBuffer[], int lnShorts, uint16_t *lShorts, bool littleEndian) { +//read array of unsigned shorts US http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_6.2.html + if ((lnShorts < 1) || (lByteLength != (lnShorts * 2))) return; + memcpy(&lShorts[0], (uint16_t *)&lBuffer[0], lByteLength); + bool swap = (littleEndian != littleEndianPlatform()); + if (swap) + nifti_swap_2bytes(lnShorts, &lShorts[0]); +} //dcmMultiShorts() + void dcmMultiFloat (int lByteLength, char lBuffer[], int lnFloats, float *lFloats) { //warning: lFloats indexed from 1! will fill lFloats[1]..[nFloats] if ((lnFloats < 1) || (lByteLength < 1)) return; @@ -2559,6 +2569,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kRadionuclidePositronFraction 0x0018+(0x1076<< 16 ) #define kGantryTilt 0x0018+(0x1120 << 16 ) #define kXRayExposure 0x0018+(0x1152 << 16 ) +#define kAcquisitionMatrix 0x0018+(0x1310 << 16 ) //US #define kFlipAngle 0x0018+(0x1314 << 16 ) #define kInPlanePhaseEncodingDirection 0x0018+(0x1312<< 16 ) //CS #define kPatientOrient 0x0018+(0x5100<< 16 ) //0018,5100. patient orientation - 'HFS' @@ -2567,6 +2578,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kDiffusionDirectionGEX 0x0019+(0x10BB<< 16 ) //DS #define kDiffusionDirectionGEY 0x0019+(0x10BC<< 16 ) //DS #define kDiffusionDirectionGEZ 0x0019+(0x10BD<< 16 ) //DS +#define kBandwidthPerPixelPhaseEncode 0x0019+(0x1028<< 16 ) //FD #define kStudyID 0x0020+(0x0010 << 16 ) #define kSeriesNum 0x0020+(0x0011 << 16 ) #define kAcquNum 0x0020+(0x0012 << 16 ) @@ -2599,6 +2611,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru #define kRealWorldSlope 0x0040+uint32_t(0x9225 << 16 ) //IS dicm2nii's SlopInt_6_9 #define kDiffusionBFactorGE 0x0043+(0x1039 << 16 ) //IS dicm2nii's SlopInt_6_9 #define kCoilSiemens 0x0051+(0x100F << 16 ) +#define kImaPATModeText 0x0051+(0x1011 << 16 ) #define kLocationsInAcquisition 0x0054+(0x0081 << 16 ) #define kDoseCalibrationFactor 0x0054+(0x1322<< 16 ) #define kIconImageSequence 0x0088+(0x0200 << 16 ) @@ -2656,7 +2669,6 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru bool is2005140FSQwarned = false; //for buggy Philips bool isAtFirstPatientPosition = false; //for 3d and 4d files: flag is true for slices at same position as first slice bool isMosaic = false; - int phaseEncodingSteps = 0; int patientPositionNum = 0; int sqDepth = 0; float patientPosition[4] = {NAN, NAN, NAN, NAN}; //used to compute slice direction for Philips 4D @@ -2753,7 +2765,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru case kTransferSyntax: { char transferSyntax[kDICOMStr]; dcmStr (lLength, &buffer[lPos], transferSyntax); - //printMessage("transfer syntax '%s'\n", transferSyntax); + //printMessage("%d transfer syntax>>> '%s'\n", compressFlag, transferSyntax); if (strcmp(transferSyntax, "1.2.840.10008.1.2.1") == 0) ; //default isExplicitVR=true; //d.isLittleEndian=true else if (strcmp(transferSyntax, "1.2.840.10008.1.2.4.50") == 0) { @@ -2900,6 +2912,9 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru d.CSA.numDti = 1; } break; + case kBandwidthPerPixelPhaseEncode: + d.bandwidthPerPixelPhaseEncode = dcmFloatDouble(lLength, &buffer[lPos],d.isLittleEndian); + break; case kStudyInstanceUID : dcmStr (lLength, &buffer[lPos], d.studyInstanceUID); break; @@ -3008,11 +3023,21 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru zSpacing = dcmStrFloat(lLength, &buffer[lPos]); break; case kPhaseEncodingSteps : - phaseEncodingSteps = dcmStrInt(lLength, &buffer[lPos]); + d.phaseEncodingSteps = dcmStrInt(lLength, &buffer[lPos]); break; case kEchoTrainLength : d.echoTrainLength = dcmStrInt(lLength, &buffer[lPos]); - // printf(">>>>>>>>>>>>>>>> %d", d.echoTrainLength); + break; + case kAcquisitionMatrix : + if (lLength == 8) { + uint16_t acquisitionMatrix[4]; + dcmMultiShorts(lLength, &buffer[lPos], 4, &acquisitionMatrix[0],d.isLittleEndian); //slice position + //phaseEncodingLines stored in either image columns or rows + if (acquisitionMatrix[3] > 0) + d.phaseEncodingLines = acquisitionMatrix[3]; + if (acquisitionMatrix[2] > 0) + d.phaseEncodingLines = acquisitionMatrix[2]; + } break; case kFlipAngle : d.flipAngle = dcmStrFloat(lLength, &buffer[lPos]); @@ -3062,6 +3087,15 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru d.coilNum = 0; } break; } + case kImaPATModeText : { //e.g. Siemens iPAT x2 listed as "p2" + char accelStr[kDICOMStr]; + dcmStr (lLength, &buffer[lPos], accelStr); + char *ptr; + dcmStrDigitsOnly(accelStr); + d.accelFactPE = (float)strtof(accelStr, &ptr); + if (*ptr != '\0') + d.accelFactPE = 0.0; + break; } case kLocationsInAcquisition : d.locationsInAcquisition = dcmInt(lLength,&buffer[lPos],d.isLittleEndian); break; @@ -3292,7 +3326,6 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru } //switch/case for groupElement } //if nest - //#ifdef MY_DEBUG if (isVerbose > 1) { if ((lLength > 12) && (lLength < 128)) { //if length is greater than 8 bytes (+4 hdr) the data must be a string [or image data] char tagStr[kDICOMStr]; @@ -3311,7 +3344,6 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru } else printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tnest=%d\n", groupElement & 65535,groupElement>>16, lLength, lPos, nest); } //printMessage(" tag=%04x,%04x length=%u pos=%ld %c%c nest=%d\n", groupElement & 65535,groupElement>>16, lLength, lPos,vr[0], vr[1], nest); - //#endif lPos = lPos + (lLength); //printMessage("%d\n",d.imageStart); } //while d.imageStart == 0 @@ -3376,8 +3408,8 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru printError("Unable to decode %d-bit images with Transfer Syntax 1.2.840.10008.1.2.4.51, decompress with dcmdjpg\n", d.bitsAllocated); d.isValid = false; } - if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (isMosaic) && (d.CSA.mosaicSlices < 1) && (phaseEncodingSteps > 0) && ((d.xyzDim[1] % phaseEncodingSteps) == 0) && ((d.xyzDim[2] % phaseEncodingSteps) == 0) ) { - d.CSA.mosaicSlices = (d.xyzDim[1] / phaseEncodingSteps) * (d.xyzDim[2] / phaseEncodingSteps); + if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (isMosaic) && (d.CSA.mosaicSlices < 1) && (d.phaseEncodingSteps > 0) && ((d.xyzDim[1] % d.phaseEncodingSteps) == 0) && ((d.xyzDim[2] % d.phaseEncodingSteps) == 0) ) { + d.CSA.mosaicSlices = (d.xyzDim[1] / d.phaseEncodingSteps) * (d.xyzDim[2] / d.phaseEncodingSteps); printWarning("Mosaic inferred without CSA header (check number of slices and spatial orientation)\n"); } if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.CSA.dtiV[1] < -1.0) && (d.CSA.dtiV[2] < -1.0) && (d.CSA.dtiV[3] < -1.0)) diff --git a/console/nii_dicom.h b/console/nii_dicom.h index f6d3bb90..01f0ebe0 100644 --- a/console/nii_dicom.h +++ b/console/nii_dicom.h @@ -99,12 +99,12 @@ static const int kCompress50 = 3; //obsolete JPEG lossy struct TDICOMdata { long seriesNum; int xyzDim[5]; - int echoTrainLength, patientPositionNumPhilips, coilNum, echoNum, sliceOrient,numberOfDynamicScans, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,patientPositionSequentialRepeats,locationsInAcquisition, compressionScheme; - float flipAngle, fieldStrength, TE, TI, TR, intenScale, intenIntercept, intenScalePhilips, gantryTilt, lastScanLoc, angulation[4]; + int phaseEncodingLines, phaseEncodingSteps, echoTrainLength, patientPositionNumPhilips, coilNum, echoNum, sliceOrient,numberOfDynamicScans, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,patientPositionSequentialRepeats,locationsInAcquisition, compressionScheme; + float accelFactPE, flipAngle, fieldStrength, TE, TI, TR, intenScale, intenIntercept, intenScalePhilips, gantryTilt, lastScanLoc, angulation[4]; float orient[7], patientPosition[4], patientPositionLast[4], xyzMM[4], stackOffcentre[4]; float radionuclidePositronFraction, radionuclideTotalDose, radionuclideHalfLife, doseCalibrationFactor; //PET ISOTOPE MODULE ATTRIBUTES (C.8-57) float ecat_isotope_halflife, ecat_dosage; - double dateTime, acquisitionTime, acquisitionDate; + double dateTime, acquisitionTime, acquisitionDate, bandwidthPerPixelPhaseEncode; bool isXRay, isMultiEcho, isSlicesSpatiallySequentialPhilips, isNonImage, isValid, is3DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase,isHasMagnitude,isHasMixed, isFloat, isResampled; char phaseEncodingRC; char softwareVersions[kDICOMStr], deviceSerialNumber[kDICOMStr], institutionAddress[kDICOMStr], institutionName[kDICOMStr], referringPhysicianName[kDICOMStr], seriesInstanceUID[kDICOMStr], studyInstanceUID[kDICOMStr], bodyPartExamined[kDICOMStr], procedureStepDescription[kDICOMStr], imageType[kDICOMStr], manufacturersModelName[kDICOMStr], patientID[kDICOMStr], patientOrient[kDICOMStr], patientName[kDICOMStr],seriesDescription[kDICOMStr], studyID[kDICOMStr], sequenceName[kDICOMStr], protocolName[kDICOMStr],sequenceVariant[kDICOMStr],scanningSequence[kDICOMStr], birthDate[kDICOMStr], gender[kDICOMStr], age[kDICOMStr], studyDate[kDICOMStr],studyTime[kDICOMStr], imageComments[kDICOMStr]; diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index b94b87a8..d01cd0d1 100755 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -332,6 +332,11 @@ bool isDerived(struct TDICOMdata d) { return true; } +#ifdef myReadAsciiCsa +//read from the ASCII portion of the Siemens CSA series header +// this is not recommended: poorly documented +// it is better to stick to the binary portion of the Siemens CSA image header + #if defined(_WIN64) || defined(_WIN32) //https://opensource.apple.com/source/Libc/Libc-1044.1.2/string/FreeBSD/memmem.c /*- @@ -475,6 +480,8 @@ int siemensEchoEPIFactor(const char * filename, int csaOffset, int csaLength, i return ret; } +#endif //myReadAsciiCsa + void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, struct TDTI4D *dti4D, struct nifti_1_header *h, const char * filename) { //https://docs.google.com/document/d/1HFUkAEE-pB-angVcYe6pf_-fVf4sCpOHKesUvfb8Grc/edit# // Generate Brain Imaging Data Structure (BIDS) info @@ -518,19 +525,20 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, //if (strlen(d.patientID) > 0) // fprintf(fp, "\t\"PatientID\": \"%s\",\n", d.patientID ); } - //printMessage("-->%d %d %s\n", d.CSA.SeriesHeader_offset, d.CSA.SeriesHeader_length, filename); + #ifdef myReadAsciiCsa if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.CSA.SeriesHeader_offset > 0) && (d.CSA.SeriesHeader_length > 0) && (strlen(d.scanningSequence) > 1) && (d.scanningSequence[0] == 'E') && (d.scanningSequence[1] == 'P')) { //for EPI scans only int echoSpacing, echoTrainDuration, epiFactor; epiFactor = siemensEchoEPIFactor(filename, d.CSA.SeriesHeader_offset, d.CSA.SeriesHeader_length, &echoSpacing, &echoTrainDuration); //printMessage("ES %d ETD %d EPI %d\n", echoSpacing, echoTrainDuration, epiFactor); if (echoSpacing > 0) - fprintf(fp, "\t\"EchoSpacing\": %d,\n", echoSpacing); + fprintf(fp, "\t\"EchoSpacing\": %g,\n", echoSpacing / 1000000.0); //usec -> sec if (echoTrainDuration > 0) - fprintf(fp, "\t\"EchoTrainDuration\": %d,\n", echoTrainDuration); + fprintf(fp, "\t\"EchoTrainDuration\": %g,\n", echoTrainDuration / 1000000.0); //usec -> sec if (epiFactor > 0) fprintf(fp, "\t\"EPIFactor\": %d,\n", epiFactor); } + #endif if (d.echoTrainLength > 1) //>1 as for Siemens EPI this is 1, Siemens uses EPI factor http://mriquestions.com/echo-planar-imaging.html fprintf(fp, "\t\"EchoTrainLength\": %d,\n", d.echoTrainLength); if (d.echoNum > 1) @@ -624,18 +632,38 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, if (d.TI > 0.0) fprintf(fp, "\t\"InversionTime\": %g,\n", d.TI / 1000.0 ); if (d.ecat_isotope_halflife > 0.0) fprintf(fp, "\t\"IsotopeHalfLife\": %g,\n", d.ecat_isotope_halflife); if (d.ecat_dosage > 0.0) fprintf(fp, "\t\"Dosage\": %g,\n", d.ecat_dosage); - //fprintf(fp, "\t\"XXXX\": %g,\n", d.CSA.bandwidthPerPixelPhaseEncode ); - if ((d.CSA.bandwidthPerPixelPhaseEncode > 0.0) && (h->dim[2] > 0) && (h->dim[1] > 0)) { - float dwellTime = 0.0f; + double bandwidthPerPixelPhaseEncode = d.bandwidthPerPixelPhaseEncode; + int phaseEncodingLines = d.phaseEncodingLines; + if ((phaseEncodingLines == 0) && (h->dim[2] > 0) && (h->dim[1] > 0)) { if (h->dim[2] == h->dim[2]) //phase encoding does not matter - dwellTime = 1.0/d.CSA.bandwidthPerPixelPhaseEncode/h->dim[2]; + phaseEncodingLines = h->dim[2]; else if (d.phaseEncodingRC =='R') - dwellTime = 1.0/d.CSA.bandwidthPerPixelPhaseEncode/h->dim[2]; + phaseEncodingLines = h->dim[2]; else if (d.phaseEncodingRC =='C') - dwellTime = 1.0/d.CSA.bandwidthPerPixelPhaseEncode/h->dim[1]; - if (dwellTime != 0.0f) //as long as phase encoding = R or C or does not matter - fprintf(fp, "\t\"EffectiveEchoSpacing\": %g,\n", dwellTime ); - } + phaseEncodingLines = h->dim[1]; + } + if (bandwidthPerPixelPhaseEncode == 0.0) + bandwidthPerPixelPhaseEncode = d.CSA.bandwidthPerPixelPhaseEncode; + if (phaseEncodingLines > 0.0) fprintf(fp, "\t\"PhaseEncodingLines\": %d,\n", phaseEncodingLines ); + if (bandwidthPerPixelPhaseEncode > 0.0) + fprintf(fp, "\t\"BandwidthPerPixelPhaseEncode\": %g,\n", bandwidthPerPixelPhaseEncode ); + double effectiveEchoSpacing = 0.0; + if ((phaseEncodingLines > 0) && (bandwidthPerPixelPhaseEncode > 0.0)) + effectiveEchoSpacing = 1.0 / (bandwidthPerPixelPhaseEncode * phaseEncodingLines) ; + if (effectiveEchoSpacing > 0.0) + fprintf(fp, "\t\"EffectiveEchoSpacing\": %g,\n", effectiveEchoSpacing); + //FSL definition is start of first line until start of last line, so n-1 unless accelerated in-plane acquisition + // to check: partial Fourier, iPAT, etc. + int fencePost = 1; + if (d.accelFactPE > 1.0) + fencePost = (int)round(d.accelFactPE); //e.g. if 64 lines with iPAT=2, we want time from start of first until start of 62nd effective line + if ((d.phaseEncodingSteps > 1) && (effectiveEchoSpacing > 0.0)) + fprintf(fp, "\t\"TotalReadoutDuration\": %g,\n", effectiveEchoSpacing * ((float)d.phaseEncodingSteps - fencePost)); + if (d.accelFactPE > 1.0) { + fprintf(fp, "\t\"AccelFactPE\": %g,\n", d.accelFactPE); + if (effectiveEchoSpacing > 0.0) + fprintf(fp, "\t\"TrueEchoSpacing\": %g,\n", effectiveEchoSpacing * d.accelFactPE); + } bool first = 1; if (dti4D->S[0].sliceTiming >= 0.0) { fprintf(fp, "\t\"SliceTiming\": [\n");