subcluster.cpp

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

   double subsum;
   ClassData *Data;

   /* set data pointer */
   Data = &(Sig->classData);

   /* compute likelihoods */
   likelihood = 0;
   for(s=0; s<Data->npixels; s++)
   {

     for(i=0; i<Sig->nsubclasses; i++)
     {
       tmp = loglike(Data->x[s],&(Sig->subSig[i]),nbands);
       Data->p[s][i] = tmp;
       if(i==0) maxlike = tmp;
       if(tmp>maxlike) maxlike = tmp;
     }

     subsum = 0;
     for(i=0; i<Sig->nsubclasses; i++)
     {
       tmp = exp( Data->p[s][i]-maxlike )*Sig->subSig[i].pi;
       subsum += tmp;
       Data->p[s][i] = tmp;
     }
     likelihood += log(subsum) + maxlike;

     for(i=0; i<Sig->nsubclasses; i++)
       Data->p[s][i] /= subsum;
   }

   return(likelihood);
}


static void reduce_order(
    ClassSig *Sig,
    int nbands,
    int *min_ii,
    int *min_jj)
{
    int i,j;
    int min_i,min_j;
    double dist;
    double min_dist;
    SubSig *SubSig1,*SubSig2;

    static int first=1;
    SigSet S;
    static ClassSig *Sig3;
    static SubSig *SubSig3;

    /* allocate scratch space first time subroutine is called */
    if(first)
    {
      I_InitSigSet (&S);
      I_SigSetNBands (&S, nbands);
      Sig3 = I_NewClassSig(&S);
      I_NewSubSig (&S, Sig3);
      SubSig3 = Sig3->subSig;
      first = 0;
    }

    if(Sig->nsubclasses>1)
    {
      /* find the closest subclasses */
      for(i=0; i<Sig->nsubclasses-1; i++)
      for(j=i+1; j<Sig->nsubclasses; j++)
      {
        dist = distance(&(Sig->subSig[i]),&(Sig->subSig[j]),nbands);
        if((i==0)&&(j==1))
        {
          min_dist = dist;
          min_i = i;
          min_j = j;
        }
        if(dist<min_dist)
        {
          min_dist = dist;
          min_i = i;
          min_j = j;
        }
      }

      /* Save result for output */
      *min_ii = min_i; *min_jj = min_j;

      /* Combine Subclasses */
      SubSig1 = &(Sig->subSig[min_i]);
      SubSig2 = &(Sig->subSig[min_j]);
      add_SubSigs(SubSig1,SubSig2,SubSig3,nbands);
      copy_SubSig(SubSig3,SubSig1,nbands);

      /* remove extra subclass */
      for(i=min_j; i<Sig->nsubclasses-1; i++)
        copy_SubSig(&(Sig->subSig[i+1]),&(Sig->subSig[i]),nbands);

      /* Remove last Subclass */
      /* (Sig->nsubclasses)--; */
      I_DeallocSubSig(Sig); 

      /* Rerun compute_constants */
      compute_constants(Sig,nbands);
      normalize_pi(Sig);
    }
}


static double loglike(double *x, SubSig *SubSig, int nbands)
{
    int b1,b2;
    double diff1,diff2;
    double sum;

    sum = 0;
    for(b1=0; b1<nbands; b1++) 
    for(b2=0; b2<nbands; b2++)
    {
      diff1 = x[b1]-SubSig->means[b1];
      diff2 = x[b2]-SubSig->means[b2];
      sum += diff1*diff2*SubSig->Rinv[b1][b2];
    }

    sum = -0.5*sum + SubSig->cnst;
    return(sum);
}


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

    static int first=1;
    SigSet S;
    static ClassSig *Sig3;
    static SubSig *SubSig3;


    /* allocate scratch space first time subroutine is called */
    if(first)
    {
      I_InitSigSet (&S);
      I_SigSetNBands (&S, nbands);
      Sig3 = I_NewClassSig(&S);
      I_NewSubSig (&S, Sig3);
      SubSig3 = Sig3->subSig;
      first = 0;
    }

    /* form SubSig3 by adding SubSig1 and SubSig2 */
    add_SubSigs(SubSig1,SubSig2,SubSig3,nbands);

    /* compute constant for SubSig3 */
    compute_constants(Sig3,nbands);

    /* compute distance */
    dist = SubSig1->N*SubSig1->cnst + SubSig2->N*SubSig2->cnst
           - SubSig3->N*SubSig3->cnst;

    return(dist);
}


/**********************************************************/
/* invert matrix and compute Sig->subSig[i].cnst          */
/**********************************************************/
static void compute_constants(ClassSig *Sig, int nbands)
{
   int i;
   int b1,b2;
   double det_man;
   int    det_exp;

   static int first=1;
   static int *indx;
   static double **y;
   static double *col;


   /* allocate memory first time subroutine is called */
   if(first)
   {
     indx = G_alloc_ivector(nbands);
     y = G_alloc_matrix(nbands,nbands); 
     col = G_alloc_vector(nbands);
     first = 0;
   }

   /* invert matrix and compute constant for each subclass */
   for(i=0; i<Sig->nsubclasses; i++) {
     for(b1=0; b1<nbands; b1++)
     for(b2=0; b2<nbands; b2++)
       Sig->subSig[i].Rinv[b1][b2] = Sig->subSig[i].R[b1][b2];

     clust_invert(Sig->subSig[i].Rinv,nbands,&det_man,&det_exp,indx,y,col);

     Sig->subSig[i].cnst = (-nbands/2.0)*log(2*PI) 
                            - 0.5*log(det_man) - 0.5*det_exp*log(10.0);
   } 
}


static void normalize_pi(ClassSig *Sig)
{
    int i;
    double sum;

    sum = 0.0;
    for(i=0; i<Sig->nsubclasses; i++) sum += Sig->subSig[i].pi;

    if(sum>0) {
      for(i=0; i<Sig->nsubclasses; i++) Sig->subSig[i].pi /= sum;
    }
    else {
      for(i=0; i<Sig->nsubclasses; i++) Sig->subSig[i].pi = 0.0;
    }
}


/*******************************************/
/* add SubSig1 and SubSig2 to form SubSig3 */
/*******************************************/
static void add_SubSigs(
    SubSig *SubSig1,
    SubSig *SubSig2,
    SubSig *SubSig3,
    int nbands)
{
    int b1,b2;
    double wt1,wt2;
    double tmp;

    wt1 = SubSig1->N/(SubSig1->N + SubSig2->N);
    wt2 = 1 - wt1;

    /* compute means */
    for(b1=0; b1<nbands; b1++)
      SubSig3->means[b1] = wt1*SubSig1->means[b1] + wt2*SubSig2->means[b1];

    /* compute covariance */
    for(b1=0; b1<nbands; b1++)
    for(b2=b1; b2<nbands; b2++)
    {
      tmp = (SubSig3->means[b1]-SubSig1->means[b1])
            *(SubSig3->means[b2]-SubSig1->means[b2]);
      SubSig3->R[b1][b2] = wt1*(SubSig1->R[b1][b2] + tmp);
      tmp = (SubSig3->means[b1]-SubSig2->means[b1])
            *(SubSig3->means[b2]-SubSig2->means[b2]);
      SubSig3->R[b1][b2] += wt2*(SubSig2->R[b1][b2] + tmp);
      SubSig3->R[b2][b1] = SubSig3->R[b1][b2];
    }

    /* compute pi and N */
    SubSig3->pi = SubSig1->pi + SubSig2->pi;
    SubSig3->N = SubSig1->N + SubSig2->N;
}


/**********************/
/* saves Sig1 to Sig2 */
/**********************/
static void save_ClassSig(
    ClassSig *Sig1,
    SigSet *S,
    int nbands)
{
    ClassSig *Sig2;

    I_InitSigSet (S);
    I_SigSetNBands (S, nbands);
    Sig2 = I_NewClassSig(S);
    while(Sig2->nsubclasses<Sig1->nsubclasses) I_NewSubSig(S, Sig2);
    copy_ClassSig(Sig1,Sig2,nbands);
}


/*********************/
/* copy Sig1 to Sig2 */
/*********************/
static void copy_ClassSig(
    ClassSig *Sig1,
    ClassSig *Sig2,
    int nbands)
{
    int i;

    Sig2->classnum = Sig1->classnum;
    Sig2->title = Sig1->title;
    Sig2->used = Sig1->used;
    Sig2->type = Sig1->type;
    Sig2->nsubclasses = Sig1->nsubclasses;
    for(i=0; i<Sig1->nsubclasses; i++)
      copy_SubSig(&(Sig1->subSig[i]),&(Sig2->subSig[i]),nbands);
}


/***************************/
/* copy SubSig1 to SubSig2 */
/***************************/
static void copy_SubSig(
    SubSig *SubSig1,
    SubSig *SubSig2,
    int nbands)
{
    int b1,b2;

    SubSig2->N = SubSig1->N;
    SubSig2->pi = SubSig1->pi;
    SubSig2->cnst = SubSig1->cnst;
    SubSig2->used = SubSig1->used;

    for(b1=0; b1<nbands; b1++)
      SubSig2->means[b1] = SubSig1->means[b1];

    for(b1=0; b1<nbands; b1++)
    for(b2=0; b2<nbands; b2++)
    {
      SubSig2->R[b1][b2] = SubSig1->R[b1][b2];
      SubSig2->Rinv[b1][b2] = SubSig1->Rinv[b1][b2];
    }
}

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

    for(b1=0; b1<nbands; b1++)
    for(b2=0; b2<nbands; b2++)
      if(b1!=b2) R[b1][b2] = 0;
}

#if 0
static void list_Sig(ClassSig *Sig, int nbands)
{
    int i,j,k;

    for(i=0; i<Sig->nsubclasses; i++)
    {
      print

f("Subclass %d: pi = %f, ",i,Sig->subSig[i].pi);
      printf("cnst = %f\n",Sig->subSig[i].cnst);
      for(j=0; j<nbands; j++)
      {
        printf("%f;    ",Sig->subSig[i].means[j]);
        for(k=0; k<nbands; k++)
          printf("%f ",Sig->subSig[i].R[j][k]);
        printf("\n");
      }
      printf("\n");
    }
}

static void print_class(ClassSig *Sig, char *fname)
{
    FILE *fp;
    int s,i;


    if((fp=fopen(fname,"w"))==NULL) {
      fprintf(stderr,"can't open data file\n"); exit(1);}

    for(s=0; s<Sig->classData.npixels; s++)
    {
      /* fprintf(fp,"Pixel number %d:  ", s); */
      for(i=0; i<Sig->nsubclasses; i++)
        fprintf(fp,"%f  ",Sig->classData.p[s][i]);
      fprintf(fp,"\n");
    }

    fclose(fp);

}
#endif