Merge tag 'v1.0.20180614' into debian

* tag 'v1.0.20180614': If integer scaling is used, append " isN" to NIfTI header (rordenlab#198) Check yaml-cpp version before build. Reseting defaults does not ignore prior arguments in call rordenlab@2cd185f Add '-g i' option (rordenlab#152) Refine lossless 16-bit scaling (rordenlab#198) integer scaling option (rordenlab#198)
neurodebian · Jun 18, 2018 · 749a0d7 · 749a0d7
2 parents 5fc6578 + e77bbbd
commit 749a0d7
Show file tree

Hide file tree

Showing 6 changed files with 141 additions and 9 deletions.
diff --git a/SuperBuild/SuperBuild.cmake b/SuperBuild/SuperBuild.cmake
@@ -28,7 +28,7 @@ if(USE_STATIC_RUNTIME)
         mark_as_advanced(STATIC_LIBCXX)
         if(NOT STATIC_LIBCXX)
             unset(STATIC_LIBCXX CACHE)
-            # only on some Centos/Redhat systems
+            # Only on some Centos/Redhat systems
             message(FATAL_ERROR
                 "\"USE_STATIC_RUNTIME\" set to ON but \"libstdcxx.a\" not found! \
                  \"yum install libstdc++-static\" to resolve the error.")
@@ -91,7 +91,8 @@ if(BATCH_VERSION)
             pkg_check_modules(YAML-CPP yaml-cpp)
         endif()
 
-        if(YAML-CPP_FOUND)
+        # Build from github if not found or version < 0.5.3
+        if(YAML-CPP_FOUND AND NOT (YAML-CPP_VERSION VERSION_LESS "0.5.3"))
             set(YAML-CPP_DIR ${YAML-CPP_LIBDIR}/cmake/yaml-cpp CACHE PATH "Path to yaml-cpp configuration file"  FORCE)
             message("--     Using yaml-cpp library from ${YAML-CPP_DIR}")
         else()

diff --git a/VERSIONS.md b/VERSIONS.md
@@ -8,6 +8,7 @@
  - Support for Philips Private RLE (1.3.46.670589.33.1.4.1) transfer syntax.
  - Optional support for JPEG-LS (1.2.840.10008.1.2.4.80/1.2.840.10008.1.2.4.81) transfer syntaxes (using [CharLS](https://github.com/team-charls/charls)). Requires c++14.
  - [Improved GE support](https://github.com/rordenlab/dcm2niix/issues/163)
+ - Optional [lossless integer scaling](https://github.com/rordenlab/dcm2niix/issues/198) for INT16 and UINT16 DICOM images that only use a fraction of the possible range.
 
 15-Dec-2017
  - Support [Siemens XA10 images](https://github.com/rordenlab/dcm2niix/pull/145).

diff --git a/console/main_console.cpp b/console/main_console.cpp
@@ -86,9 +86,12 @@ void showHelp(const char * argv[], struct TDCMopts opts) {
      #define kQstr ""
     #endif
     printf("  -f : filename (%%a=antenna  (coil) number, %%c=comments, %%d=description, %%e=echo number, %%f=folder name, %%i=ID of patient, %%j=seriesInstanceUID, %%k=studyInstanceUID, %%m=manufacturer, %%n=name of patient, %%p=protocol,%s %%s=series number, %%t=time, %%u=acquisition number, %%v=vendor, %%x=study ID; %%z=sequence name; default '%s')\n", kQstr, opts.filename);
-    printf("  -g : generate defaults file (y/n/o [o=only: reset and write defaults], default n)\n");
+    printf("  -g : generate defaults file (y/n/o/i [o=only: reset and write defaults; i=ignore: reset defaults], default n)\n");
     printf("  -h : show help\n");
     printf("  -i : ignore derived, localizer and 2D images (y/n, default n)\n");
+    char max16Ch = 'n';
+    if (opts.isMaximize16BitRange) max16Ch = 'y';
+    printf("  -l : losslessly scale 16-bit integers to use dynamic range (y/n, default %c)\n", max16Ch);
     printf("  -m : merge 2D slices from same series regardless of study time, echo, coil, orientation, etc. (y/n, default n)\n");
     printf("  -n : only convert this series number - can be used up to %i times (default convert all)\n", MAX_NUM_SERIES);
     printf("  -o : output directory (omit to save to input folder)\n");
@@ -268,6 +271,12 @@ int main(int argc, const char * argv[])
                 if (invalidParam(i, argv)) return 0;
                 if ((argv[i][0] == 'y') || (argv[i][0] == 'Y'))
                     isSaveIni = true;
+                if (((argv[i][0] == 'i') || (argv[i][0] == 'I')) && (!isResetDefaults)) {
+                    isResetDefaults = true;
+                    printf("Defaults reset\n");
+                    setDefaultOpts(&opts, argv);
+                    i = 0; //re-read all settings for this pass, e.g. "dcm2niix -f %p_%s -d o" should save filename as "%p_%s"
+                }
                 if (((argv[i][0] == 'o') || (argv[i][0] == 'O')) && (!isResetDefaults)) {
                 	//reset defaults - do not read, but do write defaults
                     isSaveIni = true;
@@ -284,6 +293,13 @@ int main(int argc, const char * argv[])
                     opts.isIgnoreDerivedAnd2D = false;
                 else
                     opts.isIgnoreDerivedAnd2D = true;
+            } else if ((argv[i][1] == 'l') && ((i+1) < argc)) {
+                i++;
+                if (invalidParam(i, argv)) return 0;
+                if ((argv[i][0] == 'n') || (argv[i][0] == 'N')  || (argv[i][0] == '0'))
+                    opts.isMaximize16BitRange = false;
+                else
+                    opts.isMaximize16BitRange = true;
             } else if ((argv[i][1] == 'm') && ((i+1) < argc)) {
                 i++;
                 if (invalidParam(i, argv)) return 0;

diff --git a/console/nii_dicom.h b/console/nii_dicom.h
@@ -43,7 +43,7 @@ extern "C" {
 	#define kCCsuf " CompilerNA" //unknown compiler!
 #endif
 
-#define kDCMvers "v1.0.20180606" kJP2suf kLSsuf kCCsuf
+#define kDCMvers "v1.0.20180614" kJP2suf kLSsuf kCCsuf
 
 static const int kMaxEPI3D = 1024; //maximum number of EPI images in Siemens Mosaic
 static const int kMaxDTI4D = 18000; //maximum number of DTI directions for 4D (Philips) images, also maximum number of 3D slices for Philips 3D and 4D images

diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp
@@ -2124,6 +2124,109 @@ int nii_saveNII3D(char * niiFilename, struct nifti_1_header hdr, unsigned char*
     return EXIT_SUCCESS;
 }// nii_saveNII3D()
 
+/*
+//this version can convert INT16->UINT16
+// some were concerned about this https://github.com/rordenlab/dcm2niix/issues/198
+void nii_scale16bitSigned(unsigned char *img, struct nifti_1_header *hdr){
+    if (hdr->datatype != DT_INT16) return;
+    int dim3to7 = 1;
+    for (int i = 3; i < 8; i++)
+        if (hdr->dim[i] > 1) dim3to7 = dim3to7 * hdr->dim[i];
+    int nVox = hdr->dim[1]*hdr->dim[2]* dim3to7;
+    if (nVox < 1) return;
+    int16_t * img16 = (int16_t*) img;
+    int16_t max16 = img16[0];
+    int16_t min16 = max16;
+    //clock_t start = clock();
+    for (int i=0; i < nVox; i++) {
+        if (img16[i] < min16)
+            min16 = img16[i];
+        if (img16[i] > max16)
+            max16 = img16[i];
+    }
+    int kMx = 32000; //actually 32767 - maybe a bit of padding for interpolation ringing
+    bool isConvertToUint16 = true; //if false output is always same as input: INT16, if true and no negative values output will be UINT16
+    if ((isConvertToUint16) && (min16 >= 0))
+    	kMx = 64000;
+    int scale = kMx / (int)max16;
+	if (abs(min16) > max16)
+		scale = kMx / (int)abs(min16);
+	if (scale < 2) return; //already uses dynamic range
+    hdr->scl_slope = hdr->scl_slope/ scale;
+	if ((isConvertToUint16) && (min16 >= 0)) { //only positive values: save as UINT16 0..65535
+		hdr->datatype = DT_UINT16;
+		uint16_t * uimg16 = (uint16_t*) img;
+    	for (int i=0; i < nVox; i++)
+    		uimg16[i] = (int)img16[i] * scale;
+	} else {//includes negative values: save as INT16 -32768..32768
+		for (int i=0; i < nVox; i++)
+    		img16[i] = img16[i] * scale;
+	}
+    printMessage("Maximizing 16-bit range: raw %d..%d\n", min16, max16);
+}*/
+
+void nii_storeIntegerScaleFactor(int scale, struct nifti_1_header *hdr) {
+//appends NIfTI header description field with " isN" where N is integer scaling
+	char newstr[256];
+	sprintf(newstr, " is%d", scale);
+	if ((strlen(newstr)+strlen(hdr->descrip)) < 80)
+		strcat (hdr->descrip,newstr);
+}
+void nii_scale16bitSigned(unsigned char *img, struct nifti_1_header *hdr) {
+//lossless scaling of INT16 data: e.g. input with range -100...3200 and scl_slope=1
+//  will be stored as -1000...32000 with scl_slope 0.1
+    if (hdr->datatype != DT_INT16) return;
+    int dim3to7 = 1;
+    for (int i = 3; i < 8; i++)
+        if (hdr->dim[i] > 1) dim3to7 = dim3to7 * hdr->dim[i];
+    int nVox = hdr->dim[1]*hdr->dim[2]* dim3to7;
+    if (nVox < 1) return;
+    int16_t * img16 = (int16_t*) img;
+    int16_t max16 = img16[0];
+    int16_t min16 = max16;
+    for (int i=0; i < nVox; i++) {
+        if (img16[i] < min16)
+            min16 = img16[i];
+        if (img16[i] > max16)
+            max16 = img16[i];
+    }
+    int kMx = 32000; //actually 32767 - maybe a bit of padding for interpolation ringing
+    int scale = kMx / (int)max16;
+	if (abs(min16) > max16)
+		scale = kMx / (int)abs(min16);
+	if (scale < 2) return; //already uses dynamic range
+    hdr->scl_slope = hdr->scl_slope/ scale;
+	for (int i=0; i < nVox; i++)
+    	img16[i] = img16[i] * scale;
+    printMessage("Maximizing 16-bit range: raw %d..%d is%d\n", min16, max16, scale);
+    nii_storeIntegerScaleFactor(scale, hdr);
+}
+
+
+void nii_scale16bitUnsigned(unsigned char *img, struct nifti_1_header *hdr){
+//lossless scaling of UINT16 data: e.g. input with range 0...3200 and scl_slope=1
+//  will be stored as 0...64000 with scl_slope 0.05
+    if (hdr->datatype != DT_UINT16) return;
+    int dim3to7 = 1;
+    for (int i = 3; i < 8; i++)
+        if (hdr->dim[i] > 1) dim3to7 = dim3to7 * hdr->dim[i];
+    int nVox = hdr->dim[1]*hdr->dim[2]* dim3to7;
+    if (nVox < 1) return;
+    uint16_t * img16 = (uint16_t*) img;
+    uint16_t max16 = img16[0];
+    for (int i=0; i < nVox; i++)
+        if (img16[i] > max16)
+            max16 = img16[i];
+    int kMx = 64000; //actually 65535 - maybe a bit of padding for interpolation ringing
+    int scale = kMx / (int)max16;
+	if (scale < 2) return; //already uses dynamic range
+	hdr->scl_slope = hdr->scl_slope/ scale;
+	for (int i=0; i < nVox; i++)
+    	img16[i] = img16[i] * scale;
+    printMessage("Maximizing 16-bit range: raw max %d is%d\n", max16, scale);
+    nii_storeIntegerScaleFactor(scale, hdr);
+}
+
 void nii_check16bitUnsigned(unsigned char *img, struct nifti_1_header *hdr){
     //default NIfTI 16-bit is signed, set to unusual 16-bit unsigned if required...
     if (hdr->datatype != DT_UINT16) return;
@@ -2889,9 +2992,14 @@ int saveDcm2NiiCore(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dc
 	nii_SaveText(pathoutname, dcmList[dcmSort[0].indx], opts, &hdr0, nameList->str[indx]);
 	int numADC = 0;
     int * volOrderIndex = nii_SaveDTI(pathoutname,nConvert, dcmSort, dcmList, opts, sliceDir, dti4D, &numADC);
-    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, opts.isVerbose);
+    if ((opts.isMaximize16BitRange) && (hdr0.datatype == DT_INT16)) {
+    	nii_scale16bitSigned(imgM, &hdr0); //allow INT16 to use full dynamic range
+    } else if ((opts.isMaximize16BitRange) && (hdr0.datatype == DT_UINT16) &&  (!dcmList[dcmSort[0].indx].isSigned)) {
+    	nii_scale16bitUnsigned(imgM, &hdr0); //allow UINT16 to use full dynamic range
+    } else if ((!opts.isMaximize16BitRange) && (hdr0.datatype == DT_UINT16) &&  (!dcmList[dcmSort[0].indx].isSigned))
+    	nii_check16bitUnsigned(imgM, &hdr0); //save UINT16 as INT16 if we can do this losslessly
+    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]);
     //~ if (!dcmList[dcmSort[0].indx].isSlicesSpatiallySequentialPhilips)
     //~ 	nii_reorderSlices(imgM, &hdr0, dti4D);
     if (hdr0.dim[3] < 2)
@@ -3752,6 +3860,7 @@ void setDefaultOpts (struct TDCMopts *opts, const char * argv[]) { //either "set
     #else
     	opts->gzLevel = MZ_DEFAULT_LEVEL; //-1;
     #endif
+    opts->isMaximize16BitRange = false; //e.g. if INT16 image has range 0..500 scale to be 0..50000 with hdr.scl_slope =  hdr.scl_slope * 0.01
     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;
@@ -3784,9 +3893,9 @@ void saveIniFile (struct TDCMopts opts) {
 	}
 	printMessage("Saving defaults to registry\n");
 	DWORD dwValue    = opts.isGz;
-	//RegSetValueEx(hKey,"isGZ", 0, REG_DWORD,reinterpret_cast<BYTE *>(&dwValue),sizeof(dwValue));
-	//RegSetValueExA(hKey, "isGZ", 0, REG_DWORD, (LPDWORD)&dwValue, sizeof(dwValue));
 	RegSetValueExA(hKey, "isGZ", 0, REG_DWORD, reinterpret_cast<BYTE *>(&dwValue), sizeof(dwValue));
+	dwValue    = opts.isMaximize16BitRange;
+	RegSetValueExA(hKey, "isMaximize16BitRange", 0, REG_DWORD, reinterpret_cast<BYTE *>(&dwValue), sizeof(dwValue));
 	RegSetValueExA(hKey,"filename",0, REG_SZ,(LPBYTE)opts.filename, strlen(opts.filename)+1);
 	RegCloseKey(hKey);
 } //saveIniFile()
@@ -3807,6 +3916,8 @@ void readIniFile (struct TDCMopts *opts, const char * argv[]) {
     //if(RegQueryValueExA(hKey,"isGZ", 0, (LPDWORD )&dwDataType, (&dwValue), &vSize) == ERROR_SUCCESS)
     if(RegQueryValueExA(hKey,"isGZ", 0, (LPDWORD )&dwDataType, reinterpret_cast<BYTE *>(&dwValue), &vSize) == ERROR_SUCCESS)
     	opts->isGz = dwValue;
+    if(RegQueryValueExA(hKey,"isMaximize16BitRange", 0, (LPDWORD )&dwDataType, reinterpret_cast<BYTE *>(&dwValue), &vSize) == ERROR_SUCCESS)
+    	opts->isMaximize16BitRange = dwValue;
     vSize = 512;
     char buffer[512];
     if(RegQueryValueExA(hKey,"filename", 0,NULL,(LPBYTE)buffer,&vSize ) == ERROR_SUCCESS )
@@ -3830,6 +3941,8 @@ void readIniFile (struct TDCMopts *opts, const char * argv[]) {
         //printMessage(">%s<->'%s'\n",Setting,Value);
         if ( strcmp(Setting,"isGZ") == 0 )
             opts->isGz = atoi(Value);
+        if ( strcmp(Setting,"isMaximize16BitRange") == 0 )
+            opts->isMaximize16BitRange = atoi(Value);
         else if ( strcmp(Setting,"isBIDS") == 0 )
             opts->isCreateBIDS = atoi(Value);
         else if ( strcmp(Setting,"filename") == 0 )
@@ -3844,6 +3957,7 @@ void saveIniFile (struct TDCMopts opts) {
     if (fp == NULL) return;
     printMessage("Saving defaults file %s\n", opts.optsname);
     fprintf(fp, "isGZ=%d\n", opts.isGz);
+    fprintf(fp, "isMaximize16BitRange=%d\n", opts.isMaximize16BitRange);
     fprintf(fp, "isBIDS=%d\n", opts.isCreateBIDS);
     fprintf(fp, "filename=%s\n", opts.filename);
     fclose(fp);

diff --git a/console/nii_dicom_batch.h b/console/nii_dicom_batch.h
@@ -25,7 +25,7 @@ extern "C" {
 #define MAX_NUM_SERIES 16
 
     struct TDCMopts {
-        bool isSave3D,isGz, isFlipY,  isCreateBIDS, isSortDTIbyBVal, isAnonymizeBIDS, isOnlyBIDS, isCreateText, isIgnoreDerivedAnd2D, isPhilipsFloatNotDisplayScaling, isTiltCorrect, isRGBplanar, isOnlySingleFile, isForceStackSameSeries, isCrop;
+        bool isMaximize16BitRange, isSave3D,isGz, isFlipY,  isCreateBIDS, isSortDTIbyBVal, isAnonymizeBIDS, isOnlyBIDS, isCreateText, isIgnoreDerivedAnd2D, isPhilipsFloatNotDisplayScaling, isTiltCorrect, isRGBplanar, isOnlySingleFile, isForceStackSameSeries, isCrop;
         int isVerbose, compressFlag, dirSearchDepth, gzLevel; //support for compressed data 0=none,
         char filename[512], outdir[512], indir[512], pigzname[512], optsname[512], indirParent[512], imageComments[24];
         float seriesNumber[MAX_NUM_SERIES];