[GRASS-SVN] r61547 - in grass-addons/grass7/raster3d: . r3.gradient
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Aug 5 19:46:37 PDT 2014
Author: annakrat
Date: 2014-08-05 19:46:37 -0700 (Tue, 05 Aug 2014)
New Revision: 61547
Added:
grass-addons/grass7/raster3d/r3.gradient/
grass-addons/grass7/raster3d/r3.gradient/Makefile
grass-addons/grass7/raster3d/r3.gradient/main.c
grass-addons/grass7/raster3d/r3.gradient/r3.gradient.html
grass-addons/grass7/raster3d/r3.gradient/r3gradient_structs.h
Modified:
grass-addons/grass7/raster3d/Makefile
Log:
r3.gradient added to addons
Modified: grass-addons/grass7/raster3d/Makefile
===================================================================
--- grass-addons/grass7/raster3d/Makefile 2014-08-05 22:54:57 UTC (rev 61546)
+++ grass-addons/grass7/raster3d/Makefile 2014-08-06 02:46:37 UTC (rev 61547)
@@ -1,7 +1,8 @@
MODULE_TOPDIR = ..
SUBDIRS = \
- r3.flow
+ r3.flow \
+ r3.gradient
include $(MODULE_TOPDIR)/include/Make/Dir.make
Added: grass-addons/grass7/raster3d/r3.gradient/Makefile
===================================================================
--- grass-addons/grass7/raster3d/r3.gradient/Makefile (rev 0)
+++ grass-addons/grass7/raster3d/r3.gradient/Makefile 2014-08-06 02:46:37 UTC (rev 61547)
@@ -0,0 +1,11 @@
+MODULE_TOPDIR = ../..
+
+PGM=r3.gradient
+
+LIBES = $(RASTER3DLIB) $(RASTERLIB) $(GISLIB) $(OMPLIB)
+DEPENDENCIES = $(RASTER3DDEP) $(RASTERDEP) $(GISDEP)
+EXTRA_CFLAGS = $(OMPCFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/grass7/raster3d/r3.gradient/main.c
===================================================================
--- grass-addons/grass7/raster3d/r3.gradient/main.c (rev 0)
+++ grass-addons/grass7/raster3d/r3.gradient/main.c 2014-08-06 02:46:37 UTC (rev 61547)
@@ -0,0 +1,299 @@
+/****************************************************************************
+ *
+ * MODULE: r3.gradient
+ * AUTHOR(S): Anna Petrasova kratochanna <at> gmail <dot> com
+ * PURPOSE: Computes gradient of a 3D raster map
+ * COPYRIGHT: (C) 2014 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
+
+#include <grass/raster3d.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "r3gradient_structs.h"
+
+int main(int argc, char *argv[])
+{
+ struct Option *input_opt, *output_opt, *block_opt, *process_opt;
+ struct GModule *module;
+ RASTER3D_Region region;
+ RASTER3D_Map *input;
+ RASTER3D_Map *output[3];
+ struct Gradient_block *blocks;
+ int block_x, block_y, block_z;
+ int index_x, index_y, index_z;
+ int n_x, n_y, n_z;
+ int start_x, start_y, start_z;
+ int i, max_i, k, j, N;
+ double step[3];
+ int *bl_indices;
+ int *bl_overlap;
+ int r, c, d;
+ DCELL value;
+
+ module = G_define_module();
+ G_add_keyword(_("raster3d"));
+ G_add_keyword(_("voxel"));
+ G_add_keyword(_("gradient"));
+ module->description =
+ _("Computes gradient of a 3D raster map "
+ "and outputs gradient components as three 3D raster maps.");
+
+ input_opt = G_define_standard_option(G_OPT_R3_INPUT);
+
+ /* TODO: define G_OPT_R3_OUTPUTS or use separate options for each map? */
+ output_opt = G_define_standard_option(G_OPT_R3_OUTPUT);
+ output_opt->multiple = YES;
+ output_opt->description = _("Name for output 3D raster map(s)");
+
+ block_opt = G_define_option();
+ block_opt->key = "block_size";
+ block_opt->multiple = TRUE;
+ block_opt->answer = "30,30,20"; /* based on testing */
+ block_opt->key_desc = "size_x,size_y,size_z";
+ block_opt->description = _("Size of blocks");
+
+ process_opt = G_define_option();
+ process_opt->key = "nprocs";
+ process_opt->type = TYPE_INTEGER;
+ process_opt->required = NO;
+ process_opt->description = _("Number of parallel processes");
+ process_opt->options = "1-100";
+ process_opt->answer = "1";
+
+ G_gisinit(argv[0]);
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ N = atoi(process_opt->answer);
+#if defined(_OPENMP)
+ omp_set_num_threads(N);
+#endif
+
+ Rast3d_init_defaults();
+ Rast3d_get_window(®ion);
+
+ block_x = atoi(block_opt->answers[0]);
+ block_y = atoi(block_opt->answers[1]);
+ block_z = atoi(block_opt->answers[2]);
+
+ if (block_x < 3 || block_y < 3 || block_z < 3)
+ G_warning("block size too small, set to 3");
+
+ block_x = block_x < 3 ? 3 : block_x;
+ block_y = block_y < 3 ? 3 : block_y;
+ block_z = block_z < 3 ? 3 : block_z;
+ block_x = block_x > region.cols ? region.cols : block_x;
+ block_y = block_y > region.rows ? region.rows : block_y;
+ block_z = block_y > region.depths ? region.depths : block_z;
+
+ step[0] = region.ew_res;
+ step[1] = region.ns_res;
+ step[2] = region.tb_res;
+
+ input = Rast3d_open_cell_old(input_opt->answer,
+ G_find_raster3d(input_opt->answer, ""),
+ ®ion, RASTER3D_TILE_SAME_AS_FILE,
+ RASTER3D_USE_CACHE_DEFAULT);
+ if (!input)
+ Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
+ input_opt->answer);
+
+ for (i = 0; i < 3; i++) {
+ output[i] =
+ Rast3d_open_new_opt_tile_size(output_opt->answers[i],
+ RASTER3D_USE_CACHE_DEFAULT,
+ ®ion, DCELL_TYPE, 32);
+ if (!output[i]) {
+ Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
+ output_opt->answers[i]);
+ }
+
+ }
+
+ blocks = G_calloc(N, sizeof(struct Gradient_block));
+ if (!blocks)
+ G_fatal_error(_("Failed to allocate memory for blocks"));
+ for (i = 0; i < N; i++) {
+ blocks[i].input.array = G_malloc(((block_x + 2) * (block_y + 2)
+ * (block_z + 2)) * sizeof(DCELL));
+ blocks[i].dx.array = G_malloc(((block_x + 2) * (block_y + 2)
+ * (block_z + 2)) * sizeof(DCELL));
+ blocks[i].dy.array = G_malloc(((block_x + 2) * (block_y + 2)
+ * (block_z + 2)) * sizeof(DCELL));
+ blocks[i].dz.array = G_malloc(((block_x + 2) * (block_y + 2)
+ * (block_z + 2)) * sizeof(DCELL));
+ }
+
+ bl_indices = G_calloc(N * 3, sizeof(int));
+ bl_overlap = G_calloc(N * 6, sizeof(int));
+
+ max_i = (int)ceil(region.cols / (float)block_x) *
+ (int)ceil(region.rows / (float)block_y) *
+ (int)ceil(region.depths / (float)block_z);
+ i = j = 0;
+ index_z = 0;
+
+ /* loop through the blocks */
+ while (index_z < region.depths) {
+ index_y = 0;
+ while (index_y < region.rows) {
+ index_x = 0;
+ while (index_x < region.cols) {
+ G_percent(i, max_i, 1);
+ /* generally overlap is 1 on both sides */
+ bl_overlap[j * 6 + 0] = bl_overlap[j * 6 + 2] =
+ bl_overlap[j * 6 + 4] = 1;
+ bl_overlap[j * 6 + 1] = bl_overlap[j * 6 + 3] =
+ bl_overlap[j * 6 + 5] = 1;
+
+ /* compute the starting index of the block and its size */
+ start_x = fmax(index_x - 1, 0);
+ n_x = fmin(index_x + block_x, region.cols - 1) - start_x + 1;
+ start_y = fmax(index_y - 1, 0);
+ n_y = fmin(index_y + block_y, region.rows - 1) - start_y + 1;
+ start_z = fmax(index_z - 1, 0);
+ n_z =
+ fmin(index_z + block_z, region.depths - 1) - start_z + 1;
+
+ /* adjust offset on edges */
+ /* start offset */
+ if (index_x == 0)
+ bl_overlap[j * 6 + 0] = 0;
+ if (index_y == 0)
+ bl_overlap[j * 6 + 2] = 0;
+ if (index_z == 0)
+ bl_overlap[j * 6 + 4] = 0;
+ /* end offset */
+ if (index_x + block_x >= region.cols)
+ bl_overlap[j * 6 + 1] = 0;
+ if (index_y + block_y >= region.rows)
+ bl_overlap[j * 6 + 3] = 0;
+ if (index_z + block_z >= region.depths)
+ bl_overlap[j * 6 + 5] = 0;
+ /* adjust offset when the end block would be too small */
+ if (n_x <= 2) {
+ start_x -= 1;
+ n_x += 1;
+ bl_overlap[j * 6 + 0] = 2;
+ }
+ if (n_y <= 2) {
+ start_y -= 1;
+ n_y += 1;
+ bl_overlap[j * 6 + 2] = 2;
+ }
+ if (n_z <= 2) {
+ start_z -= 1;
+ n_z += 1;
+ bl_overlap[j * 6 + 4] = 2;
+ }
+ /* store indices for later writing */
+ bl_indices[j * 3 + 0] = index_x;
+ bl_indices[j * 3 + 1] = index_y;
+ bl_indices[j * 3 + 2] = index_z;
+
+ blocks[j].input.sx = n_x;
+ blocks[j].input.sy = n_y;
+ blocks[j].input.sz = n_z;
+ blocks[j].dx.sx = blocks[j].dy.sx = blocks[j].dz.sx = n_x;
+ blocks[j].dx.sy = blocks[j].dy.sy = blocks[j].dz.sy = n_y;
+ blocks[j].dx.sz = blocks[j].dy.sz = blocks[j].dz.sz = n_z;
+
+ /* read */
+ Rast3d_get_block(input, start_x, start_y, start_z,
+ n_x, n_y, n_z, blocks[j].input.array,
+ DCELL_TYPE);
+ if ((j + 1) == N || i == max_i - 1) {
+
+ /* compute gradient */
+ #pragma omp parallel for schedule (static) private (k)
+ for (k = 0; k <= j; k++) {
+ Rast3d_gradient_double(&(blocks[k].input), step,
+ &(blocks[k].dx), &(blocks[k].dy),
+ &(blocks[k].dz));
+ }
+
+ /* write */
+ for (k = 0; k <= j; k++) {
+ for (c = 0;c < blocks[k].input.sx - bl_overlap[k * 6 + 0] -
+ bl_overlap[k * 6 + 1]; c++) {
+ for (r = 0; r < blocks[k].input.sy - bl_overlap[k * 6 + 2] -
+ bl_overlap[k * 6 + 3]; r++) {
+ for (d = 0; d < blocks[k].input.sz - bl_overlap[k * 6 + 4] -
+ bl_overlap[k * 6 + 5]; d++) {
+ value = RASTER3D_ARRAY_ACCESS(&(blocks[k].dx),
+ c + bl_overlap[k * 6 + 0],
+ r + bl_overlap[k * 6 + 2],
+ d + bl_overlap[k * 6 + 4]);
+ Rast3d_put_value(output[0],
+ c + bl_indices[k * 3 + 0],
+ r + bl_indices[k * 3 + 1],
+ d + bl_indices[k * 3 + 2],
+ &value, DCELL_TYPE);
+
+ value = RASTER3D_ARRAY_ACCESS(&(blocks[k].dy),
+ c + bl_overlap[k * 6 + 0],
+ r + bl_overlap[k * 6 + 2],
+ d + bl_overlap[k * 6 + 4]);
+ Rast3d_put_value(output[1],
+ c + bl_indices[k * 3 + 0],
+ r + bl_indices[k * 3 + 1],
+ d + bl_indices[k * 3 + 2],
+ &value, DCELL_TYPE);
+
+ value = RASTER3D_ARRAY_ACCESS(&(blocks[k].dz),
+ c + bl_overlap[k * 6 + 0],
+ r + bl_overlap[k * 6 + 2],
+ d + bl_overlap[k * 6 + 4]);
+ Rast3d_put_value(output[2],
+ c + bl_indices[k * 3 + 0],
+ r + bl_indices[k * 3 + 1],
+ d + bl_indices[k * 3 + 2],
+ &value, DCELL_TYPE);
+ }
+ }
+ }
+ }
+ j = -1;
+ }
+ i++;
+ j++;
+ index_x += block_x;
+ }
+ index_y += block_y;
+ }
+ index_z += block_z;
+ }
+ G_percent(1, 1, 1);
+ for (i = 0; i < N; i++) {
+ G_free(blocks[i].input.array);
+ G_free(blocks[i].dx.array);
+ G_free(blocks[i].dy.array);
+ G_free(blocks[i].dz.array);
+ }
+ G_free(blocks);
+ G_free(bl_indices);
+ G_free(bl_overlap);
+
+ G_message(_("Writing gradient 3D raster maps..."));
+ G_percent(0, 3, 1);
+ Rast3d_close(output[0]);
+ G_percent(1, 3, 1);
+ Rast3d_close(output[1]);
+ G_percent(2, 3, 1);
+ Rast3d_close(output[2]);
+ G_percent(1, 1, 1);
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass-addons/grass7/raster3d/r3.gradient/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster3d/r3.gradient/r3.gradient.html
===================================================================
--- grass-addons/grass7/raster3d/r3.gradient/r3.gradient.html (rev 0)
+++ grass-addons/grass7/raster3d/r3.gradient/r3.gradient.html 2014-08-06 02:46:37 UTC (rev 61547)
@@ -0,0 +1,23 @@
+<h2>DESCRIPTION</h2>
+Module <em>r3.gradient</em> computes gradient from a 3D raster map.
+Results are three 3D raster maps describing the x, y, z components of the computed gradient field.
+
+<h2>EXAMPLES</h2>
+<div class="code"><pre>
+# create a 3D raster
+g.region s=0 n=100 w=0 e=100 b=0 t=100
+r3.mapcalc "test_gradient = sqrt(row()*row() +col()*col()+ depth()*depth())"
+
+# compute gradient
+r3.gradient input=test_gradient output=grad_x,grad_y,grad_z
+</pre></div>
+
+<h2>AUTHORS</h2>
+Anna Petrasova, <a href="http://geospatial.ncsu.edu/osgeorel/">NCSU OSGeoREL</a>, developed during GSoC 2014.
+
+<h2>SEE ALSO</h2>
+<em>
+<a href="http://grass.osgeo.org/grass71/manuals/r.slope.aspect.html">r.flow</a>
+</em>
+
+<p><i>Last changed: $Date$</i>
Property changes on: grass-addons/grass7/raster3d/r3.gradient/r3.gradient.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster3d/r3.gradient/r3gradient_structs.h
===================================================================
--- grass-addons/grass7/raster3d/r3.gradient/r3gradient_structs.h (rev 0)
+++ grass-addons/grass7/raster3d/r3.gradient/r3gradient_structs.h 2014-08-06 02:46:37 UTC (rev 61547)
@@ -0,0 +1,13 @@
+#ifndef R3GRADIENT_STRUCTS_H
+#define R3GRADIENT_STRUCTS_H
+
+#include <grass/raster3d.h>
+
+struct Gradient_block {
+ RASTER3D_Array_double input;
+ RASTER3D_Array_double dx;
+ RASTER3D_Array_double dy;
+ RASTER3D_Array_double dz;
+};
+
+#endif // R3GRADIENT_STRUCTS_H
Property changes on: grass-addons/grass7/raster3d/r3.gradient/r3gradient_structs.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list