1
2 /****************************************************************************
3 *
4 * MODULE: i.ortho.transform (cloned from m.transform nee g.transform)
5 * AUTHOR(S): Brian J. Buckley
6 * Glynn Clements
7 * Hamish Bowman
8 * Markus Metz
9 * PURPOSE: Utility to compute transformation based upon GCPs and
10 * output error measurements
11 * COPYRIGHT: (C) 2006-2013 by the GRASS Development Team
12 *
13 * This program is free software under the GNU General Public
14 * License (>=v2). Read the file COPYING that comes with GRASS
15 * for details.
16 *
17 *****************************************************************************/
18
19 #include <unistd.h>
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include <string.h>
23 #include <math.h>
24
25 #include <grass/gis.h>
26 #include <grass/glocale.h>
27 #include <grass/imagery.h>
28 #include <grass/glocale.h>
29 #include "orthophoto.h"
30
31 struct Max
32 {
33 int idx;
34 double val;
35 };
36
37 struct Stats
38 {
39 struct Max x, y, g;
40 double sum2, rms;
41 };
42
43 static int summary;
44 static int forward;
45 static char **columns;
46 static int need_fwd;
47 static int need_rev;
48 static int need_fd;
49 static int need_rd;
50 static char *coord_file;
51
52 struct Ortho_Image_Group group;
53
54 static struct Ortho_Control_Points *points;
55
56 static int count;
57 static struct Stats fwd, rev;
58
59 /* target fns taken from i.ortho.rectify */
60 static int which_env = -1; /* 0 = cur, 1 = target */
61
select_current_env(void)62 int select_current_env(void)
63 {
64 if (which_env < 0) {
65 G_create_alt_env();
66 which_env = 0;
67 }
68 if (which_env != 0) {
69 G_switch_env();
70 which_env = 0;
71 }
72
73 return 0;
74 }
75
select_target_env(void)76 int select_target_env(void)
77 {
78 if (which_env < 0) {
79 G_create_alt_env();
80 which_env = 1;
81 }
82 if (which_env != 1) {
83 G_switch_env();
84 which_env = 1;
85 }
86
87 return 0;
88 }
89
get_target(void)90 static int get_target(void)
91 {
92 char location[GMAPSET_MAX];
93 char mapset[GMAPSET_MAX];
94 char buf[1024];
95 int stat;
96
97 if (!I_get_target(group.name, location, mapset)) {
98 sprintf(buf, _("Target information for group <%s> missing"), group.name);
99 goto error;
100 }
101
102 sprintf(buf, "%s/%s", G_gisdbase(), location);
103 if (access(buf, 0) != 0) {
104 sprintf(buf, _("Target location <%s> not found"), location);
105 goto error;
106 }
107 select_target_env();
108 G_setenv_nogisrc("LOCATION_NAME", location);
109 stat = G_mapset_permissions(mapset);
110 if (stat > 0) {
111 G_setenv_nogisrc("MAPSET", mapset);
112 select_current_env();
113 return 1;
114 }
115 sprintf(buf, _("Mapset <%s> in target location <%s> - "), mapset, location);
116 strcat(buf, stat == 0 ? _("permission denied") : _("not found"));
117 error:
118 strcat(buf, "\n");
119 strcat(buf, _("Please run i.target for group "));
120 strcat(buf, group.name);
121 G_fatal_error("%s", buf);
122 return 1; /* never reached */
123 }
124
update_max(struct Max * m,int n,double k)125 static void update_max(struct Max *m, int n, double k)
126 {
127 if (k > m->val) {
128 m->idx = n;
129 m->val = k;
130 }
131 }
132
update_stats(struct Stats * st,int n,double dx,double dy,double dg,double d2)133 static void update_stats(struct Stats *st, int n, double dx, double dy,
134 double dg, double d2)
135 {
136 update_max(&st->x, n, dx);
137 update_max(&st->y, n, dy);
138 update_max(&st->g, n, dg);
139 st->sum2 += d2;
140 }
141
diagonal(double * dg,double * d2,double dx,double dy)142 static void diagonal(double *dg, double *d2, double dx, double dy)
143 {
144 *d2 = dx * dx + dy * dy;
145 *dg = sqrt(*d2);
146 }
147
compute_transformation(void)148 static void compute_transformation(void)
149 {
150 int n, i, status;
151 double e0, e1, e2, n0, n1, n2, z1, z2;
152 struct Ortho_Control_Points temp_points;
153
154 /* compute photo <-> image equations */
155 group.ref_equation_stat = I_compute_ref_equations(&group.photo_points,
156 group.E12, group.N12,
157 group.E21, group.N21);
158
159 if (group.ref_equation_stat <= 0)
160 G_fatal_error(_("Error conducting transform (%d)"),
161 group.ref_equation_stat);
162
163 /* compute target <-> photo equations */
164
165 /* alloc and fill temp control points */
166 temp_points.count = 0;
167 temp_points.status = NULL;
168 temp_points.e1 = NULL;
169 temp_points.n1 = NULL;
170 temp_points.z1 = NULL;
171 temp_points.e2 = NULL;
172 temp_points.n2 = NULL;
173 temp_points.z2 = NULL;
174
175 /* e0, n0, equal photo coordinates not image coords */
176 for (i = 0; i < group.control_points.count; i++) {
177 status = group.control_points.status[i];
178 e1 = group.control_points.e1[i];
179 n1 = group.control_points.n1[i];
180 z1 = group.control_points.z1[i];
181 e2 = group.control_points.e2[i];
182 n2 = group.control_points.n2[i];
183 z2 = group.control_points.z2[i];
184
185 /* image to photo transformation */
186 I_georef(e1, n1, &e0, &n0, group.E12, group.N12, 1);
187 I_new_con_point(&temp_points, e0, n0, z1, e2, n2, z2, status);
188 }
189
190
191 group.con_equation_stat = I_compute_ortho_equations(&temp_points,
192 &group.camera_ref,
193 &group.camera_exp,
194 &group.XC, &group.YC,
195 &group.ZC,
196 &group.omega,
197 &group.phi,
198 &group.kappa,
199 &group.M,
200 &group.MI);
201
202 if (group.con_equation_stat <= 0)
203 G_fatal_error(_("Error conducting transform (%d)"),
204 group.con_equation_stat);
205
206 count = 0;
207
208 for (n = 0; n < points->count; n++) {
209 double etmp, ntmp;
210 double fx, fy, fd, fd2;
211 double rx, ry, rd, rd2;
212
213 if (points->status[n] <= 0)
214 continue;
215
216 count++;
217
218 fd = fd2 = rd = rd2 = 0;
219
220 if (need_fwd) {
221 /* image -> photo -> target */
222
223 /* image coordinates ex, nx to photo coordinates ex1, nx1 */
224 I_georef(points->e1[n], points->n1[n], &etmp, &ntmp, group.E12, group.N12, 1);
225
226 /* photo coordinates ex1, nx1 to target coordinates e1, n1 */
227 I_inverse_ortho_ref(etmp, ntmp, points->z1[n], &e2, &n2, &z2,
228 &group.camera_ref,
229 group.XC, group.YC, group.ZC, group.MI);
230
231 fx = fabs(e2 - points->e2[n]);
232 fy = fabs(n2 - points->n2[n]);
233
234 if (need_fd)
235 diagonal(&fd, &fd2, fx, fy);
236
237 if (summary)
238 update_stats(&fwd, n, fx, fy, fd, fd2);
239 }
240
241 if (need_rev) {
242 /* target -> photo -> image */
243
244 /* target coordinates e1, n1 to photo coordinates ex1, nx1 */
245 I_ortho_ref(points->e2[n], points->n2[n], points->z2[n],
246 &etmp, &ntmp, &z2, &group.camera_ref,
247 group.XC, group.YC, group.ZC, group.M);
248
249 /* photo coordinates ex1, nx1 to image coordinates ex, nx */
250 I_georef(etmp, ntmp, &e1, &n1, group.E21, group.N21, 1);
251
252 rx = fabs(e1 - points->e1[n]);
253 ry = fabs(n1 - points->n1[n]);
254
255 if (need_rd)
256 diagonal(&rd, &rd2, rx, ry);
257
258 if (summary)
259 update_stats(&rev, n, rx, ry, rd, rd2);
260 }
261
262 if (!columns[0])
263 continue;
264
265 if (coord_file)
266 continue;
267
268 for (i = 0;; i++) {
269 const char *col = columns[i];
270
271 if (!col)
272 break;
273
274 if (strcmp("idx", col) == 0)
275 printf(" %d", n);
276 if (strcmp("src", col) == 0)
277 printf(" %f %f", points->e1[n], points->n1[n]);
278 if (strcmp("dst", col) == 0)
279 printf(" %f %f", points->e2[n], points->n2[n]);
280 if (strcmp("fwd", col) == 0)
281 printf(" %f %f", e2, n2);
282 if (strcmp("rev", col) == 0)
283 printf(" %f %f", e1, n1);
284 if (strcmp("fxy", col) == 0)
285 printf(" %f %f", fx, fy);
286 if (strcmp("rxy", col) == 0)
287 printf(" %f %f", rx, ry);
288 if (strcmp("fd", col) == 0)
289 printf(" %f", fd);
290 if (strcmp("rd", col) == 0)
291 printf(" %f", rd);
292 }
293
294 printf("\n");
295 }
296
297 if (summary && count > 0) {
298 fwd.rms = sqrt(fwd.sum2 / count);
299 rev.rms = sqrt(rev.sum2 / count);
300 }
301 }
302
do_max(char name,const struct Max * m)303 static void do_max(char name, const struct Max *m)
304 {
305 printf("%c[%d] = %.2f\n", name, m->idx, m->val);
306 }
307
do_stats(const char * name,const struct Stats * st)308 static void do_stats(const char *name, const struct Stats *st)
309 {
310 printf("%s:\n", name);
311 do_max('x', &st->x);
312 do_max('y', &st->y);
313 do_max('g', &st->g);
314 printf("RMS = %.2f\n", st->rms);
315 }
316
analyze(void)317 static void analyze(void)
318 {
319 if (group.ref_equation_stat == -1)
320 G_warning(_("Poorly placed image to photo control points"));
321 else if (group.con_equation_stat == -1)
322 G_warning(_("Poorly placed image to target control points"));
323 else if (group.ref_equation_stat == -2 || group.con_equation_stat == -2)
324 G_fatal_error(_("Insufficient memory"));
325 else if (group.ref_equation_stat < 0 || group.con_equation_stat < 0)
326 G_fatal_error(_("Parameter error"));
327 else if (group.ref_equation_stat == 0 || group.con_equation_stat == 0)
328 G_fatal_error(_("No active control points"));
329 else if (summary) {
330 printf("Number of active points: %d\n", count);
331 do_stats("Forward", &fwd);
332 do_stats("Reverse", &rev);
333 }
334 }
335
parse_format(void)336 static void parse_format(void)
337 {
338 int i;
339
340 if (summary) {
341 need_fwd = need_rev = need_fd = need_rd = 1;
342 return;
343 }
344
345 if (!columns)
346 return;
347
348 for (i = 0;; i++) {
349 const char *col = columns[i];
350
351 if (!col)
352 break;
353
354 if (strcmp("fwd", col) == 0)
355 need_fwd = 1;
356 if (strcmp("fxy", col) == 0)
357 need_fwd = 1;
358 if (strcmp("fd", col) == 0)
359 need_fwd = need_fd = 1;
360 if (strcmp("rev", col) == 0)
361 need_rev = 1;
362 if (strcmp("rxy", col) == 0)
363 need_rev = 1;
364 if (strcmp("rd", col) == 0)
365 need_rev = need_rd = 1;
366 }
367 }
368
dump_cooefs(void)369 static void dump_cooefs(void)
370 {
371 int i;
372
373 for (i = 0; i < 3; i++)
374 fprintf(stdout, "E%d=%.15g\n", i, forward ? group.E12[i] : group.E21[i]);
375
376 for (i = 0; i < 3; i++)
377 fprintf(stdout, "N%d=%.15g\n", i, forward ? group.N12[i] : group.N21[i]);
378
379 /* print ortho transformation matrix ? */
380 }
381
xform_value(double east,double north,double height)382 static void xform_value(double east, double north, double height)
383 {
384 double e1, n1, z1, xe, xn, xz;
385
386 if (forward) {
387 /* image -> photo -> target */
388 I_georef(east, north, &e1, &n1, group.E12, group.N12, 1);
389 z1 = height;
390 I_inverse_ortho_ref(e1, n1, z1, &xe, &xn, &xz, &group.camera_ref,
391 group.XC, group.YC, group.ZC, group.MI);
392 xz = z1;
393 }
394 else {
395 /* target -> photo -> image */
396 I_ortho_ref(east, north, height, &e1, &n1, &z1, &group.camera_ref,
397 group.XC, group.YC, group.ZC, group.M);
398
399 I_georef(e1, n1, &xe, &xn, group.E21, group.N21, 1);
400 xz = 0.;
401 }
402
403 fprintf(stdout, "%.15g %.15g %.15g\n", xe, xn, xz);
404 }
405
do_pt_xforms(void)406 static void do_pt_xforms(void)
407 {
408 double easting, northing, height;
409 int ret;
410 FILE *fp;
411
412 if (strcmp(coord_file, "-") == 0)
413 fp = stdin;
414 else {
415 fp = fopen(coord_file, "r");
416 if (!fp)
417 G_fatal_error(_("Unable to open file <%s>"), coord_file);
418 }
419
420 for (;;) {
421 char buf[1024];
422
423 if (!G_getl2(buf, sizeof(buf), fp))
424 break;
425
426 if ((buf[0] == '#') || (buf[0] == '\0'))
427 continue;
428
429 /* ? sscanf(buf, "%s %s", &east_str, &north_str)
430 ? G_scan_easting(,,-1)
431 ? G_scan_northing(,,-1) */
432 /* ? muliple delims with sscanf(buf, "%[ ,|\t]", &dummy) ? */
433
434 ret = sscanf(buf, "%lf %lf %lf", &easting, &northing, &height);
435 if (ret != 3)
436 G_fatal_error(_("Invalid coordinates: [%s]"), buf);
437
438 xform_value(easting, northing, height);
439 }
440
441 if (fp != stdin)
442 fclose(fp);
443 }
444
main(int argc,char ** argv)445 int main(int argc, char **argv)
446 {
447 struct Option *grp, *fmt, *xfm_pts;
448 struct Flag *sum, *rev_flag, *dump_flag, *pan_flag;
449 struct GModule *module;
450 char *desc;
451 char *camera;
452
453 G_gisinit(argv[0]);
454
455 /* Get Args */
456 module = G_define_module();
457 G_add_keyword(_("imagery"));
458 G_add_keyword(_("orthorectify"));
459 G_add_keyword(_("transformation"));
460 G_add_keyword(_("GCP"));
461 module->description =
462 _("Computes a coordinate transformation based on the control points.");
463
464 grp = G_define_standard_option(G_OPT_I_GROUP);
465
466 fmt = G_define_option();
467 fmt->key = "format";
468 fmt->type = TYPE_STRING;
469 fmt->required = NO;
470 fmt->multiple = YES;
471 fmt->options = "idx,src,dst,fwd,rev,fxy,rxy,fd,rd";
472 desc = NULL;
473 G_asprintf(&desc,
474 "idx;%s;src;%s;dst;%s;fwd;%s;rev;%s;fxy;%s;rxy;%s;fd;%s;rd;%s",
475 _("point index"),
476 _("source coordinates"),
477 _("destination coordinates"),
478 _("forward coordinates (destination)"),
479 _("reverse coordinates (source)"),
480 _("forward coordinates difference (destination)"),
481 _("reverse coordinates difference (source)"),
482 _("forward error (destination)"),
483 _("reverse error (source)"));
484 fmt->descriptions = desc;
485 fmt->answer = "fd,rd";
486 fmt->description = _("Output format");
487
488 sum = G_define_flag();
489 sum->key = 's';
490 sum->description = _("Display summary information");
491
492 xfm_pts = G_define_standard_option(G_OPT_F_INPUT);
493 xfm_pts->key = "coords";
494 xfm_pts->required = NO;
495 xfm_pts->label =
496 _("File containing coordinates to transform (\"-\" to read from stdin)");
497 xfm_pts->description = _("Local x,y,z coordinates to target east,north,height");
498
499 rev_flag = G_define_flag();
500 rev_flag->key = 'r';
501 rev_flag->label = _("Reverse transform of coords file or coeff. dump");
502 rev_flag->description = _("Target east,north,height coordinates to local x,y,z");
503
504 dump_flag = G_define_flag();
505 dump_flag->key = 'x';
506 dump_flag->description = _("Display transform matrix coefficients");
507
508 pan_flag = G_define_flag();
509 pan_flag->key = 'p';
510 pan_flag->description = _("Enable panorama camera correction");
511
512 if (G_parser(argc, argv))
513 exit(EXIT_FAILURE);
514
515
516 G_strip(grp->answer);
517 strcpy(group.name, grp->answer);
518
519 summary = !!sum->answer;
520 columns = fmt->answers;
521 forward = !rev_flag->answer;
522 coord_file = xfm_pts->answer;
523
524 if (pan_flag->answer)
525 I_ortho_panorama();
526
527 if (!I_get_ref_points(group.name, &group.photo_points)) {
528 G_fatal_error(_("Can not read reference points for group <%s>"),
529 group.name);
530 }
531 if (!I_get_con_points(group.name, &group.control_points)) {
532 group.control_points.count = 0;
533 }
534
535 camera = (char *)G_malloc(GNAME_MAX * sizeof(char));
536 if (!I_get_group_camera(group.name, camera))
537 G_fatal_error(_("No camera reference file selected for group <%s>"),
538 group.name);
539
540 if (!I_get_cam_info(camera, &group.camera_ref))
541 G_fatal_error(_("Bad format in camera file for group <%s>"),
542 group.name);
543
544 /* get initial camera exposure station, if any */
545 if (I_find_initial(group.name)) {
546 if (!I_get_init_info(group.name, &group.camera_exp))
547 G_warning(_("Bad format in initial exposure station file for group <%s>"),
548 group.name);
549 }
550
551 /* get the target */
552 get_target();
553
554 points = &group.control_points;
555
556 parse_format();
557
558 /* I_compute_ortho_equations() must be run in the target location */
559 select_target_env();
560 compute_transformation();
561 select_current_env();
562
563 analyze();
564
565 if (dump_flag->answer)
566 dump_cooefs();
567
568 if (coord_file)
569 do_pt_xforms();
570
571 return 0;
572 }
573