Skip to content

Commit

Permalink
Improved UIH support (#225)
Browse files Browse the repository at this point in the history
  • Loading branch information
@neurolabusc
neurolabusc committed Sep 30, 2018
1 parent 773016b commit 5af76a9
Show file tree
Hide file tree
Showing 3 changed files with 22 additions and 41 deletions.
2 changes: 1 addition & 1 deletion GE/README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ dcm2niix attempts to convert GE DICOM format images to NIfTI. The current genera

## Diffusion Tensor Notes

As noted by Jaemin Shin (GE), the GE convention for reported diffusion gradient direction has always been in “MR physics” logical coordinate, i.e Freq (X), Phase (Y), Slice (Z). Note that this is neither “with reference to the scanner bore” (like Siemens or Philips) nor “with reference to the imaging plane” (as expected by FSL tools). This is the main source of confusion. This explains why the dcm2niix function geCorrectBvecs() checks whether the DICOM tag In-plane Phase Encoding Direction (0018,1312) is 'ROW' or 'COL'. In addition, it will generate the warning 'reorienting for ROW phase-encoding untested' if you acquire DTI data with the phase encoding in the ROW direction. If you do test this feature, please report your findings as a Github issue.
The [NA-MIC Wiki]](https://www.na-mic.org/wiki/NAMIC_Wiki:DTI:DICOM_for_DWI_and_DTI#Private_vendor:_GE) provides a nice description of the GE diffusion tags. In brief, the B-value is stored as the first element in the array of 0043,1039. The DICOM elements 0019,10bb, 0019,10bc and 0019,10bd provide the gradient direction relative to the frequency, phase and slice. As noted by Jaemin Shin (GE), the GE convention for reported diffusion gradient direction has always been in “MR physics” logical coordinate, i.e Freq (X), Phase (Y), Slice (Z). Note that this is neither “with reference to the scanner bore” (like Siemens or Philips) nor “with reference to the imaging plane” (as expected by FSL tools). This is the main source of confusion. This explains why the dcm2niix function geCorrectBvecs() checks whether the DICOM tag In-plane Phase Encoding Direction (0018,1312) is 'ROW' or 'COL'. In addition, it will generate the warning 'reorienting for ROW phase-encoding untested' if you acquire DTI data with the phase encoding in the ROW direction. If you do test this feature, please report your findings as a Github issue. Assuming you have COL phase encoding, dcm2niix should provide [FSL format](http://justinblaber.org/brief-introduction-to-dwmri/) [bvec files](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FDT/FAQ#What_conventions_do_the_bvecs_use.3F).

## Slice Timing

Expand Down
6 changes: 3 additions & 3 deletions console/nii_dicom.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2292,7 +2292,8 @@ unsigned char * nii_demosaic(unsigned char* inImg, struct nifti_1_header *hdr, i
int nRow = nCol;
//n.b. Siemens store 20 images as 5x5 grid, UIH as 5rows, 4 Col https://github.com/rordenlab/dcm2niix/issues/225
if (isUIH)
nRow = ceil(nMosaicSlices/nCol);
nRow = ceil((float)nMosaicSlices/(float)nCol);
//printf("%d = %dx%d\n", nMosaicSlices, nCol, nRow);
int colBytes = hdr->dim[1]/nCol * hdr->bitpix/8;
int lineBytes = hdr->dim[1] * hdr->bitpix/8;
int rowBytes = hdr->dim[1] * hdr->dim[2]/nRow * hdr->bitpix/8;
Expand Down Expand Up @@ -4518,7 +4519,7 @@ double TE = 0.0; //most recent echo time recorded
if (d.manufacturer != kMANUFACTURER_UIH) break;
//UIH slice timing
d.CSA.sliceTiming[acquisitionTimesUIH] = d.acquisitionTime;
//printf("%g\n", d.CSA.sliceTiming[acquisitionTimesUIH]);
//printf("UIHsliceTime\t%s\t%0.8f\n", acquisitionTimeTxt, d.CSA.sliceTiming[acquisitionTimesUIH]);
acquisitionTimesUIH ++;
break;
case kContentTime :
Expand Down Expand Up @@ -4769,7 +4770,6 @@ double TE = 0.0; //most recent echo time recorded
case kPEDirectionDisplayedUIH :
if (d.manufacturer != kMANUFACTURER_UIH) break;
dcmStr (lLength, &buffer[lPos], d.phaseEncodingDirectionDisplayedUIH);
printf("---> %s\n", d.phaseEncodingDirectionDisplayedUIH);
break;
case kDiffusion_bValueUIH : {
if (d.manufacturer != kMANUFACTURER_UIH) break;
Expand Down
55 changes: 18 additions & 37 deletions console/nii_dicom_batch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -283,10 +283,8 @@ void siemensPhilipsCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI
+ (vx[i].V[2]*vx[i].V[2])
+ (vx[i].V[3]*vx[i].V[3]));
if ((vx[i].V[0] <= FLT_EPSILON)|| (vLen <= FLT_EPSILON) ) { //bvalue=0
if (vx[i].V[0] > FLT_EPSILON)
printWarning("Volume %d appears to be an ADC map (non-zero b-value with zero vector length)\n", i);
//for (int v= 0; v < 4; v++)
// vx[i].V[v] =0.0f;
if (vx[i].V[0] > 5.0) //Philip stores n.b. UIH B=1.25126 Vec=0,0,0 while Philips stored isotropic images
printWarning("Volume %d appears to be derived image ADC/Isotropic (non-zero b-value with zero vector length)\n", i);
continue; //do not normalize or reorient b0 vectors
}//if bvalue=0
vec3 bvecs_old =setVec3(vx[i].V[1],vx[i].V[2],vx[i].V[3]);
Expand Down Expand Up @@ -399,28 +397,6 @@ const void * memmem(const char *l, size_t l_len, const char *s, size_t s_len) {
//n.b. memchr returns "const void *" not "void *" for Windows C++ https://msdn.microsoft.com/en-us/library/d7zdhf37.aspx
#endif //for systems without memmem

/*int readKeyX(const char * key, char * buffer, int remLength) { //look for text key in binary data stream, return subsequent integer value
int ret = 0;
char *keyPos = (char *)memmem(buffer, remLength, key, strlen(key));
printWarning("<><><>\n");
if (!keyPos) return 0;
printWarning("<><> %d\n", strlen(keyPos));
int i = (int)strlen(key);
int numDigits = 0;
while( ( i< remLength) && (numDigits >= 0) ) {
printMessage("%c", keyPos[i]);
if( keyPos[i] >= '0' && keyPos[i] <= '9' ) {
ret = (10 * ret) + keyPos[i] - '0';
numDigits ++;
} else if (numDigits > 0)
numDigits = -1;
i++;
}
printWarning("---> %d\n", ret);
return ret;
} //readKey()
*/

int readKey(const char * key, char * buffer, int remLength) { //look for text key in binary data stream, return subsequent integer value
int ret = 0;
char *keyPos = (char *)memmem(buffer, remLength, key, strlen(key));
Expand Down Expand Up @@ -2783,6 +2759,7 @@ int nii_saveCrop(char * niiFilename, struct nifti_1_header hdr, unsigned char* i
}// nii_saveCrop()

float dicomTimeToSec (float dicomTime) {
//convert HHMMSS to seconds, 135300.024 -> 135259.731 are 0.293 sec apart
char acqTimeBuf[64];
snprintf(acqTimeBuf, sizeof acqTimeBuf, "%+013.5f", (double)dicomTime);
int ahour,amin;
Expand Down Expand Up @@ -2813,25 +2790,29 @@ void checkDateTimeOrder(struct TDICOMdata * d, struct TDICOMdata * d1) {
void checkSliceTiming(struct TDICOMdata * d, struct TDICOMdata * d1) {
//detect images with slice timing errors. https://github.com/rordenlab/dcm2niix/issues/126
if ((d->TR < 0.0) || (d->CSA.sliceTiming[0] < 0.0)) return; //no slice timing
int nSlices = 0;
while ((nSlices < kMaxEPI3D) && (d->CSA.sliceTiming[nSlices] >= 0.0))
nSlices++;
if (nSlices < 1) return;
if (d->manufacturer == kMANUFACTURER_UIH) {//convert HHMMSS to Sec
for (int i = 0; i < nSlices; i++)
d->CSA.sliceTiming[i] = dicomTimeToSec(d->CSA.sliceTiming[i]);
float minT = d->CSA.sliceTiming[0];
for (int i = 0; i < nSlices; i++)
if (d->CSA.sliceTiming[i] < minT) minT = d->CSA.sliceTiming[i];
for (int i = 0; i < nSlices; i++)
d->CSA.sliceTiming[i] = d->CSA.sliceTiming[i] - minT;
}
float minT = d->CSA.sliceTiming[0];
float maxT = minT;
for (int i = 0; i < kMaxEPI3D; i++) {
if (d->CSA.sliceTiming[i] < 0.0) break;
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;
if (d->manufacturer == kMANUFACTURER_UIH) //convert HHMMSS to Sec
for (int i = 0; i < kMaxEPI3D; i++)
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

0 comments on commit 5af76a9

Please sign in to comment.