[med-svn] [dcm2niix] 01/04: New upstream version 1.0.20170923
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Mon Oct 9 17:49:43 UTC 2017
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 8a48fba94c1d3e816bccee519eae9a46c59273af
Author: Ghislain Antony Vaillant <ghisvail at gmail.com>
Date: Mon Oct 9 18:40:16 2017 +0100
New upstream version 1.0.20170923
---
README.md | 11 +-
console/CMakeLists.txt | 8 +-
console/jpg_0XC3.cpp | 14 +-
console/main_console.cpp | 31 +-
console/nifti1_io_core.cpp | 6 +-
console/nii_dicom.cpp | 154 ++++---
console/nii_dicom.h | 14 +-
console/nii_dicom_batch.cpp | 948 ++++++++++++++++++++++++++++++--------------
console/nii_dicom_batch.h | 7 +-
console/nii_foreign.cpp | 5 +-
console/nii_ortho.cpp | 8 +-
console/ucm.cmake | 10 +-
old_nii_SaveBIDS.cpp | 274 +++++++++++++
13 files changed, 1095 insertions(+), 395 deletions(-)
diff --git a/README.md b/README.md
index fdbcb8d..3411e16 100644
--- a/README.md
+++ b/README.md
@@ -15,6 +15,12 @@ This software should run on macOS, Linux and Windows without requiring any other
## Versions
+23-Sept-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).
+ - BIDS tag [TotalReadoutTime](https://github.com/rordenlab/dcm2niix/issues/130) handles partial fourier, Phase Resolution, etc (Michael Harms).
+ - Additional [json fields](https://github.com/rordenlab/dcm2niix/issues/127).
+
18-Aug-2017
- Better BVec extraction for [PAR/REC 4.1](https://www.nitrc.org/forum/forum.php?thread_id=8387&forum_id=4703).
- Support for [Segami Cerescan volumes](https://www.nitrc.org/forum/forum.php?thread_id=8076&forum_id=4703).
@@ -194,9 +200,12 @@ This requires a compiler that supports c++11.
## Links
- - [Dcm2Bids](https://github.com/cbedetti/Dcm2Bids) uses dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets.
+ - [Dcm2Bids](https://github.com/cbedetti/Dcm2Bids) uses dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets.
+ - [bidskit](https://github.com/jmtyszka/bidskit) uses dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets.
+ - [DAC2BIDS](https://github.com/dangom/dac2bids) uses dcm2niibatch to create [BIDS](http://bids.neuroimaging.io/) datasets.
- [heudiconv](https://github.com/nipy/heudiconv) can use dcm2niix to create [BIDS](http://bids.neuroimaging.io/) datasets.
- [nipype](https://github.com/nipy/nipype) can use dcm2niix to convert images.
+ - [pydcm2niix is a Python module for working with dcm2niix](https://github.com/jstutters/pydcm2niix).
- [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.
diff --git a/console/CMakeLists.txt b/console/CMakeLists.txt
index dcd78d3..01dea3c 100644
--- a/console/CMakeLists.txt
+++ b/console/CMakeLists.txt
@@ -25,9 +25,11 @@ if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -Wl,-dead_strip")
elseif("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
# using GCC
- add_definitions(-Wno-unused-result)
- if (${CMAKE_CXX_COMPILER_VERSION} VERSION_GREATER 4.7)
- add_definitions(-fno-diagnostics-show-caret)
+ if (${CMAKE_CXX_COMPILER_VERSION} VERSION_GREATER 4.4.7)
+ add_definitions(-Wno-unused-result) # available since GCC 4.5
+ endif()
+ if (${CMAKE_CXX_COMPILER_VERSION} VERSION_GREATER 4.7.4)
+ add_definitions(-fno-diagnostics-show-caret) # available since GCC 4.8
endif()
elseif("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
# using Intel C++
diff --git a/console/jpg_0XC3.cpp b/console/jpg_0XC3.cpp
index c8ca375..e45906a 100644
--- a/console/jpg_0XC3.cpp
+++ b/console/jpg_0XC3.cpp
@@ -125,7 +125,7 @@ unsigned char * decode_JPEG_SOF_0XC3 (const char *fn, int skipBytes, bool verbo
unsigned char *lRawRA = (unsigned char*) malloc(lRawSz);
size_t lSz = fread(lRawRA, 1, lRawSz, reader);
fclose(reader);
- if ((lSz < lRawSz) || (lRawRA[0] != 0xFF) || (lRawRA[1] != 0xD8) || (lRawRA[2] != 0xFF)) {
+ if ((lSz < (size_t)lRawSz) || (lRawRA[0] != 0xFF) || (lRawRA[1] != 0xD8) || (lRawRA[2] != 0xFF)) {
abortGoto("JPEG signature 0xFFD8FF not found at offset %d of %s\n", skipBytes, fn);//signature failure http://en.wikipedia.org/wiki/List_of_file_signatures
}
if (verbose)
@@ -133,9 +133,13 @@ unsigned char * decode_JPEG_SOF_0XC3 (const char *fn, int skipBytes, bool verbo
//next: read header
long lRawPos = 2; //Skip initial 0xFFD8, begin with third byte
//long lRawPos = 0; //Skip initial 0xFFD8, begin with third byte
- unsigned char btS1, btS2, SOSss, SOSse, SOSahal, SOSpttrans, btMarkerType, SOSns = 0x00; //tag
- uint8_t SOFnf, SOFprecision;
- uint16_t SOFydim, SOFxdim; //, lRestartSegmentSz;
+ unsigned char btS1, btS2, SOSse, SOSahal, btMarkerType, SOSns = 0x00; //tag
+ unsigned char SOSpttrans = 0;
+ unsigned char SOSss = 0;
+ uint8_t SOFnf = 0;
+ uint8_t SOFprecision = 0;
+ uint16_t SOFydim = 0;
+ uint16_t SOFxdim = 0;
// long SOSarrayPos; //SOFarrayPos
int lnHufTables = 0;
const int kmaxFrames = 4;
@@ -186,7 +190,7 @@ unsigned char * decode_JPEG_SOF_0XC3 (const char *fn, int skipBytes, bool verbo
DHTnLi = DHTnLi + l[lFrameCount].DHTliRA[lInc];
if (l[lFrameCount].DHTliRA[lInc] != 0) l[lFrameCount].MaxHufSi = lInc;
if (verbose) printMessage("DHT has %d combinations with %d bits\n", l[lFrameCount].DHTliRA[lInc], lInc);
-
+
}
if (DHTnLi > 17) {
abortGoto("Huffman table corrupted.\n");
diff --git a/console/main_console.cpp b/console/main_console.cpp
index 4c4d55b..3608038 100644
--- a/console/main_console.cpp
+++ b/console/main_console.cpp
@@ -77,6 +77,9 @@ void showHelp(const char * argv[], struct TDCMopts opts) {
printf(" -b : BIDS sidecar (y/n/o(o=only: no NIfTI), default %c)\n", bidsCh);
if (opts.isAnonymizeBIDS) bidsCh = 'y'; else bidsCh = 'n';
printf(" -ba : anonymize BIDS (y/n, default %c)\n", bidsCh);
+ printf(" -c : comment stored as NIfTI aux_file (up to 24 characters)\n");
+ if (opts.isSortDTIbyBVal) bidsCh = 'y'; else bidsCh = 'n';
+ printf(" -d : diffusion volumes sorted by b-value (y/n, default %c)\n", bidsCh);
#ifdef mySegmentByAcq
#define kQstr " %%q=sequence number,"
#else
@@ -96,14 +99,14 @@ void showHelp(const char * argv[], struct TDCMopts opts) {
if (opts.isGz) gzCh = 'y';
#ifdef myDisableZLib
if (strlen(opts.pigzname) > 0)
- printf(" -z : gz compress images (y/n, default %c)\n", gzCh);
+ printf(" -z : gz compress images (y/n/3, default %c) [y=pigz, n=no, 3=no,3D]\n", gzCh);
else
- printf(" -z : gz compress images (y/n, default %c) [REQUIRES pigz]\n", gzCh);
+ printf(" -z : gz compress images (y/n/3, default %c) [y=pigz(MISSING!), n=no, 3=no,3D]\n", gzCh);
#else
#ifdef myDisableMiniZ
- printf(" -z : gz compress images (y/i/n, default %c) [y=pigz, i=internal:zlib, n=no]\n", gzCh);
+ printf(" -z : gz compress images (y/i/n/3, default %c) [y=pigz, i=internal:zlib, n=no, 3=no,3D]\n", gzCh);
#else
- printf(" -z : gz compress images (y/i/n, default %c) [y=pigz, i=internal, n=no]\n", gzCh);
+ printf(" -z : gz compress images (y/i/n/3, default %c) [y=pigz, i=internal, n=no, 3=no,3D]\n", gzCh);
#endif
#endif
@@ -111,6 +114,7 @@ void showHelp(const char * argv[], struct TDCMopts opts) {
printf(" Defaults stored in Windows registry\n");
printf(" Examples :\n");
printf(" %s c:\\DICOM\\dir\n", cstr);
+ printf(" %s -c \"my comment\" c:\\DICOM\\dir\n", cstr);
printf(" %s -o c:\\out\\dir c:\\DICOM\\dir\n", cstr);
printf(" %s -f mystudy%%s c:\\DICOM\\dir\n", cstr);
printf(" %s -o \"c:\\dir with spaces\\dir\" c:\\dicomdir\n", cstr);
@@ -118,6 +122,7 @@ void showHelp(const char * argv[], struct TDCMopts opts) {
printf(" Defaults file : %s\n", opts.optsname);
printf(" Examples :\n");
printf(" %s /Users/chris/dir\n", cstr);
+ printf(" %s -c \"my comment\" /Users/chris/dir\n", cstr);
printf(" %s -o /users/cr/outdir/ -z y ~/dicomdir\n", cstr);
printf(" %s -f %%p_%%s -b y -ba n ~/dicomdir\n", cstr);
printf(" %s -f mystudy%%s ~/dicomdir\n", cstr);
@@ -131,7 +136,8 @@ int invalidParam(int i, const char * argv[]) {
|| (argv[i][0] == 'o') || (argv[i][0] == 'O')
|| (argv[i][0] == 'h') || (argv[i][0] == 'H')
|| (argv[i][0] == 'i') || (argv[i][0] == 'I')
- || (argv[i][0] == '0') || (argv[i][0] == '1') || (argv[i][0] == '2'))
+ || (argv[i][0] == '0') || (argv[i][0] == '1')
+ || (argv[i][0] == '2') || (argv[i][0] == '3') )
return 0;
//if (argv[i][0] != '-') return 0;
@@ -192,6 +198,16 @@ int main(int argc, const char * argv[])
opts.isAnonymizeBIDS = true;
} else
printf("Error: Unknown command line argument: '%s'\n", argv[i]);
+ } else if ((argv[i][1] == 'c') && ((i+1) < argc)) {
+ i++;
+ snprintf(opts.imageComments,24,"%s",argv[i]);
+ } else if ((argv[i][1] == 'd') && ((i+1) < argc)) {
+ i++;
+ if (invalidParam(i, argv)) return 0;
+ if ((argv[i][0] == 'n') || (argv[i][0] == 'N') || (argv[i][0] == '0'))
+ opts.isSortDTIbyBVal = false;
+ else
+ opts.isSortDTIbyBVal = true;
} else if ((argv[i][1] == 'i') && ((i+1) < argc)) {
i++;
if (invalidParam(i, argv)) return 0;
@@ -247,7 +263,10 @@ int main(int argc, const char * argv[])
} else if ((argv[i][1] == 'z') && ((i+1) < argc)) {
i++;
if (invalidParam(i, argv)) return 0;
- if ((argv[i][0] == 'i') || (argv[i][0] == 'I') ) {
+ if ((argv[i][0] == '3') ) {
+ opts.isGz = false; //uncompressed 3D
+ opts.isSave3D = true;
+ } else if ((argv[i][0] == 'i') || (argv[i][0] == 'I') ) {
opts.isGz = true; //force use of internal compression instead of pigz
strcpy(opts.pigzname,"");
} else if ((argv[i][0] == 'n') || (argv[i][0] == 'N') || (argv[i][0] == '0'))
diff --git a/console/nifti1_io_core.cpp b/console/nifti1_io_core.cpp
index 34efd04..45f34f3 100644
--- a/console/nifti1_io_core.cpp
+++ b/console/nifti1_io_core.cpp
@@ -86,19 +86,19 @@ int isSameFloat (float a, float b) {
ivec3 setiVec3(int x, int y, int z)
{
- ivec3 v = {x, y, z};
+ ivec3 v = {{x, y, z}};
return v;
}
vec3 setVec3(float x, float y, float z)
{
- vec3 v = {x, y, z};
+ vec3 v = {{x, y, z}};
return v;
}
vec4 setVec4(float x, float y, float z)
{
- vec4 v= {x, y, z, 1};
+ vec4 v= {{x, y, z, 1}};
return v;
}
diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp
index 92fdd72..af3fbc2 100644
--- a/console/nii_dicom.cpp
+++ b/console/nii_dicom.cpp
@@ -35,8 +35,8 @@
#include "nifti1_io_core.h"
#ifdef HAVE_R
-#undef isnan
-#define isnan ISNAN
+ #undef isnan
+ #define isnan ISNAN
#endif
#ifndef myDisableClassicJPEG
@@ -53,7 +53,7 @@
#include "openjpeg.h"
#ifdef myEnableJasper
-ERROR: YOU CAN NOT COMPILE WITH myEnableJasper AND NOT myDisableOpenJPEG OPTIONS SET SIMULTANEOUSLY
+ ERROR: YOU CAN NOT COMPILE WITH myEnableJasper AND NOT myDisableOpenJPEG OPTIONS SET SIMULTANEOUSLY
#endif
unsigned char * imagetoimg(opj_image_t * image)
@@ -335,7 +335,7 @@ int verify_slice_dir (struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1_
flip = ((sliceV.v[0]+sliceV.v[1]+sliceV.v[2]) < 0);
//printMessage("verify slice dir %g %g %g\n",sliceV.v[0],sliceV.v[1],sliceV.v[2]);
if (isVerbose) { //1st pass only
- if (!d.isNonImage) //do not warn user if image is derived
+ if (!d.isDerived) //do not warn user if image is derived
printWarning("Unable to determine slice direction: please check whether slices are flipped\n");
else
printWarning("Unable to determine slice direction: please check whether slices are flipped (derived image)\n");
@@ -345,16 +345,6 @@ int verify_slice_dir (struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1_
for (int i = 0; i < 4; i++)
R->m[i][2] = -R->m[i][2];
}
- /*if (d.manufacturer == kMANUFACTURER_SEGAMI) { //set origin as center of volume
- for (int i = 0; i < 3; i++)
- R->m[i][3] = 0; //remove old offset
- vec4 originVx = setVec4( (h->dim[1]+1.0f)/2.0f, (h->dim[2]+1.0f)/2.0f, (h->dim[3]+1.0f)/2.0f);
- printMessage("origin (vx) %g %g %g\n",originVx.v[0],originVx.v[1],originVx.v[2]);
- vec4 originMm = nifti_vect44mat44_mul(originVx, *R);
- printMessage("origin (mm) %g %g %g\n",originMm.v[0],originMm.v[1],originMm.v[2]);
- for (int i = 0; i < 3; i++)
- R->m[i][3] = originMm.v[i]; //remove old offset
- }*/
if (flip)
iSL = -iSL;
#ifdef MY_DEBUG
@@ -539,7 +529,7 @@ mat44 set_nii_header(struct TDICOMdata d) {
mat44 set_nii_header_x(struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1_header *h, int* sliceDir, int isVerbose) {
*sliceDir = 0;
mat44 Q44 = nifti_dicom2mat(d.orient, d.patientPosition, d.xyzMM);
- if (d.manufacturer == kMANUFACTURER_SEGAMI) {
+ if (d.isSegamiOasis == true) {
//Segami reconstructions appear to disregard DICOM spatial parameters: assume center of volume is isocenter and no table tilt
// Consider sample image with d.orient (0020,0037) = -1 0 0; 0 1 0: this suggests image RAI (L->R, P->A, S->I) but the vendors viewing software suggests LPS
//Perhaps we should ignore 0020,0037 and 0020,0032 as they are hidden in sequence 0054,0022, but in this case no positioning is provided
@@ -608,7 +598,7 @@ int headerDcm2NiiSForm(struct TDICOMdata d, struct TDICOMdata d2, struct nifti_
//we will have to guess, assume axial acquisition saved in standard Siemens style?
d.orient[1] = 1.0f; d.orient[2] = 0.0f; d.orient[3] = 0.0f;
d.orient[1] = 0.0f; d.orient[2] = 1.0f; d.orient[3] = 0.0f;
- if ((d.isNonImage) || ((d.bitsAllocated == 8) && (d.samplesPerPixel == 3) && (d.manufacturer == kMANUFACTURER_SIEMENS))) {
+ if ((d.isDerived) || ((d.bitsAllocated == 8) && (d.samplesPerPixel == 3) && (d.manufacturer == kMANUFACTURER_SIEMENS))) {
printMessage("Unable to determine spatial orientation: 0020,0037 missing (probably not a problem: derived image)\n");
} else {
printMessage("Unable to determine spatial orientation: 0020,0037 missing!\n");
@@ -629,16 +619,6 @@ int headerDcm2Nii2(struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1_hea
sprintf(dtxt, ";phase=%d", d.CSA.phaseEncodingDirectionPositive);
strcat(txt,dtxt);
}
- if ((d.CSA.bandwidthPerPixelPhaseEncode > 0) && ((d.phaseEncodingRC =='C') || (d.phaseEncodingRC =='R'))) {
- float dwellTime = 0;
- if (d.phaseEncodingRC =='C')
- dwellTime = 1000/d.CSA.bandwidthPerPixelPhaseEncode/h->dim[2];
- else
- dwellTime = 1000/d.CSA.bandwidthPerPixelPhaseEncode/h->dim[1];
- char dtxt[1024] = {""};
- sprintf(dtxt, ";dwell=%.3f", dwellTime);
- strcat(txt,dtxt);
- }
//from dicm2nii 20151117 InPlanePhaseEncodingDirection
if (d.phaseEncodingRC =='R')
h->dim_info = (3 << 4) + (1 << 2) + 2;
@@ -696,6 +676,8 @@ struct TDICOMdata clear_dicom_data() {
strcpy(d.institutionAddress, "");
strcpy(d.deviceSerialNumber, "");
strcpy(d.softwareVersions, "");
+ strcpy(d.stationName, "");
+ strcpy(d.scanOptions, "");
strcpy(d.seriesInstanceUID, "");
strcpy(d.studyID, "");
strcpy(d.studyInstanceUID, "");
@@ -720,6 +702,8 @@ struct TDICOMdata clear_dicom_data() {
d.numberOfDynamicScans = 0;
d.echoNum = 1;
d.echoTrainLength = 0;
+ d.phaseFieldofView = 0.0;
+ d.dwellTime = 0;
d.phaseEncodingSteps = 0;
d.coilNum = 1;
d.accelFactPE = 0.0;
@@ -741,7 +725,8 @@ struct TDICOMdata clear_dicom_data() {
d.imageStart = 0;
d.is3DAcq = false; //e.g. MP-RAGE, SPACE, TFE
d.isSlicesSpatiallySequentialPhilips = true; //Philips can save slices in random order, e.g. 4,5,6,1,2,3
- d.isNonImage = false; //0008,0008 = DERIVED,CSAPARALLEL,POSDISP
+ d.isDerived = false; //0008,0008 = DERIVED,CSAPARALLEL,POSDISP
+ d.isSegamiOasis = false; //these images do not store spatial coordinates
d.bitsAllocated = 16;//bits
d.bitsStored = 0;
d.samplesPerPixel = 1;
@@ -773,6 +758,21 @@ struct TDICOMdata clear_dicom_data() {
return d;
} //clear_dicom_data()
+void dcmStrDigitsOnlyKey(char key, char* lStr) {
+ //e.g. string "p2s3" returns 2 if key=="p" and 3 if key=="s"
+ size_t len = strlen(lStr);
+ if (len < 1) return;
+ bool isKey = false;
+ for (int i = 0; i < (int) len; i++) {
+ if (!isdigit(lStr[i]) ) {
+ isKey = (lStr[i] == key);
+ lStr[i] = ' ';
+
+ } else if (!isKey)
+ lStr[i] = ' ';
+ }
+} //dcmStrDigitsOnlyKey()
+
void dcmStrDigitsOnly(char* lStr) {
//e.g. change "H11" to " 11"
size_t len = strlen(lStr);
@@ -938,8 +938,6 @@ int dcmStrManufacturer (int lByteLength, unsigned char lBuffer[]) {//read float
ret = kMANUFACTURER_PHILIPS;
if ((toupper(cString[0])== 'T') && (toupper(cString[1])== 'O'))
ret = kMANUFACTURER_TOSHIBA;
- if ((toupper(cString[0])== 'S') && (toupper(cString[1])== 'E'))
- ret = kMANUFACTURER_SEGAMI;
//#ifdef _MSC_VER
free(cString);
//#endif
@@ -1007,7 +1005,8 @@ bool csaIsPhaseMap (unsigned char buff[], int nItems) {
return false;
} //csaIsPhaseMap()
-int readCSAImageHeader(unsigned char *buff, int lLength, struct TCSAdata *CSA, int isVerbose, struct TDTI4D *dti4D) {
+//int readCSAImageHeader(unsigned char *buff, int lLength, struct TCSAdata *CSA, int isVerbose, struct TDTI4D *dti4D) {
+int readCSAImageHeader(unsigned char *buff, int lLength, struct TCSAdata *CSA, int isVerbose) {
//see also http://afni.nimh.nih.gov/pub/dist/src/siemens_dicom_csa.c
//printMessage("%c%c%c%c\n",buff[0],buff[1],buff[2],buff[3]);
if (lLength < 36) return EXIT_FAILURE;
@@ -1687,7 +1686,8 @@ size_t nii_ImgBytes(struct nifti_1_header hdr) {
return imgsz;
} //nii_ImgBytes()
-unsigned char * nii_demosaic(unsigned char* inImg, struct nifti_1_header *hdr, int nMosaicSlices, int ProtocolSliceNumber1) {
+//unsigned char * nii_demosaic(unsigned char* inImg, struct nifti_1_header *hdr, int nMosaicSlices, int ProtocolSliceNumber1) {
+unsigned char * nii_demosaic(unsigned char* inImg, struct nifti_1_header *hdr, int nMosaicSlices) {
//demosaic http://nipy.org/nibabel/dicom/dicom_mosaic.html
if (nMosaicSlices < 2) return inImg;
//Byte inImg[ [img length] ];
@@ -2075,8 +2075,7 @@ unsigned char * nii_XYTZ_XYZT(unsigned char* bImg, struct nifti_1_header *hdr, i
//memcpy(&tempImg[0], &bImg[0], imgSz);
uint64_t origPos = 0;
uint64_t Pos = 0; //
-
- for (int t = 0; t < typeRepeats; t++) { //for each volume
+ for (int t = 0; t < (int)typeRepeats; t++) { //for each volume
for (int s = 0; s < seqRepeats; s++) {
origPos = (t*typeBytes) +s*sliceBytes;
for (int z = 0; z < hdr->dim[3]; z++) { //for each slice
@@ -2256,7 +2255,7 @@ TJPEG * decode_JPEG_SOF_0XC3_stack (const char *fn, int skipBytes, bool isVerbo
unsigned char *lRawRA = (unsigned char*) malloc(lRawSz);
size_t lSz = fread(lRawRA, 1, lRawSz, reader);
fclose(reader);
- if (lSz < lRawSz) {
+ if (lSz < (size_t)lRawSz) {
printError("Unable to read %s\n", fn);
abortGoto(); //read failure
}
@@ -2346,7 +2345,8 @@ unsigned char * nii_loadImgJPEGC3(char* imgname, struct nifti_1_header hdr, stru
#ifdef myTurboJPEG //if turboJPEG instead of nanoJPEG for classic JPEG decompression
-unsigned char * nii_loadImgJPEG50(char* imgname, struct nifti_1_header hdr, struct TDICOMdata dcm) {
+//unsigned char * nii_loadImgJPEG50(char* imgname, struct nifti_1_header hdr, struct TDICOMdata dcm) {
+unsigned char * nii_loadImgJPEG50(char* imgname, struct TDICOMdata dcm) {
//decode classic JPEG using nanoJPEG
//printMessage("50 offset %d\n", dcm.imageStart);
if ((dcm.samplesPerPixel != 1) && (dcm.samplesPerPixel != 3)) {
@@ -2396,7 +2396,8 @@ unsigned char * nii_loadImgJPEG50(char* imgname, struct nifti_1_header hdr, stru
#else //if turboJPEG else use nanojpeg...
-unsigned char * nii_loadImgJPEG50(char* imgname, struct nifti_1_header hdr, struct TDICOMdata dcm) {
+//unsigned char * nii_loadImgJPEG50(char* imgname, struct nifti_1_header hdr, struct TDICOMdata dcm) {
+unsigned char * nii_loadImgJPEG50(char* imgname, struct TDICOMdata dcm) {
//decode classic JPEG using nanoJPEG
//printMessage("50 offset %d\n", dcm.imageStart);
if( access(imgname, F_OK ) == -1 ) {
@@ -2459,7 +2460,7 @@ unsigned char * nii_loadImgRLE(char* imgname, struct nifti_1_header hdr, struct
}
fseek(file, 0, SEEK_END);
long fileLen=ftell(file);
- if (fileLen < (dcm.imageBytes+dcm.imageStart)) {
+ if ((fileLen < 1) || (dcm.imageBytes < 1) || (fileLen < (dcm.imageBytes+dcm.imageStart))) {
printMessage("File not large enough to store image data: %s\n", imgname);
fclose(file);
return NULL;
@@ -2469,7 +2470,7 @@ unsigned char * nii_loadImgRLE(char* imgname, struct nifti_1_header hdr, struct
unsigned char *cImg = (unsigned char *)malloc(dcm.imageBytes); //compressed input
size_t sz = fread(cImg, 1, dcm.imageBytes, file);
fclose(file);
- if (sz < dcm.imageBytes) {
+ if (sz < (size_t)dcm.imageBytes) {
printError("Only loaded %zu of %d bytes for %s\n", sz, dcm.imageBytes, imgname);
free(cImg);
return NULL;
@@ -2478,17 +2479,17 @@ unsigned char * nii_loadImgRLE(char* imgname, struct nifti_1_header hdr, struct
bool swap = (dcm.isLittleEndian != littleEndianPlatform());
int bytesPerSample = dcm.samplesPerPixel * (dcm.bitsAllocated / 8);
uint32_t bytesPerSampleRLE = rleInt(0, cImg, swap);
- if (bytesPerSampleRLE != (bytesPerSample)) {
+ if ((bytesPerSample < 0) || (bytesPerSampleRLE != (uint32_t)bytesPerSample)) {
printError("RLE header corrupted %d != %d\n", bytesPerSampleRLE, bytesPerSample);
free(cImg);
return NULL;
}
unsigned char *bImg = (unsigned char *)malloc(imgsz); //binary output
- for (int i = 0; i < imgsz; i++)
+ for (size_t i = 0; i < imgsz; i++)
bImg[i] = 0;
for (int i = 0; i < bytesPerSample; i++) {
uint32_t offset = rleInt(i+1, cImg, swap);
- if (offset > dcm.imageBytes) {
+ if ((dcm.imageBytes < 0) || (offset > (uint32_t)dcm.imageBytes)) {
printError("RLE header error\n");
free(cImg);
free(bImg);
@@ -2532,6 +2533,13 @@ unsigned char * nii_loadImgRLE(char* imgname, struct nifti_1_header hdr, struct
return bImg;
} // nii_loadImgRLE()
+#ifdef myDisableOpenJPEG
+ #ifndef myEnableJasper
+ //avoid compiler warning, see https://stackoverflow.com/questions/3599160/unused-parameter-warnings-in-c
+ #define UNUSED(x) (void)(x)
+ #endif
+#endif
+
unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct TDICOMdata dcm, bool iVaries, int compressFlag, int isVerbose) {
//provided with a filename (imgname) and DICOM header (dcm), creates NIfTI header (hdr) and img
if (headerDcm2Nii(dcm, hdr, true) == EXIT_FAILURE) return NULL; //TOFU
@@ -2541,7 +2549,8 @@ unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct
printMessage("Software not compiled to decompress classic JPEG DICOM images\n");
return NULL;
#else
- img = nii_loadImgJPEG50(imgname, *hdr, dcm);
+ //img = nii_loadImgJPEG50(imgname, *hdr, dcm);
+ img = nii_loadImgJPEG50(imgname, dcm);
if (hdr->datatype ==DT_RGB24) //convert to planar
img = nii_rgb2planar(img, hdr, dcm.isPlanarRGB);//do this BEFORE Y-Flip, or RGB order can be flipped
#endif
@@ -2562,6 +2571,8 @@ unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct
if ((dcm.compressionScheme == kCompressYes) && (compressFlag != kCompressNone) )
img = nii_loadImgCoreJasper(imgname, *hdr, dcm, compressFlag);
else
+ #else
+ UNUSED(compressFlag); //avoid compiler -Wunused-parameter warning: compressFlag required when myEnableJasper or not myDisableOpenJPEG
#endif
#endif
if (dcm.compressionScheme == kCompressYes) {
@@ -2576,7 +2587,7 @@ unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct
img = nii_rgb2planar(img, hdr, dcm.isPlanarRGB);//do this BEFORE Y-Flip, or RGB order can be flipped
dcm.isPlanarRGB = true;
if (dcm.CSA.mosaicSlices > 1) {
- img = nii_demosaic(img, hdr, dcm.CSA.mosaicSlices, dcm.CSA.protocolSliceNumber1);
+ img = nii_demosaic(img, hdr, dcm.CSA.mosaicSlices); //, dcm.CSA.protocolSliceNumber1);
/* we will do this in nii_dicom_batch #ifdef obsolete_mosaic_flip
img = nii_flipImgY(img, hdr);
#endif*/
@@ -2663,7 +2674,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
#else
size_t MaxBufferSz = 1000000; //ideally size of DICOM header, but this varies from 2D to 4D files
#endif
- if (MaxBufferSz > fileLen)
+ if (MaxBufferSz > (size_t)fileLen)
MaxBufferSz = fileLen;
//printf("%d -> %d\n", MaxBufferSz, fileLen);
long lFileOffset = 0;
@@ -2689,6 +2700,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 kSourceApplicationEntityTitle 0x0002+(0x0016 << 16 )
//#define kSpecificCharacterSet 0x0008+(0x0005 << 16 ) //someday we should handle foreign characters...
#define kImageTypeTag 0x0008+(0x0008 << 16 )
#define kStudyDate 0x0008+(0x0020 << 16 )
@@ -2701,15 +2713,18 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
#define kInstitutionName 0x0008+(0x0080 << 16 )
#define kInstitutionAddress 0x0008+(0x0081 << 16 )
#define kReferringPhysicianName 0x0008+(0x0090 << 16 )
+#define kStationName 0x0008+(0x1010 << 16 )
#define kSeriesDescription 0x0008+(0x103E << 16 ) // '0008' '103E' 'LO' 'SeriesDescription'
#define kManufacturersModelName 0x0008+(0x1090 << 16 )
#define kDerivationDescription 0x0008+(0x2111 << 16 )
#define kComplexImageComponent (uint32_t) 0x0008+(0x9208 << 16 )//'0008' '9208' 'CS' 'ComplexImageComponent'
#define kPatientName 0x0010+(0x0010 << 16 )
#define kPatientID 0x0010+(0x0020 << 16 )
+#define kAnatomicalOrientationType 0x0010+(0x2210 << 16 )
#define kBodyPartExamined 0x0018+(0x0015 << 16)
#define kScanningSequence 0x0018+(0x0020 << 16)
#define kSequenceVariant 0x0018+(0x0021 << 16)
+#define kScanOptions 0x0018+(0x0022 << 16)
#define kMRAcquisitionType 0x0018+(0x0023 << 16)
#define kSequenceName 0x0018+(0x0024 << 16)
#define kZThick 0x0018+(0x0050 << 16 )
@@ -2721,6 +2736,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
#define kZSpacing 0x0018+(0x0088 << 16 ) //'DS' 'SpacingBetweenSlices'
#define kPhaseEncodingSteps 0x0018+(0x0089 << 16 ) //'IS'
#define kEchoTrainLength 0x0018+(0x0091 << 16 ) //IS
+#define kPhaseFieldofView 0x0018+(0x0094 << 16 ) //'DS'
#define kDeviceSerialNumber 0x0018+(0x1000 << 16 ) //LO
#define kSoftwareVersions 0x0018+(0x1020 << 16 ) //LO
#define kProtocolName 0x0018+(0x1030<< 16 )
@@ -2734,6 +2750,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
#define kInPlanePhaseEncodingDirection 0x0018+(0x1312<< 16 ) //CS
#define kPatientOrient 0x0018+(0x5100<< 16 ) //0018,5100. patient orientation - 'HFS'
//#define kDiffusionBFactorSiemens 0x0019+(0x100C<< 16 ) // 0019;000C;SIEMENS MR HEADER ;B_value
+#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
@@ -2775,6 +2792,9 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
#define kCoilSiemens 0x0051+(0x100F << 16 )
#define kImaPATModeText 0x0051+(0x1011 << 16 )
#define kLocationsInAcquisition 0x0054+(0x0081 << 16 )
+//ftp://dicom.nema.org/MEDICAL/dicom/2014c/output/chtml/part03/sect_C.8.9.4.html
+//If ImageType is REPROJECTION we slice direction is reversed - need example to test
+// #define kSeriesType 0x0054+(0x1000 << 16 )
#define kDoseCalibrationFactor 0x0054+(0x1322<< 16 )
#define kIconImageSequence 0x0088+(0x0200 << 16 )
#define kDiffusionBFactor 0x2001+(0x1003 << 16 )// FL
@@ -2800,9 +2820,9 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
#define kImageStart 0x7FE0+(0x0010 << 16 )
#define kImageStartFloat 0x7FE0+(0x0008 << 16 )
#define kImageStartDouble 0x7FE0+(0x0009 << 16 )
-#define kNest 0xFFFE +(0xE000 << 16 ) //Item follows SQ
-#define kUnnest 0xFFFE +(0xE00D << 16 ) //ItemDelimitationItem [length defined] http://www.dabsoft.ch/dicom/5/7.5/
-#define kUnnest2 0xFFFE +(0xE0DD << 16 )//SequenceDelimitationItem [length undefined]
+uint32_t kNest = 0xFFFE +(0xE000 << 16 ); //#define kNest 0xFFFE +(0xE000 << 16 ) //Item follows SQ
+uint32_t kUnnest = 0xFFFE +(0xE00D << 16 ); //#define kUnnest 0xFFFE +(0xE00D << 16 ) //ItemDelimitationItem [length defined] http://www.dabsoft.ch/dicom/5/7.5/
+uint32_t kUnnest2 = 0xFFFE +(0xE0DD << 16 ); //#define kUnnest2 0xFFFE +(0xE0DD << 16 )//SequenceDelimitationItem [length undefined]
int nest = 0;
double zSpacing = -1.0l; //includes slice thickness plus gap
int locationsInAcquisitionGE = 0;
@@ -2840,10 +2860,9 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
float patientPositionStartPhilips[4] = {NAN, NAN, NAN, NAN};
while ((d.imageStart == 0) && ((lPos+8+lFileOffset) < fileLen)) {
#ifndef myLoadWholeFileToReadHeader //read one segment at a time
- if ((lPos + 128) > MaxBufferSz) { //avoid overreading the file
+ if ((size_t)(lPos + 128) > MaxBufferSz) { //avoid overreading the file
lFileOffset = lFileOffset + lPos;
- if ((lFileOffset+MaxBufferSz) > fileLen)
- MaxBufferSz = fileLen - lFileOffset;
+ if ((lFileOffset+MaxBufferSz) > (size_t)fileLen) MaxBufferSz = fileLen - lFileOffset;
fseek(file, lFileOffset, SEEK_SET);
size_t sz = fread(buffer, 1, MaxBufferSz, file);
if (sz < MaxBufferSz) {
@@ -2869,6 +2888,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
} //transfer syntax requests switching VR after group 0001
//uint32_t group = (groupElement & 0xFFFF);
lPos += 4;
+ if ((groupElement == kUnnest) || (groupElement == kUnnest2)) isIconImageSequence = false;
if (((groupElement == kNest) || (groupElement == kUnnest) || (groupElement == kUnnest2)) && (!isEncapsulatedData)) {
//if (((groupElement == kNest) || (groupElement == kUnnest) || (groupElement == kUnnest2)) ) {
vr[0] = 'N';
@@ -2988,6 +3008,13 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
d.imageStart = 1;//abort as invalid (imageStart MUST be >128)
}
break;} //{} provide scope for variable 'transferSyntax
+ case kSourceApplicationEntityTitle: {
+ char saeTxt[kDICOMStr];
+ dcmStr (lLength, &buffer[lPos], saeTxt);
+ int slen = (int) strlen(saeTxt);
+ if((slen < 5) || (strstr(saeTxt, "oasis") == NULL) ) break;
+ d.isSegamiOasis = true;
+ break; }
case kImageTypeTag:
dcmStr (lLength, &buffer[lPos], d.imageType);
int slen;
@@ -2998,7 +3025,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
//isNonImage 0008,0008 = DERIVED,CSAPARALLEL,POSDISP
// attempt to detect non-images, see https://github.com/scitran/data/blob/a516fdc39d75a6e4ac75d0e179e18f3a5fc3c0af/scitran/data/medimg/dcm/mr/siemens.py
if((slen > 6) && (strstr(d.imageType, "DERIVED") != NULL) )
- d.isNonImage = true;
+ d.isDerived = true;
//if((slen > 4) && (strstr(typestr, "DIS2D") != NULL) )
// d.isNonImage = true;
break;
@@ -3057,9 +3084,19 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
case kPatientName :
dcmStr (lLength, &buffer[lPos], d.patientName);
break;
+ case kAnatomicalOrientationType: {
+ char aotTxt[kDICOMStr]; //ftp://dicom.nema.org/MEDICAL/dicom/2015b/output/chtml/part03/sect_C.7.6.2.html#sect_C.7.6.2.1.1
+ dcmStr (lLength, &buffer[lPos], aotTxt);
+ int slen = (int) strlen(aotTxt);
+ if((slen < 9) || (strstr(aotTxt, "QUADRUPED") == NULL) ) break;
+ printError("Anatomical Orientation Type (0010,2210) is QUADRUPED: rotate coordinates accordingly");
+ break; }
case kPatientID :
dcmStr (lLength, &buffer[lPos], d.patientID);
break;
+ case kStationName :
+ dcmStr (lLength, &buffer[lPos], d.stationName);
+ break;
case kSeriesDescription: {
dcmStr (lLength, &buffer[lPos], d.seriesDescription);
break; }
@@ -3088,6 +3125,9 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
case kPatientOrient :
dcmStr (lLength, &buffer[lPos], d.patientOrient);
break;
+ case kDwellTime :
+ d.dwellTime = dcmStrInt(lLength, &buffer[lPos]);
+ break;
case kLastScanLoc :
d.lastScanLoc = dcmStrFloat(lLength, &buffer[lPos]);
break;
@@ -3230,6 +3270,9 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
case kEchoTrainLength :
d.echoTrainLength = dcmStrInt(lLength, &buffer[lPos]);
break;
+ case kPhaseFieldofView :
+ d.phaseFieldofView = dcmStrFloat(lLength, &buffer[lPos]);
+ break;
case kAcquisitionMatrix :
if (lLength == 8) {
uint16_t acquisitionMatrix[4];
@@ -3293,7 +3336,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
char accelStr[kDICOMStr];
dcmStr (lLength, &buffer[lPos], accelStr);
char *ptr;
- dcmStrDigitsOnly(accelStr);
+ dcmStrDigitsOnlyKey('p', accelStr); //e.g. if "p2s4" return "2", if "s4" return ""
d.accelFactPE = (float)strtof(accelStr, &ptr);
if (*ptr != '\0')
d.accelFactPE = 0.0;
@@ -3329,6 +3372,9 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
dcmStr (lLength, &buffer[lPos], d.sequenceVariant);
break;
}
+ case kScanOptions:
+ dcmStr (lLength, &buffer[lPos], d.scanOptions);
+ break;
case kSequenceName : {
//if (strlen(d.protocolName) < 1) //precedence given to kProtocolName and kProtocolNameGE
dcmStr (lLength, &buffer[lPos], d.sequenceName);
@@ -3455,7 +3501,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
printMessage("Skipping DICOM (audio not image) '%s'\n", fname);
break;
case kCSAImageHeaderInfo:
- readCSAImageHeader(&buffer[lPos], lLength, &d.CSA, isVerbose, dti4D);
+ readCSAImageHeader(&buffer[lPos], lLength, &d.CSA, isVerbose); //, dti4D);
d.isHasPhase = d.CSA.isPhaseMap;
break;
//case kObjectGraphics:
@@ -3539,7 +3585,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
tagStr[0] = 'X'; //avoid compiler warning: orientStr filled by dcmStr
dcmStr (lLength, &buffer[lPos], tagStr);
if (strlen(tagStr) > 1) {
- for (int pos = 0; pos<strlen(tagStr); pos ++)
+ for (size_t pos = 0; pos<strlen(tagStr); pos ++)
if ((tagStr[pos] == '<') || (tagStr[pos] == '>') || (tagStr[pos] == ':')
|| (tagStr[pos] == '"') || (tagStr[pos] == '\\') || (tagStr[pos] == '/')
|| (tagStr[pos] == '^') || (tagStr[pos] < 33)
diff --git a/console/nii_dicom.h b/console/nii_dicom.h
index fcdefc5..f1f9e5c 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.20170818" kDCMsuf kCCsuf
+#define kDCMvers "v1.0.20170923" 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
@@ -48,7 +48,6 @@ static const int kMaxDTI4D = 4096; //maximum number of DTI directions for 4D (Ph
#define kMANUFACTURER_GE 2
#define kMANUFACTURER_PHILIPS 3
#define kMANUFACTURER_TOSHIBA 4
-#define kMANUFACTURER_SEGAMI 5
//note: note a complete modality list, e.g. XA,PX, etc
#define kMODALITY_UNKNOWN 0
@@ -58,7 +57,6 @@ static const int kMaxDTI4D = 4096; //maximum number of DTI directions for 4D (Ph
#define kMODALITY_PT 4
#define kMODALITY_US 5
-
#define kEXIT_NO_VALID_FILES_FOUND 2
static const int kSliceOrientUnknown = 0;
static const int kSliceOrientTra = 1;
@@ -113,15 +111,15 @@ static const int kCompressRLE = 4; //run length encoding
struct TDICOMdata {
long seriesNum;
int xyzDim[5];
- int modality, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, echoTrainLength, patientPositionNumPhilips, coilNum, echoNum, sliceOrient,numberOfDynamicScans, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,patientPositionSequentialRepeats,locationsInAcquisition, compressionScheme;
- float accelFactPE, flipAngle, fieldStrength, TE, TI, TR, intenScale, intenIntercept, intenScalePhilips, gantryTilt, lastScanLoc, angulation[4];
+ int modality, dwellTime, effectiveEchoSpacingGE, phaseEncodingLines, phaseEncodingSteps, echoTrainLength, patientPositionNumPhilips, coilNum, echoNum, sliceOrient,numberOfDynamicScans, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,patientPositionSequentialRepeats,locationsInAcquisition, compressionScheme;
+ float phaseFieldofView, accelFactPE, flipAngle, fieldStrength, TE, TI, TR, intenScale, intenIntercept, intenScalePhilips, gantryTilt, lastScanLoc, angulation[4];
float orient[7], patientPosition[4], patientPositionLast[4], xyzMM[4], stackOffcentre[4];
float radionuclidePositronFraction, radionuclideTotalDose, radionuclideHalfLife, doseCalibrationFactor; //PET ISOTOPE MODULE ATTRIBUTES (C.8-57)
float ecat_isotope_halflife, ecat_dosage;
double dateTime, acquisitionTime, acquisitionDate, bandwidthPerPixelPhaseEncode;
- bool isXRay, isMultiEcho, isSlicesSpatiallySequentialPhilips, isNonImage, isValid, is3DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase,isHasMagnitude,isHasMixed, isFloat, isResampled;
+ bool isSegamiOasis, isDerived, isXRay, isMultiEcho, isSlicesSpatiallySequentialPhilips, isValid, is3DAcq, isExplicitVR, isLittleEndian, isPlanarRGB, isSigned, isHasPhase,isHasMagnitude,isHasMixed, isFloat, isResampled;
char phaseEncodingRC;
- char softwareVersions[kDICOMStr], deviceSerialNumber[kDICOMStr], institutionAddress[kDICOMStr], institutionName[kDICOMStr], referringPhysicianName[kDICOMStr], seriesInstanceUID[kDICOMStr], studyInstanceUID[kDICOMStr], bodyPartExamined[kDICOMStr], procedureStepDescription[kDICOMStr], imageType[kDICOMStr], manufacturersModelName[kDICOMStr], patientID[kDICOMStr], patientOrient[kDICOMStr], patientName[kDICOMStr],seriesDescription[kDICOMStr], studyID[kDICOMStr], sequenceName[kDICOMStr [...]
+ 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], manufacturersModelName[kDICOMStr], patientID[kDICOMStr], patientOrient[kDICOMStr], patientName[kDICOMStr],seriesDescription[kDICOM [...]
struct TCSAdata CSA;
};
@@ -139,9 +137,7 @@ static const int kCompressRLE = 4; //run length encoding
void setQSForm(struct nifti_1_header *h, mat44 Q44i, bool isVerbose);
int headerDcm2Nii2(struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1_header *h, int isVerbose);
int headerDcm2Nii(struct TDICOMdata d, struct nifti_1_header *h, bool isComputeSForm) ;
- //unsigned char * nii_loadImgX(char* imgname, struct nifti_1_header *hdr, struct TDICOMdata dcm, bool iVaries);
unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct TDICOMdata dcm, bool iVaries, int compressFlag, int isVerbose);
- //int foo (float vx);
#ifdef __cplusplus
}
#endif
diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp
old mode 100755
new mode 100644
index 484150d..55bfe23
--- a/console/nii_dicom_batch.cpp
+++ b/console/nii_dicom_batch.cpp
@@ -78,6 +78,10 @@ struct TSearchList {
char **str;
};
+#ifndef PATH_MAX
+ #define PATH_MAX 4096
+#endif
+
void dropFilenameFromPath(char *path) { //
const char *dirPath = strrchr(path, '/'); //UNIX
if (dirPath == 0)
@@ -151,21 +155,36 @@ bool is_exe(const char* path) { //requires #include <sys/stat.h>
//return (S_ISREG(buf.st_mode) && (buf.st_mode & 0111) );
} //is_exe()
-int is_dir(const char *pathname, int follow_link) {
- struct stat s;
- if ((NULL == pathname) || (0 == strlen(pathname)))
- return 0;
- int err = stat(pathname, &s);
- if(-1 == err)
- return 0; /* does not exist */
- else {
- if(S_ISDIR(s.st_mode)) {
- return 1; /* it's a dir */
- } else {
- return 0;/* exists but is no dir */
- }
- }
-}// is_dir()
+#if defined(_WIN64) || defined(_WIN32)
+ //Windows does not support lstat
+ int is_dir(const char *pathname, int follow_link) {
+ struct stat s;
+ if ((NULL == pathname) || (0 == strlen(pathname)))
+ return 0;
+ int err = stat(pathname, &s);
+ if(-1 == err)
+ return 0; // does not exist
+ else {
+ if(S_ISDIR(s.st_mode)) {
+ return 1; // it's a dir
+ } else {
+ return 0;// exists but is no dir
+ }
+ }
+ }// is_dir()
+#else //if windows else Unix
+ int is_dir(const char *pathname, int follow_link)
+ {
+ struct stat s;
+ int retval;
+ if ((NULL == pathname) || (0 == strlen(pathname)))
+ return 0; // does not exist
+ retval = follow_link ? stat(pathname, &s) : lstat(pathname, &s);
+ if ((-1 != retval) && (S_ISDIR(s.st_mode)))
+ return 1; // it's a dir
+ return 0; // exists but is no dir
+ }// is_dir()
+#endif
void geCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI *vx){
//0018,1312 phase encoding is either in row or column direction
@@ -203,8 +222,10 @@ void geCorrectBvecs(struct TDICOMdata *d, int sliceDir, struct TDTI *vx){
flp = setiVec3(0, 1, 1); //CORONAL
else if (abs(sliceDir) == 3)
flp = setiVec3(0, 0, 1); //AXIAL
- else
+ else {
printMessage("Impossible GE slice orientation!");
+ flp = setiVec3(0, 0, 1); //AXIAL???
+ }
if (sliceDir < 0)
flp.v[2] = 1 - flp.v[2];
printMessage("Saving %d DTI gradients. GE Reorienting %s : please validate. isCol=%d sliceDir=%d flp=%d %d %d\n", d->CSA.numDti, d->protocolName, col, sliceDir, flp.v[0], flp.v[1],flp.v[2]);
@@ -324,13 +345,7 @@ void nii_SaveText(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts,
fclose(fp);
}// nii_SaveText()
-bool isDerived(struct TDICOMdata d) {
- #define kDerivedStr "DERIVED"
- if ((strlen(d.imageType) < strlen(kDerivedStr)) || (strstr(d.imageType, kDerivedStr) == NULL))
- return false;
- else
- return true;
-}
+#define myReadAsciiCsa
#ifdef myReadAsciiCsa
//read from the ASCII portion of the Siemens CSA series header
@@ -402,6 +417,49 @@ int readKey(const char * key, char * buffer, int remLength) { //look for text k
return ret;
} //readKey()
+float readKeyFloat(const char * key, char * buffer, int remLength) { //look for text key in binary data stream, return subsequent integer value
+ char *keyPos = (char *)memmem(buffer, remLength, key, strlen(key));
+ if (!keyPos) return 0.0;
+ char str[kDICOMStr];
+ strcpy(str, "");
+ char tmpstr[2];
+ tmpstr[1] = 0;
+ int i = (int)strlen(key);
+ while( ( i< remLength) && (keyPos[i] != 0x0A) ) {
+ if( (keyPos[i] >= '0' && keyPos[i] <= '9') || (keyPos[i] <= '.') || (keyPos[i] <= '-') ) {
+ tmpstr[0] = keyPos[i];
+ strcat (str, tmpstr);
+ }
+ i++;
+ }
+ if (strlen(str) < 1) return 0.0;
+ return atof(str);
+} //readKeyFloat()
+
+void readKeyStr(const char * key, char * buffer, int remLength, char* outStr) {
+//if key is CoilElementID.tCoilID the string 'CoilElementID.tCoilID = ""Head_32""' returns 'Head32'
+ strcpy(outStr, "");
+ char *keyPos = (char *)memmem(buffer, remLength, key, strlen(key));
+ if (!keyPos) return;
+ int i = (int)strlen(key);
+ int outLen = 0;
+ char tmpstr[2];
+ tmpstr[1] = 0;
+ bool isQuote = false;
+ while( ( i < remLength) && (keyPos[i] != 0x0A) ) {
+ if ((isQuote) && (keyPos[i] != '"') && (outLen < kDICOMStr)) {
+ tmpstr[0] = keyPos[i];
+ strcat (outStr, tmpstr);
+ outLen ++;
+ }
+ if (keyPos[i] == '"') {
+ if (outLen > 0) break;
+ isQuote = true;
+ }
+ i++;
+ }
+} //readKeyStr()
+
int phoenixOffsetCSASeriesHeader(unsigned char *buff, int lLength) {
//returns offset to ASCII Phoenix data
if (lLength < 36) return 0;
@@ -432,34 +490,41 @@ int phoenixOffsetCSASeriesHeader(unsigned char *buff, int lLength) {
return 0;
} // phoenixOffsetCSASeriesHeader()
-int siemensEchoEPIFactor(const char * filename, int csaOffset, int csaLength, int* echoSpacing, int* echoTrainDuration) {
+void siemensCsaAscii(const char * filename, int csaOffset, int csaLength, float* phaseOversampling, float* phaseResolution, int* baseResolution, int* interp, int* partialFourier, int* echoSpacing, int* parallelReductionFactorInPlane, char* coilID, char* consistencyInfo, char* coilElements, char* pulseSequenceDetails, char* fmriExternalInfo) {
//reads ASCII portion of CSASeriesHeaderInfo and returns lEchoTrainDuration or lEchoSpacing value
// returns 0 if no value found
+ *phaseOversampling = 0.0;
+ *phaseResolution = 0.0;
+ *baseResolution = 0;
+ *interp = 0;
+ *partialFourier = 0;
*echoSpacing = 0;
- *echoTrainDuration = 0;
- if ((csaOffset < 0) || (csaLength < 8)) return 0;
+ strcpy(coilID, "");
+ strcpy(consistencyInfo, "");
+ strcpy(coilElements, "");
+ strcpy(pulseSequenceDetails, "");
+ if ((csaOffset < 0) || (csaLength < 8)) return;
FILE * pFile = fopen ( filename, "rb" );
- if(pFile==NULL) return 0;
+ if(pFile==NULL) return;
fseek (pFile , 0 , SEEK_END);
long lSize = ftell (pFile);
if (lSize < (csaOffset+csaLength)) {
fclose (pFile);
- return 0;
+ return;
}
fseek(pFile, csaOffset, SEEK_SET);
char * buffer = (char*) malloc (csaLength);
- if(buffer == NULL) return 0;
+ if(buffer == NULL) return;
size_t result = fread (buffer,1,csaLength,pFile);
- if(result != csaLength) return 0;
+ if ((int)result != csaLength) return;
//next bit complicated: restrict to ASCII portion to avoid buffer overflow errors in BINARY portion
int startAscii = phoenixOffsetCSASeriesHeader((unsigned char *)buffer, csaLength);
int csaLengthTrim = csaLength;
char * bufferTrim = buffer;
- if ((startAscii > 0) && (startAscii < csaLengthTrim) ){ //ignore binary data at start
+ if ((startAscii > 0) && (startAscii < csaLengthTrim) ) { //ignore binary data at start
bufferTrim += startAscii;
csaLengthTrim -= startAscii;
}
- int ret = 0;
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) {
@@ -470,19 +535,66 @@ int siemensEchoEPIFactor(const char * filename, int csaOffset, int csaLength, i
csaLengthTrim = (int)(keyPosEnd - keyPos);
char keyStrES[] = "sFastImaging.lEchoSpacing";
*echoSpacing = readKey(keyStrES, keyPos, csaLengthTrim);
- char keyStrETD[] = "sFastImaging.lEchoTrainDuration";
- *echoTrainDuration = readKey(keyStrETD, keyPos, csaLengthTrim);
- char keyStrEF[] = "sFastImaging.lEPIFactor";
- ret = readKey(keyStrEF, keyPos, csaLengthTrim);
+ char keyStrBase[] = "sKSpace.lBaseResolution";
+ *baseResolution = readKey(keyStrBase, keyPos, csaLengthTrim);
+ char keyStrInterp[] = "sKSpace.uc2DInterpolation";
+ *interp = readKey(keyStrInterp, keyPos, csaLengthTrim);
+ char keyStrPF[] = "sKSpace.ucPhasePartialFourier";
+ *partialFourier = readKey(keyStrPF, keyPos, csaLengthTrim);
+ //char keyStrETD[] = "sFastImaging.lEchoTrainDuration";
+ //*echoTrainDuration = readKey(keyStrETD, keyPos, csaLengthTrim);
+ char keyStrAF[] = "sPat.lAccelFactPE";
+ *parallelReductionFactorInPlane = readKey(keyStrAF, keyPos, csaLengthTrim);
+ //char keyStrEF[] = "sFastImaging.lEPIFactor";
+ //ret = readKey(keyStrEF, keyPos, csaLengthTrim);
+ char keyStrCoil[] = "sCoilElementID.tCoilID";
+ readKeyStr(keyStrCoil, keyPos, csaLengthTrim, coilID);
+ char keyStrCI[] = "sProtConsistencyInfo.tMeasuredBaselineString";
+ readKeyStr(keyStrCI, keyPos, csaLengthTrim, consistencyInfo);
+ char keyStrCS[] = "sCoilSelectMeas.sCoilStringForConversion";
+ readKeyStr(keyStrCS, keyPos, csaLengthTrim, coilElements);
+ char keyStrSeq[] = "tSequenceFileName";
+ readKeyStr(keyStrSeq, keyPos, csaLengthTrim, pulseSequenceDetails);
+ char keyStrExt[] = "FmriExternalInfo";
+ readKeyStr(keyStrExt, keyPos, csaLengthTrim, fmriExternalInfo);
+ char keyStrPhase[] = "sKSpace.dPhaseResolution";
+ *phaseResolution = readKeyFloat(keyStrPhase, keyPos, csaLengthTrim);
+ char keyStrOver[] = "sKSpace.dPhaseOversamplingForDialog";
+ *phaseOversampling = readKeyFloat(keyStrOver, keyPos, csaLengthTrim);
}
fclose (pFile);
free (buffer);
- return ret;
-}
-
-#endif //myReadAsciiCsa
-
-void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, struct TDTI4D *dti4D, struct nifti_1_header *h, const char * filename) {
+ return;
+} // siemensCsaAscii()
+#endif //myReadAsciiCsa()
+
+void json_Str(FILE *fp, const char *sLabel, char *sVal) {
+ if (strlen(sVal) < 1) return;
+ //fprintf(fp, sLabel, sVal );
+ //convert \ ' " characters to _ see https://github.com/rordenlab/dcm2niix/issues/131
+ for (size_t pos = 0; pos < strlen(sVal); pos ++) {
+ if ((sVal[pos] == '\'') || (sVal[pos] == '"') || (sVal[pos] == '\\'))
+ sVal[pos] = '_';
+ }
+ fprintf(fp, sLabel, sVal );
+ /*char outname[PATH_MAX] = {""};
+ char appendChar[2] = {"\\"};
+ char passChar[2] = {"\\"};
+ for (int pos = 0; pos<strlen(sVal); pos ++) {
+ if ((sVal[pos] == '\'') || (sVal[pos] == '"') || (sVal[pos] == '\\'))
+ strcpy(outname, appendChar);
+ passChar[0] = sVal[pos];
+ strcpy(outname, passChar);
+ }
+ fprintf(fp, sLabel, outname );*/
+} //json_Str
+
+void json_Float(FILE *fp, const char *sLabel, float sVal) {
+ if (sVal <= 0.0) return;
+ fprintf(fp, sLabel, sVal );
+} //json_Float
+
+void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, struct nifti_1_header *h, const char * filename) {
//https://docs.google.com/document/d/1HFUkAEE-pB-angVcYe6pf_-fVf4sCpOHKesUvfb8Grc/edit#
// Generate Brain Imaging Data Structure (BIDS) info
// sidecar JSON file (with the same filename as the .nii.gz file, but with .json extension).
@@ -512,6 +624,7 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts,
fprintf(fp, "\t\"Modality\": \"US\",\n" );
break;
};
+ if (d.fieldStrength > 0.0) fprintf(fp, "\t\"MagneticFieldStrength\": %g,\n", d.fieldStrength );
switch (d.manufacturer) {
case kMANUFACTURER_SIEMENS:
fprintf(fp, "\t\"Manufacturer\": \"Siemens\",\n" );
@@ -525,74 +638,34 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts,
case kMANUFACTURER_TOSHIBA:
fprintf(fp, "\t\"Manufacturer\": \"Toshiba\",\n" );
break;
- case kMANUFACTURER_SEGAMI:
- fprintf(fp, "\t\"Manufacturer\": \"Segami\",\n" );
- break;
};
- fprintf(fp, "\t\"ManufacturersModelName\": \"%s\",\n", d.manufacturersModelName );
+ json_Str(fp, "\t\"ManufacturersModelName\": \"%s\",\n", d.manufacturersModelName);
+ json_Str(fp, "\t\"InstitutionName\": \"%s\",\n", d.institutionName);
+ json_Str(fp, "\t\"InstitutionAddress\": \"%s\",\n", d.institutionAddress);
+ json_Str(fp, "\t\"DeviceSerialNumber\": \"%s\",\n", d.deviceSerialNumber );
+ json_Str(fp, "\t\"StationName\": \"%s\",\n", d.stationName );
if (!opts.isAnonymizeBIDS) {
- if (strlen(d.seriesInstanceUID) > 0)
- fprintf(fp, "\t\"SeriesInstanceUID\": \"%s\",\n", d.seriesInstanceUID );
- if (strlen(d.studyInstanceUID) > 0)
- fprintf(fp, "\t\"StudyInstanceUID\": \"%s\",\n", d.studyInstanceUID );
- if (strlen(d.referringPhysicianName) > 0)
- fprintf(fp, "\t\"ReferringPhysicianName\": \"%s\",\n", d.referringPhysicianName );
- if (strlen(d.studyID) > 0)
- fprintf(fp, "\t\"StudyID\": \"%s\",\n", d.studyID );
+ json_Str(fp, "\t\"SeriesInstanceUID\": \"%s\",\n", d.seriesInstanceUID);
+ json_Str(fp, "\t\"StudyInstanceUID\": \"%s\",\n", d.studyInstanceUID);
+ json_Str(fp, "\t\"ReferringPhysicianName\": \"%s\",\n", d.referringPhysicianName);
+ json_Str(fp, "\t\"StudyID\": \"%s\",\n", d.studyID);
//Next lines directly reveal patient identity
- //if (strlen(d.patientName) > 0)
- // fprintf(fp, "\t\"PatientName\": \"%s\",\n", d.patientName );
- //if (strlen(d.patientID) > 0)
- // fprintf(fp, "\t\"PatientID\": \"%s\",\n", d.patientID );
+ // json_Str(fp, "\t\"PatientName\": \"%s\",\n", d.patientName);
+ // json_Str(fp, "\t\"PatientID\": \"%s\",\n", d.patientID);
}
- #ifdef myReadAsciiCsa
- if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.CSA.SeriesHeader_offset > 0) && (d.CSA.SeriesHeader_length > 0) &&
- (strlen(d.scanningSequence) > 1) && (d.scanningSequence[0] == 'E') && (d.scanningSequence[1] == 'P')) { //for EPI scans only
- int echoSpacing, echoTrainDuration, epiFactor;
- epiFactor = siemensEchoEPIFactor(filename, d.CSA.SeriesHeader_offset, d.CSA.SeriesHeader_length, &echoSpacing, &echoTrainDuration);
- //printMessage("ES %d ETD %d EPI %d\n", echoSpacing, echoTrainDuration, epiFactor);
- if (echoSpacing > 0)
- fprintf(fp, "\t\"EchoSpacing\": %g,\n", echoSpacing / 1000000.0); //usec -> sec
- if (echoTrainDuration > 0)
- fprintf(fp, "\t\"EchoTrainDuration\": %g,\n", echoTrainDuration / 1000000.0); //usec -> sec
- if (epiFactor > 0)
- fprintf(fp, "\t\"EPIFactor\": %d,\n", epiFactor);
- }
- #endif
- if (d.echoTrainLength > 1) //>1 as for Siemens EPI this is 1, Siemens uses EPI factor http://mriquestions.com/echo-planar-imaging.html
- fprintf(fp, "\t\"EchoTrainLength\": %d,\n", d.echoTrainLength);
- if (d.echoNum > 1)
- fprintf(fp, "\t\"EchoNumber\": %d,\n", d.echoNum);
- if (d.isNonImage) //DICOM is derived image or non-spatial file (sounds, etc)
- fprintf(fp, "\t\"RawImage\": false,\n");
- if (d.acquNum > 0)
- fprintf(fp, "\t\"AcquisitionNumber\": %d,\n", d.acquNum);
- if (strlen(d.institutionName) > 0)
- fprintf(fp, "\t\"InstitutionName\": \"%s\",\n", d.institutionName );
- if (strlen(d.institutionAddress) > 0)
- fprintf(fp, "\t\"InstitutionAddress\": \"%s\",\n", d.institutionAddress );
- if (strlen(d.deviceSerialNumber) > 0)
- fprintf(fp, "\t\"DeviceSerialNumber\": \"%s\",\n", d.deviceSerialNumber );
- if (strlen(d.softwareVersions) > 0)
- fprintf(fp, "\t\"SoftwareVersions\": \"%s\",\n", d.softwareVersions );
- if (strlen(d.procedureStepDescription) > 0)
- fprintf(fp, "\t\"ProcedureStepDescription\": \"%s\",\n", d.procedureStepDescription );
- if (strlen(d.scanningSequence) > 0)
- fprintf(fp, "\t\"ScanningSequence\": \"%s\",\n", d.scanningSequence );
- if (strlen(d.sequenceVariant) > 0)
- fprintf(fp, "\t\"SequenceVariant\": \"%s\",\n", d.sequenceVariant );
- if (strlen(d.seriesDescription) > 0)
- fprintf(fp, "\t\"SeriesDescription\": \"%s\",\n", d.seriesDescription );
- if (strlen(d.bodyPartExamined) > 0)
- fprintf(fp, "\t\"BodyPartExamined\": \"%s\",\n", d.bodyPartExamined );
- if (strlen(d.protocolName) > 0)
- fprintf(fp, "\t\"ProtocolName\": \"%s\",\n", d.protocolName );
- if (strlen(d.sequenceName) > 0)
- fprintf(fp, "\t\"SequenceName\": \"%s\",\n", d.sequenceName );
+ json_Str(fp, "\t\"BodyPartExamined\": \"%s\",\n", d.bodyPartExamined);
+ json_Str(fp, "\t\"ProcedureStepDescription\": \"%s\",\n", d.procedureStepDescription);
+ json_Str(fp, "\t\"SoftwareVersions\": \"%s\",\n", d.softwareVersions);
+ json_Str(fp, "\t\"SeriesDescription\": \"%s\",\n", d.seriesDescription);
+ json_Str(fp, "\t\"ProtocolName\": \"%s\",\n", d.protocolName);
+ json_Str(fp, "\t\"ScanningSequence\": \"%s\",\n", d.scanningSequence);
+ json_Str(fp, "\t\"SequenceVariant\": \"%s\",\n", d.sequenceVariant);
+ json_Str(fp, "\t\"ScanOptions\": \"%s\",\n", d.scanOptions);
+ json_Str(fp, "\t\"SequenceName\": \"%s\",\n", d.sequenceName);
if (strlen(d.imageType) > 0) {
fprintf(fp, "\t\"ImageType\": [\"");
bool isSep = false;
- for (int i = 0; i < strlen(d.imageType); i++) {
+ for (size_t i = 0; i < strlen(d.imageType); i++) {
if (d.imageType[i] != '_') {
if (isSep)
fprintf(fp, "\", \"");
@@ -603,6 +676,8 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts,
}
fprintf(fp, "\"],\n");
}
+ if (d.isDerived) //DICOM is derived image or non-spatial file (sounds, etc)
+ fprintf(fp, "\t\"RawImage\": false,\n");
//Chris Gorgolewski: BIDS standard specifies ISO8601 date-time format (Example: 2016-07-06T12:49:15.679688)
//Lines below directly save DICOM values
if (d.acquisitionTime > 0.0 && d.acquisitionDate > 0.0){
@@ -625,12 +700,15 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts,
fprintf(fp, "\t\"AcquisitionDateTime\": ");
fprintf(fp, (ayear >= 0 && ayear <= 9999) ? "\"%4d" : "\"%+4d", ayear);
fprintf(fp, "-%02d-%02dT%02d:%02d:%02.6f\",\n", amonth, aday, ahour, amin, asec);
-
}
} //if (count)
} //if acquisitionTime and acquisitionDate recorded
// if (d.acquisitionTime > 0.0) fprintf(fp, "\t\"AcquisitionTime\": %f,\n", d.acquisitionTime );
// if (d.acquisitionDate > 0.0) fprintf(fp, "\t\"AcquisitionDate\": %8.0f,\n", d.acquisitionDate );
+ if (d.acquNum > 0)
+ fprintf(fp, "\t\"AcquisitionNumber\": %d,\n", d.acquNum);
+ json_Str(fp, "\t\"ImageComments\": \"%s\",\n", d.imageComments);
+ json_Str(fp, "\t\"ConversionComments\": \"%s\",\n", opts.imageComments);
//if conditionals: the following values are required for DICOM MRI, but not available for CT
if ((d.intenScalePhilips != 0) || (d.manufacturer == kMANUFACTURER_PHILIPS)) { //for details, see PhilipsPrecise()
fprintf(fp, "\t\"PhilipsRescaleSlope\": %g,\n", d.intenScale );
@@ -639,64 +717,132 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts,
fprintf(fp, "\t\"UsePhilipsFloatNotDisplayScaling\": %d,\n", opts.isPhilipsFloatNotDisplayScaling);
}
//PET ISOTOPE MODULE ATTRIBUTES
- if (d.radionuclidePositronFraction > 0.0) fprintf(fp, "\t\"RadionuclidePositronFraction\": %g,\n", d.radionuclidePositronFraction );
- if (d.radionuclideTotalDose > 0.0) fprintf(fp, "\t\"RadionuclideTotalDose\": %g,\n", d.radionuclideTotalDose );
- if (d.radionuclideHalfLife > 0.0) fprintf(fp, "\t\"RadionuclideHalfLife\": %g,\n", d.radionuclideHalfLife );
- if (d.doseCalibrationFactor > 0.0) fprintf(fp, "\t\"DoseCalibrationFactor\": %g,\n", d.doseCalibrationFactor );
- //MRI parameters
- if (d.fieldStrength > 0.0) fprintf(fp, "\t\"MagneticFieldStrength\": %g,\n", d.fieldStrength );
- if (d.flipAngle > 0.0) fprintf(fp, "\t\"FlipAngle\": %g,\n", d.flipAngle );
+ json_Float(fp, "\t\"RadionuclidePositronFraction\": %g,\n", d.radionuclidePositronFraction );
+ json_Float(fp, "\t\"RadionuclideTotalDose\": %g,\n", d.radionuclideTotalDose );
+ json_Float(fp, "\t\"RadionuclideHalfLife\": %g,\n", d.radionuclideHalfLife );
+ json_Float(fp, "\t\"DoseCalibrationFactor\": %g,\n", d.doseCalibrationFactor );
+ json_Float(fp, "\t\"IsotopeHalfLife\": %g,\n", d.ecat_isotope_halflife);
+ json_Float(fp, "\t\"Dosage\": %g,\n", d.ecat_dosage);
+ //CT parameters
+ if ((d.TE > 0.0) && (d.isXRay)) fprintf(fp, "\t\"XRayExposure\": %g,\n", d.TE );
+ //MRI parameters
+ if (d.echoNum > 1) fprintf(fp, "\t\"EchoNumber\": %d,\n", d.echoNum);
if ((d.TE > 0.0) && (!d.isXRay)) fprintf(fp, "\t\"EchoTime\": %g,\n", d.TE / 1000.0 );
- if ((d.TE > 0.0) && (d.isXRay)) fprintf(fp, "\t\"XRayExposure\": %g,\n", d.TE );
- if (d.TR > 0.0) fprintf(fp, "\t\"RepetitionTime\": %g,\n", d.TR / 1000.0 );
- if (d.TI > 0.0) fprintf(fp, "\t\"InversionTime\": %g,\n", d.TI / 1000.0 );
- if (d.ecat_isotope_halflife > 0.0) fprintf(fp, "\t\"IsotopeHalfLife\": %g,\n", d.ecat_isotope_halflife);
- if (d.ecat_dosage > 0.0) fprintf(fp, "\t\"Dosage\": %g,\n", d.ecat_dosage);
- double bandwidthPerPixelPhaseEncode = d.bandwidthPerPixelPhaseEncode;
- int phaseEncodingLines = d.phaseEncodingLines;
- if ((phaseEncodingLines == 0) && (h->dim[2] > 0) && (h->dim[1] > 0)) {
+ json_Float(fp, "\t\"RepetitionTime\": %g,\n", d.TR / 1000.0 );
+ json_Float(fp, "\t\"InversionTime\": %g,\n", d.TI / 1000.0 );
+ json_Float(fp, "\t\"FlipAngle\": %g,\n", d.flipAngle );
+ float pf = 1.0f; //partial fourier
+ bool interp = false; //2D interpolation
+ float phaseOversampling = 0.0;
+ #ifdef myReadAsciiCsa
+ if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.CSA.SeriesHeader_offset > 0) && (d.CSA.SeriesHeader_length > 0)) {
+ int baseResolution, interpInt, partialFourier, echoSpacing, parallelReductionFactorInPlane;
+ float phaseResolution;
+ char fmriExternalInfo[kDICOMStr], coilID[kDICOMStr], consistencyInfo[kDICOMStr], coilElements[kDICOMStr], pulseSequenceDetails[kDICOMStr];
+ siemensCsaAscii(filename, d.CSA.SeriesHeader_offset, d.CSA.SeriesHeader_length, &phaseOversampling, &phaseResolution, &baseResolution, &interpInt, &partialFourier, &echoSpacing, ¶llelReductionFactorInPlane, coilID, consistencyInfo, coilElements, pulseSequenceDetails, fmriExternalInfo);
+ if (partialFourier > 0) {
+ //https://github.com/ismrmrd/siemens_to_ismrmrd/blob/master/parameter_maps/IsmrmrdParameterMap_Siemens_EPI_FLASHREF.xsl
+ if (partialFourier == 1) pf = 0.5; // 4/8
+ if (partialFourier == 2) pf = 0.625; // 5/8
+ if (partialFourier == 4) pf = 0.75;
+ if (partialFourier == 8) pf = 0.875;
+ fprintf(fp, "\t\"PartialFourier\": %g,\n", pf);
+ }
+ if (interpInt > 0) {
+ interp = true;
+ fprintf(fp, "\t\"Interpolation2D\": %d,\n", interp);
+ }
+ if (baseResolution > 0) fprintf(fp, "\t\"BaseResolution\": %d,\n", baseResolution );
+ json_Float(fp, "\t\"PhaseResolution\": %g,\n", phaseResolution);
+ json_Float(fp, "\t\"PhaseOversampling\": %g,\n", phaseOversampling); //usec -> sec
+ json_Float(fp, "\t\"VendorReportedEchoSpacing\": %g,\n", echoSpacing / 1000000.0); //usec -> sec
+ //ETD and epiFactor not useful/reliable https://github.com/rordenlab/dcm2niix/issues/127
+ //if (echoTrainDuration > 0) fprintf(fp, "\t\"EchoTrainDuration\": %g,\n", echoTrainDuration / 1000000.0); //usec -> sec
+ //if (epiFactor > 0) fprintf(fp, "\t\"EPIFactor\": %d,\n", epiFactor);
+ json_Str(fp, "\t\"ReceiveCoilName\": \"%s\",\n", coilID);
+ json_Str(fp, "\t\"ReceiveCoilActiveElements\": \"%s\",\n", coilElements);
+ json_Str(fp, "\t\"PulseSequenceDetails\": \"%s\",\n", pulseSequenceDetails);
+ json_Str(fp, "\t\"FmriExternalInfo\": \"%s\",\n", fmriExternalInfo);
+ json_Str(fp, "\t\"ConsistencyInfo\": \"%s\",\n", consistencyInfo);
+ if (parallelReductionFactorInPlane > 0) {//AccelFactorPE -> phase encoding
+ if (d.accelFactPE < 1.0) { //value not found in DICOM header, but WAS found in CSA ascii
+ d.accelFactPE = parallelReductionFactorInPlane; //value found in ASCII but not in DICOM (0051,1011)
+ //fprintf(fp, "\t\"ParallelReductionFactorInPlane\": %g,\n", d.accelFactPE);
+ }
+ if (parallelReductionFactorInPlane != (int)(d.accelFactPE))
+ printWarning("ParallelReductionFactorInPlane reported in DICOM [0051,1011] (%d) does not match CSA series value %d\n", (int)(d.accelFactPE), parallelReductionFactorInPlane);
+ }
+ }
+ #endif
+ if (d.CSA.multiBandFactor > 1) //AccelFactorSlice
+ fprintf(fp, "\t\"MultibandAccelerationFactor\": %d,\n", d.CSA.multiBandFactor);
+ json_Float(fp, "\t\"PercentPhaseFOV\": %g,\n", d.phaseFieldofView);
+ if (d.echoTrainLength > 1) //>1 as for Siemens EPI this is 1, Siemens uses EPI factor http://mriquestions.com/echo-planar-imaging.html
+ fprintf(fp, "\t\"EchoTrainLength\": %d,\n", d.echoTrainLength); //0018,0091 Combination of partial fourier and in-plane parallel imaging
+ if (d.phaseEncodingSteps > 0) fprintf(fp, "\t\"PhaseEncodingSteps\": %d,\n", d.phaseEncodingSteps );
+ if (d.phaseEncodingLines > 0) fprintf(fp, "\t\"AcquisitionMatrixPE\": %d,\n", d.phaseEncodingLines );
+
+ //Compute ReconMatrixPE
+ // Actual size of the *reconstructed* data in the PE dimension, which does NOT match
+ // phaseEncodingLines in the case of interpolation or phaseResolution < 100%
+ // We'll need this for generating a value for effectiveEchoSpacing that is consistent
+ // with the *reconstructed* data.
+ int reconMatrixPE = d.phaseEncodingLines;
+ if ((h->dim[2] > 0) && (h->dim[1] > 0)) {
if (h->dim[2] == h->dim[2]) //phase encoding does not matter
- phaseEncodingLines = h->dim[2];
+ reconMatrixPE = h->dim[2];
else if (d.phaseEncodingRC =='R')
- phaseEncodingLines = h->dim[2];
+ reconMatrixPE = h->dim[2];
else if (d.phaseEncodingRC =='C')
- phaseEncodingLines = h->dim[1];
+ reconMatrixPE = h->dim[1];
}
+ if (reconMatrixPE > 0) fprintf(fp, "\t\"ReconMatrixPE\": %d,\n", reconMatrixPE );
+
+ double bandwidthPerPixelPhaseEncode = d.bandwidthPerPixelPhaseEncode;
if (bandwidthPerPixelPhaseEncode == 0.0)
bandwidthPerPixelPhaseEncode = d.CSA.bandwidthPerPixelPhaseEncode;
- if (phaseEncodingLines > 0.0) fprintf(fp, "\t\"PhaseEncodingLines\": %d,\n", phaseEncodingLines );
- if (bandwidthPerPixelPhaseEncode > 0.0)
- fprintf(fp, "\t\"BandwidthPerPixelPhaseEncode\": %g,\n", bandwidthPerPixelPhaseEncode );
+ json_Float(fp, "\t\"BandwidthPerPixelPhaseEncode\": %g,\n", bandwidthPerPixelPhaseEncode );
+ if (d.accelFactPE > 1.0) fprintf(fp, "\t\"ParallelReductionFactorInPlane\": %g,\n", d.accelFactPE);
+
+ //EffectiveEchoSpacing
+ // Siemens bandwidthPerPixelPhaseEncode already accounts for the effects of parallel imaging,
+ // interpolation, phaseOversampling, and phaseResolution, in the context of the size of the
+ // *reconstructed* data in the PE dimension
double effectiveEchoSpacing = 0.0;
- if ((phaseEncodingLines > 0) && (bandwidthPerPixelPhaseEncode > 0.0))
- effectiveEchoSpacing = 1.0 / (bandwidthPerPixelPhaseEncode * phaseEncodingLines) ;
+ if ((reconMatrixPE > 0) && (bandwidthPerPixelPhaseEncode > 0.0))
+ effectiveEchoSpacing = 1.0 / (bandwidthPerPixelPhaseEncode * reconMatrixPE);
if (d.effectiveEchoSpacingGE > 0.0)
effectiveEchoSpacing = d.effectiveEchoSpacingGE / 1000000.0;
- if (effectiveEchoSpacing > 0.0)
- fprintf(fp, "\t\"EffectiveEchoSpacing\": %g,\n", effectiveEchoSpacing);
- //FSL definition is start of first line until start of last line, so n-1 unless accelerated in-plane acquisition
- // to check: partial Fourier, iPAT, etc.
- int fencePost = 1;
- if (d.accelFactPE > 1.0)
- fencePost = (int)round(d.accelFactPE); //e.g. if 64 lines with iPAT=2, we want time from start of first until start of 62nd effective line
- if ((d.phaseEncodingSteps > 1) && (effectiveEchoSpacing > 0.0))
- fprintf(fp, "\t\"TotalReadoutTime\": %g,\n", effectiveEchoSpacing * ((float)d.phaseEncodingSteps - fencePost));
- if (d.accelFactPE > 1.0) {
- fprintf(fp, "\t\"AccelFactPE\": %g,\n", d.accelFactPE);
- if (effectiveEchoSpacing > 0.0)
- fprintf(fp, "\t\"TrueEchoSpacing\": %g,\n", effectiveEchoSpacing * d.accelFactPE);
- }
- 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 );
- }
- fprintf(fp, "\t],\n");
- }
- if (((d.phaseEncodingRC == 'R') || (d.phaseEncodingRC == 'C')) && ((d.CSA.phaseEncodingDirectionPositive == 1) || (d.CSA.phaseEncodingDirectionPositive == 0))) {
+ json_Float(fp, "\t\"EffectiveEchoSpacing\": %g,\n", effectiveEchoSpacing);
+
+ // Calculate true echo spacing (should match what Siemens reports on the console)
+ // i.e., should match "echoSpacing" extracted from the ASCII CSA header, when that exists
+ double trueESfactor = 1.0;
+ if (d.accelFactPE > 1.0) trueESfactor /= d.accelFactPE;
+ if (phaseOversampling > 0.0)
+ trueESfactor *= (1.0 + phaseOversampling);
+ float derivedEchoSpacing = 0.0;
+ derivedEchoSpacing = bandwidthPerPixelPhaseEncode * trueESfactor * reconMatrixPE;
+ if (derivedEchoSpacing != 0) derivedEchoSpacing = 1/derivedEchoSpacing;
+ json_Float(fp, "\t\"DerivedVendorReportedEchoSpacing\": %g,\n", derivedEchoSpacing);
+
+ //TotalReadOutTime: Really should be called "EffectiveReadOutTime", by analogy with "EffectiveEchoSpacing".
+ // But BIDS spec calls it "TotalReadOutTime".
+ // So, we DO NOT USE EchoTrainLength, because not trying to compute the actual (physical) readout time.
+ // Rather, the point of computing "EffectiveEchoSpacing" properly is so that this
+ // "Total(Effective)ReadOutTime" can be computed straightforwardly as the product of the
+ // EffectiveEchoSpacing and the size of the *reconstructed* matrix in the PE direction.
+ // see https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/topup/TopupUsersGuide#A--datain
+ // FSL definition is start of first line until start of last line.
+ // Other than the use of (n-1), the value is basically just 1.0/bandwidthPerPixelPhaseEncode.
+ // https://github.com/rordenlab/dcm2niix/issues/130
+ if ((reconMatrixPE > 0) && (effectiveEchoSpacing > 0.0))
+ fprintf(fp, "\t\"TotalReadoutTime\": %g,\n", effectiveEchoSpacing * (reconMatrixPE - 1.0));
+
+ if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.dwellTime > 0))
+ fprintf(fp, "\t\"DwellTime\": %g,\n", d.dwellTime * 1E-9);
+ // Phase encoding polarity
+ if (((d.phaseEncodingRC == 'R') || (d.phaseEncodingRC == 'C')) && (!d.is3DAcq) && ((d.CSA.phaseEncodingDirectionPositive == 1) || (d.CSA.phaseEncodingDirectionPositive == 0))) {
if (d.phaseEncodingRC == 'C') //Values should be "R"ow, "C"olumn or "?"Unknown
fprintf(fp, "\t\"PhaseEncodingDirection\": \"j");
else if (d.phaseEncodingRC == 'R')
@@ -706,14 +852,29 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts,
//phaseEncodingDirectionPositive has one of three values: UNKNOWN (-1), NEGATIVE (0), POSITIVE (1)
//However, DICOM and NIfTI are reversed in the j (ROW) direction
//Equivalent to dicm2nii's "if flp(iPhase), phPos = ~phPos; end"
+ //for samples see https://github.com/rordenlab/dcm2niix/issues/125
if (d.CSA.phaseEncodingDirectionPositive == -1)
fprintf(fp, "?"); //unknown
- else if ((d.CSA.phaseEncodingDirectionPositive == 1) && ((opts.isFlipY)))
+ else if ((d.CSA.phaseEncodingDirectionPositive == 0) && (d.phaseEncodingRC != 'C'))
+ fprintf(fp, "-");
+ else if ((d.phaseEncodingRC == 'C') && (d.CSA.phaseEncodingDirectionPositive == 1) && (opts.isFlipY))
fprintf(fp, "-");
- else if ((d.CSA.phaseEncodingDirectionPositive == 0) && ((!opts.isFlipY)))
+ else if ((d.phaseEncodingRC == 'C') && (d.CSA.phaseEncodingDirectionPositive == 0) && (!opts.isFlipY))
fprintf(fp, "-");
fprintf(fp, "\",\n");
} //only save PhaseEncodingDirection if BOTH direction and POLARITY are known
+ // 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 );
+ }
+ fprintf(fp, "\t],\n");
+ }
+ // Finish up with info on the conversion tool
fprintf(fp, "\t\"ConversionSoftware\": \"dcm2niix\",\n");
fprintf(fp, "\t\"ConversionSoftwareVersion\": \"%s\"\n", kDCMvers );
//fprintf(fp, "\t\"DicomConversion\": [\"dcm2niix\", \"%s\"]\n", kDCMvers );
@@ -726,34 +887,78 @@ bool isADCnotDTI(TDTI bvec) { //returns true if bval!=0 but all bvecs == 0 (Phil
((isSameFloat(bvec.V[1],0.0f)) && (isSameFloat(bvec.V[2],0.0f)) && (isSameFloat(bvec.V[3],0.0f)) ) );
}
-unsigned char * removeADC(struct nifti_1_header *hdr, unsigned char *inImg, bool * isADC) {
- int numVolOut = 0;
- int numVolIn = hdr->dim[4];
+unsigned char * removeADC(struct nifti_1_header *hdr, unsigned char *inImg, int numADC) {
+//for speed we just clip the number of volumes, the realloc routine would be nice
+// we do not want to copy input to a new smaller array since 4D DTI datasets can be huge
+// and that would require almost twice as much RAM
+ if (numADC < 1) return inImg;
+ hdr->dim[4] = hdr->dim[4] - numADC;
+ if (hdr->dim[4] < 2)
+ hdr->dim[0] = 3; //e.g. 4D 2-volume DWI+ADC becomes 3D DWI if ADC is removed
+ return inImg;
+} //removeADC()
+
+//#define naive_reorder_vols //for simple, fast re-ordering that consumes a lot of RAM
+#ifdef naive_reorder_vols
+unsigned char * reorderVolumes(struct nifti_1_header *hdr, unsigned char *inImg, int * volOrderIndex) {
+//reorder volumes to place ADC at end and (optionally) B=0 at start
+// volOrderIndex[0] reports location of desired first volume
+// naive solution creates an output buffer that doubles RAM usage (2 *numVol)
+ int numVol = hdr->dim[4];
int numVolBytes = hdr->dim[1]*hdr->dim[2]*hdr->dim[3]*(hdr->bitpix/8);
- if ((!isADC) || (numVolIn < 1) || (numVolBytes < 1)) return inImg;
- for (int i = 0; i < numVolIn; i++)
- if (!isADC[i])
- numVolOut++;
- if (numVolOut < 1) return inImg;
- //printMessage("Removing ADC maps, %d volumes reduced to %d\n", numVolIn, numVolOut);
- unsigned char *outImg = (unsigned char *)malloc(numVolBytes * numVolOut);
+ if ((!volOrderIndex) || (numVol < 1) || (numVolBytes < 1)) return inImg;
+ unsigned char *outImg = (unsigned char *)malloc(numVolBytes * numVol);
int outPos = 0;
- for (int i = 0; i < numVolIn; i++) {
- if (!isADC[i]) {
- memcpy(&outImg[outPos], &inImg[i * numVolBytes], numVolBytes); // dest, src, bytes
- outPos += numVolBytes;
- }
+ for (int i = 0; i < numVol; i++) {
+ memcpy(&outImg[outPos], &inImg[volOrderIndex[i] * numVolBytes], numVolBytes); // dest, src, bytes
+ outPos += numVolBytes;
} //for each volume
- free(isADC);
+ free(volOrderIndex);
free(inImg);
- hdr->dim[4] = numVolOut;
return outImg;
-} //removeADC()
-
-
-bool * nii_SaveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmList[], struct TDCMopts opts, int sliceDir, struct TDTI4D *dti4D) {
+} //reorderVolumes()
+#else // naive_reorder_vols
+unsigned char * reorderVolumes(struct nifti_1_header *hdr, unsigned char *inImg, int * volOrderIndex) {
+//reorder volumes to place ADC at end and (optionally) B=0 at start
+// volOrderIndex[0] reports location of desired first volume
+// complicated by fact that 4D DTI data is often huge
+// simple solutions would create an output buffer that would double RAM usage (2 *numVol)
+// here we bubble-sort volumes in place to use numVols+1 memory
+ int numVol = hdr->dim[4];
+ int numVolBytes = hdr->dim[1]*hdr->dim[2]*hdr->dim[3]*(hdr->bitpix/8);
+ int * inPos = (int *) malloc(numVol * sizeof(int));
+ for (int i = 0; i < numVol; i++)
+ inPos[i] = i;
+ unsigned char *tempVol = (unsigned char *)malloc(numVolBytes);
+ int outPos = 0;
+ for (int o = 0; o < numVol; o++) {
+ int i = inPos[volOrderIndex[o]]; //input volume
+ if (i == o) continue; //volume in correct order
+ memcpy(&tempVol[0], &inImg[o * numVolBytes], numVolBytes); //make temp
+ memcpy(&inImg[o * numVolBytes], &inImg[i * numVolBytes], numVolBytes); //copy volume to desire location dest, src, bytes
+ memcpy(&inImg[i * numVolBytes], &tempVol[0], numVolBytes); //copy unsorted volume
+ inPos[o] = i;
+ outPos += numVolBytes;
+ } //for each volume
+ free(inPos);
+ free(volOrderIndex);
+ free(tempVol);
+ return inImg;
+} //reorderVolumes()
+#endif // naive_reorder_vols
+
+float * bvals; //global variable for cmp_bvals
+int cmp_bvals(const void *a, const void *b){
+ int ia = *(int *)a;
+ int ib = *(int *)b;
+ //return bvals[ia] > bvals[ib] ? -1 : bvals[ia] < bvals[ib];
+ return bvals[ia] < bvals[ib] ? -1 : bvals[ia] > bvals[ib];
+} // cmp_bvals()
+
+int * nii_SaveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmList[], struct TDCMopts opts, int sliceDir, struct TDTI4D *dti4D, int * numADC) {
//reports non-zero if any volumes should be excluded (e.g. philip stores an ADC maps)
//to do: works with 3D mosaics and 4D files, must remove repeated volumes for 2D sequences....
+ *numADC = 0;
if (opts.isOnlyBIDS) return NULL;
uint64_t indx0 = dcmSort[0].indx; //first volume
int numDti = dcmList[indx0].CSA.numDti;
@@ -798,47 +1003,85 @@ bool * nii_SaveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],st
free(vx);
return NULL;
}
- int firstB0 = -1;
- for (int i = 0; i < numDti; i++) //check if all bvalues match first volume
- if (isSameFloat(vx[i].V[0],0) ) {
- firstB0 = i;
- break;
- }
- if (firstB0 < 0) printWarning("This diffusion series does not have a B0 (reference) volume\n");
- if (firstB0 > 0) printMessage("Note: B0 not the first volume in the series (FSL eddy reference volume is %d)\n", firstB0);
- int numADC = 0;
- bool * isADC = NULL;
- if (dcmList[dcmSort[0].indx].manufacturer == kMANUFACTURER_PHILIPS) {
- isADC = (bool *)malloc(numDti * sizeof(bool));
- int o = 0; //output index
- for (int i = 0; i < numDti; i++) {//check if all bvalues match first volume
- if (isADCnotDTI(vx[i])) {
- isADC[i] = true;
- numADC++;
- printMessage("Volume %d is not a normal DTI image (ADC?)\n", i+1);
- } else { //save output
- vx[o] = vx[i];
- isADC[i] = false;
- o++;
- }
- } //
- numDti = numDti - numADC;
- dcmList[indx0].CSA.numDti = numDti; //warning structure not changed outside scope!
- if (numADC > 0) {
- printMessage("Note: %d volumes appear to be ADC images that will be removed to allow processing\n", numADC);
- if (numDti == 0) {
- printWarning("No bvec/bval files created: no volumes with bvecs applied \n");
- free(isADC);
- free(vx);
- return NULL;
- }
+ int minB0idx = 0;
+ float minB0 = vx[0].V[0];
+ for (int i = 0; i < numDti; i++)
+ if (vx[i].V[0] < minB0) {
+ minB0 = vx[i].V[0];
+ minB0idx = i;
}
- if (numADC == 0) { //no ADC images
- free(isADC);
- isADC = NULL;
+ float maxB0 = vx[0].V[0];
+ for (int i = 0; i < numDti; i++)
+ if (vx[i].V[0] > maxB0)
+ maxB0 = vx[i].V[0];
+ //for CMRR sequences unweighted volumes are not actually B=0 but they have B near zero
+ if (minB0 > 50) printWarning("This diffusion series does not have a B0 (reference) volume\n");
+ if ((!opts.isSortDTIbyBVal) && (minB0idx > 0))
+ printMessage("Note: B0 not the first volume in the series (FSL eddy reference volume is %d)\n", minB0idx);
+
+ float kADCval = maxB0 + 1; //mark as unusual
+ *numADC = 0;
+ bvals = (float *) malloc(numDti * sizeof(float));
+ for (int i = 0; i < numDti; i++) {
+ bvals[i] = vx[i].V[0];
+ if (isADCnotDTI(vx[i])) {
+ *numADC = *numADC + 1;
+ //printMessage("Volume %d is not a normal DTI image (ADC?)\n", i+1);
+ bvals[i] = kADCval;
}
- }
- // philipsCorrectBvecs(&dcmList[indx0]); //<- replaced by unified siemensPhilips solution
+ 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);
+ }
+ //sort ALL including ADC
+ int * volOrderIndex = (int *) malloc(numDti * sizeof(int));
+ for (int i = 0; i < numDti; i++)
+ volOrderIndex[i] = i;
+ if (opts.isSortDTIbyBVal)
+ qsort(volOrderIndex, numDti, sizeof(*volOrderIndex), cmp_bvals);
+ else if (*numADC > 0) {
+ int o = 0;
+ for (int i = 0; i < numDti; i++) {
+ if (bvals[i] < kADCval) {
+ volOrderIndex[o] = i;
+ o++;
+ } //if not ADC
+ } //for each volume
+ } //if sort else if has ADC
+ free(bvals);
+ //save VX as sorted
+ TDTI * vxOrig = (TDTI *)malloc(numDti * sizeof(TDTI));
+ for (int i = 0; i < numDti; i++)
+ vxOrig[i] = vx[i];
+ //remove ADC
+ numDti = numDti - *numADC;
+ free(vx);
+ vx = (TDTI *)malloc(numDti * sizeof(TDTI));
+ for (int i = 0; i < numDti; i++)
+ vx[i] = vxOrig[volOrderIndex[i]];
+ free(vxOrig);
+ //if no ADC or sequential, the is no need to re-order volumes
+ bool isSequential = true;
+ for (int i = 1; i < (numDti + *numADC); i++)
+ if (volOrderIndex[i] <= volOrderIndex[i-1])
+ isSequential = false;
+ if (isSequential) {
+ free(volOrderIndex);
+ volOrderIndex = NULL;
+ }
+ if (!isSequential)
+ printMessage("DTI volumes re-ordered by ascending b-value\n");
+ dcmList[indx0].CSA.numDti = numDti; //warning structure not changed outside scope!
+
geCorrectBvecs(&dcmList[indx0],sliceDir, vx);
siemensPhilipsCorrectBvecs(&dcmList[indx0],sliceDir, vx);
if (!opts.isFlipY ) { //!FLIP_Y&& (dcmList[indx0].CSA.mosaicSlices < 2) mosaics are always flipped in the Y direction
@@ -876,7 +1119,7 @@ bool * nii_SaveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],st
FILE *fp = fopen(txtname, "w");
if (fp == NULL) {
free(vx);
- return isADC;
+ return volOrderIndex;
}
for (int i = 0; i < (numDti-1); i++) {
if (opts.isCreateBIDS) {
@@ -893,7 +1136,7 @@ bool * nii_SaveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],st
fp = fopen(txtname, "w");
if (fp == NULL) {
free(vx);
- return isADC;
+ return volOrderIndex;
}
for (int v = 1; v < 4; v++) {
for (int i = 0; i < (numDti-1); i++) {
@@ -908,7 +1151,7 @@ bool * nii_SaveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],st
fclose(fp);
#endif
free(vx);
- return isADC;
+ return volOrderIndex;
}// nii_SaveDTI()
float sqr(float v){
@@ -1029,8 +1272,8 @@ int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopt
if (strlen(inname) < 1) {
strcpy(inname, "T%t_N%n_S%s");
}
- int start = 0;
- int pos = 0;
+ size_t start = 0;
+ size_t pos = 0;
bool isCoilReported = false;
bool isEchoReported = false;
bool isSeriesReported = false;
@@ -1049,10 +1292,8 @@ int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopt
sprintf(newstr, "%02d", dcm.coilNum);
strcat (outname,newstr);
}
- if (f == 'C')
- strcat (outname,dcm.imageComments);
- if (f == 'D')
- strcat (outname,dcm.seriesDescription);
+ if (f == 'C') strcat (outname,dcm.imageComments);
+ if (f == 'D') strcat (outname,dcm.seriesDescription);
if (f == 'E') {
isEchoReported = true;
sprintf(newstr, "%d", dcm.echoNum);
@@ -1157,19 +1398,51 @@ int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopt
if (strlen(outname) < 1) strcpy(outname, "dcm2nii_invalidName");
if (outname[0] == '.') outname[0] = '_'; //make sure not a hidden file
//eliminate illegal characters http://msdn.microsoft.com/en-us/library/windows/desktop/aa365247(v=vs.85).aspx
- for (int pos = 0; pos<strlen(outname); pos ++)
+ #if defined(_WIN64) || defined(_WIN32) //https://stackoverflow.com/questions/1976007/what-characters-are-forbidden-in-windows-and-linux-directory-names
+ for (size_t pos = 0; pos<strlen(outname); pos ++)
if ((outname[pos] == '<') || (outname[pos] == '>') || (outname[pos] == ':')
- || (outname[pos] == '"') || (outname[pos] == '\\') || (outname[pos] == '/')
+ || (outname[pos] == '"') // || (outname[pos] == '/') || (outname[pos] == '\\')
|| (outname[pos] == '^')
|| (outname[pos] == '*') || (outname[pos] == '|') || (outname[pos] == '?'))
outname[pos] = '_';
- //printMessage("outname=*%s* %d %d\n", outname, pos,start);
+ for (int pos = 0; pos<strlen(outname); pos ++)
+ if (outname[pos] == '/')
+ outname[pos] = kPathSeparator; //for Windows, convert "/" to "\"
+ #else
+ for (size_t pos = 0; pos<strlen(outname); pos ++)
+ if (outname[pos] == ':') //not allowed by MacOS
+ outname[pos] = '_';
+ #endif
char baseoutname[2048] = {""};
strcat (baseoutname,pth);
char appendChar[2] = {"a"};
appendChar[0] = kPathSeparator;
- if (pth[strlen(pth)-1] != kPathSeparator)
+ if ((pth[strlen(pth)-1] != kPathSeparator) && (outname[0] != kPathSeparator))
strcat (baseoutname,appendChar);
+ //Allow user to specify new folders, e.g. "-f dir/%p" or "-f %s/%p/%m"
+ // These folders are created if they do not exist
+ char *sep = strchr(outname, kPathSeparator);
+ if (sep) {
+ char newdir[2048] = {""};
+ strcat (newdir,baseoutname);
+ //struct stat st = {0};
+ for (size_t pos = 0; pos< strlen(outname); pos ++) {
+ if (outname[pos] == kPathSeparator) {
+ //if (stat(newdir, &st) == -1)
+ if (!is_dir(newdir,true))
+ #if defined(_WIN64) || defined(_WIN32)
+ mkdir(newdir);
+ #else
+ mkdir(newdir, 0700);
+ #endif
+ }
+ char ch[12] = {""};
+ sprintf(ch,"%c",outname[pos]);
+ strcat (newdir,ch);
+ }
+ }
+ //printMessage("path='%s' name='%s'\n", pathoutname, outname);
+ //make sure outname is unique
strcat (baseoutname,outname);
char pathoutname[2048] = {""};
strcat (pathoutname,baseoutname);
@@ -1236,6 +1509,10 @@ unsigned long mz_crc32(unsigned long crc, const unsigned char *ptr, size_t buf_l
#define MZ_UBER_COMPRESSION 9
#endif
+#ifndef MZ_DEFAULT_LEVEL
+ #define MZ_DEFAULT_LEVEL 6
+#endif
+
void writeNiiGz (char * baseName, struct nifti_1_header hdr, unsigned char* src_buffer, unsigned long src_len, int gzLevel) {
//create gz file in RAM, save to disk http://www.zlib.net/zlib_how.html
// in general this single-threaded approach is slower than PIGZ but is useful for slow (network attached) disk drives
@@ -1253,7 +1530,7 @@ void writeNiiGz (char * baseName, struct nifti_1_header hdr, unsigned char* src
strm.opaque = Z_NULL;
strm.next_out = pCmp; // output char array
strm.avail_out = (unsigned int)cmp_len; // size of output
- int zLevel = Z_DEFAULT_COMPRESSION;
+ int zLevel = MZ_DEFAULT_LEVEL;//Z_DEFAULT_COMPRESSION;
if ((gzLevel > 0) && (gzLevel < 11))
zLevel = gzLevel;
if (zLevel > MZ_UBER_COMPRESSION)
@@ -1401,16 +1678,29 @@ int nii_saveNII(char * niiFilename, struct nifti_1_header hdr, unsigned char* im
printMessage("Error: Image size is zero bytes %s\n", niiFilename);
return EXIT_FAILURE;
}
- #define kMaxPigz 3758096384
- #ifndef myDisableZLib
- if ((opts.isGz) && (strlen(opts.pigzname) < 1) && ((imgsz+hdr.vox_offset) >= 2147483647) ) { //use internal compressor
- printWarning("Saving uncompressed data: internal compressor limited to 2Gb images.\n");
- if ((imgsz+hdr.vox_offset) < 3758096384)
- printWarning(" Hint: using external compressor (pigz) should help.\n");
- } else if ((opts.isGz) && (strlen(opts.pigzname) < 1) && ((imgsz+hdr.vox_offset) < 2147483647) ) { //use internal compressor
- writeNiiGz (niiFilename, hdr, im, imgsz, opts.gzLevel);
- return EXIT_SUCCESS;
- }
+ #ifndef myDisableGzSizeLimits
+ //see https://github.com/rordenlab/dcm2niix/issues/124
+ uint64_t kMaxPigz = 4294967264;
+ //https://stackoverflow.com/questions/5272825/detecting-64bit-compile-in-c
+ #ifndef UINTPTR_MAX
+ uint64_t kMaxGz = 2147483647;
+ #elif UINTPTR_MAX == 0xffffffff
+ uint64_t kMaxGz = 2147483647;
+ #elif UINTPTR_MAX == 0xffffffffffffffff
+ uint64_t kMaxGz = kMaxPigz;
+ #else
+ compiler error: unable to determine is 32 or 64 bit
+ #endif
+ #ifndef myDisableZLib
+ if ((opts.isGz) && (strlen(opts.pigzname) < 1) && ((imgsz+hdr.vox_offset) >= kMaxGz) ) { //use internal compressor
+ printWarning("Saving uncompressed data: internal compressor unable to process such large files.\n");
+ if ((imgsz+hdr.vox_offset) < kMaxPigz)
+ printWarning(" Hint: using external compressor (pigz) should help.\n");
+ } else if ((opts.isGz) && (strlen(opts.pigzname) < 1) && ((imgsz+hdr.vox_offset) < kMaxGz) ) { //use internal compressor
+ writeNiiGz (niiFilename, hdr, im, imgsz, opts.gzLevel);
+ return EXIT_SUCCESS;
+ }
+ #endif
#endif
char fname[2048] = {""};
strcpy (fname,niiFilename);
@@ -1423,10 +1713,12 @@ int nii_saveNII(char * niiFilename, struct nifti_1_header hdr, unsigned char* im
fwrite(&im[0], imgsz, 1, fp);
fclose(fp);
if ((opts.isGz) && (strlen(opts.pigzname) > 0) ) {
+ #ifndef myDisableGzSizeLimits
if ((imgsz+hdr.vox_offset) > kMaxPigz) {
printWarning("Saving uncompressed data: image too large for pigz.\n");
return EXIT_SUCCESS;
}
+ #endif
char command[768];
strcpy(command, "\"" );
strcat(command, opts.pigzname );
@@ -1879,12 +2171,48 @@ int nii_saveCrop(char * niiFilename, struct nifti_1_header hdr, unsigned char* i
return returnCode;
}// nii_saveCrop()
+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
+ 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 ((minT != maxT) && (maxT <= d->TR)) return; //looks fine
+ if ((minT == maxT) && (d->CSA.multiBandFactor == d->CSA.mosaicSlices)) return; //fine: all slices single excitation
+ if ((strlen(d->seriesDescription) > 0) && (strstr(d->seriesDescription, "SBRef") != NULL)) return; //fine: single-band calibration data, the slice timing WILL exceed the TR
+ //check if 2nd image has valud slice timing
+ float minT1 = d1->CSA.sliceTiming[0];
+ float maxT1 = minT1;
+ for (int i = 0; i < kMaxEPI3D; i++) {
+ if (d1->CSA.sliceTiming[i] < 0.0) break;
+ if (d1->CSA.sliceTiming[i] < minT1) minT1 = d1->CSA.sliceTiming[i];
+ if (d1->CSA.sliceTiming[i] > maxT1) maxT1 = d1->CSA.sliceTiming[i];
+ }
+ if ((minT1 == maxT1) || (maxT1 >= d->TR)) { //both first and second image corrupted
+ printWarning("CSA slice timing appears corrupted (range %g..%g, TR=%gms)\n", minT, maxT, d->TR);
+ return;
+ }
+ //1st image corrupted, but 2nd looks ok - substitute values from 2nd image
+ for (int i = 0; i < kMaxEPI3D; i++) {
+ d->CSA.sliceTiming[i] = d1->CSA.sliceTiming[i];
+ if (d1->CSA.sliceTiming[i] < 0.0) break;
+ }
+ d->CSA.multiBandFactor = d1->CSA.multiBandFactor;
+ printMessage("CSA slice timing based on 2nd volume, 1st volume corrupted (CMRR bug, range %g..%g, TR=%gms)\n", minT, maxT, d->TR);
+}//checkSliceTiming
+
int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmList[], struct TSearchList *nameList, struct TDCMopts opts, struct TDTI4D *dti4D) {
bool iVaries = intensityScaleVaries(nConvert,dcmSort,dcmList);
float *sliceMMarray = NULL; //only used if slices are not equidistant
uint64_t indx = dcmSort[0].indx;
uint64_t indx0 = dcmSort[0].indx;
- if (opts.isIgnoreDerivedAnd2D && isDerived(dcmList[indx])) {
+ uint64_t indx1 = indx0;
+ if (nConvert > 1) indx1 = dcmSort[1].indx;
+ if (opts.isIgnoreDerivedAnd2D && dcmList[indx].isDerived) {
printMessage("Ignoring derived image(s) of series %ld %s\n", dcmList[indx].seriesNum, nameList->str[indx]);
return EXIT_SUCCESS;
}
@@ -1899,6 +2227,10 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis
bool saveAs3D = dcmList[indx].isHasPhase;
struct nifti_1_header hdr0;
unsigned char * img = nii_loadImgXL(nameList->str[indx], &hdr0,dcmList[indx], iVaries, opts.compressFlag, opts.isVerbose);
+ if (strlen(opts.imageComments) > 0) {
+ for (int i = 0; i < 24; i++) hdr0.aux_file[i] = 0; //remove dcm.imageComments
+ snprintf(hdr0.aux_file,24,"%s",opts.imageComments);
+ }
if (opts.isVerbose)
printMessage("Converting %s\n",nameList->str[indx]);
if (img == NULL) return EXIT_FAILURE;
@@ -2035,7 +2367,9 @@ 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!
}
- nii_SaveBIDS(pathoutname, dcmList[dcmSort[0].indx], opts, dti4D, &hdr0, nameList->str[dcmSort[0].indx]);
+ 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]);
if (opts.isOnlyBIDS) {
//note we waste time loading every image, however this ensures hdr0 matches actual output
free(imgM);
@@ -2043,11 +2377,11 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis
}
nii_SaveText(pathoutname, dcmList[dcmSort[0].indx], opts, &hdr0, nameList->str[indx]);
- bool * isADC = nii_SaveDTI(pathoutname,nConvert, dcmSort, dcmList, opts, sliceDir, dti4D);
+ int numADC = 0;
+ int * volOrderIndex = nii_SaveDTI(pathoutname,nConvert, dcmSort, dcmList, opts, sliceDir, dti4D, &numADC);
if ((hdr0.datatype == DT_UINT16) && (!dcmList[dcmSort[0].indx].isSigned)) nii_check16bitUnsigned(imgM, &hdr0);
printMessage( "Convert %d DICOM as %s (%dx%dx%dx%d)\n", nConvert, pathoutname, hdr0.dim[1],hdr0.dim[2],hdr0.dim[3],hdr0.dim[4]);
PhilipsPrecise(&dcmList[dcmSort[0].indx], opts.isPhilipsFloatNotDisplayScaling, &hdr0);
-
if (!dcmList[dcmSort[0].indx].isSlicesSpatiallySequentialPhilips)
nii_reorderSlices(imgM, &hdr0, dti4D);
if (hdr0.dim[3] < 2)
@@ -2071,18 +2405,25 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis
if ((hdr0.dim[4] > 1) && (saveAs3D))
returnCode = nii_saveNII3D(pathoutname, hdr0, imgM,opts);
else {
- if (isADC) { //ADC maps can disrupt analysis: save a copy with the ADC map, and another without
+ if (volOrderIndex) //reorder volumes
+ imgM = reorderVolumes(&hdr0, imgM, volOrderIndex);
#ifndef HAVE_R
- char pathoutnameADC[2048] = {""};
- strcat(pathoutnameADC,pathoutname);
- strcat(pathoutnameADC,"_ADC");
- nii_saveNII(pathoutnameADC, hdr0, imgM, opts);
+ if (numADC > 0) {//ADC maps can disrupt analysis: save a copy with the ADC map, and another without
+ char pathoutnameADC[2048] = {""};
+ strcat(pathoutnameADC,pathoutname);
+ strcat(pathoutnameADC,"_ADC");
+ if (opts.isSave3D)
+ nii_saveNII3D(pathoutnameADC, hdr0, imgM, opts);
+ else
+ nii_saveNII(pathoutnameADC, hdr0, imgM, opts);
+ }
#endif
- imgM = removeADC(&hdr0, imgM, isADC);
- //hdr0.dim[4] = hdr0.dim[4]-numFinalADC;
- };
+ imgM = removeADC(&hdr0, imgM, numADC);
#ifndef HAVE_R
- returnCode = nii_saveNII(pathoutname, hdr0, imgM, opts);
+ if (opts.isSave3D)
+ returnCode = nii_saveNII3D(pathoutname, hdr0, imgM, opts);
+ else
+ returnCode = nii_saveNII(pathoutname, hdr0, imgM, opts);
#endif
}
#endif
@@ -2114,7 +2455,7 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis
#endif
free(imgM);
- return EXIT_SUCCESS;
+ return returnCode;//EXIT_SUCCESS;
}// saveDcm2Nii()
int compareTDCMsort(void const *item1, void const *item2) {
@@ -2156,11 +2497,14 @@ TWarnings setWarnings() {
return r;
}
-bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, bool isForceStackSameSeries, struct TWarnings* warnings, bool *isMultiEcho) {
+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!
+ if (d1.isDerived != d2.isDerived) return false; //do not stack raw and derived image types
+ if (d1.manufacturer != d2.manufacturer) return false; //do not stack data from different vendors
if (d1.seriesNum != d2.seriesNum) return false;
#ifdef mySegmentByAcq
if (d1.acquNum != d2.acquNum) return false;
@@ -2177,7 +2521,7 @@ bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, bool isForceStackSam
warnings->bitDepthVaries = true;
return false;
}
- if (isForceStackSameSeries) {
+ if (opts->isForceStackSameSeries) {
if ((d1.TE != d2.TE) || (d1.echoNum != d2.echoNum))
*isMultiEcho = true;
return true; //we will stack these images, even if they differ in the following attributes
@@ -2190,9 +2534,9 @@ bool isSameSet (struct TDICOMdata d1, struct TDICOMdata d2, bool isForceStackSam
}
if ((d1.TE != d2.TE) || (d1.echoNum != d2.echoNum)) {
if ((!warnings->echoVaries) && (d1.isXRay)) //for CT/XRay we check DICOM tag 0018,1152 (XRayExposure)
- printMessage("slices not stacked: X-Ray Exposure varies (%g, %g; number %d, %d)\n", d1.TE, d2.TE,d1.echoNum, d2.echoNum );
+ printMessage("slices not stacked: X-Ray Exposure varies (exposure %g, %g; number %d, %d). Use 'merge 2D slices' option to force stacking\n", d1.TE, d2.TE,d1.echoNum, d2.echoNum );
if ((!warnings->echoVaries) && (!d1.isXRay)) //for MRI
- printMessage("slices not stacked: echo varies (TE %g, %g; echo %d, %d)\n", d1.TE, d2.TE,d1.echoNum, d2.echoNum );
+ printMessage("slices not stacked: echo varies (TE %g, %g; echo %d, %d). Use 'merge 2D slices' option to force stacking\n", d1.TE, d2.TE,d1.echoNum, d2.echoNum );
warnings->echoVaries = true;
*isMultiEcho = true;
return false;
@@ -2364,7 +2708,7 @@ int convert_parRec(struct TDCMopts opts) {
if (nameList.numItems < 1)
printMessage("No valid PAR/REC files were found\n");
if (nameList.numItems > 0)
- for (int i = 0; i < nameList.numItems; i++)
+ for (int i = 0; i < (int)nameList.numItems; i++)
free(nameList.str[i]);
free(nameList.str);
return ret;
@@ -2459,7 +2803,7 @@ int nii_loadDir(struct TDCMopts* opts) {
int nConvertTotal = 0;
bool compressionWarning = false;
bool convertError = false;
- for (int i = 0; i < nDcm; i++ ) {
+ for (int i = 0; i < (int)nDcm; i++ ) {
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];
@@ -2491,7 +2835,7 @@ int nii_loadDir(struct TDCMopts* opts) {
// If the file matches an existing series, add it to the corresponding file list
for (int j = 0; j < opts->series.size(); j++) {
bool isMultiEchoUnused;
- if (isSameSet(opts->series[j].representativeData, dcmList[i], opts->isForceStackSameSeries, &warnings, &isMultiEchoUnused)) {
+ if (isSameSet(opts->series[j].representativeData, dcmList[i], opts, &warnings, &isMultiEchoUnused)) {
opts->series[j].files.push_back(nameList.str[i]);
matched = true;
break;
@@ -2510,19 +2854,19 @@ int nii_loadDir(struct TDCMopts* opts) {
} else {
#endif
//3: stack DICOMs with the same Series
- for (int i = 0; i < nDcm; i++ ) {
+ for (int i = 0; i < (int)nDcm; i++ ) {
if ((dcmList[i].converted2NII == 0) && (dcmList[i].isValid)) {
int nConvert = 0;
struct TWarnings warnings = setWarnings();
bool isMultiEcho;
- for (int j = i; j < nDcm; j++)
- if (isSameSet(dcmList[i], dcmList[j], opts->isForceStackSameSeries, &warnings, &isMultiEcho ) )
+ 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;
- for (int j = i; j < nDcm; j++)
- if (isSameSet(dcmList[i], dcmList[j], opts->isForceStackSameSeries, &warnings, &isMultiEcho)) {
+ 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;
@@ -2585,40 +2929,40 @@ int nii_loadDir(struct TDCMopts* opts) {
#if defined(_WIN64) || defined(_WIN32)
#else //UNIX
-#define PATH_MAX 1024
int findpathof(char *pth, const char *exe) {
//Find executable by searching the PATH environment variable.
// http://www.linuxquestions.org/questions/programming-9/get-full-path-of-a-command-in-c-117965/
- char *searchpath;
- char *beg, *end;
- int stop, found;
- int len;
- if (strchr(exe, '/') != NULL) {
- if (realpath(exe, pth) == NULL) return 0;
- return is_exe(pth);
- }
- searchpath = getenv("PATH");
- if (searchpath == NULL) return 0;
- if (strlen(searchpath) <= 0) return 0;
- beg = searchpath;
- stop = 0; found = 0;
- do {
- end = strchr(beg, ':');
- if (end == NULL) {
- stop = 1;
- strncpy(pth, beg, PATH_MAX);
- len = strlen(pth);
- } else {
- strncpy(pth, beg, end - beg);
- pth[end - beg] = '\0';
- len = end - beg;
- }
- if (pth[len - 1] != '/') strncat(pth, "/", 1);
- strncat(pth, exe, PATH_MAX - len);
- found = is_exe(pth);
- if (!stop) beg = end + 1;
- } while (!stop && !found);
- return found;
+ char *searchpath;
+ char *beg, *end;
+ int stop, found;
+ size_t len;
+ if (strchr(exe, '/') != NULL) {
+ if (realpath(exe, pth) == NULL) return 0;
+ return is_exe(pth);
+ }
+ searchpath = getenv("PATH");
+ if (searchpath == NULL) return 0;
+ if (strlen(searchpath) <= 0) return 0;
+ beg = searchpath;
+ stop = 0; found = 0;
+ do {
+ end = strchr(beg, ':');
+ if (end == NULL) {
+ len = strlen(beg);
+ if (len == 0) return 0;
+ strncpy(pth, beg, len);
+ stop = 1;
+ } else {
+ strncpy(pth, beg, end - beg);
+ pth[end - beg] = '\0';
+ len = end - beg;
+ }
+ if (pth[len - 1] != '/') strncat(pth, "/", 1);
+ strncat(pth, exe, PATH_MAX - len);
+ found = is_exe(pth);
+ if (!stop) beg = end + 1;
+ } while (!stop && !found);
+ return found;
}
#endif
@@ -2643,7 +2987,7 @@ void readFindPigz (struct TDCMopts *opts, const char * argv[]) {
"pigz_afni",
};
#define n_nam (sizeof (nams) / sizeof (const char *))
- for (int n = 0; n < n_nam; n++) {
+ for (int n = 0; n < (int)n_nam; n++) {
if (findpathof(str, nams[n])) {
strcpy(opts->pigzname,str);
//printMessage("Found pigz: %s\n", str);
@@ -2663,12 +3007,11 @@ void readFindPigz (struct TDCMopts *opts, const char * argv[]) {
appendChar[0] = kPathSeparator;
if (exepth[strlen(exepth)-1] != kPathSeparator) strcat (exepth,appendChar);
//see if pigz in any path
- for (int n = 0; n < n_nam; n++) {
+ for (int n = 0; n < (int)n_nam; n++) {
//printf ("%d: %s\n", i, nams[n]);
- for (int p = 0; p < n_pth; p++) {
+ for (int p = 0; p < (int)n_pth; p++) {
strcpy(str, pths[p]);
strcat(str, nams[n]);
- //printf (">>>%s\n", str);
if (is_exe(str))
goto pigzFound;
} //p
@@ -2706,17 +3049,20 @@ void setDefaultOpts (struct TDCMopts *opts, const char * argv[]) { //either "set
//printMessage("%d %s\n",opts->compressFlag, opts->compressname);
strcpy(opts->indir,"");
strcpy(opts->outdir,"");
+ strcpy(opts->imageComments,"");
opts->isOnlySingleFile = false; //convert all files in a directory, not just a single file
opts->isForceStackSameSeries = false;
opts->isIgnoreDerivedAnd2D = false;
opts->isPhilipsFloatNotDisplayScaling = true;
opts->isCrop = false;
opts->isGz = false;
- opts->gzLevel = -1;
+ opts->isSave3D = false;
+ opts->gzLevel = MZ_DEFAULT_LEVEL; //-1;
opts->isFlipY = true; //false: images in raw DICOM orientation, true: image rows flipped to cartesian coordinates
opts->isRGBplanar = false; //false for NIfTI (RGBRGB...), true for Analyze (RRR..RGGG..GBBB..B)
opts->isCreateBIDS = true;
opts->isOnlyBIDS = false;
+ opts->isSortDTIbyBVal = false;
#ifdef myNoAnonymizeBIDS
opts->isAnonymizeBIDS = false;
#else
@@ -2805,4 +3151,4 @@ void saveIniFile (struct TDCMopts opts) {
fclose(fp);
} //saveIniFile()
-#endif
+#endif
\ No newline at end of file
diff --git a/console/nii_dicom_batch.h b/console/nii_dicom_batch.h
index 533b6fa..f9f0ba9 100644
--- a/console/nii_dicom_batch.h
+++ b/console/nii_dicom_batch.h
@@ -23,9 +23,9 @@ extern "C" {
#endif
struct TDCMopts {
- bool isGz, isFlipY, isCreateBIDS, isAnonymizeBIDS, isOnlyBIDS, isCreateText, isIgnoreDerivedAnd2D, isPhilipsFloatNotDisplayScaling, isTiltCorrect, isRGBplanar, isOnlySingleFile, isForceStackSameSeries, isCrop;
+ 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];
+ char filename[512], outdir[512], indir[512], pigzname[512], optsname[512], indirParent[512], imageComments[24];
#ifdef HAVE_R
bool isScanOnly;
void *imageList;
@@ -38,7 +38,8 @@ extern "C" {
int nii_saveNII(char * niiFilename, struct nifti_1_header hdr, unsigned char* im, struct TDCMopts opts);
//void readIniFile (struct TDCMopts *opts);
int nii_loadDir (struct TDCMopts *opts);
- void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, struct TDTI4D *dti4D, struct nifti_1_header *h, const char * filename);
+ //void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, struct TDTI4D *dti4D, struct nifti_1_header *h, const char * filename);
+ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, struct nifti_1_header *h, const char * filename);
int nii_createFilename(struct TDICOMdata dcm, char * niiFilename, struct TDCMopts opts);
void nii_createDummyFilename(char * niiFilename, struct TDCMopts opts);
//void findExe(char name[512], const char * argv[]);
diff --git a/console/nii_foreign.cpp b/console/nii_foreign.cpp
index 13989be..b1666d6 100644
--- a/console/nii_foreign.cpp
+++ b/console/nii_foreign.cpp
@@ -399,8 +399,9 @@ int convert_foreign (const char *fn, struct TDCMopts opts){
int ret = nii_createFilename(dcm, niiFilename, opts);
printMessage("Saving ECAT as '%s'\n", niiFilename);
if (ret != EXIT_SUCCESS) return ret;
- struct TDTI4D dti4D;
- nii_SaveBIDS(niiFilename, dcm, opts, &dti4D, &hdr, fn);
+ //struct TDTI4D dti4D;
+ //nii_SaveBIDS(niiFilename, dcm, opts, &dti4D, &hdr, fn);
+ nii_SaveBIDS(niiFilename, dcm, opts, &hdr, fn);
ret = nii_saveNII(niiFilename, hdr, img, opts);
free(img);
return ret;
diff --git a/console/nii_ortho.cpp b/console/nii_ortho.cpp
index 6f72ce4..832f3a6 100644
--- a/console/nii_ortho.cpp
+++ b/console/nii_ortho.cpp
@@ -111,7 +111,7 @@ vec3i setOrientVec(mat33 m)
// Assumes isOrthoMat NOT computed on INVERSE, hence return INVERSE of solution...
//e.g. [-1,2,3] means reflect x axis, [2,1,3] means swap x and y dimensions
{
- vec3i ret = {0, 0, 0};
+ vec3i ret = {{0, 0, 0}};
//mat33 m = {-1,0,0, 0,1,0, 0,0,1};
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) {
@@ -223,8 +223,8 @@ unsigned char * reOrient(unsigned char* img, struct nifti_1_header *h, vec3i or
{
size_t nvox = h->dim[1] * h->dim[2] * h->dim[3];
if (nvox < 1) return img;
- vec3i outDim= {0,0,0};
- vec3i outInc= {0,0,0};
+ vec3i outDim= {{0,0,0}};
+ vec3i outInc= {{0,0,0}};
for (int i = 0; i < 3; i++) { //set dimension, pixdim and
outDim.v[i] = h->dim[abs(orientVec.v[i])];
if (abs(orientVec.v[i]) == 1) outInc.v[i] = 1;
@@ -239,7 +239,7 @@ unsigned char * reOrient(unsigned char* img, struct nifti_1_header *h, vec3i or
}
reOrientImg(img, outDim, outInc, h->bitpix / 8, nvol);
//now change the header....
- vec3 outPix= {h->pixdim[abs(orientVec.v[0])],h->pixdim[abs(orientVec.v[1])],h->pixdim[abs(orientVec.v[2])]};
+ vec3 outPix= {{h->pixdim[abs(orientVec.v[0])],h->pixdim[abs(orientVec.v[1])],h->pixdim[abs(orientVec.v[2])]}};
for (int i = 0; i < 3; i++) {
h->dim[i+1] = outDim.v[i];
h->pixdim[i+1] = outPix.v[i];
diff --git a/console/ucm.cmake b/console/ucm.cmake
index 2ab74ea..e6e7fc8 100644
--- a/console/ucm.cmake
+++ b/console/ucm.cmake
@@ -19,7 +19,7 @@ if(NOT COMMAND cotire)
include(${CMAKE_CURRENT_LIST_DIR}/../cotire/CMake/cotire.cmake OPTIONAL)
endif()
-if(COMMAND cotire AND "1.7.7" VERSION_LESS "${COTIRE_CMAKE_MODULE_VERSION}")
+if(COMMAND cotire AND "1.7.9" VERSION_LESS "${COTIRE_CMAKE_MODULE_VERSION}")
set(ucm_with_cotire 1)
else()
set(ucm_with_cotire 0)
@@ -87,7 +87,7 @@ macro(ucm_add_linker_flags)
endif()
foreach(CONFIG ${ARG_CONFIG})
- string(TOUPPER "${ARG_CONFIG}" ARG_CONFIG)
+ string(TOUPPER "${CONFIG}" CONFIG)
if(NOT ${ARG_EXE} AND NOT ${ARG_MODULE} AND NOT ${ARG_SHARED} AND NOT ${ARG_STATIC})
set(ARG_EXE 1)
@@ -210,8 +210,10 @@ macro(ucm_set_runtime)
if("${ARG_STATIC}")
foreach(flags ${flags_configs})
if(CMAKE_CXX_COMPILER_ID MATCHES "GNU")
- if(NOT ${flags} MATCHES "-static-libstdc\\+\\+")
- set(${flags} "${${flags}} -static-libstdc++")
+ if (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER 4.4.7) # option "-static-libstdc++" available since GCC 4.5
+ if(NOT ${flags} MATCHES "-static-libstdc\\+\\+")
+ set(${flags} "${${flags}} -static-libstdc++")
+ endif()
endif()
if(NOT ${flags} MATCHES "-static-libgcc")
set(${flags} "${${flags}} -static-libgcc")
diff --git a/old_nii_SaveBIDS.cpp b/old_nii_SaveBIDS.cpp
new file mode 100644
index 0000000..a055e18
--- /dev/null
+++ b/old_nii_SaveBIDS.cpp
@@ -0,0 +1,274 @@
+void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts, struct nifti_1_header *h, const char * filename) {
+//https://docs.google.com/document/d/1HFUkAEE-pB-angVcYe6pf_-fVf4sCpOHKesUvfb8Grc/edit#
+// Generate Brain Imaging Data Structure (BIDS) info
+// sidecar JSON file (with the same filename as the .nii.gz file, but with .json extension).
+// we will use %g for floats since exponents are allowed
+// we will not set the locale, so decimal separator is always a period, as required
+// https://www.ietf.org/rfc/rfc4627.txt
+ if (!opts.isCreateBIDS) return;
+ char txtname[2048] = {""};
+ strcpy (txtname,pathoutname);
+ strcat (txtname,".json");
+ FILE *fp = fopen(txtname, "w");
+ fprintf(fp, "{\n");
+ switch (d.modality) {
+ case kMODALITY_CR:
+ fprintf(fp, "\t\"Modality\": \"CR\",\n" );
+ break;
+ case kMODALITY_CT:
+ fprintf(fp, "\t\"Modality\": \"CT\",\n" );
+ break;
+ case kMODALITY_MR:
+ fprintf(fp, "\t\"Modality\": \"MR\",\n" );
+ break;
+ case kMODALITY_PT:
+ fprintf(fp, "\t\"Modality\": \"PT\",\n" );
+ break;
+ case kMODALITY_US:
+ fprintf(fp, "\t\"Modality\": \"US\",\n" );
+ break;
+ };
+ switch (d.manufacturer) {
+ case kMANUFACTURER_SIEMENS:
+ fprintf(fp, "\t\"Manufacturer\": \"Siemens\",\n" );
+ break;
+ case kMANUFACTURER_GE:
+ fprintf(fp, "\t\"Manufacturer\": \"GE\",\n" );
+ break;
+ case kMANUFACTURER_PHILIPS:
+ fprintf(fp, "\t\"Manufacturer\": \"Philips\",\n" );
+ break;
+ case kMANUFACTURER_TOSHIBA:
+ fprintf(fp, "\t\"Manufacturer\": \"Toshiba\",\n" );
+ break;
+ };
+ fprintf(fp, "\t\"ManufacturersModelName\": \"%s\",\n", d.manufacturersModelName );
+ if (!opts.isAnonymizeBIDS) {
+ if (strlen(d.seriesInstanceUID) > 0)
+ fprintf(fp, "\t\"SeriesInstanceUID\": \"%s\",\n", d.seriesInstanceUID );
+ if (strlen(d.studyInstanceUID) > 0)
+ fprintf(fp, "\t\"StudyInstanceUID\": \"%s\",\n", d.studyInstanceUID );
+ if (strlen(d.referringPhysicianName) > 0)
+ fprintf(fp, "\t\"ReferringPhysicianName\": \"%s\",\n", d.referringPhysicianName );
+ if (strlen(d.studyID) > 0)
+ fprintf(fp, "\t\"StudyID\": \"%s\",\n", d.studyID );
+ //Next lines directly reveal patient identity
+ //if (strlen(d.patientName) > 0)
+ // fprintf(fp, "\t\"PatientName\": \"%s\",\n", d.patientName );
+ //if (strlen(d.patientID) > 0)
+ // fprintf(fp, "\t\"PatientID\": \"%s\",\n", d.patientID );
+ }
+ #ifdef myReadAsciiCsa
+ if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.CSA.SeriesHeader_offset > 0) && (d.CSA.SeriesHeader_length > 0)) {
+ //&& (strlen(d.scanningSequence) > 1) && (d.scanningSequence[0] == 'E') && (d.scanningSequence[1] == 'P')) { //for EPI scans only
+ int partialFourier, echoSpacing, echoTrainDuration, epiFactor, parallelReductionFactorInPlane;
+ char fmriExternalInfo[kDICOMStr], coilID[kDICOMStr], consistencyInfo[kDICOMStr], coilElements[kDICOMStr], pulseSequenceDetails[kDICOMStr];
+ epiFactor = siemensCsaAscii(filename, d.CSA.SeriesHeader_offset, d.CSA.SeriesHeader_length, &partialFourier, &echoSpacing, &echoTrainDuration, ¶llelReductionFactorInPlane, coilID, consistencyInfo, coilElements, pulseSequenceDetails, fmriExternalInfo);
+ //printMessage("ES %d ETD %d EPI %d\n", echoSpacing, echoTrainDuration, epiFactor);
+ if (partialFourier > 0) {
+ //https://github.com/ismrmrd/siemens_to_ismrmrd/blob/master/parameter_maps/IsmrmrdParameterMap_Siemens_EPI_FLASHREF.xsl
+ float pf = 1.0f;
+ if (partialFourier == 1) pf = 0.5;
+ if (partialFourier == 2) pf = 0.75;
+ if (partialFourier == 4) pf = 0.875;
+ fprintf(fp, "\t\"PartialFourier\": %g,\n", pf);
+ }
+ if (echoSpacing > 0)
+ fprintf(fp, "\t\"EchoSpacing\": %g,\n", echoSpacing / 1000000.0); //usec -> sec
+ if (echoTrainDuration > 0)
+ fprintf(fp, "\t\"EchoTrainDuration\": %g,\n", echoTrainDuration / 1000000.0); //usec -> sec
+ if (epiFactor > 0)
+ fprintf(fp, "\t\"EPIFactor\": %d,\n", epiFactor);
+ if (strlen(coilID) > 0)
+ fprintf(fp, "\t\"ReceiveCoilName\": \"%s\",\n", coilID);
+ if (strlen(coilElements) > 0)
+ fprintf(fp, "\t\"ReceiveCoilActiveElements\": \"%s\",\n", coilElements);
+ if (strlen(pulseSequenceDetails) > 0)
+ fprintf(fp, "\t\"PulseSequenceDetails\": \"%s\",\n", pulseSequenceDetails);
+ if (strlen(fmriExternalInfo) > 0)
+ fprintf(fp, "\t\"FmriExternalInfo\": \"%s\",\n", fmriExternalInfo);
+ if (strlen(consistencyInfo) > 0)
+ fprintf(fp, "\t\"ConsistencyInfo\": \"%s\",\n", consistencyInfo);
+ if (parallelReductionFactorInPlane > 0) {//AccelFactorPE -> phase encoding
+ if (d.accelFactPE < 1.0) d.accelFactPE = parallelReductionFactorInPlane; //value found in ASCII but not in DICOM (0051,1011)
+ if (parallelReductionFactorInPlane != round(d.accelFactPE))
+ printWarning("ParallelReductionFactorInPlane reported in DICOM [0051,1011] (%g) does not match CSA series value %g\n", round(d.accelFactPE), parallelReductionFactorInPlane);
+ }
+ }
+ #endif
+ if (d.CSA.multiBandFactor > 1) //AccelFactorSlice
+ fprintf(fp, "\t\"MultibandAccelerationFactor\": %d,\n", d.CSA.multiBandFactor);
+ if (strlen(d.imageComments) > 0)
+ fprintf(fp, "\t\"ImageComments\": \"%s\",\n", d.imageComments);
+ if (strlen(opts.imageComments) > 0)
+ fprintf(fp, "\t\"ConversionComments\": \"%s\",\n", opts.imageComments);
+ if (d.echoTrainLength > 1) //>1 as for Siemens EPI this is 1, Siemens uses EPI factor http://mriquestions.com/echo-planar-imaging.html
+ fprintf(fp, "\t\"EchoTrainLength\": %d,\n", d.echoTrainLength);
+ if (d.echoNum > 1)
+ fprintf(fp, "\t\"EchoNumber\": %d,\n", d.echoNum);
+ if (d.isDerived) //DICOM is derived image or non-spatial file (sounds, etc)
+ fprintf(fp, "\t\"RawImage\": false,\n");
+ if (d.acquNum > 0)
+ fprintf(fp, "\t\"AcquisitionNumber\": %d,\n", d.acquNum);
+ if (strlen(d.institutionName) > 0)
+ fprintf(fp, "\t\"InstitutionName\": \"%s\",\n", d.institutionName );
+ if (strlen(d.institutionAddress) > 0)
+ fprintf(fp, "\t\"InstitutionAddress\": \"%s\",\n", d.institutionAddress );
+ if (strlen(d.deviceSerialNumber) > 0)
+ fprintf(fp, "\t\"DeviceSerialNumber\": \"%s\",\n", d.deviceSerialNumber );
+ if (strlen(d.stationName) > 0)
+ fprintf(fp, "\t\"StationName\": \"%s\",\n", d.stationName );
+ if (strlen(d.scanOptions) > 0)
+ fprintf(fp, "\t\"ScanOptions\": \"%s\",\n", d.scanOptions );
+ if (strlen(d.softwareVersions) > 0)
+ fprintf(fp, "\t\"SoftwareVersions\": \"%s\",\n", d.softwareVersions );
+ if (strlen(d.procedureStepDescription) > 0)
+ fprintf(fp, "\t\"ProcedureStepDescription\": \"%s\",\n", d.procedureStepDescription );
+ if (strlen(d.scanningSequence) > 0)
+ fprintf(fp, "\t\"ScanningSequence\": \"%s\",\n", d.scanningSequence );
+ if (strlen(d.sequenceVariant) > 0)
+ fprintf(fp, "\t\"SequenceVariant\": \"%s\",\n", d.sequenceVariant );
+ if (strlen(d.seriesDescription) > 0)
+ fprintf(fp, "\t\"SeriesDescription\": \"%s\",\n", d.seriesDescription );
+ if (strlen(d.bodyPartExamined) > 0)
+ fprintf(fp, "\t\"BodyPartExamined\": \"%s\",\n", d.bodyPartExamined );
+ if (strlen(d.protocolName) > 0)
+ fprintf(fp, "\t\"ProtocolName\": \"%s\",\n", d.protocolName );
+ if (strlen(d.sequenceName) > 0)
+ fprintf(fp, "\t\"SequenceName\": \"%s\",\n", d.sequenceName );
+ if (strlen(d.imageType) > 0) {
+ fprintf(fp, "\t\"ImageType\": [\"");
+ bool isSep = false;
+ for (int i = 0; i < strlen(d.imageType); i++) {
+ if (d.imageType[i] != '_') {
+ if (isSep)
+ fprintf(fp, "\", \"");
+ isSep = false;
+ fprintf(fp, "%c", d.imageType[i]);
+ } else
+ isSep = true;
+ }
+ fprintf(fp, "\"],\n");
+ }
+ //Chris Gorgolewski: BIDS standard specifies ISO8601 date-time format (Example: 2016-07-06T12:49:15.679688)
+ //Lines below directly save DICOM values
+ if (d.acquisitionTime > 0.0 && d.acquisitionDate > 0.0){
+ long acquisitionDate = d.acquisitionDate;
+ double acquisitionTime = d.acquisitionTime;
+ char acqDateTimeBuf[64];
+ //snprintf(acqDateTimeBuf, sizeof acqDateTimeBuf, "%+08ld%+08f", acquisitionDate, acquisitionTime);
+ snprintf(acqDateTimeBuf, sizeof acqDateTimeBuf, "%+08ld%+013.5f", acquisitionDate, acquisitionTime); //CR 20170404 add zero pad so 1:23am appears as +012300.00000 not +12300.00000
+ //printMessage("acquisitionDateTime %s\n",acqDateTimeBuf);
+ int ayear,amonth,aday,ahour,amin;
+ double asec;
+ int count = 0;
+ sscanf(acqDateTimeBuf, "%5d%2d%2d%3d%2d%lf%n", &ayear, &amonth, &aday, &ahour, &amin, &asec, &count); //CR 20170404 %lf not %f for double precision
+ //printf("-%02d-%02dT%02d:%02d:%02.6f\",\n", amonth, aday, ahour, amin, asec);
+ if (count) { // ISO 8601 specifies a sign must exist for distant years.
+ //report time of the day only format, https://www.cs.tut.fi/~jkorpela/iso8601.html
+ fprintf(fp, "\t\"AcquisitionTime\": \"%02d:%02d:%02.6f\",\n",ahour, amin, asec);
+ //report date and time together
+ if (!opts.isAnonymizeBIDS) {
+ fprintf(fp, "\t\"AcquisitionDateTime\": ");
+ fprintf(fp, (ayear >= 0 && ayear <= 9999) ? "\"%4d" : "\"%+4d", ayear);
+ fprintf(fp, "-%02d-%02dT%02d:%02d:%02.6f\",\n", amonth, aday, ahour, amin, asec);
+
+ }
+ } //if (count)
+ } //if acquisitionTime and acquisitionDate recorded
+ // if (d.acquisitionTime > 0.0) fprintf(fp, "\t\"AcquisitionTime\": %f,\n", d.acquisitionTime );
+ // if (d.acquisitionDate > 0.0) fprintf(fp, "\t\"AcquisitionDate\": %8.0f,\n", d.acquisitionDate );
+ //if conditionals: the following values are required for DICOM MRI, but not available for CT
+ if ((d.intenScalePhilips != 0) || (d.manufacturer == kMANUFACTURER_PHILIPS)) { //for details, see PhilipsPrecise()
+ fprintf(fp, "\t\"PhilipsRescaleSlope\": %g,\n", d.intenScale );
+ fprintf(fp, "\t\"PhilipsRescaleIntercept\": %g,\n", d.intenIntercept );
+ fprintf(fp, "\t\"PhilipsScaleSlope\": %g,\n", d.intenScalePhilips );
+ fprintf(fp, "\t\"UsePhilipsFloatNotDisplayScaling\": %d,\n", opts.isPhilipsFloatNotDisplayScaling);
+ }
+ //PET ISOTOPE MODULE ATTRIBUTES
+ if (d.radionuclidePositronFraction > 0.0) fprintf(fp, "\t\"RadionuclidePositronFraction\": %g,\n", d.radionuclidePositronFraction );
+ if (d.radionuclideTotalDose > 0.0) fprintf(fp, "\t\"RadionuclideTotalDose\": %g,\n", d.radionuclideTotalDose );
+ if (d.radionuclideHalfLife > 0.0) fprintf(fp, "\t\"RadionuclideHalfLife\": %g,\n", d.radionuclideHalfLife );
+ if (d.doseCalibrationFactor > 0.0) fprintf(fp, "\t\"DoseCalibrationFactor\": %g,\n", d.doseCalibrationFactor );
+ //MRI parameters
+ if (d.fieldStrength > 0.0) fprintf(fp, "\t\"MagneticFieldStrength\": %g,\n", d.fieldStrength );
+ if (d.flipAngle > 0.0) fprintf(fp, "\t\"FlipAngle\": %g,\n", d.flipAngle );
+ if ((d.TE > 0.0) && (!d.isXRay)) fprintf(fp, "\t\"EchoTime\": %g,\n", d.TE / 1000.0 );
+ if ((d.TE > 0.0) && (d.isXRay)) fprintf(fp, "\t\"XRayExposure\": %g,\n", d.TE );
+ if (d.TR > 0.0) fprintf(fp, "\t\"RepetitionTime\": %g,\n", d.TR / 1000.0 );
+ if (d.TI > 0.0) fprintf(fp, "\t\"InversionTime\": %g,\n", d.TI / 1000.0 );
+ if (d.ecat_isotope_halflife > 0.0) fprintf(fp, "\t\"IsotopeHalfLife\": %g,\n", d.ecat_isotope_halflife);
+ if (d.ecat_dosage > 0.0) fprintf(fp, "\t\"Dosage\": %g,\n", d.ecat_dosage);
+ double bandwidthPerPixelPhaseEncode = d.bandwidthPerPixelPhaseEncode;
+ int phaseEncodingLines = d.phaseEncodingLines;
+ if ((phaseEncodingLines == 0) && (h->dim[2] > 0) && (h->dim[1] > 0)) {
+ if (h->dim[2] == h->dim[2]) //phase encoding does not matter
+ phaseEncodingLines = h->dim[2];
+ else if (d.phaseEncodingRC =='R')
+ phaseEncodingLines = h->dim[2];
+ else if (d.phaseEncodingRC =='C')
+ phaseEncodingLines = h->dim[1];
+ }
+ if (bandwidthPerPixelPhaseEncode == 0.0)
+ bandwidthPerPixelPhaseEncode = d.CSA.bandwidthPerPixelPhaseEncode;
+ if (phaseEncodingLines > 0.0) fprintf(fp, "\t\"PhaseEncodingLines\": %d,\n", phaseEncodingLines );
+ if (bandwidthPerPixelPhaseEncode > 0.0)
+ fprintf(fp, "\t\"BandwidthPerPixelPhaseEncode\": %g,\n", bandwidthPerPixelPhaseEncode );
+ double effectiveEchoSpacing = 0.0;
+ if ((phaseEncodingLines > 0) && (bandwidthPerPixelPhaseEncode > 0.0))
+ effectiveEchoSpacing = 1.0 / (bandwidthPerPixelPhaseEncode * phaseEncodingLines) ;
+ if (d.effectiveEchoSpacingGE > 0.0)
+ effectiveEchoSpacing = d.effectiveEchoSpacingGE / 1000000.0;
+ if (effectiveEchoSpacing > 0.0)
+ fprintf(fp, "\t\"EffectiveEchoSpacing\": %g,\n", effectiveEchoSpacing);
+ //FSL definition is start of first line until start of last line, so n-1 unless accelerated in-plane acquisition
+ // to check: partial Fourier, iPAT, etc.
+ int fencePost = 1;
+ if (d.accelFactPE > 1.0)
+ fencePost = (int)round(d.accelFactPE); //e.g. if 64 lines with iPAT=2, we want time from start of first until start of 62nd effective line
+ if ((d.phaseEncodingSteps > 1) && (effectiveEchoSpacing > 0.0))
+ fprintf(fp, "\t\"TotalReadoutTime\": %g,\n", effectiveEchoSpacing * ((float)d.phaseEncodingSteps - fencePost));
+ if (d.accelFactPE > 1.0) {
+ fprintf(fp, "\t\"ParallelReductionFactorInPlane\": %g,\n", d.accelFactPE);
+ if (effectiveEchoSpacing > 0.0)
+ fprintf(fp, "\t\"TrueEchoSpacing\": %g,\n", effectiveEchoSpacing * d.accelFactPE);
+ }
+ if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (d.dwellTime > 0))
+ fprintf(fp, "\t\"DwellTime\": %g,\n", d.dwellTime * 1E-9);
+ 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 );
+ }
+ fprintf(fp, "\t],\n");
+ }
+ if (((d.phaseEncodingRC == 'R') || (d.phaseEncodingRC == 'C')) && (!d.is3DAcq) && ((d.CSA.phaseEncodingDirectionPositive == 1) || (d.CSA.phaseEncodingDirectionPositive == 0))) {
+ if (d.phaseEncodingRC == 'C') //Values should be "R"ow, "C"olumn or "?"Unknown
+ fprintf(fp, "\t\"PhaseEncodingDirection\": \"j");
+ else if (d.phaseEncodingRC == 'R')
+ fprintf(fp, "\t\"PhaseEncodingDirection\": \"i");
+ else
+ fprintf(fp, "\t\"PhaseEncodingDirection\": \"?");
+ //phaseEncodingDirectionPositive has one of three values: UNKNOWN (-1), NEGATIVE (0), POSITIVE (1)
+ //However, DICOM and NIfTI are reversed in the j (ROW) direction
+ //Equivalent to dicm2nii's "if flp(iPhase), phPos = ~phPos; end"
+ //for samples see https://github.com/rordenlab/dcm2niix/issues/125
+ if (d.CSA.phaseEncodingDirectionPositive == -1)
+ fprintf(fp, "?"); //unknown
+ else if ((d.CSA.phaseEncodingDirectionPositive == 0) && (d.phaseEncodingRC != 'C'))
+ fprintf(fp, "-");
+ else if ((d.phaseEncodingRC == 'C') && (d.CSA.phaseEncodingDirectionPositive == 1) && (opts.isFlipY))
+ fprintf(fp, "-");
+ else if ((d.phaseEncodingRC == 'C') && (d.CSA.phaseEncodingDirectionPositive == 0) && (!opts.isFlipY))
+ fprintf(fp, "-");
+ fprintf(fp, "\",\n");
+ } //only save PhaseEncodingDirection if BOTH direction and POLARITY are known
+ fprintf(fp, "\t\"ConversionSoftware\": \"dcm2niix\",\n");
+ fprintf(fp, "\t\"ConversionSoftwareVersion\": \"%s\"\n", kDCMvers );
+ //fprintf(fp, "\t\"DicomConversion\": [\"dcm2niix\", \"%s\"]\n", kDCMvers );
+ fprintf(fp, "}\n");
+ fclose(fp);
+}// nii_SaveBIDS()
--
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