cluster.c Source File

00001 #include        "sadie.h"
00002 
00003 typedef struct _MERGE {
00004     short  cls1;
00005     short  cls2;
00006     PIXEL  dist;
00007     struct _MERGE *next;
00008 } MERGE;
00009 
00010 /*-Copyright Information------------------------------------------------------*/
00011 /* Copyright (c) 1990 by the University of Arizona Digital Image Analysis Lab */
00012 /*----------------------------------------------------------------------------*/
00013 /*-General Information--------------------------------------------------------*/
00014 /*                                                                            */
00015 /*   This function determines the intrinsic groupings of pixels in feature    */
00016 /*   space, using an iterative K-means clustering algorithm. Clusters with    */
00017 /*   too few pixels are discarded, and clusters that are too close are merged */
00018 /*   in the iteration process. After the clusters have been determined, the   */
00019 /*   entire input image is classified in a final pass.                        */
00020 /*                                                                            */
00021 /*----------------------------------------------------------------------------*/
00022 /*-Background Information-----------------------------------------------------*/
00023 /*                                                                            */
00024 /*   Duda, R.D., and Hart, P.E.:                                              */
00025 /*   "Pattern Classification and Scene Analysis,"                             */
00026 /*   J. Wiley and Sons, New York, 1973.                                       */
00027 /*                                                                            */
00028 /*----------------------------------------------------------------------------*/
00029 /* Interface Information------------------------------------------------------*/
00030 void CLUSTER (
00031 IMAGE *in,      /*  I   Pointer to the input image.                           */
00032 short niter,    /*  I   Number of iterations.                                 */
00033 short nclass,   /*  I   Number of seed classes.                               */
00034 short csize,    /*  |   Minimum number of pixels in each class.               */
00035 PIXEL mthr,     /*  I   Threshold for merging clusters.                       */
00036 short iinc,     /*  I   Band increment.                                       */
00037 short jinc,     /*  I   Line increment.                                       */
00038 short kinc,     /*  I   Pixels/line increment.                                */
00039 PIXEL othr,     /*  I   Outlier threshold.                                    */
00040 IMAGE **out     /*  O   Address of a pointer to the output image.             */
00041 /*----------------------------------------------------------------------------*/
00042 ) { register short i, j, k, m, n, class, nsize, *mrgd=0;
00043     char   msg[SLEN];
00044     long   *cnt=0;
00045     double pinc, psum;
00046     PIXEL  *mean=0, *nwmn=0, *dev=0, gmin, gmax, ginc, gsum, gmig, gdiff, gval;
00047     MERGE  *head=0, *tail, *ptr, *newClass;
00048 
00049     if (TIMES) TIMING(T_CLUSTER);
00050     if (NAMES) {
00051         MESSAGE('I',"");
00052         MESSAGE('I',"CLUSTER");
00053         MESSAGE('I',"");
00054         sprintf(msg," Input image:                   %s",in->text);
00055         MESSAGE('I',msg);
00056         sprintf(msg," Number of iterations:          %d",niter);
00057         MESSAGE('I',msg);
00058         sprintf(msg," Number of seed classes:        %d",nclass);
00059         MESSAGE('I',msg);
00060         sprintf(msg," Minimum class size:            %d",csize);
00061         MESSAGE('I',msg);
00062         sprintf(msg," Class merging threshold:     %12.4e",mthr);
00063         MESSAGE('I',msg);
00064         sprintf(msg," Subsample increment: bands:    %d",iinc);
00065         MESSAGE('I',msg);
00066         sprintf(msg,"                      lines:    %d",jinc);
00067         MESSAGE('I',msg);
00068         sprintf(msg,"                      pixels:   %d",kinc);
00069         MESSAGE('I',msg);
00070         sprintf(msg," Outlier threshold:           %12.4e",othr);
00071         MESSAGE('I',msg);
00072         MESSAGE('I'," ...............");
00073     }
00074 
00075     /* check input */
00076     if (!CHECKIMG(in)) {
00077         MESSAGE('E'," Can't identify image.");
00078         goto the_end;
00079     } else if (iinc <= 0) {
00080         MESSAGE('E'," Band subsample increment must be greater than zero.");
00081         goto the_end;
00082     } else if (jinc <= 0) {
00083         MESSAGE('E'," Line subsample increment must be greater than zero.");
00084         goto the_end;
00085     } else if (kinc <= 0) {
00086         MESSAGE('E'," Pixel subsample increment must be greater than zero.");
00087         goto the_end;
00088     } else if (niter <= 0) {
00089         MESSAGE('E'," Number of iterations must be greater than zero.");
00090         goto the_end;
00091     } else if (nclass <= 0) {
00092         MESSAGE('E'," Number of classes must be greater than zero.");
00093         goto the_end;
00094     } else if (csize <= 0) {
00095         MESSAGE('E'," Minimum class size must be greater than zero.");
00096         goto the_end;
00097     } else if (mthr < 0.0) {
00098         MESSAGE('E'," Distance threshold must be greater than or equal to zero.");
00099         goto the_end;
00100     } else if (othr < 0.0) {
00101         MESSAGE('E'," Distance threshold must be greater than or equal to zero.");
00102         goto the_end;
00103     }
00104 
00105     /* create image of appropriate size */
00106     if (!CHECKIMG(*out)) GETMEM(1,in->nlin,in->npix,out);
00107     if (!*out) goto the_end;
00108 
00109     /* allocate space */
00110     nsize = nclass + 1;
00111     if (!(mrgd=(short *)malloc(nsize*sizeof(short)))) {
00112         MESSAGE('E'," Memory request failed.");
00113         goto the_end;
00114     }
00115     if (!(cnt=(long *)malloc(nsize*sizeof(long)))) {
00116         MESSAGE('E'," Memory request failed.");
00117         goto the_end;
00118     }
00119     if (!(mean=(PIXEL *)malloc(nsize*in->nbnd*sizeof(PIXEL)))) {
00120         MESSAGE('E'," Memory request failed.");
00121         goto the_end;
00122     }
00123     if (!(nwmn=(PIXEL *)malloc(nsize*in->nbnd*sizeof(PIXEL)))) {
00124         MESSAGE('E'," Memory request failed.");
00125         goto the_end;
00126     }
00127     if (!(dev=(PIXEL *)malloc(nsize*in->nbnd*sizeof(PIXEL)))) {
00128         MESSAGE('E'," Memory request failed.");
00129         goto the_end;
00130     }
00131     if (!(head=(MERGE *)malloc((((nclass-1)*nclass)/2+2)*sizeof(MERGE)))) {
00132         MESSAGE('E'," Memory request failed.");
00133         goto the_end;
00134     }
00135 
00136     /* initialize progress indicator */
00137     if (LINES  &&  !PROGRESS(psum=0.0)) goto the_end;
00138     pinc = 0.5/(double)(niter*(in->nlin/jinc + (in->nlin%jinc ? 1:0)) + in->nlin);
00139 
00140     /* initialize */
00141     for (i=0; i<in->nbnd; i++) {
00142         gmax = gmin = in->data[i][0][0];
00143         for (j=0; j<in->nlin; j++) {
00144             for (k=0; k<in->npix; k++) {
00145                 gmin = min(gmin,in->data[i][j][k]);
00146                 gmax = max(gmax,in->data[i][j][k]);
00147             }
00148         }
00149         ginc = (gmax-gmin) / (PIXEL)nclass;
00150         for (m=1; m<nsize; m++) {
00151             mean[i*nsize+m] = gmin + (PIXEL)(((double)m-0.5)*(double)ginc);
00152         }
00153     }
00154     tail = head + 1;
00155     tail->cls1 = 0;
00156     tail->cls2 = 0;
00157     tail->dist = mthr * (PIXEL)in->nbnd;
00158     tail->next = 0;
00159 
00160     /* iterate */
00161     for (m=0; m<niter; m++) {
00162 
00163         /* initialize */
00164         for (n=0; n<nsize; n++) {
00165             cnt[n]  = 0L;
00166             mrgd[n] = TRUE;
00167             for (i=0; i<in->nbnd; i++) {
00168                 nwmn[i*nsize+n] = (PIXEL)0;
00169                 dev[i*nsize+n]  = (PIXEL)0;
00170             }
00171         }
00172 
00173         /* cluster */
00174         for (class=1,j=0; j<in->nlin; j+=jinc) {
00175             for (k=0; k<in->npix; k+=kinc) {
00176                 gmin = (PIXEL)0;
00177                 for (i=0; i<in->nbnd; i+=iinc) {
00178                     gmin += (PIXEL)fabs((double)(in->data[i][j][k]-mean[i*nsize+class]));
00179                 }
00180                 for (n=1; n<=nclass; n++) {
00181                     if (n != class) {
00182                         gsum = (PIXEL)0;
00183                         for (i=0; i<in->nbnd && gmin>gsum; i+=iinc) {
00184                             gsum += (PIXEL)fabs((double)(in->data[i][j][k]-mean[i*nsize+n]));
00185                         }
00186                         if (gmin > gsum) {
00187                             gmin = gsum;
00188                             class = n;
00189                         }
00190                     }
00191                 }
00192                 (*out)->data[0][j][k] = (PIXEL)class;
00193                 for (i=0; i<in->nbnd; i+=iinc) {
00194                     nwmn[i*nsize+class] += in->data[i][j][k];
00195                 }
00196                 cnt[class] += 1L;
00197             }
00198             if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00199         }
00200 
00201         /* compute new mean vectors and average migration */
00202         gmig = (PIXEL)0;
00203         for (i=0; i<in->nbnd; i+=iinc) {
00204             for (n=1; n<=nclass; n++) {
00205                 if (cnt[n] > 0L) {
00206                     nwmn[i*nsize+n] /= (PIXEL)cnt[n];
00207                             gmig += (PIXEL)fabs((double)(nwmn[i*nsize+n]-mean[i*nsize+n]));
00208                             mean[i*nsize+n] = nwmn[i*nsize+n];
00209                 }
00210             }
00211         }
00212             gmig /= (double)(in->nbnd*nclass);
00213 
00214         /* compute new standard deviation vectors */
00215         for (j=0; j<in->nlin; j+=jinc) {
00216             for (k=0; k<in->npix; k+=kinc) {
00217                 for (i=0; i<in->nbnd; i+=iinc) {
00218                     class = (short)(*out)->data[0][j][k];
00219                     gdiff = in->data[i][j][k] - mean[i*nsize+class];
00220                     dev[i*nsize+class] += gdiff * gdiff;
00221                 }
00222             }
00223             if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00224         }
00225         for (n=1; n<=nclass; n++) {
00226             for (i=0; i<in->nbnd; i+=iinc) {
00227                 if (cnt[n] > 0L) {
00228                     dev[i*nsize+n] = (PIXEL)sqrt((double)dev[i*nsize+n]/(double)cnt[n]);
00229                 }
00230             }
00231         }
00232 
00233         /* eliminate classes that are too small */
00234         for (n=1; n<=nclass; n++) {
00235             if (cnt[n] < csize) {
00236                 for (j=n; j<nclass; j++) {
00237                     for (i=0; i<in->nbnd; i+=iinc) {
00238                         mean[i*nsize+j] = mean[i*nsize+j+1];
00239                         dev[i*nsize+j]  = dev[i*nsize+j+1];
00240                     }
00241                     cnt[j] = cnt[j+1];
00242                 }
00243                 nclass -= 1;
00244                 n -= 1;
00245             }
00246         }
00247 
00248         /* find and merge classes that are too close together */
00249         newClass = head + 2;
00250         head->next = tail;
00251         for (n=1; n<=nclass-1; n++) {
00252             for (k=n+1; k<=nclass; k++) {
00253                 gsum = (PIXEL)0;
00254                 for (i=0; i<in->nbnd; i+=iinc) {
00255                     gval = dev[i*nsize+n] + dev[i*nsize+k];
00256                     if (gval != (PIXEL)0  &&  (gval=(PIXEL)fabs((double)(mean[i*nsize+n]-mean[i*nsize+k]))/gval) < mthr) {
00257                             gsum += gval;
00258                             if (i+iinc >= in->nbnd) {
00259                             for (ptr=head; ptr->next->dist<=gsum; ptr=ptr->next);
00260                                                 newClass->cls1 = n;
00261                                                 newClass->cls2 = k;
00262                                                 newClass->dist = gsum;
00263                                                 newClass->next = ptr->next;
00264                                                 ptr->next = newClass++;
00265                         }
00266                     } else {
00267                         break;
00268                     }
00269                 }
00270             }
00271         }
00272                 for (ptr=head->next; ptr!=tail; ptr=ptr->next) {
00273                     if (mrgd[ptr->cls1] == TRUE  &&  mrgd[ptr->cls2] == TRUE) {
00274                 for (i=0; i<in->nbnd; i+=iinc) {
00275                     gval = mean[i*nsize+ptr->cls1]*cnt[ptr->cls1] + mean[i*nsize+ptr->cls2]*cnt[ptr->cls2];
00276                     mean[i*nsize+ptr->cls1] = gval / (PIXEL)(cnt[ptr->cls1]+cnt[ptr->cls2]);
00277                 }
00278                 for (j=ptr->cls2; j<=nclass-1; j++) {
00279                     for (i=0; i<in->nbnd; i++) {
00280                         mean[i*nsize+j] = mean[i*nsize+j+1];
00281                     }
00282                 }
00283                 mrgd[ptr->cls1] = FALSE;
00284                 mrgd[ptr->cls2] = FALSE;
00285                 nclass -= 1;
00286             }
00287         }
00288 
00289         /* output iteration statistics */
00290         if (m == 0) {
00291             MESSAGE('I',"");
00292             MESSAGE('I'," Mean vector migration averaged over all classes:");
00293             MESSAGE('I'," Iteration     Migration     Number of classes after merging");
00294         }
00295         sprintf(msg,"%6d%18.4e%20d",m+1,(double)gmig,nclass);
00296         MESSAGE('I',msg);
00297     }
00298 
00299     /* initialize */
00300     for (n=0; n<nsize; n++) {
00301         cnt[n] = 0L;
00302         for (i=0; i<in->nbnd; i++) {
00303             nwmn[i*nsize+n] = (PIXEL)0;
00304         }
00305     }
00306 
00307     /* classify */
00308     for (class=1,j=0; j<in->nlin; j++) {
00309         for (k=0; k<in->npix; k++) {
00310             class = max(1,class);
00311             gmin = (PIXEL)0;
00312             for (i=0; i<in->nbnd; i++) {
00313                 gmin += (PIXEL)fabs((double)(in->data[i][j][k]-mean[i*nsize+class]));
00314             }
00315             for (n=1; n<=nclass; n++) {
00316                 if (n != class) {
00317                     gsum = (PIXEL)0;
00318                     for (i=0; i<in->nbnd && gmin>gsum; i++) {
00319                         gsum += (PIXEL)fabs((double)(in->data[i][j][k]-mean[i*nsize+n]));
00320                     }
00321                     if (gmin > gsum) {
00322                         gmin = gsum;
00323                         class = n;
00324                     }
00325                 }
00326             }
00327             for (i=0; i<in->nbnd ; i++) {
00328                         if ((PIXEL)fabs((double)(in->data[i][j][k]-mean[i*nsize+class])) > othr*dev[i*nsize+class]) {
00329                             class = 0;
00330                             break;
00331                         }
00332             }
00333                     (*out)->data[0][j][k] = class;
00334             if (class != 0) {
00335                 for (i=0; i<in->nbnd; i++) {
00336                     nwmn[i*nsize+class] += in->data[i][j][k];
00337                 }
00338             }
00339             cnt[class] += 1L;
00340         }
00341         if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00342     }
00343 
00344     /* compute new mean vectors and average migration */
00345     gmig = (PIXEL)0;
00346     for (i=0; i<in->nbnd; i++) {
00347         for (n=1; n<=nclass; n++) {
00348             if (cnt[n] > 0L) {
00349                 nwmn[i*nsize+n] /= (PIXEL)cnt[n];
00350                 gmig += (PIXEL)fabs((double)(nwmn[i*nsize+n]-mean[i*nsize+n]));
00351                         mean[i*nsize+n] = nwmn[i*nsize+n];
00352             }
00353         }
00354     }
00355     gmig /= (double)(in->nbnd*nclass);
00356 
00357     /* compute new standard deviation vectors */
00358         for (i=0; i<nsize*in->nbnd; i++) {
00359             dev[i] = (PIXEL)0;
00360         }
00361     for (j=0; j<in->nlin; j++) {
00362         for (k=0; k<in->npix; k++) {
00363             for (i=0; i<in->nbnd; i++) {
00364                 class = (short)(*out)->data[0][j][k];
00365                 if (class != 0) {
00366                     gdiff = in->data[i][j][k] - mean[i*nsize+class];
00367                     dev[i*nsize+class] += gdiff * gdiff;
00368                 }
00369             }
00370         }
00371         if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00372     }
00373     for (n=1; n<=nclass; n++) {
00374         for (i=0; i<in->nbnd; i++) {
00375             if (cnt[n] > 0L) {
00376                 dev[i*nsize+n] = (PIXEL)sqrt((double)dev[i*nsize+n]/(double)cnt[n]);
00377             }
00378         }
00379     }
00380 
00381     /* output iteration statistics */
00382     sprintf(msg,"     *%18.4e%20d",(double)gmig,nclass);
00383     MESSAGE('I',msg);
00384 
00385     /* output classification statistics */
00386     MESSAGE('I',"");
00387     MESSAGE('I'," Mean and standard deviation vectors of each class:");
00388     sprintf(msg," Class  ");
00389     for (i=0; i<in->nbnd; i++) {
00390         sprintf(msg+strlen(msg),"           Band %-11d",i+1);
00391     }
00392     sprintf(msg+strlen(msg),"# Pixels     %%");
00393     MESSAGE('I',msg);
00394     sprintf(msg,"%6d",0);
00395     for (i=0; i<in->nbnd; i++) {
00396         sprintf(msg+strlen(msg),"              N/A          ","");
00397     }
00398     sprintf(msg+strlen(msg),"%10ld%7.1f   (Threshold class)",cnt[0],100.0*(double)cnt[0]/((double)in->nlin*(double)in->npix));
00399     MESSAGE('I',msg);
00400     for (n=1; n<=nclass; n++) {
00401         sprintf(msg,"%6d",n);
00402         for (i=0; i<in->nbnd; i++) {
00403             sprintf(msg+strlen(msg),"%13.4e +/-%10.4e",(double)mean[i*nsize+n],(double)dev[i*nsize+n]);
00404         }
00405         sprintf(msg+strlen(msg),"%10ld%7.1f",cnt[n],100.0*(double)cnt[n]/((double)in->nlin*(double)in->npix));
00406         MESSAGE('I',msg);
00407     }
00408 
00409     the_end:
00410     if (mrgd)  free(mrgd);
00411     if (cnt)   free(cnt);
00412     if (mean)  free(mean);
00413     if (nwmn)  free(nwmn);
00414     if (dev)   free(dev);
00415     if (head)  free(head);
00416     if (TIMES) TIMING(T_EXIT);
00417 }