GRASS Programmer's Manual  6.4.2(2012)
georef.c
Go to the documentation of this file.
00001 
00002 /****************************************************************************
00003  *
00004  * MODULE:       imagery library
00005  * AUTHOR(S):    Original author(s) name(s) unknown - written by CERL
00006  * PURPOSE:      Image processing library
00007  * COPYRIGHT:    (C) 1999, 2005 by the GRASS Development Team
00008  *
00009  *               This program is free software under the GNU General Public
00010  *               License (>=v2). Read the file COPYING that comes with GRASS
00011  *               for details.
00012  *
00013  *****************************************************************************/
00014 
00015 #include <grass/config.h>
00016 #include <grass/imagery.h>
00017 #include <signal.h>
00018 
00019 static int floating_exception;
00020 static void catch(int);
00021 static double determinant(double, double,
00022                           double, double, double, double, double, double,
00023                           double);
00024 
00025 /* find coefficients A,B,C for e2 = A + B*e1 + C*n1
00026  * also compute the reverse equations
00027  *
00028  * return 0 if no points
00029  *       -1 if not solvable
00030  *        1 if ok
00031  *
00032  * method is least squares.
00033  * the least squares problem reduces to solving the following
00034  * system of equations for A,B,C
00035  *
00036  *   s0*A + s1*B + s2*C = x0
00037  *   s1*A + s3*B + s4*C = x1
00038  *   s2*A + s4*B + s5*C = x2
00039  *
00040  * use Cramer's rule
00041  *
00042  *     | x0 s1 s2 |      | s0 x0 s2 |      | s0 s1 x0 |
00043  *     | x1 s3 s4 |      | s1 x1 s4 |      | s1 s3 x1 |
00044  *     | x2 s4 s5 |      | s2 x2 s5 |      | s2 s4 x2 |
00045  * A = ------------  B = ------------  C = ------------ 
00046  *     | s0 s1 s2 |      | s0 s1 s2 |      | s0 s1 s2 |
00047  *     | s1 s3 s4 |      | s1 s3 s4 |      | s1 s3 s4 |
00048  *     | s2 s4 s5 |      | s2 s4 s5 |      | s2 s4 s5 |
00049  *
00050  */
00051 
00052 int I_compute_georef_equations(struct Control_Points *cp,
00053                                double E12[3], double N12[3], double E21[3],
00054                                double N21[3])
00055 {
00056     RETSIGTYPE(*sigfpe) (int);
00057     double s0, s1, s2, s3, s4, s5;
00058     double x0, x1, x2;
00059     double det;
00060     int i;
00061 
00062 
00063     s0 = s1 = s2 = s3 = s4 = s5 = 0.0;
00064     for (i = 0; i < cp->count; i++) {
00065         if (cp->status[i] <= 0)
00066             continue;
00067         s0 += 1.0;
00068         s1 += cp->e1[i];
00069         s2 += cp->n1[i];
00070         s3 += cp->e1[i] * cp->e1[i];
00071         s4 += cp->e1[i] * cp->n1[i];
00072         s5 += cp->n1[i] * cp->n1[i];
00073     }
00074     if (s0 < 0.5)
00075         return 0;
00076 
00077     floating_exception = 0;
00078     sigfpe = signal(SIGFPE, catch);
00079 
00080     /* eastings */
00081     x0 = x1 = x2 = 0.0;
00082     for (i = 0; i < cp->count; i++) {
00083         if (cp->status[i] <= 0)
00084             continue;
00085         x0 += cp->e2[i];
00086         x1 += cp->e1[i] * cp->e2[i];
00087         x2 += cp->n1[i] * cp->e2[i];
00088     }
00089 
00090     det = determinant(s0, s1, s2, s1, s3, s4, s2, s4, s5);
00091     if (det == 0.0) {
00092         signal(SIGFPE, sigfpe);
00093         return -1;
00094     }
00095     E12[0] = determinant(x0, s1, s2, x1, s3, s4, x2, s4, s5) / det;
00096     E12[1] = determinant(s0, x0, s2, s1, x1, s4, s2, x2, s5) / det;
00097     E12[2] = determinant(s0, s1, x0, s1, s3, x1, s2, s4, x2) / det;
00098 
00099     /* northings */
00100     x0 = x1 = x2 = 0.0;
00101     for (i = 0; i < cp->count; i++) {
00102         if (cp->status[i] <= 0)
00103             continue;
00104         x0 += cp->n2[i];
00105         x1 += cp->e1[i] * cp->n2[i];
00106         x2 += cp->n1[i] * cp->n2[i];
00107     }
00108 
00109     det = determinant(s0, s1, s2, s1, s3, s4, s2, s4, s5);
00110     if (det == 0.0) {
00111         signal(SIGFPE, sigfpe);
00112         return -1;
00113     }
00114     N12[0] = determinant(x0, s1, s2, x1, s3, s4, x2, s4, s5) / det;
00115     N12[1] = determinant(s0, x0, s2, s1, x1, s4, s2, x2, s5) / det;
00116     N12[2] = determinant(s0, s1, x0, s1, s3, x1, s2, s4, x2) / det;
00117 
00118     /* the inverse equations */
00119 
00120     s0 = s1 = s2 = s3 = s4 = s5 = 0.0;
00121     for (i = 0; i < cp->count; i++) {
00122         if (cp->status[i] <= 0)
00123             continue;
00124         s0 += 1.0;
00125         s1 += cp->e2[i];
00126         s2 += cp->n2[i];
00127         s3 += cp->e2[i] * cp->e2[i];
00128         s4 += cp->e2[i] * cp->n2[i];
00129         s5 += cp->n2[i] * cp->n2[i];
00130     }
00131 
00132     /* eastings */
00133     x0 = x1 = x2 = 0.0;
00134     for (i = 0; i < cp->count; i++) {
00135         if (cp->status[i] <= 0)
00136             continue;
00137         x0 += cp->e1[i];
00138         x1 += cp->e2[i] * cp->e1[i];
00139         x2 += cp->n2[i] * cp->e1[i];
00140     }
00141 
00142     det = determinant(s0, s1, s2, s1, s3, s4, s2, s4, s5);
00143     if (det == 0.0) {
00144         signal(SIGFPE, sigfpe);
00145         return -1;
00146     }
00147     E21[0] = determinant(x0, s1, s2, x1, s3, s4, x2, s4, s5) / det;
00148     E21[1] = determinant(s0, x0, s2, s1, x1, s4, s2, x2, s5) / det;
00149     E21[2] = determinant(s0, s1, x0, s1, s3, x1, s2, s4, x2) / det;
00150 
00151     /* northings */
00152     x0 = x1 = x2 = 0.0;
00153     for (i = 0; i < cp->count; i++) {
00154         if (cp->status[i] <= 0)
00155             continue;
00156         x0 += cp->n1[i];
00157         x1 += cp->e2[i] * cp->n1[i];
00158         x2 += cp->n2[i] * cp->n1[i];
00159     }
00160 
00161     det = determinant(s0, s1, s2, s1, s3, s4, s2, s4, s5);
00162     if (det == 0.0) {
00163         signal(SIGFPE, sigfpe);
00164         return -1;
00165     }
00166     N21[0] = determinant(x0, s1, s2, x1, s3, s4, x2, s4, s5) / det;
00167     N21[1] = determinant(s0, x0, s2, s1, x1, s4, s2, x2, s5) / det;
00168     N21[2] = determinant(s0, s1, x0, s1, s3, x1, s2, s4, x2) / det;
00169 
00170     signal(SIGFPE, sigfpe);
00171     return floating_exception ? -1 : 1;
00172 }
00173 
00174 static double determinant(double a, double b, double c, double d, double e,
00175                           double f, double g, double h, double i)
00176 {
00177     /* compute determinant of 3x3 matrix
00178      *     | a b c |
00179      *     | d e f |
00180      *     | g h i |
00181      */
00182     return a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g);
00183 }
00184 
00185 static void catch(int n)
00186 {
00187     floating_exception = 1;
00188     signal(n, catch);
00189 }
00190 
00191 int I_georef(double e1, double n1,
00192              double *e2, double *n2, double E[3], double N[3])
00193 {
00194     *e2 = E[0] + E[1] * e1 + E[2] * n1;
00195     *n2 = N[0] + N[1] * e1 + N[2] * n1;
00196 
00197     return 0;
00198 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines