Skip to content

Commit

Permalink
Computer TotalReadoutDuration for Siemens
Browse files Browse the repository at this point in the history
  • Loading branch information
@neurolabusc
neurolabusc committed Jun 9, 2017
1 parent 0102f8b commit 518bc77
Show file tree
Hide file tree
Showing 4 changed files with 89 additions and 29 deletions.
4 changes: 2 additions & 2 deletions COMPILE.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
```

```
Expand Down
56 changes: 44 additions & 12 deletions console/nii_dicom.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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, "");
Expand All @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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'
Expand All @@ -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 )
Expand Down Expand Up @@ -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 )
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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]);
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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];
Expand All @@ -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
Expand Down Expand Up @@ -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))
Expand Down
6 changes: 3 additions & 3 deletions console/nii_dicom.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand Down
52 changes: 40 additions & 12 deletions console/nii_dicom_batch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
/*-
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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");
Expand Down

0 comments on commit 518bc77

Please sign in to comment.