/* 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; } } }