[GRASS-SVN] r61257 - grass-addons/grass7/raster3d/r3.flow
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jul 15 19:44:20 PDT 2014
Author: annakrat
Date: 2014-07-15 19:44:20 -0700 (Tue, 15 Jul 2014)
New Revision: 61257
Added:
grass-addons/grass7/raster3d/r3.flow/voxel_traversal.c
grass-addons/grass7/raster3d/r3.flow/voxel_traversal.h
Modified:
grass-addons/grass7/raster3d/r3.flow/flowline.c
grass-addons/grass7/raster3d/r3.flow/flowline.h
grass-addons/grass7/raster3d/r3.flow/main.c
grass-addons/grass7/raster3d/r3.flow/r3flow_structs.h
Log:
r3.flow: added flowaccumulation (not well tested yet)
Modified: grass-addons/grass7/raster3d/r3.flow/flowline.c
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/flowline.c 2014-07-15 13:37:55 UTC (rev 61256)
+++ grass-addons/grass7/raster3d/r3.flow/flowline.c 2014-07-16 02:44:20 UTC (rev 61257)
@@ -18,6 +18,7 @@
#include "r3flow_structs.h"
#include "integrate.h"
#include "flowline.h"
+#include "voxel_traversal.h"
/*!
\brief Computes flowline by integrating velocity field.
@@ -32,22 +33,33 @@
\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 ** velocity_field, RASTER3D_Map *flowacc,
struct Integration *integration,
struct Map_info *flowline_vec, struct line_cats *cats,
struct line_pnts *points, const int cat)
{
- int i, count;
+ 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;
+ int *trav_coords;
+ int size, trav_count;
point[0] = seed->x;
point[1] = seed->y;
point[2] = seed->z;
+ last_col = last_row = last_depth = -1;
+
+ size = 5;
+ trav_coords = G_malloc(3 * size * sizeof(int));
+
if (seed->flowline) {
/* append first point */
Vect_append_point(points, seed->x, seed->y, seed->z);
@@ -80,18 +92,45 @@
if (seed->flowline)
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);
+ coor_diff = (fabs(last_col - col) + fabs(last_row - row) + fabs(last_row - row));
+ /* if not run for the 1. time and previous and next point coordinates
+ differ by more than 1 voxel coordinate */
+ if (last_col >= 0 && 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 && 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);
+ 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);
}
- 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-15 13:37:55 UTC (rev 61256)
+++ grass-addons/grass7/raster3d/r3.flow/flowline.h 2014-07-16 02:44:20 UTC (rev 61257)
@@ -4,11 +4,13 @@
#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,
RASTER3D_Map ** velocity_field,
- struct Integration *integration,
- struct Map_info *flowline_vec, struct line_cats *cats,
- struct line_pnts *points, const int cat);
+ 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/main.c
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/main.c 2014-07-15 13:37:55 UTC (rev 61256)
+++ grass-addons/grass7/raster3d/r3.flow/main.c 2014-07-16 02:44:20 UTC (rev 61257)
@@ -71,14 +71,26 @@
}
}
+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,
+ struct Option *vector_opt, *seed_opt, *flowlines_opt, *flowacc_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 Map_info seed_Map;
@@ -118,9 +130,14 @@
flowlines_opt = G_define_standard_option(G_OPT_V_OUTPUT);
flowlines_opt->key = "flowline";
- flowlines_opt->required = YES;
+ flowlines_opt->required = NO;
flowlines_opt->description = _("Name for vector map of flow lines");
+ 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");
+
unit_opt = G_define_option();
unit_opt->key = "unit";
unit_opt->type = TYPE_STRING;
@@ -217,15 +234,30 @@
/* open raster 3D maps of velocity components */
load_input_velocity_maps(vector_opt, velocity_field, ®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);
+ }
+
/* 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);
+ 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);
+ }
- Vect_hist_command(&fl_map);
-
n_seeds = 0;
/* open vector map of seeds */
if (seed_opt->answer) {
@@ -237,10 +269,14 @@
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]);
+ 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);
@@ -269,7 +305,7 @@
}
G_percent(seed_count, n_seeds, 1);
cat = seed_count + 1;
- compute_flowline(®ion, &seed, velocity_field, &integration,
+ compute_flowline(®ion, &seed, velocity_field, flowacc, &integration,
&fl_map, fl_cats, fl_points, cat);
seed_count++;
}
@@ -278,7 +314,7 @@
Vect_destroy_cats_struct(seed_cats);
Vect_close(&seed_Map);
}
- else {
+ 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++) {
@@ -291,26 +327,35 @@
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)) {
+ if (flowacc_opt->answer)
+ seed.flowaccum = TRUE;
+
+ if (flowlines_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);
cat = seed_count + 1;
compute_flowline(®ion, &seed, velocity_field,
- &integration, &fl_map, fl_cats,
- fl_points, cat);
+ flowacc, &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);
+ if (flowlines_opt->answer) {
+ 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;
}
Modified: grass-addons/grass7/raster3d/r3.flow/r3flow_structs.h
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/r3flow_structs.h 2014-07-15 13:37:55 UTC (rev 61256)
+++ grass-addons/grass7/raster3d/r3.flow/r3flow_structs.h 2014-07-16 02:44:20 UTC (rev 61257)
@@ -1,8 +1,6 @@
#ifndef R3FLOW_STRUCTS_H
#define R3FLOW_STRUCTS_H
-#include <grass/raster3d.h>
-
struct Seed
{
double x;
Added: grass-addons/grass7/raster3d/r3.flow/voxel_traversal.c
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/voxel_traversal.c (rev 0)
+++ grass-addons/grass7/raster3d/r3.flow/voxel_traversal.c 2014-07-16 02:44:20 UTC (rev 61257)
@@ -0,0 +1,93 @@
+#include <math.h>
+#include <grass/raster3d.h>
+
+void traverse(RASTER3D_Region * region, double *start, double *end,
+ int *coordinates, int *size, int *coor_count)
+{
+ double dx, dy, dz;
+ double 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;
+ }
+ }
+ coordinates[count * 3 + 0] = x;
+ coordinates[count * 3 + 1] = region->rows - y - 1;
+ coordinates[count * 3 + 2] = z;
+ count++;
+ 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;
+
+ /* reallocation for cases when the steps would be too big */
+ if (*size <= count) {
+ *size = 2 * (*size);
+ coordinates = G_realloc(coordinates, (*size) * sizeof(int));
+ }
+ }
+ *coor_count = count;
+}
Added: grass-addons/grass7/raster3d/r3.flow/voxel_traversal.h
===================================================================
--- grass-addons/grass7/raster3d/r3.flow/voxel_traversal.h (rev 0)
+++ grass-addons/grass7/raster3d/r3.flow/voxel_traversal.h 2014-07-16 02:44:20 UTC (rev 61257)
@@ -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
More information about the grass-commit
mailing list