1 /* rectification code */
2
3 /* Modified to support Grass 5.0 fp format 11 april 2000
4 *
5 * Pierre de Mouveaux - pmx@audiovu.com
6 *
7 */
8
9 #include <stdlib.h>
10 #include <string.h>
11 #include <unistd.h>
12 #include "global.h"
13
rectify(struct Image_Group * group,char * name,char * mapset,char * result,int order,char * interp_method)14 int rectify(struct Image_Group *group, char *name, char *mapset,
15 char *result, int order, char *interp_method)
16 {
17 struct Cell_head cellhd;
18 int ncols, nrows;
19 int row, col;
20 double row_idx, col_idx;
21 int infd, outfd;
22 RASTER_MAP_TYPE map_type;
23 int cell_size;
24 void *trast, *tptr;
25 double n1, e1, nx, ex;
26 struct cache *ibuffer;
27
28 select_current_env();
29 Rast_get_cellhd(name, mapset, &cellhd);
30
31 /* open the file to be rectified
32 * set window to cellhd first to be able to read file exactly
33 */
34 Rast_set_input_window(&cellhd);
35 infd = Rast_open_old(name, mapset);
36 map_type = Rast_get_map_type(infd);
37 cell_size = Rast_cell_size(map_type);
38
39 ibuffer = readcell(infd, seg_mb_img);
40
41 Rast_close(infd); /* (pmx) 17 april 2000 */
42
43 G_message(_("Rectify <%s@%s> (location <%s>)"),
44 name, mapset, G_location());
45 select_target_env();
46 G_message(_("into <%s@%s> (location <%s>) ..."),
47 result, G_mapset(), G_location());
48
49 nrows = target_window.rows;
50 ncols = target_window.cols;
51
52 if (strcmp(interp_method, "nearest") != 0) {
53 map_type = DCELL_TYPE;
54 cell_size = Rast_cell_size(map_type);
55 }
56
57 /* open the result file into target window
58 * this open must be first since we change the window later
59 * raster maps open for writing are not affected by window changes
60 * but those open for reading are
61 */
62
63 outfd = Rast_open_new(result, map_type);
64 trast = Rast_allocate_output_buf(map_type);
65
66 for (row = 0; row < nrows; row++) {
67 n1 = target_window.north - (row + 0.5) * target_window.ns_res;
68
69 G_percent(row, nrows, 2);
70
71 Rast_set_null_value(trast, ncols, map_type);
72 tptr = trast;
73 for (col = 0; col < ncols; col++) {
74 e1 = target_window.west + (col + 0.5) * target_window.ew_res;
75
76 /* backwards transformation of target cell center */
77 if (order == 0)
78 I_georef_tps(e1, n1, &ex, &nx, group->E21_t, group->N21_t, &group->control_points, 0);
79 else
80 I_georef(e1, n1, &ex, &nx, group->E21, group->N21, order);
81
82 /* convert to row/column indices of source raster */
83 row_idx = (cellhd.north - nx) / cellhd.ns_res;
84 col_idx = (ex - cellhd.west) / cellhd.ew_res;
85
86 /* resample data point */
87 interpolate(ibuffer, tptr, map_type, &row_idx, &col_idx, &cellhd);
88
89 tptr = G_incr_void_ptr(tptr, cell_size);
90 }
91 Rast_put_row(outfd, trast, map_type);
92 }
93 G_percent(1, 1, 1);
94
95 Rast_close(outfd); /* (pmx) 17 april 2000 */
96 G_free(trast);
97
98 close(ibuffer->fd);
99 release_cache(ibuffer);
100
101 Rast_get_cellhd(result, G_mapset(), &cellhd);
102
103 if (cellhd.proj == 0) { /* x,y imagery */
104 cellhd.proj = target_window.proj;
105 cellhd.zone = target_window.zone;
106 }
107
108 if (target_window.proj != cellhd.proj) {
109 cellhd.proj = target_window.proj;
110 G_warning(_("Raster map <%s@%s>: projection don't match current settings"),
111 name, mapset);
112 }
113
114 if (target_window.zone != cellhd.zone) {
115 cellhd.zone = target_window.zone;
116 G_warning(_("Raster map <%s@%s>: zone don't match current settings"),
117 name, mapset);
118 }
119
120 select_current_env();
121
122 return 1;
123 }
124