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 "glim.h"
#include "sigproc.h"
#include "miscmaths.h"
#include "gaussComparer.h"
#include "Log.h"
#include "AutoCorrEstimator.h"
#include "paradigm.h"
#include "FilmGlsOptions.h"
#include "glimGls.h"
#include <vector>
#include <string>

#ifndef NO_NAMESPACE
using namespace NEWMAT;
using namespace SIGPROC;
using namespace TACO;
using namespace UTILS;
#endif

int main(int argc, char *argv[])
{
  try{
    rand();
    // parse command line to find out directory name for logging:
    ofstream out2;
    FilmGlsOptions& globalopts = FilmGlsOptions::getInstance();
    globalopts.parse_command_line(argc, argv, out2);

    // Setup logging:
    Log& logger = Log::getInstance();
    logger.setLogFile("glslogfile");
    logger.establishDir(globalopts.datadir);

    // parse command line again to output arguments to logfile
    globalopts.parse_command_line(argc, argv, logger.str());

    // load non-temporally filtered data
    VolumeSeries x;
    x.read(globalopts.inputfname);

    // if needed 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 later 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);

    // Sort out detrending:
    if(globalopts.detrend)
      {
	SIGPROC::Detrend(x, false);
      }

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

    vector<Matrix> dms;
    for(int i = 1; i <= numTS; i++)
      {
	dms.push_back(parad.getDesignMatrix());
      }
    
    // Setup OLS GLM for temporally filtered data:
    int numParams = parad.getDesignMatrix().Ncols();
    GlimGls glimGls(numTS, sizeTS, numParams);
    
    VolumeSeries residuals(sizeTS, numTS, x.getDims(), x.getPreThresholdPositions());
    AutoCorrEstimator acEst(residuals);

    if(!globalopts.noest)
      {
	int numiters = globalopts.numiters+1;
	int iters = 1;

	// iters==1 is for high freq removal
	if(!globalopts.highfreqremoval) 
	  iters++;

	for(; iters <= numiters; iters++)
	  { 
	    cerr << "iters = " << iters << endl;
	    // Loop through voxels:
	    cerr << "Calculating residuals..." << endl; 
	    for(int i = 1; i <= numTS; i++)
	      {						    
		glimGls.setData(x.getSeries(i), dms[i-1], i);
		residuals.getSeries(i) = glimGls.getResiduals();
	      }
	    cerr << "Completed" << endl; 

	    cerr << "Estimating residual autocorrelation..." << endl; 
	
	    // Estimate Autocorrelation:	    
	    acEst.calcRaw();

	    if(iters==1)
	      {
		// iters==1 is for high freq removal
		acEst.tukey(10);
	      }
	    else
	      {
		if(globalopts.fitAutoRegressiveModel)
		  {
		    acEst.fitAutoRegressiveModel();
		  }
		else if(globalopts.tukey)
		  {    
		    // 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; 

	    // Loop through voxels:
	    cerr << "Prewhitening..." << endl;   
	    int co = 1;
	    
	    for(int i = 1; i <= numTS; i++)
	      {					
		ColumnVector I(sizeTS);
		I = 0;
		I(1) = 1;
	    
		ColumnVector xw(sizeTS);
		ColumnVector xprew(sizeTS);
	    
		acEst.setDesignMatrix(dms[i-1]);
	    
		// Use autocorr estimate to prewhiten data:
		xprew = x.getSeries(i);
	    
		Matrix designmattw;
	    
		// iters==1 is for high freq removel
		acEst.preWhiten(xprew, xw, i, designmattw, bool(iters==1));
	    
		if(co > 1000)
		  {
		    co = 1;
		    cerr << (float)i/(float)numTS << ",";
		  }
		else
		  co++;
	   
		x.getSeries(i) = xw;
		dms[i-1] = designmattw;
	      }  
      	
	    cerr << "Completed" << endl;
	
	    // Add param number to "pe" to create filename:
	    char strc[4];
	    ostrstream osc(strc,4);
	    osc << iters - 1 << '\0';
	
	    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() + "/threshac" + strc);
	    threshac.thresholdSeries();

	    if(globalopts.verbose)
	      {	
		cerr << "Saving results... " << endl;
	    
		ColumnVector& countLargeE = acEst.getCountLargeE();
	
		logger.out(string("countLargeE") + strc, countLargeE);

		residuals.unthresholdSeries(x.getDims(),x.getPreThresholdPositions());
		residuals.writeAsFloat(logger.getDir() + "/res4d" + strc);
		residuals.thresholdSeries();
		cerr << "Completed" << endl;
	      }
	  }
      }
    else // no estimation of autocorrelations
      {
	// do nothing
      }

    // Do once more to compute real param ests:
    cerr << "Computing parameter estimates..." << endl;
    for(int i = 1; i <= numTS; i++)
      {						    
	glimGls.setData(x.getSeries(i), dms[i-1], i);
	
	if(globalopts.verbose)
	  {
	    residuals.getSeries(i) = glimGls.getResiduals();
	  }
      }
    cerr << "Completed" << endl;

    // Write out necessary data:
    cerr << "Saving results... " << endl;

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

    // write out design matrices - a volume Series for each param
    if(globalopts.verbose_dms)
      {
	VolumeSeries::Dims dims = x.getDims();
	for(int j = 1; j <= numParams; j++)
	  {	  
	    char strc[4];
	    ostrstream osc(strc,4);
	    osc << j << '\0';
	
	    VolumeSeries dmsmat(sizeTS, numTS);
	
	    for(int i = 1; i <= numTS; i++)
	      {
		dmsmat.getSeries(i) = dms[i-1].Column(j);
	      }
	    dmsmat.setDims(dims);
	    dmsmat.setPreThresholdPositions(x.getPreThresholdPositions());
	    dmsmat.unthresholdSeries();
	    dmsmat.writeAsFloat(logger.getDir() + "/dms" + strc);
	  }
	// output x
	x.unthresholdSeries();
	x.writeAsFloat(logger.getDir() + "/wx");
      }

    if(globalopts.verbose)
      {	
	VolumeSeries::Dims dims = x.getDims();

       	// Save E
	VolumeSeries& E = acEst.getE();
	dims.v = acEst.getZeroPad();
	E.setDims(dims);
	E.setPreThresholdPositions(x.getPreThresholdPositions());
	E.unthresholdSeries();	
	E.writeAsFloat(logger.getDir() + "/E");

      }

    glimGls.Save(x.getDims(), x.getPreThresholdPositions());
    cerr << "Completed" << endl;
  }  
  catch(Exception p_excp) 
    {
      cerr << p_excp.what() << endl;
    }
  catch(...) 
    {
      cerr << "Image error" << endl;
    } 
  return 0;
}