[GRASS-SVN] r68944 - grass/trunk/imagery/i.segment
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jul 11 11:52:07 PDT 2016
Author: mmetz
Date: 2016-07-11 11:52:07 -0700 (Mon, 11 Jul 2016)
New Revision: 68944
Modified:
grass/trunk/imagery/i.segment/create_isegs.c
grass/trunk/imagery/i.segment/iseg.h
grass/trunk/imagery/i.segment/main.c
grass/trunk/imagery/i.segment/ngbrtree.c
grass/trunk/imagery/i.segment/ngbrtree.h
grass/trunk/imagery/i.segment/open_files.c
grass/trunk/imagery/i.segment/parse_args.c
grass/trunk/imagery/i.segment/region_growing.c
grass/trunk/imagery/i.segment/regtree.c
grass/trunk/imagery/i.segment/write_output.c
Log:
i.segment: improvement for #3084
Modified: grass/trunk/imagery/i.segment/create_isegs.c
===================================================================
--- grass/trunk/imagery/i.segment/create_isegs.c 2016-07-10 16:39:23 UTC (rev 68943)
+++ grass/trunk/imagery/i.segment/create_isegs.c 2016-07-11 18:52:07 UTC (rev 68944)
@@ -84,6 +84,48 @@
}
}
+ if (globals->method == ORM_RG) {
+ /* renumber */
+ int i, *new_id, max_id;
+
+ G_debug(1, "Largest assigned ID: %d", globals->max_rid);
+
+ new_id = G_malloc((globals->max_rid + 1) * sizeof(int));
+
+ for (i = 0; i <= globals->max_rid; i++)
+ new_id[i] = 0;
+
+ for (row = 0; row < globals->nrows; row++) {
+ for (col = 0; col < globals->ncols; col++) {
+ Segment_get(&globals->rid_seg, &rid, row, col);
+ if (!Rast_is_c_null_value(&rid))
+ new_id[rid]++;
+ }
+ }
+
+ max_id = 0;
+ for (i = 0; i <= globals->max_rid; i++) {
+ if (new_id[i] > 0) {
+ max_id++;
+ new_id[i] = max_id;
+ }
+ }
+ globals->max_rid = max_id;
+ G_debug(1, "Largest renumbered ID: %d", globals->max_rid);
+
+ for (row = 0; row < globals->nrows; row++) {
+ for (col = 0; col < globals->ncols; col++) {
+ Segment_get(&globals->rid_seg, &rid, row, col);
+ if (!Rast_is_c_null_value(&rid)) {
+ rid = new_id[rid];
+ Segment_put(&globals->rid_seg, &rid, row, col);
+ }
+ }
+ }
+
+ G_free(new_id);
+ }
+
return successflag;
}
Modified: grass/trunk/imagery/i.segment/iseg.h
===================================================================
--- grass/trunk/imagery/i.segment/iseg.h 2016-07-10 16:39:23 UTC (rev 68943)
+++ grass/trunk/imagery/i.segment/iseg.h 2016-07-11 18:52:07 UTC (rev 68944)
@@ -13,12 +13,21 @@
*****************************************************************************/
#include <grass/segment.h>
+#include <grass/imagery.h>
#include "flag.h"
#include "regtree.h"
#include "ngbrtree.h"
/* #def _OR_SHAPE_ */
+#ifdef HAVE_LONG_LONG_INT
+#define LARGEINT long long
+#elif defined HAVE_LARGEFILES
+#define LARGEINT off_t
+#else
+#define LARGEINT long
+#endif
+
/* methods */
#define ORM_RG 1 /* region growing */
#define ORM_MS 2 /* mean shift */
@@ -41,59 +50,76 @@
/* input and output files, as well as some other processing info */
struct globals
{
- /* user parameters */
+ /* input */
char *image_group;
+ struct Ref Ref; /* group reference list */
+ DCELL *min, *max;
int weighted; /* 0 if false/not selected, so we should scale input.
* 1 if the scaling should be skipped */
+ int nbands; /* number of rasters in the image group */
+ size_t datasize; /* nbands * sizeof(double) */
+ int mb; /* amount of memory to use in MB */
+
+ char *seeds, *bounds_map; /* optional segment seeds and polygon constraints/boundaries */
+ CELL lower_bound, upper_bound;
+ const char *bounds_mapset;
+
+ /* output */
+ /* region growing */
+ char *out_name; /* name of output raster map with regions */
+ char *out_band; /* indicator for segment heterogeneity / goodness of fit */
+ /* mean_shift */
+ char *ms_suf; /* suffix to be appended to input bands */
+
+ /* general segmentation */
int method; /* Segmentation method code */
int (*method_fn)(); /* Segmentation method function */
int nn; /* number of neighbors, 4 or 8 */
double max_diff; /* max possible difference */
double alpha; /* similarity threshold */
+ int end_t; /* maximum number of iterations */
+
+ /* region growing */
int min_segment_size; /* smallest number of pixels/cells allowed in a final segment */
+ /* inactive options for region growing */
double radio_weight; /* weighing factor radiometric - shape */
double smooth_weight; /* weighing factor smoothness - compactness */
- int end_t; /* maximum number of iterations */
- int mb;
+ /* mean shift */
+ double hs, hr; /* spectral and range bandwidth */
/* region info */
int nrows, ncols;
int row_min, row_max, col_min, col_max; /* region constraints */
- long ncells, notnullcells;
+ LARGEINT ncells, notnullcells;
- char *out_name; /* name of output raster map */
- 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 */
- int nbands; /* number of rasters in the image group */
SEGMENT bands_seg, /* input group with one or more bands */
+ bands_seg2, /* copy of bands_seg for mean shift */
bounds_seg,
rid_seg;
+ SEGMENT *bands_in, *bands_out; /* pointers to input/output bands_seg, for mean shift */
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 */
+ /* maximum used region ID */
+ CELL max_rid;
+
+ /* region growing internal structure */
struct RG_TREE *reg_tree; /* search tree with region stats */
int min_reg_size; /* minimum region size */
struct reg_stats rs, rs_i, rs_k;
struct ngbr_stats ns;
- size_t datasize; /* nbands * sizeof(double) */
- int n_regions;
/* processing flags */
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? */
/* number of remaining cells to check */
- long candidate_count;
+ LARGEINT candidate_count;
/* functions */
-
void (*find_neighbors) (int, int, int[8][2]); /*parameters: row, col, neighbors */
double (*calculate_similarity) (struct ngbr_stats *,
struct ngbr_stats *,
@@ -142,7 +168,8 @@
int rclist_drop(struct rclist *, struct rc *);
void rclist_destroy(struct rclist *);
-
/* write_output.c */
-int write_output(struct globals *);
+int write_ids(struct globals *);
+int write_gof_rg(struct globals *);
+int write_bands_ms(struct globals *);
int close_files(struct globals *);
Modified: grass/trunk/imagery/i.segment/main.c
===================================================================
--- grass/trunk/imagery/i.segment/main.c 2016-07-10 16:39:23 UTC (rev 68943)
+++ grass/trunk/imagery/i.segment/main.c 2016-07-11 18:52:07 UTC (rev 68944)
@@ -20,6 +20,7 @@
*****************************************************************************/
#include <stdlib.h>
+#include <stdio.h>
#include <grass/gis.h>
#include <grass/glocale.h>
#include "iseg.h"
@@ -40,7 +41,7 @@
_("Identifies segments (objects) from imagery data.");
parse_args(argc, argv, &globals);
-
+
G_debug(1, "Main: starting open_files()");
if (open_files(&globals) != TRUE)
G_fatal_error(_("Error in reading data"));
@@ -50,13 +51,23 @@
G_fatal_error(_("Error in creating segments"));
G_debug(1, "Main: starting write_output()");
- if (write_output(&globals) != TRUE)
- G_fatal_error(_("Error in writing data"));
+ if (write_ids(&globals) != TRUE)
+ G_fatal_error(_("Error in writing IDs"));
+ if (globals.method == ORM_RG && globals.out_band) {
+ if (write_gof_rg(&globals) != TRUE)
+ G_fatal_error(_("Error in writing goodness of fit"));
+ }
+
+ if (globals.method == ORM_MS) {
+ if (write_bands_ms(&globals) != TRUE)
+ G_fatal_error(_("Error in writing new band values"));
+ }
+
G_debug(1, "Main: starting close_files()");
close_files(&globals);
- G_done_msg(_("Number of segments created: %d"), globals.n_regions);
+ G_done_msg(_("Number of segments created: %d"), globals.max_rid);
exit(EXIT_SUCCESS);
}
Modified: grass/trunk/imagery/i.segment/ngbrtree.c
===================================================================
--- grass/trunk/imagery/i.segment/ngbrtree.c 2016-07-10 16:39:23 UTC (rev 68943)
+++ grass/trunk/imagery/i.segment/ngbrtree.c 2016-07-11 18:52:07 UTC (rev 68944)
@@ -44,7 +44,13 @@
int cmp_ngbr(struct ngbr_stats *a, struct ngbr_stats *b)
{
- return (a->id - b->id);
+ if (a->id > 0 || b->id > 0)
+ return (a->id - b->id);
+
+ /* a->id == b->id == 0 */
+ if (a->row == b->row)
+ return (a->col - b->col);
+ return (a->row - b->row);
}
@@ -417,7 +423,7 @@
return NULL; /* finished traversing */
}
-/* destroy the tree */
+/* clear the tree: remove all nodes */
void nbtree_clear(struct NB_TREE *tree)
{
struct NB_NODE *it;
@@ -451,6 +457,15 @@
return;
}
+/* destroy the tree */
+void nbtree_destroy(struct NB_TREE *tree)
+{
+ nbtree_clear(tree);
+ free(tree);
+
+ return;
+}
+
/* used for debugging: check for errors in tree structure */
int nbtree_debug(struct NB_TREE *tree, struct NB_NODE *root)
{
Modified: grass/trunk/imagery/i.segment/ngbrtree.h
===================================================================
--- grass/trunk/imagery/i.segment/ngbrtree.h 2016-07-10 16:39:23 UTC (rev 68943)
+++ grass/trunk/imagery/i.segment/ngbrtree.h 2016-07-11 18:52:07 UTC (rev 68944)
@@ -110,9 +110,11 @@
/* tree functions */
struct NB_TREE *nbtree_create(int, size_t);
void nbtree_clear(struct NB_TREE *);
+void nbtree_destroy(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 *);
+int cmp_ngbr(struct ngbr_stats *, struct ngbr_stats *);
/* tree traversal functions */
int nbtree_init_trav(struct NB_TRAV *, struct NB_TREE *);
Modified: grass/trunk/imagery/i.segment/open_files.c
===================================================================
--- grass/trunk/imagery/i.segment/open_files.c 2016-07-10 16:39:23 UTC (rev 68943)
+++ grass/trunk/imagery/i.segment/open_files.c 2016-07-11 18:52:07 UTC (rev 68944)
@@ -13,13 +13,12 @@
int open_files(struct globals *globals)
{
- struct Ref Ref; /* group reference list */
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_val;
int have_bounds = 0;
- CELL s, id;
+ CELL id;
struct Range range; /* min/max values of bounds map */
struct FPRange *fp_range; /* min/max values of each input raster */
DCELL *min, *max;
@@ -36,38 +35,41 @@
/* ****** open the input rasters ******* */
- if (!I_get_group_ref(globals->image_group, &Ref))
+ if (!I_get_group_ref(globals->image_group, &globals->Ref))
G_fatal_error(_("Group <%s> not found in the current mapset"),
globals->image_group);
- if (Ref.nfiles <= 0)
+ if (globals->Ref.nfiles <= 0)
G_fatal_error(_("Group <%s> contains no raster maps"),
globals->image_group);
/* Read Imagery Group */
- in_fd = G_malloc(Ref.nfiles * sizeof(int));
- inbuf = (DCELL **) G_malloc(Ref.nfiles * sizeof(DCELL *));
- fp_range = G_malloc(Ref.nfiles * sizeof(struct FPRange));
- min = G_malloc(Ref.nfiles * sizeof(DCELL));
- max = G_malloc(Ref.nfiles * sizeof(DCELL));
+ in_fd = G_malloc(globals->Ref.nfiles * sizeof(int));
+ inbuf = (DCELL **) G_malloc(globals->Ref.nfiles * sizeof(DCELL *));
+ fp_range = G_malloc(globals->Ref.nfiles * sizeof(struct FPRange));
+ min = G_malloc(globals->Ref.nfiles * sizeof(DCELL));
+ max = G_malloc(globals->Ref.nfiles * sizeof(DCELL));
+
+ globals->min = min;
+ globals->max = max;
G_debug(1, "Opening input rasters...");
- for (n = 0; n < Ref.nfiles; n++) {
+ for (n = 0; n < globals->Ref.nfiles; n++) {
inbuf[n] = Rast_allocate_d_buf();
- in_fd[n] = Rast_open_old(Ref.file[n].name, Ref.file[n].mapset);
+ in_fd[n] = Rast_open_old(globals->Ref.file[n].name, globals->Ref.file[n].mapset);
}
/* Get min/max values of each input raster for scaling */
globals->max_diff = 0.;
- globals->nbands = Ref.nfiles;
+ globals->nbands = globals->Ref.nfiles;
- for (n = 0; n < Ref.nfiles; n++) {
+ for (n = 0; n < globals->Ref.nfiles; n++) {
/* returns -1 on error, 2 on empty range, quitting either way. */
- if (Rast_read_fp_range(Ref.file[n].name, Ref.file[n].mapset, &fp_range[n]) != 1)
+ if (Rast_read_fp_range(globals->Ref.file[n].name, globals->Ref.file[n].mapset, &fp_range[n]) != 1)
G_fatal_error(_("No min/max found in raster map <%s>"),
- Ref.file[n].name);
+ globals->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",
@@ -75,7 +77,7 @@
}
if (globals->weighted == FALSE)
- globals->max_diff = Ref.nfiles;
+ globals->max_diff = globals->Ref.nfiles;
else {
/* max difference with selected similarity method */
Ri.mean = max;
@@ -89,21 +91,21 @@
/* size of each element to be stored */
- inlen = sizeof(DCELL) * Ref.nfiles;
+ inlen = sizeof(DCELL) * globals->Ref.nfiles;
outlen = sizeof(CELL);
G_debug(1, "data element size, in: %d , out: %d ", inlen, outlen);
globals->datasize = sizeof(double) * globals->nbands;
/* count non-null cells */
- globals->notnullcells = (long)globals->nrows * globals->ncols;
+ globals->notnullcells = (LARGEINT)globals->nrows * globals->ncols;
for (row = 0; row < globals->nrows; row++) {
- for (n = 0; n < Ref.nfiles; n++) {
+ for (n = 0; n < globals->Ref.nfiles; n++) {
Rast_get_d_row(in_fd[n], inbuf[n], row);
}
for (col = 0; col < globals->ncols; col++) {
is_null = 0; /*Assume there is data */
- for (n = 0; n < Ref.nfiles; n++) {
+ for (n = 0; n < globals->Ref.nfiles; n++) {
if (Rast_is_d_null_value(&inbuf[n][col])) {
is_null = 1;
}
@@ -114,7 +116,6 @@
}
}
}
- G_verbose_message(_("Non-NULL cells: %ld"), globals->notnullcells);
if (globals->notnullcells < 2)
G_fatal_error(_("Insufficient number of non-NULL cells in current region"));
@@ -130,35 +131,45 @@
scols, inlen, nseg) != 1)
G_fatal_error("Unable to create input temporary files");
+ if (globals->method == ORM_MS) {
+ if (Segment_open
+ (&globals->bands_seg2, G_tempfile(), globals->nrows, globals->ncols, srows,
+ scols, inlen, nseg) != 1)
+ G_fatal_error("Unable to create input temporary files");
+
+ globals->bands_in = &globals->bands_seg;
+ globals->bands_out = &globals->bands_seg2;
+ }
+
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");
/* load input bands to segment structure */
- if (Ref.nfiles > 1)
+ if (globals->Ref.nfiles > 1)
G_message(_("Loading input bands..."));
else
G_message(_("Loading input band..."));
globals->bands_val = (double *)G_malloc(inlen);
globals->second_val = (double *)G_malloc(inlen);
- /* initial segment ID */
- s = 1;
+ globals->max_rid = 0;
+
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++) {
G_percent(row, globals->nrows, 4);
- for (n = 0; n < Ref.nfiles; n++) {
+ for (n = 0; n < globals->Ref.nfiles; n++) {
Rast_get_d_row(in_fd[n], inbuf[n], row);
}
for (col = 0; col < globals->ncols; col++) {
is_null = 0; /*Assume there is data */
- for (n = 0; n < Ref.nfiles; n++) {
+ for (n = 0; n < globals->Ref.nfiles; n++) {
globals->bands_val[n] = inbuf[n][col];
if (Rast_is_d_null_value(&inbuf[n][col])) {
is_null = 1;
@@ -173,13 +184,14 @@
(void *)globals->bands_val, row, col) != 1)
G_fatal_error(_("Unable to write to temporary file"));
+ if (globals->method == ORM_MS) {
+ if (Segment_put(&globals->bands_seg2,
+ (void *)globals->bands_val, row, col) != 1)
+ G_fatal_error(_("Unable to write to temporary file"));
+ }
+
+ id = 0;
if (!is_null) {
- 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;
@@ -195,11 +207,9 @@
Rast_set_c_null_value(&id, 1);
FLAG_SET(globals->null_flag, row, col);
}
- if (!globals->seeds || is_null) {
- if (Segment_put(&globals->rid_seg,
- (void *)&id, row, col) != 1)
- G_fatal_error(_("Unable to write to temporary file"));
- }
+ if (Segment_put(&globals->rid_seg,
+ (void *)&id, row, col) != 1)
+ G_fatal_error(_("Unable to write to temporary file"));
}
}
G_percent(1, 1, 1);
@@ -210,7 +220,7 @@
globals->row_max++;
globals->col_max++;
- globals->ncells = (long)(globals->row_max - globals->row_min) *
+ globals->ncells = (LARGEINT)(globals->row_max - globals->row_min) *
(globals->col_max - globals->col_min);
/* bounds/constraints */
@@ -281,7 +291,7 @@
/* Free memory */
- for (n = 0; n < Ref.nfiles; n++) {
+ for (n = 0; n < globals->Ref.nfiles; n++) {
G_free(inbuf[n]);
Rast_close(in_fd[n]);
}
@@ -290,19 +300,15 @@
globals->rs.mean = G_malloc(globals->datasize);
globals->reg_tree = rgtree_create(globals->nbands, globals->datasize);
- globals->n_regions = s - 1;
- if (globals->seeds) {
+ if (globals->method == ORM_RG && globals->seeds) {
load_seeds(globals, srows, scols, nseg);
+ G_debug(1, "Number of initial regions: %d", globals->max_rid);
}
- G_debug(1, "Number of initial regions: %d", globals->n_regions);
-
G_free(inbuf);
G_free(in_fd);
G_free(fp_range);
- G_free(min);
- G_free(max);
return TRUE;
}
@@ -313,12 +319,17 @@
int row, col;
SEGMENT seeds_seg;
CELL *seeds_buf, seeds_val;
- int seeds_fd;
- int spos, sneg, have_seeds;
+ int seeds_fd, have_seeds;
+ CELL new_id, cellmax, noid;
struct rc Ri;
G_debug(1, "load_seeds()");
-
+
+ cellmax = (1 << (sizeof(CELL) * 8 - 2)) - 1;
+ cellmax += (1 << (sizeof(CELL) * 8 - 2));
+
+ noid = 0;
+
G_message(_("Loading seeds from raster map <%s>..."), globals->seeds);
if (Segment_open
@@ -356,8 +367,7 @@
return 0;
}
- spos = 1;
- sneg = -1;
+ new_id = 0;
/* convert seeds to regions */
G_debug(1, "convert seeds to regions");
@@ -368,17 +378,18 @@
if (!(FLAG_GET(globals->null_flag, row, col)) &&
!(FLAG_GET(globals->candidate_flag, row, col))) {
- if (Rast_is_c_null_value(&(seeds_buf[col]))) {
- if (Segment_put(&globals->rid_seg, &sneg, row, col) != 1)
- G_fatal_error(_("Unable to write to temporary file"));
- sneg--;
- globals->n_regions--;
- }
- else {
+ if (!Rast_is_c_null_value(&(seeds_buf[col]))) {
+ if (new_id == cellmax)
+ G_fatal_error(_("Too many seeds: integer overflow"));
+
+ new_id++;
+
Ri.row = row;
Ri.col = col;
- read_seed(globals, &seeds_seg, &Ri, spos);
- spos++;
+ if (!read_seed(globals, &seeds_seg, &Ri, new_id)) {
+ new_id--;
+ Segment_put(&globals->rid_seg, (void *)&noid, Ri.row, Ri.col);
+ }
}
}
}
@@ -388,7 +399,7 @@
Rast_close(seeds_fd);
Segment_close(&seeds_seg);
- globals->n_regions = spos - 1;
+ globals->max_rid = new_id;
flag_clear_all(globals->candidate_flag);
@@ -499,9 +510,10 @@
else {
if (globals->rs.count > 1)
update_band_vals(Ri->row, Ri->col, &(globals->rs), globals);
+ else if (globals->rs.count == 1) {
+ return 0;
+ }
}
- if (globals->rs.count > 1)
- globals->n_regions -= (globals->rs.count - 1);
return 1;
}
@@ -511,46 +523,64 @@
{
double reg_size_mb, segs_mb;
int reg_size_count, nseg, nseg_total;
+
+ segs_mb = globals->mb;
+ if (globals->method == ORM_RG) {
- /* minimum region size to store in search tree */
- reg_size_mb = 2 * globals->datasize + /* mean, sum */
- 2 * sizeof(int) + /* id, count */
- sizeof(unsigned char) +
- 2 * sizeof(struct REG_NODE *);
- reg_size_mb /= (1024. * 1024.);
+ /* minimum region size to store in search tree */
+ reg_size_mb = 2 * globals->datasize + /* mean, sum */
+ 2 * sizeof(int) + /* id, count */
+ sizeof(unsigned char) +
+ 2 * sizeof(struct REG_NODE *);
+ reg_size_mb /= (1024. * 1024.);
- /* put aside some memory for segment structures */
- segs_mb = globals->mb * 0.1;
- if (segs_mb > 10)
- segs_mb = 10;
+ /* put aside some memory for segment structures */
+ segs_mb = globals->mb * 0.1;
+ if (segs_mb > 10)
+ segs_mb = 10;
- /* calculate number of region stats that can be kept in memory */
- reg_size_count = (globals->mb - segs_mb) / reg_size_mb;
- globals->min_reg_size = 3;
- if (reg_size_count < (double) globals->notnullcells / globals->min_reg_size) {
- globals->min_reg_size = (double) globals->notnullcells / reg_size_count;
+ /* calculate number of region stats that can be kept in memory */
+ reg_size_count = (globals->mb - segs_mb) / reg_size_mb;
+ globals->min_reg_size = 3;
+ if (reg_size_count < (double) globals->notnullcells / globals->min_reg_size) {
+ globals->min_reg_size = (double) globals->notnullcells / reg_size_count;
+ }
+ else {
+ reg_size_count = (double) globals->notnullcells / globals->min_reg_size;
+ /* recalculate segs_mb */
+ segs_mb = globals->mb - reg_size_count * reg_size_mb;
+ }
+
+ G_verbose_message(_("Regions with at least %d cells are stored in memory"),
+ globals->min_reg_size);
}
- else {
- reg_size_count = (double) globals->notnullcells / globals->min_reg_size;
- /* recalculate segs_mb */
- segs_mb = globals->mb - reg_size_count * reg_size_mb;
- }
- G_verbose_message(_("Regions with at least %d cells are stored in memory"),
- globals->min_reg_size);
-
/* calculate number of segments in memory */
if (globals->bounds_map != NULL) {
/* input bands, segment ids, bounds map */
- nseg = (1024. * 1024. * segs_mb) /
- (sizeof(DCELL) * globals->nbands * srows * scols +
- sizeof(CELL) * 4 * srows * scols);
+ if (globals->method == ORM_MS) {
+ nseg = (1024. * 1024. * segs_mb) /
+ (sizeof(DCELL) * 2 * globals->nbands * srows * scols +
+ sizeof(CELL) * 4 * srows * scols);
+ }
+ else {
+ nseg = (1024. * 1024. * segs_mb) /
+ (sizeof(DCELL) * globals->nbands * srows * scols +
+ sizeof(CELL) * 4 * srows * scols);
+ }
}
else {
/* input bands, segment ids */
- nseg = (1024. * 1024. * segs_mb) /
- (sizeof(DCELL) * globals->nbands * srows * scols +
- sizeof(CELL) * 2 * srows * scols);
+ if (globals->method == ORM_MS) {
+ nseg = (1024. * 1024. * segs_mb) /
+ (sizeof(DCELL) * 2 * globals->nbands * srows * scols +
+ sizeof(CELL) * 2 * srows * scols);
+ }
+ else {
+ nseg = (1024. * 1024. * segs_mb) /
+ (sizeof(DCELL) * globals->nbands * srows * scols +
+ sizeof(CELL) * 2 * srows * scols);
+ }
}
nseg_total = (globals->nrows / srows + (globals->nrows % srows > 0)) *
(globals->ncols / scols + (globals->ncols % scols > 0));
Modified: grass/trunk/imagery/i.segment/parse_args.c
===================================================================
--- grass/trunk/imagery/i.segment/parse_args.c 2016-07-10 16:39:23 UTC (rev 68943)
+++ grass/trunk/imagery/i.segment/parse_args.c 2016-07-11 18:52:07 UTC (rev 68944)
@@ -97,7 +97,7 @@
endt->key = "iterations";
endt->type = TYPE_INTEGER;
endt->required = NO;
- endt->answer = "20";
+ endt->answer = "50";
endt->description = _("Maximum number of iterations");
endt->guisection = _("Settings");
@@ -146,9 +146,7 @@
else
G_fatal_error("Invalid output raster name");
- /* Note: this threshold is scaled after we know more at the beginning of create_isegs() */
globals->alpha = atof(threshold->answer);
-
if (globals->alpha <= 0 || globals->alpha >= 1)
G_fatal_error(_("Threshold should be > 0 and < 1"));
@@ -241,6 +239,23 @@
globals->nrows = Rast_window_rows();
globals->ncols = Rast_window_cols();
+ if (sizeof(LARGEINT) < 8) {
+ int i;
+
+ LARGEINT intmax;
+
+ intmax = ((LARGEINT)1 << (sizeof(LARGEINT) * 8 - 2)) - 1;
+ intmax += ((LARGEINT)1 << (sizeof(LARGEINT) * 8 - 2));
+
+ globals->ncells = globals->ncols;
+ for (i = 1; i < globals->nrows; i++) {
+ if (globals->ncols > intmax - globals->ncells)
+ G_fatal_error(_("Integer overflow: too many cells in current region"));
+
+ globals->ncells += globals->ncols;
+ }
+ }
+
/* debug help */
if (outband->answer == NULL)
globals->out_band = NULL;
@@ -255,12 +270,13 @@
if (atoi(endt->answer) > 0)
globals->end_t = atoi(endt->answer);
else {
- globals->end_t = 100;
- G_warning(_("Invalid number of iterations, 100 will be used"));
+ globals->end_t = 50;
+ G_warning(_("Invalid number of iterations, %d will be used"),
+ globals->end_t);
}
}
else
- globals->end_t = 1000;
+ globals->end_t = 50;
if (mem->answer && atoi(mem->answer) > 10)
globals->mb = atoi(mem->answer);
Modified: grass/trunk/imagery/i.segment/region_growing.c
===================================================================
--- grass/trunk/imagery/i.segment/region_growing.c 2016-07-10 16:39:23 UTC (rev 68943)
+++ grass/trunk/imagery/i.segment/region_growing.c 2016-07-11 18:52:07 UTC (rev 68944)
@@ -25,6 +25,15 @@
#endif
#define MIN(a,b) ( ((a) < (b)) ? (a) : (b) )
+struct idlist
+{
+ int *ids;
+ int nids, nalloc;
+ CELL cellmax;
+};
+
+static struct idlist idlist;
+
/* internal functions */
static int merge_regions(struct ngbr_stats *, struct reg_stats *, /* Ri */
struct ngbr_stats *, struct reg_stats *, /* Rk */
@@ -45,6 +54,61 @@
static int calculate_reg_stats(int, int, struct reg_stats *,
struct globals *);
+void init_free_ids(void)
+{
+ idlist.nalloc = 10;
+ idlist.nids = 0;
+
+ idlist.ids = G_malloc(idlist.nalloc * sizeof(int));
+
+ idlist.cellmax = ((CELL)1 << (sizeof(CELL) * 8 - 2)) - 1;
+ idlist.cellmax += ((CELL)1 << (sizeof(CELL) * 8 - 2));
+
+ return;
+}
+
+void add_free_id(int id)
+{
+ if (id <= 0)
+ return;
+
+ if (idlist.nalloc <= idlist.nids) {
+ idlist.nalloc = idlist.nids + 10;
+ idlist.ids = G_realloc(idlist.ids, idlist.nalloc * sizeof(int));
+ }
+ idlist.ids[idlist.nids++] = id;
+
+ return;
+}
+
+int get_free_id(struct globals *globals)
+{
+ if (idlist.nids > 0) {
+ idlist.nids--;
+
+ return idlist.ids[idlist.nids];
+ }
+
+ if (globals->max_rid == idlist.cellmax)
+ G_fatal_error(_("Too many objects: integer overflow"));
+
+ globals->max_rid++;
+
+ return globals->max_rid;
+}
+
+void free_free_ids(void)
+{
+ if (idlist.nalloc) {
+ G_free(idlist.ids);
+ idlist.ids = NULL;
+ idlist.nalloc = 0;
+ idlist.nids = 0;
+ }
+
+ return;
+}
+
/* function used by binary tree to compare items */
static int compare_rc(const void *first, const void *second)
{
@@ -73,21 +137,39 @@
static int compare_double(double first, double second)
{
-
- /* standard comparison, gives good results */
if (first < second)
return -1;
return (first > second);
+}
- /* fuzzy comparison,
- * can give weird results if EPSILON is too large or
- * if the formula is changed because this is operating at the
- * limit of double fp precision */
- if (first < second && first + first * EPSILON < second)
- return -1;
- if (first > second && first > second + second * EPSILON)
- return 1;
- return 0;
+static int compare_sim_ngbrs(double simi, double simk, int candi, int candk,
+ struct ngbr_stats *Ri, struct ngbr_stats *Rk)
+{
+ if (simi < simk)
+ return -1;
+
+ if (simi > simk)
+ return 1;
+
+ if (Rk->count == 0 || Ri->count < Rk->count)
+ return -1;
+ if (Ri->count > Rk->count)
+ return 1;
+
+ if (candi && !candk)
+ return -1;
+
+ if (candk && !candi)
+ return 1;
+
+ if (Ri->row < Rk->row)
+ return -1;
+ if (Ri->row > Rk->row)
+ return 1;
+
+ if (Ri->col < Rk->col)
+ return -1;
+ return (Ri->col > Rk->col);
}
static int dump_Ri(struct ngbr_stats *Ri, struct reg_stats *Ri_rs, double *Ri_sim,
@@ -132,9 +214,16 @@
struct reg_stats Ri_rs, Rk_rs, Rk_bestn_rs;
double *dp;
struct NB_TREE *tmpnbtree;
+ CELL cellmax;
+ struct Cell_head cellhd;
G_verbose_message("Running region growing algorithm");
+ cellmax = ((CELL)1 << (sizeof(CELL) * 8 - 2)) - 1;
+ cellmax += ((CELL)1 << (sizeof(CELL) * 8 - 2));
+
+ init_free_ids();
+
/* init neighbor stats */
Ri.mean = G_malloc(globals->datasize);
Rk.mean = G_malloc(globals->datasize);
@@ -159,9 +248,11 @@
threshold = alpha2;
G_debug(1, "Squared threshold: %g", threshold);
- /* make the divisor a constant ? */
- divisor = globals->nrows + globals->ncols;
+ Rast_get_cellhd(globals->Ref.file[0].name, globals->Ref.file[0].mapset, &cellhd);
+ divisor = cellhd.rows + cellhd.cols;
+ /* TODO: renumber seeds */
+
while (t < globals->end_t && n_merges > 1) {
G_message(_("Processing pass %d..."), ++t);
@@ -208,11 +299,6 @@
/* get Ri's segment ID */
Segment_get(&globals->rid_seg, (void *)&Ri.id, Ri.row, Ri.col);
-
- if (Ri.id < 0)
- continue;
- if (Ri.id == 0)
- G_fatal_error("Zero segment id at row %d, col %d", Ri.row, Ri.col);
/* find segment neighbors */
/* find Ri's best neighbor, clear candidate flag */
@@ -230,15 +316,21 @@
Ri_nn = find_best_neighbor(&Ri, &Ri_rs, Ri_ngbrs,
&Rk, &Rk_rs, &Ri_similarity,
1, globals);
+
/* Rk is now complete */
G_debug(4, "Rk is now complete");
- if (Rk.id == 0) {
+ if (Rk.id < 0) {
/* this can only happen if the segment is surrounded by NULL data */
G_debug(4, "Segment had no valid neighbors");
continue;
}
+ if (compare_double(Ri_similarity, threshold) >= 0) {
+ G_debug(4, "Best neighbor is not similar enough");
+ continue;
+ }
+
if (/* !(t & 1) && */ Ri_nn == 1 &&
!(FLAG_GET(globals->candidate_flag, Rk.row, Rk.col)) &&
compare_double(Ri_similarity, threshold) == -1) {
@@ -264,16 +356,13 @@
pathflag = FALSE;
}
- /* this is slow ??? */
- if (/* t & */ 1) {
- if ((globals->nn < 8 && Rk.count <= 8) ||
- (globals->nn >= 8 && Rk.count <= globals->nn))
- candidates_only = FALSE;
- }
while (pathflag) {
pathflag = FALSE;
-
+
+ if (Rk.count <= globals->nn || Rk.count <= globals->min_segment_size)
+ candidates_only = FALSE;
+
/* optional check if Rk is candidate
* to prevent backwards merging */
if (candidates_only &&
@@ -331,11 +420,13 @@
Ri_nn -= Ri_ngbrs->count;
Ri_nn += (Rk_nn - Rk_ngbrs->count);
globals->ns.id = Rk.id;
+ globals->ns.row = Rk.row;
+ globals->ns.col = Rk.col;
nbtree_remove(Ri_ngbrs, &(globals->ns));
nbtree_init_trav(&travngbr, Rk_ngbrs);
while ((next = nbtree_traverse(&travngbr))) {
- if (!nbtree_find(Ri_ngbrs, next) && next->id != Ri.id)
+ if (!nbtree_find(Ri_ngbrs, next) && cmp_ngbr(next, &Ri) != 0)
nbtree_insert(Ri_ngbrs, next);
}
nbtree_clear(Rk_ngbrs);
@@ -358,14 +449,14 @@
search_neighbors(&Ri, &Ri_rs, Ri_ngbrs, &Ri_similarity,
&Rk, &Rk_rs, globals);
- if (Rk.id != 0 && Ri_nn > 0 &&
+ if (Rk.id >= 0 && Ri_nn > 0 &&
compare_double(Ri_similarity, threshold) == -1) {
pathflag = TRUE;
/* candidates_only:
- * FALSE: less passes, takes a bit longer, but less memory
- * TRUE: more passes, is a bit faster */
- candidates_only = FALSE;
+ * FALSE: less passes, slower, but less memory
+ * TRUE: more passes but faster */
+ candidates_only = TRUE;
}
/* else end of Ri -> Rk chain since we merged Ri and Rk
* go to next row, col */
@@ -381,16 +472,16 @@
if (Rk_nn < 2)
pathflag = FALSE;
- if (Rk.id < 1)
+ if (Rk.id < 0)
pathflag = FALSE;
- if (Rk_bestn.id == 0) {
- G_debug(4, "Rk's best neighour is zero");
+ if (Rk_bestn.id < 0) {
+ G_debug(4, "Rk's best neighour is negative");
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);
@@ -409,7 +500,7 @@
Ri_rs.id = Rk_rs.id;
Rk_rs.id = Rk_bestn_rs.id;
- Rk_bestn_rs.id = 0;
+ Rk_bestn_rs.id = -1;
Ri_rs.count = Rk_rs.count;
Rk_rs.count = Rk_bestn_rs.count;
Rk_bestn_rs.count = 0;
@@ -446,7 +537,26 @@
else
G_message(_("Segmentation converged after %d iterations"), t);
+ /* assign region IDs to remaining 0 IDs */
+ G_message(_("Assigning region IDs to remaining single-cell regions..."));
+ for (row = globals->row_min; row < globals->row_max; row++) {
+ G_percent(row - globals->row_min,
+ globals->row_max - globals->row_min, 4);
+ for (col = globals->col_min; col < globals->col_max; col++) {
+ if (!(FLAG_GET(globals->null_flag, row, col))) {
+ /* get segment id */
+ Segment_get(&globals->rid_seg, (void *) &Ri.id, row, col);
+ if (Ri.id == 0) {
+ Ri.id = get_free_id(globals);
+ Segment_put(&globals->rid_seg, (void *) &Ri.id, row, col);
+ }
+ }
+ }
+ }
+ G_percent(1, 1, 1);
+ free_free_ids();
+
/* ****************************************************************************************** */
/* final pass, ignore threshold and force a merge for small segments with their best neighbor */
/* ****************************************************************************************** */
@@ -518,7 +628,7 @@
Ri_nn = 0;
Ri_similarity = 2;
- Rk.id = 0;
+ Rk.id = -1;
if (do_merge) {
@@ -530,7 +640,7 @@
}
do_merge = 0;
- if (Rk.id != 0) {
+ if (Rk.id >= 0) {
/* merge Ri with Rk */
/* do not clear candidate flag for Rk */
merge_regions(&Ri, &Ri_rs, &Rk, &Rk_rs, 0, globals);
@@ -547,7 +657,7 @@
/* finished one pass for processing candidate pixels */
G_verbose_message("%d merges", n_merges);
}
-
+
/* free neighbor stats */
G_free(Ri.mean);
G_free(Rk.mean);
@@ -585,15 +695,14 @@
int neighbors[8][2];
struct RB_TREE *no_check_tree; /* cells already checked */
struct reg_stats *rs_found;
+ int candk, candtmp;
G_debug(4, "find_best_neighbor()");
if (Ri->id != Ri_rs->id)
G_fatal_error("Ri = %d but Ri_rs = %d", Ri->id, Ri_rs->id);
- if (Ri->id <= 0)
+ if (Ri->id < 0)
G_fatal_error("Ri is %d", Ri->id);
- if (Ri_rs->id <= 0)
- G_fatal_error("Ri_rs is %d", Ri_rs->id);
/* dynamics of the region growing algorithm
* some regions are growing fast, often surrounded by many small regions
@@ -609,8 +718,9 @@
nbtree_clear(Ri_ngbrs);
n_ngbrs = 0;
/* TODO: add size of largest region to reg_tree, use this as min */
- Rk->count = globals->ncells + 1;
- Rk->id = Rk_rs->id = 0;
+ Rk->count = Rk_rs->count = 0;
+ Rk->id = Rk_rs->id = -1;
+ candk = 0;
/* go through segment, spreading outwards from head */
rclist_init(&rilist);
@@ -659,7 +769,7 @@
(void *) &(globals->ns.id),
ngbr_rc.row, ngbr_rc.col);
- if (globals->ns.id == Ri->id) {
+ if (Ri->id > 0 && globals->ns.id == Ri->id) {
/* want to check this neighbor's neighbors */
rclist_add(&rilist, ngbr_rc.row, ngbr_rc.col);
@@ -684,10 +794,13 @@
/* globals->ns is now complete */
tempsim = (globals->calculate_similarity)(Ri, &globals->ns, globals);
+ candtmp = (FLAG_GET(globals->candidate_flag, ngbr_rc.row, ngbr_rc.col)) != 0;
- cmp = compare_double(tempsim, *sim);
+ cmp = compare_sim_ngbrs(tempsim, *sim, candtmp, candk, &globals->ns, Rk);
+
if (cmp == -1) {
*sim = tempsim;
+ candk = candtmp;
/* copy temp Rk to Rk */
Rk->row = ngbr_rc.row;
Rk->col = ngbr_rc.col;
@@ -704,28 +817,7 @@
memcpy(Rk_rs->sum, rs_found->sum,
globals->datasize);
}
- else if (cmp == 0) {
- /* resolve ties: prefer smaller regions */
- if (Rk->count > globals->ns.count) {
- /* copy temp Rk to Rk */
- Rk->row = ngbr_rc.row;
- Rk->col = ngbr_rc.col;
-
- Rk->id = rs_found->id;
- Rk->count = rs_found->count;
- memcpy(Rk->mean, rs_found->mean,
- globals->datasize);
-
- Rk_rs->id = Rk->id;
- Rk_rs->count = Rk->count;
- memcpy(Rk_rs->mean, rs_found->mean,
- globals->datasize);
- memcpy(Rk_rs->sum, rs_found->sum,
- globals->datasize);
- }
- }
-
n_ngbrs++;
nbtree_insert(Ri_ngbrs, &globals->ns);
}
@@ -738,10 +830,12 @@
/* clean up */
rbtree_destroy(no_check_tree);
+ rclist_destroy(&rilist);
return n_ngbrs;
}
+#ifdef _OR_SHAPE_
double calculate_shape(struct reg_stats *rsi, struct reg_stats *rsk,
int nshared, struct globals *globals)
{
@@ -817,8 +911,8 @@
return globals->smooth_weight * smooth + (1 - globals->smooth_weight) * compact;
}
+#endif
-
static int search_neighbors(struct ngbr_stats *Ri,
struct reg_stats *Ri_rs,
struct NB_TREE *Ri_ngbrs,
@@ -830,7 +924,7 @@
double tempsim, *dp;
struct NB_TRAV travngbr;
struct ngbr_stats *next;
- int cmp;
+ int cmp, candk, candtmp;
G_debug(4, "search_neighbors");
@@ -842,36 +936,29 @@
G_fatal_error("Ri_rs is %d", Ri_rs->id);
nbtree_init_trav(&travngbr, Ri_ngbrs);
- Rk->count = globals->ncells + 1;
- Rk->id = Rk_rs->id = 0;
+ Rk->count = 0;
+ Rk->id = Rk_rs->id = -1;
+ candk = 0;
while ((next = nbtree_traverse(&travngbr))) {
tempsim = (globals->calculate_similarity)(Ri, next, globals);
+ candtmp = (FLAG_GET(globals->candidate_flag, next->row, next->col)) != 0;
- cmp = compare_double(tempsim, *sim);
+ cmp = compare_sim_ngbrs(tempsim, *sim, candtmp, candk, next, Rk);
+
if (cmp == -1) {
*sim = tempsim;
+ candk = candtmp;
dp = Rk->mean;
*Rk = *next;
Rk->mean = dp;
memcpy(Rk->mean, next->mean, globals->datasize);
}
- else if (cmp == 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);
- }
- }
}
Rk_rs->id = Rk->id;
- if (Rk->id != 0) {
+ if (Rk->id >= 0) {
fetch_reg_stats(Rk->row, Rk->col, Rk_rs, globals);
}
@@ -899,6 +986,10 @@
G_fatal_error(_("Region ids are different"));
}
+ if (rs->id < 1) {
+ G_fatal_error(_("Region id %d is invalid"), rs->id);
+ }
+
if (rs->count == 1) {
G_warning(_("Region consists of only one cell, nothing to update"));
return rs->count;
@@ -1032,9 +1123,10 @@
G_debug(4, "merge_regions");
/* Ri ID must always be positive */
- if (Ri_rs->id < 1)
- G_fatal_error("Ri id is not positive: %d", Ri_rs->id);
- /* if Rk ID is negative (no seed), Rk count must be 1 */
+ if (Ri_rs->id < 1 && Ri_rs->count > 1)
+ G_fatal_error("Ri id is not positive: %d, but count is > 1: %d",
+ Ri_rs->id, Ri_rs->count);
+ /* if Rk ID is zero (no seed), Rk count must be 1 */
if (Rk_rs->id < 1 && Rk_rs->count > 1)
G_fatal_error("Rk id is not positive: %d, but count is > 1: %d",
Rk_rs->id, Rk_rs->count);
@@ -1058,6 +1150,12 @@
} while (n--);
if (Ri->count >= Rk->count) {
+
+ if (Ri->id == 0) {
+ Ri->id = get_free_id(globals);
+ Ri_rs->id = Ri->id;
+ Segment_put(&globals->rid_seg, (void *) &Ri->id, Ri->row, Ri->col);
+ }
if (Rk->count >= globals->min_reg_size) {
if (rgtree_find(globals->reg_tree, Rk_rs) == NULL)
@@ -1065,6 +1163,7 @@
/* remove from tree */
rgtree_remove(globals->reg_tree, Rk_rs);
}
+ add_free_id(Rk->id);
}
else {
@@ -1074,6 +1173,7 @@
/* remove from tree */
rgtree_remove(globals->reg_tree, Ri_rs);
}
+ add_free_id(Ri->id);
/* magic switch */
Ri_rs->id = Rk->id;
@@ -1093,7 +1193,19 @@
Ri->count = Ri_rs->count;
memcpy(Ri->mean, Ri_rs->mean, globals->datasize);
- if (Ri->id == Ri_rs->id) {
+ if (Rk->id == 0) {
+ /* the actual merge: change region id */
+ Segment_put(&globals->rid_seg, (void *) &Ri->id, Rk->row, Rk->col);
+
+ if (do_cand) {
+ 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--;
+ }
+ }
+ }
+ else if (Ri->id == Ri_rs->id) {
/* Ri is already updated, including candidate flags
* need to clear candidate flag for Rk and set new id */
@@ -1111,7 +1223,8 @@
}
rclist_init(&rlist);
- rclist_add(&rlist, Rk->row, Rk->col);
+ if (Rk->count > 1)
+ rclist_add(&rlist, Rk->row, Rk->col);
while (rclist_drop(&rlist, &next)) {
@@ -1137,9 +1250,9 @@
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);
+ ngbr_rc.row, ngbr_rc.col);
- if (R_id == Rk->id) {
+ if (Rk->id > 0 && R_id == Rk->id) {
/* the actual merge: change region id */
Segment_put(&globals->rid_seg, (void *) &Ri->id, ngbr_rc.row, ngbr_rc.col);
@@ -1166,7 +1279,8 @@
Segment_put(&globals->rid_seg, (void *) &Rk->id, Ri->row, Ri->col);
rclist_init(&rlist);
- rclist_add(&rlist, Ri->row, Ri->col);
+ if (Ri->count > 1)
+ rclist_add(&rlist, Ri->row, Ri->col);
while (rclist_drop(&rlist, &next)) {
@@ -1188,7 +1302,7 @@
Segment_get(&globals->rid_seg, (void *) &R_id, ngbr_rc.row, ngbr_rc.col);
- if (R_id == Ri->id) {
+ if (Ri->id > 0 && R_id == Ri->id) {
/* the actual merge: change region id */
Segment_put(&globals->rid_seg, (void *) &Rk->id, ngbr_rc.row, ngbr_rc.col);
@@ -1205,12 +1319,9 @@
if (Ri->id != Rk->id)
G_fatal_error("Ri ID should be set to Rk ID");
}
-
- if (Rk->id > 0)
- globals->n_regions--;
/* disable Rk */
- Rk->id = Rk_rs->id = 0;
+ Rk->id = Rk_rs->id = -1;
Rk->count = Rk_rs->count = 0;
/* update Ri */
@@ -1232,14 +1343,11 @@
G_debug(4, "set_candidate_flag");
- if (!(FLAG_GET(globals->candidate_flag, head->row, head->col)) != value) {
+ if ((!(FLAG_GET(globals->candidate_flag, head->row, head->col))) != value) {
G_warning(_("Candidate flag is already %s"), value ? _("set") : _("unset"));
return FALSE;
}
- rclist_init(&rlist);
- rclist_add(&rlist, head->row, head->col);
-
/* (un)set candidate flag */
if (value == TRUE) {
FLAG_SET(globals->candidate_flag, head->row, head->col);
@@ -1249,7 +1357,13 @@
FLAG_UNSET(globals->candidate_flag, head->row, head->col);
globals->candidate_count--;
}
+
+ if (head->id == 0)
+ return TRUE;
+ rclist_init(&rlist);
+ rclist_add(&rlist, head->row, head->col);
+
while (rclist_drop(&rlist, &next)) {
globals->find_neighbors(next.row, next.col, neighbors);
@@ -1267,7 +1381,7 @@
if (!(FLAG_GET(globals->null_flag, ngbr_rc.row, ngbr_rc.col))) {
- if (!(FLAG_GET(globals->candidate_flag, ngbr_rc.row, ngbr_rc.col)) == value) {
+ if ((!(FLAG_GET(globals->candidate_flag, ngbr_rc.row, ngbr_rc.col))) == value) {
Segment_get(&globals->rid_seg, (void *) &R_id, ngbr_rc.row, ngbr_rc.col);
@@ -1290,6 +1404,7 @@
}
} while (n--);
}
+ rclist_destroy(&rlist);
return TRUE;
}
@@ -1298,11 +1413,11 @@
struct globals *globals)
{
struct reg_stats *rs_found;
-
- if (rs->id <= 0)
+
+ if (rs->id < 0)
G_fatal_error("fetch_reg_stats(): invalid region id %d", rs->id);
- if ((rs_found = rgtree_find(globals->reg_tree, rs)) != NULL) {
+ if (rs->id > 0 && (rs_found = rgtree_find(globals->reg_tree, rs)) != NULL) {
memcpy(rs->mean, rs_found->mean, globals->datasize);
memcpy(rs->sum, rs_found->sum, globals->datasize);
@@ -1323,7 +1438,7 @@
G_debug(4, "calculate_reg_stats()");
- if (rs->id <= 0)
+ if (rs->id < 0)
G_fatal_error("Invalid region id %d", rs->id);
Segment_get(&globals->bands_seg, (void *)globals->bands_val,
@@ -1331,6 +1446,12 @@
rs->count = 1;
memcpy(rs->sum, globals->bands_val, globals->datasize);
+ if (rs->id == 0) {
+ memcpy(rs->mean, rs->sum, globals->datasize);
+
+ return 1;
+ }
+
if (globals->min_reg_size < 3)
ret = 1;
else if (globals->min_reg_size == 3) {
Modified: grass/trunk/imagery/i.segment/regtree.c
===================================================================
--- grass/trunk/imagery/i.segment/regtree.c 2016-07-10 16:39:23 UTC (rev 68943)
+++ grass/trunk/imagery/i.segment/regtree.c 2016-07-11 18:52:07 UTC (rev 68944)
@@ -77,6 +77,7 @@
int rgtree_insert(struct RG_TREE *tree, struct reg_stats *data)
{
assert(tree && data);
+ assert(data->id > 0);
if (tree->root == NULL) {
/* create a new root node for tree */
Modified: grass/trunk/imagery/i.segment/write_output.c
===================================================================
--- grass/trunk/imagery/i.segment/write_output.c 2016-07-10 16:39:23 UTC (rev 68943)
+++ grass/trunk/imagery/i.segment/write_output.c 2016-07-11 18:52:07 UTC (rev 68944)
@@ -1,6 +1,3 @@
-/* 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>
@@ -9,11 +6,7 @@
#include <grass/glocale.h>
#include "iseg.h"
-/* 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 write_ids(struct globals *globals)
{
int out_fd, row, col, maxid;
CELL *outbuf, rid;
@@ -64,167 +57,230 @@
Rast_command_history(&hist);
Rast_write_history(globals->out_name, &hist);
- /* write goodness of fit */
- if (globals->out_band) {
- int mean_fd;
- FCELL *meanbuf;
- double thresh, maxdev, sim, mingood;
- struct ngbr_stats Ri, Rk;
- struct Ref Ref; /* group reference list */
- DCELL **inbuf; /* buffers to store lines from each of the imagery group rasters */
- int n, *in_fd;
- struct FPRange *fp_range; /* min/max values of each input raster */
- DCELL *min, *max;
+ /* free memory */
+ Rast_free_colors(&colors);
- mean_fd = Rast_open_new(globals->out_band, FCELL_TYPE);
- meanbuf = Rast_allocate_f_buf();
+ return TRUE;
+}
- /* 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) */
+/* write goodness of fit */
+int write_gof_rg(struct globals *globals)
+{
+ int row, col;
+ int mean_fd;
+ CELL rid;
+ FCELL *meanbuf;
+ double thresh, maxdev, sim, mingood;
+ struct ngbr_stats Ri, Rk;
+ DCELL **inbuf; /* buffers to store lines from each of the imagery group rasters */
+ int n, *in_fd;
+ struct FPRange *fp_range; /* min/max values of each input raster */
+ struct Colors colors;
+ struct History hist;
+ DCELL *min, *max;
- thresh = globals->alpha * globals->alpha * globals->max_diff;
- maxdev = globals->max_diff * (1 - globals->alpha * globals->alpha);
- mingood = 1;
+ mean_fd = Rast_open_new(globals->out_band, FCELL_TYPE);
+ meanbuf = Rast_allocate_f_buf();
- /* open input bands */
- if (!I_get_group_ref(globals->image_group, &Ref))
- G_fatal_error(_("Group <%s> not found in the current mapset"),
- globals->image_group);
- if (Ref.nfiles <= 0)
- G_fatal_error(_("Group <%s> contains no raster maps"),
- globals->image_group);
+ /* 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) */
- in_fd = G_malloc(Ref.nfiles * sizeof(int));
- inbuf = (DCELL **) G_malloc(Ref.nfiles * sizeof(DCELL *));
- fp_range = G_malloc(Ref.nfiles * sizeof(struct FPRange));
- min = G_malloc(Ref.nfiles * sizeof(DCELL));
- max = G_malloc(Ref.nfiles * sizeof(DCELL));
+ thresh = globals->alpha * globals->alpha * globals->max_diff;
+ maxdev = globals->max_diff * (1 - globals->alpha * globals->alpha);
+ mingood = 1;
- 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);
+ /* open input bands */
+ in_fd = G_malloc(globals->Ref.nfiles * sizeof(int));
+ inbuf = (DCELL **) G_malloc(globals->Ref.nfiles * sizeof(DCELL *));
+ fp_range = G_malloc(globals->Ref.nfiles * sizeof(struct FPRange));
+ min = G_malloc(globals->Ref.nfiles * sizeof(DCELL));
+ max = G_malloc(globals->Ref.nfiles * sizeof(DCELL));
- /* returns -1 on error, 2 on empty range, quitting 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]);
+ G_debug(1, "Opening input rasters...");
+ for (n = 0; n < globals->Ref.nfiles; n++) {
+ inbuf[n] = Rast_allocate_d_buf();
+ in_fd[n] = Rast_open_old(globals->Ref.file[n].name, globals->Ref.file[n].mapset);
- G_debug(1, "Range for layer %d: min = %f, max = %f",
- n, min[n], max[n]);
- }
+ /* returns -1 on error, 2 on empty range, quitting either way. */
+ if (Rast_read_fp_range(globals->Ref.file[n].name, globals->Ref.file[n].mapset, &fp_range[n]) != 1)
+ G_fatal_error(_("No min/max found in raster map <%s>"),
+ globals->Ref.file[n].name);
+ Rast_get_fp_range_min_max(&(fp_range[n]), &min[n], &max[n]);
- G_message(_("Writing out goodness of fit"));
- for (row = 0; row < globals->nrows; row++) {
+ G_debug(1, "Range for layer %d: min = %f, max = %f",
+ n, min[n], max[n]);
+ }
- G_percent(row, globals->nrows, 9);
+ G_message(_("Writing out goodness of fit"));
+ for (row = 0; row < globals->nrows; row++) {
- Rast_set_f_null_value(meanbuf, globals->ncols);
+ G_percent(row, globals->nrows, 9);
- for (n = 0; n < Ref.nfiles; n++) {
- Rast_get_d_row(in_fd[n], inbuf[n], row);
- }
+ Rast_set_f_null_value(meanbuf, globals->ncols);
- for (col = 0; col < globals->ncols; col++) {
+ for (n = 0; n < globals->Ref.nfiles; n++) {
+ Rast_get_d_row(in_fd[n], inbuf[n], row);
+ }
- if (!(FLAG_GET(globals->null_flag, row, col))) {
+ 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) {
- Segment_get(&globals->rid_seg, (void *) &rid, row, col);
+ Ri.row = Rk.row = row;
+ Ri.col = Rk.col = col;
- if (rid > 0) {
-
- Ri.row = Rk.row = row;
- Ri.col = Rk.col = col;
+ /* get values for Ri = this region */
+ globals->rs.id = rid;
+ fetch_reg_stats(row, col, &globals->rs, globals);
+ Ri.mean = globals->rs.mean;
+ Ri.count = globals->rs.count;
- /* get values for Ri = this region */
- globals->rs.id = rid;
- fetch_reg_stats(row, col, &globals->rs, globals);
- Ri.mean = globals->rs.mean;
- Ri.count = globals->rs.count;
+ sim = 0.;
+ /* region consists of more than one cell */
+ if (Ri.count > 1) {
- sim = 0.;
- /* region consists of more than one cell */
- if (Ri.count > 1) {
+ /* get values for Rk = this cell */
+ for (n = 0; n < globals->Ref.nfiles; n++) {
+ if (globals->weighted == FALSE)
+ /* scaled version */
+ globals->second_val[n] = (inbuf[n][col] - min[n]) / (max[n] - min[n]);
+ else
+ globals->second_val[n] = inbuf[n][col];
+ }
- /* get values for Rk = this cell */
- for (n = 0; n < Ref.nfiles; n++) {
- if (globals->weighted == FALSE)
- /* scaled version */
- globals->second_val[n] = (inbuf[n][col] - min[n]) / (max[n] - min[n]);
- else
- globals->second_val[n] = inbuf[n][col];
- }
+ Rk.mean = globals->second_val;
- Rk.mean = globals->second_val;
-
- /* calculate similarity */
- sim = (*globals->calculate_similarity) (&Ri, &Rk, globals);
- }
-
- if (0) {
- if (sim < thresh)
- meanbuf[col] = 1;
- else {
- sim = 1. - (sim - thresh) / maxdev;
- meanbuf[col] = sim;
- if (mingood > sim)
- mingood = sim;
- }
- }
+ /* calculate similarity */
+ sim = (*globals->calculate_similarity) (&Ri, &Rk, globals);
+ }
+
+ if (0) {
+ if (sim < thresh)
+ meanbuf[col] = 1;
else {
- sim = 1 - sim;
+ sim = 1. - (sim - thresh) / maxdev;
meanbuf[col] = sim;
if (mingood > sim)
mingood = sim;
}
}
+ else {
+ sim = 1 - sim;
+ meanbuf[col] = sim;
+ if (mingood > sim)
+ mingood = sim;
+ }
}
}
- Rast_put_row(mean_fd, meanbuf, FCELL_TYPE);
}
+ Rast_put_row(mean_fd, meanbuf, FCELL_TYPE);
+ }
- Rast_close(mean_fd);
+ Rast_close(mean_fd);
- Rast_init_colors(&colors);
- Rast_make_grey_scale_fp_colors(&colors, mingood, 1);
- Rast_write_colors(globals->out_band, G_mapset(), &colors);
+ Rast_init_colors(&colors);
+ Rast_make_grey_scale_fp_colors(&colors, mingood, 1);
+ Rast_write_colors(globals->out_band, G_mapset(), &colors);
- Rast_short_history(globals->out_band, "raster", &hist);
- Rast_command_history(&hist);
- Rast_write_history(globals->out_band, &hist);
+ Rast_short_history(globals->out_band, "raster", &hist);
+ Rast_command_history(&hist);
+ Rast_write_history(globals->out_band, &hist);
- G_free(meanbuf);
+ G_free(meanbuf);
- G_debug(1, "Closing input rasters...");
- for (n = 0; n < Ref.nfiles; n++) {
- Rast_close(in_fd[n]);
- G_free(inbuf[n]);
+ G_debug(1, "Closing input rasters...");
+ for (n = 0; n < globals->Ref.nfiles; n++) {
+ Rast_close(in_fd[n]);
+ G_free(inbuf[n]);
+ }
+
+ G_free(inbuf);
+ G_free(in_fd);
+ G_free(fp_range);
+ G_free(min);
+ G_free(max);
+
+ return TRUE;
+}
+
+int write_bands_ms(struct globals *globals)
+{
+ int *out_fd, row, col, n;
+ DCELL **outbuf;
+ char **name;
+ struct Colors colors;
+ struct History hist;
+ struct ngbr_stats Rout;
+
+ out_fd = G_malloc(sizeof(int) * globals->nbands);
+ name = G_malloc(sizeof(char *) * globals->nbands);
+ outbuf = G_malloc(sizeof(DCELL) * globals->nbands);
+ for (n = 0; n < globals->nbands; n++) {
+ outbuf[n] = Rast_allocate_d_buf();
+ G_asprintf(&name[n], "%s%s", globals->Ref.file[n].name, globals->ms_suf);
+ out_fd[n] = Rast_open_new(name[n], DCELL_TYPE);
+ }
+
+ Rout.mean = G_malloc(globals->datasize);
+
+
+ for (row = 0; row < globals->nrows; row++) {
+
+ G_percent(row, globals->nrows, 9);
+
+ for (n = 0; n < globals->nbands; n++)
+ Rast_set_d_null_value(outbuf[n], globals->ncols);
+ for (col = 0; col < globals->ncols; col++) {
+
+ if (!(FLAG_GET(globals->null_flag, row, col))) {
+ Segment_get(globals->bands_out, (void *) Rout.mean, row, col);
+
+ for (n = 0; n < globals->nbands; n++) {
+ outbuf[n][col] = Rout.mean[n];
+ if (globals->weighted == FALSE)
+ outbuf[n][col] = Rout.mean[n] * (globals->max[n] - globals->min[n]) + globals->min[n];
+ }
+ }
}
- G_free(inbuf);
- G_free(in_fd);
- G_free(fp_range);
- G_free(min);
- G_free(max);
+ for (n = 0; n < globals->nbands; n++)
+ Rast_put_row(out_fd[n], outbuf[n], DCELL_TYPE);
}
- /* free memory */
- Rast_free_colors(&colors);
+ for (n = 0; n < globals->nbands; n++) {
+ Rast_close(out_fd[n]);
+ Rast_read_colors(globals->Ref.file[n].name, globals->Ref.file[n].mapset, &colors);
+ Rast_write_colors(name[n], G_mapset(), &colors);
+
+ Rast_short_history(name[n], "raster", &hist);
+ Rast_command_history(&hist);
+ Rast_write_history(name[n], &hist);
+ }
+
+ /* free */
+
return TRUE;
}
int close_files(struct globals *globals)
{
+ G_debug(1, "Closing input rasters...");
+
/* close segmentation files and output raster */
G_debug(1, "closing files");
Segment_close(&globals->bands_seg);
+ if (globals->method == ORM_MS)
+ Segment_close(&globals->bands_seg2);
if (globals->bounds_map)
Segment_close(&globals->bounds_seg);
More information about the grass-commit
mailing list