Update explicit naming of DICOMs (#252), GE file naming changes (#476)

rordenlab · Jan 29, 2021 · 275969c · 275969c
1 parent e32fb1b
commit 275969c
Show file tree

Hide file tree

Showing 7 changed files with 97 additions and 103 deletions.
diff --git a/GE/README.md b/GE/README.md
@@ -34,10 +34,26 @@ Some sequences allow the user to interpolate images in plane (e.g. saving a 2D 6
 
 ## Total Readout Time
 
-One often wants to determine [echo spacing, bandwidth](https://support.brainvoyager.com/brainvoyager/functional-analysis-preparation/29-pre-processing/78-epi-distortion-correction-echo-spacing-and-bandwidth) and total read-out time for EPI data so they can be undistorted. Speifically, we are interested in FSL's definition of total read-out time, which may differ from the actual read-out time. FSL expects “the time from the middle of the first each to the middle of the last echo, as it would have been had partial k-space not been used”. So total read-out time is influenced by parallel acceleration factor, bandwidth, number of EPI lines, but not partial Fourier. For GE data we can use the Acquisition Matrix (0018,1310) in the phase-encoding direction, the in-plane acceleration ASSET R factor (0043,1083) and the Effective Echo Spacing (0043,102C).
+One often wants to determine [echo spacing, bandwidth](https://support.brainvoyager.com/brainvoyager/functional-analysis-preparation/29-pre-processing/78-epi-distortion-correction-echo-spacing-and-bandwidth) and total read-out time for EPI data so they can be undistorted. Speifically, we are interested in FSL's definition of total read-out time, which may differ from the actual read-out time. FSL expects “the time from the middle of the first each to the middle of the last echo, as it would have been had partial k-space not been used”. So total read-out time is influenced by parallel acceleration factor, bandwidth, number of EPI lines, but not partial Fourier. For GE data we can use the Acquisition Matrix (0018,1310) in the phase-encoding direction, the in-plane acceleration ASSET R factor (the reciprocal of this is stored as the first element of 0043,1083) and the Effective Echo Spacing (0043,102C). While GE does not tell us the [partial Fourier fraction](https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html), is does reveal if it is present with the ScanOptions (0018,1022) reporting [PFF](http://dicomlookup.com/lookup.asp?sw=Ttable&q=C.8-4) (in my experience, GE does not populate [(0018,9081)](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,9081))). While partial Fourier does not impact FSL's totalReadoutTime directly, it can interact with the number of lines acquired when combined with parallel imaging (the `Round_factor`).
+
+The formula for FSL's definition of TotalReadoutTime (in seconds) is:
+
+```
+TotalReadoutTime = ( ( ceil ((1/Round_factor) * PE_AcquisitionMatrix / Asset_R_factor ) * Round_factor) - 1 ] * EchoSpacing * 0.000001
+```
+
+Consider an example:
+```
+(0018,1310) US 128\0\0\128                                 #   8, 4 AcquisitionMatrix
+(0018,0022) CS [SAT_GEMS\MP_GEMS\EPI_GEMS\ACC_GEMS\PFF\FS] #  42, 6 ScanOptions
+(0043,102c) SS 636                                         #   2, 1 EchoSpacing
+(0043,1083) DS [0.666667\1]                                #  10, 2 Acceleration
+```
+
+From this we can derive:
 
 ```
-TotalReadoutTime = ( ( ceil (0.5 * PE_AcquisitionMatrix / Asset_R_factor ) * 2) - 1 ] * EchoSpacing
+ASSET= 1.5 PE_AcquisitionMatrix= 128 EchoSpacing= 636 Round_Factor= 4 TotalReadoutTime= 0.055332
 ```
 
 ## Image Acceleration

diff --git a/console/main_console.cpp b/console/main_console.cpp
@@ -102,7 +102,8 @@ void showHelp(const char * argv[], struct TDCMopts opts) {
     printf("  -p : Philips precise float (not display) scaling (y/n, default y)\n");
     printf("  -r : rename instead of convert DICOMs (y/n, default n)\n");
     printf("  -s : single file mode, do not convert other images in folder (y/n, default n)\n");
-    printf("  -t : text notes includes private patient details (y/n, default n)\n");
+	//text notes replaced with BIDS: this function is deprecated
+    //printf("  -t : text notes includes private patient details (y/n, default n)\n");
     #if !defined(_WIN64) && !defined(_WIN32) //shell script for Unix only
 	printf("  -u : up-to-date check\n");
 	#endif

diff --git a/console/nifti1_io_core.cpp b/console/nifti1_io_core.cpp
@@ -192,6 +192,18 @@ vec3 nifti_vect33_norm (vec3 v) { //normalize vector length
     return vO;
 }
 
+vec4 nifti_vect44_norm (vec4 v) { //normalize vector length
+    vec4 vO = v;
+    float vLen = sqrt( (v.v[0]*v.v[0])
+                      + (v.v[1]*v.v[1])
+                      + (v.v[2]*v.v[2])
+                      + (v.v[3]*v.v[3]));
+    if (vLen <= FLT_EPSILON) return vO; //avoid divide by zero
+    for (int i = 0; i < 4; i++)
+        vO.v[i] = v.v[i]/vLen;
+    return vO;
+}
+
 vec3 nifti_vect33mat33_mul(vec3 v, mat33 m ) { //multiply vector * 3x3matrix
     vec3 vO;
     for (int i=0; i<3; i++) { //multiply Pcrs * m

diff --git a/console/nifti1_io_core.h b/console/nifti1_io_core.h
@@ -68,6 +68,7 @@ mat44 nifti_mat44_inverse( mat44 R );
 mat44 nifti_mat44_mul( mat44 A , mat44 B );
 vec3 crossProduct(vec3 u, vec3 v);
 vec3 nifti_vect33_norm (vec3 v);
+vec4 nifti_vect44_norm (vec4 v);
 vec3 nifti_vect33mat33_mul(vec3 v, mat33 m );
 ivec3 setiVec3(int x, int y, int z);
 vec3 setVec3(float x, float y, float z);

diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp
@@ -5322,7 +5322,7 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 );
             	break;
             }
             case kProtocolName : {
-                dcmStr(lLength, &buffer[lPos], d.protocolName); //see also kSequenceName
+                dcmStr(lLength, &buffer[lPos], d.protocolName);
                 break; }
             case kPatientOrient :
                 dcmStr(lLength, &buffer[lPos], d.patientOrient);
@@ -6027,7 +6027,7 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 );
 				//Canon/Toshiba omits 0018,9018, but reports [SE\EP];[GR\EP] for 0018,0020
 				//GE omits 0018,9018, but reports [EP\GR];[EP\SE] for 0018,0020
 				//Philips reports 0018,9018, but reports [SE];[GR] for 0018,0020
-                if ((lLength > 1) && (strstr(d.scanningSequence, "IR") != NULL))
+				if ((lLength > 1) && (strstr(d.scanningSequence, "IR") != NULL))
                 	d.isIR = true;
                 if ((lLength > 1) && (strstr(d.scanningSequence, "EP") != NULL))
                 	d.isEPI = true;
@@ -6042,6 +6042,8 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 );
             }
             case kScanOptions:
             	dcmStr(lLength, &buffer[lPos], d.scanOptions);
+				if ((lLength > 1) && (strstr(d.scanOptions, "PFF") != NULL))
+                	d.isPartialFourier = true; //e.g. GE does not populate (0018,9081)			
             	break;
             case kSequenceName : {
                 dcmStr(lLength, &buffer[lPos], d.sequenceName);
@@ -6865,13 +6867,14 @@ uint32_t kSequenceDelimitationItemTag = 0xFFFE +(0xE0DD << 16 );
     	printWarning("0008,0008=MOSAIC but number of slices not specified: %s\n", fname);
     if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.CSA.dtiV[1] < -1.0) && (d.CSA.dtiV[2] < -1.0) && (d.CSA.dtiV[3] < -1.0))
     	d.CSA.dtiV[0] = 0; //SiemensTrio-Syngo2004A reports B=0 images as having impossible b-vectors.
-    if ((d.manufacturer == kMANUFACTURER_GE) && (d.modality == kMODALITY_MR) && (strlen(d.seriesDescription) > 1)) {//GE uses a generic session name here: do not overwrite kProtocolNameGE
-		//issue 473: swap protocolName and seriesDescription
+	/*
+    if ((d.manufacturer == kMANUFACTURER_GE) && (d.modality == kMODALITY_MR) && (strlen(d.seriesDescription) > 1)) {
+		//issue 473, 476: swap protocolName and seriesDescription
 		char tempTxt[kDICOMStr] = "";
 		strcpy(tempTxt, d.protocolName);
 		strcpy(d.protocolName, d.seriesDescription);
 		strcpy(d.seriesDescription, tempTxt);
-	}
+	}*/
     if ((strlen(d.protocolName) < 1) && (strlen(d.seriesDescription) > 1))
 		strcpy(d.protocolName, d.seriesDescription);
 	if ((strlen(d.protocolName) > 1) && (isMoCo))

diff --git a/console/nii_dicom.h b/console/nii_dicom.h
@@ -50,7 +50,7 @@ extern "C" {
     #define kCPUsuf " " //unknown CPU
 #endif
 
-#define kDCMdate "v1.0.20210124"
+#define kDCMdate "v1.0.20210129"
 #define kDCMvers kDCMdate " " kJP2suf kLSsuf kCCsuf kCPUsuf
 
 static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic

diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp
@@ -57,11 +57,9 @@
 #endif
 
 #if defined(_WIN64) || defined(_WIN32)
-	#define myTextFileInputLists //comment out to disable this feature: https://github.com/rordenlab/dcm2niix/issues/288
 	const char kPathSeparator ='\\';
 	const char kFileSep[2] = "\\";
 #else
-	#define myTextFileInputLists
 	const char kPathSeparator ='/';
 	const char kFileSep[2] = "/";
 #endif
@@ -6249,77 +6247,6 @@ int singleDICOM(struct TDCMopts* opts, char *fname) {
     return ret;
 }// singleDICOM()
 
-#ifdef myTextFileInputLists //https://github.com/rordenlab/dcm2niix/issues/288
-int textDICOM(struct TDCMopts* opts, char *fname) {
-	//check input file
-    FILE *fp = fopen(fname, "r");
-    if (fp == NULL)
-#ifdef USING_R
-        return EXIT_FAILURE;
-#else
-    	exit(EXIT_FAILURE);
-#endif
-	int nConvert = 0;
-    char dcmname[2048];
-    while (fgets(dcmname, sizeof(dcmname), fp)) {
-		int sz = strlen(dcmname);
-		if (sz > 0 && dcmname[sz-1] == '\n') dcmname[sz-1] = 0; //Unix LF
-		if (sz > 1 && dcmname[sz-2] == '\r') dcmname[sz-2] = 0; //Windows CR/LF
-		//if (isDICOMfile(dcmname) == 0) { //<- this will reject DICOM metadata not wrapped with a header
-        if ((!is_fileexists(dcmname)) || (!is_fileNotDir(dcmname)) ) { //<-this will accept meta data
-        	fclose(fp);
-        	printError("Problem with file '%s'\n", dcmname);
-        	return EXIT_FAILURE;
-    	}
-    	//printf("%s\n", dcmname);
-		nConvert ++;
-    }
-    fclose(fp);
-    if (nConvert < 1) {
-    	printError("No DICOM files found '%s'\n", dcmname);
-    	return EXIT_FAILURE;
-    }
-    printMessage("Found %d DICOM file(s)\n", nConvert);
-    #ifndef USING_R
-    fflush(stdout); //show immediately if run from MRIcroGL GUI
-    #endif
-    TDCMsort * dcmSort = (TDCMsort *)malloc(nConvert * sizeof(TDCMsort));
-	struct TDICOMdata *dcmList  = (struct TDICOMdata *)malloc(nConvert * sizeof(struct  TDICOMdata));
-    struct TDTI4D *dti4D  = (struct TDTI4D *)malloc(sizeof(struct  TDTI4D));
-    struct TSearchList nameList;
-    nameList.maxItems = nConvert; // larger requires more memory, smaller more passes
-    nameList.str = (char **) malloc((nameList.maxItems+1) * sizeof(char *)); //reserve one pointer (32 or 64 bits) per potential file
-    nameList.numItems = 0;
-	nConvert = 0;
-	fp = fopen(fname, "r");
-    while (fgets(dcmname, sizeof(dcmname), fp)) {
-		int sz = strlen(dcmname);
-		if (sz > 0 && dcmname[sz-1] == '\n') dcmname[sz-1] = 0; //Unix LF
-		if (sz > 1 && dcmname[sz-2] == '\r') dcmname[sz-2] = 0; //Windows CR/LF
-		nameList.str[nameList.numItems]  = (char *)malloc(strlen(dcmname)+1);
-    	strcpy(nameList.str[nameList.numItems],dcmname);
-    	nameList.numItems++;
-		dcmList[nConvert] = readDICOMv(nameList.str[nConvert], opts->isVerbose, opts->compressFlag, dti4D); //ignore compile warning - memory only freed on first of 2 passes
-		fillTDCMsort(dcmSort[nConvert], nConvert, dcmList[nConvert]);
-		nConvert ++;
-    }
-    fclose(fp);
-    qsort(dcmSort, nConvert, sizeof(struct TDCMsort), compareTDCMsort); //sort based on series and image numbers....
-	int ret = saveDcm2Nii(nConvert, dcmSort, dcmList, &nameList, *opts, dti4D);
-    free(dcmSort);
-    free(dcmList);
-	free(dti4D);
-    freeNameList(nameList);
-    return ret;
-}//textDICOM()
-
-#else //ifdef myTextFileInputLists
-int textDICOM(struct TDCMopts* opts, char *fname) {
-	printError("Unable to parse txt files: re-compile with 'myTextFileInputLists' (see issue 288)");
-	return EXIT_FAILURE;
-}
-#endif
-
 size_t fileBytes(const char * fname) {
     FILE *fp = fopen(fname, "rb");
 	if (!fp)  return 0;
@@ -6684,25 +6611,57 @@ int nii_loadDirCore(char *indir, struct TDCMopts* opts) {
     #ifdef myTimer
     clock_t start = clock();
     #endif
-	//1: find filenames of dicom files: up to two passes if we found more files than we allocated memory
-    for (int i = 0; i < 2; i++ ) {
-        nameList.str = (char **) malloc((nameList.maxItems+1) * sizeof(char *)); //reserve one pointer (32 or 64 bits) per potential file
-        nameList.numItems = 0;
-        searchDirForDICOM(indir, &nameList, opts->dirSearchDepth, 0, opts);
-        if (nameList.numItems <= nameList.maxItems)
-            break;
-        freeNameList(nameList);
-        nameList.maxItems = nameList.numItems+1;
-        //printMessage("Second pass required, found %ld images\n", nameList.numItems);
-    }
-    if (nameList.numItems < 1) {
-        if (opts->dirSearchDepth > 0)
-        	printError("Unable to find any DICOM images in %s (or subfolders %d deep)\n", indir, opts->dirSearchDepth);
-        else //keep silent for dirSearchDepth = 0 - presumably searching multiple folders
-        	{};
-        free(nameList.str); //ignore compile warning - memory only freed on first of 2 passes
-        return kEXIT_NO_VALID_FILES_FOUND;
-    }
+    if ((is_fileNotDir(opts->indir)) && isExt(opts->indir, ".txt") ){
+		nameList.str = (char **) malloc((nameList.maxItems+1) * sizeof(char *)); //reserve one pointer (32 or 64 bits) per potential file
+		nameList.numItems = 0;
+		FILE *fp = fopen(opts->indir, "r"); //textDICOM
+		if (fp == NULL)
+			return EXIT_FAILURE;
+		char dcmname[2048];
+		while (fgets(dcmname, sizeof(dcmname), fp)) {
+			int sz = (int)strlen(dcmname);
+			if (sz > 0 && dcmname[sz-1] == '\n') dcmname[sz-1] = 0; //Unix LF
+			if (sz > 1 && dcmname[sz-2] == '\r') dcmname[sz-2] = 0; //Windows CR/LF
+			if ((!is_fileexists(dcmname)) || (!is_fileNotDir(dcmname)) ) { //<-this will accept meta data
+				fclose(fp);
+				printError("Problem with file '%s'\n", dcmname);
+				return EXIT_FAILURE;
+			}
+			if (nameList.numItems < nameList.maxItems) {
+				nameList.str[nameList.numItems]  = (char *)malloc(strlen(dcmname)+1);
+				strcpy(nameList.str[nameList.numItems],dcmname);
+			}
+			nameList.numItems ++;
+		}
+		fclose(fp);
+		if (nameList.numItems >= nameList.maxItems) {
+			printError("Too many file names in '%s'\n", opts->indir);
+			return EXIT_FAILURE;
+		}
+		if (nameList.numItems < 1)
+			return kEXIT_NO_VALID_FILES_FOUND;
+		printMessage("Found %lu files in '%s'\n", nameList.numItems, opts->indir);
+    } else {
+		//1: find filenames of dicom files: up to two passes if we found more files than we allocated memory
+		for (int i = 0; i < 2; i++ ) {
+			nameList.str = (char **) malloc((nameList.maxItems+1) * sizeof(char *)); //reserve one pointer (32 or 64 bits) per potential file
+			nameList.numItems = 0;
+			searchDirForDICOM(indir, &nameList, opts->dirSearchDepth, 0, opts);
+			if (nameList.numItems <= nameList.maxItems)
+				break;
+			freeNameList(nameList);
+			nameList.maxItems = nameList.numItems+1;
+			//printMessage("Second pass required, found %ld images\n", nameList.numItems);
+		}
+		if (nameList.numItems < 1) {
+			if (opts->dirSearchDepth > 0)
+				printError("Unable to find any DICOM images in %s (or subfolders %d deep)\n", indir, opts->dirSearchDepth);
+			else //keep silent for dirSearchDepth = 0 - presumably searching multiple folders
+				{};
+			free(nameList.str); //ignore compile warning - memory only freed on first of 2 passes
+			return kEXIT_NO_VALID_FILES_FOUND;
+		}
+	}
     size_t nDcm = nameList.numItems;
     printMessage( "Found %lu DICOM file(s)\n", nameList.numItems); //includes images and other non-image DICOMs
     #ifdef myTimer
@@ -7043,8 +7002,10 @@ int nii_loadDir(struct TDCMopts* opts) {
             return convert_parRec(pname, *opts);
         };
     }
-    if (isFile && (opts->isOnlySingleFile) && isExt(indir, ".txt") )
-    	return textDICOM(opts, indir);
+	if (isFile && (opts->isOnlySingleFile) && isExt(indir, ".txt") ) {
+		strcpy(opts->indir,indir);
+		return nii_loadDirCore(opts->indir, opts);
+	}
 	if (opts->isRenameNotConvert) {
 		int nConvert = searchDirRenameDICOM(opts->indir, opts->dirSearchDepth, 0, opts);
 		if (nConvert < 0) return kEXIT_RENAME_ERROR;