[GRASS-SVN] r43480 - grass-addons/grass7/raster/r.clump2
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Sep 16 12:06:56 EDT 2010
Author: mmetz
Date: 2010-09-16 16:06:56 +0000 (Thu, 16 Sep 2010)
New Revision: 43480
Modified:
grass-addons/grass7/raster/r.clump2/main.c
grass-addons/grass7/raster/r.clump2/pq.c
grass-addons/grass7/raster/r.clump2/r.clump2.html
Log:
add start point option
Modified: grass-addons/grass7/raster/r.clump2/main.c
===================================================================
--- grass-addons/grass7/raster/r.clump2/main.c 2010-09-16 13:00:46 UTC (rev 43479)
+++ grass-addons/grass7/raster/r.clump2/main.c 2010-09-16 16:06:56 UTC (rev 43480)
@@ -38,7 +38,8 @@
struct Option *opt_in;
struct Option *opt_out;
struct Option *opt_title;
- struct Flag *sides_flag;
+ struct Option *opt_coord;
+ struct Flag *flag_sides;
CELL clump_no, *out_buf;
struct Cell_head window, cellhd;
int nrows, ncols, r, c;
@@ -48,7 +49,7 @@
FLAG *inlist;
int ramseg;
long index, index_nbr;
- int sides, s, r_nbr, c_nbr;
+ int sides, s, r_nbr, c_nbr, have_seeds = 0, *id_map = NULL;
int nextdr[8] = { 1, -1, 0, 0, -1, 1, 1, -1 };
int nextdc[8] = { 0, 0, -1, 1, 1, -1, 1, -1 };
@@ -73,12 +74,19 @@
opt_title->type = TYPE_STRING;
opt_title->required = NO;
opt_title->description = _("Title");
+
+ opt_coord = G_define_option();
+ opt_coord->key = "coordinate";
+ opt_coord->type = TYPE_STRING;
+ opt_coord->key_desc = "x,y";
+ opt_coord->multiple = YES;
+ opt_coord->description =
+ _("Map grid coordinates of starting points (E,N)");
- sides_flag = G_define_flag();
- sides_flag->key = 'e';
- sides_flag->description = _("Ignore diagonal cells");
+ flag_sides = G_define_flag();
+ flag_sides->key = 'e';
+ flag_sides->description = _("Ignore diagonal cells");
-
/* parse options */
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -91,7 +99,7 @@
if (out_fd < 0)
G_fatal_error(_("Unable to create raster map <%s>"), opt_out->answer);
- if (sides_flag->answer)
+ if (flag_sides->answer)
sides = 4;
else
sides = 8;
@@ -148,9 +156,6 @@
Rast_close(in_fd);
G_free(in_buf);
- /* determine clumps */
- G_message(_("determine clumps ..."));
-
ncells = nrows * ncols;
counter = 1;
@@ -158,12 +163,60 @@
init_pq(nrows * ncols * 0.02);
clump_no = 1;
- index = SEG_INDEX(ramseg, 0, 0);
- if (clump_id[index] != 0)
- clump_id[index] = clump_no++;
- add_pnt(index);
- FLAG_SET(inlist, 0, 0);
+ if (opt_coord->answers) {
+ int col, row, i;
+ double east, north;
+ char **answers = opt_coord->answers;
+ struct Cell_head window;
+
+ /* get start point coordinates */
+ G_message(_("get start point coordinates ..."));
+ G_get_window(&window);
+
+ for (; *answers != NULL; answers += 2) {
+ if (!G_scan_easting(*answers, &east, G_projection()))
+ G_fatal_error(_("Illegal x coordinate <%s>"), *answers);
+ if (*(answers + 1) == NULL)
+ G_fatal_error(_("Missing y corrdinate"));
+ if (!G_scan_northing(*(answers + 1), &north, G_projection()))
+ G_fatal_error(_("Illegal y coordinate <%s>"), *(answers + 1));
+
+ if (east < window.west || east > window.east ||
+ north < window.south || north > window.north) {
+ G_warning(_("Warning, ignoring point outside window: %.4f,%.4f"),
+ east, north);
+ continue;
+ }
+ row = (window.north - north) / window.ns_res;
+ col = (east - window.west) / window.ew_res;
+
+ index = SEG_INDEX(ramseg, row, col);
+ if (clump_id[index] != 0) {
+ /* add seed points with decreasing clump ID: first point will seed first clump */
+ clump_id[index] = clump_no++;
+ add_pnt(index);
+ FLAG_SET(inlist, row, col);
+ }
+ }
+
+ id_map = G_malloc(clump_no * sizeof(int));
+ for (i = 0; i < clump_no; i++)
+ id_map[i] = i;
+
+ have_seeds = 1;
+ }
+ else {
+ index = SEG_INDEX(ramseg, 0, 0);
+ if (clump_id[index] != 0)
+ clump_id[index] = clump_no++;
+ add_pnt(index);
+ FLAG_SET(inlist, 0, 0);
+ }
+
+ /* determine clumps */
+ G_message(_("determine clumps ..."));
+
while (pqsize > 0) {
int start_new = 1;
@@ -192,27 +245,53 @@
}
}
else if (clump_id[index] == 0 && clump_id[index_nbr] == 0) {
- add_pnt(index_nbr);
- FLAG_SET(inlist, r_nbr, c_nbr);
+ if (!have_seeds) {
+ add_pnt(index_nbr);
+ FLAG_SET(inlist, r_nbr, c_nbr);
+ }
}
else if (clump_id[index] != 0) {
if (clump_id[index_nbr] != 0)
clump_id[index_nbr] = clump_id[index];
- add_pnt(index_nbr);
- FLAG_SET(inlist, r_nbr, c_nbr);
+ if (!have_seeds ||
+ (have_seeds && clump_id[index_nbr] != 0)) {
+ add_pnt(index_nbr);
+ FLAG_SET(inlist, r_nbr, c_nbr);
+ }
}
}
+ else if (have_seeds) {
+ /* safety check */
+ index_nbr = SEG_INDEX(ramseg, r_nbr, c_nbr);
+ if (clump_id[index] != clump_id[index_nbr]) {
+ G_warning(_("%d. and %d. start point are in the same clump, ignoring %d. start point"),
+ clump_id[index], clump_id[index_nbr], clump_id[index_nbr]);
+ id_map[clump_id[index_nbr]] = clump_id[index];
+ clump_id[index_nbr] = clump_id[index];
+ }
+ }
}
}
}
free_pq();
- if (counter < ncells)
- G_warning("missed some cells!");
-
- flag_destroy(inlist);
-
+ if (counter < ncells) {
+ if (!have_seeds)
+ G_warning("missed some cells!");
+ else
+ G_percent(ncells, ncells, 1);
+ }
+
+ if (have_seeds && clump_no > 1) {
+ int i;
+
+ for (i = 1; i < clump_no; i++) {
+ if (id_map[i] - id_map[i - 1] > 1)
+ id_map[i] = id_map[i - 1] + 1;
+ }
+ }
+
/* write output */
G_message(_("write output map ..."));
out_buf = Rast_allocate_buf(CELL_TYPE);
@@ -220,6 +299,13 @@
G_percent(r, nrows, 2);
for (c = 0; c < ncols; c++) {
index = SEG_INDEX(ramseg, r, c);
+ if (have_seeds) {
+ if ((FLAG_GET(inlist, r, c)) == 0) {
+ clump_id[index] = 0;
+ }
+ else
+ clump_id[index] = id_map[clump_id[index]];
+ }
if (clump_id[index] == 0)
Rast_set_c_null_value(&out_buf[c], 1);
else
@@ -232,7 +318,11 @@
G_free(out_buf);
Rast_close(out_fd);
G_free(clump_id);
+ flag_destroy(inlist);
+ if (id_map)
+ G_free(id_map);
+
G_debug(1, "Creating support files...");
/* build title */
Modified: grass-addons/grass7/raster/r.clump2/pq.c
===================================================================
--- grass-addons/grass7/raster/r.clump2/pq.c 2010-09-16 13:00:46 UTC (rev 43479)
+++ grass-addons/grass7/raster/r.clump2/pq.c 2010-09-16 16:06:56 UTC (rev 43480)
@@ -27,6 +27,15 @@
static int heap_step;
static long *heap_index;
+int cmp_clump(long cid1, long cid2)
+{
+ if (!cid1)
+ return 0;
+ if (!cid2)
+ return 1;
+ return (cid1 < cid2);
+}
+
int init_pq(int step)
{
pqsize = 0;
@@ -57,7 +66,7 @@
parent = GET_PARENT(child);
/* child is larger */
- if (clump_id[child_pnt] > clump_id[heap_index[parent]]) {
+ if (cmp_clump(clump_id[child_pnt], clump_id[heap_index[parent]])) {
/* push parent point down */
heap_index[child] = heap_index[parent];
child = parent;
@@ -120,8 +129,8 @@
i = child + 3;
/* get largest child */
while (childr < i && childr <= pqsize) {
- if (clump_id[heap_index[childr]] >
- clump_id[heap_index[child]]) {
+ if (cmp_clump(clump_id[heap_index[childr]],
+ clump_id[heap_index[child]])) {
child = childr;
}
childr++;
Modified: grass-addons/grass7/raster/r.clump2/r.clump2.html
===================================================================
--- grass-addons/grass7/raster/r.clump2/r.clump2.html 2010-09-16 13:00:46 UTC (rev 43479)
+++ grass-addons/grass7/raster/r.clump2/r.clump2.html 2010-09-16 16:06:56 UTC (rev 43480)
@@ -3,7 +3,8 @@
<em>r.clump2</em> finds all areas of contiguous cell category values in the
input raster map layer <em>name.</em> It assigns a unique category value
to each such area ("clump") in the resulting output raster map layer
-<em>name.</em>
+<em>name.</em> Optionally <em>r.clump2</em> uses start points and finds
+only those clumps where the start points are falling into.
<h4>Differences to r.clump</h4>
More information about the grass-commit
mailing list