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

Candidate for Newimage film stable

parent 68bb6176
No related branches found
No related tags found
No related merge requests found
......@@ -70,7 +70,6 @@ namespace FILM {
void ContrastMgr::run()
{
Tracer_Plus ts("ContrastMgr::run");
Load();
// Loop through tcontrasts:
......@@ -116,31 +115,22 @@ namespace FILM {
Tracer_Plus ts("ContrastMgr::Load");
// Need to read in b, sigmaSquareds, corrections and dof
Log& logger = LogSingleton::getInstance();
// Load contrasts:
parad.load("", ContrastMgrOptions::getInstance().contrastfname, ContrastMgrOptions::getInstance().fcontrastfname, false, 0);
numParams = parad.getTContrasts().Ncols();
if(ContrastMgrOptions::getInstance().verbose)
{
cerr << "T Contrasts:" << endl << parad.getTContrasts();
cerr << "F Contrasts:" << endl << parad.getFContrasts();
}
volume4D<float> input_data;
read_volume4D(input_data,logger.getDir() + "/sigmasquareds",vinfo);
mask=input_data[0];
mask.binarise(0.0,exclusive);
//input_data*=mask;
sigmaSquareds=input_data.matrix(mask);
numTS = sigmaSquareds.Ncols();
mask.binarise(0.0,mask.max()+1,exclusive);
sigmaSquareds=input_data.matrix(mask).AsColumn();
numTS = sigmaSquareds.Nrows();
// b:
ColumnVector peVol;
b.ReSize(numTS, numParams);
......@@ -150,7 +140,7 @@ namespace FILM {
ostringstream osc;
osc << i;
read_volume4D(input_data,logger.getDir() + "/pe" + osc.str().c_str());
peVol=input_data.matrix(mask);
peVol=input_data.matrix(mask).AsColumn();
b.Column(i) = peVol;
}
......@@ -160,7 +150,7 @@ namespace FILM {
if(!in) //avw format
{
read_volume4D(input_data,logger.getDir() + "/dof");
dof=input_data.matrix(mask);
dof=input_data.matrix(mask).AsColumn();
}
else
{
......@@ -201,14 +191,13 @@ namespace FILM {
ostringstream osc;
osc << suffix << c_counter;
volume4D<float> output;
output.set_intent(NIFTI_INTENT_FTEST,0,0,0);
output.setmatrix(fstat);
output.setmatrix(fstat.AsRow(),mask);
FslSetCalMinMax(&vinfo,output.min(),output.max());
output.set_intent(NIFTI_INTENT_FTEST,0,0,0);
save_volume4D(output,logger.getDir() + "/fstat" + osc.str().c_str(),vinfo);
output.set_intent(NIFTI_INTENT_ZSCORE,0,0,0);
output.setmatrix(zstat);
output.setmatrix(zstat.AsRow(),mask);
FslSetCalMinMax(&vinfo,output.min(),output.max());
output.set_intent(NIFTI_INTENT_ZSCORE,0,0,0);
save_volume4D(output,logger.getDir() + "/zfstat" + osc.str().c_str(),vinfo);
}
......@@ -216,32 +205,31 @@ namespace FILM {
{
Tracer_Plus ts("ContrastMgr::SaveTContrast");
Log& logger = LogSingleton::getInstance();
// prepare contrast number:
ostringstream osc;
osc << suffix << c_counter;
volume4D<float> output;
output.set_intent(NIFTI_INTENT_NONE,0,0,0);
output.setmatrix(neff);
output.setmatrix(neff.AsRow(),mask);
FslSetCalMinMax(&vinfo,output.min(),output.max());
output.set_intent(NIFTI_INTENT_NONE,0,0,0);
save_volume4D(output,logger.getDir() + "/neff" + osc.str().c_str(),vinfo);
output.set_intent(NIFTI_INTENT_ESTIMATE,0,0,0);
output.setmatrix(cb);
output.setmatrix(cb.AsRow(),mask);
FslSetCalMinMax(&vinfo,output.min(),output.max());
save_volume4D(output,logger.getDir() + "/cope" + osc.str().c_str(),vinfo);
output.set_intent(NIFTI_INTENT_ESTIMATE,0,0,0);
output.setmatrix(varcb);
save_volume4D(output,logger.getDir() + "/cope" + osc.str().c_str(),vinfo);
output.setmatrix(varcb.AsRow(),mask);
FslSetCalMinMax(&vinfo,output.min(),output.max());
output.set_intent(NIFTI_INTENT_ESTIMATE,0,0,0);
save_volume4D(output,logger.getDir() + "/varcope" + osc.str().c_str(),vinfo);
output.set_intent(NIFTI_INTENT_TTEST,0,0,0);
output.setmatrix(tstat);
output.setmatrix(tstat.AsRow(),mask);
FslSetCalMinMax(&vinfo,output.min(),output.max());
output.set_intent(NIFTI_INTENT_TTEST,0,0,0);
cerr << output.intent_code() << endl;
save_volume4D(output,logger.getDir() + "/tstat" + osc.str().c_str(),vinfo);
output.set_intent(NIFTI_INTENT_ZSCORE,0,0,0);
output.setmatrix(zstat);
output.setmatrix(zstat.AsRow(),mask);
FslSetCalMinMax(&vinfo,output.min(),output.max());
output.set_intent(NIFTI_INTENT_ZSCORE,0,0,0);
save_volume4D(output,logger.getDir() + "/zstat" + osc.str().c_str(),vinfo);
}
......@@ -252,7 +240,6 @@ namespace FILM {
// puts ColumnVector of length p*p from correction
// into Matrix corr which is p*p:
corr.ReSize(numParams, numParams);
for (int i = 1; i <= numParams; i++)
{
for (int j = 1; j <= numParams; j++)
......
......@@ -9,7 +9,7 @@ LIBS = -lnewimage -lmiscmaths -lutils -lm -lnewmat -lfslio -lniftiio -lznz -lpro
XFILES = film_gls contrast_mgr ttoz ftoz ttologp
OBJS = ContrastMgrOptions.o ContrastMgr.o BandPassOptions.o gaussComparer.o AutoCorrEstimator.o glm.o paradigm.o glim.o FilmOlsOptions.o FilmGlsOptionsRes.o FilmGlsOptions.o glimGls.o
OBJS = ContrastMgrOptions.o ContrastMgr.o BandPassOptions.o gaussComparer.o AutoCorrEstimator.o glm.o paradigm.o FilmOlsOptions.o FilmGlsOptionsRes.o FilmGlsOptions.o glimGls.o
#OPTFLAGS =
......
......@@ -61,7 +61,6 @@ int main(int argc, char *argv[])
int numTS = datam.Ncols();
ColumnVector epivol = reference.matrix(mask).t();
//epivol = datam.Row(int(sizeTS/2)).AsColumn();
// Load paradigm:
Paradigm parad;
......@@ -240,7 +239,6 @@ int main(int argc, char *argv[])
Matrix& threshacm = acEst.getEstimates();
int cutoff = sizeTS/2;
if(globalopts.tukey) cutoff = globalopts.tukeysize;
//cutoff = MISCMATHS::Max(1,cutoff);
threshacm = threshacm.Rows(1,MISCMATHS::Max(1,cutoff));
input_data.setmatrix(threshacm,mask);
......
......@@ -26,7 +26,6 @@ public:
string fsfname;
string doffile1;
string doffile2;
string zscoresfname;
public:
globaloptions();
......
......@@ -11,11 +11,10 @@
#include "glimGls.h"
#include "miscmaths/miscmaths.h"
#include "utils/log.h"
#include "miscmaths/volume.h"
#include "miscmaths/volumeseries.h"
using namespace MISCMATHS;
using namespace Utilities;
using namespace NEWIMAGE;
namespace FILM {
......@@ -61,50 +60,7 @@ namespace FILM {
SetCorrection(inv_xx, ind);
}
void GlimGls::Save(const VolumeInfo& pvolinfo, const ColumnVector& prethreshpos)
{
// Need to save b, sigmaSquareds, corrections and dof
Log& logger = LogSingleton::getInstance();
VolumeInfo volinfo = pvolinfo;
// b:
Volume peVol;
for(int i = 1; i <= numParams; i++)
{
peVol = b.Row(i).AsColumn();
volinfo.intent_code = NIFTI_INTENT_ESTIMATE;
peVol.setInfo(volinfo);
peVol.setPreThresholdPositions(prethreshpos);
peVol.unthreshold();
// Add param number to "pe" to create filename:
ostringstream osc;
osc << i;
peVol.writeAsFloat(logger.getDir() + "/pe" + osc.str().c_str());
}
// sigmaSquareds:
volinfo.intent_code = NIFTI_INTENT_ESTIMATE;
/*sigmaSquareds.setInfo(volinfo);
sigmaSquareds.setPreThresholdPositions(prethreshpos);
sigmaSquareds.unthreshold();
sigmaSquareds.writeAsFloat(logger.getDir() + "/sigmasquareds");
// dof:
ColumnVector dofVec(1);
dofVec = dof;
write_ascii_matrix(logger.appendDir("dof"), dofVec);
// corrections (are the pes correlation matrix (x.t()*x).i() reshapen to a vector):
VolumeInfo newvolinfo = volinfo;
newvolinfo.v = numParams*numParams;
newvolinfo.intent_code = NIFTI_INTENT_NONE;
corrections.writeThresholdedSeriesAsFloat(newvolinfo,prethreshpos,logger.getDir() + "/corrections");*/
}
void GlimGls::Save(volumeinfo& vinfo, const volume<float>& mask, const float reftdim=1.0)
void GlimGls::Save(volumeinfo vinfo, const volume<float>& mask, const float reftdim=1.0)
{
// Need to save b, sigmaSquareds, corrections and dof
Log& logger = LogSingleton::getInstance();
......
......@@ -14,13 +14,9 @@
#define WANT_STREAM
#define WANT_MATH
#include "miscmaths/volume.h"
#include "miscmaths/volumeseries.h"
#include "newimage/newimageall.h"
using namespace NEWMAT;
using namespace MISCMATHS;
using namespace NEWIMAGE;
namespace FILM {
......@@ -32,8 +28,7 @@ namespace FILM {
GlimGls(const int pnumTS, const int psizeTS, const int pnumParams);
void setData(const ColumnVector& p_y, const Matrix& p_x, const int ind);
void Save(const VolumeInfo& volinfo, const ColumnVector& prethreshpos);
void Save(volumeinfo& vinfo, const volume<float>& mask,const float reftdim);
void Save(NEWIMAGE::volumeinfo vinfo, const NEWIMAGE::volume<float>& mask,const float reftdim);
ColumnVector& getResiduals() { return r; }
void CleanUp();
......
......@@ -133,9 +133,9 @@ int main(int argc,char *argv[])
ColumnVector ps(numTS);
T2z::ComputeZStats(vars, cbs, globalopts.dof, ps);
input.set_intent(NIFTI_INTENT_LOGPVAL,0,0,0);
input.setmatrix(ps.AsRow());
FslSetCalMinMax(&vinfo,input.min(),input.max());
input.set_intent(NIFTI_INTENT_LOGPVAL,0,0,0);
save_volume4D(input,globalopts.logpfname.c_str(),vinfo);
}
......
......@@ -137,9 +137,9 @@ int main(int argc,char *argv[])
ColumnVector zs(numTS);
T2z::ComputeZStats(varsm, cbsm, dofsm, zs);
input.set_intent(NIFTI_INTENT_ZSCORE,0,0,0);
input.setmatrix(zs.AsRow());
FslSetCalMinMax(&vinfo,input.min(),input.max());
input.set_intent(NIFTI_INTENT_ZSCORE,0,0,0);
save_volume4D(input,globalopts.zscoresfname,vinfo);
}
catch(Exception p_excp)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment