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 "VolumeSeries.h"
#include "Volume.h"
#include "miscmaths.h"
#include "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;

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

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

    // load data
    VolumeSeries x;
    x.read(globalopts.inputfname);

    // if needed for SUSAN spatial smoothing, output the 12th volume for use later
    Volume epivol;
    if(globalopts.smoothACEst)
      {
	epivol = x.getVolume(12).AsColumn();
	epivol.setDims(x.getDims());
	
	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 sizeTS = x.getNumVolumes();
    int numTS = x.getNumSeries();
    
    // Load paradigm: 
    Paradigm parad;
    parad.load(globalopts.paradigmfname, "", "", false, sizeTS);

    if(globalopts.verbose)
      {
	logger.out("Gc", parad.getDesignMatrix());
      }
    
    // Setup GLM:
    int numParams = parad.getDesignMatrix().Ncols();
    GlimGls glimGls(numTS, sizeTS, numParams);
    
    // Residuals container:
    VolumeSeries residuals(sizeTS, numTS, x.getDims(), x.getPreThresholdPositions());

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

    if(!globalopts.noest)
      {
	cerr << "Calculating residuals..." << endl; 
	for(int i = 1; i <= numTS; i++)
	  {						    
	    glimGls.setData(x.getSeries(i), parad.getDesignMatrix(), i);
	    residuals.getSeries(i) = glimGls.getResiduals();
	  }
	cerr << "Completed" << endl; 
	
	cerr << "Estimating residual autocorrelation..." << endl; 
	
	acEst.calcRaw();
	
	if(globalopts.fitAutoRegressiveModel)
	  {
	    acEst.fitAutoRegressiveModel();
	  }
	else if(globalopts.tukey)
	  {    
	    // SUSAN smooth raw estimates:
	    if(globalopts.smoothACEst)
	      {
		acEst.spatiallySmooth(logger.getDir() + "/" + globalopts.epifname, epivol, globalopts.ms, globalopts.epifname, globalopts.susanpath, globalopts.epith);		    
	      }
	    if(globalopts.tukeysize == 0)
	      globalopts.tukeysize = (int)(2*sqrt(sizeTS))/2;
		
	    acEst.tukey(globalopts.tukeysize);
	  }
	else if(globalopts.multitaper)
	  {
	    acEst.multitaper(globalopts.multitapersize);
	  }
	else if(globalopts.pava)
	  {
	    // Smooth raw estimates:
	    if(globalopts.smoothACEst)
	      { 
		acEst.spatiallySmooth(logger.getDir() + "/" + globalopts.epifname, epivol, globalopts.ms, globalopts.epifname, globalopts.susanpath, globalopts.epith);		    
	      }
	    
	    acEst.pava();
	  }
	    
      }
    cerr << "Completed" << endl; 

    cerr << "Prewhitening and Computing PEs..." << endl;   
    int co = 1;
    for(int i = 1; i <= numTS; i++)
      {						
	ColumnVector xw(sizeTS);
	ColumnVector xprew(sizeTS);
	    
	acEst.setDesignMatrix(parad.getDesignMatrix());
	
	// Use autocorr estimate to prewhiten data:
	xprew = x.getSeries(i);	
	Matrix designmattw;
	acEst.preWhiten(xprew, xw, i, designmattw);
	    
	if(co > 1000)
	  {
	    co = 1;
	    cerr << (float)i/(float)numTS << ",";
	  }
	else
	  co++;
	   
	x.getSeries(i) = xw;
    
	glimGls.setData(x.getSeries(i), designmattw, 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;
    threshac = threshac.Rows(1,cutoff);
    VolumeSeries::Dims dims = x.getDims();
    dims.v = cutoff;
    threshac.unthresholdSeries(dims,x.getPreThresholdPositions());
    threshac.writeAsFloat(logger.getDir() + "/threshac1");
    threshac.thresholdSeries();
    threshac.CleanUp();

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

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

    cerr << "Completed" << endl;
  }  
  catch(Exception p_excp) 
    {
      cerr << p_excp.what() << endl;
    }
  catch(...) 
    {
      cerr << "Image error" << endl;
    } 
  return 0;
}