1
2 /****************************************************************************
3 *
4 * MODULE: r.resamp.interp
5 * AUTHOR(S): Glynn Clements <glynn gclements.plus.com> (original contributor),
6 * Paul Kelly <paul-grass stjohnspoint.co.uk>,
7 * Dylan Beaudette, Hamish Bowman <hamish_b yahoo.com>,
8 * Markus Metz (lanczos)
9 * PURPOSE:
10 * COPYRIGHT: (C) 2006-2007 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 /* TODO: add fallback methods, see e.g. r.proj */
19
20 #include <stdlib.h>
21 #include <string.h>
22 #include <math.h>
23 #include <grass/gis.h>
24 #include <grass/raster.h>
25 #include <grass/glocale.h>
26
27 static int neighbors;
28 static DCELL *bufs[5];
29 static int cur_row;
30
read_rows(int infile,int row)31 static void read_rows(int infile, int row)
32 {
33 int end_row = cur_row + neighbors;
34 int first_row = row;
35 int last_row = row + neighbors;
36 int offset = first_row - cur_row;
37 int keep = end_row - first_row;
38 int i;
39
40 if (end_row >= last_row)
41 return;
42
43 if (keep > 0 && offset > 0)
44 for (i = 0; i < keep; i++) {
45 DCELL *tmp = bufs[i];
46
47 bufs[i] = bufs[i + offset];
48 bufs[i + offset] = tmp;
49 }
50
51 if (keep < 0)
52 keep = 0;
53
54 for (i = keep; i < neighbors; i++)
55 Rast_get_d_row(infile, bufs[i], first_row + i);
56
57 cur_row = first_row;
58 }
59
main(int argc,char * argv[])60 int main(int argc, char *argv[])
61 {
62 struct GModule *module;
63 struct Option *rastin, *rastout, *method;
64 struct History history;
65 char title[64];
66 char buf_nsres[100], buf_ewres[100];
67 struct Colors colors;
68 int infile, outfile;
69 DCELL *outbuf;
70 int row, col;
71 struct Cell_head dst_w, src_w;
72
73 G_gisinit(argv[0]);
74
75 module = G_define_module();
76 G_add_keyword(_("raster"));
77 G_add_keyword(_("resample"));
78 G_add_keyword(_("interpolation"));
79 G_add_keyword(_("nearest neighbor"));
80 G_add_keyword(_("bilinear"));
81 G_add_keyword(_("bicubic"));
82 G_add_keyword(_("lanczos"));
83 module->description =
84 _("Resamples raster map to a finer grid using interpolation.");
85
86 rastin = G_define_standard_option(G_OPT_R_INPUT);
87 rastout = G_define_standard_option(G_OPT_R_OUTPUT);
88
89 method = G_define_standard_option(G_OPT_R_INTERP_TYPE);
90 method->options = "nearest,bilinear,bicubic,lanczos";
91 method->answer = "bilinear";
92 method->guisection = _("Method");
93
94 if (G_parser(argc, argv))
95 exit(EXIT_FAILURE);
96
97 if (G_strcasecmp(method->answer, "nearest") == 0)
98 neighbors = 1;
99 else if (G_strcasecmp(method->answer, "bilinear") == 0)
100 neighbors = 2;
101 else if (G_strcasecmp(method->answer, "bicubic") == 0)
102 neighbors = 4;
103 else if (G_strcasecmp(method->answer, "lanczos") == 0)
104 neighbors = 5;
105 else
106 G_fatal_error(_("Invalid method: %s"), method->answer);
107
108 G_get_set_window(&dst_w);
109
110 /* set window to old map */
111 Rast_get_cellhd(rastin->answer, "", &src_w);
112
113 if (G_projection() == PROJECTION_LL) {
114 /* try to shift source window to overlap with destination window */
115 while (src_w.west >= dst_w.east && src_w.east - 360.0 > dst_w.west) {
116 src_w.east -= 360.0;
117 src_w.west -= 360.0;
118 }
119 while (src_w.east <= dst_w.west && src_w.west + 360.0 < dst_w.east) {
120 src_w.east += 360.0;
121 src_w.west += 360.0;
122 }
123 }
124
125 /* enlarge source window */
126 {
127 double north = Rast_row_to_northing(0.5, &dst_w);
128 double south = Rast_row_to_northing(dst_w.rows - 0.5, &dst_w);
129 int r0 = (int)floor(Rast_northing_to_row(north, &src_w) - 0.5) - 2;
130 int r1 = (int)floor(Rast_northing_to_row(south, &src_w) - 0.5) + 3;
131 double west = Rast_col_to_easting(0.5, &dst_w);
132 double east = Rast_col_to_easting(dst_w.cols - 0.5, &dst_w);
133 /* do not use Rast_easting_to_col() because it does ll wrap */
134 /*
135 int c0 = (int)floor(Rast_easting_to_col(west, &src_w) - 0.5) - 2;
136 int c1 = (int)floor(Rast_easting_to_col(east, &src_w) - 0.5) + 3;
137 */
138 int c0 = (int)floor(((west - src_w.west) / src_w.ew_res) - 0.5) - 2;
139 int c1 = (int)floor(((east - src_w.west) / src_w.ew_res) - 0.5) + 3;
140
141 src_w.south -= src_w.ns_res * (r1 - src_w.rows);
142 src_w.north += src_w.ns_res * (-r0);
143 src_w.west -= src_w.ew_res * (-c0);
144 src_w.east += src_w.ew_res * (c1 - src_w.cols);
145 src_w.rows = r1 - r0;
146 src_w.cols = c1 - c0;
147 }
148
149 Rast_set_input_window(&src_w);
150
151 /* allocate buffers for input rows */
152 for (row = 0; row < neighbors; row++)
153 bufs[row] = Rast_allocate_d_input_buf();
154
155 cur_row = -100;
156
157 /* open old map */
158 infile = Rast_open_old(rastin->answer, "");
159
160 /* reset window to current region */
161 Rast_set_output_window(&dst_w);
162
163 outbuf = Rast_allocate_d_output_buf();
164
165 /* open new map */
166 outfile = Rast_open_new(rastout->answer, DCELL_TYPE);
167
168 switch (neighbors) {
169 case 1: /* nearest */
170 for (row = 0; row < dst_w.rows; row++) {
171 double north = Rast_row_to_northing(row + 0.5, &dst_w);
172 double maprow_f = Rast_northing_to_row(north, &src_w) - 0.5;
173 int maprow0 = (int)floor(maprow_f + 0.5);
174
175 G_percent(row, dst_w.rows, 2);
176
177 read_rows(infile, maprow0);
178
179 for (col = 0; col < dst_w.cols; col++) {
180 double east = Rast_col_to_easting(col + 0.5, &dst_w);
181 double mapcol_f = Rast_easting_to_col(east, &src_w) - 0.5;
182 int mapcol0 = (int)floor(mapcol_f + 0.5);
183
184 double c = bufs[0][mapcol0];
185
186 if (Rast_is_d_null_value(&c)) {
187 Rast_set_d_null_value(&outbuf[col], 1);
188 }
189 else {
190 outbuf[col] = c;
191 }
192 }
193
194 Rast_put_d_row(outfile, outbuf);
195 }
196 break;
197
198 case 2: /* bilinear */
199 for (row = 0; row < dst_w.rows; row++) {
200 double north = Rast_row_to_northing(row + 0.5, &dst_w);
201 double maprow_f = Rast_northing_to_row(north, &src_w) - 0.5;
202 int maprow0 = (int)floor(maprow_f);
203 double v = maprow_f - maprow0;
204
205 G_percent(row, dst_w.rows, 2);
206
207 read_rows(infile, maprow0);
208
209 for (col = 0; col < dst_w.cols; col++) {
210 double east = Rast_col_to_easting(col + 0.5, &dst_w);
211 double mapcol_f = Rast_easting_to_col(east, &src_w) - 0.5;
212 int mapcol0 = (int)floor(mapcol_f);
213 int mapcol1 = mapcol0 + 1;
214 double u = mapcol_f - mapcol0;
215
216 double c00 = bufs[0][mapcol0];
217 double c01 = bufs[0][mapcol1];
218 double c10 = bufs[1][mapcol0];
219 double c11 = bufs[1][mapcol1];
220
221 if (Rast_is_d_null_value(&c00) ||
222 Rast_is_d_null_value(&c01) ||
223 Rast_is_d_null_value(&c10) || Rast_is_d_null_value(&c11)) {
224 Rast_set_d_null_value(&outbuf[col], 1);
225 }
226 else {
227 outbuf[col] = Rast_interp_bilinear(u, v, c00, c01, c10, c11);
228 }
229 }
230
231 Rast_put_d_row(outfile, outbuf);
232 }
233 break;
234
235 case 4: /* bicubic */
236 for (row = 0; row < dst_w.rows; row++) {
237 double north = Rast_row_to_northing(row + 0.5, &dst_w);
238 double maprow_f = Rast_northing_to_row(north, &src_w) - 0.5;
239 int maprow1 = (int)floor(maprow_f);
240 int maprow0 = maprow1 - 1;
241 double v = maprow_f - maprow1;
242
243 G_percent(row, dst_w.rows, 2);
244
245 read_rows(infile, maprow0);
246
247 for (col = 0; col < dst_w.cols; col++) {
248 double east = Rast_col_to_easting(col + 0.5, &dst_w);
249 double mapcol_f = Rast_easting_to_col(east, &src_w) - 0.5;
250 int mapcol1 = (int)floor(mapcol_f);
251 int mapcol0 = mapcol1 - 1;
252 int mapcol2 = mapcol1 + 1;
253 int mapcol3 = mapcol1 + 2;
254 double u = mapcol_f - mapcol1;
255
256 double c00 = bufs[0][mapcol0];
257 double c01 = bufs[0][mapcol1];
258 double c02 = bufs[0][mapcol2];
259 double c03 = bufs[0][mapcol3];
260
261 double c10 = bufs[1][mapcol0];
262 double c11 = bufs[1][mapcol1];
263 double c12 = bufs[1][mapcol2];
264 double c13 = bufs[1][mapcol3];
265
266 double c20 = bufs[2][mapcol0];
267 double c21 = bufs[2][mapcol1];
268 double c22 = bufs[2][mapcol2];
269 double c23 = bufs[2][mapcol3];
270
271 double c30 = bufs[3][mapcol0];
272 double c31 = bufs[3][mapcol1];
273 double c32 = bufs[3][mapcol2];
274 double c33 = bufs[3][mapcol3];
275
276 if (Rast_is_d_null_value(&c00) ||
277 Rast_is_d_null_value(&c01) ||
278 Rast_is_d_null_value(&c02) ||
279 Rast_is_d_null_value(&c03) ||
280 Rast_is_d_null_value(&c10) ||
281 Rast_is_d_null_value(&c11) ||
282 Rast_is_d_null_value(&c12) ||
283 Rast_is_d_null_value(&c13) ||
284 Rast_is_d_null_value(&c20) ||
285 Rast_is_d_null_value(&c21) ||
286 Rast_is_d_null_value(&c22) ||
287 Rast_is_d_null_value(&c23) ||
288 Rast_is_d_null_value(&c30) ||
289 Rast_is_d_null_value(&c31) ||
290 Rast_is_d_null_value(&c32) || Rast_is_d_null_value(&c33)) {
291 Rast_set_d_null_value(&outbuf[col], 1);
292 }
293 else {
294 outbuf[col] = Rast_interp_bicubic(u, v,
295 c00, c01, c02, c03,
296 c10, c11, c12, c13,
297 c20, c21, c22, c23,
298 c30, c31, c32, c33);
299 }
300 }
301
302 Rast_put_d_row(outfile, outbuf);
303 }
304 break;
305
306 case 5: /* lanczos */
307 for (row = 0; row < dst_w.rows; row++) {
308 double north = Rast_row_to_northing(row + 0.5, &dst_w);
309 double maprow_f = Rast_northing_to_row(north, &src_w) - 0.5;
310 int maprow1 = (int)floor(maprow_f + 0.5);
311 int maprow0 = maprow1 - 2;
312 double v = maprow_f - maprow1;
313
314 G_percent(row, dst_w.rows, 2);
315
316 read_rows(infile, maprow0);
317
318 for (col = 0; col < dst_w.cols; col++) {
319 double east = Rast_col_to_easting(col + 0.5, &dst_w);
320 double mapcol_f = Rast_easting_to_col(east, &src_w) - 0.5;
321 int mapcol2 = (int)floor(mapcol_f + 0.5);
322 int mapcol0 = mapcol2 - 2;
323 int mapcol4 = mapcol2 + 2;
324 double u = mapcol_f - mapcol2;
325 double c[25];
326 int ci = 0, i, j, do_lanczos = 1;
327
328 for (i = 0; i < 5; i++) {
329 for (j = mapcol0; j <= mapcol4; j++) {
330 c[ci] = bufs[i][j];
331 if (Rast_is_d_null_value(&(c[ci]))) {
332 Rast_set_d_null_value(&outbuf[col], 1);
333 do_lanczos = 0;
334 break;
335 }
336 ci++;
337 }
338 if (!do_lanczos)
339 break;
340 }
341
342 if (do_lanczos) {
343 outbuf[col] = Rast_interp_lanczos(u, v, c);
344 }
345 }
346
347 Rast_put_d_row(outfile, outbuf);
348 }
349 break;
350 }
351
352 G_percent(dst_w.rows, dst_w.rows, 2);
353
354 Rast_close(infile);
355 Rast_close(outfile);
356
357
358 /* record map metadata/history info */
359 sprintf(title, "Resample by %s interpolation", method->answer);
360 Rast_put_cell_title(rastout->answer, title);
361
362 Rast_short_history(rastout->answer, "raster", &history);
363 Rast_set_history(&history, HIST_DATSRC_1, rastin->answer);
364 G_format_resolution(src_w.ns_res, buf_nsres, src_w.proj);
365 G_format_resolution(src_w.ew_res, buf_ewres, src_w.proj);
366 Rast_format_history(&history, HIST_DATSRC_2,
367 "Source map NS res: %s EW res: %s",
368 buf_nsres, buf_ewres);
369 Rast_command_history(&history);
370 Rast_write_history(rastout->answer, &history);
371
372 /* copy color table from source map */
373 if (Rast_read_colors(rastin->answer, "", &colors) < 0)
374 G_fatal_error(_("Unable to read color table for %s"), rastin->answer);
375 Rast_mark_colors_as_fp(&colors);
376 Rast_write_colors(rastout->answer, G_mapset(), &colors);
377
378 return (EXIT_SUCCESS);
379 }
380