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