From 61d6cfff3f8a1905b091b1047b87ad8a02ca3e6e Mon Sep 17 00:00:00 2001
From: Matthew Webster <mwebster@fmrib.ox.ac.uk>
Date: Tue, 11 Jun 2019 15:25:37 +0000
Subject: [PATCH] read4D to avoid calling legacy read
---
fslchfiletype.cc | 2 +-
fslcomplex.cc | 37 ++++++++--------
fslstats.cc | 107 +++++++++++++++++++++++------------------------
3 files changed, 72 insertions(+), 74 deletions(-)
diff --git a/fslchfiletype.cc b/fslchfiletype.cc
index ed79241..1fb07f9 100644
--- a/fslchfiletype.cc
+++ b/fslchfiletype.cc
@@ -26,7 +26,7 @@ template <class T>
int fmrib_main(int argc, char *argv[])
{
volume<T> input;
- read_volume(input,string(argv[2]));
+ read_volume4D(input,string(argv[2]));
int format(-1);
const vector<pair<string,int> > formats=NEWIMAGE::fileFormats::allFormats();
for ( std::vector<pair<string,int> >::const_iterator it=formats.begin(); it != formats.end(); it++ )
diff --git a/fslcomplex.cc b/fslcomplex.cc
index 1731849..f3b07f5 100644
--- a/fslcomplex.cc
+++ b/fslcomplex.cc
@@ -13,19 +13,19 @@ using namespace NEWIMAGE;
void print_usage(const string& progname) {
cout << "Usage: " << progname << " <-outputtype> <input> <output> [startvol [endvol]]" << endl << endl;
- cout << " " << progname << " -realabs complexvol absvol"
+ cout << " " << progname << " -realabs complexvol absvol"
<< " [startvol [endvol]]" << endl;
- cout << " " << progname << " -realphase complexvol phasevol"
+ cout << " " << progname << " -realphase complexvol phasevol"
<< " [startvol [endvol]]" << endl;
- cout << " " << progname << " -realpolar complexvol absvol phasevol"
+ cout << " " << progname << " -realpolar complexvol absvol phasevol"
<< " [startvol [endvol]]" << endl;
cout << " " << progname << " -realcartesian complexvol realvol"
<< " imagvol [startvol [endvol]]" << endl;
- cout << " " << progname << " -complex realvol imagvol"
+ cout << " " << progname << " -complex realvol imagvol"
<< " complexvol [startvol [endvol]]" << endl;
- cout << " " << progname << " -complexpolar absvol phasevol"
+ cout << " " << progname << " -complexpolar absvol phasevol"
<< " complexvol [startvol [endvol]]" << endl;
- cout << " " << progname << " -complexsplit source dest"
+ cout << " " << progname << " -complexsplit source dest"
<< " [startvol [endvol]]" << endl;
cout << " " << progname << " -complexmerge source1 source2 dest" <<endl;
cout << " " << progname << " -copyonly source dest" << endl;
@@ -42,8 +42,8 @@ void fix_start_and_end(int& start, int& end, int mint, int maxt)
}
-void realpolar(const string& fin, const string& fabs, const string& fphase,
- int start, int end)
+void realpolar(const string& fin, const string& fabs, const string& fphase,
+ int start, int end)
{
volume<float> vreal, vimag;
read_complexvolume(vreal, vimag, fin);
@@ -68,7 +68,7 @@ void realpolar(const string& fin, const string& fabs, const string& fphase,
}
}
-void realcartesian(const string& fin, const string& freal, const string& fimag,
+void realcartesian(const string& fin, const string& freal, const string& fimag,
int start, int end)
{
volume<float> vreal, vimag;
@@ -95,12 +95,12 @@ void realcartesian(const string& fin, const string& freal, const string& fimag,
}
-void complexsave(const string& freal, const string& fimag,
+void complexsave(const string& freal, const string& fimag,
const string& fcomplex, int start, int end)
{
volume<float> vreal, vimag;
- read_volume(vreal, freal);
- read_volume(vimag, fimag);
+ read_volume4D(vreal, freal);
+ read_volume4D(vimag, fimag);
fix_start_and_end(start,end,Max(vreal.mint(),vimag.mint()),
Min(vreal.maxt(),vimag.maxt()));
@@ -117,12 +117,12 @@ void complexsave(const string& freal, const string& fimag,
}
-void complexpolar(const string& fabsvol, const string& fphasevol,
+void complexpolar(const string& fabsvol, const string& fphasevol,
const string& fcomplex, int start, int end)
{
volume<float> vabs, vphase;
- read_volume(vabs, fabsvol);
- read_volume(vphase, fphasevol);
+ read_volume4D(vabs, fabsvol);
+ read_volume4D(vphase, fphasevol);
fix_start_and_end(start,end,Max(vabs.mint(),vphase.mint()),
Min(vabs.maxt(),vphase.maxt()));
@@ -138,7 +138,7 @@ void complexpolar(const string& fabsvol, const string& fphasevol,
}
}
-void complexsplit(const string& fsource, const string& fdest,
+void complexsplit(const string& fsource, const string& fdest,
int start, int end)
{
volume4D<float> vreal, vimag;
@@ -204,9 +204,9 @@ int main(int argc,char *argv[])
Tracer tr("main");
string progname=argv[0];
- if (argc<4) {
+ if (argc<4) {
print_usage(progname);
- return -1;
+ return -1;
}
string arg = argv[1];
@@ -259,4 +259,3 @@ int main(int argc,char *argv[])
return 0;
}
-
diff --git a/fslstats.cc b/fslstats.cc
index 7196159..5fc26e3 100644
--- a/fslstats.cc
+++ b/fslstats.cc
@@ -15,7 +15,7 @@
using namespace NEWIMAGE;
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 -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;
@@ -42,7 +42,7 @@ int print_usage(const string& progname) {
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> : 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 << "Note - thresholds are not inclusive ie lthresh<allowed<uthresh" << endl;
return 1;
@@ -64,12 +64,12 @@ 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++ )
+ 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)
+ if (totn>0)
total/=totn;
return total;
}
@@ -78,7 +78,7 @@ 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++ )
+ 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);
@@ -95,31 +95,31 @@ double nonzerostddev(const volume4D<float>& vol)
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& lthr, const float& uthr)
{
mask = binarise(input,lthr,uthr,exclusive);
return generateNonZeroMask(mask,masknz,input);
}
-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)
{
- cout.setf(ios::dec);
- cout.setf(ios::fixed, ios::floatfield);
- cout.setf(ios::left, ios::adjustfield);
- cout.precision(6);
+ 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_volume(inputMaster,argv[1]);
- else
- read_volume(vol,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_volume(indexMask,indexMaskName);
+ 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++ ) {
@@ -127,7 +127,7 @@ int fmrib_main_float(int argc, char* argv[],const bool timeseriesMode, const str
vol=inputMaster[timepoint];
volume<float> mask, masknz;
float lthr(-numeric_limits<float>::max());
- float uthr(numeric_limits<float>::max());
+ float uthr(numeric_limits<float>::max());
if ( indexMask.nvoxels() ) {
if ( indexMask.dimensionality() > 3 ) {
cerr << "Index mask must be 3D" << endl;
@@ -165,27 +165,27 @@ int fmrib_main_float(int argc, char* argv[],const bool timeseriesMode, const str
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)) << " " <<
+ cout << MISCMATHS::round(coord(1)) << " " <<
MISCMATHS::round(coord(2)) << " " << MISCMATHS::round(coord(3)) << " ";
- } else if (sarg=="-x") {
+ } 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)) << " " <<
+ 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;
+ 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++)
+
+ 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;
@@ -215,8 +215,8 @@ int fmrib_main_float(int argc, char* argv[],const bool timeseriesMode, const str
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;
+ generate_masks(mask,masknz,vol,lthr,uthr);
+ vol*=mask;
}
ColumnVector hist;
int nbins=1000;
@@ -225,24 +225,24 @@ int fmrib_main_float(int argc, char* argv[],const bool timeseriesMode, const str
double ntot = hist.Sum();
for (int j=1; j<=nbins; j++) {
if (hist(j)>0) {
- entropy -= (hist(j)/ntot) * log(hist(j)/ntot);
+ entropy -= (hist(j)/ntot) * log(hist(j)/ntot);
}
}
entropy /= log((double) nbins);
cout << entropy << " ";
- } else if (sarg=="-E") {
+ } 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;
+ 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 -= (hist(j)/ntot) * log(hist(j)/ntot);
}
}
entropy /= log((double) nbins);
@@ -260,14 +260,14 @@ int fmrib_main_float(int argc, char* argv[],const bool timeseriesMode, const str
exit(3);
}
if (timeseriesMode && mask.tsize() != 1 ) { mask2=mask[timepoint]; }
- if ( mask.tsize() != vol.tsize() && mask.tsize() != 1)
+ 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) {
+ if (mask.tsize()!=1) {
generateNonZeroMask(mask,masknz,vol);
- vol*=mask;
+ vol*=mask;
}
else {
generateNonZeroMask(mask[0],masknz,vol);
@@ -296,9 +296,9 @@ int fmrib_main_float(int argc, char* argv[],const bool timeseriesMode, const str
image2.addvolume(image2[image2.maxt()]);
}
}
- // now substract the new image from the base image
+ // now substract the new image from the base image
vol -= image2;
- } else if (sarg=="-l") {
+ } else if (sarg=="-l") {
narg++;
if (narg<argc) {
lthr = atof(argv[narg]);
@@ -324,7 +324,7 @@ int fmrib_main_float(int argc, char* argv[],const bool timeseriesMode, const str
if (mask.nvoxels()>0) {
long int nvox = mask.sum();
if (mask.tsize() == 1) nvox = nvox * vol.tsize();
- cout << (long int) nvox << " "
+ cout << (long int) nvox << " "
<< nvox * vol.xdim() * vol.ydim() * vol.zdim() << " ";
} else {
cout << (long int) vol.nvoxels() * vol.tsize() << " "
@@ -332,12 +332,12 @@ int fmrib_main_float(int argc, char* argv[],const bool timeseriesMode, const str
}
} else if (sarg=="-V") {
if (masknz.nvoxels()>0) {
- cout << (long int) masknz.sum() << " "
+ cout << (long int) masknz.sum() << " "
<< masknz.sum() * vol.xdim() * vol.ydim() * vol.zdim() << " ";
} else {
long int nzvox;
nzvox = nonzerocount(vol);
- cout << nzvox << " "
+ cout << nzvox << " "
<< nzvox * vol.xdim() * vol.ydim() * vol.zdim() << " ";
}
} else if (sarg=="-D") {
@@ -362,7 +362,7 @@ int fmrib_main_float(int argc, char* argv[],const bool timeseriesMode, const str
// 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;
+ cog = vol.newimagevox2mm_mat() * cog;
cout << cog(1) << " " << cog(2) << " " << cog(3) << " " ;
} else if (sarg=="-C") {
ColumnVector cog(4);
@@ -386,7 +386,7 @@ int fmrib_main_float(int argc, char* argv[],const bool timeseriesMode, const str
}
if (mask.nvoxels()>0) cout << vol.percentile((float) n/100.0, mask) << " ";
else cout << vol.percentile((float) n/100.0) << " ";
- } else if (sarg=="-P") {
+ } else if (sarg=="-P") {
float n;
narg++;
if (narg<argc) {
@@ -400,8 +400,8 @@ int fmrib_main_float(int argc, char* argv[],const bool timeseriesMode, const str
exit(1);
}
if (mask.nvoxels()<1) {
- generate_masks(mask,masknz,vol,lthr,uthr);
- vol*=mask;
+ generate_masks(mask,masknz,vol,lthr,uthr);
+ vol*=mask;
}
cout << vol.percentile((float) n/100.0,masknz) << " ";
} else if (sarg=="-h") {
@@ -462,12 +462,12 @@ int fmrib_main_float(int argc, char* argv[],const bool timeseriesMode, const str
cerr << "Unrecognised option: " << sarg << endl;
exit(3);
}
-
+
narg++;
}
}
cout << endl;
-
+
}
return 0;
}
@@ -493,7 +493,7 @@ int main(int argc,char *argv[])
argc--;
}
if (argc < 3 )
- return print_usage(progname);
+ return print_usage(progname);
try {
return fmrib_main_float(argc,argv,timeseriesMode,indexMask);
} catch(std::exception &e) {
@@ -503,6 +503,5 @@ int main(int argc,char *argv[])
}
return -1;
-
-}
+}
--
GitLab