transform.c

Go to the documentation of this file.
00001 
00020 /****************************************************************
00021 note: uses sqrt() from math library
00022 *****************************************************************
00023 Points from one system may be converted into the second by
00024 use of one of the two equation routines.
00025 
00026 transform_a_into_b (ax,ay,bx,by)
00027 
00028     double ax,ay;            input point from system a
00029     double *bx,*by;          resultant point in system b
00030 
00031 transform_b_into_a (bx,by,ax,ay)
00032 
00033     double bx,by;            input point from system b
00034     double *ax,*ay;          resultant point in system a
00035 *****************************************************************
00036 Residual analysis on the equation can be run to test how well
00037 the equations work.  Either test how well b is predicted by a
00038 or vice versa.
00039 
00040 residuals_a_predicts_b (ax,ay,bx,by,use,n,residuals,rms)
00041 residuals_b_predicts_a (ax,ay,bx,by,use,n,residuals,rms)
00042 
00043     double ax[], ay[];       coordinate from system a
00044     double bx[], by[];       coordinate from system b
00045     char   use[];            use point flags
00046     int n;                   number of points in ax,ay,bx,by
00047     double residual[]        residual error for each point
00048     double *rms;             overall root mean square error
00049 ****************************************************************/
00050 
00051 #include <stdio.h>
00052 #include <math.h>
00053 #include <grass/libtrans.h>
00054 
00055 /* the coefficients */
00056 static double A0, A1, A2, A3, A4, A5;
00057 static double B0, B1, B2, B3, B4, B5;
00058 
00059 /* function prototypes */
00060 static int resid(double *, double *, double *, double *, int *, int, double *,
00061                  double *, int);
00062 
00063 
00090 int compute_transformation_coef(double ax[], double ay[], double bx[],
00091                                 double by[], int *use, int n)
00092 {
00093     int i;
00094     int j;
00095     int count;
00096     double aa[3];
00097     double aar[3];
00098     double bb[3];
00099     double bbr[3];
00100 
00101     double cc[3][3];
00102     double x;
00103 
00104     count = 0;
00105     for (i = 0; i < n; i++)
00106         if (use[i])
00107             count++;
00108     if (count < 4)
00109         return -2;              /* must have at least 4 points */
00110 
00111     for (i = 0; i < 3; i++) {
00112         aa[i] = bb[i] = 0.0;
00113 
00114         for (j = 0; j < 3; j++)
00115             cc[i][j] = 0.0;
00116     }
00117 
00118     for (i = 0; i < n; i++) {
00119         if (!use[i])
00120             continue;           /* skip this point */
00121         cc[0][0] += 1;
00122         cc[0][1] += bx[i];
00123         cc[0][2] += by[i];
00124 
00125         cc[1][1] += bx[i] * bx[i];
00126         cc[1][2] += bx[i] * by[i];
00127         cc[2][2] += by[i] * by[i];
00128 
00129         aa[0] += ay[i];
00130         aa[1] += ay[i] * bx[i];
00131         aa[2] += ay[i] * by[i];
00132 
00133         bb[0] += ax[i];
00134         bb[1] += ax[i] * bx[i];
00135         bb[2] += ax[i] * by[i];
00136     }
00137 
00138     cc[1][0] = cc[0][1];
00139     cc[2][0] = cc[0][2];
00140     cc[2][1] = cc[1][2];
00141 
00142     /* aa and bb are solved */
00143     if (inverse(cc) < 0)
00144         return (-1);
00145     if (m_mult(cc, aa, aar) < 0 || m_mult(cc, bb, bbr) < 0)
00146         return (-1);
00147 
00148     /* the equation coefficients */
00149     B0 = aar[0];
00150     B1 = aar[1];
00151     B2 = aar[2];
00152 
00153     B3 = bbr[0];
00154     B4 = bbr[1];
00155     B5 = bbr[2];
00156 
00157     /* the inverse equation */
00158     x = B2 * B4 - B1 * B5;
00159 
00160     if (!x)
00161         return (-1);
00162 
00163     A0 = (B1 * B3 - B0 * B4) / x;
00164     A1 = -B1 / x;
00165     A2 = B4 / x;
00166     A3 = (B0 * B5 - B2 * B3) / x;
00167     A4 = B2 / x;
00168     A5 = -B5 / x;
00169 
00170     return 1;
00171 }
00172 
00173 
00174 int transform_a_into_b(double ax, double ay, double *bx, double *by)
00175 {
00176     *by = A0 + A1 * ax + A2 * ay;
00177     *bx = A3 + A4 * ax + A5 * ay;
00178 
00179     return 0;
00180 }
00181 
00182 
00183 int transform_b_into_a(double bx, double by, double *ax, double *ay)
00184 {
00185     *ay = B0 + B1 * bx + B2 * by;
00186     *ax = B3 + B4 * bx + B5 * by;
00187 
00188     return 0;
00189 }
00190 
00191 /**************************************************************
00192 These routines are internal to this source code
00193 
00194 solve (a, b)
00195     double a[3][3]
00196     double b[3]
00197 
00198     equation solver used by compute_transformation_coef()
00199 **************************************************************/
00200 
00201 /*  #define abs(xx) (xx >= 0 ? xx : -xx)  */
00202 /*      #define N 3  */
00203 
00204 
00205 int residuals_a_predicts_b(double ax[], double ay[], double bx[], double by[],
00206                            int use[], int n, double residuals[], double *rms)
00207 {
00208     resid(ax, ay, bx, by, use, n, residuals, rms, 1);
00209 
00210     return 0;
00211 }
00212 
00213 
00214 int residuals_b_predicts_a(double ax[], double ay[], double bx[], double by[],
00215                            int use[], int n, double residuals[], double *rms)
00216 {
00217     resid(ax, ay, bx, by, use, n, residuals, rms, 0);
00218 
00219     return 0;
00220 }
00221 
00222 
00231 int print_transform_matrix(void)
00232 {
00233     fprintf(stdout, "\nTransformation Matrix\n");
00234     fprintf(stdout, "| xoff a b |\n");
00235     fprintf(stdout, "| yoff d e |\n");
00236     fprintf(stdout, "-------------------------------------------\n");
00237     fprintf(stdout, "%f %f %f \n", -B3, B2, -B5);
00238     fprintf(stdout, "%f %f %f \n", -B0, -B1, B4);
00239     fprintf(stdout, "-------------------------------------------\n");
00240 
00241     return 1;
00242 }
00243 
00244 
00245 static int resid(double ax[], double ay[], double bx[], double by[],
00246                  int use[], int n, double residuals[], double *rms, int atob)
00247 {
00248     double x, y;
00249     int i;
00250     int count;
00251     double sum;
00252     double delta;
00253     double dx, dy;
00254 
00255     count = 0;
00256     sum = 0.0;
00257     for (i = 0; i < n; i++) {
00258         if (!use[i])
00259             continue;
00260 
00261         count++;
00262         if (atob) {
00263             transform_a_into_b(ax[i], ay[i], &x, &y);
00264             dx = x - bx[i];
00265             dy = y - by[i];
00266         }
00267         else {
00268             transform_b_into_a(bx[i], by[i], &x, &y);
00269             dx = x - ax[i];
00270             dy = y - ay[i];
00271         }
00272 
00273         delta = dx * dx + dy * dy;
00274         residuals[i] = sqrt(delta);
00275         sum += delta;
00276     }
00277     *rms = sqrt(sum / count);
00278 
00279     return 0;
00280 }

Generated on Thu Jul 16 13:21:18 2009 for GRASS Programmer's Manual by  doxygen 1.5.6