mindist.c Source File

00001 #include        "sadie.h"
00002 
00003 /*-Copyright Information------------------------------------------------------*/
00004 /* Copyright (c) 1992 by the University of Arizona Digital Image Analysis Lab */
00005 /*----------------------------------------------------------------------------*/
00006 /*-General Information--------------------------------------------------------*/
00007 /*                                                                            */
00008 /*   This function performs a minimum distance classification on an image.    */ 
00009 /*                                                                            */
00010 /*   The discriminant function for classification of a pixel to a class is:   */ 
00011 /*                                                                            */ 
00012 /*   City-block distance classifier:                                          */ 
00013 /*      DC(X) = |X - M|                                                       */ 
00014 /*                                                                            */ 
00015 /*   Euclidean distance classifier:                                           */ 
00016 /*      DE(X) = (X - M)**t * (X - M)                                          */ 
00017 /*                                                                            */ 
00018 /*   Mahalanobis distance classifier:                                         */
00019 /*      DM(X) = (X - M)**t * C**-1 * (X - M)                                  */ 
00020 /*                                                                            */ 
00021 /*   where X is the feature vector, M is the class mean vector, and C**-1     */ 
00022 /*   is the inverse of the covariance matrix (equal for all classes).         */ 
00023 /*                                                                            */ 
00024 /*----------------------------------------------------------------------------*/ 
00025 /*-Interface Information------------------------------------------------------*/ 
00026 void MINDIST (
00027 IMAGE  *in,    /*  I  Pointer to the input image.                             */ 
00028 short  option, /*  I  Switch, indicating the type of classifier:              */
00029                /*     CITYBLOCK     -   DC(X)                                 */
00030                /*     EUCLIDEAN     -   DE(X)                                 */
00031                /*     MAHALANOBIS   -   DM(X)                                 */
00032 short  nlev,   /*  I  Number of classes                                       */ 
00033                /*     (and number of graylevels in the output image).         */ 
00034 double *mean,  /*  I  Pointer to the matrix[nlev][nbnd] of mean graylevels.   */ 
00035 double *cov,   /*  I  Pointer to the covariance matrix[nbnd][nbnd].           */ 
00036 IMAGE  **out   /*  O  Address of a pointer to the output image.               */ 
00037 /*--------------------------------------------------------------------------- */ 
00038 ) { register short i, j, k, m, n, l1, l2, ll, nl, np;                              
00039     char   msg[SLEN];                                                                
00040     PIXEL  *dist=0, mindist;                      
00041     double *inv=0;                                                          
00042                                                                                 
00043     if (TIMES) TIMING(T_MINDIST);
00044     if (NAMES) {                                                                    
00045         MESSAGE('I',"");                                                           
00046         MESSAGE('I',"MINDIST");                                                   
00047         MESSAGE('I',"");                                                           
00048         sprintf(msg," Input image:          %s",in->text);                          
00049         MESSAGE('I',msg);                                                          
00050         if (option == CITYBLOCK) {
00051             MESSAGE('I'," Classifier:           City-block");
00052         } else if (option == EUCLIDEAN) {
00053             MESSAGE('I'," Classifier:           Euclidean");
00054         } else if (option == MAHALANOBIS) {
00055             MESSAGE('I'," Classifier:           Mahalanobis");
00056         }
00057         sprintf(msg," Number of classes:    %d",nlev);                              
00058         MESSAGE('I',msg);
00059         sprintf(msg," Mean graylevels:");
00060         MESSAGE('I',msg);
00061         for (n=j=0; j<nlev; j++) {
00062             sprintf(msg,"  Class%3d:          ",j+1);
00063             for (k=0,m=strlen(msg); k<in->nbnd; k++,m=strlen(msg)) {
00064                 sprintf(msg+m,"%12.4e",mean[n++]);
00065             }
00066             MESSAGE('I',msg);
00067         }
00068         MESSAGE('I'," Covariance matrix:");
00069         for (n=j=0; j<in->nbnd; j++) {
00070                 sprintf(msg,"  Band%4d:          ",j+1);
00071             for (k=0,m=strlen(msg); k<in->nbnd; k++,m=strlen(msg)) {
00072                 sprintf(msg+m,"%12.4e",cov[n++]);
00073             }
00074             MESSAGE('I',msg);
00075         }
00076         MESSAGE('I'," ...............");
00077     }                                                                           
00078                                                                                 
00079     /* check input  */                 
00080     if (!CHECKIMG(in)) {                                                             
00081         MESSAGE('E',"Can't identify image.");                                     
00082         goto the_end;                                                             
00083     } else if (option != CITYBLOCK  &&  option != EUCLIDEAN  &&  option != MAHALANOBIS) {
00084         MESSAGE('E'," Option must be one of CITYBLOCK, EUCLIDEAN, or MAHALANOBIS.");
00085         goto the_end;
00086     } else if (nlev <= 0) {                                                            
00087         MESSAGE('E',"Number of classes must be greater than zero.");               
00088         goto the_end;                                                              
00089     }                                                                           
00090 
00091     /* create image of appropriate size */
00092     if (!CHECKIMG(*out)) GETMEM(1,in->nlin,in->npix,out);                          
00093     if (!*out) goto the_end;
00094 
00095 
00096     /* allocate space for  distance array */   
00097     if (!(dist=(PIXEL *)malloc(nlev*sizeof(PIXEL)))) {                                                                  
00098             MESSAGE('E',"Memory request failed.");                         
00099             goto the_end;                                                             
00100     }                                                                          
00101 
00102 
00103     /* City-block distance classifier */           
00104     if (option == CITYBLOCK) {                                    
00105         for (nl=0; nl<in->nlin; nl++) {                                                 
00106             for (np=0; np<in->npix; np++) {                                                
00107                 for (k=0; k<nlev; k++) {                                                     
00108                     dist[k] = 0.0;                                                        
00109                     for (i=0; i<in->nbnd; i++) {                                               
00110                         l1 = k*in->nbnd + i;                                                      
00111                         dist[k] += fabs(in->data[i][nl][np]-mean[l1]);                    
00112                     }                                                                     
00113                 }                                                                       
00114                 mindist = dist[0];                                                       
00115                 for (i=0; i<nlev; i++) {                                                     
00116                     if (dist[i] <= mindist) {                                                  
00117                         (*out)->data[0][nl][np] = (PIXEL)(i+1);
00118                     }                              
00119                 }
00120             }
00121         }
00122 
00123     /* Euclidean distance classifier */              
00124     } else if (option == EUCLIDEAN) {                                    
00125         for (nl=0; nl<in->nlin; nl++) {                                                 
00126             for (np=0; np<in->npix; np++) {                                                
00127                 for (k=0; k<nlev; k++) {                                                     
00128                     dist[k] = 0.0;                                                         
00129                     for (i=0; i<in->nbnd; i++) {                                               
00130                         l1 = k*in->nbnd + i;                                                      
00131                         dist[k] += (in->data[i][nl][np]-mean[l1])*(in->data[i][nl][np]-mean[l1]);                                                                            
00132                     }                                                                     
00133                 }                                                                       
00134                 mindist = dist[0];                                                        
00135                 for (i=0; i<nlev; i++) {                                                     
00136                     if (dist[i] <= mindist) {                                                    
00137                         (*out)->data[0][nl][np] = (PIXEL)(i+1);                                
00138                     }
00139                 }                                                                       
00140             }                                                                         
00141         }                                                                          
00142 
00143     /* Mahalanobis distance classifier */                
00144     } else if (option == MAHALANOBIS) {                                          
00145                                                                                 
00146         /* allocate space for inverse of covariance matrix and distance array */   
00147         if (!(inv=(double *)malloc(in->nbnd*in->nbnd*sizeof(double)))) {                                                                     
00148                 MESSAGE('E',"Memory request failed.");                            
00149                 goto the_end;                                                             
00150         }                                                                          
00151         if (!(dist=(PIXEL *)malloc(nlev*sizeof(PIXEL)))) {                                                                  
00152                 MESSAGE('E',"Memory request failed.");                         
00153                 goto the_end;                                                             
00154         }                                                                          
00155         MATRIX_INVERT(cov,in->nbnd,inv);
00156                                                   
00157         for (nl=0; nl<in->nlin; nl++) {                                                  
00158             for (np=0; np<in->npix; np++) {                                                
00159                 for (k=0; k<nlev; k++) {                                                     
00160                     dist[k] = 0.0;                                                        
00161                     for (i=0; i<in->nbnd; i++) {                                               
00162                         for (j=0; j<in->nbnd; j++) {                                              
00163                                 l1 = k*in->nbnd + i;                                                      
00164                                 l2 = k*in->nbnd + j;                                                      
00165                             ll = i*in->nbnd + j;                                                  
00166                             dist[k]+=(in->data[i][nl][np]-mean[l1])*inv[ll]*(in->data[j][nl][np] - mean[l2]);                                                              
00167                         }
00168                     }                                                                    
00169                 }                                                                       
00170                 mindist = dist[0];                                                       
00171                 for (i=0; i<nlev; i++) {                                                     
00172                     if (dist[i] <= mindist) {                                                
00173                         (*out)->data[0][nl][np] = (PIXEL)(i+1);
00174                     }
00175                 }                                                                       
00176             }                                                                          
00177         }                                                                            
00178     }
00179 
00180     the_end:
00181     if (inv)   free(inv);
00182     if (dist)  free(dist);
00183     if (TIMES) TIMING(T_EXIT);
00184 }