/**************************************************************************** * * MODULE: r.object.geometry * * AUTHOR(S): Moritz Lennert * Markus Metz * * PURPOSE: Fetch object geometry parameters. * * COPYRIGHT: (C) 2016-2021 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 #include #include #include #include #include #include /* 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, x, y; int min_row, max_row, min_col, max_col; /* bounding box */ int num; } *obj_geos; double unit_area; int n_objects; int planimetric, compute_areas; struct Cell_head cellhd; 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 must be of type CELL")); if (opt_out->answer != NULL && strcmp(opt_out->answer, "-") != 0) { if (!(out_fp = fopen(opt_out->answer, "w"))) G_fatal_error(_("Unable to open file <%s> for writing"), opt_out->answer); } else { out_fp = stdout; } Rast_get_window(&cellhd); nrows = cellhd.rows; ncols = cellhd.cols; /* 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); /* REMARK: The following is only true if object ids are continuously numbered */ 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; obj_geos[i].x = 0; obj_geos[i].y = 0; obj_geos[i].num = 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; G_begin_distance_calculations(); G_message(_("Calculating statistics")); 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; obj_geos[cur - min].num += 1; } else { obj_geos[cur - min].area += 1; } obj_geos[cur - min].x += Rast_col_to_easting(col - 0.5, &cellhd); obj_geos[cur - min].y += Rast_row_to_northing(row + 0.5, &cellhd); 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; } 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 += 1; if (!top_null) obj_geos[top - min].perimeter += 1; } } 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; } } } /* last col, right borders */ if (flag_m->answer) { double perimeter; perimeter = G_distance(cellhd.east, Rast_row_to_northing(row, &cellhd), cellhd.east, Rast_row_to_northing(row + 1, &cellhd)); if (!cur_null) obj_geos[cur - min].perimeter += perimeter; } else { if (!cur_null) obj_geos[cur - 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 += 1; } } G_percent(1, 1, 1); Rast_close(in_fd); G_free(cur_in); G_free(prev_in); G_message(_("Writing output")); /* print table */ fprintf(out_fp, "cat%s", sep); fprintf(out_fp, "area%s", sep); fprintf(out_fp, "perimeter%s", sep); fprintf(out_fp, "compact_square%s", sep); fprintf(out_fp, "compact_circle%s", sep); fprintf(out_fp, "fd%s", sep); fprintf(out_fp, "mean_x%s", sep); fprintf(out_fp, "mean_y"); fprintf(out_fp, "\n"); /* print table body */ for (i = 0; i < n_objects; i++) { G_percent(i, n_objects - 1, 1); /* skip empty objects */ if (obj_geos[i].area == 0) continue; fprintf(out_fp, "%d%s", min + i, sep); fprintf(out_fp, "%f%s", obj_geos[i].area, sep); fprintf(out_fp, "%f%s", obj_geos[i].perimeter, sep); fprintf(out_fp, "%f%s", 4 * sqrt(obj_geos[i].area) / obj_geos[i].perimeter, sep); fprintf(out_fp, "%f%s", obj_geos[i].perimeter / (2 * sqrt(M_PI * obj_geos[i].area)), sep); /* log 1 = 0, so avoid that by always adding 0.001 to the area: */ fprintf(out_fp, "%f%s", 2 * log(obj_geos[i].perimeter) / log(obj_geos[i].area + 0.001), sep); if (!flag_m->answer) obj_geos[i].num = obj_geos[i].area; fprintf(out_fp, "%f%s", obj_geos[i].x / obj_geos[i].num, sep); fprintf(out_fp, "%f", obj_geos[i].y / obj_geos[i].num); /* object id: i + min */ /* TODO */ /* 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 */ /* variance of X and Y to approximate bounding ellipsoid */ fprintf(out_fp, "\n"); } if (out_fp != stdout) fclose(out_fp); exit(EXIT_SUCCESS); }