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

svn_grass at osgeo.org svn_grass at osgeo.org
Fri Oct 26 07:03:34 PDT 2012


Author: momsen
Date: 2012-10-26 07:03:33 -0700 (Fri, 26 Oct 2012)
New Revision: 53561

Added:
   grass-addons/grass7/imagery/i.segment/estimate_threshold.c
Modified:
   grass-addons/grass7/imagery/i.segment/create_isegs.c
   grass-addons/grass7/imagery/i.segment/i.segment.html
   grass-addons/grass7/imagery/i.segment/iseg.h
   grass-addons/grass7/imagery/i.segment/main.c
   grass-addons/grass7/imagery/i.segment/open_files.c
   grass-addons/grass7/imagery/i.segment/parse_args.c
Log:
added shape parameters and threshold estimating option

Modified: grass-addons/grass7/imagery/i.segment/create_isegs.c
===================================================================
--- grass-addons/grass7/imagery/i.segment/create_isegs.c	2012-10-26 13:37:15 UTC (rev 53560)
+++ grass-addons/grass7/imagery/i.segment/create_isegs.c	2012-10-26 14:03:33 UTC (rev 53561)
@@ -4,7 +4,7 @@
 
 #include <stdlib.h>
 #include <float.h>		/* for DBL_MAX */
-#include <math.h>		/* for fabs() */
+#include <math.h>		/* for fabs() and sqrt() */
 #include <grass/gis.h>
 #include <grass/glocale.h>
 #include <grass/raster.h>
@@ -18,10 +18,19 @@
 #endif
 
 /* This will do a typical rowmajor processing of the image(s).  
- * Commenting it out will switch to z-order processing.
- * Z-order does NOT support the bounds option, and assumes a somewhat square rectangle (otherwise the power of 2 square could overrun the long long maximum). */
-#define ROWMAJOR		/* note: gives nesting error with INDENT program since this changes the for loops.  comment out when running INDENT. */
+ * Z-order was implemented in Revision 53236, but since removed.
+ * It did not show a speed increase on small rasters.  It could
+ * be considered again after other processing steps are sped up. */
 
+/* is there a better way to do this? */
+#ifndef max
+#define max(a,b) ( ((a)>(b)) ? (a) : (b) )
+#endif
+#ifndef min
+#define min(a,b) ( ((a)<(b)) ? (a) : (b) )
+#endif
+
+
 int create_isegs(struct files *files, struct functions *functions)
 {
     int lower_bound, upper_bound, row, col;
@@ -159,11 +168,6 @@
 	*current, *newpixel, *Ri_bestn;
     int Ri_count, Rk_count;	/* number of pixels/cells in Ri and Rk */
 
-#ifndef ROWMAJOR
-    unsigned long int z, end_z;	/* only used for z-order */
-    int j;
-#endif
-
 #ifdef VCLOSE
     struct pixels *Rclose_head, *Rc_head, *Rc_tail, *Rcn_head;
     int Rc_count;
@@ -212,6 +216,7 @@
 
     /* do while loop until no merges are made, or until t reaches maximum number of iterations */
     do {
+
 #ifdef PROFILE
 	pass_start = clock();
 #endif
@@ -230,301 +235,260 @@
 	/*process candidate pixels for this iteration */
 
 	/*check each pixel, start the processing only if it is a candidate pixel */
-#ifdef ROWMAJOR
+
 	for (row = files->minrow; row < files->maxrow; row++) {
 	    for (col = files->mincol; col < files->maxcol; col++) {
-#else
-	{
-	    /* z-order traversal: */
-	    /* need to get a "square" power of 2 around our processing area */
 
-	    /*largest dimension: */
-	    if (files->nrows > files->ncols)
-		end_z = files->nrows;
-	    else
-		end_z = files->ncols;
+		if (FLAG_GET(files->candidate_flag, row, col) &&
+		    FLAG_GET(files->seeds_flag, row, col)) {
 
-	    /* largest power of 2: */
-	    end_z--;		/* in case we are already a power of two. */
-	    end_z = (end_z >> 1) | end_z;
-	    end_z = (end_z >> 2) | end_z;
-	    end_z = (end_z >> 4) | end_z;
-	    end_z = (end_z >> 8) | end_z;
-	    end_z = (end_z >> 16) | end_z;
-	    end_z = (end_z >> 32) | end_z;	/* only for 64-bit architecture TODO, would this mess things up on 32? */
-	    /*TODO does this need to repeat more since z long unsigned??? */
-	    end_z++;
-
-	    /*squared: */
-	    end_z *= end_z;
-
-	    for (z = 0; z < end_z; z++) {
-		row = col = 0;
-		/*bit wise construct row and col from i */
-		for (j = 8 * sizeof(long int) - 1; j > 1; j--) {
-		    row = row | (1 & (z >> j));
-		    row = row << 1;
-		    j--;
-		    col = col | (1 & (z >> j));
-		    col = col << 1;
-		}
-		row = row | (1 & (z >> j));
-		j--;
-		col = col | (1 & (z >> j));
-		if (row >= files->nrows || col >= files->ncols)
-		    continue;	/* slight speed enhancement, if z-order is helpful, figure out how to skip to the next z that is in the processing area. */
-	    }
-#endif
-
-	    if (FLAG_GET(files->candidate_flag, row, col) &&
-		FLAG_GET(files->seeds_flag, row, col)) {
-
-		/*free memory for linked lists */
-		my_dispose_list(files->token, &Ri_head);
-		my_dispose_list(files->token, &Rk_head);
-		my_dispose_list(files->token, &Rin_head);
-		my_dispose_list(files->token, &Rkn_head);
+		    /*free memory for linked lists */
+		    my_dispose_list(files->token, &Ri_head);
+		    my_dispose_list(files->token, &Rk_head);
+		    my_dispose_list(files->token, &Rin_head);
+		    my_dispose_list(files->token, &Rkn_head);
 #ifdef VCLOSE
-		my_dispose_list(files->token, &Rclose_head);
+		    my_dispose_list(files->token, &Rclose_head);
 #endif
-		Rk_count = 0;
+		    Rk_count = 0;
 
-		/* First pixel in Ri is current row/col pixel.  We may add more later if it is part of a segment */
-		Ri_count = 1;
-		newpixel = (struct pixels *)link_new(files->token);
-		newpixel->next = NULL;
-		newpixel->row = row;
-		newpixel->col = col;
-		Ri_head = newpixel;
+		    /* First pixel in Ri is current row/col pixel.  We may add more later if it is part of a segment */
+		    Ri_count = 1;
+		    newpixel = (struct pixels *)link_new(files->token);
+		    newpixel->next = NULL;
+		    newpixel->row = row;
+		    newpixel->col = col;
+		    Ri_head = newpixel;
 
-		pathflag = TRUE;
+		    pathflag = TRUE;
 
-		while (pathflag == TRUE) {	/*if don't find mutual neighbors on first try, will use Rk as next Ri. */
-		    G_debug(4, "Next starting pixel: row, %d, col, %d",
-			    Ri_head->row, Ri_head->col);
+		    while (pathflag == TRUE) {	/*if don't find mutual neighbors on first try, will use Rk as next Ri. */
+			G_debug(4, "Next starting pixel: row, %d, col, %d",
+				Ri_head->row, Ri_head->col);
 
-		    pathflag = FALSE;
+			pathflag = FALSE;
 
-		    /* find segment neighbors, if we don't already have them */
-		    if (Rin_head == NULL) {
+			/* find segment neighbors, if we don't already have them */
+			if (Rin_head == NULL) {
 #ifdef PROFILE
-			fn_start = clock();
+			    fn_start = clock();
 #endif
 
-			find_segment_neighbors
-			    (&Ri_head, &Rin_head, &Ri_count, files,
-			     functions);
+			    find_segment_neighbors
+				(&Ri_head, &Rin_head, &Ri_count, files,
+				 functions);
 #ifdef PROFILE
-			fn_end = clock();
-			fn_lap =
-			    ((double)(fn_end - fn_start)) / CLOCKS_PER_SEC;
-			fn_accum += fn_lap;
-			fprintf(stdout, "fsn(Ri): %g\t", fn_lap);
+			    fn_end = clock();
+			    fn_lap =
+				((double)(fn_end - fn_start)) /
+				CLOCKS_PER_SEC;
+			    fn_accum += fn_lap;
+			    fprintf(stdout, "fsn(Ri): %g\t", fn_lap);
 #endif
-		    }
+			}
 
-		    if (Rin_head != NULL) {	/*found neighbors, find best neighbor then see if is mutually best neighbor */
+			if (Rin_head != NULL) {	/*found neighbors, find best neighbor then see if is mutually best neighbor */
 
-			/* ********  find Ri's most similar neighbor  ******** */
-			Ri_bestn = NULL;
-			Ri_similarity = threshold + 1;	/* set current similarity to max value */
-			segment_get(&files->bands_seg, (void *)files->bands_val, Ri_head->row, Ri_head->col);	/* current segment values */
+			    /* ********  find Ri's most similar neighbor  ******** */
+			    Ri_bestn = NULL;
+			    Ri_similarity = threshold + 1;	/* set current similarity to max value */
+			    segment_get(&files->bands_seg, (void *)files->bands_val, Ri_head->row, Ri_head->col);	/* current segment values */
 
-			/* for each of Ri's neighbors */
-			for (current = Rin_head; current != NULL;
-			     current = current->next) {
-			    tempsim = (*functions->calculate_similarity)
-				(Ri_head, current, files, functions);
+			    /* for each of Ri's neighbors */
+			    for (current = Rin_head; current != NULL;
+				 current = current->next) {
+				tempsim = (*functions->calculate_similarity)
+				    (Ri_head, current, files, functions);
 
 #ifdef VCLOSE
-			    /* if very close, will merge, but continue checking other neighbors */
-			    if (tempsim < functions->very_close * threshold) {
-				/* add to Rclose list */
+				/* if very close, will merge, but continue checking other neighbors */
+				if (tempsim <
+				    functions->very_close * threshold) {
+				    /* add to Rclose list */
+				    newpixel = (struct pixels *)
+					link_new(files->token);
+				    newpixel->next = Rclose_head;
+				    newpixel->row = current->row;
+				    newpixel->col = current->col;
+				    Rclose_head = newpixel;
+				}
+				/* If "sort of" close, merge only if it is the mutually most similar */
+				else
+#endif
+				if (tempsim < Ri_similarity) {
+				    Ri_similarity = tempsim;
+				    Ri_bestn = current;
+				}
+			    }	/* finished similiarity check for all neighbors */
+
+			    /* *** merge all the "very close" pixels/segments *** */
+			    /* doing this after checking all Rin, so we don't change the bands_val between similarity comparisons
+			     * ... but that leaves the possibility that we have the wrong best Neighbor after doing these merges... 
+			     * but it seems we can't put this merge after the Rk/Rkn portion of the loop, because we are changing the available neighbors
+			     * ...maybe this extra "very close" idea has to be done completely differently or dropped???  */
+#ifdef VCLOSE
+			    for (current = Rclose_head; current != NULL;
+				 current = current->next) {
+				my_dispose_list(files->token, &Rc_head);
+				my_dispose_list(files->token, &Rcn_head);
+
+				/* get membership of neighbor segment */
+				Rc_count = 1;
 				newpixel =
 				    (struct pixels *)link_new(files->token);
-				newpixel->next = Rclose_head;
+				newpixel->next = NULL;
 				newpixel->row = current->row;
 				newpixel->col = current->col;
-				Rclose_head = newpixel;
-			    }
-			    /* If "sort of" close, merge only if it is the mutually most similar */
-			    else
-#endif
-			    if (tempsim < Ri_similarity) {
-				Ri_similarity = tempsim;
-				Ri_bestn = current;
-			    }
-			}	/* finished similiarity check for all neighbors */
+				Rc_head = Rc_tail = newpixel;
+				find_segment_neighbors(&Rc_head, &Rcn_head, &Rc_count, files, functions);	/* just to get members, not looking at neighbors now */
+				merge_values(Ri_head, Rc_head, Ri_count,
+					     Rc_count, files);
 
-			/* *** merge all the "very close" pixels/segments *** */
-			/* doing this after checking all Rin, so we don't change the bands_val between similarity comparisons
-			 * ... but that leaves the possibility that we have the wrong best Neighbor after doing these merges... 
-			 * but it seems we can't put this merge after the Rk/Rkn portion of the loop, because we are changing the available neighbors
-			 * ...maybe this extra "very close" idea has to be done completely differently or dropped???  */
-#ifdef VCLOSE
-			for (current = Rclose_head; current != NULL;
-			     current = current->next) {
-			    my_dispose_list(files->token, &Rc_head);
-			    my_dispose_list(files->token, &Rcn_head);
+				/* Add Rc pixels to Ri */
+				Rc_tail->next = Ri_head;
+				Ri_head = Rc_head;
 
-			    /* get membership of neighbor segment */
-			    Rc_count = 1;
-			    newpixel =
-				(struct pixels *)link_new(files->token);
-			    newpixel->next = NULL;
-			    newpixel->row = current->row;
-			    newpixel->col = current->col;
-			    Rc_head = Rc_tail = newpixel;
-			    find_segment_neighbors(&Rc_head, &Rcn_head, &Rc_count, files, functions);	/* just to get members, not looking at neighbors now */
-			    merge_values(Ri_head, Rc_head, Ri_count,
-					 Rc_count, files);
+				/*to consider, recurse?  Check all Rcn neighbors if they are very close? */
 
-			    /* Add Rc pixels to Ri */
-			    Rc_tail->next = Ri_head;
-			    Ri_head = Rc_head;
-
-			    /*to consider, recurse?  Check all Rcn neighbors if they are very close? */
-
-			    Rc_head = NULL;
-			    my_dispose_list(files->token, &Rcn_head);
-			}
-			my_dispose_list(files->token, &Rclose_head);
+				Rc_head = NULL;
+				my_dispose_list(files->token, &Rcn_head);
+			    }
+			    my_dispose_list(files->token, &Rclose_head);
 #endif
 
-			/* check if we have a bestn that is valid to use to look at Rk */
-			if (Ri_bestn != NULL) {
-			    if ((functions->limited == TRUE) && !
-				(FLAG_GET
-				 (files->candidate_flag,
-				  Ri_bestn->row, Ri_bestn->col))) {
-				/* this check is important:
-				 * best neighbor is not a valid candidate, was already merged earlier in this time step */
-				Ri_bestn = NULL;
+			    /* check if we have a bestn that is valid to use to look at Rk */
+			    if (Ri_bestn != NULL) {
+				if ((functions->limited == TRUE) && !
+				    (FLAG_GET
+				     (files->candidate_flag,
+				      Ri_bestn->row, Ri_bestn->col))) {
+				    /* this check is important:
+				     * best neighbor is not a valid candidate, was already merged earlier in this time step */
+				    Ri_bestn = NULL;
+				}
 			    }
-			}
-			if (Ri_bestn != NULL && Ri_similarity < threshold) {	/* small TODO: should this be < or <= for threshold? */
-			    /* Rk starts from Ri's best neighbor */
-			    if (Rk_head) {
-				G_warning(_("Rk_head is not NULL!"));
-				my_dispose_list(files->token, &Rk_head);
-			    }
-			    if (Rkn_head) {
-				G_warning(_("Rkn_head is not NULL!"));
-				my_dispose_list(files->token, &Rkn_head);
-			    }
-			    Rk_count = 1;
-			    newpixel =
-				(struct pixels *)link_new(files->token);
-			    newpixel->next = NULL;
-			    newpixel->row = Ri_bestn->row;
-			    newpixel->col = Ri_bestn->col;
-			    Rk_head = newpixel;
+			    if (Ri_bestn != NULL && Ri_similarity < threshold) {	/* small TODO: should this be < or <= for threshold? */
+				/* Rk starts from Ri's best neighbor */
+				if (Rk_head) {
+				    G_warning(_("Rk_head is not NULL!"));
+				    my_dispose_list(files->token, &Rk_head);
+				}
+				if (Rkn_head) {
+				    G_warning(_("Rkn_head is not NULL!"));
+				    my_dispose_list(files->token, &Rkn_head);
+				}
+				Rk_count = 1;
+				newpixel =
+				    (struct pixels *)link_new(files->token);
+				newpixel->next = NULL;
+				newpixel->row = Ri_bestn->row;
+				newpixel->col = Ri_bestn->col;
+				newpixel->count_shared =
+				    Ri_bestn->count_shared;
+				Rk_head = newpixel;
 
-			    find_segment_neighbors(&Rk_head, &Rkn_head,
-						   &Rk_count, files,
-						   functions);
+				find_segment_neighbors(&Rk_head, &Rkn_head,
+						       &Rk_count, files,
+						       functions);
 
-			    /* ********  find Rk's most similar neighbor  ******** */
-			    Rk_similarity = Ri_similarity;	/*Ri gets first priority - ties won't change anything, so we'll accept Ri and Rk as mutually best neighbors */
-			    segment_get(&files->bands_seg, (void *)files->bands_val, Rk_head->row, Rk_head->col);	/* current segment values */
+				/* ********  find Rk's most similar neighbor  ******** */
+				Rk_similarity = Ri_similarity;	/*Ri gets first priority - ties won't change anything, so we'll accept Ri and Rk as mutually best neighbors */
+				segment_get(&files->bands_seg, (void *)files->bands_val, Rk_head->row, Rk_head->col);	/* current segment values */
 
-			    /* check similarity for each of Rk's neighbors */
-			    for (current = Rkn_head; current != NULL;
-				 current = current->next) {
-				tempsim =
-				    functions->calculate_similarity
-				    (Rk_head, current, files, functions);
+				/* check similarity for each of Rk's neighbors */
+				for (current = Rkn_head; current != NULL;
+				     current = current->next) {
+				    tempsim =
+					functions->calculate_similarity
+					(Rk_head, current, files, functions);
 
-				if (tempsim < Rk_similarity) {
-				    Rk_similarity = tempsim;
-				    break;	/* exit for Rk's neighbors loop here, we know that Ri and Rk aren't mutually best neighbors */
-				}
-			    }	/* have checked all of Rk's neighbors */
+				    if (tempsim < Rk_similarity) {
+					Rk_similarity = tempsim;
+					break;	/* exit for Rk's neighbors loop here, we know that Ri and Rk aren't mutually best neighbors */
+				    }
+				}	/* have checked all of Rk's neighbors */
 
-			    if (Rk_similarity == Ri_similarity) {	/* mutually most similar neighbors */
+				if (Rk_similarity == Ri_similarity) {	/* mutually most similar neighbors */
 
 #ifdef PROFILE
-				merge_start = clock();
+				    merge_start = clock();
 #endif
-				merge_values(Ri_head, Rk_head, Ri_count,
-					     Rk_count, files);
+				    merge_values(Ri_head, Rk_head, Ri_count,
+						 Rk_count, files);
 #ifdef PROFILE
-				merge_end = clock();
-				merge_lap =
-				    ((double)(merge_end - merge_start)) /
-				    CLOCKS_PER_SEC;
-				merge_accum += merge_lap;
-				fprintf(stdout, "merge time: %g\n",
-					merge_lap);
+				    merge_end = clock();
+				    merge_lap =
+					((double)(merge_end - merge_start)) /
+					CLOCKS_PER_SEC;
+				    merge_accum += merge_lap;
+				    fprintf(stdout, "merge time: %g\n",
+					    merge_lap);
 #endif
-				endflag = FALSE;	/* we've made at least one merge, so want another t iteration */
-			    }
-			    else {	/* they weren't mutually best neighbors */
+				    endflag = FALSE;	/* we've made at least one merge, so want another t iteration */
+				}
+				else {	/* they weren't mutually best neighbors */
 
-				/* checked Ri once, didn't find a mutually best neighbor, so remove all members of Ri from candidate pixels for this iteration */
+				    /* checked Ri once, didn't find a mutually best neighbor, so remove all members of Ri from candidate pixels for this iteration */
+				    set_candidate_flag(Ri_head, FALSE, files);
+
+				    if (FLAG_GET
+					(files->candidate_flag, Rk_head->row,
+					 Rk_head->col) &&
+					FLAG_GET(files->seeds_flag,
+						 Rk_head->row, Rk_head->col))
+					pathflag = TRUE;
+				}
+			    }	/* end if (Ri_bestn != NULL && Ri_similarity < threshold) */
+			    else {
+				/* no valid best neighbor for this Ri
+				 * exclude this Ri from further comparisons 
+				 * because we checked already Ri for a mutually best neighbor with all valid candidates
+				 * thus Ri can not be the mutually best neighbor later on during this pass
+				 * unfortunately this does happen sometimes */
 				set_candidate_flag(Ri_head, FALSE, files);
+			    }
 
-				if (FLAG_GET
-				    (files->candidate_flag, Rk_head->row,
-				     Rk_head->col) &&
-				    FLAG_GET(files->seeds_flag, Rk_head->row,
-					     Rk_head->col))
-				    pathflag = TRUE;
-			    }
-			}	/* end if (Ri_bestn != NULL && Ri_similarity < threshold) */
-			else {
-			    /* no valid best neighbor for this Ri
-			     * exclude this Ri from further comparisons 
-			     * because we checked already Ri for a mutually best neighbor with all valid candidates
-			     * thus Ri can not be the mutually best neighbor later on during this pass
-			     * unfortunately this does happen sometimes */
+			}	/* end if(Rin_head != NULL) */
+			else {	/* Ri didn't have a neighbor */
+			    G_debug(4, "Segment had no neighbors");
 			    set_candidate_flag(Ri_head, FALSE, files);
 			}
 
-		    }		/* end if(Rin_head != NULL) */
-		    else {	/* Ri didn't have a neighbor */
-			G_debug(4, "Segment had no neighbors");
-			set_candidate_flag(Ri_head, FALSE, files);
-		    }
+			if (pathflag) {	/*initialize Ri, Rin, using Rk as Ri and Rkn as Rin. */
 
-		    if (pathflag) {	/*initialize Ri, Rin, Rk, Rin using Rk as Ri. */
-			/* For the next iteration, lets start with Rk as the focus segment */
-			if (functions->path == TRUE) {
-			    Ri_count = Rk_count;
-			    Rk_count = 0;
-			    my_dispose_list(files->token, &Ri_head);
-			    Ri_head = Rk_head;
-			    Rk_head = NULL;
-			    if (Rkn_head != NULL) {
-				my_dispose_list(files->token, &Rin_head);
-				Rin_head = Rkn_head;
-				Rkn_head = NULL;
+			    /* For the next iteration, lets start with Rk as the focus segment */
+			    if (functions->path == TRUE) {
+				Ri_count = Rk_count;
+				Rk_count = 0;
+				my_dispose_list(files->token, &Ri_head);
+				Ri_head = Rk_head;
+				Rk_head = NULL;
+				if (Rkn_head != NULL) {
+				    my_dispose_list(files->token, &Rin_head);
+				    Rin_head = Rkn_head;
+				    Rkn_head = NULL;
+				}
+				else
+				    my_dispose_list(files->token, &Rin_head);
 			    }
 			    else
-				my_dispose_list(files->token, &Rin_head);
+				pathflag = FALSE;
+
 			}
-			else
-			    pathflag = FALSE;
-
-		    }
-
-		}		/*end pathflag do loop */
-	    }		/*end if pixel is candidate and seed pixel */
-	}			/*next column (or next z) */
-#ifdef ROWMAJOR
-    }				/*next row */
-#endif
+		    }		/*end pathflag do loop */
+		}		/*end if pixel is candidate and seed pixel */
+	    }			/*next column */
+	}			/*next row */
 #ifdef PROFILE
-    pass_end = clock();
-    fprintf(stdout, "pass %d took: %g\n", t,
-	    ((double)(pass_end - pass_start)) / CLOCKS_PER_SEC);
+	pass_end = clock();
+	fprintf(stdout, "pass %d took: %g\n", t,
+		((double)(pass_end - pass_start)) / CLOCKS_PER_SEC);
 #endif
 
-
-    /* finished one iteration over entire raster */
-    t++;
+	/* finished one iteration over entire raster */
+	t++;
     }
-    while (t <= functions->end_t && endflag == FALSE) ;	/*end t loop, either reached max iterations or didn't merge any segments */
+    while (t <= functions->end_t && endflag == FALSE);	/*end t loop, either reached max iterations or didn't merge any segments */
 
     if (t == 2 && files->bounds_map == NULL)
 	G_warning(_("No segments were created. Verify threshold and region settings."));
@@ -648,9 +612,9 @@
 		    }		/* end if neighbors found */
 		    else {	/* no neighbors were found */
 			if (files->bounds_map == NULL)
-			G_warning
-			    (_("no neighbors found, this means only one segment was created."));
-			
+			    G_warning
+				(_("no neighbors found, this means only one segment was created."));
+
 			set_candidate_flag(Ri_head, FALSE, files);
 		    }
 		}		/* end if pixel is candidate pixel */
@@ -668,46 +632,104 @@
 	    ((double)(end - start) / CLOCKS_PER_SEC));
 #endif
     return TRUE;
-    }
+}
 
     /* perimeter todo, for now will return borderPixels instead of passing a pointer, I saw mentioned that each parameter slows down the function call? */
     /* perimeter todo, My first impression is that the borderPixels count is ONLY needed for the case of initial seeds, and not used later on.  Another reason to split the function... */
-    int find_segment_neighbors(struct pixels **R_head,
-			       struct pixels **neighbors_head, int *seg_count,
-			       struct files *files,
-			       struct functions *functions)
-    {
-	int n, borderPixels, current_seg_ID, R_iseg = -1;
-	struct pixels *newpixel, *current, *to_check, tree_pix;	/* need to check the pixel neighbors of to_check */
-	int pixel_neighbors[8][2];
-	struct RB_TREE *no_check_tree;	/* pixels that should no longer be checked on this current find_neighbors() run */
-	struct RB_TREE *known_iseg;
+int find_segment_neighbors(struct pixels **R_head,
+			   struct pixels **neighbors_head, int *seg_count,
+			   struct files *files, struct functions *functions)
+{
+    int n, borderPixels, current_seg_ID, temp_ID, R_iseg = -1;	/* borderPixels is just used for open_files to get a starting perimeter with seeded segmentation. */
+    struct pixels *newpixel, *current, *to_check, tree_pix, *pixel_iter;	/* need to check the pixel neighbors of to_check */
+    int pixel_neighbors[8][2];
+    struct RB_TREE *no_check_tree;	/* pixels that should no longer be checked on this current find_neighbors() run */
+    struct RB_TREE *known_iseg;
 
 #ifdef DEBUG
-	struct RB_TRAV trav;
+    struct RB_TRAV trav;
 #endif
 
-	/* speed enhancement, any time savings to move any variables to files (mem allocation once in open_files) */
+    /* speed enhancement, any time savings to move any variables to files (mem allocation once in open_files) */
 
-	/* neighbor list will be a listing of pixels that are neighbors, but will be limited to just one pixel from each neighboring segment.
-	 * */
+    /* neighbor list will be a listing of pixels that are neighbors, but will be limited to just one pixel from each neighboring segment.
+     * */
 
-	/* parameter: R, current segment membership, could be single pixel (incomplete list) or list of pixels.
-	 * parameter: neighbors/Rin/Rik, neighbor pixels, could have a list already, or could be empty
-	 * functions->num_pn  int, 4 or 8, for number of pixel neighbors 
-	 * */
+    /* parameter: R, current segment membership, could be single pixel (incomplete list) or list of pixels.
+     * parameter: neighbors/Rin/Rik, neighbor pixels, could have a list already, or could be empty
+     * functions->num_pn  int, 4 or 8, for number of pixel neighbors 
+     * */
 
 
-	/* *** initialize data *** */
-	borderPixels = 0;
+    /* *** initialize data *** */
+    borderPixels = 0;
 
-	segment_get(&files->iseg_seg, &R_iseg, (*R_head)->row,
-		    (*R_head)->col);
+    segment_get(&files->iseg_seg, &R_iseg, (*R_head)->row, (*R_head)->col);
 
-	if (R_iseg == 0) {	/* if seeds were provided, this is just a single non-seed pixel, only return neighbors that are segments or seeds */
+    if (R_iseg == 0) {		/* if seeds were provided, this is just a single non-seed pixel, only return neighbors that are segments or seeds */
 
-	    functions->find_pixel_neighbors((*R_head)->row, (*R_head)->col,
+	functions->find_pixel_neighbors((*R_head)->row, (*R_head)->col,
+					pixel_neighbors, files);
+	for (n = 0; n < functions->num_pn; n++) {
+
+	    /* skip pixel if out of computational area or null */
+	    if (pixel_neighbors[n][0] < files->minrow ||
+		pixel_neighbors[n][0] >= files->maxrow ||
+		pixel_neighbors[n][1] < files->mincol ||
+		pixel_neighbors[n][1] >= files->maxcol ||
+		FLAG_GET(files->null_flag, pixel_neighbors[n][0],
+			 pixel_neighbors[n][1])
+		)
+		continue;
+
+	    segment_get(&files->iseg_seg, &current_seg_ID,
+			pixel_neighbors[n][0], pixel_neighbors[n][1]);
+
+	    if (current_seg_ID > 0) {
+		newpixel = (struct pixels *)link_new(files->token);
+		newpixel->next = *neighbors_head;	/*point the new pixel to the current first pixel */
+		newpixel->row = pixel_neighbors[n][0];
+		newpixel->col = pixel_neighbors[n][1];
+		*neighbors_head = newpixel;	/*change the first pixel to be the new pixel. */
+	    }
+	    borderPixels++;	/* increment for all non null pixels *//* TODO perimeter: OK to ignore NULL cells? */
+	}
+
+    }
+    else {			/*normal processing, look for all adjacent pixels or segments */
+	no_check_tree = rbtree_create(compare_pixels, sizeof(struct pixels));
+	known_iseg = rbtree_create(compare_ids, sizeof(int));
+	to_check = NULL;
+
+	/* Copy R in to_check and no_check data structures (don't expand them if we find them again) */
+
+	for (current = *R_head; current != NULL; current = current->next) {
+	    /* put in to_check linked list */
+	    newpixel = (struct pixels *)link_new(files->token);
+	    newpixel->next = to_check;	/*point the new pixel to the current first pixel */
+	    newpixel->row = current->row;
+	    newpixel->col = current->col;
+	    to_check = newpixel;	/*change the first pixel to be the new pixel. */
+
+	    /* put in no_check tree */
+	    tree_pix.row = current->row;
+	    tree_pix.col = current->col;
+	    if (rbtree_insert(no_check_tree, &tree_pix) == 0)	/* don't check it again */
+		G_warning(_("could not insert data into tree, out of memory?"));
+	}
+
+	while (to_check != NULL) {	/* removing from to_check list as we go, NOT iterating over the list. */
+
+	    functions->find_pixel_neighbors(to_check->row,
+					    to_check->col,
 					    pixel_neighbors, files);
+
+	    /* Done using this to_check pixels coords, remove from list */
+	    current = to_check;	/* temporary store the old head */
+	    to_check = to_check->next;	/*head now points to the next element in the list */
+	    link_dispose(files->token, (VOID_T *) current);
+
+	    /* for each pixel neighbors, check if they should be processed, check segment ID, and add to appropriate lists */
 	    for (n = 0; n < functions->num_pn; n++) {
 
 		/* skip pixel if out of computational area or null */
@@ -720,238 +742,254 @@
 		    )
 		    continue;
 
-		segment_get(&files->iseg_seg, &current_seg_ID,
-			    pixel_neighbors[n][0], pixel_neighbors[n][1]);
+		tree_pix.row = pixel_neighbors[n][0];
+		tree_pix.col = pixel_neighbors[n][1];
 
-		if (current_seg_ID > 0) {
-		    newpixel = (struct pixels *)link_new(files->token);
-		    newpixel->next = *neighbors_head;	/*point the new pixel to the current first pixel */
-		    newpixel->row = pixel_neighbors[n][0];
-		    newpixel->col = pixel_neighbors[n][1];
-		    *neighbors_head = newpixel;	/*change the first pixel to be the new pixel. */
-		}
-		borderPixels++;	/* increment for all non null pixels *//* TODO perimeter: OK to ignore NULL cells? */
-	    }
+		if (rbtree_find(no_check_tree, &tree_pix) == FALSE) {	/* want to check this neighbor */
+		    segment_get(&files->iseg_seg, &current_seg_ID,
+				pixel_neighbors[n][0], pixel_neighbors[n][1]);
 
-	}
-	else {			/*normal processing, look for all adjacent pixels or segments */
-	    no_check_tree =
-		rbtree_create(compare_pixels, sizeof(struct pixels));
-	    known_iseg = rbtree_create(compare_ids, sizeof(int));
-	    to_check = NULL;
+		    rbtree_insert(no_check_tree, &tree_pix);	/* don't check it again */
 
-	    /* Copy R in to_check and no_check data structures (don't expand them if we find them again) */
+		    if (current_seg_ID == R_iseg) {	/* pixel is member of current segment, add to R */
+			/* put pixel_neighbor[n] in Ri */
+			newpixel = (struct pixels *)link_new(files->token);
+			newpixel->next = *R_head;	/*point the new pixel to the current first pixel */
+			newpixel->row = pixel_neighbors[n][0];
+			newpixel->col = pixel_neighbors[n][1];
+			newpixel->count_shared = (*R_head)->count_shared;	/* Needed for Rk - to have the head remember how many pixels were shared with Ri. Might be a better way to do this... */
+			*R_head = newpixel;	/*change the first pixel to be the new pixel. */
+			*seg_count = *seg_count + 1;	/* zero index... Ri[0] had first pixel and set count =1.  increment after save data. */
 
-	    for (current = *R_head; current != NULL; current = current->next) {
-		/* put in to_check linked list */
-		newpixel = (struct pixels *)link_new(files->token);
-		newpixel->next = to_check;	/*point the new pixel to the current first pixel */
-		newpixel->row = current->row;
-		newpixel->col = current->col;
-		to_check = newpixel;	/*change the first pixel to be the new pixel. */
+			/* put pixel_neighbor[n] in to_check -- want to check this pixels neighbors */
+			newpixel = (struct pixels *)link_new(files->token);
+			newpixel->next = to_check;	/*point the new pixel to the current first pixel */
+			newpixel->row = pixel_neighbors[n][0];
+			newpixel->col = pixel_neighbors[n][1];
+			to_check = newpixel;	/*change the first pixel to be the new pixel. */
 
-		/* put in no_check tree */
-		tree_pix.row = current->row;
-		tree_pix.col = current->col;
-		if (rbtree_insert(no_check_tree, &tree_pix) == 0)	/* don't check it again */
-		    G_warning(_("could not insert data into tree, out of memory?"));
-	    }
+		    }
+		    else {	/* segment id's were different */
+			borderPixels++;	/* increment for all non null pixels that are non in no-check or R_iseg TODO perimeter: move this to include pixels in no-check ??? */
 
-	    while (to_check != NULL) {	/* removing from to_check list as we go, NOT iterating over the list. */
+			if (!rbtree_find(known_iseg, &current_seg_ID)) {	/* we don't have any neighbors yet from this segment */
+			    if (current_seg_ID != 0)
+				/* with seeds, non seed pixels are defaulted to zero.  Should we use null instead?? then could skip this check?  Or we couldn't insert it??? */
+				/* add to known neighbors list */
+				rbtree_insert(known_iseg, &current_seg_ID);
 
-		functions->find_pixel_neighbors(to_check->row,
-						to_check->col,
-						pixel_neighbors, files);
+			    /* put pixel_neighbor[n] in Rin */
+			    newpixel =
+				(struct pixels *)link_new(files->token);
+			    newpixel->next = *neighbors_head;	/*point the new pixel to the current first pixel */
+			    newpixel->row = pixel_neighbors[n][0];
+			    newpixel->col = pixel_neighbors[n][1];
+			    newpixel->count_shared = 1;	/*for shared perimeter */
+			    *neighbors_head = newpixel;	/*change the first pixel to be the new pixel. */
+			}
+			else {	/* todo perimeter we need to keep track of (and return!) a total count of neighbors pixels for each neighbor segment, to update the perimeter value in the similarity calculation. */
+			    /* todo perimeter: need to initalize this somewhere!!! */
+			    /* todo perimeter... need to find pixel with same segment ID....  countShared++;
+			     * hmmm,  Should we change the known_iseg tree to sort on segment ID...need to think of fast way to return this count?  with pixel?  or with something else? */
 
-		/* Done using this to_check pixels coords, remove from list */
-		current = to_check;	/* temporary store the old head */
-		to_check = to_check->next;	/*head now points to the next element in the list */
-		link_dispose(files->token, (VOID_T *) current);
+			    /* need to increment the count of border pixels for later know the shared border for calculating the new perimeter. */
+			    if (functions->radio_weight < 1) {	/* TODO: any reason to calculate this if we have a weight of 0 for the shape features? */
 
-		/* for each pixel neighbors, check if they should be processed, check segment ID, and add to appropriate lists */
-		for (n = 0; n < functions->num_pn; n++) {
+				for (pixel_iter = *neighbors_head;
+				     pixel_iter != NULL;
+				     pixel_iter = pixel_iter->next) {
 
-		    /* skip pixel if out of computational area or null */
-		    if (pixel_neighbors[n][0] < files->minrow ||
-			pixel_neighbors[n][0] >= files->maxrow ||
-			pixel_neighbors[n][1] < files->mincol ||
-			pixel_neighbors[n][1] >= files->maxcol ||
-			FLAG_GET(files->null_flag, pixel_neighbors[n][0],
-				 pixel_neighbors[n][1])
-			)
-			continue;
+				    segment_get(&files->iseg_seg,
+						&temp_ID, pixel_iter->row,
+						pixel_iter->col);
 
-		    tree_pix.row = pixel_neighbors[n][0];
-		    tree_pix.col = pixel_neighbors[n][1];
 
-		    if (rbtree_find(no_check_tree, &tree_pix) == FALSE) {	/* want to check this neighbor */
-			segment_get(&files->iseg_seg, &current_seg_ID,
-				    pixel_neighbors[n][0],
-				    pixel_neighbors[n][1]);
+				    if (temp_ID == current_seg_ID) {
+					pixel_iter->count_shared++;
+					break;
+				    }
+				}
+			    }
+			}
+		    }
 
-			rbtree_insert(no_check_tree, &tree_pix);	/* don't check it again */
 
-			if (current_seg_ID == R_iseg) {	/* pixel is member of current segment, add to R */
-			    /* put pixel_neighbor[n] in Ri */
-			    newpixel =
-				(struct pixels *)link_new(files->token);
-			    newpixel->next = *R_head;	/*point the new pixel to the current first pixel */
-			    newpixel->row = pixel_neighbors[n][0];
-			    newpixel->col = pixel_neighbors[n][1];
-			    *R_head = newpixel;	/*change the first pixel to be the new pixel. */
-			    *seg_count = *seg_count + 1;	/* zero index... Ri[0] had first pixel and set count =1.  increment after save data. */
+		}		/*end if for pixel_neighbor was in "don't check" list */
+		/* even if no_check tree, if we are using shape measurements we need to count if there is a shared border. */
+		else if (functions->radio_weight < 1) {
+		    segment_get(&files->iseg_seg, &current_seg_ID,
+				pixel_neighbors[n][0], pixel_neighbors[n][1]);
+		    if (current_seg_ID != R_iseg) {
+			for (pixel_iter = *neighbors_head;
+			     pixel_iter != NULL;
+			     pixel_iter = pixel_iter->next) {
 
-			    /* put pixel_neighbor[n] in to_check -- want to check this pixels neighbors */
-			    newpixel =
-				(struct pixels *)link_new(files->token);
-			    newpixel->next = to_check;	/*point the new pixel to the current first pixel */
-			    newpixel->row = pixel_neighbors[n][0];
-			    newpixel->col = pixel_neighbors[n][1];
-			    to_check = newpixel;	/*change the first pixel to be the new pixel. */
+			    segment_get(&files->iseg_seg,
+					&temp_ID, pixel_iter->row,
+					pixel_iter->col);
 
-			}
-			else {	/* segment id's were different */
-			    borderPixels++;	/* increment for all non null pixels that are non in no-check or R_iseg TODO perimeter: move this to include pixels in no-check ??? */
 
-			    if (!rbtree_find(known_iseg, &current_seg_ID)) {	/* we don't have any neighbors yet from this segment */
-				if (current_seg_ID != 0)
-				    /* with seeds, non seed pixels are defaulted to zero.  Should we use null instead?? then could skip this check?  Or we couldn't insert it??? */
-				    /* add to known neighbors list */
-				    rbtree_insert(known_iseg,
-						  &current_seg_ID);
-
-				/* put pixel_neighbor[n] in Rin */
-				newpixel =
-				    (struct pixels *)link_new(files->token);
-				newpixel->next = *neighbors_head;	/*point the new pixel to the current first pixel */
-				newpixel->row = pixel_neighbors[n][0];
-				newpixel->col = pixel_neighbors[n][1];
-				*neighbors_head = newpixel;	/*change the first pixel to be the new pixel. */
+			    if (temp_ID == current_seg_ID) {
+				pixel_iter->count_shared++;
+				break;
 			    }
-			    else {	/* todo perimeter we need to keep track of (and return!) a total count of neighbors pixels for each neighbor segment, to update the perimeter value in the similarity calculation. */
-				/* todo perimeter: need to initalize this somewhere!!! */
-				/* todo perimeter... need to find pixel with same segment ID....  countShared++;
-				 * hmmm,  Should we change the known_iseg tree to sort on segment ID...need to think of fast way to return this count?  with pixel?  or with something else? */
-			    }
 			}
+		    }
+		}
 
+	    }			/* end for loop - next pixel neighbor */
+	}			/* end while to_check has more elements */
 
-		    }		/*end if for pixel_neighbor was in "don't check" list */
-		}		/* end for loop - next pixel neighbor */
-	    }			/* end while to_check has more elements */
-
-	    /* clean up */
-	    rbtree_destroy(no_check_tree);
-	    rbtree_destroy(known_iseg);
-	}
-	return borderPixels;
+	/* clean up */
+	rbtree_destroy(no_check_tree);
+	rbtree_destroy(known_iseg);
     }
+    return borderPixels;
+}
 
-    int find_four_pixel_neighbors(int p_row, int p_col,
-				  int pixel_neighbors[8][2],
-				  struct files *files)
-    {
-	/* Note: this will return neighbors outside of the raster boundary.
-	 * Check in the calling routine if the pixel should be processed.
-	 */
+int find_four_pixel_neighbors(int p_row, int p_col,
+			      int pixel_neighbors[8][2], struct files *files)
+{
+    /* Note: this will return neighbors outside of the raster boundary.
+     * Check in the calling routine if the pixel should be processed.
+     */
 
-	/* north */
-	pixel_neighbors[0][1] = p_col;
-	pixel_neighbors[0][0] = p_row - 1;
+    /* north */
+    pixel_neighbors[0][1] = p_col;
+    pixel_neighbors[0][0] = p_row - 1;
 
-	/* east */
-	pixel_neighbors[1][0] = p_row;
-	pixel_neighbors[1][1] = p_col + 1;
+    /* east */
+    pixel_neighbors[1][0] = p_row;
+    pixel_neighbors[1][1] = p_col + 1;
 
-	/* south */
-	pixel_neighbors[2][1] = p_col;
-	pixel_neighbors[2][0] = p_row + 1;
+    /* south */
+    pixel_neighbors[2][1] = p_col;
+    pixel_neighbors[2][0] = p_row + 1;
 
-	/* west */
-	pixel_neighbors[3][0] = p_row;
-	pixel_neighbors[3][1] = p_col - 1;
+    /* west */
+    pixel_neighbors[3][0] = p_row;
+    pixel_neighbors[3][1] = p_col - 1;
 
-	return TRUE;
-    }
+    return TRUE;
+}
 
-    int find_eight_pixel_neighbors(int p_row, int p_col,
-				   int pixel_neighbors[8][2],
-				   struct files *files)
-    {
-	/* get the 4 orthogonal neighbors: */
-	find_four_pixel_neighbors(p_row, p_col, pixel_neighbors, files);
+int find_eight_pixel_neighbors(int p_row, int p_col,
+			       int pixel_neighbors[8][2], struct files *files)
+{
+    /* get the 4 orthogonal neighbors: */
+    find_four_pixel_neighbors(p_row, p_col, pixel_neighbors, files);
 
-	/* and then the diagonals: */
+    /* and then the diagonals: */
 
-	/* north west */
-	pixel_neighbors[4][0] = p_row - 1;
-	pixel_neighbors[4][1] = p_col - 1;
+    /* north west */
+    pixel_neighbors[4][0] = p_row - 1;
+    pixel_neighbors[4][1] = p_col - 1;
 
-	/* north east */
-	pixel_neighbors[5][0] = p_row - 1;
-	pixel_neighbors[5][1] = p_col + 1;
+    /* north east */
+    pixel_neighbors[5][0] = p_row - 1;
+    pixel_neighbors[5][1] = p_col + 1;
 
-	/* south east */
-	pixel_neighbors[6][0] = p_row + 1;
-	pixel_neighbors[6][1] = p_col + 1;
+    /* south east */
+    pixel_neighbors[6][0] = p_row + 1;
+    pixel_neighbors[6][1] = p_col + 1;
 
-	/* south west */
-	pixel_neighbors[7][0] = p_row + 1;
-	pixel_neighbors[7][1] = p_col - 1;
+    /* south west */
+    pixel_neighbors[7][0] = p_row + 1;
+    pixel_neighbors[7][1] = p_col - 1;
 
-	return TRUE;
-    }
+    return TRUE;
+}
 
     /* similarity / distance functions between two points based on their input raster values */
     /* assumes first point values already saved in files->bands_seg */
     /* speed enhancement: segment_get was already done for a[] values in the main function.  Could remove a[] from these parameters, reducing number of parameters in function call could provide a speed improvement. */
 
-    double calculate_euclidean_similarity(struct pixels *a, struct pixels *b,
-					  struct files *files,
-					  struct functions *functions)
-    {
-	double val = 0;
-	int n;
+double calculate_euclidean_similarity(struct pixels *a, struct pixels *b,
+				      struct files *files,
+				      struct functions *functions)
+{
+    double val = 0;
+    double smooth, compact, shape, PL;
+    int n;
 
-	/* get values for pixel b */
-	segment_get(&files->bands_seg, (void *)files->second_val, b->row,
-		    b->col);
+    /* get values for pixel b */
+    segment_get(&files->bands_seg, (void *)files->second_val, b->row, b->col);
 
-	/* euclidean distance, sum the square differences for each dimension */
-	for (n = 0; n < files->nbands; n++) {
-	    val =
-		val + (files->bands_val[n] -
-		       files->second_val[n]) * (files->bands_val[n] -
-						files->second_val[n]);
-	}
+    /* euclidean distance, sum the square differences for each dimension */
+    for (n = 0; n < files->nbands; n++) {
+	val =
+	    val + (files->bands_val[n] -
+		   files->second_val[n]) * (files->bands_val[n] -
+					    files->second_val[n]);
+    }
 
-	/* use squared distance, save the calculation time. We squared the similarity threshold earlier to allow for this. */
-	/* val = sqrt(val); */
+    /* use squared distance, save the calculation time. We squared the similarity threshold earlier to allow for this. */
+    /* val = sqrt(val); */
 
-	return val;
+    if (functions->radio_weight < 1) {
+	/*weight the raster data */
+	val = val * functions->radio_weight;
 
+	/*remaining part of similarity is based on the current shape of the object. */
+	/*I assume the idea is to add to the similarity information about the shape of the new segment if the two candidates were to be merged. */
+
+	PL = files->bands_val[files->nbands + 1] +
+	    files->second_val[files->nbands + 1] - b->count_shared;
+
+	/* compact = PL/sqrt(Npx) */
+
+	compact =
+	    PL / sqrt(files->bands_val[files->nbands] +
+		      files->second_val[files->nbands]);
+
+	/* smooth = PL/Pbbox */
+
+	smooth = PL /
+	    (2 *
+	     (max
+	      (files->bands_val[files->nbands + 2],
+	       files->second_val[files->nbands + 2])
+	      - min(files->bands_val[files->nbands + 3],
+		    files->second_val[files->nbands + 3]))
+	     +
+	     2 *
+	     (max
+	      (files->bands_val[files->nbands + 4],
+	       files->second_val[files->nbands + 4])
+	      - min(files->bands_val[files->nbands + 5],
+		    files->second_val[files->nbands + 5])));
+
+	shape =
+	    functions->smooth_weight * smooth + (1 - functions->smooth_weight)
+	    * compact;
+
+
+	val = val + (1 - functions->radio_weight) * shape;
     }
 
-    double calculate_manhattan_similarity(struct pixels *a, struct pixels *b,
-					  struct files *files,
-					  struct functions *functions)
-    {
-	double val = 0;
-	int n;
+    return val;
 
-	/* get values for pixel b */
-	segment_get(&files->bands_seg, (void *)files->second_val, b->row,
-		    b->col);
+}
 
-	/* Manhattan distance, sum the absolute difference between values for each dimension */
-	for (n = 0; n < files->nbands; n++) {
-	    val += fabs(files->bands_val[n] - files->second_val[n]);	/* speed enhancement: is fabs() is the "fast" way for absolute value calculations? */
-	}
+double calculate_manhattan_similarity(struct pixels *a, struct pixels *b,
+				      struct files *files,
+				      struct functions *functions)
+{
+    double val = 0;
+    int n;
 
-	return val;
+    /* get values for pixel b */
+    segment_get(&files->bands_seg, (void *)files->second_val, b->row, b->col);
 
+    /* Manhattan distance, sum the absolute difference between values for each dimension */
+    for (n = 0; n < files->nbands; n++) {
+	val += fabs(files->bands_val[n] - files->second_val[n]);	/* speed enhancement: is fabs() is the "fast" way for absolute value calculations? */
     }
 
+    return val;
+
+}
+
     /* TODO: add shape parameter...
      * 
      In the eCognition literature, we find that the key factor in the
@@ -972,192 +1010,197 @@
      object, and Pbbox the perimeter of the bounding box of the object.
      */
 
-    int merge_values(struct pixels *Ri_head, struct pixels *Rk_head,
-		     int Ri_count, int Rk_count, struct files *files)
-    {
-	int n, Ri_iseg, Rk_iseg;
-	struct pixels *current;
+int merge_values(struct pixels *Ri_head, struct pixels *Rk_head,
+		 int Ri_count, int Rk_count, struct files *files)
+{
+    int n, Ri_iseg, Rk_iseg;
+    struct pixels *current;
 
-	/*get input values *//*speed enhancement: Confirm if we can assume we already have bands_val for Ri, so don't need to segment_get() again?  note...current very_close implementation requires getting this value again... */
-	segment_get(&files->bands_seg, (void *)files->bands_val, Ri_head->row,
-		    Ri_head->col);
-	segment_get(&files->bands_seg, (void *)files->second_val,
-		    Rk_head->row, Rk_head->col);
+    /*get input values *//*speed enhancement: Confirm if we can assume we already have bands_val for Ri, so don't need to segment_get() again?  note...current very_close implementation requires getting this value again... */
+    segment_get(&files->bands_seg, (void *)files->bands_val, Ri_head->row,
+		Ri_head->col);
+    segment_get(&files->bands_seg, (void *)files->second_val,
+		Rk_head->row, Rk_head->col);
 
-	segment_get(&files->iseg_seg, &Rk_iseg, Rk_head->row, Rk_head->col);
-	segment_get(&files->iseg_seg, &Ri_iseg, Ri_head->row, Ri_head->col);
+    segment_get(&files->iseg_seg, &Rk_iseg, Rk_head->row, Rk_head->col);
+    segment_get(&files->iseg_seg, &Ri_iseg, Ri_head->row, Ri_head->col);
 
-	for (n = 0; n < files->nbands; n++) {
-	    files->bands_val[n] =
-		(files->bands_val[n] * Ri_count +
-		 files->second_val[n] * Rk_count) / (Ri_count + Rk_count);
-	}
+    for (n = 0; n < files->nbands; n++) {
+	files->bands_val[n] =
+	    (files->bands_val[n] * Ri_count +
+	     files->second_val[n] * Rk_count) / (Ri_count + Rk_count);
+    }
 
-	/* update segment number and candidate flag ==0 */
+    /* update shape parameters */
+
+    files->bands_val[n] += files->second_val[n];	/* area */
+    files->bands_val[n + 1] = files->bands_val[n + 1] + files->second_val[n + 1] - 2 * Rk_head->count_shared;	/* Perimeter Length */
+    files->bands_val[n + 2] = max(files->bands_val[n + 2], files->second_val[n + 2]);	/*max col */
+    files->bands_val[n + 3] = min(files->bands_val[n + 3], files->second_val[n + 3]);	/*min col */
+    files->bands_val[n + 4] = max(files->bands_val[n + 4], files->second_val[n + 4]);	/*max row */
+    files->bands_val[n + 5] = min(files->bands_val[n + 5], files->second_val[n + 5]);	/*min row */
+
+    /* update segment number and candidate flag ==0 */
 #ifdef SIGNPOST
-	fprintf(stdout,
-		"merging Ri (pixel count): %d (%d) with Rk (count): %d (%d).\n",
-		Ri_iseg, Ri_count, Rk_iseg, Rk_count);
+    fprintf(stdout,
+	    "merging Ri (pixel count): %d (%d) with Rk (count): %d (%d).\n",
+	    Ri_iseg, Ri_count, Rk_iseg, Rk_count);
 #endif
 
-	/* for each member of Ri and Rk, write new average bands values and segment values */
-	for (current = Ri_head; current != NULL; current = current->next) {
-	    segment_put(&files->bands_seg, (void *)files->bands_val,
-			current->row, current->col);
-	    FLAG_UNSET(files->candidate_flag, current->row, current->col);	/*candidate pixel flag, only one merge allowed per t iteration */
-	}
-	for (current = Rk_head; current != NULL; current = current->next) {
-	    segment_put(&files->bands_seg, (void *)files->bands_val,
-			current->row, current->col);
-	    segment_put(&files->iseg_seg, &Ri_iseg, current->row,
-			current->col);
-	    FLAG_UNSET(files->candidate_flag, current->row, current->col);
-	}
+    /* for each member of Ri and Rk, write new average bands values and segment values */
+    for (current = Ri_head; current != NULL; current = current->next) {
+	segment_put(&files->bands_seg, (void *)files->bands_val,
+		    current->row, current->col);
+	FLAG_UNSET(files->candidate_flag, current->row, current->col);	/*candidate pixel flag, only one merge allowed per t iteration */
+    }
+    for (current = Rk_head; current != NULL; current = current->next) {
+	segment_put(&files->bands_seg, (void *)files->bands_val,
+		    current->row, current->col);
+	segment_put(&files->iseg_seg, &Ri_iseg, current->row, current->col);
+	FLAG_UNSET(files->candidate_flag, current->row, current->col);
+    }
 
-	/* merged two segments, decrement count if Rk was an actual segment (not a non-seed pixel) */
-	if (Rk_iseg > 0)
-	    files->nsegs--;
+    /* merged two segments, decrement count if Rk was an actual segment (not a non-seed pixel) */
+    if (Rk_iseg > 0)
+	files->nsegs--;
 
-	return TRUE;
-    }
+    return TRUE;
+}
 
     /* calculates and stores the mean value for all pixels in a list, assuming they are all in the same segment */
-    int merge_pixels(struct pixels *R_head, int borderPixels,
-		     struct files *files)
-    {
-	int n, count = 0;
-	struct pixels *current;
+int merge_pixels(struct pixels *R_head, int borderPixels, struct files *files)
+{
+    int n, count = 0;
+    struct pixels *current;
 
-	/* Note: using files->bands_val for current pixel values, and files->second_val for the accumulated value */
+    /* Note: using files->bands_val for current pixel values, and files->second_val for the accumulated value */
 
-	/* initialize second_val */
-	for (n = 0; n < files->nbands; n++) {
-	    files->second_val[n] = 0;
-	}
+    /* initialize second_val */
+    for (n = 0; n < files->nbands; n++) {
+	files->second_val[n] = 0;
+    }
 
-	if (R_head->next != NULL) {
-	    /* total up bands values for all pixels */
-	    for (current = R_head; current != NULL; current = current->next) {
-		segment_get(&files->bands_seg, (void *)files->bands_val,
-			    current->row, current->col);
-		for (n = 0; n < files->nbands; n++) {
-		    files->second_val[n] += files->bands_val[n];
-		}
-		count++;
-	    }
-
-	    /* calculate the mean */
+    if (R_head->next != NULL) {
+	/* total up bands values for all pixels */
+	for (current = R_head; current != NULL; current = current->next) {
+	    segment_get(&files->bands_seg, (void *)files->bands_val,
+			current->row, current->col);
 	    for (n = 0; n < files->nbands; n++) {
-		files->second_val[n] = files->second_val[n] / count;
+		files->second_val[n] += files->bands_val[n];
 	    }
+	    count++;
+	}
 
-	    /* add in the shape values */
-	    files->bands_val[files->nbands] = count;	/* area (Num Pixels) */
-	    files->bands_val[files->nbands + 1] = borderPixels;	/* Perimeter Length *//* todo perimeter, not exact for edges...close enough for now? */
+	/* calculate the mean */
+	for (n = 0; n < files->nbands; n++) {
+	    files->second_val[n] = files->second_val[n] / count;
+	}
 
-	    /* save the results */
-	    for (current = R_head; current != NULL; current = current->next) {
-		segment_put(&files->bands_seg, (void *)files->second_val,
-			    current->row, current->col);
-	    }
+	/* add in the shape values */
+	files->bands_val[files->nbands] = count;	/* area (Num Pixels) */
+	files->bands_val[files->nbands + 1] = borderPixels;	/* Perimeter Length *//* todo perimeter, not exact for edges...close enough for now? */
 
+	/* save the results */
+	for (current = R_head; current != NULL; current = current->next) {
+	    segment_put(&files->bands_seg, (void *)files->second_val,
+			current->row, current->col);
 	}
 
-	return TRUE;
     }
 
+    return TRUE;
+}
+
     /* besides setting flag, also increments how many pixels remain to be processed */
-    int set_candidate_flag(struct pixels *head, int value,
-			   struct files *files)
-    {
-	/* head is linked list of pixels, value is new value of flag */
-	struct pixels *current;
+int set_candidate_flag(struct pixels *head, int value, struct files *files)
+{
+    /* head is linked list of pixels, value is new value of flag */
+    struct pixels *current;
 
-	for (current = head; current != NULL; current = current->next) {
+    for (current = head; current != NULL; current = current->next) {
 
 
-	    if (value == FALSE) {
-		FLAG_UNSET(files->candidate_flag, current->row, current->col);
-	    }
-	    else if (value == TRUE) {
-		FLAG_SET(files->candidate_flag, current->row, current->col);
-	    }
-	    else
-		G_fatal_error
-		    (_("programming bug, helper function called with invalid argument"));
+	if (value == FALSE) {
+	    FLAG_UNSET(files->candidate_flag, current->row, current->col);
 	}
-	return TRUE;
+	else if (value == TRUE) {
+	    FLAG_SET(files->candidate_flag, current->row, current->col);
+	}
+	else
+	    G_fatal_error
+		(_("programming bug, helper function called with invalid argument"));
     }
+    return TRUE;
+}
 
     /* let memory manager know space is available again and reset head to NULL */
-    int my_dispose_list(struct link_head *token, struct pixels **head)
-    {
-	struct pixels *current;
+int my_dispose_list(struct link_head *token, struct pixels **head)
+{
+    struct pixels *current;
 
-	while ((*head) != NULL) {
-	    current = *head;	/* remember "old" head */
-	    *head = (*head)->next;	/* move head to next pixel */
-	    link_dispose(token, (VOID_T *) current);	/* remove "old" head */
-	}
-
-	return TRUE;
+    while ((*head) != NULL) {
+	current = *head;	/* remember "old" head */
+	*head = (*head)->next;	/* move head to next pixel */
+	link_dispose(token, (VOID_T *) current);	/* remove "old" head */
     }
 
+    return TRUE;
+}
+
     /* functions used by binary tree to compare items */
 
     /* speed enhancement: Maybe changing this would be an improvement? "static" was used in break_polygons.c  extern was suggested in docs.  */
 
-    int compare_ids(const void *first, const void *second)
-    {
-	int *a = (int *)first, *b = (int *)second;
+int compare_ids(const void *first, const void *second)
+{
+    int *a = (int *)first, *b = (int *)second;
 
-	if (*a < *b)
-	    return -1;
-	else if (*a > *b)
-	    return 1;
-	else if (*a == *b)
-	    return 0;
+    if (*a < *b)
+	return -1;
+    else if (*a > *b)
+	return 1;
+    else if (*a == *b)
+	return 0;
 
 
-	G_warning(_("find neighbors: Bug in binary tree!"));
-	return 1;
+    G_warning(_("find neighbors: Bug in binary tree!"));
+    return 1;
 
-    }
+}
 
-    int compare_pixels(const void *first, const void *second)
-    {
-	struct pixels *a = (struct pixels *)first, *b =
-	    (struct pixels *)second;
+int compare_pixels(const void *first, const void *second)
+{
+    struct pixels *a = (struct pixels *)first, *b = (struct pixels *)second;
 
-	if (a->row < b->row)
+    if (a->row < b->row)
+	return -1;
+    else if (a->row > b->row)
+	return 1;
+    else {
+	/* same row */
+	if (a->col < b->col)
 	    return -1;
-	else if (a->row > b->row)
+	else if (a->col > b->col)
 	    return 1;
-	else {
-	    /* same row */
-	    if (a->col < b->col)
-		return -1;
-	    else if (a->col > b->col)
-		return 1;
-	}
-	/* same row and col */
-	return 0;
     }
+    /* same row and col */
+    return 0;
+}
 
     /* Set candidate flag to true/1 or false/0 for all pixels in current processing area
      * checks for NULL flag and if it is in current "polygon" if a bounds map is given */
-    int set_all_candidate_flags(struct files *files)
-    {
-	int row, col;
+int set_all_candidate_flags(struct files *files)
+{
+    int row, col;
 
-	for (row = files->minrow; row < files->maxrow; row++) {
-	    for (col = files->mincol; col < files->maxcol; col++) {
-		if (!(FLAG_GET(files->null_flag, row, col))) {
-		    FLAG_SET(files->candidate_flag, row, col);
-		}
-		else
-		    FLAG_UNSET(files->candidate_flag, row, col);
+    for (row = files->minrow; row < files->maxrow; row++) {
+	for (col = files->mincol; col < files->maxcol; col++) {
+	    if (!(FLAG_GET(files->null_flag, row, col))) {
+		FLAG_SET(files->candidate_flag, row, col);
 	    }
+	    else
+		FLAG_UNSET(files->candidate_flag, row, col);
 	}
-	return TRUE;
     }
+    return TRUE;
+}

Added: grass-addons/grass7/imagery/i.segment/estimate_threshold.c
===================================================================
--- grass-addons/grass7/imagery/i.segment/estimate_threshold.c	                        (rev 0)
+++ grass-addons/grass7/imagery/i.segment/estimate_threshold.c	2012-10-26 14:03:33 UTC (rev 53561)
@@ -0,0 +1,94 @@
+/* Purpose: open input files and suggest a reasonable threshold */
+
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/imagery.h>
+#include "iseg.h"
+
+int estimate_threshold(char *image_group)
+{
+    double min, max, est_t;
+
+    /* check if the input image group is valid. */
+    check_group(image_group);
+
+    /* read the raster files to find the minimum and maximum values */
+    read_range(&min, &max, image_group);
+
+    /* perform the calculations to estimate the threshold */
+    est_t = calc_t(min, max);
+
+    /* set the output message and finish */
+    G_done_msg(_("Suggested threshold (if using -w flag and radioweight=1) is: <%g>"),
+	       est_t);
+    return TRUE;
+}
+
+/* function to validate that the user input is readable and contains files */
+int check_group(char *image_group)
+{
+    struct Ref Ref;		/* group reference list */
+
+    /* check that the input image group can be found */
+    if (!I_get_group_ref(image_group, &Ref))
+	G_fatal_error(_("Unable to read REF file for group <%s>"),
+		      image_group);
+
+    /* check that the input image group contains rasters */
+    if (Ref.nfiles <= 0)
+	G_fatal_error(_("Group <%s> contains no raster maps"), image_group);
+
+    return TRUE;
+}
+
+/* function to find the min and max values in all rasters of an image group */
+int read_range(double *min, double *max, char *image_group)
+{
+    int n;
+    struct FPRange fp_range;	/* range for a raster */
+    struct Ref Ref;		/* group reference list */
+    double candidate_min, candidate_max;	/* for min/max in each raster */
+
+    /* read the image group */
+    I_get_group_ref(image_group, &Ref);
+
+    /* initialize min/max with the first raster min/max */
+    if (Rast_read_fp_range(Ref.file[0].name, Ref.file[0].mapset, &fp_range) != 1)	/* returns -1 on error, 2 on empty range, quiting either way. */
+	G_fatal_error(_("No min/max found in raster map <%s>"),
+		      Ref.file[0].name);
+    else
+	Rast_get_fp_range_min_max(&fp_range, min, max);
+
+    /* check the min/max for any remaining rasters */
+    for (n = 1; n < Ref.nfiles; n++) {
+	if (Rast_read_fp_range(Ref.file[n].name, Ref.file[n].mapset, &fp_range) != 1) {	/* returns -1 on error, 2 on empty range, quiting either way. */
+	    G_fatal_error(_("No min/max found in raster map <%s>"),
+			  Ref.file[n].name);
+
+	    Rast_get_fp_range_min_max(&fp_range, &candidate_min,
+				      &candidate_max);
+
+	    if (candidate_min < *min)
+		*min = candidate_min;
+	    if (candidate_max > *max)
+		*max = candidate_max;
+	}
+    }
+
+    return TRUE;
+}
+
+/* function to calculate a suggested threshold based on the min and max values of the rasters */
+double calc_t(double min, double max)
+{
+    double t, fraction;
+
+    /* Empirical testing indicated 1 to 5% of the differences was a good starting point. */
+    fraction = 0.03;
+
+    /* TODO: allow the community to test this estimate, the formula can be updated based on their advice. */
+    t = fraction * (max - min);
+
+    return t;
+}

Modified: grass-addons/grass7/imagery/i.segment/i.segment.html
===================================================================
--- grass-addons/grass7/imagery/i.segment/i.segment.html	2012-10-26 13:37:15 UTC (rev 53560)
+++ grass-addons/grass7/imagery/i.segment/i.segment.html	2012-10-26 14:03:33 UTC (rev 53561)
@@ -2,7 +2,8 @@
 Image segmentation is the process of grouping similar pixels into 
 unique segments. Boundary and region based algorithms are described in 
 the literature, currently a region growing and merging algorithm is 
-implemented.  Each object found during the segmentation process is 
+implemented.  Each grouping (usually refered to as objects or segments) 
+found during the segmentation process is 
 given a unique ID and is a collection of contiguous pixels meeting some 
 criteria.  (Note the contrast with image classification, where 
 continuity and spatial characteristics are not important, but rather 
@@ -25,11 +26,14 @@
 made during a complete pass.
 
 <h3>Similarity and Threshold</h3>
-Similarity between objects is used to determine what objects are 
-merged.  The current implementation uses only the radiometric distance 
-between the two objects, but eventually region properties will be 
-included as well.  Thus a lower value for the similarity is a closer 
-match, with a similarity score of zero for identical pixels.
+The similarity between segments and unmerged pixels is used to determine 
+which are merged.  The Euclidean version uses the radiometric distance 
+between the two segments and also the shape characteristics.  
+The Manhatten calculations currently only uses only the radiometric distance 
+between the two segments, but eventually shape characteristics will be 
+included as well.  NOTE: Closer/smaller distances mean a lower value for 
+the similarity indicates a closer match, with a similarity score of zero 
+for identical pixels.
 <p>
 During normal processing, merges are only allowed when the similarity 
 between two segments is lower then the calculated threshold value.  
@@ -46,7 +50,23 @@
 The threshold will be multiplied by the number of rasters included in 
 the image group.  This will allow the same threshold to achieve similar 
 segmentation results when the number of rasters in the image group varies.
+<p>
+The -t flag will estimate the threshold, it is calculated at 3% of the
+range of data in the imagery group.  Initial empirical tests indicate
+threshold values of 1% to 5% are reasonable places to start.
 
+<h4>Calculation Formulas</h4>
+Both Euclidean and Manhattan distances use the normal definition, considering
+each raster in the image group as a dimension.
+
+Furthermore, the Euclidean calculation also takes into account the shape
+characteristics of the segments.  The normal distances are multiplied by 
+the input radiometric weight.  Next an additional contribution is added:
+(1-radioweight) * {smoothness * smoothness weight + compactness * (1-smoothness weight)},
+where compactness = the Perimeter Length / sqrt( Area ) and smoothness = Perimeter Length / 
+the Bounding Box.  The perimeter length is estimated as the number of pixel sides
+the segment has.
+
 <h3>Seeds</h3>
 The seeds map can be used to provide either seed pixels (random or 
 selected points from which to start the segmentation process) or seed 
@@ -133,17 +153,15 @@
           threshold=10 method=region_growing minsize=5 endt=5000
 </pre></div>
 <p>
-Processing the entire ortho image (over 9 million cells) took about a day.
+Processing the entire ortho image (over 9 million pixels) took about a day.
 
 <h2>TODO</h2>
 <h3>Functionality</h3>
 <ul>
-<li>Add shape characteristics (smoothness, compactness) to the 
-similarity measurement.  This will add two input parameters (weight of 
-radiometric to shape, and weight of compactness to smoothness.)</li>
+<li>Further testing of the shape characteristics (smoothness, compactness),
+if it looks good it should be added to the Manhatten option.
+<b>in progress</b></li>
 <li>Malahanobis distance for the similarity calculation.</li>
-<li>Estimating the threshold value. (1 to 5% of the max difference gave 
-me subjectively good results.)</li>
 </ul>
 <h3>Use of Segmentation Results</h3>
 <ul>
@@ -180,7 +198,7 @@
 <p>
 Information about
 <a href="http://grass.osgeo.org/wiki/Image_classification">classification 
-in GRASS GIS</a> is at available on the wiki.
+in GRASS GIS</a> is also available on the wiki.
 
 <h2>SEE ALSO</h2>
 <em>

Modified: grass-addons/grass7/imagery/i.segment/iseg.h
===================================================================
--- grass-addons/grass7/imagery/i.segment/iseg.h	2012-10-26 13:37:15 UTC (rev 53560)
+++ grass-addons/grass7/imagery/i.segment/iseg.h	2012-10-26 14:03:33 UTC (rev 53561)
@@ -30,7 +30,7 @@
     struct pixels *next;
     int row;
     int col;
-    int countShared;		/* todo perimeter: will hold the count how many pixels are shared on the Border Between Ri and Rk.  Not used for all pixels... see if this is an OK way to do this... */
+    int count_shared;		/* todo perimeter: will hold the count how many pixels are shared on the Border Between Ri and Rk.  Not used for all pixels... see if this is an OK way to do this... */
 };
 
 /* input and output files, as well as some other processing info */
@@ -50,6 +50,7 @@
 
     /* file processing */
     /* bands_seg is initialized with the input raster valuess, then is updated with current mean values for the segment. */
+    /* for now, also storing: Area, Perimeter, maxcol, mincol, maxrow, minrow */
     int nbands;			/* number of rasters in the image group */
     SEGMENT bands_seg, bounds_seg;	/* bands is for input, normal application is landsat bands, but other input can be included in the group. */
     double *bands_val;		/* array, to hold all input values at one pixel */
@@ -77,7 +78,7 @@
     int num_pn;			/* number of pixel neighbors  int, 4 or 8. */
     float threshold;		/* similarity threshold */
     int min_segment_size;	/* smallest number of pixels/cells allowed in a final segment */
-
+    float radio_weight, smooth_weight;	/* radiometric (bands) vs. shape and smoothness vs. compactness */
     /* Some function pointers to set in parse_args() */
     int (*find_pixel_neighbors) (int, int, int[8][2], struct files *);	/*parameters: row, col, pixel_neighbors */
     double (*calculate_similarity) (struct pixels *, struct pixels *, struct files *, struct functions *);	/*parameters: two pixels (each with row,col) to compare */
@@ -87,6 +88,7 @@
 
     int path;			/* flag if we are using Rk as next Ri for non-mutually best neighbor. */
     int limited;		/* flag if we are limiting merges to one per pass */
+    int estimate_threshold;	/* flag if we just want to estimate a suggested threshold value and exit. */
 
     /* todo: is there a fast way (and valid from an algorithm standpoint) to merge all neighbors that are within some small % of the treshold?
      * There is some code using "very_close" that is excluded with IFDEF
@@ -99,6 +101,11 @@
 #endif
 };
 
+/* main.c */
+int estimate_threshold(char *);
+int check_group(char *);
+int read_range(double *, double *, char *);
+double calc_t(double, double);
 
 /* parse_args.c */
 /* gets input from user, validates, and sets up functions */

Modified: grass-addons/grass7/imagery/i.segment/main.c
===================================================================
--- grass-addons/grass7/imagery/i.segment/main.c	2012-10-26 13:37:15 UTC (rev 53560)
+++ grass-addons/grass7/imagery/i.segment/main.c	2012-10-26 14:03:33 UTC (rev 53561)
@@ -15,12 +15,14 @@
  *               Library for the data files/tiling, so "iseg" (image segmentation)
  *               will be used to refer to the image segmentation.
  * 
+ * 				 First developed for GSoC 2012 with mentor: Markus Metz
  * 
  *****************************************************************************/
 
 #include <stdlib.h>
 #include <grass/gis.h>
 #include <grass/glocale.h>
+#include <grass/imagery.h>
 #include "iseg.h"
 
 int main(int argc, char *argv[])
@@ -40,22 +42,29 @@
     if (parse_args(argc, argv, &files, &functions) != TRUE)
 	G_fatal_error(_("Error in parse_args()"));
 
-    G_debug(1, "Main: starting open_files()");
-    if (open_files(&files, &functions) != TRUE)
-	G_fatal_error(_("Error in open_files()"));
+    /* check if we are doing normal processing or if the estimate threshold and exit flag has been selected */
 
-    G_debug(1, "Main: starting create_isegs()");
-    if (create_isegs(&files, &functions) != TRUE)
-	G_fatal_error(_("Error in create_isegs()"));
+    if (functions.estimate_threshold == FALSE) {
 
-    G_debug(1, "Main: starting write_output()");
-    if (write_output(&files) != TRUE)
-	G_fatal_error(_("Error in write_output()"));
+	G_debug(1, "Main: starting open_files()");
+	if (open_files(&files, &functions) != TRUE)
+	    G_fatal_error(_("Error in open_files()"));
 
-    G_debug(1, "Main: starting close_files()");
-    close_files(&files);
+	G_debug(1, "Main: starting create_isegs()");
+	if (create_isegs(&files, &functions) != TRUE)
+	    G_fatal_error(_("Error in create_isegs()"));
 
-    G_done_msg("Number of segments created: %d", files.nsegs);
+	G_debug(1, "Main: starting write_output()");
+	if (write_output(&files) != TRUE)
+	    G_fatal_error(_("Error in write_output()"));
 
+	G_debug(1, "Main: starting close_files()");
+	close_files(&files);
+
+	G_done_msg(_("Number of segments created: <%d>"), files.nsegs);
+    }
+    else {
+	estimate_threshold(files.image_group);
+    }
     exit(EXIT_SUCCESS);
 }

Modified: grass-addons/grass7/imagery/i.segment/open_files.c
===================================================================
--- grass-addons/grass7/imagery/i.segment/open_files.c	2012-10-26 13:37:15 UTC (rev 53560)
+++ grass-addons/grass7/imagery/i.segment/open_files.c	2012-10-26 14:03:33 UTC (rev 53561)
@@ -116,8 +116,10 @@
 
     /* size of each element to be stored */
 
-    inlen = sizeof(double) * Ref.nfiles;
+    /* save for bands, plus area, perimeter, bounding box. */
 
+    inlen = sizeof(double) * (Ref.nfiles + 6);
+
     /* when fine tuning, should be a power of 2 and not larger than 256 for speed reasons */
     srows = 64;
     scols = 64;
@@ -151,8 +153,11 @@
 
     /* save the area and perimeter as well */
     /* perimeter todo: currently saving this with the input DCELL values.  Better to have a second segment structure to save as integers ??? */
-    inlen = inlen + sizeof(double) * 2;
+    /* along with P and A, also saving the bounding box - min/max row/col */
 
+    /* Why was this being reset here??? commented out... 
+       inlen = inlen + sizeof(double) * 6; */
+
     files->bands_val = (double *)G_malloc(inlen);
     files->second_val = (double *)G_malloc(inlen);
 
@@ -212,10 +217,17 @@
 		else
 		    files->bands_val[n] = (inbuf[n][col] - min[n]) / (max[n] - min[n]);	/* scaled */
 	    }
+
+	    /* besides the user input rasters, also save the shape parameters */
+
 	    files->bands_val[Ref.nfiles] = 1;	/* area (just using the number of pixels) */
 	    files->bands_val[Ref.nfiles + 1] = 4;	/* Perimeter Length *//* todo perimeter, not exact for edges...close enough for now? */
+	    files->bands_val[Ref.nfiles + 2] = col;	/*max col */
+	    files->bands_val[Ref.nfiles + 3] = col;	/*min col */
+	    files->bands_val[Ref.nfiles + 4] = row;	/*max row */
+	    files->bands_val[Ref.nfiles + 5] = row;	/*min row */
+
 	    segment_put(&files->bands_seg, (void *)files->bands_val, row, col);	/* store input bands */
-
 	    if (null_check != -1) {	/*good pixel */
 		FLAG_UNSET(files->null_flag, row, col);	/*flag */
 		if (files->seeds_map != NULL) {

Modified: grass-addons/grass7/imagery/i.segment/parse_args.c
===================================================================
--- grass-addons/grass7/imagery/i.segment/parse_args.c	2012-10-26 13:37:15 UTC (rev 53560)
+++ grass-addons/grass7/imagery/i.segment/parse_args.c	2012-10-26 14:03:33 UTC (rev 53561)
@@ -11,7 +11,8 @@
 	       struct functions *functions)
 {
     struct Option *group, *seeds, *bounds, *output, *method, *similarity, *threshold, *min_segment_size, *endt;	/* Establish an Option pointer for each option */
-    struct Flag *diagonal, *weighted, *limited;	/* Establish a Flag pointer for each option */
+    struct Option *radio_weight, *smooth_weight;
+    struct Flag *estimate_threshold, *diagonal, *weighted, *limited;	/* Establish a Flag pointer for each option */
     struct Option *outband;	/* optional saving of segment data, until a seperate module is written */
 
 #ifdef VCLOSE
@@ -43,7 +44,7 @@
     similarity->required = YES;
     similarity->answer = "euclidean";
     similarity->options = "euclidean, manhattan";
-    similarity->description = _("Distance calculation method.");
+    similarity->description = _("Similarity calculation method.");
 
     min_segment_size = G_define_option();
     min_segment_size->key = "minsize";
@@ -67,8 +68,33 @@
 	_("Neighbors similarity lower then this fraction of the threshold will be merged without regard to any other processing rules.");
 #endif
 
+    /* for the weights of bands values vs. shape parameters, and weight of smoothness vs. compactness */
+
+    radio_weight = G_define_option();
+    radio_weight->key = "radioweight";
+    radio_weight->type = TYPE_DOUBLE;
+    radio_weight->required = YES;
+    radio_weight->answer = "0.9";
+    radio_weight->options = "0-1";
+    radio_weight->label =
+	_("Importance of radiometric (input raseters) values relative to shape");
+
+    smooth_weight = G_define_option();
+    smooth_weight->key = "smoothweight";
+    smooth_weight->type = TYPE_DOUBLE;
+    smooth_weight->required = YES;
+    smooth_weight->answer = "0.5";
+    smooth_weight->options = "0-1";
+    smooth_weight->label =
+	_("Importance of smoothness relative to compactness");
+
     /* optional parameters */
 
+    estimate_threshold = G_define_flag();
+    estimate_threshold->key = 't';
+    estimate_threshold->description =
+	_("Estimate a threshold based on input image group and exit.");
+
     diagonal = G_define_flag();
     diagonal->key = 'd';
     diagonal->description =
@@ -135,6 +161,12 @@
 
     files->image_group = group->answer;
 
+    functions->estimate_threshold = estimate_threshold->answer;
+
+    /* if we are just estimating a threshold, skip remaining input validation */
+    if (functions->estimate_threshold == TRUE)
+	return TRUE;
+
     if (G_legal_filename(output->answer) == TRUE)
 	files->out_name = output->answer;	/* name of output (segment ID) raster map */
     else
@@ -166,6 +198,9 @@
 
     functions->min_segment_size = atoi(min_segment_size->answer);
 
+    functions->radio_weight = atof(radio_weight->answer);
+    functions->smooth_weight = atof(smooth_weight->answer);
+
     if (diagonal->answer == FALSE) {
 	functions->find_pixel_neighbors = &find_four_pixel_neighbors;
 	functions->num_pn = 4;



More information about the grass-commit mailing list