subcluster.cpp

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

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



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

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

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

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

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

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

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

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

static void normalize_pi(ClassSig *Sig);

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

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

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

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

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

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


int subcluster(
    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;
    ClassSig *Sig;
    static SigSet Smin;

    status = 0;

    /* set class pointer */
    Sig = &(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 = 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(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(
    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(ClassSig *Sig, int nbands, double Rmin, int option)
{
     int i;
     int s;
     int b1,b2;
     double diff1,diff2;
     ClassData *Data;

     /* set data pointer */
     Data = &(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);
     }

     /* Normalize probabilities for subclusters */
     normalize_pi(Sig);

     /* Compute constants */
     compute_constants(Sig,nbands);
     normalize_pi(Sig);
}


static double regroup(ClassSig *Sig, int nbands)
{
   int s;
   int i;
   double tmp;
   double maxlike;
   double likelihood;