GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-112dd97adf
georef.c
Go to the documentation of this file.
1 /****************************************************************************
2  *
3  * MODULE: imagery library
4  * AUTHOR(S): Original author(s) name(s) unknown - written by CERL
5  * Written By: Brian J. Buckley
6  *
7  * At: The Center for Remote Sensing
8  * Michigan State University
9  *
10  * PURPOSE: Image processing library
11  * COPYRIGHT: (C) 1999, 2005 by the GRASS Development Team
12  *
13  * This program is free software under the GNU General Public
14  * License (>=v2). Read the file COPYING that comes with GRASS
15  * for details.
16  *
17  *****************************************************************************/
18 
19 /*
20  * Written: 12/19/91
21  *
22  * Last Update: 12/26/91 Brian J. Buckley
23  * Last Update: 1/24/92 Brian J. Buckley
24  * Added printout of trnfile. Triggered by BDEBUG.
25  * Last Update: 1/27/92 Brian J. Buckley
26  * Fixed bug so that only the active control points were used.
27  *
28  */
29 
30 #include <stdlib.h>
31 #include <math.h>
32 #include <grass/gis.h>
33 #include <grass/imagery.h>
34 #include <signal.h>
35 
36 /* STRUCTURE FOR USE INTERNALLY WITH THESE FUNCTIONS. THESE FUNCTIONS EXPECT
37  SQUARE MATRICES SO ONLY ONE VARIABLE IS GIVEN (N) FOR THE MATRIX SIZE */
38 
39 struct MATRIX {
40  int n; /* SIZE OF THIS MATRIX (N x N) */
41  double *v;
42 };
43 
44 /* CALCULATE OFFSET INTO ARRAY BASED ON R/C */
45 
46 #define M(row, col) m->v[(((row)-1) * (m->n)) + (col)-1]
47 
48 #define MSUCCESS 1 /* SUCCESS */
49 #define MNPTERR 0 /* NOT ENOUGH POINTS */
50 #define MUNSOLVABLE -1 /* NOT SOLVABLE */
51 #define MMEMERR -2 /* NOT ENOUGH MEMORY */
52 #define MPARMERR -3 /* PARAMETER ERROR */
53 #define MINTERR -4 /* INTERNAL ERROR */
54 
55 #define MAXORDER 3 /* HIGHEST SUPPORTED ORDER OF TRANSFORMATION */
56 
57 /***********************************************************************
58 
59  FUNCTION PROTOTYPES FOR STATIC (INTERNAL) FUNCTIONS
60 
61 ************************************************************************/
62 
63 static int calccoef(struct Control_Points *, double *, double *, int);
64 static int calcls(struct Control_Points *, struct MATRIX *, double *, double *,
65  double *, double *);
66 static int exactdet(struct Control_Points *, struct MATRIX *, double *,
67  double *, double *, double *);
68 static int solvemat(struct MATRIX *, double *, double *, double *, double *);
69 static double term(int, double, double);
70 
71 /***********************************************************************
72 
73  TRANSFORM A SINGLE COORDINATE PAIR.
74 
75 ************************************************************************/
76 
77 int I_georef(double e1, /* EASTING TO BE TRANSFORMED */
78  double n1, /* NORTHING TO BE TRANSFORMED */
79  double *e, /* EASTING, TRANSFORMED */
80  double *n, /* NORTHING, TRANSFORMED */
81  double E[], /* EASTING COEFFICIENTS */
82  double N[], /* NORTHING COEFFICIENTS */
83  int order /* ORDER OF TRANSFORMATION TO BE PERFORMED, MUST MATCH
84  THE ORDER USED TO CALCULATE THE COEFFICIENTS */
85 )
86 {
87  double e3, e2n, en2, n3, e2, en, n2;
88 
89  switch (order) {
90  case 1:
91  *e = E[0] + E[1] * e1 + E[2] * n1;
92  *n = N[0] + N[1] * e1 + N[2] * n1;
93  break;
94 
95  case 2:
96  e2 = e1 * e1;
97  n2 = n1 * n1;
98  en = e1 * n1;
99 
100  *e = E[0] + E[1] * e1 + E[2] * n1 + E[3] * e2 + E[4] * en + E[5] * n2;
101  *n = N[0] + N[1] * e1 + N[2] * n1 + N[3] * e2 + N[4] * en + N[5] * n2;
102  break;
103 
104  case 3:
105  e2 = e1 * e1;
106  en = e1 * n1;
107  n2 = n1 * n1;
108  e3 = e1 * e2;
109  e2n = e2 * n1;
110  en2 = e1 * n2;
111  n3 = n1 * n2;
112 
113  *e = E[0] + E[1] * e1 + E[2] * n1 + E[3] * e2 + E[4] * en + E[5] * n2 +
114  E[6] * e3 + E[7] * e2n + E[8] * en2 + E[9] * n3;
115  *n = N[0] + N[1] * e1 + N[2] * n1 + N[3] * e2 + N[4] * en + N[5] * n2 +
116  N[6] * e3 + N[7] * e2n + N[8] * en2 + N[9] * n3;
117  break;
118 
119  default:
120  return MPARMERR;
121  }
122 
123  return MSUCCESS;
124 }
125 
126 /***********************************************************************
127 
128  COMPUTE THE FORWARD AND BACKWARD GEOREFFERENCING COEFFICIENTS
129  BASED ON A SET OF CONTROL POINTS
130 
131 ************************************************************************/
132 
133 int I_compute_georef_equations(struct Control_Points *cp, double E12[],
134  double N12[], double E21[], double N21[],
135  int order)
136 {
137  double *tempptr;
138  int status;
139 
140  if (order < 1 || order > MAXORDER)
141  return MPARMERR;
142 
143  /* CALCULATE THE FORWARD TRANSFORMATION COEFFICIENTS */
144 
145  status = calccoef(cp, E12, N12, order);
146 
147  if (status != MSUCCESS)
148  return status;
149 
150  /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS */
151 
152  tempptr = cp->e1;
153  cp->e1 = cp->e2;
154  cp->e2 = tempptr;
155  tempptr = cp->n1;
156  cp->n1 = cp->n2;
157  cp->n2 = tempptr;
158 
159  /* CALCULATE THE BACKWARD TRANSFORMATION COEFFICIENTS */
160 
161  status = calccoef(cp, E21, N21, order);
162 
163  /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS BACK */
164 
165  tempptr = cp->e1;
166  cp->e1 = cp->e2;
167  cp->e2 = tempptr;
168  tempptr = cp->n1;
169  cp->n1 = cp->n2;
170  cp->n2 = tempptr;
171 
172  return status;
173 }
174 
175 /***********************************************************************
176 
177  COMPUTE THE GEOREFFERENCING COEFFICIENTS
178  BASED ON A SET OF CONTROL POINTS
179 
180 ************************************************************************/
181 
182 static int calccoef(struct Control_Points *cp, double E[], double N[],
183  int order)
184 {
185  struct MATRIX m;
186  double *a;
187  double *b;
188  int numactive; /* NUMBER OF ACTIVE CONTROL POINTS */
189  int status, i;
190 
191  /* CALCULATE THE NUMBER OF VALID CONTROL POINTS */
192 
193  for (i = numactive = 0; i < cp->count; i++) {
194  if (cp->status[i] > 0)
195  numactive++;
196  }
197 
198  /* CALCULATE THE MINIMUM NUMBER OF CONTROL POINTS NEEDED TO DETERMINE
199  A TRANSFORMATION OF THIS ORDER */
200 
201  m.n = ((order + 1) * (order + 2)) / 2;
202 
203  if (numactive < m.n)
204  return MNPTERR;
205 
206  /* INITIALIZE MATRIX */
207 
208  m.v = G_calloc(m.n * m.n, sizeof(double));
209  a = G_calloc(m.n, sizeof(double));
210  b = G_calloc(m.n, sizeof(double));
211 
212  if (numactive == m.n)
213  status = exactdet(cp, &m, a, b, E, N);
214  else
215  status = calcls(cp, &m, a, b, E, N);
216 
217  G_free(m.v);
218  G_free(a);
219  G_free(b);
220 
221  return status;
222 }
223 
224 /***********************************************************************
225 
226  CALCULATE THE TRANSFORMATION COEFFICIENTS WITH EXACTLY THE MINIMUM
227  NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION.
228 
229 ************************************************************************/
230 
231 static int exactdet(struct Control_Points *cp, struct MATRIX *m, double a[],
232  double b[], double E[], /* EASTING COEFFICIENTS */
233  double N[] /* NORTHING COEFFICIENTS */
234 )
235 {
236  int pntnow, currow, j;
237 
238  currow = 1;
239  for (pntnow = 0; pntnow < cp->count; pntnow++) {
240  if (cp->status[pntnow] > 0) {
241  /* POPULATE MATRIX M */
242 
243  for (j = 1; j <= m->n; j++)
244  M(currow, j) = term(j, cp->e1[pntnow], cp->n1[pntnow]);
245 
246  /* POPULATE MATRIX A AND B */
247 
248  a[currow - 1] = cp->e2[pntnow];
249  b[currow - 1] = cp->n2[pntnow];
250 
251  currow++;
252  }
253  }
254 
255  if (currow - 1 != m->n)
256  return MINTERR;
257 
258  return solvemat(m, a, b, E, N);
259 }
260 
261 /***********************************************************************
262 
263  CALCULATE THE TRANSFORMATION COEFFICIENTS WITH MORE THAN THE MINIMUM
264  NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION. THIS
265  ROUTINE USES THE LEAST SQUARES METHOD TO COMPUTE THE COEFFICIENTS.
266 
267 ************************************************************************/
268 
269 static int calcls(struct Control_Points *cp, struct MATRIX *m, double a[],
270  double b[], double E[], /* EASTING COEFFICIENTS */
271  double N[] /* NORTHING COEFFICIENTS */
272 )
273 {
274  int i, j, n, numactive = 0;
275 
276  /* INITIALIZE THE UPPER HALF OF THE MATRIX AND THE TWO COLUMN VECTORS */
277 
278  for (i = 1; i <= m->n; i++) {
279  for (j = i; j <= m->n; j++)
280  M(i, j) = 0.0;
281  a[i - 1] = b[i - 1] = 0.0;
282  }
283 
284  /* SUM THE UPPER HALF OF THE MATRIX AND THE COLUMN VECTORS ACCORDING TO
285  THE LEAST SQUARES METHOD OF SOLVING OVER DETERMINED SYSTEMS */
286 
287  for (n = 0; n < cp->count; n++) {
288  if (cp->status[n] > 0) {
289  numactive++;
290  for (i = 1; i <= m->n; i++) {
291  for (j = i; j <= m->n; j++)
292  M(i, j) += term(i, cp->e1[n], cp->n1[n]) *
293  term(j, cp->e1[n], cp->n1[n]);
294 
295  a[i - 1] += cp->e2[n] * term(i, cp->e1[n], cp->n1[n]);
296  b[i - 1] += cp->n2[n] * term(i, cp->e1[n], cp->n1[n]);
297  }
298  }
299  }
300 
301  if (numactive <= m->n)
302  return MINTERR;
303 
304  /* TRANSPOSE VALUES IN UPPER HALF OF M TO OTHER HALF */
305 
306  for (i = 2; i <= m->n; i++)
307  for (j = 1; j < i; j++)
308  M(i, j) = M(j, i);
309 
310  return solvemat(m, a, b, E, N);
311 }
312 
313 /***********************************************************************
314 
315  CALCULATE THE X/Y TERM BASED ON THE TERM NUMBER
316 
317  ORDER\TERM 1 2 3 4 5 6 7 8 9 10
318  1 e0n0 e1n0 e0n1
319  2 e0n0 e1n0 e0n1 e2n0 e1n1 e0n2
320  3 e0n0 e1n0 e0n1 e2n0 e1n1 e0n2 e3n0 e2n1 e1n2 e0n3
321 
322 ************************************************************************/
323 
324 static double term(int term, double e, double n)
325 {
326  switch (term) {
327  case 1:
328  return 1.0;
329  case 2:
330  return e;
331  case 3:
332  return n;
333  case 4:
334  return e * e;
335  case 5:
336  return e * n;
337  case 6:
338  return n * n;
339  case 7:
340  return e * e * e;
341  case 8:
342  return e * e * n;
343  case 9:
344  return e * n * n;
345  case 10:
346  return n * n * n;
347  }
348 
349  return 0.0;
350 }
351 
352 /***********************************************************************
353 
354  SOLVE FOR THE 'E' AND 'N' COEFFICIENTS BY USING A SOMEWHAT MODIFIED
355  GAUSSIAN ELIMINATION METHOD.
356 
357  | M11 M12 ... M1n | | E0 | | a0 |
358  | M21 M22 ... M2n | | E1 | = | a1 |
359  | . . . . | | . | | . |
360  | Mn1 Mn2 ... Mnn | | En-1 | | an-1 |
361 
362  and
363 
364  | M11 M12 ... M1n | | N0 | | b0 |
365  | M21 M22 ... M2n | | N1 | = | b1 |
366  | . . . . | | . | | . |
367  | Mn1 Mn2 ... Mnn | | Nn-1 | | bn-1 |
368 
369 ************************************************************************/
370 
371 static int solvemat(struct MATRIX *m, double a[], double b[], double E[],
372  double N[])
373 {
374  int i, j, i2, j2, imark;
375  double factor, temp;
376  double pivot; /* ACTUAL VALUE OF THE LARGEST PIVOT CANDIDATE */
377 
378  for (i = 1; i <= m->n; i++) {
379  j = i;
380 
381  /* find row with largest magnitude value for pivot value */
382 
383  pivot = M(i, j);
384  imark = i;
385  for (i2 = i + 1; i2 <= m->n; i2++) {
386  temp = fabs(M(i2, j));
387  if (temp > fabs(pivot)) {
388  pivot = M(i2, j);
389  imark = i2;
390  }
391  }
392 
393  /* if the pivot is very small then the points are nearly co-linear */
394  /* co-linear points result in an undefined matrix, and nearly */
395  /* co-linear points results in a solution with rounding error */
396 
397  if (pivot == 0.0)
398  return MUNSOLVABLE;
399 
400  /* if row with highest pivot is not the current row, switch them */
401 
402  if (imark != i) {
403  for (j2 = 1; j2 <= m->n; j2++) {
404  temp = M(imark, j2);
405  M(imark, j2) = M(i, j2);
406  M(i, j2) = temp;
407  }
408 
409  temp = a[imark - 1];
410  a[imark - 1] = a[i - 1];
411  a[i - 1] = temp;
412 
413  temp = b[imark - 1];
414  b[imark - 1] = b[i - 1];
415  b[i - 1] = temp;
416  }
417 
418  /* compute zeros above and below the pivot, and compute
419  values for the rest of the row as well */
420 
421  for (i2 = 1; i2 <= m->n; i2++) {
422  if (i2 != i) {
423  factor = M(i2, j) / pivot;
424  for (j2 = j; j2 <= m->n; j2++)
425  M(i2, j2) -= factor * M(i, j2);
426  a[i2 - 1] -= factor * a[i - 1];
427  b[i2 - 1] -= factor * b[i - 1];
428  }
429  }
430  }
431 
432  /* SINCE ALL OTHER VALUES IN THE MATRIX ARE ZERO NOW, CALCULATE THE
433  COEFFICIENTS BY DIVIDING THE COLUMN VECTORS BY THE DIAGONAL VALUES. */
434 
435  for (i = 1; i <= m->n; i++) {
436  E[i - 1] = a[i - 1] / M(i, i);
437  N[i - 1] = b[i - 1] / M(i, i);
438  }
439 
440  return MSUCCESS;
441 }
int order(int i_x, int i_y, int yNum)
Definition: InterpSpline.c:53
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
#define G_calloc(m, n)
Definition: defs/gis.h:95
#define N
Definition: e_intersect.c:926
int I_georef(double e1, double n1, double *e, double *n, double E[], double N[], int order)
Definition: georef.c:77
#define MAXORDER
Definition: georef.c:55
#define MNPTERR
Definition: georef.c:49
int I_compute_georef_equations(struct Control_Points *cp, double E12[], double N12[], double E21[], double N21[], int order)
Definition: georef.c:133
#define MINTERR
Definition: georef.c:53
#define MSUCCESS
Definition: georef.c:48
#define MPARMERR
Definition: georef.c:52
#define MUNSOLVABLE
Definition: georef.c:50
#define M(row, col)
Definition: georef.c:46
double b
Definition: r_raster.c:39
int * status
Definition: imagery.h:44
double * e2
Definition: imagery.h:41
double * e1
Definition: imagery.h:38
double * n2
Definition: imagery.h:42
double * n1
Definition: imagery.h:39