GRASS Programmer's Manual  6.5.svn(2014)-r66266
georef.c
Go to the documentation of this file.
1
2 /****************************************************************************
3  *
4  * MODULE: imagery library
5  * AUTHOR(S): Original author(s) name(s) unknown - written by CERL
6  * PURPOSE: Image processing library
7  * COPYRIGHT: (C) 1999, 2005 by the GRASS Development Team
8  *
9  * This program is free software under the GNU General Public
10  * License (>=v2). Read the file COPYING that comes with GRASS
11  * for details.
12  *
13  *****************************************************************************/
14
15 #include <grass/config.h>
16 #include <grass/imagery.h>
17 #include <signal.h>
18
19 static int floating_exception;
20 static void catch(int);
21 static double determinant(double, double,
22  double, double, double, double, double, double,
23  double);
24
25 /* find coefficients A,B,C for e2 = A + B*e1 + C*n1
26  * also compute the reverse equations
27  *
28  * return 0 if no points
29  * -1 if not solvable
30  * 1 if ok
31  *
32  * method is least squares.
33  * the least squares problem reduces to solving the following
34  * system of equations for A,B,C
35  *
36  * s0*A + s1*B + s2*C = x0
37  * s1*A + s3*B + s4*C = x1
38  * s2*A + s4*B + s5*C = x2
39  *
40  * use Cramer's rule
41  *
42  * | x0 s1 s2 | | s0 x0 s2 | | s0 s1 x0 |
43  * | x1 s3 s4 | | s1 x1 s4 | | s1 s3 x1 |
44  * | x2 s4 s5 | | s2 x2 s5 | | s2 s4 x2 |
45  * A = ------------ B = ------------ C = ------------
46  * | s0 s1 s2 | | s0 s1 s2 | | s0 s1 s2 |
47  * | s1 s3 s4 | | s1 s3 s4 | | s1 s3 s4 |
48  * | s2 s4 s5 | | s2 s4 s5 | | s2 s4 s5 |
49  *
50  */
51
52 int I_compute_georef_equations(struct Control_Points *cp,
53  double E12[3], double N12[3], double E21[3],
54  double N21[3])
55 {
56  RETSIGTYPE(*sigfpe) (int);
57  double s0, s1, s2, s3, s4, s5;
58  double x0, x1, x2;
59  double det;
60  int i;
61
62
63  s0 = s1 = s2 = s3 = s4 = s5 = 0.0;
64  for (i = 0; i < cp->count; i++) {
65  if (cp->status[i] <= 0)
66  continue;
67  s0 += 1.0;
68  s1 += cp->e1[i];
69  s2 += cp->n1[i];
70  s3 += cp->e1[i] * cp->e1[i];
71  s4 += cp->e1[i] * cp->n1[i];
72  s5 += cp->n1[i] * cp->n1[i];
73  }
74  if (s0 < 0.5)
75  return 0;
76
77  floating_exception = 0;
78  sigfpe = signal(SIGFPE, catch);
79
80  /* eastings */
81  x0 = x1 = x2 = 0.0;
82  for (i = 0; i < cp->count; i++) {
83  if (cp->status[i] <= 0)
84  continue;
85  x0 += cp->e2[i];
86  x1 += cp->e1[i] * cp->e2[i];
87  x2 += cp->n1[i] * cp->e2[i];
88  }
89
90  det = determinant(s0, s1, s2, s1, s3, s4, s2, s4, s5);
91  if (det == 0.0) {
92  signal(SIGFPE, sigfpe);
93  return -1;
94  }
95  E12[0] = determinant(x0, s1, s2, x1, s3, s4, x2, s4, s5) / det;
96  E12[1] = determinant(s0, x0, s2, s1, x1, s4, s2, x2, s5) / det;
97  E12[2] = determinant(s0, s1, x0, s1, s3, x1, s2, s4, x2) / det;
98
99  /* northings */
100  x0 = x1 = x2 = 0.0;
101  for (i = 0; i < cp->count; i++) {
102  if (cp->status[i] <= 0)
103  continue;
104  x0 += cp->n2[i];
105  x1 += cp->e1[i] * cp->n2[i];
106  x2 += cp->n1[i] * cp->n2[i];
107  }
108
109  det = determinant(s0, s1, s2, s1, s3, s4, s2, s4, s5);
110  if (det == 0.0) {
111  signal(SIGFPE, sigfpe);
112  return -1;
113  }
114  N12[0] = determinant(x0, s1, s2, x1, s3, s4, x2, s4, s5) / det;
115  N12[1] = determinant(s0, x0, s2, s1, x1, s4, s2, x2, s5) / det;
116  N12[2] = determinant(s0, s1, x0, s1, s3, x1, s2, s4, x2) / det;
117
118  /* the inverse equations */
119
120  s0 = s1 = s2 = s3 = s4 = s5 = 0.0;
121  for (i = 0; i < cp->count; i++) {
122  if (cp->status[i] <= 0)
123  continue;
124  s0 += 1.0;
125  s1 += cp->e2[i];
126  s2 += cp->n2[i];
127  s3 += cp->e2[i] * cp->e2[i];
128  s4 += cp->e2[i] * cp->n2[i];
129  s5 += cp->n2[i] * cp->n2[i];
130  }
131
132  /* eastings */
133  x0 = x1 = x2 = 0.0;
134  for (i = 0; i < cp->count; i++) {
135  if (cp->status[i] <= 0)
136  continue;
137  x0 += cp->e1[i];
138  x1 += cp->e2[i] * cp->e1[i];
139  x2 += cp->n2[i] * cp->e1[i];
140  }
141
142  det = determinant(s0, s1, s2, s1, s3, s4, s2, s4, s5);
143  if (det == 0.0) {
144  signal(SIGFPE, sigfpe);
145  return -1;
146  }
147  E21[0] = determinant(x0, s1, s2, x1, s3, s4, x2, s4, s5) / det;
148  E21[1] = determinant(s0, x0, s2, s1, x1, s4, s2, x2, s5) / det;
149  E21[2] = determinant(s0, s1, x0, s1, s3, x1, s2, s4, x2) / det;
150
151  /* northings */
152  x0 = x1 = x2 = 0.0;
153  for (i = 0; i < cp->count; i++) {
154  if (cp->status[i] <= 0)
155  continue;
156  x0 += cp->n1[i];
157  x1 += cp->e2[i] * cp->n1[i];
158  x2 += cp->n2[i] * cp->n1[i];
159  }
160
161  det = determinant(s0, s1, s2, s1, s3, s4, s2, s4, s5);
162  if (det == 0.0) {
163  signal(SIGFPE, sigfpe);
164  return -1;
165  }
166  N21[0] = determinant(x0, s1, s2, x1, s3, s4, x2, s4, s5) / det;
167  N21[1] = determinant(s0, x0, s2, s1, x1, s4, s2, x2, s5) / det;
168  N21[2] = determinant(s0, s1, x0, s1, s3, x1, s2, s4, x2) / det;
169
170  signal(SIGFPE, sigfpe);
171  return floating_exception ? -1 : 1;
172 }
173
174 static double determinant(double a, double b, double c, double d, double e,
175  double f, double g, double h, double i)
176 {
177  /* compute determinant of 3x3 matrix
178  * | a b c |
179  * | d e f |
180  * | g h i |
181  */
182  return a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g);
183 }
184
185 static void catch(int n)
186 {
187  floating_exception = 1;
188  signal(n, catch);
189 }
190
191 int I_georef(double e1, double n1,
192  double *e2, double *n2, double E[3], double N[3])
193 {
194  *e2 = E[0] + E[1] * e1 + E[2] * n1;
195  *n2 = N[0] + N[1] * e1 + N[2] * n1;
196
197  return 0;
198 }
float b
Definition: named_colr.c:8
float g
Definition: named_colr.c:8
int
Definition: g3dcolor.c:48
int I_georef(double e1, double n1, double *e2, double *n2, double E[3], double N[3])
Definition: georef.c:191
#define N
Definition: inverse.c:8
int I_compute_georef_equations(struct Control_Points *cp, double E12[3], double N12[3], double E21[3], double N21[3])
Definition: georef.c:52
Definition: spawn.c:69
int n
Definition: dataquad.c:291
tuple h
panel.defaultSize = wx.CheckBox(panel, id = wx.ID_ANY, label = _(&quot;Use default size&quot;)) panel...