Main Page   Data Structures   File List   Data Fields   Globals  

geomcoef.c

Go to the documentation of this 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
00023 GEOMCOEF (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   )
00039 {
00040   register short i, j, k;
00041   char msg[SLEN];
00042   double *a = 0, *ata = 0, *atainv = 0, *xt = 0, *yt = 0, err;
00043 
00044   if (TIMES)
00045     TIMING (T_GEOMCOEF);
00046   if (NAMES)
00047     {
00048       MESSAGE ('I', "");
00049       MESSAGE ('I', "GEOMCOEF");
00050       MESSAGE ('I', "");
00051       sprintf (msg, " Number of control points:  %d", npts);
00052       MESSAGE ('I', msg);
00053       MESSAGE ('I', " Control points:");
00054       MESSAGE ('I',
00055                " Point #        x           y              x'          y'");
00056       for (i = 0; i < npts; i++)
00057         {
00058 #ifndef MAC
00059           sprintf (msg, " %4d%15.4e%12.4e%15.4e%12.4e", i + 1, x[i], y[i],
00060                    xp[i], yp[i]);
00061 #else
00062           sprintf (msg, " %4d%15.4e%12.4e%15.4e%12.4e", i + 1, x[i] + 1.0,
00063                    y[i] + 1.0, xp[i] + 1.0, yp[i] + 1.0);
00064 #endif
00065           MESSAGE ('I', msg);
00066         }
00067       sprintf (msg, " Number of polynomial terms:     %d", nterms);
00068       MESSAGE ('I', msg);
00069       MESSAGE ('I', " ...............");
00070     }
00071 
00072   /* check input */
00073   if (npts < nterms)
00074     {
00075       MESSAGE ('E',
00076                " Number of points must be equal to or greater than number of polynomial terms.");
00077       goto the_end;
00078     }
00079   else if (MINTERMS > nterms)
00080     {
00081       sprintf (msg, " Number of terms must be equal to or greater than %d.",
00082                MINTERMS);
00083       MESSAGE ('E', msg);
00084       goto the_end;
00085     }
00086   else if (nterms > MAXTERMS)
00087     {
00088       sprintf (msg, " Number of terms must be less than or equal to %d.",
00089                MAXTERMS);
00090       MESSAGE ('E', msg);
00091       goto the_end;
00092     }
00093 
00094   /* allocate matrices and vectors */
00095   if (!(a = (double *) malloc (npts * nterms * sizeof (double))))
00096     {
00097       MESSAGE ('E', " Memory request failed.");
00098       goto the_end;
00099     }
00100   if (!(ata = (double *) malloc (nterms * nterms * sizeof (double))))
00101     {
00102       MESSAGE ('E', " Memory request failed.");
00103       goto the_end;
00104     }
00105   if (!(atainv = (double *) malloc (nterms * nterms * sizeof (double))))
00106     {
00107       MESSAGE ('E', " Memory request failed.");
00108       goto the_end;
00109     }
00110   if (!(xt = (double *) malloc (npts * sizeof (double))))
00111     {
00112       MESSAGE ('E', " Memory request failed.");
00113       goto the_end;
00114     }
00115   if (!(yt = (double *) malloc (npts * sizeof (double))))
00116     {
00117       MESSAGE ('E', " Memory request failed.");
00118       goto the_end;
00119     }
00120 
00121   /* compute matrix a */
00122   for (i = 0; i < npts; i++)
00123     {
00124       switch (nterms)
00125         {
00126         case 6:
00127           a[i * nterms + 5] = yp[i] * yp[i];
00128         case 5:
00129           a[i * nterms + 4] = xp[i] * xp[i];
00130         case 4:
00131           a[i * nterms + 3] = xp[i] * yp[i];
00132         case 3:
00133           a[i * nterms + 2] = yp[i];
00134           a[i * nterms + 1] = xp[i];
00135           a[i * nterms + 0] = 1.0;
00136         }
00137     }
00138 
00139   /* compute and invert matrix ata */
00140   for (i = 0; i < nterms; i++)
00141     {
00142       for (j = 0; j < nterms; j++)
00143         {
00144           for (ata[i * nterms + j] = 0.0, k = 0; k < npts; k++)
00145             {
00146               ata[i * nterms + j] += a[k * nterms + i] * a[k * nterms + j];
00147             }
00148         }
00149     }
00150   MATRIX_INVERT (ata, nterms, atainv);
00151 
00152   /* compute transformation coefficients */
00153   for (i = 0; i < nterms; i++)
00154     {
00155       for (xt[i] = yt[i] = 0.0, j = 0; j < npts; j++)
00156         {
00157           xt[i] += a[j * nterms + i] * x[j];
00158           yt[i] += a[j * nterms + i] * y[j];
00159         }
00160     }
00161   for (i = 0; i < nterms; i++)
00162     {
00163       for (xc[i] = yc[i] = 0.0, j = 0; j < nterms; j++)
00164         {
00165           xc[i] += atainv[i * nterms + j] * xt[j];
00166           yc[i] += atainv[i * nterms + j] * yt[j];
00167         }
00168     }
00169 
00170   /* output transformation coefficients */
00171   MESSAGE ('I', "");
00172   MESSAGE ('I', " Transformation Coefficients:");
00173   sprintf (msg, "    x  =   %12.4e                 y  =   %12.4e", xc[0],
00174            yc[0]);
00175   MESSAGE ('I', msg);
00176   sprintf (msg, "          +%12.4e * x'                  +%12.4e * x'", xc[1],
00177            yc[1]);
00178   MESSAGE ('I', msg);
00179   sprintf (msg, "          +%12.4e * y'                  +%12.4e * y'", xc[2],
00180            yc[2]);
00181   MESSAGE ('I', msg);
00182   if (nterms > 3)
00183     {
00184       sprintf (msg,
00185                "          +%12.4e * x' * y'             +%12.4e * x' * y'",
00186                xc[3], yc[3]);
00187       MESSAGE ('I', msg);
00188     }
00189   if (nterms > 4)
00190     {
00191       sprintf (msg,
00192                "          +%12.4e * x' * x'             +%12.4e * x' * x'",
00193                xc[4], yc[4]);
00194       MESSAGE ('I', msg);
00195     }
00196   if (nterms > 5)
00197     {
00198       sprintf (msg,
00199                "          +%12.4e * y' * y'             +%12.4e * y' * y'",
00200                xc[5], yc[5]);
00201       MESSAGE ('I', msg);
00202     }
00203 
00204   /* compute and output control point errors */
00205   MESSAGE ('I', "");
00206   MESSAGE ('I', " Goodness-of-fit at Control Points:");
00207   MESSAGE ('I',
00208            "                    given                    computed                   difference");
00209   MESSAGE ('I',
00210            " Point #        x           y              x           y              dx          dy");
00211   for (err = 0.0, i = 0; i < npts; i++)
00212     {
00213       xt[i] = yt[i] = 0.0;
00214       switch (nterms)
00215         {
00216         case 6:
00217           xt[i] += xc[5] * yp[i] * yp[i];
00218           yt[i] += yc[5] * yp[i] * yp[i];
00219         case 5:
00220           xt[i] += xc[4] * xp[i] * xp[i];
00221           yt[i] += yc[4] * xp[i] * xp[i];
00222         case 4:
00223           xt[i] += xc[3] * xp[i] * yp[i];
00224           yt[i] += yc[3] * xp[i] * yp[i];
00225         case 3:
00226           xt[i] += xc[0] + xc[1] * xp[i] + xc[2] * yp[i];
00227           yt[i] += yc[0] + yc[1] * xp[i] + yc[2] * yp[i];
00228         }
00229       err +=
00230         (xt[i] - x[i]) * (xt[i] - x[i]) + (yt[i] - y[i]) * (yt[i] - y[i]);
00231 #ifndef MAC
00232       sprintf (msg, " %4d%15.4e%12.4e%15.4e%12.4e%15.4e%12.4e", i + 1, x[i],
00233                y[i], xt[i], yt[i], x[i] - xt[i], y[i] - yt[i]);
00234 #else
00235       sprintf (msg, " %4d%15.4e%12.4e%15.4e%12.4e%15.4e%12.4e", i + 1,
00236                x[i] + 1.0, y[i] + 1.0, xt[i] + 1.0, yt[i] + 1.0, x[i] - xt[i],
00237                y[i] - yt[i]);
00238 #endif
00239       MESSAGE ('I', msg);
00240     }
00241   MESSAGE ('I', "");
00242   sprintf (msg, " Root-Mean-Squared Error at Control Points = %10.4e",
00243            sqrt (err / (double) npts));
00244   MESSAGE ('I', msg);
00245 
00246 the_end:
00247   if (a)
00248     free (a);
00249   if (ata)
00250     free (ata);
00251   if (atainv)
00252     free (atainv);
00253   if (xt)
00254     free (xt);
00255   if (yt)
00256     free (yt);
00257   if (TIMES)
00258     TIMING (T_EXIT);
00259 }

Generated on Sun May 18 15:36:10 2003 for tclSadie by doxygen1.3