clust.cpp

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

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


double AverageVariance(ClassSig *Sig, int nbands);


int Clust(int argc, char *argv[]) 
{
    FILE *fp,*info_fp;

    int i,j,k;
    int init_num_of_subclasses,max_num;
    int vector_dimension;
    int nclasses;
    int num_of_samples;
    char option1[16];
    int  option2;
    double Rmin;
    char fname[512];
    SigSet S;
    ClassSig *Sig;

    /* set level of diagnostic printing */
    //clusterMessageVerboseLevel = 2;

    /* print usage message if arguments are not valid */
    if((argc!=4) && (argc!=5) && (argc!=6)) {
      fprintf(stderr,"\n\nUsage: %s #_subclasses info_file output_params [option1 option2]\n\n",argv[0] );
      fprintf(stderr,"    #_subclasses - initial number of clusters for each class\n\n");
      fprintf(stderr,"    info_file - name of file which contains the following information:\n");
      fprintf(stderr,"      <# of classes>\n");
      fprintf(stderr,"      <data vector length>\n");
      fprintf(stderr,"      <class 1 data file name> <# of data vectors in class 1>\n");
      fprintf(stderr,"      <class 2 data file name> <# of data vectors in class 2>\n");
      fprintf(stderr,"                     .                        .\n");
      fprintf(stderr,"                     .                        .\n");
      fprintf(stderr,"                     .                        .\n");
      fprintf(stderr,"      <last class data file name> <# of data vectors in last class>\n\n");
      fprintf(stderr,"    output_params - name of file containing output clustering");
      fprintf(stderr," parameters\n\n");
      fprintf(stderr,"    option1 - (optional) controls clustering model\n");
      fprintf(stderr,"      full - (default) use full convariance matrices\n");
      fprintf(stderr,"      diag - use diagonal convariance matrices\n\n");
      fprintf(stderr,"    option2 - (optional) controls number of clusters\n");
      fprintf(stderr,"      0 - (default) estimate number of clusters\n");
      fprintf(stderr,"      n - use n clusters in mixture model with n<#_subclasses\n");
      exit(1);
    } 

    /* read number of initial subclasses to use */
    sscanf(argv[1],"%d",&init_num_of_subclasses);

/* Set option 1 */
    if(argc==4) {
      sprintf(option1,"full") ;
    }

    if((argc==5) || (argc==6)) {
      /* set default option 1 */
      sscanf(argv[4], "%s", option1) ;
      if((strcmp(option1, "full")!=0) && (strcmp(option1, "diag")!=0)) {
        fprintf(stderr,"\nInvalid option1: %s\n\n",option1);
        fprintf(stderr,"There are 2 valid assumptions:\n");
        fprintf(stderr,"    full - default option which allows full convariance matrices\n");
        fprintf(stderr,"    diag - use diagonal convariance matrices\n\n");
        exit(1);
      }
    }
          
/* Set option 2 */
    if((argc==4) || (argc==5)) {
      option2 = 0;
    }
          
    if(argc==6) {
      /* set default option 2 */
      sscanf(argv[5], "%d", &option2) ;
      if( (option2<0) || (option2>init_num_of_subclasses) ) {
        fprintf(stderr,"\nInvalid option2: %d \n\n",option2);
        fprintf(stderr,"There are 2 valid assumptions:\n");
        fprintf(stderr,"      0 - (default) estimate number of clusters\n");
        fprintf(stderr,"      n - use n clusters in mixture model with n<#_subclasses\n\n");
        exit(1);
      }
    }


    /* open information file */
    if((info_fp=fopen(argv[2],"r"))==NULL) {
      fprintf(stderr,"can't open information file");
      exit(1);
    }

    /* read number of classes from info file */
    fscanf(info_fp,"%d\n",&nclasses);

    /* read vector dimension from info file */
    fscanf(info_fp,"%d\n",&vector_dimension);


    /* 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 = I_NewClassSig(&S);
      I_SetClassTitle (Sig, "test class signature");
      for(i=0; i<init_num_of_subclasses; i++)
        I_NewSubSig (&S, Sig);
    }

    /* Read data for each class */
    for(k=0; k<nclasses; k++) {
      /* read class k data file name */
      fscanf(info_fp,"%s",fname);

      /* read number of samples for class k */
      fscanf(info_fp,"%d\n",&num_of_samples);

      Sig = &(S.classSig[k]);

      I_AllocClassData (&S, Sig, num_of_samples);

      /* Read Data */
      if((fp=fopen(fname,"r"))==NULL) {
        fprintf(stderr,"can't open data file %s", fname);
        exit(1);
      }

      for(i=0; i<Sig->classData.npixels; i++) {
        for(j=0; j<vector_dimension; j++) {
          fscanf(fp,"%lf",&(Sig->classData.x[i][j]) );
        }
        fscanf(fp,"\n");
      }
      fclose(fp);

      /* 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];
      }
    }
    fclose(info_fp);


    /* Compute the average variance over all classes */
    Rmin = 0;
    for(k=0; k<nclasses; k++) {
      Sig = &(S.classSig[k]);
      Rmin += AverageVariance(Sig, vector_dimension);
    }
    Rmin = Rmin/(COVAR_DYNAMIC_RANGE*nclasses);

    /* Perform clustering for each class */
    for(k=0; k<nclasses; k++) {

      Sig = &(S.classSig[k]);

      //if(1<=clusterMessageVerboseLevel) {
      //  fprintf(stdout,"Start clustering class %d\n\n",k);
      //}

      if(strcmp(option1, "diag")==0) {
        /* assume covariance matrices to be diagonal */
        subcluster(&S,k,option2,(int)CLUSTER_DIAG,Rmin,&max_num);

      } else {
        /* no assumption for covariance matrices */
        subcluster(&S,k,option2,(int)CLUSTER_FULL,Rmin,&max_num);
      }

      //if(2<=clusterMessageVerboseLevel) {
      //  fprintf(stdout,"Maximum number of subclasses = %d\n",max_num);
      //}

      I_DeallocClassData(&S, Sig);
    }

    /* Write out result to output parameter file */
    if((fp=fopen(argv[3],"w"))==NULL) {
      fprintf(stderr,"can't open parameter output file");
      exit(1);
    }
    I_WriteSigSet(fp, &S);
    fclose(fp);

    /* De-allocate cluster signature memory */
    I_DeallocSigSet(&S); 

    return(0);
}


double AverageVariance(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);
}