i.segment: improvement for https://trac.osgeo.org/grass/ticket/3084

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

imagery/i.segment/create_isegs.c

@@ -84,6 +84,48 @@ int create_isegs(struct globals *globals)
 	}
     }
 
+    if (globals->method == ORM_RG) {
+	/* renumber */
+	int i, *new_id, max_id;
+
+	G_debug(1, "Largest assigned ID: %d", globals->max_rid);
+
+	new_id = G_malloc((globals->max_rid + 1) * sizeof(int));
+	
+	for (i = 0; i <= globals->max_rid; i++)
+	    new_id[i] = 0;
+
+	for (row = 0; row < globals->nrows; row++) {
+	    for (col = 0; col < globals->ncols; col++) {
+		Segment_get(&globals->rid_seg, &rid, row, col);
+		if (!Rast_is_c_null_value(&rid))
+		    new_id[rid]++;
+	    }
+	}
+
+	max_id = 0;
+	for (i = 0; i <= globals->max_rid; i++) {
+	    if (new_id[i] > 0) {
+		max_id++;
+		new_id[i] = max_id;
+	    }
+	}
+	globals->max_rid = max_id;
+	G_debug(1, "Largest renumbered ID: %d", globals->max_rid);
+	
+	for (row = 0; row < globals->nrows; row++) {
+	    for (col = 0; col < globals->ncols; col++) {
+		Segment_get(&globals->rid_seg, &rid, row, col);
+		if (!Rast_is_c_null_value(&rid)) {
+		    rid = new_id[rid];
+		    Segment_put(&globals->rid_seg, &rid, row, col);
+		}
+	    }
+	}
+
+	G_free(new_id);
+    }
+
     return successflag;
 }
 

imagery/i.segment/iseg.h

@@ -13,12 +13,21 @@
  *****************************************************************************/
 
 #include <grass/segment.h>
+#include <grass/imagery.h>
 #include "flag.h"
 #include "regtree.h"
 #include "ngbrtree.h"
 
 /* #def _OR_SHAPE_ */
 
+#ifdef HAVE_LONG_LONG_INT
+#define LARGEINT long long
+#elif defined HAVE_LARGEFILES
+#define LARGEINT off_t
+#else
+#define LARGEINT long
+#endif
+
 /* methods */
 #define ORM_RG 1	/* region growing */
 #define ORM_MS 2	/* mean shift */
@@ -41,59 +50,76 @@ struct rclist
 /* input and output files, as well as some other processing info */
 struct globals
 {
-    /* user parameters */
+    /* input */
     char *image_group;
+    struct Ref Ref;		/* group reference list */
+    DCELL *min, *max;
     int weighted;		/* 0 if false/not selected, so we should scale input.
 				 * 1 if the scaling should be skipped */
+    int nbands;			/* number of rasters in the image group */
+    size_t datasize;		/* nbands * sizeof(double) */
+    int mb;			/* amount of memory to use in MB */
+
+    char *seeds, *bounds_map;	/* optional segment seeds and polygon constraints/boundaries */
+    CELL lower_bound, upper_bound;
+    const char *bounds_mapset;
+
+    /* output */
+    /* region growing */
+    char *out_name;		/* name of output raster map with regions */
+    char *out_band;		/* indicator for segment heterogeneity / goodness of fit */
+    /* mean_shift */
+    char *ms_suf;		/* suffix to be appended to input bands */
+
+    /* general segmentation */
     int method;			/* Segmentation method code */
     int (*method_fn)();		/* Segmentation method function */
     int nn;			/* number of neighbors, 4 or 8 */
     double max_diff;		/* max possible difference */
     double alpha;		/* similarity threshold */
+    int end_t;			/* maximum number of iterations */
+
+    /* region growing */
     int min_segment_size;	/* smallest number of pixels/cells allowed in a final segment */
 
+    /* inactive options for region growing */
     double radio_weight;	/* weighing factor radiometric - shape */
     double smooth_weight;       /* weighing factor smoothness - compactness */
 
-    int end_t;			/* maximum number of iterations */
-    int mb;
+    /* mean shift */
+    double hs, hr;		/* spectral and range bandwidth */
 
     /* region info */
     int nrows, ncols;
     int row_min, row_max, col_min, col_max; /* region constraints */
-    long ncells, notnullcells;
-
-    char *out_name;		/* name of output raster map */
-    char *seeds, *bounds_map;	/* optional segment seeds and polygon constraints/boundaries */
-    CELL lower_bound, upper_bound;
-    const char *bounds_mapset;
-    char *out_band;		/* indicator for segment heterogeneity */
+    LARGEINT ncells, notnullcells;
 
     /* file processing */
-    int nbands;			/* number of rasters in the image group */
     SEGMENT bands_seg, 	        /* input group with one or more bands */
+            bands_seg2,		/* copy of bands_seg for mean shift */
             bounds_seg,
 	    rid_seg;
+    SEGMENT *bands_in, *bands_out; /* pointers to input/output bands_seg, for mean shift */
     DCELL *bands_min, *bands_max;
     DCELL *bands_val;		/* array to hold all input values for one cell */
     DCELL *second_val;		/* array to hold all input values for another cell */
 
-    /* results */
+    /* maximum used region ID */
+    CELL max_rid;
+
+    /* region growing internal structure */
     struct RG_TREE *reg_tree;   /* search tree with region stats */
     int min_reg_size;		/* minimum region size */
     struct reg_stats rs, rs_i, rs_k;
     struct ngbr_stats ns;
-    size_t datasize;		/* nbands * sizeof(double) */
-    int n_regions;
 
     /* processing flags */
     FLAG *candidate_flag, *null_flag;	/*TODO, need some way to remember MASK/NULL values.  Was using -1, 0, 1 in int array.  Better to use 2 FLAG structures, better readibility? */
 
     /* number of remaining cells to check */
-    long candidate_count;
+    LARGEINT candidate_count;
 
     /* functions */
-	
     void (*find_neighbors) (int, int, int[8][2]);	/*parameters: row, col, neighbors */
     double (*calculate_similarity) (struct ngbr_stats *,
                                     struct ngbr_stats *,
@@ -142,7 +168,8 @@ void rclist_add(struct rclist *, int, int);
 int rclist_drop(struct rclist *, struct rc *);
 void rclist_destroy(struct rclist *);
 
-
 /* write_output.c */
-int write_output(struct globals *);
+int write_ids(struct globals *);
+int write_gof_rg(struct globals *);
+int write_bands_ms(struct globals *);
 int close_files(struct globals *);

imagery/i.segment/main.c

@@ -20,6 +20,7 @@
  *****************************************************************************/
 
 #include <stdlib.h>
+#include <stdio.h>
 #include <grass/gis.h>
 #include <grass/glocale.h>
 #include "iseg.h"
@@ -40,7 +41,7 @@ int main(int argc, char *argv[])
 	_("Identifies segments (objects) from imagery data.");
 
     parse_args(argc, argv, &globals);
-	
+
     G_debug(1, "Main: starting open_files()");
     if (open_files(&globals) != TRUE)
 	G_fatal_error(_("Error in reading data"));
@@ -50,13 +51,23 @@ int main(int argc, char *argv[])
 	G_fatal_error(_("Error in creating segments"));
 
     G_debug(1, "Main: starting write_output()");
-    if (write_output(&globals) != TRUE)
-	G_fatal_error(_("Error in writing data"));
+    if (write_ids(&globals) != TRUE)
+	G_fatal_error(_("Error in writing IDs"));
+
+    if (globals.method == ORM_RG && globals.out_band) {
+	if (write_gof_rg(&globals) != TRUE)
+	    G_fatal_error(_("Error in writing goodness of fit"));
+    }
+
+    if (globals.method == ORM_MS) {
+	if (write_bands_ms(&globals) != TRUE)
+	    G_fatal_error(_("Error in writing new band values"));
+    }
 
     G_debug(1, "Main: starting close_files()");
     close_files(&globals);
 
-    G_done_msg(_("Number of segments created: %d"), globals.n_regions);
+    G_done_msg(_("Number of segments created: %d"), globals.max_rid);
 
     exit(EXIT_SUCCESS);
 }

imagery/i.segment/ngbrtree.c

@@ -44,7 +44,13 @@ static int is_red(struct NB_NODE *);
 
 int cmp_ngbr(struct ngbr_stats *a, struct ngbr_stats *b)
 {
-    return (a->id - b->id);
+    if (a->id > 0 || b->id > 0)
+	return (a->id - b->id);
+
+    /* a->id == b->id == 0 */
+    if (a->row == b->row)
+	return (a->col - b->col);
+    return (a->row - b->row);
 }
 
 
@@ -417,7 +423,7 @@ static struct ngbr_stats *nbtree_next(struct NB_TRAV *trav)
 	return NULL;		/* finished traversing */
 }
 
-/* destroy the tree */
+/* clear the tree: remove all nodes */
 void nbtree_clear(struct NB_TREE *tree)
 {
     struct NB_NODE *it;
@@ -451,6 +457,15 @@ void nbtree_clear(struct NB_TREE *tree)
     return;
 }
 
+/* destroy the tree */
+void nbtree_destroy(struct NB_TREE *tree)
+{
+    nbtree_clear(tree);
+    free(tree);
+    
+    return;
+}
+
 /* used for debugging: check for errors in tree structure */
 int nbtree_debug(struct NB_TREE *tree, struct NB_NODE *root)
 {

imagery/i.segment/ngbrtree.h

@@ -110,9 +110,11 @@ struct NB_TRAV
 /* tree functions */
 struct NB_TREE *nbtree_create(int, size_t);
 void nbtree_clear(struct NB_TREE *);
+void nbtree_destroy(struct NB_TREE *);
 int nbtree_insert(struct NB_TREE *, struct ngbr_stats *);
 int nbtree_remove(struct NB_TREE *, struct ngbr_stats *);
 struct ngbr_stats *nbtree_find(struct NB_TREE *, struct ngbr_stats *);
+int cmp_ngbr(struct ngbr_stats *, struct ngbr_stats *);
 
 /* tree traversal functions */
 int nbtree_init_trav(struct NB_TRAV *, struct NB_TREE *);

imagery/i.segment/open_files.c

@@ -13,13 +13,12 @@ static int manage_memory(int, int, struct globals *);
 
 int open_files(struct globals *globals)
 {
-    struct Ref Ref;		/* group reference list */
     int *in_fd, bounds_fd, is_null;
     int n, row, col, srows, scols, inlen, outlen, nseg;
     DCELL **inbuf;		/* buffers to store lines from each of the imagery group rasters */
     CELL *boundsbuf, bounds_val;
     int have_bounds = 0;
-    CELL s, id;
+    CELL id;
     struct Range range;	/* min/max values of bounds map */
     struct FPRange *fp_range;	/* min/max values of each input raster */
     DCELL *min, *max;
@@ -36,38 +35,41 @@ int open_files(struct globals *globals)
 
     /* ****** open the input rasters ******* */
 
-    if (!I_get_group_ref(globals->image_group, &Ref))
+    if (!I_get_group_ref(globals->image_group, &globals->Ref))
 	G_fatal_error(_("Group <%s> not found in the current mapset"),
 		      globals->image_group);
 
-    if (Ref.nfiles <= 0)
+    if (globals->Ref.nfiles <= 0)
 	G_fatal_error(_("Group <%s> contains no raster maps"),
 		      globals->image_group);
 
     /* Read Imagery Group */
 
-    in_fd = G_malloc(Ref.nfiles * sizeof(int));
-    inbuf = (DCELL **) G_malloc(Ref.nfiles * sizeof(DCELL *));
-    fp_range = G_malloc(Ref.nfiles * sizeof(struct FPRange));
-    min = G_malloc(Ref.nfiles * sizeof(DCELL));
-    max = G_malloc(Ref.nfiles * sizeof(DCELL));
+    in_fd = G_malloc(globals->Ref.nfiles * sizeof(int));
+    inbuf = (DCELL **) G_malloc(globals->Ref.nfiles * sizeof(DCELL *));
+    fp_range = G_malloc(globals->Ref.nfiles * sizeof(struct FPRange));
+    min = G_malloc(globals->Ref.nfiles * sizeof(DCELL));
+    max = G_malloc(globals->Ref.nfiles * sizeof(DCELL));
+    
+    globals->min = min;
+    globals->max = max;
 
     G_debug(1, "Opening input rasters...");
-    for (n = 0; n < Ref.nfiles; n++) {
+    for (n = 0; n < globals->Ref.nfiles; n++) {
 	inbuf[n] = Rast_allocate_d_buf();
-	in_fd[n] = Rast_open_old(Ref.file[n].name, Ref.file[n].mapset);
+	in_fd[n] = Rast_open_old(globals->Ref.file[n].name, globals->Ref.file[n].mapset);
     }
 
     /* Get min/max values of each input raster for scaling */
 
     globals->max_diff = 0.;
-    globals->nbands = Ref.nfiles;
+    globals->nbands = globals->Ref.nfiles;
 
-    for (n = 0; n < Ref.nfiles; n++) {
+    for (n = 0; n < globals->Ref.nfiles; n++) {
 	/* returns -1 on error, 2 on empty range, quitting either way. */
-	if (Rast_read_fp_range(Ref.file[n].name, Ref.file[n].mapset, &fp_range[n]) != 1)
+	if (Rast_read_fp_range(globals->Ref.file[n].name, globals->Ref.file[n].mapset, &fp_range[n]) != 1)
 	    G_fatal_error(_("No min/max found in raster map <%s>"),
-			  Ref.file[n].name);
+			  globals->Ref.file[n].name);
 	Rast_get_fp_range_min_max(&(fp_range[n]), &min[n], &max[n]);
 
 	G_debug(1, "Range for layer %d: min = %f, max = %f",
@@ -75,7 +77,7 @@ int open_files(struct globals *globals)
 	
     }
     if (globals->weighted == FALSE)
-	globals->max_diff = Ref.nfiles;
+	globals->max_diff = globals->Ref.nfiles;
     else {
 	/* max difference with selected similarity method */
 	Ri.mean = max;
@@ -89,21 +91,21 @@ int open_files(struct globals *globals)
 
     /* size of each element to be stored */
 
-    inlen = sizeof(DCELL) * Ref.nfiles;
+    inlen = sizeof(DCELL) * globals->Ref.nfiles;
     outlen = sizeof(CELL);
     G_debug(1, "data element size, in: %d , out: %d ", inlen, outlen);
     globals->datasize = sizeof(double) * globals->nbands;
 
     /* count non-null cells */
-    globals->notnullcells = (long)globals->nrows * globals->ncols;
+    globals->notnullcells = (LARGEINT)globals->nrows * globals->ncols;
     for (row = 0; row < globals->nrows; row++) {
-	for (n = 0; n < Ref.nfiles; n++) {
+	for (n = 0; n < globals->Ref.nfiles; n++) {
 	    Rast_get_d_row(in_fd[n], inbuf[n], row);
 	}
 	for (col = 0; col < globals->ncols; col++) {
 
 	    is_null = 0;	/*Assume there is data */
-	    for (n = 0; n < Ref.nfiles; n++) {
+	    for (n = 0; n < globals->Ref.nfiles; n++) {
 		if (Rast_is_d_null_value(&inbuf[n][col])) {
 		    is_null = 1;
 		}
@@ -114,7 +116,6 @@ int open_files(struct globals *globals)
 	    }
 	}
     }
-    G_verbose_message(_("Non-NULL cells: %ld"), globals->notnullcells);
     if (globals->notnullcells < 2)
 	G_fatal_error(_("Insufficient number of non-NULL cells in current region"));
 
@@ -130,21 +131,31 @@ int open_files(struct globals *globals)
 	 scols, inlen, nseg) != 1)
 	G_fatal_error("Unable to create input temporary files");
 
+    if (globals->method == ORM_MS) {
+	if (Segment_open
+	    (&globals->bands_seg2, G_tempfile(), globals->nrows, globals->ncols, srows,
+	     scols, inlen, nseg) != 1)
+	    G_fatal_error("Unable to create input temporary files");
+	
+	globals->bands_in = &globals->bands_seg;
+	globals->bands_out = &globals->bands_seg2;
+    }
+
     if (Segment_open
 	(&globals->rid_seg, G_tempfile(), globals->nrows, globals->ncols, srows,
 	 scols, outlen, nseg * 2) != 1)
 	G_fatal_error("Unable to create input temporary files");
 
     /* load input bands to segment structure */
-    if (Ref.nfiles > 1)
+    if (globals->Ref.nfiles > 1)
 	G_message(_("Loading input bands..."));
     else
 	G_message(_("Loading input band..."));
 
     globals->bands_val = (double *)G_malloc(inlen);
     globals->second_val = (double *)G_malloc(inlen);
-    /* initial segment ID */
-    s = 1;
+
+    globals->max_rid = 0;
 
     globals->row_min = globals->nrows;
     globals->row_max = 0;
@@ -152,13 +163,13 @@ int open_files(struct globals *globals)
     globals->col_max = 0;
     for (row = 0; row < globals->nrows; row++) {
 	G_percent(row, globals->nrows, 4);
-	for (n = 0; n < Ref.nfiles; n++) {
+	for (n = 0; n < globals->Ref.nfiles; n++) {
 	    Rast_get_d_row(in_fd[n], inbuf[n], row);
 	}
 	for (col = 0; col < globals->ncols; col++) {
 
 	    is_null = 0;	/*Assume there is data */
-	    for (n = 0; n < Ref.nfiles; n++) {
+	    for (n = 0; n < globals->Ref.nfiles; n++) {
 		globals->bands_val[n] = inbuf[n][col];
 		if (Rast_is_d_null_value(&inbuf[n][col])) {
 		    is_null = 1;
@@ -173,13 +184,14 @@ int open_files(struct globals *globals)
 	                    (void *)globals->bands_val, row, col) != 1)
 		G_fatal_error(_("Unable to write to temporary file"));
 
-	    if (!is_null) {
-		if (!globals->seeds) {
-		    /* sequentially number all cells with a unique segment ID */
-		    id = s;
-		    s++;
-		}
+	    if (globals->method == ORM_MS) {
+		if (Segment_put(&globals->bands_seg2,
+				(void *)globals->bands_val, row, col) != 1)
+		    G_fatal_error(_("Unable to write to temporary file"));
+	    }
 
+	    id = 0;
+	    if (!is_null) {
 		/* get min/max row/col to narrow the processing window */
 		if (globals->row_min > row)
 		    globals->row_min = row;
@@ -195,11 +207,9 @@ int open_files(struct globals *globals)
 		Rast_set_c_null_value(&id, 1);
 		FLAG_SET(globals->null_flag, row, col);
 	    }
-	    if (!globals->seeds || is_null) {
-		if (Segment_put(&globals->rid_seg,
-		                (void *)&id, row, col) != 1)
-		    G_fatal_error(_("Unable to write to temporary file"));
-	    }
+	    if (Segment_put(&globals->rid_seg,
+			    (void *)&id, row, col) != 1)
+		G_fatal_error(_("Unable to write to temporary file"));
 	}
     }
     G_percent(1, 1, 1);
@@ -210,7 +220,7 @@ int open_files(struct globals *globals)
     
     globals->row_max++;
     globals->col_max++;
-    globals->ncells = (long)(globals->row_max - globals->row_min) *
+    globals->ncells = (LARGEINT)(globals->row_max - globals->row_min) *
 			    (globals->col_max - globals->col_min);
 
     /* bounds/constraints */
@@ -281,7 +291,7 @@ int open_files(struct globals *globals)
 
     /* Free memory */
 
-    for (n = 0; n < Ref.nfiles; n++) {
+    for (n = 0; n < globals->Ref.nfiles; n++) {
 	G_free(inbuf[n]);
 	Rast_close(in_fd[n]);
     }
@@ -290,19 +300,15 @@ int open_files(struct globals *globals)
     globals->rs.mean = G_malloc(globals->datasize);
 
     globals->reg_tree = rgtree_create(globals->nbands, globals->datasize);
-    globals->n_regions = s - 1;
 
-    if (globals->seeds) {
+    if (globals->method == ORM_RG && globals->seeds) {
 	load_seeds(globals, srows, scols, nseg);
+	G_debug(1, "Number of initial regions: %d", globals->max_rid);
     }
 
-    G_debug(1, "Number of initial regions: %d", globals->n_regions);
-
     G_free(inbuf);
     G_free(in_fd);
     G_free(fp_range);
-    G_free(min);
-    G_free(max);
 
     return TRUE;
 }
@@ -313,12 +319,17 @@ static int load_seeds(struct globals *globals, int srows, int scols, int nseg)
     int row, col;
     SEGMENT seeds_seg;
     CELL *seeds_buf, seeds_val;
-    int seeds_fd;
-    int spos, sneg, have_seeds;
+    int seeds_fd, have_seeds;
+    CELL new_id, cellmax, noid;
     struct rc Ri;
 
     G_debug(1, "load_seeds()");
-    
+
+    cellmax = (1 << (sizeof(CELL) * 8 - 2)) - 1;
+    cellmax += (1 << (sizeof(CELL) * 8 - 2));
+
+    noid = 0;
+
     G_message(_("Loading seeds from raster map <%s>..."), globals->seeds);
 
     if (Segment_open
@@ -356,8 +367,7 @@ static int load_seeds(struct globals *globals, int srows, int scols, int nseg)
 	return 0;
     }
 
-    spos = 1;
-    sneg = -1;
+    new_id = 0;
 
     /* convert seeds to regions */
     G_debug(1, "convert seeds to regions");
@@ -368,17 +378,18 @@ static int load_seeds(struct globals *globals, int srows, int scols, int nseg)
 	    if (!(FLAG_GET(globals->null_flag, row, col)) && 
 	        !(FLAG_GET(globals->candidate_flag, row, col))) {
 
-		if (Rast_is_c_null_value(&(seeds_buf[col]))) {
-		    if (Segment_put(&globals->rid_seg, &sneg, row, col) != 1)
-			G_fatal_error(_("Unable to write to temporary file"));
-		    sneg--;
-		    globals->n_regions--;
-		}
-		else {
+		if (!Rast_is_c_null_value(&(seeds_buf[col]))) {
+		    if (new_id == cellmax)
+			G_fatal_error(_("Too many seeds: integer overflow"));
+
+		    new_id++;
+
 		    Ri.row = row;
 		    Ri.col = col;
-		    read_seed(globals, &seeds_seg, &Ri, spos);
-		    spos++;
+		    if (!read_seed(globals, &seeds_seg, &Ri, new_id)) {
+			new_id--;
+			Segment_put(&globals->rid_seg, (void *)&noid, Ri.row, Ri.col);
+		    }
 		}
 	    }
 	}
@@ -388,7 +399,7 @@ static int load_seeds(struct globals *globals, int srows, int scols, int nseg)
     Rast_close(seeds_fd);
     Segment_close(&seeds_seg);
 
-    globals->n_regions = spos - 1;
+    globals->max_rid = new_id;
     
     flag_clear_all(globals->candidate_flag);
     
@@ -499,9 +510,10 @@ static int read_seed(struct globals *globals, SEGMENT *seeds_seg, struct rc *Ri,
     else {
 	if (globals->rs.count > 1)
 	    update_band_vals(Ri->row, Ri->col, &(globals->rs), globals);
+	else if (globals->rs.count == 1) {
+	    return 0;
+	}
     }
-    if (globals->rs.count > 1)
-	globals->n_regions -= (globals->rs.count - 1);
 
     return 1;
 }
@@ -511,46 +523,64 @@ static int manage_memory(int srows, int scols, struct globals *globals)
 {
     double reg_size_mb, segs_mb;
     int reg_size_count, nseg, nseg_total;
+    
+    segs_mb = globals->mb;
+    if (globals->method == ORM_RG) {
+
+	/* minimum region size to store in search tree */
+	reg_size_mb = 2 * globals->datasize +     /* mean, sum */
+		      2 * sizeof(int) +           /* id, count */
+		      sizeof(unsigned char) + 
+		      2 * sizeof(struct REG_NODE *);
+	reg_size_mb /= (1024. * 1024.);
+
+	/* put aside some memory for segment structures */
+	segs_mb = globals->mb * 0.1;
+	if (segs_mb > 10)
+	    segs_mb = 10;
+
+	/* calculate number of region stats that can be kept in memory */
+	reg_size_count = (globals->mb - segs_mb) / reg_size_mb;
+	globals->min_reg_size = 3;
+	if (reg_size_count < (double) globals->notnullcells / globals->min_reg_size) {
+	    globals->min_reg_size = (double) globals->notnullcells / reg_size_count;
+	}
+	else {
+	    reg_size_count = (double) globals->notnullcells / globals->min_reg_size;
+	    /* recalculate segs_mb */
+	    segs_mb = globals->mb - reg_size_count * reg_size_mb;
+	}
 
-    /* minimum region size to store in search tree */
-    reg_size_mb = 2 * globals->datasize +     /* mean, sum */
-                  2 * sizeof(int) +           /* id, count */
-		  sizeof(unsigned char) + 
-		  2 * sizeof(struct REG_NODE *);
-    reg_size_mb /= (1024. * 1024.);
-
-    /* put aside some memory for segment structures */
-    segs_mb = globals->mb * 0.1;
-    if (segs_mb > 10)
-	segs_mb = 10;
-
-    /* calculate number of region stats that can be kept in memory */
-    reg_size_count = (globals->mb - segs_mb) / reg_size_mb;
-    globals->min_reg_size = 3;
-    if (reg_size_count < (double) globals->notnullcells / globals->min_reg_size) {
-	globals->min_reg_size = (double) globals->notnullcells / reg_size_count;
+	G_verbose_message(_("Regions with at least %d cells are stored in memory"),
+			  globals->min_reg_size);
     }
-    else {
-	reg_size_count = (double) globals->notnullcells / globals->min_reg_size;
-	/* recalculate segs_mb */
-	segs_mb = globals->mb - reg_size_count * reg_size_mb;
-    }
-
-    G_verbose_message(_("Regions with at least %d cells are stored in memory"),
-                      globals->min_reg_size);
 
     /* calculate number of segments in memory */
     if (globals->bounds_map != NULL) {
 	/* input bands, segment ids, bounds map */
-	nseg = (1024. * 1024. * segs_mb) /
-	       (sizeof(DCELL) * globals->nbands * srows * scols + 
-		sizeof(CELL) * 4 * srows * scols);
+	if (globals->method == ORM_MS) {
+	    nseg = (1024. * 1024. * segs_mb) /
+		   (sizeof(DCELL) * 2 * globals->nbands * srows * scols + 
+		    sizeof(CELL) * 4 * srows * scols);
+	}
+	else {
+	    nseg = (1024. * 1024. * segs_mb) /
+		   (sizeof(DCELL) * globals->nbands * srows * scols + 
+		    sizeof(CELL) * 4 * srows * scols);
+	}
     }
     else {
 	/* input bands, segment ids */
-	nseg = (1024. * 1024. * segs_mb) /
-	       (sizeof(DCELL) * globals->nbands * srows * scols + 
-		sizeof(CELL) * 2 * srows * scols);
+	if (globals->method == ORM_MS) {
+	    nseg = (1024. * 1024. * segs_mb) /
+		   (sizeof(DCELL) * 2 * globals->nbands * srows * scols + 
+		    sizeof(CELL) * 2 * srows * scols);
+	}
+	else {
+	    nseg = (1024. * 1024. * segs_mb) /
+		   (sizeof(DCELL) * globals->nbands * srows * scols + 
+		    sizeof(CELL) * 2 * srows * scols);
+	}
     }
     nseg_total = (globals->nrows / srows + (globals->nrows % srows > 0)) *
                  (globals->ncols / scols + (globals->ncols % scols > 0));

imagery/i.segment/parse_args.c

@@ -97,7 +97,7 @@ int parse_args(int argc, char *argv[], struct globals *globals)
     endt->key = "iterations";
     endt->type = TYPE_INTEGER;
     endt->required = NO;
-    endt->answer = "20";
+    endt->answer = "50";
     endt->description = _("Maximum number of iterations");
     endt->guisection = _("Settings");
 
@@ -146,9 +146,7 @@ int parse_args(int argc, char *argv[], struct globals *globals)
     else
 	G_fatal_error("Invalid output raster name");
 
-    /* Note: this threshold is scaled after we know more at the beginning of create_isegs() */
     globals->alpha = atof(threshold->answer);
-
     if (globals->alpha <= 0 || globals->alpha >= 1)
 	G_fatal_error(_("Threshold should be > 0 and < 1"));
 
@@ -241,6 +239,23 @@ int parse_args(int argc, char *argv[], struct globals *globals)
     globals->nrows = Rast_window_rows();
     globals->ncols = Rast_window_cols();
 
+    if (sizeof(LARGEINT) < 8) {
+	int i;
+
+	LARGEINT intmax;
+
+	intmax = ((LARGEINT)1 << (sizeof(LARGEINT) * 8 - 2)) - 1;
+	intmax += ((LARGEINT)1 << (sizeof(LARGEINT) * 8 - 2));
+
+	globals->ncells = globals->ncols;
+	for (i = 1; i < globals->nrows; i++) {
+	    if (globals->ncols > intmax - globals->ncells)
+		G_fatal_error(_("Integer overflow: too many cells in current region"));
+
+	    globals->ncells += globals->ncols;
+	}
+    }
+
     /* debug help */
     if (outband->answer == NULL)
 	globals->out_band = NULL;
@@ -255,12 +270,13 @@ int parse_args(int argc, char *argv[], struct globals *globals)
 	if (atoi(endt->answer) > 0)
 	    globals->end_t = atoi(endt->answer);
 	else {
-	    globals->end_t = 100;
-	    G_warning(_("Invalid number of iterations, 100 will be used"));
+	    globals->end_t = 50;
+	    G_warning(_("Invalid number of iterations, %d will be used"),
+	              globals->end_t);
 	}
     }
     else
-	globals->end_t = 1000;
+	globals->end_t = 50;
 
     if (mem->answer && atoi(mem->answer) > 10)
 	globals->mb = atoi(mem->answer);

imagery/i.segment/region_growing.c

@@ -25,6 +25,15 @@
 #endif
 #define MIN(a,b) ( ((a) < (b)) ? (a) : (b) )
 
+struct idlist
+{
+    int *ids;
+    int nids, nalloc;
+    CELL cellmax;
+};
+
+static struct idlist idlist;
+
 /* internal functions */
 static int merge_regions(struct ngbr_stats *, struct reg_stats *, /* Ri */
                          struct ngbr_stats *, struct reg_stats *, /* Rk */
@@ -45,6 +54,61 @@ static int find_best_neighbor(struct ngbr_stats *, struct reg_stats *,
 static int calculate_reg_stats(int, int, struct reg_stats *, 
                          struct globals *);
 
+void init_free_ids(void)
+{
+    idlist.nalloc = 10;
+    idlist.nids = 0;
+    
+    idlist.ids = G_malloc(idlist.nalloc * sizeof(int));
+
+    idlist.cellmax = ((CELL)1 << (sizeof(CELL) * 8 - 2)) - 1;
+    idlist.cellmax += ((CELL)1 << (sizeof(CELL) * 8 - 2));
+
+    return;
+}
+
+void add_free_id(int id)
+{
+    if (id <= 0)
+	return;
+
+    if (idlist.nalloc <= idlist.nids) {
+	idlist.nalloc = idlist.nids + 10;
+	idlist.ids = G_realloc(idlist.ids, idlist.nalloc * sizeof(int));
+    }
+    idlist.ids[idlist.nids++] = id;
+
+    return;
+}
+
+int get_free_id(struct globals *globals)
+{
+    if (idlist.nids > 0) {
+	idlist.nids--;
+
+	return idlist.ids[idlist.nids];
+    }
+
+    if (globals->max_rid == idlist.cellmax)
+	G_fatal_error(_("Too many objects: integer overflow"));
+
+    globals->max_rid++;
+    
+    return globals->max_rid;
+}
+
+void free_free_ids(void)
+{
+    if (idlist.nalloc) {
+	G_free(idlist.ids);
+	idlist.ids = NULL;
+	idlist.nalloc = 0;
+	idlist.nids = 0;
+    }
+    
+    return;
+}
+
 /* function used by binary tree to compare items */
 static int compare_rc(const void *first, const void *second)
 {
@@ -73,21 +137,39 @@ static int compare_ints(const void *first, const void *second)
 
 static int compare_double(double first, double second)
 {
-
-    /* standard comparison, gives good results */
     if (first < second)
 	return -1;
     return (first > second);
+}
 
-    /* fuzzy comparison, 
-     * can give weird results if EPSILON is too large or 
-     * if the formula is changed because this is operating at the 
-     * limit of double fp precision */
-    if (first < second && first + first * EPSILON < second)
-	    return -1;
-    if (first > second && first > second + second * EPSILON)
-	    return 1;
-    return 0;
+static int compare_sim_ngbrs(double simi, double simk, int candi, int candk,
+                             struct ngbr_stats *Ri, struct ngbr_stats *Rk)
+{
+    if (simi < simk)
+	return -1;
+
+    if (simi > simk)
+	return 1;
+
+    if (Rk->count == 0 || Ri->count < Rk->count)
+	return -1;
+    if (Ri->count > Rk->count)
+	return 1;
+
+    if (candi && !candk)
+	return -1;
+    
+    if (candk && !candi)
+	return 1;
+
+    if (Ri->row < Rk->row)
+	return -1;
+    if (Ri->row > Rk->row)
+	return 1;
+
+    if (Ri->col < Rk->col)
+	return -1;
+    return (Ri->col > Rk->col);
 }
 
 static int dump_Ri(struct ngbr_stats *Ri, struct reg_stats *Ri_rs, double *Ri_sim,
@@ -132,9 +214,16 @@ int region_growing(struct globals *globals)
     struct reg_stats Ri_rs, Rk_rs, Rk_bestn_rs;
     double *dp;
     struct NB_TREE *tmpnbtree;
+    CELL cellmax;
+    struct Cell_head cellhd;
 
     G_verbose_message("Running region growing algorithm");
 
+    cellmax = ((CELL)1 << (sizeof(CELL) * 8 - 2)) - 1;
+    cellmax += ((CELL)1 << (sizeof(CELL) * 8 - 2));
+
+    init_free_ids();
+
     /* init neighbor stats */
     Ri.mean = G_malloc(globals->datasize);
     Rk.mean = G_malloc(globals->datasize);
@@ -159,8 +248,10 @@ int region_growing(struct globals *globals)
     threshold = alpha2;
     G_debug(1, "Squared threshold: %g", threshold);
 
-    /* make the divisor a constant ? */
-    divisor = globals->nrows + globals->ncols;
+    Rast_get_cellhd(globals->Ref.file[0].name, globals->Ref.file[0].mapset, &cellhd);
+    divisor = cellhd.rows + cellhd.cols;
+
+    /* TODO: renumber seeds */
 
     while (t < globals->end_t && n_merges > 1) {
 
@@ -208,11 +299,6 @@ int region_growing(struct globals *globals)
 
 		/* get Ri's segment ID */
 		Segment_get(&globals->rid_seg, (void *)&Ri.id, Ri.row, Ri.col);
-		
-		if (Ri.id < 0)
-		    continue;
-		if (Ri.id == 0)
-		    G_fatal_error("Zero segment id at row %d, col %d", Ri.row, Ri.col);
 
 		/* find segment neighbors */
 		/* find Ri's best neighbor, clear candidate flag */
@@ -230,15 +316,21 @@ int region_growing(struct globals *globals)
 		Ri_nn = find_best_neighbor(&Ri, &Ri_rs, Ri_ngbrs,
 					   &Rk, &Rk_rs, &Ri_similarity,
 					   1, globals);
+
 		/* Rk is now complete */
 		G_debug(4, "Rk is now complete");
 
-		if (Rk.id == 0) {
+		if (Rk.id < 0) {
 		    /* this can only happen if the segment is surrounded by NULL data */
 		    G_debug(4, "Segment had no valid neighbors");
 		    continue;
 		}
 
+		if (compare_double(Ri_similarity, threshold) >= 0) {
+		    G_debug(4, "Best neighbor is not similar enough");
+		    continue;
+		}
+
 		if (/* !(t & 1) && */ Ri_nn == 1 &&
 		    !(FLAG_GET(globals->candidate_flag, Rk.row, Rk.col)) &&
 		    compare_double(Ri_similarity, threshold) == -1) {
@@ -264,16 +356,13 @@ int region_growing(struct globals *globals)
 
 		    pathflag = FALSE;
 		}
-		/* this is slow ??? */
-		if (/* t & */ 1) {
-		    if ((globals->nn < 8 && Rk.count <= 8) || 
-		        (globals->nn >= 8 && Rk.count <= globals->nn))
-		    candidates_only = FALSE;
-		}
 
 		while (pathflag) {
 		    pathflag = FALSE;
-		    
+
+		    if (Rk.count <= globals->nn || Rk.count <= globals->min_segment_size)
+			candidates_only = FALSE;
+
 		    /* optional check if Rk is candidate
 		     * to prevent backwards merging */
 		    if (candidates_only && 
@@ -331,11 +420,13 @@ int region_growing(struct globals *globals)
 			    Ri_nn -= Ri_ngbrs->count;
 			    Ri_nn += (Rk_nn - Rk_ngbrs->count);
 			    globals->ns.id = Rk.id;
+			    globals->ns.row = Rk.row;
+			    globals->ns.col = Rk.col;
 			    nbtree_remove(Ri_ngbrs, &(globals->ns));
 
 			    nbtree_init_trav(&travngbr, Rk_ngbrs);
 			    while ((next = nbtree_traverse(&travngbr))) {
-				if (!nbtree_find(Ri_ngbrs, next) && next->id != Ri.id)
+				if (!nbtree_find(Ri_ngbrs, next) && cmp_ngbr(next, &Ri) != 0)
 				    nbtree_insert(Ri_ngbrs, next);
 			    }
 			    nbtree_clear(Rk_ngbrs);
@@ -358,14 +449,14 @@ int region_growing(struct globals *globals)
 			    search_neighbors(&Ri, &Ri_rs, Ri_ngbrs, &Ri_similarity,
 					     &Rk, &Rk_rs, globals);
 
-			    if (Rk.id != 0 && Ri_nn > 0 && 
+			    if (Rk.id >= 0 && Ri_nn > 0 &&
 			        compare_double(Ri_similarity, threshold) == -1) {
 
 				pathflag = TRUE;
 				/* candidates_only:
-				 * FALSE: less passes, takes a bit longer, but less memory
-				 * TRUE: more passes, is a bit faster */
-				candidates_only = FALSE;
+				 * FALSE: less passes, slower, but less memory
+				 * TRUE: more passes but faster */
+				candidates_only = TRUE;
 			    }
 			    /* else end of Ri -> Rk chain since we merged Ri and Rk
 			     * go to next row, col */
@@ -381,16 +472,16 @@ int region_growing(struct globals *globals)
 			    if (Rk_nn < 2)
 				pathflag = FALSE;
 
-			    if (Rk.id < 1)
+			    if (Rk.id < 0)
 				pathflag = FALSE;
 
-			    if (Rk_bestn.id == 0) {
-				G_debug(4, "Rk's best neighour is zero");
+			    if (Rk_bestn.id < 0) {
+				G_debug(4, "Rk's best neighour is negative");
 				pathflag = FALSE;
 			    }
 
 			    if (pathflag) {
-				
+
 				/* clear candidate flag for Rk */
 				if (FLAG_GET(globals->candidate_flag, Rk.row, Rk.col)) {
 				    set_candidate_flag(&Rk, FALSE, globals);
@@ -409,7 +500,7 @@ int region_growing(struct globals *globals)
 
 				Ri_rs.id = Rk_rs.id;
 				Rk_rs.id = Rk_bestn_rs.id;
-				Rk_bestn_rs.id = 0;
+				Rk_bestn_rs.id = -1;
 				Ri_rs.count = Rk_rs.count;
 				Rk_rs.count = Rk_bestn_rs.count;
 				Rk_bestn_rs.count = 0;
@@ -446,6 +537,25 @@ int region_growing(struct globals *globals)
     else
 	G_message(_("Segmentation converged after %d iterations"), t);
 
+    /* assign region IDs to remaining 0 IDs */
+    G_message(_("Assigning region IDs to remaining single-cell regions..."));
+    for (row = globals->row_min; row < globals->row_max; row++) {
+	G_percent(row - globals->row_min,
+		  globals->row_max - globals->row_min, 4);
+	for (col = globals->col_min; col < globals->col_max; col++) {
+	    if (!(FLAG_GET(globals->null_flag, row, col))) {
+		/* get segment id */
+		Segment_get(&globals->rid_seg, (void *) &Ri.id, row, col);
+		if (Ri.id == 0) {
+		    Ri.id = get_free_id(globals);
+		    Segment_put(&globals->rid_seg, (void *) &Ri.id, row, col);
+		}
+	    }
+	}
+    }
+    G_percent(1, 1, 1);
+
+    free_free_ids();
 
     /* ****************************************************************************************** */
     /* final pass, ignore threshold and force a merge for small segments with their best neighbor */
@@ -518,7 +628,7 @@ int region_growing(struct globals *globals)
 		    Ri_nn = 0;
 		    Ri_similarity = 2;
 		    
-		    Rk.id = 0;
+		    Rk.id = -1;
 
 		    if (do_merge) {
 
@@ -530,7 +640,7 @@ int region_growing(struct globals *globals)
 		    }
 		    do_merge = 0;
 
-		    if (Rk.id != 0) {
+		    if (Rk.id >= 0) {
 			/* merge Ri with Rk */
 			/* do not clear candidate flag for Rk */
 			merge_regions(&Ri, &Ri_rs, &Rk, &Rk_rs, 0, globals);
@@ -547,7 +657,7 @@ int region_growing(struct globals *globals)
 	/* finished one pass for processing candidate pixels */
 	G_verbose_message("%d merges", n_merges);
     }
-    
+
     /* free neighbor stats */
     G_free(Ri.mean);
     G_free(Rk.mean);
@@ -585,15 +695,14 @@ static int find_best_neighbor(struct ngbr_stats *Ri,
     int neighbors[8][2];
     struct RB_TREE *no_check_tree;	/* cells already checked */
     struct reg_stats *rs_found;
+    int candk, candtmp;
 
     G_debug(4, "find_best_neighbor()");
 
     if (Ri->id != Ri_rs->id)
 	G_fatal_error("Ri = %d but Ri_rs = %d", Ri->id, Ri_rs->id);
-    if (Ri->id <= 0)
+    if (Ri->id < 0)
 	G_fatal_error("Ri is %d", Ri->id);
-    if (Ri_rs->id <= 0)
-	G_fatal_error("Ri_rs is %d", Ri_rs->id);
 
     /* dynamics of the region growing algorithm
      * some regions are growing fast, often surrounded by many small regions
@@ -609,8 +718,9 @@ static int find_best_neighbor(struct ngbr_stats *Ri,
     nbtree_clear(Ri_ngbrs);
     n_ngbrs = 0;
     /* TODO: add size of largest region to reg_tree, use this as min */
-    Rk->count = globals->ncells + 1;
-    Rk->id = Rk_rs->id = 0;
+    Rk->count = Rk_rs->count = 0;
+    Rk->id = Rk_rs->id = -1;
+    candk = 0;
 
     /* go through segment, spreading outwards from head */
     rclist_init(&rilist);
@@ -659,7 +769,7 @@ static int find_best_neighbor(struct ngbr_stats *Ri,
 				    (void *) &(globals->ns.id),
 				    ngbr_rc.row, ngbr_rc.col);
 
-			if (globals->ns.id == Ri->id) {
+			if (Ri->id > 0 && globals->ns.id == Ri->id) {
 
 			    /* want to check this neighbor's neighbors */
 			    rclist_add(&rilist, ngbr_rc.row, ngbr_rc.col);
@@ -684,10 +794,13 @@ static int find_best_neighbor(struct ngbr_stats *Ri,
 				/* globals->ns is now complete */
 
 				tempsim = (globals->calculate_similarity)(Ri, &globals->ns, globals);
+				candtmp = (FLAG_GET(globals->candidate_flag, ngbr_rc.row, ngbr_rc.col)) != 0;
+
+				cmp = compare_sim_ngbrs(tempsim, *sim, candtmp, candk, &globals->ns, Rk);
 
-				cmp = compare_double(tempsim, *sim);
 				if (cmp == -1) {
 				    *sim = tempsim;
+				    candk = candtmp;
 				    /* copy temp Rk to Rk */
 				    Rk->row = ngbr_rc.row;
 				    Rk->col = ngbr_rc.col;
@@ -704,27 +817,6 @@ static int find_best_neighbor(struct ngbr_stats *Ri,
 				    memcpy(Rk_rs->sum, rs_found->sum,
 				           globals->datasize);
 				}
-				else if (cmp == 0) {
-				    /* resolve ties: prefer smaller regions */
-
-				    if (Rk->count > globals->ns.count) {
-					/* copy temp Rk to Rk */
-					Rk->row = ngbr_rc.row;
-					Rk->col = ngbr_rc.col;
-
-					Rk->id = rs_found->id;
-					Rk->count = rs_found->count;
-					memcpy(Rk->mean, rs_found->mean,
-					       globals->datasize);
-
-					Rk_rs->id = Rk->id;
-					Rk_rs->count = Rk->count;
-					memcpy(Rk_rs->mean, rs_found->mean,
-					       globals->datasize);
-					memcpy(Rk_rs->sum, rs_found->sum,
-					       globals->datasize);
-				    }
-				}
 
 				n_ngbrs++;
 				nbtree_insert(Ri_ngbrs, &globals->ns);
@@ -738,10 +830,12 @@ static int find_best_neighbor(struct ngbr_stats *Ri,
 
     /* clean up */
     rbtree_destroy(no_check_tree);
+    rclist_destroy(&rilist);
 
     return n_ngbrs;
 }
 
+#ifdef _OR_SHAPE_
 double calculate_shape(struct reg_stats *rsi, struct reg_stats *rsk,
                        int nshared, struct globals *globals)
 {
@@ -817,7 +911,7 @@ double calculate_shape(struct reg_stats *rsi, struct reg_stats *rsk,
 
     return globals->smooth_weight * smooth + (1 - globals->smooth_weight) * compact;
 }
-
+#endif
 
 static int search_neighbors(struct ngbr_stats *Ri,
 			    struct reg_stats *Ri_rs,
@@ -830,7 +924,7 @@ static int search_neighbors(struct ngbr_stats *Ri,
     double tempsim, *dp;
     struct NB_TRAV travngbr;
     struct ngbr_stats *next;
-    int cmp;
+    int cmp, candk, candtmp;
 
     G_debug(4, "search_neighbors");
 
@@ -842,36 +936,29 @@ static int search_neighbors(struct ngbr_stats *Ri,
 	G_fatal_error("Ri_rs is %d", Ri_rs->id);
 
     nbtree_init_trav(&travngbr, Ri_ngbrs);
-    Rk->count = globals->ncells + 1;
-    Rk->id = Rk_rs->id = 0;
+    Rk->count = 0;
+    Rk->id = Rk_rs->id = -1;
+    candk = 0;
 
     while ((next = nbtree_traverse(&travngbr))) {
 	tempsim = (globals->calculate_similarity)(Ri, next, globals);
+	candtmp = (FLAG_GET(globals->candidate_flag, next->row, next->col)) != 0;
+
+	cmp = compare_sim_ngbrs(tempsim, *sim, candtmp, candk, next, Rk);
 
-	cmp = compare_double(tempsim, *sim);
 	if (cmp == -1) {
 	    *sim = tempsim;
+	    candk = candtmp;
 
 	    dp = Rk->mean;
 	    *Rk = *next;
 	    Rk->mean = dp;
 	    memcpy(Rk->mean, next->mean, globals->datasize);
 	}
-	else if (cmp == 0) {
-	    /* resolve ties, prefer smaller regions */
-	    G_debug(4, "resolve ties");
-
-	    if (Rk->count > next->count) {
-		dp = Rk->mean;
-		*Rk = *next;
-		Rk->mean = dp;
-		memcpy(Rk->mean, next->mean, globals->datasize);
-	    }
-	}
     }
     Rk_rs->id = Rk->id;
 
-    if (Rk->id != 0) {
+    if (Rk->id >= 0) {
 	fetch_reg_stats(Rk->row, Rk->col, Rk_rs, globals);
     }
 
@@ -899,6 +986,10 @@ int update_band_vals(int row, int col, struct reg_stats *rs,
 	G_fatal_error(_("Region ids are different"));
     }
 
+    if (rs->id < 1) {
+	G_fatal_error(_("Region id %d is invalid"), rs->id);
+    }
+
     if (rs->count == 1) {
 	G_warning(_("Region consists of only one cell, nothing to update"));
 	return rs->count;
@@ -1032,9 +1123,10 @@ static int merge_regions(struct ngbr_stats *Ri, struct reg_stats *Ri_rs,
     G_debug(4, "merge_regions");
 
     /* Ri ID must always be positive */
-    if (Ri_rs->id < 1)
-	G_fatal_error("Ri id is not positive: %d", Ri_rs->id);
-    /* if Rk ID is negative (no seed), Rk count must be 1  */
+    if (Ri_rs->id < 1 && Ri_rs->count > 1)
+	G_fatal_error("Ri id is not positive: %d, but count is > 1: %d",
+	              Ri_rs->id, Ri_rs->count);
+    /* if Rk ID is zero (no seed), Rk count must be 1  */
     if (Rk_rs->id < 1 && Rk_rs->count > 1)
 	G_fatal_error("Rk id is not positive: %d, but count is > 1: %d",
 	              Rk_rs->id, Rk_rs->count);
@@ -1058,6 +1150,12 @@ static int merge_regions(struct ngbr_stats *Ri, struct reg_stats *Ri_rs,
     } while (n--);
 
     if (Ri->count >= Rk->count) {
+	
+	if (Ri->id == 0) {
+	    Ri->id = get_free_id(globals);
+	    Ri_rs->id = Ri->id;
+	    Segment_put(&globals->rid_seg, (void *) &Ri->id, Ri->row, Ri->col);
+	}
 
 	if (Rk->count >= globals->min_reg_size) {
 	    if (rgtree_find(globals->reg_tree, Rk_rs) == NULL)
@@ -1065,6 +1163,7 @@ static int merge_regions(struct ngbr_stats *Ri, struct reg_stats *Ri_rs,
 	    /* remove from tree */
 	    rgtree_remove(globals->reg_tree, Rk_rs);
 	}
+	add_free_id(Rk->id);
     }
     else {
 
@@ -1074,6 +1173,7 @@ static int merge_regions(struct ngbr_stats *Ri, struct reg_stats *Ri_rs,
 	    /* remove from tree */
 	    rgtree_remove(globals->reg_tree, Ri_rs);
 	}
+	add_free_id(Ri->id);
 
 	/* magic switch */
 	Ri_rs->id = Rk->id;
@@ -1093,7 +1193,19 @@ static int merge_regions(struct ngbr_stats *Ri, struct reg_stats *Ri_rs,
     Ri->count = Ri_rs->count;
     memcpy(Ri->mean, Ri_rs->mean, globals->datasize);
 
-    if (Ri->id == Ri_rs->id) {
+    if (Rk->id == 0) {
+	/* the actual merge: change region id */
+	Segment_put(&globals->rid_seg, (void *) &Ri->id, Rk->row, Rk->col);
+
+	if (do_cand) {
+	    if (FLAG_GET(globals->candidate_flag, Rk->row, Rk->col)) {
+		/* clear candidate flag */
+		FLAG_UNSET(globals->candidate_flag, Rk->row, Rk->col);
+		globals->candidate_count--;
+	    }
+	}
+    }
+    else if (Ri->id == Ri_rs->id) {
 	/* Ri is already updated, including candidate flags
 	 * need to clear candidate flag for Rk and set new id */
 	 
@@ -1111,7 +1223,8 @@ static int merge_regions(struct ngbr_stats *Ri, struct reg_stats *Ri_rs,
 	}
 
 	rclist_init(&rlist);
-	rclist_add(&rlist, Rk->row, Rk->col);
+	if (Rk->count > 1)
+	    rclist_add(&rlist, Rk->row, Rk->col);
 
 	while (rclist_drop(&rlist, &next)) {
 
@@ -1137,9 +1250,9 @@ static int merge_regions(struct ngbr_stats *Ri, struct reg_stats *Ri_rs,
 		    if (!(FLAG_GET(globals->null_flag, ngbr_rc.row, ngbr_rc.col))) {
 
 			Segment_get(&globals->rid_seg, (void *) &R_id,
-			            ngbr_rc.row, ngbr_rc.col);
+				    ngbr_rc.row, ngbr_rc.col);
 
-			if (R_id == Rk->id) {
+			if (Rk->id > 0 && R_id == Rk->id) {
 			    /* the actual merge: change region id */
 			    Segment_put(&globals->rid_seg, (void *) &Ri->id, ngbr_rc.row, ngbr_rc.col);
 
@@ -1166,7 +1279,8 @@ static int merge_regions(struct ngbr_stats *Ri, struct reg_stats *Ri_rs,
 	Segment_put(&globals->rid_seg, (void *) &Rk->id, Ri->row, Ri->col);
 
 	rclist_init(&rlist);
-	rclist_add(&rlist, Ri->row, Ri->col);
+	if (Ri->count > 1)
+	    rclist_add(&rlist, Ri->row, Ri->col);
 
 	while (rclist_drop(&rlist, &next)) {
 
@@ -1188,7 +1302,7 @@ static int merge_regions(struct ngbr_stats *Ri, struct reg_stats *Ri_rs,
 
 			Segment_get(&globals->rid_seg, (void *) &R_id, ngbr_rc.row, ngbr_rc.col);
 
-			if (R_id == Ri->id) {
+			if (Ri->id > 0 && R_id == Ri->id) {
 			    /* the actual merge: change region id */
 			    Segment_put(&globals->rid_seg, (void *) &Rk->id, ngbr_rc.row, ngbr_rc.col);
 
@@ -1205,12 +1319,9 @@ static int merge_regions(struct ngbr_stats *Ri, struct reg_stats *Ri_rs,
 	if (Ri->id != Rk->id)
 	    G_fatal_error("Ri ID should be set to Rk ID");
     }
-    
-    if (Rk->id > 0)
-	globals->n_regions--;
 
     /* disable Rk */
-    Rk->id = Rk_rs->id = 0;
+    Rk->id = Rk_rs->id = -1;
     Rk->count = Rk_rs->count = 0;
     
     /* update Ri */
@@ -1232,14 +1343,11 @@ static int set_candidate_flag(struct ngbr_stats *head, int value, struct globals
 
     G_debug(4, "set_candidate_flag");
 
-    if (!(FLAG_GET(globals->candidate_flag, head->row, head->col)) != value) {
+    if ((!(FLAG_GET(globals->candidate_flag, head->row, head->col))) != value) {
 	G_warning(_("Candidate flag is already %s"), value ? _("set") : _("unset"));
 	return FALSE;
     }
 
-    rclist_init(&rlist);
-    rclist_add(&rlist, head->row, head->col);
-
     /* (un)set candidate flag */
     if (value == TRUE) {
 	FLAG_SET(globals->candidate_flag, head->row, head->col);
@@ -1249,6 +1357,12 @@ static int set_candidate_flag(struct ngbr_stats *head, int value, struct globals
 	FLAG_UNSET(globals->candidate_flag, head->row, head->col);
 	globals->candidate_count--;
     }
+    
+    if (head->id == 0)
+	return TRUE;
+
+    rclist_init(&rlist);
+    rclist_add(&rlist, head->row, head->col);
 
     while (rclist_drop(&rlist, &next)) {
 
@@ -1267,7 +1381,7 @@ static int set_candidate_flag(struct ngbr_stats *head, int value, struct globals
 
 		if (!(FLAG_GET(globals->null_flag, ngbr_rc.row, ngbr_rc.col))) {
 
-		    if (!(FLAG_GET(globals->candidate_flag, ngbr_rc.row, ngbr_rc.col)) == value) {
+		    if ((!(FLAG_GET(globals->candidate_flag, ngbr_rc.row, ngbr_rc.col))) == value) {
 
 			Segment_get(&globals->rid_seg, (void *) &R_id, ngbr_rc.row, ngbr_rc.col);
 
@@ -1290,6 +1404,7 @@ static int set_candidate_flag(struct ngbr_stats *head, int value, struct globals
 	    }
 	} while (n--);
     }
+    rclist_destroy(&rlist);
 
     return TRUE;
 }
@@ -1298,11 +1413,11 @@ int fetch_reg_stats(int row, int col, struct reg_stats *rs,
                            struct globals *globals)
 {
     struct reg_stats *rs_found;
-    
-    if (rs->id <= 0)
+
+    if (rs->id < 0)
 	G_fatal_error("fetch_reg_stats(): invalid region id %d", rs->id);
 
-    if ((rs_found = rgtree_find(globals->reg_tree, rs)) != NULL) {
+    if (rs->id > 0 && (rs_found = rgtree_find(globals->reg_tree, rs)) != NULL) {
 
 	memcpy(rs->mean, rs_found->mean, globals->datasize);
 	memcpy(rs->sum, rs_found->sum, globals->datasize);
@@ -1323,7 +1438,7 @@ static int calculate_reg_stats(int row, int col, struct reg_stats *rs,
 
     G_debug(4, "calculate_reg_stats()");
 
-    if (rs->id <= 0)
+    if (rs->id < 0)
 	G_fatal_error("Invalid region id %d", rs->id);
 
     Segment_get(&globals->bands_seg, (void *)globals->bands_val,
@@ -1331,6 +1446,12 @@ static int calculate_reg_stats(int row, int col, struct reg_stats *rs,
     rs->count = 1;
     memcpy(rs->sum, globals->bands_val, globals->datasize);
 
+    if (rs->id == 0) {
+	memcpy(rs->mean, rs->sum, globals->datasize);
+
+	return 1;
+    }
+
     if (globals->min_reg_size < 3)
 	ret = 1;
     else if (globals->min_reg_size == 3) {

imagery/i.segment/regtree.c

@@ -77,6 +77,7 @@ struct RG_TREE *rgtree_create(int nbands, size_t rb_datasize)
 int rgtree_insert(struct RG_TREE *tree, struct reg_stats *data)
 {
     assert(tree && data);
+    assert(data->id > 0);
 
     if (tree->root == NULL) {
 	/* create a new root node for tree */

imagery/i.segment/write_output.c

@@ -1,6 +1,3 @@
-/* transfer the segmented regions from the segmented data file to a raster file */
-/* close_files() function is at bottom */
-
 #include <stdlib.h>
 #include <grass/gis.h>
 #include <grass/raster.h>
@@ -9,11 +6,7 @@
 #include <grass/glocale.h>
 #include "iseg.h"
 
-/* TODO some time delay here with meanbuf, etc being processed.  I only put if statements on the actual files
- * to try and keep the code more clear.  Need to see if this raster makes stats processing easier?  Or IFDEF it out?
- */
-
-int write_output(struct globals *globals)
+int write_ids(struct globals *globals)
 {
     int out_fd, row, col, maxid;
     CELL *outbuf, rid;
@@ -64,167 +57,230 @@ int write_output(struct globals *globals)
     Rast_command_history(&hist);
     Rast_write_history(globals->out_name, &hist);
 
-    /* write goodness of fit */
-    if (globals->out_band) {
-	int mean_fd;
-	FCELL *meanbuf;
-	double thresh, maxdev, sim, mingood;
-	struct ngbr_stats Ri, Rk;
-	struct Ref Ref;		/* group reference list */
-	DCELL **inbuf;		/* buffers to store lines from each of the imagery group rasters */
-	int n, *in_fd;
-	struct FPRange *fp_range;	/* min/max values of each input raster */
-	DCELL *min, *max;
-
-	mean_fd = Rast_open_new(globals->out_band, FCELL_TYPE);
-	meanbuf = Rast_allocate_f_buf();
-
-	/* goodness of fit for each cell: 1 = good fit, 0 = bad fit */
-	/* similarity of each cell to region mean
-	 * max possible difference: globals->threshold
-	 * if similarity < globals->alpha * globals->alpha * globals->threshold
-	 * 1
-	 * else 
-	 * (similarity - globals->alpha * globals->alpha * globals->threshold) /
-	 * (globals->threshold * (1 - globals->alpha * globals->alpha) */
-
-	thresh = globals->alpha * globals->alpha * globals->max_diff;
-	maxdev = globals->max_diff * (1 - globals->alpha * globals->alpha);
-	mingood = 1;
-
-	/* open input bands */
-	if (!I_get_group_ref(globals->image_group, &Ref))
-	    G_fatal_error(_("Group <%s> not found in the current mapset"),
-			  globals->image_group);
-	if (Ref.nfiles <= 0)
-	    G_fatal_error(_("Group <%s> contains no raster maps"),
-			  globals->image_group);
-
-	in_fd = G_malloc(Ref.nfiles * sizeof(int));
-	inbuf = (DCELL **) G_malloc(Ref.nfiles * sizeof(DCELL *));
-	fp_range = G_malloc(Ref.nfiles * sizeof(struct FPRange));
-	min = G_malloc(Ref.nfiles * sizeof(DCELL));
-	max = G_malloc(Ref.nfiles * sizeof(DCELL));
-
-	G_debug(1, "Opening input rasters...");
-	for (n = 0; n < Ref.nfiles; n++) {
-	    inbuf[n] = Rast_allocate_d_buf();
-	    in_fd[n] = Rast_open_old(Ref.file[n].name, Ref.file[n].mapset);
-
-	    /* returns -1 on error, 2 on empty range, quitting either way. */
-	    if (Rast_read_fp_range(Ref.file[n].name, Ref.file[n].mapset, &fp_range[n]) != 1)
-		G_fatal_error(_("No min/max found in raster map <%s>"),
-			      Ref.file[n].name);
-	    Rast_get_fp_range_min_max(&(fp_range[n]), &min[n], &max[n]);
-
-	    G_debug(1, "Range for layer %d: min = %f, max = %f",
-			n, min[n], max[n]);
-	}
+    /* free memory */
+    Rast_free_colors(&colors);
+
+    return TRUE;
+}
 
-	G_message(_("Writing out goodness of fit"));
-	for (row = 0; row < globals->nrows; row++) {
+/* write goodness of fit */
+int write_gof_rg(struct globals *globals)
+{
+    int row, col;
+    int mean_fd;
+    CELL rid;
+    FCELL *meanbuf;
+    double thresh, maxdev, sim, mingood;
+    struct ngbr_stats Ri, Rk;
+    DCELL **inbuf;		/* buffers to store lines from each of the imagery group rasters */
+    int n, *in_fd;
+    struct FPRange *fp_range;	/* min/max values of each input raster */
+    struct Colors colors;
+    struct History hist;
+    DCELL *min, *max;
+
+    mean_fd = Rast_open_new(globals->out_band, FCELL_TYPE);
+    meanbuf = Rast_allocate_f_buf();
+
+    /* goodness of fit for each cell: 1 = good fit, 0 = bad fit */
+    /* similarity of each cell to region mean
+     * max possible difference: globals->threshold
+     * if similarity < globals->alpha * globals->alpha * globals->threshold
+     * 1
+     * else 
+     * (similarity - globals->alpha * globals->alpha * globals->threshold) /
+     * (globals->threshold * (1 - globals->alpha * globals->alpha) */
+
+    thresh = globals->alpha * globals->alpha * globals->max_diff;
+    maxdev = globals->max_diff * (1 - globals->alpha * globals->alpha);
+    mingood = 1;
+
+    /* open input bands */
+    in_fd = G_malloc(globals->Ref.nfiles * sizeof(int));
+    inbuf = (DCELL **) G_malloc(globals->Ref.nfiles * sizeof(DCELL *));
+    fp_range = G_malloc(globals->Ref.nfiles * sizeof(struct FPRange));
+    min = G_malloc(globals->Ref.nfiles * sizeof(DCELL));
+    max = G_malloc(globals->Ref.nfiles * sizeof(DCELL));
+
+    G_debug(1, "Opening input rasters...");
+    for (n = 0; n < globals->Ref.nfiles; n++) {
+	inbuf[n] = Rast_allocate_d_buf();
+	in_fd[n] = Rast_open_old(globals->Ref.file[n].name, globals->Ref.file[n].mapset);
+
+	/* returns -1 on error, 2 on empty range, quitting either way. */
+	if (Rast_read_fp_range(globals->Ref.file[n].name, globals->Ref.file[n].mapset, &fp_range[n]) != 1)
+	    G_fatal_error(_("No min/max found in raster map <%s>"),
+			  globals->Ref.file[n].name);
+	Rast_get_fp_range_min_max(&(fp_range[n]), &min[n], &max[n]);
+
+	G_debug(1, "Range for layer %d: min = %f, max = %f",
+		    n, min[n], max[n]);
+    }
 
-	    G_percent(row, globals->nrows, 9);
+    G_message(_("Writing out goodness of fit"));
+    for (row = 0; row < globals->nrows; row++) {
 
-	    Rast_set_f_null_value(meanbuf, globals->ncols);
+	G_percent(row, globals->nrows, 9);
 
-	    for (n = 0; n < Ref.nfiles; n++) {
-		Rast_get_d_row(in_fd[n], inbuf[n], row);
-	    }
+	Rast_set_f_null_value(meanbuf, globals->ncols);
 
-	    for (col = 0; col < globals->ncols; col++) {
+	for (n = 0; n < globals->Ref.nfiles; n++) {
+	    Rast_get_d_row(in_fd[n], inbuf[n], row);
+	}
 
-		if (!(FLAG_GET(globals->null_flag, row, col))) {
+	for (col = 0; col < globals->ncols; col++) {
+
+	    if (!(FLAG_GET(globals->null_flag, row, col))) {
+		
+		Segment_get(&globals->rid_seg, (void *) &rid, row, col);
+
+		if (rid > 0) {
 		    
-		    Segment_get(&globals->rid_seg, (void *) &rid, row, col);
-
-		    if (rid > 0) {
-			
-			Ri.row = Rk.row = row;
-			Ri.col = Rk.col = col;
-
-			/* get values for Ri = this region */
-			globals->rs.id = rid;
-			fetch_reg_stats(row, col, &globals->rs, globals);
-			Ri.mean = globals->rs.mean;
-			Ri.count = globals->rs.count; 
-
-			sim = 0.;
-			/* region consists of more than one cell */
-			if (Ri.count > 1) {
-
-			    /* get values for Rk = this cell */
-			    for (n = 0; n < Ref.nfiles; n++) {
-				if (globals->weighted == FALSE)
-				    /* scaled version */
-				    globals->second_val[n] = (inbuf[n][col] - min[n]) / (max[n] - min[n]);
-				else
-				    globals->second_val[n] = inbuf[n][col];
-			    }
-
-			    Rk.mean = globals->second_val;
-
-			    /* calculate similarity */
-			    sim = (*globals->calculate_similarity) (&Ri, &Rk, globals);
-			}
-			
-			if (0) {
-			    if (sim < thresh)
-				meanbuf[col] = 1;
-			    else {
-				sim = 1. - (sim - thresh) / maxdev;
-				meanbuf[col] = sim;
-				if (mingood > sim)
-				    mingood = sim;
-			    }
+		    Ri.row = Rk.row = row;
+		    Ri.col = Rk.col = col;
+
+		    /* get values for Ri = this region */
+		    globals->rs.id = rid;
+		    fetch_reg_stats(row, col, &globals->rs, globals);
+		    Ri.mean = globals->rs.mean;
+		    Ri.count = globals->rs.count; 
+
+		    sim = 0.;
+		    /* region consists of more than one cell */
+		    if (Ri.count > 1) {
+
+			/* get values for Rk = this cell */
+			for (n = 0; n < globals->Ref.nfiles; n++) {
+			    if (globals->weighted == FALSE)
+				/* scaled version */
+				globals->second_val[n] = (inbuf[n][col] - min[n]) / (max[n] - min[n]);
+			    else
+				globals->second_val[n] = inbuf[n][col];
 			}
+
+			Rk.mean = globals->second_val;
+
+			/* calculate similarity */
+			sim = (*globals->calculate_similarity) (&Ri, &Rk, globals);
+		    }
+		    
+		    if (0) {
+			if (sim < thresh)
+			    meanbuf[col] = 1;
 			else {
-			    sim = 1 - sim;
+			    sim = 1. - (sim - thresh) / maxdev;
 			    meanbuf[col] = sim;
 			    if (mingood > sim)
 				mingood = sim;
 			}
 		    }
+		    else {
+			sim = 1 - sim;
+			meanbuf[col] = sim;
+			if (mingood > sim)
+			    mingood = sim;
+		    }
 		}
 	    }
-	    Rast_put_row(mean_fd, meanbuf, FCELL_TYPE);
 	}
+	Rast_put_row(mean_fd, meanbuf, FCELL_TYPE);
+    }
 
-	Rast_close(mean_fd);
+    Rast_close(mean_fd);
 
-	Rast_init_colors(&colors);
-	Rast_make_grey_scale_fp_colors(&colors, mingood, 1);
-	Rast_write_colors(globals->out_band, G_mapset(), &colors);
+    Rast_init_colors(&colors);
+    Rast_make_grey_scale_fp_colors(&colors, mingood, 1);
+    Rast_write_colors(globals->out_band, G_mapset(), &colors);
 
-	Rast_short_history(globals->out_band, "raster", &hist);
-	Rast_command_history(&hist);
-	Rast_write_history(globals->out_band, &hist);
+    Rast_short_history(globals->out_band, "raster", &hist);
+    Rast_command_history(&hist);
+    Rast_write_history(globals->out_band, &hist);
+
+    G_free(meanbuf);
+
+    G_debug(1, "Closing input rasters...");
+    for (n = 0; n < globals->Ref.nfiles; n++) {
+	Rast_close(in_fd[n]);
+	G_free(inbuf[n]);
+    }
+
+    G_free(inbuf);
+    G_free(in_fd);
+    G_free(fp_range);
+    G_free(min);
+    G_free(max);
+
+    return TRUE;
+}
 
-	G_free(meanbuf);
+int write_bands_ms(struct globals *globals)
+{
+    int *out_fd, row, col, n;
+    DCELL **outbuf;
+    char **name;
+    struct Colors colors;
+    struct History hist;
+    struct ngbr_stats Rout;
+
+    out_fd = G_malloc(sizeof(int) * globals->nbands);
+    name = G_malloc(sizeof(char *) * globals->nbands);
+    outbuf = G_malloc(sizeof(DCELL) * globals->nbands);
+    for (n = 0; n < globals->nbands; n++) {
+	outbuf[n] = Rast_allocate_d_buf();
+	G_asprintf(&name[n], "%s%s", globals->Ref.file[n].name, globals->ms_suf);
+	out_fd[n] = Rast_open_new(name[n], DCELL_TYPE);
+    }
+    
+    Rout.mean = G_malloc(globals->datasize);
+
+
+    for (row = 0; row < globals->nrows; row++) {
 
-	G_debug(1, "Closing input rasters...");
-	for (n = 0; n < Ref.nfiles; n++) {
-	    Rast_close(in_fd[n]);
-	    G_free(inbuf[n]);
+	G_percent(row, globals->nrows, 9);
+
+	for (n = 0; n < globals->nbands; n++)
+	    Rast_set_d_null_value(outbuf[n], globals->ncols);
+	for (col = 0; col < globals->ncols; col++) {
+
+	    if (!(FLAG_GET(globals->null_flag, row, col))) {
+		Segment_get(globals->bands_out, (void *) Rout.mean, row, col);
+		
+		for (n = 0; n < globals->nbands; n++) {
+		    outbuf[n][col] = Rout.mean[n];
+		    if (globals->weighted == FALSE)
+			outbuf[n][col] = Rout.mean[n] * (globals->max[n] - globals->min[n]) + globals->min[n];
+		}
+	    }
 	}
-	G_free(inbuf);
-	G_free(in_fd);
-	G_free(fp_range);
-	G_free(min);
-	G_free(max);
+	for (n = 0; n < globals->nbands; n++)
+	    Rast_put_row(out_fd[n], outbuf[n], DCELL_TYPE);
     }
 
-    /* free memory */
-    Rast_free_colors(&colors);
+    for (n = 0; n < globals->nbands; n++) {
+	Rast_close(out_fd[n]);
+
+	Rast_read_colors(globals->Ref.file[n].name, globals->Ref.file[n].mapset, &colors);
+	Rast_write_colors(name[n], G_mapset(), &colors);
+
+	Rast_short_history(name[n], "raster", &hist);
+	Rast_command_history(&hist);
+	Rast_write_history(name[n], &hist);
+    }
+
+    /* free */
 
     return TRUE;
 }
 
 int close_files(struct globals *globals)
 {
+    G_debug(1, "Closing input rasters...");
+
     /* close segmentation files and output raster */
     G_debug(1, "closing files");
     Segment_close(&globals->bands_seg);
+    if (globals->method == ORM_MS)
+	Segment_close(&globals->bands_seg2);
     if (globals->bounds_map)
 	Segment_close(&globals->bounds_seg);