00001 #include "sadie.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 void
00023 GEOMCOEF (double *x,
00024
00025 double *y,
00026
00027 double *xp,
00028
00029 double *yp,
00030
00031 short npts,
00032 short nterms,
00033 double *xc,
00034
00035 double *yc
00036
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
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
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
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
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
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
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
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 }