diff --git a/console/main_console.cpp b/console/main_console.cpp index 152eb41e..adee3e1b 100644 --- a/console/main_console.cpp +++ b/console/main_console.cpp @@ -280,7 +280,11 @@ int main(int argc, const char * argv[]) } else if ((argv[i][1] == 'i') && ((i+1) < argc)) { i++; if (invalidParam(i, argv)) return 0; - if ((argv[i][0] == 'n') || (argv[i][0] == 'N') || (argv[i][0] == '0')) + if ((strlen(argv[i-1]) > 2) && (argv[i-1][2] == 's')) {//"-is y" + printf("Warning: integer scaling is an experimental feature\n"); + if ((argv[i][0] == 'y') || (argv[i][0] == 'Y')) + opts.isMaximize16BitRange = true; + } else if ((argv[i][0] == 'n') || (argv[i][0] == 'N') || (argv[i][0] == '0')) opts.isIgnoreDerivedAnd2D = false; else opts.isIgnoreDerivedAnd2D = true; diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp index 0064df7f..25b2ab73 100644 --- a/console/nii_dicom_batch.cpp +++ b/console/nii_dicom_batch.cpp @@ -2124,6 +2124,66 @@ int nii_saveNII3D(char * niiFilename, struct nifti_1_header hdr, unsigned char* return EXIT_SUCCESS; }// nii_saveNII3D() +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_scale16bitUnsigned(unsigned char *img, struct nifti_1_header *hdr){ + 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]; + //clock_t start = clock(); + 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\n", max16); +} + 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 +2949,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 +3817,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; diff --git a/console/nii_dicom_batch.h b/console/nii_dicom_batch.h index 89209127..663d3363 100644 --- 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];