move i.segment.xl to trunk

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

imagery/i.segment/Makefile

@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.segment.xl
+
+LIBES = $(IMAGERYLIB) $(RASTERLIB) $(SEGMENTLIB) $(GISLIB) $(BTREE2LIB)
+DEPENDENCIES = $(IMAGERYDEP) $(RASTERDEP) $(SEGMENTDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

imagery/i.segment/flag.c

@@ -0,0 +1,60 @@
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "flag.h"
+
+FLAG *flag_create(int nrows, int ncols)
+{
+    unsigned char *temp;
+
+    FLAG *new_flag;
+
+    register int i;
+
+    new_flag = (FLAG *) G_malloc(sizeof(FLAG));
+    new_flag->nrows = nrows;
+    new_flag->ncols = ncols;
+    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);
+}
+
+int flag_destroy(FLAG * flags)
+{
+    G_free(flags->array[0]);
+    G_free(flags->array);
+    G_free(flags);
+
+    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;
+}

imagery/i.segment/flag.h

@@ -0,0 +1,64 @@
+#ifndef __FLAG_H__
+#define __FLAG_H__
+
+
+/* flag.[ch] is a set of routines which will set up an array of bits
+ ** that allow the programmer to "flag" cells in a raster map.
+ **
+ ** FLAG *
+ ** flag_create(nrows,ncols)
+ ** int nrows, ncols;
+ **     opens the structure flag.  
+ **     The flag structure will be a two dimensional array of bits the
+ **     size of nrows by ncols.  Will initalize flags to zero (unset).
+ **
+ ** flag_destroy(flags)
+ ** FLAG *flags;
+ **     closes flags and gives the memory back to the system.
+ **
+ ** flag_clear_all(flags)
+ ** FLAG *flags;
+ **     sets all values in flags to zero.
+ **
+ ** flag_unset(flags, row, col)
+ ** FLAG *flags;
+ ** int row, col;
+ **     sets the value of (row, col) in flags to zero.
+ **
+ ** flag_set(flags, row, col)
+ ** FLAG *flags;
+ ** int row, col;
+ **     will set the value of (row, col) in flags to one.
+ **
+ ** int
+ ** flag_get(flags, row, col)
+ ** FLAG *flags;
+ ** int row, col;
+ **     returns the value in flags that is at (row, col).
+ **
+ ** idea by Michael Shapiro
+ ** code by Chuck Ehlschlaeger
+ ** April 03, 1989
+ */
+#define FLAG struct _flagsss_
+FLAG {
+    int nrows, ncols, leng;
+    unsigned char **array;
+};
+
+#define FLAG_UNSET(flags,row,col) \
+	(flags)->array[(row)][(col)>>3] &= ~(1<<((col) & 7))
+
+#define FLAG_SET(flags,row,col) \
+	(flags)->array[(row)][(col)>>3] |= (1<<((col) & 7))
+
+#define FLAG_GET(flags,row,col) \
+	(flags)->array[(row)][(col)>>3] & (1<<((col) & 7))
+
+/* flag.c */
+FLAG *flag_create(int, int);
+int flag_destroy(FLAG *);
+int flag_clear_all(FLAG *);
+
+
+#endif /* __FLAG_H__ */

imagery/i.segment/i.segment.xl.html

@@ -0,0 +1,179 @@
+<h2>DESCRIPTION</h2>
+Image segmentation or object recognition is the process of grouping 
+similar pixels into unique segments, also refered to as objects. 
+Boundary and region based algorithms are described in the literature, 
+currently a region growing and merging algorithm is implemented. Each 
+object found during the segmentation process is given a unique ID and 
+is a collection of contiguous pixels meeting some criteria. Note the 
+contrast with image classification where all pixels similar to each 
+other are assigned to the same class and do not need to be contiguous. 
+The image segmentation results can be useful on their own, or used as a 
+preprocessing step for image classification. The segmentation 
+preprocessing step can reduce noise and speed up the classification.
+
+<H2>NOTES</h2>
+
+<h3>Region Growing and Merging</h3>
+This segmentation algorithm sequentially examines all current 
+segments in the map. The similarity between the current segment and 
+each of its neighbors is calculated according to the given distance 
+formula. Segments will be merged if they meet a number of criteria, 
+including: 1. The pair is mutually most similar to each other (the 
+similarity distance will be smaller than to any other neighbor), and 
+2. The similarity must be lower than the input threshold. The process 
+is repeated until no merges are made during a complete pass.
+
+<h3>Similarity and Threshold</h3>
+The similarity between segments and unmerged objects is used to 
+determine which objects are merged. Smaller distance values indicate a 
+closer match, with a similarity score of zero for identical pixels.
+<p>
+During normal processing, merges are only allowed when the 
+similarity between two segments is lower than the givem 
+threshold value. During the final pass, however, if a minimum 
+segment size of 2 or larger is given with the <em>minsize</em> 
+parameter, segments with a smaller pixel count will be merged with 
+their most similar neighbor even if the similarity is greater than 
+the threshold.
+<p>
+The threshold should be set by the user between 0 and 1.0. A threshold 
+of 0 would allow only identical valued pixels to be merged, while a 
+threshold of 1 would allow everything to be merged. Initial empirical 
+tests indicate threshold values of 0.01 to 0.05 are reasonable values 
+to start.
+
+<h4>Calculation Formulas</h4>
+Both Euclidean and Manhattan distances use the normal definition, 
+considering each raster in the image group as a dimension.
+
+In future , the distance calculation will also take into account the 
+shape characteristics of the segments. The normal distances are then
+multiplied by the input radiometric weight. Next an additional 
+contribution is added: (1-radioweight) * {smoothness * smoothness 
+weight + compactness * (1-smoothness weight)}, where compactness = 
+the Perimeter Length / sqrt( Area ) and smoothness = Perimeter 
+Length / the Bounding Box. The perimeter length is estimated as the 
+number of pixel sides the segment has.
+
+<h3>Seeds</h3>
+The seeds map can be used to provide either seed pixels (random or 
+selected points from which to start the segmentation process) or 
+seed segments (results of previous segmentations or 
+classifications). The different approaches are automatically 
+detected by the program: any pixels that have identical seed values 
+and are contiguous will be assigned a unique segment ID.
+<p>
+It is expected that the <em>minsize</em> will be set to 1 if a seed 
+map is used, but the program will allow other values to be used.  If 
+both options are used, the final iteration that ignores the 
+threshold also will ignore the seed map and force merges for all 
+pixels (not just segments that have grown/merged from the seeds).
+
+<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 IDs.  If the initial map has a larger 
+number of non-null pixels, there are two workarounds:
+<p>
+1.  Use starting seed pixels. (Maximum 2 billion pixels can be 
+labeled with positive integers.)
+<p>
+2.  Use starting seed segments. (By initial classification or other 
+methods.)
+
+<h3>Boundary Constraints</h3>
+Boundary constraints limit the adjacency of pixels and segments.  
+Each unique value present in the <em>bounds</em> raster are 
+considered as a MASK. Thus no segments in the final segmentated map 
+will cross a boundary, even if their spectral data is very similar.
+
+<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>
+<div class="code"><pre>
+i.group group=ortho_group input=ortho_2001_t792_1m@PERMANENT
+</pre></div>
+
+<p>Because the segmentation process is computationally expensive, 
+start with a small processing area to confirm if the segmentation 
+results meet your requirements.  Some manual adjustment of the 
+threshold may be required. <br>
+
+<div class="code"><pre>
+g.region rast=ortho_2001_t792_1m@PERMANENT
+</pre></div>
+
+Try out a first threshold and check the results.<br>
+<div class="code"><pre>
+i.segment -w group=ortho_group output=ortho_segs threshold=0.04 \
+          method=region_growing 
+</pre></div>
+<p></p>
+From a visual inspection, it seems this results in oversegmentation.  
+Increasing the threshold: <br>
+<div class="code"><pre>
+i.segment -w group=ortho_group output=ortho_segs \
+          threshold=0.1 method=region_growing
+</pre></div>
+<p></p>
+This looks better.  There is some noise in the image, lets next force 
+all segments smaller than 5 pixels to be merged into their most similar 
+neighbor (even if they are less similar then required by our 
+threshold):<br>
+<div class="code"><pre>
+i.segment -w --overwrite group=ortho_group output=ortho_segs \
+          threshold=0.1 method=region_growing minsize=5
+</pre></div>
+<p></p>
+Processing the entire ortho image with nearly 10 million pixels took about 
+15 minutes.
+
+<h2>TODO</h2>
+<h3>Functionality</h3>
+<ul>
+<li>Further testing of the shape characteristics (smoothness, 
+compactness), if it looks good it should be added.
+<b>in progress</b></li>
+<li>Malahanobis distance for the similarity calculation.</li>
+</ul>
+<h3>Use of Segmentation Results</h3>
+<ul>
+<li>Improve the optional output from this module, or better yet, add a 
+module for <em>i.segment.metrics</em>.</li>
+<li>Providing updates to i.maxlik to ensure the segmentation output can 
+be used as input for the existing classification functionality.</li>
+<li>Integration/workflow for <em>r.fuzzy</em>.</li>
+</ul>
+<h3>Speed</h3>
+<ul>
+<li>See create_isegs.c</li>
+</ul>
+<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<br>
+Markus Metz (GSoC Mentor)
+

imagery/i.segment/iseg.h

@@ -0,0 +1,134 @@
+
+/****************************************************************************
+ *
+ * MODULE:       i.segment
+ * AUTHOR(S):    Eric Momsen <eric.momsen at gmail com>
+ * PURPOSE:      structure definition and function listing
+ * COPYRIGHT:    (C) 2012 by Eric Momsen, and the GRASS Development Team
+ *
+ *               This program is free software under the GNU General Public
+ *               License (>=v2). Read the COPYING file that comes with GRASS
+ *               for details.
+ *
+ *****************************************************************************/
+
+#include <grass/segment.h>
+#include "flag.h"
+#include "regtree.h"
+#include "ngbrtree.h"
+
+/* #def _OR_SHAPE_ */
+
+
+/* row/col list */
+struct rc
+{
+    struct rc *next;
+    int row;
+    int col;
+};
+
+struct rclist
+{
+    struct rc *tail, *head;
+};
+
+/* input and output files, as well as some other processing info */
+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 method;			/* Segmentation method */
+    int nn;			/* number of neighbors, 4 or 8 */
+    double max_diff;		/* max possible difference */
+    double alpha;		/* similarity threshold */
+    int min_segment_size;	/* smallest number of pixels/cells allowed in a final segment */
+
+    double radio_weight;	/* weighing factor radiometric - shape */
+    double smooth_weight;       /* weighing factor smoothness - compactness */
+
+    int end_t;			/* maximum number of iterations */
+    int mb;
+
+    /* region info */
+    int nrows, ncols;
+    int row_min, row_max, col_min, col_max; /* region constraints */
+    int ncells;
+
+    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 */
+            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 */
+    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 */
+    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 */
+};
+
+
+/* parse_args.c */
+/* gets input from user, validates, and sets up functions */
+int parse_args(int, char *[], struct globals *);
+
+/* open_files.c */
+int open_files(struct globals *);
+
+/* create_isegs.c */
+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 *);
+double calculate_manhattan_similarity(struct ngbr_stats *, 
+                                      struct ngbr_stats *,
+				      struct globals *);
+double calculate_shape(struct reg_stats *, struct reg_stats *,
+                       int, struct globals *);
+int fetch_reg_stats(int , int , struct reg_stats *, 
+                           struct globals *);
+
+/* void calculate_reg_stats(int, int, struct reg_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 globals *);
+int close_files(struct globals *);

imagery/i.segment/main.c

@@ -0,0 +1,60 @@
+
+/****************************************************************************
+ *
+ * MODULE:       i.segment
+ * AUTHOR(S):    Eric Momsen <eric.momsen at gmail com>
+ * PURPOSE:      Segments an image group.
+ * COPYRIGHT:    (C) 2012 by Eric Momsen, and the GRASS Development Team
+ *
+ *               This program is free software under the GNU General Public
+ *               License (>=v2). Read the COPYING file that comes with GRASS
+ *               for details.
+ * 
+ *
+ *               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.
+ * 
+ * 
+ *****************************************************************************/
+
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "iseg.h"
+
+int main(int argc, char *argv[])
+{
+    struct globals globals;		/* input and output file descriptors, data structure, buffers */
+    struct GModule *module;
+
+    G_gisinit(argv[0]);
+
+    module = G_define_module();
+    G_add_keyword(_("imagery"));
+    G_add_keyword(_("segmentation"));
+    module->description =
+	_("Outputs a single segmented map (raster) based on input values in an image group.");
+
+    if (parse_args(argc, argv, &globals) != TRUE)
+	G_fatal_error(_("Error in parse_args()"));
+
+    G_debug(1, "Main: starting open_files()");
+    if (open_files(&globals) != TRUE)
+	G_fatal_error(_("Error in open_files()"));
+
+    G_debug(1, "Main: starting create_isegs()");
+    if (create_isegs(&globals) != TRUE)
+	G_fatal_error(_("Error in create_isegs()"));
+
+    G_debug(1, "Main: starting write_output()");
+    if (write_output(&globals) != TRUE)
+	G_fatal_error(_("Error in write_output()"));
+
+    G_debug(1, "Main: starting close_files()");
+    close_files(&globals);
+
+    G_done_msg(_("Number of segments created: %d"), globals.n_regions);
+
+    exit(EXIT_SUCCESS);
+}

imagery/i.segment/ngbrtree.c

@@ -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);
+}

imagery/i.segment/ngbrtree.h

@@ -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

imagery/i.segment/open_files.c

@@ -0,0 +1,524 @@
+/* PURPOSE:      opening input rasters and creating segmentation files */
+
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/imagery.h>
+#include <grass/segment.h>	/* segmentation library */
+#include "iseg.h"
+
+static int load_seeds(struct globals *, int, int, int);
+static int read_seed(struct globals *, SEGMENT *, struct rc *, int);
+static int manage_memory(int, int, struct globals *);
+
+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_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 ngbr_stats Ri, Rk;
+
+    /*allocate memory for flags */
+    globals->null_flag = flag_create(globals->nrows, globals->ncols);
+    globals->candidate_flag = flag_create(globals->nrows, globals->ncols);
+
+    G_debug(1, "Checking image group...");
+
+    /* ****** open the input rasters ******* */
+
+    if (!I_get_group_ref(globals->image_group, &Ref))
+	G_fatal_error(_("Unable to read REF file for group <%s>"),
+		      globals->image_group);
+
+    if (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));
+
+    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);
+    }
+
+    /* Get min/max values of each input raster for scaling */
+
+    globals->max_diff = 0.;
+
+    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]);
+
+	G_debug(1, "Range for layer %d: min = %f, max = %f",
+		    n, min[n], max[n]);
+	
+    }
+    if (globals->weighted == FALSE)
+	globals->max_diff = Ref.nfiles;
+    else {
+	/* max difference with selected similarity method */
+	Ri.mean = max;
+	Rk.mean = min;
+	globals->max_diff = 1;
+	globals->max_diff = (*globals->calculate_similarity) (&Ri, &Rk, globals);
+    }
+
+    /* ********** find out file segmentation size ************ */
+    G_debug(1, "Calculate temp file sizes...");
+
+    globals->nbands = Ref.nfiles;
+
+    /* size of each element to be stored */
+
+    inlen = sizeof(DCELL) * Ref.nfiles;
+    outlen = sizeof(CELL);
+    G_debug(1, "data element size, in: %d , out: %d ", inlen, outlen);
+    globals->datasize = sizeof(double) * globals->nbands;
+
+    /* segment lib segment size */
+    srows = 64;
+    scols = 64;
+
+    nseg = manage_memory(srows, scols, globals);
+
+    /* create segment structures */
+    if (segment_open
+	(&globals->bands_seg, G_tempfile(), globals->nrows, globals->ncols, srows,
+	 scols, inlen, nseg) != 1)
+	G_fatal_error("Unable to create input temporary files");
+
+    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 */
+    G_debug(1, "Reading input rasters into segmentation data files...");
+
+    globals->bands_val = (double *)G_malloc(inlen);
+    globals->second_val = (double *)G_malloc(inlen);
+    /* initial segment ID */
+    s = 1;
+
+    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++) {
+
+	    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
+	    (&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();
+
+	for (row = 0; row < globals->nrows; row++) {
+	    Rast_get_c_row(bounds_fd, boundsbuf, row);
+	    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);
+
+	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.");
+    }
+
+    /* other info */
+    globals->candidate_count = 0;	/* counter for remaining candidate pixels */
+
+    /* Free memory */
+
+    for (n = 0; n < Ref.nfiles; n++) {
+	G_free(inbuf[n]);
+	Rast_close(in_fd[n]);
+    }
+
+    globals->rs.sum = G_malloc(globals->datasize);
+    globals->rs.mean = G_malloc(globals->datasize);
+
+    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()");
+    
+    G_message(_("Loading seeds from '%s'"), globals->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 {
+		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);
+	}
+    }
+
+    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++;
+		}
+	    }
+	}
+    }
+
+    G_free(seeds_buf);
+    Rast_close(seeds_fd);
+    segment_close(&seeds_seg);
+
+    globals->n_regions = spos - 1;
+    
+    flag_clear_all(globals->candidate_flag);
+    
+    return 1;
+}
+
+
+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];
+
+    G_debug(4, "read_seed()");
+
+    /* 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);
+
+    /* initialize region stats */
+    globals->rs.count = 1;
+
+    globals->rs.id = new_id;
+    segment_get(&globals->bands_seg, (void *)globals->bands_val,
+		Ri->row, Ri->col);
+
+    for (i = 0; i < globals->nbands; i++) {
+	globals->rs.sum[i] = globals->bands_val[i];
+	globals->rs.mean[i] = globals->bands_val[i];
+    }
+
+    /* 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;
+	    }
+
+	    if (FLAG_GET(globals->null_flag, ngbr_rc.row, ngbr_rc.col)) {
+		continue;
+	    }
+
+	    if (FLAG_GET(globals->candidate_flag, ngbr_rc.row, ngbr_rc.col)) {
+		continue;
+	    }
+
+	    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 (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 >= globals->min_reg_size) {
+	for (i = 0; i < globals->nbands; i++)
+	    globals->rs.mean[i] = globals->rs.sum[i] / globals->rs.count;
+
+	rgtree_insert(globals->reg_tree, &(globals->rs));
+    }
+
+    return 1;
+}
+
+
+static int manage_memory(int srows, int scols, struct globals *globals)
+{
+    double reg_size_mb, segs_mb;
+    int reg_size_count, nseg, nseg_total;
+
+    /* 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;
+
+    /* 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 = 4;
+    if (reg_size_count < (double) globals->nrows * globals->ncols / globals->min_reg_size) {
+	globals->min_reg_size = (double) globals->nrows * globals->ncols / reg_size_count;
+    }
+    else {
+	reg_size_count = (double) globals->nrows * globals->ncols / globals->min_reg_size;
+	/* recalculate segs_mb */
+	segs_mb = globals->mb - reg_size_count * reg_size_mb;
+    }
+
+    G_verbose_message(_("Stats for 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);
+    }
+    else {
+	/* input bands, segment ids, bounds map */
+	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));
+
+    if (nseg > nseg_total)
+	nseg = nseg_total;
+    
+    G_debug(1, "current region:  %d rows, %d cols", globals->nrows, globals->ncols);
+    G_debug(1, "segmented to tiles with size:  %d rows, %d cols", srows,
+	    scols);
+    G_verbose_message(_("Number of segments in memory: %d of %d total"),
+                      nseg, nseg_total);
+    
+    return nseg;
+}

imagery/i.segment/parse_args.c

@@ -0,0 +1,259 @@
+/* PURPOSE:      Parse and validate the input */
+
+#include <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/raster.h>
+#include "iseg.h"
+
+int parse_args(int argc, char *argv[], struct globals *globals)
+{
+    struct Option *group, *seeds, *bounds, *output,
+                  *method, *similarity, *threshold, *min_segment_size,
+#ifdef _OR_SHAPE_
+		  *shape_weight, *smooth_weight,
+#endif
+		   *mem;
+    struct Flag *diagonal, *weighted;
+    struct Option *outband, *endt;
+
+    /* required parameters */
+    group = G_define_standard_option(G_OPT_I_GROUP);
+
+    output = G_define_standard_option(G_OPT_R_OUTPUT);
+
+    threshold = G_define_option();
+    threshold->key = "threshold";
+    threshold->type = TYPE_DOUBLE;
+    threshold->required = YES;
+    threshold->label = _("Difference threshold between 0 and 1.");
+    threshold->description = _("Threshold = 0 merges only identical segments; threshold = 1 merges all.");
+
+    /* optional parameters */
+
+    method = G_define_option();
+    method->key = "method";
+    method->type = TYPE_STRING;
+    method->required = NO;
+    method->answer = "region_growing";
+    method->description = _("Segmentation method.");
+    method->guisection = _("Settings");
+
+    similarity = G_define_option();
+    similarity->key = "similarity";
+    similarity->type = TYPE_STRING;
+    similarity->required = YES;
+    similarity->answer = "euclidean";
+    similarity->options = "euclidean, manhattan";
+    similarity->description = _("Similarity calculation method.");
+    similarity->guisection = _("Settings");
+
+    min_segment_size = G_define_option();
+    min_segment_size->key = "min";
+    min_segment_size->type = TYPE_INTEGER;
+    min_segment_size->required = NO;
+    min_segment_size->answer = "1";
+    min_segment_size->options = "1-100000";
+    min_segment_size->label = _("Minimum number of cells in a segment.");
+    min_segment_size->description =
+	_("The final step will merge small segments with their best neighbor.");
+    min_segment_size->guisection = _("Settings");
+
+#ifdef _OR_SHAPE_
+    radio_weight = G_define_option();
+    radio_weight->key = "radio_weight";
+    radio_weight->type = TYPE_DOUBLE;
+    radio_weight->required = YES;
+    radio_weight->answer = "1";
+    radio_weight->options = "0-1";
+    radio_weight->label =
+	_("Importance of radiometric (input raster) values relative to shape.");
+    radio_weight->guisection = _("Settings");
+
+    smooth_weight = G_define_option();
+    smooth_weight->key = "smooth_weight";
+    smooth_weight->type = TYPE_DOUBLE;
+    smooth_weight->required = YES;
+    smooth_weight->answer = "0.5";
+    smooth_weight->options = "0-1";
+    smooth_weight->label =
+	_("Importance of smoothness relative to compactness.");
+    smooth_weight->guisection = _("Settings");
+#endif
+
+    /* 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->description = _("Optional raster map with starting seeds.");
+
+    /* Polygon constraints. */
+    bounds = G_define_standard_option(G_OPT_R_INPUT);
+    bounds->key = "bounds";
+    bounds->required = NO;
+    bounds->label = _("Optional bounding/constraining raster map");
+    bounds->description =
+	_("Must be integer values, each area will be segmented independent of the others.");
+
+    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 = "iterations";
+    endt->type = TYPE_INTEGER;
+    endt->required = NO;
+    endt->answer = "100";
+    endt->description = _("Maximum number of iterations.");
+
+    outband = G_define_standard_option(G_OPT_R_OUTPUT);
+    outband->key = "goodness";
+    outband->required = NO;
+    outband->description =
+	_("Goodness of fit estimate.");
+
+    diagonal = G_define_flag();
+    diagonal->key = 'd';
+    diagonal->description =
+	_("Use 8 neighbors (3x3 neighborhood) instead of the default 4 neighbors for each pixel.");
+
+    weighted = G_define_flag();
+    weighted->key = 'w';
+    weighted->description =
+	_("Weighted input, don't perform the default scaling of input maps.");
+
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    /* Check and save parameters */
+
+    globals->image_group = group->answer;
+
+    if (G_legal_filename(output->answer) == TRUE)
+	globals->out_name = output->answer;
+    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"));
+
+    /* segmentation methods:  1 = region growing */
+    if (strcmp(method->answer, "region_growing") == 0)
+	globals->method = 1;
+    else
+	G_fatal_error("Couldn't assign segmentation method.");
+
+    G_debug(1, "segmentation method: %d", globals->method);
+
+    /* distance methods for similarity measurement */
+    if (strcmp(similarity->answer, "euclidean") == 0)
+	globals->calculate_similarity = calculate_euclidean_similarity;
+    else if (strcmp(similarity->answer, "manhattan") == 0)
+	globals->calculate_similarity = calculate_manhattan_similarity;
+    else
+	G_fatal_error(_("Invalid similarity method."));
+
+#ifdef _OR_SHAPE_
+    /* consider shape */
+    globals->radio_weight = atof(radio_weight->answer);
+    if (globals->radio_weight <= 0)
+	G_fatal_error(_("Option '%s' must be > 0"), radio_weight->key);
+    if (globals->radio_weight > 1)
+	G_fatal_error(_("Option '%s' must be <= 1"), radio_weight->key);
+    globals->smooth_weight = atof(smooth_weight->answer);
+    if (globals->smooth_weight < 0)
+	G_fatal_error(_("Option '%s' must be >= 0"), smooth_weight->key);
+    if (globals->smooth_weight > 1)
+	G_fatal_error(_("Option '%s' must be <= 1"), smooth_weight->key);
+#else
+    globals->radio_weight = 1;
+    globals->smooth_weight = 0.5;
+#endif
+
+    globals->min_segment_size = atoi(min_segment_size->answer);
+
+    if (diagonal->answer == FALSE) {
+	globals->find_neighbors = find_four_neighbors;
+	globals->nn = 4;
+	G_debug(1, "four pixel neighborhood");
+    }
+    else if (diagonal->answer == TRUE) {
+	globals->find_neighbors = find_eight_neighbors;
+	globals->nn = 8;
+	G_debug(1, "eight (3x3) pixel neighborhood");
+    }
+
+    /* default/0 for performing the scaling
+     * selected/1 if scaling should be skipped. */
+    globals->weighted = weighted->answer;
+
+    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(globals->seeds, "") !=
+	    CELL_TYPE) {
+	    G_fatal_error(_("Seeeds map must be CELL type (integers)"));
+	}
+    }
+
+    if (bounds->answer == NULL) {
+	globals->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(globals->bounds_map, globals->bounds_mapset) !=
+	    CELL_TYPE) {
+	    G_fatal_error(_("Segmentation constraint map must be CELL type (integers)"));
+	}
+    }
+
+    /* other data */
+    globals->nrows = Rast_window_rows();
+    globals->ncols = Rast_window_cols();
+
+    /* debug help */
+    if (outband->answer == NULL)
+	globals->out_band = NULL;
+    else {
+	if (G_legal_filename(outband->answer) == TRUE)
+	    globals->out_band = outband->answer;
+	else
+	    G_fatal_error("Invalid output raster name for goodness of fit.");
+    }
+
+    if (endt->answer) {
+	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."));
+	}
+    }
+    else
+	globals->end_t = 1000;
+
+    if (mem->answer && atoi(mem->answer) > 10)
+	globals->mb = atoi(mem->answer);
+    else {
+	globals->mb = 300;
+	G_warning(_("Invalid number of MB, 300 will be used."));
+    }
+
+    return TRUE;
+}

imagery/i.segment/rclist.c

@@ -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;
+    }
+
+    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;
+}
+

imagery/i.segment/regtree.c

@@ -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);
+}

imagery/i.segment/regtree.h

@@ -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

imagery/i.segment/write_output.c

@@ -0,0 +1,179 @@
+/* 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"
+
+/* 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, row, col;
+    CELL *outbuf, rid;
+    struct Colors colors;
+    struct History hist;
+
+    outbuf = Rast_allocate_c_buf();
+
+    G_debug(1, "preparing output raster");
+    /* open output raster map */
+    out_fd = Rast_open_new(globals->out_name, CELL_TYPE);
+
+    G_debug(1, "start data transfer from segmentation file to raster");
+
+    G_message(_("Writing out segment IDs"));
+    for (row = 0; row < globals->nrows; row++) {
+
+	G_percent(row, globals->nrows, 9);
+
+	Rast_set_c_null_value(outbuf, 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;
+		}
+	    }
+	}
+	Rast_put_row(out_fd, outbuf, CELL_TYPE);
+    }
+
+    /* close and save segment id file */
+    Rast_close(out_fd);
+
+    /* set colors */
+    Rast_init_colors(&colors);
+    Rast_make_random_colors(&colors, 1, globals->nrows * globals->ncols);
+    Rast_write_colors(globals->out_name, G_mapset(), &colors);
+
+    Rast_short_history(globals->out_name, "raster", &hist);
+    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;
+
+	mean_fd = Rast_open_new(globals->out_band, FCELL_TYPE);
+	meanbuf = Rast_allocate_f_buf();
+
+	/* 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) */
+
+	thresh = globals->alpha * globals->alpha * globals->max_diff;
+	maxdev = globals->max_diff * (1 - globals->alpha * globals->alpha);
+	mingood = 1;
+
+	G_message(_("Writing out goodness of fit"));
+	for (row = 0; row < globals->nrows; row++) {
+
+	    G_percent(row, globals->nrows, 9);
+
+	    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) {
+
+			/* get values for Ri = this region */
+			globals->rs.id = rid;
+			fetch_reg_stats(Ri.row, Ri.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) {
+
+			    /* 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);
+			}
+			
+			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->max_diff;
+			    meanbuf[col] = sim;
+			    if (mingood > sim)
+				mingood = sim;
+			}
+		    }
+		}
+	    }
+	    Rast_put_row(mean_fd, meanbuf, FCELL_TYPE);
+	}
+
+	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_short_history(globals->out_band, "raster", &hist);
+	Rast_command_history(&hist);
+	Rast_write_history(globals->out_band, &hist);
+
+	G_free(meanbuf);
+    }
+
+    /* free memory */
+    G_free(outbuf);
+    Rast_free_colors(&colors);
+
+    return TRUE;
+}
+
+int close_files(struct globals *globals)
+{
+    /* close segmentation files and output raster */
+    G_debug(1, "closing files");
+    segment_close(&globals->bands_seg);
+    if (globals->bounds_map)
+	segment_close(&globals->bounds_seg);
+
+    G_free(globals->bands_val);
+    G_free(globals->second_val);
+
+    segment_close(&globals->rid_seg);
+
+    flag_destroy(globals->null_flag);
+    flag_destroy(globals->candidate_flag);
+
+    /* anything else left to clean up? */
+
+    return TRUE;
+}