paradigm.cc

/*  paradigm.cc

    Mark Woolrich, FMRIB Image Analysis Group

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

/*  CCOPYRIGHT  */

#include "paradigm.h"
#include "utils/log.h"
#include "miscmaths/miscmaths.h"
#include <fstream>

using namespace FILM;

namespace FILM {

  void Paradigm::read_vest_waveform(string p_fname, Matrix& p_mat)
  {
    ifstream in;
    in.open(p_fname.c_str(), ios::in);
    
    if(!in)
      {
	cerr << "Unable to open " << p_fname << endl;
	throw;
      }
    
    int numWaves = 0;
    int numPoints = 0;
  
    string str;

    while(true)
      {
	in >> str;
	if(str == "/Matrix")
	  break;
	else if(str == "/NumWaves")
	  {
	    in >> numWaves;
	  }
	else if(str == "/NumPoints" || str == "/NumContrasts")
	  {
	    in >> numPoints;
	  }
      }

    if(numWaves != 0)
      {
	p_mat.ReSize(numPoints, numWaves);
      }
    else
      {
	numWaves = p_mat.Ncols();
	numPoints = p_mat.Nrows();
      }
      
    for(int i = 1; i <= numPoints; i++)
      {
	for(int j = 1; j <= numWaves; j++)    
	  {
	    in >> p_mat(i,j);
	  }
      }

    in.close();
  }
  
  void Paradigm::load(const string& p_paradfname, const string& p_tcontrastfname, const string& p_fcontrastfname, bool p_blockdesign, int p_sizets)
    {
      if(p_paradfname != "")
	read_vest_waveform(p_paradfname, designMatrix);

      if(p_tcontrastfname != "")
	read_vest_waveform(p_tcontrastfname, tcontrasts);       

      if(p_fcontrastfname != "")
	read_vest_waveform(p_fcontrastfname, fcontrasts);   
       
      // Check time series match up with design
      if(p_paradfname != "" && designMatrix.Nrows() != p_sizets)
	{
	  cerr << "Num scans = "<< p_sizets << ", design matrix rows = " << designMatrix.Nrows() << endl;
	  throw Exception("size of design matrix does not match number of scans\n");
	}
      
      // Check contrasts match
      if(p_tcontrastfname != "" &&  p_fcontrastfname != "" && tcontrasts.Nrows() != fcontrasts.Ncols())
	{ 
	  cerr << "tcontrasts.Nrows()  = " << tcontrasts.Nrows() << endl;
	  cerr << "fcontrasts.Ncols()  = " << fcontrasts.Ncols() << endl;
	  throw Exception("size of tcontrasts does not match fcontrasts\n");
	}

      // Check contrast matches design matrix
      if(p_paradfname != "" &&  p_tcontrastfname != "" && designMatrix.Ncols() != tcontrasts.Ncols())
	 { 
	  cerr << "Num tcontrast cols  = " << tcontrasts.Ncols() << ", design matrix cols = " 
	    << designMatrix.Ncols() << endl;
	  throw Exception("size of design matrix does not match t contrasts\n");
	 }

      if(p_blockdesign)
	{
	  // Need to adjust amplitude from (-1,1) to (0,1)
	  designMatrix = (designMatrix + 1)/2;
	}
    }
  
}