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 <strstream>
#define WANT_STREAM
#define WANT_MATH

#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;

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
    VolumeSeries x;
    x.read(globalopts.inputfname);

    int sizeTS = x.getNumVolumes();

    // if needed for SUSAN spatial smoothing, output the 12th volume for use later
    Volume epivol;
    if(globalopts.smoothACEst)
      {
	epivol = x.getVolume(int(sizeTS/2)).AsColumn();
	epivol.setInfo(x.getInfo());
	
	epivol.writeAsInt(logger.getDir() + "/" + globalopts.epifname);
      }

    // This also removes the mean from each of the time series:
    x.thresholdSeries(globalopts.thresh, true);

    // if needed for SUSAN spatial smoothing, also threshold the epi volume
    if(globalopts.smoothACEst)
      {
	epivol.setPreThresholdPositions(x.getPreThresholdPositions());
	epivol.threshold();
      }
    
    int numTS = x.getNumSeries();

    // 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_ascii_matrix(logger.appendDir("Gc"), parad.getDesignMatrix());
      }
    
    // Setup GLM:
    int numParams = parad.getDesignMatrix().Ncols();
    GlimGls glimGls(numTS, sizeTS, numParams);

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

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

    if(!globalopts.noest)
      {
	cout << "Calculating residuals..." << endl; 
	for(int i = 1; i <= numTS; i++)
	  {						    
	    glimGls.setData(x.getSeries(i), parad.getDesignMatrix(), i);
	    residuals.setSeries(glimGls.getResiduals(),i);
	  }
	cout << "Completed" << endl; 
	
	cout << "Estimating residual autocorrelation..." << endl; 
		
	if(globalopts.fitAutoRegressiveModel)
	  {
	    acEst.fitAutoRegressiveModel();
	  }
	else if(globalopts.tukey)
	  {    
	    if(globalopts.tukeysize == 0)
	      globalopts.tukeysize = (int)(2*sqrt(sizeTS))/2;

	    acEst.calcRaw();

	    // SUSAN smooth raw estimates:
	    if(globalopts.smoothACEst)
	      {
		acEst.spatiallySmooth(logger.getDir() + "/" + globalopts.epifname, epivol, globalopts.ms, globalopts.epifname, globalopts.susanpath, 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)
	      { 
		acEst.spatiallySmooth(logger.getDir() + "/" + globalopts.epifname, epivol, globalopts.ms, globalopts.epifname, globalopts.susanpath, globalopts.epith);		    
	      }
	    
	    acEst.pava();
	  }
	    
      }
    cout << "Completed" << endl; 

    cout << "Prewhitening and Computing PEs..." << endl;
    cout << "Percentage done:" << endl;
    int co = 1;
    
    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 = x.getSeries(i);	
	    Matrix designmattw;
	    acEst.preWhiten(xprew, xw, i, designmattw);
	    
	    if ( (100.0*i)/numTS > co )
	      {
		cout << co << ",";
		cout.flush();
		co++;
	      }
	    
	    x.setSeries(xw,i);
	    
	    glimGls.setData(x.getSeries(i), designmattw, i);
	  }
	else
	  {
	    if ( (100.0*i)/numTS > co )
	      {
		cout << co << ",";
		cout.flush();
		co++;
	      }
	    
	    glimGls.setData(x.getSeries(i), parad.getDesignMatrix(), i);
	    residuals.setSeries(glimGls.getResiduals(),i);
	  }
      }

    cerr << "Completed" << endl;

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

    if(globalopts.verbose)
      {
	// output whitened x
	x.unthresholdSeries();
	x.writeAsFloat(logger.getDir() + "/wx");
      }
    x.CleanUp();

    // Write out threshac:
    VolumeSeries& threshac = acEst.getEstimates();
    int cutoff = sizeTS/2;
    if(globalopts.tukey)
      cutoff = globalopts.tukeysize;
    cutoff = MISCMATHS::Max(3,cutoff);
    threshac = threshac.Rows(1,cutoff);
    VolumeInfo volinfo = x.getInfo();
    volinfo.v = cutoff;
    threshac.unthresholdSeries(volinfo,x.getPreThresholdPositions());
    threshac.writeAsFloat(logger.getDir() + "/threshac1");
    threshac.thresholdSeries();
    threshac.CleanUp();

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

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

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

  return 0;
}