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