GRASS GIS 8 Programmer's Manual  8.2.1dev(2022)-f67d8e515
georef_tps.c
Go to the documentation of this file.
1 
2 /****************************************************************************
3  *
4  * MODULE: imagery library
5  * AUTHOR(S): Markus Metz
6  *
7  * PURPOSE: Image processing library
8  * COPYRIGHT: (C) 2013 by the GRASS Development Team
9  *
10  * This program is free software under the GNU General Public
11  * License (>=v2). Read the file COPYING that comes with GRASS
12  * for details.
13  *
14  *****************************************************************************/
15 
16 #include <stdlib.h>
17 #include <math.h>
18 #include <grass/gis.h>
19 #include <grass/imagery.h>
20 #include <grass/glocale.h>
21 #include <signal.h>
22 
23 /* STRUCTURE FOR USE INTERNALLY WITH THESE FUNCTIONS. THESE FUNCTIONS EXPECT
24  SQUARE MATRICES SO ONLY ONE VARIABLE IS GIVEN (N) FOR THE MATRIX SIZE */
25 
26 struct MATRIX
27 {
28  int n; /* SIZE OF THIS MATRIX (N x N) */
29  double *v;
30 };
31 
32 /* CALCULATE OFFSET INTO ARRAY BASED ON R/C */
33 
34 #define M(row,col) m->v[(((row)-1)*(m->n))+(col)-1]
35 
36 #define MSUCCESS 1 /* SUCCESS */
37 #define MNPTERR 0 /* NOT ENOUGH POINTS */
38 #define MUNSOLVABLE -1 /* NOT SOLVABLE */
39 #define MMEMERR -2 /* NOT ENOUGH MEMORY */
40 #define MPARMERR -3 /* PARAMETER ERROR */
41 #define MINTERR -4 /* INTERNAL ERROR */
42 
43 #define MAXORDER 3 /* HIGHEST SUPPORTED ORDER OF TRANSFORMATION */
44 
45 /***********************************************************************
46 
47  FUNCTION PROTOTYPES FOR STATIC (INTERNAL) FUNCTIONS
48 
49 ************************************************************************/
50 
51 static int calccoef(struct Control_Points *, double **, double **);
52 static int calcls(struct Control_Points *, struct MATRIX *, double *,
53  double *, double *, double *);
54 
55 static double tps_base_func(const double x1, const double y1,
56  const double x2, const double y2);
57 static int solvemat(struct MATRIX *, double *, double *, double *, double *);
58 
59 /***********************************************************************
60 
61  TRANSFORM A SINGLE COORDINATE PAIR.
62 
63 ************************************************************************/
64 
65 int I_georef_tps(double e1, /* EASTING TO BE TRANSFORMED */
66  double n1, /* NORTHING TO BE TRANSFORMED */
67  double *e, /* EASTING, TRANSFORMED */
68  double *n, /* NORTHING, TRANSFORMED */
69  double *E, /* EASTING COEFFICIENTS */
70  double *N, /* NORTHING COEFFICIENTS */
71  struct Control_Points *cp,
72  int fwd
73  )
74 {
75  int i, j;
76  double dist, *pe, *pn;
77 
78  if (fwd) {
79  pe = cp->e1;
80  pn = cp->n1;
81  }
82  else {
83  pe = cp->e2;
84  pn = cp->n2;
85  }
86 
87  /* global affine (1st order poly) */
88  *e = E[0] + e1 * E[1] + n1 * E[2];
89  *n = N[0] + e1 * N[1] + n1 * N[2];
90 
91 
92  for (i = 0, j = 0; i < cp->count; i++) {
93  if (cp->status[i] > 0) {
94 
95  dist = tps_base_func(e1, n1, pe[i], pn[i]);
96 
97  *e += E[j + 3] * dist;
98  *n += N[j + 3] * dist;
99  j++;
100  }
101  }
102 
103  return MSUCCESS;
104 }
105 
106 /***********************************************************************
107 
108  COMPUTE THE FORWARD AND BACKWARD GEOREFFERENCING COEFFICIENTS
109  BASED ON A SET OF CONTROL POINTS
110 
111 ************************************************************************/
112 
114  double **E12tps, double **N12tps,
115  double **E21tps, double **N21tps)
116 {
117  double *tempptr;
118  int numactive; /* NUMBER OF ACTIVE CONTROL POINTS */
119  int status, i;
120  double xmax, xmin, ymax, ymin;
121  double delx, dely;
122  double xx, yy;
123  double sumx, sumy, sumx2, sumy2, sumxy;
124  double SSxx, SSyy, SSxy;
125 
126  /* CALCULATE THE NUMBER OF VALID CONTROL POINTS */
127 
128  for (i = numactive = 0; i < cp->count; i++) {
129  if (cp->status[i] > 0)
130  numactive++;
131  }
132 
133  if (numactive < 3)
134  return MNPTERR;
135 
136  if (numactive > 100000) /* arbitrary, admittedly */
137  return MNPTERR;
138 
139  xmin = xmax = cp->e1[0];
140  ymin = ymax = cp->n1[0];
141 
142  sumx = sumy = sumx2 = sumy2 = sumxy = 0.0;
143 
144  for (i = 0; i < cp->count; i++ ) {
145  if (cp->status[i] > 0) {
146  xx = cp->e1[i];
147  yy = cp->n1[i];
148 
149  xmax = MAX(xmax, xx);
150  xmin = MIN(xmin, xx);
151  ymax = MAX(ymax, yy);
152  ymin = MIN(ymin, yy);
153 
154  sumx += xx;
155  sumx2 += xx * xx;
156  sumy += yy;
157  sumy2 += yy * yy;
158  sumxy += xx * yy;
159  }
160  }
161 
162  delx = xmax - xmin;
163  dely = ymax - ymin;
164 
165  SSxx = sumx2 - sumx * sumx / numactive;
166  SSyy = sumy2 - sumy * sumy / numactive;
167  SSxy = sumxy - sumx * sumy / numactive;
168 
169  if (delx < 0.001 * dely || dely < 0.001 * delx ||
170  fabs(SSxy * SSxy / (SSxx * SSyy)) > 0.99) {
171  /* points are colinear */
172  return MUNSOLVABLE;
173  }
174 
175  xmin = xmax = cp->e2[0];
176  ymin = ymax = cp->n2[0];
177 
178  sumx = sumy = sumx2 = sumy2 = sumxy = 0.0;
179  for (i = 0; i < cp->count; i++ ) {
180  if (cp->status[i] > 0) {
181  xx = cp->e2[i];
182  yy = cp->n2[i];
183 
184  xmax = MAX(xmax, xx);
185  xmin = MIN(xmin, xx);
186  ymax = MAX(ymax, yy);
187  ymin = MIN(ymin, yy);
188 
189  sumx += xx;
190  sumx2 += xx * xx;
191  sumy += yy;
192  sumy2 += yy * yy;
193  sumxy += xx * yy;
194  }
195  }
196 
197  delx = xmax - xmin;
198  dely = ymax - ymin;
199 
200  SSxx = sumx2 - sumx * sumx / numactive;
201  SSyy = sumy2 - sumy * sumy / numactive;
202  SSxy = sumxy - sumx * sumy / numactive;
203 
204  if (delx < 0.001 * dely || dely < 0.001 * delx ||
205  fabs(SSxy * SSxy / (SSxx * SSyy)) > 0.99) {
206  /* points are colinear */
207  return MUNSOLVABLE;
208  }
209 
210  /* CALCULATE THE FORWARD TRANSFORMATION COEFFICIENTS */
211 
212  G_message(_("Calculating forward transformation coefficients"));
213  status = calccoef(cp, E12tps, N12tps);
214 
215  if (status != MSUCCESS)
216  return status;
217 
218  /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS */
219 
220  tempptr = cp->e1;
221  cp->e1 = cp->e2;
222  cp->e2 = tempptr;
223  tempptr = cp->n1;
224  cp->n1 = cp->n2;
225  cp->n2 = tempptr;
226 
227  /* CALCULATE THE BACKWARD TRANSFORMATION COEFFICIENTS */
228 
229  G_message(_("Calculating backward transformation coefficients"));
230  status = calccoef(cp, E21tps, N21tps);
231 
232  /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS BACK */
233 
234  tempptr = cp->e1;
235  cp->e1 = cp->e2;
236  cp->e2 = tempptr;
237  tempptr = cp->n1;
238  cp->n1 = cp->n2;
239  cp->n2 = tempptr;
240 
241  return status;
242 }
243 
244 /***********************************************************************
245 
246  COMPUTE THE GEOREFFERENCING COEFFICIENTS
247  BASED ON A SET OF CONTROL POINTS
248 
249 ************************************************************************/
250 
251 static int calccoef(struct Control_Points *cp, double **E, double **N)
252 {
253  struct MATRIX m;
254  double *a;
255  double *b;
256  int numactive; /* NUMBER OF ACTIVE CONTROL POINTS */
257  int status, i;
258 
259  /* CALCULATE THE NUMBER OF VALID CONTROL POINTS */
260 
261  for (i = numactive = 0; i < cp->count; i++) {
262  if (cp->status[i] > 0)
263  numactive++;
264  }
265 
266  /* INITIALIZE MATRIX */
267 
268  m.n = numactive + 3;
269 
270  m.v = G_calloc(m.n * m.n, sizeof(double));
271  if (m.v == NULL)
272  G_fatal_error(_("%s: out of memory"), "I_compute_georef_equations_tps()");
273  a = G_calloc(m.n, sizeof(double));
274  if (a == NULL)
275  G_fatal_error(_("%s: out of memory"), "I_compute_georef_equations_tps()");
276  b = G_calloc(m.n, sizeof(double));
277  if (b == NULL)
278  G_fatal_error(_("%s: out of memory"), "I_compute_georef_equations_tps()");
279 
280  /* equation coefficients */
281  *E = G_calloc(m.n, sizeof(double));
282  if (*E == NULL)
283  G_fatal_error(_("%s: out of memory"), "I_compute_georef_equations_tps()");
284  *N = G_calloc(m.n, sizeof(double));
285  if (*N == NULL)
286  G_fatal_error(_("%s: out of memory"), "I_compute_georef_equations_tps()");
287 
288  status = calcls(cp, &m, a, b, *E, *N);
289 
290  G_free(m.v);
291  G_free(a);
292  G_free(b);
293 
294  return status;
295 }
296 
297 
298 /***********************************************************************
299 
300  CALCULATE THE TRANSFORMATION COEFFICIENTS FOR THIN PLATE SPLINE
301  INTERPOLATION.
302  THIS ROUTINE USES THE LEAST SQUARES METHOD TO COMPUTE THE COEFFICIENTS.
303 
304 ************************************************************************/
305 
306 static int calcls(struct Control_Points *cp, struct MATRIX *m,
307  double a[], double b[],
308  double E[], /* EASTING COEFFICIENTS */
309  double N[] /* NORTHING COEFFICIENTS */
310  )
311 {
312  int i, j, n, o, numactive = 0;
313  double dist = 0.0, dx, dy, regularization;
314 
315  /* INITIALIZE THE MATRIX AND THE TWO COLUMN VECTORS */
316 
317  for (i = 1; i <= m->n; i++) {
318  for (j = i; j <= m->n; j++) {
319  M(i, j) = 0.0;
320  if (i != j)
321  M(j, i) = 0.0;
322  }
323  a[i - 1] = b[i - 1] = 0.0;
324  }
325 
326  /* SUM THE UPPER HALF OF THE MATRIX AND THE COLUMN VECTORS ACCORDING TO
327  THE LEAST SQUARES METHOD OF SOLVING OVER DETERMINED SYSTEMS */
328 
329  for (n = 0; n < cp->count; n++) {
330  if (cp->status[n] > 0) {
331 
332  a[numactive + 3] = cp->e2[n];
333  b[numactive + 3] = cp->n2[n];
334 
335  numactive++;
336  M(1, numactive + 3) = 1.0;
337  M(2, numactive + 3) = cp->e1[n];
338  M(3, numactive + 3) = cp->n1[n];
339 
340  M(numactive + 3, 1) = 1.0;
341  M(numactive + 3, 2) = cp->e1[n];
342  M(numactive + 3, 3) = cp->n1[n];
343  }
344  }
345 
346  if (numactive < m->n - 3)
347  return MINTERR;
348 
349  i = 0;
350  for (n = 0; n < cp->count; n++) {
351  if (cp->status[n] > 0) {
352  i++;
353 
354  j = 0;
355  for (o = 0; o <= n; o++) {
356  if (cp->status[o] > 0) {
357  j++;
358  M(i + 3, j + 3) = tps_base_func(cp->e1[n], cp->n1[n],
359  cp->e1[o], cp->n1[o]);
360 
361  if (i != j)
362  M(j + 3, i + 3) = M(i + 3, j + 3);
363 
364  dx = cp->e1[n] - cp->e1[o];
365  dy = cp->n1[n] - cp->n1[o];
366  dist += sqrt(dx * dx + dy * dy);
367  }
368  }
369  }
370  }
371 
372  /* regularization */
373  dist /= (numactive * numactive);
374  regularization = 0.01 * dist * dist;
375 
376  /* set diagonal to regularization, but not the first 3x3 (global affine) */
377 
378  return solvemat(m, a, b, E, N);
379 }
380 
381 
382 /***********************************************************************
383 
384  SOLVE FOR THE 'E' AND 'N' COEFFICIENTS BY USING A SOMEWHAT MODIFIED
385  GAUSSIAN ELIMINATION METHOD.
386 
387  | M11 M12 ... M1n | | E0 | | a0 |
388  | M21 M22 ... M2n | | E1 | = | a1 |
389  | . . . . | | . | | . |
390  | Mn1 Mn2 ... Mnn | | En-1 | | an-1 |
391 
392  and
393 
394  | M11 M12 ... M1n | | N0 | | b0 |
395  | M21 M22 ... M2n | | N1 | = | b1 |
396  | . . . . | | . | | . |
397  | Mn1 Mn2 ... Mnn | | Nn-1 | | bn-1 |
398 
399 ************************************************************************/
400 
401 static int solvemat(struct MATRIX *m, double a[], double b[], double E[],
402  double N[])
403 {
404  int i, j, i2, j2, imark;
405  double factor, temp;
406  double pivot; /* ACTUAL VALUE OF THE LARGEST PIVOT CANDIDATE */
407 
408  for (i = 1; i <= m->n; i++) {
409  G_percent(i - 1, m->n, 4);
410  j = i;
411 
412  /* find row with largest magnitude value for pivot value */
413 
414  pivot = M(i, j);
415  imark = i;
416  for (i2 = i + 1; i2 <= m->n; i2++) {
417  temp = fabs(M(i2, j));
418  if (temp > fabs(pivot)) {
419  pivot = M(i2, j);
420  imark = i2;
421  }
422  }
423 
424  /* if the pivot is very small then the points are nearly co-linear */
425  /* co-linear points result in an undefined matrix, and nearly */
426  /* co-linear points results in a solution with rounding error */
427 
428  if (pivot == 0.0)
429  return MUNSOLVABLE;
430 
431  /* if row with highest pivot is not the current row, switch them */
432 
433  if (imark != i) {
434  for (j2 = 1; j2 <= m->n; j2++) {
435  temp = M(imark, j2);
436  M(imark, j2) = M(i, j2);
437  M(i, j2) = temp;
438  }
439 
440  temp = a[imark - 1];
441  a[imark - 1] = a[i - 1];
442  a[i - 1] = temp;
443 
444  temp = b[imark - 1];
445  b[imark - 1] = b[i - 1];
446  b[i - 1] = temp;
447  }
448 
449  /* compute zeros above and below the pivot, and compute
450  values for the rest of the row as well */
451 
452  for (i2 = 1; i2 <= m->n; i2++) {
453  if (i2 != i) {
454  factor = M(i2, j) / pivot;
455  for (j2 = j; j2 <= m->n; j2++)
456  M(i2, j2) -= factor * M(i, j2);
457  a[i2 - 1] -= factor * a[i - 1];
458  b[i2 - 1] -= factor * b[i - 1];
459  }
460  }
461  }
462  G_percent(1, 1, 1);
463 
464  /* SINCE ALL OTHER VALUES IN THE MATRIX ARE ZERO NOW, CALCULATE THE
465  COEFFICIENTS BY DIVIDING THE COLUMN VECTORS BY THE DIAGONAL VALUES. */
466 
467  for (i = 1; i <= m->n; i++) {
468  E[i - 1] = a[i - 1] / M(i, i);
469  N[i - 1] = b[i - 1] / M(i, i);
470  }
471 
472  return MSUCCESS;
473 }
474 
475 static double tps_base_func(const double x1, const double y1,
476  const double x2, const double y2)
477 {
478  /* official: r * r * log(r) */
479  double dist;
480 
481  if ((x1 == x2) && (y1 == y2))
482  return 0.0;
483 
484  dist = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1);
485 
486  return dist * log(dist) * 0.5;
487 }
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
int I_compute_georef_equations_tps(struct Control_Points *cp, double **E12tps, double **N12tps, double **E21tps, double **N21tps)
Definition: georef_tps.c:113
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:149
#define MSUCCESS
Definition: georef_tps.c:36
#define N
Definition: e_intersect.c:917
#define NULL
Definition: ccmath.h:32
double * n2
Definition: imagery.h:47
#define G_calloc(m, n)
Definition: defs/gis.h:113
void G_message(const char *,...) __attribute__((format(printf
double * n1
Definition: imagery.h:44
int I_georef_tps(double e1, double n1, double *e, double *n, double *E, double *N, struct Control_Points *cp, int fwd)
Definition: georef_tps.c:65
double b
Definition: r_raster.c:39
#define MIN(a, b)
Definition: gis.h:140
#define MUNSOLVABLE
Definition: georef_tps.c:38
#define MAX(a, b)
Definition: gis.h:135
int * status
Definition: imagery.h:49
#define M(row, col)
Definition: georef_tps.c:34
void G_percent(long, long, int)
Print percent complete messages.
Definition: percent.c:62
#define MINTERR
Definition: georef_tps.c:41
double * e2
Definition: imagery.h:46
#define _(str)
Definition: glocale.h:10
double * e1
Definition: imagery.h:43
#define MNPTERR
Definition: georef_tps.c:37