[GRASS-SVN] r68789 - in sandbox: . moritz moritz/r.object.geometry
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jun 28 01:40:38 PDT 2016
Author: mlennert
Date: 2016-06-28 01:40:38 -0700 (Tue, 28 Jun 2016)
New Revision: 68789
Added:
sandbox/moritz/
sandbox/moritz/r.object.geometry/
sandbox/moritz/r.object.geometry/Makefile
sandbox/moritz/r.object.geometry/main.c
sandbox/moritz/r.object.geometry/r.object.geometry.html
Log:
New module: r.object.geometry
Added: sandbox/moritz/r.object.geometry/Makefile
===================================================================
--- sandbox/moritz/r.object.geometry/Makefile (rev 0)
+++ sandbox/moritz/r.object.geometry/Makefile 2016-06-28 08:40:38 UTC (rev 68789)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.object.geometry
+
+LIBES = $(RASTERLIB) $(GISLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: sandbox/moritz/r.object.geometry/main.c
===================================================================
--- sandbox/moritz/r.object.geometry/main.c (rev 0)
+++ sandbox/moritz/r.object.geometry/main.c 2016-06-28 08:40:38 UTC (rev 68789)
@@ -0,0 +1,316 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.object.geometry
+ *
+ * AUTHOR(S): Moritz Lennert
+ * Markus Metz
+ *
+ * PURPOSE: Fetch object geometry parameters.
+ *
+ * COPYRIGHT: (C) 2016 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ ***************************************************************************/
+
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+/* compare two cell values
+ * return 0 if equal, 1 if different */
+static int cmp_cells(CELL a, CELL b, int a_null, int b_null)
+{
+ return (a_null + b_null == 1 || (a_null + b_null == 0 && a != b));
+}
+
+int main(int argc, char *argv[])
+{
+ int row, col, nrows, ncols;
+
+ struct Range range;
+ CELL min, max;
+ int in_fd;
+ int i;
+ struct GModule *module;
+ struct Option *opt_in;
+ struct Option *opt_out;
+ struct Option *opt_sep;
+ struct Flag *flag_m;
+ char *sep;
+ FILE *out_fp;
+ CELL *prev_in, *cur_in, *temp_in;
+ CELL cur, top, left;
+ int cur_null, top_null, left_null;
+ int len;
+ struct obj_geo
+ {
+ double area, perimeter;
+ int min_row, max_row, min_col, max_col; /* bounding box */
+ } *obj_geos;
+ double unit_area;
+ int n_objects;
+ int planimetric, compute_areas;
+ struct Cell_head cellhd;
+ double ew_factor;
+
+ G_gisinit(argv[0]);
+
+ /* Define the different options */
+
+ module = G_define_module();
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("statistics"));
+ G_add_keyword(_("reclass"));
+ G_add_keyword(_("clumps"));
+ module->description =
+ _("Calculates geometry parameters for raster objects.");
+
+ opt_in = G_define_standard_option(G_OPT_R_INPUT);
+
+ opt_out = G_define_standard_option(G_OPT_F_OUTPUT);
+ opt_out->required = NO;
+
+ opt_sep = G_define_standard_option(G_OPT_F_SEP);
+
+ flag_m = G_define_flag();
+ flag_m->key = 'm';
+ flag_m->label = _("Use meters as units instead of cells");
+
+ /* parse options */
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ sep = G_option_to_separator(opt_sep);
+ in_fd = Rast_open_old(opt_in->answer, "");
+
+ if (Rast_get_map_type(in_fd) != CELL_TYPE)
+ G_fatal_error(_("Input raster mus be of type CELL"));
+
+ out_fp = stdout;
+ /* TODO: open optional output file */
+
+ Rast_get_window(&cellhd);
+ nrows = cellhd.rows;
+ ncols = cellhd.cols;
+ /* for unit = cells */
+ ew_factor = cellhd.ew_res / cellhd.ns_res;
+
+ /* allocate CELL buffers two columns larger than current window */
+ len = (ncols + 2) * sizeof(CELL);
+ prev_in = (CELL *) G_malloc(len);
+ cur_in = (CELL *) G_malloc(len);
+
+ /* fake a previous row which is all NULL */
+ Rast_set_c_null_value(prev_in, ncols + 2);
+
+ /* set left and right edge to NULL */
+ Rast_set_c_null_value(&cur_in[0], 1);
+ Rast_set_c_null_value(&cur_in[ncols + 1], 1);
+
+ Rast_read_range(opt_in->answer, "", &range);
+ Rast_get_range_min_max(&range, &min, &max);
+ if (Rast_is_c_null_value(&min) || Rast_is_c_null_value(&max))
+ G_fatal_error(_("Empty input map <%s>"), opt_in->answer);
+
+ n_objects = max - min + 1;
+ obj_geos = G_malloc(n_objects * sizeof(struct obj_geo));
+
+ for (i = 0; i < n_objects; i++) {
+ obj_geos[i].area = 0;
+ obj_geos[i].perimeter = 0;
+ obj_geos[i].min_row = nrows;
+ obj_geos[i].max_row = 0;
+ obj_geos[i].min_col = ncols;
+ obj_geos[i].max_col = 0;
+ }
+
+ unit_area = 0.0;
+ if (flag_m->answer) {
+ switch (G_begin_cell_area_calculations()) {
+ case 0: /* areas don't make sense, but ignore this for now */
+ case 1:
+ planimetric = 1;
+ unit_area = G_area_of_cell_at_row(0);
+ break;
+ default:
+ planimetric = 0;
+ break;
+ }
+ }
+ compute_areas = flag_m->answer && !planimetric;
+
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows + 1, 2);
+
+ Rast_get_c_row(in_fd, cur_in + 1, row);
+ cur = cur_in[0];
+ cur_null = 1;
+ if (compute_areas)
+ unit_area = G_area_of_cell_at_row(row);
+ for (col = 1; col <= ncols; col++) {
+ left = cur;
+ cur = cur_in[col];
+ top = prev_in[col];
+
+ left_null = cur_null;
+ cur_null = Rast_is_c_null_value(&cur);
+ top_null = Rast_is_c_null_value(&top);
+
+ if (!cur_null) {
+ if (flag_m->answer) {
+ obj_geos[cur - min].area += unit_area;
+ } else {
+ obj_geos[cur - min].area += 1 * ew_factor;
+ }
+ if (obj_geos[cur - min].min_row > row)
+ obj_geos[cur - min].min_row = row;
+ if (obj_geos[cur - min].max_row < row + 1)
+ obj_geos[cur - min].max_row = row + 1;
+ if (obj_geos[cur - min].min_col > col)
+ obj_geos[cur - min].min_col = col;
+ if (obj_geos[cur - min].max_col < col + 1)
+ obj_geos[cur - min].max_col = col + 1;
+ }
+
+ /* TODO: convert perimeter to meters for -m flag */
+ if (cmp_cells(cur, top, cur_null, top_null)) {
+ if (flag_m->answer) {
+ double perimeter;
+
+ /* could be optimized */
+ perimeter = G_distance(cellhd.west + col * cellhd.ew_res,
+ Rast_row_to_northing(row, &cellhd),
+ cellhd.west + (col + 1) * cellhd.ew_res,
+ Rast_row_to_northing(row, &cellhd));
+
+ if (!cur_null)
+ obj_geos[cur - min].perimeter += perimeter;
+ if (!top_null)
+ obj_geos[top - min].perimeter += perimeter;
+ }
+ else {
+ if (!cur_null)
+ obj_geos[cur - min].perimeter += ew_factor;
+ if (!top_null)
+ obj_geos[top - min].perimeter += ew_factor;
+ }
+ }
+ if (cmp_cells(cur, left, cur_null, left_null)) {
+ if (flag_m->answer) {
+ double perimeter;
+
+ /* could be optimized */
+ perimeter = G_distance(cellhd.west + col * cellhd.ew_res,
+ Rast_row_to_northing(row, &cellhd),
+ cellhd.west + (col) * cellhd.ew_res,
+ Rast_row_to_northing(row + 1, &cellhd));
+
+ if (!cur_null)
+ obj_geos[cur - min].perimeter += perimeter;
+ if (!left_null)
+ obj_geos[left - min].perimeter += perimeter;
+ }
+ else {
+ if (!cur_null)
+ obj_geos[cur - min].perimeter += 1;
+ if (!left_null)
+ obj_geos[left - min].perimeter += 1;
+ }
+ }
+ }
+
+ /* switch the buffers so that the current buffer becomes the previous */
+ temp_in = cur_in;
+ cur_in = prev_in;
+ prev_in = temp_in;
+ }
+
+ /* last row, bottom borders */
+ G_percent(row, nrows + 1, 2);
+ for (col = 1; col <= ncols; col++) {
+ top = prev_in[col];
+ top_null = Rast_is_c_null_value(&top);
+
+ if (flag_m->answer) {
+ double perimeter;
+
+ /* could be optimized */
+ perimeter = G_distance(cellhd.west + col * cellhd.ew_res,
+ Rast_row_to_northing(row, &cellhd),
+ cellhd.west + (col + 1) * cellhd.ew_res,
+ Rast_row_to_northing(row, &cellhd));
+
+ if (!top_null)
+ obj_geos[top - min].perimeter += perimeter;
+ }
+ else {
+ if (!top_null)
+ obj_geos[top - min].perimeter += ew_factor;
+ }
+ }
+ G_percent(1, 1, 1);
+
+ Rast_close(in_fd);
+ G_free(cur_in);
+ G_free(prev_in);
+
+ /* print table */
+ /* TODO: use separator */
+ /* print table header (column names) */
+ fprintf(out_fp, "\n");
+
+ /* print table body */
+ for (i = 0; i < n_objects; i++) {
+ /* skip emtpy objects */
+ if (obj_geos[i].area == 0)
+ continue;
+
+ fprintf(out_fp, "%d%s", i, sep);
+ fprintf(out_fp, "%f%s", obj_geos[i].area, sep);
+ fprintf(out_fp, "%f", obj_geos[i].perimeter);
+ /* TODO: convert bounding box area and perimeter to (square) meters for -m flag */
+
+
+ /* object id: i + min */
+
+ /* area */
+
+ /* perimeter */
+
+ /* compactness */
+ /* area / perimeter -> more compact objects have higher compactness value
+ * most compact object in a grid is a square with side length r
+ * area = r^2, perimeter = 4 * r = 4 * sqrt(area)
+ * compactness = area / (4 * sqrt(area)) = sqrt(area) / 4
+ * -> compactness is dependent of area size
+ * -> normalize to square to make compactness of objects of
+ * different size comparable:
+ * compactness = (area / perimeter) / (sqrt(area) / 4)
+ * compactness = 4 * sqrt(area) / perimeter
+ * -> compactness is in the range 0 < compactness <= 1
+ */
+
+ /* fractal dimension */
+
+ /* smoothness */
+ /* perimeter of bounding box / perimeter -> smoother objects have a higher smoothness value
+ * smoothness is in the range 0 < smoothness <= 1 */
+
+ /* bounding box perimeter */
+
+ /* bounding box size */
+
+ fprintf(out_fp, "\n");
+ }
+ if (out_fp != stdout)
+ fclose(out_fp);
+
+ exit(EXIT_SUCCESS);
+}
Added: sandbox/moritz/r.object.geometry/r.object.geometry.html
===================================================================
--- sandbox/moritz/r.object.geometry/r.object.geometry.html (rev 0)
+++ sandbox/moritz/r.object.geometry/r.object.geometry.html 2016-06-28 08:40:38 UTC (rev 68789)
@@ -0,0 +1,25 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.object.geometry</em> does ...
+
+
+<h2>NOTES</h2>
+
+
+<h2>EXAMPLES</h2>
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="i.segment.html">i.segment</a>,
+<a href="r.clump.html">r.clump</a>,
+</em>
+
+<h2>AUTHOR</h2>
+
+Moritz Lennert<br>
+Markus Metz (diagonal clump tracing)
+
+<p>
+<i>Last changed: $Date: 2016-05-31 00:01:21 +0200 (Tue, 31 May 2016) $</i>
More information about the grass-commit
mailing list