move r3.flow and r3.gradient to trunk

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@61630 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster3d/Makefile

@@ -2,6 +2,8 @@ MODULE_TOPDIR = ..
 
 SUBDIRS = \
 	r3.cross.rast \
+	r3.flow \
+	r3.gradient \
 	r3.gwflow \
 	r3.in.ascii \
 	r3.in.bin \

raster3d/r3.flow/Makefile

@@ -0,0 +1,15 @@
+MODULE_TOPDIR = ../..
+
+PROGRAMS=r3.flow test.r3flow 
+
+LIBES = $(RASTER3DLIB) $(RASTERLIB) $(GISLIB) $(VECTORLIB) $(DBMILIB)
+DEPENDENCIES = $(RASTER3DDEP) $(RASTERDEP) $(GISDEP) $(VECTORDEP) $(DBMIDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+r3_flow_OBJS = main.o flowline.o integrate.o interpolate.o voxel_traversal.o
+test_r3flow_OBJS = test_main.o flowline.o integrate.o interpolate.o voxel_traversal.o
+
+include $(MODULE_TOPDIR)/include/Make/Multi.make
+
+default: multi

raster3d/r3.flow/flowline.c

@@ -0,0 +1,214 @@
+/*!
+   \file flowline.c
+
+   \brief Generates flowlines as vector lines
+
+   (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.
+
+   \author Anna Petrasova
+ */
+#include <grass/raster3d.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+#include <grass/dbmi.h>
+
+#include "r3flow_structs.h"
+#include "integrate.h"
+#include "flowline.h"
+#include "voxel_traversal.h"
+
+static void write_segment(struct Map_info *flowline_vec,
+			  struct line_pnts *points, struct line_cats *cats,
+			  const double *point, int *cat)
+{
+    Vect_append_point(points, point[0], point[1], point[2]);
+
+    Vect_cat_set(cats, 1, *cat);
+    (*cat)++;
+    Vect_write_line(flowline_vec, GV_LINE, points, cats);
+    Vect_reset_line(points);
+    Vect_reset_cats(cats);
+    Vect_append_point(points, point[0], point[1], point[2]);
+}
+
+static void write_segment_db(struct field_info *finfo, dbDriver * driver,
+			     dbString * sql, const double velocity,
+			     double scalar_value, int write_scalar,
+			     const int cat)
+{
+    char buf[200];
+
+    if (write_scalar)
+	sprintf(buf, "insert into %s values ( %d, %e, %e )", finfo->table,
+		cat, scalar_value, velocity);
+    else
+	sprintf(buf, "insert into %s values ( %d, %e)", finfo->table, cat,
+		velocity);
+
+    db_set_string(sql, buf);
+
+    if (db_execute_immediate(driver, sql) != DB_OK) {
+	G_fatal_error(_("Unable to insert new record: '%s'"),
+		      db_get_string(sql));
+    }
+}
+
+static double get_scalar_value(RASTER3D_Region * region,
+			       struct Gradient_info *gradient_info,
+			       double north, double east, double top)
+{
+    int col, row, depth;
+    double val;
+
+    Rast3d_location2coord(region, north, east, top, &col, &row, &depth);
+    Rast3d_get_value(gradient_info->scalar_map, col, row, depth, &val,
+		     DCELL_TYPE);
+
+    return val;
+}
+
+/*!
+   \brief Computes flowline by integrating velocity field.
+
+   \param region pointer to current 3D region
+   \param seed starting seed (point)
+   \param velocity_field pointer to array of 3 3D raster maps
+   \param integration pointer to integration struct
+   \param flowline_vec pointer to Map_info struct of flowline vector
+   \param cats pointer to line_cats struct of flowline vector
+   \param points pointer to line_pnts struct of flowline vector
+   \param[in,out] cat starting category of the newly created flow line
+   \param if_table TRUE if attribute table should be created and filled
+ */
+void compute_flowline(RASTER3D_Region * region, const struct Seed *seed,
+		      struct Gradient_info *gradient_info,
+		      RASTER3D_Map * flowacc, struct Integration *integration,
+		      struct Map_info *flowline_vec, struct line_cats *cats,
+		      struct line_pnts *points, int *cat, int if_table,
+		      struct field_info *finfo, dbDriver * driver)
+{
+    int i, j, count;
+    double delta_t;
+    double velocity_norm;
+    double point[3], new_point[3];
+    double vel_x, vel_y, vel_z;
+    double min_step, max_step;
+    int col, row, depth;
+    int last_col, last_row, last_depth;
+    int coor_diff;
+    FCELL value;
+    double scalar_value;
+    int *trav_coords;
+    int size, trav_count;
+    dbString sql;
+    double velocity;
+
+    point[0] = seed->x;
+    point[1] = seed->y;
+    point[2] = seed->z;
+
+    last_col = last_row = last_depth = -1;
+
+    size = 5;
+    scalar_value = 0;
+    trav_coords = G_malloc(3 * size * sizeof(int));
+
+    if (seed->flowline) {
+	/* append first point */
+	Vect_append_point(points, seed->x, seed->y, seed->z);
+	db_init_string(&sql);
+    }
+    count = 1;
+    while (count <= integration->limit) {
+	if (get_velocity(region, gradient_info, point[0], point[1], point[2],
+			 &vel_x, &vel_y, &vel_z) < 0)
+	    break;		/* outside region */
+	velocity_norm = norm(vel_x, vel_y, vel_z);
+
+	if (velocity_norm <= VELOCITY_EPSILON)
+	    break;		/* zero velocity means end of propagation */
+
+	/* convert to time */
+	delta_t = get_time_step(integration->unit, integration->step,
+				velocity_norm, integration->cell_size);
+
+	/* integrate */
+	min_step = get_time_step("cell", integration->min_step, velocity_norm,
+				 integration->cell_size);
+	max_step = get_time_step("cell", integration->max_step, velocity_norm,
+				 integration->cell_size);
+	delta_t *= (integration->actual_direction == FLOWDIR_UP ? 1 : -1);
+	if (rk45_integrate_next
+	    (region, gradient_info, point, new_point,
+	     &delta_t, &velocity, min_step, max_step,
+	     integration->max_error) < 0)
+	    break;
+
+	if (seed->flowline) {
+	    if (if_table) {
+		write_segment(flowline_vec, points, cats, new_point, cat);
+		if (gradient_info->compute_gradient)
+		    scalar_value = get_scalar_value(region, gradient_info,
+						    point[1], point[0],
+						    point[2]);
+		write_segment_db(finfo, driver, &sql, velocity, scalar_value,
+				 gradient_info->compute_gradient, *cat - 1);
+	    }
+	    else
+		Vect_append_point(points, point[0], point[1], point[2]);
+	}
+	if (seed->flowaccum) {
+	    Rast3d_location2coord(region, new_point[1], new_point[0],
+				  new_point[2], &col, &row, &depth);
+	    if (!(last_col == col && last_row == row && last_depth == depth)) {
+		value = Rast3d_get_float(flowacc, col, row, depth);
+		Rast3d_put_float(flowacc, col, row, depth, value + 1);
+		if (last_col >= 0) {
+		    coor_diff = (abs(last_col - col) + abs(last_row - row) +
+				 abs(last_depth - depth));
+		    /* if not run for the 1. time and previous and next point coordinates
+		       differ by more than 1 voxel coordinate */
+		    if (coor_diff > 1) {
+			traverse(region, point, new_point, trav_coords, &size,
+				 &trav_count);
+			for (j = 0; j < trav_count; j++) {
+			    value =
+				Rast3d_get_float(flowacc,
+						 trav_coords[3 * j + 0],
+						 trav_coords[3 * j + 1],
+						 trav_coords[3 * j + 2]);
+			    Rast3d_put_float(flowacc, trav_coords[3 * j + 0],
+					     trav_coords[3 * j + 1],
+					     trav_coords[3 * j + 2],
+					     value + 1);
+			}
+		    }
+		}
+		last_col = col;
+		last_row = row;
+		last_depth = depth;
+	    }
+	}
+	for (i = 0; i < 3; i++) {
+	    point[i] = new_point[i];
+	}
+	count++;
+
+    }
+    if (seed->flowline) {
+	if (points->n_points > 1) {
+	    Vect_cat_set(cats, 1, *cat);
+	    (*cat)++;
+	    Vect_write_line(flowline_vec, GV_LINE, points, cats);
+	    G_debug(1, "Flowline ended after %d steps", count - 1);
+	}
+	Vect_reset_line(points);
+	Vect_reset_cats(cats);
+	db_free_string(&sql);
+    }
+    G_free(trav_coords);
+}

raster3d/r3.flow/flowline.h

@@ -0,0 +1,18 @@
+#ifndef FLOWLINE_H
+#define FLOWLINE_H
+
+#include <grass/raster3d.h>
+#include <grass/vector.h>
+
+#include "r3flow_structs.h"
+
+static const double VELOCITY_EPSILON = 1e-8;
+
+void compute_flowline(RASTER3D_Region * region, const struct Seed *seed,
+		      struct Gradient_info *gradient_info,
+		      RASTER3D_Map * flowacc,
+		      struct Integration *integration,
+		      struct Map_info *flowline_vec, struct line_cats *cats,
+		      struct line_pnts *points, int *cat, int if_table,
+		      struct field_info *finfo, dbDriver *driver);
+#endif // FLOWLINE_H

raster3d/r3.flow/integrate.c

@@ -0,0 +1,227 @@
+/*!
+   \file integrate.c
+
+   \brief Flowline integration using Runge-Kutta 45 method
+   with adaptive step size control. Implementation based
+   on VTK class vtkRungeKutta45.
+
+   (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.
+
+   \author Anna Petrasova
+ */
+
+#include <math.h>
+#include <float.h>
+#include <string.h>
+#include <grass/raster3d.h>
+
+#include "interpolate.h"
+#include "integrate.h"
+#include "r3flow_structs.h"
+
+/*!
+   \brief Computes vector norm.
+
+   \param x,y,z vector components
+
+   \return norm
+ */
+double norm(const double x, const double y, const double z)
+{
+    return sqrt(x * x + y * y + z * z);
+}
+
+/*!
+   \brief Converts integration step to time
+
+   \param unit type of unit
+   \param step step value
+   \param velocity velocity value
+   \param cell_size size of voxel diagonal
+
+   \return time step
+ */
+double get_time_step(const char *unit, const double step,
+		     const double velocity, const double cell_size)
+{
+    if (strcmp(unit, "time") == 0)
+	return step;
+    else if (strcmp(unit, "length") == 0)
+	return step / velocity;
+    else			/* cell */
+	return (step * cell_size) / velocity;
+}
+
+
+int get_velocity(RASTER3D_Region * region, struct Gradient_info *gradient_info,
+		 const double x, const double y, const double z,
+		 double *vel_x, double *vel_y, double *vel_z)
+{
+    if (gradient_info->compute_gradient)
+        return get_gradient(region, gradient_info, y, x, z, vel_x, vel_y, vel_z);
+
+    return interpolate_velocity(region, gradient_info->velocity_maps, y, x, z,
+				vel_x, vel_y, vel_z);
+}
+
+/*!
+   \brief Runge-Kutta45 (inner steps) 
+
+   \param region pointer to current 3D region
+   \param velocity_field pointer to array of 3 3D raster maps (velocity field)
+   \param point pointer to array of geographic coordinates of starting point
+   \param[out] next_point pointer to array of geographic coordinates of integrated point
+   \param delta_t integration time step
+   \param error pointer to array of error values
+
+   \return 0 success
+   \return -1 out of region or null values
+ */
+static int rk45_next(RASTER3D_Region * region, struct Gradient_info *gradient_info,
+		     const double *point, double *next_point,
+		     const double delta_t, double *velocity, double *error)
+{
+    double tmp1[6][3];		/* 3 is 3 dimensions, 6 is the number of k's */
+    double tmp_point[3];
+    double vel_x, vel_y, vel_z;
+    double sum_tmp;
+    int i, j, k;
+    double vel_sq;
+
+    if (get_velocity(region, gradient_info, point[0], point[1], point[2],
+		     &vel_x, &vel_y, &vel_z) < 0)
+	return -1;
+
+    tmp1[0][0] = vel_x;
+    tmp1[0][1] = vel_y;
+    tmp1[0][2] = vel_z;
+
+    /* compute k's */
+    for (i = 1; i < 6; i++) {
+	for (j = 0; j < 3; j++) {	/* for each coordinate */
+	    sum_tmp = 0;
+	    for (k = 0; k < i; k++) {	/* k1; k1, k2; ... */
+		sum_tmp += B[i - 1][k] * tmp1[k][j];
+	    }
+	    tmp_point[j] = point[j] + delta_t * sum_tmp;
+	}
+	if (get_velocity
+	    (region, gradient_info, tmp_point[0], tmp_point[1], tmp_point[2],
+	     &vel_x, &vel_y, &vel_z) < 0)
+	    return -1;
+
+	tmp1[i][0] = vel_x;
+	tmp1[i][1] = vel_y;
+	tmp1[i][2] = vel_z;
+    }
+
+    vel_sq = 0;
+    /* compute next point */
+    for (j = 0; j < 3; j++) {
+	sum_tmp = 0;
+	for (i = 0; i < 6; i++) {
+	    sum_tmp += C[i] * tmp1[i][j];
+	}
+	next_point[j] = point[j] + delta_t * sum_tmp;
+	vel_sq += sum_tmp * sum_tmp;
+    }
+    *velocity = sqrt(vel_sq);
+
+    if (!Rast3d_is_valid_location
+	(region, next_point[1], next_point[0], next_point[2]))
+	return -1;
+
+    /* compute error vector */
+    for (j = 0; j < 3; j++) {
+	sum_tmp = 0;
+	for (i = 0; i < 6; i++) {
+	    sum_tmp += DC[i] * tmp1[i][j];
+	}
+	error[j] = delta_t * sum_tmp;
+    }
+
+    return 0;
+}
+
+/*!
+   \brief Runge-Kutta45 with step size control 
+
+   \param region pointer to current 3D region
+   \param velocity_field pointer to array of 3 3D raster maps (velocity field)
+   \param point pointer to array of geographic coordinates of starting point
+   \param[out] next_point pointer to array of geographic coordinates of integrated point
+   \param[in,out] delta_t integration time step
+   \param min_step minimum time step
+   \param max_step maximum time step
+
+   \return 0 success
+   \return -1 out of region or null values
+ */
+int rk45_integrate_next(RASTER3D_Region * region,
+			struct Gradient_info *gradient_info, const double *point,
+			double *next_point, double *delta_t, double *velocity,
+			const double min_step, const double max_step,
+			const double max_error)
+{
+    double estimated_error;
+    double error_ratio;
+    double error[3];
+    double tmp, tmp2;
+    int do_break;
+
+    estimated_error = DBL_MAX;
+
+    /* check if min_step < delta_t < max_step */
+    if (fabs(*delta_t) < min_step)
+	*delta_t = min_step * (*delta_t > 0 ? 1 : -1);
+    else if (fabs(*delta_t) > max_step)
+	*delta_t = max_step * (*delta_t > 0 ? 1 : -1);
+
+    /* try to iteratively decrease error to less than max error */
+    while (estimated_error > max_error) {
+	/* compute next point and get estimated error */
+	if (rk45_next
+	    (region, gradient_info, point, next_point, *delta_t,
+	     velocity, error) == 0)
+	    estimated_error = norm(error[0], error[1], error[2]);
+	else
+	    return -1;
+
+	/* compute new step size (empirically) */
+	error_ratio = estimated_error / max_error;
+	if (error_ratio == 0.0)
+	    tmp = *delta_t > 0 ? min_step : -min_step;
+	else if (error_ratio > 1)
+	    tmp = 0.9 * *delta_t * pow(error_ratio, -0.25);
+	else
+	    tmp = 0.9 * *delta_t * pow(error_ratio, -0.2);
+	tmp2 = fabs(tmp);
+
+	do_break = FALSE;
+
+	/* adjust new step size to be within min max limits */
+	if (tmp2 > max_step) {
+	    *delta_t = max_step * (*delta_t > 0 ? 1 : -1);
+	    do_break = TRUE;
+	}
+	else if (tmp2 < min_step) {
+	    *delta_t = min_step * (*delta_t > 0 ? 1 : -1);
+	    do_break = TRUE;
+	}
+	else
+	    *delta_t = tmp;
+
+	/* break when the adjustment was needed (not sure why) */
+	if (do_break) {
+	    if (rk45_next(region, gradient_info, point, next_point,
+			  *delta_t, velocity, error) < 0)
+		return -1;
+	    break;
+	}
+    }
+    return 0;
+}

raster3d/r3.flow/integrate.h

@@ -0,0 +1,35 @@
+#ifndef INTEGRATE_H
+#define INTEGRATE_H
+
+#include <grass/raster3d.h>
+
+#include "r3flow_structs.h"
+
+/* Cash-Karp parameters */
+static const double B[5][5] = { {1. / 5, 0, 0, 0, 0},
+{3. / 40, 9. / 40, 0, 0, 0},
+{3. / 10, -9. / 10, 6. / 5, 0, 0},
+{-11. / 54, 5. / 2, -70. / 27, 35. / 27, 0},
+{1631. / 55296, 175. / 512, 575. / 13824,
+ 44275. / 110592, 253. / 4096}
+};
+static const double C[6] =
+    { 37. / 378, 0, 250. / 621, 125. / 594, 0, 512. / 1771 };
+static const double DC[6] = { 37. / 378 - 2825. / 27648, 0,
+    250. / 621 - 18575. / 48384, 125. / 594 - 13525. / 55296,
+    -277. / 14336, 512. / 1771 - 1. / 4
+};
+
+double norm(const double x, const double y, const double z);
+int get_velocity(RASTER3D_Region * region, struct Gradient_info *gradient_info,
+		 const double x, const double y, const double z,
+		 double *vel_x, double *vel_y, double *vel_z);
+double get_time_step(const char *unit, const double step,
+		     const double velocity, const double cell_size);
+int rk45_integrate_next(RASTER3D_Region * region,
+			struct Gradient_info *gradient_info, const double *point,
+			double *next_point, double *delta_t,
+			double *velocity, const double min_step,
+			const double max_step, const double max_error);
+
+#endif // INTEGRATE_H

raster3d/r3.flow/interpolate.c

@@ -0,0 +1,337 @@
+/*!
+   \file interpolate.c
+
+   \brief Trilinear interpolation
+
+   (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.
+
+   \author Anna Petrasova
+ */
+
+#include <grass/gis.h>
+#include <grass/raster3d.h>
+
+#include "r3flow_structs.h"
+#include "interpolate.h"
+
+/*!
+   \brief Finds 8 nearest voxels from a point.
+
+   \param region pointer to current 3D region
+   \param north,east,top geographic coordinates
+   \param[out] x,y,z pointer to indices of neighbouring voxels
+ */
+static void find_nearest_voxels(RASTER3D_Region * region,
+				const double north, const double east,
+				const double top, int *x, int *y, int *z)
+{
+    double n_minus, e_minus, t_minus;
+    double n_plus, e_plus, t_plus;
+
+    n_minus = north - region->ns_res / 2;
+    n_plus = north + region->ns_res / 2;
+    e_minus = east - region->ew_res / 2;
+    e_plus = east + region->ew_res / 2;
+    t_minus = top - region->tb_res / 2;
+    t_plus = top + region->tb_res / 2;
+
+    Rast3d_location2coord(region, n_minus, e_minus, t_minus, x++, y++, z++);
+    Rast3d_location2coord(region, n_minus, e_plus, t_minus, x++, y++, z++);
+    Rast3d_location2coord(region, n_plus, e_minus, t_minus, x++, y++, z++);
+    Rast3d_location2coord(region, n_plus, e_plus, t_minus, x++, y++, z++);
+    Rast3d_location2coord(region, n_minus, e_minus, t_plus, x++, y++, z++);
+    Rast3d_location2coord(region, n_minus, e_plus, t_plus, x++, y++, z++);
+    Rast3d_location2coord(region, n_plus, e_minus, t_plus, x++, y++, z++);
+    Rast3d_location2coord(region, n_plus, e_plus, t_plus, x, y, z);
+}
+
+/*!
+   \brief Trilinear interpolation
+
+   Computation is based on the sum of values of nearest voxels
+   weighted by the relative distance of a point
+   to the center of the nearest voxels.
+
+   \param array_values pointer to values of 8 (neigboring) voxels
+   \param x,y,z relative coordinates (0 - 1)
+   \param[out] interpolated pointer to the array (of size 3) of interpolated values
+ */
+static void trilinear_interpolation(const double *array_values,
+				    const double x, const double y,
+				    const double z, double *interpolated)
+{
+    int i, j;
+    double rx, ry, rz, value;
+    double weights[8];
+
+    rx = 1 - x;
+    ry = 1 - y;
+    rz = 1 - z;
+    weights[0] = rx * ry * rz;
+    weights[1] = x * ry * rz;
+    weights[2] = rx * y * rz;
+    weights[3] = x * y * rz;
+    weights[4] = rx * ry * z;
+    weights[5] = x * ry * z;
+    weights[6] = rx * y * z;
+    weights[7] = x * y * z;
+
+
+    /* weighted sum of surrounding values */
+    for (i = 0; i < 3; i++) {
+	value = 0;
+	for (j = 0; j < 8; j++) {
+	    value += weights[j] * array_values[i * 8 + j];
+	}
+	interpolated[i] = value;
+    }
+}
+
+/*!
+   \brief Converts geographic to relative coordinates
+
+   Converts geographic to relative coordinates
+   which are needed for trilinear interpolation.
+
+
+   \param region pointer to current 3D region
+   \param north,east,top geographic coordinates
+   \param[out] x,y,z pointer to relative coordinates (0 - 1)
+ */
+static void get_relative_coords_for_interp(RASTER3D_Region * region,
+					   const double north,
+					   const double east,
+					   const double top, double *x,
+					   double *y, double *z)
+{
+    int col, row, depth;
+    double temp;
+
+    Rast3d_location2coord(region, north, east, top, &col, &row, &depth);
+
+    /* x */
+    temp = east - region->west - col * region->ew_res;
+    *x = (temp > region->ew_res / 2 ?
+	  temp - region->ew_res / 2 : temp + region->ew_res / 2)
+	/ region->ew_res;
+    /* y */
+    temp = north - region->south - (region->rows - row - 1) * region->ns_res;
+    *y = (temp > region->ns_res / 2 ?
+	  temp - region->ns_res / 2 : temp + region->ns_res / 2)
+	/ region->ns_res;
+    /* z */
+    temp = top - region->bottom - depth * region->tb_res;
+    *z = (temp > region->tb_res / 2 ?
+	  temp - region->tb_res / 2 : temp + region->tb_res / 2)
+	/ region->tb_res;
+}
+
+/*!
+   \brief Interpolates velocity at a given point.
+
+   \param region pointer to current 3D region
+   \param map pointer to array of 3 3D raster maps (velocity field)
+   \param north,east,top geographic coordinates
+   \param[out] vel_x,vel_y,vel_z interpolated velocity
+
+   \return 0 success
+   \return -1 out of region
+ */
+int interpolate_velocity(RASTER3D_Region * region, RASTER3D_Map ** map,
+			 const double north, const double east,
+			 const double top, double *vel_x, double *vel_y,
+			 double *vel_z)
+{
+    int i, j;
+    double values[24];		/* 3 x 8, 3 dimensions, 8 neighbors */
+    double value;
+    double interpolated[3];
+    int x[8], y[8], z[8];
+    double rel_x, rel_y, rel_z;
+    int type;
+
+    /* check if we are out of region, any array should work */
+    if (!Rast3d_is_valid_location(region, north, east, top))
+	return -1;
+
+    find_nearest_voxels(region, north, east, top, x, y, z);
+    /* get values of the nearest cells */
+    for (i = 0; i < 3; i++) {
+	type = Rast3d_tile_type_map(map[i]);
+	for (j = 0; j < 8; j++) {
+
+	    Rast3d_get_value_region(map[i], x[j], y[j], z[j], &value, type);
+	    if (Rast_is_null_value(&value, type))
+		values[i * 8 + j] = 0;
+	    else
+		values[i * 8 + j] = value;
+	}
+    }
+
+    /* compute weights */
+    get_relative_coords_for_interp(region, north, east, top,
+				   &rel_x, &rel_y, &rel_z);
+
+    trilinear_interpolation(values, rel_x, rel_y, rel_z, interpolated);
+    *vel_x = interpolated[0];
+    *vel_y = interpolated[1];
+    *vel_z = interpolated[2];
+
+    return 0;
+}
+/*!
+   \brief Computes gradient for a point.
+
+   \param region pointer to current 3D region
+   \param gradient_info struct which remembers values
+          related to gradient computation to avoid computation every time
+   \param north,east,top geographic coordinates
+   \param[out] vel_x,vel_y,vel_z interpolated gradient components
+
+   \return 0 success
+   \return -1 out of region
+ */
+int get_gradient(RASTER3D_Region * region,
+		 struct Gradient_info *gradient_info, const double north,
+		 const double east, const double top, double *vel_x, double *vel_y,
+		 double *vel_z)
+{
+
+    int i, d, r, c, count;
+    int near_x[8], near_y[8], near_z[8];
+    int minx, maxx, miny, maxy, minz, maxz;
+    int xshift, yshift, zshift;
+    double rel_x, rel_y, rel_z;
+    double scalar_map_array[64];
+    double grad_x_map_array[64], grad_y_map_array[64], grad_z_map_array[64];
+    RASTER3D_Array_double array;
+    RASTER3D_Array_double grad_x, grad_y, grad_z;
+    RASTER3D_Array_double *grad_xyz[3];
+    double step[3];
+    double interpolated[3];
+
+    step[0] = region->ew_res;
+    step[1] = region->ns_res;
+    step[2] = region->tb_res;
+
+    array.array = scalar_map_array;
+    array.sx = array.sy = array.sz = 4;
+    grad_x.array = grad_x_map_array;
+    grad_y.array = grad_y_map_array;
+    grad_z.array = grad_z_map_array;
+    grad_x.sx = grad_x.sy = grad_x.sz = 4;
+    grad_y.sx = grad_y.sy = grad_y.sz = 4;
+    grad_z.sx = grad_z.sy = grad_z.sz = 4;
+
+    find_nearest_voxels(region, north, east, top, near_x, near_y, near_z);
+    minx = near_x[0];
+    maxx = near_x[7];
+    miny = near_y[7];
+    maxy = near_y[0];
+    minz = near_z[0];
+    maxz = near_z[7];
+
+    /* position of x, y, z neighboring voxels */
+    if (!gradient_info->initialized ||
+	(gradient_info->neighbors_pos[0] != minx ||
+	 gradient_info->neighbors_pos[1] != miny ||
+	 gradient_info->neighbors_pos[2] != minz)) {
+
+	gradient_info->neighbors_pos[0] = minx;
+	gradient_info->neighbors_pos[1] = miny;
+	gradient_info->neighbors_pos[2] = minz;
+	gradient_info->initialized = TRUE;
+
+	/* just to be sure, we check that at least one voxel is inside */
+	if (maxx < 0 || minx >= region->cols ||
+	    maxy < 0 || miny >= region->rows ||
+	    maxz < 0 || minz >= region->depths)
+	    return -1;
+
+	/* these if's are here to handle edge cases
+	   min is changed to represent the min coords of the 4x4x4 array
+	   from which the gradient will be computed
+	   shift is relative position of the neighbors within this 4x4x4 array */
+	if (minx == 0 || minx == -1) {
+	    xshift = minx;
+	    minx = 0;
+	}
+	else if (maxx >= region->cols - 1) {
+	    minx = maxx < region->cols ? maxx - 3 : maxx - 4;
+	    xshift = maxx < region->cols ? 2 : 3;
+	}
+	else {
+	    minx -= 1;
+	    xshift = 1;
+	}
+
+	if (miny == 0 || miny == -1) {
+	    yshift = miny;
+	    miny = 0;
+	}
+	else if (maxy >= region->rows - 1) {
+	    miny = maxy < region->rows ? maxy - 3 : maxy - 4;
+	    yshift = maxy < region->rows ? 2 : 3;
+	}
+	else {
+	    miny -= 1;
+	    yshift = 1;
+	}
+
+	if (minz == 0 || minz == -1) {
+	    zshift = minz;
+	    minz = 0;
+	}
+	else if (maxz >= region->depths - 1) {
+	    minz = maxz < region->depths ? maxz - 3 : maxz - 4;
+	    zshift = maxz < region->depths ? 2 : 3;
+	}
+	else {
+	    minz -= 1;
+	    zshift = 1;
+	}
+
+	/* get the 4x4x4 block of the array */
+	Rast3d_get_block(gradient_info->scalar_map, minx, miny, minz,
+			 4, 4, 4, array.array, DCELL_TYPE);
+	Rast3d_gradient_double(&array, step, &grad_x, &grad_y, &grad_z);
+	grad_xyz[0] = &grad_x;
+	grad_xyz[1] = &grad_y;
+	grad_xyz[2] = &grad_z;
+	/* go through x, y, z and all 8 neighbors and store their value
+	   if the voxel is outside, add 0 (weight) */
+	for (i = 0; i < 3; i++) {
+	    count = 0;
+	    for (d = 0; d < 2; d++)
+		for (r = 1; r > -1; r--)
+		    for (c = 0; c < 2; c++) {
+			if (d + zshift < 0 || d + zshift > 3 ||
+			    r + yshift < 0 || r + yshift > 3 ||
+			    c + xshift < 0 || c + xshift > 3)
+			    gradient_info->neighbors_values[i * 8 + count] =
+				0;
+			else
+			    gradient_info->neighbors_values[i * 8 + count] =
+				RASTER3D_ARRAY_ACCESS(grad_xyz[i], c + xshift,
+						      r + yshift, d + zshift);
+			count++;
+		    }
+	}
+    }
+    get_relative_coords_for_interp(region, north, east, top, &rel_x, &rel_y, &rel_z);
+    trilinear_interpolation(gradient_info->neighbors_values,
+			    rel_x, rel_y, rel_z, interpolated);
+
+    *vel_x = interpolated[0];
+    *vel_y = interpolated[1];
+    *vel_z = interpolated[2];
+
+    return 0;
+
+
+
+}

raster3d/r3.flow/interpolate.h

@@ -0,0 +1,16 @@
+#ifndef INTERPOLATE_H
+#define INTERPOLATE_H
+
+#include <grass/raster3d.h>
+
+#include "r3flow_structs.h"
+
+int interpolate_velocity(RASTER3D_Region * region, RASTER3D_Map ** map,
+			 const double north, const double east,
+			 const double top, double *vel_x, double *vel_y,
+			 double *vel_z);
+int get_gradient(RASTER3D_Region * region,
+		 struct Gradient_info *gradient_info, const double north,
+		 const double east, const double top, double *vel_x, double *vel_y,
+		 double *vel_z);
+#endif // INTERPOLATE_H

raster3d/r3.flow/main.c

@@ -0,0 +1,512 @@
+
+/****************************************************************************
+ * 
+ * MODULE:       r3.flow     
+ * AUTHOR(S):    Anna Petrasova kratochanna <at> gmail <dot> com
+ * PURPOSE:      Computes 3D flow lines and flow accumulation based on 3D raster map(s)
+ * 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 <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/raster3d.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/vector.h>
+#include <grass/dbmi.h>
+#include <grass/glocale.h>
+
+#include "r3flow_structs.h"
+#include "flowline.h"
+
+static void create_table(struct Map_info *flowline_vec,
+			 struct field_info **f_info, dbDriver ** driver,
+			 int write_scalar)
+{
+    dbString sql;
+    char buf[200];
+    dbDriver *drvr;
+    struct field_info *fi;
+
+    db_init_string(&sql);
+    fi = Vect_default_field_info(flowline_vec, 1, NULL, GV_1TABLE);
+    *f_info = fi;
+    Vect_map_add_dblink(flowline_vec, 1, NULL, fi->table, GV_KEY_COLUMN,
+			fi->database, fi->driver);
+    drvr = db_start_driver_open_database(fi->driver,
+					 Vect_subst_var(fi->database,
+							flowline_vec));
+    if (drvr == NULL) {
+	G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
+		      Vect_subst_var(fi->database, flowline_vec), fi->driver);
+    }
+    *driver = drvr;
+    if (write_scalar)
+	sprintf(buf, "create table %s (cat integer, input double precision, "
+		"velocity double precision)", fi->table);
+    else
+	sprintf(buf,
+		"create table %s (cat integer, velocity double precision)",
+		fi->table);
+    db_set_string(&sql, buf);
+
+    db_begin_transaction(drvr);
+    /* Create table */
+    if (db_execute_immediate(drvr, &sql) != DB_OK) {
+	G_fatal_error(_("Unable to create table: %s"), db_get_string(&sql));
+    }
+    if (db_create_index2(drvr, fi->table, fi->key) != DB_OK)
+	G_warning(_("Unable to create index for table <%s>, key <%s>"),
+		  fi->table, fi->key);
+    /* Grant */
+    if (db_grant_on_table
+	(drvr, fi->table, DB_PRIV_SELECT, DB_GROUP | DB_PUBLIC) != DB_OK) {
+	G_fatal_error(_("Unable to grant privileges on table <%s>"),
+		      fi->table);
+    }
+}
+
+static void check_vector_input_maps(struct Option *vector_opt,
+				    struct Option *seed_opt)
+{
+    int i;
+
+    /* Check for velocity components maps. */
+    if (vector_opt->answers != NULL) {
+	for (i = 0; i < 3; i++) {
+	    if (vector_opt->answers[i] != NULL) {
+		if (NULL == G_find_raster3d(vector_opt->answers[i], ""))
+		    Rast3d_fatal_error(_("3D raster map <%s> not found"),
+				       vector_opt->answers[i]);
+	    }
+	    else {
+		Rast3d_fatal_error(_("Please provide three 3D raster maps"));
+	    }
+	}
+    }
+
+    if (seed_opt->answer != NULL) {
+	if (NULL == G_find_vector2(seed_opt->answer, ""))
+	    G_fatal_error(_("Vector seed map <%s> not found"),
+			  seed_opt->answer);
+    }
+
+}
+
+static void load_input_raster3d_maps(struct Option *scalar_opt,
+				     struct Option *vector_opt,
+				     struct Gradient_info *gradient_info,
+				     RASTER3D_Region * region)
+{
+    int i;
+
+    if (scalar_opt->answer) {
+	gradient_info->scalar_map =
+	    Rast3d_open_cell_old(scalar_opt->answer,
+				 G_find_raster3d(scalar_opt->answer, ""),
+				 region, RASTER3D_TILE_SAME_AS_FILE,
+				 RASTER3D_USE_CACHE_DEFAULT);
+	if (!gradient_info->scalar_map)
+	    Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
+			       scalar_opt->answer);
+	gradient_info->compute_gradient = TRUE;
+    }
+    else {
+	for (i = 0; i < 3; i++) {
+	    gradient_info->velocity_maps[i] =
+		Rast3d_open_cell_old(vector_opt->answers[i],
+				     G_find_raster3d(vector_opt->answers[i],
+						     ""), region,
+				     RASTER3D_TILE_SAME_AS_FILE,
+				     RASTER3D_USE_CACHE_DEFAULT);
+
+	    if (!gradient_info->velocity_maps[i])
+		Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
+				   vector_opt->answers[i]);
+	}
+	gradient_info->compute_gradient = FALSE;
+    }
+}
+
+static void init_flowaccum(RASTER3D_Region * region, RASTER3D_Map * flowacc)
+{
+    int c, r, d;
+
+    for (d = 0; d < region->depths; d++)
+	for (r = 0; r < region->rows; r++)
+	    for (c = 0; c < region->cols; c++)
+		if (Rast3d_put_float(flowacc, c, r, d, 0) != 1)
+		    Rast3d_fatal_error(_("init_flowaccum: error in Rast3d_put_float"));
+}
+
+int main(int argc, char *argv[])
+{
+    struct Option *vector_opt, *seed_opt, *flowlines_opt, *flowacc_opt,
+	*scalar_opt, *unit_opt, *step_opt, *limit_opt, *skip_opt, *dir_opt,
+	*error_opt;
+    struct Flag *table_fl;
+    struct GModule *module;
+    RASTER3D_Region region;
+    RASTER3D_Map *flowacc;
+    struct Integration integration;
+    struct Seed seed;
+    struct Gradient_info gradient_info;
+    struct Map_info seed_Map;
+    struct line_pnts *seed_points;
+    struct line_cats *seed_cats;
+    struct Map_info fl_map;
+    struct line_cats *fl_cats;	/* for flowlines */
+    struct line_pnts *fl_points;	/* for flowlines */
+    struct field_info *finfo;
+    dbDriver *driver;
+    int cat;			/* cat of flowlines */
+    int if_table;
+    int i, r, c, d;
+    char *desc;
+    int n_seeds, seed_count, ltype;
+    int skip[3];
+
+    G_gisinit(argv[0]);
+    module = G_define_module();
+    G_add_keyword(_("raster3d"));
+    G_add_keyword(_("voxel"));
+    G_add_keyword(_("flowline"));
+    module->description =
+	_("Computes 3D flow lines and 3D flow accumulation.");
+
+
+    scalar_opt = G_define_standard_option(G_OPT_R3_INPUT);
+    scalar_opt->required = NO;
+    scalar_opt->guisection = _("Input");
+
+    vector_opt = G_define_standard_option(G_OPT_R3_INPUTS);
+    vector_opt->key = "vector_field";
+    vector_opt->required = NO;
+    vector_opt->description = _("Names of three 3D raster maps describing "
+				"x, y, z components of vector field");
+    vector_opt->guisection = _("Input");
+
+    seed_opt = G_define_standard_option(G_OPT_V_INPUT);
+    seed_opt->required = NO;
+    seed_opt->key = "seed_points";
+    seed_opt->description = _("If no map is provided, "
+			      "flow lines are generated "
+			      "from each cell of the input 3D raster");
+    seed_opt->label = _("Name of vector map with points "
+			"from which flow lines are generated");
+    seed_opt->guisection = _("Input");
+
+    flowlines_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+    flowlines_opt->key = "flowline";
+    flowlines_opt->required = NO;
+    flowlines_opt->description = _("Name for vector map of flow lines");
+    flowlines_opt->guisection = _("Output");
+
+    flowacc_opt = G_define_standard_option(G_OPT_R3_OUTPUT);
+    flowacc_opt->key = "flowaccumulation";
+    flowacc_opt->required = NO;
+    flowacc_opt->description =
+	_("Name for output flowaccumulation 3D raster");
+    flowacc_opt->guisection = _("Output");
+
+    unit_opt = G_define_option();
+    unit_opt->key = "unit";
+    unit_opt->type = TYPE_STRING;
+    unit_opt->required = NO;
+    unit_opt->answer = "cell";
+    unit_opt->options = "time,length,cell";
+    desc = NULL;
+    G_asprintf(&desc,
+	       "time;%s;"
+	       "length;%s;"
+	       "cell;%s",
+	       _("elapsed time"),
+	       _("length in map units"), _("length in cells (voxels)"));
+    unit_opt->descriptions = desc;
+    unit_opt->label = _("Unit of integration step");
+    unit_opt->description = _("Default unit is cell");
+    unit_opt->guisection = _("Integration");
+
+    step_opt = G_define_option();
+    step_opt->key = "step";
+    step_opt->type = TYPE_DOUBLE;
+    step_opt->required = NO;
+    step_opt->answer = "0.25";
+    step_opt->label = _("Integration step in selected unit");
+    step_opt->description = _("Default step is 0.25 cell");
+    step_opt->guisection = _("Integration");
+
+    limit_opt = G_define_option();
+    limit_opt->key = "limit";
+    limit_opt->type = TYPE_INTEGER;
+    limit_opt->required = NO;
+    limit_opt->answer = "2000";
+    limit_opt->description = _("Maximum number of steps");
+    limit_opt->guisection = _("Integration");
+
+    error_opt = G_define_option();
+    error_opt->key = "max_error";
+    error_opt->type = TYPE_DOUBLE;
+    error_opt->required = NO;
+    error_opt->answer = "1e-5";
+    error_opt->label = _("Maximum error of integration");
+    error_opt->description = _("Influences step, increase maximum error "
+			       "to allow bigger steps");
+    error_opt->guisection = _("Integration");
+
+    skip_opt = G_define_option();
+    skip_opt->key = "skip";
+    skip_opt->type = TYPE_INTEGER;
+    skip_opt->required = NO;
+    skip_opt->multiple = YES;
+    skip_opt->description =
+	_("Number of cells between flow lines in x, y and z direction");
+
+    dir_opt = G_define_option();
+    dir_opt->key = "direction";
+    dir_opt->type = TYPE_STRING;
+    dir_opt->required = NO;
+    dir_opt->multiple = NO;
+    dir_opt->options = "up,down,both";
+    dir_opt->answer = "down";
+    dir_opt->description = _("Compute flowlines upstream, "
+			     "downstream or in both direction.");
+
+    table_fl = G_define_flag();
+    table_fl->key = 'a';
+    table_fl->description = _("Create and fill attribute table");
+
+    G_option_required(scalar_opt, vector_opt, NULL);
+    G_option_exclusive(scalar_opt, vector_opt, NULL);
+    G_option_required(flowlines_opt, flowacc_opt, NULL);
+    G_option_requires(seed_opt, flowlines_opt, NULL);
+    G_option_requires(table_fl, flowlines_opt);
+
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    driver = NULL;
+    finfo = NULL;
+
+    if_table = table_fl->answer ? TRUE : FALSE;
+
+    check_vector_input_maps(vector_opt, seed_opt);
+
+    Rast3d_init_defaults();
+    Rast3d_get_window(&region);
+
+    /* set up integration variables */
+    if (step_opt->answer) {
+	integration.step = atof(step_opt->answer);
+	integration.unit = unit_opt->answer;
+    }
+    else {
+	integration.unit = "cell";
+	integration.step = 0.25;
+    }
+    integration.max_error = atof(error_opt->answer);
+    integration.max_step = 5 * integration.step;
+    integration.min_step = integration.step / 5;
+    integration.limit = atof(limit_opt->answer);
+    if (strcmp(dir_opt->answer, "up") == 0)
+	integration.direction_type = FLOWDIR_UP;
+    else if (strcmp(dir_opt->answer, "down") == 0)
+	integration.direction_type = FLOWDIR_DOWN;
+    else
+	integration.direction_type = FLOWDIR_BOTH;
+
+
+    /* cell size is the diagonal */
+    integration.cell_size = sqrt(region.ns_res * region.ns_res +
+				 region.ew_res * region.ew_res +
+				 region.tb_res * region.tb_res);
+
+    /* set default skip if needed */
+    if (skip_opt->answers) {
+	for (i = 0; i < 3; i++) {
+	    if (skip_opt->answers[i] != NULL) {
+		skip[i] = atoi(skip_opt->answers[i]);
+	    }
+	    else {
+		G_fatal_error(_("Please provide 3 integer values for skip option."));
+	    }
+	}
+    }
+    else {
+	skip[0] = fmax(1, region.cols / 10);
+	skip[1] = fmax(1, region.rows / 10);
+	skip[2] = fmax(1, region.depths / 10);
+
+    }
+
+    /* open raster 3D maps of velocity components */
+    gradient_info.initialized = FALSE;
+    load_input_raster3d_maps(scalar_opt, vector_opt, &gradient_info, &region);
+
+
+    /* open new 3D raster map of flowacumulation */
+    if (flowacc_opt->answer) {
+	flowacc = Rast3d_open_new_opt_tile_size(flowacc_opt->answer,
+						RASTER3D_USE_CACHE_DEFAULT,
+						&region, FCELL_TYPE, 32);
+
+
+	if (!flowacc)
+	    Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
+			       flowacc_opt->answer);
+	init_flowaccum(&region, flowacc);
+    }
+
+    /* open new vector map of flowlines */
+    if (flowlines_opt->answer) {
+	fl_cats = Vect_new_cats_struct();
+	fl_points = Vect_new_line_struct();
+	if (Vect_open_new(&fl_map, flowlines_opt->answer, TRUE) < 0)
+	    G_fatal_error(_("Unable to create vector map <%s>"),
+			  flowlines_opt->answer);
+
+	Vect_hist_command(&fl_map);
+
+	if (if_table) {
+	    create_table(&fl_map, &finfo, &driver,
+			 gradient_info.compute_gradient);
+	}
+    }
+
+    n_seeds = 0;
+    /* open vector map of seeds */
+    if (seed_opt->answer) {
+	if (Vect_open_old2(&seed_Map, seed_opt->answer, "", "1") < 0)
+	    G_fatal_error(_("Unable to open vector map <%s>"),
+			  seed_opt->answer);
+	if (!Vect_is_3d(&seed_Map))
+	    G_fatal_error(_("Vector map <%s> is not 3D"), seed_opt->answer);
+
+	n_seeds = Vect_get_num_primitives(&seed_Map, GV_POINT);
+    }
+    if (flowacc_opt->answer || (!seed_opt->answer && flowlines_opt->answer)) {
+	if (flowacc_opt->answer)
+	    n_seeds += region.cols * region.rows * region.depths;
+	else {
+	    n_seeds += ceil(region.cols / (double)skip[0]) *
+		ceil(region.rows / (double)skip[1]) *
+		ceil(region.depths / (double)skip[2]);
+	}
+    }
+    G_debug(1, "Number of seeds is %d", n_seeds);
+
+    seed_count = 0;
+    cat = 1;
+    if (seed_opt->answer) {
+
+	seed_points = Vect_new_line_struct();
+	seed_cats = Vect_new_cats_struct();
+
+	/* compute flowlines from vector seed map */
+	while (TRUE) {
+	    ltype = Vect_read_next_line(&seed_Map, seed_points, seed_cats);
+	    if (ltype == -1) {
+		Vect_close(&seed_Map);
+		G_fatal_error(_("Error during reading seed vector map"));
+	    }
+	    else if (ltype == -2) {
+		break;
+	    }
+	    else if (ltype == GV_POINT) {
+		seed.x = seed_points->x[0];
+		seed.y = seed_points->y[0];
+		seed.z = seed_points->z[0];
+		seed.flowline = TRUE;
+		seed.flowaccum = FALSE;
+	    }
+	    G_percent(seed_count, n_seeds, 1);
+	    if (integration.direction_type == FLOWDIR_UP ||
+		integration.direction_type == FLOWDIR_BOTH) {
+		integration.actual_direction = FLOWDIR_UP;
+		compute_flowline(&region, &seed, &gradient_info, flowacc,
+				 &integration, &fl_map, fl_cats, fl_points,
+				 &cat, if_table, finfo, driver);
+	    }
+	    if (integration.direction_type == FLOWDIR_DOWN ||
+		integration.direction_type == FLOWDIR_BOTH) {
+		integration.actual_direction = FLOWDIR_DOWN;
+		compute_flowline(&region, &seed, &gradient_info, flowacc,
+				 &integration, &fl_map, fl_cats, fl_points,
+				 &cat, if_table, finfo, driver);
+	    }
+	    seed_count++;
+	}
+
+	Vect_destroy_line_struct(seed_points);
+	Vect_destroy_cats_struct(seed_cats);
+	Vect_close(&seed_Map);
+    }
+    if (flowacc_opt->answer || (!seed_opt->answer && flowlines_opt->answer)) {
+	/* compute flowlines from points on grid */
+	for (r = region.rows; r > 0; r--) {
+	    for (c = 0; c < region.cols; c++) {
+		for (d = 0; d < region.depths; d++) {
+		    seed.x =
+			region.west + c * region.ew_res + region.ew_res / 2;
+		    seed.y =
+			region.south + r * region.ns_res - region.ns_res / 2;
+		    seed.z =
+			region.bottom + d * region.tb_res + region.tb_res / 2;
+		    seed.flowline = FALSE;
+		    seed.flowaccum = FALSE;
+		    if (flowacc_opt->answer)
+			seed.flowaccum = TRUE;
+
+		    if (flowlines_opt->answer && !seed_opt->answer &&
+		       (c % skip[0] == 0) && (r % skip[1] == 0) && (d % skip[2] == 0))
+			seed.flowline = TRUE;
+
+		    if (seed.flowaccum || seed.flowline) {
+			G_percent(seed_count, n_seeds, 1);
+
+			if (integration.direction_type == FLOWDIR_UP ||
+			    integration.direction_type == FLOWDIR_BOTH) {
+			    integration.actual_direction = FLOWDIR_UP;
+			    compute_flowline(&region, &seed, &gradient_info,
+					     flowacc, &integration, &fl_map,
+					     fl_cats, fl_points, &cat,
+					     if_table, finfo, driver);
+			}
+			if (integration.direction_type == FLOWDIR_DOWN ||
+			    integration.direction_type == FLOWDIR_BOTH) {
+			    integration.actual_direction = FLOWDIR_DOWN;
+			    compute_flowline(&region, &seed, &gradient_info,
+					     flowacc, &integration, &fl_map,
+					     fl_cats, fl_points, &cat,
+					     if_table, finfo, driver);
+			}
+			seed_count++;
+		    }
+		}
+	    }
+	}
+    }
+    G_percent(1, 1, 1);
+    if (flowlines_opt->answer) {
+	if (if_table) {
+	    db_commit_transaction(driver);
+	    db_close_database_shutdown_driver(driver);
+	}
+	Vect_destroy_line_struct(fl_points);
+	Vect_destroy_cats_struct(fl_cats);
+	Vect_build(&fl_map);
+	Vect_close(&fl_map);
+    }
+
+    if (flowacc_opt->answer)
+	Rast3d_close(flowacc);
+
+
+    return EXIT_SUCCESS;
+}

raster3d/r3.flow/r3.flow.html

@@ -0,0 +1,96 @@
+<h2>DESCRIPTION</h2>
+
+Module <em>r3.flow</em> computes 3D flow lines and 3D flow accumulation.
+It accepts either three 3D raster maps representing the vector field or one 3D raster map.
+In case of one map, it computes on-the-fly gradient field.
+
+<h3>Flow lines</h3>
+
+Flow lines are computed either from points (seeds) provided in <b>seed_points</b> vector map,
+or if there are no seeds, it creates seeds in a regular grid in the center of voxels (3D raster cells).
+Parameter <b>skip</b> then controls the step between the regularly distributed seeds.
+If skip is not provided, r3.flow decides optimal skip for each dimension based on current 3D region
+as one tenth of the number of columns, rows, and depths.
+Flow lines can be computed in upstream direction (in the direction of gradient or vector field),
+in downstream direction or in both directions.
+
+<h3>Flow accumulation</h3>
+Flow accumulation is computed as the number of flow lines traversing each voxel.
+Since the flow lines are computed for each voxel, the flow accumulation computation
+can be more demanding.
+Parameter skip does not influence the flow accumulation computation, parameter direction does.
+
+<h3>Flow line integration</h3>
+Flow line integration can be influenced by several parameters.
+Option <b>step</b> controls the integration step and influences the precision and computational time.
+The unit of step can be defined either in terms of the size of the voxel (3D raster cell),
+length in map units, or as elapsed time. 
+Option <b>limit</b> specifies the maximum number of steps of each flow line.
+
+<h3>Attributes</h3>
+Without using flag <b>a</b>, no attribute table is created and each flow line
+is represented by one vector line with one category. With <b>a</b> flag, an attribute table is created
+and each category (record) represents one segment of a flowline, so that attributes
+specific for segments can be written. In case of <b>vector_field</b> input, only velocity is written,
+in case of <b>input</b> option, also values of the input 3D raster are written.
+Note that using <b>a</b> flag results in longer computation time, so consider increasing
+<b>step</b> and <b>max_error</b> parameter.
+
+<h2>NOTES</h2>
+r3.flow uses Runge-Kutta with adaptive step size
+(<a href="http://en.wikipedia.org/wiki/Cash-Karp_method">Cash-Karp method</a>).
+
+<h2>EXAMPLES</h2>
+First we create input data using
+<a href="http://grass.osgeo.org/grass71/manuals/r3.gwflow.html#example-1">example 1</a> from
+<a href="http://grass.osgeo.org/grass71/manuals/r3.gwflow.html">r3.gwflow manual page</a>:
+
+<div class="code"><pre>
+# set the region accordingly
+g.region res=25 res3=25 t=100 b=0 n=1000 s=0 w=0 e=1000 -p3
+
+# now create the input raster maps for a confined aquifer
+r3.mapcalc expression="phead = if(row() == 1 && depth() == 4, 50, 40)"
+r3.mapcalc expression="status = if(row() == 1 && depth() == 4, 2, 1)"
+r3.mapcalc expression="well = if(row() == 20 && col() == 20 && depth() == 2, -0.25, 0)"
+r3.mapcalc expression="hydcond = 0.00025"
+r3.mapcalc expression="syield = 0.0001"
+r.mapcalc  expression="recharge = 0.0"
+
+r3.gwflow solver=cg phead=phead status=status hc_x=hydcond hc_y=hydcond  \
+   hc_z=hydcond q=well s=syield r=recharge output=gwresult dt=8640000 vx=vx vy=vy vz=vz budget=budget
+</pre></div>
+
+Then we compute flow lines in both directions and downstream flowaccumulation.
+
+<div class="code"><pre>
+r3.flow vector_field=vx,vy,vz flowline=gw_flowlines skip=5,5,2 direction=both
+r3.flow vector_field=vx,vy,vz flowaccumulation=gw_flowacc
+</pre></div>
+
+We can visualize the result in 3D view:
+<p>
+<img src="r3flow_flowlines.png">
+
+<p>
+We can store velocity values (and values of the input 3D raster map if we use option <b>input</b>) for each segment of flow line
+in an attribute table.
+<div class="code"><pre>
+r3.flow -a vector_field=vx,vy,vz flowline=gw_flowlines skip=5,5,2 direction=both
+v.colors map=flowlines_color@user1 use=attr column=velocity color=bcyr
+</pre></div>
+
+Again, we visualize the result in 3D view and we check 'use color for thematic rendering' on 3D view vector page.
+<p>
+<img src="r3flow_flowlines_color.png">
+
+<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.flow.html">r.flow</a>,
+<a href="http://grass.osgeo.org/grass71/manuals/r3.gwflow.html">r3.gwflow</a>
+</em>
+
+<p><i>Last changed: $Date$</i>

raster3d/r3.flow/r3flow_structs.h

@@ -0,0 +1,40 @@
+#ifndef R3FLOW_STRUCTS_H
+#define R3FLOW_STRUCTS_H
+
+#include <grass/raster3d.h>
+
+struct Seed
+{
+    double x;
+    double y;
+    double z;
+    int flowline;
+    int flowaccum;
+};
+
+enum flowdir {FLOWDIR_UP, FLOWDIR_DOWN, FLOWDIR_BOTH};
+
+struct Integration
+{
+    enum flowdir direction_type;
+    enum flowdir actual_direction;
+    char *unit;
+    double step;
+    double cell_size;
+    int limit;
+    double max_error;
+    double max_step;
+    double min_step;
+};
+
+struct Gradient_info
+{
+    int compute_gradient;
+    RASTER3D_Map *velocity_maps[3];
+    RASTER3D_Map *scalar_map;
+    double neighbors_values[24];
+    int neighbors_pos[3];
+    int initialized;
+};
+
+#endif // R3FLOW_STRUCTS_H

raster3d/r3.flow/test_main.c

@@ -0,0 +1,105 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/raster3d.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+
+#include "r3flow_structs.h"
+#include "flowline.h"
+#include "interpolate.h"
+
+static void test_interpolation(RASTER3D_Region * region,
+			       RASTER3D_Map ** input_maps, double north,
+			       double east, double top)
+{
+    double interpolated[3];
+
+    if (interpolate_velocity(region, input_maps, north, east, top,
+			     &interpolated[0], &interpolated[1],
+			     &interpolated[2]) < 0) {
+	fprintf(stdout, "return=-1\n");
+    }
+    else
+	fprintf(stdout, "return=0\nvalues=%.10f,%.10f,%.10f\n",
+		interpolated[0], interpolated[1], interpolated[2]);
+
+}
+
+int main(int argc, char *argv[])
+{
+    int i;
+    struct GModule *module;
+    struct Option *test_opt, *coordinates_opt, *input_opt;
+    RASTER3D_Region region;
+    RASTER3D_Map *input_3drasters[3];
+    double coordinates[3];
+
+    G_gisinit(argv[0]);
+    module = G_define_module();
+    G_add_keyword(_("raster3d"));
+    G_add_keyword(_("voxel"));
+    G_add_keyword(_("flowline"));
+    module->description = _("Testing flow lines.");
+
+    test_opt = G_define_option();
+    test_opt->key = "test";
+    test_opt->required = YES;
+    test_opt->type = TYPE_STRING;
+    test_opt->options = "interpolation,gradient";
+    test_opt->description = "Select what is tested";
+
+    coordinates_opt = G_define_option();
+    coordinates_opt->key = "coordinates";
+    coordinates_opt->required = NO;
+    coordinates_opt->type = TYPE_DOUBLE;
+    coordinates_opt->multiple = YES;
+    coordinates_opt->description = "x,y,z coordinates";
+
+    input_opt = G_define_standard_option(G_OPT_R3_INPUTS);
+    input_opt->required = NO;
+
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    Rast3d_init_defaults();
+    Rast3d_get_window(&region);
+
+    if (strcmp(test_opt->answer, "interpolation") == 0) {
+
+
+	if (input_opt->answers) {
+	    for (i = 0; i < 3; i++) {
+		input_3drasters[i] =
+		    Rast3d_open_cell_old(input_opt->answers[i],
+					 G_find_raster3d(input_opt->
+							 answers[i], ""),
+					 &region, RASTER3D_TILE_SAME_AS_FILE,
+					 RASTER3D_USE_CACHE_DEFAULT);
+		if (input_3drasters[i] == NULL)
+		    Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
+				       input_opt->answers[i]);
+	    }
+	}
+	else
+	    G_fatal_error("No input map for interpolation test");
+	if (coordinates_opt->answers) {
+	    for (i = 0; i < 3; i++) {
+		if (coordinates_opt->answers[i]) {
+		    coordinates[i] = atof(coordinates_opt->answers[i]);
+		}
+		else
+		    G_fatal_error("Provide 3 coordinates");
+	    }
+	}
+	else
+	    G_fatal_error("No coordinates for interpolation test");
+	test_interpolation(&region, input_3drasters, coordinates[1],
+			   coordinates[0], coordinates[2]);
+    }
+
+    return EXIT_SUCCESS;
+}

raster3d/r3.flow/testsuite/data/flowline.ascii

@@ -0,0 +1,220 @@
+ORGANIZATION: 
+DIGIT DATE:   
+DIGIT NAME:   anna
+MAP NAME:     
+MAP DATE:     Thu Aug  7 14:07:09 2014
+MAP SCALE:    1
+OTHER INFO:   
+ZONE:         0
+MAP THRESH:   0.000000
+VERTI:
+L  29 1
+ 84.806094    35.197336    39.43322    
+ 84.806094    35.197336    39.43322    
+ 81.397053    37.087846    37.746534   
+ 77.990596    38.976923    36.061127   
+ 74.586931    40.864452    34.377101   
+ 71.186292    42.750303    32.694572   
+ 67.78894     44.634331    31.013669   
+ 64.395174    46.516369    29.334541   
+ 61.005337    48.39623     27.657357   
+ 57.619819    50.273695    25.982309   
+ 54.239073    52.148514    24.309623   
+ 50.863626    54.020394    22.639558   
+ 47.494096    55.888993    20.972421   
+ 44.131215    57.753905    19.308574   
+ 40.775859    59.614643    17.64845    
+ 37.429084    61.470623    15.992571   
+ 34.092185    63.321127    14.341579   
+ 30.766765    65.165264    12.696266   
+ 27.454848    67.001913    11.057634   
+ 24.159033    68.829633    9.426969    
+ 20.882729    70.646532    7.805957    
+ 17.63053     72.450064    6.196872    
+ 16.2203      73.232118    5.499135    
+ 15.524693    73.617873    5.15497     
+ 14.844693    73.994972    4.816822    
+ 13.957167    74.487156    4.359276    
+ 13.294472    74.854659    4.003445    
+ 12.641464    75.219149    3.639976    
+ 10.056633    76.803483    2.055642    
+ 1     1         
+L  11 1
+ 14.429159    56.869515    22.42339    
+ 14.429159    56.869515    22.42339    
+ 12.601355    59.515898    19.843526   
+ 10.782742    62.148975    17.276635   
+ 8.975302     64.765875    14.725514   
+ 7.181717     67.362715    12.193949   
+ 5.405736     69.934067    9.687231    
+ 5.028815     70.479791    9.155225    
+ 4.652154     71.016388    8.632116    
+ 2.867253     73.302256    6.403706    
+ 2.265763     73.985122    5.738005    
+ 1     2         
+L  3 1
+ 29.06094     78.060291    39.317079   
+ 29.06094     78.060291    39.317079   
+ 26.687827    78.782477    36.229394   
+ 1     3         
+L  15 1
+ 64.955248    50.766146    12.023014   
+ 64.955248    50.766146    12.023014   
+ 61.256917    52.57599     11.123057   
+ 57.564137    54.383119    10.22445    
+ 53.877539    56.187221    9.327348    
+ 50.197872    57.987932    8.431932    
+ 46.526023    59.784817    7.538419    
+ 42.863063    61.577353    6.647069    
+ 39.210293    63.364901    5.758198    
+ 38.018565    63.948095    5.468202    
+ 36.824201    64.532578    5.177563    
+ 36.071128    64.901108    4.994309    
+ 32.597809    66.600839    4.110713    
+ 29.181582    68.272631    3.158125    
+ 25.83767     69.909033    2.128626    
+ 1     4         
+L  25 1
+ 75.476779    18.36009     29.362122   
+ 75.476779    18.36009     29.362122   
+ 72.370137    20.932588    28.035642   
+ 69.265932    23.503067    26.710203   
+ 66.164365    26.071363    25.38589    
+ 63.065658    28.63729     24.062799   
+ 59.970067    31.200637    22.741038   
+ 56.877879    33.761166    21.42073    
+ 53.789424    36.318604    20.102016   
+ 50.705085    38.872634    18.785059   
+ 47.625301    41.422891    17.470047   
+ 44.550588    43.96895     16.157201   
+ 41.481554    46.510306    14.846779   
+ 38.41892     49.046362    13.53909    
+ 35.363553    51.576401    12.234504   
+ 32.316507    54.09955     10.93347    
+ 29.279077    56.614736    9.636543    
+ 26.252882    59.120619    8.344413    
+ 23.23998     61.615494    7.057958    
+ 20.243039    64.097152    5.778319    
+ 19.425572    64.774067    5.429276    
+ 18.800082    65.292013    5.162202    
+ 18.176085    65.808722    4.89522     
+ 15.400616    68.106988    3.619661    
+ 12.732901    70.316026    2.218227    
+ 1     5         
+L  27 1
+ 88.842317    31.231087    34.335553   
+ 88.842317    31.231087    34.335553   
+ 85.377082    33.216565    32.883043   
+ 81.91433     35.200622    31.431573   
+ 78.454253    37.183145    29.981225   
+ 74.99707     39.16401     28.532089   
+ 71.543022    41.143079    27.084269   
+ 68.092387    43.120192    25.637878   
+ 64.645478    45.095171    24.193049   
+ 61.202652    47.067809    22.749932   
+ 57.764323    49.037872    21.3087     
+ 54.330967    51.005085    19.869552   
+ 50.903141    52.969129    18.432722   
+ 47.4815      54.92963     16.998485   
+ 44.06682     56.886142    15.567166   
+ 40.660031    58.838132    14.139154   
+ 37.262263    60.784955    12.714924   
+ 33.874899    62.725815    11.295054   
+ 30.499664    64.659727    9.880269    
+ 27.138741    66.585438    8.471483    
+ 23.794948    68.501334    7.069877    
+ 20.472006    70.405283    5.677011    
+ 19.348268    71.049152    5.205978    
+ 18.653231    71.447388    4.914227    
+ 15.59501     73.199659    3.529718    
+ 12.680501    74.869588    2.003854    
+ 12.040711    75.238889    1.635283    
+ 1     6         
+L  12 1
+ 41.138221    64.704133    25.85159    
+ 41.138221    64.704133    25.85159    
+ 37.878495    66.138062    23.671884   
+ 34.62772     67.568054    21.498164   
+ 31.387297    68.993493    19.331366   
+ 28.158978    70.413607    17.172662   
+ 24.94499     71.827417    15.023541   
+ 21.748226    73.23365     12.885936   
+ 18.572528    74.630617    10.762419   
+ 17.773874    74.981939    10.228377   
+ 16.750722    75.441692    9.544218    
+ 13.85757     76.869017    7.609633    
+ 1     7         
+L  10 1
+ 14.157683    45.265562    10.630492   
+ 14.157683    45.265562    10.630492   
+ 12.464634    48.777068    9.249157    
+ 10.77904     52.273115    7.873903    
+ 9.102396     55.750597    6.505953    
+ 7.436691     59.205391    5.146927    
+ 7.088389     59.927794    4.861766    
+ 5.550659     63.11716     3.501194    
+ 5.214971     63.813401    3.177855    
+ 4.881477     64.503578    2.846669    
+ 1     8         
+L  30 1
+ 69.363152    2.929942     39.636635   
+ 69.363152    2.929942     39.636635   
+ 68.798791    3.539399     39.297875   
+ 68.116225    4.282763     38.888164   
+ 67.559167    4.894601     38.553789   
+ 67.016431    5.493706     38.22801    
+ 64.370003    8.415368     36.639487   
+ 61.725453    11.334957    35.05209    
+ 59.082925    14.252313    33.465907   
+ 56.442583    17.167257    31.881036   
+ 53.80461     20.079585    30.297588   
+ 51.169212    22.989069    28.715685   
+ 48.536626    25.89545     27.13547    
+ 45.90712     28.798431    25.557103   
+ 43.281002    31.697671    23.98077    
+ 40.65863     34.592775    22.406686   
+ 38.040419    37.483285    20.8351     
+ 35.42686     40.36866     19.266306   
+ 32.818531    43.248261    17.700651   
+ 30.216125    46.121322    16.138551   
+ 27.620482    48.986918    14.580511   
+ 25.032628    51.843914    13.027147   
+ 22.453842    54.6909      11.479225   
+ 19.885737    57.526093    9.937715    
+ 17.33039     60.347202    8.403862    
+ 14.790534    63.151208    6.879309    
+ 12.269864    65.934034    5.366271    
+ 11.73894     66.520174    5.047583    
+ 11.214399    67.099269    4.729052    
+ 8.923633     69.62828     3.205931    
+ 1     9         
+L  28 1
+ 95.880286    0.802102     28.572064   
+ 95.880286    0.802102     28.572064   
+ 93.148886    2.998018     27.663079   
+ 92.408842    3.603505     27.416799   
+ 91.636987    4.239939     27.159932   
+ 90.785863    4.947655     26.876686   
+ 90.140973    5.486583     26.662072   
+ 75.332663    17.862527    21.733998   
+ 72.17706     20.499801    20.683841   
+ 69.024022    23.13493     19.634538   
+ 65.873765    25.767736    18.586161   
+ 62.726529    28.398017    17.538788   
+ 59.582587    31.025545    16.492512   
+ 56.442251    33.65006     15.447436   
+ 53.305878    36.271261    14.403679   
+ 50.173883    38.888805    13.361379   
+ 47.046745    41.502289    12.320695   
+ 43.925029    44.111242    11.281816   
+ 40.809402    46.715106    10.244963   
+ 37.70066     49.313216    9.210401    
+ 34.599764    51.904768    8.17845     
+ 31.507889    54.488782    7.149502    
+ 28.426489    57.064041    6.124039    
+ 25.357395    59.629015    5.102672    
+ 24.603013    60.259484    4.85051     
+ 21.709165    62.677998    3.821816    
+ 18.883632    65.039418    2.703682    
+ 16.147183    67.326387    1.486815    
+ 1     10        

raster3d/r3.flow/testsuite/r3flow_test.py

@@ -0,0 +1,96 @@
+# -*- coding: utf-8 -*-
+"""
+Test of r3.flow
+
+@author Anna Petrasova
+"""
+import os
+from grass.gunittest import TestCase, test
+
+seeds = """
+84.80609404|35.19733594|39.43321996
+14.42915927|56.86951467|22.42338987
+29.06094033|78.06029074|39.31707858
+64.95524796|50.76614609|12.02301418
+75.47677891|18.36008965|29.362122
+88.84231714|31.23108675|34.33555293
+41.13822083|64.70413255|25.85158957
+14.15768282|45.26556161|10.63049231
+69.36315244|2.92994235|39.63663467
+95.88028571|0.80210167|28.57206445
+"""
+
+flowaccum = """
+n=480
+null_cells=0
+cells=480
+min=0
+max=91
+range=91
+mean=6.38333333333333
+mean_of_abs=6.38333333333333
+stddev=11.3685848821312
+variance=129.244722222222
+coeff_var=178.097935490306
+sum=3064
+"""
+
+class FlowlineTest(TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        """Use temporary region settings"""
+        cls.use_temp_region()
+        cls.runModule("g.region", res=10, res3=10, n=80, s=0, w=0, e=120, b=0, t=50)
+        cls.runModule("r3.mapcalc", expression="map_1 = 100")
+        cls.runModule("r3.mapcalc", expression="map_2 = -20")
+        cls.runModule("r3.mapcalc", expression="map_3 = 0.01")
+        cls.runModule("r3.mapcalc", expression="map_4 = col() + row() + depth()")
+        cls.runModule("r3.mapcalc", expression="map_5 = col() * col() + row() * row() + depth() * depth()")
+        cls.runModule('v.in.ascii', input='-', output='test_seeds', z=3, flags='zt',
+                      stdin=seeds)
+
+    @classmethod
+    def tearDownClass(cls):
+        """!Remove the temporary region"""
+        cls.del_temp_region()
+        cls.runModule('g.remove', rast3d=['map_1', 'map_2', 'map_3', 'map_4', 'map_5', 'test_flowaccum'])
+        cls.runModule('g.remove', vect=['test_flowline', 'test_seeds'])
+        os.remove('./data/flowline_tmp.ascii')
+
+    def test_interpolation(self):
+        self.assertModuleKeyValue('test.r3flow', test='interpolation',
+                                  coordinates=[100, 55, 11], input=['map_1', 'map_2', 'map_3'],
+                                  reference={'return': 0, 'values': [100, -20, 0.01]},
+                                  precision=1e-10, sep='=')
+        self.assertModuleKeyValue('test.r3flow', test='interpolation',
+                                  coordinates=[5, 5, 5], input=['map_1', 'map_2', 'map_3'],
+                                  reference={'return': 0, 'values': [100, -20, 0.01]},
+                                  precision=1e-10, sep='=')
+        self.assertModuleKeyValue('test.r3flow', test='interpolation',
+                                  coordinates=[10, 10, 60], input=['map_1', 'map_2', 'map_3'],
+                                  reference={'return': -1},
+                                  precision=1e-10, sep='=')
+        self.assertModuleKeyValue('test.r3flow', test='interpolation',
+                                  coordinates=[25, 69, 17], input=['map_4', 'map_4', 'map_4'],
+                                  reference={'return': 0, 'values': [7.8, 7.8, 7.8]},
+                                  precision=1e-10, sep='=')
+        self.assertModuleKeyValue('test.r3flow', test='interpolation',
+                                  coordinates=[81, 30, 25], input=['map_4', 'map_4', 'map_4'],
+                                  reference={'return': 0, 'values': [18.1, 18.1, 18.1]},
+                                  precision=1e-10, sep='=')
+
+    def test_flowlines(self):
+        self.assertModule('r3.flow', input='map_5', flowline='test_flowline',
+                          seed_points='test_seeds', flowaccumulation='test_flowaccum',
+                          direction='down')
+        self.runModule('v.out.ascii', input='test_flowline',
+                       format='standard', output='./data/flowline_tmp.ascii',
+                       dp=6)
+        self.assertVectorAsciiEqualsVectorAscii(actual='./data/flowline_tmp.ascii',
+                                                reference='./data/flowline.ascii')
+        self.assertRaster3dFitsUnivar('test_flowaccum', reference=flowaccum, precision=1e-6)
+
+
+if __name__ == '__main__':
+    test()

raster3d/r3.flow/voxel_traversal.c

@@ -0,0 +1,96 @@
+#include <math.h>
+#include <grass/raster3d.h>
+
+#include "voxel_traversal.h"
+
+void traverse(RASTER3D_Region * region, double *start, double *end,
+	      int *coordinates, int *size, int *coor_count)
+{
+    double dx, dy, dz;
+    int step_x, step_y, step_z;
+    int x, y, z;
+    int x_end, y_end, z_end;
+    double t_delta_x, t_delta_y, t_delta_z;
+    double t_max_x, t_max_y, t_max_z;
+    double xtmp, ytmp, ztmp;
+    int count;
+
+    /* initialize */
+    dx = end[0] - start[0];
+    dy = end[1] - start[1];
+    dz = end[2] - start[2];
+
+    step_x = start[0] < end[0] ? 1 : -1;
+    step_y = start[1] < end[1] ? 1 : -1;
+    step_z = start[2] < end[2] ? 1 : -1;
+
+    x = (int)floor((start[0] - region->west) / region->ew_res);
+    y = (int)floor((start[1] - region->south) / region->ns_res);
+    z = (int)floor((start[2] - region->bottom) / region->tb_res);
+    x_end = (int)floor((end[0] - region->west) / region->ew_res);
+    y_end = (int)floor((end[1] - region->south) / region->ns_res);
+    z_end = (int)floor((end[2] - region->bottom) / region->tb_res);
+
+    t_delta_x = fabs(region->ew_res / (dx != 0 ? dx : 1e6));
+    t_delta_y = fabs(region->ns_res / (dy != 0 ? dy : 1e6));
+    t_delta_z = fabs(region->tb_res / (dz != 0 ? dz : 1e6));
+
+    xtmp = (start[0] - region->west) / region->ew_res;
+    ytmp = (start[1] - region->south) / region->ns_res;
+    ztmp = (start[2] - region->bottom) / region->tb_res;
+
+    if (step_x > 0)
+	t_max_x = t_delta_x * (1.0 - (xtmp - floor(xtmp)));
+    else
+	t_max_x = t_delta_x * (xtmp - floor(xtmp));
+    if (step_y > 0)
+	t_max_y = t_delta_y * (1.0 - (ytmp - floor(ytmp)));
+    else
+	t_max_y = t_delta_y * (ytmp - floor(ytmp));
+    if (step_z > 0)
+	t_max_z = t_delta_z * (1.0 - (ztmp - floor(ztmp)));
+    else
+	t_max_z = t_delta_z * (ztmp - floor(ztmp));
+
+    count = 0;
+    while (TRUE) {
+	if (t_max_x < t_max_y) {
+	    if (t_max_x < t_max_z) {
+		t_max_x = t_max_x + t_delta_x;
+		x = x + step_x;
+	    }
+	    else {
+		t_max_z = t_max_z + t_delta_z;
+		z = z + step_z;
+	    }
+	}
+	else {
+	    if (t_max_y < t_max_z) {
+		t_max_y = t_max_y + t_delta_y;
+		y = y + step_y;
+	    }
+	    else {
+		t_max_z = t_max_z + t_delta_z;
+		z = z + step_z;
+	    }
+	}
+	if ((x == x_end && y == y_end && z == z_end) ||
+	    /* just to make sure it breaks */
+	    (step_x * (x - x_end) > 0 || step_y * (y - y_end) > 0 ||
+	     step_z * (z - z_end) > 0))
+
+	    break;
+
+	coordinates[count * 3 + 0] = x;
+	coordinates[count * 3 + 1] = region->rows - y - 1;
+	coordinates[count * 3 + 2] = z;
+	count++;
+
+	/* reallocation for cases when the steps would be too big */
+	if (*size <= count) {
+	    *size = 2 * (*size);
+	    coordinates = G_realloc(coordinates, (*size) * 3 * sizeof(int));
+	}
+    }
+    *coor_count = count;
+}

raster3d/r3.flow/voxel_traversal.h

@@ -0,0 +1,9 @@
+#ifndef VOXEL_TRAVERSAL_H
+#define VOXEL_TRAVERSAL_H
+
+#include <grass/raster3d.h>
+
+void traverse(RASTER3D_Region *region, double *start, double *end,
+              int *coordinates, int *size, int *coor_count);
+
+#endif // VOXEL_TRAVERSAL_H

raster3d/r3.gradient/Makefile

@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM=r3.gradient
+
+LIBES = $(RASTER3DLIB) $(RASTERLIB) $(GISLIB)
+DEPENDENCIES = $(RASTER3DDEP) $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

raster3d/r3.gradient/main.c

@@ -0,0 +1,302 @@
+/****************************************************************************
+ * 
+ * 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>
+
+#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;
+    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->key_desc = "grad_x,grad_y,grad_z";
+    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");
+
+    /* disabled - was there for openMP
+    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 = 1;
+    /* disabled - was there for openMP
+    N = atoi(process_opt->answer);
+#if defined(_OPENMP)
+    omp_set_num_threads(N);
+#endif
+    */
+
+    Rast3d_init_defaults();
+    Rast3d_get_window(&region);
+
+    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, ""),
+				 &region, 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,
+					  &region, 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 */
+		    /* disabled openMP #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);
+}

raster3d/r3.gradient/r3.gradient.html

@@ -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>

raster3d/r3.gradient/r3gradient_structs.h

@@ -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

raster3d/r3.gradient/testsuite/data/test_map_1

@@ -0,0 +1,61 @@
+version: grass7
+order: nsbt
+north: 100.000000
+south: 0.000000
+east: 120.000000
+west: 0.000000
+top: 50.000000
+bottom: 0.000000
+rows: 10
+cols: 12
+levels: 5
+1.73205081 2.44948974 3.31662479 4.24264069 5.19615242 6.16441400 7.14142843 8.12403840 9.11043358 10.09950494 11.09053651 12.08304597 
+2.44948974 3.00000000 3.74165739 4.58257569 5.47722558 6.40312424 7.34846923 8.30662386 9.27361850 10.24695077 11.22497216 12.20655562 
+3.31662479 3.74165739 4.35889894 5.09901951 5.91607978 6.78232998 7.68114575 8.60232527 9.53939201 10.48808848 11.44552314 12.40967365 
+4.24264069 4.58257569 5.09901951 5.74456265 6.48074070 7.28010989 8.12403840 9.00000000 9.89949494 10.81665383 11.74734012 12.68857754 
+5.19615242 5.47722558 5.91607978 6.48074070 7.14142843 7.87400787 8.66025404 9.48683298 10.34408043 11.22497216 12.12435565 13.03840481 
+6.16441400 6.40312424 6.78232998 7.28010989 7.87400787 8.54400375 9.27361850 10.04987562 10.86278049 11.70469991 12.56980509 13.45362405 
+7.14142843 7.34846923 7.68114575 8.12403840 8.66025404 9.27361850 9.94987437 10.67707825 11.44552314 12.24744871 13.07669683 13.92838828 
+8.12403840 8.30662386 8.60232527 9.00000000 9.48683298 10.04987562 10.67707825 11.35781669 12.08304597 12.84523258 13.63818170 14.45683229 
+9.11043358 9.27361850 9.53939201 9.89949494 10.34408043 10.86278049 11.44552314 12.08304597 12.76714533 13.49073756 14.24780685 15.03329638 
+10.09950494 10.24695077 10.48808848 10.81665383 11.22497216 11.70469991 12.24744871 12.84523258 13.49073756 14.17744688 14.89966443 15.65247584 
+2.44948974 3.00000000 3.74165739 4.58257569 5.47722558 6.40312424 7.34846923 8.30662386 9.27361850 10.24695077 11.22497216 12.20655562 
+3.00000000 3.46410162 4.12310563 4.89897949 5.74456265 6.63324958 7.54983444 8.48528137 9.43398113 10.39230485 11.35781669 12.32882801 
+3.74165739 4.12310563 4.69041576 5.38516481 6.16441400 7.00000000 7.87400787 8.77496439 9.69535971 10.63014581 11.57583690 12.52996409 
+4.58257569 4.89897949 5.38516481 6.00000000 6.70820393 7.48331477 8.30662386 9.16515139 10.04987562 10.95445115 11.87434209 12.80624847 
+5.47722558 5.74456265 6.16441400 6.70820393 7.34846923 8.06225775 8.83176087 9.64365076 10.48808848 11.35781669 12.24744871 13.15294644 
+6.40312424 6.63324958 7.00000000 7.48331477 8.06225775 8.71779789 9.43398113 10.19803903 11.00000000 11.83215957 12.68857754 13.56465997 
+7.34846923 7.54983444 7.87400787 8.30662386 8.83176087 9.43398113 10.09950494 10.81665383 11.57583690 12.36931688 13.19090596 14.03566885 
+8.30662386 8.48528137 8.77496439 9.16515139 9.64365076 10.19803903 10.81665383 11.48912529 12.20655562 12.96148140 13.74772708 14.56021978 
+9.27361850 9.43398113 9.69535971 10.04987562 10.48808848 11.00000000 11.57583690 12.20655562 12.88409873 13.60147051 14.35270009 15.13274595 
+10.24695077 10.39230485 10.63014581 10.95445115 11.35781669 11.83215957 12.36931688 12.96148140 13.60147051 14.28285686 15.00000000 15.74801575 
+3.31662479 3.74165739 4.35889894 5.09901951 5.91607978 6.78232998 7.68114575 8.60232527 9.53939201 10.48808848 11.44552314 12.40967365 
+3.74165739 4.12310563 4.69041576 5.38516481 6.16441400 7.00000000 7.87400787 8.77496439 9.69535971 10.63014581 11.57583690 12.52996409 
+4.35889894 4.69041576 5.19615242 5.83095189 6.55743852 7.34846923 8.18535277 9.05538514 9.94987437 10.86278049 11.78982612 12.72792206 
+5.09901951 5.38516481 5.83095189 6.40312424 7.07106781 7.81024968 8.60232527 9.43398113 10.29563014 11.18033989 12.08304597 13.00000000 
+5.91607978 6.16441400 6.55743852 7.07106781 7.68114575 8.36660027 9.11043358 9.89949494 10.72380529 11.57583690 12.44989960 13.34166406 
+6.78232998 7.00000000 7.34846923 7.81024968 8.36660027 9.00000000 9.69535971 10.44030651 11.22497216 12.04159458 12.88409873 13.74772708 
+7.68114575 7.87400787 8.18535277 8.60232527 9.11043358 9.69535971 10.34408043 11.04536102 11.78982612 12.56980509 13.37908816 14.21267040 
+8.60232527 8.77496439 9.05538514 9.43398113 9.89949494 10.44030651 11.04536102 11.70469991 12.40967365 13.15294644 13.92838828 14.73091986 
+9.53939201 9.69535971 9.94987437 10.29563014 10.72380529 11.22497216 11.78982612 12.40967365 13.07669683 13.78404875 14.52583905 15.29705854 
+10.48808848 10.63014581 10.86278049 11.18033989 11.57583690 12.04159458 12.56980509 13.15294644 13.78404875 14.45683229 15.16575089 15.90597372 
+4.24264069 4.58257569 5.09901951 5.74456265 6.48074070 7.28010989 8.12403840 9.00000000 9.89949494 10.81665383 11.74734012 12.68857754 
+4.58257569 4.89897949 5.38516481 6.00000000 6.70820393 7.48331477 8.30662386 9.16515139 10.04987562 10.95445115 11.87434209 12.80624847 
+5.09901951 5.38516481 5.83095189 6.40312424 7.07106781 7.81024968 8.60232527 9.43398113 10.29563014 11.18033989 12.08304597 13.00000000 
+5.74456265 6.00000000 6.40312424 6.92820323 7.54983444 8.24621125 9.00000000 9.79795897 10.63014581 11.48912529 12.36931688 13.26649916 
+6.48074070 6.70820393 7.07106781 7.54983444 8.12403840 8.77496439 9.48683298 10.24695077 11.04536102 11.87434209 12.72792206 13.60147051 
+7.28010989 7.48331477 7.81024968 8.24621125 8.77496439 9.38083152 10.04987562 10.77032961 11.53256259 12.32882801 13.15294644 14.00000000 
+8.12403840 8.30662386 8.60232527 9.00000000 9.48683298 10.04987562 10.67707825 11.35781669 12.08304597 12.84523258 13.63818170 14.45683229 
+9.00000000 9.16515139 9.43398113 9.79795897 10.24695077 10.77032961 11.35781669 12.00000000 12.68857754 13.41640786 14.17744688 14.96662955 
+9.89949494 10.04987562 10.29563014 10.63014581 11.04536102 11.53256259 12.08304597 12.68857754 13.34166406 14.03566885 14.76482306 15.52417470 
+10.81665383 10.95445115 11.18033989 11.48912529 11.87434209 12.32882801 12.84523258 13.41640786 14.03566885 14.69693846 15.39480432 16.12451550 
+5.19615242 5.47722558 5.91607978 6.48074070 7.14142843 7.87400787 8.66025404 9.48683298 10.34408043 11.22497216 12.12435565 13.03840481 
+5.47722558 5.74456265 6.16441400 6.70820393 7.34846923 8.06225775 8.83176087 9.64365076 10.48808848 11.35781669 12.24744871 13.15294644 
+5.91607978 6.16441400 6.55743852 7.07106781 7.68114575 8.36660027 9.11043358 9.89949494 10.72380529 11.57583690 12.44989960 13.34166406 
+6.48074070 6.70820393 7.07106781 7.54983444 8.12403840 8.77496439 9.48683298 10.24695077 11.04536102 11.87434209 12.72792206 13.60147051 
+7.14142843 7.34846923 7.68114575 8.12403840 8.66025404 9.27361850 9.94987437 10.67707825 11.44552314 12.24744871 13.07669683 13.92838828 
+7.87400787 8.06225775 8.36660027 8.77496439 9.27361850 9.84885780 10.48808848 11.18033989 11.91637529 12.68857754 13.49073756 14.31782106 
+8.66025404 8.83176087 9.11043358 9.48683298 9.94987437 10.48808848 11.09053651 11.74734012 12.44989960 13.19090596 13.96424004 14.76482306 
+9.48683298 9.64365076 9.89949494 10.24695077 10.67707825 11.18033989 11.74734012 12.36931688 13.03840481 13.74772708 14.49137675 15.26433752 
+10.34408043 10.48808848 10.72380529 11.04536102 11.44552314 11.91637529 12.44989960 13.03840481 13.67479433 14.35270009 15.06651917 15.81138830 
+11.22497216 11.35781669 11.57583690 11.87434209 12.24744871 12.68857754 13.19090596 13.74772708 14.35270009 15.00000000 15.68438714 16.40121947 

raster3d/r3.gradient/testsuite/data/test_map_2

@@ -0,0 +1,36 @@
+version: grass7
+order: nsbt
+north: 5.000000
+south: 0.000000
+east: 5.000000
+west: 0.000000
+top: 5.000000
+bottom: 0.000000
+rows: 5
+cols: 5
+levels: 5
+3.00000000 6.00000000 11.00000000 18.00000000 27.00000000 
+6.00000000 9.00000000 14.00000000 * 30.00000000 
+11.00000000 14.00000000 19.00000000 26.00000000 35.00000000 
+18.00000000 * 26.00000000 33.00000000 42.00000000 
+27.00000000 30.00000000 35.00000000 42.00000000 51.00000000 
+6.00000000 9.00000000 14.00000000 * 30.00000000 
+9.00000000 12.00000000 17.00000000 24.00000000 33.00000000 
+14.00000000 17.00000000 22.00000000 29.00000000 * 
+* 24.00000000 29.00000000 36.00000000 45.00000000 
+30.00000000 33.00000000 * 45.00000000 54.00000000 
+11.00000000 14.00000000 19.00000000 26.00000000 35.00000000 
+14.00000000 17.00000000 22.00000000 29.00000000 * 
+19.00000000 22.00000000 27.00000000 * 43.00000000 
+26.00000000 29.00000000 * 41.00000000 50.00000000 
+35.00000000 * 43.00000000 50.00000000 * 
+18.00000000 * 26.00000000 33.00000000 42.00000000 
+* 24.00000000 29.00000000 36.00000000 45.00000000 
+26.00000000 29.00000000 * 41.00000000 50.00000000 
+33.00000000 36.00000000 41.00000000 48.00000000 57.00000000 
+42.00000000 45.00000000 50.00000000 57.00000000 66.00000000 
+27.00000000 30.00000000 35.00000000 42.00000000 51.00000000 
+30.00000000 33.00000000 * 45.00000000 54.00000000 
+35.00000000 * 43.00000000 50.00000000 * 
+42.00000000 45.00000000 50.00000000 57.00000000 66.00000000 
+51.00000000 54.00000000 * 66.00000000 75.00000000 

raster3d/r3.gradient/testsuite/r3gradient_test.py

@@ -0,0 +1,157 @@
+# -*- coding: utf-8 -*-
+"""
+Test of r3.gradient
+
+@author Anna Petrasova
+"""
+from grass.gunittest import TestCase, test
+
+
+r3univar_test_grad_x = """
+n=600
+null_cells=0
+cells=600
+min=0.00902566899999995
+max=0.0993248405000001
+range=0.0902991715000001
+mean=0.0641879624599999
+mean_of_abs=0.0641879624599999
+stddev=0.0243482677445681
+variance=0.000592838142161176
+coeff_var=37.9327631091908
+sum=38.512777476
+"""
+
+
+r3univar_test_grad_y = """
+n=600
+null_cells=0
+cells=600
+min=-0.0990409449999998
+max=-0.00774536350000012
+range=0.0912955814999997
+mean=-0.0563959154616667
+mean_of_abs=0.0563959154616667
+stddev=0.0244377519801364
+variance=0.000597203721842658
+coeff_var=-43.3324856597942
+sum=-33.837549277"""
+
+
+r3univar_test_grad_z = """
+n=600
+null_cells=0
+cells=600
+min=0.00643308800000026
+max=0.0967259644999999
+range=0.0902928764999997
+mean=0.0336457494116667
+mean_of_abs=0.0336457494116667
+stddev=0.0186882020765464
+variance=0.000349248896853835
+coeff_var=55.5440208743464
+sum=20.187449647
+"""
+
+r3univar_test_nulls_grad_x = """
+n=107
+null_cells=18
+cells=125
+min=0
+max=10
+range=10
+mean=3.70093457943925
+mean_of_abs=3.70093457943925
+stddev=3.6357902977452
+variance=13.2189710891781
+coeff_var=98.2397883481656
+sum=396
+"""
+
+r3univar_test_nulls_grad_y = """
+n=107
+null_cells=18
+cells=125
+min=-10
+max=0
+range=10
+mean=-3.70093457943925
+mean_of_abs=3.70093457943925
+stddev=3.6357902977452
+variance=13.2189710891781
+coeff_var=-98.2397883481656
+sum=-396
+"""
+
+r3univar_test_nulls_grad_z = """
+n=107
+null_cells=18
+cells=125
+min=0
+max=10
+range=10
+mean=3.70093457943925
+mean_of_abs=3.70093457943925
+stddev=3.6357902977452
+variance=13.2189710891781
+coeff_var=98.2397883481656
+sum=396
+"""
+
+
+class GradientTest(TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        """Use temporary region settings"""
+        cls.use_temp_region()
+        cls.runModule('g.region', res3=10, n=100, s=0, w=0, e=120, b=0, t=50)
+        cls.runModule('r3.in.ascii', input='data/test_map_1', output='test_map_1_ref')
+        cls.runModule('g.region', res3=1, n=5, s=0, w=0, e=5, b=0, t=5)
+        cls.runModule('r3.in.ascii', input='data/test_map_2', output='test_map_2_ref')
+
+    @classmethod
+    def tearDownClass(cls):
+        """!Remove the temporary region"""
+        cls.del_temp_region()
+        cls.runModule('g.remove', rast3d=['test_map_1_ref', 'test_map_2_ref', 'test_grad_x',
+                                          'test_grad_y', 'test_grad_z', 'test_null_grad_x',
+                                          'test_null_grad_y', 'test_null_grad_z'])
+
+    def test_gradient_runs(self):
+        self.runModule('g.region', res3=10, n=100, s=0, w=0, e=120, b=0, t=50)
+        self.assertModuleFail('r3.gradient', input='test_map_1_ref',
+                              output=['test_grad_x', 'test_grad_y'], overwrite=True)
+        self.assertModule('r3.gradient', input='test_map_1_ref',
+                          output=['test_grad_x', 'test_grad_y', 'test_grad_z'], overwrite=True)
+
+    def test_gradient(self):
+        self.runModule('g.region', res3=10, n=100, s=0, w=0, e=120, b=0, t=50)
+        self.runModule('r3.gradient', input='test_map_1_ref',
+                       output=['test_grad_x', 'test_grad_y', 'test_grad_z'], overwrite=True)
+        self.assertRaster3dFitsUnivar(raster='test_grad_x', reference=r3univar_test_grad_x, precision=1e-8)
+        self.assertRaster3dFitsUnivar(raster='test_grad_y', reference=r3univar_test_grad_y, precision=1e-8)
+        self.assertRaster3dFitsUnivar(raster='test_grad_z', reference=r3univar_test_grad_z, precision=1e-8)
+
+    def test_gradient_block(self):
+        self.runModule('g.region', res3=10, n=100, s=0, w=0, e=120, b=0, t=50)
+        self.assertModule('r3.gradient', input='test_map_1_ref', block_size=[200, 2, 50],
+                          output=['test_grad_x', 'test_grad_y', 'test_grad_z'], overwrite=True)
+        self.assertRaster3dFitsUnivar(raster='test_grad_x', reference=r3univar_test_grad_x, precision=1e-8)
+        self.assertRaster3dFitsUnivar(raster='test_grad_y', reference=r3univar_test_grad_y, precision=1e-8)
+        self.assertRaster3dFitsUnivar(raster='test_grad_z', reference=r3univar_test_grad_z, precision=1e-8)
+
+    def test_gradient_nulls(self):
+        self.runModule('g.region', res3=1, n=5, s=0, w=0, e=5, b=0, t=5)
+        self.assertModule('r3.gradient', input='test_map_2_ref', block_size=[200, 2, 50],
+                          output=['test_null_grad_x', 'test_null_grad_y', 'test_null_grad_z'])
+        self.assertRaster3dFitsUnivar(raster='test_null_grad_x',
+                                      reference=r3univar_test_nulls_grad_x, precision=1e-8)
+        self.assertRaster3dFitsUnivar(raster='test_null_grad_y',
+                                      reference=r3univar_test_nulls_grad_y, precision=1e-8)
+        self.assertRaster3dFitsUnivar(raster='test_null_grad_z',
+                                      reference=r3univar_test_nulls_grad_z, precision=1e-8)
+
+
+if __name__ == '__main__':
+    test()