GRASS Programmer's Manual
6.4.2(2012)
|
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 }