[GRASS-SVN] r67210 - in grass/trunk/raster/r.in.lidar: . testsuite
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Dec 17 19:44:10 PST 2015
Author: wenzeslaus
Date: 2015-12-17 19:44:10 -0800 (Thu, 17 Dec 2015)
New Revision: 67210
Modified:
grass/trunk/raster/r.in.lidar/Makefile
grass/trunk/raster/r.in.lidar/local_proto.h
grass/trunk/raster/r.in.lidar/main.c
grass/trunk/raster/r.in.lidar/point_binning.c
grass/trunk/raster/r.in.lidar/point_binning.h
grass/trunk/raster/r.in.lidar/testsuite/test_base_resolution.sh
Log:
r.in.lidar: vector output with grid-based decimation [news]
Creates one point per raster cell. Point gets
XYZ which are mean of all points in the cell.
No categories are created or preserved as it
would be unclear what to take (we cannot average
retrun number).
Some improvements are needed for history and done
message code. Vector writting could be more separated.
Modified: grass/trunk/raster/r.in.lidar/Makefile
===================================================================
--- grass/trunk/raster/r.in.lidar/Makefile 2015-12-17 23:19:23 UTC (rev 67209)
+++ grass/trunk/raster/r.in.lidar/Makefile 2015-12-18 03:44:10 UTC (rev 67210)
@@ -2,15 +2,14 @@
PGM = r.in.lidar
-LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB) $(GPROJLIB) $(LASLIBS) $(SEGMENTLIB)
-DEPENDENCIES = $(RASTERDEP) $(GISDEP) $(SEGMENTDEP)
+LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB) $(GPROJLIB) $(LASLIBS) $(SEGMENTLIB) $(VECTORLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP) $(SEGMENTDEP) $(VECTORDEP)
-EXTRA_INC = $(PROJINC) $(LASINC)
-EXTRA_CFLAGS = $(GDALCFLAGS)
+EXTRA_INC = $(LASINC) $(VECT_INC) $(PROJINC)
+EXTRA_CFLAGS = $(VECT_CFLAGS) $(GDALCFLAGS)
include $(MODULE_TOPDIR)/include/Make/Module.make
ifneq ($(USE_LIBLAS),)
default: cmd
endif
-
Modified: grass/trunk/raster/r.in.lidar/local_proto.h
===================================================================
--- grass/trunk/raster/r.in.lidar/local_proto.h 2015-12-17 23:19:23 UTC (rev 67209)
+++ grass/trunk/raster/r.in.lidar/local_proto.h 2015-12-18 03:44:10 UTC (rev 67210)
@@ -27,6 +27,7 @@
#define BUFFSIZE 256
+#define METHOD_NONE 0
#define METHOD_N 1
#define METHOD_MIN 2
#define METHOD_MAX 3
@@ -89,4 +90,21 @@
void string_list_from_one_item(struct StringList *string_list, char *item);
void string_list_free(struct StringList *string_list);
+/* forward declarations */
+struct Map_info;
+struct line_pnts;
+struct line_cats;
+
+struct VectorWriter
+{
+ struct Map_info *info;
+ struct line_pnts *points;
+ struct line_cats *cats;
+#ifdef HAVE_LONG_LONG_INT
+ unsigned long long count;
+#else
+ unsigned long count;
+#endif
+};
+
#endif /* __LOCAL_PROTO_H__ */
Modified: grass/trunk/raster/r.in.lidar/main.c
===================================================================
--- grass/trunk/raster/r.in.lidar/main.c 2015-12-17 23:19:23 UTC (rev 67209)
+++ grass/trunk/raster/r.in.lidar/main.c 2015-12-18 03:44:10 UTC (rev 67210)
@@ -28,6 +28,7 @@
#include <grass/segment.h>
#include <grass/gprojects.h>
#include <grass/glocale.h>
+#include <grass/vector.h>
#include <liblas/capi/liblas.h>
#include "local_proto.h"
@@ -79,8 +80,10 @@
struct Option *method_opt, *base_raster_opt, *zrange_opt, *zscale_opt;
struct Option *trim_opt, *pth_opt, *res_opt;
struct Option *file_list_opt;
+ struct Option *voutput_opt;
struct Flag *print_flag, *scan_flag, *shell_style, *over_flag, *extents_flag, *intens_flag;
struct Flag *base_rast_res_flag;
+ struct Flag *notopo_flag;
/* LAS */
LASReaderH LAS_reader;
@@ -110,6 +113,7 @@
input_opt->guisection = _("Input");
output_opt = G_define_standard_option(G_OPT_R_OUTPUT);
+ output_opt->required = NO;
output_opt->guisection = _("Output");
file_list_opt = G_define_standard_option(G_OPT_F_INPUT);
@@ -119,6 +123,14 @@
file_list_opt->required = NO;
file_list_opt->guisection = _("Input");
+ voutput_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+ voutput_opt->key = "vector_output";
+ voutput_opt->required = NO;
+ voutput_opt->label = _("Grid-decimated point cloud");
+ voutput_opt->description = _("Grid-decimated point cloud with"
+ " XYZ coordinates which are mean for all points in a raster cell");
+ voutput_opt->guisection = _("Output");
+
method_opt = G_define_option();
method_opt->key = "method";
method_opt->type = TYPE_STRING;
@@ -254,6 +266,10 @@
base_rast_res_flag->description =
_("Use base raster actual resolution instead of computational region");
+ notopo_flag = G_define_standard_flag(G_FLG_V_TOPO);
+
+ G_option_required(output_opt, voutput_opt, NULL);
+
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -269,7 +285,10 @@
}
/* parse input values */
- outmap = output_opt->answer;
+ if (output_opt->answer)
+ outmap = output_opt->answer;
+ else
+ outmap = NULL;
if (shell_style->answer && !scan_flag->answer) {
scan_flag->answer = 1; /* pointer not int, so set = shell_style->answer ? */
@@ -378,7 +397,8 @@
}
}
- point_binning_set(&point_binning, method_opt->answer, pth_opt->answer, trim_opt->answer);
+ point_binning_set(&point_binning, method_opt->answer, pth_opt->answer,
+ trim_opt->answer, voutput_opt->answer ? TRUE : FALSE);
base_array = NULL;
@@ -466,11 +486,26 @@
}
/* open output map */
- out_fd = Rast_open_new(outmap, rtype);
+ if (outmap)
+ out_fd = Rast_open_new(outmap, rtype);
+ else
+ out_fd = 0; /* TODO: is this correct? */
/* allocate memory for a single row of output data */
raster_row = Rast_allocate_output_buf(rtype);
+ struct Map_info voutput;
+ struct VectorWriter vector_writer;
+ if (voutput_opt->answer) {
+ if (Vect_open_new(&voutput, voutput_opt->answer, 1) < 0)
+ G_fatal_error(_("Unable to create vector map <%s>"), voutput_opt->answer);
+ vector_writer.info = &voutput;
+ vector_writer.points = Vect_new_line_struct();
+ vector_writer.cats = Vect_new_cats_struct();
+ vector_writer.count = 0;
+ Vect_hist_command(&voutput);
+ }
+
G_message(_("Reading data ..."));
count_total = line_total = 0;
@@ -590,7 +625,8 @@
count++;
/* G_debug(5, "x: %f, y: %f, z: %f", x, y, z); */
- update_value(&point_binning, &bin_index_nodes, cols, arr_row, arr_col, rtype, z);
+ update_value(&point_binning, &bin_index_nodes, cols,
+ arr_row, arr_col, rtype, x, y, z);
} /* while !EOF of one input file */
/* close input LAS file */
LASReader_Destroy(LAS_reader);
@@ -604,9 +640,12 @@
/* calc stats and output */
G_message(_("Writing to map ..."));
for (row = 0; row < rows; row++) {
- write_values(&point_binning, &bin_index_nodes, raster_row, row, cols, rtype);
+ /* potentially vector writing can be independent on the binning */
+ write_values(&point_binning, &bin_index_nodes, raster_row, row,
+ cols, rtype, &vector_writer);
/* write out line of raster data */
- Rast_put_row(out_fd, raster_row, rtype);
+ if (outmap)
+ Rast_put_row(out_fd, raster_row, rtype);
}
/* free memory */
@@ -621,17 +660,29 @@
G_free(raster_row);
/* close raster file & write history */
- Rast_close(out_fd);
+ if (outmap)
+ Rast_close(out_fd);
- sprintf(title, "Raw x,y,z data binned into a raster grid by cell %s",
- method_opt->answer);
- Rast_put_cell_title(outmap, title);
+ if (outmap) {
+ sprintf(title, "Raw x,y,z data binned into a raster grid by cell %s",
+ method_opt->answer);
+ Rast_put_cell_title(outmap, title);
- Rast_short_history(outmap, "raster", &history);
- Rast_command_history(&history);
- Rast_set_history(&history, HIST_DATSRC_1, infile);
- Rast_write_history(outmap, &history);
+ Rast_short_history(outmap, "raster", &history);
+ Rast_command_history(&history);
+ Rast_set_history(&history, HIST_DATSRC_1, infile);
+ Rast_write_history(outmap, &history);
+ }
+ /* close output vector map */
+ if (voutput_opt->answer) {
+ if (!notopo_flag->answer)
+ Vect_build(vector_writer.info);
+ Vect_close(vector_writer.info);
+ Vect_destroy_line_struct(vector_writer.points);
+ Vect_destroy_cats_struct(vector_writer.cats);
+ }
+
if (infiles.num_items > 1) {
sprintf(buff, _("Raster map <%s> created."
" %lu points from %d files found in region."),
Modified: grass/trunk/raster/r.in.lidar/point_binning.c
===================================================================
--- grass/trunk/raster/r.in.lidar/point_binning.c 2015-12-17 23:19:23 UTC (rev 67209)
+++ grass/trunk/raster/r.in.lidar/point_binning.c 2015-12-18 03:44:10 UTC (rev 67210)
@@ -18,6 +18,7 @@
#include <grass/gis.h>
#include <grass/glocale.h>
#include <grass/raster.h>
+#include <grass/vector.h>
#include "point_binning.h"
#include "local_proto.h"
@@ -113,7 +114,7 @@
}
void point_binning_set(struct PointBinning *point_binning, char *method,
- char *percentile, char *trim)
+ char *percentile, char *trim, int bin_coordinates)
{
/* figure out what maps we need in memory */
@@ -131,12 +132,14 @@
skewness n array index to linked list
trimmean n array index to linked list
*/
+ point_binning->method = METHOD_NONE;
point_binning->bin_n = FALSE;
point_binning->bin_min = FALSE;
point_binning->bin_max = FALSE;
point_binning->bin_sum = FALSE;
point_binning->bin_sumsq = FALSE;
point_binning->bin_index = FALSE;
+ point_binning->bin_coordinates = FALSE;
point_binning->n_array = NULL;
point_binning->min_array = NULL;
@@ -144,6 +147,8 @@
point_binning->sum_array = NULL;
point_binning->sumsq_array = NULL;
point_binning->index_array = NULL;
+ point_binning->x_array = NULL;
+ point_binning->y_array = NULL;
if (strcmp(method, "n") == 0) {
point_binning->method = METHOD_N;
@@ -213,6 +218,13 @@
point_binning->method = METHOD_TRIMMEAN;
point_binning->bin_index = TRUE;
}
+ if (bin_coordinates) {
+ /* x, y */
+ point_binning->bin_coordinates = TRUE;
+ /* z and n */
+ point_binning->bin_sum = TRUE;
+ point_binning->bin_n = TRUE;
+ }
}
@@ -251,6 +263,12 @@
if (point_binning->bin_index)
point_binning->index_array =
G_calloc((size_t) rows * (cols + 1), Rast_cell_size(CELL_TYPE));
+ if (point_binning->bin_coordinates) {
+ point_binning->x_array =
+ G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
+ point_binning->y_array =
+ G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
+ }
/* TODO: perhaps none of them needs to be freed */
/* and then free it again */
@@ -266,6 +284,10 @@
G_free(point_binning->sumsq_array);
if (point_binning->bin_index)
G_free(point_binning->index_array);
+ if (point_binning->bin_coordinates) {
+ G_free(point_binning->x_array);
+ G_free(point_binning->y_array);
+ }
}
@@ -308,6 +330,15 @@
G_calloc((size_t) rows * (cols + 1), Rast_cell_size(CELL_TYPE));
blank_array(point_binning->index_array, rows, cols, CELL_TYPE, -1); /* fill with NULLs */
}
+ if (point_binning->bin_coordinates) {
+ G_debug(2, "allocating x_array and y_array");
+ point_binning->x_array =
+ G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
+ blank_array(point_binning->x_array, rows, cols, rtype, 0);
+ point_binning->y_array =
+ G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
+ blank_array(point_binning->y_array, rows, cols, rtype, 0);
+ }
}
void point_binning_free(struct PointBinning *point_binning,
@@ -330,6 +361,10 @@
bin_index_nodes->max_nodes = 0;
bin_index_nodes->nodes = NULL;
}
+ if (point_binning->bin_coordinates) {
+ G_free(point_binning->x_array);
+ G_free(point_binning->y_array);
+ }
}
void write_variance(void *raster_row, void *n_array, void *sum_array,
@@ -617,7 +652,8 @@
void write_values(struct PointBinning *point_binning,
struct BinIndex *bin_index_nodes, void *raster_row, int row,
- int cols, RASTER_MAP_TYPE rtype)
+ int cols, RASTER_MAP_TYPE rtype,
+ struct VectorWriter *vector_writer)
{
void *ptr = NULL;
int col;
@@ -707,13 +743,75 @@
break;
default:
- G_fatal_error("?");
+ G_debug(2, "No method selected");
}
+ if (point_binning->bin_coordinates) {
+ for (col = 0; col < cols; col++) {
+ size_t offset = (row * cols + col) * Rast_cell_size(rtype);
+ size_t n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
+ int n = Rast_get_c_value(point_binning->n_array + n_offset,
+ CELL_TYPE);
+
+ if (n == 0)
+ continue;
+
+ double sum_x =
+ Rast_get_d_value(point_binning->x_array + offset, rtype);
+ double sum_y =
+ Rast_get_d_value(point_binning->y_array + offset, rtype);
+ /* TODO: we do this also in mean writing */
+ double sum_z =
+ Rast_get_d_value(point_binning->sum_array + offset, rtype);
+
+ /* We are not writing any categories. They are not needed
+ * and potentially it is too much trouble to do it and it is
+ * unclear what to write. Not writing them is also a little
+ * bit faster. */
+ Vect_append_point(vector_writer->points, sum_x, sum_y, sum_z / n);
+ Vect_write_line(vector_writer->info, GV_POINT,
+ vector_writer->points, vector_writer->cats);
+ Vect_reset_line(vector_writer->points);
+ vector_writer->count++;
+ }
+ }
}
+/* TODO: duplication with support.c, refactoring needed */
+static void *get_cell_ptr(void *array, int cols, int row, int col,
+ RASTER_MAP_TYPE map_type)
+{
+ return G_incr_void_ptr(array,
+ ((row * (size_t) cols) +
+ col) * Rast_cell_size(map_type));
+}
+
+int update_val(void *array, int cols, int row, int col,
+ RASTER_MAP_TYPE map_type, double value)
+{
+ void *ptr = get_cell_ptr(array, cols, row, col, map_type);
+
+ Rast_set_d_value(ptr, value, map_type);
+ return 0;
+}
+
+int update_moving_mean(void *array, int cols, int row, int col,
+ RASTER_MAP_TYPE rtype, double value, int n)
+{
+ /* for xy we do this check twice */
+ if (n != 0) {
+ double m_v;
+
+ row_array_get_value_row_col(array, row, col, cols, rtype, &m_v);
+ value = m_v + (value - m_v) / n;
+ }
+ /* else we just write the initial value */
+ return update_val(array, cols, row, col, rtype, value);;
+}
+
void update_value(struct PointBinning *point_binning,
struct BinIndex *bin_index_nodes, int cols, int arr_row,
- int arr_col, RASTER_MAP_TYPE rtype, double z)
+ int arr_col, RASTER_MAP_TYPE rtype, double x, double y,
+ double z)
{
if (point_binning->bin_n)
update_n(point_binning->n_array, cols, arr_row, arr_col);
@@ -732,4 +830,15 @@
if (point_binning->bin_index)
update_bin_index(bin_index_nodes, point_binning->index_array, cols,
arr_row, arr_col, rtype, z);
+ if (point_binning->bin_coordinates) {
+ /* this assumes that n is already computed for this xyz */
+ void *ptr = get_cell_ptr(point_binning->n_array, cols, arr_row,
+ arr_col, CELL_TYPE);
+ int n = Rast_get_c_value(ptr, CELL_TYPE);
+
+ update_moving_mean(point_binning->x_array, cols, arr_row, arr_col,
+ rtype, x, n);
+ update_moving_mean(point_binning->y_array, cols, arr_row, arr_col,
+ rtype, y, n);
+ }
}
Modified: grass/trunk/raster/r.in.lidar/point_binning.h
===================================================================
--- grass/trunk/raster/r.in.lidar/point_binning.h 2015-12-17 23:19:23 UTC (rev 67209)
+++ grass/trunk/raster/r.in.lidar/point_binning.h 2015-12-18 03:44:10 UTC (rev 67210)
@@ -17,6 +17,8 @@
#include <grass/raster.h>
+/* forward declaration */
+struct Map_info;
struct node
{
@@ -41,6 +43,7 @@
int bin_sum;
int bin_sumsq;
int bin_index;
+ int bin_coordinates;
void *n_array;
void *min_array;
@@ -48,6 +51,8 @@
void *sum_array;
void *sumsq_array;
void *index_array;
+ void *x_array;
+ void *y_array;
int pth;
double trim;
@@ -58,7 +63,7 @@
int cols, RASTER_MAP_TYPE rtype);
void point_binning_set(struct PointBinning *point_binning, char *method,
- char *percentile, char *trim);
+ char *percentile, char *trim, int coordinates);
void point_binning_allocate(struct PointBinning *point_binning, int rows,
int cols, RASTER_MAP_TYPE rtype);
@@ -84,12 +89,17 @@
void *index_array, int row, int cols,
RASTER_MAP_TYPE rtype, double trim);
+/* forward declaration */
+struct VectorWriter;
+
void write_values(struct PointBinning *point_binning,
struct BinIndex *bin_index_nodes, void *raster_row, int row,
- int cols, RASTER_MAP_TYPE rtype);
+ int cols, RASTER_MAP_TYPE rtype,
+ struct VectorWriter *vector_writer);
void update_value(struct PointBinning *point_binning,
struct BinIndex *bin_index_nodes, int cols, int arr_row,
- int arr_col, RASTER_MAP_TYPE rtype, double z);
+ int arr_col, RASTER_MAP_TYPE rtype, double x, double y,
+ double z);
#endif /* __POINT_BINNING_H__ */
Modified: grass/trunk/raster/r.in.lidar/testsuite/test_base_resolution.sh
===================================================================
--- grass/trunk/raster/r.in.lidar/testsuite/test_base_resolution.sh 2015-12-17 23:19:23 UTC (rev 67209)
+++ grass/trunk/raster/r.in.lidar/testsuite/test_base_resolution.sh 2015-12-18 03:44:10 UTC (rev 67210)
@@ -30,7 +30,7 @@
r.univar ${basename}with_base -g \
| grep -ve "=0$" | grep -ve "=-nan$" | grep -e "=[^12356789]$"
-echo "Test successful"
-echo "When running manually maps can be now removed with:"
-echo " g.remove type=rast pattern='test_rinlidar_*' -f"
-echo "However, the region was changed to whatever the test needed."
+echo "Test successful
+When running manually maps can be now removed with:
+ g.remove type=rast pattern='test_rinlidar_*' -f
+However, the region was changed to whatever the test needed."
More information about the grass-commit
mailing list