1
2 /****************************************************************************
3 *
4 * MODULE: r3.neighbors
5 *
6 * AUTHOR(S): Original author
7 * Soeren Gebbert soerengebbert <at> googlemail <dot> co
8 * with code from r.series and r.neighbors for parameter menu handling
9 *
10 * PURPOSE: Makes each voxel value a function of the values assigned to the voxels
11 * around it, and stores new voxel values in an output 3D raster map
12 *
13 * COPYRIGHT: (C) 2013 by the GRASS Development Team
14 *
15 * This program is free software under the GNU General Public
16 * License (>=v2). Read the file COPYING that comes with GRASS
17 * for details.
18 *
19 *****************************************************************************/
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include <string.h>
23 #include <grass/raster3d.h>
24 #include <grass/stats.h>
25 #include <grass/gis.h>
26 #include <grass/raster.h>
27 #include <grass/glocale.h>
28
29 int nx, ny, nz; /* Number of cells in x, y and z direction */
30 int x_dist, y_dist, z_dist; /* Distance of cells from the center
31 to the edge of the moving window */
32 int x_size, y_size, z_size; /* The size of the moving window in x,
33 y and z direction */
34 int size; /* The maximum size of the value buffer */
35
36 struct menu
37 {
38 stat_func *method; /* routine to compute new value */
39 char *name; /* method name */
40 char *text; /* menu display - full description */
41 } menu[] = {
42 {
43 c_ave, "average", "average value"}, {
44 c_median, "median", "median value"}, {
45 c_mode, "mode", "most frequently occurring value"}, {
46 c_min, "minimum", "lowest value"}, {
47 c_max, "maximum", "highest value"}, {
48 c_range, "range", "range value"}, {
49 c_stddev, "stddev", "standard deviation"}, {
50 c_sum, "sum", "sum of values"}, {
51 c_count, "count", "count of non-NULL values"}, {
52 c_var, "variance", "statistical variance"}, {
53 c_divr, "diversity", "number of different values"}, {
54 c_intr, "interspersion", "number of values different than center value"},
55 {
56 c_quart1, "quart1", "first quartile"}, {
57 c_quart3, "quart3", "third quartile"}, {
58 c_perc90, "perc90", "ninetieth percentile"}, {
59 c_quant, "quantile", "arbitrary quantile"}, {
60 NULL, NULL, NULL}
61 };
62
63 /* ************************************************************************* */
64
build_method_list(void)65 static char *build_method_list(void)
66 {
67 char *buf = G_malloc(1024);
68 char *p = buf;
69 int i;
70
71 for (i = 0; menu[i].name; i++) {
72 char *q;
73
74 if (i)
75 *p++ = ',';
76 for (q = menu[i].name; *q; p++, q++)
77 *p = *q;
78 }
79 *p = '\0';
80
81 return buf;
82 }
83
84 /* ************************************************************************* */
85
find_method(const char * method_name)86 static int find_method(const char *method_name)
87 {
88 int i;
89
90 for (i = 0; menu[i].name; i++)
91 if (strcmp(menu[i].name, method_name) == 0)
92 return i;
93
94 G_fatal_error(_("Unknown method <%s>"), method_name);
95
96 return -1;
97 }
98
99 /* ************************************************************************* */
100
101 typedef struct
102 {
103 struct Option *input, *output, *window, *method, *quantile;
104 } paramType;
105
106 paramType param;
107
set_params()108 static void set_params()
109 {
110 param.input = G_define_standard_option(G_OPT_R3_INPUT);
111
112 param.output = G_define_standard_option(G_OPT_R3_OUTPUT);
113
114 param.method = G_define_option();
115 param.method->key = "method";
116 param.method->type = TYPE_STRING;
117 param.method->required = YES;
118 param.method->options = build_method_list();
119 param.method->description = _("Aggregate operation");
120 param.method->multiple = NO;
121
122 param.quantile = G_define_option();
123 param.quantile->key = "quantile";
124 param.quantile->type = TYPE_DOUBLE;
125 param.quantile->required = NO;
126 param.quantile->description =
127 _("Quantile to calculate for method=quantile");
128 param.quantile->options = "0.0-1.0";
129 param.quantile->multiple = NO;
130
131 param.window = G_define_option();
132 param.window->key = "window";
133 param.window->type = TYPE_INTEGER;
134 param.window->required = YES;
135 param.window->key_desc = "x,y,z";
136 param.window->description =
137 _("The size of the window in x, y and z direction,"
138 " values must be odd integer numbers, eg: 3,3,3");
139 }
140
141 /* ************************************************************************* */
142
gather_values(RASTER3D_Map * map,DCELL * buff,int x,int y,int z)143 static int gather_values(RASTER3D_Map * map, DCELL * buff, int x,
144 int y, int z)
145 {
146 int i, j, k, l;
147 DCELL value;
148
149 int start_z = z - z_dist;
150 int start_y = y - y_dist;
151 int start_x = x - x_dist;
152 int end_z = start_z + z_size;
153 int end_y = start_y + y_size;
154 int end_x = start_x + x_size;
155
156 if (start_z < 0)
157 start_z = 0;
158
159 if (start_y < 0)
160 start_y = 0;
161
162 if (start_x < 0)
163 start_x = 0;
164
165 if (end_z > nz)
166 end_z = nz;
167
168 if (end_y > ny)
169 end_y = ny;
170
171 if (end_x > nx)
172 end_x = nx;
173
174 l = 0;
175
176 for (i = start_z; i < end_z; i++) {
177 for (j = start_y; j < end_y; j++) {
178 for (k = start_x; k < end_x; k++) {
179 value = (DCELL) Rast3d_get_double(map, k, j, i);
180
181 if (Rast_is_d_null_value(&value))
182 continue;
183
184 buff[l] = value;
185 l++;
186 }
187 }
188 }
189 return l;
190 }
191
192 /* ************************************************************************* */
193
main(int argc,char ** argv)194 int main(int argc, char **argv)
195 {
196 RASTER3D_Map *input;
197 RASTER3D_Map *output;
198 RASTER3D_Region region;
199 struct GModule *module;
200 stat_func *method_fn;
201 double quantile;
202 int x, y, z;
203
204 /* Initialize GRASS */
205 G_gisinit(argv[0]);
206
207 module = G_define_module();
208 G_add_keyword(_("raster3d"));
209 G_add_keyword(_("algebra"));
210 G_add_keyword(_("voxel"));
211 G_add_keyword(_("statistics"));
212 G_add_keyword(_("aggregation"));
213 G_add_keyword(_("neighbor"));
214 G_add_keyword(_("focal statistics"));
215 G_add_keyword(_("filter"));
216
217 module->description =
218 _("Makes each voxel value a "
219 "function of the values assigned to the voxels "
220 "around it, and stores new voxel values in an output 3D raster map");
221
222 /* Get parameters from user */
223 set_params();
224
225 if (G_parser(argc, argv))
226 exit(EXIT_FAILURE);
227
228 if (NULL == G_find_raster3d(param.input->answer, ""))
229 Rast3d_fatal_error(_("3D raster map <%s> not found"),
230 param.input->answer);
231
232 Rast3d_init_defaults();
233 Rast3d_get_window(®ion);
234
235 nx = region.cols;
236 ny = region.rows;
237 nz = region.depths;
238
239 /* Size fo the moving window */
240 x_size = atoi(param.window->answers[0]);
241 y_size = atoi(param.window->answers[1]);
242 z_size = atoi(param.window->answers[2]);
243
244 /* Distances in all directions */
245 x_dist = x_size / 2;
246 y_dist = y_size / 2;
247 z_dist = z_size / 2;
248
249 /* Maximum size of the buffer */
250 size = x_size * y_size * z_size;
251
252 /* Set the computation method */
253 method_fn = menu[find_method(param.method->answer)].method;
254
255 if (param.quantile->answer)
256 quantile = atof(param.quantile->answer);
257 else
258 quantile = 0.0;
259
260 input = Rast3d_open_cell_old(param.input->answer,
261 G_find_raster3d(param.input->answer, ""),
262 ®ion, RASTER3D_TILE_SAME_AS_FILE,
263 RASTER3D_USE_CACHE_DEFAULT);
264
265 if (input == NULL)
266 Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
267 param.input->answer);
268
269 output =
270 Rast3d_open_new_opt_tile_size(param.output->answer,
271 RASTER3D_USE_CACHE_X, ®ion,
272 DCELL_TYPE, 32);
273
274 if (output == NULL)
275 Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
276 param.output->answer);
277
278 Rast3d_min_unlocked(output, RASTER3D_USE_CACHE_X);
279 Rast3d_autolock_on(output);
280 Rast3d_unlock_all(output);
281
282 DCELL *buff = NULL, value;
283
284 buff = (DCELL *) calloc(size, sizeof(DCELL));
285
286 if (buff == NULL)
287 Rast3d_fatal_error(_("Unable to allocate buffer"));
288
289 for (z = 0; z < nz; z++) {
290 G_percent(z, nz, 1);
291 for (y = 0; y < ny; y++) {
292 for (x = 0; x < nx; x++) {
293 /* Gather values in moving window */
294 int num = gather_values(input, buff, x, y, z);
295
296 /* Compute the resulting value */
297 if (num > 0)
298 (*method_fn) (&value, buff, num, &quantile);
299 else
300 Rast_set_d_null_value(&value, 1);
301 /* Write the value */
302 Rast3d_put_double(output, x, y, z, value);
303 }
304 }
305 }
306 G_percent(z, nz, 1);
307
308 free(buff);
309
310 if (!Rast3d_flush_all_tiles(output))
311 G_fatal_error(_("Error flushing tiles"));
312
313 Rast3d_autolock_off(output);
314 Rast3d_unlock_all(output);
315
316 Rast3d_close(input);
317 Rast3d_close(output);
318
319 return 0;
320 }
321