gaussian.cpp

混合高斯模型的C++程序

#include "Gaussian.h"

#include "gmm/alloc_util.h"
#include "gmm/clust_io.h"
#include "gmm/clust_util.h"
#include "gmm/subcluster.h"
#include "gmm/classify_util.h"

CGaussian::CGaussian(void)
{
	Initialize();
}


CGaussian::~CGaussian(void)
{
	Clear();
}


// 从文件中导入GMM模型，path为文件目录，返回是否装载成功
bool CGaussian::LoadModel(const char *path)
{
	Clear();
	FILE *fp = fopen(path,"r");
	if(fp == NULL) {
		fprintf(stderr, "\nError: Can't open parameter file %s", path); 
		return false;
	}
	I_ReadSigSet(fp,&m_model);
	fclose(fp);
	return true;
}


// 将GMM模型导出文件，path为文件目录，返回是否保存成功
bool CGaussian::SaveModel(const char *path)
{
	FILE *fp = fopen(path,"w");
	if(fp == NULL) {
		fprintf(stderr,"can't open parameter output file");
		return false;
	}
	I_WriteSigSet(fp, &m_model);
	fclose(fp);
	return true;
}


// 设置GMM的高斯数目，clusterClass为聚类后的数目，若值为0则表示自动确定，maxClass为初始化时的最大类别数
void CGaussian::SetClass(int clusterClass, int maxClass)
{
	m_clusterClass = clusterClass;
	m_maxClass = maxClass;
}


// 设置协方差矩阵的类型，isFull为false表示为对角矩阵，true表示为全矩阵
void CGaussian::SetCovarianceType(bool isFull)
{
	m_isFull = isFull;
}


// 设置GMM的数据向量的维数，dimension为维数
void CGaussian::SetDimension(int dimension)
{
	m_dimension = dimension;
}


// 根据给定的数据训练GMM，pSample为数据块的指针，number为数据的数目
void CGaussian::TrainSample(double **pSample, int number)
{
	I_SigSetNBands(&m_model, m_dimension);
	I_SetSigTitle(&m_model, "GMM model");

	ClassSig *sig = I_NewClassSig(&m_model);
	I_SetClassTitle(sig, "GMM class");
	for(int index = 0; index < m_maxClass; ++ index)
		I_NewSubSig(&m_model, sig);

	I_AllocClassData(&m_model, sig, number);
	for(int n = 0; n < number; ++ n) {
		for(int d = 0; d < m_dimension; ++ d) {
			sig->classData.x[n][d] = pSample[n][d];
		}
	}

	sig->classData.SummedWeights = 0.0;
	for(int n = 0; n < sig->classData.npixels; ++ n) {
		sig->classData.w[n] = 1.0;
		sig->classData.SummedWeights += sig->classData.w[n];
	}

	double *mean = G_alloc_vector(m_dimension);
	double **covariance = G_alloc_matrix(m_dimension, m_dimension);
	for(int d = 0; d < m_dimension; ++ d) {
		mean[d] = 0;
		for(int n = 0; n < sig->classData.npixels; ++ n) {
			mean[d] += (sig->classData.x[n][d]) * (sig->classData.w[n]);
		}
		mean[d] /= sig->classData.SummedWeights;
	}
	for(int d = 0; d < m_dimension; ++ d) {
		covariance[d][d] = 0;
		for(int n = 0; n < sig->classData.npixels; ++ n) {
			covariance[d][d] += (sig->classData.x[n][d]) * (sig->classData.x[n][d]) * (sig->classData.w[n]);
		}
		covariance[d][d] /= sig->classData.SummedWeights;
		covariance[d][d] -= mean[d] * mean[d];
	}

	double rmin = 0.0;
	for(int d = 0; d < m_dimension; ++ d)
		rmin += covariance[d][d];
	rmin = rmin / m_dimension;

	G_free_vector(mean);
	G_free_matrix(covariance);

	rmin /= COVAR_DYNAMIC_RANGE;
	int maxNumber = 0;
	if(m_isFull)
		subcluster(&m_model, 0, m_clusterClass,(int)CLUSTER_FULL, rmin, &maxNumber);
	else
		subcluster(&m_model, 0, m_clusterClass,(int)CLUSTER_DIAG, rmin, &maxNumber);
	I_DeallocClassData(&m_model, sig);

	m_clusterClass = sig->nsubclasses;
	m_classInitialize = true;
}


// 获得训练后的GMM的高斯维数
int CGaussian::GetClass()
{
	return m_clusterClass;
}


// 获得数据向量的维数
int CGaussian::GetDimension()
{
	return m_dimension;
}


// 获得数据出现的概率的log值
double CGaussian::GetLogProbability(double *pSample)
{
	if(m_classInitialize) {
		ClassLogLikelihood_init(&m_model);
		m_classInitialize = false;
	}
	double *logPropability = G_alloc_vector(m_model.nclasses);
	ClassLogLikelihood(pSample, logPropability, &m_model);
	double answer = logPropability[0];
	G_free_vector(logPropability);
	return answer;
}


// 获得各维高斯的权重
double CGaussian::GetPi(int index)
{
	return m_model.classSig[0].subSig[index].pi;
}


// 获得各维高斯的均值
double *CGaussian::GetMean(int index)
{
	return m_model.classSig[0].subSig[index].means;
}


// 获得各维高斯的协方差矩阵
double **CGaussian::GetCovariance(int index)
{
	return m_model.classSig[0].subSig[index].R;
}


// 设置各维高斯的权重
void CGaussian::SetPi(int index, double pi)
{
	m_model.classSig[0].subSig[index].pi = pi;
}


// 设置各维高斯的均值
void CGaussian::SetMean(int index, double *mean)
{
	double *curMean = GetMean(index);
	memcpy(curMean, mean, m_dimension * sizeof(double));
}


// 设置各维高斯的协方差矩阵
void CGaussian::SetCovariance(int index, double **covariance)
{
	double **curCovariance = GetCovariance(index);
	for(int d = 0; d < m_dimension; ++ d)
		memcpy(curCovariance[d], covariance[d], m_dimension * sizeof(double));
	m_classInitialize = true;
}


// 初始化
void CGaussian::Initialize()
{
	I_InitSigSet(&m_model);
	m_clusterClass = 0, m_maxClass = 20;
	m_isFull = true;
	m_dimension = 1;
	m_classInitialize = true;;
}


// 清空
void CGaussian::Clear()
{
	I_DeallocSigSet(&m_model);
}