subcluster.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.
*/

#include "stdafx.h"

#include "clust_defs.h"
#include "subcluster.h"
#include "alloc_util.h"
#include "clust_util.h"


int clusterMessageVerboseLevel;

static void seed(struct ClassSig *Sig, int nbands, double Rmin, int option);

static double refine_clusters(
    struct ClassSig *Sig, 
    int nbands, 
    double Rmin, 
    int option);

static void reestimate(struct ClassSig *Sig, int nbands, double Rmin, int option);

static double regroup(struct ClassSig *Sig, int nbands);

static void reduce_order(
    struct ClassSig *Sig,
    int nbands,
    int *min_ii,
    int *min_jj);

static double loglike(
    double *x, 
    struct SubSig *SubSig, 
    int nbands);

static double distance(
    struct SubSig *SubSig1,
    struct SubSig *SubSig2,
    int nbands);

static void compute_constants(struct ClassSig *Sig, int nbands);

static void normalize_pi(struct ClassSig *Sig);

static void add_SubSigs(
    struct SubSig *SubSig1,
    struct SubSig *SubSig2,
    struct SubSig *SubSig3,
    int nbands);

static void save_ClassSig(
    struct ClassSig *Sig1,
    struct SigSet *S,
    int nbands);

static void copy_ClassSig(
    struct ClassSig *Sig1,
    struct ClassSig *Sig2,
    int nbands);

static void copy_SubSig(
    struct SubSig *SubSig1,
    struct SubSig *SubSig2,
    int nbands);

static void DiagonalizeMatrix(double **R, int nbands);

#if 0
static void list_Sig(struct ClassSig *Sig, int nbands);
static void print_class(struct ClassSig *Sig, char *fname);
#endif


int subcluster(
    struct SigSet *S, /* Input: structure contataining input data */ 
    int Class_Index,  /* Input: index corresponding to class to be processed */
    int desired_num,  /* Input: desired number of subclusters. */
                      /*      0=>ignore this input. */
    int option,       /* Input: type of clustering to use */
                      /*      option=1=CLUSTER_FULL=>full covariance matrix */
                      /*      option=0=CLUSTER_DIAG=>diagonal covariance matrix */
    double Rmin,      /* Minimum value for diagonal elements of convariance */
    int *Max_num)     /* Output: maximum number of allowed subclusters */
{
    int nparams_clust;
    int ndata_points;
    int min_i,min_j;
    int status;
    int nbands;
    double rissanen;
    double min_riss;
    struct ClassSig *Sig;
    static struct SigSet Smin;

    status = 0;

    /* set class pointer */
    Sig = (struct ClassSig *)&(S->ClassSig[Class_Index]);

    /* set number of bands */
    nbands = S->nbands;

    /* compute number of parameters per cluster */
    nparams_clust = 1+nbands+0.5*(nbands+1)*nbands;
    if(option==CLUSTER_DIAG) nparams_clust = 1+nbands+nbands;

    /* compute number of data points */
    ndata_points = Sig->ClassData.npixels*nbands;

    /* compute maximum number of subclasses */
    ndata_points = Sig->ClassData.npixels*nbands;
    *Max_num = (ndata_points+1)/nparams_clust - 1;

    /* check for too many subclasses */
    if(Sig->nsubclasses > (*Max_num/2) )
    {
      Sig->nsubclasses = *Max_num/2;
      fprintf(stderr,"Too many subclasses for class index %d\n",Class_Index);
      fprintf(stderr,"         number of subclasses set to %d\n\n",Sig->nsubclasses);
      status = -2;
    }

    /* initialize clustering */
    seed(Sig,nbands,Rmin,option);

    /* EM algorithm */
    min_riss = refine_clusters(Sig,nbands,Rmin,option);

    if(2<=clusterMessageVerboseLevel) {
      fprintf(stdout,"Subclasses = %d; Rissanen = %f; \n",Sig->nsubclasses,min_riss);
    }

    /* Save contents of Class Signature to Smin */
    save_ClassSig(Sig,&Smin,nbands);

    if(desired_num==0) {
      while(Sig->nsubclasses>1) {
        reduce_order(Sig,nbands,&min_i,&min_j);

        if(2<=clusterMessageVerboseLevel) {
          fprintf(stdout,"Combining Subclasses (%d,%d)\n",min_i,min_j);
        }

        rissanen = refine_clusters(Sig,nbands,Rmin,option);

        if(2<=clusterMessageVerboseLevel) {
          fprintf(stdout,"Subclasses = %d; Rissanen = %f; \n",
                       Sig->nsubclasses, rissanen);
        }

        if(rissanen<min_riss)
        {
          min_riss = rissanen;

          /* Delete old Smin, and save new Smin */
          I_DeallocSigSet(&Smin);
          save_ClassSig(Sig,&Smin,nbands);
        }
      }
    }
    else {
      while( (Sig->nsubclasses>desired_num)&&(Sig->nsubclasses>0) ) {
        reduce_order(Sig,nbands,&min_i,&min_j);

        if(2<=clusterMessageVerboseLevel) {
          fprintf(stdout,"Combining Subclasses (%d,%d)\n",min_i,min_j);
        }
 
        rissanen = refine_clusters(Sig,nbands,Rmin,option);

        if(2<=clusterMessageVerboseLevel) {
          fprintf(stdout,"Subclasses = %d; Rissanen = %f; \n",
                       Sig->nsubclasses,rissanen);
        }

        /* Delete old Smin, and save new Smin */
        I_DeallocSigSet(&Smin);
        save_ClassSig(Sig,&Smin,nbands);
      }  
    }

    /* Deallocate memory for class, and replace with solution */
    while(Sig->nsubclasses>0) I_DeallocSubSig(Sig);
    Sig->SubSig = (struct ClassSig::SubSig *)Smin.ClassSig[0].SubSig;
    Sig->nsubclasses = Smin.ClassSig[0].nsubclasses;

    /* return warning status */
    return(status);
}


/******************************************************************/
/* Computes initial values for parameters of Gaussian Mixture     */
/* model. The subroutine returns the minimum allowed value for    */
/* the diagonal entries of the convariance matrix of each class.  */
/*****************************************************************/
static void seed(struct ClassSig *Sig, int nbands, double Rmin, int option)
{
     int     i,b1,b2;
     double  period;
     double  *mean,**R;

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

     /* If diagonal clustering is desired, then diagonalize matrix */
     if(option==CLUSTER_DIAG) DiagonalizeMatrix(R,nbands);

     /* Compute the sampling period for seeding */
     if(Sig->nsubclasses>1) {
       period = (Sig->ClassData.npixels-1)/(Sig->nsubclasses-1.0);
     }
     else period =0;


     /* Seed the means and set the covarience components */
     for(i=0; i<Sig->nsubclasses; i++) {
       for(b1=0; b1<nbands; b1++) {
         Sig->SubSig[i].means[b1] = Sig->ClassData.x[(int)(i*period)][b1];
       }

       for(b1=0; b1<nbands; b1++)
       for(b2=0; b2<nbands; b2++) {
         Sig->SubSig[i].R[b1][b2] = R[b1][b2];
       }
       for(b1=0; b1<nbands; b1++) {
         Sig->SubSig[i].R[b1][b1] += Rmin;
       }
       Sig->SubSig[i].pi = 1.0/Sig->nsubclasses;
     }

     G_free_vector(mean);
     G_free_matrix(R);

     compute_constants(Sig,nbands);
     normalize_pi(Sig);
}


/*****************************************************************/
/* Computes ML clustering of data using Gaussian Mixture model.  */
/* Returns the values of the Rissen constant for the clustering. */
/*****************************************************************/
static double refine_clusters(
    struct ClassSig *Sig, 
    int nbands, 
    double Rmin, 
    int option)
{
     int nparams_clust;
     int num_params;
     int ndata_points;
     int repeat;
     double rissanen_const;
     double change,ll_new,ll_old;
     double epsilon;

     /* compute number of parameters per cluster */
     nparams_clust = 1+nbands+0.5*(nbands+1)*nbands;
     if(option==CLUSTER_DIAG) nparams_clust = 1+nbands+nbands;

     /* compute number of data points */
     ndata_points = Sig->ClassData.npixels*nbands;

     /* compute epsilon */
     epsilon = nparams_clust*log((double)ndata_points);
     epsilon *= 0.01;

     /* Perform initial regrouping */
     ll_new = regroup(Sig,nbands);

     /* Perform EM algorithm */
     change = 2*epsilon;
     do {
       ll_old = ll_new;
       reestimate(Sig,nbands,Rmin,option);

       ll_new = regroup(Sig,nbands);
       change = ll_new-ll_old;
       repeat = change>epsilon;
     } while(repeat);

     /* compute Rissanens expression */
     if(Sig->nsubclasses>0) {
       num_params = Sig->nsubclasses*nparams_clust - 1;
       rissanen_const = -ll_new + 0.5*num_params*log((double)ndata_points);
       return(rissanen_const);
     }
     else {
       return((double)0);
     }
}


static void reestimate(struct ClassSig *Sig, int nbands, double Rmin, int option)
{
     int i;
     int s;
     int b1,b2;
     double diff1,diff2;
     struct ClassData *Data;

     /* set data pointer */
     Data = (struct ClassData *)&(Sig->ClassData);
     
     /* Compute N */
     for(i=0; i<Sig->nsubclasses; i++) 
     {
       Sig->SubSig[i].N = 0;
       for(s=0; s<Data->npixels; s++)
         Sig->SubSig[i].N += (Data->p[s][i])*(Data->w[s]);
       Sig->SubSig[i].pi = Sig->SubSig[i].N;
     }

     /* Compute means and variances for each subcluster */
     for(i=0; i<Sig->nsubclasses; i++) 
     {
       /* Compute mean */
       for(b1=0; b1<nbands; b1++) 
       {
         Sig->SubSig[i].means[b1] = 0;
         for(s=0; s<Data->npixels; s++)
           Sig->SubSig[i].means[b1] += Data->p[s][i]*Data->x[s][b1]*Data->w[s];
         Sig->SubSig[i].means[b1] /= Sig->SubSig[i].N;
       }
	
       /* Compute R */
       for(b1=0; b1<nbands; b1++) 
       for(b2=b1; b2<nbands; b2++)
       {
         Sig->SubSig[i].R[b1][b2] = 0;
         for(s=0; s<Data->npixels; s++)
         {
           diff1 = Data->x[s][b1] - Sig->SubSig[i].means[b1];
           diff2 = Data->x[s][b2] - Sig->SubSig[i].means[b2];
           Sig->SubSig[i].R[b1][b2] += Data->p[s][i]*diff1*diff2*Data->w[s];
         }
         Sig->SubSig[i].R[b1][b2] /= Sig->SubSig[i].N;
         Sig->SubSig[i].R[b2][b1] = Sig->SubSig[i].R[b1][b2];
       }

       /* Regularize matrix */
       for(b1=0; b1<nbands; b1++) {
         Sig->SubSig[i].R[b1][b1] += Rmin;
       }

       if(option==CLUSTER_DIAG) DiagonalizeMatrix(Sig->SubSig[i].R,nbands);
     }