← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to raster/r.resamp.interp/main.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
raster/r.resamp.interp/main.c
4
4
 * MODULE:       r.resamp.interp
5
5
 * AUTHOR(S):    Glynn Clements <glynn gclements.plus.com> (original contributor),
6
6
 *               Paul Kelly <paul-grass stjohnspoint.co.uk>,
7
 
 *               Dylan Beaudette, Hamish Bowman <hamish_b yahoo.com>
 
7
 *               Dylan Beaudette, Hamish Bowman <hamish_b yahoo.com>,
 
8
 *               Markus Metz (lanczos)
8
9
 * PURPOSE:      
9
10
 * COPYRIGHT:    (C) 2006-2007 by the GRASS Development Team
10
11
 *
13
14
 *               for details.
14
15
 *
15
16
 *****************************************************************************/
 
17
 
 
18
/* TODO: add fallback methods, see e.g. r.proj */
 
19
 
16
20
#include <stdlib.h>
17
21
#include <string.h>
18
22
#include <math.h>
19
23
#include <grass/gis.h>
 
24
#include <grass/raster.h>
20
25
#include <grass/glocale.h>
21
26
 
22
27
static int neighbors;
23
 
static DCELL *bufs[4];
 
28
static DCELL *bufs[5];
24
29
static int cur_row;
25
30
 
26
31
static void read_rows(int infile, int row)
47
52
        keep = 0;
48
53
 
49
54
    for (i = keep; i < neighbors; i++)
50
 
        G_get_d_raster_row(infile, bufs[i], first_row + i);
 
55
        Rast_get_d_row(infile, bufs[i], first_row + i);
51
56
 
52
57
    cur_row = first_row;
53
58
}
60
65
    char title[64];
61
66
    char buf_nsres[100], buf_ewres[100];
62
67
    struct Colors colors;
63
 
    char *inmap;
64
68
    int infile, outfile;
65
69
    DCELL *outbuf;
66
70
    int row, col;
69
73
    G_gisinit(argv[0]);
70
74
 
71
75
    module = G_define_module();
72
 
    module->keywords = _("raster, resample");
 
76
    G_add_keyword(_("raster"));
 
77
    G_add_keyword(_("resample"));
 
78
    G_add_keyword(_("interpolation"));
73
79
    module->description =
74
 
        _("Resamples raster map layers to a finer grid using interpolation.");
 
80
        _("Resamples raster map to a finer grid using interpolation.");
75
81
 
76
82
    rastin = G_define_standard_option(G_OPT_R_INPUT);
77
83
    rastout = G_define_standard_option(G_OPT_R_OUTPUT);
78
84
 
79
 
    method = G_define_option();
80
 
    method->key = "method";
81
 
    method->type = TYPE_STRING;
82
 
    method->required = NO;
83
 
    method->description = _("Interpolation method");
84
 
    method->options = "nearest,bilinear,bicubic";
 
85
    method = G_define_standard_option(G_OPT_R_INTERP_TYPE);
 
86
    method->options = "nearest,bilinear,bicubic,lanczos";
85
87
    method->answer = "bilinear";
86
 
 
 
88
    method->guisection = _("Method");
 
89
    
87
90
    if (G_parser(argc, argv))
88
91
        exit(EXIT_FAILURE);
89
92
 
93
96
        neighbors = 2;
94
97
    else if (G_strcasecmp(method->answer, "bicubic") == 0)
95
98
        neighbors = 4;
 
99
    else if (G_strcasecmp(method->answer, "lanczos") == 0)
 
100
        neighbors = 5;
96
101
    else
97
102
        G_fatal_error(_("Invalid method: %s"), method->answer);
98
103
 
99
104
    G_get_set_window(&dst_w);
100
105
 
101
 
    inmap = G_find_cell2(rastin->answer, "");
102
 
    if (!inmap)
103
 
        G_fatal_error(_("Raster map <%s> not found"), rastin->answer);
104
 
 
105
106
    /* set window to old map */
106
 
    G_get_cellhd(rastin->answer, inmap, &src_w);
 
107
    Rast_get_cellhd(rastin->answer, "", &src_w);
107
108
 
108
109
    /* enlarge source window */
109
110
    {
110
 
        double north = G_row_to_northing(0.5, &dst_w);
111
 
        double south = G_row_to_northing(dst_w.rows - 0.5, &dst_w);
112
 
        int r0 = (int)floor(G_northing_to_row(north, &src_w) - 0.5) - 1;
113
 
        int r1 = (int)floor(G_northing_to_row(south, &src_w) - 0.5) + 3;
114
 
        double west = G_col_to_easting(0.5, &dst_w);
115
 
        double east = G_col_to_easting(dst_w.cols - 0.5, &dst_w);
116
 
        int c0 = (int)floor(G_easting_to_col(west, &src_w) - 0.5) - 1;
117
 
        int c1 = (int)floor(G_easting_to_col(east, &src_w) - 0.5) + 3;
 
111
        double north = Rast_row_to_northing(0.5, &dst_w);
 
112
        double south = Rast_row_to_northing(dst_w.rows - 0.5, &dst_w);
 
113
        int r0 = (int)floor(Rast_northing_to_row(north, &src_w) - 0.5) - 2;
 
114
        int r1 = (int)floor(Rast_northing_to_row(south, &src_w) - 0.5) + 3;
 
115
        double west = Rast_col_to_easting(0.5, &dst_w);
 
116
        double east = Rast_col_to_easting(dst_w.cols - 0.5, &dst_w);
 
117
        int c0 = (int)floor(Rast_easting_to_col(west, &src_w) - 0.5) - 2;
 
118
        int c1 = (int)floor(Rast_easting_to_col(east, &src_w) - 0.5) + 3;
118
119
 
119
120
        src_w.south -= src_w.ns_res * (r1 - src_w.rows);
120
121
        src_w.north += src_w.ns_res * (-r0);
124
125
        src_w.cols = c1 - c0;
125
126
    }
126
127
 
127
 
    G_set_window(&src_w);
 
128
    Rast_set_input_window(&src_w);
128
129
 
129
130
    /* allocate buffers for input rows */
130
131
    for (row = 0; row < neighbors; row++)
131
 
        bufs[row] = G_allocate_d_raster_buf();
 
132
        bufs[row] = Rast_allocate_d_input_buf();
132
133
 
133
134
    cur_row = -100;
134
135
 
135
136
    /* open old map */
136
 
    infile = G_open_cell_old(rastin->answer, inmap);
137
 
    if (infile < 0)
138
 
        G_fatal_error(_("Unable to open raster map <%s>"), rastin->answer);
 
137
    infile = Rast_open_old(rastin->answer, "");
139
138
 
140
139
    /* reset window to current region */
141
 
    G_set_window(&dst_w);
 
140
    Rast_set_output_window(&dst_w);
142
141
 
143
 
    outbuf = G_allocate_d_raster_buf();
 
142
    outbuf = Rast_allocate_d_output_buf();
144
143
 
145
144
    /* open new map */
146
 
    outfile = G_open_raster_new(rastout->answer, DCELL_TYPE);
147
 
    if (outfile < 0)
148
 
        G_fatal_error(_("Unable to create raster map <%s>"), rastout->answer);
149
 
 
150
 
    G_suppress_warnings(1);
151
 
    /* otherwise get complaints about window changes */
 
145
    outfile = Rast_open_new(rastout->answer, DCELL_TYPE);
152
146
 
153
147
    switch (neighbors) {
154
148
    case 1:                     /* nearest */
155
149
        for (row = 0; row < dst_w.rows; row++) {
156
 
            double north = G_row_to_northing(row + 0.5, &dst_w);
157
 
            double maprow_f = G_northing_to_row(north, &src_w) - 0.5;
 
150
            double north = Rast_row_to_northing(row + 0.5, &dst_w);
 
151
            double maprow_f = Rast_northing_to_row(north, &src_w) - 0.5;
158
152
            int maprow0 = (int)floor(maprow_f + 0.5);
159
153
 
160
154
            G_percent(row, dst_w.rows, 2);
161
155
 
162
 
            G_set_window(&src_w);
163
156
            read_rows(infile, maprow0);
164
157
 
165
158
            for (col = 0; col < dst_w.cols; col++) {
166
 
                double east = G_col_to_easting(col + 0.5, &dst_w);
167
 
                double mapcol_f = G_easting_to_col(east, &src_w) - 0.5;
 
159
                double east = Rast_col_to_easting(col + 0.5, &dst_w);
 
160
                double mapcol_f = Rast_easting_to_col(east, &src_w) - 0.5;
168
161
                int mapcol0 = (int)floor(mapcol_f + 0.5);
169
162
 
170
163
                double c = bufs[0][mapcol0];
171
164
 
172
 
                if (G_is_d_null_value(&c)) {
173
 
                    G_set_d_null_value(&outbuf[col], 1);
 
165
                if (Rast_is_d_null_value(&c)) {
 
166
                    Rast_set_d_null_value(&outbuf[col], 1);
174
167
                }
175
168
                else {
176
169
                    outbuf[col] = c;
177
170
                }
178
171
            }
179
172
 
180
 
            G_set_window(&dst_w);
181
 
            G_put_d_raster_row(outfile, outbuf);
 
173
            Rast_put_d_row(outfile, outbuf);
182
174
        }
183
175
        break;
184
176
 
185
177
    case 2:                     /* bilinear */
186
178
        for (row = 0; row < dst_w.rows; row++) {
187
 
            double north = G_row_to_northing(row + 0.5, &dst_w);
188
 
            double maprow_f = G_northing_to_row(north, &src_w) - 0.5;
 
179
            double north = Rast_row_to_northing(row + 0.5, &dst_w);
 
180
            double maprow_f = Rast_northing_to_row(north, &src_w) - 0.5;
189
181
            int maprow0 = (int)floor(maprow_f);
190
182
            double v = maprow_f - maprow0;
191
183
 
192
184
            G_percent(row, dst_w.rows, 2);
193
185
 
194
 
            G_set_window(&src_w);
195
186
            read_rows(infile, maprow0);
196
187
 
197
188
            for (col = 0; col < dst_w.cols; col++) {
198
 
                double east = G_col_to_easting(col + 0.5, &dst_w);
199
 
                double mapcol_f = G_easting_to_col(east, &src_w) - 0.5;
 
189
                double east = Rast_col_to_easting(col + 0.5, &dst_w);
 
190
                double mapcol_f = Rast_easting_to_col(east, &src_w) - 0.5;
200
191
                int mapcol0 = (int)floor(mapcol_f);
201
192
                int mapcol1 = mapcol0 + 1;
202
193
                double u = mapcol_f - mapcol0;
206
197
                double c10 = bufs[1][mapcol0];
207
198
                double c11 = bufs[1][mapcol1];
208
199
 
209
 
                if (G_is_d_null_value(&c00) ||
210
 
                    G_is_d_null_value(&c01) ||
211
 
                    G_is_d_null_value(&c10) || G_is_d_null_value(&c11)) {
212
 
                    G_set_d_null_value(&outbuf[col], 1);
 
200
                if (Rast_is_d_null_value(&c00) ||
 
201
                    Rast_is_d_null_value(&c01) ||
 
202
                    Rast_is_d_null_value(&c10) || Rast_is_d_null_value(&c11)) {
 
203
                    Rast_set_d_null_value(&outbuf[col], 1);
213
204
                }
214
205
                else {
215
 
                    outbuf[col] = G_interp_bilinear(u, v, c00, c01, c10, c11);
 
206
                    outbuf[col] = Rast_interp_bilinear(u, v, c00, c01, c10, c11);
216
207
                }
217
208
            }
218
209
 
219
 
            G_set_window(&dst_w);
220
 
            G_put_d_raster_row(outfile, outbuf);
 
210
            Rast_put_d_row(outfile, outbuf);
221
211
        }
222
212
        break;
223
213
 
224
214
    case 4:                     /* bicubic */
225
215
        for (row = 0; row < dst_w.rows; row++) {
226
 
            double north = G_row_to_northing(row + 0.5, &dst_w);
227
 
            double maprow_f = G_northing_to_row(north, &src_w) - 0.5;
 
216
            double north = Rast_row_to_northing(row + 0.5, &dst_w);
 
217
            double maprow_f = Rast_northing_to_row(north, &src_w) - 0.5;
228
218
            int maprow1 = (int)floor(maprow_f);
229
219
            int maprow0 = maprow1 - 1;
230
220
            double v = maprow_f - maprow1;
231
221
 
232
222
            G_percent(row, dst_w.rows, 2);
233
223
 
234
 
            G_set_window(&src_w);
235
224
            read_rows(infile, maprow0);
236
225
 
237
226
            for (col = 0; col < dst_w.cols; col++) {
238
 
                double east = G_col_to_easting(col + 0.5, &dst_w);
239
 
                double mapcol_f = G_easting_to_col(east, &src_w) - 0.5;
 
227
                double east = Rast_col_to_easting(col + 0.5, &dst_w);
 
228
                double mapcol_f = Rast_easting_to_col(east, &src_w) - 0.5;
240
229
                int mapcol1 = (int)floor(mapcol_f);
241
230
                int mapcol0 = mapcol1 - 1;
242
231
                int mapcol2 = mapcol1 + 1;
263
252
                double c32 = bufs[3][mapcol2];
264
253
                double c33 = bufs[3][mapcol3];
265
254
 
266
 
                if (G_is_d_null_value(&c00) ||
267
 
                    G_is_d_null_value(&c01) ||
268
 
                    G_is_d_null_value(&c02) ||
269
 
                    G_is_d_null_value(&c03) ||
270
 
                    G_is_d_null_value(&c10) ||
271
 
                    G_is_d_null_value(&c11) ||
272
 
                    G_is_d_null_value(&c12) ||
273
 
                    G_is_d_null_value(&c13) ||
274
 
                    G_is_d_null_value(&c20) ||
275
 
                    G_is_d_null_value(&c21) ||
276
 
                    G_is_d_null_value(&c22) ||
277
 
                    G_is_d_null_value(&c23) ||
278
 
                    G_is_d_null_value(&c30) ||
279
 
                    G_is_d_null_value(&c31) ||
280
 
                    G_is_d_null_value(&c32) || G_is_d_null_value(&c33)) {
281
 
                    G_set_d_null_value(&outbuf[col], 1);
 
255
                if (Rast_is_d_null_value(&c00) ||
 
256
                    Rast_is_d_null_value(&c01) ||
 
257
                    Rast_is_d_null_value(&c02) ||
 
258
                    Rast_is_d_null_value(&c03) ||
 
259
                    Rast_is_d_null_value(&c10) ||
 
260
                    Rast_is_d_null_value(&c11) ||
 
261
                    Rast_is_d_null_value(&c12) ||
 
262
                    Rast_is_d_null_value(&c13) ||
 
263
                    Rast_is_d_null_value(&c20) ||
 
264
                    Rast_is_d_null_value(&c21) ||
 
265
                    Rast_is_d_null_value(&c22) ||
 
266
                    Rast_is_d_null_value(&c23) ||
 
267
                    Rast_is_d_null_value(&c30) ||
 
268
                    Rast_is_d_null_value(&c31) ||
 
269
                    Rast_is_d_null_value(&c32) || Rast_is_d_null_value(&c33)) {
 
270
                    Rast_set_d_null_value(&outbuf[col], 1);
282
271
                }
283
272
                else {
284
 
                    outbuf[col] = G_interp_bicubic(u, v,
 
273
                    outbuf[col] = Rast_interp_bicubic(u, v,
285
274
                                                   c00, c01, c02, c03,
286
275
                                                   c10, c11, c12, c13,
287
276
                                                   c20, c21, c22, c23,
289
278
                }
290
279
            }
291
280
 
292
 
            G_set_window(&dst_w);
293
 
            G_put_d_raster_row(outfile, outbuf);
 
281
            Rast_put_d_row(outfile, outbuf);
 
282
        }
 
283
        break;
 
284
 
 
285
    case 5:                     /* lanczos */
 
286
        for (row = 0; row < dst_w.rows; row++) {
 
287
            double north = Rast_row_to_northing(row + 0.5, &dst_w);
 
288
            double maprow_f = Rast_northing_to_row(north, &src_w) - 0.5;
 
289
            int maprow1 = (int)floor(maprow_f + 0.5);
 
290
            int maprow0 = maprow1 - 2;
 
291
            double v = maprow_f - maprow1;
 
292
 
 
293
            G_percent(row, dst_w.rows, 2);
 
294
 
 
295
            read_rows(infile, maprow0);
 
296
 
 
297
            for (col = 0; col < dst_w.cols; col++) {
 
298
                double east = Rast_col_to_easting(col + 0.5, &dst_w);
 
299
                double mapcol_f = Rast_easting_to_col(east, &src_w) - 0.5;
 
300
                int mapcol2 = (int)floor(mapcol_f + 0.5);
 
301
                int mapcol0 = mapcol2 - 2;
 
302
                int mapcol4 = mapcol2 + 2;
 
303
                double u = mapcol_f - mapcol2;
 
304
                double c[25];
 
305
                int ci = 0, i, j, do_lanczos = 1;
 
306
 
 
307
                for (i = 0; i < 5; i++) {
 
308
                    for (j = mapcol0; j <= mapcol4; j++) {
 
309
                        c[ci] = bufs[i][j];
 
310
                        if (Rast_is_d_null_value(&(c[ci]))) {
 
311
                            Rast_set_d_null_value(&outbuf[col], 1);
 
312
                            do_lanczos = 0;
 
313
                            break;
 
314
                        }
 
315
                        ci++;
 
316
                    }
 
317
                    if (!do_lanczos)
 
318
                        break;
 
319
                }
 
320
 
 
321
                if (do_lanczos) {
 
322
                    outbuf[col] = Rast_interp_lanczos(u, v, c);
 
323
                }
 
324
            }
 
325
 
 
326
            Rast_put_d_row(outfile, outbuf);
294
327
        }
295
328
        break;
296
329
    }
297
330
 
298
331
    G_percent(dst_w.rows, dst_w.rows, 2);
299
332
 
300
 
    G_close_cell(infile);
301
 
    G_close_cell(outfile);
 
333
    Rast_close(infile);
 
334
    Rast_close(outfile);
302
335
 
303
336
 
304
337
    /* record map metadata/history info */
305
338
    sprintf(title, "Resample by %s interpolation", method->answer);
306
 
    G_put_cell_title(rastout->answer, title);
 
339
    Rast_put_cell_title(rastout->answer, title);
307
340
 
308
 
    G_short_history(rastout->answer, "raster", &history);
309
 
    strncpy(history.datsrc_1, rastin->answer, RECORD_LEN);
310
 
    history.datsrc_1[RECORD_LEN - 1] = '\0';    /* strncpy() doesn't null terminate if maxfill */
 
341
    Rast_short_history(rastout->answer, "raster", &history);
 
342
    Rast_set_history(&history, HIST_DATSRC_1, rastin->answer);
311
343
    G_format_resolution(src_w.ns_res, buf_nsres, src_w.proj);
312
344
    G_format_resolution(src_w.ew_res, buf_ewres, src_w.proj);
313
 
    sprintf(history.datsrc_2, "Source map NS res: %s   EW res: %s", buf_nsres,
314
 
            buf_ewres);
315
 
    G_command_history(&history);
316
 
    G_write_history(rastout->answer, &history);
 
345
    Rast_format_history(&history, HIST_DATSRC_2,
 
346
                        "Source map NS res: %s   EW res: %s",
 
347
                        buf_nsres, buf_ewres);
 
348
    Rast_command_history(&history);
 
349
    Rast_write_history(rastout->answer, &history);
317
350
 
318
351
    /* copy color table from source map */
319
 
    if (G_read_colors(rastin->answer, inmap, &colors) < 0)
 
352
    if (Rast_read_colors(rastin->answer, "", &colors) < 0)
320
353
        G_fatal_error(_("Unable to read color table for %s"), rastin->answer);
321
 
    G_mark_colors_as_fp(&colors);
322
 
    if (G_write_colors(rastout->answer, G_mapset(), &colors) < 0)
323
 
        G_fatal_error(_("Unable to write color table for %s"),
324
 
                      rastout->answer);
 
354
    Rast_mark_colors_as_fp(&colors);
 
355
    Rast_write_colors(rastout->answer, G_mapset(), &colors);
325
356
 
326
357
    return (EXIT_SUCCESS);
327
358
}