[GRASS-SVN] r61276 - grass-addons/grass7/raster3d/r3.flow
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Jul 18 13:30:23 PDT 2014
Author: annakrat
Date: 2014-07-18 13:30:23 -0700 (Fri, 18 Jul 2014)
New Revision: 61276
Modified:
grass-addons/grass7/raster3d/r3.flow/flowline.c
grass-addons/grass7/raster3d/r3.flow/flowline.h
grass-addons/grass7/raster3d/r3.flow/gradient.c
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/r3flow_structs.h
Log:
r3.flow: added flowaccumulation and on-the-fly gradient computation
Modified: grass-addons/grass7/raster3d/r3.flow/flowline.c
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/flowline.c 2014-07-18 09:18:29 UTC (rev 61275)
+++ grass-addons/grass7/raster3d/r3.flow/flowline.c 2014-07-18 20:30:23 UTC (rev 61276)
@@ -10,7 +10,7 @@
for details.
\author Anna Petrasova
-*/
+ */
#include <grass/raster3d.h>
#include <grass/vector.h>
#include <grass/glocale.h>
@@ -31,10 +31,10 @@
\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, RASTER3D_Map *flowacc,
- struct Integration *integration,
+ struct Gradient_info *gradient_info,
+ RASTER3D_Map * flowacc, struct Integration *integration,
struct Map_info *flowline_vec, struct line_cats *cats,
struct line_pnts *points, const int cat)
{
@@ -66,7 +66,7 @@
}
count = 1;
while (count <= integration->limit) {
- if (get_velocity(region, velocity_field, point[0], point[1], point[2],
+ 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);
@@ -85,7 +85,7 @@
integration->cell_size);
delta_t *= integration->direction;
if (rk45_integrate_next
- (region, velocity_field, point, new_point,
+ (region, gradient_info, point, new_point,
&delta_t, min_step, max_step) < 0)
break;
@@ -93,34 +93,37 @@
Vect_append_point(points, new_point[0], new_point[1],
new_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;
- }
- }
+ 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];
}
@@ -128,12 +131,12 @@
}
if (seed->flowline) {
- if (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);
+ if (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);
}
}
Modified: grass-addons/grass7/raster3d/r3.flow/flowline.h
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/flowline.h 2014-07-18 09:18:29 UTC (rev 61275)
+++ grass-addons/grass7/raster3d/r3.flow/flowline.h 2014-07-18 20:30:23 UTC (rev 61276)
@@ -9,8 +9,9 @@
static const double VELOCITY_EPSILON = 1e-8;
void compute_flowline(RASTER3D_Region * region, const struct Seed *seed,
- RASTER3D_Map ** velocity_field,
- RASTER3D_Map *flowacc,
- struct Integration *integration, struct Map_info *flowline_vec,
- struct line_cats *cats, struct line_pnts *points, const int cat);
+ struct Gradient_info *gradient_info,
+ RASTER3D_Map * flowacc,
+ struct Integration *integration,
+ struct Map_info *flowline_vec, struct line_cats *cats,
+ struct line_pnts *points, const int cat);
#endif // FLOWLINE_H
Modified: grass-addons/grass7/raster3d/r3.flow/gradient.c
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/gradient.c 2014-07-18 09:18:29 UTC (rev 61275)
+++ grass-addons/grass7/raster3d/r3.flow/gradient.c 2014-07-18 20:30:23 UTC (rev 61276)
@@ -3,10 +3,10 @@
\brief Gradient computation
- Gradient computation (second order approximation)
- using central differencing scheme (plus forward and backward
- difference of second order approx.)
-
+ Gradient computation (second order approximation)
+ using central differencing scheme (plus forward and backward
+ difference of second order approx.)
+
(C) 2014 by the GRASS Development Team
This program is free software under the GNU General Public
@@ -35,7 +35,7 @@
4 * ACCESS(array, array->sx - 2, row, depth) +
ACCESS(array, array->sx - 3, row, depth)) / (2 * step[0]);
- for (col = 0; col < array->sx; col++) {
+ for (col = 1; col < array->sx - 1; col++) {
ACCESS(grad_x, col, row, depth) =
(ACCESS(array, col + 1, row, depth) -
ACCESS(array, col - 1, row, depth)) / (2 * step[0]);
@@ -45,19 +45,20 @@
for (depth = 0; depth < array->sz; depth++) {
for (col = 0; col < array->sx; col++) {
ACCESS(grad_y, col, 0, depth) =
- (-3 * ACCESS(array, col, 0, depth) +
- 4 * ACCESS(array, col, 1, depth) -
- ACCESS(array, col, 2, depth)) / (2 * step[1]);
+ -(-3 * ACCESS(array, col, 0, depth) +
+ 4 * ACCESS(array, col, 1, depth) -
+ ACCESS(array, col, 2, depth)) / (2 * step[1]);
ACCESS(grad_y, col, array->sy - 1, depth) =
- (3 * ACCESS(array, col, array->sy - 1, depth) -
- 4 * ACCESS(array, col, array->sy - 2, depth) +
- ACCESS(array, col, array->sy - 3, depth)) / (2 * step[1]);
+ -(3 * ACCESS(array, col, array->sy - 1, depth) -
+ 4 * ACCESS(array, col, array->sy - 2, depth) +
+ ACCESS(array, col, array->sy - 3, depth)) / (2 * step[1]);
- for (row = 0; row < array->sy; row++) {
+ for (row = 1; row < array->sy - 1; row++) {
+ /* is minus here? */
ACCESS(grad_y, col, row, depth) =
- (ACCESS(array, col, row + 1, depth) -
- ACCESS(array, col, row - 1, depth)) / (2 * step[1]);
+ -(ACCESS(array, col, row + 1, depth) -
+ ACCESS(array, col, row - 1, depth)) / (2 * step[1]);
}
}
}
@@ -73,11 +74,10 @@
4 * ACCESS(array, col, row, array->sz - 2) +
ACCESS(array, col, row, array->sz - 3)) / (2 * step[2]);
- for (depth = 0; depth < array->sz; depth++) {
- /* is minus here? */
- ACCESS(grad_y, col, row, depth) =
- -(ACCESS(array, col, row, depth + 1) -
- ACCESS(array, col, row, depth - 1)) / (2 * step[2]);
+ for (depth = 1; depth < array->sz - 1; depth++) {
+ ACCESS(grad_z, col, row, depth) =
+ (ACCESS(array, col, row, depth + 1) -
+ ACCESS(array, col, row, depth - 1)) / (2 * step[2]);
}
}
}
Modified: grass-addons/grass7/raster3d/r3.flow/integrate.c
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/integrate.c 2014-07-18 09:18:29 UTC (rev 61275)
+++ grass-addons/grass7/raster3d/r3.flow/integrate.c 2014-07-18 20:30:23 UTC (rev 61276)
@@ -21,6 +21,7 @@
#include "interpolate.h"
#include "integrate.h"
+#include "r3flow_structs.h"
/*!
\brief Computes vector norm.
@@ -56,11 +57,14 @@
}
-int get_velocity(RASTER3D_Region * region, RASTER3D_Map ** velocity_field,
+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)
{
- return interpolate_velocity(region, velocity_field, y, x, z,
+ if (gradient_info->compute_gradient)
+ return get_gradient(region, gradient_info, x, y, z, vel_x, vel_y, vel_z);
+
+ return interpolate_velocity(region, gradient_info->velocity_maps, y, x, z,
vel_x, vel_y, vel_z);
}
@@ -77,7 +81,7 @@
\return 0 success
\return -1 out of region or null values
*/
-static int rk45_next(RASTER3D_Region * region, RASTER3D_Map ** velocity_field,
+static int rk45_next(RASTER3D_Region * region, struct Gradient_info *gradient_info,
const double *point, double *next_point,
const double delta_t, double *error)
{
@@ -87,7 +91,7 @@
double sum_tmp;
int i, j, k;
- if (get_velocity(region, velocity_field, point[0], point[1], point[2],
+ if (get_velocity(region, gradient_info, point[0], point[1], point[2],
&vel_x, &vel_y, &vel_z) < 0)
return -1;
@@ -105,7 +109,7 @@
tmp_point[j] = point[j] + delta_t * sum_tmp;
}
if (get_velocity
- (region, velocity_field, tmp_point[0], tmp_point[1], tmp_point[2],
+ (region, gradient_info, tmp_point[0], tmp_point[1], tmp_point[2],
&vel_x, &vel_y, &vel_z) < 0)
return -1;
@@ -153,7 +157,7 @@
\return -1 out of region or null values
*/
int rk45_integrate_next(RASTER3D_Region * region,
- RASTER3D_Map ** velocity_field, const double *point,
+ struct Gradient_info *gradient_info, const double *point,
double *next_point, double *delta_t,
const double min_step, const double max_step)
{
@@ -175,7 +179,7 @@
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)
+ (region, gradient_info, point, next_point, *delta_t, error) == 0)
estimated_error = norm(error[0], error[1], error[2]);
else
return -1;
@@ -206,7 +210,7 @@
/* break when the adjustment was needed (not sure why) */
if (do_break) {
- if (rk45_next(region, velocity_field, point, next_point,
+ if (rk45_next(region, gradient_info, point, next_point,
*delta_t, error) < 0)
return -1;
break;
Modified: grass-addons/grass7/raster3d/r3.flow/integrate.h
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/integrate.h 2014-07-18 09:18:29 UTC (rev 61275)
+++ grass-addons/grass7/raster3d/r3.flow/integrate.h 2014-07-18 20:30:23 UTC (rev 61276)
@@ -3,6 +3,8 @@
#include <grass/raster3d.h>
+#include "r3flow_structs.h"
+
static const double MAX_ERROR = 1.0e-6;
static const double MIN_STEP = 0.01;
static const double MAX_STEP = 1;
@@ -23,13 +25,13 @@
};
double norm(const double x, const double y, const double z);
-int get_velocity(RASTER3D_Region * region, RASTER3D_Map ** velocity_field,
+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,
- RASTER3D_Map ** velocity_field, const double *point,
+ struct Gradient_info *gradient_info, const double *point,
double *next_point, double *delta_t,
const double min_step, const double max_step);
Modified: grass-addons/grass7/raster3d/r3.flow/interpolate.c
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/interpolate.c 2014-07-18 09:18:29 UTC (rev 61275)
+++ grass-addons/grass7/raster3d/r3.flow/interpolate.c 2014-07-18 20:30:23 UTC (rev 61276)
@@ -12,8 +12,12 @@
\author Anna Petrasova
*/
+#include <grass/gis.h>
#include <grass/raster3d.h>
+#include "r3flow_structs.h"
+#include "gradient.h"
+
/*!
\brief Finds 8 nearest voxels from a point.
@@ -170,7 +174,7 @@
/* compute weights */
get_relative_coords_for_interp(region, north, east, top,
- &rel_x, &rel_y, &rel_z);
+ &rel_x, &rel_y, &rel_z);
trilinear_interpolation(values, rel_x, rel_y, rel_z, interpolated);
*vel_x = interpolated[0];
@@ -179,3 +183,144 @@
return 0;
}
+
+int get_gradient(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)
+{
+
+ 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];
+ struct Array array;
+ struct Array grad_x, grad_y, grad_z;
+ struct Array *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, y, x, z, 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);
+ gradient(&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] =
+ ACCESS(grad_xyz[i], c + xshift, r + yshift,
+ d + zshift);
+ count++;
+ }
+ }
+ }
+ get_relative_coords_for_interp(region, y, x, z, &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;
+
+
+
+}
Modified: grass-addons/grass7/raster3d/r3.flow/interpolate.h
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/interpolate.h 2014-07-18 09:18:29 UTC (rev 61275)
+++ grass-addons/grass7/raster3d/r3.flow/interpolate.h 2014-07-18 20:30:23 UTC (rev 61276)
@@ -3,8 +3,14 @@
#include <grass/raster3d.h>
-int interpolate_velocity(RASTER3D_Region *region, RASTER3D_Map **map,
+#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 x,
+ const double y, const double z, double *vel_x, double *vel_y,
+ double *vel_z);
#endif // INTERPOLATE_H
Modified: grass-addons/grass7/raster3d/r3.flow/main.c
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/main.c 2014-07-18 09:18:29 UTC (rev 61275)
+++ grass-addons/grass7/raster3d/r3.flow/main.c 2014-07-18 20:30:23 UTC (rev 61276)
@@ -52,47 +52,63 @@
}
-static void load_input_velocity_maps(struct Option *vector_opt,
- RASTER3D_Map ** velocity_field,
+static void load_input_raster3d_maps(struct Option *scalar_opt,
+ struct Option *vector_opt,
+ struct Gradient_info *gradient_info,
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], ""),
+ 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 (!velocity_field[i])
+ if (!gradient_info->scalar_map)
Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
- vector_opt->answers[i]);
+ 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) {
-
+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"));
+ 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,
- *unit_opt, *step_opt, *limit_opt, *skip_opt;
+ *scalar_opt, *unit_opt, *step_opt, *limit_opt, *skip_opt;
struct Flag *up_flag;
struct GModule *module;
RASTER3D_Region region;
- RASTER3D_Map *velocity_field[3];
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;
@@ -113,9 +129,12 @@
module->description = _("Compute flow lines.");
+ scalar_opt = G_define_standard_option(G_OPT_R3_INPUT);
+ scalar_opt->required = NO;
+
vector_opt = G_define_standard_option(G_OPT_R3_INPUTS);
vector_opt->key = "vector_field";
- vector_opt->required = YES;
+ vector_opt->required = NO;
vector_opt->description = _("Names of three 3D raster maps describing "
"x, y, z components of vector field");
@@ -136,7 +155,8 @@
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->description =
+ _("Name for output flowaccumulation 3D raster");
unit_opt = G_define_option();
unit_opt->key = "unit";
@@ -191,6 +211,13 @@
check_vector_input_maps(vector_opt, seed_opt);
+ if (scalar_opt->answer || vector_opt->answers) {
+ if (scalar_opt->answer && vector_opt->answers)
+ G_fatal_error(_("Options 'input' and 'vector_field' are mutually exclusive."));
+ }
+ else
+ G_fatal_error(_("Use one of options 'input' and 'vector_field'."));
+
Rast3d_init_defaults();
Rast3d_get_window(®ion);
@@ -219,8 +246,7 @@
skip[i] = atoi(skip_opt->answers[i]);
}
else {
- G_fatal_error
- ("Please provide 3 integer values for skip option.");
+ G_fatal_error(_("Please provide 3 integer values for skip option."));
}
}
}
@@ -232,30 +258,32 @@
}
/* open raster 3D maps of velocity components */
- load_input_velocity_maps(vector_opt, velocity_field, ®ion);
+ gradient_info.initialized = FALSE;
+ load_input_raster3d_maps(scalar_opt, vector_opt, &gradient_info, ®ion);
+
/* 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, ®ion,
- FCELL_TYPE, 32);
-
-
- if (!flowacc)
- Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
- flowacc_opt->answer);
- init_flowaccum(®ion, flowacc);
+ flowacc = Rast3d_open_new_opt_tile_size(flowacc_opt->answer,
+ RASTER3D_USE_CACHE_DEFAULT,
+ ®ion, FCELL_TYPE, 32);
+
+
+ if (!flowacc)
+ Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
+ flowacc_opt->answer);
+ init_flowaccum(®ion, 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);
+ 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;
@@ -270,13 +298,13 @@
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]);
- }
+ 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);
@@ -305,8 +333,8 @@
}
G_percent(seed_count, n_seeds, 1);
cat = seed_count + 1;
- compute_flowline(®ion, &seed, velocity_field, flowacc, &integration,
- &fl_map, fl_cats, fl_points, cat);
+ compute_flowline(®ion, &seed, &gradient_info, flowacc,
+ &integration, &fl_map, fl_cats, fl_points, cat);
seed_count++;
}
@@ -329,15 +357,15 @@
seed.flowaccum = FALSE;
if (flowacc_opt->answer)
seed.flowaccum = TRUE;
-
+
if (flowlines_opt->answer && (c % skip[0] == 0) &&
- (r % skip[1] == 0) && (d % skip[2] == 0))
+ (r % skip[1] == 0) && (d % skip[2] == 0))
seed.flowline = TRUE;
-
+
if (seed.flowaccum || seed.flowline) {
G_percent(seed_count, n_seeds, 1);
cat = seed_count + 1;
- compute_flowline(®ion, &seed, velocity_field,
+ compute_flowline(®ion, &seed, &gradient_info,
flowacc, &integration, &fl_map,
fl_cats, fl_points, cat);
seed_count++;
@@ -347,14 +375,14 @@
}
}
if (flowlines_opt->answer) {
- Vect_destroy_line_struct(fl_points);
- Vect_destroy_cats_struct(fl_cats);
- Vect_build(&fl_map);
- Vect_close(&fl_map);
+ 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);
+ Rast3d_close(flowacc);
return EXIT_SUCCESS;
Modified: grass-addons/grass7/raster3d/r3.flow/r3flow_structs.h
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/r3flow_structs.h 2014-07-18 09:18:29 UTC (rev 61275)
+++ grass-addons/grass7/raster3d/r3.flow/r3flow_structs.h 2014-07-18 20:30:23 UTC (rev 61276)
@@ -1,6 +1,8 @@
#ifndef R3FLOW_STRUCTS_H
#define R3FLOW_STRUCTS_H
+#include <grass/raster3d.h>
+
struct Seed
{
double x;
@@ -27,6 +29,16 @@
int sz;
};
+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;
+};
+
#define ACCESS(arr, x, y, z) ((arr)->array[(arr)->sx * (arr)->sy * (z) + (arr)->sx * (y) + (x)])
#endif // R3FLOW_STRUCTS_H
More information about the grass-commit
mailing list