melgmix.cc

/*  MELODIC - Multivariate exploratory linear optimized decomposition into 
              independent components
    
    melgmix.cc - Gaussian Mixture Model

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

/*  CCOPYRIGHT */

#include "newimage/newimageall.h"
//#include "melmmopts.h"
#include "melgmix.h"
//#include "melmm.h"
#include "utils/log.h"
#include "miscmaths/miscprob.h"
#include <time.h>
#include "libvis/miscplot.h"
#include "libvis/miscpic.h"

using namespace Utilities;
using namespace NEWIMAGE;
using namespace MISCPLOT;
using namespace MISCPIC;

string float2str(float f,int width, int prec, bool scientif){
 	ostringstream os;
  int redw = int(std::abs(std::log10(std::abs(f))))+1;
  if(width>0)
    os.width(width);
  if(scientif)
    os.setf(ios::scientific);
  os.precision(redw+std::abs(prec));
  os.setf(ios::internal, ios::adjustfield);
  os << f;
  return os.str();
} 
  
namespace Melodic{
  
  void MelGMix::setup(const RowVector& dat, 
		const string dirname,
		int cnum, volume<float> themask, 
		volume<float> themean, 
		int num_mix, float eps, bool fixit){
    	cnumber = cnum;
    	Mask = themask;
    	Mean = themean;
    	prefix = string("IC_")+num2str(cnum);
    	
    	fitted = false;
    	nummix = num_mix;
    	numdata = dat.Ncols();
    
    	//normalise data 
    	datamean = mean(dat,2).AsScalar();
    	datastdev= stdev(dat,2).AsScalar();
    	data=(dat - datamean)/datastdev;

			dbgmsg(" mapdat; mean: " << datamean << " std: " <<datastdev << endl);

    	props=zeros(1,nummix);
    	vars=zeros(1,nummix);
    	means=zeros(1,nummix);
    	Params=zeros(1,nummix);
    	logprobY = 1.0;

    	props = std::pow(float(nummix),float(-1.0));