[GRASS-SVN] r52939 - grass-addons/grass7/imagery/i.segment.xl
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Aug 27 12:08:41 PDT 2012
Author: mmetz
Date: 2012-08-27 12:08:40 -0700 (Mon, 27 Aug 2012)
New Revision: 52939
Added:
grass-addons/grass7/imagery/i.segment.xl/i.segment.xl.html
grass-addons/grass7/imagery/i.segment.xl/ngbrtree.c
grass-addons/grass7/imagery/i.segment.xl/ngbrtree.h
grass-addons/grass7/imagery/i.segment.xl/rclist.c
grass-addons/grass7/imagery/i.segment.xl/regtree.c
grass-addons/grass7/imagery/i.segment.xl/regtree.h
Removed:
grass-addons/grass7/imagery/i.segment.xl/i.segment.html
Modified:
grass-addons/grass7/imagery/i.segment.xl/Makefile
grass-addons/grass7/imagery/i.segment.xl/create_isegs.c
grass-addons/grass7/imagery/i.segment.xl/flag.c
grass-addons/grass7/imagery/i.segment.xl/flag.h
grass-addons/grass7/imagery/i.segment.xl/iseg.h
grass-addons/grass7/imagery/i.segment.xl/main.c
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:
xl version of i.segment
Modified: grass-addons/grass7/imagery/i.segment.xl/Makefile
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/Makefile 2012-08-27 18:43:23 UTC (rev 52938)
+++ grass-addons/grass7/imagery/i.segment.xl/Makefile 2012-08-27 19:08:40 UTC (rev 52939)
@@ -1,9 +1,9 @@
MODULE_TOPDIR = ../..
-PGM = i.segment
+PGM = i.segment4
-LIBES = $(IMAGERYLIB) $(RASTERLIB) $(SEGMENTLIB) $(GISLIB) $(LINKMLIB) $(BTREE2LIB)
-DEPENDENCIES = $(IMAGERYDEP) $(RASTERDEP) $(SEGMENTDEP) $(GISDEP) $(LINKMDEP)
+LIBES = $(IMAGERYLIB) $(RASTERLIB) $(SEGMENTLIB) $(GISLIB) $(BTREE2LIB)
+DEPENDENCIES = $(IMAGERYDEP) $(RASTERDEP) $(SEGMENTDEP) $(GISDEP)
include $(MODULE_TOPDIR)/include/Make/Module.make
Modified: grass-addons/grass7/imagery/i.segment.xl/create_isegs.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/create_isegs.c 2012-08-27 18:43:23 UTC (rev 52938)
+++ grass-addons/grass7/imagery/i.segment.xl/create_isegs.c 2012-08-27 19:08:40 UTC (rev 52939)
@@ -3,1159 +3,1123 @@
/* Currently only region growing is implemented */
#include <stdlib.h>
-#include <float.h> /* for DBL_MAX */
-#include <math.h> /* for fabs() */
+#include <string.h>
+#include <float.h>
+#include <math.h>
+#include <time.h>
#include <grass/gis.h>
#include <grass/glocale.h>
#include <grass/raster.h>
#include <grass/segment.h> /* segmentation library */
-#include <grass/linkm.h> /* memory manager for linked lists */
#include <grass/rbtree.h> /* Red Black Tree library functions */
#include "iseg.h"
-#ifdef PROFILE
-#include <time.h>
-#endif
+#define EPSILON 1.0e-9
-/* 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. */
+int debug_level;
-int create_isegs(struct files *files, struct functions *functions)
+/* internal functions */
+static int merge_regions(struct ngbr_stats *, /* Ri */
+ struct reg_stats *,
+ struct ngbr_stats *, /* Rk */
+ struct reg_stats *,
+ struct globals *);
+static int search_neighbors(struct ngbr_stats *, /* Ri */
+ struct NB_TREE *, /* Ri's neighbors */
+ double *, /* highest similarity */
+ struct ngbr_stats *, /* Ri's best neighbor */
+ struct globals *);
+static int set_candidate_flag(struct ngbr_stats *, int , struct globals *);
+static int find_best_neighbor(struct ngbr_stats *, struct NB_TREE *,
+ struct reg_stats **, struct ngbr_stats *,
+ double *, int, struct globals *);
+
+/* function used by binary tree to compare items */
+static int compare_rc(const void *first, const void *second)
{
- int lower_bound, upper_bound, row, col;
- int successflag = 1;
- struct Range range;
+ struct rc *a = (struct rc *)first, *b = (struct rc *)second;
- /* Modify the threshold for easier similarity comparisons.
- * For Euclidean, square the threshold so we don't need to calculate the root value in the distance calculation.
- * In either case, multiplie by the number of input bands, so the same threshold will achieve similar thresholds even for different numbers of input bands.*/
- if (functions->calculate_similarity == calculate_euclidean_similarity)
- functions->threshold =
- functions->threshold * functions->threshold * files->nbands;
- else
- functions->threshold = functions->threshold * files->nbands;
+ if (a->row < b->row)
+ return -1;
+ if (a->row > b->row)
+ return 1;
+ /* same row */
+ if (a->col < b->col)
+ return -1;
+ if (a->col > b->col)
+ return 1;
+ /* same row and col */
+ return 0;
+}
- /* set parameters for outer processing loop for polygon constraints */
- if (files->bounds_map == NULL) { /*normal processing */
- lower_bound = upper_bound = 0; /* so run the segmentation algorithm just one time */
- }
- else {
- if (Rast_read_range(files->bounds_map, files->bounds_mapset, &range) != 1) { /* returns -1 on error, 2 on empty range, quiting either way. */
- G_fatal_error(_("No min/max found in boundary raster map <%s>"),
- files->bounds_map);
- }
- Rast_get_range_min_max(&range, &lower_bound, &upper_bound);
- /* speed enhancement, when processing with bounds: get the unique values.
- * As is, we will iterate at one time over the entire raster for each integer between the upper and lower bound, even if no regions exist with that value. */
- }
+static int compare_ints(const void *first, const void *second)
+{
+ int a = *(int *)first, b = *(int *)second;
- /* processing loop for polygon/boundary constraints */
- if (files->bounds_map != NULL)
- G_message(_("Starting image segmentation within boundary constraints, the percent complete is based the range of values in the boundary constraints map"));
- for (files->current_bound = lower_bound;
- files->current_bound <= upper_bound; files->current_bound++) {
+ return (a < b ? -1 : (a > b));
+}
- if (files->bounds_map != NULL)
- G_percent(files->current_bound - lower_bound,
- upper_bound - lower_bound, 1);
+static int compare_double(double first, double second)
+{
+ if (first > second)
+ return (first > (second * (1 + EPSILON)));
+ if (first < second)
+ return (-1 * (second > (first * (1 + EPSILON))));
+ return 0;
+}
- /* *** check the processing window *** */
+int create_isegs(struct globals *globals)
+{
+ int row, col;
+ int successflag = 1;
+ int have_bound;
+ CELL current_bound, bounds_val;
- /* set boundaries at "opposite" end, change until reach lowest/highest */
- files->minrow = files->nrows;
- files->mincol = files->ncols;
- files->maxrow = files->maxcol = 0;
+ G_debug(1, "Threshold: %g", globals->threshold);
+ G_debug(1, "segmentation method: %d", globals->method);
- if (files->bounds_map == NULL) {
- /* check the NULL flag to see where the first/last row/col of real data are, and reduce the processing window.
- * This could help (a little?) if a MASK is used that removes a large border portion of the map. */
- for (row = 0; row < files->nrows; row++) {
- for (col = 0; col < files->ncols; col++) {
+ if (globals->bounds_map == NULL) {
+ /* just one time through loop */
+ successflag = region_growing(globals);
+ }
+ else {
+ /* outer processing loop for polygon constraints */
+ for (current_bound = globals->lower_bound;
+ current_bound <= globals->upper_bound; current_bound++) {
- if (!(FLAG_GET(files->null_flag, row, col))) {
+ G_debug(1, "current_bound = %d", current_bound);
- if (files->minrow > row)
- files->minrow = row;
- if (files->maxrow < row)
- files->maxrow = row;
- if (files->mincol > col)
- files->mincol = col;
- if (files->maxcol < col)
- files->maxcol = col;
- }
- }
- }
- }
- else {
- for (row = 0; row < files->nrows; row++) {
- for (col = 0; col < files->ncols; col++) {
+ have_bound = 0;
- segment_get(&files->bounds_seg, &files->bounds_val, row,
- col);
- if (files->bounds_val == files->current_bound &&
- !(FLAG_GET(files->orig_null_flag, row, col))) {
- FLAG_UNSET(files->null_flag, row, col);
+ /* get min/max row/col to narrow the processing window */
+ globals->row_min = globals->nrows;
+ globals->row_max = 0;
+ globals->col_min = globals->ncols;
+ globals->col_max = 0;
+ for (row = 0; row < globals->nrows; row++) {
+ for (col = 0; col < globals->ncols; col++) {
+ segment_get(&globals->bounds_seg, &bounds_val,
+ row, col);
- if (files->minrow > row)
- files->minrow = row;
- if (files->maxrow < row)
- files->maxrow = row;
- if (files->mincol > col)
- files->mincol = col;
- if (files->maxcol < col)
- files->maxcol = col;
+ if (bounds_val == current_bound) {
+ have_bound = 1;
+ FLAG_UNSET(globals->null_flag, row, col);
+
+ if (globals->row_min > row)
+ globals->row_min = row;
+ if (globals->row_max < row)
+ globals->row_max = row;
+ if (globals->col_min > col)
+ globals->col_min = col;
+ if (globals->col_max < col)
+ globals->col_max = col;
}
- else /* pixel is outside the current boundary or was null in the input bands */
- FLAG_SET(files->null_flag, row, col);
+ else
+ FLAG_SET(globals->null_flag, row, col);
}
}
+ globals->row_max++;
+ globals->col_max++;
- } /* end of else, set up for bounded segmentation */
-
- /* consider maxrow/maxcol as nrow, ncol, loops will have: row < maxrow
- * so need to increment by one */
- files->maxrow++;
- files->maxcol++;
-
- /* run the segmentation algorithm */
-
- if (functions->method == 1) {
- successflag = region_growing(files, functions);
- }
-
- /* check if something went wrong */
- if (successflag == FALSE)
- G_fatal_error(_("Error during segmentation"));
-
- } /* end outer loop for processing bounds constraints (will just run once if not provided) */
-
- /* reset null flag to the original if we have boundary constraints */
- if (files->bounds_map != NULL) {
- for (row = 0; row < files->nrows; row++) {
- for (col = 0; col < files->ncols; col++) {
- if (FLAG_GET(files->orig_null_flag, row, col))
- FLAG_SET(files->null_flag, row, col);
- else
- FLAG_UNSET(files->null_flag, row, col);
- }
- }
+ if (have_bound)
+ successflag = region_growing(globals);
+ } /* end outer loop for processing polygons */
}
-
return successflag;
}
-int region_growing(struct files *files, struct functions *functions)
+
+int region_growing(struct globals *globals)
{
int row, col, t;
- double threshold, Ri_similarity, Rk_similarity, tempsim;
- int endflag; /* =TRUE if there were no merges on that processing iteration */
- int pathflag; /* =TRUE if we didn't find mutual neighbors, and should continue with Rk */
- struct pixels *Ri_head, *Rk_head, *Rin_head, *Rkn_head,
- *current, *newpixel, *Ri_bestn;
- int Ri_count, Rk_count; /* number of pixels/cells in Ri and Rk */
+ double threshold, adjthresh, Ri_similarity, Rk_similarity;
+ double alpha2, divisor; /* threshold parameters */
+ int n_merges, do_merge; /* number of merges on that iteration */
+ int pathflag; /* =1 if we didn't find mutually best neighbors, continue with Rk */
+ int candidates_only;
+ struct ngbr_stats Ri,
+ Rk,
+ Rk_bestn, /* Rk's best neighbor */
+ *next;
+ int Ri_nn, Rk_nn; /* number of neighbors for Ri/Rk */
+ struct NB_TREE *Ri_ngbrs, *Rk_ngbrs;
+ struct NB_TRAV travngbr;
+ struct reg_stats *Ri_rs, *Rk_rs, *Rk_bestn_rs;
+ int verbose = (G_verbose() >= 3);
+ double *dp;
+ struct NB_TREE *tmpnbtree;
-#ifndef ROWMAJOR
- unsigned long int z, end_z; /* only used for z-order */
- int j;
-#endif
+ G_verbose_message("Running region growing algorithm");
+
+ Ri.mean = G_malloc(globals->datasize);
+ Rk.mean = G_malloc(globals->datasize);
+ Rk_bestn.mean = G_malloc(globals->datasize);
+
+ Ri_ngbrs = nbtree_create(globals->nbands, globals->datasize);
+ Rk_ngbrs = nbtree_create(globals->nbands, globals->datasize);
-#ifdef VCLOSE
- struct pixels *Rclose_head, *Rc_head, *Rc_tail, *Rcn_head;
- int Rc_count;
-#endif
+ t = 0;
+ n_merges = 1;
-#ifdef PROFILE
- clock_t start, end;
- clock_t merge_start, merge_end;
- double merge_accum, merge_lap;
- clock_t fn_start, fn_end;
- double fn_accum, fn_lap;
- clock_t pass_start, pass_end;
+ /* threshold calculation */
+ alpha2 = globals->alpha * globals->alpha;
+ /* make the divisor a constant ? */
+ divisor = globals->nrows + globals->ncols;
+
+ debug_level = 5;
+
+ while (t < globals->end_t && n_merges > 0) {
- merge_accum = fn_accum = 0;
- start = clock();
-#endif
- /* files->token has the "link_head" for linkm: linked list memory allocation.
- *
- * 4 linked lists of pixels:
- * Ri = current focus segment
- * Rk = Ri's most similar neighbor
- * Rkn = Rk's neighbors
- * Rin = Ri's neigbors
- * */
- if (files->bounds_map == NULL)
- G_message(_("Running region growing algorithm, the percent completed is based on %d max iterations, but the process will end earlier if no further merges can be made."),
- functions->end_t);
+ /* optional threshold as a function of t. */
+ threshold = alpha2 * globals->threshold;
- t = 1;
+ t++;
+ G_debug(3, "####### Starting outer do loop! t = %d #######", t);
+ if (verbose)
+ G_verbose_message(_("Pass %d:"), t);
+ else
+ G_percent(t, globals->end_t, 1);
- /*set linked list head pointers to null. */
- Ri_head = NULL;
- Rk_head = NULL;
- Rin_head = NULL;
- Rkn_head = NULL;
- Ri_bestn = NULL;
-#ifdef VCLOSE
- Rclose_head = NULL;
- Rc_head = NULL;
- Rcn_head = NULL;
-#endif
+ if (t >= 50)
+ debug_level = 5;
- /* One paper mentioned gradually lowering the threshold at each iteration.
- * if this is implemented, move this assignment inside the do loop and make it a function of t. */
- threshold = functions->threshold;
+ n_merges = 0;
+ globals->candidate_count = 0;
+ flag_clear_all(globals->candidate_flag);
- /* do while loop until no merges are made, or until t reaches maximum number of iterations */
- do {
-#ifdef PROFILE
- pass_start = clock();
-#endif
-#ifdef SIGNPOST
- fprintf(stdout, "pass %d\n", t);
-#endif
- G_debug(3, "####### Starting outer do loop! t = %d #######", t);
- if (files->bounds_map == NULL)
- G_percent(t, functions->end_t, 1);
+ /* 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++;
+ }
+ }
+ }
- endflag = TRUE;
+ G_debug(4, "Starting to process %d candidate cells",
+ globals->candidate_count);
- /* Set candidate flag to true/1 for all pixels */
- set_all_candidate_flags(files);
+ /*process candidate cells */
+ for (row = globals->row_min; row < globals->row_max; row++) {
+ if (verbose)
+ G_percent(row, globals->row_max, 4);
+ for (col = globals->col_min; col < globals->col_max; col++) {
+ if (!(FLAG_GET(globals->candidate_flag, row, col)))
+ continue;
- /*process candidate pixels for this iteration */
+ pathflag = TRUE;
+ candidates_only = TRUE;
- /*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 */
+ nbtree_clear(Ri_ngbrs);
+ nbtree_clear(Rk_ngbrs);
- /*largest dimension: */
- if (files->nrows > files->ncols)
- end_z = files->nrows;
- else
- end_z = files->ncols;
+ G_debug(4, "Next starting cell: row, %d, col, %d",
+ 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++;
+ /* First cell in Ri is current row/col */
+ Ri.row = row;
+ Ri.col = col;
- /*squared: */
- end_z *= end_z;
+ /* get Ri's segment ID */
+ segment_get(&globals->rid_seg, (void *)&Ri.id, Ri.row, Ri.col);
+
+ if (Ri.id < 0)
+ continue;
- 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;
+ /* find segment neighbors */
+ /* find Ri's best neighbor, clear candidate flag */
+ Ri_similarity = globals->threshold + 1;
+
+ globals->rs.id = Ri.id;
+ if ((Ri_rs = rgtree_find(globals->reg_tree, &(globals->rs))) != NULL) {
+ memcpy(Ri.mean, Ri_rs->mean, globals->datasize);
+ Ri.count = Ri_rs->count;
}
- 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
+ else {
+ segment_get(&globals->bands_seg, (void *)Ri.mean,
+ Ri.row, Ri.col);
+ Ri.count = 1;
+ }
+ /* Ri is now complete */
- if (FLAG_GET(files->candidate_flag, row, col) &&
- FLAG_GET(files->seeds_flag, row, col)) {
+ Rk_rs = NULL;
+ Ri_nn = find_best_neighbor(&Ri, Ri_ngbrs, &Rk_rs,
+ &Rk, &Ri_similarity, 1,
+ globals);
+ /* Rk is now complete */
- /*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);
-#endif
- Rk_count = 0;
+ if (Ri_nn == 0) {
+ /* this can only happen if only one segment is left */
+ G_debug(4, "Segment had no valid neighbors");
+ pathflag = FALSE;
+ Ri.count = 0;
+ }
+ if (!(t & 1) && Ri_nn == 1 &&
+ !(FLAG_GET(globals->candidate_flag, Rk.row, Rk.col)) &&
+ compare_double(Ri_similarity, threshold) == -1) {
+ /* this is slow ??? */
+ int smaller = Rk.count;
- /* 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;
+ if (Ri.count < Rk.count)
+ smaller = Ri.count;
- pathflag = TRUE;
+ adjthresh = pow(alpha2, 1. + (double) smaller / divisor) *
+ globals->threshold;
- 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);
+ if (compare_double(Ri_similarity, adjthresh) == -1) {
+ G_debug(4, "Ri nn == 1");
+ if (Rk.count < 2)
+ G_fatal_error("Rk count too low");
+ if (Rk.count < Ri.count)
+ G_debug(1, "Rk count lower than Ri count");
+ merge_regions(&Ri, Ri_rs, &Rk, Rk_rs, globals);
+ n_merges++;
+ }
+
pathflag = FALSE;
+ }
+ /* this is slow ??? */
+ if (0 && (t & 1) && Rk.count <= globals->nn)
+ candidates_only = FALSE;
+
+ while (pathflag) {
+ pathflag = FALSE;
+
+ /* optional check if Rk is candidate
+ * to prevent backwards merging */
+ if (candidates_only &&
+ !(FLAG_GET(globals->candidate_flag, Rk.row, Rk.col))) {
- /* find segment neighbors, if we don't already have them */
- if (Rin_head == NULL) {
-#ifdef PROFILE
- fn_start = clock();
-#endif
-
- 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);
-#endif
+ Ri_similarity = globals->threshold + 1;
}
- if (Rin_head != NULL) { /*found neighbors, find best neighbor then see if is mutually best neighbor */
+ candidates_only = TRUE;
- /* ******** 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 */
+ if (compare_double(Ri_similarity, threshold) == -1) {
+ do_merge = 1;
- /* for each of Ri's neighbors */
- for (current = Rin_head; current != NULL;
- current = current->next) {
- tempsim = (*functions->calculate_similarity)
- (Ri_head, current, files, functions);
+ /* we'll have the neighbor pixel to start with. */
+ G_debug(4, "Working with Rk");
-#ifdef VCLOSE
- /* 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 */
+ /* find Rk's best neighbor, do not clear candidate flag */
+ Rk_similarity = globals->threshold + 1;
+ Rk_bestn_rs = NULL;
+ Rk_nn = find_best_neighbor(&Rk, Rk_ngbrs,
+ &Rk_bestn_rs,
+ &Rk_bestn,
+ &Rk_similarity, 0,
+ globals);
- /* *** 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);
+ /* not mutually best neighbors */
+ if (Rk_similarity != Ri_similarity) {
+ do_merge = 0;
+ }
+ /* Ri has only one neighbor, merge */
+ if (Ri_nn == 1 && Rk_nn > 1)
+ do_merge = 1;
- /* 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);
+ /* adjust threshold */
+ if (do_merge) {
+ int smaller = Rk.count;
- /* Add Rc pixels to Ri */
- Rc_tail->next = Ri_head;
- Ri_head = Rc_head;
+ if (Ri.count < Rk.count)
+ smaller = Ri.count;
- /*to consider, recurse? Check all Rcn neighbors if they are very close? */
+ adjthresh = pow(alpha2, 1. + (double) smaller / divisor) *
+ globals->threshold;
- Rc_head = NULL;
- my_dispose_list(files->token, &Rcn_head);
+ if (compare_double(Ri_similarity, adjthresh) == 1)
+ do_merge = 0;
}
- 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;
- }
- }
- 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 (do_merge) {
+
+ G_debug(4, "merge neighbor trees");
- find_segment_neighbors(&Rk_head, &Rkn_head,
- &Rk_count, files,
- functions);
+ Ri_nn -= Ri_ngbrs->count;
+ Ri_nn += (Rk_nn - Rk_ngbrs->count);
+ globals->ns.id = Rk.id;
+ nbtree_remove(Ri_ngbrs, &(globals->ns));
- /* ******** 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 */
+ nbtree_init_trav(&travngbr, Rk_ngbrs);
+ while ((next = nbtree_traverse(&travngbr))) {
+ if (!nbtree_find(Ri_ngbrs, next) && next->id != Ri.id)
+ nbtree_insert(Ri_ngbrs, next);
+ }
+ nbtree_clear(Rk_ngbrs);
+ Ri_nn += Ri_ngbrs->count;
- /* 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);
+ merge_regions(&Ri, Ri_rs, &Rk, Rk_rs, globals);
+ /* Ri is now updated, Rk no longer usable */
- 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 */
+ /* made a merge, need another iteration */
+ n_merges++;
- if (Rk_similarity == Ri_similarity) { /* mutually most similar neighbors */
+ Ri_similarity = globals->threshold + 1;
+ Rk_similarity = globals->threshold + 1;
-#ifdef PROFILE
- merge_start = clock();
-#endif
- 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);
-#endif
- endflag = FALSE; /* we've made at least one merge, so want another t iteration */
- }
- else { /* they weren't mutually best neighbors */
+ /* we have checked the neighbors of Ri, Rk already
+ * use faster version of finding the best neighbor
+ */
+
+ /* use neighbor tree to find Ri's new best neighbor */
+ search_neighbors(&Ri, Ri_ngbrs, &Ri_similarity,
+ &Rk, globals);
- /* 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 (Ri_nn > 0 && compare_double(Ri_similarity, threshold) == -1) {
- 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);
- }
+ globals->rs.id = Ri.id;
+ Ri_rs = rgtree_find(globals->reg_tree, &(globals->rs));
+ globals->rs.id = Rk.id;
+ Rk_rs = rgtree_find(globals->reg_tree, &(globals->rs));
- } /* 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, 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;
+ pathflag = TRUE;
+ /* candidates_only:
+ * FALSE: less passes, takes a bit longer, but less memory
+ * TRUE: more passes, is a bit faster */
+ candidates_only = FALSE;
}
- else
- my_dispose_list(files->token, &Rin_head);
+ /* else end of Ri -> Rk chain since we merged Ri and Rk
+ * go to next row, col */
}
- else
- pathflag = FALSE;
+ else {
+ if (compare_double(Rk_similarity, threshold) == -1) {
+ pathflag = TRUE;
+ }
+ /* test this */
+ if (!(FLAG_GET(globals->candidate_flag, Rk.row, Rk.col)))
+ pathflag = FALSE;
- }
+ if (Rk_nn < 2)
+ pathflag = FALSE;
- } /*end pathflag do loop */
- } /*end if pixel is candidate and seed pixel */
- } /*next column (or next z) */
-#ifdef ROWMAJOR
- } /*next row */
-#endif
-#ifdef PROFILE
- pass_end = clock();
- fprintf(stdout, "pass %d took: %g\n", t,
- ((double)(pass_end - pass_start)) / CLOCKS_PER_SEC);
-#endif
+ if (Rk.id < 1)
+ pathflag = FALSE;
+ if (pathflag) {
+
+ /* clear candidate flag for Rk */
+ if (FLAG_GET(globals->candidate_flag, Rk.row, Rk.col)) {
+ set_candidate_flag(&Rk, FALSE, globals);
+ }
- /* 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 */
+ /* Use Rk as next Ri:
+ * this is the eCognition technique. */
+ G_debug(4, "do ecog");
+ Ri_nn = Rk_nn;
+ Ri_similarity = Rk_similarity;
- if (t == 2 && files->bounds_map == NULL)
- G_warning(_("No segments were created. Verify threshold and region settings."));
- /* future enhancement, remove the "&& bounds_map == NULL" if we check for unique bounds values. */
+ dp = Ri.mean;
+ Ri = Rk;
+ Rk = Rk_bestn;
+ Rk_bestn.mean = dp;
+
+ Ri_rs = Rk_rs;
+ Rk_rs = Rk_bestn_rs;
+ Rk_bestn_rs = NULL;
- if (endflag == FALSE)
- G_message(_("Merging processes stopped due to reaching max iteration limit, more merges may be possible"));
+ tmpnbtree = Ri_ngbrs;
+ Ri_ngbrs = Rk_ngbrs;
+ Rk_ngbrs = tmpnbtree;
+ nbtree_clear(Rk_ngbrs);
+ }
+ }
+ } /* end if < threshold */
+ } /* end pathflag */
+ } /* next col */
+ } /* next row */
+ /* finished one pass for processing candidate pixels */
+ G_verbose_message("%d merges", n_merges);
+ G_debug(4, "Finished pass %d", t);
+ }
- /* ****************************************************************************************** */
- /* speed enhancement: after < 3 (?) merges are made in one pass, switch processing modes.
- * It seems a significant portion of the time is spent merging small pixels into the largest
- * segments. So consider allowing multiple merges after finding one large Ri.
- *
- * Maybe related to this, and/or integrated into the first loops: If a segment has only one neighbor
- * go ahead and merge it if similarity is < threshold.
- *
- * ****************************************************************************************** */
+ /*end t loop *//*TODO, should there be a max t that it can iterate for? Include t in G_message? */
+ if (n_merges > 0)
+ G_message(_("Segmentation processes stopped at %d due to reaching max iteration limit, more merges may be possible"), t);
+ else
+ G_message(_("Segmentation converged after %d iterations."), t);
-
/* ****************************************************************************************** */
/* final pass, ignore threshold and force a merge for small segments with their best neighbor */
/* ****************************************************************************************** */
+
+ if (globals->min_segment_size > 1) {
+ G_message(_("Merging segments smaller than %d cells"), globals->min_segment_size);
+ /* optional threshold as a function of t. */
+ threshold = globals->alpha * globals->alpha * globals->threshold;
- if (functions->min_segment_size > 1 && t > 2) { /* NOTE: added t > 2, it doesn't make sense to force merges if no merges were made on the original pass. Something should be adjusted first */
+ flag_clear_all(globals->candidate_flag);
- if (files->bounds_map == NULL) {
- G_message
- (_("Final iteration, forcing merges for small segments, percent complete based on rows."));
+ /* 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++;
+ }
+ }
}
- /* for the final forced merge, the candidate flag is just to keep track if we have confirmed if:
- * a. the segment size is >= to the minimum allowed size or
- * b. we have merged it with its best neighbor
- */
- set_all_candidate_flags(files);
+ G_debug(4, "Starting to process %d candidate cells",
+ globals->candidate_count);
- for (row = files->minrow; row < files->maxrow; row++) {
- if (files->bounds_map == NULL)
- G_percent(row, files->maxrow - 1, 1);
- for (col = files->mincol; col < files->maxcol; col++) {
+ /* process candidate cells */
+ for (row = globals->row_min; row < globals->row_max; row++) {
+ G_percent(row, globals->row_max, 9);
+ for (col = globals->col_min; col < globals->col_max; col++) {
+ int do_merge = 1;
+
+ if (!(FLAG_GET(globals->candidate_flag, row, col)))
+ continue;
+
+ nbtree_clear(Ri_ngbrs);
+ nbtree_clear(Rk_ngbrs);
- if (FLAG_GET(files->candidate_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);
- Rk_count = 0;
+ Ri.row = row;
+ Ri.col = col;
- /* 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 = Ri_head;
- newpixel->row = row;
- newpixel->col = col;
- Ri_head = newpixel;
+ /* get segment id */
+ segment_get(&globals->rid_seg, (void *) &Ri.id, row, col);
+
+ if (Ri.id < 0)
+ continue;
- G_debug(4, "Next starting pixel: row, %d, col, %d",
- Ri_head->row, Ri_head->col);
+ while (do_merge) {
- /* find segment neighbors */
- find_segment_neighbors
- (&Ri_head, &Rin_head, &Ri_count, files, functions);
+ /* get segment size */
+ globals->rs.id = Ri.id;
+ Ri_rs = rgtree_find(globals->reg_tree, &(globals->rs));
- if (Rin_head != NULL) { /*found neighbors */
- if (Ri_count >= functions->min_segment_size) /* don't force a merge */
- set_candidate_flag(Ri_head, FALSE, files);
+ Ri.count = 1;
+ do_merge = 0;
+ if (Ri_rs != NULL) {
+ Ri.count = Ri_rs->count;
+ memcpy(Ri.mean, Ri_rs->mean, globals->datasize);
+ }
- else { /* Merge with most similar neighbor */
+ /* merge all smaller than min size */
+ if (Ri.count < globals->min_segment_size)
+ do_merge = 1;
- /* find Ri's most similar neighbor */
- Ri_bestn = NULL;
- Ri_similarity = DBL_MAX; /* set current similarity to max value */
- segment_get(&files->bands_seg, (void *)files->bands_val, Ri_head->row, Ri_head->col); /* current segment values */
+ Ri_nn = 0;
+ Ri_similarity = globals->threshold + 1;
- /* for each of Ri's neighbors */
- for (current = Rin_head; current != NULL;
- current = current->next) {
- tempsim = (*functions->calculate_similarity)
- (Ri_head, current, files, functions);
+ if (do_merge) {
- if (tempsim < Ri_similarity) {
- Ri_similarity = tempsim;
- Ri_bestn = current;
- }
- }
- if (Ri_bestn != NULL) {
+ segment_get(&globals->bands_seg, (void *)Ri.mean,
+ Ri.row, Ri.col);
- /* we'll have the neighbor pixel to start with. */
- 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;
+ /* find Ri's best neighbor, clear candidate flag */
+ Ri_nn = find_best_neighbor(&Ri, Ri_ngbrs,
+ &Rk_rs, &Rk,
+ &Ri_similarity, 1,
+ globals);
+ }
- /* get the full pixel/cell membership list for Rk *//* speed enhancement: a seperate function for this, since we don't need the neighbors */
- find_segment_neighbors(&Rk_head, &Rkn_head,
- &Rk_count, files,
- functions);
+ if (do_merge) {
- merge_values(Ri_head, Rk_head, Ri_count,
- Rk_count, files);
+ nbtree_clear(Ri_ngbrs);
+
+ /* merge Ri with Rk */
+ merge_regions(&Ri, Ri_rs, &Rk, Rk_rs, globals);
- /* merge_values sets Ri and Rk candidate flag to FALSE. Put Rk back to TRUE if the size is too small. */
- if (Ri_count + Rk_count <
- functions->min_segment_size)
- set_candidate_flag(Rk_head, TRUE, files);
- } /* end if best neighbor != null */
- else
- G_warning
- (_("No best neighbor found in final merge for small segment, this shouldn't happen!"));
-
-
- } /* end else - pixel count was below minimum allowed */
- } /* end if neighbors found */
- else { /* no neighbors were found */
- G_warning
- (_("no neighbors found, this means only one segment was created."));
- set_candidate_flag(Ri_head, FALSE, files);
+ if (Ri_nn <= 0 || Ri.count >= globals->min_segment_size)
+ do_merge = 0;
}
- } /* end if pixel is candidate pixel */
- } /* next column */
- } /* next row */
- t++; /* to count one more "iteration" */
- } /* end if for force merge */
- else if (t > 2 && files->bounds_map == NULL)
- G_verbose_message(_("Number of passes completed: %d"), t - 1);
-#ifdef PROFILE
- end = clock();
- fprintf(stdout, "time spent merging: %g\n", merge_accum);
- fprintf(stdout, "time spent finding neighbors: %g\n", fn_accum);
- fprintf(stdout, "total time: %g\n",
- ((double)(end - start) / CLOCKS_PER_SEC));
-#endif
- return TRUE;
+ }
+ }
+ }
}
+
+ G_free(Ri.mean);
+ G_free(Rk.mean);
+ G_free(Rk_bestn.mean);
+
+ nbtree_clear(Ri_ngbrs);
+ nbtree_clear(Rk_ngbrs);
+ free(Ri_ngbrs);
+ free(Rk_ngbrs);
- /* 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;
+ return TRUE;
+}
-#ifdef DEBUG
- struct RB_TRAV trav;
-#endif
+static int find_best_neighbor(struct ngbr_stats *Ri,
+ struct NB_TREE *Ri_ngbrs,
+ struct reg_stats **Rk_rs,
+ struct ngbr_stats *Rk,
+ double *sim, int clear_cand,
+ struct globals *globals)
+{
+ int n, n_ngbrs, no_check;
+ struct rc ngbr_rc, next;
+ struct rclist rilist;
+ double tempsim, *dp;
+ int neighbors[8][2];
+ struct RB_TREE *no_check_tree; /* cells already checked */
+ struct reg_stats *rs_found;
- /* speed enhancement, any time savings to move any variables to files (mem allocation once in open_files) */
+ G_debug(4, "find_best_neighbor()");
- /* neighbor list will be a listing of pixels that are neighbors, but will be limited to just one pixel from each neighboring segment.
- * */
+ /* dynamics of the region growing algorithm
+ * some regions are growing fast, often surrounded by many small regions
+ * not all regions are equally growing, some will only grow at a later stage ? */
- /* 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 *** */
+ no_check_tree = rbtree_create(compare_rc, sizeof(struct rc));
+ ngbr_rc.row = Ri->row;
+ ngbr_rc.col = Ri->col;
+ rbtree_insert(no_check_tree, &ngbr_rc);
- /* *** initialize data *** */
- borderPixels = 0;
+ Ri->count = 1;
- segment_get(&files->iseg_seg, &R_iseg, (*R_head)->row,
- (*R_head)->col);
+ n_ngbrs = 0;
+ /* TODO: add size of largest region to reg_tree, use this as min */
+ Rk->count = globals->ncells;
- if (R_iseg == 0) { /* if seeds were provided, this is just a single non-seed pixel, only return neighbors that are segments or seeds */
+ /* go through segment, spreading outwards from head */
+ rclist_init(&rilist);
+ rclist_add(&rilist, Ri->row, Ri->col);
- functions->find_pixel_neighbors((*R_head)->row, (*R_head)->col,
- pixel_neighbors, files);
- for (n = 0; n < functions->num_pn; n++) {
+ while (rclist_drop(&rilist, &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;
+ /* remove from candidates */
+ if (clear_cand)
+ FLAG_UNSET(globals->candidate_flag, next.row, next.col);
- segment_get(&files->iseg_seg, ¤t_seg_ID,
- pixel_neighbors[n][0], pixel_neighbors[n][1]);
+ G_debug(5, "find_pixel_neighbors for row: %d , col %d",
+ next.row, next.col);
- 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? */
- }
+ globals->find_neighbors(next.row, next.col, neighbors);
+
+ n = globals->nn - 1;
+ do {
- }
- 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;
+ globals->ns.row = ngbr_rc.row = neighbors[n][0];
+ globals->ns.col = ngbr_rc.col = neighbors[n][1];
- /* Copy R in to_check and no_check data structures (don't expand them if we find them again) */
+ no_check = (ngbr_rc.row < globals->row_min ||
+ ngbr_rc.row >= globals->row_max ||
+ ngbr_rc.col < globals->col_min ||
+ ngbr_rc.col >= globals->col_max);
- 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. */
+ n_ngbrs += no_check;
- /* 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?"));
- }
+ if (!no_check) {
- while (to_check != NULL) { /* removing from to_check list as we go, NOT iterating over the list. */
+ no_check = ((FLAG_GET(globals->null_flag, ngbr_rc.row,
+ ngbr_rc.col)) != 0);
+ n_ngbrs += no_check;
- functions->find_pixel_neighbors(to_check->row,
- to_check->col,
- pixel_neighbors, files);
+ if (!no_check) {
- /* 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);
+ if (!rbtree_find(no_check_tree, &ngbr_rc)) {
- /* 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++) {
+ /* not yet checked, don't check it again */
+ rbtree_insert(no_check_tree, &ngbr_rc);
- /* 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;
+ /* get Rk ID */
+ segment_get(&globals->rid_seg,
+ (void *) &(globals->ns.id),
+ ngbr_rc.row, ngbr_rc.col);
- tree_pix.row = pixel_neighbors[n][0];
- tree_pix.col = pixel_neighbors[n][1];
+ if (globals->ns.id == Ri->id) {
- if (rbtree_find(no_check_tree, &tree_pix) == FALSE) { /* want to check this neighbor */
- segment_get(&files->iseg_seg, ¤t_seg_ID,
- pixel_neighbors[n][0],
- pixel_neighbors[n][1]);
+ /* want to check this neighbor's neighbors */
+ rclist_add(&rilist, ngbr_rc.row, ngbr_rc.col);
+ Ri->count++;
+ }
+ else {
- rbtree_insert(no_check_tree, &tree_pix); /* don't check it again */
+ /* new neighbor ? */
+ if (nbtree_find(Ri_ngbrs, &globals->ns) == NULL) {
- 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. */
+ /* get values for Rk */
+ globals->rs.id = globals->ns.id;
+ rs_found = rgtree_find(globals->reg_tree,
+ &(globals->rs));
+ if (rs_found != NULL) {
+ globals->ns.mean = rs_found->mean;
+ globals->ns.count = rs_found->count;
+ }
+ else {
+ segment_get(&globals->bands_seg,
+ (void *)globals->bands_val,
+ ngbr_rc.row, ngbr_rc.col);
- /* 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. */
+ globals->ns.mean = globals->bands_val;
+ globals->ns.count = 1;
+ }
+ /* next Rk = globals->ns is now complete */
- }
- 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 ??? */
+ tempsim = (globals->calculate_similarity)(Ri, &globals->ns, globals);
- if (!rbtree_find(known_iseg, ¤t_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,
- ¤t_seg_ID);
+ if (compare_double(tempsim, *sim) == -1) {
+ *sim = tempsim;
+ /* copy temp Rk to Rk */
+ dp = Rk->mean;
+ *Rk = globals->ns;
+ Rk->mean = dp;
+ memcpy(Rk->mean, globals->ns.mean,
+ globals->datasize);
+ *Rk_rs = rs_found;
+ }
+ else if (compare_double(tempsim, *sim) == 0) {
+ /* resolve ties: prefer smaller regions */
- /* 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 (Rk->count > globals->ns.count) {
+ /* copy temp Rk to Rk */
+ dp = Rk->mean;
+ *Rk = globals->ns;
+ Rk->mean = dp;
+ memcpy(Rk->mean,
+ globals->ns.mean,
+ globals->datasize);
+ *Rk_rs = rs_found;
+ }
+ }
+
+ n_ngbrs++;
+ nbtree_insert(Ri_ngbrs, &globals->ns);
}
- 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? */
- }
}
+ }
+ }
+ }
+ } while (n--); /* end do loop - next neighbor */
+ } /* while there are cells to check */
+ /* clean up */
+ rbtree_destroy(no_check_tree);
- } /*end if for pixel_neighbor was in "don't check" list */
- } /* end for loop - next pixel neighbor */
- } /* end while to_check has more elements */
+ return n_ngbrs;
+}
- /* clean up */
- rbtree_destroy(no_check_tree);
- rbtree_destroy(known_iseg);
- }
- return borderPixels;
- }
+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;
- 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.
- */
+ /* east */
+ neighbors[1][0] = p_row;
+ neighbors[1][1] = p_col + 1;
- /* north */
- pixel_neighbors[0][1] = p_col;
- pixel_neighbors[0][0] = p_row - 1;
+ /* south */
+ neighbors[2][0] = p_row + 1;
+ neighbors[2][1] = p_col;
- /* east */
- pixel_neighbors[1][0] = p_row;
- pixel_neighbors[1][1] = p_col + 1;
+ /* west */
+ neighbors[3][0] = p_row;
+ neighbors[3][1] = p_col - 1;
- /* south */
- pixel_neighbors[2][1] = p_col;
- pixel_neighbors[2][0] = p_row + 1;
+ return;
+}
- /* west */
- pixel_neighbors[3][0] = p_row;
- pixel_neighbors[3][1] = p_col - 1;
+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);
- return TRUE;
- }
+ /* get the 4 diagonal neighbors */
+ /* north-west */
+ neighbors[4][0] = p_row - 1;
+ neighbors[4][1] = p_col - 1;
- 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);
+ /* north-east */
+ neighbors[5][0] = p_row - 1;
+ neighbors[5][1] = p_col + 1;
- /* and then the diagonals: */
+ /* south-west */
+ neighbors[6][0] = p_row + 1;
+ neighbors[6][1] = p_col - 1;
- /* north west */
- pixel_neighbors[4][0] = p_row - 1;
- pixel_neighbors[4][1] = p_col - 1;
+ /* south-east */
+ neighbors[7][0] = p_row + 1;
+ neighbors[7][1] = p_col + 1;
- /* north east */
- pixel_neighbors[5][0] = p_row - 1;
- pixel_neighbors[5][1] = p_col + 1;
+ return;
+}
- /* south east */
- pixel_neighbors[6][0] = p_row + 1;
- pixel_neighbors[6][1] = p_col + 1;
+/* 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;
- /* south west */
- pixel_neighbors[7][0] = p_row + 1;
- pixel_neighbors[7][1] = p_col - 1;
+ /* squared euclidean distance, sum the square differences for each dimension */
+ do {
+ diff = Ri->mean[n] - Rk->mean[n];
+
+ val += diff * diff;
+ } while (n--);
- return TRUE;
- }
+ return val;
+}
- /* 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;
+static int search_neighbors(struct ngbr_stats *Ri,
+ struct NB_TREE *Ri_ngbrs,
+ double *sim,
+ struct ngbr_stats *Rk,
+ struct globals *globals)
+{
+ double tempsim, *dp;
+ struct NB_TRAV travngbr;
+ struct ngbr_stats *next;
- /* get values for pixel b */
- segment_get(&files->bands_seg, (void *)files->second_val, b->row,
- b->col);
+ G_debug(4, "search_neighbors");
- /* 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]);
- }
+ nbtree_init_trav(&travngbr, Ri_ngbrs);
+ Rk->count = globals->ncells;
- /* use squared distance, save the calculation time. We squared the similarity threshold earlier to allow for this. */
- /* val = sqrt(val); */
+ while ((next = nbtree_traverse(&travngbr))) {
+ tempsim = (globals->calculate_similarity)(Ri, next, globals);
- return val;
+ if (compare_double(tempsim, *sim) == -1) {
+ *sim = tempsim;
+ dp = Rk->mean;
+ *Rk = *next;
+ Rk->mean = dp;
+ memcpy(Rk->mean, next->mean, globals->datasize);
+ }
+ else if (compare_double(tempsim, *sim) == 0) {
+ /* resolve ties, prefer smaller regions */
+ G_debug(4, "resolve ties");
+ if (Rk->count > next->count) {
+ dp = Rk->mean;
+ *Rk = *next;
+ Rk->mean = dp;
+ memcpy(Rk->mean, next->mean, globals->datasize);
+ }
+ }
}
- double calculate_manhattan_similarity(struct pixels *a, struct pixels *b,
- struct files *files,
- struct functions *functions)
- {
- double val = 0;
- int n;
+ return 1;
+}
- /* get values for pixel b */
- segment_get(&files->bands_seg, (void *)files->second_val, b->row,
- b->col);
+static int merge_regions(struct ngbr_stats *Ri, struct reg_stats *Ri_rs,
+ struct ngbr_stats *Rk, struct reg_stats *Rk_rs,
+ struct globals *globals)
+{
+ int n, do_cand;
+ int Ri_id, Rk_id, larger_id, R_id;
+ struct rc next, ngbr_rc;
+ struct rclist rlist;
+ int neighbors[8][2];
+ struct reg_stats *new_rs, old_rs;
- /* 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? */
- }
+ G_debug(4, "merge_regions");
- return val;
+ Ri_id = Ri->id;
+ Rk_id = Rk->id;
- }
+ /* update segment number and clear candidate flag */
+
+ /* cases
+ * Ri, Rk are single cells, not in tree
+ * Ri, Rk are both in the tree
+ * Ri is in the tree, Rk is not
+ * Rk is in the tree, Ri is not
+ */
- /* TODO: add shape parameter...
- *
- In the eCognition literature, we find that the key factor in the
- multi-scale segmentation algorithm used by Definiens is the scale
- factor f:
+ /* for each member of Ri and Rk, write new average bands values and segment values */
- 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
+ if (Ri->count >= Rk->count) {
+ larger_id = Ri_id;
- 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, and Pbbox the perimeter of the bounding box of the object.
- */
+ if (Ri_rs) {
+ new_rs = Ri_rs;
+ }
+ else {
+ new_rs = &(globals->rs);
+ new_rs->id = Ri_id;
+ new_rs->count = 1;
+
+ memcpy(new_rs->sum, Ri->mean, globals->datasize);
+ memcpy(new_rs->mean, Ri->mean, globals->datasize);
+ }
- 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;
+ /* add Rk */
+ if (Rk_rs) {
+ new_rs->count += Rk_rs->count;
+ n = globals->nbands - 1;
+ do {
+ new_rs->sum[n] += Rk_rs->sum[n];
+ } while (n--);
+ }
+ else {
+ new_rs->count++;
+ n = globals->nbands - 1;
+ do {
+ new_rs->sum[n] += Rk->mean[n];
+ } while (n--);
+ }
- /*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);
+ n = globals->nbands - 1;
+ do {
+ new_rs->mean[n] = new_rs->sum[n] / new_rs->count;
+ } while (n--);
+
+ memcpy(Ri->mean, new_rs->mean, globals->datasize);
+ Ri->count = new_rs->count;
- 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);
+ if (!Ri_rs) {
+ /* add to tree */
+ rgtree_insert(globals->reg_tree, new_rs);
+ }
+ if (Rk_rs) {
+ /* remove from tree */
+ old_rs.id = Rk_id;
+ rgtree_remove(globals->reg_tree, &old_rs);
+ }
+ }
+ else {
+ larger_id = Rk_id;
+ Ri->id = Rk_id;
- 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);
+ if (Rk_rs) {
+ new_rs = Rk_rs;
}
+ else {
+ new_rs = &(globals->rs);
+ new_rs->id = Rk_id;
+ new_rs->count = 1;
- /* 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);
-#endif
+ memcpy(new_rs->sum, Rk->mean, globals->datasize);
+ memcpy(new_rs->mean, Rk->mean, globals->datasize);
+ }
- /* 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 */
+ /* add Ri */
+ if (Ri_rs) {
+ new_rs->count += Ri_rs->count;
+ n = globals->nbands - 1;
+ do {
+ new_rs->sum[n] += Ri_rs->sum[n];
+ } while (n--);
}
- 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);
+ else {
+ new_rs->count++;
+ n = globals->nbands - 1;
+ do {
+ new_rs->sum[n] += Ri->mean[n];
+ } while (n--);
}
- /* merged two segments, decrement count if Rk was an actual segment (not a non-seed pixel) */
- if (Rk_iseg > 0)
- files->nsegs--;
+ n = globals->nbands - 1;
+ do {
+ new_rs->mean[n] = new_rs->sum[n] / new_rs->count;
+ } while (n--);
- return TRUE;
+ memcpy(Ri->mean, new_rs->mean, globals->datasize);
+ Ri->count = new_rs->count;
+
+ if (!Rk_rs) {
+ /* add to tree */
+ rgtree_insert(globals->reg_tree, new_rs);
+ }
+ if (Ri_rs) {
+ /* remove from tree */
+ old_rs.id = Ri_id;
+ rgtree_remove(globals->reg_tree, &old_rs);
+ }
}
- /* 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;
+ if (larger_id == Ri_id) {
+ /* Ri is already updated, including candidate flags
+ * need to clear candidate flag for Rk and set new id */
+
+ /* the actual merge: change region id */
+ segment_put(&globals->rid_seg, (void *) &Ri_id, Rk->row, Rk->col);
- /* 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;
+ do_cand = 0;
+ if (FLAG_GET(globals->candidate_flag, Rk->row, Rk->col)) {
+ /* clear candidate flag */
+ FLAG_UNSET(globals->candidate_flag, Rk->row, Rk->col);
+ globals->candidate_count--;
+ do_cand = 1;
}
- if (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++;
- }
+ rclist_init(&rlist);
+ rclist_add(&rlist, Rk->row, Rk->col);
- /* calculate the mean */
- for (n = 0; n < files->nbands; n++) {
- files->second_val[n] = files->second_val[n] / count;
+ while (rclist_drop(&rlist, &next)) {
+
+ if (do_cand) {
+ /* clear candidate flag */
+ FLAG_UNSET(globals->candidate_flag, next.row, next.col);
+ globals->candidate_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? */
+ globals->find_neighbors(next.row, next.col, neighbors);
+
+ n = globals->nn - 1;
+ do {
- /* 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);
- }
+ ngbr_rc.row = neighbors[n][0];
+ ngbr_rc.col = neighbors[n][1];
- }
+ if (ngbr_rc.row >= globals->row_min &&
+ ngbr_rc.row < globals->row_max &&
+ ngbr_rc.col >= globals->col_min &&
+ ngbr_rc.col < globals->col_max) {
- return TRUE;
+ if (!(FLAG_GET(globals->null_flag, ngbr_rc.row, ngbr_rc.col))) {
+
+ segment_get(&globals->rid_seg, (void *) &R_id, ngbr_rc.row, ngbr_rc.col);
+
+ if (R_id == Rk_id) {
+ /* the actual merge: change region id */
+ segment_put(&globals->rid_seg, (void *) &Ri_id, ngbr_rc.row, ngbr_rc.col);
+
+ /* want to check this neighbor's neighbors */
+ rclist_add(&rlist, ngbr_rc.row, ngbr_rc.col);
+ }
+ }
+ }
+ } while (n--);
+ }
}
+ else {
+ /* larger_id == Rk_id */
- /* 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;
+ /* clear candidate flag for Rk */
+ if (FLAG_GET(globals->candidate_flag, Rk->row, Rk->col)) {
+ set_candidate_flag(Rk, FALSE, globals);
+ }
- for (current = head; current != NULL; current = current->next) {
+ /* update region id for Ri */
+ /* the actual merge: change region id */
+ segment_put(&globals->rid_seg, (void *) &Rk_id, Ri->row, Ri->col);
- 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"));
+ rclist_init(&rlist);
+ rclist_add(&rlist, Ri->row, Ri->col);
+
+ while (rclist_drop(&rlist, &next)) {
+
+ globals->find_neighbors(next.row, next.col, neighbors);
+
+ n = globals->nn - 1;
+ do {
+
+ ngbr_rc.row = neighbors[n][0];
+ ngbr_rc.col = neighbors[n][1];
+
+ if (ngbr_rc.row >= globals->row_min &&
+ ngbr_rc.row < globals->row_max &&
+ ngbr_rc.col >= globals->col_min &&
+ ngbr_rc.col < globals->col_max) {
+
+
+ if (!(FLAG_GET(globals->null_flag, ngbr_rc.row, ngbr_rc.col))) {
+
+ segment_get(&globals->rid_seg, (void *) &R_id, ngbr_rc.row, ngbr_rc.col);
+
+ if (R_id == Ri_id) {
+ /* the actual merge: change region id */
+ segment_put(&globals->rid_seg, (void *) &Rk_id, ngbr_rc.row, ngbr_rc.col);
+
+ /* want to check this neighbor's neighbors */
+ rclist_add(&rlist, ngbr_rc.row, ngbr_rc.col);
+ }
+ }
+ }
+ } while (n--);
}
- return TRUE;
}
+
+ if (Rk_id > 0)
+ globals->n_regions--;
- /* 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;
+ 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 */
- }
+static int set_candidate_flag(struct ngbr_stats *head, int value, struct globals *globals)
+{
+ int n, R_id;
+ struct rc next, ngbr_rc;
+ struct rclist rlist;
+ int neighbors[8][2];
- return TRUE;
+ G_debug(4, "set_candidate_flag");
+
+ if (!(FLAG_GET(globals->candidate_flag, head->row, head->col)) != value) {
+ G_warning(_("Candidate flag is already %s"), value ? _("set") : _("unset"));
+ return FALSE;
}
- /* functions used by binary tree to compare items */
+ rclist_init(&rlist);
+ rclist_add(&rlist, head->row, head->col);
- /* speed enhancement: Maybe changing this would be an improvement? "static" was used in break_polygons.c extern was suggested in docs. */
+ /* (un)set candidate flag */
+ if (value == TRUE) {
+ FLAG_SET(globals->candidate_flag, head->row, head->col);
+ globals->candidate_count++;
+ }
+ else {
+ FLAG_UNSET(globals->candidate_flag, head->row, head->col);
+ globals->candidate_count--;
+ }
- int compare_ids(const void *first, const void *second)
- {
- int *a = (int *)first, *b = (int *)second;
+ while (rclist_drop(&rlist, &next)) {
- if (*a < *b)
- return -1;
- else if (*a > *b)
- return 1;
- else if (*a == *b)
- return 0;
+ globals->find_neighbors(next.row, next.col, neighbors);
+
+ n = globals->nn - 1;
+ do {
+ ngbr_rc.row = neighbors[n][0];
+ ngbr_rc.col = neighbors[n][1];
- G_warning(_("find neighbors: Bug in binary tree!"));
- return 1;
+ if (ngbr_rc.row >= globals->row_min &&
+ ngbr_rc.row < globals->row_max &&
+ ngbr_rc.col >= globals->col_min &&
+ ngbr_rc.col < globals->col_max) {
- }
+ if (!(FLAG_GET(globals->null_flag, ngbr_rc.row, ngbr_rc.col))) {
- int compare_pixels(const void *first, const void *second)
- {
- struct pixels *a = (struct pixels *)first, *b =
- (struct pixels *)second;
+ if (!(FLAG_GET(globals->candidate_flag, ngbr_rc.row, ngbr_rc.col)) == value) {
- 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->col > b->col)
- return 1;
- }
- /* same row and col */
- return 0;
- }
+ segment_get(&globals->rid_seg, (void *) &R_id, ngbr_rc.row, ngbr_rc.col);
- /* 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;
+ if (R_id == head->id) {
+ /* want to check this neighbor's neighbors */
+ rclist_add(&rlist, ngbr_rc.row, ngbr_rc.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);
+ /* (un)set candidate flag */
+ if (value == TRUE) {
+ FLAG_SET(globals->candidate_flag, ngbr_rc.row, ngbr_rc.col);
+ globals->candidate_count++;
+ }
+ else {
+ FLAG_UNSET(globals->candidate_flag, ngbr_rc.row, ngbr_rc.col);
+ globals->candidate_count--;
+ }
+ }
+ }
}
- else
- FLAG_UNSET(files->candidate_flag, row, col);
}
- }
- return TRUE;
+ } while (n--);
}
+
+ return TRUE;
+}
Modified: grass-addons/grass7/imagery/i.segment.xl/flag.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/flag.c 2012-08-27 18:43:23 UTC (rev 52938)
+++ grass-addons/grass7/imagery/i.segment.xl/flag.c 2012-08-27 19:08:40 UTC (rev 52939)
@@ -1,23 +1,13 @@
#include <grass/gis.h>
+#include <grass/glocale.h>
#include "flag.h"
-int flag_clear_all(FLAG * flags)
-{
- register int r, c;
-
- for (r = 0; r < flags->nrows; r++) {
- for (c = 0; c < flags->leng; c++) {
- flags->array[r][c] = 0;
- }
- }
-
- return 0;
-}
-
FLAG *flag_create(int nrows, int ncols)
{
unsigned char *temp;
+
FLAG *new_flag;
+
register int i;
new_flag = (FLAG *) G_malloc(sizeof(FLAG));
@@ -26,15 +16,23 @@
new_flag->leng = (ncols + 7) / 8;
new_flag->array =
(unsigned char **)G_malloc(nrows * sizeof(unsigned char *));
+
+ if (!new_flag->array)
+ G_fatal_error(_("Out of memory!"));
temp =
(unsigned char *)G_malloc(nrows * new_flag->leng *
sizeof(unsigned char));
+ if (!temp)
+ G_fatal_error(_("Out of memory!"));
+
for (i = 0; i < nrows; i++) {
new_flag->array[i] = temp;
temp += new_flag->leng;
}
+
+ flag_clear_all(new_flag);
return (new_flag);
}
@@ -47,3 +45,16 @@
return 0;
}
+
+int flag_clear_all(FLAG * flags)
+{
+ register int r, c;
+
+ for (r = 0; r < flags->nrows; r++) {
+ for (c = 0; c < flags->leng; c++) {
+ flags->array[r][c] = 0;
+ }
+ }
+
+ return 0;
+}
Modified: grass-addons/grass7/imagery/i.segment.xl/flag.h
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/flag.h 2012-08-27 18:43:23 UTC (rev 52938)
+++ grass-addons/grass7/imagery/i.segment.xl/flag.h 2012-08-27 19:08:40 UTC (rev 52939)
@@ -20,8 +20,6 @@
** FLAG *flags;
** sets all values in flags to zero.
**
- * following 3 were changed to macros, same usage
- *
** flag_unset(flags, row, col)
** FLAG *flags;
** int row, col;
@@ -58,8 +56,9 @@
(flags)->array[(row)][(col)>>3] & (1<<((col) & 7))
/* flag.c */
-int flag_clear_all(FLAG *);
FLAG *flag_create(int, int);
int flag_destroy(FLAG *);
+int flag_clear_all(FLAG *);
+
#endif /* __FLAG_H__ */
Deleted: grass-addons/grass7/imagery/i.segment.xl/i.segment.html
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/i.segment.html 2012-08-27 18:43:23 UTC (rev 52938)
+++ grass-addons/grass7/imagery/i.segment.xl/i.segment.html 2012-08-27 19:08:40 UTC (rev 52939)
@@ -1,98 +0,0 @@
-<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.
-
-<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.
-
-<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.
-<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>
-<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.
-<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>
-<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:
-<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>
-<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.
-
-<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.
-
-<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>
-
-<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>
-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>
-<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>
-<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>
-<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.
-
-<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>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>
-</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>.
-<p>
-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>
-<a href="i.group.html">i.group</a>,
-<a href="i.maxlik.html">i.maxlik</a>,
-<a href="r.fuzzy">r.fuzzy</a>,
-<a href="i.smap.html">i.smap</a>,
-<a href="r.seg.html">r.seg</a> (Addon)
-</em>
-
-<h2>AUTHORS</h2>
-Eric Momsen - North Dakota State University
Copied: grass-addons/grass7/imagery/i.segment.xl/i.segment.xl.html (from rev 52938, grass-addons/grass7/imagery/i.segment.xl/i.segment.html)
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/i.segment.xl.html (rev 0)
+++ grass-addons/grass7/imagery/i.segment.xl/i.segment.xl.html 2012-08-27 19:08:40 UTC (rev 52939)
@@ -0,0 +1,98 @@
+<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.
+
+<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.
+
+<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.
+<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>
+<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.
+<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>
+<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:
+<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>
+<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.
+
+<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.
+
+<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>
+
+<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>
+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>
+<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>
+<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>
+<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.
+
+<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>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>
+</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>.
+<p>
+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>
+<a href="i.group.html">i.group</a>,
+<a href="i.maxlik.html">i.maxlik</a>,
+<a href="r.fuzzy">r.fuzzy</a>,
+<a href="i.smap.html">i.smap</a>,
+<a href="r.seg.html">r.seg</a> (Addon)
+</em>
+
+<h2>AUTHORS</h2>
+Eric Momsen - North Dakota State University
Modified: grass-addons/grass7/imagery/i.segment.xl/iseg.h
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/iseg.h 2012-08-27 18:43:23 UTC (rev 52938)
+++ grass-addons/grass7/imagery/i.segment.xl/iseg.h 2012-08-27 19:08:40 UTC (rev 52939)
@@ -13,119 +13,103 @@
*****************************************************************************/
#include <grass/segment.h>
-#include <grass/linkm.h>
#include "flag.h"
+#include "regtree.h"
+#include "ngbrtree.h"
-/* PROFILE will add some rough time checks for finding neighbors, merging, and pass times. */
-/* #define PROFILE */
-
-/* SIGNPOST will add some fprintf statments to indicate what segments are being merged.
- * other diagnostics could be included here for during further development and speed improvements.
- */
-/* #define SIGNPOST */
-
-/* pixel stack */
-struct pixels
+/* row/col list */
+struct rc
{
- struct pixels *next;
+ struct rc *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... */
};
+struct rclist
+{
+ struct rc *tail, *head;
+};
+
/* input and output files, as well as some other processing info */
-struct files
+struct globals
{
/* user parameters */
char *image_group;
- int weighted; /* 0 if false/not selected, so we should scale input. 1 if the scaling should be skipped */
+ int weighted; /* 0 if false/not selected, so we should scale input.
+ * 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;
+ int min_segment_size; /* smallest number of pixels/cells allowed in a final segment */
+ int end_t; /* maximum number of iterations */
+ int mb;
/* region info */
int nrows, ncols;
+ int row_min, row_max, col_min, col_max; /* region constraints */
+ int ncells;
- /* files */
char *out_name; /* name of output raster map */
- const char *seeds_map, *seeds_mapset, *bounds_map, *bounds_mapset; /* optional segment seeds and polygon constraints/boundaries */
- char *out_band; /* for segment average values */
+ char *seeds, *bounds_map; /* optional segment seeds and polygon constraints/boundaries */
+ CELL lower_bound, upper_bound;
+ const char *bounds_mapset;
+ char *out_band; /* indicator for segment heterogeneity */
/* file processing */
- /* bands_seg is initialized with the input raster valuess, then is updated with current mean values for the segment. */
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 */
- double *second_val; /* to hold values at second point for similarity comparison */
- int bounds_val, current_bound;
- int minrow, maxrow, mincol, maxcol;
+ SEGMENT bands_seg, /* input group with one or more bands */
+ bounds_seg,
+ rid_seg;
+ DCELL *bands_min, *bands_max;
+ DCELL *bands_val; /* array to hold all input values for one cell */
+ DCELL *second_val; /* array to hold all input values for another cell */
/* results */
- SEGMENT iseg_seg; /* segment ID assignment */
- int nsegs; /* number of segments */
+ struct RG_TREE *reg_tree; /* search tree with region stats */
+ struct reg_stats rs;
+ struct ngbr_stats ns;
+ size_t datasize; /* nbands * sizeof(double) */
+ int n_regions;
/* processing flags */
- /* candidate flag for if a cell/segment has already been merged in that pass. */
- /* seeds flag for if a cell/segment is a seed (can be Ri to start a merge). All cells are valid seeds if a starting seeds map is not supplied. */
- FLAG *candidate_flag, *null_flag, *orig_null_flag, *seeds_flag;
+ FLAG *candidate_flag, *null_flag; /*TODO, need some way to remember MASK/NULL values. Was using -1, 0, 1 in int array. Better to use 2 FLAG structures, better readibility? */
- /* memory management, linked lists */
- struct link_head *token; /* for linkm.h linked list memory management. */
+ /* number of remaining cells to check */
+ int candidate_count;
+ /* functions */
+
+ void (*find_neighbors) (int, int, int[8][2]); /*parameters: row, col, neighbors */
+ double (*calculate_similarity) (struct ngbr_stats *,
+ struct ngbr_stats *,
+ struct globals *); /*parameters: two regions to compare */
};
-struct functions
-{
- int method; /* Segmentation method */
- 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 */
- /* 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 */
-
- /* max number of iterations/passes */
- int end_t;
-
- 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 */
-
- /* 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
- * The goal is to speed processing, since in the end many of these very similar neighbors will be merged.
- * But the problem is that the find_segment_neighbors() function only returns a single pixel, not the entire segment membership.
- * The commented out code actually slowed down processing times in the first tries. */
-#ifdef VCLOSE
- double very_close; /* segments with very_close similarity will be merged without changing or checking the candidate flag.
- * The algorithm will continue looking for the "most similar" neighbor that isn't "very close". */
-#endif
-};
-
-
/* parse_args.c */
/* gets input from user, validates, and sets up functions */
-int parse_args(int, char *[], struct files *, struct functions *);
+int parse_args(int, char *[], struct globals *);
/* open_files.c */
-int open_files(struct files *, struct functions *);
+int open_files(struct globals *);
/* create_isegs.c */
-int create_isegs(struct files *, struct functions *);
-int region_growing(struct files *, struct functions *);
-int find_segment_neighbors(struct pixels **, struct pixels **, int *,
- struct files *, struct functions *);
-int set_candidate_flag(struct pixels *, int, struct files *);
-int merge_values(struct pixels *, struct pixels *, int, int, struct files *);
-int merge_pixels(struct pixels *, int, struct files *);
-int find_four_pixel_neighbors(int, int, int[][2], struct files *);
-int find_eight_pixel_neighbors(int, int, int[8][2], struct files *);
-double calculate_euclidean_similarity(struct pixels *, struct pixels *,
- struct files *, struct functions *);
-double calculate_manhattan_similarity(struct pixels *, struct pixels *,
- struct files *, struct functions *);
-int my_dispose_list(struct link_head *, struct pixels **);
-int compare_ids(const void *, const void *);
-int compare_pixels(const void *, const void *);
-int set_all_candidate_flags(struct files *);
+int create_isegs(struct globals *);
+int region_growing(struct globals *);
+void find_four_neighbors(int, int, int[][2]);
+void find_eight_neighbors(int, int, int[8][2]);
+double calculate_euclidean_similarity(struct ngbr_stats *,
+ struct ngbr_stats *, struct globals *);
+
+/* rclist.c */
+void rclist_init(struct rclist *);
+void rclist_add(struct rclist *, int, int);
+int rclist_drop(struct rclist *, struct rc *);
+void rclist_destroy(struct rclist *);
+
+
/* write_output.c */
-int write_output(struct files *);
-int close_files(struct files *);
+int write_output(struct globals *);
+int close_files(struct globals *);
Modified: grass-addons/grass7/imagery/i.segment.xl/main.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/main.c 2012-08-27 18:43:23 UTC (rev 52938)
+++ grass-addons/grass7/imagery/i.segment.xl/main.c 2012-08-27 19:08:40 UTC (rev 52939)
@@ -11,8 +11,8 @@
* for details.
*
*
- * NOTE: the names "segment" and "SEG" are already used by the Segmentation
- * Library for the data files/tiling, so "iseg" (image segmentation)
+ * NOTE: the word "segment" is already used by the Segmentation
+ * Library for the data files/tiling, so iseg (image segmentation)
* will be used to refer to the image segmentation.
*
*
@@ -25,8 +25,7 @@
int main(int argc, char *argv[])
{
- struct files files; /* input and output file descriptors, data structure, buffers */
- struct functions functions; /* function pointers and parameters for the calculations */
+ struct globals globals; /* input and output file descriptors, data structure, buffers */
struct GModule *module;
G_gisinit(argv[0]);
@@ -37,25 +36,25 @@
module->description =
_("Outputs a single segmented map (raster) based on input values in an image group.");
- if (parse_args(argc, argv, &files, &functions) != TRUE)
+ if (parse_args(argc, argv, &globals) != TRUE)
G_fatal_error(_("Error in parse_args()"));
G_debug(1, "Main: starting open_files()");
- if (open_files(&files, &functions) != TRUE)
+ if (open_files(&globals) != TRUE)
G_fatal_error(_("Error in open_files()"));
G_debug(1, "Main: starting create_isegs()");
- if (create_isegs(&files, &functions) != TRUE)
+ if (create_isegs(&globals) != TRUE)
G_fatal_error(_("Error in create_isegs()"));
G_debug(1, "Main: starting write_output()");
- if (write_output(&files) != TRUE)
+ if (write_output(&globals) != TRUE)
G_fatal_error(_("Error in write_output()"));
G_debug(1, "Main: starting close_files()");
- close_files(&files);
+ close_files(&globals);
- G_done_msg("Number of segments created: %d", files.nsegs);
+ G_done_msg(_("Number of segments created: %d"), globals.n_regions);
exit(EXIT_SUCCESS);
}
Added: grass-addons/grass7/imagery/i.segment.xl/ngbrtree.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/ngbrtree.c (rev 0)
+++ grass-addons/grass7/imagery/i.segment.xl/ngbrtree.c 2012-08-27 19:08:40 UTC (rev 52939)
@@ -0,0 +1,563 @@
+/*!
+ * \file rbtree.c
+ *
+ * \brief binary search tree
+ *
+ * Generic balanced binary search tree (Red Black Tree) implementation
+ *
+ * (C) 2009 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public License
+ * (>=v2). Read the file COPYING that comes with GRASS for details.
+ *
+ * \author Original author Julienne Walker 2003, 2008
+ * GRASS implementation Markus Metz, 2009
+ */
+
+/* balanced binary search tree implementation
+ *
+ * this one is a Red Black Tree, no parent pointers, no threads
+ * The core code comes from Julienne Walker's tutorials on binary search trees
+ * original license: public domain
+ * http://eternallyconfuzzled.com/tuts/datastructures/jsw_tut_rbtree.aspx
+ * some ideas come from libavl (GPL >= 2)
+ *
+ * Red Black Trees are used to maintain a data structure with
+ * search, insertion and deletion in O(log N) time
+ */
+
+#include <assert.h>
+#include <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "ngbrtree.h"
+
+/* internal functions */
+static struct NB_NODE *nbtree_single(struct NB_NODE *, int);
+static struct NB_NODE *nbtree_double(struct NB_NODE *, int);
+static struct ngbr_stats *nbtree_first(struct NB_TRAV *);
+static struct ngbr_stats *nbtree_next(struct NB_TRAV *);
+static struct NB_NODE *nbtree_make_node(size_t, struct ngbr_stats *);
+static int is_red(struct NB_NODE *);
+
+
+int cmp_ngbr(struct ngbr_stats *a, struct ngbr_stats *b)
+{
+ return (a->id < b->id ? -1 : (a->id > b->id));
+}
+
+
+/* create new tree and initialize
+ * returns pointer to new tree, NULL for memory allocation error
+ */
+struct NB_TREE *nbtree_create(int nbands, size_t rb_datasize)
+{
+ struct NB_TREE *tree = (struct NB_TREE *)malloc(sizeof(struct NB_TREE));
+
+ if (tree == NULL) {
+ G_warning("RB tree: Out of memory!");
+ return NULL;
+ }
+
+ tree->datasize = rb_datasize;
+ tree->count = 0;
+ tree->nbands = nbands;
+ tree->root = NULL;
+
+ return tree;
+}
+
+/* add an item to a tree
+ * non-recursive top-down insertion
+ * the algorithm does not allow duplicates and also does not warn about a duplicate
+ * returns 1 on success, 0 on failure
+ */
+int nbtree_insert(struct NB_TREE *tree, struct ngbr_stats *data)
+{
+ assert(tree && data);
+
+ if (tree->root == NULL) {
+ /* create a new root node for tree */
+ tree->root = nbtree_make_node(tree->datasize, data);
+ if (tree->root == NULL)
+ return 0;
+ }
+ else {
+ struct NB_NODE head = { 0, {0, 0}, {0, 0, 0, 0, 0} }; /* False tree root */
+ struct NB_NODE *g, *t; /* Grandparent & parent */
+ struct NB_NODE *p, *q; /* Iterator & parent */
+ int dir = 0, last = 0;
+
+ /* Set up helpers */
+ t = &head;
+ g = p = NULL;
+ q = t->link[1] = tree->root;
+
+ /* Search down the tree */
+ for (;;) {
+ if (q == NULL) {
+ /* Insert new node at the bottom */
+ p->link[dir] = q = nbtree_make_node(tree->datasize, data);
+ if (q == NULL)
+ return 0;
+ }
+ else if (is_red(q->link[0]) && is_red(q->link[1])) {
+ /* Color flip */
+ q->red = 1;
+ q->link[0]->red = 0;
+ q->link[1]->red = 0;
+ }
+
+ /* Fix red violation */
+ if (is_red(q) && is_red(p)) {
+ int dir2 = t->link[1] == g;
+
+ if (q == p->link[last])
+ t->link[dir2] = nbtree_single(g, !last);
+ else
+ t->link[dir2] = nbtree_double(g, !last);
+ }
+
+ last = dir;
+ dir = cmp_ngbr(&(q->data), data);
+
+ /* Stop if found. This check also disallows duplicates in the tree */
+ if (dir == 0)
+ break;
+
+ dir = dir < 0;
+
+ /* Move the helpers down */
+ if (g != NULL)
+ t = g;
+
+ g = p, p = q;
+ q = q->link[dir];
+ }
+
+ /* Update root */
+ tree->root = head.link[1];
+ }
+
+ /* Make root black */
+ tree->root->red = 0;
+
+ tree->count++;
+
+ return 1;
+}
+
+/* remove an item from a tree that matches given data
+ * non-recursive top-down removal
+ * returns 1 on successful removal
+ * returns 0 if data item was not found
+ */
+int nbtree_remove(struct NB_TREE *tree, struct ngbr_stats *data)
+{
+ struct NB_NODE head = { 0, {0, 0}, {0, 0, 0, 0, 0} }; /* False tree root */
+ struct NB_NODE *q, *p, *g; /* Helpers */
+ struct NB_NODE *f = NULL; /* Found item */
+ int dir = 1, removed = 0;
+
+ assert(tree && data);
+
+ if (tree->root == NULL) {
+ return 0; /* empty tree, nothing to remove */
+ }
+
+ /* Set up helpers */
+ q = &head;
+ g = p = NULL;
+ q->link[1] = tree->root;
+
+ /* Search and push a red down */
+ while (q->link[dir] != NULL) {
+ int last = dir;
+
+ /* Update helpers */
+ g = p, p = q;
+ q = q->link[dir];
+ dir = cmp_ngbr(&(q->data), data);
+
+ /* Save found node */
+ if (dir == 0)
+ f = q;
+
+ dir = dir < 0;
+
+ /* Push the red node down */
+ if (!is_red(q) && !is_red(q->link[dir])) {
+ if (is_red(q->link[!dir]))
+ p = p->link[last] = nbtree_single(q, dir);
+ else if (!is_red(q->link[!dir])) {
+ struct NB_NODE *s = p->link[!last];
+
+ if (s != NULL) {
+ if (!is_red(s->link[!last]) && !is_red(s->link[last])) {
+ /* Color flip */
+ p->red = 0;
+ s->red = 1;
+ q->red = 1;
+ }
+ else {
+ int dir2 = g->link[1] == p;
+
+ if (is_red(s->link[last]))
+ g->link[dir2] = nbtree_double(p, last);
+ else if (is_red(s->link[!last]))
+ g->link[dir2] = nbtree_single(p, last);
+
+ /* Ensure correct coloring */
+ q->red = g->link[dir2]->red = 1;
+ g->link[dir2]->link[0]->red = 0;
+ g->link[dir2]->link[1]->red = 0;
+ }
+ }
+ }
+ }
+ }
+
+ /* Replace and remove if found */
+ if (f != NULL) {
+ f->data.id = q->data.id;
+ f->data.row = q->data.row;
+ f->data.col = q->data.col;
+ f->data.count = q->data.count;
+ memcpy(f->data.mean, q->data.mean, tree->datasize);
+ p->link[p->link[1] == q] = q->link[q->link[0] == NULL];
+
+ free(q->data.mean);
+ free(q);
+ q = NULL;
+ tree->count--;
+ removed = 1;
+ }
+ else
+ G_debug(2, "RB tree: data not found in search tree");
+
+ /* Update root and make it black */
+ tree->root = head.link[1];
+ if (tree->root != NULL)
+ tree->root->red = 0;
+
+ return removed;
+}
+
+/* find data item in tree
+ * returns pointer to data item if found else NULL
+ */
+struct ngbr_stats *nbtree_find(struct NB_TREE *tree, struct ngbr_stats *data)
+{
+ struct NB_NODE *curr_node = tree->root;
+ int cmp;
+
+ assert(tree && data);
+
+ while (curr_node != NULL) {
+ cmp = cmp_ngbr(&(curr_node->data), data);
+ if (cmp == 0)
+ return &curr_node->data; /* found */
+
+ curr_node = curr_node->link[cmp < 0];
+ }
+ return NULL;
+}
+
+/* initialize tree traversal
+ * (re-)sets trav structure
+ * returns 0
+ */
+int nbtree_init_trav(struct NB_TRAV *trav, struct NB_TREE *tree)
+{
+ assert(trav && tree);
+
+ trav->tree = tree;
+ trav->curr_node = tree->root;
+ trav->first = 1;
+ trav->top = 0;
+
+ return 0;
+}
+
+/* traverse the tree in ascending order
+ * useful to get all items in the tree non-recursively
+ * struct NB_TRAV *trav needs to be initialized first
+ * returns pointer to data, NULL when finished
+ */
+struct ngbr_stats *nbtree_traverse(struct NB_TRAV *trav)
+{
+ assert(trav);
+
+ if (trav->curr_node == NULL) {
+ if (trav->first)
+ G_debug(1, "RB tree: empty tree");
+ else
+ G_debug(1, "RB tree: finished traversing");
+
+ return NULL;
+ }
+
+ if (!trav->first)
+ return nbtree_next(trav);
+
+ trav->first = 0;
+ return nbtree_first(trav);
+}
+
+/* find start point to traverse the tree in ascending order
+ * useful to get a selection of items in the tree
+ * magnitudes faster than traversing the whole tree
+ * may return first item that's smaller or first item that's larger
+ * struct NB_TRAV *trav needs to be initialized first
+ * returns pointer to data, NULL when finished
+ */
+struct ngbr_stats *nbtree_traverse_start(struct NB_TRAV *trav, struct ngbr_stats *data)
+{
+ int dir = 0;
+
+ assert(trav && data);
+
+ if (trav->curr_node == NULL) {
+ if (trav->first)
+ G_warning("RB tree: empty tree");
+ else
+ G_warning("RB tree: finished traversing");
+
+ return NULL;
+ }
+
+ if (!trav->first)
+ return nbtree_next(trav);
+
+ /* else first time, get start node */
+
+ trav->first = 0;
+ trav->top = 0;
+
+ while (trav->curr_node != NULL) {
+ dir = cmp_ngbr(&(trav->curr_node->data), data);
+ /* exact match, great! */
+ if (dir == 0)
+ return &(trav->curr_node->data);
+ else {
+ dir = dir < 0;
+ /* end of branch, also reached if
+ * smallest item is larger than search template or
+ * largest item is smaller than search template */
+ if (trav->curr_node->link[dir] == NULL)
+ return &(trav->curr_node->data);
+
+ trav->up[trav->top++] = trav->curr_node;
+ trav->curr_node = trav->curr_node->link[dir];
+ }
+ }
+
+ return NULL; /* should not happen */
+}
+
+/* two functions needed to fully traverse the tree: initialize and continue
+ * useful to get all items in the tree non-recursively
+ * this one here uses a stack
+ * parent pointers or threads would also be possible
+ * but these would need to be added to NB_NODE
+ * -> more memory needed for standard operations
+ */
+
+/* start traversing the tree
+ * returns pointer to smallest data item
+ */
+static struct ngbr_stats *nbtree_first(struct NB_TRAV *trav)
+{
+ /* get smallest item */
+ while (trav->curr_node->link[0] != NULL) {
+ trav->up[trav->top++] = trav->curr_node;
+ trav->curr_node = trav->curr_node->link[0];
+ }
+
+ return &(trav->curr_node->data); /* return smallest item */
+}
+
+/* continue traversing the tree in ascending order
+ * returns pointer to data item, NULL when finished
+ */
+static struct ngbr_stats *nbtree_next(struct NB_TRAV *trav)
+{
+ if (trav->curr_node->link[1] != NULL) {
+ /* something on the right side: larger item */
+ trav->up[trav->top++] = trav->curr_node;
+ trav->curr_node = trav->curr_node->link[1];
+
+ /* go down, find smallest item in this branch */
+ while (trav->curr_node->link[0] != NULL) {
+ trav->up[trav->top++] = trav->curr_node;
+ trav->curr_node = trav->curr_node->link[0];
+ }
+ }
+ else {
+ /* at smallest item in this branch, go back up */
+ struct NB_NODE *last;
+
+ do {
+ if (trav->top == 0) {
+ trav->curr_node = NULL;
+ break;
+ }
+ last = trav->curr_node;
+ trav->curr_node = trav->up[--trav->top];
+ } while (last == trav->curr_node->link[1]);
+ }
+
+ if (trav->curr_node != NULL) {
+ return &(trav->curr_node->data);
+ }
+ else
+ return NULL; /* finished traversing */
+}
+
+/* destroy the tree */
+void nbtree_clear(struct NB_TREE *tree)
+{
+ struct NB_NODE *it;
+ struct NB_NODE *save = tree->root;
+
+ /*
+ Rotate away the left links so that
+ we can treat this like the destruction
+ of a linked list
+ */
+ while((it = save) != NULL) {
+ if (it->link[0] == NULL) {
+ /* No left links, just kill the node and move on */
+ save = it->link[1];
+ free(it->data.mean);
+ free(it);
+ it = NULL;
+ }
+ else {
+ /* Rotate away the left link and check again */
+ save = it->link[0];
+ it->link[0] = save->link[1];
+ save->link[1] = it;
+ }
+ }
+
+ tree->root = NULL;
+
+ tree->count = 0;
+
+ return;
+}
+
+/* used for debugging: check for errors in tree structure */
+int nbtree_debug(struct NB_TREE *tree, struct NB_NODE *root)
+{
+ int lh, rh;
+
+ if (root == NULL)
+ return 1;
+ else {
+ struct NB_NODE *ln = root->link[0];
+ struct NB_NODE *rn = root->link[1];
+ int lcmp = 0, rcmp = 0;
+
+ /* Consecutive red links */
+ if (is_red(root)) {
+ if (is_red(ln) || is_red(rn)) {
+ G_warning("Red Black Tree debugging: Red violation");
+ return 0;
+ }
+ }
+
+ lh = nbtree_debug(tree, ln);
+ rh = nbtree_debug(tree, rn);
+
+ if (ln) {
+ lcmp = cmp_ngbr(&(ln->data), &(root->data));
+ }
+
+ if (rn) {
+ rcmp = cmp_ngbr(&(rn->data), &(root->data));
+ }
+
+ /* Invalid binary search tree:
+ * left node >= parent or right node <= parent */
+ if ((ln != NULL && lcmp > -1)
+ || (rn != NULL && rcmp < 1)) {
+ G_warning("Red Black Tree debugging: Binary tree violation");
+ return 0;
+ }
+
+ /* Black height mismatch */
+ if (lh != 0 && rh != 0 && lh != rh) {
+ G_warning("Red Black Tree debugging: Black violation");
+ return 0;
+ }
+
+ /* Only count black links */
+ if (lh != 0 && rh != 0)
+ return is_red(root) ? lh : lh + 1;
+ else
+ return 0;
+ }
+}
+
+/*******************************************************
+ * *
+ * internal functions for Red Black Tree maintenance *
+ * *
+ *******************************************************/
+
+/* add a new node to the tree */
+static struct NB_NODE *nbtree_make_node(size_t datasize, struct ngbr_stats *data)
+{
+ struct NB_NODE *new_node = (struct NB_NODE *)malloc(sizeof(struct NB_NODE));
+
+ if (new_node == NULL)
+ G_fatal_error("RB Search Tree: Out of memory!");
+
+ if ((new_node->data.mean = malloc(datasize)) == NULL)
+ G_fatal_error("RB Search Tree: Out of memory!");
+
+ new_node->data.id = data->id;
+ new_node->data.row = data->row;
+ new_node->data.col = data->col;
+ new_node->data.count = data->count;
+ memcpy(new_node->data.mean, data->mean, datasize);
+
+ new_node->red = 1; /* 1 is red, 0 is black */
+ new_node->link[0] = NULL;
+ new_node->link[1] = NULL;
+
+ return new_node;
+}
+
+/* check for red violation */
+static int is_red(struct NB_NODE *root)
+{
+ if (root)
+ return root->red == 1;
+
+ return 0;
+}
+
+/* single rotation */
+static struct NB_NODE *nbtree_single(struct NB_NODE *root, int dir)
+{
+ struct NB_NODE *newroot = root->link[!dir];
+
+ root->link[!dir] = newroot->link[dir];
+ newroot->link[dir] = root;
+
+ root->red = 1;
+ newroot->red = 0;
+
+ return newroot;
+}
+
+/* double rotation */
+static struct NB_NODE *nbtree_double(struct NB_NODE *root, int dir)
+{
+ root->link[!dir] = nbtree_single(root->link[!dir], !dir);
+ return nbtree_single(root, dir);
+}
Property changes on: grass-addons/grass7/imagery/i.segment.xl/ngbrtree.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/imagery/i.segment.xl/ngbrtree.h
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/ngbrtree.h (rev 0)
+++ grass-addons/grass7/imagery/i.segment.xl/ngbrtree.h 2012-08-27 19:08:40 UTC (rev 52939)
@@ -0,0 +1,125 @@
+/*************************************************************
+ * USAGE *
+ *************************************************************
+ *
+ * NOTE: duplicates are not supported
+ *
+ * custom compare function
+ * extern int my_compare_fn(const void *, const void *);
+ * int my_compare_fn(const void *a, const void *b) {
+ * if ((mydatastruct *) a < (mydatastruct *) b)
+ * return -1;
+ * else if ((mydatastruct *) a > (mydatastruct *) b)
+ * return 1;
+ * else if ((mydatastruct *) a == (mydatastruct *) b)
+ * return 0;
+ * }
+ *
+ * create and initialize tree:
+ * struct NB_TREE *mytree = nbtree_create(my_compare_fn, item_size);
+ *
+ * insert items to tree:
+ * struct mydatastruct data = <some data>;
+ * if (nbtree_insert(mytree, &data) == 0)
+ * G_warning("could not insert data");
+ *
+ * find item in tree:
+ * struct mydatastruct data = <some data>;
+ * if (nbtree_find(mytree, &data) == 0)
+ * G_message("data not found");
+ *
+ * delete item from tree:
+ * struct mydatastruct data = <some data>;
+ * if (nbtree_remove(mytree, &data) == 0)
+ * G_warning("could not find data in tree");
+ *
+ * traverse tree (get all items in tree in ascending order):
+ * struct NB_TRAV trav;
+ * nbtree_init_trav(&trav, tree);
+ * while ((data = nbtree_traverse(&trav)) != NULL) {
+ * if (my_compare_fn(data, threshold_data) == 0) break;
+ * <do something with data>;
+ * }
+ *
+ * get a selection of items: all data > data1 and < data2
+ * start in tree where data is last smaller or first larger compared to data1
+ * struct NB_TRAV trav;
+ * nbtree_init_trav(&trav, tree);
+ * data = nbtree_traverse_start(&trav, &data1);
+ * <do something with data>;
+ * while ((data = nbtree_traverse(&trav)) != NULL) {
+ * if (data > data2) break;
+ * <do something with data>;
+ * }
+ *
+ * destroy tree:
+ * nbtree_destroy(mytree);
+ *
+ * debug the whole tree with
+ * nbtree_debug(mytree, mytree->root);
+ *
+ *************************************************************/
+
+#ifndef NBTREE_H
+#define NBTREE_H
+
+#include <stddef.h>
+
+/* maximum RB Tree height */
+#define NBTREE_MAX_HEIGHT 64 /* should be more than enough */
+
+/* routine to compare data items
+ * return -1 if rb_a < rb_b
+ * return 0 if rb_a == rb_b
+ * return 1 if rb_a > rb_b
+ */
+struct ngbr_stats
+{
+ int id; /* region ID */
+ int row, col; /* row, col of one cell in this region */
+ int count; /* number cells in this region */
+ double *mean; /* mean for each band = sum[b] / count */
+};
+
+
+struct NB_NODE
+{
+ unsigned char red; /* 0 = black, 1 = red */
+ struct NB_NODE *link[2]; /* link to children: link[0] for smaller, link[1] for larger */
+ struct ngbr_stats data; /* any kind of data */
+};
+
+struct NB_TREE
+{
+ struct NB_NODE *root; /* root node */
+ size_t datasize; /* item size */
+ size_t count; /* number of items in tree. */
+ int nbands; /* number of bands */
+};
+
+struct NB_TRAV
+{
+ struct NB_TREE *tree; /* tree being traversed */
+ struct NB_NODE *curr_node; /* current node */
+ struct NB_NODE *up[NBTREE_MAX_HEIGHT]; /* stack of parent nodes */
+ int top; /* index for stack */
+ int first; /* little helper flag */
+};
+
+
+/* tree functions */
+struct NB_TREE *nbtree_create(int, size_t);
+void nbtree_clear(struct NB_TREE *);
+int nbtree_insert(struct NB_TREE *, struct ngbr_stats *);
+int nbtree_remove(struct NB_TREE *, struct ngbr_stats *);
+struct ngbr_stats *nbtree_find(struct NB_TREE *, struct ngbr_stats *);
+
+/* tree traversal functions */
+int nbtree_init_trav(struct NB_TRAV *, struct NB_TREE *);
+struct ngbr_stats *nbtree_traverse(struct NB_TRAV *);
+struct ngbr_stats *nbtree_traverse_start(struct NB_TRAV *, struct ngbr_stats *);
+
+/* debug tree from given node downwards */
+int nbtree_debug(struct NB_TREE *, struct NB_NODE *);
+
+#endif
Property changes on: grass-addons/grass7/imagery/i.segment.xl/ngbrtree.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Modified: grass-addons/grass7/imagery/i.segment.xl/open_files.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/open_files.c 2012-08-27 18:43:23 UTC (rev 52938)
+++ grass-addons/grass7/imagery/i.segment.xl/open_files.c 2012-08-27 19:08:40 UTC (rev 52939)
@@ -1,6 +1,5 @@
/* PURPOSE: opening input rasters and creating segmentation files */
-#include <limits.h> /* for INT_MAX */
#include <stdlib.h>
#include <grass/gis.h>
#include <grass/glocale.h>
@@ -8,65 +7,38 @@
#include <grass/segment.h> /* segmentation library */
#include "iseg.h"
-int open_files(struct files *files, struct functions *functions)
+static int load_seeds(struct globals *, int, int, int);
+static int read_seed(struct globals *, SEGMENT *, struct rc *, int);
+
+int open_files(struct globals *globals)
{
struct Ref Ref; /* group reference list */
- int *in_fd, seeds_fd, bounds_fd, null_check, out_fd, mean_fd;
- int n, s, row, col, srows, scols, inlen, nseg, borderPixels;
- DCELL **inbuf; /* buffer array, to store lines from each of the imagery group rasters */
- CELL *boundsbuf, *seedsbuf;
- void *ptr; /* for iterating seedsbuf */
- size_t ptrsize;
- struct FPRange *fp_range; /* for getting min/max values on each input raster */
+ int *in_fd, bounds_fd, is_null;
+ int n, row, col, srows, scols, inlen, outlen, nseg;
+ DCELL **inbuf; /* buffers to store lines from each of the imagery group rasters */
+ CELL *boundsbuf, bounds_null, bounds_val;
+ int have_bounds = 0;
+ CELL s, id;
+ struct Range range; /* min/max values of bounds map */
+ struct FPRange *fp_range; /* min/max values of each input raster */
DCELL *min, *max;
- struct Range range; /* for seeds range */
- int seeds_min, seeds_max;
+ struct ngbr_stats Ri, Rk;
-
- /* for merging seed values */
- struct pixels *R_head, *Rn_head, *newpixel, *current;
- int R_count;
-
- G_verbose_message(_("Opening files and initializing"));
-
- /* confirm output maps can be opened (don't want to do all this work for nothing!) */
- out_fd = Rast_open_new(files->out_name, CELL_TYPE);
- if (out_fd < 0)
- G_fatal_error(_("Could not open output raster for writing segment ID's"));
- else
- Rast_unopen(out_fd);
-
- if (files->out_band != NULL) {
- mean_fd = Rast_open_new(files->out_band, DCELL_TYPE);
- if (mean_fd < 0)
- G_fatal_error(_("Could not open output raster for writing mean segment values"));
- else
- Rast_unopen(mean_fd);
- }
-
/*allocate memory for flags */
- files->null_flag = flag_create(files->nrows, files->ncols);
- flag_clear_all(files->null_flag);
- files->candidate_flag = flag_create(files->nrows, files->ncols);
+ globals->null_flag = flag_create(globals->nrows, globals->ncols);
+ globals->candidate_flag = flag_create(globals->nrows, globals->ncols);
- if (files->bounds_map != NULL)
- files->orig_null_flag = flag_create(files->nrows, files->ncols);
+ G_debug(1, "Checking image group...");
- files->seeds_flag = flag_create(files->nrows, files->ncols);
- flag_clear_all(files->seeds_flag);
-
- /* references for segmentation library: i.cost r.watershed/seg and http://grass.osgeo.org/programming7/segmentlib.html */
-
/* ****** open the input rasters ******* */
- /* Note: I confirmed, the API does not check this: */
- if (!I_get_group_ref(files->image_group, &Ref))
+ if (!I_get_group_ref(globals->image_group, &Ref))
G_fatal_error(_("Unable to read REF file for group <%s>"),
- files->image_group);
+ globals->image_group);
if (Ref.nfiles <= 0)
G_fatal_error(_("Group <%s> contains no raster maps"),
- files->image_group);
+ globals->image_group);
/* Read Imagery Group */
@@ -76,272 +48,421 @@
min = G_malloc(Ref.nfiles * sizeof(DCELL));
max = G_malloc(Ref.nfiles * sizeof(DCELL));
+ G_debug(1, "Opening input rasters...");
for (n = 0; n < Ref.nfiles; n++) {
inbuf[n] = Rast_allocate_d_buf();
in_fd[n] = Rast_open_old(Ref.file[n].name, Ref.file[n].mapset);
- if (in_fd[n] < 0)
- G_fatal_error(_("Error opening %s@%s"), Ref.file[n].name,
- Ref.file[n].mapset);
}
- /* open seeds raster and confirm all positive integers were given */
- if (files->seeds_map != NULL) {
- seeds_fd = Rast_open_old(files->seeds_map, "");
- seedsbuf = Rast_allocate_c_buf();
- ptrsize = sizeof(CELL);
+ /* Get min/max values of each input raster for scaling */
- if (Rast_read_range(files->seeds_map, files->seeds_mapset, &range) != 1) { /* returns -1 on error, 2 on empty range, quiting either way. */
- G_fatal_error(_("No min/max found in seeds raster map <%s>"),
- files->seeds_map);
- }
- Rast_get_range_min_max(&range, &seeds_min, &seeds_max);
- if (seeds_min < 0)
- G_fatal_error(_("Seeds raster should have postive integers for starting seeds, and zero or NULL for all other pixels."));
- }
+ globals->threshold = 0.;
- /* Get min/max values of each input raster for scaling */
+ for (n = 0; n < Ref.nfiles; n++) {
+ /* returns -1 on error, 2 on empty range, quiting either way. */
+ if (Rast_read_fp_range(Ref.file[n].name, Ref.file[n].mapset, &fp_range[n]) != 1)
+ G_fatal_error(_("No min/max found in raster map <%s>"),
+ Ref.file[n].name);
+ Rast_get_fp_range_min_max(&(fp_range[n]), &min[n], &max[n]);
- if (files->weighted == FALSE) { /*default, we will scale */
- for (n = 0; n < Ref.nfiles; n++) {
- if (Rast_read_fp_range(Ref.file[n].name, Ref.file[n].mapset, &fp_range[n]) != 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[n]), &min[n], &max[n]);
- }
+ G_debug(1, "Range for layer %d: min = %f, max = %f",
+ n, min[n], max[n]);
+
}
+ if (globals->weighted == FALSE)
+ globals->threshold = Ref.nfiles;
+ else {
+ /* max difference with selected similarity method */
+ Ri.mean = max;
+ Rk.mean = min;
+ globals->threshold = (*globals->calculate_similarity) (&Ri, &Rk, globals);
+ }
/* ********** find out file segmentation size ************ */
+ G_debug(1, "Calculate temp file sizes...");
- files->nbands = Ref.nfiles;
+ globals->nbands = Ref.nfiles;
/* size of each element to be stored */
- inlen = sizeof(double) * Ref.nfiles;
+ inlen = sizeof(DCELL) * Ref.nfiles;
+ outlen = sizeof(CELL);
- /* when fine tuning, should be a power of 2 and not larger than 256 for speed reasons */
+ /* segment lib segment size */
srows = 64;
scols = 64;
- /* RAM enhancement: have user input limit and make calculations for this, reference i.cost and i.watershed
- * One segment tile is tile_mb = (nbands * sizeof(double) + sizeof(CELL) * srows * scols / (1024*1024)
- * (check if sizeof(CELL) was from when iseg would be included, or the extra overhead?)
- * If user inputs total RAM available, need to subtract the size of the flags, and the size of the linked lists.
- * I'm not sure how to estimate the size of the linked lists, it is allowed to grow when needed. Assume one segment
- * could be 50% of the map? Or more? So ll_mb = sizeof(pixels) * nrows * ncols * 0.5 / (1024*1024)
- * then split the remaining RAM between bands_seg and iseg_seg. */
+ /* calculate number of segments in memory */
+ nseg = (1024. * 1024. * globals->mb) /
+ (sizeof(DCELL) * Ref.nfiles * srows * scols +
+ + sizeof(CELL) * 2 * srows * scols);
- nseg = 16;
-
-
- /* ******* create temporary segmentation files ********* */
- /* Initalize access to database and create temporary files */
-
- G_debug(1, "Image size: %d rows, %d cols", files->nrows, files->ncols);
- G_debug(1, "Segmented to tiles with size: %d rows, %d cols", srows,
+ /* create segment structures */
+ G_debug(1, "image size: %d rows, %d cols", globals->nrows, globals->ncols);
+ G_debug(1, "degmented to tiles with size: %d rows, %d cols", srows,
scols);
- G_debug(1, "Data element size, in: %d", inlen);
- G_debug(1, "number of segments to have in memory: %d", nseg);
+ G_debug(1, "data element size, in: %d , out: %d ", inlen, outlen);
+ G_debug(1, "number of segments in memory: %d", nseg);
if (segment_open
- (&files->bands_seg, G_tempfile(), files->nrows, files->ncols, srows,
- scols, inlen, nseg) != TRUE)
+ (&globals->bands_seg, G_tempfile(), globals->nrows, globals->ncols, srows,
+ scols, inlen, nseg) != 1)
G_fatal_error("Unable to create input temporary files");
- /* ******* remaining memory allocation ********* */
+ if (segment_open
+ (&globals->rid_seg, G_tempfile(), globals->nrows, globals->ncols, srows,
+ scols, outlen, nseg * 2) != 1)
+ G_fatal_error("Unable to create input temporary files");
- /* 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;
+ /* load input bands to segment structure */
+ G_debug(1, "Reading input rasters into segmentation data files...");
- files->bands_val = (double *)G_malloc(inlen);
- files->second_val = (double *)G_malloc(inlen);
+ globals->bands_val = (double *)G_malloc(inlen);
+ globals->second_val = (double *)G_malloc(inlen);
+ /* initial segment ID */
+ s = 1;
- if (segment_open
- (&files->iseg_seg, G_tempfile(), files->nrows, files->ncols, srows,
- scols, sizeof(int), nseg) != TRUE)
- G_fatal_error(_("Unable to allocate memory for initial segment ID's"));
+ globals->row_min = globals->nrows;
+ globals->row_max = 0;
+ globals->col_min = globals->ncols;
+ globals->col_max = 0;
+ for (row = 0; row < globals->nrows; row++) {
+ for (n = 0; n < Ref.nfiles; n++) {
+ Rast_get_d_row(in_fd[n], inbuf[n], row);
+ }
+ for (col = 0; col < globals->ncols; col++) {
- /* bounds/constraints (start with processing constraints to get any possible NULL values) */
- if (files->bounds_map != NULL) {
+ is_null = 0; /*Assume there is data */
+ for (n = 0; n < Ref.nfiles; n++) {
+ if (Rast_is_d_null_value(&inbuf[n][col])) {
+ is_null = 1;
+ globals->bands_val[n] = inbuf[n][col];
+ }
+ else {
+ if (globals->weighted == TRUE)
+ globals->bands_val[n] = inbuf[n][col];
+ else
+ /* scaled version */
+ globals->bands_val[n] = (inbuf[n][col] - min[n]) / (max[n] - min[n]);
+ }
+ }
+ segment_put(&globals->bands_seg, (void *)globals->bands_val, row, col);
+
+ if (!is_null) {
+
+ FLAG_UNSET(globals->null_flag, row, col);
+
+ if (!globals->seeds) {
+ /* sequentially number all cells with a unique segment ID */
+ id = s;
+ s++;
+ }
+
+ /* get min/max row/col to narrow the processing window */
+
+ if (globals->row_min > row)
+ globals->row_min = row;
+ if (globals->row_max < row)
+ globals->row_max = row;
+ if (globals->col_min > col)
+ globals->col_min = col;
+ if (globals->col_max < col)
+ globals->col_max = col;
+ }
+ else {
+ /* all input bands NULL */
+ Rast_set_c_null_value(&id, 1);
+ /*Rast_set_c_null_value(&(globals->rid[row][col]), 1);*/
+ FLAG_SET(globals->null_flag, row, col);
+ }
+ if (!globals->seeds)
+ segment_put(&globals->rid_seg, (void *)&id, row, col);
+ }
+ }
+ G_debug(1, "nrows: %d, min row: %d, max row %d", globals->nrows, globals->row_min, globals->row_max);
+ G_debug(1, "ncols: %d, min col: %d, max col %d", globals->ncols, globals->col_min, globals->col_max);
+
+ globals->row_max++;
+ globals->col_max++;
+ globals->ncells = (globals->row_max - globals->row_min) * (globals->col_max - globals->col_min);
+
+ /* bounds/constraints */
+
+ Rast_set_c_null_value(&globals->upper_bound, 1);
+ Rast_set_c_null_value(&globals->lower_bound, 1);
+
+ if (globals->bounds_map != NULL) {
if (segment_open
- (&files->bounds_seg, G_tempfile(), files->nrows, files->ncols,
- srows, scols, sizeof(int), nseg) != TRUE)
- G_fatal_error(_("Unable to create bounds temporary files"));
+ (&globals->bounds_seg, G_tempfile(), globals->nrows, globals->ncols,
+ srows, scols, sizeof(CELL), nseg) != TRUE)
+ G_fatal_error("Unable to create bounds temporary files");
+ if (Rast_read_range(globals->bounds_map, globals->bounds_mapset, &range) != 1)
+ G_fatal_error(_("No min/max found in raster map <%s>"),
+ globals->bounds_map);
+ Rast_get_range_min_max(&range, &globals->upper_bound,
+ &globals->lower_bound);
+
+ if (Rast_is_c_null_value(&globals->upper_bound) ||
+ Rast_is_c_null_value(&globals->lower_bound)) {
+
+ G_fatal_error(_("No min/max found in raster map <%s>"),
+ globals->bounds_map);
+ }
+
+ bounds_null = globals->upper_bound = globals->lower_bound + 1;
+
+ bounds_fd = Rast_open_old(globals->bounds_map, globals->bounds_mapset);
boundsbuf = Rast_allocate_c_buf();
- bounds_fd = Rast_open_old(files->bounds_map, files->bounds_mapset);
- for (row = 0; row < files->nrows; row++) {
+ for (row = 0; row < globals->nrows; row++) {
Rast_get_c_row(bounds_fd, boundsbuf, row);
- for (col = 0; col < files->ncols; col++) {
- files->bounds_val = boundsbuf[col];
- segment_put(&files->bounds_seg, &files->bounds_val, row, col);
- if (Rast_is_c_null_value(&boundsbuf[col]) == TRUE) {
- FLAG_SET(files->null_flag, row, col);
+ for (col = 0; col < globals->ncols; col++) {
+ if (FLAG_GET(globals->null_flag, row, col)) {
+ Rast_set_c_null_value(&bounds_val, 1);
}
+ else {
+ if (Rast_is_c_null_value(&boundsbuf[col]))
+ boundsbuf[col] = bounds_null;
+ else {
+ have_bounds = 1;
+ if (globals->lower_bound > boundsbuf[col])
+ globals->lower_bound = boundsbuf[col];
+ }
+ }
+ segment_put(&globals->bounds_seg, &bounds_val, row, col);
}
}
Rast_close(bounds_fd);
G_free(boundsbuf);
- } /* end bounds/constraints opening */
+ if (!have_bounds) {
+ G_warning(_("There are no boundary constraints in '%s'"), globals->bounds_map);
+ Rast_set_c_null_value(&globals->upper_bound, 1);
+ Rast_set_c_null_value(&globals->lower_bound, 1);
+ segment_close(&globals->bounds_seg);
+ globals->bounds_map = NULL;
+ globals->bounds_mapset = NULL;
+ }
+ }
+ else {
+ G_debug(1, "no boundary constraint supplied.");
+ }
- /* ******** load input bands to segment structure and fill initial seg ID's ******** */
- s = 0; /* initial segment ID will be 1 */
+ /* other info */
+ globals->candidate_count = 0; /* counter for remaining candidate pixels */
- for (row = 0; row < files->nrows; row++) {
+ /* Free memory */
- /* read in rows of data (each input band from the imagery group and the optional seeds map) */
- for (n = 0; n < Ref.nfiles; n++) {
- Rast_get_d_row(in_fd[n], inbuf[n], row);
- }
- if (files->seeds_map != NULL) {
- Rast_get_c_row(seeds_fd, seedsbuf, row);
- ptr = seedsbuf;
- }
+ for (n = 0; n < Ref.nfiles; n++) {
+ G_free(inbuf[n]);
+ Rast_close(in_fd[n]);
+ }
- for (col = 0; col < files->ncols; col++) {
- if (FLAG_GET(files->null_flag, row, col))
- continue;
- null_check = 1; /*Assume there is data */
- for (n = 0; n < Ref.nfiles; n++) {
- if (Rast_is_d_null_value(&inbuf[n][col]))
- null_check = -1;
- if (files->weighted == TRUE)
- files->bands_val[n] = inbuf[n][col]; /*unscaled */
- else
- files->bands_val[n] = (inbuf[n][col] - min[n]) / (max[n] - min[n]); /* scaled */
- }
- 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? */
- segment_put(&files->bands_seg, (void *)files->bands_val, row, col); /* store input bands */
+ globals->datasize = sizeof(double) * globals->nbands;
+ globals->rs.sum = G_malloc(globals->datasize);
+ globals->rs.mean = G_malloc(globals->datasize);
- if (null_check != -1) { /*good pixel */
- FLAG_UNSET(files->null_flag, row, col); /*flag */
- if (files->seeds_map != NULL) {
- if (Rast_is_c_null_value(ptr) == TRUE) {
- /* when using iseg_seg the segmentation file is already initialized to zero. Just initialize seeds_flag: */
- FLAG_UNSET(files->seeds_flag, row, col);
- }
- else {
- FLAG_SET(files->seeds_flag, row, col); /* RAM enhancement, but it might cost speed. Could look for seg ID > 0 instead of using seed_flag. */
- /* seed value is starting segment ID. */
- segment_put(&files->iseg_seg, ptr, row, col);
- }
- ptr = G_incr_void_ptr(ptr, ptrsize);
- }
- else { /* no seeds provided */
- s++; /* sequentially number all pixels with their own segment ID */
- if (s < INT_MAX) { /* Check that the starting seeds aren't too large. */
- segment_put(&files->iseg_seg, &s, row, col); /*starting segment number */
- FLAG_SET(files->seeds_flag, row, col); /*all pixels are seeds */
- }
- else
- G_fatal_error(_("Exceeded integer storage limit, too many initial pixels."));
- }
+ globals->reg_tree = rgtree_create(globals->nbands, globals->datasize);
+ globals->n_regions = s - 1;
+
+ if (globals->seeds) {
+ load_seeds(globals, srows, scols, nseg);
+ }
+
+ G_free(inbuf);
+ G_free(in_fd);
+ G_free(fp_range);
+ G_free(min);
+ G_free(max);
+ /* Need to clean up anything else? */
+
+ return TRUE;
+}
+
+
+static int load_seeds(struct globals *globals, int srows, int scols, int nseg)
+{
+ int row, col;
+ SEGMENT seeds_seg;
+ CELL *seeds_buf, seeds_val;
+ int seeds_fd;
+ int spos, sneg, have_seeds;
+ struct rc Ri;
+
+ G_debug(1, "load_seeds()");
+
+ if (segment_open
+ (&seeds_seg, G_tempfile(), globals->nrows, globals->ncols,
+ srows, scols, sizeof(CELL), nseg) != TRUE)
+ G_fatal_error("Unable to create bounds temporary files");
+
+ seeds_fd = Rast_open_old(globals->seeds, "");
+ seeds_buf = Rast_allocate_c_buf();
+
+ have_seeds = 0;
+
+ /* load seeds map to segment structure */
+ for (row = 0; row < globals->nrows; row++) {
+ Rast_get_c_row(seeds_fd, seeds_buf, row);
+ for (col = 0; col < globals->ncols; col++) {
+ if (FLAG_GET(globals->null_flag, row, col)) {
+ Rast_set_c_null_value(&seeds_val, 1);
}
- else { /*don't use this pixel */
- FLAG_SET(files->null_flag, row, col); /*flag */
+ else {
+ seeds_val = seeds_buf[col];
+ if (!Rast_is_c_null_value(&seeds_buf[col]))
+ have_seeds = 1;
}
+ segment_put(&seeds_seg, &seeds_val, row, col);
}
}
- /* keep original copy of null flag if we have boundary constraints */
- if (files->bounds_map != NULL) {
- for (row = 0; row < files->nrows; row++) {
- for (col = 0; col < files->ncols; col++) {
- if (FLAG_GET(files->null_flag, row, col))
- FLAG_SET(files->orig_null_flag, row, col);
- else
- FLAG_UNSET(files->orig_null_flag, row, col);
+ if (!have_seeds) {
+ G_warning(_("No seeds found in '%s'!"), globals->seeds);
+ G_free(seeds_buf);
+ Rast_close(seeds_fd);
+ segment_close(&seeds_seg);
+ return 0;
+ }
+
+ spos = 1;
+ sneg = -1;
+
+ /* convert seeds to regions */
+ G_debug(1, "convert seeds to regions");
+ Rast_set_c_null_value(&seeds_val, 1);
+ for (row = 0; row < globals->nrows; row++) {
+ Rast_get_c_row(seeds_fd, seeds_buf, row);
+ for (col = 0; col < globals->ncols; col++) {
+ if (FLAG_GET(globals->null_flag, row, col)) {
+ segment_put(&globals->rid_seg, &seeds_val, row, col);
}
+ else if (!(FLAG_GET(globals->candidate_flag, row, col))) {
+ if (Rast_is_c_null_value(&(seeds_buf[col])) || seeds_buf[col] < 1) {
+ segment_put(&globals->rid_seg, &sneg, row, col);
+ sneg--;
+ }
+ else {
+ Ri.row = row;
+ Ri.col = col;
+ read_seed(globals, &seeds_seg, &Ri, spos);
+ spos++;
+ }
+ }
}
- } /* end: if (files->bounds_map != NULL) */
+ }
- /* linked list memory management linkm */
- link_set_chunk_size(1000); /* TODO RAM: fine tune this number */
- files->token = link_init(sizeof(struct pixels));
+ G_free(seeds_buf);
+ Rast_close(seeds_fd);
+ segment_close(&seeds_seg);
- /* if we have seeds that are segments (not pixels) we need to update the bands_seg */
- /* also renumber the segment ID's in case they were classified (duplicating numbers) instead of output from i.segment. */
- if (files->seeds_map != NULL) {
+ globals->n_regions = spos - 1;
+
+ flag_clear_all(globals->candidate_flag);
+
+ return 1;
+}
- /*initialization */
- files->minrow = files->mincol = 0;
- files->maxrow = files->nrows;
- files->maxcol = files->ncols;
- R_count = 1;
- R_head = NULL;
- Rn_head = NULL;
- newpixel = NULL;
- current = NULL;
- set_all_candidate_flags(files);
- for (row = 0; row < files->nrows; row++) {
- G_percent(row, files->nrows, 1); /* I think this is the longest part of open_files() - not entirely accurate for the actual %, but will give the user something to see. */
- for (col = 0; col < files->ncols; col++) {
- if (!(FLAG_GET(files->candidate_flag, row, col)) ||
- FLAG_GET(files->null_flag, row, col))
- continue;
- /*start R_head */
- newpixel = (struct pixels *)link_new(files->token);
- newpixel->next = NULL;
- newpixel->row = row;
- newpixel->col = col;
- R_head = newpixel;
- /* get pixel list, possible initialization speed enhancement: could use a custom (shorter) function, some results from find_segment_neighbors are not used here */
- /* bug todo: There is a small chance that a renumbered segment matches and borders an original segment. This would be a good reason to write a custom function - use the candidate flag to see if the pixel was already processed. */
- borderPixels =
- find_segment_neighbors(&R_head, &Rn_head, &R_count, files,
- functions);
+static int read_seed(struct globals *globals, SEGMENT *seeds_seg, struct rc *Ri, int new_id)
+{
+ int n, i, Ri_id, Rk_id;
+ struct rc ngbr_rc, next;
+ struct rclist rilist;
+ int neighbors[8][2];
- /* update the segment ID */
+ G_debug(4, "read_seed()");
- s++;
- if (s == INT_MAX) /* Check that the starting seeds aren't too large. */
- G_fatal_error(_("Exceeded integer storage limit, too many initial pixels."));
+ /* get Ri's segment ID from input seeds */
+ segment_get(seeds_seg, &Ri_id, Ri->row, Ri->col);
+
+ /* set new segment id */
+ segment_put(&globals->rid_seg, &new_id, Ri->row, Ri->col);
+ /* set candidate flag */
+ FLAG_SET(globals->candidate_flag, Ri->row, Ri->col);
- for (current = R_head; current != NULL;
- current = current->next) {
- segment_put(&files->iseg_seg, &s, current->row,
- current->col);
- FLAG_UNSET(files->candidate_flag, current->row,
- current->col);
- }
+ /* initialize region stats */
+ globals->rs.count = 1;
- /*merge pixels (updates the bands_seg) */
- merge_pixels(R_head, borderPixels, files);
+ globals->rs.id = new_id;
+ segment_get(&globals->bands_seg, (void *)globals->bands_val,
+ Ri->row, Ri->col);
- /*todo calculate perimeter (?and area?) here? */
+ for (i = 0; i < globals->nbands; i++) {
+ globals->rs.sum[i] = globals->bands_val[i];
+ globals->rs.mean[i] = globals->bands_val[i];
+ }
- /*clean up */
- my_dispose_list(files->token, &R_head);
- my_dispose_list(files->token, &Rn_head);
- R_count = 1;
+ /* go through seed, spreading outwards from head */
+ rclist_init(&rilist);
+ rclist_add(&rilist, Ri->row, Ri->col);
+
+ while (rclist_drop(&rilist, &next)) {
+
+ G_debug(5, "find_pixel_neighbors for row: %d , col %d",
+ next.row, next.col);
+
+ globals->find_neighbors(next.row, next.col, neighbors);
+
+ for (n = 0; n < globals->nn; n++) {
+
+ ngbr_rc.row = neighbors[n][0];
+ ngbr_rc.col = neighbors[n][1];
+
+ if (ngbr_rc.row < 0 || ngbr_rc.row >= globals->nrows ||
+ ngbr_rc.col < 0 || ngbr_rc.col >= globals->ncols) {
+ continue;
}
- }
- }
- files->nsegs = s;
+ if (FLAG_GET(globals->null_flag, ngbr_rc.row, ngbr_rc.col)) {
+ continue;
+ }
- /* Free memory */
+ if (FLAG_GET(globals->candidate_flag, ngbr_rc.row, ngbr_rc.col)) {
+ continue;
+ }
- for (n = 0; n < Ref.nfiles; n++) {
- G_free(inbuf[n]);
- Rast_close(in_fd[n]);
+ segment_get(seeds_seg, (void *) &Rk_id, ngbr_rc.row, ngbr_rc.col);
+
+ G_debug(5, "Rk ID = %d Ri ID = %d", Rk_id, Ri_id);
+
+ if (Rk_id != Ri_id) {
+ continue;
+ }
+
+ /* set segment id */
+ segment_put(&globals->rid_seg, &new_id, ngbr_rc.row, ngbr_rc.col);
+
+ /* set candidate flag */
+ FLAG_SET(globals->candidate_flag, ngbr_rc.row, ngbr_rc.col);
+
+ /* add to list of cells to check */
+ rclist_add(&rilist, ngbr_rc.row, ngbr_rc.col);
+
+ /* update region stats */
+ segment_get(&globals->bands_seg, (void *)globals->bands_val,
+ ngbr_rc.row, ngbr_rc.col);
+
+ for (i = 0; i < globals->nbands; i++) {
+ globals->rs.sum[i] += globals->bands_val[i];
+ }
+ globals->rs.count++;
+ }
}
- if (files->seeds_map != NULL) {
- Rast_close(seeds_fd);
- G_free(seedsbuf);
+ if (rgtree_find(globals->reg_tree, &(globals->rs)) != NULL) {
+ G_fatal_error(_("Segment %d is already registered!"), new_id);
}
+
+ /* insert into region tree */
+ if (globals->rs.count > 1) {
+ for (i = 0; i < globals->nbands; i++)
+ globals->rs.mean[i] = globals->rs.sum[i] / globals->rs.count;
- G_free(inbuf);
- G_free(in_fd);
- G_free(fp_range);
- G_free(min);
- G_free(max);
+ rgtree_insert(globals->reg_tree, &(globals->rs));
+ }
- return TRUE;
+ return 1;
}
Modified: grass-addons/grass7/imagery/i.segment.xl/parse_args.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/parse_args.c 2012-08-27 18:43:23 UTC (rev 52938)
+++ grass-addons/grass7/imagery/i.segment.xl/parse_args.c 2012-08-27 19:08:40 UTC (rev 52939)
@@ -7,66 +7,45 @@
#include <grass/raster.h>
#include "iseg.h"
-int parse_args(int argc, char *argv[], struct files *files,
- struct functions *functions)
+int parse_args(int argc, char *argv[], struct globals *globals)
{
- 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 *outband; /* optional saving of segment data, until a seperate module is written */
+ struct Option *group, *seeds, *bounds, *output,
+ *method, *threshold, *min_segment_size, *mem;
+ struct Flag *diagonal, *weighted;
+ struct Option *outband, *endt;
-#ifdef VCLOSE
- struct Option *very_close;
-#endif
-
/* required parameters */
- group = G_define_standard_option(G_OPT_I_GROUP); /* enhancement: consider giving the option to process just one raster directly, without creating an image group. */
+ /* 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->description = _("Similarity threshold.");
+ threshold->label = _("Difference threshold between 0 and 1.");
+ threshold->description = _("Threshold = 0 would merge only identical segments; threshold = 1 would merge all.");
method = G_define_option();
method->key = "method";
method->type = TYPE_STRING;
method->required = YES;
method->answer = "region_growing";
- method->options = "region_growing";
method->description = _("Segmentation method.");
- similarity = G_define_option();
- similarity->key = "similarity";
- similarity->type = TYPE_STRING;
- similarity->required = YES;
- similarity->answer = "euclidean";
- similarity->options = "euclidean, manhattan";
- similarity->description = _("Distance calculation method.");
-
min_segment_size = G_define_option();
- min_segment_size->key = "minsize";
+ min_segment_size->key = "min";
min_segment_size->type = TYPE_INTEGER;
min_segment_size->required = YES;
- min_segment_size->answer = "1"; /* default: no merges, a minimum of 1 pixel is allowed in a segment. */
+ 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 iteration will ignore the threshold for any segments with fewer pixels.");
+ _("The final step will merge small segments with their best neighbor.");
-#ifdef VCLOSE
- very_close = G_define_option();
- very_close->key = "very_close";
- very_close->type = TYPE_DOUBLE;
- very_close->required = YES;
- very_close->answer = "0";
- very_close->options = "0-1";
- very_close->label = _("Fraction of the threshold.");
- very_close->description =
- _("Neighbors similarity lower then this fraction of the threshold will be merged without regard to any other processing rules.");
-#endif
-
/* optional parameters */
diagonal = G_define_flag();
@@ -79,13 +58,12 @@
weighted->description =
_("Weighted input, don't perform the default scaling of input maps.");
- /* Raster for initial segment seeds *//* future enhancement: allow vector points/centroids for seed input. */
+ /* Using raster for seeds
+ * Low priority TODO: allow vector points/centroids seed input. */
seeds = G_define_standard_option(G_OPT_R_INPUT);
seeds->key = "seeds";
seeds->required = NO;
- seeds->label = _("Optional raster map with starting seeds.");
- seeds->description =
- _("Pixel values with positive integers are used as starting seeds.");
+ seeds->description = _("Optional raster map with starting seeds.");
/* Polygon constraints. */
bounds = G_define_standard_option(G_OPT_R_INPUT);
@@ -93,150 +71,132 @@
bounds->required = NO;
bounds->label = _("Optional bounding/constraining raster map");
bounds->description =
- _("Pixels with the same integer value will be segmented independent of the others.");
+ _("Must be integer values, each area will be segmented independent of the others.");
- /* other parameters */
+ mem = G_define_option();
+ mem->key = "memory";
+ mem->type = TYPE_INTEGER;
+ mem->required = NO;
+ mem->answer = "300";
+ mem->description = _("Memory in MB.");
+
+ /* TODO input for distance function */
+
+ /* debug parameters */
endt = G_define_option();
- endt->key = "endt";
+ endt->key = "iterations";
endt->type = TYPE_INTEGER;
endt->required = NO;
endt->answer = "1000";
- endt->description =
- _("Maximum number of passes (time steps) to complete.");
+ endt->description = _("Maximum number of iterations.");
- /* Leaving path flag out of user options, will hardcode TRUE for this option. This does change the resulting segments, but I don't see any
- * reason why one version is more valid. It reduced the processing time by 50% when segmenting the ortho image.
- * Code in the segmenation algorithm remains, in case more validation of this option should be done. */
- /*
- path = G_define_flag();
- path->key = 'p';
- path->description =
- _("temporary option, pathflag, select to use Rk as next Ri if not mutually best neighbors.");
- */
- limited = G_define_flag();
- limited->key = 'l';
- limited->description =
- _("segments are limited to be included in only one merge per pass");
-
outband = G_define_standard_option(G_OPT_R_OUTPUT);
- outband->key = "final_mean";
+ outband->key = "goodness";
outband->required = NO;
outband->description =
- _("Save the final mean values for the first band in the imagery group.");
- /* enhancement: save mean values for all bands. Multiple rasters or switch to polygons? */
+ _("debug - save band mean, currently implemented for only 1 band.");
-
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
- /* Validation */
-
/* Check and save parameters */
- files->image_group = group->answer;
+ globals->image_group = group->answer;
if (G_legal_filename(output->answer) == TRUE)
- files->out_name = output->answer; /* name of output (segment ID) raster map */
+ globals->out_name = output->answer;
else
- G_fatal_error(_("Invalid output raster name."));
+ G_fatal_error("Invalid output raster name.");
- functions->threshold = atof(threshold->answer); /* Note: this threshold is scaled after we know more at the beginning of create_isegs() */
+ /* Note: this threshold is scaled after we know more at the beginning of create_isegs() */
+ globals->alpha = atof(threshold->answer);
- if (weighted->answer == FALSE &&
- (functions->threshold <= 0 || functions->threshold >= 1))
+ if (globals->threshold <= 0 || globals->threshold >= 1)
G_fatal_error(_("threshold should be >= 0 and <= 1"));
- /* segmentation methods: 1 = region growing */
+ /* segmentation methods: 1 = region growing */
if (strncmp(method->answer, "region_growing", 10) == 0)
- functions->method = 1;
+ globals->method = 1;
else
- G_fatal_error(_("Couldn't assign segmentation method.")); /*shouldn't be able to get here */
+ G_fatal_error("Couldn't assign segmentation method.");
- /* distance methods for similarity measurement */
- if (strncmp(similarity->answer, "euclidean", 5) == 0)
- functions->calculate_similarity = &calculate_euclidean_similarity;
- else if (strncmp(similarity->answer, "manhattan", 5) == 0)
- functions->calculate_similarity = &calculate_manhattan_similarity;
- else
- G_fatal_error(_("Couldn't assign similarity method.")); /*shouldn't be able to get here */
+ G_debug(1, "segmentation method: %d", globals->method);
-#ifdef VCLOSE
- functions->very_close = atof(very_close->answer);
-#endif
+ globals->min_segment_size = atoi(min_segment_size->answer);
- functions->min_segment_size = atoi(min_segment_size->answer);
-
if (diagonal->answer == FALSE) {
- functions->find_pixel_neighbors = &find_four_pixel_neighbors;
- functions->num_pn = 4;
+ globals->find_neighbors = find_four_neighbors;
+ globals->nn = 4;
+ G_debug(1, "four pixel neighborhood");
}
else if (diagonal->answer == TRUE) {
- functions->find_pixel_neighbors = &find_eight_pixel_neighbors;
- functions->num_pn = 8;
+ globals->find_neighbors = find_eight_neighbors;
+ globals->nn = 8;
+ G_debug(1, "eight (3x3) pixel neighborhood");
}
- /* speed enhancement: Check if function pointer or IF statement is faster */
- /* default/0 for performing the scaling, but selected/1 if user has weighted values so scaling should be skipped. */
- files->weighted = weighted->answer;
+ /* default/0 for performing the scaling
+ * selected/1 if scaling should be skipped. */
+ globals->weighted = weighted->answer;
- /* check if optional seeds map was given, if not, use all pixels as starting seeds. */
- if (seeds->answer == NULL) {
- files->seeds_map = NULL;
- }
- else { /* seeds provided, check if valid map */
- files->seeds_map = seeds->answer;
- if ((files->seeds_mapset =
- G_find_raster2(files->seeds_map, "")) == NULL) {
- G_fatal_error(_("Starting seeds map not found."));
+ /* 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) {
+ G_fatal_error(_("Seeds map not found."));
}
- if (Rast_map_type(files->seeds_map, files->seeds_mapset) != CELL_TYPE) {
- G_fatal_error(_("Starting seeds map must be CELL type (integers)"));
+ if (Rast_map_type(globals->seeds, "") !=
+ CELL_TYPE) {
+ G_fatal_error(_("Seeeds map must be CELL type (integers)"));
}
}
- /* check if optional processing boundaries were provided */
- if (bounds->answer == NULL) { /*no processing constraints */
- files->bounds_map = NULL;
+ if (bounds->answer == NULL) {
+ globals->bounds_map = NULL;
}
- else { /* processing constraints given, check if valid map */
- files->bounds_map = bounds->answer;
- if ((files->bounds_mapset =
- G_find_raster2(files->bounds_map, "")) == NULL) {
+ else {
+ globals->bounds_map = bounds->answer;
+ if ((globals->bounds_mapset = G_find_raster(globals->bounds_map, "")) == NULL) {
G_fatal_error(_("Segmentation constraint/boundary map not found."));
}
- if (Rast_map_type(files->bounds_map, files->bounds_mapset) !=
+ if (Rast_map_type(globals->bounds_map, globals->bounds_mapset) !=
CELL_TYPE) {
G_fatal_error(_("Segmentation constraint map must be CELL type (integers)"));
}
}
/* other data */
- files->nrows = Rast_window_rows();
- files->ncols = Rast_window_cols();
+ globals->nrows = Rast_window_rows();
+ globals->ncols = Rast_window_cols();
- if (endt->answer != NULL && atoi(endt->answer) >= 0)
- functions->end_t = atoi(endt->answer);
+ /* debug help */
+ if (outband->answer == NULL)
+ globals->out_band = NULL;
else {
- functions->end_t = 1000;
- G_warning(_("invalid number of iterations, 1000 will be used."));
+ if (G_legal_filename(outband->answer) == TRUE)
+ globals->out_band = outband->answer;
+ else
+ G_fatal_error("Invalid output raster name for goodness of fit.");
}
- /* other parameters */
+ if (endt->answer) {
+ if (atoi(endt->answer) > 0)
+ globals->end_t = atoi(endt->answer);
+ else {
+ globals->end_t = 1000;
+ G_warning(_("Invalid number of iterations, 1000 will be used."));
+ }
+ }
+ else
+ globals->end_t = 1000;
- /* default/0 for no pathflag, but selected/1 to use Rk as next Ri if not mutually best neighbors. */
- /* functions->path = path->answer; */
- functions->path = TRUE;
- /* see notes above about pathflag. */
-
- functions->limited = limited->answer;
-
- if (outband->answer == NULL)
- files->out_band = NULL;
+ if (mem->answer && atoi(mem->answer) > 10)
+ globals->mb = atoi(mem->answer);
else {
- if (G_legal_filename(outband->answer) == TRUE)
- files->out_band = outband->answer; /* name of current means */
- else
- G_fatal_error(_("Invalid output raster name for means."));
+ globals->mb = 300;
+ G_warning(_("Invalid number of MB, 300 will be used."));
}
return TRUE;
Added: grass-addons/grass7/imagery/i.segment.xl/rclist.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/rclist.c (rev 0)
+++ grass-addons/grass7/imagery/i.segment.xl/rclist.c 2012-08-27 19:08:40 UTC (rev 52939)
@@ -0,0 +1,68 @@
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "iseg.h"
+
+void rclist_init(struct rclist *list)
+{
+ list->head = list->tail = NULL;
+
+ return;
+}
+
+void rclist_add(struct rclist *list, int row, int col)
+{
+ struct rc *new = G_malloc(sizeof(struct rc));
+
+ if (!new)
+ G_fatal_error(_("rclist out of memory"));
+
+ new->next = NULL;
+ new->row = row;
+ new->col = col;
+
+ if (list->head) {
+ list->head->next = new;
+ list->head = list->head->next;
+ }
+ else {
+ list->head = list->tail = new;
+ }
+
+ return;
+}
+
+/* return 1 if an element was dropped
+ * return 0 if list is empty
+ */
+int rclist_drop(struct rclist *list, struct rc *rc)
+{
+ if (list->tail) {
+ struct rc *next = list->tail->next;
+
+ rc->row = list->tail->row;
+ rc->col = list->tail->col;
+ G_free(list->tail);
+ list->tail = next;
+ if (!list->tail)
+ list->head = NULL;
+
+ return 1;
+ }
+ else
+ return 0;
+}
+
+void rclist_destroy(struct rclist *list)
+{
+ struct rc *next = list->tail;
+
+ while (next) {
+ next = next->next;
+ G_free(list->tail);
+ list->tail = next;
+ }
+ list->head = NULL;
+
+ return;
+}
+
Property changes on: grass-addons/grass7/imagery/i.segment.xl/rclist.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/imagery/i.segment.xl/regtree.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/regtree.c (rev 0)
+++ grass-addons/grass7/imagery/i.segment.xl/regtree.c 2012-08-27 19:08:40 UTC (rev 52939)
@@ -0,0 +1,584 @@
+/*!
+ * \file rbtree.c
+ *
+ * \brief binary search tree
+ *
+ * Generic balanced binary search tree (Red Black Tree) implementation
+ *
+ * (C) 2009 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public License
+ * (>=v2). Read the file COPYING that comes with GRASS for details.
+ *
+ * \author Original author Julienne Walker 2003, 2008
+ * GRASS implementation Markus Metz, 2009
+ */
+
+/* balanced binary search tree implementation
+ *
+ * this one is a Red Black Tree, no parent pointers, no threads
+ * The core code comes from Julienne Walker's tutorials on binary search trees
+ * original license: public domain
+ * http://eternallyconfuzzled.com/tuts/datastructures/jsw_tut_rbtree.aspx
+ * some ideas come from libavl (GPL >= 2)
+ *
+ * Red Black Trees are used to maintain a data structure with
+ * search, insertion and deletion in O(log N) time
+ */
+
+#include <assert.h>
+#include <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "regtree.h"
+
+/* internal functions */
+static struct RG_NODE *rgtree_single(struct RG_NODE *, int);
+static struct RG_NODE *rgtree_double(struct RG_NODE *, int);
+static struct reg_stats *rgtree_first(struct RG_TRAV *);
+static struct reg_stats *rgtree_next(struct RG_TRAV *);
+static struct RG_NODE *rgtree_make_node(size_t, struct reg_stats *);
+static int is_red(struct RG_NODE *);
+
+
+int compare_regstat(struct reg_stats *a, struct reg_stats *b)
+{
+ return (a->id < b->id ? -1 : (a->id > b->id));
+}
+
+
+/* create new tree and initialize
+ * returns pointer to new tree, NULL for memory allocation error
+ */
+struct RG_TREE *rgtree_create(int nbands, size_t rb_datasize)
+{
+ struct RG_TREE *tree = (struct RG_TREE *)malloc(sizeof(struct RG_TREE));
+
+ if (tree == NULL) {
+ G_warning("RB tree: Out of memory!");
+ return NULL;
+ }
+
+ tree->datasize = rb_datasize;
+ tree->cmp = compare_regstat;
+ tree->count = 0;
+ tree->nbands = nbands;
+ tree->root = NULL;
+
+ return tree;
+}
+
+/* add an item to a tree
+ * non-recursive top-down insertion
+ * the algorithm does not allow duplicates and also does not warn about a duplicate
+ * returns 1 on success, 0 on failure
+ */
+int rgtree_insert(struct RG_TREE *tree, struct reg_stats *data)
+{
+ assert(tree && data);
+
+ if (tree->root == NULL) {
+ /* create a new root node for tree */
+ tree->root = rgtree_make_node(tree->datasize, data);
+ if (tree->root == NULL)
+ return 0;
+ }
+ else {
+ struct RG_NODE head = { 0, {0, 0}, {0, 0, 0, 0} }; /* False tree root */
+ struct RG_NODE *g, *t; /* Grandparent & parent */
+ struct RG_NODE *p, *q; /* Iterator & parent */
+ int dir = 0, last = 0;
+
+ /* Set up helpers */
+ t = &head;
+ g = p = NULL;
+ q = t->link[1] = tree->root;
+
+ /* Search down the tree */
+ for (;;) {
+ if (q == NULL) {
+ /* Insert new node at the bottom */
+ p->link[dir] = q = rgtree_make_node(tree->datasize, data);
+ if (q == NULL)
+ return 0;
+ }
+ else if (is_red(q->link[0]) && is_red(q->link[1])) {
+ /* Color flip */
+ q->red = 1;
+ q->link[0]->red = 0;
+ q->link[1]->red = 0;
+ }
+
+ /* Fix red violation */
+ if (is_red(q) && is_red(p)) {
+ int dir2 = t->link[1] == g;
+
+ if (q == p->link[last])
+ t->link[dir2] = rgtree_single(g, !last);
+ else
+ t->link[dir2] = rgtree_double(g, !last);
+ }
+
+ last = dir;
+ dir = tree->cmp(&(q->data), data);
+
+ /* Stop if found. This check also disallows duplicates in the tree */
+ if (dir == 0)
+ break;
+
+ dir = dir < 0;
+
+ /* Move the helpers down */
+ if (g != NULL)
+ t = g;
+
+ g = p, p = q;
+ q = q->link[dir];
+ }
+
+ /* Update root */
+ tree->root = head.link[1];
+ }
+
+ /* Make root black */
+ tree->root->red = 0;
+
+ tree->count++;
+
+ return 1;
+}
+
+/* remove an item from a tree that matches given data
+ * non-recursive top-down removal
+ * returns 1 on successful removal
+ * returns 0 if data item was not found
+ */
+int rgtree_remove(struct RG_TREE *tree, struct reg_stats *data)
+{
+ struct RG_NODE head = { 0, {0, 0}, {0, 0, 0, 0} }; /* False tree root */
+ struct RG_NODE *q, *p, *g; /* Helpers */
+ struct RG_NODE *f = NULL; /* Found item */
+ int dir = 1, removed = 0;
+
+ assert(tree && data);
+
+ if (tree->root == NULL) {
+ return 0; /* empty tree, nothing to remove */
+ }
+
+ /* Set up helpers */
+ q = &head;
+ g = p = NULL;
+ q->link[1] = tree->root;
+
+ /* Search and push a red down */
+ while (q->link[dir] != NULL) {
+ int last = dir;
+
+ /* Update helpers */
+ g = p, p = q;
+ q = q->link[dir];
+ dir = tree->cmp(&(q->data), data);
+
+ /* Save found node */
+ if (dir == 0)
+ f = q;
+
+ dir = dir < 0;
+
+ /* Push the red node down */
+ if (!is_red(q) && !is_red(q->link[dir])) {
+ if (is_red(q->link[!dir]))
+ p = p->link[last] = rgtree_single(q, dir);
+ else if (!is_red(q->link[!dir])) {
+ struct RG_NODE *s = p->link[!last];
+
+ if (s != NULL) {
+ if (!is_red(s->link[!last]) && !is_red(s->link[last])) {
+ /* Color flip */
+ p->red = 0;
+ s->red = 1;
+ q->red = 1;
+ }
+ else {
+ int dir2 = g->link[1] == p;
+
+ if (is_red(s->link[last]))
+ g->link[dir2] = rgtree_double(p, last);
+ else if (is_red(s->link[!last]))
+ g->link[dir2] = rgtree_single(p, last);
+
+ /* Ensure correct coloring */
+ q->red = g->link[dir2]->red = 1;
+ g->link[dir2]->link[0]->red = 0;
+ g->link[dir2]->link[1]->red = 0;
+ }
+ }
+ }
+ }
+ }
+
+ /* Replace and remove if found */
+ if (f != NULL) {
+ f->data.id = q->data.id;
+ f->data.count = q->data.count;
+ memcpy(f->data.sum, q->data.sum, tree->datasize);
+ memcpy(f->data.mean, q->data.mean, tree->datasize);
+ /* unused:
+ memcpy(f->data.min, q->data.min, tree->datasize);
+ memcpy(f->data.max, q->data.max, tree->datasize);
+ */
+ p->link[p->link[1] == q] = q->link[q->link[0] == NULL];
+
+ free(q->data.sum);
+ free(q->data.mean);
+ /* unused:
+ free(q->data.min);
+ free(q->data.max);
+ */
+ free(q);
+ q = NULL;
+ tree->count--;
+ removed = 1;
+ }
+ else
+ G_debug(2, "RB tree: data not found in search tree");
+
+ /* Update root and make it black */
+ tree->root = head.link[1];
+ if (tree->root != NULL)
+ tree->root->red = 0;
+
+ return removed;
+}
+
+/* find data item in tree
+ * returns pointer to data item if found else NULL
+ */
+struct reg_stats *rgtree_find(struct RG_TREE *tree, struct reg_stats *data)
+{
+ struct RG_NODE *curr_node = tree->root;
+ int cmp;
+
+ assert(tree && data);
+
+ while (curr_node != NULL) {
+ cmp = tree->cmp(&(curr_node->data), data);
+ if (cmp == 0)
+ return &curr_node->data; /* found */
+
+ curr_node = curr_node->link[cmp < 0];
+ }
+ return NULL;
+}
+
+/* initialize tree traversal
+ * (re-)sets trav structure
+ * returns 0
+ */
+int rgtree_init_trav(struct RG_TRAV *trav, struct RG_TREE *tree)
+{
+ assert(trav && tree);
+
+ trav->tree = tree;
+ trav->curr_node = tree->root;
+ trav->first = 1;
+ trav->top = 0;
+
+ return 0;
+}
+
+/* traverse the tree in ascending order
+ * useful to get all items in the tree non-recursively
+ * struct RG_TRAV *trav needs to be initialized first
+ * returns pointer to data, NULL when finished
+ */
+struct reg_stats *rgtree_traverse(struct RG_TRAV *trav)
+{
+ assert(trav);
+
+ if (trav->curr_node == NULL) {
+ if (trav->first)
+ G_debug(1, "RB tree: empty tree");
+ else
+ G_debug(1, "RB tree: finished traversing");
+
+ return NULL;
+ }
+
+ if (!trav->first)
+ return rgtree_next(trav);
+ else {
+ trav->first = 0;
+ return rgtree_first(trav);
+ }
+}
+
+/* find start point to traverse the tree in ascending order
+ * useful to get a selection of items in the tree
+ * magnitudes faster than traversing the whole tree
+ * may return first item that's smaller or first item that's larger
+ * struct RG_TRAV *trav needs to be initialized first
+ * returns pointer to data, NULL when finished
+ */
+struct reg_stats *rgtree_traverse_start(struct RG_TRAV *trav, struct reg_stats *data)
+{
+ int dir = 0;
+
+ assert(trav && data);
+
+ if (trav->curr_node == NULL) {
+ if (trav->first)
+ G_warning("RB tree: empty tree");
+ else
+ G_warning("RB tree: finished traversing");
+
+ return NULL;
+ }
+
+ if (!trav->first)
+ return rgtree_next(trav);
+
+ /* else first time, get start node */
+
+ trav->first = 0;
+ trav->top = 0;
+
+ while (trav->curr_node != NULL) {
+ dir = trav->tree->cmp(&(trav->curr_node->data), data);
+ /* exact match, great! */
+ if (dir == 0)
+ return &(trav->curr_node->data);
+ else {
+ dir = dir < 0;
+ /* end of branch, also reached if
+ * smallest item is larger than search template or
+ * largest item is smaller than search template */
+ if (trav->curr_node->link[dir] == NULL)
+ return &(trav->curr_node->data);
+
+ trav->up[trav->top++] = trav->curr_node;
+ trav->curr_node = trav->curr_node->link[dir];
+ }
+ }
+
+ return NULL; /* should not happen */
+}
+
+/* two functions needed to fully traverse the tree: initialize and continue
+ * useful to get all items in the tree non-recursively
+ * this one here uses a stack
+ * parent pointers or threads would also be possible
+ * but these would need to be added to RG_NODE
+ * -> more memory needed for standard operations
+ */
+
+/* start traversing the tree
+ * returns pointer to smallest data item
+ */
+static struct reg_stats *rgtree_first(struct RG_TRAV *trav)
+{
+ /* get smallest item */
+ while (trav->curr_node->link[0] != NULL) {
+ trav->up[trav->top++] = trav->curr_node;
+ trav->curr_node = trav->curr_node->link[0];
+ }
+
+ return &(trav->curr_node->data); /* return smallest item */
+}
+
+/* continue traversing the tree in ascending order
+ * returns pointer to data item, NULL when finished
+ */
+static struct reg_stats *rgtree_next(struct RG_TRAV *trav)
+{
+ if (trav->curr_node->link[1] != NULL) {
+ /* something on the right side: larger item */
+ trav->up[trav->top++] = trav->curr_node;
+ trav->curr_node = trav->curr_node->link[1];
+
+ /* go down, find smallest item in this branch */
+ while (trav->curr_node->link[0] != NULL) {
+ trav->up[trav->top++] = trav->curr_node;
+ trav->curr_node = trav->curr_node->link[0];
+ }
+ }
+ else {
+ /* at smallest item in this branch, go back up */
+ struct RG_NODE *last;
+
+ do {
+ if (trav->top == 0) {
+ trav->curr_node = NULL;
+ break;
+ }
+ last = trav->curr_node;
+ trav->curr_node = trav->up[--trav->top];
+ } while (last == trav->curr_node->link[1]);
+ }
+
+ if (trav->curr_node != NULL) {
+ return &(trav->curr_node->data);
+ }
+ else
+ return NULL; /* finished traversing */
+}
+
+/* destroy the tree */
+void rgtree_destroy(struct RG_TREE *tree)
+{
+ struct RG_NODE *it;
+ struct RG_NODE *save = tree->root;
+
+ /*
+ Rotate away the left links so that
+ we can treat this like the destruction
+ of a linked list
+ */
+ while((it = save) != NULL) {
+ if (it->link[0] == NULL) {
+ /* No left links, just kill the node and move on */
+ save = it->link[1];
+ free(it->data.sum);
+ free(it->data.mean);
+ free(it);
+ it = NULL;
+ }
+ else {
+ /* Rotate away the left link and check again */
+ save = it->link[0];
+ it->link[0] = save->link[1];
+ save->link[1] = it;
+ }
+ }
+ free(tree);
+ tree = NULL;
+
+ return;
+}
+
+/* used for debugging: check for errors in tree structure */
+int rgtree_debug(struct RG_TREE *tree, struct RG_NODE *root)
+{
+ int lh, rh;
+
+ if (root == NULL)
+ return 1;
+ else {
+ struct RG_NODE *ln = root->link[0];
+ struct RG_NODE *rn = root->link[1];
+ int lcmp = 0, rcmp = 0;
+
+ /* Consecutive red links */
+ if (is_red(root)) {
+ if (is_red(ln) || is_red(rn)) {
+ G_warning("Red Black Tree debugging: Red violation");
+ return 0;
+ }
+ }
+
+ lh = rgtree_debug(tree, ln);
+ rh = rgtree_debug(tree, rn);
+
+ if (ln) {
+ lcmp = tree->cmp(&(ln->data), &(root->data));
+ }
+
+ if (rn) {
+ rcmp = tree->cmp(&(rn->data), &(root->data));
+ }
+
+ /* Invalid binary search tree:
+ * left node >= parent or right node <= parent */
+ if ((ln != NULL && lcmp > -1)
+ || (rn != NULL && rcmp < 1)) {
+ G_warning("Red Black Tree debugging: Binary tree violation");
+ return 0;
+ }
+
+ /* Black height mismatch */
+ if (lh != 0 && rh != 0 && lh != rh) {
+ G_warning("Red Black Tree debugging: Black violation");
+ return 0;
+ }
+
+ /* Only count black links */
+ if (lh != 0 && rh != 0)
+ return is_red(root) ? lh : lh + 1;
+ else
+ return 0;
+ }
+}
+
+/*******************************************************
+ * *
+ * internal functions for Red Black Tree maintenance *
+ * *
+ *******************************************************/
+
+/* add a new node to the tree */
+static struct RG_NODE *rgtree_make_node(size_t datasize, struct reg_stats *data)
+{
+ struct RG_NODE *new_node = (struct RG_NODE *)malloc(sizeof(*new_node));
+
+ if (new_node == NULL)
+ G_fatal_error("RB Search Tree: Out of memory!");
+
+ if ((new_node->data.sum = malloc(datasize)) == NULL)
+ G_fatal_error("RB Search Tree: Out of memory!");
+ if ((new_node->data.mean = malloc(datasize)) == NULL)
+ G_fatal_error("RB Search Tree: Out of memory!");
+ /* unused:
+ if ((new_node->data.min = malloc(datasize)) == NULL)
+ G_fatal_error("RB Search Tree: Out of memory!");
+ if ((new_node->data.max = malloc(datasize)) == NULL)
+ G_fatal_error("RB Search Tree: Out of memory!");
+ */
+
+
+ new_node->data.id = data->id;
+ new_node->data.count = data->count;
+ memcpy(new_node->data.sum, data->sum, datasize);
+ memcpy(new_node->data.mean, data->mean, datasize);
+ /* unused
+ memcpy(new_node->data.min, data->min, datasize);
+ memcpy(new_node->data.max, data->max, datasize);
+ */
+
+ new_node->red = 1; /* 1 is red, 0 is black */
+ new_node->link[0] = NULL;
+ new_node->link[1] = NULL;
+
+ return new_node;
+}
+
+/* check for red violation */
+static int is_red(struct RG_NODE *root)
+{
+ if (root)
+ return root->red == 1;
+
+ return 0;
+}
+
+/* single rotation */
+static struct RG_NODE *rgtree_single(struct RG_NODE *root, int dir)
+{
+ struct RG_NODE *newroot = root->link[!dir];
+
+ root->link[!dir] = newroot->link[dir];
+ newroot->link[dir] = root;
+
+ root->red = 1;
+ newroot->red = 0;
+
+ return newroot;
+}
+
+/* double rotation */
+static struct RG_NODE *rgtree_double(struct RG_NODE *root, int dir)
+{
+ root->link[!dir] = rgtree_single(root->link[!dir], !dir);
+ return rgtree_single(root, dir);
+}
Property changes on: grass-addons/grass7/imagery/i.segment.xl/regtree.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/imagery/i.segment.xl/regtree.h
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/regtree.h (rev 0)
+++ grass-addons/grass7/imagery/i.segment.xl/regtree.h 2012-08-27 19:08:40 UTC (rev 52939)
@@ -0,0 +1,136 @@
+/*************************************************************
+ * USAGE *
+ *************************************************************
+ *
+ * NOTE: duplicates are not supported
+ *
+ * custom compare function
+ * extern int my_compare_fn(const void *, const void *);
+ * int my_compare_fn(const void *a, const void *b) {
+ * if ((mydatastruct *) a < (mydatastruct *) b)
+ * return -1;
+ * else if ((mydatastruct *) a > (mydatastruct *) b)
+ * return 1;
+ * else if ((mydatastruct *) a == (mydatastruct *) b)
+ * return 0;
+ * }
+ *
+ * create and initialize tree:
+ * struct RG_TREE *mytree = rgtree_create(my_compare_fn, item_size);
+ *
+ * insert items to tree:
+ * struct mydatastruct data = <some data>;
+ * if (rgtree_insert(mytree, &data) == 0)
+ * G_warning("could not insert data");
+ *
+ * find item in tree:
+ * struct mydatastruct data = <some data>;
+ * if (rgtree_find(mytree, &data) == 0)
+ * G_message("data not found");
+ *
+ * delete item from tree:
+ * struct mydatastruct data = <some data>;
+ * if (rgtree_remove(mytree, &data) == 0)
+ * G_warning("could not find data in tree");
+ *
+ * traverse tree (get all items in tree in ascending order):
+ * struct RG_TRAV trav;
+ * rgtree_init_trav(&trav, tree);
+ * while ((data = rgtree_traverse(&trav)) != NULL) {
+ * if (my_compare_fn(data, threshold_data) == 0) break;
+ * <do something with data>;
+ * }
+ *
+ * get a selection of items: all data > data1 and < data2
+ * start in tree where data is last smaller or first larger compared to data1
+ * struct RG_TRAV trav;
+ * rgtree_init_trav(&trav, tree);
+ * data = rgtree_traverse_start(&trav, &data1);
+ * <do something with data>;
+ * while ((data = rgtree_traverse(&trav)) != NULL) {
+ * if (data > data2) break;
+ * <do something with data>;
+ * }
+ *
+ * destroy tree:
+ * rgtree_destroy(mytree);
+ *
+ * debug the whole tree with
+ * rgtree_debug(mytree, mytree->root);
+ *
+ *************************************************************/
+
+#ifndef REGTREE_H
+#define REGTREE_H
+
+#include <stddef.h>
+
+/* maximum RB Tree height */
+#define REGTREE_MAX_HEIGHT 64 /* should be more than enough */
+
+/* routine to compare data items
+ * return -1 if rb_a < rb_b
+ * return 0 if rb_a == rb_b
+ * return 1 if rb_a > rb_b
+ */
+struct reg_stats
+{
+ int id; /* region ID */
+ int count; /* number of cells of this region */
+ double *sum, /* sum for each band */
+ *mean; /* mean for each band = sum[b] / count */
+
+ /* unused; stddev and thus sumsq may be needed for
+ * eCognition-like multi-scale segmentation */
+ /*
+ *min,
+ *max,
+ *sumsq,
+ *stddev;
+ */
+};
+
+typedef int rg_compare_fn(struct reg_stats *rb_a, struct reg_stats *rb_b);
+
+struct RG_NODE
+{
+ unsigned char red; /* 0 = black, 1 = red */
+ struct RG_NODE *link[2]; /* link to children: link[0] for smaller, link[1] for larger */
+ struct reg_stats data; /* any kind of data */
+};
+
+struct RG_TREE
+{
+ struct RG_NODE *root; /* root node */
+ size_t datasize; /* item size */
+ size_t count; /* number of items in tree. */
+ int nbands; /* number of bands */
+ rg_compare_fn *cmp; /* function to compare data */
+};
+
+struct RG_TRAV
+{
+ struct RG_TREE *tree; /* tree being traversed */
+ struct RG_NODE *curr_node; /* current node */
+ struct RG_NODE *up[REGTREE_MAX_HEIGHT]; /* stack of parent nodes */
+ int top; /* index for stack */
+ int first; /* little helper flag */
+};
+
+
+/* tree functions */
+struct RG_TREE *rgtree_create(int, size_t);
+void rgtree_destroy(struct RG_TREE *);
+int rgtree_insert(struct RG_TREE *, struct reg_stats *);
+int rgtree_remove(struct RG_TREE *, struct reg_stats *);
+struct reg_stats *rgtree_find(struct RG_TREE *, struct reg_stats *);
+
+/* tree traversal functions */
+int rgtree_init_trav(struct RG_TRAV *, struct RG_TREE *);
+struct reg_stats *rgtree_traverse(struct RG_TRAV *);
+struct reg_stats *rgtree_traverse_start(struct RG_TRAV *, struct reg_stats *);
+
+/* debug tree from given node downwards */
+int rgtree_debug(struct RG_TREE *, struct RG_NODE *);
+
+#endif
Property changes on: grass-addons/grass7/imagery/i.segment.xl/regtree.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Modified: grass-addons/grass7/imagery/i.segment.xl/write_output.c
===================================================================
--- grass-addons/grass7/imagery/i.segment.xl/write_output.c 2012-08-27 18:43:23 UTC (rev 52938)
+++ grass-addons/grass7/imagery/i.segment.xl/write_output.c 2012-08-27 19:08:40 UTC (rev 52939)
@@ -1,68 +1,135 @@
-/* write_output(): transfer the segmented regions from the segmented data file to a raster file */
-/* close_files(): close SEG files and free memory */
+/* transfer the segmented regions from the segmented data file to a raster file */
+/* close_files() function is at bottom */
#include <stdlib.h>
#include <grass/gis.h>
#include <grass/raster.h>
#include <grass/segment.h> /* segmentation library */
+#include <grass/glocale.h>
#include "iseg.h"
-int write_output(struct files *files)
+/* TODO some time delay here with meanbuf, etc being processed. I only put if statements on the actual files
+ * to try and keep the code more clear. Need to see if this raster makes stats processing easier? Or IFDEF it out?
+ */
+
+int write_output(struct globals *globals)
{
- int out_fd, mean_fd, row, col; /* mean_fd for validiating/debug of means */
- CELL *outbuf;
- DCELL *meanbuf;
+ int out_fd, mean_fd, row, col;
+ CELL *outbuf, rid;
+ FCELL *meanbuf;
struct Colors colors;
- struct History history;
+ double thresh, maxdev, sim, mingood;
+ struct reg_stats *rs_found;
+ struct ngbr_stats Ri, Rk;
+ outbuf = Rast_allocate_c_buf();
- outbuf = Rast_allocate_c_buf(); /* hold one row of data to put into raster */
- meanbuf = Rast_allocate_d_buf();
+ G_debug(1, "preparing output raster");
+ /* open output raster map */
+ out_fd = Rast_open_new(globals->out_name, CELL_TYPE);
+ if (globals->out_band) {
+ mean_fd = Rast_open_new(globals->out_band, FCELL_TYPE);
+ meanbuf = Rast_allocate_f_buf();
+ }
+ else {
+ mean_fd = -1;
+ meanbuf = NULL;
+ }
- /* force all data to disk */
- segment_flush(&files->bands_seg);
- segment_flush(&files->iseg_seg);
+ /* goodness of fit for each cell: 1 = good fit, 0 = bad fit */
+ /* similarity of each cell to region mean
+ * max possible difference: globals->threshold
+ * if similarity < globals->alpha * globals->alpha * globals->threshold
+ * 1
+ * else
+ * (similarity - globals->alpha * globals->alpha * globals->threshold) /
+ * (globals->threshold * (1 - globals->alpha * globals->alpha) */
- /* open output raster map */
- out_fd = Rast_open_new(files->out_name, CELL_TYPE);
- if (files->out_band != NULL)
- mean_fd = Rast_open_new(files->out_band, DCELL_TYPE);
+ thresh = globals->alpha * globals->alpha * globals->threshold;
+ maxdev = globals->threshold * (1 - globals->alpha * globals->alpha);
+ mingood = 1;
- /* transfer data from segmentation file to raster */
- for (row = 0; row < files->nrows; row++) {
- Rast_set_c_null_value(outbuf, files->ncols); /*set buffer to NULLs, only write those that weren't originally masked */
- Rast_set_d_null_value(meanbuf, files->ncols);
- for (col = 0; col < files->ncols; col++) {
- segment_get(&files->bands_seg, (void *)files->bands_val, row,
- col);
- if (!(FLAG_GET(files->null_flag, row, col))) {
- segment_get(&files->iseg_seg, &(outbuf[col]), row, col);
- meanbuf[col] = files->bands_val[0];
+ G_debug(1, "start data transfer from segmentation file to raster");
+
+ G_message(_("Writing output"));
+ for (row = 0; row < globals->nrows; row++) {
+
+ G_percent(row, globals->nrows, 9);
+
+ Rast_set_c_null_value(outbuf, globals->ncols);
+ if (globals->out_band)
+ Rast_set_f_null_value(meanbuf, globals->ncols);
+ for (col = 0; col < globals->ncols; col++) {
+
+ if (!(FLAG_GET(globals->null_flag, row, col))) {
+ segment_get(&globals->rid_seg, (void *) &rid, row, col);
+
+ if (rid > 0) {
+ outbuf[col] = rid;
+
+ if (globals->out_band) {
+
+ /* get values for Ri = larger region */
+ globals->rs.id = rid;
+ rs_found = rgtree_find(globals->reg_tree, &(globals->rs));
+
+ if (rs_found != NULL) {
+ Ri.mean = rs_found->mean;
+
+ /* get values for Rk = this cell */
+ segment_get(&globals->bands_seg,
+ (void *)globals->second_val, row, col);
+
+ Rk.mean = globals->second_val;
+
+ /* calculate similarity */
+ sim = (*globals->calculate_similarity) (&Ri, &Rk, globals);
+ }
+ else
+ /* region consists of only one cell */
+ sim = 0.;
+
+ if (0) {
+ if (sim < thresh)
+ meanbuf[col] = 1;
+ else {
+ sim = 1. - (sim - thresh) / maxdev;
+ meanbuf[col] = sim;
+ if (mingood > sim)
+ mingood = sim;
+ }
+ }
+ else {
+ sim = 1 - sim / globals->threshold;
+ meanbuf[col] = sim;
+ if (mingood > sim)
+ mingood = sim;
+ }
+ }
+ }
}
}
Rast_put_row(out_fd, outbuf, CELL_TYPE);
- if (files->out_band != NULL)
- Rast_put_row(mean_fd, meanbuf, DCELL_TYPE);
-
- G_percent(row, files->nrows, 1);
+ if (globals->out_band)
+ Rast_put_row(mean_fd, meanbuf, FCELL_TYPE);
}
/* close and save file */
Rast_close(out_fd);
- if (files->out_band != NULL)
+ if (globals->out_band)
Rast_close(mean_fd);
/* set colors */
Rast_init_colors(&colors);
- Rast_make_random_colors(&colors, 1, files->nrows * files->ncols);
- Rast_write_colors(files->out_name, G_mapset(), &colors);
+ Rast_make_random_colors(&colors, 1, globals->nrows * globals->ncols);
+ Rast_write_colors(globals->out_name, G_mapset(), &colors);
+
+ if (globals->out_band) {
+ Rast_init_colors(&colors);
+ Rast_make_grey_scale_fp_colors(&colors, mingood, 1);
+ Rast_write_colors(globals->out_band, G_mapset(), &colors);
+ }
- /* add command line to history */
- /* see i.pca as an example of setting custom info */
- Rast_short_history(files->out_name, "raster", &history);
- Rast_command_history(&history);
- Rast_write_history(files->out_name, &history);
-
/* free memory */
G_free(outbuf);
G_free(meanbuf);
@@ -71,24 +138,29 @@
return TRUE;
}
-int close_files(struct files *files)
+int close_files(struct globals *globals)
{
-
/* close segmentation files and output raster */
- segment_close(&files->bands_seg);
- if (files->bounds_map != NULL)
- segment_close(&files->bounds_seg);
+ G_debug(1, "closing files");
+ segment_close(&globals->bands_seg);
+ if (globals->bounds_map)
+ segment_close(&globals->bounds_seg);
- G_free(files->bands_val);
- G_free(files->second_val);
+ G_free(globals->bands_val);
+ G_free(globals->second_val);
- segment_close(&files->iseg_seg);
+ segment_close(&globals->rid_seg);
- flag_destroy(files->null_flag);
- flag_destroy(files->candidate_flag);
- flag_destroy(files->seeds_flag);
+ /*
+ for (i = 0; i < globals->nrows; i++)
+ G_free(globals->rid[i]);
+ G_free(globals->rid);
+ */
- link_cleanup(files->token);
+ flag_destroy(globals->null_flag);
+ flag_destroy(globals->candidate_flag);
+ /* anything else left to clean up? */
+
return TRUE;
}
More information about the grass-commit
mailing list