34 #define M(row,col) m->v[(((row)-1)*(m->n))+(col)-1] 38 #define MUNSOLVABLE -1 51 static int calccoef(
struct Control_Points *,
double **,
double **);
52 static int calcls(
struct Control_Points *,
struct MATRIX *,
double *,
53 double *,
double *,
double *);
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 *);
76 double dist, *pe, *pn;
88 *e = E[0] + e1 * E[1] + n1 * E[2];
89 *n = N[0] + e1 * N[1] + n1 * N[2];
92 for (i = 0, j = 0; i < cp->
count; i++) {
95 dist = tps_base_func(e1, n1, pe[i], pn[i]);
97 *e += E[j + 3] * dist;
98 *n += N[j + 3] * dist;
114 double **E12tps,
double **N12tps,
115 double **E21tps,
double **N21tps)
120 double xmax, xmin, ymax, ymin;
123 double sumx, sumy, sumx2, sumy2, sumxy;
124 double SSxx, SSyy, SSxy;
128 for (i = numactive = 0; i < cp->
count; i++) {
136 if (numactive > 100000)
139 xmin = xmax = cp->
e1[0];
140 ymin = ymax = cp->
n1[0];
142 sumx = sumy = sumx2 = sumy2 = sumxy = 0.0;
144 for (i = 0; i < cp->
count; i++ ) {
149 xmax =
MAX(xmax, xx);
150 xmin =
MIN(xmin, xx);
151 ymax =
MAX(ymax, yy);
152 ymin =
MIN(ymin, yy);
165 SSxx = sumx2 - sumx * sumx / numactive;
166 SSyy = sumy2 - sumy * sumy / numactive;
167 SSxy = sumxy - sumx * sumy / numactive;
169 if (delx < 0.001 * dely || dely < 0.001 * delx ||
170 fabs(SSxy * SSxy / (SSxx * SSyy)) > 0.99) {
175 xmin = xmax = cp->
e2[0];
176 ymin = ymax = cp->
n2[0];
178 sumx = sumy = sumx2 = sumy2 = sumxy = 0.0;
179 for (i = 0; i < cp->
count; i++ ) {
184 xmax =
MAX(xmax, xx);
185 xmin =
MIN(xmin, xx);
186 ymax =
MAX(ymax, yy);
187 ymin =
MIN(ymin, yy);
200 SSxx = sumx2 - sumx * sumx / numactive;
201 SSyy = sumy2 - sumy * sumy / numactive;
202 SSxy = sumxy - sumx * sumy / numactive;
204 if (delx < 0.001 * dely || dely < 0.001 * delx ||
205 fabs(SSxy * SSxy / (SSxx * SSyy)) > 0.99) {
212 G_message(
_(
"Calculating forward transformation coefficients"));
213 status = calccoef(cp, E12tps, N12tps);
229 G_message(
_(
"Calculating backward transformation coefficients"));
230 status = calccoef(cp, E21tps, N21tps);
261 for (i = numactive = 0; i < cp->
count; i++) {
270 m.v =
G_calloc(m.n * m.n,
sizeof(
double));
272 G_fatal_error(
_(
"%s: out of memory"),
"I_compute_georef_equations_tps()");
275 G_fatal_error(
_(
"%s: out of memory"),
"I_compute_georef_equations_tps()");
278 G_fatal_error(
_(
"%s: out of memory"),
"I_compute_georef_equations_tps()");
283 G_fatal_error(
_(
"%s: out of memory"),
"I_compute_georef_equations_tps()");
286 G_fatal_error(
_(
"%s: out of memory"),
"I_compute_georef_equations_tps()");
288 status = calcls(cp, &m, a, b, *E, *N);
307 double a[],
double b[],
312 int i, j, n, o, numactive = 0;
313 double dist = 0.0, dx, dy, regularization;
317 for (i = 1; i <= m->n; i++) {
318 for (j = i; j <= m->n; j++) {
323 a[i - 1] =
b[i - 1] = 0.0;
329 for (n = 0; n < cp->
count; n++) {
332 a[numactive + 3] = cp->
e2[n];
333 b[numactive + 3] = cp->
n2[n];
336 M(1, numactive + 3) = 1.0;
337 M(2, numactive + 3) = cp->
e1[n];
338 M(3, numactive + 3) = cp->
n1[n];
340 M(numactive + 3, 1) = 1.0;
341 M(numactive + 3, 2) = cp->
e1[n];
342 M(numactive + 3, 3) = cp->
n1[n];
346 if (numactive < m->n - 3)
350 for (n = 0; n < cp->
count; n++) {
355 for (o = 0; o <= n; o++) {
358 M(i + 3, j + 3) = tps_base_func(cp->
e1[n], cp->
n1[n],
359 cp->
e1[o], cp->
n1[o]);
362 M(j + 3, i + 3) =
M(i + 3, j + 3);
364 dx = cp->
e1[n] - cp->
e1[o];
365 dy = cp->
n1[n] - cp->
n1[o];
366 dist += sqrt(dx * dx + dy * dy);
373 dist /= (numactive * numactive);
374 regularization = 0.01 * dist * dist;
378 return solvemat(m, a,
b, E, N);
401 static int solvemat(
struct MATRIX *m,
double a[],
double b[],
double E[],
404 int i, j, i2, j2, imark;
408 for (i = 1; i <= m->n; i++) {
416 for (i2 = i + 1; i2 <= m->n; i2++) {
417 temp = fabs(
M(i2, j));
418 if (temp > fabs(pivot)) {
434 for (j2 = 1; j2 <= m->n; j2++) {
436 M(imark, j2) =
M(i, j2);
441 a[imark - 1] = a[i - 1];
445 b[imark - 1] =
b[i - 1];
452 for (i2 = 1; i2 <= m->n; i2++) {
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];
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);
475 static double tps_base_func(
const double x1,
const double y1,
476 const double x2,
const double y2)
481 if ((x1 == x2) && (y1 == y2))
484 dist = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1);
486 return dist * log(dist) * 0.5;
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)
void G_free(void *)
Free allocated memory.
void G_message(const char *,...) __attribute__((format(printf
int I_georef_tps(double e1, double n1, double *e, double *n, double *E, double *N, struct Control_Points *cp, int fwd)
void G_percent(long, long, int)
Print percent complete messages.