film_gls.cc

/*  film_gls.cc

    Mark Woolrich, FMRIB Image Analysis Group

    Copyright (C) 1999-2000 University of Oxford  */

/*  CCOPYRIGHT  */

#include <iostream>
#include <fstream>
#include <sstream>
#define WANT_STREAM
#define WANT_MATH

#include "newimage/newimageall.h"
#include "newmatap.h"
#include "newmatio.h"
#include "miscmaths/volumeseries.h"
#include "miscmaths/volume.h"
#include "miscmaths/miscmaths.h"
#include "utils/log.h"
#include "AutoCorrEstimator.h"
#include "paradigm.h"
#include "FilmGlsOptions.h"
#include "glimGls.h"
#include <vector>
#include <string>

using namespace NEWMAT;
using namespace FILM;
using namespace Utilities;
using namespace MISCMATHS;
using namespace NEWIMAGE;

int main(int argc, char *argv[])
{
  try{
    
    // Setup logging:
    Log& logger = LogSingleton::getInstance();

    // parse command line
    FilmGlsOptions& globalopts = FilmGlsOptions::getInstance();
    globalopts.parse_command_line(argc, argv, logger);

    // load data
    volume4D<float> input_data;
    volumeinfo vinfo;
    read_volume4D(input_data,globalopts.inputfname,vinfo);
    int sizeTS = input_data.tsize();

    Matrix datam;
    VolumeSeries x,y;
    y.read(globalopts.inputfname);
    x.setInfo(y.getInfo());
    x=input_data.matrix();
    x.thresholdSeries(globalopts.thresh, true);

    //save_volume_dtype(input_data[int(sizeTS/2)-1],logger.getDir() + "/" + globalopts.epifname,DT_SIGNED_INT );
    save_volume(input_data[int(sizeTS/2)-1],logger.getDir() + "/" + globalopts.epifname);

    volume<float> mask(input_data[0]);
    volume4D<float> epivol2(mask.xsize(),mask.ysize(),mask.zsize(),1);
    epivol2[0]=input_data[int(sizeTS/2)-1];
    mask=0;
    for(int t=0;t<input_data.tsize();t++) mask+=input_data[t];
    mask/=input_data.tsize();
    for(int z=0;z<input_data.zsize();z++)
      for(int y=0;y<input_data.ysize();y++)          
	for(int x=0;x<input_data.xsize();x++)
	  if ( mask(x,y,z) > globalopts.thresh )
	    for(int t=0;t<input_data.tsize();t++) input_data(x,y,z,t)-= mask(x,y,z);
    mask.binarise(globalopts.thresh,mask.max()+1,exclusive);
    input_data*=mask;
    datam=input_data.matrix(mask);
   
    int numTS = mask.sum();
    
    // Load paradigm:
    Paradigm parad;
    if(!globalopts.ac_only)
      {	
	parad.load(globalopts.paradigmfname, "", "", false, sizeTS);
      }
    else
      {
	// set design matrix to be one ev with all ones:
	Matrix mat(sizeTS,1);
	mat = 1;
	parad.setDesignMatrix(mat);
      }

    if(globalopts.verbose)
      {
	write_vest(logger.appendDir("Gc"), parad.getDesignMatrix());
      }
    

    OUT(parad.getDesignMatrix().Nrows());
    OUT(parad.getDesignMatrix().Ncols());
    OUT(sizeTS);
    OUT(numTS);

    // Setup GLM:
    int numParams = parad.getDesignMatrix().Ncols();
    GlimGls glimGls(numTS, sizeTS, numParams);

    // Residuals container:
    VolumeSeries residuals(sizeTS, numTS, x.getInfo(), x.getPreThresholdPositions());
    Matrix residualsm = datam;

    // Setup autocorrelation estimator:
    AutoCorrEstimator acEst(residuals);

    acEst.mask=mask;

    if(!globalopts.noest)
      {
	cout << "Calculating residuals..." << endl; 
	for(int i = 1; i <= numTS; i++)
	  {						    
            glimGls.setData(datam.Column(i), parad.getDesignMatrix(), i);
	    residuals.setSeries(glimGls.getResiduals(),i);
	  }
	cout << "Completed" << endl; 
	
	cout << "Estimating residual autocorrelation..." << endl; 
		
	if(globalopts.fitAutoRegressiveModel)
	  {
	    volume4D<float> beta;
	    beta.setmatrix(acEst.fitAutoRegressiveModel(),mask);
	    copybasicproperties(input_data,beta);
	    FslSetCalMinMax(&vinfo,beta.min(),beta.max());
	    save_volume4D(beta,LogSingleton::getInstance().getDir() + "/betas",vinfo);
	  }
	else if(globalopts.tukey)
	  {    
	    if(globalopts.tukeysize == 0)
	      globalopts.tukeysize = (int)(2*sqrt(sizeTS))/2;

	    acEst.calcRaw();

	    // SUSAN smooth raw estimates:
	    if(globalopts.smoothACEst)
	      {
		ColumnVector epivol = epivol2.matrix(mask).t();
		acEst.spatiallySmooth(logger.getDir() + "/" + globalopts.epifname, epivol, globalopts.ms, globalopts.epifname, globalopts.epith, globalopts.tukeysize);	
	      }	    
		
	    acEst.tukey(globalopts.tukeysize);
	  }
	else if(globalopts.multitaper)
	  {
	    acEst.calcRaw();
	    acEst.multitaper(int(globalopts.multitapersize));
	  }
	else if(globalopts.pava)
	  {
	    acEst.calcRaw();

	    // Smooth raw estimates:
	    if(globalopts.smoothACEst)
	      { 
		ColumnVector epivol = epivol2.matrix(mask).t();
		acEst.spatiallySmooth(logger.getDir() + "/" + globalopts.epifname, epivol, globalopts.ms, globalopts.epifname, globalopts.epith);
	      }
	    
	    acEst.pava();
	  }
	    
      }
    cout << "Completed" << endl; 

    cout << "Prewhitening and Computing PEs..." << endl;
    cout << "Percentage done:" << endl;
    int co = 1;

    Matrix mean_prewhitened_dm(parad.getDesignMatrix().Nrows(),parad.getDesignMatrix().Ncols());
    mean_prewhitened_dm=0;

    for(int i = 1; i <= numTS; i++)
      {						
	ColumnVector xw(sizeTS);
	ColumnVector xprew(sizeTS);
	   
	if(!globalopts.noest)
	  {
	    acEst.setDesignMatrix(parad.getDesignMatrix());
	    
	    // Use autocorr estimate to prewhiten data:
	    xprew = datam.Column(i); 
	    Matrix designmattw;
	    acEst.preWhiten(xprew, xw, i, designmattw);
	    
	    if ( (100.0*i)/numTS > co )
	      {
		cout << co << ",";
		cout.flush();
		co++;
	      }
	    
            datam.Column(i)=xw;
	    glimGls.setData(datam.Column(i), designmattw, i);
	    residuals.setSeries(glimGls.getResiduals(),i);

	    if(globalopts.output_pwdata || globalopts.verbose)
	      {
		mean_prewhitened_dm=mean_prewhitened_dm+designmattw;		
	      }
	  }
	else
	  {
	    if ( (100.0*i)/numTS > co )
	      {
		cout << co << ",";
		cout.flush();
		co++;
	      }
	    glimGls.setData(datam.Column(i), parad.getDesignMatrix(), i);
	    residuals.setSeries(glimGls.getResiduals(),i);

	    if(globalopts.output_pwdata || globalopts.verbose)
	      {
		mean_prewhitened_dm=mean_prewhitened_dm+parad.getDesignMatrix();		
	      }
	  }
      }

    if(globalopts.output_pwdata || globalopts.verbose)
      {
	mean_prewhitened_dm=mean_prewhitened_dm/numTS;
      }

    cerr << "Completed" << endl;

    cerr << "Saving results... " << endl;

    // write residuals
    residuals.writeThresholdedSeriesAsFloat(x.getInfo(),x.getPreThresholdPositions(),logger.getDir() + "/res4d");
    residuals.CleanUp();

    if(globalopts.output_pwdata || globalopts.verbose)
      {
	// Write out whitened data
        input_data.setmatrix(datam,mask);
	FslSetCalMinMax(&vinfo,input_data.min(),input_data.max());
        save_volume4D(input_data,logger.getDir() + "/prewhitened_data",vinfo);
	// Write out whitened design matrix
	write_vest(logger.appendDir("mean_prewhitened_dm.mat"), mean_prewhitened_dm);
		
      }

    x.CleanUp();

    // Write out threshac:
    Matrix& threshacm = acEst.getEstimates();
    int cutoff = sizeTS/2;
    if(globalopts.tukey)
      cutoff = globalopts.tukeysize;
    cutoff = MISCMATHS::Max(1,cutoff);
    threshacm = threshacm.Rows(1,cutoff); 
    volume4D<float> threshac;
    threshac.setmatrix(threshacm,mask);
    copybasicproperties(input_data,threshac);
    threshac.set_intent(NIFTI_INTENT_ESTIMATE,0,0,0);
    FslSetCalMinMax(&vinfo,threshac.min(),threshac.max());
    save_volume4D(threshac,logger.getDir() + "/threshac1",vinfo);
    threshacm.CleanUp();

    // save gls results:
    glimGls.Save(x.getInfo(), x.getPreThresholdPositions());
    glimGls.CleanUp();

    cerr << "Completed" << endl;
  }  
  catch(Exception p_excp) 
    {
      cerr << p_excp.what() << endl;
    }
  catch(...)
    {
      cerr << "Uncaught exception!" << endl;
    }

  return 0;
}