From 84f782e8d23cb95569ac24098515949809d4d855 Mon Sep 17 00:00:00 2001
From: Taylor Hanayik <hanayik@gmail.com>
Date: Fri, 16 Jun 2023 11:11:33 +0100
Subject: [PATCH] add --json option and associated logic
---
.gitignore | 1 +
fslstats.cc | 1263 +++++++++++++++++++++++++++++++++++----------------
2 files changed, 867 insertions(+), 397 deletions(-)
create mode 100644 .gitignore
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..600d2d3
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+.vscode
\ No newline at end of file
diff --git a/fslstats.cc b/fslstats.cc
index 03514a9..83a9873 100644
--- a/fslstats.cc
+++ b/fslstats.cc
@@ -15,16 +15,18 @@
#include "newimage/newimageall.h"
#include "newimage/costfns.h"
-
using namespace std;
using namespace NEWMAT;
using namespace NEWIMAGE;
-int print_usage(const string& progname) {
- cout << "Usage: fslstats [preoptions] <input> [options]" << endl << endl;
+int print_usage(const string &progname)
+{
+ cout << "Usage: fslstats [preoptions] <input> [options]" << endl
+ << endl;
cout << "preoption -t will give a separate output line for each 3D volume of a 4D timeseries" << endl;
cout << "preoption -K < indexMask > will generate seperate n submasks from indexMask, for indexvalues 1..n where n is the maximum index value in indexMask, and generate statistics for each submask" << endl;
- cout << "Note - options are applied in order, e.g. -M -l 10 -M will report the non-zero mean, apply a threshold and then report the new nonzero mean" << endl << endl;
+ cout << "Note - options are applied in order, e.g. -M -l 10 -M will report the non-zero mean, apply a threshold and then report the new nonzero mean" << endl
+ << endl;
cout << "-l <lthresh> : set lower threshold" << endl;
cout << "-u <uthresh> : set upper threshold" << endl;
cout << "-r : output <robust min intensity> <robust max intensity>" << endl;
@@ -44,12 +46,16 @@ int print_usage(const string& progname) {
cout << "-C : output centre-of-gravity (cog) in voxel coordinates" << endl;
cout << "-p <n> : output nth percentile (n between 0 and 100)" << endl;
cout << "-P <n> : output nth percentile (for nonzero voxels)" << endl;
- cout << "-a : use absolute values of all image intensities"<< endl;
+ cout << "-a : use absolute values of all image intensities" << endl;
cout << "-n : treat NaN or Inf as zero for subsequent stats" << endl;
cout << "-k <mask> : use the specified image (filename) for masking - overrides lower and upper thresholds" << endl;
cout << "-d <image> : take the difference between the base image and the image specified here" << endl;
cout << "-h <nbins> : output a histogram (for the thresholded/masked voxels only) with nbins" << endl;
- cout << "-H <nbins> <min> <max> : output a histogram (for the thresholded/masked voxels only) with nbins and histogram limits of min and max" << endl << endl;
+ cout << "-H <nbins> <min> <max> : output a histogram (for the thresholded/masked voxels only) with nbins and histogram limits of min and max" << endl;
+ // add --json output format option
+ // keys will correspond to the options above
+ cout << "--json : output in JSON format" << endl
+ << endl;
cout << "Note - thresholds are not inclusive ie lthresh<allowed<uthresh" << endl;
return 1;
}
@@ -57,461 +63,924 @@ int print_usage(const string& progname) {
// Some specialised nonzero functions just for speedup
// (it avoids generating masks when not absolutely necessary)
-long nonzerocount(const volume4D<float>& vol)
+long nonzerocount(const volume4D<float> &vol)
{
- long totn=0;
- for(volume4D<float>::fast_const_iterator it=vol.fbegin(), end=vol.fend(); it<end; it++ )
- if ( (*it) != 0 )
+ long totn = 0;
+ for (volume4D<float>::fast_const_iterator it = vol.fbegin(), end = vol.fend(); it < end; it++)
+ if ((*it) != 0)
totn++;
return totn;
}
-double nonzeromean(const volume4D<float>& vol)
+double nonzeromean(const volume4D<float> &vol)
{
- double total=0.0;
- long int totn=0;
- for(volume4D<float>::fast_const_iterator it=vol.fbegin(), end=vol.fend(); it<end; it++ )
- if ( (*it) != 0 ) {
- total+=(double)(*it);
+ double total = 0.0;
+ long int totn = 0;
+ for (volume4D<float>::fast_const_iterator it = vol.fbegin(), end = vol.fend(); it < end; it++)
+ if ((*it) != 0)
+ {
+ total += (double)(*it);
totn++;
}
- if (totn>0)
- total/=totn;
+ if (totn > 0)
+ total /= totn;
return total;
}
-double nonzerostddev(const volume4D<float>& vol)
+double nonzerostddev(const volume4D<float> &vol)
{
- double totv=0.0, totvv=0.0;
- long int totn=0;
- for(volume4D<float>::fast_const_iterator it=vol.fbegin(), end=vol.fend(); it<end; it++ )
- if ( (*it) != 0 ) {
- totvv+=(double)((*it)*(*it));
- totv+=(double)(*it);
+ double totv = 0.0, totvv = 0.0;
+ long int totn = 0;
+ for (volume4D<float>::fast_const_iterator it = vol.fbegin(), end = vol.fend(); it < end; it++)
+ if ((*it) != 0)
+ {
+ totvv += (double)((*it) * (*it));
+ totv += (double)(*it);
totn++;
}
- double var=0.0;
- if (totn>1) {
- double meanval = totv/totn;
- var = (totvv - totn*meanval*meanval)/(totn-1);
+ double var = 0.0;
+ if (totn > 1)
+ {
+ double meanval = totv / totn;
+ var = (totvv - totn * meanval * meanval) / (totn - 1);
}
return std::sqrt(var);
}
-template<class M>
+template <class M>
int generateNonZeroMask(const M &mask, volume4D<float> &masknz, const volume4D<float> &input)
{
- masknz = binarise(input,0.0f,0.0f,inclusive,true)*mask;
+ masknz = binarise(input, 0.0f, 0.0f, inclusive, true) * mask;
return 0;
}
-int generate_masks(volume4D<float>& mask, volume4D<float>& masknz, const volume4D<float>& input, const float& lthr, const float& uthr)
+int generate_masks(volume4D<float> &mask, volume4D<float> &masknz, const volume4D<float> &input, const float <hr, const float &uthr)
+{
+ mask = binarise(input, lthr, uthr, exclusive);
+ return generateNonZeroMask(mask, masknz, input);
+}
+
+// print function that can handle outputting in JSON format or not
+void print_value(const string &key, const string &value, string &json, const bool jsonOutput)
{
- mask = binarise(input,lthr,uthr,exclusive);
- return generateNonZeroMask(mask,masknz,input);
+ // append the key and value to the json string
+ if (jsonOutput)
+ {
+ json.append("\"" + key + "\": \"" + value + "\", ");
+ }
+ else
+ {
+ cout << value << " ";
+ }
}
-int fmrib_main_float(int argc, char* argv[],const bool timeseriesMode, const string& indexMaskName)
+
+int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const string &indexMaskName)
{
+ // determine if the user wants to output in JSON format
+ bool jsonOutput = false;
+ for (int i = 1; i < argc; i++)
+ {
+ if (strcmp(argv[i], "--json") == 0)
+ {
+ jsonOutput = true;
+ break;
+ }
+ }
+ // if JSON output is requested, then we need to store the output as a string
+ // and then output it at the end via cout
+ // don't forget to close the json at the end with a "}"
+ string json = "";
+
cout.setf(ios::dec);
cout.setf(ios::fixed, ios::floatfield);
cout.setf(ios::left, ios::adjustfield);
cout.precision(6);
volume<float> vol, inputMaster;
volume<int> indexMask;
- if ( timeseriesMode || indexMaskName != "" )
- read_volume4D(inputMaster,argv[1]);
+ if (timeseriesMode || indexMaskName != "")
+ read_volume4D(inputMaster, argv[1]);
else
- read_volume4D(vol,argv[1]);
- volume<float> & indexMaster = (timeseriesMode ) ? vol : inputMaster;
- if ( indexMaskName != "" )
- read_volume4D(indexMask,indexMaskName);
+ read_volume4D(vol, argv[1]);
+ volume<float> &indexMaster = (timeseriesMode) ? vol : inputMaster;
+ if (indexMaskName != "")
+ read_volume4D(indexMask, indexMaskName);
int nTimepoints((timeseriesMode) ? inputMaster.tsize() : 1), nIndices((indexMaskName != "") ? indexMask.max() : 1);
- for ( int timepoint=0; timepoint < nTimepoints ; timepoint++ ) {
- for ( int index=1; index <= nIndices; index++ ) {
- if ( timeseriesMode )
- vol=inputMaster[timepoint];
- volume<float> mask, masknz;
- float lthr(-numeric_limits<float>::max());
- float uthr(numeric_limits<float>::max());
- if ( indexMask.nvoxels() ) {
- if ( indexMask.dimensionality() > 3 ) {
- cerr << "Index mask must be 3D" << endl;
- return 3;
- }
- copyconvert(indexMask,mask);
- mask.binarise(index-1,index+1,exclusive);
- vol=indexMaster*mask;
- if (mask.max() < 1) {
- cout << "missing label: " << index << " in mask" << endl;
- continue;
- }
- generateNonZeroMask(mask,masknz,vol);
- }
- int narg(2);
-
- while (narg<argc) {
- string sarg(argv[narg]);
- if (sarg=="-n") {
- for (int t=0; t<vol.tsize(); t++)
- for (int z=0; z<vol.zsize(); z++)
- for (int y=0; y<vol.ysize(); y++)
- for (int x=0; x<vol.xsize(); x++)
- if (!isfinite((double)vol(x,y,z,t)))
- vol(x,y,z,t)=0;
- } else if (sarg=="-m") {
- if (mask.nvoxels()>0) cout << vol.mean(mask) << " ";
- else cout << vol.mean() << " ";
- } else if (sarg=="-M") {
- if (masknz.nvoxels()>0) cout << vol.mean(masknz) << " ";
- else {
- double nzmean=0;
- nzmean = nonzeromean(vol);
- cout << nzmean << " ";
- }
- } else if (sarg=="-X") {
- ColumnVector coord(4);
- vector<int64_t> iCoords;
- vol.min(mask,iCoords);
- coord << (Real)iCoords[7] << (Real)iCoords[8] << (Real)iCoords[9] << 1.0;
- coord = vol.niftivox2newimagevox_mat().i() * coord;
- cout << MISCMATHS::round(coord(1)) << " " <<
- MISCMATHS::round(coord(2)) << " " << MISCMATHS::round(coord(3)) << " ";
- } else if (sarg=="-x") {
- ColumnVector coord(4);
- vector<int64_t> iCoords;
- vol.min(mask,iCoords);
- coord << (Real)iCoords[0] << (Real)iCoords[1] << (Real)iCoords[2] << 1.0;
- coord = vol.niftivox2newimagevox_mat().i() * coord;
- cout << MISCMATHS::round(coord(1)) << " " <<
- MISCMATHS::round(coord(2)) << " " << MISCMATHS::round(coord(3)) << " ";
- } else if (sarg=="-w") {
- if (masknz.nvoxels()<1) { //Need to generate non-zeromask
- generate_masks(mask,masknz,vol,lthr,uthr);
- vol*=mask;
- }
- int xmin=masknz.xsize()-1,xmax(0),ymin=masknz.ysize()-1,ymax(0),zmin=masknz.zsize()-1,zmax(0),tmin=masknz.tsize()-1,tmax(0);
- for(int t=0;t<masknz.tsize();t++)
- for(int z=0;z<masknz.zsize();z++)
- for(int y=0;y<masknz.ysize();y++)
- for(int x=0;x<masknz.xsize();x++)
- if (masknz(x,y,z,t)>0.5) {
- if (x<xmin) xmin=x;
- if (x>xmax) xmax=x;
- if (y<ymin) ymin=y;
- if (y>ymax) ymax=y;
- if (z<zmin) zmin=z;
- if (z>zmax) zmax=z;
- if (t<tmin) tmin=t;
- if (t>tmax) tmax=t;
- }
- // change voxel coords from newimage to nifti convention for output
- ColumnVector v(4);
- v << xmin << ymin << zmin << 1.0;
- v = masknz.niftivox2newimagevox_mat().i() * v;
- xmin = MISCMATHS::round(v(1));
- ymin = MISCMATHS::round(v(2));
- zmin = MISCMATHS::round(v(3));
- v << xmax << ymax << zmax << 1.0;
- v = masknz.niftivox2newimagevox_mat().i() * v;
- xmax = MISCMATHS::round(v(1));
- ymax = MISCMATHS::round(v(2));
- zmax = MISCMATHS::round(v(3));
- if (xmin>xmax) { int tmp=xmax; xmax=xmin; xmin=tmp; }
- if (ymin>ymax) { int tmp=ymax; ymax=ymin; ymin=tmp; }
- if (zmin>zmax) { int tmp=zmax; zmax=zmin; zmin=tmp; }
- // now output nifti coords
- cout << xmin << " " << 1+xmax-xmin << " " << ymin << " " << 1+ymax-ymin << " " << zmin << " " << 1+zmax-zmin << " " << tmin << " " << 1+tmax-tmin << " ";
- } else if (sarg=="-e") {
- if (mask.nvoxels()<1) {
- generate_masks(mask,masknz,vol,lthr,uthr);
- vol*=mask;
- }
- ColumnVector hist;
- int nbins=1000;
- double entropy=0;
- hist = vol.histogram(nbins,mask);
- double ntot = hist.Sum();
- for (int j=1; j<=nbins; j++) {
- if (hist(j)>0) {
- entropy -= (hist(j)/ntot) * log(hist(j)/ntot);
- }
- }
- entropy /= log((double) nbins);
- cout << entropy << " ";
- } else if (sarg=="-E") {
- ColumnVector hist;
- int nbins=1000;
- double entropy=0;
- if (mask.nvoxels()<1) {
- generate_masks(mask,masknz,vol,lthr,uthr);
- vol*=mask;
- }
- hist = vol.histogram(nbins,masknz);
- double ntot = hist.Sum();
- for (int j=1; j<=nbins; j++) {
- if (hist(j)>0) {
- entropy -= (hist(j)/ntot) * log(hist(j)/ntot);
- }
- }
- entropy /= log((double) nbins);
- cout << entropy << " ";
- } else if (sarg=="-k") {
- narg++;
- volume<float> mask2;
- if (narg>=argc) {
- cerr << "Must specify an argument to -k" << endl;
- exit(2);
- }
- read_volume4D(mask,argv[narg]);
- if ( !samesize(mask,vol,3) || mask.tsize() > vol.tsize() ) {
- cerr << "Mask and image must be the same size" << endl;
- exit(3);
- }
- if (timeseriesMode && mask.tsize() != 1 ) { mask2=mask[timepoint]; }
- if ( mask.tsize() != vol.tsize() && mask.tsize() != 1)
- while (mask.tsize() < vol.tsize() ) { // copy the last max volume until the correct size is reached
- mask.addvolume(mask[mask.tsize()-1]);
- }
- mask.binarise(0.5);
- if (mask.tsize()!=1) {
- generateNonZeroMask(mask,masknz,vol);
- vol*=mask;
- }
- else {
- generateNonZeroMask(mask[0],masknz,vol);
- vol*=mask[0];
- }
- } else if (sarg=="-d") {
- narg++;
- if (narg>=argc) {
- cerr << "Must specify an argument to -d" << endl;
- exit(2);
- }
- volume4D<float> image2,image3;
- read_volume4D(image2,argv[narg]);
- if (!samesize(image2[0],vol[0])) {
- cerr << "Image used with -d must be the same size as the base image" << endl;
- exit(3);
- }
- if (timeseriesMode) { image3=image2[timepoint]; }
- if ( image2.tsize() > vol.tsize() ) {
- cerr << "Image must be the same size as the base image" << endl;
- exit(3);
- }
- if ( image2.tsize() != vol.tsize() && image2.tsize() != 1) {
- // copy the last max volume until the correct size is reached
- while (image2.tsize() < vol.tsize() ) {
- image2.addvolume(image2[image2.maxt()]);
- }
- }
- // now substract the new image from the base image
- vol -= image2;
- } else if (sarg=="-l") {
- narg++;
- if (narg<argc) {
- lthr = atof(argv[narg]);
- } else {
- cerr << "Must specify an argument to -l" << endl;
- exit(2);
- }
- generate_masks(mask,masknz,vol,lthr,uthr);
- vol*=mask;
- } else if (sarg=="-u") {
- narg++;
- if (narg<argc) {
- uthr = atof(argv[narg]);
- } else {
- cerr << "Must specify an argument to -u" << endl;
- exit(2);
- }
- generate_masks(mask,masknz,vol,lthr,uthr);
- vol*=mask;
- } else if (sarg=="-a") {
- vol = abs(vol);
- } else if (sarg=="-v") {
- if (mask.nvoxels()>0) {
- long int nvox = mask.sum();
- if (mask.tsize() == 1) nvox = nvox * vol.tsize();
- cout << (long int) nvox << " "
- << nvox * vol.xdim() * vol.ydim() * vol.zdim() << " ";
- } else {
- cout << (long int) vol.nvoxels() * vol.tsize() << " "
- << vol.nvoxels() * vol.tsize() * vol.xdim() * vol.ydim() * vol.zdim() << " ";
- }
- } else if (sarg=="-V") {
- if (masknz.nvoxels()>0) {
- cout << (long int) masknz.sum() << " "
- << masknz.sum() * vol.xdim() * vol.ydim() * vol.zdim() << " ";
- } else {
- long int nzvox;
- nzvox = nonzerocount(vol);
- cout << nzvox << " "
- << nzvox * vol.xdim() * vol.ydim() * vol.zdim() << " ";
- }
- } else if (sarg=="-D") {
- // hidden debug option!
- cout << vol.sum() << " ";
- } else if (sarg=="-s") {
- if (mask.nvoxels()>0) cout << vol.stddev(mask) << " ";
- else cout << vol.stddev() << " ";
- } else if (sarg=="-S") {
- if (masknz.nvoxels()>0) {
- cout << vol.stddev(masknz) << " ";
- } else {
- cout << nonzerostddev(vol) << " ";
- }
- } else if (sarg=="-r") {
- vector<float> limits(vol.robustlimits(mask));
- cout << limits[0] << " " << limits[1] << " ";
- } else if (sarg=="-R") {
- cout << vol.min(mask) << " " << vol.max(mask) << " ";
- } else if (sarg=="-c") {
- ColumnVector cog(4);
- // convert from fsl mm to voxel to nifti sform coord
- cog.SubMatrix(1,3,1,1) = vol.cog();
- cog(4) = 1.0;
- cog = vol.newimagevox2mm_mat() * cog;
- cout << cog(1) << " " << cog(2) << " " << cog(3) << " " ;
- } else if (sarg=="-C") {
- ColumnVector cog(4);
- // convert from fsl mm to fsl voxel coord to nifti voxel coord
- cog.SubMatrix(1,3,1,1) = vol.cog();
- cog(4) = 1.0;
- cog = vol.niftivox2newimagevox_mat().i() * cog;
- cout << cog(1) << " " << cog(2) << " " << cog(3) << " " ;
- } else if (sarg=="-p") {
- float n;
- narg++;
- if (narg<argc) {
- n = atof(argv[narg]);
- } else {
- cerr << "Must specify an argument to -p" << endl;
- exit(2);
- }
- if ( (n<0) || (n>100) ) {
- cerr << "Percentile must be between 0 and 100" << endl;
- exit(1);
- }
- if (mask.nvoxels()>0) cout << vol.percentile((float) n/100.0, mask) << " ";
- else cout << vol.percentile((float) n/100.0) << " ";
- } else if (sarg=="-P") {
- float n;
- narg++;
- if (narg<argc) {
- n = atof(argv[narg]);
- } else {
- cerr << "Must specify an argument to -P" << endl;
- exit(2);
- }
- if ( (n<0) || (n>100) ) {
- cerr << "Percentile must be between 0 and 100" << endl;
- exit(1);
- }
- if (mask.nvoxels()<1) {
- generate_masks(mask,masknz,vol,lthr,uthr);
- vol*=mask;
- }
- cout << vol.percentile((float) n/100.0,masknz) << " ";
- } else if (sarg=="-h") {
- float n;
- narg++;
- if (narg<argc) {
- n = atof(argv[narg]);
- } else {
- cerr << "Must specify the number of bins" << endl;
- exit(2);
- }
- int nbins = (int) n;
- if (nbins<1) {
- cerr << "Must specify at least 1 bin" << endl;
- exit(1);
- }
- if (mask.nvoxels()>0) {
- cout << vol.histogram(nbins,vol.min(),vol.max(),mask) << " ";
+ for (int timepoint = 0; timepoint < nTimepoints; timepoint++)
+ {
+ if (jsonOutput)
+ {
+ if (timepoint == 0){
+ json.append("{\"volumes\": [");
} else {
- cout << vol.histogram(nbins,vol.min(),vol.max()) << " ";
+ json.append(", ");
}
- } else if (sarg=="-H") {
- float n;
- narg++;
- if (narg<argc) {
- n = atof(argv[narg]);
- } else {
- cerr << "Must specify the number of bins" << endl;
- exit(2);
+ }
+ for (int index = 1; index <= nIndices; index++)
+ {
+ // if json mode is enabled
+ if (jsonOutput)
+ {
+ if (index == 1){
+ json.append("{\"indices\": [{");
+ } else {
+ json.append(", {");
+ }
}
- int nbins = (int) n;
- if (nbins<1) {
- cerr << "Must specify at least 1 bin" << endl;
- exit(1);
+ if (timeseriesMode)
+ vol = inputMaster[timepoint];
+ volume<float> mask, masknz;
+ float lthr(-numeric_limits<float>::max());
+ float uthr(numeric_limits<float>::max());
+ if (indexMask.nvoxels())
+ {
+ if (indexMask.dimensionality() > 3)
+ {
+ cerr << "Index mask must be 3D" << endl;
+ return 3;
+ }
+ copyconvert(indexMask, mask);
+ mask.binarise(index - 1, index + 1, exclusive);
+ vol = indexMaster * mask;
+ if (mask.max() < 1)
+ {
+ cout << "missing label: " << index << " in mask" << endl;
+ continue;
+ }
+ generateNonZeroMask(mask, masknz, vol);
}
- float min=0;
- narg++;
- if (narg<argc) {
- min = atof(argv[narg]);
- } else {
- cerr << "Must specify the histogram minimum intensity" << endl;
- exit(2);
+ int narg(2);
+ bool last = false;
+ while (narg < argc)
+ {
+
+ last = !(narg < argc - 1);
+ string sarg(argv[narg]);
+ if (sarg == "-n")
+ {
+ for (int t = 0; t < vol.tsize(); t++)
+ for (int z = 0; z < vol.zsize(); z++)
+ for (int y = 0; y < vol.ysize(); y++)
+ for (int x = 0; x < vol.xsize(); x++)
+ if (!isfinite((double)vol(x, y, z, t)))
+ vol(x, y, z, t) = 0;
+ }
+ else if (sarg == "-m")
+ {
+ if (mask.nvoxels() > 0){
+ string val = to_string(vol.mean(mask));
+ print_value(sarg, val, json, jsonOutput);
+ } else {
+ string val = to_string(vol.mean());
+ print_value(sarg, val, json, jsonOutput);
+ }
+ }
+ else if (sarg == "-M")
+ {
+ if (masknz.nvoxels() > 0)
+ print_value(sarg, to_string(vol.mean(masknz)), json, jsonOutput);
+ else
+ {
+ double nzmean = 0;
+ nzmean = nonzeromean(vol);
+ string val = to_string(nzmean);
+ print_value(sarg, val, json, jsonOutput);
+ }
+ }
+ else if (sarg == "-X")
+ {
+ ColumnVector coord(4);
+ vector<int64_t> iCoords;
+ vol.min(mask, iCoords);
+ coord << (Real)iCoords[7] << (Real)iCoords[8] << (Real)iCoords[9] << 1.0;
+ coord = vol.niftivox2newimagevox_mat().i() * coord;
+ //cout << MISCMATHS::round(coord(1)) << " " << MISCMATHS::round(coord(2)) << " " << MISCMATHS::round(coord(3)) << " ";
+ auto coord1 = MISCMATHS::round(coord(1));
+ auto coord2 = MISCMATHS::round(coord(2));
+ auto coord3 = MISCMATHS::round(coord(3));
+ string coordString = "";
+ coordString.append(to_string(coord1));
+ coordString.append(" ");
+ coordString.append(to_string(coord2));
+ coordString.append(" ");
+ coordString.append(to_string(coord3));
+ print_value(
+ sarg,
+ coordString,
+ json,
+ jsonOutput);
+
+ }
+ else if (sarg == "-x")
+ {
+ ColumnVector coord(4);
+ vector<int64_t> iCoords;
+ vol.min(mask, iCoords);
+ coord << (Real)iCoords[0] << (Real)iCoords[1] << (Real)iCoords[2] << 1.0;
+ coord = vol.niftivox2newimagevox_mat().i() * coord;
+ //cout << MISCMATHS::round(coord(1)) << " " << MISCMATHS::round(coord(2)) << " " << MISCMATHS::round(coord(3)) << " ";
+ auto coord1 = MISCMATHS::round(coord(1));
+ auto coord2 = MISCMATHS::round(coord(2));
+ auto coord3 = MISCMATHS::round(coord(3));
+ string coordString = "";
+ coordString.append(to_string(coord1));
+ coordString.append(" ");
+ coordString.append(to_string(coord2));
+ coordString.append(" ");
+ coordString.append(to_string(coord3));
+ print_value(
+ sarg,
+ coordString,
+ json,
+ jsonOutput);
+ }
+ else if (sarg == "-w")
+ {
+ if (masknz.nvoxels() < 1)
+ { // Need to generate non-zeromask
+ generate_masks(mask, masknz, vol, lthr, uthr);
+ vol *= mask;
+ }
+ int xmin = masknz.xsize() - 1, xmax(0), ymin = masknz.ysize() - 1, ymax(0), zmin = masknz.zsize() - 1, zmax(0), tmin = masknz.tsize() - 1, tmax(0);
+
+ for (int t = 0; t < masknz.tsize(); t++)
+ for (int z = 0; z < masknz.zsize(); z++)
+ for (int y = 0; y < masknz.ysize(); y++)
+ for (int x = 0; x < masknz.xsize(); x++)
+ if (masknz(x, y, z, t) > 0.5)
+ {
+ if (x < xmin)
+ xmin = x;
+ if (x > xmax)
+ xmax = x;
+ if (y < ymin)
+ ymin = y;
+ if (y > ymax)
+ ymax = y;
+ if (z < zmin)
+ zmin = z;
+ if (z > zmax)
+ zmax = z;
+ if (t < tmin)
+ tmin = t;
+ if (t > tmax)
+ tmax = t;
+ }
+ // change voxel coords from newimage to nifti convention for output
+ ColumnVector v(4);
+ v << xmin << ymin << zmin << 1.0;
+ v = masknz.niftivox2newimagevox_mat().i() * v;
+ xmin = MISCMATHS::round(v(1));
+ ymin = MISCMATHS::round(v(2));
+ zmin = MISCMATHS::round(v(3));
+ v << xmax << ymax << zmax << 1.0;
+ v = masknz.niftivox2newimagevox_mat().i() * v;
+ xmax = MISCMATHS::round(v(1));
+ ymax = MISCMATHS::round(v(2));
+ zmax = MISCMATHS::round(v(3));
+ if (xmin > xmax)
+ {
+ int tmp = xmax;
+ xmax = xmin;
+ xmin = tmp;
+ }
+ if (ymin > ymax)
+ {
+ int tmp = ymax;
+ ymax = ymin;
+ ymin = tmp;
+ }
+ if (zmin > zmax)
+ {
+ int tmp = zmax;
+ zmax = zmin;
+ zmin = tmp;
+ }
+ // now output nifti coords
+ string coords = "";
+ coords.append(to_string(xmin));
+ coords.append(" ");
+ coords.append(to_string(1 + xmax - xmin));
+ coords.append(" ");
+ coords.append(to_string(ymin));
+ coords.append(" ");
+ coords.append(to_string(1 + ymax - ymin));
+ coords.append(" ");
+ coords.append(to_string(zmin));
+ coords.append(" ");
+ coords.append(to_string(1 + zmax - zmin));
+ coords.append(" ");
+ coords.append(to_string(tmin));
+ coords.append(" ");
+ coords.append(to_string(1 + tmax - tmin));
+ print_value(
+ sarg,
+ coords,
+ json,
+ jsonOutput);
+ }
+ else if (sarg == "-e")
+ {
+ if (mask.nvoxels() < 1)
+ {
+ generate_masks(mask, masknz, vol, lthr, uthr);
+ vol *= mask;
+ }
+ ColumnVector hist;
+ int nbins = 1000;
+ double entropy = 0;
+ hist = vol.histogram(nbins, mask);
+ double ntot = hist.Sum();
+ for (int j = 1; j <= nbins; j++)
+ {
+ if (hist(j) > 0)
+ {
+ entropy -= (hist(j) / ntot) * log(hist(j) / ntot);
+ }
+ }
+ entropy /= log((double)nbins);
+ print_value(
+ sarg,
+ to_string(entropy),
+ json,
+ jsonOutput
+ );
+ }
+ else if (sarg == "-E")
+ {
+ ColumnVector hist;
+ int nbins = 1000;
+ double entropy = 0;
+ if (mask.nvoxels() < 1)
+ {
+ generate_masks(mask, masknz, vol, lthr, uthr);
+ vol *= mask;
+ }
+ hist = vol.histogram(nbins, masknz);
+ double ntot = hist.Sum();
+ for (int j = 1; j <= nbins; j++)
+ {
+ if (hist(j) > 0)
+ {
+ entropy -= (hist(j) / ntot) * log(hist(j) / ntot);
+ }
+ }
+ entropy /= log((double)nbins);
+ print_value(
+ sarg,
+ to_string(entropy),
+ json,
+ jsonOutput
+ );
+ }
+ else if (sarg == "-k")
+ {
+ narg++;
+ volume<float> mask2;
+ if (narg >= argc)
+ {
+ cerr << "Must specify an argument to -k" << endl;
+ exit(2);
+ }
+ read_volume4D(mask, argv[narg]);
+ if (!samesize(mask, vol, 3) || mask.tsize() > vol.tsize())
+ {
+ cerr << "Mask and image must be the same size" << endl;
+ exit(3);
+ }
+ if (timeseriesMode && mask.tsize() != 1)
+ {
+ mask2 = mask[timepoint];
+ }
+ if (mask.tsize() != vol.tsize() && mask.tsize() != 1)
+ while (mask.tsize() < vol.tsize())
+ { // copy the last max volume until the correct size is reached
+ mask.addvolume(mask[mask.tsize() - 1]);
+ }
+ mask.binarise(0.5);
+ if (mask.tsize() != 1)
+ {
+ generateNonZeroMask(mask, masknz, vol);
+ vol *= mask;
+ }
+ else
+ {
+ generateNonZeroMask(mask[0], masknz, vol);
+ vol *= mask[0];
+ }
+ }
+ else if (sarg == "-d")
+ {
+ narg++;
+ if (narg >= argc)
+ {
+ cerr << "Must specify an argument to -d" << endl;
+ exit(2);
+ }
+ volume4D<float> image2, image3;
+ read_volume4D(image2, argv[narg]);
+ if (!samesize(image2[0], vol[0]))
+ {
+ cerr << "Image used with -d must be the same size as the base image" << endl;
+ exit(3);
+ }
+ if (timeseriesMode)
+ {
+ image3 = image2[timepoint];
+ }
+ if (image2.tsize() > vol.tsize())
+ {
+ cerr << "Image must be the same size as the base image" << endl;
+ exit(3);
+ }
+ if (image2.tsize() != vol.tsize() && image2.tsize() != 1)
+ {
+ // copy the last max volume until the correct size is reached
+ while (image2.tsize() < vol.tsize())
+ {
+ image2.addvolume(image2[image2.maxt()]);
+ }
+ }
+ // now substract the new image from the base image
+ vol -= image2;
+ }
+ else if (sarg == "-l")
+ {
+ narg++;
+ if (narg < argc)
+ {
+ lthr = atof(argv[narg]);
+ }
+ else
+ {
+ cerr << "Must specify an argument to -l" << endl;
+ exit(2);
+ }
+ generate_masks(mask, masknz, vol, lthr, uthr);
+ vol *= mask;
+ }
+ else if (sarg == "-u")
+ {
+ narg++;
+ if (narg < argc)
+ {
+ uthr = atof(argv[narg]);
+ }
+ else
+ {
+ cerr << "Must specify an argument to -u" << endl;
+ exit(2);
+ }
+ generate_masks(mask, masknz, vol, lthr, uthr);
+ vol *= mask;
+ }
+ else if (sarg == "-a")
+ {
+ vol = abs(vol);
+ }
+ else if (sarg == "-v")
+ {
+ if (mask.nvoxels() > 0)
+ {
+ long int nvox = mask.sum();
+ if (mask.tsize() == 1)
+ nvox = nvox * vol.tsize();
+
+ string nvox_str = "";
+ nvox_str.append(to_string(nvox));
+ nvox_str.append(" ");
+ nvox_str.append(to_string(nvox * vol.xdim() * vol.ydim() * vol.zdim()));
+ print_value(
+ sarg,
+ nvox_str,
+ json,
+ jsonOutput
+ );
+ }
+ else
+ {
+ string nvox_str = "";
+ nvox_str.append(to_string((long int)vol.nvoxels() * vol.tsize()));
+ nvox_str.append(" ");
+ nvox_str.append(to_string(vol.nvoxels() * vol.tsize() * vol.xdim() * vol.ydim() * vol.zdim()));
+ print_value(
+ sarg,
+ nvox_str,
+ json,
+ jsonOutput);
+ }
+ }
+ else if (sarg == "-V")
+ {
+ if (masknz.nvoxels() > 0)
+ {
+ string nvox_str = "";
+ nvox_str.append(to_string((long int)masknz.sum()));
+ nvox_str.append(" ");
+ nvox_str.append(to_string(masknz.sum() * vol.xdim() * vol.ydim() * vol.zdim()));
+ print_value(
+ sarg,
+ nvox_str,
+ json,
+ jsonOutput);
+ }
+ else
+ {
+ long int nzvox;
+ nzvox = nonzerocount(vol);
+ string nzvox_str = "";
+ nzvox_str.append(to_string(nzvox));
+ nzvox_str.append(" ");
+ nzvox_str.append(to_string(nzvox * vol.xdim() * vol.ydim() * vol.zdim()));
+ print_value(
+ sarg,
+ nzvox_str,
+ json,
+ jsonOutput);
+ }
+ }
+ else if (sarg == "-D")
+ {
+ // hidden debug option!
+ print_value(
+ sarg,
+ to_string(vol.sum()),
+ json,
+ jsonOutput
+ );
+ }
+ else if (sarg == "-s")
+ {
+ if (mask.nvoxels() > 0)
+ print_value(
+ sarg,
+ to_string(vol.stddev(mask)),
+ json,
+ jsonOutput
+ );
+ else
+ print_value(
+ sarg,
+ to_string(vol.stddev()),
+ json,
+ jsonOutput
+ );
+ }
+ else if (sarg == "-S")
+ {
+ if (masknz.nvoxels() > 0)
+ {
+ print_value(
+ sarg,
+ to_string(vol.stddev(masknz)),
+ json,
+ jsonOutput
+ );
+ }
+ else
+ {
+ print_value(
+ sarg,
+ to_string(nonzerostddev(vol)),
+ json,
+ jsonOutput
+ );
+ }
+ }
+ else if (sarg == "-r")
+ {
+ vector<float> limits(vol.robustlimits(mask));
+ print_value(
+ sarg,
+ to_string(limits[0]).append(" ").append(to_string(limits[1])),
+ json,
+ jsonOutput
+ );
+ }
+ else if (sarg == "-R")
+ {
+ print_value(
+ sarg,
+ to_string(vol.min(mask)).append(" ").append(to_string(vol.max(mask))),
+ json,
+ jsonOutput
+ );
+ }
+ else if (sarg == "-c")
+ {
+ ColumnVector cog(4);
+ cog.SubMatrix(1, 3, 1, 1) = vol.cog();
+ cog(4) = 1.0;
+ cog = vol.newimagevox2mm_mat() * cog;
+ print_value(
+ sarg,
+ to_string(cog(1)).append(" ").append(to_string(cog(2))).append(" ").append(to_string(cog(3))),
+ json,
+ jsonOutput
+ );
+ }
+ else if (sarg == "-C")
+ {
+ ColumnVector cog(4);
+ // convert from fsl mm to fsl voxel coord to nifti voxel coord
+ cog.SubMatrix(1, 3, 1, 1) = vol.cog();
+ cog(4) = 1.0;
+ cog = vol.niftivox2newimagevox_mat().i() * cog;
+ print_value(
+ sarg,
+ to_string(cog(1)).append(" ").append(to_string(cog(2))).append(" ").append(to_string(cog(3))),
+ json,
+ jsonOutput
+ );
+ }
+ else if (sarg == "-p")
+ {
+ float n;
+ narg++;
+ if (narg < argc)
+ {
+ n = atof(argv[narg]);
+ }
+ else
+ {
+ cerr << "Must specify an argument to -p" << endl;
+ exit(2);
+ }
+ if ((n < 0) || (n > 100))
+ {
+ cerr << "Percentile must be between 0 and 100" << endl;
+ exit(1);
+ }
+ if (mask.nvoxels() > 0)
+ print_value(
+ sarg,
+ to_string(vol.percentile((float)n / 100.0, mask)),
+ json,
+ jsonOutput
+ );
+ else
+ print_value(
+ sarg,
+ to_string(vol.percentile((float)n / 100.0)),
+ json,
+ jsonOutput
+ );
+ }
+ else if (sarg == "-P")
+ {
+ float n;
+ narg++;
+ if (narg < argc)
+ {
+ n = atof(argv[narg]);
+ }
+ else
+ {
+ cerr << "Must specify an argument to -P" << endl;
+ exit(2);
+ }
+ if ((n < 0) || (n > 100))
+ {
+ cerr << "Percentile must be between 0 and 100" << endl;
+ exit(1);
+ }
+ if (mask.nvoxels() < 1)
+ {
+ generate_masks(mask, masknz, vol, lthr, uthr);
+ vol *= mask;
+ }
+ print_value(
+ sarg,
+ to_string(vol.percentile((float)n / 100.0, masknz)),
+ json,
+ jsonOutput
+ );
+ }
+ else if (sarg == "-h")
+ {
+ float n;
+ narg++;
+ if (narg < argc)
+ {
+ n = atof(argv[narg]);
+ }
+ else
+ {
+ cerr << "Must specify the number of bins" << endl;
+ exit(2);
+ }
+ int nbins = (int)n;
+ if (nbins < 1)
+ {
+ cerr << "Must specify at least 1 bin" << endl;
+ exit(1);
+ }
+ if (mask.nvoxels() > 0)
+ {
+ auto hist = vol.histogram(nbins, vol.min(), vol.max(), mask);
+ // loop over histogram values and and build a string to print
+ string hist_str = "";
+ for (int i = 0; i < hist.Nrows(); i++)
+ {
+ hist_str += to_string(hist(i+1));
+ if (i < hist.Nrows() - 1)
+ {
+ hist_str += " ";
+ }
+ }
+ print_value(
+ sarg,
+ hist_str,
+ json,
+ jsonOutput
+ );
+ }
+ else
+ {
+ auto hist = vol.histogram(nbins, vol.min(), vol.max());
+ string hist_str = "";
+ for (int i = 0; i < hist.Nrows(); i++)
+ {
+ hist_str += to_string(hist(i+1));
+ if (i < hist.Nrows() - 1)
+ {
+ hist_str += " ";
+ }
+ }
+ print_value(
+ sarg,
+ hist_str,
+ json,
+ jsonOutput
+ );
+ }
+ }
+ else if (sarg == "-H")
+ {
+ float n;
+ narg++;
+ if (narg < argc)
+ {
+ n = atof(argv[narg]);
+ }
+ else
+ {
+ cerr << "Must specify the number of bins" << endl;
+ exit(2);
+ }
+ int nbins = (int)n;
+ if (nbins < 1)
+ {
+ cerr << "Must specify at least 1 bin" << endl;
+ exit(1);
+ }
+ float min = 0;
+ narg++;
+ if (narg < argc)
+ {
+ min = atof(argv[narg]);
+ }
+ else
+ {
+ cerr << "Must specify the histogram minimum intensity" << endl;
+ exit(2);
+ }
+ float max = 0;
+ narg++;
+ if (narg < argc)
+ {
+ max = atof(argv[narg]);
+ }
+ else
+ {
+ cerr << "Must specify the histogram maximum intensity" << endl;
+ exit(2);
+ }
+ if (mask.nvoxels() > 0)
+ {
+ auto hist = vol.histogram(nbins, min, max, mask);
+ string hist_str = "";
+ for (int i = 0; i < hist.Nrows(); i++)
+ {
+ hist_str += to_string(hist(i+1));
+ if (i < hist.Nrows() - 1)
+ {
+ hist_str += " ";
+ }
+ }
+ print_value(
+ sarg,
+ hist_str,
+ json,
+ jsonOutput
+ );
+ }
+ else
+ {
+ auto hist = vol.histogram(nbins, min, max);
+ string hist_str = "";
+ for (int i = 0; i < hist.Nrows(); i++)
+ {
+ hist_str += to_string(hist(i+1));
+ if (i < hist.Nrows() - 1)
+ {
+ hist_str += " ";
+ }
+ }
+ print_value(
+ sarg,
+ hist_str,
+ json,
+ jsonOutput
+ );
+ }
+ }
+ else
+ {
+ // only print cerr if not the --json option
+ if (sarg != "--json")
+ {
+ cerr << "Unrecognised option: " << sarg << endl;
+ exit(3);
+ } else if (sarg == "--json" && last) {
+ // remove the trailing comma and space
+ json.pop_back();
+ json.pop_back();
+ }
+ }
+
+ narg++;
}
- float max=0;
- narg++;
- if (narg<argc) {
- max = atof(argv[narg]);
- } else {
- cerr << "Must specify the histogram maximum intensity" << endl;
- exit(2);
+ // if json mode is enabled
+ if (jsonOutput)
+ {
+ // close the json object
+ json.append("}");
}
- if (mask.nvoxels()>0) {
- cout << vol.histogram(nbins,min,max,mask) << " ";
- } else {
- cout << vol.histogram(nbins,min,max) << " ";
+ if (index < nIndices && !jsonOutput)
+ {
+ cout << endl;
}
- } else {
- cerr << "Unrecognised option: " << sarg << endl;
- exit(3);
+ // cout << endl;
+ } // end nIndices
+ // if json mode is enabled
+ if (jsonOutput)
+ {
+ // close the json object
+ json.append("]}");
}
-
- narg++;
- }
- cout << endl;
- }
-
+ if (timepoint < nTimepoints - 1 && !jsonOutput)
+ {
+ cout << endl;
+ }
+ // cout << endl;
+ } // end timepoints
+ if (jsonOutput)
+ {
+ // close the json object
+ json.append("]}");
+ cout << json;
+ } else {
+ cout << endl;
}
return 0;
}
-
-
-int main(int argc,char *argv[])
+int main(int argc, char *argv[])
{
Tracer tr("main");
string progname(argv[0]);
bool timeseriesMode(false);
string indexMask("");
- while ( argc > 2 && ( string(argv[1])=="-t" || string(argv[1]) =="-K" ) ) {
- if ( string(argv[1])=="-t" )
- timeseriesMode=true;
- if ( string(argv[1])=="-K" ) {
- indexMask=string(argv[2]);
+ while (argc > 2 && (string(argv[1]) == "-t" || string(argv[1]) == "-K"))
+ {
+ if (string(argv[1]) == "-t")
+ timeseriesMode = true;
+ if (string(argv[1]) == "-K")
+ {
+ indexMask = string(argv[2]);
argv++;
argc--;
}
argv++;
argc--;
}
- if (argc < 3 )
+ if (argc < 3)
return print_usage(progname);
- try {
- return fmrib_main_float(argc,argv,timeseriesMode,indexMask);
- } catch(std::exception &e) {
+ try
+ {
+ return fmrib_main_float(argc, argv, timeseriesMode, indexMask);
+ }
+ catch (std::exception &e)
+ {
cerr << e.what() << endl;
- } catch (...) {
+ }
+ catch (...)
+ {
// do nothing - just exit without garbage message
}
return -1;
-
}
--
GitLab