← Back to branch summary

~ubuntu-branches/ubuntu/precise/grass/precise

« back to all changes in this revision

Viewing changes to imagery/i.ortho.photo/i.photo.rectify/main.c

  • Committer: Bazaar Package Importer
  • Author(s): Francesco Paolo Lovergine
  • Date: 2011-04-13 17:08:41 UTC
  • mfrom: (8.1.7 sid)
  • Revision ID: james.westby@ubuntu.com-20110413170841-ss1t9bic0d0uq0gz
Tags: 6.4.1-1
* New upstream version.
* Now build-dep on libjpeg-dev and current libreadline6-dev.
* Removed patch swig: obsolete.
* Policy bumped to 3.9.2, without changes.

Show diffs side-by-side

added added

removed removed

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