Skip to content

Commit

Permalink
Improve PET and PET BIDS (#886)
Browse files Browse the repository at this point in the history
  • Loading branch information
@neurolabusc
neurolabusc committed Oct 31, 2024
1 parent ead9ed1 commit b3605c4
Show file tree
Hide file tree
Showing 4 changed files with 79 additions and 9 deletions.
9 changes: 9 additions & 0 deletions BIDS/README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,6 +248,15 @@ Data unique to [GE](https://github.com/rordenlab/dcm2niix/tree/master/GE). Deter
| TablePosition | mm | The 3rd value of DICOM tag 0043,10B2 - the value of DICOM tag 0019,107F | B |
| DeepLearningFactor | | DICOM tag 0043,10CA | D |

### Manufacturer GE (Positron Emission Tomography)

[GE Private Tags](https://www.gehealthcare.com/-/jssmedia/b86f641ae0bd4a7b919c79215e5c01e7)

| Field | Unit | Comments | Defined By |
|------------------------------------|------|-----------------------|------------|
| InjectedVolume | | DICOM tag 0009,103A | B |
| ReconFilterSize | | DICOM tag 0009,108F | B |

### Manufacturer Philips

Data unique to Philips, including [custom intensity scaling](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3998685/)
Expand Down
15 changes: 15 additions & 0 deletions console/nii_dicom.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -865,6 +865,8 @@ struct TDICOMdata clear_dicom_data() {
d.radionuclidePositronFraction = 0.0;
d.radionuclideHalfLife = 0.0;
d.doseCalibrationFactor = 0.0;
d.injectedVolume = 0.0;
d.reconFilterSize = NAN;
d.ecat_isotope_halflife = 0.0;
d.frameDuration = -1.0;
d.frameReferenceTime = -1.0;
Expand Down Expand Up @@ -4376,6 +4378,8 @@ struct TDICOMdata readDICOMx(char *fname, struct TDCMprefs *prefs, struct TDTI4D
#define kReferencedImageEvidenceSQ (uint32_t)0x0008 + (0x9092 << 16)
#define kComplexImageComponent (uint32_t)0x0008 + (0x9208 << 16) //'0008' '9208' 'CS' 'ComplexImageComponent'
#define kAcquisitionContrast (uint32_t)0x0008 + (0x9209 << 16) //'0008' '9209' 'CS' 'AcquisitionContrast'
#define kInjectedVolumeGE 0x0009 + (0x103A << 16) // FL
#define kReconFilterSizeGE 0x0009 + (0x108F << 16) // FL bp_filter_cutoff
#define kIconSQ 0x0009 + (0x1110 << 16)
#define kPatientName 0x0010 + (0x0010 << 16)
#define kPatientID 0x0010 + (0x0020 << 16)
Expand Down Expand Up @@ -4465,6 +4469,7 @@ struct TDICOMdata readDICOMx(char *fname, struct TDCMprefs *prefs, struct TDTI4D
// #define kFrameAcquisitionDuration 0x0018+uint32_t(0x9220 << 16 ) //FD
#define kArterialSpinLabelingContrast 0x0018 + uint32_t(0x9250 << 16) // CS
#define kASLPulseTrainDuration 0x0018 + uint32_t(0x9258 << 16) // UL
//TODO ASL LabelingOrientation 0018,9255, VascularCrushing 0x0018,9259 CS, VascularCrushingVENC 0018,925A
#define kDiffusionBValueXX 0x0018 + uint32_t(0x9602 << 16) // FD
// #define kDiffusionBValueXY 0x0018 + uint32_t(0x9603 << 16) //FD
// #define kDiffusionBValueXZ 0x0018 + uint32_t(0x9604 << 16) //FD
Expand Down Expand Up @@ -5664,6 +5669,16 @@ struct TDICOMdata readDICOMx(char *fname, struct TDCMprefs *prefs, struct TDTI4D
if (strlen(d.studyTime) < 2)
dcmStr(lLength, &buffer[lPos], d.studyTime);
break;
case kInjectedVolumeGE:
if (d.manufacturer != kMANUFACTURER_GE)
break;
d.injectedVolume = dcmFloat(lLength, &buffer[lPos], d.isLittleEndian);
break;
case kReconFilterSizeGE:
if (d.manufacturer != kMANUFACTURER_GE)
break;
d.reconFilterSize = dcmFloat(lLength, &buffer[lPos], d.isLittleEndian);
break;
case kIconSQ: { // issue760 GEIIS icon strikes again
if ((vr[0] != 'S') || (vr[1] != 'Q') || (lLength != 8))
break; // only risk kludge for explicit VR with length
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 @@ -254,7 +254,7 @@ struct TDICOMdata {
int postLabelDelay, shimGradientX, shimGradientY, shimGradientZ, phaseNumber, spoiling, mtState, partialFourierDirection, interp3D, aslFlags, durationLabelPulseGE, epiVersionGE, internalepiVersionGE, maxEchoNumGE, rawDataRunNumber, numberOfTR, numberOfImagesInGridUIH, numberOfDiffusionT2GE, numberOfDiffusionDirectionGE, tensorFileGE, diffCyclingModeGE, phaseEncodingGE, protocolBlockStartGE, protocolBlockLengthGE, modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, frequencyEncodingSteps, phaseEncodingStepsOutOfPlane, echoTrainLength, echoNum, sliceOrient, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel, locationsInAcquisition, locationsInAcquisitionConflict, compressionScheme;
float compressedSensingFactor, xRayTubeCurrent, exposureTimeMs, numberOfExcitations, numberOfArms, numberOfPointsPerArm, groupDelay, decayFactor, scatterFraction, percentSampling, waterFatShift, numberOfAverages, patientSize, patientWeight, zSpacing, zThick, pixelBandwidth, SAR, phaseFieldofView, accelFactPE, accelFactOOP, flipAngle, fieldStrength, TE, TI, TR, intenScale, intenIntercept, intenScalePhilips, gantryTilt, lastScanLoc, angulation[4], velocityEncodeScaleGE;
float orient[7], patientPosition[4], patientPositionLast[4], xyzMM[4], stackOffcentre[4];
float rtia_timerGE, radionuclidePositronFraction, radionuclideTotalDose, radionuclideHalfLife, doseCalibrationFactor; // PET ISOTOPE MODULE ATTRIBUTES (C.8-57)
float rtia_timerGE, radionuclidePositronFraction, radionuclideTotalDose, radionuclideHalfLife, doseCalibrationFactor, injectedVolume, reconFilterSize; // PET ISOTOPE MODULE ATTRIBUTES (C.8-57)
float frameReferenceTime, frameDuration, ecat_isotope_halflife, ecat_dosage;
float pixelPaddingValue; // used for both FloatPixelPaddingValue (0028, 0122) and PixelPaddingValue (0028, 0120); NaN if not present.
double imagingFrequency, acquisitionDuration, triggerDelayTime, RWVScale, RWVIntercept, dateTime, acquisitionTime, acquisitionDate, bandwidthPerPixelPhaseEncode;
Expand Down
62 changes: 54 additions & 8 deletions console/nii_dicom_batch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -793,6 +793,7 @@ void siemensCsaAscii(const char *filename, TCsaAscii *csaAscii, int csaOffset, i
csaAscii->lInvContrasts = readKey(keyStrNumInv, keyPos, csaLengthTrim);
char keyStrNumEcho[] = "lContrasts";
csaAscii->lContrasts = readKey(keyStrNumEcho, keyPos, csaLengthTrim);
// TODO: read sAsl.ulSuppressionMode for required BackgroundSuppression
char keyStrDS[] = "sDiffusion.dsScheme";
csaAscii->difBipolar = readKey(keyStrDS, keyPos, csaLengthTrim);
if (csaAscii->difBipolar == 0) {
Expand Down Expand Up @@ -1207,7 +1208,7 @@ void print_FloatNotNan(const char *sLabel, int iVal, float sVal) {
printMessage(sLabel, iVal, sVal);
} // json_Float

void json_Float(FILE *fp, const char *sLabel, float sVal) {
void json_Float(FILE *fp, const char *sLabel, double sVal) {
if (!isfinite(sVal)) { // isfinite() defined in C99
// https://github.com/bids-standard/bids-2-devel/issues/12
printWarning(sLabel, sVal);
Expand Down Expand Up @@ -1565,14 +1566,18 @@ tse3d: T2*/
}
// https://bids-specification--622.org.readthedocs.build/en/622/04-modality-specific-files/01-magnetic-resonance-imaging-data.html#case-3-direct-field-mapping
if ((d.isRealIsPhaseMapHz) && ((d.manufacturer == kMANUFACTURER_GE) || (d.isHasReal)))
fprintf(fp, "\t\"Units\": \"Hz\",\n"); //
fprintf(fp, "\t\"Units\": \"Hz\",\n");
// PET ISOTOPE MODULE ATTRIBUTES
json_Str(fp, "\t\"TracerRadionuclide\": \"%s\",\n", d.tracerRadionuclide);
json_Str(fp, "\t\"Radiopharmaceutical\": \"%s\",\n", d.radiopharmaceutical);
json_Float(fp, "\t\"RadionuclidePositronFraction\": %g,\n", d.radionuclidePositronFraction);
json_Float(fp, "\t\"InjectedRadioactivity\": %g,\n", d.radionuclideTotalDose); // renamed https://bids-specification.readthedocs.io/en/stable/glossary.html#objects.metadata.InjectedRadioactivity
json_Float(fp, "\t\"InjectedRadioactivity\": %.8g,\n", d.radionuclideTotalDose); // renamed https://bids-specification.readthedocs.io/en/stable/glossary.html#objects.metadata.InjectedRadioactivity
if (d.radionuclideTotalDose > 0.0)
fprintf(fp, "\t\"InjectedRadioactivityUnits\": \"MBq\",\n");
json_Float(fp, "\t\"InjectedVolume\": %g,\n", d.injectedVolume);

json_Float(fp, "\t\"RadionuclideHalfLife\": %g,\n", d.radionuclideHalfLife);
json_Float(fp, "\t\"DoseCalibrationFactor\": %g,\n", d.doseCalibrationFactor);
json_Float(fp, "\t\"DoseCalibrationFactor\": %.8g,\n", d.doseCalibrationFactor);
json_Float(fp, "\t\"IsotopeHalfLife\": %g,\n", d.ecat_isotope_halflife);
json_Float(fp, "\t\"Dosage\": %g,\n", d.ecat_dosage);
json_Str(fp, "\t\"ConvolutionKernel\": \"%s\",\n", d.convolutionKernel);
Expand Down Expand Up @@ -1671,7 +1676,6 @@ tse3d: T2*/
// printf("::::%s ->'%s' : s%d i%d\n", d.reconstructionMethod, reconMethodName, subsets, iterations);
// END issue 802
json_Float(fp, "\t\"ScatterFraction\": %g,\n", d.scatterFraction);

if (dti4D->decayFactor[0] >= 0.0) { // see BEP009 PET https://docs.google.com/document/d/1mqMLnxVdLwZjDd4ZiWFqjEAmOmfcModA_R535v3eQs0
fprintf(fp, "\t\"DecayFactor\": [\n");
for (int i = 0; i < h->dim[4]; i++) {
Expand All @@ -1682,6 +1686,8 @@ tse3d: T2*/
fprintf(fp, "\t\t%g", dti4D->decayFactor[i]);
}
fprintf(fp, "\t],\n");
} else if (d.decayFactor > 0) { //single volume
fprintf(fp, "\t\"DecayFactor\": [\n\t\t%g\t],\n", d.decayFactor);
}
if ((h->dim[4] > 1) && (dti4D->volumeOnsetTime[0] >= 0.0)) { // see BEP009 PET https://docs.google.com/document/d/1mqMLnxVdLwZjDd4ZiWFqjEAmOmfcModA_R535v3eQs0
fprintf(fp, "\t\"FrameTimesStart\": [\n");
Expand All @@ -1693,8 +1699,11 @@ tse3d: T2*/
fprintf(fp, "\t\t%g", dti4D->volumeOnsetTime[i]);
}
fprintf(fp, "\t],\n");
} else if ((d.decayFactor > 0) && (h->dim[4] == 1)) { //single volume
fprintf(fp, "\t\"FrameTimesStart\": [\n\t\t0\t],\n");
}
if ((h->dim[4] > 1) && (dti4D->frameDuration[0] >= 0.0)) { // see BEP009 PET https://docs.google.com/document/d/1mqMLnxVdLwZjDd4ZiWFqjEAmOmfcModA_R535v3eQs0

if ((h->dim[4] > 0) && (dti4D->frameDuration[0] >= 0.0)) { // see BEP009 PET https://docs.google.com/document/d/1mqMLnxVdLwZjDd4ZiWFqjEAmOmfcModA_R535v3eQs0
fprintf(fp, "\t\"FrameDuration\": [\n");
for (int i = 0; i < h->dim[4]; i++) {
if (i != 0)
Expand Down Expand Up @@ -1725,6 +1734,19 @@ tse3d: T2*/
fprintf(fp, "\t],\n");
}
}
json_FloatNotNan(fp, "\t\"ReconFilterSize\": %g,\n", d.reconFilterSize);
if ((d.modality == kMODALITY_PT) && (d.acquisitionTime > 0.0)) {
// to do: should we use post_inj_datetime 0009,103d for GE?
// are we really sure that acquisitionTime is timezero?
// Convert float to integer for easier manipulation
int time = (int)d.acquisitionTime;
// Extract hours, minutes, and seconds
int hours = time / 10000;
int minutes = (time / 100) % 100;
int seconds = time % 100;
// Format into HH:MM:SS and store in timeString
fprintf(fp, "\t\"TimeZero\": \"%02d:%02d:%02d\",\n", hours, minutes, seconds);
}
// CT parameters
json_Float(fp, "\t\"ExposureTime\": %g,\n", d.exposureTimeMs / 1000.0);
json_Float(fp, "\t\"XRayTubeCurrent\": %g,\n", d.xRayTubeCurrent);
Expand Down Expand Up @@ -1860,6 +1882,7 @@ tse3d: T2*/
*/
bool isPCASL = false;
bool isPASL = false;
int nPLD = 0;
// ASL specific tags - 2D pCASL Danny J.J. Wang http://www.loft-lab.org
// ASL specific tags - 2D PASL Siemens Product XA30, n.b pasl/casl/pcasl set using 0018,9250
if (strstr(pulseSequenceDetails, "ep2d_asl")) {
Expand Down Expand Up @@ -1896,19 +1919,22 @@ tse3d: T2*/
// ASL specific tags - 2D PASL Siemens Product
if (strstr(pulseSequenceDetails, "ep2d_pasl")) {
isPASL = true;
fprintf(fp, "\t\"PASLType\": \"PICORE\",\n");
json_FloatNotNan(fp, "\t\"BolusDuration\": %g,\n", csaAscii.alTI[0] * (1.0 / 1000000.0)); // usec->sec
json_FloatNotNan(fp, "\t\"InversionTime\": %g,\n", csaAscii.alTI[2] * (1.0 / 1000000.0)); // usec -> sec
}
// ASL specific tags - 3D PASL Siemens Product http://adni.loni.usc.edu/wp-content/uploads/2010/05/ADNI3_Basic_Siemens_Skyra_E11.pdf
if (strstr(pulseSequenceDetails, "tgse_pasl")) {
isPASL = true;
fprintf(fp, "\t\"PASLType\": \"FAIR QII\",\n");
json_FloatNotNan(fp, "\t\"BolusDuration\": %g,\n", csaAscii.alTI[0] * (1.0 / 1000000.0)); // usec->sec
json_FloatNotNan(fp, "\t\"InversionTime\": %g,\n", csaAscii.alTI[2] * (1.0 / 1000000.0)); // usec -> sec
// json_FloatNotNan(fp, "\t\"SaturationStopTime\": %g,\n", csaAscii.alTI[2] * (1.0/1000.0));
}
// PASL http://www.pubmed.com/11746944 http://www.pubmed.com/21606572
if (strstr(pulseSequenceDetails, "ep2d_fairest")) {
isPASL = true;
fprintf(fp, "\t\"PASLType\": \"FAIR\",\n");
json_FloatNotNan(fp, "\t\"PostInversionDelay\": %g,\n", csaAscii.adFree[2] * (1.0 / 1000.0)); // usec->sec
json_FloatNotNan(fp, "\t\"PostLabelingDelay\": %g,\n", csaAscii.adFree[4] * (1.0 / 1000.0)); // usec -> sec
}
Expand All @@ -1928,7 +1954,6 @@ tse3d: T2*/
json_Float(fp, "\t\"TagDuration\": %g,\n", csaAscii.alFree[9] / 1000.0); // ms -> sec
json_Float(fp, "\t\"MaximumT1Opt\": %g,\n", csaAscii.alFree[10] / 1000.0); // ms -> sec
// report post label delay
int nPLD = 0;
bool isValid = true; // detect gaps in PLD array: If user sets PLD1=250, PLD2=0 PLD3=375 only PLD1 was acquired
for (int k = 11; k < 31; k++) {
if ((isnan(csaAscii.alFree[k])) || (csaAscii.alFree[k] <= 0.0))
Expand Down Expand Up @@ -1971,7 +1996,6 @@ tse3d: T2*/
json_Float(fp, "\t\"Tag1\": %g,\n", csaAscii.alFree[11] / 1000.0); // DelayTimeInTR usec -> sec
json_Float(fp, "\t\"Tag2\": %g,\n", csaAscii.alFree[12] / 1000.0); // DelayTimeInTR usec -> sec
json_Float(fp, "\t\"Tag3\": %g,\n", csaAscii.alFree[13] / 1000.0); // DelayTimeInTR usec -> sec
int nPLD = 0;
bool isValid = true; // detect gaps in PLD array: If user sets PLD1=250, PLD2=0 PLD3=375 only PLD1 was acquired
for (int k = 30; k < 38; k++) {
if ((isnan(csaAscii.alFree[k])) || (csaAscii.alFree[k] <= 0.0))
Expand Down Expand Up @@ -2012,6 +2036,28 @@ tse3d: T2*/
// https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html#common-metadata-fields-applicable-to-both-pcasl-and-pasl
if (((isPASL) || (isPCASL)) && (csaAscii.interp <= 0))
fprintf(fp, "\t\"AcquisitionVoxelSize\": [\n\t\t%g,\n\t\t%g,\n\t\t%g\t],\n", d.xyzMM[1], d.xyzMM[2], d.zThick);
int maxEchoNum = csaAscii.lContrasts; // this stores number of echoes, but maybe other contrasts (PLD)
if (maxEchoNum < 1)
maxEchoNum = 1;
if (nPLD < 1)
nPLD = 1;
int nContrasts = nPLD * maxEchoNum;
if (isPASL || isPCASL) {
int totalAcquiredPairs = 0;
if ((h->dim[4] % (nContrasts * 2)) == 0) {
// n.b. M0Type requires information about other series
// either "Separate" or "Absent"
// fprintf(fp, "\t\"M0Type\": \"Separate\",\n");
int totalAcquiredPairs = h->dim[4] / nContrasts;
} else if (((h->dim[4] - 1) % (nContrasts * 2)) == 0) {
fprintf(fp, "\t\"M0Type\": \"Included\",\n");
int totalAcquiredPairs = (h->dim[4] - 1) / nContrasts;
} else {
printWarning("Unable to determine M0Type\n");
}
if (totalAcquiredPairs > 0)
fprintf(fp, "\t\"TotalAcquiredPairs\": %d,\n", totalAcquiredPairs);
}
// lund free waveform sequence, see https://github.com/filip-szczepankiewicz/fwf_sequence_tools
if ((strstr(pulseSequenceDetails, "ep2d_diff_fwf") != 0) || (strstr(pulseSequenceDetails, "ep2d_diff_sms_fwf_simple") != 0)) {
for (int i = 0; i < kMaxWipFree; i++) {
Expand Down

0 comments on commit b3605c4

Please sign in to comment.