Skip to content
Snippets Groups Projects
Commit 61d6cfff authored by Matthew Webster's avatar Matthew Webster
Browse files

read4D to avoid calling legacy read

parent d4bd099c
Branches cvsHEAD
Tags cvsFINAL
No related merge requests found
Hide whitespace changes
Inline Side-by-side
fslchfiletype.cc
+ 1
− 1
View file @ 61d6cfff
......@@ -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++ )
......
fslcomplex.cc
+ 18
− 19
View file @ 61d6cfff
......@@ -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;
}
fslstats.cc
+ 53
− 54
View file @ 61d6cfff
......@@ -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;
}
}
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