i.segment: more refactoring

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

imagery/i.segment/create_isegs.c

@@ -23,7 +23,7 @@ int create_isegs(struct globals *globals)
 
     if (globals->bounds_map == NULL) {
 	/* just one time through loop */
-	successflag = globals->method(globals);
+	successflag = globals->method_fn(globals);
     }
     else {
 	/* outer processing loop for polygon constraints */
@@ -70,7 +70,7 @@ int create_isegs(struct globals *globals)
 	    globals->col_max++;
 
 	    if (have_bound)
-		successflag = globals->method(globals);
+		successflag = globals->method_fn(globals);
 	}    /* end outer loop for processing polygons */
 
 	/* restore NULL flag */
@@ -82,8 +82,114 @@ int create_isegs(struct globals *globals)
 		    FLAG_SET(globals->null_flag, row, col);
 	    }
 	}
-
     }
 
     return successflag;
 }
+
+void find_four_neighbors(int p_row, int p_col,
+			        int neighbors[8][2])
+{
+    /* north */
+    neighbors[0][0] = p_row - 1;
+    neighbors[0][1] = p_col;
+
+    /* east */
+    neighbors[1][0] = p_row;
+    neighbors[1][1] = p_col + 1;
+
+    /* south */
+    neighbors[2][0] = p_row + 1;
+    neighbors[2][1] = p_col;
+
+    /* west */
+    neighbors[3][0] = p_row;
+    neighbors[3][1] = p_col - 1;
+
+    return;
+}
+
+void find_eight_neighbors(int p_row, int p_col,
+			         int neighbors[8][2])
+{
+    /* get the 4 orthogonal neighbors */
+    find_four_neighbors(p_row, p_col, neighbors);
+
+    /* get the 4 diagonal neighbors */
+    /* north-west */
+    neighbors[4][0] = p_row - 1;
+    neighbors[4][1] = p_col - 1;
+
+    /* north-east */
+    neighbors[5][0] = p_row - 1;
+    neighbors[5][1] = p_col + 1;
+
+    /* south-west */
+    neighbors[6][0] = p_row + 1;
+    neighbors[6][1] = p_col - 1;
+
+    /* south-east */
+    neighbors[7][0] = p_row + 1;
+    neighbors[7][1] = p_col + 1;
+
+    return;
+}
+
+/* similarity / distance between two points based on their input raster values */
+/* assumes first point values already saved in files->bands_seg - only run Segment_get once for that value... */
+/* TODO: Segment_get already happened for a[] values in the main function.  Could remove a[] from these parameters */
+double calculate_euclidean_similarity(struct ngbr_stats *Ri,
+                                      struct ngbr_stats *Rk,
+				      struct globals *globals)
+{
+    double val = 0., diff;
+    int n = globals->nbands - 1;
+
+    /* squared euclidean distance, sum the square differences for each dimension */
+    do {
+	diff = Ri->mean[n] - Rk->mean[n];
+	    
+	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 manhattan distance, sum the 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;
+}

imagery/i.segment/iseg.h

@@ -19,6 +19,10 @@
 
 /* #def _OR_SHAPE_ */
 
+/* methods */
+#define ORM_RG 1	/* region growing */
+#define ORM_MS 1	/* mean shift */
+#define ORM_WS 1	/* watershed */
 
 /* row/col list */
 struct rc
@@ -41,7 +45,8 @@ struct globals
     char *image_group;
     int weighted;		/* 0 if false/not selected, so we should scale input.
 				 * 1 if the scaling should be skipped */
-    int (*method)();		/* Segmentation method function */
+    int method;			/* Segmentation method code */
+    int (*method_fn)();		/* Segmentation method function */
     int nn;			/* number of neighbors, 4 or 8 */
     double max_diff;		/* max possible difference */
     double alpha;		/* similarity threshold */

imagery/i.segment/mean_shift.c

@@ -13,10 +13,102 @@
 #include <grass/rbtree.h>	/* Red Black Tree library functions */
 #include "iseg.h"
 
+/* standard gauss funtion:
+ * a * exp(-(x - m)^2 / (2 * stddev^2)
+ * a is not needed because the sum of weights is calculated for each 
+ * sampling window
+ * (x - m)^2 is the squared difference = diff2
+ * stddev^2 is the variance
+ * this code can be further simplified, e.g. by supplying 2 * var instead 
+ * of var
+ * 
+ * the standard deviation is the bandwidth
+ * */
+ 
+static double gauss_kernel(double diff2, double var)
+{
+    return exp(-diff2 / (2 * var));
+}
+
+
 int mean_shift(struct globals *globals)
 {
+    int row, col, t;
+    int n_changes;
+    double alpha2;
+    struct ngbr_stats Rin, Rout;
+    double diff2;
+
     G_fatal_error(_("Mean shift is not yet implemented"));
-    
     return FALSE;
+
+
+    Rin.mean = G_malloc(globals->datasize);
+    Rout.mean = G_malloc(globals->datasize);
+
+    /* TODO: need another segment structure holding output
+     * for mean shift, output of the previous iteration becomes input of the next iteration
+     * initially, input and output are original band values */
+
+    alpha2 = globals->alpha * globals->alpha;
+    
+    t = 0;
+    n_changes = 1;
+    while (t < globals->end_t && n_changes > 0) {
+
+	G_message(_("Processing pass %d..."), ++t);
+
+	n_changes = 0;
+	globals->candidate_count = 0;
+	flag_clear_all(globals->candidate_flag);
+
+	/* Set candidate flag to true/1 for all non-NULL cells */
+	for (row = globals->row_min; row < globals->row_max; row++) {
+	    for (col = globals->col_min; col < globals->col_max; col++) {
+		if (!(FLAG_GET(globals->null_flag, row, col))) {
+		    
+		    FLAG_SET(globals->candidate_flag, row, col);
+		    globals->candidate_count++;
+		}
+	    }
+	}
+
+	G_debug(4, "Starting to process %ld candidate cells",
+		globals->candidate_count);
+
+	/*process candidate cells */
+	G_percent_reset();
+	for (row = globals->row_min; row < globals->row_max; row++) {
+	    G_percent(row - globals->row_min,
+	              globals->row_max - globals->row_min, 4);
+	    for (col = globals->col_min; col < globals->col_max; col++) {
+		if (!(FLAG_GET(globals->candidate_flag, row, col)))
+		    continue;
+
+		/* get current band values */
+		Segment_get(&globals->bands_seg, (void *)Rin.mean,
+			    row, col);
+
+		/* adapt initial spatial and range bandwiths */
+
+		/* calculate new band values */
+		
+		/* if the squared difference between old and new band values 
+		 * is larger than alpha2, then increase n_changes */
+		
+		diff2 = (globals->calculate_similarity)(&Rin, &Rout, globals);
+		if (diff2 > alpha2)
+		    n_changes++;
+	    }
+	}
+    }
+    if (n_changes > 1)
+	G_message(_("Mean shift stopped at %d due to reaching max iteration limit, more changes may be possible"), t);
+    else
+	G_message(_("Mean shift converged after %d iterations"), t);
+    
+    /* identify connected components */
+    
+    return TRUE;
 }
 

imagery/i.segment/parse_args.c

@@ -152,13 +152,19 @@ int parse_args(int argc, char *argv[], struct globals *globals)
     if (globals->alpha <= 0 || globals->alpha >= 1)
 	G_fatal_error(_("Threshold should be > 0 and < 1"));
 
-    /* segmentation methods:  1 = region growing */
-    if (strcmp(method->answer, "region_growing") == 0)
-	globals->method = region_growing;
-    else if (strcmp(method->answer, "mean_shift") == 0)
-	globals->method = mean_shift;
-    else if (strcmp(method->answer, "watershed") == 0)
-	globals->method = watershed;
+    /* segmentation methods */
+    if (strcmp(method->answer, "region_growing") == 0) {
+	globals->method = ORM_RG;
+	globals->method_fn = region_growing;
+    }
+    else if (strcmp(method->answer, "mean_shift") == 0) {
+	globals->method = ORM_MS;
+	globals->method_fn = mean_shift;
+    }
+    else if (strcmp(method->answer, "watershed") == 0) {
+	globals->method = ORM_WS;
+	globals->method_fn = watershed;
+    }
     else
 	G_fatal_error(_("Unable to assign segmentation method"));
 

imagery/i.segment/region_growing.c

@@ -742,113 +742,6 @@ static int find_best_neighbor(struct ngbr_stats *Ri,
     return n_ngbrs;
 }
 
-void find_four_neighbors(int p_row, int p_col,
-			        int neighbors[8][2])
-{
-    /* north */
-    neighbors[0][0] = p_row - 1;
-    neighbors[0][1] = p_col;
-
-    /* east */
-    neighbors[1][0] = p_row;
-    neighbors[1][1] = p_col + 1;
-
-    /* south */
-    neighbors[2][0] = p_row + 1;
-    neighbors[2][1] = p_col;
-
-    /* west */
-    neighbors[3][0] = p_row;
-    neighbors[3][1] = p_col - 1;
-
-    return;
-}
-
-void find_eight_neighbors(int p_row, int p_col,
-			         int neighbors[8][2])
-{
-    /* get the 4 orthogonal neighbors */
-    find_four_neighbors(p_row, p_col, neighbors);
-
-    /* get the 4 diagonal neighbors */
-    /* north-west */
-    neighbors[4][0] = p_row - 1;
-    neighbors[4][1] = p_col - 1;
-
-    /* north-east */
-    neighbors[5][0] = p_row - 1;
-    neighbors[5][1] = p_col + 1;
-
-    /* south-west */
-    neighbors[6][0] = p_row + 1;
-    neighbors[6][1] = p_col - 1;
-
-    /* south-east */
-    neighbors[7][0] = p_row + 1;
-    neighbors[7][1] = p_col + 1;
-
-    return;
-}
-
-/* similarity / distance between two points based on their input raster values */
-/* assumes first point values already saved in files->bands_seg - only run Segment_get once for that value... */
-/* TODO: Segment_get already happened for a[] values in the main function.  Could remove a[] from these parameters */
-double calculate_euclidean_similarity(struct ngbr_stats *Ri,
-                                      struct ngbr_stats *Rk,
-				      struct globals *globals)
-{
-    double val = 0., diff;
-    int n = globals->nbands - 1;
-
-    /* squared euclidean distance, sum the square differences for each dimension */
-    do {
-	diff = Ri->mean[n] - Rk->mean[n];
-	    
-	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 manhattan distance, sum the 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)
 {