1
2 /****************************************************************************
3 *
4 * MODULE: m.transform (nee g.transform)
5 * AUTHOR(S): Brian J. Buckley
6 * Glynn Clements
7 * Hamish Bowman
8 * PURPOSE: Utility to compute transformation based upon GCPs and
9 * output error measurements
10 * COPYRIGHT: (C) 2006-2010 by the GRASS Development Team
11 *
12 * This program is free software under the GNU General Public
13 * License (>=v2). Read the file COPYING that comes with GRASS
14 * for details.
15 *
16 *****************************************************************************/
17
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <string.h>
21 #include <math.h>
22
23 #include <grass/gis.h>
24 #include <grass/glocale.h>
25 #include <grass/imagery.h>
26 #include <grass/glocale.h>
27
28 struct Max
29 {
30 int idx;
31 double val;
32 };
33
34 struct Stats
35 {
36 struct Max x, y, g;
37 double sum2, rms;
38 };
39
40 static char *name;
41 static int order;
42 static int summary;
43 static int forward;
44 static char **columns;
45 static int need_fwd;
46 static int need_rev;
47 static int need_fd;
48 static int need_rd;
49 static char *coord_file;
50
51 static double E12[10], N12[10], E21[10], N21[10];
52
53 static struct Control_Points points;
54
55 static int equation_stat;
56
57 static int count;
58 static struct Stats fwd, rev;
59
update_max(struct Max * m,int n,double k)60 static void update_max(struct Max *m, int n, double k)
61 {
62 if (k > m->val) {
63 m->idx = n;
64 m->val = k;
65 }
66 }
67
update_stats(struct Stats * st,int n,double dx,double dy,double dg,double d2)68 static void update_stats(struct Stats *st, int n, double dx, double dy,
69 double dg, double d2)
70 {
71 update_max(&st->x, n, dx);
72 update_max(&st->y, n, dy);
73 update_max(&st->g, n, dg);
74 st->sum2 += d2;
75 }
76
diagonal(double * dg,double * d2,double dx,double dy)77 static void diagonal(double *dg, double *d2, double dx, double dy)
78 {
79 *d2 = dx * dx + dy * dy;
80 *dg = sqrt(*d2);
81 }
82
compute_transformation(void)83 static void compute_transformation(void)
84 {
85 static const int order_pnts[3] = { 3, 6, 10 };
86 int n, i;
87
88 equation_stat =
89 I_compute_georef_equations(&points, E12, N12, E21, N21, order);
90
91 if (equation_stat == 0)
92 G_fatal_error(_("Not enough points, %d are required"),
93 order_pnts[order - 1]);
94
95 if (equation_stat <= 0)
96 G_fatal_error(_("Error conducting transform (%d)"), equation_stat);
97
98 count = 0;
99
100 for (n = 0; n < points.count; n++) {
101 double e1, n1, e2, n2;
102 double fx, fy, fd, fd2;
103 double rx, ry, rd, rd2;
104
105 if (points.status[n] <= 0)
106 continue;
107
108 count++;
109
110 if (need_fwd) {
111 I_georef(points.e1[n], points.n1[n], &e2, &n2, E12, N12, order);
112
113 fx = fabs(e2 - points.e2[n]);
114 fy = fabs(n2 - points.n2[n]);
115
116 if (need_fd)
117 diagonal(&fd, &fd2, fx, fy);
118
119 if (summary)
120 update_stats(&fwd, n, fx, fy, fd, fd2);
121 }
122
123 if (need_rev) {
124 I_georef(points.e2[n], points.n2[n], &e1, &n1, E21, N21, order);
125
126 rx = fabs(e1 - points.e1[n]);
127 ry = fabs(n1 - points.n1[n]);
128
129 if (need_rd)
130 diagonal(&rd, &rd2, rx, ry);
131
132 if (summary)
133 update_stats(&rev, n, rx, ry, rd, rd2);
134 }
135
136 if (!columns[0])
137 continue;
138
139 if (coord_file)
140 continue;
141
142 for (i = 0;; i++) {
143 const char *col = columns[i];
144
145 if (!col)
146 break;
147
148 if (strcmp("idx", col) == 0)
149 printf(" %d", n);
150 if (strcmp("src", col) == 0)
151 printf(" %f %f", points.e1[n], points.n1[n]);
152 if (strcmp("dst", col) == 0)
153 printf(" %f %f", points.e2[n], points.n2[n]);
154 if (strcmp("fwd", col) == 0)
155 printf(" %f %f", e2, n2);
156 if (strcmp("rev", col) == 0)
157 printf(" %f %f", e1, n1);
158 if (strcmp("fxy", col) == 0)
159 printf(" %f %f", fx, fy);
160 if (strcmp("rxy", col) == 0)
161 printf(" %f %f", rx, ry);
162 if (strcmp("fd", col) == 0)
163 printf(" %f", fd);
164 if (strcmp("rd", col) == 0)
165 printf(" %f", rd);
166 }
167
168 printf("\n");
169 }
170
171 if (summary && count > 0) {
172 fwd.rms = sqrt(fwd.sum2 / count);
173 rev.rms = sqrt(rev.sum2 / count);
174 }
175 }
176
do_max(char name,const struct Max * m)177 static void do_max(char name, const struct Max *m)
178 {
179 printf("%c[%d] = %.2f\n", name, m->idx, m->val);
180 }
181
do_stats(const char * name,const struct Stats * st)182 static void do_stats(const char *name, const struct Stats *st)
183 {
184 printf("%s:\n", name);
185 do_max('x', &st->x);
186 do_max('y', &st->y);
187 do_max('g', &st->g);
188 printf("RMS = %.2f\n", st->rms);
189 }
190
analyze(void)191 static void analyze(void)
192 {
193 if (equation_stat == -1)
194 G_warning(_("Poorly placed control points"));
195 else if (equation_stat == -2)
196 G_fatal_error(_("Insufficient memory"));
197 else if (equation_stat < 0)
198 G_fatal_error(_("Parameter error"));
199 else if (equation_stat == 0)
200 G_fatal_error(_("No active control points"));
201 else if (summary) {
202 printf("Number of active points: %d\n", count);
203 do_stats("Forward", &fwd);
204 do_stats("Reverse", &rev);
205 }
206 }
207
parse_format(void)208 static void parse_format(void)
209 {
210 int i;
211
212 if (summary) {
213 need_fwd = need_rev = need_fd = need_rd = 1;
214 return;
215 }
216
217 if (!columns)
218 return;
219
220 for (i = 0;; i++) {
221 const char *col = columns[i];
222
223 if (!col)
224 break;
225
226 if (strcmp("fwd", col) == 0)
227 need_fwd = 1;
228 if (strcmp("fxy", col) == 0)
229 need_fwd = 1;
230 if (strcmp("fd", col) == 0)
231 need_fwd = need_fd = 1;
232 if (strcmp("rev", col) == 0)
233 need_rev = 1;
234 if (strcmp("rxy", col) == 0)
235 need_rev = 1;
236 if (strcmp("rd", col) == 0)
237 need_rev = need_rd = 1;
238 }
239 }
240
dump_cooefs(void)241 static void dump_cooefs(void)
242 {
243 int i;
244 static const int order_pnts[3] = { 3, 6, 10 };
245
246 for (i = 0; i < order_pnts[order - 1]; i++)
247 fprintf(stdout, "E%d=%.15g\n", i, forward ? E12[i] : E21[i]);
248
249 for (i = 0; i < order_pnts[order - 1]; i++)
250 fprintf(stdout, "N%d=%.15g\n", i, forward ? N12[i] : N21[i]);
251 }
252
xform_value(double east,double north)253 static void xform_value(double east, double north)
254 {
255 double xe, xn;
256
257 if(forward)
258 I_georef(east, north, &xe, &xn, E12, N12, order);
259 else
260 I_georef(east, north, &xe, &xn, E21, N21, order);
261
262 fprintf(stdout, "%.15g %.15g\n", xe, xn);
263 }
264
do_pt_xforms(void)265 static void do_pt_xforms(void)
266 {
267 double easting, northing;
268 int ret;
269 FILE *fp;
270
271 if (strcmp(coord_file, "-") == 0)
272 fp = stdin;
273 else {
274 fp = fopen(coord_file, "r");
275 if (!fp)
276 G_fatal_error(_("Unable to open file <%s>"), coord_file);
277 }
278
279 for (;;) {
280 char buf[64];
281
282 if (!G_getl2(buf, sizeof(buf), fp))
283 break;
284
285 if ((buf[0] == '#') || (buf[0] == '\0'))
286 continue;
287
288 /* ? sscanf(buf, "%s %s", &east_str, &north_str)
289 ? G_scan_easting(,,-1)
290 ? G_scan_northing(,,-1) */
291 /* ? muliple delims with sscanf(buf, "%[ ,|\t]", &dummy) ? */
292
293 ret = sscanf(buf, "%lf %lf", &easting, &northing);
294 if (ret != 2)
295 G_fatal_error(_("Invalid coordinates: [%s]"), buf);
296
297 xform_value(easting, northing);
298 }
299
300 if (fp != stdin)
301 fclose(fp);
302 }
303
304
main(int argc,char ** argv)305 int main(int argc, char **argv)
306 {
307 struct Option *grp, *val, *fmt, *xfm_pts;
308 struct Flag *sum, *rev_flag, *dump_flag;
309 struct GModule *module;
310 char *desc;
311
312 G_gisinit(argv[0]);
313
314 /* Get Args */
315 module = G_define_module();
316 G_add_keyword(_("miscellaneous"));
317 G_add_keyword(_("transformation"));
318 G_add_keyword("GCP");
319 module->description =
320 _("Computes a coordinate transformation based on the control points.");
321
322 grp = G_define_standard_option(G_OPT_I_GROUP);
323
324 val = G_define_option();
325 val->key = "order";
326 val->type = TYPE_INTEGER;
327 val->required = YES;
328 val->options = "1-3";
329 val->answer = "1";
330 val->description = _("Rectification polynomial order");
331
332 fmt = G_define_option();
333 fmt->key = "format";
334 fmt->type = TYPE_STRING;
335 fmt->required = NO;
336 fmt->multiple = YES;
337 fmt->options = "idx,src,dst,fwd,rev,fxy,rxy,fd,rd";
338 desc = NULL;
339 G_asprintf(&desc,
340 "idx;%s;src;%s;dst;%s;fwd;%s;rev;%s;fxy;%s;rxy;%s;fd;%s;rd;%s",
341 _("point index"),
342 _("source coordinates"),
343 _("destination coordinates"),
344 _("forward coordinates (destination)"),
345 _("reverse coordinates (source)"),
346 _("forward coordinates difference (destination)"),
347 _("reverse coordinates difference (source)"),
348 _("forward error (destination)"),
349 _("reverse error (source)"));
350 fmt->descriptions = desc;
351 fmt->answer = "fd,rd";
352 fmt->description = _("Output format");
353
354 sum = G_define_flag();
355 sum->key = 's';
356 sum->description = _("Display summary information");
357
358 xfm_pts = G_define_standard_option(G_OPT_F_INPUT);
359 xfm_pts->required = NO;
360 xfm_pts->label =
361 _("File containing coordinates to transform (\"-\" to read from stdin)");
362 xfm_pts->description = _("Local x,y coordinates to target east,north");
363
364 rev_flag = G_define_flag();
365 rev_flag->key = 'r';
366 rev_flag->label = _("Reverse transform of coords file or coeff. dump");
367 rev_flag->description = _("Target east,north coordinates to local x,y");
368
369 dump_flag = G_define_flag();
370 dump_flag->key = 'x';
371 dump_flag->description = _("Display transform matrix coefficients");
372
373 if (G_parser(argc, argv))
374 exit(EXIT_FAILURE);
375
376
377 name = grp->answer;
378 order = atoi(val->answer);
379 summary = !!sum->answer;
380 columns = fmt->answers;
381 forward = !rev_flag->answer;
382 coord_file = xfm_pts->answer;
383
384 I_get_control_points(name, &points);
385
386 parse_format();
387
388 compute_transformation();
389
390 I_put_control_points(name, &points);
391
392 analyze();
393
394 if(dump_flag->answer)
395 dump_cooefs();
396
397 if(coord_file)
398 do_pt_xforms();
399
400 return 0;
401 }
402