[med-svn] [dcm2niix] 02/07: New upstream version 1.0.20171215
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Sun Feb 4 14:36:17 UTC 2018
This is an automated email from the git hooks/post-receive script.
ghisvail-guest pushed a commit to branch debian/master
in repository dcm2niix.
commit da6a7e864029464a35ffeb262451ede3e657c59b
Author: Ghislain Antony Vaillant <ghisvail at gmail.com>
Date: Sun Feb 4 14:10:46 2018 +0000
New upstream version 1.0.20171215
---
README.md | 5 +-
VERSIONS.md | 9 +
console/main_console.cpp | 69 +++++-
console/makefile | 11 +-
console/nii_dicom.cpp | 534 +++++++++++++++++++++++++++++++++++++-------
console/nii_dicom.h | 8 +-
console/nii_dicom_batch.cpp | 187 ++++++++++++----
console/nii_dicom_batch.h | 3 +
8 files changed, 695 insertions(+), 131 deletions(-)
diff --git a/README.md b/README.md
index 97eb97b..d1f42ec 100644
--- a/README.md
+++ b/README.md
@@ -39,7 +39,7 @@ mkdir build && cd build
cmake ..
make
```
-`dcm2niix` will be created in the `bin` subfolder. To install on the system run `make install`.
+`dcm2niix` will be created in the `bin` subfolder. To install on the system run `make install` instead of `make` - this will copy the executable to your path so you do not have to provide the full path to the executable.
**optional building with OpenJPEG:**
@@ -70,7 +70,8 @@ If you have any problems with the cmake build script described above or want to
- [dcm2niir](https://github.com/muschellij2/dcm2niir) R wrapper for dcm2niix/dcm2nii.
- [divest](https://github.com/jonclayden/divest) R interface to dcm2niix.
- [sci-tran dcm2niix](https://github.com/scitran-apps/dcm2niix) docker.
- - [neuro_docker](https://github.com/Neurita/neuro_docker) includes dcm2niix as part of a provides a single, static Dockerfile.
+ - [neuro_docker](https://github.com/Neurita/neuro_docker) includes dcm2niix as part of a single, static Dockerfile.
+ - [https://github.com/lalet/boutiques-dcm2niix boutiques-dcm2niix] is a dockerfile for installing and validating dcm2niix.
- [neurodocker](https://github.com/kaczmarj/neurodocker) generates [custom](https://github.com/rordenlab/dcm2niix/issues/138) Dockerfiles given specific versions of neuroimaging software.
- [dcm2niix_afni](https://afni.nimh.nih.gov/pub/dist/doc/program_help/dcm2niix_afni.html) is a version of dcm2niix included with the [AFNI](https://afni.nimh.nih.gov/) distribution.
- [MRIcroGL](https://github.com/neurolabusc/MRIcroGL) is available for MacOS, Linux and Windows and provides a graphical interface for dcm2niix. You can get compiled copies from the [MRIcroGL NITRC web site](https://www.nitrc.org/projects/mricrogl/).
\ No newline at end of file
diff --git a/VERSIONS.md b/VERSIONS.md
index 43ef2d4..0d8b611 100644
--- a/VERSIONS.md
+++ b/VERSIONS.md
@@ -1,5 +1,14 @@
## Versions
+15-Dec-2017
+ - Support [Siemens XA10 images](https://github.com/rordenlab/dcm2niix/pull/145).
+ - [Ability to select specific series to convert](https://github.com/rordenlab/dcm2niix/pull/146).
+
+4-Dec-2017
+ - Handle implicit VR DICOMs where [critical values nested in sequence groups (SQ)](https://github.com/rordenlab/dcm2niix/commit/7f5649c6fe6ed366d07776aa54397b50f6641aff)
+ - Better support for [PAR/REC files with segmented 3D EPI](https://github.com/rordenlab/dcm2niix/commit/66cdf2dcc60d55a6ef37f5a6db8d500d3eeb7c88).
+ - Allow Protocol Name to be [empty](https://github.com/rordenlab/dcm2niix/commit/94f3129898ba83bf310c9ff28e994f29feb13068).
+
17-Oct-2017
- Swap [phase-encoding direction polarity](https://github.com/rordenlab/dcm2niix/issues/125) for Siemens images where PE is in the Column direction.
- Sort diffusion volumes by [B-value amplitude](https://www.nitrc.org/forum/forum.php?thread_id=8396&forum_id=4703) (use "-d n"/"-d y" to turn the feature off/on).
diff --git a/console/main_console.cpp b/console/main_console.cpp
index 3608038..a784d2d 100644
--- a/console/main_console.cpp
+++ b/console/main_console.cpp
@@ -89,11 +89,15 @@ void showHelp(const char * argv[], struct TDCMopts opts) {
printf(" -h : show help\n");
printf(" -i : ignore derived, localizer and 2D images (y/n, default n)\n");
printf(" -m : merge 2D slices from same series regardless of study time, echo, coil, orientation, etc. (y/n, default n)\n");
+ printf(" -n : only convert this series number - can be used up to %i times (default convert all)\n", MAX_NUM_SERIES);
printf(" -o : output directory (omit to save to input folder)\n");
printf(" -p : Philips precise float (not display) scaling (y/n, default y)\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");
- printf(" -v : verbose (n/y or 0/1/2 [no, yes, logorrheic], default 0)\n");
+ #if !defined(_WIN64) && !defined(_WIN32) //shell script for Unix only
+ printf(" -u : up-to-date check\n");
+ #endif
+ printf(" -v : verbose (n/y or 0/1/2 [no, yes, logorrheic], default 0)\n");
printf(" -x : crop (y/n, default n)\n");
char gzCh = 'n';
if (opts.isGz) gzCh = 'y';
@@ -145,6 +149,53 @@ int invalidParam(int i, const char * argv[]) {
return 1;
}
+#if !defined(_WIN64) && !defined(_WIN32) //shell script for Unix only
+
+int checkUpToDate() {
+ #define URL "/rordenlab/dcm2niix/releases/"
+ #define APIURL "\"https://api.github.com/repos" URL "latest\""
+ #define HTMURL "https://github.com" URL
+ #define SHELLSCRIPT "#!/usr/bin/env bash\n curl --silent " APIURL " | grep '\"tag_name\":' | sed -E 's/.*\"([^\"]+)\".*/\\1/'"
+ //check first 13 characters, e.g. "v1.0.20171204"
+ #define versionChars 13
+ FILE *pipe = popen(SHELLSCRIPT, "r");
+ char ch, gitvers[versionChars+1];
+ int n = 0;
+ int nMatch = 0;
+ while ((ch = fgetc(pipe)) != EOF) {
+ if (n < versionChars) {
+ gitvers[n] = ch;
+ if (gitvers[n] == kDCMvers[n])
+ nMatch ++;
+ n ++;
+ }
+ }
+ pclose(pipe);
+ gitvers[n] = 0; //null terminate
+ if (n < 1) { //script reported nothing
+ printf("Error: unable to check version with script:\n %s\n", SHELLSCRIPT);
+ return 3; //different from EXIT_SUCCESS (0) and EXIT_FAILURE (1)
+ }
+ if (nMatch == versionChars) { //versions match
+ printf("Good news: Your version is up to date: %s\n", gitvers);
+ return EXIT_SUCCESS;
+ }
+ //report error
+ char myvers[versionChars+1];
+ for (int i = 0; i < versionChars; i++) myvers[i] = kDCMvers[i];
+ myvers[versionChars] = 0; //null terminate
+ int myv = atoi(myvers + 5); //skip "v1.0."
+ int gitv = atoi(gitvers + 5); //skip "v1.0."
+ if (myv > gitv) {
+ printf("Warning: your version ('%s') more recent than stable release ('%s')\n %s\n", myvers, gitvers, HTMURL);
+ return 2; //different from EXIT_SUCCESS (0) and EXIT_FAILURE (1)
+ }
+ printf("Error: your version ('%s') is not the latest release ('%s')\n %s\n", myvers, gitvers, HTMURL);
+ return EXIT_FAILURE;
+} //checkUpToDate()
+
+#endif //shell script for UNIX only
+
//#define mydebugtest
int main(int argc, const char * argv[])
{
@@ -244,6 +295,10 @@ int main(int argc, const char * argv[])
opts.isCreateText = false;
else
opts.isCreateText = true;
+ #if !defined(_WIN64) && !defined(_WIN32) //shell script for Unix only
+ } else if (argv[i][1] == 'u') {
+ return checkUpToDate();
+ #endif
} else if ((argv[i][1] == 'v') && ((i+1) < argc)) {
i++;
if (invalidParam(i, argv)) return 0;
@@ -279,6 +334,18 @@ int main(int argc, const char * argv[])
} else if ((argv[i][1] == 'o') && ((i+1) < argc)) {
i++;
strcpy(opts.outdir,argv[i]);
+ } else if ((argv[i][1] == 'n') && ((i+1) < argc)) {
+ i++;
+ int seriesNumber = atoi(argv[i]);
+ if (seriesNumber < 0)
+ opts.numSeries = -1; //report series: convert none
+ else if ((opts.numSeries >= 0) && (opts.numSeries < MAX_NUM_SERIES)) {
+ opts.seriesNumber[opts.numSeries] = seriesNumber;
+ opts.numSeries += 1;
+ }
+ else {
+ printf("Warning: too many series specified, ignoring -n %s\n", argv[i]);
+ }
} else
printf(" Error: invalid option '%s %s'\n", argv[i], argv[i+1]);;
lastCommandArg = i;
diff --git a/console/makefile b/console/makefile
index 3fe5ea4..4b0c688 100755
--- a/console/makefile
+++ b/console/makefile
@@ -1,9 +1,14 @@
-CFLAGS=-s
+# Regular use
+CFLAGS=-s -O3
+
+# Debugging
+#CFLAGS=-g
+
ifneq ($(OS),Windows_NT)
OS = $(shell uname)
ifeq "$(OS)" "Darwin"
- CFLAGS=-dead_strip
+ CFLAGS=-dead_strip -O3
endif
endif
all:
- g++ $(CFLAGS) -O3 -I. main_console.cpp nii_foreign.cpp nii_dicom.cpp jpg_0XC3.cpp ujpeg.cpp nifti1_io_core.cpp nii_ortho.cpp nii_dicom_batch.cpp -o dcm2niix -DmyDisableOpenJPEG
+ g++ $(CFLAGS) -I. main_console.cpp nii_foreign.cpp nii_dicom.cpp jpg_0XC3.cpp ujpeg.cpp nifti1_io_core.cpp nii_ortho.cpp nii_dicom_batch.cpp -o dcm2niix -DmyDisableOpenJPEG
diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp
index a4f7451..7ad2926 100644
--- a/console/nii_dicom.cpp
+++ b/console/nii_dicom.cpp
@@ -657,6 +657,8 @@ struct TDICOMdata clear_dicom_data() {
d.angulation[i] = 0.0f;
d.xyzMM[i] = 1;
}
+ for (int i=0; i < MAX_NUMBER_OF_DIMENSIONS; ++i)
+ d.dimensionIndexValues[i] = 0;
d.CSA.sliceTiming[0] = -1.0f; //impossible value denotes not known
d.CSA.numDti = 0;
for (int i=0; i < 5; i++)
@@ -840,7 +842,7 @@ void dcmStr(int lLength, unsigned char lBuffer[], char* lOut, bool isStrLarge =
if (c == 255) cString[i] = 'y';
}
for (int i = 0; i < lLength; i++)
- if ((cString[i]<1) || (cString[i]==' ') || (cString[i]==',') || (cString[i]=='^') || (cString[i]=='/') || (cString[i]=='\\') || (cString[i]=='%') || (cString[i]=='*')) cString[i] = '_';
+ if ((cString[i]<1) || (cString[i]==' ') || (cString[i]==',') || (cString[i]=='^') || (cString[i]=='/') || (cString[i]=='\\') || (cString[i]=='%') || (cString[i]=='*') || (cString[i] == 9) || (cString[i] == 10) || (cString[i] == 11) || (cString[i] == 13)) cString[i] = '_';
int len = 1;
for (int i = 1; i < lLength; i++) { //remove repeated "_"
if ((cString[i-1]!='_') || (cString[i]!='_')) {
@@ -888,7 +890,8 @@ float dcmFloat(int lByteLength, unsigned char lBuffer[], bool littleEndian) {//r
return swapVal;
} //dcmFloat()
-double dcmFloatDouble(int lByteLength, unsigned char lBuffer[], bool littleEndian) {//read binary 32-bit float
+double dcmFloatDouble(const size_t lByteLength, const unsigned char lBuffer[],
+ const bool littleEndian) {//read binary 32-bit float
//http://stackoverflow.com/questions/2782725/converting-float-values-from-big-endian-to-little-endian
bool swap = (littleEndian != littleEndianPlatform());
double retVal = 0.0f;
@@ -921,14 +924,14 @@ int dcmInt (int lByteLength, unsigned char lBuffer[], bool littleEndian) { //rea
return lBuffer[3]+(lBuffer[2]<<8)+(lBuffer[1]<<16)+(lBuffer[0]<<24); //shortint vs word?
} //dcmInt()
-int dcmStrInt (int lByteLength, unsigned char lBuffer[]) {//read float stored as a string
+int dcmStrInt (const int lByteLength, const unsigned char lBuffer[]) {//read float stored as a string
//#ifdef _MSC_VER
char * cString = (char *)malloc(sizeof(char) * (lByteLength + 1));
//#else
// char cString[lByteLength + 1];
//#endif
cString[lByteLength] =0;
- memcpy(cString, (char*)&lBuffer[0], lByteLength);
+ memcpy(cString, (const unsigned char*)(&lBuffer[0]), lByteLength);
//printMessage(" --> *%s* %s%s\n",cString, &lBuffer[0],&lBuffer[1]);
int ret = atoi(cString);
//#ifdef _MSC_VER
@@ -937,7 +940,7 @@ int dcmStrInt (int lByteLength, unsigned char lBuffer[]) {//read float stored as
return ret;
} //dcmStrInt()
-int dcmStrManufacturer (int lByteLength, unsigned char lBuffer[]) {//read float stored as a string
+int dcmStrManufacturer (const int lByteLength, unsigned char lBuffer[]) {//read float stored as a string
if (lByteLength < 2) return kMANUFACTURER_UNKNOWN;
//#ifdef _MSC_VER
char * cString = (char *)malloc(sizeof(char) * (lByteLength + 1));
@@ -1158,6 +1161,7 @@ int readCSAImageHeader(unsigned char *buff, int lLength, struct TCSAdata *CSA, i
}
} //if at least 1 item
}// for lT 1..lnTag
+ if (CSA->protocolSliceNumber1 > 1) CSA->sliceOrder = NIFTI_SLICE_UNKNOWN;
return EXIT_SUCCESS;
} // readCSAImageHeader()
@@ -1170,6 +1174,16 @@ void dcmMultiShorts (int lByteLength, unsigned char lBuffer[], int lnShorts, uin
nifti_swap_2bytes(lnShorts, &lShorts[0]);
} //dcmMultiShorts()
+void dcmMultiLongs (int lByteLength, unsigned char lBuffer[], int lnLongs, uint32_t *lLongs, bool littleEndian) {
+ //read array of unsigned longs UL http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_6.2.html
+ if((lnLongs < 1) || (lByteLength != (lnLongs * 4)))
+ return;
+ memcpy(&lLongs[0], (uint32_t *)&lBuffer[0], lByteLength);
+ bool swap = (littleEndian != littleEndianPlatform());
+ if (swap)
+ nifti_swap_4bytes(lnLongs, &lLongs[0]);
+} //dcmMultiLongs()
+
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;
@@ -1204,7 +1218,7 @@ void dcmMultiFloat (int lByteLength, char lBuffer[], int lnFloats, float *lFloat
//#endif
} //dcmMultiFloat()
-float dcmStrFloat (int lByteLength, unsigned char lBuffer[]) { //read float stored as a string
+float dcmStrFloat (const int lByteLength, const unsigned char lBuffer[]) { //read float stored as a string
//#ifdef _MSC_VER
char * cString = (char *)malloc(sizeof(char) * (lByteLength + 1));
//#else
@@ -1375,19 +1389,29 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D
#define kYmm 29
#define kTEcho 30
#define kDynTime 31
-#define kGradientNumber 42
#define kbval 33
+#define kInversionDelayMs 40
+#define kGradientNumber 42
#define kv1 47
#define kv2 45
#define kv3 46
#define kASL 48
char buff[LINESZ];
+ int sliceNumberMrPhilipsB2[kMaxDTI4D], sliceNumberMrPhilipsVol2[kMaxDTI4D];
+ int patientPositionNumPhilipsB2 = 0;
+ int patientPositionNumPhilipsVol2 = 0;
//float intenScalePhilips = 0.0f;
float maxBValue = 0.0f;
float maxDynTime = 0.0f;
float minDynTime = 999999.0f;
+ int minDyn = 32767;
+ int maxDyn = 0;
bool ADCwarning = false;
+ int prevDyn = -1;
+ bool dynNotAscending = false;
int parVers = 0;
+ int maxEcho = 1;
+ int maxCardiac = 1;
int nCols = 26;
int slice = 0;
//int prevSliceIndex = 0; //index of prior slice: detect if images are not in order
@@ -1426,7 +1450,6 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D
if (buff[0] == '.') { //tag
char Comment[8][50];
sscanf(buff, ". %s %s %s %s %s %s %s %s\n", Comment[0], Comment[1],Comment[2], Comment[3], Comment[4], Comment[5], Comment[6], Comment[7]);
-
if ((strcmp(Comment[0], "Acquisition") == 0) && (strcmp(Comment[1], "nr") == 0)) {
d.acquNum = atoi( Comment[3]);
d.seriesNum = d.acquNum;
@@ -1490,6 +1513,8 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D
if (slice == 1) {
//for (int i = 0; i < nCols; i++)
// cols1[i] = cols[i]; //store first slice to see if dimensions or intensity scale varies between slices
+ //for (int i = 0; i < nCols; i++)
+ // printMessage("%d %g\n",i, cols[i]); //store first slice to see if dimensions or intensity scale varies between slices
d.xyzDim[1] = (int) cols[kXdim];
d.xyzDim[2] = (int) cols[kYdim];
d.xyzMM[1] = cols[kXmm];
@@ -1503,6 +1528,7 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D
d.angulation[3] = cols[kAngulationFHs];
d.sliceOrient = (int) cols[kSliceOrients];
d.TE = cols[kTEcho];
+ d.TI = cols[kInversionDelayMs];
d.bitsAllocated = (int) cols[kBitsPerVoxel];
d.bitsStored = (int) cols[kBitsPerVoxel];
d.intenIntercept = cols[kRI];
@@ -1523,8 +1549,22 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D
}
if (cols[kImageType] == 0) d.isHasMagnitude = true;
if (cols[kImageType] != 0) d.isHasPhase = true;
+ if (cols[kDyn] > maxDyn) maxDyn = (int) cols[kDyn];
+ if (cols[kDyn] < minDyn) minDyn = (int) cols[kDyn];
+ if (cols[kDyn] < prevDyn) dynNotAscending = true;
+ prevDyn = cols[kDyn];
if (cols[kDynTime] > maxDynTime) maxDynTime = cols[kDynTime];
if (cols[kDynTime] < minDynTime) minDynTime = cols[kDynTime];
+ if (cols[kEcho] > maxEcho) maxEcho = cols[kEcho];
+ if (cols[kCardiac] > maxCardiac) maxCardiac = cols[kCardiac];
+ if ((cols[kEcho] == 1) && (cols[kDyn] == 1) && (patientPositionNumPhilipsB2 < kMaxDTI4D) && (cols[kCardiac] == 1) && (cols[kGradientNumber] == 2)) {
+ sliceNumberMrPhilipsB2[patientPositionNumPhilipsB2] = round(cols[kSlice]);
+ patientPositionNumPhilipsB2++;
+ }
+ if ((cols[kEcho] == 1) && (cols[kDyn] == 2) && (patientPositionNumPhilipsVol2 < kMaxDTI4D) && (cols[kCardiac] == 1) && (cols[kGradientNumber] == 1)) {
+ sliceNumberMrPhilipsVol2[patientPositionNumPhilipsVol2] = round(cols[kSlice]);
+ patientPositionNumPhilipsVol2++;
+ }
if ((cols[kEcho] == 1) && (cols[kDyn] == 1) && (cols[kCardiac] == 1) && (cols[kGradientNumber] == 1)) {
if (cols[kSlice] == 1) {
d.patientPosition[1] = cols[kPositionRL];
@@ -1548,14 +1588,16 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D
//it seems like the ADC is always saved as the final volume, so this solution SHOULD be foolproof.
ADCwarning = true;
}*/
- if (cols[kSlice] == 1) { //only first slice
+ //666: test (cols[kDyn] == 1)
+ //(cols[kImageType] == 0) means magnitude scan
+ if ((cols[kImageType] == 0) && (cols[kDyn] == 1) && (cols[kEcho] == 1) && (cols[kCardiac] == 1) && (cols[kSlice] == 1)) { //only first slice
d.CSA.numDti++;
int dir = d.CSA.numDti;
if (dir <= kMaxDTI4D) {
if (isVerbose ) {
if (d.CSA.numDti == 1) printMessage("n\tdir\tbValue\tV1\tV2\tV3\n");
printMessage("%d\t%g\t%g\t%g\t%g\t%g\n", dir-1, cols[kGradientNumber], cols[kbval], cols[kv1], cols[kv2], cols[kv3]);
- }
+ }
dti4D->S[dir-1].V[0] = cols[kbval];
dti4D->S[dir-1].V[1] = cols[kv1];
dti4D->S[dir-1].V[2] = cols[kv2];
@@ -1568,15 +1610,31 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D
//printMessage("%f %f %lu\n",cols[9],cols[kGradientNumber], strlen(buff))
p = fgets (buff, LINESZ, fp);//get next line
}
-
free (cols);
fclose (fp);
d.manufacturer = kMANUFACTURER_PHILIPS;
d.isValid = true;
d.isSigned = true;
+ if ((patientPositionNumPhilipsB2 > 1) && (patientPositionNumPhilipsVol2 < 1)) {
+ patientPositionNumPhilipsVol2 = patientPositionNumPhilipsB2;
+ for (int s = 0; s < patientPositionNumPhilipsVol2; s++)
+ sliceNumberMrPhilipsVol2[s] = sliceNumberMrPhilipsB2[s];
+ }
+ if ((patientPositionNumPhilipsVol2 > 1) && (d.patientPositionNumPhilips == patientPositionNumPhilipsVol2)) {
+ bool isSliceOrderConsistent = true;
+ for (int s = 0; s < patientPositionNumPhilipsVol2; s++)
+ if (sliceNumberMrPhilipsVol2[s] != dti4D->S[s].sliceNumberMrPhilips) isSliceOrderConsistent = false;
+ if (!isSliceOrderConsistent)
+ printError("PAR file order of slices varies between volumes (hint: use dicm2nii)\n");
+ d.isValid = false;
+ }
+ if (dynNotAscending) {
+ printError("PAR file volumes not saved in ascending order (hint: use dicm2nii)\n");
+ d.isValid = false;
+ }
if ((slice % d.xyzDim[3]) != 0) {
- printError("Total number of slices (%d) not divisible by slices per 3D volume (%d) [acquisition aborted]. Try nii_rescue_par to fix this: %s\n", slice, d.xyzDim[3], parname);
- d.isValid = true;
+ printError("Total number of slices (%d) not divisible by slices per 3D volume (%d) [acquisition aborted]. Try dicm2nii or nii_rescue_par to fix this: %s\n", slice, d.xyzDim[3], parname);
+ d.isValid = false;
}
d.xyzDim[4] = slice/d.xyzDim[3];
d.locationsInAcquisition = d.xyzDim[3];
@@ -1665,8 +1723,10 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D
printError("Unable to convert DTI [increase kMaxDTI4D]\n");
d.CSA.numDti = 0;
};
- if ((maxBValue <= 0.0f) && (maxDynTime > minDynTime) && (d.CSA.numDti > 1)) {
- float TRms = 1000.0f * (maxDynTime - minDynTime) / (float)(d.CSA.numDti-1);
+ if ((maxBValue <= 0.0f) && (maxDyn > minDyn) && (maxDynTime > minDynTime)) { //use max vs min Dyn instead of && (d.CSA.numDti > 1)
+ int numDyn = maxDyn - minDyn;
+ float TRms = 1000.0f * (maxDynTime - minDynTime) / (float)numDyn;
+ //float TRms = 1000.0f * (maxDynTime - minDynTime) / (float)(d.CSA.numDti-1);
if (fabs(TRms - d.TR) > 0.005f)
printWarning("Reported TR=%gms, measured TR=%gms (prospect. motion corr.?)\n", d.TR, TRms);
d.TR = TRms;
@@ -1684,6 +1744,7 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D
if ((!v1varies) || (!v2varies) || (!v3varies))
printError("Bizarre b-vectors %s\n", parname);
}
+ if ((maxEcho > 1) || (maxCardiac > 1)) printWarning("Multiple Echo (%d) or Cardiac (%d). Segment output, e.g. nii_segment4d('img.nii', %d)\n", maxEcho, maxCardiac, maxEcho*maxCardiac);
return d;
} //nii_readParRec()
@@ -1843,6 +1904,7 @@ unsigned char * nii_reorderSlices(unsigned char* bImg, struct nifti_1_header *h,
memcpy(&srcImg[0], &bImg[volStart], volBytes); //dest, src, size
for (int z = 0; z < h->dim[3]; z++) { //for each slice
int src = dti4D->S[z].sliceNumberMrPhilips - 1; //-1 as Philips indexes slices from 1 not 0
+ //printMessage("Reordering volume %d slice %d\n", v, dti4D->S[z].sliceNumberMrPhilips);
if ((src < 0) || (src >= h->dim[3])) continue;
memcpy(&bImg[volStart+(src*sliceBytes)], &srcImg[z*sliceBytes], sliceBytes); //dest, src, size
}
@@ -2621,6 +2683,17 @@ unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct
return img;
} //nii_loadImgXL()
+int isSQ(uint32_t groupElement) { //Detect sequence VR ("SQ") for implicit tags
+ static const int array_size = 34;
+ uint32_t array[array_size] = {0x2005+(uint32_t(0x140F)<<16), 0x0008+(uint32_t(0x1111)<<16), 0x0008+(uint32_t(0x1115)<<16), 0x0008+(uint32_t(0x1140)<<16), 0x0008+(uint32_t(0x1199)<<16), 0x0008+(uint32_t(0x2218)<<16), 0x0008+(uint32_t(0x9092)<<16), 0x0018+(uint32_t(0x9006)<<16), 0x0018+(uint32_t(0x9042)<<16), 0x0018+(uint32_t(0x9045)<<16), 0x0018+(uint32_t(0x9049)<<16), 0x0018+(uint32_t(0x9112)<<16), 0x0018+(uint32_t(0x9114)<<16), 0x0018+(uint32_t(0x9115)<<16), 0x0018+(uint32_t(0x9119) [...]
+ for (int i = 0; i < array_size; i++) {
+ //if (array[i] == groupElement) printMessage(" implicitSQ %04x,%04x\n", groupElement & 65535,groupElement>>16);
+ if (array[i] == groupElement)
+ return 1;
+ }
+ return 0;
+} //isSQ()
+
int isDICOMfile(const char * fname) { //0=NotDICOM, 1=DICOM, 2=Maybe(not Part 10 compliant)
FILE *fp = fopen(fname, "rb");
if (!fp) return 0;
@@ -2642,6 +2715,196 @@ int isDICOMfile(const char * fname) { //0=NotDICOM, 1=DICOM, 2=Maybe(not Part 10
return 0;
} //isDICOMfile()
+//START RIR 12/2017 Robert I. Reid
+
+// Gathering spot for all the info needed to get the b value and direction
+// for a volume.
+struct TVolumeDiffusion {
+ struct TDICOMdata* pdd; // The multivolume
+ struct TDTI4D* pdti4D; // permanent records.
+
+ uint8_t manufacturer; // kMANUFACTURER_UNKNOWN, kMANUFACTURER_SIEMENS, etc.
+
+ //void set_manufacturer(const uint8_t m) {manufacturer = m; update();} // unnecessary
+
+ // Everything after this in the structure would be private if it were a C++
+ // class, but it has been rewritten as a struct for C compatibility. I am
+ // using _ as a hint of that, although _ for privacy is not really a
+ // universal convention in C. Privacy is desired because immediately
+ // any of these are updated _update_tvd() should be called.
+
+ bool _isAtFirstPatientPosition; // Limit b vals and vecs to 1 per volume.
+
+ //float bVal0018_9087; // kDiffusion_b_value, always present in Philips/Siemens.
+ //float bVal2001_1003; // kDiffusionBFactor
+ // float dirRL2005_10b0; // kDiffusionDirectionRL
+ // float dirAP2005_10b1; // kDiffusionDirectionAP
+ // float dirFH2005_10b2; // kDiffusionDirectionFH
+
+ // Philips diffusion scans tend to have a "trace" (average of the diffusion
+ // weighted volumes) volume tacked on, usually but not always at the end,
+ // so b is > 0, but the direction is meaningless. Most software versions
+ // explicitly set the direction to 0, but version 3 does not, making (0x18,
+ // 0x9075) necessary.
+ bool _isPhilipsNonDirectional;
+
+ //char _directionality0018_9075[16]; // DiffusionDirectionality, not in Philips 2.6.
+ // float _orientation0018_9089[3]; // kDiffusionOrientation, always
+ // // present in Philips/Siemens for
+ // // volumes with a direction.
+ //char _seq0018_9117[64]; // MRDiffusionSequence, not in Philips 2.6.
+
+ float _dtiV[4];
+ //uint16_t numDti;
+};
+struct TVolumeDiffusion initTVolumeDiffusion(struct TDICOMdata* ptdd, struct TDTI4D* dti4D);
+void clear_volume(struct TVolumeDiffusion* ptvd); // Blank the volume-specific members or set them to impossible values.
+void set_directionality0018_9075(struct TVolumeDiffusion* ptvd, unsigned char* inbuf);
+void set_orientation0018_9089(struct TVolumeDiffusion* ptvd, int lLength, unsigned char* inbuf,
+ bool isLittleEndian);
+void set_isAtFirstPatientPosition_tvd(struct TVolumeDiffusion* ptvd, bool iafpp);
+void set_bValGE(struct TVolumeDiffusion* ptvd, int lLength, unsigned char* inbuf);
+void set_diffusion_directionGE(struct TVolumeDiffusion* ptvd, int lLength, unsigned char* inbuf, int axis);
+void set_bVal(struct TVolumeDiffusion* ptvd, float b);
+void _update_tvd(struct TVolumeDiffusion* ptvd);
+
+struct TVolumeDiffusion initTVolumeDiffusion(struct TDICOMdata* ptdd, struct TDTI4D* dti4D) {
+ struct TVolumeDiffusion tvd;
+ tvd.pdd = ptdd;
+ tvd.pdti4D = dti4D;
+ clear_volume(&tvd);
+ return tvd;
+} //initTVolumeDiffusion()
+
+void clear_volume(struct TVolumeDiffusion* ptvd) {
+ ptvd->_isAtFirstPatientPosition = false;
+ ptvd->manufacturer = kMANUFACTURER_UNKNOWN;
+ //bVal0018_9087 = -1;
+ //ptvd->_directionality0018_9075[0] = 0;
+ //ptvd->seq0018_9117[0] = 0;
+ //bVal2001_1003 = -1;
+ // dirRL2005_10b0 = 2;
+ // dirAP2005_10b1 = 2;
+ // dirFH2005_10b2 = 2;
+ ptvd->_isPhilipsNonDirectional = false;
+ ptvd->_dtiV[0] = -1;
+ for(int i = 1; i < 4; ++i)
+ ptvd->_dtiV[i] = 2;
+ //numDti = 0;
+}//clear_volume()
+
+void set_directionality0018_9075(struct TVolumeDiffusion* ptvd, unsigned char* inbuf) {
+ if(strncmp(( char*)(inbuf), "DIRECTIONAL", 11) && // strncmp = 0 for ==.
+ strncmp(( char*)(inbuf), "BMATRIX", 7)){ // Siemens XA10
+ ptvd->_isPhilipsNonDirectional = true;
+ // Explicitly set the direction to 0 now, because there may
+ // not be a 0018,9089 for this frame.
+ for(int i = 1; i < 4; ++i) // 1-3 is intentional.
+ ptvd->_dtiV[i] = 0.0;
+ }
+ else{
+ ptvd->_isPhilipsNonDirectional = false;
+ // Wait for 0018,9089 to get the direction.
+ }
+
+ _update_tvd(ptvd);
+} //set_directionality0018_9075()
+
+void set_bValGE(struct TVolumeDiffusion* ptvd, int lLength, unsigned char* inbuf) {
+ //int dcmStrInt (int lByteLength, unsigned char lBuffer[]) {//read float stored as a string
+dcmStrInt(lLength, inbuf);
+
+ ptvd->_dtiV[0] = dcmStrInt(lLength, inbuf);
+ //dd.CSA.numDti = 1; // Always true for GE.
+
+ _update_tvd(ptvd);
+} //set_bValGE()
+
+// axis: 0 -> x, 1 -> y , 2 -> z
+void set_diffusion_directionGE(struct TVolumeDiffusion* ptvd, int lLength, unsigned char* inbuf, const int axis){
+ ptvd->_dtiV[axis + 1] = dcmStrFloat(lLength, inbuf);
+
+ _update_tvd(ptvd);
+}//set_diffusion_directionGE()
+
+void dcmMultiFloatDouble (size_t lByteLength, unsigned char lBuffer[], size_t lnFloats, float *lFloats, bool isLittleEndian) {
+ size_t floatlen = lByteLength / lnFloats;
+
+ for(size_t i = 0; i < lnFloats; ++i)
+ lFloats[i] = dcmFloatDouble((int)floatlen, lBuffer + i * floatlen, isLittleEndian);
+} //dcmMultiFloatDouble()
+
+void set_orientation0018_9089(struct TVolumeDiffusion* ptvd, int lLength, unsigned char* inbuf, bool isLittleEndian) {
+ if(ptvd->_isPhilipsNonDirectional){
+ for(int i = 1; i < 4; ++i) // Deliberately ignore inbuf; it might be nonsense.
+ ptvd->_dtiV[i] = 0.0;
+ }
+ else
+ dcmMultiFloatDouble(lLength, inbuf, 3, ptvd->_dtiV + 1, isLittleEndian);
+ _update_tvd(ptvd);
+}//set_orientation0018_9089()
+
+void set_bVal(struct TVolumeDiffusion* ptvd, const float b) {
+ ptvd->_dtiV[0] = b;
+ _update_tvd(ptvd);
+}//set_bVal()
+
+void set_isAtFirstPatientPosition_tvd(struct TVolumeDiffusion* ptvd, const bool iafpp) {
+ ptvd->_isAtFirstPatientPosition = iafpp;
+ _update_tvd(ptvd);
+}//set_isAtFirstPatientPosition_tvd()
+
+// Update the diffusion info in dd and *pdti4D for a volume once all the
+// diffusion info for that volume has been read into pvd.
+//
+// Note that depending on the scanner software the diffusion info can arrive in
+// different tags, in different orders (because of enclosing sequence tags),
+// and the values in some tags may be invalid, or may be essential, depending
+// on the presence of other tags. Thus it is best to gather all the diffusion
+// info for a volume (frame) before taking action on it.
+//
+// On the other hand, dd and *pdti4D need to be updated as soon as the
+// diffusion info is ready, before diffusion info for the next volume is read
+// in.
+void _update_tvd(struct TVolumeDiffusion* ptvd) {
+ // Figure out if we have both the b value and direction (if any) for this
+ // volume, and if isFirstPosition.
+
+ // // GE (software version 27) is liable to NOT include kDiffusion_b_value for the
+ // // slice if it is 0, but should still have kDiffusionBFactor, which comes
+ // // after PatientPosition.
+ // if(isAtFirstPatientPosition && manufacturer == kMANUFACTURER_GE && dtiV[0] < 0)
+ // dtiV[0] = 0; // Implied 0.
+
+ bool isReady = (ptvd->_isAtFirstPatientPosition && (ptvd->_dtiV[0] >= 0));
+ if(isReady){
+ for(int i = 1; i < 4; ++i){
+ if(ptvd->_dtiV[i] > 1){
+ isReady = false;
+ break;
+ }
+ }
+ }
+ if(!isReady)
+ return;
+
+ // If still here, update dd and *pdti4D.
+ ptvd->pdd->CSA.numDti++;
+ if (ptvd->pdd->CSA.numDti == 2) { // First time we know that this is a 4D DTI dataset;
+ for(int i = 0; i < 4; ++i) // Start *pdti4D before ptvd->pdd->CSA.dtiV
+ ptvd->pdti4D->S[0].V[i] = ptvd->pdd->CSA.dtiV[i]; // is updated.
+ }
+ for(int i = 0; i < 4; ++i) // Update pdd
+ ptvd->pdd->CSA.dtiV[i] = ptvd->_dtiV[i];
+ if((ptvd->pdd->CSA.numDti > 1) && (ptvd->pdd->CSA.numDti < kMaxDTI4D)){ // Update *pdti4D
+ //d.dti4D = (TDTI *)malloc(kMaxDTI4D * sizeof(TDTI));
+ for(int i = 0; i < 4; ++i)
+ ptvd->pdti4D->S[ptvd->pdd->CSA.numDti - 1].V[i] = ptvd->_dtiV[i];
+ }
+ clear_volume(ptvd); // clear the slate for the next volume.
+}//_update_tvd()
+//END RIR
+
struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, struct TDTI4D *dti4D) {
struct TDICOMdata d = clear_dicom_data();
strcpy(d.protocolName, ""); //erase dummy with empty
@@ -2649,6 +2912,9 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
strcpy(d.seriesDescription, ""); //erase dummy with empty
strcpy(d.sequenceName, ""); //erase dummy with empty
//do not read folders - code specific to GCC (LLVM/Clang seems to recognize a small file size)
+
+ struct TVolumeDiffusion volDiffusion = initTVolumeDiffusion(&d, dti4D);
+
struct stat s;
if( stat(fname,&s) == 0 ) {
if( !(s.st_mode & S_IFREG) ){
@@ -2716,6 +2982,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
#define kUnused 0x0001+(0x0001 << 16 )
#define kStart 0x0002+(0x0000 << 16 )
#define kTransferSyntax 0x0002+(0x0010 << 16)
+//#define kImplementationVersionName 0x0002+(0x0013 << 16)
#define kSourceApplicationEntityTitle 0x0002+(0x0016 << 16 )
//#define kSpecificCharacterSet 0x0008+(0x0005 << 16 ) //someday we should handle foreign characters...
#define kImageTypeTag 0x0008+(0x0008 << 16 )
@@ -2772,12 +3039,19 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
#define kInPlanePhaseEncodingDirection 0x0018+(0x1312<< 16 ) //CS
#define kSAR 0x0018+(0x1316 << 16 ) //'DS' 'SAR'
#define kPatientOrient 0x0018+(0x5100<< 16 ) //0018,5100. patient orientation - 'HFS'
+#define kDiffusionDirectionality 0x0018+uint32_t(0x9075<< 16 ) // NONE, ISOTROPIC, or DIRECTIONAL
//#define kDiffusionBFactorSiemens 0x0019+(0x100C<< 16 ) // 0019;000C;SIEMENS MR HEADER ;B_value
+#define kDiffusion_bValue 0x0018+uint32_t(0x9087<< 16 ) // FD
+#define kDiffusionOrientation 0x0018+uint32_t(0x9089<< 16 ) // FD, seen in enhanced
+ // DICOM from Philips 5.*
+ // and Siemens XA10.
#define kDwellTime 0x0019+(0x1018<< 16 ) //IS in NSec, see https://github.com/rordenlab/dcm2niix/issues/127
#define kLastScanLoc 0x0019+(0x101B<< 16 )
#define kDiffusionDirectionGEX 0x0019+(0x10BB<< 16 ) //DS
#define kDiffusionDirectionGEY 0x0019+(0x10BC<< 16 ) //DS
#define kDiffusionDirectionGEZ 0x0019+(0x10BD<< 16 ) //DS
+#define kSharedFunctionalGroupsSequence 0x5200+uint32_t(0x9229<< 16 ) // SQ
+#define kPerFrameFunctionalGroupsSequence 0x5200+uint32_t(0x9230<< 16 ) // SQ
#define kBandwidthPerPixelPhaseEncode 0x0019+(0x1028<< 16 ) //FD
#define kStudyID 0x0020+(0x0010 << 16 )
#define kSeriesNum 0x0020+(0x0011 << 16 )
@@ -2785,11 +3059,12 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
#define kImageNum 0x0020+(0x0013 << 16 )
#define kStudyInstanceUID 0x0020+(0x000D << 16 )
#define kSeriesInstanceUID 0x0020+(0x000E << 16 )
-#define kPatientPosition 0x0020+(0x0032 << 16 )
+#define kPatientPosition 0x0020+(0x0032 << 16 ) // Actually ImagePositionPatient!
#define kOrientationACR 0x0020+(0x0035 << 16 )
#define kOrientation 0x0020+(0x0037 << 16 )
#define kImagesInAcquisition 0x0020+(0x1002 << 16 ) //IS
#define kImageComments 0x0020+(0x4000<< 16 )// '0020' '4000' 'LT' 'ImageComments'
+#define kDimensionIndexValues 0x0020+uint32_t(0x9157<< 16 ) // UL n-dimensional index of frame.
#define kLocationsInAcquisitionGE 0x0021+(0x104F<< 16 )// 'SS' 'LocationsInAcquisitionGE'
#define kSamplesPerPixel 0x0028+(0x0002 << 16 )
#define kPhotometricInterpretation 0x0028+(0x0004 << 16 )
@@ -2838,7 +3113,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
#define kDiffusionDirectionRL 0x2005+(0x10B0 << 16)
#define kDiffusionDirectionAP 0x2005+(0x10B1 << 16)
#define kDiffusionDirectionFH 0x2005+(0x10B2 << 16)
-#define k2005140F 0x2005+(0x140F << 16)
+#define kPrivatePerFrameSq 0x2005+(0x140F << 16)
#define kWaveformSq 0x5400+(0x0100 << 16)
#define kImageStart 0x7FE0+(0x0010 << 16 )
#define kImageStartFloat 0x7FE0+(0x0008 << 16 )
@@ -2864,6 +3139,7 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
//float intenScalePhilips = 0.0;
char acquisitionDateTimeTxt[kDICOMStr] = "";
bool isEncapsulatedData = false;
+ int multiBandFactor = 0;
int encapsulatedDataFragments = 0;
int encapsulatedDataFragmentStart = 0; //position of first FFFE,E000 for compressed images
int encapsulatedDataImageStart = 0; //position of 7FE0,0010 for compressed images (where actual image start should be start of first fragment)
@@ -2881,6 +3157,7 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
float patientPosition[4] = {NAN, NAN, NAN, NAN}; //used to compute slice direction for Philips 4D
float patientPositionEndPhilips[4] = {NAN, NAN, NAN, NAN};
float patientPositionStartPhilips[4] = {NAN, NAN, NAN, NAN};
+
while ((d.imageStart == 0) && ((lPos+8+lFileOffset) < fileLen)) {
#ifndef myLoadWholeFileToReadHeader //read one segment at a time
if ((size_t)(lPos + 128) > MaxBufferSz) { //avoid overreading the file
@@ -2940,7 +3217,8 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
lPos = lPos + 4; //skip 4 byte length
} else if ((buffer[lPos] == 'S') && (buffer[lPos+1] == 'Q')) {
lLength = 8; //Sequence Tag
- is2005140FSQ = (groupElement == k2005140F);
+ //printMessage(" !!!SQ\t%04x,%04x\n", groupElement & 65535,groupElement>>16);
+ is2005140FSQ = (groupElement == kPrivatePerFrameSq);
} else { //explicit VR with 16-bit length
if ((d.isLittleEndian) )
lLength = buffer[lPos+2] | (buffer[lPos+3] << 8);
@@ -2956,13 +3234,19 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
else
lLength = buffer[lPos+3] | (buffer[lPos+2] << 8) | (buffer[lPos+1] << 16) | (buffer[lPos] << 24);
lPos += 4; //we have loaded the 32-bit length
+ if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isSQ(groupElement))) { //https://github.com/rordenlab/dcm2niix/issues/144
+ vr[0] = 'S';
+ vr[1] = 'Q';
+ lLength = 8; //Sequence Tag
+ is2005140FSQ = (groupElement == kPrivatePerFrameSq);
+ //if (is2005140FSQ) printMessage(" !!!SQ\t%04x,%04x\n", groupElement & 65535,groupElement>>16);
+ }
} //if explicit else implicit VR
if (lLength == 0xFFFFFFFF) {
lLength = 8; //SQ (Sequences) use 0xFFFFFFFF [4294967295] to denote unknown length
vr[0] = 'S';
vr[1] = 'Q';
}
-
if ((vr[0] == 'S') && (vr[1] == 'Q')) {
sqDepth++;
//printMessage("SQstart %d\n", sqDepth);
@@ -2986,6 +3270,12 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
if (true) { //(nest <= 0) { //some Philips images have different 0020,0013
//verbose reporting :
// printMessage("Pos %ld GroupElement %#08x,%#08x Length %d isLittle %d\n", lPos, (groupElement & 0xFFFF), (groupElement >> 16), lLength, d.isLittleEndian);
+
+ // // Debugging
+ // int groupItem = groupElement >> 16;
+ // int groupGroup = groupElement - (groupItem << 16);
+ // printMessage("groupElement: (%04X, %04X)\n", groupGroup, groupItem);
+
switch ( groupElement ) {
case kTransferSyntax: {
char transferSyntax[kDICOMStr];
@@ -3031,6 +3321,13 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
d.imageStart = 1;//abort as invalid (imageStart MUST be >128)
}
break;} //{} provide scope for variable 'transferSyntax
+ /*case kImplementationVersionName: {
+ char impTxt[kDICOMStr];
+ dcmStr (lLength, &buffer[lPos], impTxt);
+ int slen = (int) strlen(impTxt);
+ if((slen < 6) || (strstr(impTxt, "OSIRIX") == NULL) ) break;
+ printError("OSIRIX Detected\n");
+ break; }*/
case kSourceApplicationEntityTitle: {
char saeTxt[kDICOMStr];
dcmStr (lLength, &buffer[lPos], saeTxt);
@@ -3080,6 +3377,7 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
break;
case kManufacturer:
d.manufacturer = dcmStrManufacturer (lLength, &buffer[lPos]);
+ volDiffusion.manufacturer = d.manufacturer;
break;
case kInstitutionName:
dcmStr(lLength, &buffer[lPos], d.institutionName);
@@ -3163,6 +3461,9 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
case kPatientOrient :
dcmStr (lLength, &buffer[lPos], d.patientOrient);
break;
+ case kDiffusionDirectionality : // 0018, 9075
+ set_directionality0018_9075(&volDiffusion, (&buffer[lPos]));
+ break;
case kDwellTime :
d.dwellTime = dcmStrInt(lLength, &buffer[lPos]);
break;
@@ -3175,27 +3476,27 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
break;*/
case kDiffusionDirectionGEX :
- if (d.manufacturer == kMANUFACTURER_GE) d.CSA.dtiV[1] = dcmStrFloat(lLength, &buffer[lPos]);
+ if (d.manufacturer == kMANUFACTURER_GE)
+ set_diffusion_directionGE(&volDiffusion, lLength, (&buffer[lPos]), 0);
break;
case kDiffusionDirectionGEY :
- if (d.manufacturer == kMANUFACTURER_GE) d.CSA.dtiV[2] = dcmStrFloat(lLength, &buffer[lPos]);
+ if (d.manufacturer == kMANUFACTURER_GE)
+ set_diffusion_directionGE(&volDiffusion, lLength, (&buffer[lPos]), 1);
break;
case kDiffusionDirectionGEZ :
- if (d.manufacturer == kMANUFACTURER_GE) {
- d.CSA.dtiV[3] = dcmStrFloat(lLength, &buffer[lPos]);
- d.CSA.numDti = 1;
- }
+ if (d.manufacturer == kMANUFACTURER_GE)
+ set_diffusion_directionGE(&volDiffusion, lLength, (&buffer[lPos]), 2);
break;
- case kBandwidthPerPixelPhaseEncode:
- d.bandwidthPerPixelPhaseEncode = dcmFloatDouble(lLength, &buffer[lPos],d.isLittleEndian);
- break;
- case kStudyInstanceUID :
+ case kBandwidthPerPixelPhaseEncode:
+ d.bandwidthPerPixelPhaseEncode = dcmFloatDouble(lLength, &buffer[lPos],d.isLittleEndian);
+ break;
+ case kStudyInstanceUID : // 0020, 000D
dcmStr (lLength, &buffer[lPos], d.studyInstanceUID);
break;
- case kSeriesInstanceUID :
+ case kSeriesInstanceUID : // 0020, 000E
dcmStr (lLength, &buffer[lPos], d.seriesInstanceUID);
break;
- case kPatientPosition :
+ case kPatientPosition : // 0020, 0032, ImagePositionPatient
if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (is2005140FSQ)) {
if (!is2005140FSQwarned)
printWarning("Philips R3.2.2 can report different positions for the same slice. Attempting patch.\n");
@@ -3220,6 +3521,7 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
d.patientPositionSequentialRepeats = patientPositionNum-1;
} //if different position from 1st slice in file
} //if not first slice in file
+ set_isAtFirstPatientPosition_tvd(&volDiffusion, isAtFirstPatientPosition);
if (isVerbose > 1)
printMessage(" Patient Position 0020,0032 (#,@,X,Y,Z)\t%d\t%ld\t%g\t%g\t%g\n", patientPositionNum, lPos, patientPosition[1], patientPosition[2], patientPosition[3]);
@@ -3244,6 +3546,20 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
//int dx = 3;
if (d.imageNum <= 1) d.imageNum = dcmStrInt(lLength, &buffer[lPos]); //Philips renames each image as image1 in 2001,9000
break;
+ case kDimensionIndexValues: // kImageNum is not enough for 4D series from Philips 5.*.
+ { // { necessary for initializing ndim.
+ uint8_t ndim = lLength / 4;
+
+ if(ndim > MAX_NUMBER_OF_DIMENSIONS){
+ // Ideally degenerate axes would be cleverly handled.
+ // They are commonly seen in NIfTIs, but perhaps not in DICOM.
+ printError("%d is too many dimensions. Only up to %d are supported\n", ndim,
+ MAX_NUMBER_OF_DIMENSIONS);
+ ndim = MAX_NUMBER_OF_DIMENSIONS; // Truncate
+ }
+ dcmMultiLongs(4 * ndim, &buffer[lPos], ndim, d.dimensionIndexValues, d.isLittleEndian);
+ }
+ break;
case kPhotometricInterpretation:
char interp[kDICOMStr];
dcmStr (lLength, &buffer[lPos], interp);
@@ -3385,6 +3701,13 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
d.accelFactPE = (float)strtof(accelStr, &ptr);
if (*ptr != '\0')
d.accelFactPE = 0.0;
+ //between slice accel
+ dcmStr (lLength, &buffer[lPos], accelStr);
+ dcmStrDigitsOnlyKey('s', accelStr); //e.g. if "p2s4" return "4", if "p2" return ""
+ multiBandFactor = (int)strtol(accelStr, &ptr, 10);
+ if (*ptr != '\0')
+ multiBandFactor = 0.0;
+ //printMessage("p%gs%d\n", d.accelFactPE, multiBandFactor);
break; }
case kLocationsInAcquisition :
d.locationsInAcquisition = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
@@ -3457,23 +3780,70 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
else
d.sliceOrient = kSliceOrientTra; //transverse (axial)
break; }
- case kDiffusionBFactor:
- if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) {
- d.CSA.numDti++; //increment with BFactor: on Philips slices with B=0 have B-factor but no diffusion directions
- if (d.CSA.numDti == 2) { //First time we know that this is a 4D DTI dataset
- //d.dti4D = (TDTI *)malloc(kMaxDTI4D * sizeof(TDTI));
- dti4D->S[0].V[0] = d.CSA.dtiV[0];
- dti4D->S[0].V[1] = d.CSA.dtiV[1];
- dti4D->S[0].V[2] = d.CSA.dtiV[2];
- dti4D->S[0].V[3] = d.CSA.dtiV[3];
- }
- d.CSA.dtiV[0] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);
- if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D))
- dti4D->S[d.CSA.numDti-1].V[0] = d.CSA.dtiV[0];
- /*if ((d.CSA.numDti > 0) && (d.CSA.numDti <= kMaxDTIv))
- d.CSA.dtiV[d.CSA.numDti-1][0] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);*/
- }
- break;
+ // case kDiffusionBFactor: // 2001,1003
+ // if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) {
+ // d.CSA.numDti++; //increment with BFactor: on Philips slices with B=0 have B-factor but no diffusion directions
+ // if (d.CSA.numDti == 2) { //First time we know that this is a 4D DTI dataset
+ // //d.dti4D = (TDTI *)malloc(kMaxDTI4D * sizeof(TDTI));
+ // dti4D->S[0].V[0] = d.CSA.dtiV[0];
+ // dti4D->S[0].V[1] = d.CSA.dtiV[1];
+ // dti4D->S[0].V[2] = d.CSA.dtiV[2];
+ // dti4D->S[0].V[3] = d.CSA.dtiV[3];
+ // }
+ // d.CSA.dtiV[0] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);
+ // if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D))
+ // dti4D->S[d.CSA.numDti-1].V[0] = d.CSA.dtiV[0];
+ // /*if ((d.CSA.numDti > 0) && (d.CSA.numDti <= kMaxDTIv))
+ // d.CSA.dtiV[d.CSA.numDti-1][0] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);*/
+ // }
+ // break;
+ case kDiffusion_bValue: // 0018, 9087
+ // Note that this is ahead of kPatientPosition (0020,0032), so
+ // isAtFirstPatientPosition is not necessarily set yet.
+ // Philips uses this tag too, at least as of 5.1, but they also
+ // use kDiffusionBFactor (see above), and we do not want to
+ // double count. More importantly, with Philips this tag
+ // (sometimes?) gets repeated in a nested sequence with the
+ // value *unset*!
+ // GE started using this tag in 27, and annoyingly, NOT including
+ // the b value if it is 0 for the slice.
+ if((d.manufacturer != kMANUFACTURER_PHILIPS) || !is2005140FSQ){
+ // d.CSA.numDti++;
+
+ // if (d.CSA.numDti == 2) { //First time we know that this is a 4D DTI dataset
+ // //d.dti4D = (TDTI *)malloc(kMaxDTI4D * sizeof(TDTI));
+ // dti4D->S[0].V[0] = d.CSA.dtiV[0];
+ // dti4D->S[0].V[1] = d.CSA.dtiV[1];
+ // dti4D->S[0].V[2] = d.CSA.dtiV[2];
+ // dti4D->S[0].V[3] = d.CSA.dtiV[3];
+ // }
+ //d.CSA.dtiV[0] = dcmFloatDouble(lLength, &buffer[lPos], d.isLittleEndian);
+ set_bVal(&volDiffusion, dcmFloatDouble(lLength, &buffer[lPos], d.isLittleEndian));
+
+ // if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D))
+ // dti4D->S[d.CSA.numDti-1].V[0] = d.CSA.dtiV[0];
+ }
+ break;
+ case kDiffusionOrientation: // 0018, 9089
+ // Note that this is ahead of kPatientPosition (0020,0032), so
+ // isAtFirstPatientPosition is not necessarily set yet.
+ // Philips uses this tag too, at least as of 5.1, but they also
+ // use kDiffusionDirectionRL, etc., and we do not want to double
+ // count. More importantly, with Philips this tag (sometimes?)
+ // gets repeated in a nested sequence with the value *unset*!
+ // if (((d.manufacturer == kMANUFACTURER_SIEMENS) ||
+ // ((d.manufacturer == kMANUFACTURER_PHILIPS) && !is2005140FSQ)) &&
+ // (isAtFirstPatientPosition || isnan(d.patientPosition[1])))
+ if((d.manufacturer == kMANUFACTURER_SIEMENS) ||
+ ((d.manufacturer == kMANUFACTURER_PHILIPS) && !is2005140FSQ))
+ set_orientation0018_9089(&volDiffusion, lLength, &buffer[lPos], d.isLittleEndian);
+ break;
+ // case kSharedFunctionalGroupsSequence:
+ // if ((d.manufacturer == kMANUFACTURER_SIEMENS) && isAtFirstPatientPosition) {
+ // break; // For now - need to figure out how to get the nested
+ // // part of buffer[lPos].
+ // }
+ // break;
case kSliceNumberMrPhilips :
if (d.manufacturer != kMANUFACTURER_PHILIPS)
break;
@@ -3512,35 +3882,35 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
break;
locationsInAcquisitionPhilips = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
break;
- case kDiffusionDirectionRL:
- if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) {
- d.CSA.dtiV[1] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);
- if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D))
- dti4D->S[d.CSA.numDti-1].V[1] = d.CSA.dtiV[1];
- }
- /*if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition) && (d.CSA.numDti > 0) && (d.CSA.numDti <= kMaxDTIv))
- d.CSA.dtiV[d.CSA.numDti-1][1] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);*/
- break;
- case kDiffusionDirectionAP:
- if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) {
- d.CSA.dtiV[2] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);
- if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D))
- dti4D->S[d.CSA.numDti-1].V[2] = d.CSA.dtiV[2];
- }
- /*if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition) && (d.CSA.numDti > 0) && (d.CSA.numDti <= kMaxDTIv))
- d.CSA.dtiV[d.CSA.numDti-1][2] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);*/
- break;
- case kDiffusionDirectionFH:
- if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) {
- d.CSA.dtiV[3] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);
- if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D))
- dti4D->S[d.CSA.numDti-1].V[3] = d.CSA.dtiV[3];
- //printMessage("dti XYZ %g %g %g\n",d.CSA.dtiV[1],d.CSA.dtiV[2],d.CSA.dtiV[3]);
- }
- /*if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition) && (d.CSA.numDti > 0) && (d.CSA.numDti <= kMaxDTIv))
- d.CSA.dtiV[d.CSA.numDti-1][3] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);*/
- //http://www.na-mic.org/Wiki/index.php/NAMIC_Wiki:DTI:DICOM_for_DWI_and_DTI
- break;
+ // case kDiffusionDirectionRL:
+ // if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) {
+ // d.CSA.dtiV[1] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);
+ // if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D))
+ // dti4D->S[d.CSA.numDti-1].V[1] = d.CSA.dtiV[1];
+ // }
+ // /*if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition) && (d.CSA.numDti > 0) && (d.CSA.numDti <= kMaxDTIv))
+ // d.CSA.dtiV[d.CSA.numDti-1][1] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);*/
+ // break;
+ // case kDiffusionDirectionAP:
+ // if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) {
+ // d.CSA.dtiV[2] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);
+ // if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D))
+ // dti4D->S[d.CSA.numDti-1].V[2] = d.CSA.dtiV[2];
+ // }
+ // /*if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition) && (d.CSA.numDti > 0) && (d.CSA.numDti <= kMaxDTIv))
+ // d.CSA.dtiV[d.CSA.numDti-1][2] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);*/
+ // break;
+ // case kDiffusionDirectionFH:
+ // if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition)) {
+ // d.CSA.dtiV[3] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);
+ // if ((d.CSA.numDti > 1) && (d.CSA.numDti < kMaxDTI4D))
+ // dti4D->S[d.CSA.numDti-1].V[3] = d.CSA.dtiV[3];
+ // //printMessage("dti XYZ %g %g %g\n",d.CSA.dtiV[1],d.CSA.dtiV[2],d.CSA.dtiV[3]);
+ // }
+ // /*if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (isAtFirstPatientPosition) && (d.CSA.numDti > 0) && (d.CSA.numDti <= kMaxDTIv))
+ // d.CSA.dtiV[d.CSA.numDti-1][3] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);*/
+ // //http://www.na-mic.org/Wiki/index.php/NAMIC_Wiki:DTI:DICOM_for_DWI_and_DTI
+ // break;
case kWaveformSq:
d.imageStart = 1; //abort!!!
printMessage("Skipping DICOM (audio not image) '%s'\n", fname);
@@ -3570,7 +3940,8 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
if (d.manufacturer == kMANUFACTURER_GE) d.effectiveEchoSpacingGE = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
break;
case kDiffusionBFactorGE :
- if (d.manufacturer == kMANUFACTURER_GE) d.CSA.dtiV[0] = dcmStrInt(lLength, &buffer[lPos]);
+ if (d.manufacturer == kMANUFACTURER_GE)
+ set_bValGE(&volDiffusion, lLength, &buffer[lPos]);
break;
case kGeiisFlag:
if ((lLength > 4) && (buffer[lPos]=='G') && (buffer[lPos+1]=='E') && (buffer[lPos+2]=='I') && (buffer[lPos+3]=='I')) {
@@ -3641,6 +4012,7 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
//printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tnest=%d\t%s\n", groupElement & 65535,groupElement>>16, lLength, lPos, nest, tagStr);
} else
printMessage(" Tag\t%04x,%04x\tSize=%u\tOffset=%ld\tnest=%d\n", groupElement & 65535,groupElement>>16, lLength, lFileOffset+lPos, nest);
+ //if (d.isExplicitVR) printMessage(" VR=%c%c\n", vr[0], vr[1]);
} //printMessage(" tag=%04x,%04x length=%u pos=%ld %c%c nest=%d\n", groupElement & 65535,groupElement>>16, lLength, lPos,vr[0], vr[1], nest);
lPos = lPos + (lLength);
//printMessage("%d\n",d.imageStart);
@@ -3727,8 +4099,8 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
strcpy(d.protocolName, d.seriesDescription);
if ((strlen(d.protocolName) < 1) && (strlen(d.seriesDescription) > 1))
strcpy(d.protocolName, d.seriesDescription);
- if ((strlen(d.protocolName) < 1) && (strlen(d.sequenceName) > 1))
- strcpy(d.protocolName, d.sequenceName);
+ if ((strlen(d.protocolName) < 1) && (strlen(d.sequenceName) > 1) && (d.manufacturer != kMANUFACTURER_SIEMENS))
+ strcpy(d.protocolName, d.sequenceName); //protocolName (0018,1030) optional, sequence name (0018,0024) is not a good substitute for Siemens as it can vary per volume: *ep_b0 *ep_b1000#1, *ep_b1000#2, etc https://www.nitrc.org/forum/forum.php?thread_id=8771&forum_id=4703
// if (!isOrient) {
// if (d.isNonImage)
// printWarning("Spatial orientation ambiguous (tag 0020,0037 not found) [probably not important: derived image]: %s\n", fname);
@@ -3760,6 +4132,8 @@ uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD
printError("Unable to convert DTI [recompile with increased kMaxDTI4D] detected=%d, max = %d\n", d.CSA.numDti, kMaxDTI4D);
d.CSA.numDti = 0;
}
+ if (multiBandFactor > d.CSA.multiBandFactor)
+ d.CSA.multiBandFactor = multiBandFactor; //SMS reported in 0051,1011 but not CSA header
//d.isValid = false; //debug only - will not create output!
#ifndef myLoadWholeFileToReadHeader
fclose(file);
diff --git a/console/nii_dicom.h b/console/nii_dicom.h
index 9b32399..fb513e0 100644
--- a/console/nii_dicom.h
+++ b/console/nii_dicom.h
@@ -38,7 +38,7 @@ extern "C" {
#define kCCsuf " CompilerNA" //unknown compiler!
#endif
-#define kDCMvers "v1.0.20171017" kDCMsuf kCCsuf
+#define kDCMvers "v1.0.20171215" kDCMsuf kCCsuf
static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic
static const int kMaxDTI4D = 4096; //maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images
@@ -70,6 +70,10 @@ static const int kCompressC3 = 2; //obsolete JPEG lossless
static const int kCompress50 = 3; //obsolete JPEG lossy
static const int kCompressRLE = 4; //run length encoding
+// Maximum number of dimensions for .dimensionIndexValues, i.e. possibly the
+// number of axes in the output .nii.
+static const uint8_t MAX_NUMBER_OF_DIMENSIONS = 8;
+
struct TDTI {
float V[4];
int sliceNumberMrPhilips;
@@ -77,6 +81,7 @@ static const int kCompressRLE = 4; //run length encoding
struct TDTI4D {
struct TDTI S[kMaxDTI4D];
};
+
#ifdef _MSC_VER //Microsoft nomenclature for packed structures is different...
#pragma pack(2)
typedef struct {
@@ -121,6 +126,7 @@ static const int kCompressRLE = 4; //run length encoding
//char mrAcquisitionType[kDICOMStr]
char scanOptions[kDICOMStr], stationName[kDICOMStr], softwareVersions[kDICOMStr], deviceSerialNumber[kDICOMStr], institutionAddress[kDICOMStr], institutionName[kDICOMStr], referringPhysicianName[kDICOMStr], seriesInstanceUID[kDICOMStr], studyInstanceUID[kDICOMStr], bodyPartExamined[kDICOMStr], procedureStepDescription[kDICOMStr], imageType[kDICOMStr], institutionalDepartmentName[kDICOMStr], manufacturersModelName[kDICOMStr], patientID[kDICOMStr], patientOrient[kDICOMStr], patient [...]
char imageComments[kDICOMStrLarge];
+ uint32_t dimensionIndexValues[MAX_NUMBER_OF_DIMENSIONS];
struct TCSAdata CSA;
bool isSegamiOasis, isDerived, isXRay, isMultiEcho, isSlicesSpatiallySequentialPhilips, isValid, is3DAcq, is2DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase,isHasMagnitude,isHasMixed, isFloat, isResampled;
char phaseEncodingRC, patientSex;
diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp
index e7ab697..e7df4d5 100644
--- a/console/nii_dicom_batch.cpp
+++ b/console/nii_dicom_batch.cpp
@@ -70,7 +70,8 @@
#endif
struct TDCMsort {
- uint64_t indx, img;
+ uint64_t indx, img;
+ uint32_t dimensionIndexValues[MAX_NUMBER_OF_DIMENSIONS];
};
struct TSearchList {
@@ -237,9 +238,9 @@ void geCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI *vx){
+ (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
- for (int v= 0; v < 4; v++)
- vx[i].V[v] =0.0f;
- continue; //do not normalize or reorient b0 vectors
+ for (int v= 1; v < 4; v++)
+ vx[i].V[v] = 0.0f;
+ continue; //do not normalize or reorient 0 vectors
}
if (!col) { //rows need to be swizzled
//see Stanford dataset Ax_DWI_Tetrahedral_7 unable to resolve between possible solutions
@@ -535,6 +536,7 @@ void siemensCsaAscii(const char * filename, int csaOffset, int csaLength, float
char keyStr[] = "### ASCCONV BEGIN"; //skip to start of ASCII often "### ASCCONV BEGIN ###" but also "### ASCCONV BEGIN object=MrProtDataImpl at MrProtocolData"
char *keyPos = (char *)memmem(bufferTrim, csaLengthTrim, keyStr, strlen(keyStr));
if (keyPos) {
+ //We could detect multi-echo MPRAGE here, e.g. "lContrasts = 4"- but ideally we want an earlier detection
csaLengthTrim -= (keyPos-bufferTrim);
char keyStrEnd[] = "### ASCCONV END";
char *keyPosEnd = (char *)memmem(keyPos, csaLengthTrim, keyStrEnd, strlen(keyStrEnd));
@@ -936,11 +938,28 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts,
// Slice Timing
if (d.CSA.sliceTiming[0] >= 0.0) {
fprintf(fp, "\t\"SliceTiming\": [\n");
- for (int i = 0; i < kMaxEPI3D; i++) {
- if (d.CSA.sliceTiming[i] < 0.0) break;
- if (i != 0)
- fprintf(fp, ",\n");
- fprintf(fp, "\t\t%g", d.CSA.sliceTiming[i] / 1000.0 );
+ if (d.CSA.protocolSliceNumber1 > 1) {
+ //https://github.com/rordenlab/dcm2niix/issues/40
+ //equivalent to dicm2nii "s.SliceTiming = s.SliceTiming(end:-1:1);"
+ int mx = 0;
+ for (int i = 0; i < kMaxEPI3D; i++) {
+ if (d.CSA.sliceTiming[i] < 0.0) break;
+ mx++;
+ }
+ mx--;
+ for (int i = mx; i >= 0; i--) {
+ if (d.CSA.sliceTiming[i] < 0.0) break;
+ if (i != mx)
+ fprintf(fp, ",\n");
+ fprintf(fp, "\t\t%g", d.CSA.sliceTiming[i] / 1000.0 );
+ }
+ } else {
+ for (int i = 0; i < kMaxEPI3D; i++) {
+ if (d.CSA.sliceTiming[i] < 0.0) break;
+ if (i != 0)
+ fprintf(fp, ",\n");
+ fprintf(fp, "\t\t%g", d.CSA.sliceTiming[i] / 1000.0 );
+ }
}
fprintf(fp, "\t],\n");
}
@@ -1112,15 +1131,35 @@ int * nii_SaveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],str
bvals[i] = bvals[i] + (0.5 * i/numDti); //add a small bias so ties are kept in sequential order
}
if (*numADC > 0) {
- if ((numDti - *numADC) < 2) {
- if (!dcmList[indx0].isDerived) //no need to warn if images are derived Trace/ND pair
- printWarning("No bvec/bval files created: only single value after ADC excluded\n");
- *numADC = 0;
- free(bvals);
- free(vx);
- return NULL;
- }
- printMessage("Note: %d volumes appear to be ADC images that will be removed to allow processing\n", *numADC);
+ // DWIs (i.e. short diffusion scans with too few directions to
+ // calculate tensors...they typically acquire b=0 + 3 b > 0 so
+ // the isotropic trace or MD can be calculated) often come as
+ // b=0 and trace pairs, with the b=0 and trace in either order,
+ // and often as "ORIGINAL", even though the trace is not.
+ // The bval file is needed for downstream processing to know
+ // * which is the b=0 and which is the trace, and
+ // * what b is for the trace,
+ // so dcm2niix should *always* write the bval and bvec files,
+ // AND include the b for the trace for DWIs.
+ // One hackish way to accomplish that is to set *numADC = 0
+ // when *numADC == 1 && numDti == 2.
+ // - Rob Reid, 2017-11-29.
+ if ((*numADC == 1) && ((numDti - *numADC) < 2)){
+ *numADC = 0;
+ printMessage("Note: this appears to be a b=0+trace DWI; ADC/trace removal has been disabled.\n");
+ }
+ else{
+ if ((numDti - *numADC) < 2) {
+ if (!dcmList[indx0].isDerived) //no need to warn if images are derived Trace/ND pair
+ printWarning("No bvec/bval files created: only single value after ADC excluded\n");
+ *numADC = 0;
+ free(bvals);
+ free(vx);
+ return NULL;
+ }
+ printMessage("Note: %d volumes appear to be ADC or trace images that will be removed to allow processing\n",
+ *numADC);
+ }
}
//sort ALL including ADC
int * volOrderIndex = (int *) malloc(numDti * sizeof(int));
@@ -2299,7 +2338,7 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis
printMessage("Ignoring derived image(s) of series %ld %s\n", dcmList[indx].seriesNum, nameList->str[indx]);
return EXIT_SUCCESS;
}
- if ((opts.isIgnoreDerivedAnd2D) && (strcmp(dcmList[indx].sequenceName, "_fl2d1")== 0)) {
+ if ((opts.isIgnoreDerivedAnd2D) && ((strcmp(dcmList[indx].sequenceName, "_fl3d1_ns")== 0) || (strcmp(dcmList[indx].sequenceName, "_fl2d1")== 0)) ) {
printMessage("Ignoring localizer (sequence %s) of series %ld %s\n", dcmList[indx].sequenceName, dcmList[indx].seriesNum, nameList->str[indx]);
return EXIT_SUCCESS;
}
@@ -2438,12 +2477,21 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis
free(imgM);
return EXIT_FAILURE;
}
+ // Prevent these DICOM files from being reused.
+ for(int i = 0; i < nConvert; ++i)
+ dcmList[dcmSort[i].indx].converted2NII = 1;
+ if (opts.numSeries < 0) { //report series number but do not convert
+ printMessage("\t%d\t%s\n", dcmList[dcmSort[0].indx].seriesNum, pathoutname);
+ printMessage(" %s\n",nameList->str[dcmSort[0].indx]);
+ return EXIT_SUCCESS;
+ }
checkSliceTiming(&dcmList[indx0], &dcmList[indx1]);
int sliceDir = 0;
if (hdr0.dim[3] > 1)sliceDir = headerDcm2Nii2(dcmList[dcmSort[0].indx],dcmList[dcmSort[nConvert-1].indx] , &hdr0, true);
//UNCOMMENT NEXT TWO LINES TO RE-ORDER MOSAIC WHERE CSA's protocolSliceNumber does not start with 1
if (dcmList[dcmSort[0].indx].CSA.protocolSliceNumber1 > 1) {
- printWarning("Weird CSA 'ProtocolSliceNumber' (%d): SPATIAL, SLICE-ORDER AND DTI TRANSFORMS UNTESTED\n", dcmList[dcmSort[0].indx].CSA.protocolSliceNumber1);
+ printWarning("Weird CSA 'ProtocolSliceNumber' (System/Miscellaneous/ImageNumbering reversed): VALIDATE SLICETIMING AND BVECS\n");
+ //https://www.healthcare.siemens.com/siemens_hwem-hwem_ssxa_websites-context-root/wcm/idc/groups/public/@global/@imaging/@mri/documents/download/mdaz/nzmy/~edisp/mri_60_graessner-01646277.pdf
//see https://github.com/neurolabusc/dcm2niix/issues/40
sliceDir = -1; //not sure how to handle negative determinants?
}
@@ -2451,6 +2499,18 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis
imgM = nii_flipZ(imgM, &hdr0);
sliceDir = abs(sliceDir); //change this, we have flipped the image so GE DTI bvecs no longer need to be flipped!
}
+ // skip converting if user has specified one or more series, but has not specified this one
+ if (opts.numSeries > 0) {
+ int i = 0;
+ for (; i < opts.numSeries; i++) {
+ if (opts.seriesNumber[i] == dcmList[dcmSort[0].indx].seriesNum) {
+ break;
+ }
+ }
+ if (i == opts.numSeries) {
+ return EXIT_SUCCESS;
+ }
+ }
//move before headerDcm2Nii2 checkSliceTiming(&dcmList[indx0], &dcmList[indx1]);
//nii_SaveBIDS(pathoutname, dcmList[dcmSort[0].indx], opts, dti4D, &hdr0, nameList->str[dcmSort[0].indx]);
nii_SaveBIDS(pathoutname, dcmList[dcmSort[0].indx], opts, &hdr0, nameList->str[dcmSort[0].indx]);
@@ -2459,7 +2519,6 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis
free(imgM);
return EXIT_SUCCESS;
}
-
nii_SaveText(pathoutname, dcmList[dcmSort[0].indx], opts, &hdr0, nameList->str[indx]);
int numADC = 0;
int * volOrderIndex = nii_SaveDTI(pathoutname,nConvert, dcmSort, dcmList, opts, sliceDir, dti4D, &numADC);
@@ -2530,27 +2589,56 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis
}
if ((opts.isCrop) && (dcmList[indx0].is3DAcq) && (hdr0.dim[3] > 1) && (hdr0.dim[0] < 4))//for T1 scan: && (dcmList[indx0].TE < 25)
returnCode = nii_saveCrop(pathoutname, hdr0, imgM,opts); //n.b. must be run AFTER nii_setOrtho()!
-
#ifdef HAVE_R
// Note that for R, only one image should be created per series
// Hence the logical OR here
if (returnCode == EXIT_SUCCESS || nii_saveNII(pathoutname,hdr0,imgM,opts) == EXIT_SUCCESS)
nii_saveAttributes(dcmList[dcmSort[0].indx], hdr0, opts);
#endif
-
free(imgM);
return returnCode;//EXIT_SUCCESS;
}// saveDcm2Nii()
+void fillTDCMsort(struct TDCMsort& tdcmref, const uint64_t indx, const struct TDICOMdata& dcmdata){
+ // Copy the relevant parts of dcmdata to tdcmref.
+ tdcmref.indx = indx;
+ tdcmref.img = ((uint64_t)dcmdata.seriesNum << 32) + dcmdata.imageNum;
+ for(int i = 0; i < MAX_NUMBER_OF_DIMENSIONS; ++i)
+ tdcmref.dimensionIndexValues[i] = dcmdata.dimensionIndexValues[i];
+} // fillTDCMsort()
+
int compareTDCMsort(void const *item1, void const *item2) {
//for quicksort http://blog.ablepear.com/2011/11/objective-c-tuesdays-sorting-arrays.html
struct TDCMsort const *dcm1 = (const struct TDCMsort *)item1;
struct TDCMsort const *dcm2 = (const struct TDCMsort *)item2;
+
+ int retval = 0; // tie
+
if (dcm1->img < dcm2->img)
- return -1;
+ retval = -1;
else if (dcm1->img > dcm2->img)
- return 1;
- return 0; //tie
+ retval = 1;
+
+ if(retval == 0){
+ // Check the dimensionIndexValues (useful for enhanced DICOM 4D series).
+ // ->img is basically behaving as a (seriesNum, imageNum) sort key
+ // concatenated into a (large) integer for qsort. That is unwieldy when
+ // dimensionIndexValues need to be compared, because the existence of
+ // uint128_t, uint256_t, etc. is not guaranteed. This sorts by
+ // (seriesNum, ImageNum, div[0], div[1], ...), or if you think of it as a
+ // number, the dimensionIndexValues come after the decimal point.
+ for(int i=0; i < MAX_NUMBER_OF_DIMENSIONS; ++i){
+ if(dcm1->dimensionIndexValues[i] < dcm2->dimensionIndexValues[i]){
+ retval = -1;
+ break;
+ }
+ else if(dcm1->dimensionIndexValues[i] > dcm2->dimensionIndexValues[i]){
+ retval = 1;
+ break;
+ }
+ }
+ }
+ return retval;
} //compareTDCMsort()
int isSameFloatGE (float a, float b) {
@@ -2566,7 +2654,7 @@ int isSameFloatDouble (double a, double b) {
}
struct TWarnings { //generate a warning only once per set
- bool acqNumVaries, bitDepthVaries, dateTimeVaries, echoVaries, coilVaries, nameVaries, orientVaries;
+ bool acqNumVaries, bitDepthVaries, dateTimeVaries, echoVaries, coilVaries, nameVaries, nameEmpty, orientVaries;
};
TWarnings setWarnings() {
@@ -2577,13 +2665,13 @@ TWarnings setWarnings() {
r.echoVaries = false;
r.coilVaries = false;
r.nameVaries = false;
+ r.nameEmpty = false;
r.orientVaries = false;
return r;
}
bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, struct TDCMopts* opts, struct TWarnings* warnings, bool *isMultiEcho) {
//returns true if d1 and d2 should be stacked together as a single output
- *isMultiEcho = false;
if (!d1.isValid) return false;
if (!d2.isValid) return false;
if (d1.modality != d2.modality) return false; //do not stack MR and CT data!
@@ -2607,7 +2695,7 @@ bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, struct TDCMopts* opt
}
if (!isSameFloatDouble(d1.dateTime, d2.dateTime)) { //beware, some vendors incorrectly store Image Time (0008,0033) as Study Time (0008,0030).
if (!warnings->dateTimeVaries)
- printMessage("slices not stacked: Study Data/Time (0008,0020 / 0008,0030) varies %12.12f ~= %12.12f\n", d1.dateTime, d2.dateTime);
+ printMessage("slices not stacked: Study Date/Time (0008,0020 / 0008,0030) varies %12.12f ~= %12.12f\n", d1.dateTime, d2.dateTime);
warnings->dateTimeVaries = true;
return false;
}
@@ -2631,9 +2719,13 @@ bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, struct TDCMopts* opt
warnings->coilVaries = true;
return false;
}
- if ((strcmp(d1.protocolName, d2.protocolName) != 0)) {
- if ((!warnings->nameVaries))
- printMessage("slices not stacked: protocol name varies\n");
+ if ((strlen(d1.protocolName) < 1) && (strlen(d2.protocolName) < 1)) {
+ if (!warnings->nameEmpty)
+ printWarning("Empty protocol name(s) (0018,1030)\n");
+ warnings->nameEmpty = true;
+ } else if ((strcmp(d1.protocolName, d2.protocolName) != 0)) {
+ if (!warnings->nameVaries)
+ printMessage("slices not stacked: protocol name varies '%s' != '%s'\n", d1.protocolName, d2.protocolName);
warnings->nameVaries = true;
return false;
}
@@ -2666,10 +2758,9 @@ int singleDICOM(struct TDCMopts* opts, char *fname) {
strcpy(nameList.str[nameList.numItems],filename);
nameList.numItems++;
struct TDCMsort dcmSort[1];
- dcmSort[0].indx = 0;
- dcmSort[0].img = ((uint64_t)dcmList[0].seriesNum << 32) + dcmList[0].imageNum;
dcmList[0].converted2NII = 1;
dcmList[0] = readDICOMv(nameList.str[0], opts->isVerbose, opts->compressFlag, &dti4D); //ignore compile warning - memory only freed on first of 2 passes
+ fillTDCMsort(dcmSort[0], 0, dcmList[0]);
return saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, &dti4D);
}// singleDICOM()
@@ -2736,11 +2827,13 @@ int removeDuplicates(int nConvert, struct TDCMsort dcmSort[]){
if (nConvert < 2) return nConvert;
int nDuplicates = 0;
for (int i = 1; i < nConvert; i++) {
- if (dcmSort[i].img == dcmSort[i-1].img) {
+ if (compareTDCMsort(&dcmSort[i], &dcmSort[i-1]) == 0) {
nDuplicates ++;
} else {
dcmSort[i-nDuplicates].img = dcmSort[i].img;
dcmSort[i-nDuplicates].indx = dcmSort[i].indx;
+ for(int j = 0; j < MAX_NUMBER_OF_DIMENSIONS; ++j)
+ dcmSort[i - nDuplicates].dimensionIndexValues[j] = dcmSort[i].dimensionIndexValues[j];
}
}
if (nDuplicates > 0)
@@ -2753,18 +2846,20 @@ int removeDuplicatesVerbose(int nConvert, struct TDCMsort dcmSort[], struct TSea
if (nConvert < 2) return nConvert;
int nDuplicates = 0;
for (int i = 1; i < nConvert; i++) {
- if (dcmSort[i].img == dcmSort[i-1].img) {
+ if (compareTDCMsort(&dcmSort[i], &dcmSort[i-1]) == 0) {
printMessage("\t%s\t=\t%s\n",nameList->str[dcmSort[i-1].indx],nameList->str[dcmSort[i].indx]);
nDuplicates ++;
} else {
dcmSort[i-nDuplicates].img = dcmSort[i].img;
dcmSort[i-nDuplicates].indx = dcmSort[i].indx;
+ for(int j = 0; j < MAX_NUMBER_OF_DIMENSIONS; ++j)
+ dcmSort[i - nDuplicates].dimensionIndexValues[j] = dcmSort[i].dimensionIndexValues[j];
}
}
if (nDuplicates > 0)
printMessage("Some images have identical time, series, acquisition and image values. Duplicates removed.\n");
return nConvert - nDuplicates;
-}// removeDuplicates()
+}// removeDuplicatesVerbose()
bool isExt (char *file_name, const char* ext) {
char *p_extension;
@@ -2777,6 +2872,7 @@ bool isExt (char *file_name, const char* ext) {
int convert_parRec(struct TDCMopts opts) {
//sample dataset from Ed Gronenschild <ed.gronenschild at maastrichtuniversity.nl>
struct TSearchList nameList;
+ int ret = EXIT_FAILURE;
nameList.numItems = 1;
nameList.maxItems = 1;
nameList.str = (char **) malloc((nameList.maxItems+1) * sizeof(char *)); //we reserve one pointer (32 or 64 bits) per potential file
@@ -2787,7 +2883,8 @@ int convert_parRec(struct TDCMopts opts) {
dcmList[0] = nii_readParRec(nameList.str[0], opts.isVerbose, &dti4D);
struct TDCMsort dcmSort[1];
dcmSort[0].indx = 0;
- int ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, opts, &dti4D);
+ if (dcmList[0].isValid)
+ ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, opts, &dti4D);
free(dcmList);//if (nConvertTotal == 0)
if (nameList.numItems < 1)
printMessage("No valid PAR/REC files were found\n");
@@ -2891,8 +2988,8 @@ int nii_loadDir(struct TDCMopts* opts) {
dcmList[i] = readDICOMv(nameList.str[i], opts->isVerbose, opts->compressFlag, &dti4D); //ignore compile warning - memory only freed on first of 2 passes
if ((dcmList[i].isValid) &&((dcmList[i].patientPositionNumPhilips > 1) || (dcmList[i].CSA.numDti > 1))) { //4D dataset: dti4D arrays require huge amounts of RAM - write this immediately
struct TDCMsort dcmSort[1];
- dcmSort[0].indx = i;
- dcmSort[0].img = ((uint64_t)dcmList[i].seriesNum << 32) + dcmList[i].imageNum;
+
+ fillTDCMsort(dcmSort[0], i, dcmList[i]);
dcmList[i].converted2NII = 1;
int ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, &dti4D);
if (ret == EXIT_SUCCESS)
@@ -2942,24 +3039,24 @@ int nii_loadDir(struct TDCMopts* opts) {
if ((dcmList[i].converted2NII == 0) && (dcmList[i].isValid)) {
int nConvert = 0;
struct TWarnings warnings = setWarnings();
- bool isMultiEcho;
+ bool isMultiEcho = false;
for (int j = i; j < (int)nDcm; j++)
if (isSameSet(dcmList[i], dcmList[j], opts, &warnings, &isMultiEcho ) )
nConvert++;
if (nConvert < 1) nConvert = 1; //prevents compiler warning for next line: never executed since j=i always causes nConvert ++
TDCMsort * dcmSort = (TDCMsort *)malloc(nConvert * sizeof(TDCMsort));
nConvert = 0;
+ isMultiEcho = false;
for (int j = i; j < (int)nDcm; j++)
if (isSameSet(dcmList[i], dcmList[j], opts, &warnings, &isMultiEcho)) {
- dcmSort[nConvert].indx = j;
- dcmSort[nConvert].img = ((uint64_t)dcmList[j].seriesNum << 32) + dcmList[j].imageNum;
- dcmList[j].converted2NII = 1;
+ fillTDCMsort(dcmSort[nConvert], j, dcmList[j]);
nConvert++;
} else {
dcmList[i].isMultiEcho = isMultiEcho;
dcmList[j].isMultiEcho = isMultiEcho;
}
qsort(dcmSort, nConvert, sizeof(struct TDCMsort), compareTDCMsort); //sort based on series and image numbers....
+ //dcmList[dcmSort[0].indx].isMultiEcho = isMultiEcho;
if (opts->isVerbose)
nConvert = removeDuplicatesVerbose(nConvert, dcmSort, &nameList);
else
@@ -3159,6 +3256,8 @@ void setDefaultOpts (struct TDCMopts *opts, const char * argv[]) { //either "set
opts->isVerbose = false;
#endif
opts->isTiltCorrect = true;
+ opts->numSeries = 0;
+ memset(opts->seriesNumber, 0, sizeof(opts->seriesNumber));
strcpy(opts->filename,"%f_%p_%t_%s");
} // setDefaultOpts()
diff --git a/console/nii_dicom_batch.h b/console/nii_dicom_batch.h
index f9f0ba9..ced863e 100644
--- a/console/nii_dicom_batch.h
+++ b/console/nii_dicom_batch.h
@@ -22,10 +22,13 @@ extern "C" {
};
#endif
+#define MAX_NUM_SERIES 16
+
struct TDCMopts {
bool isSave3D,isGz, isFlipY, isCreateBIDS, isSortDTIbyBVal, isAnonymizeBIDS, isOnlyBIDS, isCreateText, isIgnoreDerivedAnd2D, isPhilipsFloatNotDisplayScaling, isTiltCorrect, isRGBplanar, isOnlySingleFile, isForceStackSameSeries, isCrop;
int isVerbose, compressFlag, gzLevel; //support for compressed data 0=none,
char filename[512], outdir[512], indir[512], pigzname[512], optsname[512], indirParent[512], imageComments[24];
+ long seriesNumber[MAX_NUM_SERIES], numSeries;
#ifdef HAVE_R
bool isScanOnly;
void *imageList;
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/dcm2niix.git
More information about the debian-med-commit
mailing list