[GRASS-SVN] r53276 - grass-addons/grass7/imagery/i.segment.xl

svn_grass at osgeo.org svn_grass at osgeo.org
Thu Sep 27 08:44:17 PDT 2012


Author: mmetz
Date: 2012-09-27 08:44:17 -0700 (Thu, 27 Sep 2012)
New Revision: 53276

Modified:
   grass-addons/grass7/imagery/i.segment.xl/create_isegs.c
   grass-addons/grass7/imagery/i.segment.xl/iseg.h
   grass-addons/grass7/imagery/i.segment.xl/open_files.c
   grass-addons/grass7/imagery/i.segment.xl/rclist.c
   grass-addons/grass7/imagery/i.segment.xl/write_output.c
Log:
i.segment.xl: add control for memory consumption

Modified: grass-addons/grass7/imagery/i.segment.xl/create_isegs.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/create_isegs.c	2012-09-27 11:33:21 UTC (rev 53275)
+++ grass-addons/grass7/imagery/i.segment.xl/create_isegs.c	2012-09-27 15:44:17 UTC (rev 53276)
@@ -14,23 +14,27 @@
 #include <grass/rbtree.h>	/* Red Black Tree library functions */
 #include "iseg.h"
 
-#define EPSILON 1.0e-16
+#define EPSILON 1.0e-12
 
 /* internal functions */
-static int merge_regions(struct ngbr_stats *,    /* Ri */
-                         struct reg_stats *,
-                         struct ngbr_stats *,    /* Rk */
-			 struct reg_stats *,
+static int merge_regions(struct ngbr_stats *, struct reg_stats *, /* Ri */
+                         struct ngbr_stats *, struct reg_stats *, /* Rk */
+			 int,
 			 struct globals *);
 static int search_neighbors(struct ngbr_stats *,    /* Ri */
+			    struct reg_stats *,
                             struct NB_TREE *,       /* Ri's neighbors */ 
 			    double *,               /* highest similarity */ 
 			    struct ngbr_stats *,    /* Ri's best neighbor */
+			    struct reg_stats *,
 		            struct globals *);
 static int set_candidate_flag(struct ngbr_stats *, int , struct globals *);
-static int find_best_neighbor(struct ngbr_stats *, struct NB_TREE *,
-			      struct reg_stats **, struct ngbr_stats *,
-			      double *, int, struct globals *);
+static int find_best_neighbor(struct ngbr_stats *, struct reg_stats *,
+			      struct NB_TREE *, struct ngbr_stats *,
+			      struct reg_stats *, double *, int,
+			      struct globals *);
+static int calculate_reg_stats(int, int, struct reg_stats *, 
+                         struct globals *);
 
 /* function used by binary tree to compare items */
 static int compare_rc(const void *first, const void *second)
@@ -60,6 +64,7 @@
 
 static int compare_double(double first, double second)
 {
+
     /* standard comparison, gives good results */
     if (first < second)
 	return -1;
@@ -69,11 +74,35 @@
      * 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)
-	return ((second > (first + first * EPSILON)) ? -1 : 0);
-    return (first > (second + second * EPSILON));
+    if (first < second && first + first * EPSILON < second)
+	    return -1;
+    if (first > second && first > second + second * EPSILON)
+	    return 1;
+    return 0;
 }
 
+static int dump_Ri(struct ngbr_stats *Ri, struct reg_stats *Ri_rs, double *Ri_sim,
+                   double *Rk_sim, int *Ri_nn, int *Rk_nn, struct globals *globals)
+{
+    int i;
+
+    G_debug(0, "Ri, Ri_rs ID: %d, %d", Ri->id, Ri_rs->id);
+    G_debug(0, "Ri, Ri_rs count: %d, %d", Ri->count, Ri_rs->count);
+
+    for (i = 0; i < globals->nbands; i++) {
+	G_debug(0, "Ri, Ri_rs mean %d: %g, %g", i, Ri->mean[i], Ri_rs->mean[i]);
+	G_debug(0, "Ri_rs sum %d: %g", i, Ri_rs->sum[i]);
+    }
+    G_debug(0, "Ri nn: %d", *Ri_nn);
+    if (Rk_nn)
+	G_debug(0, "Rk nn: %d", *Rk_nn);
+    G_debug(0, "Ri similarity: %g", *Ri_sim);
+    if (Rk_sim)
+	G_debug(0, "Rk similarity: %g", *Rk_sim);
+
+    return 1;
+}
+
 int create_isegs(struct globals *globals)
 {
     int row, col;
@@ -152,40 +181,44 @@
     int Ri_nn, Rk_nn; /* number of neighbors for Ri/Rk */
     struct NB_TREE *Ri_ngbrs, *Rk_ngbrs;
     struct NB_TRAV travngbr;
-    struct reg_stats *Ri_rs, *Rk_rs, *Rk_bestn_rs;
-    int verbose = (G_verbose() >= 3);
+    /* not all regions are in the tree, but we always need reg_stats for Ri and Rk */
+    struct reg_stats Ri_rs, Rk_rs, Rk_bestn_rs;
     double *dp;
     struct NB_TREE *tmpnbtree;
 
     G_verbose_message("Running region growing algorithm");
-    
+
+    /* init neighbor stats */
     Ri.mean = G_malloc(globals->datasize);
     Rk.mean = G_malloc(globals->datasize);
     Rk_bestn.mean = G_malloc(globals->datasize);
-    
+
     Ri_ngbrs = nbtree_create(globals->nbands, globals->datasize);
     Rk_ngbrs = nbtree_create(globals->nbands, globals->datasize);
 
+    /* init region stats */
+    Ri_rs.mean = G_malloc(globals->datasize);
+    Ri_rs.sum = G_malloc(globals->datasize);
+    Rk_rs.mean = G_malloc(globals->datasize);
+    Rk_rs.sum = G_malloc(globals->datasize);
+    Rk_bestn_rs.mean = G_malloc(globals->datasize);
+    Rk_bestn_rs.sum = G_malloc(globals->datasize);
+    
     t = 0;
     n_merges = 1;
 
     /* threshold calculation */
     alpha2 = globals->alpha * globals->alpha;
+    threshold = alpha2 * globals->threshold;
+    G_debug(1, "Squared threshold: %g", threshold);
+
     /* make the divisor a constant ? */
     divisor = globals->nrows + globals->ncols;
 
-    while (t < globals->end_t && n_merges > 0) {
+    while (t++ < globals->end_t && n_merges > 0) {
 
-	/* optional threshold as a function of t. does not make sense */
-	threshold = alpha2 * globals->threshold;
+	G_message(_("Pass %d:"), t);
 
-	t++;
-	G_debug(3, "#######   Starting outer do loop! t = %d    #######", t);
-	if (verbose)
-	    G_verbose_message(_("Pass %d:"), t);
-	else
-	    G_percent(t, globals->end_t, 1);
-
 	n_merges = 0;
 	globals->candidate_count = 0;
 	flag_clear_all(globals->candidate_flag);
@@ -206,8 +239,7 @@
 
 	/*process candidate cells */
 	for (row = globals->row_min; row < globals->row_max; row++) {
-	    if (verbose)
-		G_percent(row, globals->row_max, 4);
+	    G_percent(row, globals->row_max, 4);
 	    for (col = globals->col_min; col < globals->col_max; col++) {
 		if (!(FLAG_GET(globals->candidate_flag, row, col)))
 		    continue;
@@ -235,23 +267,20 @@
 		/* find Ri's best neighbor, clear candidate flag */
 		Ri_similarity = globals->threshold + 1;
 
-		globals->rs.id = Ri.id;
-		if ((Ri_rs = rgtree_find(globals->reg_tree, &(globals->rs))) != NULL) {
-		    memcpy(Ri.mean, Ri_rs->mean, globals->datasize);
-		    Ri.count = Ri_rs->count;
-		}
-		else {
-		    segment_get(&globals->bands_seg, (void *)Ri.mean,
-				Ri.row, Ri.col);
-		    Ri.count = 1;
-		}
+		Ri_rs.id = Ri.id;
+		fetch_reg_stats(Ri.row, Ri.col, &Ri_rs, globals);
+		memcpy(Ri.mean, Ri_rs.mean, globals->datasize);
+		Ri.count = Ri_rs.count;
+
 		/* Ri is now complete */
+		G_debug(4, "Ri is now complete");
 
-		Rk_rs = NULL;
-		Ri_nn = find_best_neighbor(&Ri, Ri_ngbrs, &Rk_rs,
-					   &Rk, &Ri_similarity, 1,
-					   globals);
+		/* find best neighbor, get region stats for Rk */
+		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 (Ri_nn == 0) {
 		    /* this can only happen if only one segment is left */
@@ -259,7 +288,8 @@
 		    pathflag = FALSE;
 		    Ri.count = 0;
 		}
-		if (!(t & 1) && Ri_nn == 1 &&
+
+		if (/* !(t & 1) && */ Ri_nn == 1 &&
 		    !(FLAG_GET(globals->candidate_flag, Rk.row, Rk.col)) &&
 		    compare_double(Ri_similarity, threshold) == -1) {
 		    /* this is slow ??? */
@@ -276,21 +306,21 @@
 			if (Rk.count < 2)
 			    G_fatal_error("Rk count too low");
 			if (Rk.count < Ri.count)
-			    G_debug(1, "Rk count lower than Ri count");
+			    G_debug(4, "Rk count lower than Ri count");
 
-			merge_regions(&Ri, Ri_rs, &Rk, Rk_rs, globals);
+			merge_regions(&Ri, &Ri_rs, &Rk, &Rk_rs, 1, globals);
 			n_merges++;
 		    }
 
 		    pathflag = FALSE;
 		}
 		/* this is slow ??? */
-		if (t & 1) {
+		if (/* t & */ 1) {
 		    if ((globals->nn < 8 && Rk.count <= 8) || 
 		        (globals->nn >= 8 && Rk.count <= globals->nn))
 		    candidates_only = FALSE;
 		}
-		
+
 		while (pathflag) {
 		    pathflag = FALSE;
 		    
@@ -311,17 +341,23 @@
 			G_debug(4, "Working with Rk");
 
 			/* find Rk's best neighbor, do not clear candidate flag */
-			Rk_similarity = globals->threshold + 1;
-			Rk_bestn_rs = NULL;
-			Rk_nn = find_best_neighbor(&Rk, Rk_ngbrs,
+			/* Rk_similarity = globals->threshold + 1; */
+			Rk_similarity = Ri_similarity;
+			Rk_bestn_rs.count = 0;
+			/* Rk_rs is already complete */
+			Rk_nn = find_best_neighbor(&Rk, &Rk_rs, Rk_ngbrs,
+						   &Rk_bestn,
 						   &Rk_bestn_rs,
-						   &Rk_bestn,
 						   &Rk_similarity, 0,
 						   globals);
 
 			/* not mutually best neighbors */
 			if (Rk_similarity != Ri_similarity) {
-			    do_merge = 0;
+			    /* important for fp precision limit
+			     * because region stats may be calculated 
+			     * in two slightly different ways */
+			    if (fabs(Ri_similarity / Rk_similarity  - 1) > EPSILON)
+				do_merge = 0;
 			}
 			/* Ri has only one neighbor, merge */
 			if (Ri_nn == 1 && Rk_nn > 1)
@@ -337,8 +373,10 @@
 			    adjthresh = pow(alpha2, 1. + (double) smaller / divisor) *
 					globals->threshold;
 
-			    if (compare_double(Ri_similarity, adjthresh) == 1)
-				do_merge = 0;
+			    do_merge = 0;
+			    if (compare_double(Ri_similarity, adjthresh) == -1) {
+				do_merge = 1;
+			    }
 			}
 
 			if (do_merge) {
@@ -358,7 +396,7 @@
 			    nbtree_clear(Rk_ngbrs);
 			    Ri_nn += Ri_ngbrs->count;
 
-			    merge_regions(&Ri, Ri_rs, &Rk, Rk_rs, globals);
+			    merge_regions(&Ri, &Ri_rs, &Rk, &Rk_rs, 1, globals);
 			    /* Ri is now updated, Rk no longer usable */
 
 			    /* made a merge, need another iteration */
@@ -372,16 +410,11 @@
 			     */
 			    
 			    /* use neighbor tree to find Ri's new best neighbor */
-			    search_neighbors(&Ri, Ri_ngbrs, &Ri_similarity,
-					     &Rk, globals);
+			    search_neighbors(&Ri, &Ri_rs, Ri_ngbrs, &Ri_similarity,
+					     &Rk, &Rk_rs, globals);
 
 			    if (Ri_nn > 0 && compare_double(Ri_similarity, threshold) == -1) {
 
-				globals->rs.id = Ri.id;
-				Ri_rs = rgtree_find(globals->reg_tree, &(globals->rs));
-				globals->rs.id = Rk.id;
-				Rk_rs = rgtree_find(globals->reg_tree, &(globals->rs));
-
 				pathflag = TRUE;
 				/* candidates_only:
 				 * FALSE: less passes, takes a bit longer, but less memory
@@ -395,7 +428,7 @@
 			    if (compare_double(Rk_similarity, threshold) == -1) {
 				pathflag = TRUE;
 			    }
-			    /* test this */
+			    /* test this: can it cause an infinite loop ? */
 			    if (!(FLAG_GET(globals->candidate_flag, Rk.row, Rk.col)))
 				pathflag = FALSE;
 
@@ -422,11 +455,22 @@
 				Ri = Rk;
 				Rk = Rk_bestn;
 				Rk_bestn.mean = dp;
-				
-				Ri_rs = Rk_rs;
-				Rk_rs = Rk_bestn_rs;
-				Rk_bestn_rs = NULL;
 
+				Ri_rs.id = Rk_rs.id;
+				Rk_rs.id = Rk_bestn_rs.id;
+				Rk_bestn_rs.id = 0;
+				Ri_rs.count = Rk_rs.count;
+				Rk_rs.count = Rk_bestn_rs.count;
+				Rk_bestn_rs.count = 0;
+				dp = Ri_rs.mean;
+				Ri_rs.mean = Rk_rs.mean;
+				Rk_rs.mean = Rk_bestn_rs.mean;
+				Rk_bestn_rs.mean = dp;
+				dp = Ri_rs.sum;
+				Ri_rs.sum = Rk_rs.sum;
+				Rk_rs.sum = Rk_bestn_rs.sum;
+				Rk_bestn_rs.sum = dp;
+
 				tmpnbtree = Ri_ngbrs;
 				Ri_ngbrs = Rk_ngbrs;
 				Rk_ngbrs = tmpnbtree;
@@ -458,10 +502,11 @@
     if (globals->min_segment_size > 1) {
 	G_message(_("Merging segments smaller than %d cells"), globals->min_segment_size);
 
-	/* optional threshold as a function of t. */
 	threshold = globals->alpha * globals->alpha * globals->threshold;
 
 	flag_clear_all(globals->candidate_flag);
+	
+	n_merges = 0;
 
 	/* Set candidate flag to true/1 for all non-NULL cells */
 	for (row = globals->row_min; row < globals->row_max; row++) {
@@ -498,18 +543,17 @@
 		if (Ri.id < 0)		
 		    continue;
 
+		Ri_rs.id = Ri.id;
+
+		/* get segment size */
+
+		fetch_reg_stats(Ri.row, Ri.col, &Ri_rs, globals);
+		memcpy(Ri.mean, Ri_rs.mean, globals->datasize);
+		Ri.count = Ri_rs.count;
+
 		while (do_merge) {
 
-		    /* get segment size */
-		    globals->rs.id = Ri.id;
-		    Ri_rs = rgtree_find(globals->reg_tree, &(globals->rs));
-
-		    Ri.count = 1;
 		    do_merge = 0;
-		    if (Ri_rs != NULL) {
-			Ri.count = Ri_rs->count;
-			memcpy(Ri.mean, Ri_rs->mean, globals->datasize);
-		    }
 
 		    /* merge all smaller than min size */
 		    if (Ri.count < globals->min_segment_size)
@@ -524,8 +568,8 @@
 				    Ri.row, Ri.col);
 
 			/* find Ri's best neighbor, clear candidate flag */
-			Ri_nn = find_best_neighbor(&Ri, Ri_ngbrs,
-						   &Rk_rs, &Rk, 
+			Ri_nn = find_best_neighbor(&Ri, &Ri_rs, Ri_ngbrs,
+						   &Rk, &Rk_rs,
 						   &Ri_similarity, 1,
 						   globals);
 		    }
@@ -535,7 +579,9 @@
 			nbtree_clear(Ri_ngbrs);
 			
 			/* merge Ri with Rk */
-			merge_regions(&Ri, Ri_rs, &Rk, Rk_rs, globals);
+			/* do not clear candidate flag for Rk */
+			merge_regions(&Ri, &Ri_rs, &Rk, &Rk_rs, 0, globals);
+			n_merges++;
 
 			if (Ri_nn <= 0 || Ri.count >= globals->min_segment_size)
 			    do_merge = 0;
@@ -544,8 +590,12 @@
 	    }
 	}
 	G_percent(1, 1, 1);
+
+	/* 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);
     G_free(Rk_bestn.mean);
@@ -555,20 +605,29 @@
     free(Ri_ngbrs);
     free(Rk_ngbrs);
 
+    /* free region stats */
+    G_free(Ri_rs.mean);
+    G_free(Ri_rs.sum);
+    G_free(Rk_rs.mean);
+    G_free(Rk_rs.sum);
+    G_free(Rk_bestn_rs.mean);
+    G_free(Rk_bestn_rs.sum);
+
     return TRUE;
 }
 
 static int find_best_neighbor(struct ngbr_stats *Ri,
-		       struct NB_TREE *Ri_ngbrs, 
-		       struct reg_stats **Rk_rs,
-		       struct ngbr_stats *Rk, 
-		       double *sim, int clear_cand,
-		       struct globals *globals)
+			      struct reg_stats *Ri_rs,
+			      struct NB_TREE *Ri_ngbrs, 
+			      struct ngbr_stats *Rk, 
+			      struct reg_stats *Rk_rs,
+			      double *sim, int clear_cand,
+			      struct globals *globals)
 {
-    int n, n_ngbrs, no_check;
+    int n, n_ngbrs, no_check, cmp;
     struct rc ngbr_rc, next;
     struct rclist rilist;
-    double tempsim, *dp;
+    double tempsim;
     int neighbors[8][2];
     struct RB_TREE *no_check_tree;	/* cells already checked */
     struct reg_stats *rs_found;
@@ -586,18 +645,17 @@
     ngbr_rc.col = Ri->col;
     rbtree_insert(no_check_tree, &ngbr_rc);
 
-    Ri->count = 1;
-
     n_ngbrs = 0;
     /* TODO: add size of largest region to reg_tree, use this as min */
     Rk->count = globals->ncells;
 
     /* go through segment, spreading outwards from head */
     rclist_init(&rilist);
-    rclist_add(&rilist, Ri->row, Ri->col);
 
-    while (rclist_drop(&rilist, &next)) {
-
+    /* check neighbors of start cell */
+    next.row = Ri->row;
+    next.col = Ri->col;
+    do {
 	/* remove from candidates */
 	if (clear_cand)
 	    FLAG_UNSET(globals->candidate_flag, next.row, next.col);
@@ -633,7 +691,7 @@
 			/* not yet checked, don't check it again */
 			rbtree_insert(no_check_tree, &ngbr_rc);
 
-			/* get Rk ID */
+			/* get neighbor ID */
 			segment_get(&globals->rid_seg,
 				    (void *) &(globals->ns.id),
 				    ngbr_rc.row, ngbr_rc.col);
@@ -642,7 +700,6 @@
 
 			    /* want to check this neighbor's neighbors */
 			    rclist_add(&rilist, ngbr_rc.row, ngbr_rc.col);
-			    Ri->count++;
 			}
 			else {
 
@@ -653,44 +710,56 @@
 				globals->rs.id = globals->ns.id;
 				rs_found = rgtree_find(globals->reg_tree,
 						       &(globals->rs));
-				if (rs_found != NULL) {
-				    globals->ns.mean = rs_found->mean;
-				    globals->ns.count = rs_found->count;
+				if (!rs_found) {
+				    /* region stats are not in search tree */
+				    rs_found = &(globals->rs);
+				    calculate_reg_stats(ngbr_rc.row, ngbr_rc.col,
+							rs_found, globals);
 				}
-				else {
-				    segment_get(&globals->bands_seg,
-						(void *)globals->bands_val,
-						ngbr_rc.row, ngbr_rc.col);
+				globals->ns.mean = rs_found->mean;
+				globals->ns.count = rs_found->count;
+				/* globals->ns is now complete */
 
-				    globals->ns.mean = globals->bands_val;
-				    globals->ns.count = 1;
-				}
-				/* next Rk = globals->ns is now complete */
-
 				tempsim = (globals->calculate_similarity)(Ri, &globals->ns, globals);
 
-				if (compare_double(tempsim, *sim) == -1) {
+				cmp = compare_double(tempsim, *sim);
+				if (cmp == -1) {
 				    *sim = tempsim;
 				    /* copy temp Rk to Rk */
-				    dp = Rk->mean;
-				    *Rk = globals->ns;
-				    Rk->mean = dp;
-				    memcpy(Rk->mean, globals->ns.mean,
+				    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 = rs_found;
+
+				    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);
 				}
-				else if (compare_double(tempsim, *sim) == 0) {
+				else if (cmp == 0) {
 				    /* resolve ties: prefer smaller regions */
 
 				    if (Rk->count > globals->ns.count) {
 					/* copy temp Rk to Rk */
-					dp = Rk->mean;
-					*Rk = globals->ns;
-					Rk->mean = dp;
-					memcpy(Rk->mean,
-					       globals->ns.mean,
+					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 = rs_found;
+
+					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);
 				    }
 				}
 
@@ -702,7 +771,7 @@
 		}
 	    }
 	} while (n--);    /* end do loop - next neighbor */
-    }    /* while there are cells to check */
+    } while (rclist_drop(&rilist, &next));   /* while there are cells to check */
 
     /* clean up */
     rbtree_destroy(no_check_tree);
@@ -780,31 +849,36 @@
 
 
 static int search_neighbors(struct ngbr_stats *Ri,
+			    struct reg_stats *Ri_rs,
                             struct NB_TREE *Ri_ngbrs, 
 			    double *sim,
 			    struct ngbr_stats *Rk,
+			    struct reg_stats *Rk_rs,
 		            struct globals *globals)
 {
     double tempsim, *dp;
     struct NB_TRAV travngbr;
     struct ngbr_stats *next;
+    int cmp, i;
 
     G_debug(4, "search_neighbors");
 
     nbtree_init_trav(&travngbr, Ri_ngbrs);
-    Rk->count = globals->ncells;
+    Rk->count = globals->ncells + 1;
 
     while ((next = nbtree_traverse(&travngbr))) {
 	tempsim = (globals->calculate_similarity)(Ri, next, globals);
 
-	if (compare_double(tempsim, *sim) == -1) {
+	cmp = compare_double(tempsim, *sim);
+	if (cmp == -1) {
 	    *sim = tempsim;
+
 	    dp = Rk->mean;
 	    *Rk = *next;
 	    Rk->mean = dp;
 	    memcpy(Rk->mean, next->mean, globals->datasize);
 	}
-	else if (compare_double(tempsim, *sim) == 0) {
+	else if (cmp == 0) {
 	    /* resolve ties, prefer smaller regions */
 	    G_debug(4, "resolve ties");
 
@@ -816,150 +890,115 @@
 	    }
 	}
     }
+    Rk_rs->id = Rk->id;
 
+    /* faster, but with fp error:
+     * calculate sum from mean and count */
+    /*
+    Rk_rs->count = Rk->count;
+    memcpy(Rk_rs->mean, Rk->mean, globals->datasize);
+    i = globals->nbands - 1;
+    do {
+	Rk_rs->sum[i] = Rk_rs->mean[i] * Rk_rs->count;
+    } while (i--);
+    */
+
+    /* a bit slower but correct: */
+    fetch_reg_stats(Rk->row, Rk->col, Rk_rs, globals);
+
     return 1;
 }
 
 static int merge_regions(struct ngbr_stats *Ri, struct reg_stats *Ri_rs,
 		         struct ngbr_stats *Rk, struct reg_stats *Rk_rs,
-		         struct globals *globals)
+		         int do_cand, struct globals *globals)
 {
-    int n, do_cand;
-    int Ri_id, Rk_id, larger_id, R_id;
+    int n;
+    int R_id;
     struct rc next, ngbr_rc;
     struct rclist rlist;
     int neighbors[8][2];
-    struct reg_stats *new_rs, old_rs;
+    struct reg_stats *new_rs;
 
     G_debug(4, "merge_regions");
 
-    Ri_id = Ri->id;
-    Rk_id = Rk->id;
+    /* Ri ID must always be positive */
+    if (Ri_rs->id < 1)
+	G_fatal_error("Ri id is negative: %d", Ri_rs->id);
+    /* if Rk ID is negative (no seed), Rk count must be 1  */
+    if (Rk_rs->id < 1 && Rk_rs->count > 1)
+	G_fatal_error("Rk id is negative: %d, but count is > 1: %d",
+	              Rk_rs->id, Rk_rs->count);
 
-    /* update segment number and clear candidate flag */
+    /* update segment id and clear candidate flag */
     
     /* cases
-     * Ri, Rk are single cells, not in tree
+     * Ri, Rk are not in the tree
      * Ri, Rk are both in the tree
      * Ri is in the tree, Rk is not
      * Rk is in the tree, Ri is not
      */
 
-    /* for each member of Ri and Rk, write new average bands values and segment values */
+    /* Ri_rs, Rk_rs must always be set */
+    /* add Rk */
+    Ri_rs->count += Rk_rs->count;
+    n = globals->nbands - 1;
+    do {
+	Ri_rs->sum[n] += Rk_rs->sum[n];
+	Ri_rs->mean[n] = Ri_rs->sum[n] / Ri_rs->count;
+    } while (n--);
 
     if (Ri->count >= Rk->count) {
-	larger_id = Ri_id;
 
-	if (Ri_rs) {
-	    new_rs = Ri_rs;
-	}
-	else {
-	    new_rs = &(globals->rs);
-	    new_rs->id = Ri_id;
-	    new_rs->count = 1;
-	    
-	    memcpy(new_rs->sum, Ri->mean, globals->datasize);
-	    memcpy(new_rs->mean, Ri->mean, globals->datasize);
-	}
-
-	/* add Rk */
-	if (Rk_rs) {
-	    new_rs->count += Rk_rs->count;
-	    n = globals->nbands - 1;
-	    do {
-		new_rs->sum[n] += Rk_rs->sum[n];
-	    } while (n--);
-	}
-	else {
-	    new_rs->count++;
-	    n = globals->nbands - 1;
-	    do {
-		new_rs->sum[n] += Rk->mean[n];
-	    } while (n--);
-	}
-
-	n = globals->nbands - 1;
-	do {
-	    new_rs->mean[n] = new_rs->sum[n] / new_rs->count;
-	} while (n--);
-	
-	memcpy(Ri->mean, new_rs->mean, globals->datasize);
-	Ri->count = new_rs->count;
-
-	if (!Ri_rs) {
-	    /* add to tree */
-	    rgtree_insert(globals->reg_tree, new_rs);
-	}
-	if (Rk_rs) {
+	if (Rk->count >= globals->min_reg_size) {
+	    if (rgtree_find(globals->reg_tree, Rk_rs) == NULL)
+		G_fatal_error("merge regions: Rk should be in tree");
 	    /* remove from tree */
-	    old_rs.id = Rk_id;
-	    rgtree_remove(globals->reg_tree, &old_rs);
+	    rgtree_remove(globals->reg_tree, Rk_rs);
 	}
     }
     else {
-	larger_id = Rk_id;
-	Ri->id = Rk_id;
 
-	if (Rk_rs) {
-	    new_rs = Rk_rs;
+	if (Ri->count >= globals->min_reg_size) {
+	    if (rgtree_find(globals->reg_tree, Ri_rs) == NULL)
+		G_fatal_error("merge regions: Ri should be in tree");
+	    /* remove from tree */
+	    rgtree_remove(globals->reg_tree, Ri_rs);
 	}
-	else {
-	    new_rs = &(globals->rs);
-	    new_rs->id = Rk_id;
-	    new_rs->count = 1;
 
-	    memcpy(new_rs->sum, Rk->mean, globals->datasize);
-	    memcpy(new_rs->mean, Rk->mean, globals->datasize);
-	}
+	/* magic switch */
+	Ri_rs->id = Rk->id;
+    }
 
-	/* add Ri */
-	if (Ri_rs) {
-	    new_rs->count += Ri_rs->count;
-	    n = globals->nbands - 1;
-	    do {
-		new_rs->sum[n] += Ri_rs->sum[n];
-	    } while (n--);
-	}
-	else {
-	    new_rs->count++;
-	    n = globals->nbands - 1;
-	    do {
-		new_rs->sum[n] += Ri->mean[n];
-	    } while (n--);
-	}
+    if ((new_rs = rgtree_find(globals->reg_tree, Ri_rs)) != NULL) {
+	/* update stats for tree item */
+	new_rs->count = Ri_rs->count;
+	memcpy(new_rs->mean, Ri_rs->mean, globals->datasize);
+	memcpy(new_rs->sum, Ri_rs->sum, globals->datasize);
+    }
+    else if (Ri_rs->count >= globals->min_reg_size) {
+	/* add to tree */
+	rgtree_insert(globals->reg_tree, Ri_rs);
+    }
 
-	n = globals->nbands - 1;
-	do {
-	    new_rs->mean[n] = new_rs->sum[n] / new_rs->count;
-	} while (n--);
+    Ri->count = Ri_rs->count;
+    memcpy(Ri->mean, Ri_rs->mean, globals->datasize);
 
-	memcpy(Ri->mean, new_rs->mean, globals->datasize);
-	Ri->count = new_rs->count;
-
-	if (!Rk_rs) {
-	    /* add to tree */
-	    rgtree_insert(globals->reg_tree, new_rs);
-	}
-	if (Ri_rs) {
-	    /* remove from tree */
-	    old_rs.id = Ri_id;
-	    rgtree_remove(globals->reg_tree, &old_rs);
-	}
-    }
-
-    if (larger_id == Ri_id) {
+    if (Ri->id == Ri_rs->id) {
 	/* Ri is already updated, including candidate flags
 	 * need to clear candidate flag for Rk and set new id */
 	 
 	/* the actual merge: change region id */
-	segment_put(&globals->rid_seg, (void *) &Ri_id, Rk->row, Rk->col);
+	segment_put(&globals->rid_seg, (void *) &Ri->id, Rk->row, Rk->col);
 
-	do_cand = 0;
-	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--;
-	    do_cand = 1;
+	if (do_cand) {
+	    do_cand = 0;
+	    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--;
+		do_cand = 1;
+	    }
 	}
 
 	rclist_init(&rlist);
@@ -988,11 +1027,12 @@
 
 		    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);
+			segment_get(&globals->rid_seg, (void *) &R_id,
+			            ngbr_rc.row, ngbr_rc.col);
 
-			if (R_id == Rk_id) {
+			if (R_id == Rk->id) {
 			    /* the actual merge: change region id */
-			    segment_put(&globals->rid_seg, (void *) &Ri_id, ngbr_rc.row, ngbr_rc.col);
+			    segment_put(&globals->rid_seg, (void *) &Ri->id, ngbr_rc.row, ngbr_rc.col);
 
 			    /* want to check this neighbor's neighbors */
 			    rclist_add(&rlist, ngbr_rc.row, ngbr_rc.col);
@@ -1003,17 +1043,17 @@
 	}
     }
     else {
-	/* larger_id == Rk_id */
+	/* Rk was larger than Ri */
 
 	/* clear candidate flag for Rk */
-	if (FLAG_GET(globals->candidate_flag, Rk->row, Rk->col)) {
+	if (do_cand && FLAG_GET(globals->candidate_flag, Rk->row, Rk->col)) {
 	    set_candidate_flag(Rk, FALSE, globals);
 	}
 
 	/* update region id for Ri */
 
 	/* the actual merge: change region id */
-	segment_put(&globals->rid_seg, (void *) &Rk_id, Ri->row, Ri->col);
+	segment_put(&globals->rid_seg, (void *) &Rk->id, Ri->row, Ri->col);
 
 	rclist_init(&rlist);
 	rclist_add(&rlist, Ri->row, Ri->col);
@@ -1038,9 +1078,9 @@
 
 			segment_get(&globals->rid_seg, (void *) &R_id, ngbr_rc.row, ngbr_rc.col);
 
-			if (R_id == Ri_id) {
+			if (R_id == Ri->id) {
 			    /* the actual merge: change region id */
-			    segment_put(&globals->rid_seg, (void *) &Rk_id, ngbr_rc.row, ngbr_rc.col);
+			    segment_put(&globals->rid_seg, (void *) &Rk->id, ngbr_rc.row, ngbr_rc.col);
 
 			    /* want to check this neighbor's neighbors */
 			    rclist_add(&rlist, ngbr_rc.row, ngbr_rc.col);
@@ -1049,11 +1089,18 @@
 		}
 	    } while (n--);
 	}
+	Ri->id = Ri_rs->id;   /* == Rk->id */
+	if (Ri->id != Rk->id)
+	    G_fatal_error("Ri ID should be set to Rk ID");
     }
     
-    if (Rk_id > 0)
+    if (Rk->id > 0)
 	globals->n_regions--;
 
+    /* disable Rk */
+    Rk->id = Rk_rs->id = 0;
+    Rk->count = Rk_rs->count = 0;
+
     return TRUE;
 }
 
@@ -1127,3 +1174,168 @@
 
     return TRUE;
 }
+
+int fetch_reg_stats(int row, int col, struct reg_stats *rs, 
+                           struct globals *globals)
+{
+    struct reg_stats *rs_found;
+
+    if ((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);
+	rs->count = rs_found->count;
+
+	return 1;
+    }
+
+    calculate_reg_stats(row, col, rs, globals);
+
+    return 2;
+}
+
+static int calculate_reg_stats(int row, int col, struct reg_stats *rs, 
+                         struct globals *globals)
+{
+    G_debug(4, "calculate_reg_stats()");
+
+    segment_get(&globals->bands_seg, (void *)globals->bands_val,
+		row, col);
+    rs->count = 1;
+    memcpy(rs->mean, globals->bands_val, globals->datasize);
+    memcpy(rs->sum, globals->bands_val, globals->datasize);
+
+    if (globals->min_reg_size < 3)
+	return 1;
+
+    if (globals->min_reg_size == 3) {
+	int n, i, rid;
+	struct rc ngbr_rc;
+	int neighbors[8][2];
+
+	globals->find_neighbors(row, col, neighbors);
+
+	for (n = 0; n < globals->nn; n++) {
+
+	    ngbr_rc.row = neighbors[n][0];
+	    ngbr_rc.col = neighbors[n][1];
+
+	    if (ngbr_rc.row < 0 || ngbr_rc.row >= globals->nrows ||
+		ngbr_rc.col < 0 || ngbr_rc.col >= globals->ncols) {
+		continue;
+	    }
+
+	    if ((FLAG_GET(globals->null_flag, ngbr_rc.row, ngbr_rc.col)) == 0) {
+
+		segment_get(&globals->rid_seg, (void *) &rid,
+			    ngbr_rc.row, ngbr_rc.col);
+		
+		if (rid == rs->id) {
+
+		    /* update region stats */
+		    segment_get(&globals->bands_seg, (void *)globals->bands_val,
+				ngbr_rc.row, ngbr_rc.col);
+
+		    i = globals->nbands - 1;
+		    do {
+			rs->sum[i] += globals->bands_val[i];
+		    } while (i--);
+		    rs->count++;
+
+		    /* only one other neighbor can have the same ID */
+		    break;
+		}
+	    }
+	}
+
+	/* band mean */
+	i = globals->nbands - 1;
+	do {
+	    rs->mean[i] = rs->sum[i] / rs->count;
+	} while (i--);
+	
+	return 2;
+    }
+
+    if (globals->min_reg_size > 3) {
+	/* rs->id must be set */
+	struct RB_TREE *rc_check_tree;	/* cells already checked */
+	int n, i, rid;
+	struct rc ngbr_rc, next;
+	struct rclist rilist;
+	int neighbors[8][2];
+	int no_check;
+	
+	/* go through region, spreading outwards from head */
+	rclist_init(&rilist);
+
+	rc_check_tree = rbtree_create(compare_rc, sizeof(struct rc));
+	ngbr_rc.row = row;
+	ngbr_rc.col = col;
+	rbtree_insert(rc_check_tree, &ngbr_rc);
+
+	next.row = row;
+	next.col = col;
+	do {
+	    G_debug(5, "find_pixel_neighbors for row: %d , col %d",
+		    next.row, next.col);
+
+	    globals->find_neighbors(next.row, next.col, neighbors);
+
+	    n = globals->nn - 1;
+	    n = 0;
+	    do {
+
+		ngbr_rc.row = neighbors[n][0];
+		ngbr_rc.col = neighbors[n][1];
+
+		no_check = (ngbr_rc.row < 0 || ngbr_rc.row >= globals->nrows ||
+		    ngbr_rc.col < 0 || ngbr_rc.col >= globals->ncols);
+
+		if (!no_check) {
+		    if ((FLAG_GET(globals->null_flag, ngbr_rc.row, ngbr_rc.col)) == 0) {
+		    
+			/* already checked ? */
+			if (!rbtree_find(rc_check_tree, &ngbr_rc)) {
+
+			    /* not yet checked, don't check it again */
+			    rbtree_insert(rc_check_tree, &ngbr_rc);
+
+			    segment_get(&globals->rid_seg, (void *) &rid,
+					ngbr_rc.row, ngbr_rc.col);
+			    
+			    if (rid == rs->id) {
+
+				/* want to check this neighbor's neighbors */
+				rclist_add(&rilist, ngbr_rc.row, ngbr_rc.col);
+
+				/* update region stats */
+				segment_get(&globals->bands_seg,
+					    (void *)globals->bands_val,
+					    ngbr_rc.row, ngbr_rc.col);
+
+				i = globals->nbands - 1;
+				do {
+				    rs->sum[i] += globals->bands_val[i];
+				} while (i--);
+				rs->count++;
+			    }
+			}
+		    }
+		}
+	    } while (n++ < globals->nn - 1); /* (n--); */
+	} while (rclist_drop(&rilist, &next));
+	/* band mean */
+	i = globals->nbands - 1;
+	do {
+	    rs->mean[i] = rs->sum[i] / rs->count;
+	} while (i--);
+
+	/* clean up */
+	rbtree_destroy(rc_check_tree);
+	
+	return 3;
+    }
+    
+    return 0;
+}

Modified: grass-addons/grass7/imagery/i.segment.xl/iseg.h
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/iseg.h	2012-09-27 11:33:21 UTC (rev 53275)
+++ grass-addons/grass7/imagery/i.segment.xl/iseg.h	2012-09-27 15:44:17 UTC (rev 53276)
@@ -67,7 +67,8 @@
 
     /* results */
     struct RG_TREE *reg_tree;   /* search tree with region stats */
-    struct reg_stats rs;
+    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;
@@ -100,9 +101,15 @@
 void find_four_neighbors(int, int, int[][2]);
 void find_eight_neighbors(int, int, int[8][2]);
 double calculate_euclidean_similarity(struct ngbr_stats *, 
-                                      struct ngbr_stats *, struct globals *);
+                                      struct ngbr_stats *,
+				      struct globals *);
+int fetch_reg_stats(int , int , struct reg_stats *, 
+                           struct globals *);
 
+/* void calculate_reg_stats(int, int, struct reg_stats *, 
+                         struct globals *); */
 
+
 /* rclist.c */
 void rclist_init(struct rclist *);
 void rclist_add(struct rclist *, int, int);

Modified: grass-addons/grass7/imagery/i.segment.xl/open_files.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/open_files.c	2012-09-27 11:33:21 UTC (rev 53275)
+++ grass-addons/grass7/imagery/i.segment.xl/open_files.c	2012-09-27 15:44:17 UTC (rev 53276)
@@ -9,6 +9,7 @@
 
 static int load_seeds(struct globals *, int, int, int);
 static int read_seed(struct globals *, SEGMENT *, struct rc *, int);
+static int manage_memory(int, int, struct globals *);
 
 int open_files(struct globals *globals)
 {
@@ -87,23 +88,16 @@
 
     inlen = sizeof(DCELL) * Ref.nfiles;
     outlen = sizeof(CELL);
+    G_debug(1, "data element size, in: %d , out: %d ", inlen, outlen);
+    globals->datasize = sizeof(double) * globals->nbands;
 
     /* segment lib segment size */
     srows = 64;
     scols = 64;
 
-    /* calculate number of segments in memory */
-    nseg = (1024. * 1024. * globals->mb) /
-           (sizeof(DCELL) * Ref.nfiles * srows * scols + 
-	    + sizeof(CELL) * 2 * srows * scols);
+    nseg = manage_memory(srows, scols, globals);
 
     /* create segment structures */
-    G_debug(1, "image size:  %d rows, %d cols", globals->nrows, globals->ncols);
-    G_debug(1, "degmented to tiles with size:  %d rows, %d cols", srows,
-	    scols);
-    G_debug(1, "data element size, in: %d , out: %d ", inlen, outlen);
-    G_debug(1, "number of segments in memory: %d", nseg);
-
     if (segment_open
 	(&globals->bands_seg, G_tempfile(), globals->nrows, globals->ncols, srows,
 	 scols, inlen, nseg) != 1)
@@ -259,7 +253,6 @@
 	Rast_close(in_fd[n]);
     }
 
-    globals->datasize = sizeof(double) * globals->nbands;
     globals->rs.sum = G_malloc(globals->datasize);
     globals->rs.mean = G_malloc(globals->datasize);
 
@@ -291,6 +284,8 @@
     struct rc Ri;
 
     G_debug(1, "load_seeds()");
+    
+    G_message(_("Loading seeds from '%s'"), globals->seeds);
 
     if (segment_open
 	(&seeds_seg, G_tempfile(), globals->nrows, globals->ncols,
@@ -353,7 +348,6 @@
 	}
     }
 
-
     G_free(seeds_buf);
     Rast_close(seeds_fd);
     segment_close(&seeds_seg);
@@ -457,7 +451,7 @@
     }
     
     /* insert into region tree */
-    if (globals->rs.count > 1) {
+    if (globals->rs.count >= globals->min_reg_size) {
 	for (i = 0; i < globals->nbands; i++)
 	    globals->rs.mean[i] = globals->rs.sum[i] / globals->rs.count;
 
@@ -466,3 +460,64 @@
 
     return 1;
 }
+
+
+static int manage_memory(int srows, int scols, struct globals *globals)
+{
+    double reg_size_mb, segs_mb;
+    int reg_size_count, nseg, nseg_total;
+
+    /* 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 = 4;
+    if (reg_size_count < (double) globals->nrows * globals->ncols / globals->min_reg_size) {
+	globals->min_reg_size = (double) globals->nrows * globals->ncols / reg_size_count;
+    }
+    else {
+	reg_size_count = (double) globals->nrows * globals->ncols / globals->min_reg_size;
+	/* recalculate segs_mb */
+	segs_mb = globals->mb - reg_size_count * reg_size_mb;
+    }
+
+    G_verbose_message(_("Stats for 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);
+    }
+    else {
+	/* input bands, segment ids, bounds map */
+	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));
+
+    if (nseg > nseg_total)
+	nseg = nseg_total;
+    
+    G_debug(1, "current region:  %d rows, %d cols", globals->nrows, globals->ncols);
+    G_debug(1, "segmented to tiles with size:  %d rows, %d cols", srows,
+	    scols);
+    G_verbose_message(_("Number of segments in memory: %d of %d total"),
+                      nseg, nseg_total);
+    
+    return nseg;
+}

Modified: grass-addons/grass7/imagery/i.segment.xl/rclist.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/rclist.c	2012-09-27 11:33:21 UTC (rev 53275)
+++ grass-addons/grass7/imagery/i.segment.xl/rclist.c	2012-09-27 15:44:17 UTC (rev 53276)
@@ -48,8 +48,8 @@
 
 	return 1;
     }
-    else
-	return 0;
+
+    return 0;
 }
 
 void rclist_destroy(struct rclist *list)

Modified: grass-addons/grass7/imagery/i.segment.xl/write_output.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/write_output.c	2012-09-27 11:33:21 UTC (rev 53275)
+++ grass-addons/grass7/imagery/i.segment.xl/write_output.c	2012-09-27 15:44:17 UTC (rev 53276)
@@ -14,51 +14,25 @@
 
 int write_output(struct globals *globals)
 {
-    int out_fd, mean_fd, row, col;
+    int out_fd, row, col;
     CELL *outbuf, rid;
-    FCELL *meanbuf;
     struct Colors colors;
-    double thresh, maxdev, sim, mingood;
-    struct reg_stats *rs_found;
-    struct ngbr_stats Ri, Rk;
+    struct History hist;
 
     outbuf = Rast_allocate_c_buf();
 
     G_debug(1, "preparing output raster");
     /* open output raster map */
     out_fd = Rast_open_new(globals->out_name, CELL_TYPE);
-    if (globals->out_band) {
-	mean_fd = Rast_open_new(globals->out_band, FCELL_TYPE);
-	meanbuf = Rast_allocate_f_buf();
-    }
-    else {
-	mean_fd = -1;
-	meanbuf = NULL;
-    }
 
-    /* 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->threshold;
-    maxdev = globals->threshold * (1 - globals->alpha * globals->alpha);
-    mingood = 1;
-
     G_debug(1, "start data transfer from segmentation file to raster");
 
-    G_message(_("Writing output"));
+    G_message(_("Writing out segment IDs"));
     for (row = 0; row < globals->nrows; row++) {
 
 	G_percent(row, globals->nrows, 9);
 
 	Rast_set_c_null_value(outbuf, globals->ncols);
-	if (globals->out_band)
-	    Rast_set_f_null_value(meanbuf, globals->ncols);
 	for (col = 0; col < globals->ncols; col++) {
 
 	    if (!(FLAG_GET(globals->null_flag, row, col))) {
@@ -66,15 +40,70 @@
 
 		if (rid > 0) {
 		    outbuf[col] = rid;
+		}
+	    }
+	}
+	Rast_put_row(out_fd, outbuf, CELL_TYPE);
+    }
 
-		    if (globals->out_band) {
+    /* close and save segment id file */
+    Rast_close(out_fd);
 
-			/* get values for Ri = larger region */
+    /* set colors */
+    Rast_init_colors(&colors);
+    Rast_make_random_colors(&colors, 1, globals->nrows * globals->ncols);
+    Rast_write_colors(globals->out_name, G_mapset(), &colors);
+
+    Rast_short_history(globals->out_name, "raster", &hist);
+    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;
+
+	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->threshold;
+	maxdev = globals->threshold * (1 - globals->alpha * globals->alpha);
+	mingood = 1;
+
+	G_message(_("Writing out goodness of fit"));
+	for (row = 0; row < globals->nrows; row++) {
+
+	    G_percent(row, globals->nrows, 9);
+
+	    Rast_set_f_null_value(meanbuf, globals->ncols);
+
+	    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) {
+
+			/* get values for Ri = this region */
 			globals->rs.id = rid;
-			rs_found = rgtree_find(globals->reg_tree, &(globals->rs));
+			fetch_reg_stats(Ri.row, Ri.col, &globals->rs, globals);
+			Ri.mean = globals->rs.mean;
+			Ri.count = globals->rs.count; 
 
-			if (rs_found != NULL) {
-			    Ri.mean = rs_found->mean;
+			sim = 0.;
+			/* region consists of more than one cell */
+			if (Ri.count > 1) {
 
 			    /* get values for Rk = this cell */
 			    segment_get(&globals->bands_seg,
@@ -85,9 +114,6 @@
 			    /* calculate similarity */
 			    sim = (*globals->calculate_similarity) (&Ri, &Rk, globals);
 			}
-			else
-			    /* region consists of only one cell */
-			    sim = 0.;
 			
 			if (0) {
 			    if (sim < thresh)
@@ -108,31 +134,24 @@
 		    }
 		}
 	    }
+	    Rast_put_row(mean_fd, meanbuf, FCELL_TYPE);
 	}
-	Rast_put_row(out_fd, outbuf, CELL_TYPE);
-	if (globals->out_band)
-	    Rast_put_row(mean_fd, meanbuf, FCELL_TYPE);
-    }
 
-    /* close and save file */
-    Rast_close(out_fd);
-    if (globals->out_band)
 	Rast_close(mean_fd);
 
-    /* set colors */
-    Rast_init_colors(&colors);
-    Rast_make_random_colors(&colors, 1, globals->nrows * globals->ncols);
-    Rast_write_colors(globals->out_name, G_mapset(), &colors);
-    
-    if (globals->out_band) {
 	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);
+
+	G_free(meanbuf);
     }
 
     /* free memory */
     G_free(outbuf);
-    G_free(meanbuf);
     Rast_free_colors(&colors);
 
     return TRUE;
@@ -151,12 +170,6 @@
 
     segment_close(&globals->rid_seg);
 
-    /*
-    for (i = 0; i < globals->nrows; i++)
-	G_free(globals->rid[i]);
-    G_free(globals->rid);
-    */
-
     flag_destroy(globals->null_flag);
     flag_destroy(globals->candidate_flag);
 



More information about the grass-commit mailing list