geomcoef.c Source File

00001 #include        "sadie.h"
00002 
00003 /*-Copyright Information------------------------------------------------------*/
00004 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00005 /*----------------------------------------------------------------------------*/
00006 /*-General Information--------------------------------------------------------*/
00007 /*                                                                            */
00008 /*   This function calculates a least-squares-fit polynomial to a set of      */
00009 /*   spatial control point pairs from two images, a "reference" image with    */
00010 /*   coordinate system (x',y') and a "transformation" image with coordinate   */
00011 /*   system (x,y) to be warped. The polynomial equations are:                 */
00012 /*                                                                            */
00013 /*      x  =    a1                      y  =    b1                            */
00014 /*            + a2 * x'                       + b2 * x'                       */
00015 /*            + a3 * y'                       + b3 * y'                       */
00016 /*            + a4 * x' * y'                  + b4 * x' * y'                  */
00017 /*            + a5 * x' * x'                  + b5 * x' * x'                  */
00018 /*            + a6 * y' * y'                  + b6 * y' * y'                  */
00019 /*                                                                            */
00020 /*----------------------------------------------------------------------------*/
00021 /*-Interface Information------------------------------------------------------*/
00022 void GEOMCOEF (
00023 double *x,      /*  I   Pointer to an array[npts], containing x-coordinates   */
00024                 /*      of the control points in the "reference" image.       */
00025 double *y,      /*  I   Pointer to an array[npts], containing y-coordinates   */
00026                 /*      of the control points in the "reference" image.       */
00027 double *xp,     /*  I   Pointer to an array[npts], containing x-coordinates   */
00028                 /*      of the control points in the "transformation" image.  */
00029 double *yp,     /*  I   Pointer to an array[npts], containing y-coordinates   */
00030                 /*      of the control points in the "transformation" image.  */
00031 short  npts,    /*  I   Number of control point pairs.                        */
00032 short  nterms,  /*  I   Number of polynomial terms.                           */
00033 double *xc,     /*  O   Pointer to an array[nterms], containing the           */
00034                 /*      x-coefficients of the coordinate transformation.      */
00035 double *yc      /*  O   Pointer to an array[nterms], containing               */
00036                 /*      y-coefficients of the coordinate transformation.      */
00037 /*----------------------------------------------------------------------------*/
00038 ) { register short i, j, k;
00039     char   msg[SLEN];
00040     double *a=0, *ata=0, *atainv=0, *xt=0, *yt=0, err;
00041 
00042     if (TIMES) TIMING(T_GEOMCOEF);
00043     if (NAMES) {
00044         MESSAGE('I',"");
00045         MESSAGE('I',"GEOMCOEF");
00046         MESSAGE('I',"");
00047         sprintf(msg," Number of control points:  %d",npts);
00048         MESSAGE('I',msg);
00049         MESSAGE('I'," Control points:");
00050         MESSAGE('I'," Point #        x           y              x'          y'");
00051         for (i=0; i<npts; i++) {
00052 #ifndef MAC
00053             sprintf(msg," %4d%15.4e%12.4e%15.4e%12.4e",i+1,x[i],y[i],xp[i],yp[i]);
00054 #else
00055             sprintf(msg," %4d%15.4e%12.4e%15.4e%12.4e",i+1,x[i]+1.0,y[i]+1.0,xp[i]+1.0,yp[i]+1.0);
00056 #endif
00057             MESSAGE('I',msg);
00058         }
00059         sprintf(msg," Number of polynomial terms:     %d",nterms);
00060         MESSAGE('I',msg);
00061         MESSAGE('I'," ...............");
00062     }
00063 
00064     /* check input */
00065     if (npts < nterms) {
00066         MESSAGE('E'," Number of points must be equal to or greater than number of polynomial terms.");
00067         goto the_end;
00068     } else if (MINTERMS > nterms) {
00069         sprintf(msg," Number of terms must be equal to or greater than %d.",MINTERMS);
00070         MESSAGE('E',msg);
00071         goto the_end;
00072     } else if (nterms > MAXTERMS) {
00073         sprintf(msg," Number of terms must be less than or equal to %d.",MAXTERMS);
00074         MESSAGE('E',msg);
00075         goto the_end;
00076     }
00077 
00078     /* allocate matrices and vectors */
00079     if (!(a=(double *)malloc(npts*nterms*sizeof(double)))) {
00080         MESSAGE('E'," Memory request failed.");
00081         goto the_end;
00082     }
00083     if (!(ata=(double *)malloc(nterms*nterms*sizeof(double)))) {
00084         MESSAGE('E'," Memory request failed.");
00085         goto the_end;
00086     }
00087     if (!(atainv=(double *)malloc(nterms*nterms*sizeof(double)))) {
00088         MESSAGE('E'," Memory request failed.");
00089         goto the_end;
00090     }
00091     if (!(xt=(double *)malloc(npts*sizeof(double)))) {
00092         MESSAGE('E'," Memory request failed.");
00093         goto the_end;
00094     }
00095     if (!(yt=(double *)malloc(npts*sizeof(double)))) {
00096         MESSAGE('E'," Memory request failed.");
00097         goto the_end;
00098     }
00099 
00100     /* compute matrix a */
00101     for (i=0; i<npts; i++) {
00102         switch (nterms) {
00103         case 6:
00104             a[i*nterms+5] = yp[i]*yp[i];
00105         case 5:
00106             a[i*nterms+4] = xp[i]*xp[i];
00107         case 4:
00108             a[i*nterms+3] = xp[i]*yp[i];
00109         case 3:
00110             a[i*nterms+2] = yp[i];
00111             a[i*nterms+1] = xp[i];
00112             a[i*nterms+0] = 1.0;
00113         }
00114     }
00115 
00116     /* compute and invert matrix ata */
00117     for (i=0; i<nterms; i++) {
00118         for (j=0; j<nterms; j++) {
00119             for (ata[i*nterms+j]=0.0,k=0; k<npts; k++) {
00120                 ata[i*nterms+j] += a[k*nterms+i]*a[k*nterms+j];
00121             }
00122         }
00123     }
00124     MATRIX_INVERT(ata,nterms,atainv);
00125 
00126     /* compute transformation coefficients */
00127     for (i=0; i<nterms; i++) {
00128         for (xt[i]=yt[i]=0.0,j=0; j<npts; j++) {
00129             xt[i] += a[j*nterms+i]*x[j];
00130             yt[i] += a[j*nterms+i]*y[j];
00131         }
00132     }
00133     for (i=0; i<nterms; i++) {
00134         for (xc[i]=yc[i]=0.0,j=0; j<nterms; j++) {
00135             xc[i] += atainv[i*nterms+j]*xt[j];
00136             yc[i] += atainv[i*nterms+j]*yt[j];
00137         }
00138     }
00139 
00140     /* output transformation coefficients */
00141     MESSAGE('I',"");
00142     MESSAGE('I'," Transformation Coefficients:");
00143     sprintf(msg,"    x  =   %12.4e                 y  =   %12.4e",xc[0],yc[0]);
00144     MESSAGE('I',msg);
00145     sprintf(msg,"          +%12.4e * x'                  +%12.4e * x'",xc[1],yc[1]);
00146     MESSAGE('I',msg);
00147     sprintf(msg,"          +%12.4e * y'                  +%12.4e * y'",xc[2],yc[2]);
00148     MESSAGE('I',msg);
00149     if (nterms > 3) {
00150         sprintf(msg,"          +%12.4e * x' * y'             +%12.4e * x' * y'",xc[3],yc[3]);
00151         MESSAGE('I',msg);
00152     }
00153     if (nterms > 4) {
00154         sprintf(msg,"          +%12.4e * x' * x'             +%12.4e * x' * x'",xc[4],yc[4]);
00155         MESSAGE('I',msg);
00156     }
00157     if (nterms > 5) {
00158         sprintf(msg,"          +%12.4e * y' * y'             +%12.4e * y' * y'",xc[5],yc[5]);
00159         MESSAGE('I',msg);
00160     }
00161 
00162     /* compute and output control point errors */
00163     MESSAGE('I',"");
00164     MESSAGE('I'," Goodness-of-fit at Control Points:");
00165     MESSAGE('I',"                    given                    computed                   difference");
00166     MESSAGE('I'," Point #        x           y              x           y              dx          dy");
00167     for (err=0.0,i=0; i<npts; i++) {
00168         xt[i] = yt[i] = 0.0;
00169         switch (nterms) {
00170         case 6:
00171             xt[i] += xc[5]*yp[i]*yp[i];
00172             yt[i] += yc[5]*yp[i]*yp[i];
00173         case 5:
00174             xt[i] += xc[4]*xp[i]*xp[i];
00175             yt[i] += yc[4]*xp[i]*xp[i];
00176         case 4:
00177             xt[i] += xc[3]*xp[i]*yp[i];
00178             yt[i] += yc[3]*xp[i]*yp[i];
00179         case 3:
00180             xt[i] += xc[0] + xc[1]*xp[i] + xc[2]*yp[i];
00181             yt[i] += yc[0] + yc[1]*xp[i] + yc[2]*yp[i];
00182         }
00183         err += (xt[i]-x[i])*(xt[i]-x[i]) + (yt[i]-y[i])*(yt[i]-y[i]);
00184 #ifndef MAC
00185         sprintf(msg," %4d%15.4e%12.4e%15.4e%12.4e%15.4e%12.4e",i+1,x[i],y[i],xt[i],yt[i],x[i]-xt[i],y[i]-yt[i]);
00186 #else
00187         sprintf(msg," %4d%15.4e%12.4e%15.4e%12.4e%15.4e%12.4e",i+1,x[i]+1.0,y[i]+1.0,xt[i]+1.0,yt[i]+1.0,x[i]-xt[i],y[i]-yt[i]);
00188 #endif
00189         MESSAGE('I',msg);
00190     }
00191     MESSAGE('I',"");
00192     sprintf(msg," Root-Mean-Squared Error at Control Points = %10.4e",sqrt(err/(double)npts));
00193     MESSAGE('I',msg);
00194 
00195     the_end:
00196     if (a)      free(a);
00197     if (ata)    free(ata);
00198     if (atainv) free(atainv);
00199     if (xt)     free(xt);
00200     if (yt)     free(yt);
00201     if (TIMES)  TIMING(T_EXIT);
00202 }