Skip to content
Snippets Groups Projects
Commit 84f782e8 authored by Taylor Hanayik's avatar Taylor Hanayik
Browse files

add --json option and associated logic

parent 722037d7
No related branches found
No related tags found
1 merge request!22add --json option and associated logic
Pipeline #18996 skipped
Stage: fsl-ci-build
This commit is part of merge request !22. Comments created here will be created in the context of that merge request.
Hide whitespace changes
Inline Side-by-side
.gitignore 0 → 100644
+ 1
− 0
View file @ 84f782e8
.vscode
\ No newline at end of file
fslstats.cc
+ 866
− 397
View file @ 84f782e8
...@@ -15,16 +15,18 @@ ...@@ -15,16 +15,18 @@
#include "newimage/newimageall.h" #include "newimage/newimageall.h"
#include "newimage/costfns.h" #include "newimage/costfns.h"
using namespace std; using namespace std;
using namespace NEWMAT; using namespace NEWMAT;
using namespace NEWIMAGE; using namespace NEWIMAGE;
int print_usage(const string& progname) { int print_usage(const string &progname)
cout << "Usage: fslstats [preoptions] <input> [options]" << endl << endl; {
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 -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 << "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 << "-l <lthresh> : set lower threshold" << endl;
cout << "-u <uthresh> : set upper threshold" << endl; cout << "-u <uthresh> : set upper threshold" << endl;
cout << "-r : output <robust min intensity> <robust max intensity>" << endl; cout << "-r : output <robust min intensity> <robust max intensity>" << endl;
...@@ -44,12 +46,16 @@ int print_usage(const string& progname) { ...@@ -44,12 +46,16 @@ int print_usage(const string& progname) {
cout << "-C : output centre-of-gravity (cog) in voxel coordinates" << endl; 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 (n between 0 and 100)" << endl;
cout << "-P <n> : output nth percentile (for nonzero voxels)" << 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 << "-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 << "-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 << "-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> : 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; cout << "Note - thresholds are not inclusive ie lthresh<allowed<uthresh" << endl;
return 1; return 1;
} }
...@@ -57,461 +63,924 @@ int print_usage(const string& progname) { ...@@ -57,461 +63,924 @@ int print_usage(const string& progname) {
// Some specialised nonzero functions just for speedup // Some specialised nonzero functions just for speedup
// (it avoids generating masks when not absolutely necessary) // (it avoids generating masks when not absolutely necessary)
long nonzerocount(const volume4D<float>& vol) long nonzerocount(const volume4D<float> &vol)
{ {
long totn=0; long totn = 0;
for(volume4D<float>::fast_const_iterator it=vol.fbegin(), end=vol.fend(); it<end; it++ ) for (volume4D<float>::fast_const_iterator it = vol.fbegin(), end = vol.fend(); it < end; it++)
if ( (*it) != 0 ) if ((*it) != 0)
totn++; totn++;
return totn; return totn;
} }
double nonzeromean(const volume4D<float>& vol) double nonzeromean(const volume4D<float> &vol)
{ {
double total=0.0; double total = 0.0;
long int totn=0; long int totn = 0;
for(volume4D<float>::fast_const_iterator it=vol.fbegin(), end=vol.fend(); it<end; it++ ) for (volume4D<float>::fast_const_iterator it = vol.fbegin(), end = vol.fend(); it < end; it++)
if ( (*it) != 0 ) { if ((*it) != 0)
total+=(double)(*it); {
total += (double)(*it);
totn++; totn++;
} }
if (totn>0) if (totn > 0)
total/=totn; total /= totn;
return total; return total;
} }
double nonzerostddev(const volume4D<float>& vol) double nonzerostddev(const volume4D<float> &vol)
{ {
double totv=0.0, totvv=0.0; double totv = 0.0, totvv = 0.0;
long int totn=0; long int totn = 0;
for(volume4D<float>::fast_const_iterator it=vol.fbegin(), end=vol.fend(); it<end; it++ ) for (volume4D<float>::fast_const_iterator it = vol.fbegin(), end = vol.fend(); it < end; it++)
if ( (*it) != 0 ) { if ((*it) != 0)
totvv+=(double)((*it)*(*it)); {
totv+=(double)(*it); totvv += (double)((*it) * (*it));
totv += (double)(*it);
totn++; totn++;
} }
double var=0.0; double var = 0.0;
if (totn>1) { if (totn > 1)
double meanval = totv/totn; {
var = (totvv - totn*meanval*meanval)/(totn-1); double meanval = totv / totn;
var = (totvv - totn * meanval * meanval) / (totn - 1);
} }
return std::sqrt(var); return std::sqrt(var);
} }
template<class M> template <class M>
int generateNonZeroMask(const M &mask, volume4D<float> &masknz, const volume4D<float> &input) 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; 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 &lthr, 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); // append the key and value to the json string
return generateNonZeroMask(mask,masknz,input); 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::dec);
cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::fixed, ios::floatfield);
cout.setf(ios::left, ios::adjustfield); cout.setf(ios::left, ios::adjustfield);
cout.precision(6); cout.precision(6);
volume<float> vol, inputMaster; volume<float> vol, inputMaster;
volume<int> indexMask; volume<int> indexMask;
if ( timeseriesMode || indexMaskName != "" ) if (timeseriesMode || indexMaskName != "")
read_volume4D(inputMaster,argv[1]); read_volume4D(inputMaster, argv[1]);
else else
read_volume4D(vol,argv[1]); read_volume4D(vol, argv[1]);
volume<float> & indexMaster = (timeseriesMode ) ? vol : inputMaster; volume<float> &indexMaster = (timeseriesMode) ? vol : inputMaster;
if ( indexMaskName != "" ) if (indexMaskName != "")
read_volume4D(indexMask,indexMaskName); read_volume4D(indexMask, indexMaskName);
int nTimepoints((timeseriesMode) ? inputMaster.tsize() : 1), nIndices((indexMaskName != "") ? indexMask.max() : 1); 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 timepoint = 0; timepoint < nTimepoints; timepoint++)
for(int z=0;z<masknz.zsize();z++) {
for(int y=0;y<masknz.ysize();y++) if (jsonOutput)
for(int x=0;x<masknz.xsize();x++) {
if (masknz(x,y,z,t)>0.5) { if (timepoint == 0){
if (x<xmin) xmin=x; json.append("{\"volumes\": [");
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) << " ";
} else { } else {
cout << vol.histogram(nbins,vol.min(),vol.max()) << " "; json.append(", ");
} }
} else if (sarg=="-H") { }
float n; for (int index = 1; index <= nIndices; index++)
narg++; {
if (narg<argc) { // if json mode is enabled
n = atof(argv[narg]); if (jsonOutput)
} else { {
cerr << "Must specify the number of bins" << endl; if (index == 1){
exit(2); json.append("{\"indices\": [{");
} else {
json.append(", {");
}
} }
int nbins = (int) n; if (timeseriesMode)
if (nbins<1) { vol = inputMaster[timepoint];
cerr << "Must specify at least 1 bin" << endl; volume<float> mask, masknz;
exit(1); 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; int narg(2);
narg++; bool last = false;
if (narg<argc) { while (narg < argc)
min = atof(argv[narg]); {
} else {
cerr << "Must specify the histogram minimum intensity" << endl; last = !(narg < argc - 1);
exit(2); 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; // if json mode is enabled
narg++; if (jsonOutput)
if (narg<argc) { {
max = atof(argv[narg]); // close the json object
} else { json.append("}");
cerr << "Must specify the histogram maximum intensity" << endl;
exit(2);
} }
if (mask.nvoxels()>0) { if (index < nIndices && !jsonOutput)
cout << vol.histogram(nbins,min,max,mask) << " "; {
} else { cout << endl;
cout << vol.histogram(nbins,min,max) << " ";
} }
} else { // cout << endl;
cerr << "Unrecognised option: " << sarg << endl; } // end nIndices
exit(3); // if json mode is enabled
if (jsonOutput)
{
// close the json object
json.append("]}");
} }
if (timepoint < nTimepoints - 1 && !jsonOutput)
narg++; {
} cout << endl;
cout << endl; }
} // cout << endl;
} // end timepoints
if (jsonOutput)
{
// close the json object
json.append("]}");
cout << json;
} else {
cout << endl;
} }
return 0; return 0;
} }
int main(int argc, char *argv[])
int main(int argc,char *argv[])
{ {
Tracer tr("main"); Tracer tr("main");
string progname(argv[0]); string progname(argv[0]);
bool timeseriesMode(false); bool timeseriesMode(false);
string indexMask(""); string indexMask("");
while ( argc > 2 && ( string(argv[1])=="-t" || string(argv[1]) =="-K" ) ) { while (argc > 2 && (string(argv[1]) == "-t" || string(argv[1]) == "-K"))
if ( string(argv[1])=="-t" ) {
timeseriesMode=true; if (string(argv[1]) == "-t")
if ( string(argv[1])=="-K" ) { timeseriesMode = true;
indexMask=string(argv[2]); if (string(argv[1]) == "-K")
{
indexMask = string(argv[2]);
argv++; argv++;
argc--; argc--;
} }
argv++; argv++;
argc--; argc--;
} }
if (argc < 3 ) if (argc < 3)
return print_usage(progname); return print_usage(progname);
try { try
return fmrib_main_float(argc,argv,timeseriesMode,indexMask); {
} catch(std::exception &e) { return fmrib_main_float(argc, argv, timeseriesMode, indexMask);
}
catch (std::exception &e)
{
cerr << e.what() << endl; cerr << e.what() << endl;
} catch (...) { }
catch (...)
{
// do nothing - just exit without garbage message // do nothing - just exit without garbage message
} }
return -1; return -1;
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment