[GRASS-SVN] r61235 - in grass-addons/grass7/raster3d: . r3.flow
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jul 10 11:51:58 PDT 2014
Author: annakrat
Date: 2014-07-10 11:51:58 -0700 (Thu, 10 Jul 2014)
New Revision: 61235
Added:
grass-addons/grass7/raster3d/Makefile
grass-addons/grass7/raster3d/r3.flow/
grass-addons/grass7/raster3d/r3.flow/Makefile
grass-addons/grass7/raster3d/r3.flow/flowline.c
grass-addons/grass7/raster3d/r3.flow/flowline.h
grass-addons/grass7/raster3d/r3.flow/integrate.c
grass-addons/grass7/raster3d/r3.flow/integrate.h
grass-addons/grass7/raster3d/r3.flow/interpolate.c
grass-addons/grass7/raster3d/r3.flow/interpolate.h
grass-addons/grass7/raster3d/r3.flow/main.c
grass-addons/grass7/raster3d/r3.flow/r3.flow.html
grass-addons/grass7/raster3d/r3.flow/r3flow_structs.h
Log:
r3.flow: computing 3D flowlines, not all functionality ported yet from Python version in sandbox
Added: grass-addons/grass7/raster3d/Makefile
===================================================================
--- grass-addons/grass7/raster3d/Makefile (rev 0)
+++ grass-addons/grass7/raster3d/Makefile 2014-07-10 18:51:58 UTC (rev 61235)
@@ -0,0 +1,9 @@
+MODULE_TOPDIR = ..
+
+SUBDIRS = \
+ r3.flow
+
+include $(MODULE_TOPDIR)/include/Make/Dir.make
+
+default: subdirs
+
Property changes on: grass-addons/grass7/raster3d/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster3d/r3.flow/Makefile
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/Makefile (rev 0)
+++ grass-addons/grass7/raster3d/r3.flow/Makefile 2014-07-10 18:51:58 UTC (rev 61235)
@@ -0,0 +1,12 @@
+MODULE_TOPDIR = ../..
+
+PGM=r3.flow
+
+LIBES = $(RASTER3DLIB) $(RASTERLIB) $(GISLIB) $(VECTORLIB)
+DEPENDENCIES = $(RASTER3DDEP) $(RASTERDEP) $(GISDEP) $(VECTORDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass-addons/grass7/raster3d/r3.flow/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster3d/r3.flow/flowline.c
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/flowline.c (rev 0)
+++ grass-addons/grass7/raster3d/r3.flow/flowline.c 2014-07-10 18:51:58 UTC (rev 61235)
@@ -0,0 +1,97 @@
+/*!
+ \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 "r3flow_structs.h"
+#include "integrate.h"
+#include "flowline.h"
+
+/*!
+ \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 cat category of the newly created flow line
+*/
+void compute_flowline(RASTER3D_Region * region, const struct Seed *seed,
+ RASTER3D_Map ** velocity_field,
+ struct Integration *integration,
+ struct Map_info *flowline_vec, struct line_cats *cats,
+ struct line_pnts *points, const int cat)
+{
+ int i, 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;
+
+ point[0] = seed->x;
+ point[1] = seed->y;
+ point[2] = seed->z;
+
+ if (seed->flowline) {
+ /* append first point */
+ Vect_append_point(points, seed->x, seed->y, seed->z);
+ }
+ count = 1;
+ while (count <= integration->limit) {
+ if (get_velocity(region, velocity_field, 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", MIN_STEP, velocity_norm,
+ integration->cell_size);
+ max_step = get_time_step("cell", MAX_STEP, velocity_norm,
+ integration->cell_size);
+ delta_t *= integration->direction;
+ if (rk45_integrate_next
+ (region, velocity_field, point, new_point,
+ &delta_t, min_step, max_step) < 0)
+ break;
+
+ if (seed->flowline)
+ Vect_append_point(points, new_point[0], new_point[1],
+ new_point[2]);
+
+ for (i = 0; i < 3; i++) {
+ point[i] = new_point[i];
+ }
+ count++;
+
+ }
+ if (seed->flowline && points->n_points > 1) {
+ Vect_cat_set(cats, 1, 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);
+}
Property changes on: grass-addons/grass7/raster3d/r3.flow/flowline.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster3d/r3.flow/flowline.h
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/flowline.h (rev 0)
+++ grass-addons/grass7/raster3d/r3.flow/flowline.h 2014-07-10 18:51:58 UTC (rev 61235)
@@ -0,0 +1,14 @@
+#ifndef FLOWLINE_H
+#define FLOWLINE_H
+
+#include <grass/raster3d.h>
+#include <grass/vector.h>
+
+static const double VELOCITY_EPSILON = 1e-8;
+
+void compute_flowline(RASTER3D_Region * region, const struct Seed *seed,
+ RASTER3D_Map ** velocity_field,
+ struct Integration *integration,
+ struct Map_info *flowline_vec, struct line_cats *cats,
+ struct line_pnts *points, const int cat);
+#endif // FLOWLINE_H
Property changes on: grass-addons/grass7/raster3d/r3.flow/flowline.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster3d/r3.flow/integrate.c
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/integrate.c (rev 0)
+++ grass-addons/grass7/raster3d/r3.flow/integrate.c 2014-07-10 18:51:58 UTC (rev 61235)
@@ -0,0 +1,216 @@
+/*!
+ \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"
+
+/*!
+ \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, RASTER3D_Map ** velocity_field,
+ const double x, const double y, const double z,
+ double *vel_x, double *vel_y, double *vel_z)
+{
+ return interpolate_velocity(region, velocity_field, 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, RASTER3D_Map ** velocity_field,
+ const double *point, double *next_point,
+ const double delta_t, 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;
+
+ if (get_velocity(region, velocity_field, 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, velocity_field, 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;
+ }
+
+ /* 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;
+ }
+ 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,
+ RASTER3D_Map ** velocity_field, const double *point,
+ double *next_point, double *delta_t,
+ const double min_step, const double max_step)
+{
+ 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, velocity_field, point, next_point, *delta_t, 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, velocity_field, point, next_point,
+ *delta_t, error) < 0)
+ return -1;
+ break;
+ }
+ }
+ return 0;
+}
Property changes on: grass-addons/grass7/raster3d/r3.flow/integrate.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster3d/r3.flow/integrate.h
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/integrate.h (rev 0)
+++ grass-addons/grass7/raster3d/r3.flow/integrate.h 2014-07-10 18:51:58 UTC (rev 61235)
@@ -0,0 +1,36 @@
+#ifndef INTEGRATE_H
+#define INTEGRATE_H
+
+#include <grass/raster3d.h>
+
+static const double MAX_ERROR = 1.0e-6;
+static const double MIN_STEP = 0.01;
+static const double MAX_STEP = 1;
+
+/* 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, RASTER3D_Map ** velocity_field,
+ 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,
+ RASTER3D_Map ** velocity_field, const double *point,
+ double *next_point, double *delta_t,
+ const double min_step, const double max_step);
+
+#endif // INTEGRATE_H
Property changes on: grass-addons/grass7/raster3d/r3.flow/integrate.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster3d/r3.flow/interpolate.c
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/interpolate.c (rev 0)
+++ grass-addons/grass7/raster3d/r3.flow/interpolate.c 2014-07-10 18:51:58 UTC (rev 61235)
@@ -0,0 +1,181 @@
+/*!
+ \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/raster3d.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;
+}
Property changes on: grass-addons/grass7/raster3d/r3.flow/interpolate.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster3d/r3.flow/interpolate.h
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/interpolate.h (rev 0)
+++ grass-addons/grass7/raster3d/r3.flow/interpolate.h 2014-07-10 18:51:58 UTC (rev 61235)
@@ -0,0 +1,10 @@
+#ifndef INTERPOLATE_H
+#define INTERPOLATE_H
+
+#include <grass/raster3d.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);
+#endif // INTERPOLATE_H
Property changes on: grass-addons/grass7/raster3d/r3.flow/interpolate.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster3d/r3.flow/main.c
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/main.c (rev 0)
+++ grass-addons/grass7/raster3d/r3.flow/main.c 2014-07-10 18:51:58 UTC (rev 61235)
@@ -0,0 +1,316 @@
+
+/****************************************************************************
+ *
+ * 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/glocale.h>
+
+#include "r3flow_structs.h"
+#include "flowline.h"
+
+
+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_velocity_maps(struct Option *vector_opt,
+ RASTER3D_Map ** velocity_field,
+ RASTER3D_Region * region)
+{
+ int i;
+
+ for (i = 0; i < 3; i++) {
+ velocity_field[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 (!velocity_field[i])
+ Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
+ vector_opt->answers[i]);
+ }
+}
+
+int main(int argc, char *argv[])
+{
+ struct Option *vector_opt, *seed_opt, *flowlines_opt,
+ *unit_opt, *step_opt, *limit_opt, *skip_opt;
+ struct Flag *up_flag;
+ struct GModule *module;
+ RASTER3D_Region region;
+ RASTER3D_Map *velocity_field[3];
+ struct Integration integration;
+ struct Seed seed;
+ 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 */
+ int cat; /* cat of flowlines */
+ 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 = _("Compute flow lines.");
+
+
+ vector_opt = G_define_standard_option(G_OPT_R3_INPUTS);
+ vector_opt->key = "vector_field";
+ vector_opt->required = YES;
+ vector_opt->description = _("Names of three 3D raster maps describing "
+ "x, y, z components of vector field");
+
+ 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");
+
+ flowlines_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+ flowlines_opt->key = "flowline";
+ flowlines_opt->required = YES;
+ flowlines_opt->description = _("Name for vector map of flow lines");
+
+ 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");
+
+ 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");
+
+ up_flag = G_define_flag();
+ up_flag->key = 'u';
+ up_flag->description = _("Compute upstream flowlines "
+ "instead of default downstream flowlines");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ check_vector_input_maps(vector_opt, seed_opt);
+
+ Rast3d_init_defaults();
+ Rast3d_get_window(®ion);
+
+ /* 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.limit = atof(limit_opt->answer);
+ integration.direction = up_flag->answer ? 1 : -1;
+
+
+ /* 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 */
+ load_input_velocity_maps(vector_opt, velocity_field, ®ion);
+
+ /* open new vector map of flowlines */
+ 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);
+
+ 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);
+ }
+ 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;
+ 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);
+ cat = seed_count + 1;
+ compute_flowline(®ion, &seed, velocity_field, &integration,
+ &fl_map, fl_cats, fl_points, cat);
+ seed_count++;
+ }
+
+ Vect_destroy_line_struct(seed_points);
+ Vect_destroy_cats_struct(seed_cats);
+ Vect_close(&seed_Map);
+ }
+ else {
+ /* 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 ((c % skip[0] == 0) && (r % skip[1] == 0) &&
+ (d % skip[2] == 0)) {
+ seed.flowline = TRUE;
+ G_percent(seed_count, n_seeds, 1);
+ cat = seed_count + 1;
+ compute_flowline(®ion, &seed, velocity_field,
+ &integration, &fl_map, fl_cats,
+ fl_points, cat);
+ seed_count++;
+ }
+ }
+ }
+ }
+ }
+ Vect_destroy_line_struct(fl_points);
+ Vect_destroy_cats_struct(fl_cats);
+ Vect_build(&fl_map);
+ Vect_close(&fl_map);
+
+
+ return EXIT_SUCCESS;
+}
Property changes on: grass-addons/grass7/raster3d/r3.flow/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster3d/r3.flow/r3.flow.html
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/r3.flow.html (rev 0)
+++ grass-addons/grass7/raster3d/r3.flow/r3.flow.html 2014-07-10 18:51:58 UTC (rev 61235)
@@ -0,0 +1,21 @@
+<h2>DESCRIPTION</h2>
+
+<em>r3.flow</em>
+
+
+<h2>NOTES</h2>
+
+
+
+<h2>EXAMPLES</h2>
+
+<h2>AUTHORS</h2>
+Anna Petrasova
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="g.region.html">g.region</a>
+</em>
+
+<p><i>Last changed: $Date$</i>
Property changes on: grass-addons/grass7/raster3d/r3.flow/r3.flow.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster3d/r3.flow/r3flow_structs.h
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/r3flow_structs.h (rev 0)
+++ grass-addons/grass7/raster3d/r3.flow/r3flow_structs.h 2014-07-10 18:51:58 UTC (rev 61235)
@@ -0,0 +1,24 @@
+#ifndef R3FLOW_STRUCTS_H
+#define R3FLOW_STRUCTS_H
+
+#include <grass/raster3d.h>
+
+struct Seed
+{
+ double x;
+ double y;
+ double z;
+ int flowline;
+ int flowaccum;
+};
+
+struct Integration
+{
+ int direction;
+ char *unit;
+ double step;
+ double cell_size;
+ int limit;
+};
+
+#endif // R3FLOW_STRUCTS_H
Property changes on: grass-addons/grass7/raster3d/r3.flow/r3flow_structs.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list