1
2 /****************************************************************************
3 *
4 * MODULE: i.ortho.rectify (former i.photo.rectify)
5 * AUTHOR(S): Mike Baba, DBA Systems, Inc. (original contributor)
6 * Markus Neteler <neteler itc.it>,
7 * Bernhard Reiter <bernhard intevation.de>,
8 * Glynn Clements <glynn gclements.plus.com>,
9 * Hamish Bowman <hamish_b yahoo.com>,
10 * Markus Metz
11 *
12 * PURPOSE: Rectifies an image by using the image to photo coordinate
13 * and photo to target transformation matrices
14 * COPYRIGHT: (C) 1999-2010 by the GRASS Development Team
15 *
16 * This program is free software under the GNU General Public
17 * License (>=v2). Read the file COPYING that comes with GRASS
18 * for details.
19 *
20 *****************************************************************************/
21
22 #include <stdlib.h>
23 #include <string.h>
24 #include "global.h"
25
26 int seg_mb_img, seg_mb_elev;
27
28 char *elev_name;
29 char *elev_mapset;
30
31 func interpolate;
32
33 struct Cell_head target_window;
34
35 void err_exit(struct Ref *, char *, char *);
36
37 /* modify this table to add new methods */
38 struct menu menu[] = {
39 {p_nearest, "nearest", "nearest neighbor"},
40 {p_bilinear, "linear", "linear interpolation"},
41 {p_cubic, "cubic", "cubic convolution"},
42 {p_lanczos, "lanczos", "lanczos filter"},
43 {p_bilinear_f, "linear_f", "linear interpolation with fallback"},
44 {p_cubic_f, "cubic_f", "cubic convolution with fallback"},
45 {p_lanczos_f, "lanczos_f", "lanczos filter with fallback"},
46 {NULL, NULL, NULL}
47 };
48
49 static char *make_ipol_list(void);
50
main(int argc,char * argv[])51 int main(int argc, char *argv[])
52 {
53 char extension[INAME_LEN];
54 char *ipolname; /* name of interpolation method */
55 int method;
56 char *seg_mb;
57 int i, m, k = 0;
58 int got_file = 0, target_overwrite = 0;
59 char *overstr;
60
61 struct Ortho_Image_Group group;
62 int *ref_list;
63 int n;
64 char *camera;
65 char tl[100];
66 char math_exp[100];
67 char units[100];
68 char nd[100];
69
70 struct Cell_head cellhd, elevhd;
71
72 struct Option *grp, /* imagery group */
73 *ifile, /* input files */
74 *ext, /* extension */
75 *tres, /* target resolution */
76 *mem, /* amount of memory for cache */
77 *angle, /* camera angle relative to ground surface */
78 *interpol; /* interpolation method:
79 nearest neighbor, bilinear, cubic */
80
81 struct Flag *c, *a, *pan_flag;
82 struct GModule *module;
83
84
85 G_gisinit(argv[0]);
86
87 module = G_define_module();
88 G_add_keyword(_("imagery"));
89 G_add_keyword(_("orthorectify"));
90 module->description =
91 _("Orthorectifies an image by using the image to photo coordinate transformation matrix.");
92
93 grp = G_define_standard_option(G_OPT_I_GROUP);
94
95 ifile = G_define_standard_option(G_OPT_R_INPUTS);
96 ifile->required = NO;
97
98 ext = G_define_option();
99 ext->key = "extension";
100 ext->type = TYPE_STRING;
101 ext->required = YES;
102 ext->multiple = NO;
103 ext->description = _("Output raster map(s) suffix");
104
105 tres = G_define_option();
106 tres->key = "resolution";
107 tres->type = TYPE_DOUBLE;
108 tres->required = NO;
109 tres->description = _("Target resolution (ignored if -c flag used)");
110
111 mem = G_define_standard_option(G_OPT_MEMORYMB);
112
113 ipolname = make_ipol_list();
114
115 interpol = G_define_option();
116 interpol->key = "method";
117 interpol->type = TYPE_STRING;
118 interpol->required = NO;
119 interpol->answer = "nearest";
120 interpol->options = ipolname;
121 interpol->description = _("Interpolation method to use");
122
123 angle = G_define_standard_option(G_OPT_R_OUTPUT);
124 angle->key = "angle";
125 angle->required = NO;
126 angle->description = _("Raster map with camera angle relative to ground surface");
127
128 c = G_define_flag();
129 c->key = 'c';
130 c->description =
131 _("Use current region settings in target location (def.=calculate smallest area)");
132
133 a = G_define_flag();
134 a->key = 'a';
135 a->description = _("Rectify all raster maps in group");
136
137 pan_flag = G_define_flag();
138 pan_flag->key = 'p';
139 pan_flag->description = _("Enable panorama camera correction");
140
141 if (G_parser(argc, argv))
142 exit(EXIT_FAILURE);
143
144 /* get the method */
145 for (method = 0; (ipolname = menu[method].name); method++)
146 if (strcmp(ipolname, interpol->answer) == 0)
147 break;
148
149 if (!ipolname)
150 G_fatal_error(_("<%s=%s> unknown %s"),
151 interpol->key, interpol->answer, interpol->key);
152 interpolate = menu[method].method;
153
154 G_strip(grp->answer);
155 strcpy(group.name, grp->answer);
156 strcpy(extension, ext->answer);
157
158 seg_mb = NULL;
159 if (mem->answer) {
160 if (atoi(mem->answer) > 0)
161 seg_mb = mem->answer;
162 }
163
164 if (!ifile->answers)
165 a->answer = 1; /* force all */
166
167 /* Find out how many files on command line */
168 if (!a->answer) {
169 for (k = 0; ifile->answers[k]; k++);
170 }
171
172 camera = (char *)G_malloc(GNAME_MAX * sizeof(char));
173 elev_name = (char *)G_malloc(GNAME_MAX * sizeof(char));
174 elev_mapset = (char *)G_malloc(GMAPSET_MAX * sizeof(char));
175
176 /* find group */
177 if (!I_find_group(group.name)) {
178 G_fatal_error(_("Group <%s> not found"), group.name);
179 }
180
181 /* determine the number of files in this group */
182 if (!I_get_group_ref(group.name, &group.group_ref)) {
183 G_warning(_("Location: %s"), G_location());
184 G_warning(_("Mapset: %s"), G_mapset());
185 G_fatal_error(_("Could not read REF file for group <%s>"),
186 group.name);
187 }
188
189 if (group.group_ref.nfiles <= 0) {
190 G_important_message(_("Group <%s> contains no raster maps; run i.group"),
191 grp->answer);
192 exit(EXIT_SUCCESS);
193 }
194
195 ref_list = (int *)G_malloc(group.group_ref.nfiles * sizeof(int));
196
197 if (a->answer) {
198 for (n = 0; n < group.group_ref.nfiles; n++) {
199 ref_list[n] = 1;
200 }
201 }
202 else {
203 char xname[GNAME_MAX], xmapset[GMAPSET_MAX], *name, *mapset;
204
205 for (n = 0; n < group.group_ref.nfiles; n++)
206 ref_list[n] = 0;
207
208 for (m = 0; m < k; m++) {
209 got_file = 0;
210 if (G_name_is_fully_qualified(ifile->answers[m], xname, xmapset)) {
211 name = xname;
212 mapset = xmapset;
213 }
214 else {
215 name = ifile->answers[m];
216 mapset = NULL;
217 }
218
219 got_file = 0;
220 for (n = 0; n < group.group_ref.nfiles; n++) {
221 if (mapset) {
222 if (strcmp(name, group.group_ref.file[n].name) == 0 &&
223 strcmp(mapset, group.group_ref.file[n].mapset) == 0) {
224 got_file = 1;
225 ref_list[n] = 1;
226 break;
227 }
228 }
229 else {
230 if (strcmp(name, group.group_ref.file[n].name) == 0) {
231 got_file = 1;
232 ref_list[n] = 1;
233 break;
234 }
235 }
236 }
237 if (got_file == 0)
238 err_exit(&group.group_ref, ifile->answers[m], group.name);
239 }
240 }
241
242 /** look for camera info for this block **/
243 if (!I_get_group_camera(group.name, camera))
244 G_fatal_error(_("No camera reference file selected for group <%s>"),
245 group.name);
246
247 if (!I_get_cam_info(camera, &group.camera_ref))
248 G_fatal_error(_("Bad format in camera file for group <%s>"),
249 group.name);
250
251 /* get initial camera exposure station, if any */
252 if (I_find_initial(group.name)) {
253 if (!I_get_init_info(group.name, &group.camera_exp))
254 G_warning(_("Bad format in initial exposure station file for group <%s>"),
255 group.name);
256 }
257
258 /* panorama camera correction */
259 if (pan_flag->answer)
260 I_ortho_panorama();
261
262 /* get the target */
263 get_target(group.name);
264
265 /* read the reference points for the group, compute image-to-photo trans. */
266 get_ref_points(&group);
267
268 /* read the control points for the group, convert to photo coords. */
269 get_conz_points(&group);
270
271 /* Check the GRASS_OVERWRITE environment variable */
272 if ((overstr = getenv("GRASS_OVERWRITE"))) /* OK ? */
273 target_overwrite = atoi(overstr);
274
275 if (!target_overwrite) {
276 /* check if output exists in target location/mapset */
277 char result[GNAME_MAX];
278
279 select_target_env();
280 for (i = 0; i < group.group_ref.nfiles; i++) {
281 if (!ref_list[i])
282 continue;
283
284 strcpy(result, group.group_ref.file[i].name);
285 strcat(result, extension);
286
287 if (G_legal_filename(result) < 0)
288 G_fatal_error(_("Extension <%s> is illegal"), extension);
289
290 if (G_find_raster2(result, G_mapset())) {
291 G_warning(_("The following raster map already exists in"));
292 G_warning(_("target LOCATION %s, MAPSET %s:"),
293 G_location(), G_mapset());
294 G_warning("<%s>", result);
295 G_fatal_error(_("Orthorectification cancelled."));
296 }
297 }
298 if (angle->answer) {
299 if (G_find_raster2(angle->answer, G_mapset())) {
300 G_warning(_("The following raster map already exists in"));
301 G_warning(_("target LOCATION %s, MAPSET %s:"),
302 G_location(), G_mapset());
303 G_warning("<%s>", angle->answer);
304 G_fatal_error(_("Orthorectification cancelled."));
305 }
306 }
307
308 select_current_env();
309 }
310 else
311 G_debug(1, "Overwriting OK");
312
313 /* do not use current region in target location */
314 if (!c->answer) {
315 double res = -1;
316
317 if (tres->answer) {
318 if (!((res = atof(tres->answer)) > 0))
319 G_warning(_("Target resolution must be > 0, ignored"));
320 }
321 /* get reference window from imagery group */
322 get_ref_window(&group.group_ref, ref_list, &cellhd);
323 georef_window(&group, &cellhd, &target_window, res);
324 }
325
326 G_verbose_message(_("Using region: N=%f S=%f, E=%f W=%f"), target_window.north,
327 target_window.south, target_window.east, target_window.west);
328
329 G_debug(1, "Looking for elevation file in group: <%s>", group.name);
330
331 /* get the block elevation layer raster map in target location */
332 if (!I_get_group_elev(group.name, elev_name, elev_mapset, tl,
333 math_exp, units, nd))
334 G_fatal_error(_("No target elevation model selected for group <%s>"),
335 group.name);
336
337 G_debug(1, "Block elevation: <%s> in <%s>", elev_name, elev_mapset);
338
339 /* get the elevation layer header in target location */
340 select_target_env();
341 Rast_get_cellhd(elev_name, elev_mapset, &elevhd);
342 select_current_env();
343
344 /* determine memory for elevation and imagery */
345 seg_mb_img = seg_mb_elev = -1;
346 if (seg_mb) {
347 int max_rows, max_cols;
348 int nx, ny;
349 double max_mb_img, max_mb_elev;
350 int seg_mb_total;
351
352 max_rows = max_cols = 0;
353 for (i = 0; i < group.group_ref.nfiles; i++) {
354 if (ref_list[i]) {
355 Rast_get_cellhd(group.group_ref.file[i].name,
356 group.group_ref.file[i].mapset, &cellhd);
357 if (max_rows < cellhd.rows)
358 max_rows = cellhd.rows;
359 if (max_cols < cellhd.cols)
360 max_cols = cellhd.cols;
361 }
362 }
363
364 ny = (max_rows + BDIM - 1) / BDIM;
365 nx = (max_cols + BDIM - 1) / BDIM;
366
367 max_mb_img = ((double)nx * ny * sizeof(block)) / (1<<20);
368
369 ny = (target_window.rows + BDIM - 1) / BDIM;
370 nx = (target_window.cols + BDIM - 1) / BDIM;
371
372 max_mb_elev = ((double)nx * ny * sizeof(block)) / (1<<20);
373
374 if ((seg_mb_total = atoi(seg_mb)) > 0) {
375 seg_mb_elev = max_mb_elev * (seg_mb_total / (max_mb_img + max_mb_elev)) + 0.5;
376 seg_mb_img = max_mb_img * (seg_mb_total / (max_mb_img + max_mb_elev)) + 0.5;
377 }
378 }
379
380 /* go do it */
381 exec_rectify(&group, ref_list, extension, interpol->answer, angle->answer);
382
383 G_done_msg(" ");
384
385 exit(EXIT_SUCCESS);
386 }
387
388
err_exit(struct Ref * ref,char * file,char * grp)389 void err_exit(struct Ref *ref, char *file, char *grp)
390 {
391 int n;
392
393 G_warning(_("Input raster map <%s> does not exist in group <%s>."),
394 file, grp);
395 G_message(_("Try:"));
396
397 for (n = 0; n < ref->nfiles; n++)
398 G_message("%s@%s", ref->file[n].name, ref->file[n].mapset);
399
400 G_fatal_error(_("Exit!"));
401 }
402
make_ipol_list(void)403 static char *make_ipol_list(void)
404 {
405 int size = 0;
406 int i;
407 char *buf;
408
409 for (i = 0; menu[i].name; i++)
410 size += strlen(menu[i].name) + 1;
411
412 buf = G_malloc(size);
413 *buf = '\0';
414
415 for (i = 0; menu[i].name; i++) {
416 if (i)
417 strcat(buf, ",");
418 strcat(buf, menu[i].name);
419 }
420
421 return buf;
422 }
423