[med-svn] [dcm2niix] 01/05: New upstream version 1.0.20170818
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Wed Aug 23 18:38:39 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 9b44e12b9fd7b59b10484881c38174ec202a3160
Author: Ghislain Antony Vaillant <ghisvail at gmail.com>
Date: Wed Aug 23 19:29:31 2017 +0100
New upstream version 1.0.20170818
---
README.md | 23 ++
console/main_console.cpp | 53 +++-
console/nii_dicom.cpp | 430 +++++++++++++++++++++---------
console/nii_dicom.h | 20 +-
console/nii_dicom_batch.cpp | 222 ++++++++++++----
console/nii_dicom_batch.h | 2 +-
console/nii_ostu_ml.cpp | 628 --------------------------------------------
console/nii_ostu_ml.h | 19 --
console/print.h | 9 +-
9 files changed, 558 insertions(+), 848 deletions(-)
diff --git a/README.md b/README.md
index 43bf91f..fdbcb8d 100644
--- a/README.md
+++ b/README.md
@@ -9,8 +9,19 @@ dcm2niix is a designed to convert neuroimaging data from the DICOM format to the
This software is open source. The bulk of the code is covered by the BSD license. Some units are either public domain (nifti*.*, miniz.c) or use the MIT license (ujpeg.cpp). See the license.txt file for more details.
+## Dependencies
+
+This software should run on macOS, Linux and Windows without requiring any other software. However, if you use dcm2niix to create gz-compressed images it will be faster if you have [pigz](https://github.com/madler/pigz) installed. You can get a version of both dcm2niix and pigz compiled for your operating system by downloading [MRIcroGL](https://www.nitrc.org/projects/mricrogl/).
+
## Versions
+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).
+
+24-July-2017
+ - Compiles with recent releases of [OpenJPEG](https://github.com/neurolabusc/dcm_qa/issues/5#issuecomment-317443179) for JPEG2000 support.
+
23-June-2017
- [Ensure slice timing always encoded for Siemens EPI](https://github.com/neurolabusc/dcm_qa/issues/4#issuecomment-310707906)
- [Integrates validation](https://github.com/neurolabusc/dcm_qa)
@@ -180,3 +191,15 @@ This requires a compiler that supports c++11.
### Building the command line version without cmake
[See the COMPILE.md file for details on manual compilation](./COMPILE.md).
+
+## Links
+
+ - [Dcm2Bids](https://github.com/cbedetti/Dcm2Bids) uses dcm2niix 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.
+ - [dcm2niir](https://github.com/muschellij2/dcm2niir) R wrapper for dcm2niix/dcm2nii.
+ - [divest](https://github.com/jonclayden/divest) R interface to dcm2niix.
+ - [sci-tran dcm2niix](https://github.com/scitran-apps/dcm2niix) docker.
+ - [neuro_docker](https://github.com/Neurita/neuro_docker) includes dcm2niix.
+ - [dcm2niix_afni](https://afni.nimh.nih.gov/pub/dist/doc/program_help/dcm2niix_afni.html) is a version of dcm2niix included with the [AFNI](https://afni.nimh.nih.gov/) distribution.
+ - [MRIcroGL](https://github.com/neurolabusc/MRIcroGL) is available for MacOS, Linux and Windows and provides a graphical interface for dcm2niix. You can get compiled copies from the [MRIcroGL NITRC web site](https://www.nitrc.org/projects/mricrogl/).
\ No newline at end of file
diff --git a/console/main_console.cpp b/console/main_console.cpp
index 282c461..4c4d55b 100644
--- a/console/main_console.cpp
+++ b/console/main_console.cpp
@@ -71,10 +71,10 @@ void showHelp(const char * argv[], struct TDCMopts opts) {
const char *cstr = removePath(argv[0]);
printf("usage: %s [options] <in_folder>\n", cstr);
printf(" Options :\n");
- printf(" -1..-9 : gz compression level (1=fastest, 9=smallest)\n");
+ printf(" -1..-9 : gz compression level (1=fastest..9=smallest, default %d)\n", opts.gzLevel);
char bidsCh = 'n';
if (opts.isCreateBIDS) bidsCh = 'y';
- printf(" -b : BIDS sidecar (y/n, default %c)\n", bidsCh);
+ 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);
#ifdef mySegmentByAcq
@@ -125,8 +125,21 @@ void showHelp(const char * argv[], struct TDCMopts opts) {
#endif
} //showHelp()
-//#define mydebugtest
+int invalidParam(int i, const char * argv[]) {
+ if ((argv[i][0] == 'y') || (argv[i][0] == 'Y')
+ || (argv[i][0] == 'n') || (argv[i][0] == 'N')
+ || (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'))
+ return 0;
+
+ //if (argv[i][0] != '-') return 0;
+ printf(" Error: invalid option '%s %s'\n", argv[i-1], argv[i]);
+ return 1;
+}
+//#define mydebugtest
int main(int argc, const char * argv[])
{
struct TDCMopts opts;
@@ -148,7 +161,7 @@ int main(int argc, const char * argv[])
showHelp(argv, opts);
return 0;
}
- //bool isCustomOutDir = false;
+ //for (int i = 1; i < argc; i++) { printf(" argument %d= '%s'\n", i, argv[i]);}
int i = 1;
int lastCommandArg = 0;
while (i < (argc)) { //-1 as final parameter is DICOM directory
@@ -162,33 +175,40 @@ int main(int argc, const char * argv[])
} else if ((argv[i][1] == 'b') && ((i+1) < argc)) {
if (strlen(argv[i]) < 3) { //"-b y"
i++;
+ if (invalidParam(i, argv)) return 0;
if ((argv[i][0] == 'n') || (argv[i][0] == 'N') || (argv[i][0] == '0'))
opts.isCreateBIDS = false;
- else
+ else {
opts.isCreateBIDS = true;
+ if ((argv[i][0] == 'o') || (argv[i][0] == 'O'))
+ opts.isOnlyBIDS = true;
+ }
} else if (argv[i][2] == 'a') {//"-ba y"
i++;
+ if (invalidParam(i, argv)) return 0;
if ((argv[i][0] == 'n') || (argv[i][0] == 'N') || (argv[i][0] == '0'))
opts.isAnonymizeBIDS = false;
else
opts.isAnonymizeBIDS = true;
-
} else
printf("Error: Unknown command line argument: '%s'\n", argv[i]);
} 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'))
opts.isIgnoreDerivedAnd2D = false;
else
opts.isIgnoreDerivedAnd2D = true;
} else if ((argv[i][1] == 'm') && ((i+1) < argc)) {
i++;
+ if (invalidParam(i, argv)) return 0;
if ((argv[i][0] == 'n') || (argv[i][0] == 'N') || (argv[i][0] == '0'))
opts.isForceStackSameSeries = false;
else
opts.isForceStackSameSeries = true;
} else if ((argv[i][1] == 'p') && ((i+1) < argc)) {
i++;
+ if (invalidParam(i, argv)) return 0;
if ((argv[i][0] == 'n') || (argv[i][0] == 'N') || (argv[i][0] == '0'))
opts.isPhilipsFloatNotDisplayScaling = false;
else
@@ -196,18 +216,21 @@ int main(int argc, const char * argv[])
} else if ((argv[i][1] == 's') && ((i+1) < argc)) {
i++;
+ if (invalidParam(i, argv)) return 0;
if ((argv[i][0] == 'n') || (argv[i][0] == 'N') || (argv[i][0] == '0'))
opts.isOnlySingleFile = false;
else
opts.isOnlySingleFile = true;
} else if ((argv[i][1] == 't') && ((i+1) < argc)) {
i++;
+ if (invalidParam(i, argv)) return 0;
if ((argv[i][0] == 'n') || (argv[i][0] == 'N') || (argv[i][0] == '0'))
opts.isCreateText = false;
else
opts.isCreateText = true;
} else if ((argv[i][1] == 'v') && ((i+1) < argc)) {
i++;
+ if (invalidParam(i, argv)) return 0;
if ((argv[i][0] == 'n') || (argv[i][0] == 'N') || (argv[i][0] == '0')) //0: verbose OFF
opts.isVerbose = 0;
else if ((argv[i][0] == 'h') || (argv[i][0] == 'H') || (argv[i][0] == '2')) //2: verbose HYPER
@@ -216,12 +239,14 @@ int main(int argc, const char * argv[])
opts.isVerbose = 1; //1: verbose ON
} else if ((argv[i][1] == 'x') && ((i+1) < argc)) {
i++;
+ if (invalidParam(i, argv)) return 0;
if ((argv[i][0] == 'n') || (argv[i][0] == 'N') || (argv[i][0] == '0'))
opts.isCrop = false;
else
opts.isCrop = true;
} 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') ) {
opts.isGz = true; //force use of internal compression instead of pigz
strcpy(opts.pigzname,"");
@@ -234,9 +259,9 @@ int main(int argc, const char * argv[])
strcpy(opts.filename,argv[i]);
} else if ((argv[i][1] == 'o') && ((i+1) < argc)) {
i++;
- //isCustomOutDir = true;
strcpy(opts.outdir,argv[i]);
- }
+ } else
+ printf(" Error: invalid option '%s %s'\n", argv[i], argv[i+1]);;
lastCommandArg = i;
} //if parameter is a command
i ++; //read next parameter
@@ -251,15 +276,21 @@ int main(int argc, const char * argv[])
return EXIT_SUCCESS;
}
#endif
+ #ifndef myEnableMultipleInputs
+ if ((argc-lastCommandArg-1) > 1) {
+ printf("Warning: only processing last of %d input files (recompile with 'myEnableMultipleInputs' to recursively process multiple files)\n", argc-lastCommandArg-1);
+ lastCommandArg = argc - 2;
+ }
+ #endif
#if !defined(_WIN64) && !defined(_WIN32)
double startWall = get_wall_time();
#endif
clock_t start = clock();
for (i = (lastCommandArg+1); i < argc; i++) {
strcpy(opts.indir,argv[i]); // [argc-1]
- //if (!isCustomOutDir) strcpy(opts.outdir,opts.indir);
- if (nii_loadDir(&opts) != EXIT_SUCCESS)
- return EXIT_FAILURE;
+ int ret = nii_loadDir(&opts);
+ if (ret != EXIT_SUCCESS)
+ return ret;
}
#if !defined(_WIN64) && !defined(_WIN32)
printf ("Conversion required %f seconds (%f for core code).\n",get_wall_time() - startWall, ((float)(clock()-start))/CLOCKS_PER_SEC);
diff --git a/console/nii_dicom.cpp b/console/nii_dicom.cpp
index 2a017e4..92fdd72 100644
--- a/console/nii_dicom.cpp
+++ b/console/nii_dicom.cpp
@@ -172,7 +172,7 @@ static OPJ_SIZE_T opj_skip_from_buffer(OPJ_SIZE_T p_nb_bytes, BufInfo * p_file)
//fix for https://github.com/neurolabusc/dcm_qa/issues/5
static OPJ_BOOL opj_seek_from_buffer(OPJ_SIZE_T p_nb_bytes, BufInfo * p_file) {
//printf("opj_seek_from_buffer %d + %d -> %d + %d\n", p_file->cur , p_nb_bytes, p_file->buf, p_file->len);
- if ((p_file->buf + p_nb_bytes < p_file->buf + p_file->len ) && (p_nb_bytes >= 0)){
+ if (p_nb_bytes < p_file->len ) {
p_file->cur = p_file->buf + p_nb_bytes;
return OPJ_TRUE;
}
@@ -236,7 +236,14 @@ unsigned char * nii_loadImgCoreOpenJPEG(char* imgname, struct nifti_1_header hdr
if ( !opj_setup_decoder(codec, ¶ms) ) goto cleanup2;
// Read the main header of the codestream and if necessary the JP2 boxes
if(! opj_read_header( stream, codec, &jpx)){
- printError( "OpenJPEG failed to read the header %s\n", imgname);
+ printError( "OpenJPEG failed to read the header %s (offset %d)\n", imgname, dcm.imageStart);
+ //comment these next lines to abort: include these to create zero-padded slice
+ #ifdef MY_ZEROFILLBROKENJPGS
+ //fix broken slices https://github.com/scitran-apps/dcm2niix/issues/4
+ printError( "Zero-filled slice created\n");
+ int imgbytes = (hdr.bitpix/8)*hdr.dim[1]*hdr.dim[2];
+ ret = (unsigned char*) calloc(imgbytes,1);
+ #endif
goto cleanup2;
}
// Get the decoded image
@@ -259,7 +266,6 @@ cleanup2:
#define M_PI 3.14159265358979323846
#endif
-#ifdef MY_DEBUG
float deFuzz(float v) {
if (fabs(v) < 0.00001)
return 0;
@@ -268,6 +274,7 @@ float deFuzz(float v) {
}
+#ifdef MY_DEBUG
void reportMat33(char *str, mat33 A) {
printMessage("%s = [%g %g %g ; %g %g %g; %g %g %g ]\n",str,
deFuzz(A.m[0][0]),deFuzz(A.m[0][1]),deFuzz(A.m[0][2]),
@@ -284,7 +291,7 @@ void reportMat44(char *str, mat44 A) {
}
#endif
-int verify_slice_dir (struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1_header *h, mat44 *R){
+int verify_slice_dir (struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1_header *h, mat44 *R, int isVerbose){
//returns slice direction: 1=sag,2=coronal,3=axial, -= flipped
if (h->dim[3] < 2) return 0; //don't care direction for single slice
int iSL = 1; //find Z-slice direction: row with highest magnitude of 3rd column
@@ -311,6 +318,7 @@ int verify_slice_dir (struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1_
pos = d.stackOffcentre[iSL];
if (isnan(pos) && ( !isnan(d.lastScanLoc)) )
pos = d.lastScanLoc;
+ //if (isnan(pos))
vec4 x;
x.v[0] = 0.0; x.v[1] = 0.0; x.v[2]=(float)(h->dim[3]-1.0); x.v[3] = 1.0;
vec4 pos1v = nifti_vect44mat44_mul(x, *R);
@@ -319,13 +327,34 @@ int verify_slice_dir (struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1_
if (!isnan(pos)) // we have real SliceLocation for last slice or volume center
flip = (pos > R->m[iSL-1][3]) != (pos1 > R->m[iSL-1][3]); // same direction?, note C indices from 0
else {// we do some guess work and warn user
- if (!d.isNonImage) //do not warn user if image is derived
- printWarning("Unable to determine slice direction: please check whether slices are flipped\n");
+ vec3 readV = setVec3(d.orient[1],d.orient[2],d.orient[3]);
+ vec3 phaseV = setVec3(d.orient[4],d.orient[5],d.orient[6]);
+ //printMessage("rd %g %g %g\n",readV.v[0],readV.v[1],readV.v[2]);
+ //printMessage("ph %g %g %g\n",phaseV.v[0],phaseV.v[1],phaseV.v[2]);
+ vec3 sliceV = crossProduct(readV, phaseV); //order important: this is our hail mary
+ 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
+ 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");
+ }
}
if (flip) {
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
@@ -506,18 +535,28 @@ mat44 set_nii_header(struct TDICOMdata d) {
}
#endif
-float deFuzz(float v) {
- if (fabs(v) < 0.00001)
- return 0;
- else
- return v;
-
-}
-
// This code predates Xiangrui Li's set_nii_header function
-mat44 set_nii_header_x(struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1_header *h, int* sliceDir) {
+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) {
+ //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
+ // http://www.cs.ucl.ac.uk/fileadmin/cmic/Documents/DavidAtkinson/DICOM.pdf
+ // https://www.slicer.org/wiki/Coordinate_systems
+ LOAD_MAT44(Q44, -h->pixdim[1],0,0,0, 0,-h->pixdim[2],0,0, 0,0,h->pixdim[3],0); //X and Y dimensions flipped in NIfTI (RAS) vs DICOM (LPS)
+ vec4 originVx = setVec4( (h->dim[1]+1.0f)/2.0f, (h->dim[2]+1.0f)/2.0f, (h->dim[3]+1.0f)/2.0f);
+ vec4 originMm = nifti_vect44mat44_mul(originVx, Q44);
+ for (int i = 0; i < 3; i++)
+ Q44.m[i][3] = -originMm.v[i]; //set origin to center voxel
+ if (isVerbose) {
+ //printMessage("origin (vx) %g %g %g\n",originVx.v[0],originVx.v[1],originVx.v[2]);
+ //printMessage("origin (mm) %g %g %g\n",originMm.v[0],originMm.v[1],originMm.v[2]);
+ printWarning("Segami coordinates defy DICOM convention, please check orientation\n");
+ }
+ return Q44;
+ }
if (d.CSA.mosaicSlices > 1) {
double nRowCol = ceil(sqrt((double) d.CSA.mosaicSlices));
double lFactorX = (d.xyzDim[1] -(d.xyzDim[1]/nRowCol) )/2.0;
@@ -540,7 +579,7 @@ mat44 set_nii_header_x(struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1
Q44 = nifti_mat44_mul(Q44,det);
}
} else { //not a mosaic
- *sliceDir = verify_slice_dir(d, d2, h, &Q44);
+ *sliceDir = verify_slice_dir(d, d2, h, &Q44, isVerbose);
for (int c=0; c<4; c++)// LPS to nifti RAS, xform matrix before reorient
for (int r=0; r<2; r++) //swap rows 1 & 2
Q44.m[r][c] = - Q44.m[r][c];
@@ -556,7 +595,7 @@ int headerDcm2NiiSForm(struct TDICOMdata d, struct TDICOMdata d2, struct nifti_
//returns sliceDir: 0=unknown,1=sag,2=coro,3=axial,-=reversed slices
int sliceDir = 0;
if (h->dim[3] < 2) {
- mat44 Q44 = set_nii_header_x(d, d2, h, &sliceDir);
+ mat44 Q44 = set_nii_header_x(d, d2, h, &sliceDir, isVerbose);
setQSForm(h,Q44, isVerbose);
return sliceDir; //don't care direction for single slice
}
@@ -575,12 +614,12 @@ int headerDcm2NiiSForm(struct TDICOMdata d, struct TDICOMdata d2, struct nifti_
printMessage("Unable to determine spatial orientation: 0020,0037 missing!\n");
}
}
- mat44 Q44 = set_nii_header_x(d, d2, h, &sliceDir);
+ mat44 Q44 = set_nii_header_x(d, d2, h, &sliceDir, isVerbose);
setQSForm(h,Q44, isVerbose);
return sliceDir;
} //headerDcm2NiiSForm()
-int headerDcm2Nii2(struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1_header *h) { //final pass after de-mosaic
+int headerDcm2Nii2(struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1_header *h, int isVerbose) { //final pass after de-mosaic
char txt[1024] = {""};
if (h->slice_code == NIFTI_SLICE_UNKNOWN) h->slice_code = d.CSA.sliceOrder;
if (h->slice_code == NIFTI_SLICE_UNKNOWN) h->slice_code = d2.CSA.sliceOrder; //sometimes the first slice order is screwed up https://github.com/eauerbach/CMRR-MB/issues/29
@@ -608,7 +647,7 @@ int headerDcm2Nii2(struct TDICOMdata d, struct TDICOMdata d2, struct nifti_1_hea
snprintf(h->descrip,80, "%s",txt);
if (strlen(d.imageComments) > 0)
snprintf(h->aux_file,24,"%s",d.imageComments);
- return headerDcm2NiiSForm(d,d2, h, true);
+ return headerDcm2NiiSForm(d,d2, h, isVerbose);
} //headerDcm2Nii2()
int dcmStrLen (int len) {
@@ -622,6 +661,8 @@ struct TDICOMdata clear_dicom_data() {
struct TDICOMdata d;
//d.dti4D = NULL;
d.locationsInAcquisition = 0;
+ d.modality = kMODALITY_UNKNOWN;
+ d.effectiveEchoSpacingGE = 0;
for (int i=0; i < 4; i++) {
d.CSA.dtiV[i] = 0;
d.patientPosition[i] = NAN;
@@ -897,6 +938,8 @@ 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
@@ -1253,7 +1296,7 @@ int headerDcm2Nii(struct TDICOMdata d, struct nifti_1_header *h, bool isComputeS
h->sizeof_hdr = 348; //used to signify header does not need to be byte-swapped
h->slice_code = d.CSA.sliceOrder;
if (isComputeSForm)
- headerDcm2Nii2(d, d, h);
+ headerDcm2Nii2(d, d, h, false);
return EXIT_SUCCESS;
} // headerDcm2Nii()
@@ -1356,7 +1399,7 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D
} else if (parVers < 41)
nCols = 32; //e.g PAR 4.0
else if (parVers < 42)
- nCols = 47; //e.g. PAR 4.1
+ nCols = kv1; //e.g. PAR 4.1 - last column is final diffusion b-value
else
nCols = kMaxCols; //e.g. PAR 4.2
}
@@ -1423,7 +1466,7 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D
free (cols);
return d;
}
- for (int i = 0; i < nCols; i++)
+ for (int i = 0; i <= nCols; i++)
cols[i] = strtof(p, &p); // p+1 skip comma, read a float
if ((cols[kIndex]) != slice) isIndexSequential = false; //slices 0,1,2.. should have indices 0,1,2,3...
slice ++;
@@ -1503,7 +1546,6 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D
if (cols[kbval] > maxBValue)
maxBValue = cols[kbval];
} //save DTI direction
-
}
} //if DTI directions
//printMessage("%f %f %lu\n",cols[9],cols[kGradientNumber], strlen(buff))
@@ -1612,6 +1654,19 @@ struct TDICOMdata nii_readParRec (char * parname, int isVerbose, struct TDTI4D
printWarning("Reported TR=%gms, measured TR=%gms (prospect. motion corr.?)\n", d.TR, TRms);
d.TR = TRms;
}
+ //check DTI makes sense
+ if (d.CSA.numDti > 1) {
+ bool v1varies = false;
+ bool v2varies = false;
+ bool v3varies = false;
+ for (int i = 1; i < d.CSA.numDti; i++) {
+ if (dti4D->S[0].V[1] != dti4D->S[i].V[1]) v1varies = true;
+ if (dti4D->S[0].V[2] != dti4D->S[i].V[2]) v2varies = true;
+ if (dti4D->S[0].V[3] != dti4D->S[i].V[3]) v3varies = true;
+ }
+ if ((!v1varies) || (!v2varies) || (!v3varies))
+ printError("Bizarre b-vectors %s\n", parname);
+ }
return d;
} //nii_readParRec()
@@ -1902,42 +1957,34 @@ unsigned char * nii_loadImgCore(char* imgname, struct nifti_1_header hdr, int bi
} //nii_loadImg()
unsigned char * nii_planar2rgb(unsigned char* bImg, struct nifti_1_header *hdr, int isPlanar) {
- //DICOM data saved in triples RGBRGBRGB, NIfTI RGB saved in planes RRR..RGGG..GBBBB..B
- if (bImg == NULL) return NULL;
- if (hdr->datatype != DT_RGB24) return bImg;
- if (isPlanar == 0) return bImg;//return nii_bgr2rgb(bImg,hdr);
- int dim3to7 = 1;
- for (int i = 3; i < 8; i++)
- if (hdr->dim[i] > 1) dim3to7 = dim3to7 * hdr->dim[i];
- //int sliceBytes24 = hdr->dim[1]*hdr->dim[2] * hdr->bitpix/8;
- int sliceBytes8 = hdr->dim[1]*hdr->dim[2];
- int sliceBytes24 = sliceBytes8 * 3;
-//#ifdef _MSC_VER
- unsigned char * slice24 = (unsigned char *)malloc(sizeof(unsigned char) * (sliceBytes24));
-//#else
-// unsigned char slice24[ sliceBytes24 ];
-//#endif
- int sliceOffsetRGB = 0;
- int sliceOffsetR = 0;
- int sliceOffsetG = sliceOffsetR + sliceBytes8;
- int sliceOffsetB = sliceOffsetR + 2*sliceBytes8;
- for (int sl = 0; sl < dim3to7; sl++) { //for each 2D slice
- memcpy(slice24, &bImg[sliceOffsetRGB], sliceBytes24);
- //int sliceOffsetG = sliceOffsetR + sliceBytes8;
- //int sliceOffsetB = sliceOffsetR + 2*sliceBytes8;
- int i = 0;
- for (int rgb = 0; rgb < sliceBytes8; rgb++) {
- bImg[i++] =slice24[sliceOffsetR+rgb];
- bImg[i++] =slice24[sliceOffsetG+rgb];
- bImg[i++] =slice24[sliceOffsetB+rgb];
- }
- sliceOffsetRGB += sliceBytes24;
- } //for each slice
-//#ifdef _MSC_VER
- free(slice24);
-//#endif
- return bImg;
-} //nii_rgb2Planar()
+ //DICOM data saved in triples RGBRGBRGB, NIfTI RGB saved in planes RRR..RGGG..GBBBB..B
+ if (bImg == NULL) return NULL;
+ if (hdr->datatype != DT_RGB24) return bImg;
+ if (isPlanar == 0) return bImg;
+ int dim3to7 = 1;
+ for (int i = 3; i < 8; i++)
+ if (hdr->dim[i] > 1) dim3to7 = dim3to7 * hdr->dim[i];
+ int sliceBytes8 = hdr->dim[1]*hdr->dim[2];
+ int sliceBytes24 = sliceBytes8 * 3;
+ unsigned char * slice24 = (unsigned char *)malloc(sizeof(unsigned char) * (sliceBytes24));
+ int sliceOffsetRGB = 0;
+ int sliceOffsetR = 0;
+ int sliceOffsetG = sliceOffsetR + sliceBytes8;
+ int sliceOffsetB = sliceOffsetR + 2*sliceBytes8;
+ //printMessage("planar->rgb %dx%dx%d\n", hdr->dim[1],hdr->dim[2], dim3to7);
+ int i = 0;
+ for (int sl = 0; sl < dim3to7; sl++) { //for each 2D slice
+ memcpy(slice24, &bImg[sliceOffsetRGB], sliceBytes24);
+ for (int rgb = 0; rgb < sliceBytes8; rgb++) {
+ bImg[i++] =slice24[sliceOffsetR+rgb];
+ bImg[i++] =slice24[sliceOffsetG+rgb];
+ bImg[i++] =slice24[sliceOffsetB+rgb];
+ }
+ sliceOffsetRGB += sliceBytes24;
+ } //for each slice
+ free(slice24);
+ return bImg;
+} //nii_planar2rgb()
unsigned char * nii_rgb2planar(unsigned char* bImg, struct nifti_1_header *hdr, int isPlanar) {
//DICOM data saved in triples RGBRGBRGB, Analyze RGB saved in planes RRR..RGGG..GBBBB..B
@@ -1947,14 +1994,9 @@ unsigned char * nii_rgb2planar(unsigned char* bImg, struct nifti_1_header *hdr,
int dim3to7 = 1;
for (int i = 3; i < 8; i++)
if (hdr->dim[i] > 1) dim3to7 = dim3to7 * hdr->dim[i];
- //int sliceBytes24 = hdr->dim[1]*hdr->dim[2] * hdr->bitpix/8;
int sliceBytes8 = hdr->dim[1]*hdr->dim[2];
int sliceBytes24 = sliceBytes8 * 3;
-//#ifdef _MSC_VER
unsigned char * slice24 = (unsigned char *)malloc(sizeof(unsigned char) * (sliceBytes24));
-//#else
-// unsigned char slice24[ sliceBytes24 ];
-//#endif
//printMessage("rgb->planar %dx%dx%d\n", hdr->dim[1],hdr->dim[2], dim3to7);
int sliceOffsetR = 0;
for (int sl = 0; sl < dim3to7; sl++) { //for each 2D slice
@@ -1971,9 +2013,7 @@ unsigned char * nii_rgb2planar(unsigned char* bImg, struct nifti_1_header *hdr,
}
sliceOffsetR += sliceBytes24;
} //for each slice
-//#ifdef _MSC_VER
free(slice24);
-//#endif
return bImg;
} //nii_rgb2Planar()
@@ -2242,7 +2282,7 @@ TJPEG * decode_JPEG_SOF_0XC3_stack (const char *fn, int skipBytes, bool isVerbo
}
free(lRawRA);
if (frame < frames) {
- printMessage("Only found %d of %d JPEG fragments. Please use dcmdjpeg to uncompress data.\n", frame, frames);
+ printMessage("Only found %d of %d JPEG fragments. Please use dcmdjpeg or gdcmconv to uncompress data.\n", frame, frames);
abortGoto();
}
return lOffsetRA;
@@ -2392,6 +2432,106 @@ unsigned char * nii_loadImgJPEG50(char* imgname, struct nifti_1_header hdr, stru
#endif
#endif
+uint32_t rleInt(int lIndex, unsigned char lBuffer[], bool swap) {//read binary 32-bit integer
+ uint32_t retVal = 0;
+ memcpy(&retVal, (char*)&lBuffer[lIndex * 4], 4);
+ if (!swap) return retVal;
+ uint32_t swapVal;
+ char *inInt = ( char* ) & retVal;
+ char *outInt = ( char* ) & swapVal;
+ outInt[0] = inInt[3];
+ outInt[1] = inInt[2];
+ outInt[2] = inInt[1];
+ outInt[3] = inInt[0];
+ return swapVal;
+} //rleInt()
+
+unsigned char * nii_loadImgRLE(char* imgname, struct nifti_1_header hdr, struct TDICOMdata dcm) {
+//decompress PackBits run-length encoding https://en.wikipedia.org/wiki/PackBits
+ if (dcm.imageBytes < 66 ) { //64 for header+ 2 byte minimum image
+ printError("%d is not enough bytes for RLE compression '%s'\n", dcm.imageBytes, imgname);
+ return NULL;
+ }
+ FILE *file = fopen(imgname , "rb");
+ if (!file) {
+ printError("Unable to open %s\n", imgname);
+ return NULL;
+ }
+ fseek(file, 0, SEEK_END);
+ long fileLen=ftell(file);
+ if (fileLen < (dcm.imageBytes+dcm.imageStart)) {
+ printMessage("File not large enough to store image data: %s\n", imgname);
+ fclose(file);
+ return NULL;
+ }
+ fseek(file, (long) dcm.imageStart, SEEK_SET);
+ size_t imgsz = nii_ImgBytes(hdr);
+ 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) {
+ printError("Only loaded %zu of %d bytes for %s\n", sz, dcm.imageBytes, imgname);
+ free(cImg);
+ return NULL;
+ }
+ //read header http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_G.3.html
+ bool swap = (dcm.isLittleEndian != littleEndianPlatform());
+ int bytesPerSample = dcm.samplesPerPixel * (dcm.bitsAllocated / 8);
+ uint32_t bytesPerSampleRLE = rleInt(0, cImg, swap);
+ if (bytesPerSampleRLE != (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++)
+ bImg[i] = 0;
+ for (int i = 0; i < bytesPerSample; i++) {
+ uint32_t offset = rleInt(i+1, cImg, swap);
+ if (offset > dcm.imageBytes) {
+ printError("RLE header error\n");
+ free(cImg);
+ free(bImg);
+ return NULL;
+ }
+ //save in platform's endian:
+ // The first Segment is generated by stripping off the most significant byte of each Padded Composite Pixel Code...
+ size_t vx = i;
+ if ((dcm.samplesPerPixel == 1) && (littleEndianPlatform())) //endian, except for RGB
+ vx = (bytesPerSample-1) - i;
+ while (vx < imgsz) {
+ int8_t n = cImg[offset];
+ offset++;
+ //http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_G.3.html
+ if ((n >= 0) && (n <= 127)) { //literal bytes
+ int reps = 1 + (int)n;
+ for (int r = 0; r < reps; r++) {
+ int8_t v = cImg[offset];
+ offset++;
+ if (vx >= imgsz)
+ ;//printMessage("literal overflow %d %d\n", r, reps);
+ else
+ bImg[vx] = v;
+ vx = vx + bytesPerSample;
+ }
+ } else if ((n <= -1) && (n >= -127)) { //repeated run
+ int8_t v = cImg[offset];
+ offset++;
+ int reps = -(int)n + 1;
+ for (int r = 0; r < reps; r++) {
+ if (vx >= imgsz)
+ ;//printMessage("repeat overflow %d\n", reps);
+ else
+ bImg[vx] = v;
+ vx = vx + bytesPerSample;
+ }
+ }; //n.b. we ignore -128!
+ } //while vx < imgsz
+ } //for i < bytesPerSample
+ free(cImg);
+ return bImg;
+} // nii_loadImgRLE()
+
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
@@ -2405,9 +2545,14 @@ unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct
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
- } else if (dcm.compressionScheme == kCompressC3) {
- img = nii_loadImgJPEGC3(imgname, *hdr, dcm, (isVerbose > 0));
- } else
+ } else if (dcm.compressionScheme == kCompressRLE) {
+ img = nii_loadImgRLE(imgname, *hdr, 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
+
+ } else if (dcm.compressionScheme == kCompressC3)
+ img = nii_loadImgJPEGC3(imgname, *hdr, dcm, (isVerbose > 0));
+ else
#ifndef myDisableOpenJPEG
if ( ((dcm.compressionScheme == kCompress50) || (dcm.compressionScheme == kCompressYes)) && (compressFlag != kCompressNone) )
img = nii_loadImgCoreOpenJPEG(imgname, *hdr, dcm, compressFlag);
@@ -2425,11 +2570,9 @@ unsigned char * nii_loadImgXL(char* imgname, struct nifti_1_header *hdr, struct
} else
img = nii_loadImgCore(imgname, *hdr, dcm.bitsAllocated);
if (img == NULL) return img;
- if (dcm.compressionScheme == kCompressNone) {
- if ((dcm.isLittleEndian != littleEndianPlatform()) && (hdr->bitpix > 8))
+ if ((dcm.compressionScheme == kCompressNone) && (dcm.isLittleEndian != littleEndianPlatform()) && (hdr->bitpix > 8))
img = nii_byteswap(img, hdr);
- }
- if ((dcm.compressionScheme == kCompressNone) && (hdr->datatype ==DT_RGB24)) //img = nii_planar2rgb(img, hdr, dcm.isPlanarRGB); //
+ if ((dcm.compressionScheme == kCompressNone) && (hdr->datatype ==DT_RGB24))
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) {
@@ -2522,6 +2665,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
#endif
if (MaxBufferSz > fileLen)
MaxBufferSz = fileLen;
+ //printf("%d -> %d\n", MaxBufferSz, fileLen);
long lFileOffset = 0;
fseek(file, 0, SEEK_SET);
//Allocate memory
@@ -2552,6 +2696,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
#define kAcquisitionDateTime 0x0008+(0x002A << 16 )
#define kStudyTime 0x0008+(0x0030 << 16 )
#define kAcquisitionTime 0x0008+(0x0032 << 16 )
+#define kModality 0x0008+(0x0060 << 16 ) //CS
#define kManufacturer 0x0008+(0x0070 << 16 )
#define kInstitutionName 0x0008+(0x0080 << 16 )
#define kInstitutionAddress 0x0008+(0x0081 << 16 )
@@ -2607,6 +2752,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
#define kImageComments 0x0020+(0x4000<< 16 )// '0020' '4000' 'LT' 'ImageComments'
#define kLocationsInAcquisitionGE 0x0021+(0x104F<< 16 )// 'SS' 'LocationsInAcquisitionGE'
#define kSamplesPerPixel 0x0028+(0x0002 << 16 )
+#define kPhotometricInterpretation 0x0028+(0x0004 << 16 )
#define kPlanarRGB 0x0028+(0x0006 << 16 )
#define kDim3 0x0028+(0x0008 << 16 ) //number of frames - for Philips this is Dim3*Dim4
#define kDim2 0x0028+(0x0010 << 16 )
@@ -2624,6 +2770,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
#define kProcedureStepDescription 0x0040+(0x0254 << 16 )
#define kRealWorldIntercept 0x0040+uint32_t(0x9224 << 16 ) //IS dicm2nii's SlopInt_6_9
#define kRealWorldSlope 0x0040+uint32_t(0x9225 << 16 ) //IS dicm2nii's SlopInt_6_9
+#define kEffectiveEchoSpacingGE 0x0043+(0x102C << 16 ) //SS
#define kDiffusionBFactorGE 0x0043+(0x1039 << 16 ) //IS dicm2nii's SlopInt_6_9
#define kCoilSiemens 0x0051+(0x100F << 16 )
#define kImaPATModeText 0x0051+(0x1011 << 16 )
@@ -2797,7 +2944,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
//verbose reporting :
// printMessage("Pos %ld GroupElement %#08x,%#08x Length %d isLittle %d\n", lPos, (groupElement & 0xFFFF), (groupElement >> 16), lLength, d.isLittleEndian);
switch ( groupElement ) {
- case kTransferSyntax: {
+ case kTransferSyntax: {
char transferSyntax[kDICOMStr];
dcmStr (lLength, &buffer[lPos], transferSyntax);
//printMessage("%d transfer syntax>>> '%s'\n", compressFlag, transferSyntax);
@@ -2820,20 +2967,20 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
d.compressionScheme = kCompressC3;
//printMessage("Ancient JPEG-lossless (SOF type 0xc3): please check conversion\n");
} else if (strcmp(transferSyntax, "1.2.840.10008.1.2.4.70") == 0) {
- //d.isCompressed = true;
d.compressionScheme = kCompressC3;
- //printMessage("Ancient JPEG-lossless (SOF type 0xc3): please check conversion\n");
- //d.imageStart = 1;//abort as invalid (imageStart MUST be >128)
+ } else if (strcmp(transferSyntax, "1.2.840.10008.1.2.4.80") == 0) {
+ printMessage("Unsupported transfer syntax '%s' (decode with dcmdjpls or gdcmconv)\n",transferSyntax);
+ d.imageStart = 1;//abort as invalid (imageStart MUST be >128)
} else if ((compressFlag != kCompressNone) && (strcmp(transferSyntax, "1.2.840.10008.1.2.4.90") == 0)) {
d.compressionScheme = kCompressYes;
//printMessage("JPEG2000 Lossless support is new: please validate conversion\n");
} else if ((compressFlag != kCompressNone) && (strcmp(transferSyntax, "1.2.840.10008.1.2.4.91") == 0)) {
d.compressionScheme = kCompressYes;
//printMessage("JPEG2000 support is new: please validate conversion\n");
- } else if (strcmp(transferSyntax, "1.2.840.10008.1.2.2") == 0)
+ } else if (strcmp(transferSyntax, "1.2.840.10008.1.2.5") == 0)
+ d.compressionScheme = kCompressRLE; //run length
+ else if (strcmp(transferSyntax, "1.2.840.10008.1.2.2") == 0)
isSwitchToBigEndian = true; //isExplicitVR=true;
- //else if (strcmp(transferSyntax, "1.2.840.10008.1.2.4.91") == 0)
- // ;
else if (strcmp(transferSyntax, "1.2.840.10008.1.2") == 0)
isSwitchToImplicitVR = true; //d.isLittleEndian=true
else {
@@ -2868,7 +3015,20 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
case kStudyDate:
dcmStr (lLength, &buffer[lPos], d.studyDate);
break;
- case kManufacturer:
+ case kModality:
+ if (lLength < 2) break;
+ if ((buffer[lPos]=='C') && (toupper(buffer[lPos+1]) == 'R'))
+ d.modality = kMODALITY_CR;
+ else if ((buffer[lPos]=='C') && (toupper(buffer[lPos+1]) == 'T'))
+ d.modality = kMODALITY_CT;
+ if ((buffer[lPos]=='M') && (toupper(buffer[lPos+1]) == 'R'))
+ d.modality = kMODALITY_MR;
+ if ((buffer[lPos]=='P') && (toupper(buffer[lPos+1]) == 'T'))
+ d.modality = kMODALITY_PT;
+ if ((buffer[lPos]=='U') && (toupper(buffer[lPos+1]) == 'S'))
+ d.modality = kMODALITY_US;
+ break;
+ case kManufacturer:
d.manufacturer = dcmStrManufacturer (lLength, &buffer[lPos]);
break;
case kInstitutionName:
@@ -2880,20 +3040,21 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
case kReferringPhysicianName:
dcmStr(lLength, &buffer[lPos], d.referringPhysicianName);
break;
- case kComplexImageComponent:
+ case kComplexImageComponent:
+ if (lLength < 2) break;
d.isHasPhase = (buffer[lPos]=='P') && (toupper(buffer[lPos+1]) == 'H');
d.isHasMagnitude = (buffer[lPos]=='M') && (toupper(buffer[lPos+1]) == 'A');
break;
- case kAcquisitionTime :
+ case kAcquisitionTime :
char acquisitionTimeTxt[kDICOMStr];
dcmStr (lLength, &buffer[lPos], acquisitionTimeTxt);
d.acquisitionTime = atof(acquisitionTimeTxt);
//printMessage("%s\n",acquisitionTimeTxt);
break;
- case kStudyTime :
+ case kStudyTime :
dcmStr (lLength, &buffer[lPos], d.studyTime);
break;
- case kPatientName :
+ case kPatientName :
dcmStr (lLength, &buffer[lPos], d.patientName);
break;
case kPatientID :
@@ -2924,7 +3085,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
//if ((strlen(d.protocolName) < 1) || (d.manufacturer != kMANUFACTURER_GE)) //GE uses a generic session name here: do not overwrite kProtocolNameGE
dcmStr (lLength, &buffer[lPos], d.protocolName); //see also kSequenceName
break; }
- case kPatientOrient :
+ case kPatientOrient :
dcmStr (lLength, &buffer[lPos], d.patientOrient);
break;
case kLastScanLoc :
@@ -2953,11 +3114,10 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
case kStudyInstanceUID :
dcmStr (lLength, &buffer[lPos], d.studyInstanceUID);
break;
-
case kSeriesInstanceUID :
dcmStr (lLength, &buffer[lPos], d.seriesInstanceUID);
break;
- case kPatientPosition :
+ case kPatientPosition :
if ((d.manufacturer == kMANUFACTURER_PHILIPS) && (is2005140FSQ)) {
if (!is2005140FSQwarned)
printWarning("Philips R3.2.2 can report different positions for the same slice. Attempting patch.\n");
@@ -2987,65 +3147,72 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
} //not after 2005,140F
break;
- case kInPlanePhaseEncodingDirection:
+ case kInPlanePhaseEncodingDirection:
d.phaseEncodingRC = toupper(buffer[lPos]); //first character is either 'R'ow or 'C'ol
break;
case kStudyID:
dcmStr (lLength, &buffer[lPos], d.studyID);
break;
- case kSeriesNum:
+ case kSeriesNum:
d.seriesNum = dcmStrInt(lLength, &buffer[lPos]);
break;
- case kAcquNum:
+ case kAcquNum:
d.acquNum = dcmStrInt(lLength, &buffer[lPos]);
break;
- case kImageNum:
+ case kImageNum:
//int dx = 3;
if (d.imageNum <= 1) d.imageNum = dcmStrInt(lLength, &buffer[lPos]); //Philips renames each image as image1 in 2001,9000
break;
- case kPlanarRGB:
+ case kPhotometricInterpretation:
+ char interp[kDICOMStr];
+ dcmStr (lLength, &buffer[lPos], interp);
+ if (strcmp(interp, "PALETTE_COLOR") == 0)
+ printError("Photometric Interpretation 'PALETTE COLOR' not supported\n");
+
+ break;
+ case kPlanarRGB:
d.isPlanarRGB = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
break;
- case kDim3:
+ case kDim3:
d.xyzDim[3] = dcmStrInt(lLength, &buffer[lPos]);
break;
- case kSamplesPerPixel:
+ case kSamplesPerPixel:
d.samplesPerPixel = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
break;
- case kDim2:
+ case kDim2:
d.xyzDim[2] = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
break;
- case kDim1:
+ case kDim1:
d.xyzDim[1] = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
break;
- case kXYSpacing:
+ case kXYSpacing:
dcmMultiFloat(lLength, (char*)&buffer[lPos], 2, d.xyzMM);
break;
- case kImageComments:
+ case kImageComments:
dcmStr (lLength, &buffer[lPos], d.imageComments);
break;
- case kLocationsInAcquisitionGE:
+ case kLocationsInAcquisitionGE:
locationsInAcquisitionGE = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
break;
case kDoseCalibrationFactor :
d.doseCalibrationFactor = dcmStrFloat(lLength, &buffer[lPos]);
break;
- case kBitsAllocated :
+ case kBitsAllocated :
d.bitsAllocated = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
break;
- case kBitsStored :
+ case kBitsStored :
d.bitsStored = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
break;
- case kIsSigned : //http://dicomiseasy.blogspot.com/2012/08/chapter-12-pixel-data.html
+ case kIsSigned : //http://dicomiseasy.blogspot.com/2012/08/chapter-12-pixel-data.html
d.isSigned = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
break;
- case kTR :
+ case kTR :
d.TR = dcmStrFloat(lLength, &buffer[lPos]);
break;
- case kTE :
+ case kTE :
d.TE = dcmStrFloat(lLength, &buffer[lPos]);
break;
- case kTI :
+ case kTI :
d.TI = dcmStrFloat(lLength, &buffer[lPos]);
break;
case kEchoNum :
@@ -3054,7 +3221,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
case kMagneticFieldStrength :
d.fieldStrength = dcmStrFloat(lLength, &buffer[lPos]);
break;
- case kZSpacing :
+ case kZSpacing :
zSpacing = dcmStrFloat(lLength, &buffer[lPos]);
break;
case kPhaseEncodingSteps :
@@ -3095,7 +3262,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
d.TE = dcmStrFloat(lLength, &buffer[lPos]);
}
break;
- case kSlope :
+ case kSlope :
d.intenScale = dcmStrFloat(lLength, &buffer[lPos]);
break;
case kPhilipsSlope :
@@ -3103,13 +3270,13 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
d.intenScalePhilips = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);
break;
- case kIntercept :
+ case kIntercept :
d.intenIntercept = dcmStrFloat(lLength, &buffer[lPos]);
break;
- case kZThick :
+ case kZThick :
d.xyzMM[3] = dcmStrFloat(lLength, &buffer[lPos]);
break;
- case kCoilSiemens : {
+ case kCoilSiemens : {
if (d.manufacturer == kMANUFACTURER_SIEMENS) {
//see if image from single coil "H12" or an array "HEA;HEP"
char coilStr[kDICOMStr];
@@ -3131,7 +3298,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
if (*ptr != '\0')
d.accelFactPE = 0.0;
break; }
- case kLocationsInAcquisition :
+ case kLocationsInAcquisition :
d.locationsInAcquisition = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
break;
case kIconImageSequence:
@@ -3142,7 +3309,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
printMessage("StackSliceNumber %d\n",stackSliceNumber);
break;
}*/
- case kNumberOfDynamicScans:
+ case kNumberOfDynamicScans:
d.numberOfDynamicScans = dcmStrInt(lLength, &buffer[lPos]);
break;
case kMRAcquisitionType: //detect 3D acquisition: we can reorient these without worrying about slice time correct or BVEC/BVAL orientation
@@ -3283,11 +3450,11 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
d.CSA.dtiV[d.CSA.numDti-1][3] = dcmFloat(lLength, &buffer[lPos],d.isLittleEndian);*/
//http://www.na-mic.org/Wiki/index.php/NAMIC_Wiki:DTI:DICOM_for_DWI_and_DTI
break;
- case kWaveformSq:
+ case kWaveformSq:
d.imageStart = 1; //abort!!!
printMessage("Skipping DICOM (audio not image) '%s'\n", fname);
break;
- case kCSAImageHeaderInfo:
+ case kCSAImageHeaderInfo:
readCSAImageHeader(&buffer[lPos], lLength, &d.CSA, isVerbose, dti4D);
d.isHasPhase = d.CSA.isPhaseMap;
break;
@@ -3300,14 +3467,17 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
d.CSA.SeriesHeader_offset = (int)lPos;
d.CSA.SeriesHeader_length = lLength;
break;
- case kRealWorldIntercept:
+ case kRealWorldIntercept:
if (isSameFloat(0.0, d.intenIntercept)) //give precedence to standard value
d.intenIntercept = dcmFloatDouble(lLength, &buffer[lPos],d.isLittleEndian);
break;
- case kRealWorldSlope:
+ case kRealWorldSlope:
if (isSameFloat(1.0, d.intenScale)) //give precedence to standard value
d.intenScale = dcmFloatDouble(lLength, &buffer[lPos],d.isLittleEndian);
break;
+ case kEffectiveEchoSpacingGE:
+ if (d.manufacturer == kMANUFACTURER_GE) d.effectiveEchoSpacingGE = dcmInt(lLength,&buffer[lPos],d.isLittleEndian);
+ break;
case kDiffusionBFactorGE :
if (d.manufacturer == kMANUFACTURER_GE) d.CSA.dtiV[0] = dcmStrInt(lLength, &buffer[lPos]);
break;
@@ -3322,17 +3492,17 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
case kProcedureStepDescription:
dcmStr (lLength, &buffer[lPos], d.procedureStepDescription);
break;
- case kOrientationACR : //use in emergency if kOrientation is not present!
+ case kOrientationACR : //use in emergency if kOrientation is not present!
if (!isOrient) dcmMultiFloat(lLength, (char*)&buffer[lPos], 6, d.orient);
break;
- case kOrientation :
+ case kOrientation :
dcmMultiFloat(lLength, (char*)&buffer[lPos], 6, d.orient);
isOrient = true;
break;
case kImagesInAcquisition :
imagesInAcquisition = dcmStrInt(lLength, &buffer[lPos]);
break;
- case kImageStart:
+ case kImageStart:
//if ((!geiisBug) && (!isIconImageSequence)) //do not exit for proprietary thumbnails
if ((d.compressionScheme == kCompressNone ) && (!isIconImageSequence)) //do not exit for proprietary thumbnails
d.imageStart = (int)lPos + (int)lFileOffset;
@@ -3347,13 +3517,13 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
}
isIconImageSequence = false;
break;
- case kImageStartFloat:
+ case kImageStartFloat:
d.isFloat = true;
if (!isIconImageSequence) //do not exit for proprietary thumbnails
d.imageStart = (int)lPos + (int)lFileOffset;
isIconImageSequence = false;
break;
- case kImageStartDouble:
+ case kImageStartDouble:
printWarning("Double-precision DICOM conversion untested: please provide samples to developer\n");
d.isFloat = true;
if (!isIconImageSequence) //do not exit for proprietary thumbnails
@@ -3387,14 +3557,15 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
free (buffer);
if (encapsulatedDataFragmentStart > 0) {
if (encapsulatedDataFragments > 1)
- printError(" Compressed image stored as %d fragments: decompress with Osirix, dcmdjpeg or dcmjp2k\n", encapsulatedDataFragments);
+ printError("Compressed image stored as %d fragments: decompress with gdcmconv, Osirix, dcmdjpeg or dcmjp2k\n", encapsulatedDataFragments);
else
d.imageStart = encapsulatedDataFragmentStart;
} else if ((isEncapsulatedData) && (d.imageStart < 128)) {
- printWarning(" DICOM violation (contact vendor): compressed image without image fragments, assuming image offset defined by 0x7FE0,x0010\n");
+ //http://www.dclunie.com/medical-image-faq/html/part6.html
+ //Uncompressed data (unencapsulated) is sent in DICOM as a series of raw bytes or words (little or big endian) in the Value field of the Pixel Data element (7FE0,0010). Encapsulated data on the other hand is sent not as raw bytes or words but as Fragments contained in Items that are the Value field of Pixel Data
+ printWarning("DICOM violation (contact vendor): compressed image without image fragments, assuming image offset defined by 0x7FE0,x0010: %s\n", fname);
d.imageStart = encapsulatedDataImageStart;
}
-
//Recent Philips images include DateTime (0008,002A) but not separate date and time (0008,0022 and 0008,0032)
#define kYYYYMMDDlen 8 //how many characters to encode year,month,day in "YYYYDDMM" format
if ((strlen(acquisitionDateTimeTxt) > (kYYYYMMDDlen+5)) && (!isFloatDiff(d.acquisitionTime, 0.0f)) && (!isFloatDiff(d.acquisitionDate, 0.0f)) ) {
@@ -3452,7 +3623,7 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
// d.seriesNum = d.seriesNum + (kEchoMult*d.echoNum);
if ((d.compressionScheme == kCompress50) && (d.bitsAllocated > 8) ) {
//dcmcjpg with +ee can create .51 syntax images that are 8,12,16,24-bit: we can only decode 8/24-bit
- printError("Unable to decode %d-bit images with Transfer Syntax 1.2.840.10008.1.2.4.51, decompress with dcmdjpg\n", d.bitsAllocated);
+ printError("Unable to decode %d-bit images with Transfer Syntax 1.2.840.10008.1.2.4.51, decompress with dcmdjpg or gdcmconv\n", d.bitsAllocated);
d.isValid = false;
}
if ((d.manufacturer == kMANUFACTURER_SIEMENS) && (isMosaic) && (d.CSA.mosaicSlices < 1) && (d.phaseEncodingSteps > 0) && ((d.xyzDim[1] % d.phaseEncodingSteps) == 0) && ((d.xyzDim[2] % d.phaseEncodingSteps) == 0) ) {
@@ -3484,7 +3655,8 @@ struct TDICOMdata readDICOMv(char * fname, int isVerbose, int compressFlag, stru
printMessage("Slices not spatially contiguous: please check output [new feature]\n");
}
if (isVerbose) {
- printMessage("%s\n patient position\t%g\t%g\t%g\n",fname, d.patientPosition[1],d.patientPosition[2],d.patientPosition[3]);
+ printMessage("%s\n patient position (0020,0032)\t%g\t%g\t%g\n",fname, d.patientPosition[1],d.patientPosition[2],d.patientPosition[3]);
+ printMessage("%s\n orient (0020,0037)\t%g\t%g\t%g\t%g\t%g\t%g\n",fname, d.orient[1],d.orient[2],d.orient[3], d.orient[4],d.orient[5],d.orient[6]);
printMessage(" acq %d img %d ser %ld dim %dx%dx%d mm %gx%gx%g offset %d dyn %d loc %d valid %d ph %d mag %d posReps %d nDTI %d 3d %d bits %d littleEndian %d echo %d coil %d TE %g TR %g\n",d.acquNum,d.imageNum,d.seriesNum,d.xyzDim[1],d.xyzDim[2],d.xyzDim[3],d.xyzMM[1],d.xyzMM[2],d.xyzMM[3],d.imageStart, d.numberOfDynamicScans, d.locationsInAcquisition, d.isValid, d.isHasPhase, d.isHasMagnitude,d.patientPositionSequentialRepeats, d.CSA.numDti, d.is3DAcq, d.bitsAllocated, d.isLittle [...]
if (d.CSA.dtiV[0] > 0)
printMessage(" DWI bxyz %g %g %g %g\n", d.CSA.dtiV[0], d.CSA.dtiV[1], d.CSA.dtiV[2], d.CSA.dtiV[3]);
diff --git a/console/nii_dicom.h b/console/nii_dicom.h
index 276ca4d..fcdefc5 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.20170624" kDCMsuf kCCsuf
+ #define kDCMvers "v1.0.20170818" 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,6 +48,18 @@ 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
+#define kMODALITY_CR 1
+#define kMODALITY_CT 2
+#define kMODALITY_MR 3
+#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;
static const int kSliceOrientSag = 2;
@@ -57,6 +69,8 @@ static const int kCompressNone = 0;
static const int kCompressYes = 1;
static const int kCompressC3 = 2; //obsolete JPEG lossless
static const int kCompress50 = 3; //obsolete JPEG lossy
+static const int kCompressRLE = 4; //run length encoding
+
struct TDTI {
float V[4];
int sliceNumberMrPhilips;
@@ -99,7 +113,7 @@ static const int kCompress50 = 3; //obsolete JPEG lossy
struct TDICOMdata {
long seriesNum;
int xyzDim[5];
- int phaseEncodingLines, phaseEncodingSteps, echoTrainLength, patientPositionNumPhilips, coilNum, echoNum, sliceOrient,numberOfDynamicScans, manufacturer, converted2NII, acquNum, imageNum, imageStart, imageBytes, bitsStored, bitsAllocated, samplesPerPixel,patientPositionSequentialRepeats,locationsInAcquisition, compressionScheme;
+ 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];
float orient[7], patientPosition[4], patientPositionLast[4], xyzMM[4], stackOffcentre[4];
float radionuclidePositronFraction, radionuclideTotalDose, radionuclideHalfLife, doseCalibrationFactor; //PET ISOTOPE MODULE ATTRIBUTES (C.8-57)
@@ -123,7 +137,7 @@ static const int kCompress50 = 3; //obsolete JPEG lossy
unsigned char * nii_planar2rgb(unsigned char* bImg, struct nifti_1_header *hdr, int isPlanar);
int isDICOMfile(const char * fname); //0=not DICOM, 1=DICOM, 2=NOTSURE(not part 10 compliant)
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 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);
diff --git a/console/nii_dicom_batch.cpp b/console/nii_dicom_batch.cpp
index 5453bdc..484150d 100755
--- a/console/nii_dicom_batch.cpp
+++ b/console/nii_dicom_batch.cpp
@@ -43,9 +43,9 @@
#include <stdlib.h>
#include <string.h>
#include <sys/stat.h>
-#ifdef myEnableOtsu
- #include "nii_ostu_ml.h" //provide better brain crop, but artificially reduces signal variability in air
-#endif
+//#ifdef myEnableOtsu
+// #include "nii_ostu_ml.h" //provide better brain crop, but artificially reduces signal variability in air
+//#endif
#include <time.h> // clock_t, clock, CLOCKS_PER_SEC
#include "nii_ortho.h"
#if defined(_WIN64) || defined(_WIN32)
@@ -495,6 +495,23 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts,
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" );
@@ -508,6 +525,9 @@ 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 );
if (!opts.isAnonymizeBIDS) {
@@ -650,6 +670,8 @@ void nii_SaveBIDS(char pathoutname[], struct TDICOMdata d, struct TDCMopts opts,
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
@@ -732,6 +754,7 @@ unsigned char * removeADC(struct nifti_1_header *hdr, unsigned char *inImg, bool
bool * nii_SaveDTI(char pathoutname[],int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmList[], struct TDCMopts opts, int sliceDir, struct TDTI4D *dti4D) {
//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....
+ if (opts.isOnlyBIDS) return NULL;
uint64_t indx0 = dcmSort[0].indx; //first volume
int numDti = dcmList[indx0].CSA.numDti;
if (numDti < 1) return NULL;
@@ -1371,14 +1394,20 @@ void nii_saveAttributes (struct TDICOMdata &data, struct nifti_1_header &header,
#else
int nii_saveNII(char * niiFilename, struct nifti_1_header hdr, unsigned char* im, struct TDCMopts opts) {
+ if (opts.isOnlyBIDS) return EXIT_SUCCESS;
hdr.vox_offset = 352;
size_t imgsz = nii_ImgBytes(hdr);
if (imgsz < 1) {
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
+ 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;
}
@@ -1394,6 +1423,10 @@ 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) ) {
+ if ((imgsz+hdr.vox_offset) > kMaxPigz) {
+ printWarning("Saving uncompressed data: image too large for pigz.\n");
+ return EXIT_SUCCESS;
+ }
char command[768];
strcpy(command, "\"" );
strcat(command, opts.pigzname );
@@ -1497,6 +1530,7 @@ int isSameFloatT (float a, float b, float tolerance) {
unsigned char * nii_saveNII3Dtilt(char * niiFilename, struct nifti_1_header * hdr, unsigned char* im, struct TDCMopts opts, float * sliceMMarray, float gantryTiltDeg, int manufacturer ) {
//correct for gantry tilt - http://www.mathworks.com/matlabcentral/fileexchange/24458-dicom-gantry-tilt-correction
+ if (opts.isOnlyBIDS) return im;
if (gantryTiltDeg == 0.0) return im;
struct nifti_1_header hdrIn = *hdr;
int nVox2DIn = hdrIn.dim[1]*hdrIn.dim[2];
@@ -1575,6 +1609,7 @@ unsigned char * nii_saveNII3Dtilt(char * niiFilename, struct nifti_1_header * hd
int nii_saveNII3Deq(char * niiFilename, struct nifti_1_header hdr, unsigned char* im, struct TDCMopts opts, float * sliceMMarray ) {
//convert image with unequal slice distances to equal slice distances
//sliceMMarray = 0.0 3.0 6.0 12.0 22.0 <- ascending distance from first slice
+ if (opts.isOnlyBIDS) return EXIT_SUCCESS;
int nVox2D = hdr.dim[1]*hdr.dim[2];
if ((nVox2D < 1) || (hdr.dim[0] != 3) ) return EXIT_FAILURE;
if ((hdr.datatype != DT_UINT8) && (hdr.datatype != DT_RGB24) && (hdr.datatype != DT_INT16)) {
@@ -1724,6 +1759,7 @@ void smooth1D(int num, double * im) {
int nii_saveCrop(char * niiFilename, struct nifti_1_header hdr, unsigned char* im, struct TDCMopts opts) {
//remove excess neck slices - assumes output of nii_setOrtho()
+ if (opts.isOnlyBIDS) return EXIT_SUCCESS;
int nVox2D = hdr.dim[1]*hdr.dim[2];
if ((nVox2D < 1) || (fabs(hdr.pixdim[3]) < 0.001) || (hdr.dim[0] != 3) || (hdr.dim[3] < 128)) return EXIT_FAILURE;
if ((hdr.datatype != DT_INT16) && (hdr.datatype != DT_UINT16)) {
@@ -1733,13 +1769,16 @@ int nii_saveCrop(char * niiFilename, struct nifti_1_header hdr, unsigned char* i
short * im16 = ( short*) im;
unsigned short * imu16 = (unsigned short*) im;
float kThresh = 0.09; //more than 9% of max brightness
- #ifdef myEnableOtsu
+ /*#ifdef myEnableOtsu
+ // This code removes noise "haze" yielding better volume rendering and smaller gz compressed files
+ // However, it may disrupt "differennce of gaussian" segmentation estimates
+ // Therefore feature was removed from dcm2niix, which aims for lossless conversion
kThresh = 0.0001;
if (hdr.datatype == DT_UINT16)
maskBackgroundU16 (imu16, hdr.dim[1],hdr.dim[2],hdr.dim[3], 5,2, true);
else
maskBackground16 (im16, hdr.dim[1],hdr.dim[2],hdr.dim[3], 5,2, true);
- #endif
+ #endif*/
int ventralCrop = 0;
//find max value for each slice
int slices = hdr.dim[3];
@@ -1985,7 +2024,7 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis
return EXIT_FAILURE;
}
int sliceDir = 0;
- if (hdr0.dim[3] > 1)sliceDir = headerDcm2Nii2(dcmList[dcmSort[0].indx],dcmList[dcmSort[nConvert-1].indx] , &hdr0);
+ if (hdr0.dim[3] > 1)sliceDir = headerDcm2Nii2(dcmList[dcmSort[0].indx],dcmList[dcmSort[nConvert-1].indx] , &hdr0, true);
//UNCOMMENT NEXT TWO LINES TO RE-ORDER MOSAIC WHERE CSA's protocolSliceNumber does not start with 1
if (dcmList[dcmSort[0].indx].CSA.protocolSliceNumber1 > 1) {
printWarning("Weird CSA 'ProtocolSliceNumber' (%d): SPATIAL, SLICE-ORDER AND DTI TRANSFORMS UNTESTED\n", dcmList[dcmSort[0].indx].CSA.protocolSliceNumber1);
@@ -1997,6 +2036,12 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis
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]);
+ if (opts.isOnlyBIDS) {
+ //note we waste time loading every image, however this ensures hdr0 matches actual output
+ free(imgM);
+ return EXIT_SUCCESS;
+ }
+
nii_SaveText(pathoutname, dcmList[dcmSort[0].indx], opts, &hdr0, nameList->str[indx]);
bool * isADC = nii_SaveDTI(pathoutname,nConvert, dcmSort, dcmList, opts, sliceDir, dti4D);
if ((hdr0.datatype == DT_UINT16) && (!dcmList[dcmSort[0].indx].isSigned)) nii_check16bitUnsigned(imgM, &hdr0);
@@ -2022,7 +2067,7 @@ int saveDcm2Nii(int nConvert, struct TDCMsort dcmSort[],struct TDICOMdata dcmLis
int returnCode = EXIT_FAILURE;
//printMessage(" x--> %d ----\n", nConvert);
if (! opts.isRGBplanar) //save RGB as packed RGBRGBRGB... instead of planar RRR..RGGG..GBBB..B
- imgM = nii_planar2rgb(imgM, &hdr0, true);
+ imgM = nii_planar2rgb(imgM, &hdr0, true); //NIfTI is packed while Analyze was planar
if ((hdr0.dim[4] > 1) && (saveAs3D))
returnCode = nii_saveNII3D(pathoutname, hdr0, imgM,opts);
else {
@@ -2314,7 +2359,7 @@ int convert_parRec(struct TDCMopts opts) {
dcmList[0] = nii_readParRec(nameList.str[0], opts.isVerbose, &dti4D);
struct TDCMsort dcmSort[1];
dcmSort[0].indx = 0;
- saveDcm2Nii(1, dcmSort, dcmList, &nameList, opts, &dti4D);
+ int ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, opts, &dti4D);
free(dcmList);//if (nConvertTotal == 0)
if (nameList.numItems < 1)
printMessage("No valid PAR/REC files were found\n");
@@ -2322,8 +2367,7 @@ int convert_parRec(struct TDCMopts opts) {
for (int i = 0; i < nameList.numItems; i++)
free(nameList.str[i]);
free(nameList.str);
-
- return EXIT_SUCCESS;
+ return ret;
}// convert_parRec()
void freeNameList(struct TSearchList nameList) {
@@ -2405,7 +2449,7 @@ int nii_loadDir(struct TDCMopts* opts) {
if (nameList.numItems < 1) {
printError("Unable to find any DICOM images in %s\n", opts->indir);
free(nameList.str); //ignore compile warning - memory only freed on first of 2 passes
- return EXIT_FAILURE;
+ return kEXIT_NO_VALID_FILES_FOUND;
}
size_t nDcm = nameList.numItems;
printMessage( "Found %lu DICOM image(s)\n", nameList.numItems);
@@ -2414,6 +2458,7 @@ int nii_loadDir(struct TDCMopts* opts) {
struct TDTI4D dti4D;
int nConvertTotal = 0;
bool compressionWarning = false;
+ bool convertError = false;
for (int i = 0; i < 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
@@ -2421,8 +2466,11 @@ int nii_loadDir(struct TDCMopts* opts) {
dcmSort[0].indx = i;
dcmSort[0].img = ((uint64_t)dcmList[i].seriesNum << 32) + dcmList[i].imageNum;
dcmList[i].converted2NII = 1;
- saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, &dti4D);
- nConvertTotal++;
+ int ret = saveDcm2Nii(1, dcmSort, dcmList, &nameList, *opts, &dti4D);
+ if (ret == EXIT_SUCCESS)
+ nConvertTotal++;
+ else
+ convertError = true;
}
if ((dcmList[i].compressionScheme != kCompressNone) && (!compressionWarning) && (opts->compressFlag != kCompressNone)) {
compressionWarning = true; //generate once per conversion rather than once per image
@@ -2457,7 +2505,6 @@ int nii_loadDir(struct TDCMopts* opts) {
opts->series.push_back(nextSeries);
}
}
-
// To avoid a spurious warning below
nConvertTotal = nDcm;
} else {
@@ -2489,8 +2536,11 @@ int nii_loadDir(struct TDCMopts* opts) {
nConvert = removeDuplicatesVerbose(nConvert, dcmSort, &nameList);
else
nConvert = removeDuplicates(nConvert, dcmSort);
- nConvertTotal += nConvert;
- saveDcm2Nii(nConvert, dcmSort, dcmList, &nameList, *opts, &dti4D);
+ int ret = saveDcm2Nii(nConvert, dcmSort, dcmList, &nameList, *opts, &dti4D);
+ if (ret == EXIT_SUCCESS)
+ nConvertTotal += nConvert;
+ else
+ convertError = true;
free(dcmSort);
}//convert all images of this series
}
@@ -2498,14 +2548,13 @@ int nii_loadDir(struct TDCMopts* opts) {
}
#endif
free(dcmList);
+ freeNameList(nameList);
+ if (convertError)
+ return EXIT_FAILURE; //at least one image failed to convert
if (nConvertTotal == 0) {
printMessage("No valid DICOM files were found\n");
+ return kEXIT_NO_VALID_FILES_FOUND;
}
- freeNameList(nameList);
- //if (nameList.numItems > 0)
- // for (int i = 0; i < nameList.numItems; i++)
- // free(nameList.str[i]);
- //free(nameList.str);
return EXIT_SUCCESS;
}// nii_loadDir()
@@ -2533,6 +2582,46 @@ int nii_loadDir(struct TDCMopts* opts) {
strcpy(name,""); //not found!
}*/
+#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;
+}
+#endif
+
void readFindPigz (struct TDCMopts *opts, const char * argv[]) {
#if defined(_WIN64) || defined(_WIN32)
strcpy(opts->pigzname,"pigz.exe");
@@ -2546,39 +2635,59 @@ void readFindPigz (struct TDCMopts *opts, const char * argv[]) {
} else
strcpy(opts->pigzname,".\\pigz"); //drop
#else
- strcpy(opts->pigzname,"/usr/local/bin/pigz");
- char pigz[1024];
- strcpy(pigz, opts->pigzname);
- if (!is_exe(opts->pigzname)) {
- strcpy(opts->pigzname,"/usr/bin/pigz");
- if (!is_exe(opts->pigzname)) {
- strcpy(opts->pigzname,"/usr/local/bin/pigz_mricron");
- if (argv == NULL) { //no exectuable path provided
- if (!is_exe(opts->pigzname))
- strcpy(opts->pigzname,"");
- return;
- }
- if (!is_exe(opts->pigzname)) {
- strcpy(opts->pigzname,argv[0]);
- dropFilenameFromPath(opts->pigzname);//, opts.pigzname);
- char appendChar[2] = {"a"};
- appendChar[0] = kPathSeparator;
- if (opts->pigzname[strlen(opts->pigzname)-1] != kPathSeparator) strcat (opts->pigzname,appendChar);
- strcat(opts->pigzname,"pigz_mricron");
- #if defined(_WIN64) || defined(_WIN32)
- strcat(opts->pigzname,".exe");
- #endif
- if (!is_exe(opts->pigzname)) {
- #ifdef myDisableZLib
- printMessage("Compression requires %s\n",pigz);
- #else //myUseZLib
- printMessage("Compression will be faster with %s\n",pigz);
- #endif
- strcpy(opts->pigzname,"");
- } //no pigz_mricron in exe's folder
- } //no /usr/local/pigz_mricron
- }//no /usr/bin/pigz
- } //no /usr/local/pigz
+ char str[1024];
+ //possible pigz names
+ const char * nams[] = {
+ "pigz",
+ "pigz_mricron",
+ "pigz_afni",
+ };
+ #define n_nam (sizeof (nams) / sizeof (const char *))
+ for (int n = 0; n < n_nam; n++) {
+ if (findpathof(str, nams[n])) {
+ strcpy(opts->pigzname,str);
+ //printMessage("Found pigz: %s\n", str);
+ return;
+ }
+ }
+ //possible pigz paths
+ const char * pths[] = {
+ "/usr/local/bin/",
+ "/usr/bin/",
+ };
+ #define n_pth (sizeof (pths) / sizeof (const char *))
+ char exepth[1024];
+ strcpy(exepth,argv[0]);
+ dropFilenameFromPath(exepth);//, opts.pigzname);
+ char appendChar[2] = {"a"};
+ 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++) {
+ //printf ("%d: %s\n", i, nams[n]);
+ for (int p = 0; p < n_pth; p++) {
+ strcpy(str, pths[p]);
+ strcat(str, nams[n]);
+ //printf (">>>%s\n", str);
+ if (is_exe(str))
+ goto pigzFound;
+ } //p
+ //check exepth
+ strcpy(str, exepth);
+ strcat(str, nams[n]);
+ if (is_exe(str))
+ goto pigzFound;
+ } //n
+ //Failure:
+ #ifdef myDisableZLib
+ printMessage("Compression requires 'pigz' to be installed\n");
+ #else //myUseZLib
+ printMessage("Compression will be faster with 'pigz' installed\n");
+ #endif
+ return;
+ pigzFound: //Success
+ strcpy(opts->pigzname,str);
+ //printMessage("Found pigz: %s\n", str);
#endif
} //readFindPigz()
@@ -2605,8 +2714,9 @@ void setDefaultOpts (struct TDCMopts *opts, const char * argv[]) { //either "set
opts->isGz = false;
opts->gzLevel = -1;
opts->isFlipY = true; //false: images in raw DICOM orientation, true: image rows flipped to cartesian coordinates
- opts->isRGBplanar = false;
+ opts->isRGBplanar = false; //false for NIfTI (RGBRGB...), true for Analyze (RRR..RGGG..GBBB..B)
opts->isCreateBIDS = true;
+ opts->isOnlyBIDS = false;
#ifdef myNoAnonymizeBIDS
opts->isAnonymizeBIDS = false;
#else
diff --git a/console/nii_dicom_batch.h b/console/nii_dicom_batch.h
index d15459c..533b6fa 100644
--- a/console/nii_dicom_batch.h
+++ b/console/nii_dicom_batch.h
@@ -23,7 +23,7 @@ extern "C" {
#endif
struct TDCMopts {
- bool isGz, isFlipY, isCreateBIDS, isAnonymizeBIDS, isCreateText, isIgnoreDerivedAnd2D, isPhilipsFloatNotDisplayScaling, isTiltCorrect, isRGBplanar, isOnlySingleFile, isForceStackSameSeries, isCrop;
+ bool isGz, isFlipY, isCreateBIDS, 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];
#ifdef HAVE_R
diff --git a/console/nii_ostu_ml.cpp b/console/nii_ostu_ml.cpp
deleted file mode 100644
index 799cc7a..0000000
--- a/console/nii_ostu_ml.cpp
+++ /dev/null
@@ -1,628 +0,0 @@
-
-#include "nii_ostu_ml.h"
-#include <stdint.h>
-#include <math.h>
-#include <stdlib.h> // pulls in declaration of malloc, free
-#include <string.h> // pulls in declaration for strlen.
-#include <stdio.h>
-
-//Otsu's Method [Multilevel]
-//Otsu N (1979) A threshold selection method from gray-level histograms. IEEE Trans. Sys., Man., Cyber. 9: 62-66.
-//Lookup Tables as suggested by Liao, Chen and Chung (2001) A fast algorithm for multilevel thresholding
-//note my "otsu" method is slightly faster and much simpler than otsu_ml if you only want bi-level output
-
-
-//typedef double HistoRAd[256];
-typedef double Histo2D[256][256];
-
-void otsuLUT(int *H, Histo2D result) //H is histogram 0..255
-{
- for (int u = 0; u < 256; u++)
- for (int v = 0; v < 256; v++)
- result[u][v] = 0;
- double sum,prob;
- sum = 0;
- for (int v = 0; v < 256; v++)
- sum = sum + H[v];
- if (sum <= 0) return;
- Histo2D P,S;
- P[0][0] = H[0];
- S[0][0] = H[0];
- for (int v = 1; v < 256; v++) {
- prob = H[v]/sum;
- P[0][v] = P[0][v-1]+prob;
- S[0][v] = S[0][v-1]+(v+1)*prob;
- }
- for (int u = 1; u < 256; u++) {
- for (int v = u; v < 256; v++) {
- P[u][v] = P[0][v]-P[0][u-1];
- S[u][v] = S[0][v]-S[0][u-1];
- }
- }
-
- //result is eq 29 from Liao
- for (int u = 0; u < 256; u++) {
- for (int v = u; v < 256; v++) {
- if (P[u][v] == 0) //avoid divide by zero errors...
- result[u][v] = 0;
- else
- result[u][v] = pow(S[u][v],2) /P[u][v];
- }
- }
-}
-
-
-void otsuCostFunc4(int *lHisto, int *low, int *middle, int *high)
-//call with lHisto &lo &hi
-{
- double v,max;
- Histo2D h2d;
- otsuLUT(lHisto, h2d);
- int lo = 12;
- int mid = 128;
- int hi = 192;
- max = h2d[0][lo]+h2d[lo+1][mid]+h2d[mid+1][hi]+h2d[hi+1][255];
- //exhaustively search
- for (int l = 0; l < (255-2); l++) {
- for (int m = l+1; m < (255-1); m++ ) {
- for (int h = m+1; h < (255); h++) {
- v = h2d[0][l]+h2d[l+1][m]+h2d[m+1][h]+h2d[h+1][255];
- if (v > max) {
- lo = l;
- mid = m;
- hi = h;
- max = v;
- }//new max
- }//for h -> hi
- }//for mid
- }//for l -> low
- *low = lo;
- *middle = mid;
- *high = hi;
-}//quad OtsuCostFunc4
-
-int otsuCostFunc2(int *lHisto)
-{
- double v,max;
- Histo2D h2d;
- otsuLUT(lHisto, h2d);
- //default solution - 128 is middle intensity of 0..255
- int n = 128;
- max = h2d[0][n]+h2d[n+1][255];
- int result = n;
- //exhaustively search
- for (n = 0; n < 255; n++) {
- v = h2d[0][n]+h2d[n+1][255];
- if (v > max) {
- result = n;
- max = v;
- }//max
- } //for n
- return result;
-} //bilevel otsuCostFunc2
-
-void otsuCostFunc3(int *lHisto, int *low, int *high)
-//call with lHisto &lo &hi
-{
- double v,max;
- Histo2D h2d;
- otsuLUT(lHisto, h2d);
- int lo = 85;
- int hi = 170;
- //default solution
- max = h2d[0][lo]+h2d[lo+1][hi]+h2d[hi+1][255];
- //exhaustively search
- for (int l = 0; l < (255-1); l++) {
- for (int h = l+1; h < (255); h++ ) {
- v = h2d[0][l]+h2d[l+1][h]+h2d[h+1][255];
- if (v > max) {
- lo = l;
- hi = h;
- max = v;
- } //new max
- }//for h -> hi
- }//for l -> low
- *low = lo;
- *high = hi;
-}; //trilevel OtsuCostFunc3
-
-int findOtsu2 (uint8_t *img8bit, int nVox)
-{
- int lHisto[256];
- if (nVox < 1) return 128;
- //create histogram
- for (int n = 0; n < 256; n++)
- lHisto[n] = 0;
- for (int n = 0; n < nVox; n++)
- lHisto[img8bit[n]]++;
- //now find minimum intraclass variance....
- return otsuCostFunc2(lHisto);
-}
-
-void findOtsu3 (uint8_t *img8bit, int nVox, int* lo, int* hi)
-{
- int lHisto[256];
- *lo = 85;
- *hi = 170;
- if (nVox < 1) return;
- //create histogram
- for (int n = 0; n < 256; n++)
- lHisto[n] = 0;
- for (int n = 0; n < nVox; n++)
- lHisto[img8bit[n]]++;
- otsuCostFunc3(lHisto,lo,hi);
-}
-
-void findOtsu4 (uint8_t *img8bit, int nVox, int* lo, int* mid, int* hi)
-{
- int lHisto[256];
- *lo = 64;
- *mid = 128;
- *hi = 192;
- if (nVox < 1) return;
- //create histogram
- for (int n = 0; n < 256; n++)
- lHisto[n] = 0;
- for (int n = 0; n < nVox; n++)
- lHisto[img8bit[n]]++;
- otsuCostFunc4(lHisto,lo,mid, hi);
-}
-
-void applyOtsu2 (uint8_t *img8bit, int nVox)
-{
- if (nVox < 1) return;
- int result = findOtsu2(img8bit,nVox);
- for (int n = 0; n < nVox; n++) {
- if (img8bit[n] > result)
- img8bit[n] = 255;
- else
- img8bit[n] = 0;
- }
- return;
-}//applyOtsu2
-
-void applyOtsu3 (uint8_t *img8bit, int nVox)
-{
- if (nVox < 1) return;
- int lo, hi;
- findOtsu3(img8bit,nVox, &lo, &hi);
- for (int n = 0; n < nVox; n++) {
- if (img8bit[n] <= lo)
- img8bit[n] = 0;
- else if (img8bit[n] >= hi)
- img8bit[n] = 255;
- else
- img8bit[n] = 128;
- }
- return;
-}//applyOtsu3
-
-void applyOtsu4 (uint8_t *img8bit, int nVox)
-{
- if (nVox < 1) return;
- int lo, mid, hi;
- findOtsu4(img8bit,nVox, &lo, &mid, &hi);
- for (int n = 0; n < nVox; n++) {
- if (img8bit[n] <= lo)
- img8bit[n] = 0;
- else if (img8bit[n] <= mid)
- img8bit[n] = 85;
- else if (img8bit[n] <= hi)
- img8bit[n] = 170;
- else
- img8bit[n] = 255;
- }
- return;
-}//applyOtsu4
-
-void applyOtsuBinary (uint8_t *img8bit, int nVox, int levels)
-//1=1/4, 2=1/3, 3=1/2, 4=2/3, 5=3/4
-{
- if (nVox < 1) return;
- //for (int n = 0; n < nVox; n++)
- // if (img8bit[n] < 10)
- // img8bit[n] = 0;
- if ((levels <= 1) || (levels >= 5)) {
- applyOtsu4(img8bit,nVox);
- } else if ((levels = 2) || (levels = 4))
- applyOtsu3(img8bit,nVox);
- else //level = 3
- applyOtsu2(img8bit,nVox);
- int H[256];
- if (levels <= 3) { //make dark: all except 255 equal 0
- for (int n = 0; n < 256; n++)
- H[n] = 0;
- H[255] = 255;
- } else { //make bright: all except 0 equal 255
- for (int n = 0; n < 256; n++)
- H[n] = 255;
- H[0] = 0;
- }
- for (int n = 0; n < nVox; n++)
- img8bit[n] = H[img8bit[n]];
-}
-
-void countClusterSize (uint8_t *lImg, int64_t *lClusterBuff, int lXi, int lYi, int lZi, uint8_t lClusterValue)
-{
- if ((lXi < 5) || (lYi < 5) || (lZi < 3)) return;
- int const kUndeterminedSize = -1;
- int const kCurrentSize = -2;
- //int lXY = lXi*lYi; //offset one slice
- int lXYZ =lXi*lYi*lZi;
- int64_t *lQra = (int64_t *) malloc(lXYZ * sizeof(int64_t));
- //initialize clusterbuffer: voxels either non part of clusters (cluster size 0), or part of cluster of undetermined size
- for (int v = 0; v < lXYZ; v++) {
- if (lImg[v] == lClusterValue)
- lClusterBuff[v] = kUndeterminedSize; //target voxel - will be part of a cluster of size 1..XYZ
- else
- lClusterBuff[v] = 0;//not a target, not part of a cluser, size = 0
- }//for lInc
- int const nSearch = 6;
- int64_t *searchArray = (int64_t *) malloc(nSearch * sizeof(int64_t));
- searchArray[0] = -1; //left
- searchArray[1] = 1; //right
- searchArray[2] = -lXi; //posterior
- searchArray[3] = lXi; //anterior
- searchArray[4] = -lXi*lYi; //inferior
- searchArray[5] = lXi*lYi; //superior
- for (int v = 0; v < lXYZ; v++) {
- if (lClusterBuff[v] == kUndeterminedSize) {
- int lQCurr = 0;
- int lQMax = 0;
- lQra[lQCurr] = v;
- lClusterBuff[v] = kCurrentSize;
- do {
- for (int s = 0; s < nSearch; s++) {//search each neighbor to see if it is part of this cluster
- int64_t vx = lQra[lQCurr] + searchArray[s];
- if ((vx >= 0) && (vx < lXYZ) && (lClusterBuff[vx] == kUndeterminedSize) ) {
- lClusterBuff[vx] = kCurrentSize;
- lQMax++;
- lQra[lQMax] = vx;
- }
- }
- lQCurr++;
- } while (lQMax >= lQCurr);
- //set all voxels in this cluster to reflect their size
- for (int vx = 0; vx < lQCurr; vx++)
- lClusterBuff[ lQra[vx] ] = lQCurr;
- }//if v part of new cluser
- } //for each voxel
- free(searchArray);
- // selectClusters (lClusterBuff, -1, lXi, lXYZ, lQSz, lTopSlice, lXY, lQra);
- free(lQra);
-} //CountClusterSize
-
-
- //fillBubbles works but is typically not required - increases processing time ~20%
-void fillBubbles (uint8_t *lImg, int lXi, int lYi, int lZi)
-//output voxels only zero if they are connected to edge voxels by zeros....
-{
- if ((lXi < 5) or (lYi < 5) or (lZi < 1)) return;
- int lXYZ = lXi *lYi * lZi;
- int const kUntouchedZero = 10;
- int const nSearch = 6;
- int64_t *searchArray = (int64_t *) malloc(nSearch * sizeof(int64_t));
- int64_t *lQra = (int64_t *) malloc(lXYZ * sizeof(int64_t));
- searchArray[0] = -1; //left
- searchArray[1] = 1; //right
- searchArray[2] = -lXi; //posterior
- searchArray[3] = lXi; //anterior
- searchArray[4] = -lXi*lYi; //inferior
- searchArray[5] = lXi*lYi; //superior
- for (int v = 0; v < lXYZ; v++) {//make binary
- if (lImg[v] != 0)
- lImg[v] = 255;
- else
- lImg[v] = kUntouchedZero; //zero - but not yet sure if connected to edge...
- } //for v each voxel
- for (int z = 1; z <= lZi; z++) {
- int zpos = (z-1) * lXi * lYi;
- for (int y = 1; y <= lYi; y++) {
- int ypos = (y-1) * lXi;
- for (int x = 1; x <= lXi; x++) {
- if ((x==1) || (y==1) || (z==1) || (x==lXi) || (y==lYi) || (z==lZi) ) { //voxel on edge
- int v = zpos + ypos + x -1;
- if (lImg[v] == kUntouchedZero) {
- int lQCurr = 0;
- int lQMax = 0;
- lQra[lQCurr] = v;
- lImg[v] = 0;
- do {
- for (int s = 0; s < nSearch; s++) {//search each neighbor to see if it is part of this cluster
- int64_t vx = lQra[lQCurr] + searchArray[s];
- if ((vx >= 0) && (vx < lXYZ) && (lImg[vx] == kUntouchedZero) ) {
- lImg[vx] = 0;
- lQMax++;
- lQra[lQMax] = vx;
- }
- }
- lQCurr++;
- } while (lQMax >= lQCurr);
- }//untouched zero
- } //on edge
- } //for x
- } //for y
- } //for z
- free(searchArray);
- free(lQra);
- for (int v = 0; v < lXYZ; v++) // untouched voxels not connected to any edge
- if (lImg[v] == kUntouchedZero)
- lImg[v] = 255;
-} //fillBubbles - fill holes inside object
-
-
-void preserveLargestCluster (uint8_t *lImg, int lXi, int lYi, int lZi, int lClusterValue, uint8_t ValueForSmallClusters)
-{
- if ((lXi < 5) or (lYi < 5) or (lZi < 1)) return;
- int lXYZ,lX;
- int64_t lC;
- lXYZ =lXi*lYi*lZi;
- //ensure at least some voxels exist with clusterValue
- lC = 0;
- for (lX =0; lX< lXYZ; lX++)
- if (lImg[lX] == lClusterValue)
- lC++;
- if (lC < 2)
- return;//e.g. if lC = 1 then only a single voxel, which is in fact largest cluster
- int64_t *lTemp = (int64_t *) malloc(lXYZ * sizeof(int64_t));
- for (lX = 0; lX < lXYZ; lX++) lTemp[lX] = 0; //the only purpose of this loop is to hide a compiler warning - countClusterSize will fill this array...
- countClusterSize(lImg,lTemp,lXi,lYi,lZi,lClusterValue);
- lC = 0;
- for (lX = 0; lX < lXYZ; lX++)
- if (lTemp[lX] > lC)
- lC = lTemp[lX]; //volume of largest cluster
- if (ValueForSmallClusters == 0) {
- for (lX = 0; lX < lXYZ; lX++)
- if ((lTemp[lX] >= 0) && (lTemp[lX] < lC)) //cluster, but not biggest one...
- lImg[lX] = ValueForSmallClusters;
- } else {
- for (lX = 0; lX < lXYZ; lX++)
- if ((lTemp[lX] > 0) && (lTemp[lX] < lC)) //cluster, but not biggest one...
- lImg[lX] = ValueForSmallClusters;
- }
- free(lTemp);
-}
-
-void dilate (uint8_t *lImg, int lXi, int lYi, int lZi, int lCycles, int lChange)
-//Dilates Diamonds - neighbor coefficient = 0
-//Dilate if Change=1 then all voxels where intensity <> 1 but where any neighbors = 1 will become 1
-//Erode if Change=0 then all voxels where intensity <>0 but where any neighbors = 0 will become 0
-//step is repeated for lCycles
-{
- if ((lXi < 5) || (lYi < 5) || (lZi < 1))
- return;
- int lX,lY,lZ, lXY,lXYZ,lPos,lOffset;
- lXY = lXi*lYi; //offset one slice
- lXYZ =lXY*lZi;
- uint8_t *lTemp = (uint8_t *) malloc(lXYZ);
- for (int lC = 0; lC < lCycles; lC++) {
- memcpy (lTemp, lImg, lXYZ);
- for (lZ = 0; lZ < lZi; lZ++) {
- for (lY = 0; lY < lYi; lY++) {
- lOffset = (lY*lXi) + (lZ * lXY);
- for (lX = 0; lX < lXi; lX++) {
- lPos = lOffset + lX;
- if (lTemp[lPos] != lChange) {
- if ((lX>0) && (lTemp[lPos-1] == lChange))
- lImg[lPos] = lChange;
- else if ((lX<(lXi-1)) && (lTemp[lPos+1] == lChange))
- lImg[lPos] = lChange;
- else if ((lY>0) && (lTemp[lPos-lXi] == lChange))
- lImg[lPos] = lChange;
- else if ((lY<(lYi-1)) && (lTemp[lPos+lXi] == lChange))
- lImg[lPos] = lChange;
- else if ((lZ>0) && (lTemp[lPos-lXY] == lChange))
- lImg[lPos] = lChange;
- else if ((lZ<(lZi-1)) && (lTemp[lPos+lXY] == lChange))
- lImg[lPos] = lChange;
- } //if not lChange
- } //for X
- } //for Y
- } //for Z
- } //for each cycle
- free(lTemp);
-} //dilate
-
-void dilateSphere (uint8_t *lImg, int lXi, int lYi, int lZi, float lVoxDistance, int lChange)
-//INPUT: Img is array of bytes 0..(X*Y*Z-1) that represents 3D volume, lXi,lYi,lZi are number of voxels in each dimension
-// lVoxDistance is search radius (in voxels)
-// lChange is the intensity to be changed - if background color: erosion, if foreground color: dilation
-//OUTPUT: Eroded/Dilated Img
-{
- if (lVoxDistance < 1) return;
- if (lVoxDistance == 1) { //much faster to use classic neighbor dilation
- dilate(lImg,lXi,lYi,lZi,1,lChange);
- return;
- }
- if ((lXi < 3) || (lYi < 3) || (lZi < 3))
- return;
- float lDx;
- int lXY = lXi*lYi; //voxels per slice
- int lXYZ = lXi*lYi*lZi; //voxels per volume
- //next: make 1D array of all voxels within search sphere: store offset from center
- int lDxI = trunc(lVoxDistance); //no voxel will be searched further than DxI from center
- int64_t *lSearch = (int64_t *) malloc( ((lDxI *2)+1)*((lDxI *2)+1)*((lDxI *2)+1) *sizeof(int64_t) );
- int lVoxOK = 0;
- for (int lZ = -lDxI; lZ <= lDxI; lZ++) {
- for (int lY = -lDxI; lY <= lDxI; lY++) {
- for (int lX = -lDxI; lX <= lDxI; lX++) {
- lDx = sqrt( (lX*lX)+ (lY*lY)+ (lZ*lZ) );
- if ((lDx < lVoxDistance) && (lDx > 0)) {
- lSearch[lVoxOK] = lX + (lY*lXi) + (lZ * lXY); //offset to center
- lVoxOK++;
- } //in range, not center
- } //lX
- } //lY
- }//lZ
- uint8_t *lTemp = (uint8_t *) malloc(lXYZ);
- memcpy (lTemp, lImg, lXYZ);
- for (int lVox = 0; lVox < lXYZ; lVox++)
- if (lTemp[lVox] != lChange) { //voxel not a survivor
- for (int s = 0; s < lVoxOK; s++) { //make voxel a survivor if neighbor is a survivor
- int64_t t = lVox +lSearch[s];
- if ((t >= 0) && (t < lXYZ) && (lTemp[t] == lChange)) {
- lImg[lVox] = lChange;
- break;
- } //if neighbor
- } //for s: search for neighboring survivrs
- } //if lTemp[lVox] not a survivor
- free(lTemp); //free temporary buffer
- free(lSearch); //free 1D search space
-}//proc DilateSphere
-
-void smoothFWHM2Vox (uint8_t *lImg, int lXi, int lYi, int lZi)
-{
- if ((lXi < 5) || (lYi < 5) || (lZi < 1)) return;
- int k0=240;//weight of center voxel
- int k1=120;//weight of nearest neighbors
- int k2=15;//weight of subsequent neighbors
- int kTot=k0+k1+k1+k2+k2; //weight of center plus all neighbors within 2 voxels
- int kWid = 2; //we will look +/- 2 voxels from center
- int lyPos,lPos,lWSum,lXi2,lXY,lXY2;
- lXY = lXi*lYi; //offset one slice
- int lXYZ = lXY * lZi;
- lXY2 = lXY * 2; //offset two slices
- lXi2 = lXi*2;//offset to voxel two lines above or below
- uint8_t *lTemp = (uint8_t *) malloc(lXYZ);
- memcpy (lTemp, lImg, lXYZ);
- //smooth horizontally
- for (int lZ = 0; lZ < lZi; lZ++) {
- for (int lY = 0; lY < lYi; lY++) {
- lyPos = (lY*lXi) + (lZ*lXY) ;
- for (int lX = kWid; lX < (lXi-kWid); lX++) {
- lPos = lyPos + lX;
- lWSum = lImg[lPos-2]*k2 +lImg[lPos-1]*k1 +lImg[lPos]*k0 +lImg[lPos+1]*k1 +lImg[lPos+2]*k2;
- lTemp[lPos] = lWSum / kTot;
- } //for lX
- } //for lY
- } //for lZi
- //smooth vertically
- memcpy (lImg,lTemp, lXYZ);
- for (int lZ = 0; lZ < lZi; lZ++) {
- for (int lX = 0; lX < lXi; lX++) {
- for (int lY = kWid; lY < (lYi-kWid); lY++) {
- lPos = ((lY*lXi) + lX + (lZ*lXY) ) ;
- lWSum = lTemp[lPos-lXi2]*k2+lTemp[lPos-lXi]*k1+lTemp[lPos]*k0+lTemp[lPos+lXi]*k1+lTemp[lPos+lXi2]*k2;
- lImg[lPos] = lWSum / kTot;
- }//for Y
- } //for X
- } //for Z
- //if 3rd dimension....
- if (lZi >= 5) {
- //smooth across slices
- memcpy (lTemp, lImg, lXYZ);
- for (int lZ = kWid; lZ < (lZi-kWid); lZ++) {
- for (int lY = 0; lY < lYi; lY++) {
- lyPos = (lY*lXi) + (lZ*lXY) ;
- for (int lX = 0; lX < lXi; lX++) {
- lPos = lyPos + lX;
- lWSum = lImg[lPos-lXY2]*k2+lImg[lPos-lXY]*k1+lImg[lPos]*k0+lImg[lPos+lXY]*k1+lImg[lPos+lXY2]*k2;
- lTemp[lPos] = lWSum / kTot;
- }//for lX
- }//for lY
- } //for lZi
- memcpy ( lImg, lTemp, lXYZ);
- }//if Z: at least 5 slices...
- free(lTemp);
-} //smoothFWHM2Vox
-
-void maskBackground (uint8_t *img8bit, int lXi, int lYi, int lZi, int lOtsuLevels, float lDilateVox, bool lOneContiguousObject)
-{
- if ((lXi < 3) or (lYi < 3) or (lZi < 1)) return;
- int lXYZ = lXi * lYi * lZi;
- //countSurvivors(img8bit, lXYZ);
- smoothFWHM2Vox(img8bit, lXi,lYi,lZi);
- applyOtsuBinary (img8bit, lXYZ,lOtsuLevels);
- //clip edges to prevent wrap
- int lV=0;
- for (int lZ = 1; lZ <= lZi; lZ++)
- for (int lY = 1; lY <= lYi; lY++)
- for (int lX = 1; lX <= lXi; lX++) {
- if ((lX==1) || (lX==lXi) || (lY==1) || (lY==lYi) || (lZ==1) || (lZ==lZi))
- img8bit[lV] = 0;
- lV++;
- } //for lX
- if (lOneContiguousObject) {
- preserveLargestCluster (img8bit, lXi,lYi,lZi,255,0 ); //only preserve largest single object
- if (lDilateVox > 0)
- dilateSphere (img8bit, lXi,lYi,lZi,lDilateVox,255 );
- } else {
- if (lDilateVox > 0)
- dilateSphere (img8bit, lXi,lYi,lZi,lDilateVox,255 );
- preserveLargestCluster (img8bit, lXi,lYi,lZi,0,255 ); //only erase outside air
- }
- fillBubbles(img8bit, lXi,lYi,lZi); //<- optional
- //my sense is filling air bubbles will help SPM's mixture of gaussians detect air. Without this air will have an artificially low variance
-}
-
-int compare_short (const void *a, const void *b)
-{
- const short *da = (const short *) a;
- const short *db = (const short *) b;
- return (*da > *db) - (*da < *db);
-}
-
-
-void maskBackground16 (short *img16bit, int lXi, int lYi, int lZi, int lOtsuLevels, float lDilateVox, bool lOneContiguousObject) {
- if ((lXi < 3) or (lYi < 3) or (lZi < 1)) return;
- int lXYZ = lXi * lYi * lZi;
- if (lXYZ < 100) return;
- //find image intensity range
- short * imgSort = (short *) malloc(sizeof(short)*lXYZ);
- memcpy(&imgSort[0], &img16bit[0], lXYZ * sizeof(short)); //memcpy( dest, src, bytes)
- qsort (imgSort, lXYZ, sizeof (short), compare_short);
- int pctLo = imgSort[(int)((float)lXYZ * 0.02)];
- int pctHi = imgSort[(int)((float)lXYZ * 0.98)];
- free(imgSort);
- if (pctLo >= pctHi) return; //no variability
- //re-scale data 0..255
- float scale = 255.0/(float)(pctHi - pctLo);
- uint8_t *img8bit = (uint8_t *) malloc(sizeof(uint8_t)*lXYZ);
- for (int i = 0; i < lXYZ; i++) {
- if (img16bit[i] >= pctHi)
- img8bit[i] = 255;
- else if (img16bit[i] <= pctLo)
- img8bit[i] = 0;
- else
- img8bit[i] = (int)((img16bit[i] - pctLo) * scale);
- }
- //mask 8-bit image
- maskBackground (img8bit, lXi, lYi, lZi, lOtsuLevels, lDilateVox, lOneContiguousObject);
- for (int i = 0; i < lXYZ; i++)
- if (img8bit[i] == 0)
- img16bit[i] = 0;
-}
-
-void maskBackgroundU16 (unsigned short *img16bit, int lXi, int lYi, int lZi, int lOtsuLevels, float lDilateVox, bool lOneContiguousObject) {
- if ((lXi < 3) or (lYi < 3) or (lZi < 1)) return;
- int lXYZ = lXi * lYi * lZi;
- if (lXYZ < 100) return;
- //find image intensity range
- unsigned short * imgSort = (unsigned short *) malloc(sizeof(unsigned short)*lXYZ);
- memcpy(&imgSort[0], &img16bit[0], lXYZ * sizeof(short)); //memcpy( dest, src, bytes)
- qsort (imgSort, lXYZ, sizeof (short), compare_short);
- int pctLo = imgSort[(int)((float)lXYZ * 0.02)];
- int pctHi = imgSort[(int)((float)lXYZ * 0.98)];
- free(imgSort);
- if (pctLo >= pctHi) return; //no variability
- //re-scale data 0..255
- float scale = 255.0/(float)(pctHi - pctLo);
- uint8_t *img8bit = (uint8_t *) malloc(sizeof(uint8_t)*lXYZ);
- for (int i = 0; i < lXYZ; i++) {
- if (img16bit[i] >= pctHi)
- img8bit[i] = 255;
- else if (img16bit[i] <= pctLo)
- img8bit[i] = 0;
- else
- img8bit[i] = (int)((img16bit[i] - pctLo) * scale);
- }
- //mask 8-bit image
- maskBackground (img8bit, lXi, lYi, lZi, lOtsuLevels, lDilateVox, lOneContiguousObject);
- for (int i = 0; i < lXYZ; i++)
- if (img8bit[i] == 0)
- img16bit[i] = 0;
-}
-
-
-
-
diff --git a/console/nii_ostu_ml.h b/console/nii_ostu_ml.h
deleted file mode 100644
index 589f1d4..0000000
--- a/console/nii_ostu_ml.h
+++ /dev/null
@@ -1,19 +0,0 @@
-#include <stdint.h>
-
-
-#ifndef nii_otsu_h
-#define nii_otsu_h
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
- void maskBackground ( uint8_t *img8bit, int lXi, int lYi, int lZi, int lOtsuLevels, float lDilateVox, bool lOneContiguousObject);
- void maskBackground16 (short *img16bit, int lXi, int lYi, int lZi, int lOtsuLevels, float lDilateVox, bool lOneContiguousObject);
- void maskBackgroundU16 (unsigned short *img16bit, int lXi, int lYi, int lZi, int lOtsuLevels, float lDilateVox, bool lOneContiguousObject);
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif
\ No newline at end of file
diff --git a/console/print.h b/console/print.h
index 497d3ec..fcc4498 100755
--- a/console/print.h
+++ b/console/print.h
@@ -31,15 +31,22 @@
std::cout << buf;
delete[] buf;
}
+ #define printError(...) do { printMessage("Error: "); printMessage(__VA_ARGS__);} while(0)
#else
#include<stdio.h>
#define printMessage printf
+ //#define printMessageError(...) fprintf (stderr, __VA_ARGS__)
+ #ifdef myErrorStdOut //for XCode MRIcro project, pipe errors to stdout not stderr
+ #define printError(...) do { printMessage("Error: "); printMessage(__VA_ARGS__);} while(0)
+ #else
+ #define printError(...) do { fprintf (stderr,"Error: "); fprintf (stderr, __VA_ARGS__);} while(0)
+ #endif
#endif //myUseCOut
//n.b. use ({}) for multi-line macros http://www.geeksforgeeks.org/multiline-macros-in-c/
//these next lines work on GCC but not _MSC_VER
// #define printWarning(...) ({printMessage("Warning: "); printMessage(__VA_ARGS__);})
// #define printError(...) ({ printMessage("Error: "); printMessage(__VA_ARGS__);})
#define printWarning(...) do {printMessage("Warning: "); printMessage(__VA_ARGS__);} while(0)
- #define printError(...) do { printMessage("Error: "); printMessage(__VA_ARGS__);} while(0)
+
#endif //HAVE_R
#endif //_R_PRINT_H_
--
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