← Back to branch summary

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

« back to all changes in this revision

Viewing changes to raster/r.colors/stats.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.colors/stats.c
19
19
 
20
20
#include <math.h>
21
21
#include <stdlib.h>
22
 
#include <grass/gis.h>
23
 
#include <grass/glocale.h>
24
22
#include "local_proto.h"
25
23
 
26
 
int get_stats(const char *name, const char *mapset, struct Cell_stats *statf)
27
 
{
 
24
int get_stats(struct maps_info *input_maps, struct Cell_stats *statf) {
28
25
    CELL *cell;
29
26
    int row, nrows, ncols;
30
27
    int fd;
31
 
 
32
 
    if ((fd = G_open_cell_old(name, mapset)) < 0)
33
 
        G_fatal_error(_("Unable to open raster map <%s>"),
34
 
                      G_fully_qualified_name(name, mapset));
35
 
 
36
 
    cell = G_allocate_cell_buf();
37
 
    nrows = G_window_rows();
38
 
    ncols = G_window_cols();
39
 
 
40
 
    G_init_cell_stats(statf);
41
 
    G_verbose_message(_("Reading raster map <%s>..."),
42
 
                      G_fully_qualified_name(name, mapset));
43
 
    for (row = 0; row < nrows; row++) {
44
 
        G_percent(row, nrows, 2);
45
 
        if (G_get_c_raster_row(fd, cell, row) < 0)
46
 
            G_fatal_error(_("Unable to read raster map <%s> row %d"),
47
 
                          G_fully_qualified_name(name, mapset), row);
48
 
        G_update_cell_stats(cell, ncols, statf);
 
28
    int i;
 
29
 
 
30
    Rast_init_cell_stats(statf);
 
31
 
 
32
    for(i = 0; i < input_maps->num; i++) {
 
33
                fd = Rast_open_old(input_maps->names[i], input_maps->mapsets[i]);
 
34
 
 
35
            cell = Rast_allocate_c_buf();
 
36
            nrows = Rast_window_rows();
 
37
            ncols = Rast_window_cols();
 
38
 
 
39
                G_verbose_message(_("(%i/%i) Reading raster map <%s>..."),
 
40
                                i + 1, input_maps->num,
 
41
                                G_fully_qualified_name(input_maps->names[i], input_maps->mapsets[i]));
 
42
 
 
43
                for (row = 0; row < nrows; row++) {
 
44
                        G_percent(row, nrows, 2);
 
45
                        Rast_get_c_row(fd, cell, row);
 
46
                        Rast_update_cell_stats(cell, ncols, statf);
 
47
                }
 
48
                G_percent(row, nrows, 2);
 
49
                Rast_close(fd);
 
50
            G_free(cell);
49
51
    }
50
 
    G_percent(row, nrows, 2);
51
 
    G_close_cell(fd);
52
 
    G_free(cell);
53
52
 
54
53
    return 1;
55
54
}
56
55
 
57
 
void get_fp_stats(const char *name, const char *mapset,
58
 
                  struct FP_stats *statf,
59
 
                  DCELL min, DCELL max, int geometric, int geom_abs)
60
 
{
61
 
    DCELL *dcell;
62
 
    int row, col, nrows, ncols;
 
56
void get_fp_stats(struct maps_info *input_maps,
 
57
        struct FP_stats *statf,
 
58
        DCELL min, DCELL max, int geometric, int geom_abs, int type) {
 
59
    DCELL *dcell = NULL;
 
60
    int row, col, depth, nrows, ncols, ndepths = 1;
63
61
    int fd;
64
 
 
65
 
    if ((fd = G_open_cell_old(name, mapset)) < 0)
66
 
        G_fatal_error("Unable to open raster map <%s>",
67
 
                      G_fully_qualified_name(name, mapset));
68
 
 
69
 
    dcell = G_allocate_d_raster_buf();
70
 
    nrows = G_window_rows();
71
 
    ncols = G_window_cols();
72
 
 
73
 
    statf->geometric = geometric;
74
 
    statf->geom_abs = geom_abs;
75
 
    statf->flip = 0;
76
 
 
77
 
    if (statf->geometric) {
78
 
        if (min * max < 0)
79
 
            G_fatal_error(_("Unable to use logarithmic scaling if range includes zero"));
80
 
 
81
 
        if (min < 0) {
82
 
            statf->flip = 1;
83
 
            min = -min;
84
 
            max = -max;
85
 
        }
86
 
 
87
 
        min = log(min);
88
 
        max = log(max);
89
 
    }
90
 
 
91
 
    if (statf->geom_abs) {
92
 
        double a = log(fabs(min) + 1);
93
 
        double b = log(fabs(max) + 1);
94
 
        int has_zero = min * max < 0;
95
 
        min = a < b ? a : b;
96
 
        max = a > b ? a : b;
97
 
        if (has_zero)
98
 
            min = 0;
99
 
    }
100
 
 
101
 
    statf->count = 1000;
102
 
    statf->min = min;
103
 
    statf->max = max;
104
 
    statf->stats = G_calloc(statf->count + 1, sizeof(unsigned long));
105
 
    statf->total = 0;
106
 
 
107
 
    G_verbose_message(_("Reading raster map <%s>..."),
108
 
                      G_fully_qualified_name(name, mapset));
109
 
    for (row = 0; row < nrows; row++) {
110
 
        G_percent(row, nrows, 2);
111
 
 
112
 
        if (G_get_d_raster_row(fd, dcell, row) < 0)
113
 
            G_fatal_error(_("Unable to read raster map <%s> row %d"),
114
 
                          G_fully_qualified_name(name, mapset), row);
115
 
        
116
 
        for (col = 0; col < ncols; col++) {
117
 
            DCELL x;
118
 
            int i;
119
 
 
120
 
            if (G_is_d_null_value(&dcell[col]))
121
 
                continue;
122
 
 
123
 
            x = dcell[col];
124
 
            if (statf->flip)
125
 
                x = -x;
126
 
            if (statf->geometric)
127
 
                x = log(x);
128
 
            if (statf->geom_abs)
129
 
                x = log(fabs(x) + 1);
130
 
 
131
 
            i = (int) floor(statf->count * (x - statf->min) / (statf->max - statf->min));
132
 
            statf->stats[i]++;
133
 
            statf->total++;
134
 
        }
135
 
    }
136
 
 
137
 
    G_percent(row, nrows, 2);
138
 
    G_close_cell(fd);
139
 
    G_free(dcell);
 
62
    int i;
 
63
    char *name;
 
64
    char *mapset;
 
65
    RASTER3D_Map *map3d = NULL;
 
66
 
 
67
        statf->geometric = geometric;
 
68
        statf->geom_abs = geom_abs;
 
69
        statf->flip = 0;
 
70
 
 
71
        if (statf->geometric) {
 
72
                if (min * max < 0)
 
73
                        G_fatal_error(_("Unable to use logarithmic scaling if range includes zero"));
 
74
 
 
75
                if (min < 0) {
 
76
                        statf->flip = 1;
 
77
                        min = -min;
 
78
                        max = -max;
 
79
                }
 
80
                min = log(min);
 
81
                max = log(max);
 
82
        }
 
83
 
 
84
        if (statf->geom_abs) {
 
85
                double a = log(fabs(min) + 1);
 
86
                double b = log(fabs(max) + 1);
 
87
                int has_zero = min * max < 0;
 
88
                min = a < b ? a : b;
 
89
                max = a > b ? a : b;
 
90
                if (has_zero)
 
91
                        min = 0;
 
92
        }
 
93
 
 
94
        statf->count = 1000;
 
95
        statf->min = min;
 
96
        statf->max = max;
 
97
        statf->stats = G_calloc(statf->count + 1, sizeof (unsigned long));
 
98
        statf->total = 0;
 
99
 
 
100
        /* Loop over all input maps */
 
101
        for(i = 0; i < input_maps->num; i++) {
 
102
                name = input_maps->names[i];
 
103
                mapset = input_maps->mapsets[i];
 
104
 
 
105
                if (type == RASTER_TYPE) {
 
106
                        fd = Rast_open_old(name, mapset);
 
107
                        dcell = Rast_allocate_d_buf();
 
108
                        nrows = Rast_window_rows();
 
109
                        ncols = Rast_window_cols();
 
110
                } else {
 
111
                        /* Initiate the default settings */
 
112
                        Rast3d_init_defaults();
 
113
 
 
114
                        map3d = Rast3d_open_cell_old(name, mapset, RASTER3D_DEFAULT_WINDOW,
 
115
                                        RASTER3D_TILE_SAME_AS_FILE, RASTER3D_USE_CACHE_DEFAULT);
 
116
 
 
117
                        if (map3d == NULL)
 
118
                                Rast3d_fatal_error(_("Error opening 3d raster map"));
 
119
 
 
120
                        nrows = map3d->window.rows;
 
121
                        ncols = map3d->window.cols;
 
122
                        ndepths = map3d->window.depths;
 
123
                }
 
124
 
 
125
                G_verbose_message(_("(%i/%i) Reading map <%s>..."), i, input_maps->num,
 
126
                                G_fully_qualified_name(name, mapset));
 
127
 
 
128
                for (depth = 0; depth < ndepths; depth++) {
 
129
                        for (row = 0; row < nrows; row++) {
 
130
                                G_percent(row, nrows, 2);
 
131
 
 
132
                                if (type == RASTER_TYPE)
 
133
                                        Rast_get_d_row(fd, dcell, row);
 
134
 
 
135
                                for (col = 0; col < ncols; col++) {
 
136
                                        DCELL x;
 
137
                                        int j;
 
138
 
 
139
                                        if (type == RASTER_TYPE)
 
140
                                                x = dcell[col];
 
141
                                        else
 
142
                                                x = Rast3d_get_double(map3d, col, row, depth);
 
143
 
 
144
                                        if (Rast_is_d_null_value(&x))
 
145
                                                continue;
 
146
 
 
147
                                        if (statf->flip)
 
148
                                                x = -x;
 
149
                                        if (statf->geometric)
 
150
                                                x = log(x);
 
151
                                        if (statf->geom_abs)
 
152
                                                x = log(fabs(x) + 1);
 
153
 
 
154
                                        j = (int) floor(statf->count * (x - statf->min) / (statf->max - statf->min));
 
155
                                        statf->stats[j]++;
 
156
                                        statf->total++;
 
157
                                }
 
158
                        }
 
159
                }
 
160
 
 
161
                G_percent(row, nrows, 2);
 
162
 
 
163
                if (type == RASTER_TYPE) {
 
164
                        Rast_close(fd);
 
165
                if(dcell)
 
166
                                G_free(dcell);
 
167
                } else {
 
168
                        Rast3d_close(map3d);
 
169
                }
 
170
        }
140
171
}