v.decimate: grid and count based decimation [news]

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

vector/v.decimate/Makefile

@@ -0,0 +1,14 @@
+
+MODULE_TOPDIR = ../..
+
+PGM = v.decimate
+
+LIBES = $(VECTORLIB) $(DBMILIB) $(RASTERLIB) $(GISLIB) $(BTREE2LIB) $(SEGMENTLIB)
+DEPENDENCIES = $(VECTORDEP) $(DBMIDEP) $(RASTERDEP) $(GISDEP)  $(SEGMENTDEP)
+
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

vector/v.decimate/count_decimation.c

@@ -0,0 +1,132 @@
+/****************************************************************************
+ *
+ * MODULE:       v.decimate
+ * AUTHOR(S):    Vaclav Petras
+ * PURPOSE:      Reduce the number of points in a vector map
+ * COPYRIGHT:    (C) 2015 by the GRASS Development Team
+ *
+ *               This program is free software under the GNU General Public
+ *               License (>=v2). Read the COPYING file that comes with GRASS
+ *               for details.
+ *
+ *****************************************************************************/
+
+
+/* TODO: change int */
+/* TODO: revise names */
+
+#include "count_decimation.h"
+
+#include <grass/gis.h>
+
+#include <stdlib.h>
+
+
+void count_decimation_init(struct CountDecimationControl *control,
+                           int *skip, int *preserve,
+                           int *offset, int *limit)
+{
+    control->skip_every = 0;
+    control->preserve_every = 0;
+    /* counter used by both but that's ok, skip and preserve are exclusive */
+    control->every_counter = 0;
+    control->n_count_filtered = 0;
+    control->offset_n = 0;
+    control->offset_n_counter = 0;
+    control->limit_n = 0;
+    control->limit_n_counter = 0;
+    if (skip)
+        control->skip_every = *skip;
+    if (preserve)
+        control->preserve_every = *preserve;
+    if (offset)
+        control->offset_n = *offset;
+    if (limit)
+        control->limit_n = *limit;
+}
+
+
+int count_decimation_is_valid(struct CountDecimationControl *control)
+{
+    if (control->skip_every == 1)
+        return FALSE;
+    if (control->skip_every && control->preserve_every > 1)
+        return FALSE;
+    return TRUE;
+}
+
+
+int count_decimation_is_noop(struct CountDecimationControl *control)
+{
+    if (control->skip_every < 2 && control->preserve_every < 2
+            && !control->offset_n && !control->limit_n)
+        return TRUE;
+    return FALSE;
+}
+
+void count_decimation_init_from_str(struct CountDecimationControl *control,
+                                    const char *skip, const char *preserve,
+                                    const char *offset, const char *limit)
+{
+    control->skip_every = 0;
+    control->preserve_every = 0;
+    control->every_counter = 0;
+    control->n_count_filtered = 0;
+    control->offset_n = 0;
+    control->offset_n_counter = 0;
+    control->limit_n = 0;
+    control->limit_n_counter = 0;
+    /* TODO: atoi is probably not appropriate */
+    if (skip)
+        control->skip_every = atoi(skip);
+    if (preserve)
+        control->preserve_every = atoi(preserve);
+    if (offset)
+        control->offset_n = atoi(offset);
+    if (limit)
+        control->limit_n = atoi(limit);
+}
+
+
+/* TODO: eliminate noop cases */
+int count_decimation_is_out(struct CountDecimationControl *control)
+{
+    if (control->offset_n) {
+        control->offset_n_counter++;
+        if (control->offset_n_counter < control->offset_n)
+            return TRUE;
+        else
+            control->offset_n = 0;  /* disable offset check */
+    }
+    if (control->skip_every) {
+        control->every_counter++;
+        if (control->every_counter == control->skip_every) {
+            control->n_count_filtered++;
+            control->every_counter = 0;
+            return TRUE;
+        }
+    }
+    else if (control->preserve_every) {
+        control->every_counter++;
+        if (control->every_counter == control->preserve_every) {
+            control->every_counter = 0;
+        }
+        else {
+            control->n_count_filtered++;
+            return TRUE;
+        }
+    }
+    return FALSE;
+}
+
+
+int count_decimation_is_end(struct CountDecimationControl *control)
+{
+    if (control->limit_n) {
+        control->limit_n_counter++;
+        /* this matches the last successfully imported point */
+        if (control->limit_n_counter == control->limit_n)
+            return TRUE;
+    }
+    return FALSE;
+}

vector/v.decimate/count_decimation.h

@@ -0,0 +1,43 @@
+/****************************************************************************
+ *
+ * MODULE:       v.decimate
+ * AUTHOR(S):    Vaclav Petras
+ * PURPOSE:      Reduce the number of points in a vector map
+ * COPYRIGHT:    (C) 2015 by the GRASS Development Team
+ *
+ *               This program is free software under the GNU General Public
+ *               License (>=v2). Read the COPYING file that comes with GRASS
+ *               for details.
+ *
+ *****************************************************************************/
+
+
+#ifndef GRASS_COUNT_DECIMATION_H
+#define GRASS_COUNT_DECIMATION_H
+
+/* TODO: change int to ul/ull */
+/* TODO: revise names (now partially on some vars in v.in.lidar code) */
+
+struct CountDecimationControl {
+    int offset_n;
+    int offset_n_counter;
+    int skip_every;
+    int preserve_every;
+    int every_counter;
+    int n_count_filtered;
+    int limit_n;
+    int limit_n_counter;
+};
+
+void count_decimation_init(struct CountDecimationControl *control,
+                           int *skip, int *preserve,
+                           int *offset, int *limit);
+int count_decimation_is_valid(struct CountDecimationControl *control);
+int count_decimation_is_noop(struct CountDecimationControl *control);
+void count_decimation_init_from_str(struct CountDecimationControl *control,
+                                    const char *skip, const char *preserve,
+                                    const char *offset, const char *limit);
+int count_decimation_is_out(struct CountDecimationControl *control);
+int count_decimation_is_end(struct CountDecimationControl *control);
+
+#endif /* GRASS_COUNT_DECIMATION_H */

vector/v.decimate/grid_decimation.c

@@ -0,0 +1,179 @@
+/****************************************************************************
+ *
+ * MODULE:       v.decimate
+ * AUTHOR(S):    Vaclav Petras
+ * PURPOSE:      Reduce the number of points in a vector map
+ * COPYRIGHT:    (C) 2015 by the GRASS Development Team
+ *
+ *               This program is free software under the GNU General Public
+ *               License (>=v2). Read the COPYING file that comes with GRASS
+ *               for details.
+ *
+ *****************************************************************************/
+
+
+#include "grid_decimation.h"
+
+#include <stdlib.h>
+
+
+/* max size: rows * cols < max of size_t (using 1D array) */
+void grid_decimation_create(struct GridDecimation *grid_decimation,
+                            size_t rows, size_t cols)
+{
+    grid_decimation->grid_points =
+        G_calloc(rows * cols, sizeof(struct DecimationPoint *));
+    grid_decimation->grid_sizes = G_calloc(rows * cols, sizeof(size_t));
+    grid_decimation->rows = rows;
+    grid_decimation->cols = cols;
+    grid_decimation->if_add_point = NULL;
+    grid_decimation->on_add_point = NULL;
+    grid_decimation->if_context = NULL;
+    grid_decimation->on_context = NULL;
+}
+
+
+void grid_decimation_destroy(struct GridDecimation *grid_decimation)
+{
+    /* TODO: we could also offer mode without dealloc (faster) */
+    int row, col;
+    size_t point, npoints;
+    for (row = 0; row < grid_decimation->rows; row++) {
+        for (col = 0; col < grid_decimation->cols; col++) {
+            /* TODO: make index function */
+            size_t index = row * grid_decimation->cols + col;
+            if ((npoints = grid_decimation->grid_sizes[index])) {
+                /* delete points in list */
+                for (point = 0; point < npoints; point++)
+                    G_free(grid_decimation->grid_points[index][point]);
+                /* delete list */
+                G_free(grid_decimation->grid_points[index]);
+            }
+        }
+    }
+    G_free(grid_decimation->grid_points);
+    G_free(grid_decimation->grid_sizes);
+}
+
+
+/* TODO: use Cell_head as storage? */
+void grid_decimation_create_from_region(struct GridDecimation
+                                        *grid_decimation,
+                                        struct Cell_head *region)
+{
+    grid_decimation_create(grid_decimation, region->rows, region->cols);
+    grid_decimation_set_region(grid_decimation, region->west, region->east,
+                               region->south, region->north, region->ew_res,
+                               region->ns_res);
+}
+
+
+/* TODO: change order of ns_res and ew_res to match xy */
+void grid_decimation_set_region(struct GridDecimation *grid_decimation,
+                                double minx, double maxx, double miny,
+                                double maxy, double ew_res, double ns_res)
+{
+    grid_decimation->minx = minx;
+    grid_decimation->maxx = maxx;
+    grid_decimation->miny = miny;
+    grid_decimation->maxy = maxy;
+    grid_decimation->ns_res = ns_res;
+    grid_decimation->ew_res = ew_res;
+}
+
+
+void grid_decimation_create_list_with_point(struct GridDecimation
+                                            *grid_decimation, size_t index,
+                                            struct DecimationPoint *point,
+                                            size_t npoints)
+{
+    struct DecimationPoint **point_list =
+        G_malloc(1 * sizeof(struct DecimationPoint *));
+    point_list[0] = point;
+    grid_decimation->grid_points[index] = point_list;
+    grid_decimation->grid_sizes[index] = 1;
+}
+
+
+void grid_decimation_add_point_to_list(struct GridDecimation *grid_decimation,
+                                       size_t index,
+                                       struct DecimationPoint *point,
+                                       size_t npoints)
+{
+    /* TODO: this might be too much reallocation */
+    /* TODO: line_ptns struct could be reused, it is not meant for this but it would work */
+    struct DecimationPoint **point_list =
+        G_realloc(grid_decimation->grid_points[index],
+                  (npoints + 1) * sizeof(struct DecimationPoint *));
+
+    point_list[npoints] = point;
+    grid_decimation->grid_points[index] = point_list;
+    grid_decimation->grid_sizes[index] = npoints + 1;
+}
+
+
+static size_t grid_decimation_xy_to_index(struct GridDecimation
+                                          *grid_decimation, double x,
+                                          double y)
+{
+    /* TODO: test x, y */
+    int row = (y - grid_decimation->miny) / grid_decimation->ns_res;
+    int col = (x - grid_decimation->minx) / grid_decimation->ew_res;
+
+    if (row < 0 || row > grid_decimation->rows || col < 0 ||
+        col > grid_decimation->cols) {
+        G_fatal_error
+            ("Row (%d) or column (%d) outside of range (0 - %d, 0 - %d)", row,
+             col, grid_decimation->rows, grid_decimation->cols);
+    }
+    size_t index = row * grid_decimation->cols + col;
+
+    /* TODO: are the tests really needed, especially the second one? */
+    if (row * grid_decimation->cols + col >
+        grid_decimation->rows * grid_decimation->cols) {
+        G_fatal_error("Index (%d) out of range (max: %d)",
+                      row * grid_decimation->cols + col,
+                      grid_decimation->rows * grid_decimation->cols);
+    }
+    return index;
+}
+
+
+void grid_decimation_try_add_point(struct GridDecimation *grid_decimation,
+                                   int cat, double x, double y, double z, void *point_data)
+{
+    size_t index = grid_decimation_xy_to_index(grid_decimation, x, y);
+    int npoints = grid_decimation->grid_sizes[index];
+
+    /* TODO: when max is 1, we don't have to store the point at all */
+    if (grid_decimation->max_points && grid_decimation->max_points == npoints)
+        return;
+
+    struct DecimationPoint *point = G_malloc(sizeof(struct DecimationPoint));
+
+    point->cat = cat;
+    point->x = x;
+    point->y = y;
+    point->z = z;
+
+    if (!npoints) {
+        grid_decimation_create_list_with_point(grid_decimation, index, point,
+                                               npoints);
+        if (grid_decimation->on_add_point)
+            grid_decimation->on_add_point(point, point_data, grid_decimation->on_context);
+    }
+    else {
+        if (grid_decimation->if_add_point
+            (point, point_data, grid_decimation->grid_points[index], npoints,
+             grid_decimation->if_context)) {
+            grid_decimation_add_point_to_list(grid_decimation, index, point,
+                                              npoints);
+            if (grid_decimation->on_add_point)
+                grid_decimation->on_add_point(point, point_data,
+                                              grid_decimation->on_context);
+        }
+        else {
+            G_free(point);
+        }
+    }
+}

vector/v.decimate/grid_decimation.h

@@ -0,0 +1,79 @@
+/****************************************************************************
+ *
+ * MODULE:       v.decimate
+ * AUTHOR(S):    Vaclav Petras
+ * PURPOSE:      Reduce the number of points in a vector map
+ * COPYRIGHT:    (C) 2015 by the GRASS Development Team
+ *
+ *               This program is free software under the GNU General Public
+ *               License (>=v2). Read the COPYING file that comes with GRASS
+ *               for details.
+ *
+ *****************************************************************************/
+
+
+#ifndef GRASS_GRID_DECIMATION_H
+#define GRASS_GRID_DECIMATION_H
+
+#include <grass/gis.h>
+
+struct DecimationPoint
+{
+    int cat;
+    double x;
+    double y;
+    double z;
+};
+
+
+struct GridDecimation
+{
+    struct DecimationPoint ***grid_points;
+    size_t *grid_sizes;
+    int rows;
+    int cols;
+    int max_points;
+    double minx;
+    double maxx;
+    double miny;
+    double maxy;
+    double ns_res;
+    double ew_res;
+    int (*if_add_point) (struct DecimationPoint * point, void *point_data,
+                         struct DecimationPoint ** point_list, size_t npoints,
+                         void *context);
+    void (*on_add_point) (struct DecimationPoint * point, void *point_data, void *context);
+    void *if_context;
+    void *on_context;
+};
+
+
+/* max size: rows * cols < max of size_t (using 1D array) */
+void grid_decimation_create(struct GridDecimation *grid_decimation,
+                            size_t rows, size_t cols);
+
+void grid_decimation_create_from_region(struct GridDecimation
+                                        *grid_decimation,
+                                        struct Cell_head *region);
+
+/* TODO: possible use (also?) Cell_head as input or storage */
+void grid_decimation_set_region(struct GridDecimation *grid_decimation,
+                                double minx, double maxx, double miny,
+                                double maxy, double ns_res, double ew_res);
+
+void grid_decimation_create_list_with_point(struct GridDecimation
+                                            *grid_decimation, size_t index,
+                                            struct DecimationPoint *point,
+                                            size_t npoints);
+
+void grid_decimation_add_point_to_list(struct GridDecimation *grid_decimation,
+                                       size_t index,
+                                       struct DecimationPoint *point,
+                                       size_t npoints);
+
+void grid_decimation_try_add_point(struct GridDecimation *grid_decimation,
+                                   int cat, double x, double y, double z, void *point_data);
+
+void grid_decimation_destroy(struct GridDecimation *grid_decimation);
+
+#endif /* GRASS_GRID_DECIMATION_H */

vector/v.decimate/main.c

@@ -0,0 +1,465 @@
+/*
+ ****************************************************************************
+ *
+ * MODULE:       v.decimate
+ * AUTHOR(S):    Vaclav Petras
+ * PURPOSE:      Reduce the number of points in a vector map
+ * COPYRIGHT:    (C) 2015 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.
+ *
+ *****************************************************************************/
+
+
+/* using the most-specific-first order of includes */
+#include "count_decimation.h"
+#include "grid_decimation.h"
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+
+#include <stdlib.h>
+
+
+struct DecimationContext
+{
+    int use_z;                  /*!< TRUE or FALSE */
+    double zdiff;
+    int unique_cats;            /*!< TRUE or FALSE */
+};
+
+
+static int if_add_point(struct DecimationPoint *point, void *point_data,
+                        struct DecimationPoint **point_list, size_t npoints,
+                        void *context)
+{
+    /* according to cat (which could be cluster, return or class) */
+    struct DecimationContext *dc = context;
+    double zdiff = dc->zdiff;
+    int j;
+
+    /* TODO: use something like Vect_cat_in_cat_list? */
+    for (j = 0; j < npoints; j++) {
+        if (dc->use_z && fabs(point_list[j]->z - point->z) < zdiff)
+            return FALSE;
+        if (dc->unique_cats && point_list[j]->cat == point->cat)
+            return FALSE;
+    }
+    return TRUE;
+}
+
+
+struct WriteContext
+{
+    struct Map_info *voutput;
+    struct line_pnts *line;
+    struct line_cats *cats;
+    int write_cats;
+};
+
+
+static void write_point(struct WriteContext *context, int cat, double x,
+                        double y, double z, struct line_cats *cats)
+{
+    if (Vect_append_point(context->line, x, y, z) != 1)
+        G_fatal_error
+            ("Unable to create a point in vector map (probably out of memory)");
+    struct line_cats *cats_to_write = context->cats;
+    /* only when writing cats use the ones from parameter, otherwise
+     * use the default (which is assumed to be empty) */
+    if (context->write_cats && cats)
+        cats_to_write = cats;
+    Vect_write_line(context->voutput, GV_POINT, context->line, cats_to_write);
+    Vect_reset_line(context->line);
+}
+
+
+static void on_add_point(struct DecimationPoint *point, void *point_data, void *context)
+{
+    write_point((struct WriteContext *)context, point->cat, point->x, point->y,
+                point->z, (struct line_cats *)point_data);
+}
+
+/* TODO: these have overlap with vector lib, really needed? */
+static int point_in_region_2d(struct Cell_head *region, double x, double y)
+{
+    if (x > region->east || x < region->west || y < region->south ||
+        y > region->north)
+        return FALSE;
+    return TRUE;
+}
+
+
+static int point_in_region_3d(struct Cell_head *region, double x, double y,
+                              double z)
+{
+    if (x > region->east || x < region->west || y < region->south ||
+        y > region->north || z > region->top || z < region->bottom)
+        return FALSE;
+    return TRUE;
+}
+
+/* TODO: max tries per grid cell (useful?) */
+
+
+int main(int argc, char **argv)
+{
+    struct GModule *module;
+    struct Option *map_opt, *voutput_opt;
+    struct Option *field_opt, *cats_opt;
+    struct Option *skip_opt, *preserve_opt, *offset_opt, *limit_opt;
+    struct Option *zdiff_opt, *limit_per_cell_opt, *zrange_opt;
+    struct Flag *grid_decimation_flg, *first_point_flg, *cat_in_grid_flg;
+    struct Flag *use_z_flg, *nocats_flag, *notopo_flag;
+    struct Map_info vinput, voutput;
+
+    G_gisinit(argv[0]);
+
+    module = G_define_module();
+    G_add_keyword(_("vector"));
+    G_add_keyword(_("LIDAR"));
+    G_add_keyword(_("generalization"));
+    G_add_keyword(_("decimation"));
+    G_add_keyword(_("extract"));
+    G_add_keyword(_("select"));
+    G_add_keyword(_("points"));
+    module->label = _("Decimates a point cloud");
+    module->description = _("Copies points from one vector to another"
+                            " while applying different decimations");
+
+    map_opt = G_define_standard_option(G_OPT_V_INPUT);
+
+    field_opt = G_define_standard_option(G_OPT_V_FIELD_ALL);
+    field_opt->required = NO;
+
+    voutput_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+
+    zrange_opt = G_define_option();
+    zrange_opt->key = "zrange";
+    zrange_opt->type = TYPE_DOUBLE;
+    zrange_opt->required = NO;
+    zrange_opt->key_desc = "min,max";
+    zrange_opt->description = _("Filter range for z data (min,max)");
+    zrange_opt->guisection = _("Selection");
+
+    cats_opt = G_define_standard_option(G_OPT_V_CATS);
+    cats_opt->guisection = _("Selection");
+
+    /* TODO: -r region only (now enforced for grid), spatial */
+
+    skip_opt = G_define_option();
+    skip_opt->key = "skip";
+    skip_opt->type = TYPE_INTEGER;
+    skip_opt->multiple = NO;
+    skip_opt->required = NO;
+    skip_opt->label = _("Throw away every n-th point");
+    skip_opt->description =
+        _("For example, 5 will import 80 percent of points. "
+          "If not specified, all points are copied");
+    skip_opt->guisection = _("Count");
+
+    preserve_opt = G_define_option();
+    preserve_opt->key = "preserve";
+    preserve_opt->type = TYPE_INTEGER;
+    preserve_opt->multiple = NO;
+    preserve_opt->required = NO;
+    preserve_opt->label = _("Preserve only every n-th point");
+    preserve_opt->description =
+        _("For example, 4 will import 25 percent of points. "
+          "If not specified, all points are copied");
+    preserve_opt->guisection = _("Count");
+
+    offset_opt = G_define_option();
+    offset_opt->key = "offset";
+    offset_opt->type = TYPE_INTEGER;
+    offset_opt->multiple = NO;
+    offset_opt->required = NO;
+    offset_opt->label = _("Skip first n points");
+    offset_opt->description =
+        _("Skips the given number of points at the beginning.");
+    offset_opt->guisection = _("Count");
+
+    limit_opt = G_define_option();
+    limit_opt->key = "limit";
+    limit_opt->type = TYPE_INTEGER;
+    limit_opt->multiple = NO;
+    limit_opt->required = NO;
+    limit_opt->label = _("Copy only n points");
+    limit_opt->description = _("Copies only the given number of points");
+    limit_opt->guisection = _("Count");
+
+    zdiff_opt = G_define_option();
+    zdiff_opt->key = "zdiff";
+    zdiff_opt->type = TYPE_DOUBLE;
+    zdiff_opt->required = NO;
+    zdiff_opt->label = _("Minimal difference of z values");
+    zdiff_opt->description =
+        _("Minimal difference between z values in grid-based decimation");
+    zdiff_opt->guisection = _("Grid");
+
+    limit_per_cell_opt = G_define_option();
+    limit_per_cell_opt->key = "cell_limit";
+    limit_per_cell_opt->type = TYPE_INTEGER;
+    limit_per_cell_opt->multiple = NO;
+    limit_per_cell_opt->required = NO;
+    limit_per_cell_opt->label = _("Preserve only n points per grid cell");
+    limit_per_cell_opt->description =
+        _("Preserves only the given number of points per grid cell in grid-based decimation");
+    limit_per_cell_opt->guisection = _("Grid");
+
+    grid_decimation_flg = G_define_flag();
+    grid_decimation_flg->key = 'g';
+    grid_decimation_flg->description = _("Apply grid-based decimation");
+    grid_decimation_flg->guisection = _("Grid");
+
+    first_point_flg = G_define_flag();
+    first_point_flg->key = 'f';
+    first_point_flg->description =
+        _("Use only first point in grid cell during grid-based decimation");
+    first_point_flg->guisection = _("Grid");
+
+    cat_in_grid_flg = G_define_flag();
+    cat_in_grid_flg->key = 'c';
+    cat_in_grid_flg->description = _("Only one point per cat in grid cell");
+    cat_in_grid_flg->guisection = _("Grid");
+
+    use_z_flg = G_define_flag();
+    use_z_flg->key = 'z';
+    use_z_flg->description = _("Use z in grid decimation");
+    use_z_flg->guisection = _("Grid");
+
+    nocats_flag = G_define_flag();
+    nocats_flag->key = 'x';
+    nocats_flag->label =
+        _("Store only the coordinates, throw away categories");
+    nocats_flag->description =
+        _("Do not story any categories even if they are present in input data");
+    nocats_flag->guisection = _("Speed");
+
+    notopo_flag = G_define_standard_flag(G_FLG_V_TOPO);
+    notopo_flag->guisection = _("Speed");
+
+    /* here we have different decimations but also selections/filters */
+    G_option_required(skip_opt, preserve_opt, offset_opt, limit_opt,
+                      grid_decimation_flg, zrange_opt, cats_opt, NULL);
+    /* this doesn't play well with GUI dialog unless we add defaults to options
+     * the default values would solve it but looks strange in the manual 
+     * this we use explicit check in the code */
+    /* G_option_exclusive(skip_opt, preserve_opt, NULL); */
+    G_option_requires(grid_decimation_flg, first_point_flg,
+                      limit_per_cell_opt, use_z_flg, zdiff_opt,
+                      cat_in_grid_flg, NULL);
+    G_option_requires(first_point_flg, grid_decimation_flg, NULL);
+    G_option_requires(limit_per_cell_opt, grid_decimation_flg, NULL);
+    G_option_requires(use_z_flg, grid_decimation_flg, NULL);
+    G_option_requires(zdiff_opt, grid_decimation_flg, NULL);
+    G_option_requires(cat_in_grid_flg, grid_decimation_flg, NULL);
+    G_option_exclusive(zdiff_opt, first_point_flg, limit_per_cell_opt, NULL);
+
+    if (G_parser(argc, argv))
+        exit(EXIT_FAILURE);
+
+    Vect_check_input_output_name(map_opt->answer, voutput_opt->answer,
+                                 G_FATAL_EXIT);
+
+    if (Vect_open_old2(&vinput, map_opt->answer, "", field_opt->answer) < 0)
+        G_fatal_error(_("Unable to open vector map <%s>"), map_opt->answer);
+    int layer = Vect_get_field_number(&vinput, field_opt->answer);
+
+    if (layer < 1 && (cats_opt->answer || cat_in_grid_flg->answer))
+        G_fatal_error(_("Input layer must be set to a particular layer"
+                        ", not <%s>, when using <%s> option or <-%c> flag"),
+                      field_opt->answer, cats_opt->key, cat_in_grid_flg->key);
+
+    struct cat_list *allowed_cats = NULL;
+
+    if (layer > 0)
+        allowed_cats = Vect_cats_set_constraint(&vinput, layer, NULL,
+                                                cats_opt->answer);
+
+    struct line_pnts *line = Vect_new_line_struct();
+    struct line_cats *cats = Vect_new_cats_struct();
+
+    double zrange_min, zrange_max;
+    int use_zrange = FALSE;
+
+    if (zrange_opt->answer != NULL) {
+        if (zrange_opt->answers[0] == NULL || zrange_opt->answers[1] == NULL)
+            G_fatal_error(_("Invalid zrange <%s>"), zrange_opt->answer);
+        sscanf(zrange_opt->answers[0], "%lf", &zrange_min);
+        sscanf(zrange_opt->answers[1], "%lf", &zrange_max);
+        /* for convenience, switch order to make valid input */
+        if (zrange_min > zrange_max) {
+            double tmp = zrange_max;
+
+            zrange_max = zrange_min;
+            zrange_min = tmp;
+        }
+        use_zrange = TRUE;
+    }
+
+    int use_z = FALSE;
+    double zdiff;
+
+    if (use_z_flg->answer) {
+        use_z = TRUE;
+        zdiff = atof(zdiff_opt->answer);
+    }
+
+    /* z input checks */
+    if (!Vect_is_3d(&vinput)) {
+        if (use_z)
+            G_fatal_error(_("Cannot use z for decimation, input is not 3D"));
+        if (use_zrange)
+            G_fatal_error(_("Cannot select by z range, input is not 3D"));
+    }
+
+    int do_grid_decimation = FALSE;
+
+    if (grid_decimation_flg->answer)
+        do_grid_decimation = TRUE;
+    int limit_per_cell = 0;
+
+    if (limit_per_cell_opt->answer)
+        limit_per_cell = atof(limit_per_cell_opt->answer);
+    if (first_point_flg->answer)
+        limit_per_cell = 1;
+
+    /* when user enters 1 or zero e.g. for skip, we accept it and use it
+     * but the is no advantage, however, we count it as an error when
+     * no other options are selected
+     */
+    struct CountDecimationControl count_decimation_control;
+
+    count_decimation_init_from_str(&count_decimation_control,
+                                   skip_opt->answer, preserve_opt->answer,
+                                   offset_opt->answer, limit_opt->answer);
+    if (!count_decimation_is_valid(&count_decimation_control))
+        G_fatal_error(_("Settings for count-based decimation are not valid"));
+    /* TODO: implement count_decimation_is_invalid_reason() */
+    /* the following must be in sync with required options */
+    if (count_decimation_is_noop(&count_decimation_control) &&
+        !grid_decimation_flg->answer && !zrange_opt->answer &&
+        !cats_opt->answer)
+        G_fatal_error(_("Settings for count-based decimation would cause it"
+                        " to do nothing and no other options has been set."));
+
+    struct Cell_head comp_region;
+    Rast_get_window(&comp_region);
+    if (use_zrange) {
+        comp_region.bottom = zrange_min;
+        comp_region.top = zrange_max;
+    }
+    struct GridDecimation grid_decimation;
+    struct DecimationContext decimation_context;
+    struct WriteContext write_context;
+
+    write_context.line = Vect_new_line_struct();
+    write_context.cats = Vect_new_cats_struct();
+    if (!nocats_flag->answer)
+        write_context.write_cats = TRUE;
+    else
+        write_context.write_cats = FALSE;
+    write_context.voutput = &voutput;
+    if (do_grid_decimation) {
+        grid_decimation_create_from_region(&grid_decimation, &comp_region);
+        grid_decimation.max_points = limit_per_cell;
+
+        if (cat_in_grid_flg->answer)
+            decimation_context.unique_cats = TRUE;
+        else
+            decimation_context.unique_cats = FALSE;
+        decimation_context.use_z = FALSE;
+        if (use_z) {
+            decimation_context.use_z = TRUE;
+            decimation_context.zdiff = zdiff;
+        }
+        grid_decimation.if_add_point = if_add_point;
+        grid_decimation.if_context = &decimation_context;
+
+        grid_decimation.on_add_point = on_add_point;
+        grid_decimation.on_context = &write_context;
+    }
+
+    if (Vect_open_new(&voutput, voutput_opt->answer, Vect_is_3d(&vinput)) < 0)
+        G_fatal_error(_("Unable to create vector map <%s>"),
+                      voutput_opt->answer);
+
+    /* some constraints can be set on the map */
+    Vect_set_constraint_type(&vinput, GV_POINT);
+    /* noop for layer=-1 and non-native format, skips lines without cat */
+    Vect_set_constraint_field(&vinput, layer);
+    /* TODO: replace region checks by Vect_set_constraint_region? */
+
+    int ltype;
+    int cat;
+
+    while (TRUE) {
+        ltype = Vect_read_next_line(&vinput, line, cats);
+        if (ltype == -1)
+            G_fatal_error(_("Unable to read vector map"));
+        if (ltype == -2)
+            break;              /* end of the map */
+
+        double x, y, z;
+        Vect_line_get_point(line, 0, &x, &y, &z);
+
+        /* selections/filters */
+        /* TODO: use region only when actually needed */
+        if (!use_zrange && !point_in_region_2d(&comp_region, x, y))
+            continue;
+        /* TODO: allow zrange to be used without region */
+        if (use_zrange && !point_in_region_3d(&comp_region, x, y, z))
+            continue;
+        if (allowed_cats &&
+            !Vect_cats_in_constraint(cats, layer, allowed_cats))
+            continue;
+
+        /* decimation */
+        if (count_decimation_is_out(&count_decimation_control))
+            continue;
+
+        /* TODO: test: skip points without category, unless layer=-1 */
+        /* Use cases:
+         * - all points have category (correct)
+         * - no categories for any point (correct, layer=-1 required)
+         * - some points miss category (not handled)
+         * Here we assume that only one cat has meaning for grid decimation.
+         * If no layer available, cat contains junk and shouldn't be used.
+         */
+        if (layer > 0)
+            Vect_cat_get(cats, layer, &cat);
+
+        /* using callback when using grid, direct call otherwise */
+        if (do_grid_decimation)
+            grid_decimation_try_add_point(&grid_decimation, cat, x, y, z,
+                                          cats);
+        else
+            write_point(&write_context, cat, x, y, z, cats);
+
+        if (count_decimation_is_end(&count_decimation_control))
+            break;
+    }
+
+    /* partially unnecessary as deallocated by the system */
+    Vect_destroy_line_struct(line);
+    Vect_destroy_cats_struct(cats);
+    Vect_destroy_line_struct(write_context.line);
+    Vect_destroy_cats_struct(write_context.cats);
+
+    Vect_hist_command(&voutput);
+
+    Vect_close(&vinput);
+    if (!notopo_flag->answer)
+        Vect_build(&voutput);
+    Vect_close(&voutput);
+
+    if (do_grid_decimation)
+        grid_decimation_destroy(&grid_decimation);
+
+    return EXIT_SUCCESS;
+}

vector/v.decimate/v.decimate.html

@@ -0,0 +1,138 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.decimate</em> reduces number of points in the input vector map
+and copies them over to the output vector map. Different point decimation
+techniques can be applied to reduce the number of points.
+
+<p>
+Two main decimation techniques are:
+<ul>
+<li>count-based decimation (<b>skip</b>, <b>preserve</b>, <b>offset</b>
+    and <b>limit</b> options)
+<li>grid-based decimation (<b>-g</b> flag)
+</ul>
+
+<p>
+The grid-based decimation will remove points based on:
+<ul>
+<li>similar z coordinates (<b>-z</b> flag and <b>zdiff</b> option)
+<li>same categories (<b>-c</b> flag)
+<li>count of points (<b>-f</b> flag and <b>cell_limit</b> option)
+</ul>
+
+<p>
+The grid-based decimation is currently using a 2D grid, so the points
+are placed and compared within this 2D grid. The comparison can happen
+using z coordinates or categories.
+Note that although the grid is only 2D, the module works with 3D points.
+
+<p>
+The grid-based decimation extent and resolution depend on the current
+computational region as set by <em><a href="g.region.html">g.region</a></em>.
+As a consequence, the output is limited only to computational region
+in this case.
+
+<p>
+TODO: Currently, any output is limited by the region.
+
+<p>
+The count-based decimation result highly depends on how the data are
+ordered in the input. This applies especially to <b>offset</b> and
+<b>limit</b> options where the resulting shape and densities can be
+surprising. The options <b>skip</b> and <b>preserve</b> are influenced
+by order of points in a similar way but they usually keep relative
+density of points (which may or may not be desired).
+On the other hand, the grid-based decimation will generally result in
+more even density of output points (see Figure 1).
+
+<p>
+Besides decimation, point count can be reduced by applying different
+selections or filters, these are:
+<ul>
+<li>selection by category (<b>cats</b> option)
+<li>selection by z values (<b>zrange</b> option)
+</ul>
+
+
+<h2>NOTES</h2>
+
+The grid-based decimation requires all points which will be saved in output
+to fit into the computer's memory (RAM).
+It is advantageous to have the region only in the area
+with the points, otherwise unnecessary memory is allocated.
+Higher (finer) resolutions and higher amount of preserved points
+per cell require more memory.
+The count-based decimation has no limitation regarding the available memory.
+
+<p>
+Significant speed up can be gained using <b>-b</b> flag which disables
+building of topology for the output vector map. This may limit the use
+of the vector map by some modules, but for example, this module works
+without topology as well.
+
+
+<h2>EXAMPLES</h2>
+
+Keep only every forth point, throw away the rest:
+
+<div class="code"><pre>
+v.decimate input=points_all output=points_decimated_every_4 preserve=4
+</pre></div>
+
+<p>
+Keep only points within a grid cell (given by the current computational
+region) which has unique categories (e.g. LIDAR classes):
+
+<div class="code"><pre>
+v.decimate input=points_all output=points_decimated_unique_cats layer=1 -g -c
+</pre></div>
+
+<center>
+    <img src="v_decimate_original.png">
+    <img src="v_decimate_count.png">
+    <img src="v_decimate_grid_cat.png">
+    <p><em>
+        Figure 1: Comparison of original points, decimation result
+        with every forth point preserved, and grid-based decimation
+        result with points with unique categories in each grid cell
+    </em></p>
+</center>
+
+<p>
+Keep only points with category 2 and keep only approximately 80% of the points:
+
+<div class="code"><pre>
+v.decimate input=points_all output=points_decimated_ skip=5 cats=2 layer=1
+</pre></div>
+
+<!--
+g.region map=elev_lid792_1m
+v.colors map=points_all use=cat rules=file.txt
+2 brown
+3:5 green
+6 red
+9 blue
+11 black
+13-14 gray
+d.vect map=points color=none width=1 icon=basic/point
+-->
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="v.extract.html">v.extract</a>,
+<a href="v.outlier.html">v.outlier</a>,
+<a href="v.select.html">v.select</a>,
+<a href="v.category.html">v.category</a>,
+<a href="v.build.html">v.build</a>,
+<a href="v.in.lidar.html">v.in.lidar</a>,
+<a href="g.region.html">g.region</a>
+</em>
+
+
+<h2>AUTHORS</h2>
+
+Vaclav Petras, <a href="http://geospatial.ncsu.edu/osgeorel/">NCSU OSGeoREL</a>
+
+<p><i>Last changed: $Date$</i>