meldata.cc

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

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

/*  CCOPYRIGHT  */

#include "newimage/newimageall.h"
#include "meloptions.h"
#include "meldata.h"
#include "melodic.h" 
#include "utils/log.h"
#include <time.h>
#include <algorithm>
#include "miscmaths/miscprob.h"
#include "melhlprfns.h" 

using namespace Utilities;
using namespace NEWIMAGE;
 
namespace Melodic{
// {{{ Setup

  ReturnMatrix MelodicData::process_file(string fname, int numfiles)
  {
	Matrix tmpData;
	{
    	volume4D<float> RawData;

		memmsg(" before reading file "<< fname);

    	//read data
    	message("Reading data file " << fname << "  ... ");
    	read_volume4D(RawData,fname);
    	message(" done" << endl);
		memmsg(" after reading file "<< fname);
		
		del_vols(RawData,opts.dummy.value());
    
    	Mean += meanvol(RawData)/numfiles;

		//estimate smoothness
		memmsg(" before est smoothness ");
    	if((Resels == 0)&&(!opts.filtermode))
      		Resels = est_resels(RawData,Mask);
		memmsg(" after smoothness ");
	
    	//reshape
		memmsg(" before reshape ");
    	tmpData = RawData.matrix(Mask);
		memmsg(" after reshape ");	  
	}    
        
    //convert to percent BOLD signal change
    if(opts.pbsc.value()){
      message("  Converting data to percent BOLD signal change ...");
      Matrix meanimg = convert_to_pbsc(tmpData);
      meanR = meanimg.Row(1);
      message(" done" << endl);
    } 
	else{
		if(opts.remove_meanvol.value())
		{	      
			message(string("  Removing mean image ..."));
			memmsg(" before remmean ");
      		remmean(tmpData,meanR,1);
			memmsg(" after remmean ");
      		message(" done" << endl);
		}
		else meanR=ones(1,tmpData.Ncols());
    }

	if(opts.remove_meantc.value()){
		remmean(tmpData,meanC,2);
	}
	
    //convert to power spectra
    if(opts.pspec.value()){
      message("  Converting data to powerspectra ...");
      tmpData = calc_FFT(tmpData);
      message(" done" << endl);
    }
	
    //switch dimension in case temporal ICA is required
    if(opts.temporal.value()){
      message(string("  Switching dimensions for temporal ICA") << endl);
      tmpData = tmpData.t();
      Matrix tmp;
      tmp = meanC;
      meanC = meanR.t();
      meanR = tmp.t();
      message("  Data size : " << Data.Nrows() << " x " << Data.Ncols() <<endl);
    }
      
    //variance - normalisation
	memmsg(" before VN ");
    if(opts.varnorm.value()){
      message("  Normalising by voxel-wise variance ..."); 
			if(stdDev.Storage()==0)
      	stdDev = varnorm(tmpData,std::min(30,tmpData.Nrows()-1),
					opts.vn_level.value())/numfiles;
			else 	
				stdDev += varnorm(tmpData,std::min(30,tmpData.Nrows()-1),
					opts.vn_level.value())/numfiles;
      stdDevi = pow(stdDev,-1); 
	  memmsg(" in VN ");
      message(" done" << endl);
    }

	tmpData.Release();
    return tmpData;
  }

  Matrix MelodicData::expand_mix()
  {
    Matrix out;
    out = expand_dimred(mixMatrix);
    return out;
  }

  Matrix MelodicData::expand_dimred(const Matrix& Mat)
  {
    int first, last;
    first = 1;
    last = DWM.at(0).Ncols();
    Matrix tmp = DWM.at(0) * Mat.Rows(first,last);
    for(unsigned int ctr = 1; ctr < DWM.size(); ctr++){
      first = last + 1;
      last += DWM.at(ctr).Ncols();
      tmp &= DWM.at(ctr) * Mat.Rows(first, last);
    }
    return tmp;
  }

  Matrix MelodicData::reduce_dimred(const Matrix& Mat)
  {
    int first, last;
    first = 1;
    last = WM.at(0).Ncols();
    Matrix tmp = WM.at(0) * Mat.Rows(first,last);
    for(unsigned int ctr = 1; ctr < WM.size(); ctr++){
      first = last + 1;
      last += WM.at(ctr).Ncols();
      tmp &= WM.at(ctr) * Mat.Rows(first, last);
    }
    return tmp;
  }

  void MelodicData::set_TSmode_depr()
  {
		Matrix tmp, tmpT, tmpS, tmpT2, tmpS2, tmpT3;

	    tmp = expand_dimred(mixMatrix);
	    tmpT = zeros(tmp.Nrows()/numfiles, tmp.Ncols());
	    tmpS = ones(numfiles, tmp.Ncols());

		outMsize("tmp",tmp);
		outMsize("tmpT",tmpT);
		outMsize("tmpS",tmpS);

	  	dbgmsg(string("   approach ") << opts.approach.value() << endl);	

		if(opts.approach.value()==string("tica")){
	      message("Calculating T- and S-modes " << endl);
	      explained_var = krfact(tmp,tmpT,tmpS);
	      Tmodes.clear(); Smodes.clear();
	      for(int ctr = 1; ctr <= tmp.Ncols(); ctr++){
		    tmpT3 << reshape(tmp.Column(ctr),tmpT.Nrows(),numfiles);
			outMsize("tmpT3", tmpT3);
	      	tmpT2 << tmpT.Column(ctr);
	      	tmpS2 << tmpS.Column(ctr);
			tmpT3 << SP(tmpT3,pow(ones(tmpT3.Nrows(),1)*tmpS2.t(),-1));
			if(numfiles>1)
		      tmpT2 |= tmpT3;
			if(mean(tmpS2,1).AsScalar()<0){
			  tmpT2*=-1.0;
			  tmpS2*=-1.0;
			}
	      	add_Tmodes(tmpT2);
	      	add_Smodes(tmpS2);
		  }
		}
		else{
			Tmodes.clear();
			Smodes.clear();
			for(int ctr = 1; ctr <= tmp.Ncols(); ctr++){
				tmpT3 << tmp.Column(ctr);
				add_Tmodes(tmpT3);
			}
		}

	    if(opts.approach.value()!=string("concat")){
		  //add GLM OLS fit
		  dbgmsg(string(" GLM fitting ") << endl);

		  if(Tdes.Storage()){
		    Matrix alltcs = Tmodes.at(0).Column(1);
		    for(int ctr=1; ctr < (int)Tmodes.size();ctr++)
			  alltcs|=Tmodes.at(ctr).Column(1);
		    if((alltcs.Nrows()==Tdes.Nrows())&&(Tdes.Nrows()>Tdes.Ncols()))
			  glmT.olsfit(alltcs,Tdes,Tcon);
		  }
		  if(Sdes.Storage()){
		    Matrix alltcs = Smodes.at(0);
		    for(int ctr=1; ctr < (int)Smodes.size();ctr++)
		  	  alltcs|=Smodes.at(ctr);
		    if((alltcs.Nrows()==Sdes.Nrows())&&(Sdes.Nrows()>Sdes.Ncols()&&alltcs.Nrows()>2))
			  glmS.olsfit(alltcs,Sdes,Scon);
		  }

	    }
	  //    else{
	//		dbgmsg(string(" Bypassing krfac ") << endl);
	//        add_Tmodes(tmp);
	//		add_Smodes(tmpS);
	//      }
  }

  void MelodicData::dual_regression()
  {
	dbgmsg(string("START: dual_regression") << endl);	

	Tmodes.clear();
	Smodes.clear();

	bool tmpvarnorm = opts.varnorm.value();
	opts.varnorm.set_T(false);

	Log drO;

	if(opts.dr_out.value())
		drO.makeDir(logger.appendDir("dr"),"dr.log");

	Matrix tmpData, alltcs, tmp, pinvdes;
	for(int ctr = 0; ctr < numfiles; ctr++){
		tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
		outMsize("IC",IC);
		//may want to remove the spatial means first

		
		tmp = IC*IC.t();   		outMsize("tmp",tmp);
		pinvdes = tmp.i()*IC;		  
		tmp = tmpData*pinvdes.t();      outMsize(string("tmp ")+opts.inputfname.value().at(ctr),tmp);
		if(alltcs.Storage()==0)
			alltcs=tmp;
		else
			alltcs&=tmp;
			
		// output DR
		if(opts.dr_out.value()){
			write_ascii_matrix(drO.appendDir("dr_stage1_subject"+num2str(ctr)+".txt"),tmp);
			pinvdes = tmp.t()*tmp;
			pinvdes = pinvdes.i();
			tmp = pinvdes*tmp.t()*tmpData;
			save4D(tmp,string("dr/dr_stage2_subject"+num2str(ctr)));
		}
    }

	for(int ctr = 1; ctr <= alltcs.Ncols(); ctr++){
		tmp << alltcs.Column(ctr);
		add_Tmodes(tmp);
	}

    opts.varnorm.set_T(tmpvarnorm);	
	dbgmsg(string("END: dual_regression") << endl);	
  }

  void MelodicData::set_TSmode()
  {
   	dbgmsg(string("START: set_TSmode")<< endl);	
	if(opts.dr.value())
		dual_regression();
	else
		set_TSmode_depr();
	
	dbgmsg(string("END: set_TSmode")<< endl);	
  }

  void MelodicData::setup_classic()
  {

    	Matrix alldat, tmpData;
		bool tmpvarnorm = opts.varnorm.value();

		if(numfiles > 1 && opts.joined_vn.value()){
			opts.varnorm.set_T(false);
		}
    	alldat = process_file(opts.inputfname.value().at(0), numfiles) / numfiles;
		memmsg(" after process_file ");

		if(opts.pca_dim.value() > alldat.Nrows()-2){
			cerr << "ERROR:: too many components selected \n\n";
			exit(2);
		}

 		if(opts.debug.value())
			save4D(alldat,string("preproc_dat") + num2str(1));   

		for(int ctr = 1; ctr < numfiles; ctr++){
    		tmpData = process_file(opts.inputfname.value().at(ctr), numfiles) / numfiles;
			if(opts.debug.value())
				save4D(tmpData,string("preproc_dat") + num2str(ctr+1));
			if(tmpData.Ncols() == alldat.Ncols() && tmpData.Nrows() == alldat.Nrows())
      			alldat += tmpData;	
			else{
				    if(opts.approach.value()==string("tica")){
						cerr << "ERROR:: data dimensions do not match, TICA not possible \n\n";
						exit(2); 
					}

					if(tmpData.Ncols() == alldat.Ncols()){
						int mindim = min(alldat.Nrows(),tmpData.Nrows());
						alldat = alldat.Rows(1,mindim);
						tmpData = tmpData.Rows(1,mindim);
						alldat += tmpData;	
					}				
					else 	
					message("Data dimensions do not match - ignoring "+opts.inputfname.value().at(ctr) << endl);
			}
   		}	

    	//update mask
    	if(opts.update_mask.value()){
      		message("Excluding voxels with constant value ...");
      		update_mask(Mask, alldat); 
      		message(" done" << endl);
    	}

		if((numfiles > 1 ) && opts.joined_vn.value() && tmpvarnorm){	
			//variance - normalisation
	    	message(endl<<"Normalising jointly by voxel-wise variance ..."); 
	    	stdDev = varnorm(alldat,alldat.Nrows(),3.1);
	    	stdDevi = pow(stdDev,-1); 
	    	message(" done" << endl);
		}

		if(numfiles>1)
    		message(endl << "Initial data size : "<<alldat.Nrows()<<" x "<<alldat.Ncols()<<endl<<endl);

		if(opts.debug.value())
			save4D(alldat,"alldat");
    	//estimate model order
    	Matrix tmpPPCA;
    	RowVector AdjEV, PercEV;
    	Matrix Corr, tmpE;
    	int order;

    	order = ppca_dim(remmean(alldat,2), RXweight, tmpPPCA, AdjEV, PercEV, Corr, pcaE, pcaD, Resels, opts.pca_est.value());	  
		if (opts.paradigmfname.value().length()>0)
			order += param.Ncols();

	  	if(opts.pca_dim.value() == 0){
      		opts.pca_dim.set_T(order);
			PPCA=tmpPPCA;
  		}
	  	if(opts.pca_dim.value() < 0){
      		opts.pca_dim.set_T(min(order,-1*opts.pca_dim.value()));
			PPCA=tmpPPCA;
  		}
    	order = opts.pca_dim.value();
		dbgmsg(endl << "Model order : "<<order<<endl);

		if (opts.paradigmfname.value().length()>0){
			Matrix tmpPscales;
			tmpPscales = param.t() * alldat;
			paramS = stdev(tmpPscales.t());

    		calc_white(pcaE, pcaD, order, param, paramS, whiteMatrix, dewhiteMatrix);
		}else
    		calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);

		if(opts.debug.value()){
			outMsize("pcaE",pcaE); saveascii(pcaE,"pcaE");
			outMsize("pcaD",pcaD); saveascii(pcaD,"pcaD");
			outMsize("AdjEV",AdjEV); saveascii(AdjEV,"AdjEV");
			outMsize("PercEV",PercEV); saveascii(PercEV,"PercEV");
			outMsize("tmpPPCA",tmpPPCA); saveascii(tmpPPCA,"tmpPPCA");
			outMsize("whiteMatrix",whiteMatrix); saveascii(whiteMatrix,"whiteMatrix");
			outMsize("dewhiteMatrix",dewhiteMatrix); saveascii(dewhiteMatrix,"dewhiteMatrix");
		}

		EV = AdjEV;
		EVP = PercEV;

    	if(numfiles == 1){
      		Data = alldat;
      		Matrix tmp = IdentityMatrix(Data.Nrows());
      		DWM.push_back(tmp);
      		WM.push_back(tmp);
    	} 
		else {

			dbgmsg("Multi-Subject ICA");
  			for(int ctr = 0; ctr < numfiles; ctr++){
				tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);

				if(opts.joined_vn.value() && tmpvarnorm){
					dbgmsg("tmpData normalisation"<< endl);
					dbgmsg("stdDev "  << stdDev(1,2)<< endl);
					dbgmsg("tmpData " << tmpData.SubMatrix(1,1,1,2)<< endl);
					SP3(tmpData,pow(stdDev,-1));
				}
				//  whiten (separate / joint)
				Matrix newWM,newDWM; 
				if(!opts.joined_whiten.value()){	  
					message("  Individual whitening in a " << order << " dimensional subspace " << endl);
      	  			std_pca(tmpData, RXweight, Corr, pcaE, pcaD);
	  				calc_white(pcaE, pcaD, order, newWM, newDWM);
				}else{
					if(!opts.dr_pca.value()){
						std_pca(whiteMatrix*tmpData, RXweight, Corr, pcaE, pcaD);
						calc_white(pcaE, pcaD, order, newWM, newDWM);		
						newDWM=(dewhiteMatrix*newDWM);
						newWM=(newWM*whiteMatrix);
					}
					else{
					  if(opts.debug.value())
					    dbgmsg(" --mod_pca ");
						Matrix tmp1, tmp2;
						tmp1 = whiteMatrix * alldat;
						remmean(tmp1,2);
						tmp1 *= tmpData.t();
						tmp2 = pinv(tmp1.t()).t();  
						std_pca(tmp1 * tmpData, RXweight, Corr, pcaE, pcaD);
						calc_white(pcaE, pcaD, order, newWM, newDWM);		
						newDWM=(tmp2*newDWM);
						newWM=(newWM * tmp1);
					}
				}
				DWM.push_back(newDWM);
				WM.push_back(newWM);
				tmpData = newWM * tmpData;

				//concatenate Data
				if(Data.Storage() == 0)
	  			Data = tmpData;
				else
	  			Data &= tmpData;
      	    }
        }
    	opts.varnorm.set_T(tmpvarnorm);
    
  }

  void MelodicData::setup_migp()
  {
	dbgmsg("starting MIGP");
	
	std::vector<int> myctr;
	for (int i=0; i< numfiles ; ++i) myctr.push_back(i); 

	if(opts.migp_shuffle.value()){
		message("Randomising input file order" << endl);
		std::random_shuffle ( myctr.begin(), myctr.end() );
	}
	
	Matrix tmpData;
	bool tmpvarnorm = opts.varnorm.value();

	if(numfiles > 1 && opts.joined_vn.value()){
		opts.varnorm.set_T(false);
	}
		
	for(int ctr = 0; ctr < numfiles; ctr++){
		tmpData = process_file(opts.inputfname.value().at(myctr.at(ctr)), numfiles) / numfiles;
		
		if (opts.migpN.value()==0){
			opts.migpN.set_T(2*tmpData.Nrows()-1);
		}
		if(opts.debug.value())
			save4D(tmpData,string("preproc_dat") + num2str(ctr+1));
			
		if(Data.Storage()==0)
			Data = tmpData;
		else
  			Data &= tmpData;

		outMsize("Data", Data);
		//reduce dim down to manageable level
		if(Data.Nrows() > opts.migp_factor.value()*opts.migpN.value() || ctr==numfiles-1){
			message("  Reducing data matrix to a  " << opt.migpN.value() << " dimensional subspace " << endl);
			Matrix pcaE, Corr;
			RowVector pcaD;
			std_pca(Data, RXweight, Corr, pcaE, pcaD);
		    pcaE = pcaE.Columns(pcaE.Ncols()-opts.migpN.value()+1,pcaE.Ncols());
		    Data = pcaE.t() * Data;	
		}
		outMsize("Data", Data);
		
    }

  	//update mask
    if(opts.update_mask.value()){
      message("Excluding voxels with constant value ...");
      update_mask(Mask, Data); 
      message(" done" << endl);
    }


	Matrix tmp = IdentityMatrix(Data.Nrows());
	DWM.push_back(tmp);
	WM.push_back(tmp);
   	opts.varnorm.set_T(tmpvarnorm);
  }

  void MelodicData::setup()
  { 
		
	numfiles = (int)opts.inputfname.value().size();
    setup_misc();
	if(opts.debug.value())
		memmsg(" after setup_misc ");
		
	if(opts.filtermode){ // basic setup for filtering only
		Data = process_file(opts.inputfname.value().at(0));
	}
	else{
		if((numfiles > 1) && (opts.approach.value()==string("defl") || opts.approach.value()==string("symm")))
	       opts.approach.set_T("concat");
		if(opts.migp.value())
			setup_migp();
		else
			setup_classic();   
    }

    message(endl << "  Data size : "<<Data.Nrows()<<" x "<<Data.Ncols()<<endl<<endl);
	outMsize("stdDev",stdDev);
	//meanC=mean(Data,2);
	if(opts.debug.value())
		save4D(Data,"concat_data");    
    //save the mean & mask
    save_volume(Mask,logger.appendDir("mask"));
    save_volume(Mean,logger.appendDir("mean"));
  } // void setup()
	
  void MelodicData::setup_misc()
  {

    //initialize Mean
    read_volume(Mean,opts.inputfname.value().at(0));

    //create mask
    create_mask(Mask);

	//setup background image
	if(opts.bgimage.value()>""){
		read_volume(background,opts.bgimage.value());
    	if(!samesize(Mean,background)){
        	cerr << "ERROR:: background image and data have different dimensions  \n\n";
        	exit(2);
    	}
	}else{
		background = Mean;
	}
		
    if(!samesize(Mean,Mask)){
      cerr << "ERROR:: mask and data have different dimensions  \n\n";
      exit(2);
    }

    //reset mean
    Mean *= 0;
     
    //set up weighting
    if(opts.segment.value().length()>0){
      create_RXweight();
    }

    //seed the random number generator
    double tmptime = time(NULL);
    if ( opts.seed.value() != -1 ) {
      tmptime = opts.seed.value(); 
    }
    srand((unsigned int) tmptime);

	if(opts.paradigmfname.value().length()>0){
		message("  Use columns in " << opts.paradigmfname.value() 
	      << " for PCA initialisation" <<endl);
		param = read_ascii_matrix(opts.paradigmfname.value());
		
	    Matrix tmpPU, tmpPV;
		DiagonalMatrix tmpPD;
		SVD(param, tmpPD, tmpPU, tmpPV);
		param = tmpPU;
		
		opts.pca_dim.set_T(std::max(opts.pca_dim.value(), param.Ncols()+3));		
		if(opts.debug.value()){
			outMsize("Paradigm",param); saveascii(param,"param");
		}
		//opts.guessfname.set_T(opts.paradigmfname.value());
	}

	//read in post-proc design matrices etc
	if(opts.fn_Tdesign.value().length()>0)
		Tdes = read_ascii_matrix(opts.fn_Tdesign.value());
	if(opts.fn_Sdesign.value().length()>0)
		Sdes = read_ascii_matrix(opts.fn_Sdesign.value());
	if(opts.fn_Tcon.value().length()>0)
		Tcon = read_ascii_matrix(opts.fn_Tcon.value());
	if(opts.fn_Scon.value().length()>0)
		Scon = read_ascii_matrix(opts.fn_Scon.value());
	if(opts.fn_TconF.value().length()>0)
		TconF = read_ascii_matrix(opts.fn_TconF.value());
	if(opts.fn_SconF.value().length()>0)
		SconF = read_ascii_matrix(opts.fn_SconF.value());
		
	if(numfiles>1 && Sdes.Storage() == 0){
 		Sdes = ones(numfiles,1);
		if(Scon.Storage() == 0){
			Scon = ones(1,1);
			Scon &= -1*Scon;
		}
	}
	remmean(Tdes);
  }

  void MelodicData::save()
  {   

    //check for temporal ICA
    if(opts.temporal.value()){
      message(string("temporal ICA: transform back the data ... "));
      Matrix tmpIC = mixMatrix.t();
      mixMatrix=IC.t();
      IC=tmpIC;

      unmixMatrix=pinv(mixMatrix);
      Data = Data.t();
      tmpIC = meanC;
      meanC = meanR.t();
      meanR = tmpIC.t();
      //  whiteMatrix = whiteMatrix.t;
      //  dewhiteMatrix = dewhiteMatrix.t();
      message(string("done") << endl);
      opts.temporal.set_T(false); // Do not switch again!
    }
 
    message(endl << "Writing results to : " << endl);

    //Output IC	
    if((IC.Storage()>0)&&(opts.output_origIC.value())&&(after_mm==false))
      save4D(IC,opts.outputfname.value() + "_oIC");
      
    //Output IC -- adjusted for noise	
      if(IC.Storage()>0){
				volume4D<float> tempVol;	
   
				//Matrix ICadjust;
				if(after_mm){
	  			save4D(IC,opts.outputfname.value() + "_IC");
	  			// ICadjust = IC;
				}	
				else{
					
					Matrix resids = stdev(Data - mixMatrix * IC);
					for(int ctr=1;ctr<=resids.Ncols();ctr++)
						if(resids(1,ctr) < 0.05)
							resids(1,ctr)=1;
	  //			stdNoisei = pow(stdev(Data - mixMatrix * IC)*
		//				std::sqrt((float)(Data.Nrows()-1))/
		//				std::sqrt((float)(Data.Nrows()-IC.Nrows())),-1);
	  			stdNoisei = pow(resids*
						std::sqrt((float)(Data.Nrows()-1))/
						std::sqrt((float)(Data.Nrows()-IC.Nrows())),-1);
	  	
	  			ColumnVector diagvals;
	  			diagvals=pow(diag(unmixMatrix*unmixMatrix.t()),-0.5);
	
	  			save4D(SP(IC,diagvals*stdNoisei),opts.outputfname.value() + "_IC");
				}

				if(opts.output_origIC.value())
	  			save4D(stdNoisei,string("Noise__inv"));
      }
     
    //Output T- & S-modes
 
    save_Tmodes();
    save_Smodes();

    //Output mixMatrix
    if(mixMatrix.Storage()>0){
      saveascii(expand_mix(), opts.outputfname.value() + "_mix");
      mixFFT=calc_FFT(expand_mix(), opts.logPower.value());
      saveascii(mixFFT,opts.outputfname.value() + "_FTmix");      
    }

    //Output PPCA
    if(PPCA.Storage()>0)
      saveascii(PPCA, opts.outputfname.value() + "_PPCA");
  
    //Output ICstats
    if(ICstats.Storage()>0)
      saveascii(ICstats,opts.outputfname.value() + "_ICstats"); 
      
    //Output unmixMatrix
    if(opts.output_unmix.value() && unmixMatrix.Storage()>0)
      saveascii(unmixMatrix,opts.outputfname.value() + "_unmix");

    //Output Mask
    message("  "<< logger.appendDir("mask") <<endl);

    //Output mean
    if(opts.output_mean.value() && meanC.Storage()>0 && meanR.Storage()>0){
      saveascii(meanR,opts.outputfname.value() + "_meanR");
      saveascii(meanC,opts.outputfname.value() + "_meanC");
    }

    //Output white
    if(opts.output_white.value() && whiteMatrix.Storage()>0&&
      dewhiteMatrix.Storage()>0){
      	saveascii(whiteMatrix,opts.outputfname.value() + "_white");
      	saveascii(dewhiteMatrix,opts.outputfname.value() + "_dewhite");
      	Matrix tmp;
      	tmp=calc_FFT(dewhiteMatrix, opts.logPower.value());
      	saveascii(tmp,opts.outputfname.value() + "_FTdewhite");
	  }

    //Output PCA
    if(opts.output_pca.value() && pcaD.Storage()>0&&pcaE.Storage()>0){
      saveascii(pcaE,opts.outputfname.value() + "_pcaE");
      saveascii((Matrix) diag(pcaD),opts.outputfname.value() + "_pcaD");
      
      Matrix PCAmaps;
      if(whiteMatrix.Ncols()==Data.Ncols())
				PCAmaps = dewhiteMatrix.t();
      else
				PCAmaps = whiteMatrix * Data;

      save4D(PCAmaps,opts.outputfname.value() + "_pca");
    }
 
		message("...done" << endl);
  } //void save()

  int MelodicData::remove_components()
  {  
    message("Reading " << opts.filtermix.value() << endl) 
    mixMatrix = read_ascii_matrix(opts.filtermix.value());
    if (mixMatrix.Storage()<=0) {
      cerr <<" Please specify the mixing matrix correctly" << endl;
      exit(2);
    }
    
    unmixMatrix = pinv(mixMatrix);
    IC = unmixMatrix * Data;

    string tmpstr;
    tmpstr = opts.filter.value();

    Matrix noiseMix;
    Matrix noiseIC;

    int ctr=0;    
    char *p;
    char t[1024];
    const char *discard = ", [];{(})abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ~!@#$%^&*_-=+|\':><./?";
    
    message("Filtering the data...");
    strcpy(t, tmpstr.c_str());
    p=strtok(t,discard);
    ctr = atoi(p);
    if(ctr>0 && ctr<=mixMatrix.Ncols()){
      message(" "<< ctr );
      noiseMix = mixMatrix.Column(ctr);
      noiseIC  = IC.Row(ctr).t();    
    }
    else{
      cerr << endl<< "component number "<<ctr<<" does not exist" << endl;
    }
    
    do{
      p=strtok(NULL,discard);
      if(p){
				ctr = atoi(p);
	
        if(ctr>0 && ctr<=mixMatrix.Ncols()){
	  			message(" "<<ctr);
	  			noiseMix |= mixMatrix.Column(ctr);
	  			noiseIC  |= IC.Row(ctr).t();
				}
				else{
	  			cerr << endl<< "component number "<<ctr<<" does not exist" << endl;
				}
      }
    }while(p);
    message(endl);
    Matrix newData;

		outMsize("DATA",Data);
		outMsize("IC",IC);
		outMsize("noiseIC",noiseIC);
		outMsize("noiseMix",noiseMix);
		outMsize("meanR",meanR);
		outMsize("meanC",meanC);

    newData = Data - noiseMix * noiseIC.t();

		if(meanR.Storage()>0)
    	newData = newData + ones(newData.Nrows(),1)*meanR;
    
    volume4D<float> tmp;
    read_volume4D(tmp,opts.inputfname.value().at(0)); 
    tmp.setmatrix(newData,Mask);
    save_volume4D(tmp,logger.appendDir(opts.outputfname.value() + "_ICAfiltered")); 
   
    return 0;
  } // int remove_components()

  void MelodicData::create_RXweight()
  {
    message("Reading the weights for the covariance R_X from file "<< opts.segment.value() << endl);
  
    volume4D<float> tmpRX;
    read_volume4D(tmpRX,opts.segment.value());
    RXweight = tmpRX.matrix(Mask);
  } 

  void MelodicData::est_smoothness()
  {
    if(Resels == 0){
      string SM_path = opts.binpath + "smoothest";
      string Mask_fname = logger.appendDir("mask");

      if(opts.segment.value().length()>0){
				Mask_fname =  opts.segment.value();
      } 

      // Setup external call to smoothest:
      char callSMOOTHESTstr[1000];
      ostrstream osc(callSMOOTHESTstr,1000);
      osc  << SM_path << " -d " << data_dim()
	   		<< " -r " << opts.inputfname.value().at(0) << " -m " 
	   		<< Mask_fname << " > " << logger.appendDir("smoothest") << '\0';
      
      message("  Calling Smoothest: " << callSMOOTHESTstr << endl);
      system(callSMOOTHESTstr);

      //read back the results
      ifstream in;
      string str;
      Resels = 1.0;
      
      in.open(logger.appendDir("smoothest").c_str(), ios::in);
      if(in>0){
				for(int ctr=1; ctr<7; ctr++)
					in >> str;
				in.close();
				if(str!="nan")
	  			Resels = atof(str.c_str());
      }
    }
  }

  unsigned long MelodicData::standardise(volume<float>& mask, volume4D<float>& R)
  {
    	unsigned long count = 0;
    	int M=R.tsize();
    
    	for (int z=mask.minz(); z<=mask.maxz(); z++) {
      	for (int y=mask.miny(); y<=mask.maxy(); y++) {
					for (int x=mask.minx(); x<=mask.maxx(); x++) {
	  				if( mask(x,y,z) > 0.5) {
	    				count ++;
	    				if( M > 2 ) {
	      				// For each voxel 
	      				//    calculate mean and standard deviation...
	      				double Sx = 0.0, SSx = 0.0;	      
	      				for ( int t = 0; t < M; t++ ) {
									float R_it = R(x,y,z,t);
									Sx += R_it;
									SSx += (R_it)*(R_it);
	      				}
	      
	      				float mean = Sx / M;
	      				float sdsq = (SSx - ((Sx)*(Sx) / M)) / (M - 1) ;
	      
	      				if (sdsq<=0) {
									// trap for differences between mask and invalid data
									mask(x,y,z)=0;
									count--;
	      				} else {
	      					//    ... and use them to standardise to N(0, 1).
									for ( unsigned short t = 0; t < M; t++ ) {
		  							R(x,y,z,t) = (R(x,y,z,t) - mean) / sqrt(sdsq);
									}
	      				} 
	    				}
	  				}
					}
      	}
    	}  
    	return count;
  }

  float MelodicData::est_resels(volume4D<float> R, volume<float> mask)
  {
    message("  Estimating data smoothness ... ");
    unsigned long mask_volume = standardise(mask, R);

    int dof = R.tsize();
    unsigned long N = mask_volume;

    // MJ additions to make it cope with 2D images
    bool usez = true;
    if (R.zsize() <= 1) { usez = false; }

    enum {X = 0, Y, Z};
    float SSminus[3] = {0, 0, 0}, S2[3] = {0, 0, 0};

    int zstart=1;
    if (!usez) zstart=0;
    for ( unsigned short z = zstart; z < R.zsize() ; z++ )
      for ( unsigned short y = 1; y < R.ysize() ; y++ )
				for ( unsigned short x = 1; x < R.xsize() ; x++ )
	  			// Sum over N
	  			if( (mask(x, y, z)>0.5) &&
	      		(mask(x-1, y, z)>0.5) && 
	      		(mask(x, y-1, z)>0.5) && 
	      		( (!usez) || (mask(x, y, z-1)>0.5) ) ) {
	    
	    				N++;
	  
	    				for ( unsigned short t = 0; t < R.tsize(); t++ ) {
	      				// Sum over M
	      				SSminus[X] += R(x, y, z, t) * R(x-1, y, z, t);
	      				SSminus[Y] += R(x, y, z, t) * R(x, y-1, z, t);
	      				if (usez) SSminus[Z] += R(x, y, z, t) * R(x, y, z-1, t);

	      				S2[X] += 0.5 * (R(x, y, z, t)*R(x, y, z, t) + 
									R(x-1, y, z, t)*R(x-1, y, z, t));
	      				S2[Y] += 0.5 * (R(x, y, z, t)*R(x, y, z, t) + 
									R(x, y-1, z, t)*R(x, y-1, z, t));
	      				if (usez) S2[Z] += 0.5 * (R(x, y, z, t)*R(x, y, z, t) + 
									R(x, y, z-1, t)*R(x, y, z-1, t));
	    				}
	  			}

    float norm = 1.0/(float) N;
    float v = dof;	// v - degrees of freedom (nu)  
    if(R.tsize() > 1) {
      norm = (v - 2) / ((v - 1) * N * R.tsize());
    }
 
    // for extreme smoothness 
    if (SSminus[X]>=0.99999999*S2[X]) 
      SSminus[X]=0.99999*S2[X];  
    if (SSminus[Y]>=0.99999999*S2[Y]) 
      SSminus[Y]=0.99999*S2[Y];
    if (usez) 
      if (SSminus[Z]>=0.99999999*S2[Z]) 
				SSminus[Z]=0.99999*S2[Z];
    // Convert to sigma squared
    float sigmasq[3] = {0,0,0};
    sigmasq[X] = -1.0 / (4 * log(fabs(SSminus[X]/S2[X])));
    sigmasq[Y] = -1.0 / (4 * log(fabs(SSminus[Y]/S2[Y])));
    if (usez) { sigmasq[Z] = -1.0 / (4 * log(fabs(SSminus[Z]/S2[Z]))); }
    
    // Convert to full width half maximum
    float FWHM[3] = {0,0,0};
    FWHM[X] = sqrt(8 * log(2) * sigmasq[X]);
    FWHM[Y] = sqrt(8 * log(2) * sigmasq[Y]);
    if (usez) { FWHM[Z] = sqrt(8 * log(2) * sigmasq[Z]); }
    float resels = FWHM[X] * FWHM[Y];
    if (usez) resels *= FWHM[Z];

    message(" done " <<endl);
    return resels;
  }

  void MelodicData::create_mask(volume<float>& theMask)
  {
    if(opts.use_mask.value() && opts.maskfname.value().size()>0){   // mask provided 
      read_volume(theMask,opts.maskfname.value());
      message("Mask provided : " << opts.maskfname.value()<<endl<<endl);
    }
    else{
      if(opts.perf_bet.value() && opts.use_mask.value()){ //use BET
				message("Create mask ... ");
    		//save first image
    		tmpnam(Mean_fname); // generate a tmp name
    		save_volume(Mean,Mean_fname);    

				// set up all strings
				string BET_outputfname = string(Mean_fname)+"_brain";

				string BET_path = opts.binpath + "bet";
				string BET_optarg = "-m -f 0.4"; // see man bet
				string Mask_fname = BET_outputfname+"_mask";

				// Setup external call to BET:

		//		char callBETstr[1000];
	//			ostrstream betosc(callBETstr,1000);
//				betosc  << BET_path << " " << Mean_fname << " " 
//	     		<< BET_outputfname << " " << BET_optarg << " > /dev/null " << '\0';