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

svn_grass at osgeo.org svn_grass at osgeo.org
Mon Feb 11 01:55:04 PST 2013


Author: mmetz
Date: 2013-02-11 01:55:03 -0800 (Mon, 11 Feb 2013)
New Revision: 54997

Modified:
   grass-addons/grass7/imagery/i.segment.xl/create_isegs.c
   grass-addons/grass7/imagery/i.segment.xl/i.segment.xl.html
   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/parse_args.c
   grass-addons/grass7/imagery/i.segment.xl/write_output.c
Log:
i.segment.xl: prepare for shape, update manual

Modified: grass-addons/grass7/imagery/i.segment.xl/create_isegs.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/create_isegs.c	2013-02-11 09:54:47 UTC (rev 54996)
+++ grass-addons/grass7/imagery/i.segment.xl/create_isegs.c	2013-02-11 09:55:03 UTC (rev 54997)
@@ -16,6 +16,10 @@
 
 #define EPSILON 1.0e-12
 
+#define MAX(a,b) ( ((a)>(b)) ? (a) : (b) )
+#define MIN(a,b) ( ((a)<(b)) ? (a) : (b) )
+
+
 /* internal functions */
 static int merge_regions(struct ngbr_stats *, struct reg_stats *, /* Ri */
                          struct ngbr_stats *, struct reg_stats *, /* Rk */
@@ -110,7 +114,6 @@
     int have_bound;
     CELL current_bound, bounds_val;
 
-    G_debug(1, "Threshold: %g", globals->threshold);
     G_debug(1, "segmentation method: %d", globals->method);
 
     if (globals->bounds_map == NULL) {
@@ -209,7 +212,7 @@
 
     /* threshold calculation */
     alpha2 = globals->alpha * globals->alpha;
-    threshold = alpha2 * globals->threshold;
+    threshold = alpha2;
     G_debug(1, "Squared threshold: %g", threshold);
 
     /* make the divisor a constant ? */
@@ -265,7 +268,7 @@
 
 		/* find segment neighbors */
 		/* find Ri's best neighbor, clear candidate flag */
-		Ri_similarity = globals->threshold + 1;
+		Ri_similarity = 2;
 
 		Ri_rs.id = Ri.id;
 		fetch_reg_stats(Ri.row, Ri.col, &Ri_rs, globals);
@@ -298,8 +301,7 @@
 		    if (Ri.count < Rk.count)
 			smaller = Ri.count;
 
-		    adjthresh = pow(alpha2, 1. + (double) smaller / divisor) *
-				globals->threshold;
+		    adjthresh = pow(alpha2, 1. + (double) smaller / divisor);
 
 		    if (compare_double(Ri_similarity, adjthresh) == -1) {
 			G_debug(4, "Ri nn == 1");
@@ -329,7 +331,7 @@
 		    if (candidates_only && 
 		        !(FLAG_GET(globals->candidate_flag, Rk.row, Rk.col))) {
 
-			Ri_similarity = globals->threshold + 1;
+			Ri_similarity = 2;
 		    }
 
 		    candidates_only = TRUE;
@@ -341,7 +343,6 @@
 			G_debug(4, "Working with Rk");
 
 			/* find Rk's best neighbor, do not clear candidate flag */
-			/* Rk_similarity = globals->threshold + 1; */
 			Rk_similarity = Ri_similarity;
 			Rk_bestn_rs.count = 0;
 			/* Rk_rs is already complete */
@@ -356,7 +357,7 @@
 			    /* important for fp precision limit
 			     * because region stats may be calculated 
 			     * in two slightly different ways */
-			    if (fabs(Ri_similarity / Rk_similarity  - 1) > EPSILON)
+			    if (Ri_similarity - Rk_similarity > EPSILON)
 				do_merge = 0;
 			}
 			/* Ri has only one neighbor, merge */
@@ -370,8 +371,7 @@
 			    if (Ri.count < Rk.count)
 				smaller = Ri.count;
 
-			    adjthresh = pow(alpha2, 1. + (double) smaller / divisor) *
-					globals->threshold;
+			    adjthresh = pow(alpha2, 1. + (double) smaller / divisor);
 
 			    do_merge = 0;
 			    if (compare_double(Ri_similarity, adjthresh) == -1) {
@@ -402,8 +402,8 @@
 			    /* made a merge, need another iteration */
 			    n_merges++;
 
-			    Ri_similarity = globals->threshold + 1;
-			    Rk_similarity = globals->threshold + 1;
+			    Ri_similarity = 2;
+			    Rk_similarity = 2;
 
 			    /* we have checked the neighbors of Ri, Rk already
 			     * use faster version of finding the best neighbor
@@ -502,7 +502,7 @@
     if (globals->min_segment_size > 1) {
 	G_message(_("Merging segments smaller than %d cells"), globals->min_segment_size);
 
-	threshold = globals->alpha * globals->alpha * globals->threshold;
+	threshold = globals->alpha * globals->alpha;
 
 	flag_clear_all(globals->candidate_flag);
 	
@@ -560,7 +560,7 @@
 			do_merge = 1;
 
 		    Ri_nn = 0;
-		    Ri_similarity = globals->threshold + 1;
+		    Ri_similarity = 2;
 
 		    if (do_merge) {
 
@@ -844,10 +844,125 @@
 	val += diff * diff;
     } while (n--);
 
+    /* the return value should always be in the range 0 - 1 */
+    if (val <= 0)
+	return 0.;
+
+    val /= globals->max_diff;
+
+#ifdef _OR_SHAPE_
+    if (globals->shape_weight < 1)
+	val = val * globals->shape_weight + (1 - globals->shape_weight) *
+	      calculate_shape(rsi, rsk, nshared, globals);
+#endif
+
     return val;
 }
 
+double calculate_manhattan_similarity(struct ngbr_stats *Ri,
+                                      struct ngbr_stats *Rk,
+				      struct globals *globals)
+{
+    double val = 0.;
+    int n = globals->nbands - 1;
 
+    /* squared euclidean distance, sum the square differences for each dimension */
+    do {
+	val += Ri->mean[n] - Rk->mean[n];
+    } while (n--);
+
+    /* the return value should always be in the range 0 - 1 */
+    if (val <= 0)
+	return 0.;
+
+    val /= globals->max_diff;
+
+#ifdef _OR_SHAPE_
+    if (globals->shape_weight < 1)
+	val = val * globals->shape_weight + (1 - globals->shape_weight) *
+	      calculate_shape(rsi, rsk, nshared, globals);
+#endif
+
+    return val;
+}
+
+double calculate_shape(struct reg_stats *rsi, struct reg_stats *rsk,
+                       int nshared, struct globals *globals)
+{
+     /* 
+     In the eCognition literature, we find that the key factor in the
+     multi-scale segmentation algorithm used by Definiens is the scale
+     factor f:
+
+     f = W.Hcolor + (1 - W).Hshape
+     Hcolor = sum(b = 1:nbands)(Wb.SigmaB)
+     Hshape = Ws.Hcompact + (1 - Ws).Hsmooth
+     Hcompact = PL/sqrt(Npx)
+     Hsmooth = PL/Pbbox
+
+     Where 
+     W is a user-defined weight of importance of object radiometry vs
+     shape (usually .9 vs .1)
+     Wb is the weigh given to band B
+     SigmaB is the std dev of the object for band b
+     Ws is a user-defined weight giving the importance of compactedness vs smoothness
+     PL is the perimeter lenght of the object
+     Npx the number of pixels within the object
+     Pbbox the perimeter of the bounding box of the object.
+     */
+     
+     /* here we calculate a shape index for the new object to be created
+      * the radiometric index ranges from 0 to 1, 0 = identical
+      * with the shape index we want to favour compact and smooth opjects
+      * thus the shape index should range from 0 to 1,
+      * 0 = maximum compactness and smoothness */
+
+    double smooth, compact;
+    int pl, pbbox, count;
+    double bboxdiag;
+    int pl1, pl2, count1, count2;
+    int e1, n1, s1, w1, e2, n2, s2, w2, ns_extent, ew_extent;
+
+    pl = pl1 + pl2 - nshared;
+    
+    ns_extent = MAX(n1, n2) - MIN(s1, s2);
+    ew_extent = MAX(e1, e2) - MIN(w1, w2);
+
+    pbbox = 2 * (ns_extent + ew_extent);
+
+    
+    /* Smoothness Hsmooth = PL / Pbbox
+     * the smallest possible value would be 
+     * the diagonal divided by the bbox perimeter
+     * invert such that the largest possible value would be
+     * the bbox perimeter divided by the diagonal
+     */
+
+    bboxdiag = sqrt(ns_extent * ns_extent + ew_extent * ew_extent);
+    smooth = 1. - (double)bboxdiag / pl; /* smaller -> smoother */
+    
+    count = count1 + count2;
+    
+    /* compactness Hcompact = PL / sqrt(Npx)
+     * a circle is the most compact form
+     * Npx = M_PI * r * r;
+     * r = sqrt(Npx / M_pi) 
+     * pl_circle = 2 * sqrt(count * M_PI);
+     * compact = 1 - pl_circle / (pl * sqrt(count);
+     */
+    /* compact = 1 - 2 * sqrt(M_PI) / pl; */
+
+    /* PL max = Npx */
+    /* Hcompact max = Npx / sqrt(Npx) = sqrt(Npx)
+     * Hcompact / Hcompact max = (PL / sqrt(Npx)) / sqrt(Npx)
+     *                         = PL / Npx
+     */
+    compact = (double)pl / count; /* smaller -> more compact */
+
+    return globals->smooth_weight * smooth + (1 - globals->smooth_weight) * compact;
+}
+
+
 static int search_neighbors(struct ngbr_stats *Ri,
 			    struct reg_stats *Ri_rs,
                             struct NB_TREE *Ri_ngbrs, 
@@ -859,7 +974,7 @@
     double tempsim, *dp;
     struct NB_TRAV travngbr;
     struct ngbr_stats *next;
-    int cmp, i;
+    int cmp;
 
     G_debug(4, "search_neighbors");
 

Modified: grass-addons/grass7/imagery/i.segment.xl/i.segment.xl.html
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/i.segment.xl.html	2013-02-11 09:54:47 UTC (rev 54996)
+++ grass-addons/grass7/imagery/i.segment.xl/i.segment.xl.html	2013-02-11 09:55:03 UTC (rev 54997)
@@ -1,89 +1,168 @@
 <h2>DESCRIPTION</h2>
-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 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 only the spectral similarity.)  The results can be useful on their own, or used as a preprocessing step for image classification.  The segmentation preprocessing step can reduce noise and speed up the classification.
+Image segmentation or object recognition is the process of grouping 
+similar pixels into unique segments, also refered to as objects. 
+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 given a unique ID and 
+is a collection of contiguous pixels meeting some criteria. Note the 
+contrast with image classification where all pixels similar to each 
+other are assigned to the same class and do not need to be contiguous. 
+The image segmentation results can be useful on their own, or used as a 
+preprocessing step for image classification. The segmentation 
+preprocessing step can reduce noise and speed up the classification.
 
 <H2>NOTES</h2>
 
 <h3>Region Growing and Merging</h3>
-This segmentation algorithm sequentially examines all current segments in the map.  The similarity between the current segment and each of its neighbors is calculated according to the given distance formula.  Segments will be merged if they meet a number of criteria, including: 1.  The pair is mutually most similar to each other (the similarity distance will be smaller then all other neighbors), and 2. The similarity must be lower then the input threshold.  All segments are checked once per pass.  The process is repeated until no merges are made during a complete pass.
+This segmentation algorithm sequentially examines all current 
+segments in the map. The similarity between the current segment and 
+each of its neighbors is calculated according to the given distance 
+formula. Segments will be merged if they meet a number of criteria, 
+including: 1. The pair is mutually most similar to each other (the 
+similarity distance will be smaller than to any other neighbor), and 
+2. The similarity must be lower than the input threshold. The process 
+is repeated until no merges are 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 objects is used to 
+determine which objects are merged. Smaller distance values indicate 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.  During the final pass, however, if a minimum segment size of 2 or larger is given with the <em>minsize</em> parameter, segments with a smaller pixel count will be merged with their most similar neighbor even if the similarity is greater then the threshold.
-</p><p>
-Unless the <em>-w</em> flag for weighted data is used, the threshold should be set by the user between 0 and 1.0. A threshold of 0 would allow only identical valued pixels to be merged, while a threshold of 1 would allow everything to be merged.
-</p><p>
-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>
+During normal processing, merges are only allowed when the 
+similarity between two segments is lower than the givem 
+threshold value. During the final pass, however, if a minimum 
+segment size of 2 or larger is given with the <em>minsize</em> 
+parameter, segments with a smaller pixel count will be merged with 
+their most similar neighbor even if the similarity is greater than 
+the threshold.
+<p>
+The threshold should be set by the user between 0 and 1.0. A threshold 
+of 0 would allow only identical valued pixels to be merged, while a 
+threshold of 1 would allow everything to be merged. Initial empirical 
+tests indicate threshold values of 0.01 to 0.05 are reasonable values 
+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.
+
+In future , the distance calculation will also take into account the 
+shape characteristics of the segments. The normal distances are then
+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 segments (results of previous segmentations or classifications).  The different approaches are automatically detected by the program: any pixels that have identical seed values and are contiguous will be lumped into a single segment ID.
+The seeds map can be used to provide either seed pixels (random or 
+selected points from which to start the segmentation process) or 
+seed segments (results of previous segmentations or 
+classifications). The different approaches are automatically 
+detected by the program: any pixels that have identical seed values 
+and are contiguous will be assigned a unique segment ID.
 <p>
-It is expected that the <em>minsize</em> will be set to 1 if a seed map is used, but the program will allow other values to be used.  If both options are used, the final iteration that ignores the threshold also will ignore the seed map and force merges for all pixels (not just segments that have grown/merged from the seeds).
-</p>
+It is expected that the <em>minsize</em> will be set to 1 if a seed 
+map is used, but the program will allow other values to be used.  If 
+both options are used, the final iteration that ignores the 
+threshold also will ignore the seed map and force merges for all 
+pixels (not just segments that have grown/merged from the seeds).
+
 <h3>Maximum number of starting segments</h3>
-For the region growing algorithm without starting seeds, each pixel is sequentially numbered.  The current limit with CELL storage is 2 billion starting segment ID's.  If the initial map has a larger number of non-null pixels, there are two workarounds:
+For the region growing algorithm without starting seeds, each pixel 
+is sequentially numbered.  The current limit with CELL storage is 2 
+billion starting segment IDs.  If the initial map has a larger 
+number of non-null pixels, there are two workarounds:
 <p>
-1.  Use starting seed pixels.  (Maximum 2 billion pixels can be labeled with positive integers.)
-</p><p>
-2.  Use starting seed segments.  (By initial classification or other methods.)
-</p>
+1.  Use starting seed pixels. (Maximum 2 billion pixels can be 
+labeled with positive integers.)
+<p>
+2.  Use starting seed segments. (By initial classification or other 
+methods.)
+
 <h3>Boundary Constraints</h3>
-Boundary constraints limit the adjacency of pixels and segments.  Each unique value present in the <em>bounds</em> raster are considered as a MASK.  Thus no segments in the final segmentated map will cross a boundary, even if their spectral data is very similar.
+Boundary constraints limit the adjacency of pixels and segments.  
+Each unique value present in the <em>bounds</em> raster are 
+considered as a MASK. Thus no segments in the final segmentated map 
+will cross a boundary, even if their spectral data is very similar.
 
 <h3>Minimum Segment Size</h3>
-To reduce the salt and pepper affect, a <em>minsize</em> greater than 1 will add one additional pass to the processing.  During the final pass, the threshold is ignored for any segments smaller then the set size, thus forcing very small segments to merge with their most similar neighbor.
+To reduce the salt and pepper affect, a <em>minsize</em> greater 
+than 1 will add one additional pass to the processing. During the 
+final pass, the threshold is ignored for any segments smaller then 
+the set size, thus forcing very small segments to merge with their 
+most similar neighbor.
 
-<H2>EXAMPLES</h2>
-This example uses the ortho photograph included in the NC Sample Dataset.  Set up an imagery group:<br>
-<blockquote>i.group group=ortho_group input=ortho_2001_t792_1m at PERMANENT</blockquote>
+<h2>EXAMPLES</h2>
+This example uses the ortho photograph included in the NC Sample 
+Dataset.  Set up an imagery group:<br>
+<div class="code"><pre>
+i.group group=ortho_group input=ortho_2001_t792_1m at PERMANENT
+</pre></div>
 
-<p>Because the segmentation process is computationally expensive, start with a small processing area to confirm if the segmentation results meet your requirements.  Some manual adjustment of the threshold may be required. <br>
+<p>Because the segmentation process is computationally expensive, 
+start with a small processing area to confirm if the segmentation 
+results meet your requirements.  Some manual adjustment of the 
+threshold may be required. <br>
 
-<blockquote>g.region rast=ortho_2001_t792_1m at PERMANENT n=220400 s=220200 e=639000 w=638800</blockquote>
-<p></p>
+<div class="code"><pre>
+g.region rast=ortho_2001_t792_1m at PERMANENT
+</pre></div>
+
 Try out a first threshold and check the results.<br>
-<blockquote>i.segment -w -l group=ortho_group output=ortho_segs threshold=4 method=region_growing</blockquote>
+<div class="code"><pre>
+i.segment -w group=ortho_group output=ortho_segs threshold=0.04 \
+          method=region_growing 
+</pre></div>
 <p></p>
-From a visual inspection, it seems this results in oversegmentation.  Increasing the threshold: <br>
-<blockquote>i.segment -w -l --overwrite group=ortho_group output=ortho_segs threshold=10 method=region_growing</blockquote>
+From a visual inspection, it seems this results in oversegmentation.  
+Increasing the threshold: <br>
+<div class="code"><pre>
+i.segment -w group=ortho_group output=ortho_segs \
+          threshold=0.1 method=region_growing
+</pre></div>
 <p></p>
-This looks better.  There is some noise in the image, lets next force all segments smaller then 5 pixels to be merged into their most similar neighbor (even if they are less similar then required by our threshold):<br>
-<blockquote>i.segment -w -l --overwrite group=ortho_group output=ortho_segs threshold=10 method=region_growing minsize=5</blockquote>
+This looks better.  There is some noise in the image, lets next force 
+all segments smaller than 5 pixels to be merged into their most similar 
+neighbor (even if they are less similar then required by our 
+threshold):<br>
+<div class="code"><pre>
+i.segment -w --overwrite group=ortho_group output=ortho_segs \
+          threshold=0.1 method=region_growing minsize=5
+</pre></div>
 <p></p>
-Each of these segmentation steps took less then 1 minute on a decent machine.  Now that we are satisfied with the settings, we'll process the entire image:
-<blockquote>g.region rast=ortho_2001_t792_1m at PERMANENT<br>
-i.segment -w -l --overwrite group=ortho_group output=ortho_segs threshold=10 method=region_growing minsize=5 endt=5000</blockquote>
-<p></p>
-Processing the entire ortho image (over 9 million cells) took about a day.
+Processing the entire ortho image with nearly 10 million pixels took about 
+15 minutes.
 
 <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.
+<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>
-<li>Improve the optional output from this module, or better yet, add a module for i.segment.metrics.</li>
-<li>Providing updates to i.maxlik to ensure the segmentation output can be used as input for the existing classification functionality.</li>
-<li>Integration/workflow for r.fuzzy.</li>
+<li>Improve the optional output from this module, or better yet, add a 
+module for <em>i.segment.metrics</em>.</li>
+<li>Providing updates to i.maxlik to ensure the segmentation output can 
+be used as input for the existing classification functionality.</li>
+<li>Integration/workflow for <em>r.fuzzy</em>.</li>
 </ul>
 <h3>Speed</h3>
 <ul>
 <li>See create_isegs.c</li>
 </ul>
-<h3>Memory</h3>
-<ul>
-<li>User input for how much RAM can be used.</li>
-<li>Check input map type(s), currently storing in DCELL sized SEG file, could reduce this dynamically depending on input map time. (Could only reduce to FCELL, since will be storing mean we can't use CELL. Might not be worth the added code complexity.)</li>
-</ul>
-<h2>BUGS</h2>
-If the seeds map is used to give starting seed segments, the segments are renumbered starting from 1.  There is a chance a segment could be renumbered to a seed value that has not yet been processed.  If they happen to be adjacent, they would be merged.  (Possible fixes: a.  use a processing flag to make sure the pixels hasn't been previously used, or b. use negative segment ID's as a placeholder and negate all values after the seed map has been processed.)
 <H2>REFERENCES</h2>
-This project was first developed during GSoC 2012.  Project documentation, Image Segmentation references, and other information is at the <a href="http://grass.osgeo.org/wiki/GRASS_GSoC_2012_Image_Segmentation">project wiki</a>.
+This project was first developed during GSoC 2012. Project documentation, 
+Image Segmentation references, and other information is at the 
+<a href="http://grass.osgeo.org/wiki/GRASS_GSoC_2012_Image_Segmentation">project wiki</a>.
 <p>
-Information about <a href="http://grass.osgeo.org/wiki/Image_classification">classification in GRASS</a> is at available on the wiki.
+Information about 
+<a href="http://grass.osgeo.org/wiki/Image_classification">classification in GRASS</a> 
+is at available on the wiki.
 </p>
 <h2>SEE ALSO</h2>
 <em>
@@ -95,4 +174,6 @@
 </em>
 
 <h2>AUTHORS</h2>
-Eric Momsen - North Dakota State University
+Eric Momsen - North Dakota State University<br>
+Markus Metz (GSoC Mentor)
+

Modified: grass-addons/grass7/imagery/i.segment.xl/iseg.h
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/iseg.h	2013-02-11 09:54:47 UTC (rev 54996)
+++ grass-addons/grass7/imagery/i.segment.xl/iseg.h	2013-02-11 09:55:03 UTC (rev 54997)
@@ -17,6 +17,9 @@
 #include "regtree.h"
 #include "ngbrtree.h"
 
+/* #def _OR_SHAPE_ */
+
+
 /* row/col list */
 struct rc
 {
@@ -39,9 +42,13 @@
 				 * 1 if the scaling should be skipped */
     int method;			/* Segmentation method */
     int nn;			/* number of neighbors, 4 or 8 */
-    double threshold;		/* similarity threshold */
-    double alpha;
+    double max_diff;		/* max possible difference */
+    double alpha;		/* similarity threshold */
     int min_segment_size;	/* smallest number of pixels/cells allowed in a final segment */
+
+    double radio_weight;	/* weighing factor radiometric - shape */
+    double smooth_weight;       /* weighing factor smoothness - compactness */
+
     int end_t;			/* maximum number of iterations */
     int mb;
 
@@ -103,6 +110,11 @@
 double calculate_euclidean_similarity(struct ngbr_stats *, 
                                       struct ngbr_stats *,
 				      struct globals *);
+double calculate_manhattan_similarity(struct ngbr_stats *, 
+                                      struct ngbr_stats *,
+				      struct globals *);
+double calculate_shape(struct reg_stats *, struct reg_stats *,
+                       int, struct globals *);
 int fetch_reg_stats(int , int , struct reg_stats *, 
                            struct globals *);
 

Modified: grass-addons/grass7/imagery/i.segment.xl/open_files.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/open_files.c	2013-02-11 09:54:47 UTC (rev 54996)
+++ grass-addons/grass7/imagery/i.segment.xl/open_files.c	2013-02-11 09:55:03 UTC (rev 54997)
@@ -57,7 +57,7 @@
 
     /* Get min/max values of each input raster for scaling */
 
-    globals->threshold = 0.;
+    globals->max_diff = 0.;
 
     for (n = 0; n < Ref.nfiles; n++) {
 	/* returns -1 on error, 2 on empty range, quiting either way. */
@@ -71,12 +71,13 @@
 	
     }
     if (globals->weighted == FALSE)
-	globals->threshold = Ref.nfiles;
+	globals->max_diff = Ref.nfiles;
     else {
 	/* max difference with selected similarity method */
 	Ri.mean = max;
 	Rk.mean = min;
-	globals->threshold = (*globals->calculate_similarity) (&Ri, &Rk, globals);
+	globals->max_diff = 1;
+	globals->max_diff = (*globals->calculate_similarity) (&Ri, &Rk, globals);
     }
 
     /* ********** find out file segmentation size ************ */

Modified: grass-addons/grass7/imagery/i.segment.xl/parse_args.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/parse_args.c	2013-02-11 09:54:47 UTC (rev 54996)
+++ grass-addons/grass7/imagery/i.segment.xl/parse_args.c	2013-02-11 09:55:03 UTC (rev 54997)
@@ -10,54 +10,78 @@
 int parse_args(int argc, char *argv[], struct globals *globals)
 {
     struct Option *group, *seeds, *bounds, *output,
-                  *method, *threshold, *min_segment_size, *mem;
+                  *method, *similarity, *threshold, *min_segment_size,
+#ifdef _OR_SHAPE_
+		  *shape_weight, *smooth_weight,
+#endif
+		   *mem;
     struct Flag *diagonal, *weighted;
     struct Option *outband, *endt;
 
     /* required parameters */
-    /* TODO: OK to require the user to create a group?
-     * Otherwise later add an either/or option to give just a single raster map... */
     group = G_define_standard_option(G_OPT_I_GROUP);
 
     output = G_define_standard_option(G_OPT_R_OUTPUT);
 
-    /*TODO polish: any way to recommend a threshold to the user */
     threshold = G_define_option();
     threshold->key = "threshold";
     threshold->type = TYPE_DOUBLE;
     threshold->required = YES;
     threshold->label = _("Difference threshold between 0 and 1.");
-    threshold->description = _("Threshold = 0 would merge only identical segments; threshold = 1 would merge all.");
+    threshold->description = _("Threshold = 0 merges only identical segments; threshold = 1 merges all.");
 
+    /* optional parameters */
+
     method = G_define_option();
     method->key = "method";
     method->type = TYPE_STRING;
-    method->required = YES;
+    method->required = NO;
     method->answer = "region_growing";
     method->description = _("Segmentation method.");
+    method->guisection = _("Settings");
 
+    similarity = G_define_option();
+    similarity->key = "similarity";
+    similarity->type = TYPE_STRING;
+    similarity->required = YES;
+    similarity->answer = "euclidean";
+    similarity->options = "euclidean, manhattan";
+    similarity->description = _("Similarity calculation method.");
+    similarity->guisection = _("Settings");
+
     min_segment_size = G_define_option();
     min_segment_size->key = "min";
     min_segment_size->type = TYPE_INTEGER;
-    min_segment_size->required = YES;
+    min_segment_size->required = NO;
     min_segment_size->answer = "1";
     min_segment_size->options = "1-100000";
     min_segment_size->label = _("Minimum number of cells in a segment.");
     min_segment_size->description =
 	_("The final step will merge small segments with their best neighbor.");
+    min_segment_size->guisection = _("Settings");
 
-    /* optional parameters */
+#ifdef _OR_SHAPE_
+    radio_weight = G_define_option();
+    radio_weight->key = "radio_weight";
+    radio_weight->type = TYPE_DOUBLE;
+    radio_weight->required = YES;
+    radio_weight->answer = "1";
+    radio_weight->options = "0-1";
+    radio_weight->label =
+	_("Importance of radiometric (input raster) values relative to shape.");
+    radio_weight->guisection = _("Settings");
 
-    diagonal = G_define_flag();
-    diagonal->key = 'd';
-    diagonal->description =
-	_("Use 8 neighbors (3x3 neighborhood) instead of the default 4 neighbors for each pixel.");
+    smooth_weight = G_define_option();
+    smooth_weight->key = "smooth_weight";
+    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.");
+    smooth_weight->guisection = _("Settings");
+#endif
 
-    weighted = G_define_flag();
-    weighted->key = 'w';
-    weighted->description =
-	_("Weighted input, don't perform the default scaling of input maps.");
-
     /* Using raster for seeds
      * Low priority TODO: allow vector points/centroids seed input. */
     seeds = G_define_standard_option(G_OPT_R_INPUT);
@@ -94,8 +118,18 @@
     outband->key = "goodness";
     outband->required = NO;
     outband->description =
-	_("debug - save band mean, currently implemented for only 1 band.");
+	_("Goodness of fit estimate.");
 
+    diagonal = G_define_flag();
+    diagonal->key = 'd';
+    diagonal->description =
+	_("Use 8 neighbors (3x3 neighborhood) instead of the default 4 neighbors for each pixel.");
+
+    weighted = G_define_flag();
+    weighted->key = 'w';
+    weighted->description =
+	_("Weighted input, don't perform the default scaling of input maps.");
+
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
@@ -115,13 +149,38 @@
 	G_fatal_error(_("threshold should be >= 0 and <= 1"));
 
     /* segmentation methods:  1 = region growing */
-    if (strncmp(method->answer, "region_growing", 10) == 0)
+    if (strcmp(method->answer, "region_growing") == 0)
 	globals->method = 1;
     else
 	G_fatal_error("Couldn't assign segmentation method.");
 
     G_debug(1, "segmentation method: %d", globals->method);
 
+    /* distance methods for similarity measurement */
+    if (strcmp(similarity->answer, "euclidean") == 0)
+	globals->calculate_similarity = calculate_euclidean_similarity;
+    else if (strcmp(similarity->answer, "manhattan") == 0)
+	globals->calculate_similarity = calculate_manhattan_similarity;
+    else
+	G_fatal_error(_("Invalid similarity method."));
+
+#ifdef _OR_SHAPE_
+    /* consider shape */
+    globals->radio_weight = atof(radio_weight->answer);
+    if (globals->radio_weight <= 0)
+	G_fatal_error(_("Option '%s' must be > 0"), radio_weight->key);
+    if (globals->radio_weight > 1)
+	G_fatal_error(_("Option '%s' must be <= 1"), radio_weight->key);
+    globals->smooth_weight = atof(smooth_weight->answer);
+    if (globals->smooth_weight < 0)
+	G_fatal_error(_("Option '%s' must be >= 0"), smooth_weight->key);
+    if (globals->smooth_weight > 1)
+	G_fatal_error(_("Option '%s' must be <= 1"), smooth_weight->key);
+#else
+    globals->radio_weight = 1;
+    globals->smooth_weight = 0.5;
+#endif
+
     globals->min_segment_size = atoi(min_segment_size->answer);
 
     if (diagonal->answer == FALSE) {
@@ -139,9 +198,6 @@
      * selected/1 if scaling should be skipped. */
     globals->weighted = weighted->answer;
 
-    /* TODO add user input for this */
-    globals->calculate_similarity = &calculate_euclidean_similarity;
-
     globals->seeds = seeds->answer;
     if (globals->seeds) {
 	if (G_find_raster(globals->seeds, "") == NULL) {

Modified: grass-addons/grass7/imagery/i.segment.xl/write_output.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/write_output.c	2013-02-11 09:54:47 UTC (rev 54996)
+++ grass-addons/grass7/imagery/i.segment.xl/write_output.c	2013-02-11 09:55:03 UTC (rev 54997)
@@ -77,8 +77,8 @@
 	 * (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);
+	thresh = globals->alpha * globals->alpha * globals->max_diff;
+	maxdev = globals->max_diff * (1 - globals->alpha * globals->alpha);
 	mingood = 1;
 
 	G_message(_("Writing out goodness of fit"));
@@ -126,7 +126,7 @@
 			    }
 			}
 			else {
-			    sim = 1 - sim / globals->threshold;
+			    sim = 1 - sim / globals->max_diff;
 			    meanbuf[col] = sim;
 			    if (mingood > sim)
 				mingood = sim;



More information about the grass-commit mailing list