match.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 transforms the graylevels of one image to match those      */
00009 /*   of a second image.  The first option uses a least-squares match, and     */
00010 /*   the second option uses the cumulative distribution functions of both     */
00011 /*   images.                                                                  */
00012 /*                                                                            */
00013 /*----------------------------------------------------------------------------*/
00014 /*-Background Information-----------------------------------------------------*/
00015 /*                                                                            */
00016 /*   Dallas, W.J. and Mauser, W.:                                             */
00017 /*   "Preparing Pictures for Visual Comparison."                              */
00018 /*   Applied Optics, Vol. 19, No. 21, November 1980                           */
00019 /*                                                                            */
00020 /*----------------------------------------------------------------------------*/
00021 /*-Interface Information------------------------------------------------------*/
00022 void MATCH (
00023 IMAGE *in1,      /*  I  Pointer to the transformation image.                  */
00024 IMAGE *in2,      /*  I  Pointer to the reference image.                       */
00025 short nlev,      /*  I  Number of graylevels to be used in the transformation.*/
00026 short option,    /*  I  Switch, indicating type of image matching             */
00027                  /*     LSQ   -   Images matched using gray levels            */
00028                  /*     CDF   -   Images matched using cumulative dist fncns  */
00029 IMAGE **out      /*  O  Address of a pointer to the output image.             */
00030                  /*     The input image may be overwritten (out = &in1).      */
00031 /*----------------------------------------------------------------------------*/
00032 ) { register short i, j, k, index, index2;
00033     char   msg[SLEN];
00034     long   *count=0, *hist1= 0, *hist2=0;
00035     double offset, offset2, factor, factor2,  pinc, psum, *cdf1=0, *cdf2=0;
00036     PIXEL  *table=0;
00037     char minlbl[10], maxlbl[10];
00038 
00039     if (TIMES) TIMING(T_MATCH);
00040     if (NAMES) {
00041         MESSAGE('I',"");
00042         MESSAGE('I',"MATCH");
00043         MESSAGE('I',"");
00044         sprintf(msg," Transformation input image:   %s",in1->text);
00045         MESSAGE('I',msg);
00046         sprintf(msg," Reference input image:        %s",in2->text);
00047         MESSAGE('I',msg);
00048         sprintf(msg," Lookup table entries:         %d",nlev);
00049         MESSAGE('I',msg);
00050         if (option == LSQ) {
00051             MESSAGE('I'," Least squares transformation");
00052         } else {
00053             MESSAGE('I'," CDF transformation");
00054         }
00055         MESSAGE('I'," ...............");
00056     }
00057 
00058     /* check input */
00059     if (!CHECKIMG(in1)) {
00060         MESSAGE('E'," Can't identify transformation image.");
00061         goto the_end;
00062     } else if (!CHECKIMG(in2)) {
00063         MESSAGE('E'," Can't identify reference image.");
00064         goto the_end;
00065     } else if (in1->nbnd != in2->nbnd) {
00066         MESSAGE('E'," Number of bands must be identical in both images.");
00067         goto the_end;
00068     } else if (in1->nlin != in2->nlin) {
00069         MESSAGE('E'," Number of lines must be identical in both images.");
00070         goto the_end;
00071     } else if (in1->npix != in2->npix) {
00072         MESSAGE('E'," Number of pixels/line must be identical in both images.");
00073         goto the_end;
00074     } else if (nlev <= 0) {
00075         MESSAGE('E'," Number of graylevel bins must be greater than zero.");
00076         goto the_end;
00077     } else if (option != LSQ  &&  option != CDF) {
00078         MESSAGE('E'," Option must be either LSQ or CDF.");
00079         goto the_end;
00080     }
00081 
00082     /* create image of appropriate size */
00083     if (!CHECKIMG(*out)) GETMEM(in1->nbnd,in1->nlin,in1->npix,out);
00084     if (!*out) goto the_end;
00085 
00086     /* initialize progress indicator */
00087     if (LINES  &&  !PROGRESS(psum=0.0)) goto the_end;
00088     pinc = (option == LSQ ? 0.5 : 1.0)/(double)in1->nbnd/(double)in1->nlin;
00089 
00090     /* allocate space for array */
00091     if (!(table=(PIXEL *)malloc(nlev*sizeof(PIXEL)))) {
00092         MESSAGE('E'," Memory request failed.");
00093         goto the_end;
00094     }
00095 
00096     /* initialize */
00097     RANGE(in1);
00098     offset = (double)in1->gmin;
00099     factor = (double)(nlev-1)/(double)(in1->gmax-in1->gmin);
00100 
00101     if (option == LSQ) {
00102 
00103         /* allocate space for array */
00104         if (!(count=(long *)malloc(nlev*sizeof(long)))) {
00105             MESSAGE('E'," Memory request failed.");
00106             goto the_end;
00107         }
00108 
00109         for (i=0; i<in1->nbnd; i++) {
00110 
00111             /* build look-up table */
00112             for (index=0; index<nlev; index++) {
00113                 table[index] = (PIXEL)0;
00114                 count[index] = 0L;
00115             }
00116             for (j=0; j<in1->nlin; j++) {
00117                 for (k=0; k<in1->npix; k++) {
00118                     index = (short)(((double)in1->data[i][j][k]-offset)*factor);
00119                     table[index] += in2->data[i][j][k];
00120                     count[index] += 1L;
00121                 }
00122                 if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00123             }
00124 
00125             /* average the table with the count array */
00126             for (index=0; index<nlev; index++) {
00127                 if (count[index] != 0) {
00128                     table[index] /= (double)count[index];
00129                 }
00130             }
00131 
00132             /* create output image */
00133             for (j=0; j<in1->nlin; j++) {
00134                 for (k=0; k<in1->npix; k++) {
00135                     index = (short)(((double)in1->data[i][j][k]-offset)*factor);
00136                     (*out)->data[i][j][k] = table[index];
00137                 }
00138                 if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00139             }
00140 
00141             /* output look-up table */
00142             sprintf(minlbl,"0");
00143             sprintf(maxlbl,"%d",nlev);
00144             PLOT(&table,1,nlev,minlbl,maxlbl,NORM);
00145         }
00146 
00147     } else {
00148 
00149         /* allocate space for histograms of reference and transform images and cdf's */
00150         if (!(hist1=(long *)malloc(in1->nbnd*nlev*sizeof(long)))) {
00151             MESSAGE('E'," Memory request failed.");
00152             goto the_end;
00153         }
00154         if (!(hist2=(long *)malloc(in2->nbnd*nlev*sizeof(long)))) {
00155             MESSAGE('E'," Memory request failed.");
00156             goto the_end;
00157         }
00158         if (!(cdf1=(double *)malloc(in1->nbnd*nlev*sizeof(double)))) {
00159             MESSAGE('E'," Memory request failed.");
00160             goto the_end;
00161         }
00162         if (!(cdf2=(double *)malloc(in2->nbnd*nlev*sizeof(double)))) {
00163             MESSAGE('E'," Memory request failed.");
00164             goto the_end;
00165         }
00166 
00167         /* compute cumulative distribution function for transformation image */
00168         RANGE(in1);
00169         HISTOGRM(in1,1,in1->gmin,in1->gmax,nlev,hist1);
00170         for (i=0; i<in1->nbnd; i++) {
00171             for (j=1; j<nlev; j++) {
00172                 hist1[i*nlev+j] += hist1[i*nlev+j-1];
00173             }
00174         }
00175 
00176         /* scale cumulative distribution function to [0,1] */
00177         for (i=0; i<in1->nbnd; i++) {
00178             for (j=0; j<nlev; j++) {
00179                 cdf1[i*nlev+j] = (double)(hist1[i*nlev+j]) / (double)(hist1[nlev-1]);
00180             }
00181         }
00182 
00183         /* compute cumulative distribution function for reference image */
00184         RANGE(in2);
00185         HISTOGRM(in2,1,in2->gmin,in2->gmax,nlev,hist2);
00186         for (i=0; i<in2->nbnd; i++) {
00187             for (j=1; j<nlev; j++) {
00188                 hist2[i*nlev+j] += hist2[i*nlev+j-1];
00189             }
00190         }
00191 
00192         /* scale cumulative distribution function to [0,1] */
00193         for (i=0; i<in1->nbnd; i++) {
00194             for (j=0; j<nlev; j++) {
00195                 cdf2[i*nlev+j] = (double)(hist2[i*nlev+j]) / (double)(hist2[nlev-1]);
00196             }
00197         }
00198 
00199         offset2 = (double)in2->gmin;
00200         factor2 = (double)(nlev-1)/(double)(in2->gmax-in2->gmin);
00201 
00202         /* loop through histograms making the lookup table entries table[index] = cdf2 <inverse> [cdf1[index]] */
00203         for (i=0; i<in1->nbnd; i++) {
00204             for (index=0; index<nlev; index++)  {
00205                 for (index2=0; index2<nlev; index2++)  {
00206                     if (cdf2[i*nlev+index2] == cdf1[i*nlev+index]) {
00207                         table[index] = (PIXEL)index2/(PIXEL)factor2 + offset2;
00208                         break;
00209                     } else if (cdf2[i*nlev+index2] > cdf1[i*nlev+index]) {
00210                         if (index2 == 0) {
00211                             table[index] = offset2;
00212                         } else {
00213                             if ((cdf2[i*nlev+index2] - cdf1[i*nlev+index]) <= (cdf1[i*nlev+index] - cdf2[i*nlev + index2 - 1])) {
00214                                 table[index] = (PIXEL)index2/(PIXEL)factor2 + offset2;
00215                             } else {
00216                                 table[index] = (PIXEL)(index2 - 1)/(PIXEL)factor2 + offset2;
00217                             }
00218                         }
00219                         break;
00220                     }
00221                 }
00222             }
00223 
00224             /* create output image */
00225             for (j=0; j<in1->nlin; j++) {
00226                 for (k=0; k<in1->npix; k++) {
00227                     index = (short)(((double)(in1->data[i][j][k]-offset)*factor) + 0.5);
00228                     (*out)->data[i][j][k] = table[index];
00229                 }
00230                 if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00231             }
00232 
00233             /* output the lookup table */
00234             sprintf(minlbl,"0");
00235             sprintf(maxlbl,"%d",nlev);
00236             PLOT(&table,1,nlev,minlbl,maxlbl,NORM);
00237         }
00238     }
00239 
00240     the_end:
00241     if (table) free(table);
00242     if (count) free(count);
00243     if (hist1) free(hist1);
00244     if (hist2) free(hist2);
00245     if (cdf1)  free(cdf1);
00246     if (cdf2)  free(cdf2);
00247     if (TIMES) TIMING(T_EXIT);
00248 }