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-2008 University of Oxford */

/*  Part of FSL - FMRIB's Software Library
    http://www.fmrib.ox.ac.uk/fsl
    fsl@fmrib.ox.ac.uk
    
    Developed at FMRIB (Oxford Centre for Functional Magnetic Resonance
    Imaging of the Brain), Department of Clinical Neurology, Oxford
    University, Oxford, UK
    
    
    LICENCE
    
    FMRIB Software Library, Release 4.0 (c) 2007, The University of
    Oxford (the "Software")
    
    The Software remains the property of the University of Oxford ("the
    University").
    
    The Software is distributed "AS IS" under this Licence solely for
    non-commercial use in the hope that it will be useful, but in order
    that the University as a charitable foundation protects its assets for
    the benefit of its educational and research purposes, the University
    makes clear that no condition is made or to be implied, nor is any
    warranty given or to be implied, as to the accuracy of the Software,
    or that it will be suitable for any particular purpose or for use
    under any specific conditions. Furthermore, the University disclaims
    all responsibility for the use which is made of the Software. It
    further disclaims any liability for the outcomes arising from using
    the Software.
    
    The Licensee agrees to indemnify the University and hold the
    University harmless from and against any and all claims, damages and
    liabilities asserted by third parties (including claims for
    negligence) which arise directly or indirectly from the use of the
    Software or the sale of any products based on the Software.
    
    No part of the Software may be reproduced, modified, transmitted or
    transferred in any form or by any means, electronic or mechanical,
    without the express permission of the University. The permission of
    the University is not required if the said reproduction, modification,
    transmission or transference is done without financial return, the
    conditions of this Licence are imposed upon the receiver of the
    product, and all original and amended source code is included in any
    transmitted product. You may be held legally responsible for any
    copyright infringement that is caused or encouraged by your failure to
    abide by these terms and conditions.
    
    You are not permitted under this Licence to use this Software
    commercially. Use for which any financial return is received shall be
    defined as commercial use, and includes (1) integration of all or part
    of the source code or the Software into a product for sale or license
    by or on behalf of Licensee to third parties or (2) use of the
    Software or any derivative of it for research with the final aim of
    developing software products for sale or license to a third party or
    (3) use of the Software or any derivative of it for research with the
    final aim of developing non-software products for sale or license to a
    third party, or (4) use of the Software to provide any service to an
    external organisation for which payment is received. If you are
    interested in using the Software commercially, please contact Isis
    Innovation Limited ("Isis"), the technology transfer company of the
    University, to negotiate a licence. Contact details are:
    innovation@isis.ox.ac.uk quoting reference DE/1112. */


#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 = 1);

      inline void save4D(Matrix what, string fname){
	 			volume4D<float> tempVol;
	 			tempVol.setmatrix(what,Mask);
	 			save_volume4D(tempVol,logger.appendDir(fname));
	 			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_param() {return param;} 
      inline void set_param(Matrix& Arg) {param = Arg;}	

      inline Matrix& get_paramS() {return paramS;} 
      inline void set_paramS(Matrix& Arg) {paramS = Arg;}	

      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;
				if (Tmodes.size() < 1)
	  			for (int ctr = 1; ctr <= Arg.Ncols(); ctr++){
	    			Matrix tmp = Arg.Column(ctr);
	    			add_Tmodes(tmp);
	  			}
      }

      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 volume<float>& get_bg() {
				if(opts.bgimage.value()>"")
					return background;
				else
					return Mean;
			}
      inline void set_bg(volume<float>& Arg) {background = 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) {if(EVP.Storage()==0)
																					EVP = Arg;}
      
      inline Matrix& get_EV() {return EV;}
      inline void set_EV(Matrix& Arg) {if(EV.Storage()==0)
																				 EV = Arg;}

      inline Matrix& get_PPCA() {return PPCA;}
      inline void set_PPCA(Matrix& Arg) {if(PPCA.Storage()==0)
																					 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();
	  void reregress();

      vector<Matrix> DWM, WM;
			basicGLM glmT, glmS;
			Matrix Tdes, Tcon, TconF, Sdes, Scon, SconF, param, paramS;	
			RowVector explained_var;		

    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;
      volume<float> background;

      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