melgmix.h

/*  MELODIC - Multivariate exploratory linear optimized decomposition into 
              independent components
    
    melgmix.h - class for Gaussian/Gamma Mixture Model

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

/*  CCOPYRIGHT */

#ifndef __MELGMIX_h
#define __MELGMIX_h

#include "newimage/newimageall.h"
#include "utils/log.h"
#include "melodic.h"
#include "utils/options.h"
#include "meloptions.h"
//#include "melreport.h"

using namespace Utilities;
using namespace NEWIMAGE;

namespace Melodic{
  
  class MelGMix{
    public:
 
      MelGMix(MelodicOptions &popts, Log &plogger):
				opts(popts),
				logger(plogger){}

      ~MelGMix() { 
				//mainhtml << endl << "<hr></CENTER></BODY></HTML>" << endl;
	  	}

      void save();

      void setup(const RowVector& dat, volumeinfo inf, const string dirname,
		 		int here, volume<float> themask, 
		 		volume<float> themean, int num_mix = 3, 
		 		float eps = 0.0, bool fixdim = false);
      
      void gmmfit();
      void ggmfit();

      inline void fit(string mtype = string("GGM")){
	  		mmtype = mtype;
	  		if(mmtype==string("GGM")) 
	    		this->ggmfit(); 
	  		else 
	    		this->gmmfit();

	  		//re-insert mean and stdev
	  		data = data*datastdev + datamean;
	  		//threshmaps = threshmaps*datastdev + datamean;
	  		means = means*datastdev + datamean;
	  		vars = vars*datastdev*datastdev;
			}

      inline Matrix threshold(string levels){
				return this->threshold(data, levels);
			}
      inline Matrix threshold(RowVector& levels){ 
				return this->threshold(data, levels);
			}
      Matrix threshold(const RowVector& dat, Matrix& levels);
      Matrix threshold(const RowVector& dat, string levels);

      void status(const string &txt);
 
      inline RowVector& get_means() {return means;}
      inline void set_means(RowVector& Arg) {means = Arg;}
    
      inline RowVector& get_vars() {return vars;}
      inline void set_vars(RowVector& Arg) {vars = Arg;}
      
      inline RowVector& get_pi() {return props;}
      inline void set_pi(RowVector& Arg) {props = Arg;}
      
      inline RowVector& get_data() {return data;}
      inline void set_data(RowVector& Arg) {data = Arg;}

      inline RowVector& get_prob() {return probmap;}

      inline float get_eps() {return epsilon;}
      inline void set_eps(float Arg) {epsilon = Arg;}

      inline Matrix& get_threshmaps() {return threshmaps;}
      inline void set_threshmaps(Matrix& Arg) {threshmaps = Arg;}

      inline bool isfitted(){return fitted;}

      inline int mixtures(){return nummix;}
     
      inline string get_type() { return mmtype;}
      inline void set_type(string Arg) { mmtype = Arg;}
      
      inline string get_prefix() { return prefix;}
      inline void  set_prefix(string Arg) { prefix = Arg;}

      inline RowVector get_probmap() {return probmap;}

      inline float get_offset() {return offset;}
      inline void set_offset(float Arg) {offset = Arg;}

      inline void flipres(int num){
				means = -means;
				data = -data;
				threshmaps = -threshmaps;
				if(mmtype=="GGM"){
	  			float tmp;
	  			tmp= means(2);means(2)=means(3);means(3)=tmp;
	  			tmp=vars(2);vars(2)=vars(3);vars(3)=tmp;
	  			tmp=props(2);props(2)=props(3);props(3)=tmp;
				}
      }      

      void create_rep();

      inline void add_infstr(string what){
				threshinfo.push_back(what);
      }

      inline string get_infstr(int num){
				if((threshinfo.size()<(unsigned int)(num-1))||(num<1))
	  			return string("");
				else
	  			return threshinfo[num-1];
      }

      inline int size_infstr(){
      	return threshinfo.size();
      }

      inline void clear_infstr(){
				threshinfo.clear();
      }

      inline void smooth_probs(float howmuch){
				volume4D<float> tempVol;
				tempVol.setmatrix(probmap,Mask);
        tempVol[0]= smooth(tempVol[0],howmuch);
        probmap = tempVol.matrix(Mask);
      }

      double datamean;
      double datastdev;

    private:
      MelodicOptions &opts;     
      Log &logger; //global log file

      //Log mainhtml;

      void gmmupdate();
      float gmmevidence();
      void gmmreducemm();
      void add_params(Matrix& mu, Matrix& sig, Matrix& pi, 
		    float logLH, float MDL, float Evi, bool advance = false);
      void get_params(int index, Matrix& mu, Matrix& sig, Matrix& pi, 
		    float logLH, float MDL, float Evi);

      Matrix Params;
      Matrix threshmaps;

      RowVector means;
      RowVector vars;
      RowVector props;
      RowVector data;
      RowVector probmap;

      volume<float> Mean;
      volume<float> Mask;

      float epsilon;
      float logprobY;
      float MDL;
      float Evi;
      float offset;

      int nummix;
      int numdata;
      int cnumber;

      bool fitted;
      bool fixdim;

      string prefix;
      string mmtype;
      string dirname;

      volumeinfo bginfo;
      vector<string> threshinfo;

  };
}

#endif