[GRASS-SVN] r59128 - grass/trunk/raster/r.clump
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Feb 23 14:09:37 PST 2014
Author: mmetz
Date: 2014-02-23 14:09:37 -0800 (Sun, 23 Feb 2014)
New Revision: 59128
Modified:
grass/trunk/raster/r.clump/clump.c
grass/trunk/raster/r.clump/local_proto.h
grass/trunk/raster/r.clump/main.c
Log:
r.clump: optimize, remove code duplication
Modified: grass/trunk/raster/r.clump/clump.c
===================================================================
--- grass/trunk/raster/r.clump/clump.c 2014-02-23 20:45:30 UTC (rev 59127)
+++ grass/trunk/raster/r.clump/clump.c 2014-02-23 22:09:37 UTC (rev 59128)
@@ -17,6 +17,10 @@
*
***************************************************************************/
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
#include <time.h>
#include <grass/gis.h>
#include <grass/raster.h>
@@ -25,246 +29,24 @@
#define INCR 1024
-/* 4 neighbor algorithm */
-CELL clump(int in_fd, int out_fd)
+CELL clump(int in_fd, int out_fd, int diag)
{
register int col;
register CELL *prev_clump, *cur_clump;
- register CELL *index, *index2;
+ register CELL *index, *renumber;
register int n;
CELL *prev_in, *cur_in;
CELL *temp_cell, *temp_clump, *out_cell;
- CELL X, UP, LEFT, NEW, OLD;
- CELL label;
- int nrows, ncols;
- int row;
- int len;
- int pass;
- int nalloc;
- long cur_time;
- int column;
-
- nrows = Rast_window_rows();
- ncols = Rast_window_cols();
-
- /* allocate clump index */
- nalloc = INCR;
- index = (CELL *) G_malloc(nalloc * sizeof(CELL));
- index[0] = 0;
- index2 = NULL;
-
- /* allocate CELL buffers one column larger than current window */
- len = (ncols + 1) * sizeof(CELL);
- prev_in = (CELL *) G_malloc(len);
- cur_in = (CELL *) G_malloc(len);
- prev_clump = (CELL *) G_malloc(len);
- cur_clump = (CELL *) G_malloc(len);
- out_cell = (CELL *) G_malloc(len);
-
-/******************************** PASS 1 ************************************
- * first pass thru the input simulates the clump to determine
- * the reclumping index.
- * second pass does the clumping for real
- */
- time(&cur_time);
- for (pass = 1; pass <= 2; pass++) {
- /* second pass must generate a renumbering scheme */
- if (pass == 2) {
- CELL cat;
-
- cat = 1;
- index2 = (CELL *) G_malloc((label + 1) * sizeof(CELL));
- index2[0] = 0;
- for (n = 1; n <= label; n++) {
- OLD = n;
- NEW = index[n];
- if (OLD != NEW) {
- /* find final clump id */
- while (OLD != NEW) {
- OLD = NEW;
- NEW = index[OLD];
- }
- index[n] = NEW;
- }
- else
- index2[n] = cat++;
- }
- }
-
- /* fake a previous row which is all zero */
- Rast_set_c_null_value(prev_in, ncols + 1);
- G_zero(prev_clump, len);
-
- /* create a left edge of zero */
- cur_in[0] = 0;
- cur_clump[0] = 0;
-
- /* initialize clump labels */
- label = 0;
-
- G_message(_("Pass %d..."), pass);
- for (row = 0; row < nrows; row++) {
- Rast_get_c_row(in_fd, cur_in + 1, row);
-
- G_percent(row, nrows, 4);
- X = 0;
- Rast_set_c_null_value(&X, 1);
- for (col = 1; col <= ncols; col++) {
- LEFT = X;
- X = cur_in[col];
- if (Rast_is_c_null_value(&X)) { /* don't clump NULL data */
- cur_clump[col] = 0;
- continue;
- }
-
- UP = prev_in[col];
-
- /*
- * if the cell value is different above and to the left
- * then we must start a new clump
- *
- * this new clump may eventually collide with another
- * clump and have to be merged
- */
- if (X != LEFT && X != UP) { /* start a new clump */
- label++;
- cur_clump[col] = label;
- if (pass == 1) {
- if (label >= nalloc) {
- nalloc += INCR;
- index =
- (CELL *) G_realloc(index,
- nalloc * sizeof(CELL));
- }
- index[label] = label;
- }
- continue;
- }
- if (X == LEFT && X != UP) { /* same clump as to the left */
- cur_clump[col] = cur_clump[col - 1];
- continue;
- }
- if (X == UP && X != LEFT) { /* same clump as above */
- cur_clump[col] = prev_clump[col];
- continue;
- }
-
- /*
- * at this point the cell value X is the same as LEFT and UP
- * so it should go into the same clump. It is possible for
- * the clump value on the left to differ from the clump value
- * above. If this happens we have a conflict and one of the
- * LEFT or UP needs to be reclumped
- */
- if (cur_clump[col - 1] == prev_clump[col]) { /* ok */
- cur_clump[col] = prev_clump[col];
- continue;
- }
-
- /* conflict! preserve the clump from above and change the left.
- * Must also go back to the left in the current row and to the right
- * in the previous row to change all the clump values as well.
- *
- */
-
- NEW = prev_clump[col];
- OLD = cur_clump[col - 1];
- cur_clump[col] = NEW;
-
- /* to left
- for (n = 1; n < col; n++)
- if (cur_clump[n] == OLD)
- cur_clump[n] = NEW;
- */
-
- temp_clump = cur_clump;
- n = col - 1;
- while (n-- > 0) {
- temp_clump++; /* skip left edge */
- if (*temp_clump == OLD)
- *temp_clump = NEW;
- }
-
- /* to right
- for (n = col+1; n <= ncols; n++)
- if (prev_clump[n] == OLD)
- prev_clump[n] = NEW;
- */
-
- temp_clump = prev_clump;
- temp_clump += col;
- n = ncols - col;
- while (n-- > 0) {
- temp_clump++; /* skip col */
- if (*temp_clump == OLD)
- *temp_clump = NEW;
- }
-
- /* modify the indexes
- if (pass == 1)
- for (n = 1; n <= label; n++)
- if (index[n] == OLD)
- index[n] = NEW;
- */
-
- if (pass == 1)
- index[OLD] = NEW;
- }
-
- if (pass == 2) {
- /*
- for (col = 1; col <= ncols; col++)
- out_cell[col] = index[cur_clump[col]];
-
- Rast_put_row (out_fd, out_cell+1, CELL_TYPE);
- */
- temp_clump = cur_clump;
- temp_cell = out_cell;
-
- for (column = 0; column < ncols; column++) {
- temp_clump++; /* skip left edge */
- *temp_cell = index2[index[*temp_clump]];
- if (*temp_cell == 0)
- Rast_set_null_value(temp_cell, 1, CELL_TYPE);
- temp_cell++;
- }
- Rast_put_row(out_fd, out_cell, CELL_TYPE);
- }
-
- /* switch the buffers so that the current buffer becomes the previous */
- temp_cell = cur_in;
- cur_in = prev_in;
- prev_in = temp_cell;
-
- temp_clump = cur_clump;
- cur_clump = prev_clump;
- prev_clump = temp_clump;
- }
- G_percent(1, 1, 1);
-
- print_time(&cur_time);
- }
- return 0;
-}
-
-/* 8 neighbor algorithm */
-CELL clumpd(int in_fd, int out_fd)
-{
- register int col;
- register CELL *prev_clump, *cur_clump;
- register CELL *index, *index2;
- register int n;
- CELL *prev_in, *cur_in;
- CELL *temp_cell, *temp_clump, *out_cell;
CELL X, UP, LEFT, UL, UR, NEW, OLD;
CELL label;
int nrows, ncols;
int row;
int len;
- int pass;
int nalloc;
long cur_time;
- int column;
+ char *cname;
+ int cfd, csize;
+ CELL cat;
nrows = Rast_window_rows();
ncols = Rast_window_cols();
@@ -273,7 +55,7 @@
nalloc = INCR;
index = (CELL *) G_malloc(nalloc * sizeof(CELL));
index[0] = 0;
- index2 = NULL;
+ renumber = NULL;
/* allocate CELL buffers two columns larger than current window */
len = (ncols + 2) * sizeof(CELL);
@@ -283,92 +65,84 @@
cur_clump = (CELL *) G_malloc(len);
out_cell = (CELL *) G_malloc(len);
-/******************************** PASS 1 ************************************
- * first pass thru the input simulates the clump to determine
- * the reclumping index.
- * second pass does the clumping for real
- */
+ /* temp file for initial clump IDs */
+ cname = G_tempfile();
+ if ((cfd = open(cname, O_RDWR | O_CREAT | O_EXCL, 0600)) < 0)
+ G_fatal_error(_("Unable to open temp file"));
+ csize = ncols * sizeof(CELL);
+
time(&cur_time);
- for (pass = 1; pass <= 2; pass++) {
- /* second pass must generate a renumbering scheme */
- if (pass == 2) {
- CELL cat;
- cat = 1;
- index2 = (CELL *) G_malloc((label + 1) * sizeof(CELL));
- index2[0] = 0;
- for (n = 1; n <= label; n++) {
- OLD = n;
- NEW = index[n];
- if (OLD != NEW) {
- /* find final clump id */
- while (OLD != NEW) {
- OLD = NEW;
- NEW = index[OLD];
- }
- index[n] = NEW;
- }
- else
- index2[n] = cat++;
- }
- }
+ /* fake a previous row which is all zero */
+ Rast_set_c_null_value(prev_in, ncols + 2);
+ G_zero(prev_clump, len);
- /* fake a previous row which is all zero */
- Rast_set_c_null_value(prev_in, ncols + 2);
- G_zero(prev_clump, len);
+ /* 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);
- /* 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);
+ /* set above left and right to NULL for diag == 0 */
+ Rast_set_c_null_value(&UL, 1);
+ Rast_set_c_null_value(&UR, 1);
- /* create a left and right edge of zero */
- G_zero(cur_clump, len);
+ /* create a left and right edge of zero */
+ G_zero(cur_clump, len);
- /* initialize clump labels */
- label = 0;
+ /* initialize clump labels */
+ label = 0;
- G_message(_("Pass %d..."), pass);
- for (row = 0; row < nrows; row++) {
- Rast_get_c_row(in_fd, cur_in + 1, row);
+ /****************************************************
+ * PASS 1 *
+ * pass thru the input, create initial clump labels *
+ ****************************************************/
- G_percent(row, nrows, 4);
- X = 0;
- Rast_set_c_null_value(&X, 1);
- for (col = 1; col <= ncols; col++) {
- LEFT = X;
- X = cur_in[col];
- if (Rast_is_c_null_value(&X)) { /* don't clump NULL data */
- cur_clump[col] = 0;
- continue;
- }
+ G_message(_("Pass 1..."));
+ for (row = 0; row < nrows; row++) {
+ Rast_get_c_row(in_fd, cur_in + 1, row);
- UL = prev_in[col - 1];
- UP = prev_in[col];
- UR = prev_in[col + 1];
+ G_percent(row, nrows, 4);
+ Rast_set_c_null_value(&X, 1);
+ for (col = 1; col <= ncols; col++) {
+ LEFT = X;
+ X = cur_in[col];
+ if (Rast_is_c_null_value(&X)) { /* don't clump NULL data */
+ cur_clump[col] = 0;
+ continue;
+ }
- /*
- * if the cell value is different above and to the left
- * then we must start a new clump
- *
- * this new clump may eventually collide with another
- * clump and have to be merged
- */
- if (X != LEFT && X != UP && X != UL && X != UR) { /* start a new clump */
+ /*
+ * if the cell value is different to the left and above
+ * (diagonal: and above left and above right)
+ * then we must start a new clump
+ *
+ * this new clump may eventually collide with another
+ * clump and will have to be merged
+ */
+
+ /* only one "if (diag)" for speed */
+ if (diag) {
+ temp_cell = prev_in + col - 1;
+ UL = *temp_cell++;
+ UP = *temp_cell++;
+ UR = *temp_cell;
+
+ /* start a new clump */
+ if (X != LEFT && X != UP && X != UL && X != UR) {
label++;
cur_clump[col] = label;
- if (pass == 1) {
- if (label >= nalloc) {
- nalloc += INCR;
- index =
- (CELL *) G_realloc(index,
- nalloc * sizeof(CELL));
- }
- index[label] = label;
+ if (label >= nalloc) {
+ nalloc += INCR;
+ index =
+ (CELL *) G_realloc(index,
+ nalloc * sizeof(CELL));
}
+ index[label] = label;
continue;
}
+
OLD = NEW = 0;
- if (X == LEFT) { /* same clump as to the left */
+ /* same clump as to the left */
+ if (X == LEFT) {
OLD = cur_clump[col - 1];
cur_clump[col] = OLD;
}
@@ -383,7 +157,7 @@
cur_clump[col] = OLD;
}
else {
- NEW = *temp_clump;
+ NEW = *temp_clump;
cur_clump[col] = NEW;
}
}
@@ -392,80 +166,159 @@
temp_cell--;
temp_clump--;
} while (n-- > 0);
+ }
+ else {
+ UP = prev_in[col];
- if (NEW == 0 || OLD == NEW) { /* ok */
+ /* start a new clump */
+ if (X != LEFT && X != UP) {
+ label++;
+ cur_clump[col] = label;
+ if (label >= nalloc) {
+ nalloc += INCR;
+ index =
+ (CELL *) G_realloc(index,
+ nalloc * sizeof(CELL));
+ }
+ index[label] = label;
continue;
}
- /* conflict! preserve the clump from above and change the left.
- * Must also go back to the left in the current row and to the right
- * in the previous row to change all the clump values as well.
- *
- */
-
- /* left of the current row from 1 to col - 1 */
- temp_clump = cur_clump;
- n = col - 1;
- while (n-- > 0) {
- temp_clump++; /* skip left edge */
- if (*temp_clump == OLD)
- *temp_clump = NEW;
+ OLD = NEW = 0;
+ /* same clump as to the left */
+ if (X == LEFT) {
+ OLD = cur_clump[col - 1];
+ cur_clump[col] = OLD;
}
-
- /* right of previous row from col + 1 to ncols */
- temp_clump = prev_clump;
- temp_clump += col;
- n = ncols - col;
- while (n-- > 0) {
- temp_clump++; /* skip col */
- if (*temp_clump == OLD)
- *temp_clump = NEW;
+ /* same clump as above */
+ if (X == UP) {
+ if (OLD == 0) {
+ OLD = prev_clump[col];
+ cur_clump[col] = OLD;
+ }
+ else {
+ NEW = prev_clump[col];
+ cur_clump[col] = NEW;
+ }
}
+ }
- /* modify the indexes
- if (pass == 1)
- for (n = 1; n <= label; n++)
- if (index[n] == OLD)
- index[n] = NEW;
- */
-
- if (pass == 1)
- index[OLD] = NEW;
+ if (NEW == 0 || OLD == NEW) { /* ok */
+ continue;
}
- if (pass == 2) {
- /*
- for (col = 1; col <= ncols; col++)
- out_cell[col] = index[cur_clump[col]];
+ /* conflict! preserve NEW clump ID and change OLD clump ID.
+ * Must go back to the left in the current row and to the right
+ * in the previous row to change all the clump values as well.
+ */
- Rast_put_row (out_fd, out_cell+1, CELL_TYPE);
- */
- temp_clump = cur_clump;
- temp_cell = out_cell;
+ /* left of the current row from 1 to col - 1 */
+ temp_clump = cur_clump;
+ n = col - 1;
+ while (n-- > 0) {
+ temp_clump++; /* skip left edge */
+ if (*temp_clump == OLD)
+ *temp_clump = NEW;
+ }
- for (column = 0; column < ncols; column++) {
- temp_clump++; /* skip left edge */
- *temp_cell = index2[index[*temp_clump]];
- if (*temp_cell == 0)
- Rast_set_null_value(temp_cell, 1, CELL_TYPE);
- temp_cell++;
- }
- Rast_put_row(out_fd, out_cell, CELL_TYPE);
+ /* right of previous row from col + 1 to ncols */
+ temp_clump = prev_clump;
+ temp_clump += col;
+ n = ncols - col;
+ while (n-- > 0) {
+ temp_clump++; /* skip col */
+ if (*temp_clump == OLD)
+ *temp_clump = NEW;
}
- /* switch the buffers so that the current buffer becomes the previous */
- temp_cell = cur_in;
- cur_in = prev_in;
- prev_in = temp_cell;
+ /* modify the OLD index */
+ index[OLD] = NEW;
+ }
- temp_clump = cur_clump;
- cur_clump = prev_clump;
- prev_clump = temp_clump;
+ /* write initial clump IDs */
+ /* this works also with writing out cur_clump,
+ * but only prev_clump is complete */
+ if (row > 0) {
+ if (write(cfd, prev_clump + 1, csize) != csize)
+ G_fatal_error(_("Unable to write to temp file"));
}
- G_percent(1, 1, 1);
- print_time(&cur_time);
+ /* switch the buffers so that the current buffer becomes the previous */
+ temp_cell = cur_in;
+ cur_in = prev_in;
+ prev_in = temp_cell;
+
+ temp_clump = cur_clump;
+ cur_clump = prev_clump;
+ prev_clump = temp_clump;
}
+ /* write last row with initial clump IDs */
+ if (write(cfd, prev_clump + 1, csize) != csize)
+ G_fatal_error(_("Unable to write to temp file"));
+ G_percent(1, 1, 1);
+
+ /* generate a renumbering scheme */
+ G_message(_("Generating renumbering scheme..."));
+ G_debug(0, "%d initial labels", label);
+ /* allocate final clump ID */
+ renumber = (CELL *) G_malloc((label + 1) * sizeof(CELL));
+ renumber[0] = 0;
+ cat = 1;
+ G_percent(0, label, 4);
+ for (n = 1; n <= label; n++) {
+ G_percent(n, label, 4);
+ OLD = n;
+ NEW = index[n];
+ if (OLD != NEW) {
+ renumber[n] = 0;
+ /* find valid clump ID */
+ while (OLD != NEW) {
+ OLD = NEW;
+ NEW = index[OLD];
+ }
+ index[n] = NEW;
+ }
+ else
+ /* set final clump id */
+ renumber[n] = cat++;
+ }
+
+ /* rewind temp file */
+ lseek(cfd, 0, SEEK_SET);
+
+ /****************************************************
+ * PASS 2 *
+ * apply renumbering scheme *
+ ****************************************************/
+
+ /* the input raster is no longer needed,
+ * instead we use the temp file with the initial clump labels */
+ G_message(_("Pass 2..."));
+ for (row = 0; row < nrows; row++) {
+
+ G_percent(row, nrows, 4);
+
+ if (read(cfd, cur_clump, csize) != csize)
+ G_fatal_error(_("Unable to read from temp file"));
+
+ temp_clump = cur_clump;
+ temp_cell = out_cell;
+
+ for (col = 0; col < ncols; col++) {
+ *temp_cell = renumber[index[*temp_clump]];
+ if (*temp_cell == 0)
+ Rast_set_null_value(temp_cell, 1, CELL_TYPE);
+ temp_clump++;
+ temp_cell++;
+ }
+ Rast_put_row(out_fd, out_cell, CELL_TYPE);
+ }
+ G_percent(1, 1, 1);
+ close(cfd);
+ unlink(cname);
+
+ print_time(&cur_time);
+
return 0;
}
Modified: grass/trunk/raster/r.clump/local_proto.h
===================================================================
--- grass/trunk/raster/r.clump/local_proto.h 2014-02-23 20:45:30 UTC (rev 59127)
+++ grass/trunk/raster/r.clump/local_proto.h 2014-02-23 22:09:37 UTC (rev 59128)
@@ -20,8 +20,7 @@
#define __LOCAL_PROTO_H__
/* clump.c */
-CELL clump(int, int);
-CELL clumpd(int, int);
+CELL clump(int, int, int);
int print_time(long *);
/* main.c */
Modified: grass/trunk/raster/r.clump/main.c
===================================================================
--- grass/trunk/raster/r.clump/main.c 2014-02-23 20:45:30 UTC (rev 59127)
+++ grass/trunk/raster/r.clump/main.c 2014-02-23 22:09:37 UTC (rev 59128)
@@ -80,17 +80,13 @@
out_fd = Rast_open_c_new(OUTPUT);
- if (flag_diag->answer)
- clumpd(in_fd, out_fd);
- else
- clump(in_fd, out_fd);
+ clump(in_fd, out_fd, flag_diag->answer);
G_debug(1, "Creating support files...");
Rast_close(in_fd);
Rast_close(out_fd);
-
/* build title */
if (opt_title->answer != NULL)
strcpy(title, opt_title->answer);
More information about the grass-commit
mailing list