← Back to branch summary

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

« back to all changes in this revision

Viewing changes to raster/r.resamp.filter/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.filter/main.c
 
1
 
 
2
/****************************************************************************
 
3
 *
 
4
 * MODULE:       r.resamp.filter
 
5
 * AUTHOR(S):    Glynn Clements <glynn gclements.plus.com>
 
6
 * PURPOSE:      
 
7
 * COPYRIGHT:    (C) 2010 by Glynn Clements and the GRASS Development Team
 
8
 *
 
9
 *               This program is free software under the GNU General Public
 
10
 *               License (>=v2). Read the file COPYING that comes with GRASS
 
11
 *               for details.
 
12
 *
 
13
 *****************************************************************************/
 
14
 
 
15
#include <stdlib.h>
 
16
#include <string.h>
 
17
#include <limits.h>
 
18
#include <math.h>
 
19
#include <grass/gis.h>
 
20
#include <grass/raster.h>
 
21
#include <grass/glocale.h>
 
22
 
 
23
static double f_box(double x)
 
24
{
 
25
    return (x > 1) ? 0
 
26
        : 1;
 
27
}
 
28
 
 
29
static double f_bartlett(double x)
 
30
{
 
31
    return (x > 1) ? 0
 
32
        : 1 - x;
 
33
}
 
34
 
 
35
static double f_hermite(double x)
 
36
{
 
37
    return (x > 1) ? 0
 
38
        : (2 * x - 3) * x * x + 1;
 
39
}
 
40
 
 
41
static double f_gauss(double x)
 
42
{
 
43
    return exp(-2 * x * x) * sqrt(2 / M_PI);
 
44
}
 
45
 
 
46
static double f_normal(double x)
 
47
{
 
48
    return f_gauss(x/2) / 2;
 
49
}
 
50
 
 
51
static double f_sinc(double x)
 
52
{
 
53
    return (x == 0) ? 1 : sin(M_PI * x) / (M_PI * x);
 
54
}
 
55
 
 
56
static double lanczos(double x, int a)
 
57
{
 
58
    return (x > a) ? 0
 
59
        : f_sinc(x) * f_sinc(x / a);
 
60
}
 
61
 
 
62
static double f_lanczos1(double x)
 
63
{
 
64
    return lanczos(x, 1);
 
65
}
 
66
 
 
67
static double f_lanczos2(double x)
 
68
{
 
69
    return lanczos(x, 2);
 
70
}
 
71
 
 
72
static double f_lanczos3(double x)
 
73
{
 
74
    return lanczos(x, 3);
 
75
}
 
76
 
 
77
static double f_hann(double x)
 
78
{
 
79
    return cos(M_PI * x) / 2 + 0.5;
 
80
}
 
81
 
 
82
static double f_hamming(double x)
 
83
{
 
84
    return 0.46 * cos(M_PI * x) + 0.54;
 
85
}
 
86
 
 
87
 
 
88
static double f_blackman(double x)
 
89
{
 
90
    return cos(M_PI * x) / 2 + 0.08 * cos(2 * M_PI * x) + 0.42;
 
91
}
 
92
 
 
93
struct filter_type {
 
94
    const char *name;
 
95
    double (*func)(double);
 
96
    int radius;
 
97
};
 
98
 
 
99
static const struct filter_type menu[] = {
 
100
    {"box",       f_box,       1},
 
101
    {"bartlett",  f_bartlett,  1},
 
102
    {"gauss",     f_gauss,     0},
 
103
    {"normal",    f_normal,    0},
 
104
    {"hermite",   f_hermite,   1},
 
105
    {"sinc",      f_sinc,      0},
 
106
    {"lanczos1",  f_lanczos1,  1},
 
107
    {"lanczos2",  f_lanczos2,  2},
 
108
    {"lanczos3",  f_lanczos3,  3},
 
109
    {"hann",      f_hann,      0},
 
110
    {"hamming",   f_hamming,   0},
 
111
    {"blackman",  f_blackman,  0},
 
112
    {NULL},
 
113
};
 
114
 
 
115
static char *build_filter_list(void)
 
116
{
 
117
    char *buf = G_malloc(1024);
 
118
    char *p = buf;
 
119
    int i;
 
120
 
 
121
    for (i = 0; menu[i].name; i++) {
 
122
        const char *q;
 
123
 
 
124
        if (i)
 
125
            *p++ = ',';
 
126
        for (q = menu[i].name; *q; p++, q++)
 
127
            *p = *q;
 
128
    }
 
129
    *p = '\0';
 
130
 
 
131
    return buf;
 
132
}
 
133
 
 
134
static const struct filter_type *find_method(const char *name)
 
135
{
 
136
    int i;
 
137
 
 
138
    for (i = 0; menu[i].name; i++)
 
139
        if (strcmp(menu[i].name, name) == 0)
 
140
            return &menu[i];
 
141
 
 
142
    G_fatal_error(_("Filter <%s> not found"), name);
 
143
 
 
144
    return NULL;
 
145
}
 
146
 
 
147
struct filter {
 
148
    double (*func)(double);
 
149
    double x_radius;
 
150
    double y_radius;
 
151
};
 
152
 
 
153
#define MAX_FILTERS 8
 
154
 
 
155
static int infile, outfile;
 
156
static struct filter filters[MAX_FILTERS];
 
157
static int num_filters;
 
158
static int nulls;
 
159
static struct Cell_head dst_w, src_w;
 
160
static double f_x_radius, f_y_radius;
 
161
static int row_scale, col_scale;
 
162
static DCELL *inbuf;
 
163
static DCELL *outbuf;
 
164
static DCELL **bufs;
 
165
static double *h_weights;
 
166
static double *v_weights;
 
167
static int *mapcol0, *mapcol1;
 
168
static int *maprow0, *maprow1;
 
169
 
 
170
static void make_h_weights(void)
 
171
{
 
172
    int col;
 
173
 
 
174
    h_weights = G_malloc(dst_w.cols * col_scale * sizeof(double));
 
175
    mapcol0 = G_malloc(dst_w.cols * sizeof(int));
 
176
    mapcol1 = G_malloc(dst_w.cols * sizeof(int));
 
177
 
 
178
    for (col = 0; col < dst_w.cols; col++) {
 
179
        double dx = Rast_col_to_easting(col + 0.5, &dst_w);
 
180
        double x0 = Rast_easting_to_col(dx - f_x_radius, &src_w);
 
181
        double x1 = Rast_easting_to_col(dx + f_x_radius, &src_w);
 
182
        int col0 = (int)floor(x0);
 
183
        int col1 = (int)floor(x1) + 1;
 
184
        int cols = col1 - col0;
 
185
        int j;
 
186
 
 
187
        mapcol0[col] = col0;
 
188
        mapcol1[col] = col1;
 
189
 
 
190
        for (j = 0; j < cols; j++) {
 
191
            double sx = Rast_col_to_easting(col0 + j + 0.5, &src_w);
 
192
            double r = fabs(sx - dx);
 
193
            double w = 1.0;
 
194
            int k;
 
195
 
 
196
            for (k = 0; k < num_filters; k++)
 
197
                w *= (*filters[k].func)(r / filters[k].x_radius);
 
198
 
 
199
            h_weights[col * col_scale + j] = w;
 
200
        }
 
201
 
 
202
        for (j = cols; j < col_scale; j++)
 
203
            h_weights[col * col_scale + j] = 0;
 
204
    }
 
205
}
 
206
 
 
207
static void make_v_weights(void)
 
208
{
 
209
    int row;
 
210
 
 
211
    v_weights = G_malloc(dst_w.rows * row_scale * sizeof(double));
 
212
    maprow0 = G_malloc(dst_w.rows * sizeof(int));
 
213
    maprow1 = G_malloc(dst_w.rows * sizeof(int));
 
214
 
 
215
    for (row = 0; row < dst_w.rows; row++) {
 
216
        double dy = Rast_row_to_northing(row + 0.5, &dst_w);
 
217
        double y0 = Rast_northing_to_row(dy + f_y_radius, &src_w);
 
218
        double y1 = Rast_northing_to_row(dy - f_y_radius, &src_w);
 
219
        int row0 = (int)floor(y0);
 
220
        int row1 = (int)floor(y1) + 1;
 
221
        int rows = row1 - row0;
 
222
        int i;
 
223
 
 
224
        maprow0[row] = row0;
 
225
        maprow1[row] = row1;
 
226
 
 
227
        for (i = 0; i < rows; i++) {
 
228
            double sy = Rast_row_to_northing(row0 + i + 0.5, &src_w);
 
229
            double r = fabs(sy - dy);
 
230
            double w = 1.0;
 
231
            int k;
 
232
 
 
233
            for (k = 0; k < num_filters; k++)
 
234
                w *= (*filters[k].func)(r / filters[k].y_radius);
 
235
 
 
236
            v_weights[row * row_scale + i] = w;
 
237
        }
 
238
 
 
239
        for (i = rows; i < row_scale; i++)
 
240
            v_weights[row * row_scale + i] = 0;
 
241
    }
 
242
}
 
243
 
 
244
static void h_filter(DCELL *dst, const DCELL *src)
 
245
{
 
246
    int col;
 
247
 
 
248
    for (col = 0; col < dst_w.cols; col++) {
 
249
        int col0 = mapcol0[col];
 
250
        int col1 = mapcol1[col];
 
251
        int cols = col1 - col0;
 
252
        double numer = 0.0;
 
253
        double denom = 0.0;
 
254
        int null = 0;
 
255
        int j;
 
256
 
 
257
        for (j = 0; j < cols; j++) {
 
258
            double w = h_weights[col * col_scale + j];
 
259
            const DCELL *c = &src[col0 + j];
 
260
 
 
261
            if (Rast_is_d_null_value(c)) {
 
262
                if (nulls) {
 
263
                    null = 1;
 
264
                    break;
 
265
                }
 
266
            }
 
267
            else {
 
268
                numer += w * (*c);
 
269
                denom += w;
 
270
            }
 
271
        }
 
272
 
 
273
        if (null || denom == 0)
 
274
            Rast_set_d_null_value(&dst[col], 1);
 
275
        else
 
276
            dst[col] = numer / denom;
 
277
    }
 
278
}
 
279
 
 
280
static void v_filter(DCELL *dst, DCELL **src, int row, int rows)
 
281
{
 
282
    int col;
 
283
 
 
284
    for (col = 0; col < dst_w.cols; col++) {
 
285
        double numer = 0.0;
 
286
        double denom = 0.0;
 
287
        int null = 0;
 
288
        int i;
 
289
 
 
290
        for (i = 0; i < rows; i++) {
 
291
            double w = v_weights[row * row_scale + i];
 
292
            const DCELL *c = &src[i][col];
 
293
 
 
294
            if (Rast_is_d_null_value(c)) {
 
295
                if (nulls) {
 
296
                    null = 1;
 
297
                    break;
 
298
                }
 
299
            }
 
300
            else {
 
301
                numer += w * (*c);
 
302
                denom += w;
 
303
            }
 
304
        }
 
305
 
 
306
        if (null || denom == 0)
 
307
            Rast_set_d_null_value(&dst[col], 1);
 
308
        else
 
309
            dst[col] = numer / denom;
 
310
    }
 
311
}
 
312
 
 
313
static void filter(void)
 
314
{
 
315
    int cur_row;
 
316
    int num_rows = 0;
 
317
    int row;
 
318
 
 
319
    make_h_weights();
 
320
    make_v_weights();
 
321
 
 
322
    for (row = 0; row < dst_w.rows; row++) {
 
323
        int row0 = maprow0[row];
 
324
        int row1 = maprow1[row];
 
325
        int rows = row1 - row0;
 
326
        int i;
 
327
 
 
328
        G_percent(row, dst_w.rows, 2);
 
329
 
 
330
        if (row0 >= cur_row && row0 < cur_row + num_rows) {
 
331
            int m = row0 - cur_row;
 
332
            int n = cur_row + num_rows - row0;
 
333
            int i;
 
334
 
 
335
            for (i = 0; i < n; i++) {
 
336
                DCELL *tmp = bufs[i];
 
337
                bufs[i] = bufs[m + i];
 
338
                bufs[m + i] = tmp;
 
339
            }
 
340
 
 
341
            cur_row = row0;
 
342
            num_rows = n;
 
343
        }
 
344
        else {
 
345
            cur_row = row0;
 
346
            num_rows = 0;
 
347
        }
 
348
 
 
349
        for (i = num_rows; i < rows; i++) {
 
350
            G_debug(5, "read: %p = %d", bufs[i], row0 + i);
 
351
            Rast_get_d_row(infile, inbuf, row0 + i);
 
352
            h_filter(bufs[i], inbuf);
 
353
        }
 
354
 
 
355
        num_rows = rows;
 
356
 
 
357
        v_filter(outbuf, bufs, row, rows);
 
358
 
 
359
        Rast_put_d_row(outfile, outbuf);
 
360
        G_debug(5, "write: %d", row);
 
361
    }
 
362
 
 
363
    G_percent(dst_w.rows, dst_w.rows, 2);
 
364
}
 
365
 
 
366
int main(int argc, char *argv[])
 
367
{
 
368
    struct GModule *module;
 
369
    struct
 
370
    {
 
371
        struct Option *rastin, *rastout, *method,
 
372
            *radius, *x_radius, *y_radius;
 
373
    } parm;
 
374
    struct
 
375
    {
 
376
        struct Flag *nulls;
 
377
    } flag;
 
378
    char title[64];
 
379
    int i;
 
380
 
 
381
    G_gisinit(argv[0]);
 
382
 
 
383
    module = G_define_module();
 
384
    G_add_keyword(_("raster"));
 
385
    G_add_keyword(_("resample"));
 
386
    G_add_keyword(_("kernel filter"));
 
387
    module->description =
 
388
        _("Resamples raster map layers using an analytic kernel.");
 
389
 
 
390
    parm.rastin = G_define_standard_option(G_OPT_R_INPUT);
 
391
 
 
392
    parm.rastout = G_define_standard_option(G_OPT_R_OUTPUT);
 
393
 
 
394
    parm.method = G_define_option();
 
395
    parm.method->key = "filter";
 
396
    parm.method->type = TYPE_STRING;
 
397
    parm.method->required = YES;
 
398
    parm.method->multiple = YES;
 
399
    parm.method->description = _("Filter kernel(s)");
 
400
    parm.method->options = build_filter_list();
 
401
 
 
402
    parm.radius = G_define_option();
 
403
    parm.radius->key = "radius";
 
404
    parm.radius->type = TYPE_DOUBLE;
 
405
    parm.radius->required = NO;
 
406
    parm.radius->multiple = YES;
 
407
    parm.radius->description = _("Filter radius");
 
408
 
 
409
    parm.x_radius = G_define_option();
 
410
    parm.x_radius->key = "x_radius";
 
411
    parm.x_radius->type = TYPE_DOUBLE;
 
412
    parm.x_radius->required = NO;
 
413
    parm.x_radius->multiple = YES;
 
414
    parm.x_radius->description = _("Filter radius (horizontal)");
 
415
 
 
416
    parm.y_radius = G_define_option();
 
417
    parm.y_radius->key = "y_radius";
 
418
    parm.y_radius->type = TYPE_DOUBLE;
 
419
    parm.y_radius->required = NO;
 
420
    parm.y_radius->multiple = YES;
 
421
    parm.y_radius->description = _("Filter radius (vertical)");
 
422
 
 
423
    flag.nulls = G_define_flag();
 
424
    flag.nulls->key = 'n';
 
425
    flag.nulls->description = _("Propagate NULLs");
 
426
 
 
427
    if (G_parser(argc, argv))
 
428
        exit(EXIT_FAILURE);
 
429
 
 
430
    if (parm.radius->answer) {
 
431
        if (parm.x_radius->answer || parm.y_radius->answer)
 
432
            G_fatal_error(_("%s= and %s=/%s= are mutually exclusive"),
 
433
                            parm.radius->key, parm.x_radius->key, parm.y_radius->key);
 
434
    }
 
435
    else {
 
436
        if (!parm.x_radius->answer && !parm.y_radius->answer)
 
437
            G_fatal_error(_("Either %s= or %s=/%s= required"),
 
438
                            parm.radius->key, parm.x_radius->key, parm.y_radius->key);
 
439
        if (!parm.x_radius->answer || !parm.y_radius->answer)
 
440
            G_fatal_error(_("Both %s= and %s= required"),
 
441
                            parm.x_radius->key, parm.y_radius->key);
 
442
    }
 
443
 
 
444
    nulls = flag.nulls->answer;
 
445
 
 
446
    f_x_radius = f_y_radius = 1e100;
 
447
 
 
448
    for (i = 0; ; i++) {
 
449
        const char *filter_arg = parm.method->answers[i];
 
450
        const char *x_radius_arg = parm.radius->answer
 
451
            ? parm.radius->answers[i]
 
452
            : parm.x_radius->answers[i];
 
453
        const char *y_radius_arg = parm.radius->answer
 
454
            ? parm.radius->answers[i]
 
455
            : parm.y_radius->answers[i];
 
456
        const struct filter_type *type;
 
457
        struct filter *filter;
 
458
 
 
459
        if (!filter_arg && !x_radius_arg && !y_radius_arg)
 
460
            break;
 
461
 
 
462
        if (!filter_arg || !x_radius_arg || !y_radius_arg)
 
463
            G_fatal_error(_("Differing number of values for filter= and [xy_]radius="));
 
464
 
 
465
        if (num_filters >= MAX_FILTERS)
 
466
            G_fatal_error(_("Too many filters (max: %d)"), MAX_FILTERS);
 
467
 
 
468
        filter = &filters[num_filters++];
 
469
        type = find_method(filter_arg);
 
470
        filter->func = type->func;
 
471
        filter->x_radius = fabs(atof(x_radius_arg));
 
472
        filter->y_radius = fabs(atof(y_radius_arg));
 
473
        if (type->radius) {
 
474
            double rx = type->radius * filter->x_radius;
 
475
            double ry = type->radius * filter->y_radius;
 
476
            if (rx < f_x_radius)
 
477
                f_x_radius = rx;
 
478
            if (ry < f_y_radius)
 
479
                f_y_radius = ry;
 
480
        }
 
481
    }
 
482
 
 
483
    if (f_x_radius > 1e99 || f_y_radius > 1e99)
 
484
        G_fatal_error(_("At least one filter must be finite"));
 
485
 
 
486
    G_get_set_window(&dst_w);
 
487
 
 
488
    /* set window to old map */
 
489
    Rast_get_cellhd(parm.rastin->answer, "", &src_w);
 
490
 
 
491
    /* enlarge source window */
 
492
    {
 
493
        double y0 = Rast_row_to_northing(0.5, &dst_w);
 
494
        double y1 = Rast_row_to_northing(dst_w.rows - 0.5, &dst_w);
 
495
        double x0 = Rast_col_to_easting(0.5, &dst_w);
 
496
        double x1 = Rast_col_to_easting(dst_w.cols - 0.5, &dst_w);
 
497
        int r0 = (int)floor(Rast_northing_to_row(y0 + f_y_radius, &src_w) - 0.1);
 
498
        int r1 = (int)ceil(Rast_northing_to_row(y1 - f_y_radius, &src_w) + 0.1);
 
499
        int c0 = (int)floor(Rast_easting_to_col(x0 - f_x_radius, &src_w) - 0.1);
 
500
        int c1 = (int)ceil(Rast_easting_to_col(x1 + f_x_radius, &src_w) + 0.1);
 
501
 
 
502
        src_w.south -= src_w.ns_res * (r1 - src_w.rows);
 
503
        src_w.north += src_w.ns_res * (-r0);
 
504
        src_w.west -= src_w.ew_res * (-c0);
 
505
        src_w.east += src_w.ew_res * (c1 - src_w.cols);
 
506
        src_w.rows = r1 - r0;
 
507
        src_w.cols = c1 - c0;
 
508
    }
 
509
 
 
510
    row_scale = 2 + 2 * ceil(f_y_radius / src_w.ns_res);
 
511
    col_scale = 2 + 2 * ceil(f_x_radius / src_w.ew_res);
 
512
 
 
513
    /* allocate buffers for intermediate rows */
 
514
    bufs = G_malloc(row_scale * sizeof(DCELL *));
 
515
    for (i = 0; i < row_scale; i++)
 
516
        bufs[i] = Rast_allocate_d_buf();
 
517
 
 
518
    Rast_set_input_window(&src_w);
 
519
    Rast_set_output_window(&dst_w);
 
520
 
 
521
    inbuf = Rast_allocate_d_input_buf();
 
522
    outbuf = Rast_allocate_d_output_buf();
 
523
 
 
524
    infile = Rast_open_old(parm.rastin->answer, "");
 
525
    outfile = Rast_open_new(parm.rastout->answer, DCELL_TYPE);
 
526
 
 
527
    filter();
 
528
 
 
529
    Rast_close(infile);
 
530
    Rast_close(outfile);
 
531
 
 
532
    /* record map metadata/history info */
 
533
    sprintf(title, "Filter resample by %s", parm.method->answer);
 
534
    Rast_put_cell_title(parm.rastout->answer, title);
 
535
 
 
536
    {
 
537
        struct History history;
 
538
        char buf_nsres[100], buf_ewres[100];
 
539
 
 
540
        Rast_short_history(parm.rastout->answer, "raster", &history);
 
541
        Rast_set_history(&history, HIST_DATSRC_1, parm.rastin->answer);
 
542
        G_format_resolution(src_w.ns_res, buf_nsres, src_w.proj);
 
543
        G_format_resolution(src_w.ew_res, buf_ewres, src_w.proj);
 
544
        Rast_format_history(&history, HIST_DATSRC_2,
 
545
                            "Source map NS res: %s   EW res: %s",
 
546
                            buf_nsres, buf_ewres);
 
547
        Rast_command_history(&history);
 
548
        Rast_write_history(parm.rastout->answer, &history);
 
549
    }
 
550
 
 
551
    /* copy color table from source map */
 
552
    {
 
553
        struct Colors colors;
 
554
 
 
555
        if (Rast_read_colors(parm.rastin->answer, "", &colors) < 0)
 
556
            G_fatal_error(_("Unable to read color table for %s"),
 
557
                          parm.rastin->answer);
 
558
        Rast_mark_colors_as_fp(&colors);
 
559
        Rast_write_colors(parm.rastout->answer, G_mapset(), &colors);
 
560
    }
 
561
 
 
562
    return EXIT_SUCCESS;
 
563
}