clust.cpp

高斯混合模型的C语言实现

/*
* All questions regarding the software should be addressed to
* 
*       Prof. Charles A. Bouman
*       Purdue University
*       School of Electrical and Computer Engineering
*       1285 Electrical Engineering Building
*       West Lafayette, IN 47907-1285
*       USA
*       [phone]
*       [phone] (fax)
*       email:  [e-mail]
*       http://www.ece.purdue.edu/~bouman
* 
* Copyright (c) 1995 The Board of Trustees of Purdue University.
*
* Permission to use, copy, modify, and distribute this software and its
* documentation for any purpose, without fee, and without written agreement is
* hereby granted, provided that the above copyright notice and the following
* two paragraphs appear in all copies of this software.
*
* IN NO EVENT SHALL PURDUE UNIVERSITY BE LIABLE TO ANY PARTY FOR DIRECT,
* INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
* USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF PURDUE UNIVERSITY HAS
* BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
* PURDUE UNIVERSITY SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
* PARTICULAR PURPOSE.  THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS,
* AND PURDUE UNIVERSITY HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT,
* UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
*/

/************************************************************************
* Modified for the purpose of image and video segmentation
*
*	Dr. Liu Zhi
*	School of Communication and Information Engineering
*	Shanghai University, 200072
*	P.R.China
*
************************************************************************/
#include "stdafx.h"

#include "clust.h"
#include "classify_util.h"

int GMM_Estimate(struct SigSet* S, double* sample_vector, int num_of_samples, int vector_dimension, int nclasses, int init_num_of_subclasses, int desired_num_of_subclasses, int is_cov_diag)
{
    int i,j,k;
    double Rmin;
	int max_num;
	int vector_feature_cnt;
	int flag_cov;
    struct SigSet::ClassSig *Sig;

    /* Initialize SigSet data structure */
    I_InitSigSet (S);
    I_SigSetNBands (S, vector_dimension);
    I_SetSigTitle (S, "test signature set");

    /* Allocate memory for cluster signatures */
    for(k=0; k<nclasses; k++) 
	{
      Sig = (struct SigSet::ClassSig *)I_NewClassSig(S);
      I_SetClassTitle ((struct ClassSig *)Sig, "test class signature");
      for(i=0; i<init_num_of_subclasses; i++)
        I_NewSubSig (S, (struct ClassSig *)Sig);
    }

    /* Read data for each class */
    for(k=0; k<nclasses; k++) 
	{
      Sig = &(S->ClassSig[k]);
      I_AllocClassData (S, (struct ClassSig *)Sig, num_of_samples);

	  vector_feature_cnt = 0;
      for(i=0; i<Sig->ClassData.npixels; i++) 
	  {
        for(j=0; j<vector_dimension; j++) 
		{
		  Sig->ClassData.x[i][j] = sample_vector[vector_feature_cnt++];
        }
      }

      /* Set unity weights and compute SummedWeights */
      Sig->ClassData.SummedWeights = 0.0;
      for(i=0; i<Sig->ClassData.npixels; i++) 
	  {
        Sig->ClassData.w[i] = 1.0;
        Sig->ClassData.SummedWeights += Sig->ClassData.w[i];
      }
    }

    /* Compute the average variance over all classes */
    Rmin = 0;
    for(k=0; k<nclasses; k++) 
	{
      Sig = &(S->ClassSig[k]);
      Rmin += AverageVariance((struct SigSet::ClassSig *)Sig, vector_dimension);
    }
    Rmin = Rmin/(COVAR_DYNAMIC_RANGE*nclasses);

    /* Perform clustering for each class */
    for(k=0; k<nclasses; k++) 
	{
      Sig = &(S->ClassSig[k]);

      fprintf(stdout,"Start clustering class %d\n\n",k);

	  /* Determine number of gaussians automatically? */
	  if(desired_num_of_subclasses <= 0)
		  desired_num_of_subclasses = 0;

	  /* Determine covariance matrix is diagonal? */
	  if(is_cov_diag == 1)
		  flag_cov = (int)CLUSTER_DIAG;
	  else
		  flag_cov = (int)CLUSTER_FULL;

	  /* Perform EM algorithm for this GMM model */
      subcluster(S, k, desired_num_of_subclasses, flag_cov, Rmin, &max_num);
		  
      fprintf(stdout,"Maximum number of subclasses = %d\n",max_num);

    }

   	//Get probability for each subclass
	GMM_Subclass_Prob(S);

    return 0;
}

double AverageVariance(struct SigSet::ClassSig * Sig, int nbands)
{
     int     i,b1;
     double  *mean,**R,Rmin;

     /* Compute the mean of variance for each band */
     mean = G_alloc_vector(nbands);
     R = G_alloc_matrix(nbands,nbands);

     for(b1=0; b1<nbands; b1++) {
       mean[b1] = 0.0;
       for(i=0; i<Sig->ClassData.npixels; i++) {
         mean[b1] += (Sig->ClassData.x[i][b1])*(Sig->ClassData.w[i]);
       }
       mean[b1] /= Sig->ClassData.SummedWeights;
     }

     for(b1=0; b1<nbands; b1++) {
       R[b1][b1] = 0.0;
       for(i=0; i<Sig->ClassData.npixels; i++) {
         R[b1][b1] += (Sig->ClassData.x[i][b1])*(Sig->ClassData.x[i][b1])*(Sig->ClassData.w[i]);
       }
       R[b1][b1] /= Sig->ClassData.SummedWeights;
       R[b1][b1] -= mean[b1]*mean[b1];
     }

     /* Compute average of diagonal entries */
     Rmin = 0.0;
     for(b1=0; b1<nbands; b1++) 
       Rmin += R[b1][b1];

     Rmin = Rmin/(nbands);

     G_free_vector(mean);
     G_free_matrix(R);

     return(Rmin);
}