r.clump: optimize, remove code duplication

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

raster/r.clump/clump.c

@@ -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,234 +29,11 @@
 
 #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 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 CELL *index, *renumber;
     register int n;
     CELL *prev_in, *cur_in;
     CELL *temp_cell, *temp_clump, *out_cell;
@@ -261,10 +42,11 @@ CELL clumpd(int in_fd, int out_fd)
     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 @@ CELL clumpd(int in_fd, int out_fd)
     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 @@ CELL clumpd(int in_fd, int out_fd)
     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 + 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);
+    /* 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);
 
-	/* create a left and right edge of zero */
-	G_zero(cur_clump, len);
-
-	/* 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);
+    time(&cur_time);
 
-	    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;
-		}
+    /* 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 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);
+
+    /* initialize clump labels */
+    label = 0;
+
+    /****************************************************
+     *                      PASS 1                      *
+     * pass thru the input, create initial clump labels *
+     ****************************************************/
+
+    G_message(_("Pass 1..."));
+    for (row = 0; row < nrows; row++) {
+	Rast_get_c_row(in_fd, cur_in + 1, row);
+
+	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;
+	    }
 
-		UL = prev_in[col - 1];
-		UP = prev_in[col];
-		UR = prev_in[col + 1];
-
-		/*
-		 * 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 @@ CELL clumpd(int in_fd, int out_fd)
 			    cur_clump[col] = OLD;
 			}
 			else {
-			    NEW  = *temp_clump;
+			    NEW = *temp_clump;
 			    cur_clump[col] = NEW;
 			}
 		    }
@@ -392,80 +166,159 @@ CELL clumpd(int in_fd, int out_fd)
 		    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 (NEW == 0 || OLD == NEW) {	/* ok */
+		continue;
+	    }
+
+	    /* 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.
+	     */
 
-		if (pass == 1)
-		    index[OLD] = NEW;
+	    /* 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;
 	    }
 
-	    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);
+	    /* 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;
 }
 

raster/r.clump/local_proto.h

@@ -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 */

raster/r.clump/main.c

@@ -80,17 +80,13 @@ int main(int argc, char *argv[])
 
     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);