meldata.h

/*  MELODIC - Multivariate exploratory linear optimized decomposition into 
              independent components
    
    meldata.h - data container class

    Christian F. Beckmann, FMRIB Image Analysis Group
    
    Copyright (C) 1999-2002 University of Oxford */

/*  CCOPYRIGHT */


#ifndef __MELODICDATA_h
#define __MELODICDATA_h

#include "newimage/newimageall.h"
#include "utils/log.h"
#include "meloptions.h"
#include "melhlprfns.h"

using namespace Utilities;
using namespace NEWIMAGE;

namespace Melodic{
  
  class MelodicData
    {
    public:

      //constructor
      MelodicData(MelodicOptions &popts, Log &plogger):  
	opts(popts),
	logger(plogger)
	{	
	  after_mm = false;
	  Resels = 0;
	}  
 
      void save();

      Matrix process_file(string fname, int numfiles);

      inline void save4D(Matrix what, string fname){
	 volume4D<float> tempVol;
	 tempVol.setmatrix(what,Mask);
	 save_volume4D(tempVol,logger.appendDir(fname),tempInfo);
	 message("  " << logger.appendDir(fname) << endl);
      }
      
      inline void saveascii(Matrix what, string fname){
	 write_ascii_matrix(logger.appendDir(fname),what);   
	 message("  " << logger.appendDir(fname) << endl);   
      }
 
      inline void savebinary(Matrix what, string fname){
         write_binary_matrix(what,logger.appendDir(fname));  
	 message("  " << logger.appendDir(fname) << endl);    
      }

      int  remove_components();
      void setup();
      void status(const string &txt);

      inline Matrix& get_pcaE() {return pcaE;}
      inline void set_pcaE(Matrix& Arg) {pcaE = Arg;}

      inline RowVector& get_pcaD() {return pcaD;}
      inline void set_pcaD(RowVector& Arg) {pcaD = Arg;}

      inline Matrix& get_data() {return Data;}
      inline void set_data(Matrix& Arg) {Data = Arg;}

      inline Matrix& get_IC() {return IC;}
      inline void set_IC(Matrix& Arg) {IC = Arg;}
      inline void set_IC(int ctr, Matrix& Arg) {IC.Row(ctr) = Arg;}
      
      inline vector<Matrix>& get_Smodes() {return Smodes;}
      inline Matrix& get_Smodes(int what) {return Smodes.at(what);}
      inline void add_Smodes(Matrix& Arg) {Smodes.push_back(Arg);}      
      inline void save_Smodes(){
	Matrix tmp = Smodes.at(0); 
	for(unsigned int ctr = 1; ctr < Smodes.size(); ctr++)
	  tmp |= Smodes.at(ctr);
	saveascii(tmp,opts.outputfname.value() + "_Smodes");
      }

      inline vector<Matrix>& get_Tmodes() {return Tmodes;}
      inline Matrix& get_Tmodes(int what) {return Tmodes.at(what);}
      inline void add_Tmodes(Matrix& Arg) {Tmodes.push_back(Arg);}
      inline void save_Tmodes(){
	Matrix tmp = Tmodes.at(0); 
	for(unsigned int ctr = 1; ctr < Tmodes.size(); ctr++)
	  tmp |= Tmodes.at(ctr);
	saveascii(tmp,opts.outputfname.value() + "_Tmodes");
      }

      void set_TSmode();

      inline Matrix& get_white() {return whiteMatrix;}
      inline void set_white(Matrix& Arg) {whiteMatrix = Arg;}
      
      inline Matrix& get_dewhite() {return dewhiteMatrix;}
      inline void set_dewhite(Matrix& Arg) {dewhiteMatrix = Arg;}
      
      inline Matrix& get_meanC() {return meanC;}
      inline Matrix& get_meanR() {return meanR;}

      inline Matrix& get_stdDevi() {return stdDevi;}
      inline void set_stdDevi(Matrix& Arg) {stdDevi = Arg;}
  
      inline Matrix& get_mix() {return mixMatrix;}

      inline void set_mix(Matrix& Arg) {mixMatrix = Arg;}
      Matrix expand_mix(); 
      Matrix expand_dimred(const Matrix& Mat); 
      Matrix reduce_dimred(const Matrix& Mat); 
      
      inline Matrix& get_fmix() {return mixFFT;}
      inline void set_fmix(Matrix& Arg) {mixFFT = Arg;}

      inline Matrix& get_unmix() {return unmixMatrix;}
      inline void set_unmix(Matrix& Arg) {unmixMatrix = Arg;}

      inline volume<float>& get_mask() {return Mask;}
      inline void set_mask(volume<float>& Arg) {Mask = Arg;}
  
      inline volume<float>& get_mean() {return Mean;}
      inline void set_mean(volume<float>& Arg) {Mean = Arg;}
   
      inline Matrix& get_Data() {return Data;}
      inline void set_Data(Matrix& Arg) {Data = Arg;}
    
      inline Matrix& get_RXweight() {return RXweight;}
      inline void set_RXweight(Matrix& Arg) {RXweight = Arg;}
 
      inline Matrix& get_ICstats() {return ICstats;}
      inline void set_ICstats(Matrix& Arg) {ICstats = Arg;}
     
      inline Matrix& get_EVP() {return EVP;}
      inline void set_EVP(Matrix& Arg) {EVP = Arg;}
      
      inline Matrix& get_EV() {return EV;}
      inline void set_EV(Matrix& Arg) {EV = Arg;}

      inline Matrix& get_PPCA() {return PPCA;}
      inline void set_PPCA(Matrix& Arg) {PPCA = Arg;}

      inline Matrix& get_stdNoisei() {return stdNoisei;}
      inline void set_stdNoisei(Matrix& Arg) {stdNoisei = Arg;}

      inline int data_dim() {return Data.Nrows();}
      inline int data_samples() {return Data.Ncols();}
     
      inline float get_resels() {return Resels;}
      inline void set_resels(float& Arg) {Resels = Arg;}

      inline int get_numfiles() {return numfiles;}

      inline void set_after_mm(bool val) {after_mm = val;}

      inline void flipres(int num){
	IC.Row(num) = -IC.Row(num);
	mixMatrix.Column(num) = -mixMatrix.Column(num);
	mixFFT=calc_FFT(mixMatrix);
	unmixMatrix = pinv(mixMatrix);

        if(ICstats.Storage()>0&&ICstats.Ncols()>3){
	  double tmp;
	  tmp = ICstats(num,3);
	  ICstats(num,3) = -1.0*ICstats(num,4);
	  ICstats(num,4) = -1.0*tmp;
	}
      }
      
      void sort();

      volumeinfo tempInfo;

      vector<Matrix> DWM, WM;

    private:
      MelodicOptions &opts;     
      Log &logger;       

      Matrix pcaE;
      RowVector pcaD;
      Matrix whiteMatrix;
      Matrix dewhiteMatrix;
      Matrix meanC;
      Matrix meanR;
      Matrix stdDev;
      Matrix stdDevi;
      Matrix RXweight;
      Matrix mixMatrix;
      Matrix unmixMatrix;
      Matrix mixFFT;
      Matrix IC;
      Matrix ICstats;

      vector<Matrix> Tmodes;
      vector<Matrix> Smodes;

      Matrix EVP;
      Matrix EV;
      Matrix stdNoisei;
      volume<float> Mask;
      volume<float> Mean;
      Matrix Data;
      Matrix PPCA;
      Matrix jointCC;
      
      bool after_mm;

      float Resels;
      int numfiles;

      char Mean_fname[1000];

      void setup_misc();
      void create_mask(volume<float>& theMask);
      void create_RXweight();
      void est_smoothness();

      unsigned long standardise(volume<float>& mask, 
				volume4D<float>& R);
      float est_resels(volume4D<float> R, volume<float> mask);
    };
}

#endif