i.segment: +mean shift

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

imagery/i.segment/cluster.c

@@ -0,0 +1,459 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.clump
+ *
+ * AUTHOR(S):    Michael Shapiro - CERL
+ *               Markus Metz
+ *
+ * PURPOSE:      Recategorizes data in a raster map layer by grouping cells
+ *               that form physically discrete areas into unique categories.
+ *
+ * COPYRIGHT:    (C) 2006-2016 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.
+ *
+ ***************************************************************************/
+
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <math.h>
+#include <time.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include "iseg.h"
+
+#define INCR 1024
+
+/* defeats the purpose of padding ... */
+#define CISNULL(r, c) (((c) == 0 || (c) == ncols + 1 ? 1 : \
+                       (FLAG_GET(globals->null_flag, (r), (c - 1)))))
+
+
+CELL cluster_bands(struct globals *globals)
+{
+    register int col;
+    register int i, n;
+    /* input */
+    DCELL **prev_in, **cur_in, **temp_in;
+    struct ngbr_stats Ri, Rk, Rn;
+    int nin;
+    int diag;
+    int bcol;
+    /* output */
+    CELL OLD, NEW;
+    CELL *temp_clump, out_cell;
+    CELL *prev_clump, *cur_clump;
+    CELL *index, *renumber;
+    CELL label, cellmax;
+    int nrows, ncols;
+    int row;
+    int len;
+    int nalloc;
+    char *cname;
+    int cfd, csize;
+    CELL cat;
+    int mwrow, mwrow1, mwrow2, mwnrows, mwcol, mwcol1, mwcol2, mwncols, radiusc;
+    double diff, diff2, avgdiff, ka2, hspat, hspat2;
+    double hspec, hspecad, hspec2, hspec2ad;
+    LARGEINT count;
+
+    G_message(_("%d-band clustering with threshold %g"), globals->nbands, globals->hr);
+
+    nrows = Rast_window_rows();
+    ncols = Rast_window_cols();
+
+    hspec = globals->hr;
+    hspec2 = globals->hr * globals->hr;
+    nin = globals->nbands;
+    diag = (globals->nn == 8);
+    radiusc = globals->hs;
+    mwnrows = mwncols = radiusc * 2 + 1;
+
+    /* spatial bandwidth */
+    hspat = globals->hs;
+    if (hspat < 1)
+	hspat = 1.5;
+    hspat2 = hspat * hspat;
+
+    cellmax = ((CELL)1 << (sizeof(CELL) * 8 - 2)) - 1;
+    cellmax += ((CELL)1 << (sizeof(CELL) * 8 - 2));
+
+    Ri.mean = NULL;
+    Rk.mean = NULL;
+    Rn.mean = G_malloc(globals->datasize);
+
+    /* allocate clump index */
+    /* TODO: support smallest label ID > 1 */
+    nalloc = INCR;
+    index = (CELL *) G_malloc(nalloc * sizeof(CELL));
+    index[0] = 0;
+    renumber = NULL;
+
+    /* allocate DCELL buffers two columns larger than current window */
+    prev_in = (DCELL **) G_malloc(sizeof(DCELL *) * (ncols + 2));
+    cur_in = (DCELL **) G_malloc(sizeof(DCELL *) * (ncols + 2));
+    
+    prev_in[0] = (DCELL *) G_malloc(globals->datasize * (ncols + 2) * nin);
+    cur_in[0] = (DCELL *) G_malloc(globals->datasize * (ncols + 2) * nin);
+
+    Rast_set_d_null_value(cur_in[0], (ncols + 2) * nin);
+    Rast_set_d_null_value(prev_in[0], (ncols + 2) * nin);
+
+    for (i = 1; i < ncols + 2; i++) {
+	prev_in[i] = prev_in[i - 1] + nin; 
+	cur_in[i] = cur_in[i - 1] + nin; 
+    }
+
+    /* allocate CELL buffers two columns larger than current window */
+    len = (ncols + 2) * sizeof(CELL);
+    prev_clump = (CELL *) G_malloc(len);
+    cur_clump = (CELL *) G_malloc(len);
+
+    /* temp file for initial clump IDs */
+    cname = G_tempfile();
+    if ((cfd = open(cname, O_RDWR | O_CREAT | O_EXCL, 0600)) < 0)
+	G_fatal_error(_("Unable to open temp file"));
+    csize = ncols * sizeof(CELL);
+
+    /* initialize clump labels */
+    G_zero(cur_clump, len);
+    G_zero(prev_clump, len);
+    /* TODO: support smallest label ID > 1 */
+    label = 0;
+
+    /****************************************************
+     *                      PASS 1                      *
+     * pass thru the input, create initial clump labels *
+     ****************************************************/
+
+    G_message(_("Assigning initial region IDs..."));
+    for (row = 0; row < nrows; row++) {
+	G_percent(row, nrows, 2);
+
+	for (col = 1; col <= ncols; col++) {
+
+	    /* get band values */
+	    Segment_get(globals->bands_out, (void *)cur_in[col], row, col - 1);
+	    Ri.mean = cur_in[col];
+
+	    if (CISNULL(row, col)) {
+		cur_clump[col] = 0;
+		continue;
+	    }
+	    
+	    hspec2ad = hspec2;
+
+	    if (globals->ms_adaptive) {
+		/* adapt initial range bandwith */
+
+		mwrow1 = row - radiusc;
+		mwrow2 = mwrow1 + mwnrows;
+		if (mwrow1 < 0)
+		    mwrow1 = 0;
+		if (mwrow2 > nrows)
+		    mwrow2 = nrows;
+
+		mwcol1 = col - radiusc;
+		mwcol2 = mwcol1 + mwncols;
+		if (mwcol1 < 0)
+		    mwcol1 = 0;
+		if (mwcol2 > ncols)
+		    mwcol2 = ncols;
+
+		ka2 = hspec2; 	/* OTB: conductance parameter */
+		
+		avgdiff = 0;
+		count = 0;
+		for (mwrow = mwrow1; mwrow < mwrow2; mwrow++) {
+		    for (mwcol = mwcol1; mwcol < mwcol2; mwcol++) {
+			if ((FLAG_GET(globals->null_flag, mwrow, mwcol)))
+			    continue;
+			if (mwrow == row && mwcol == col)
+			    continue;
+
+			diff = mwrow - row;
+			diff2 = diff * diff;
+			diff = mwcol - col;
+			diff2 += diff * diff;
+
+			if (diff2 <= hspat2) {
+
+			    Segment_get(globals->bands_out, (void *)Rn.mean,
+					mwrow, mwcol);
+
+			    /* get spectral distance */
+			    diff2 = (globals->calculate_similarity)(&Ri, &Rn, globals);
+
+			    avgdiff += sqrt(diff2);
+			    count++;
+			}
+		    }
+		}
+		hspec2ad = 0;
+		if (avgdiff > 0) {
+		    avgdiff /= count;
+		    hspecad = hspec;
+		    /* OTB-like, contrast enhancing */
+		    hspecad = exp(-avgdiff * avgdiff / (2 * ka2)) * avgdiff;
+		    /* preference for large regions, from Perona Malik 1990 
+		     * if the settings are right, it could be used to reduce noise */
+		    /* hspecad = 1 / (1 + (avgdiff * avgdiff / (2 * hspec2))); */
+		    hspec2ad = hspecad * hspecad;
+		    G_debug(1, "avg spectral diff: %g", avgdiff);
+		    G_debug(1, "initial hspec2: %g", hspec2);
+		    G_debug(1, "adapted hspec2: %g", hspec2ad);
+		}
+	    }
+
+	    /*
+	     * if the cell values are different to the left and above
+	     * (diagonal: and above left and above right)
+	     * then we must start a new clump
+	     *
+	     * this new clump may eventually collide with another
+	     * clump and will have to be merged
+	     */
+
+	    /* try to connect the current cell to an existing clump */
+	    OLD = NEW = 0;
+
+	    Rk.mean = cur_in[col - 1];
+
+	    /* same clump as to the left */
+	    if (!CISNULL(row, col - 1) &&
+	        globals->calculate_similarity(&Ri, &Rk, globals) <= hspec2ad) {
+		OLD = cur_clump[col] = cur_clump[col - 1];
+	    }
+
+	    if (diag) {
+		/* check above right, center, left, in that order */
+		temp_clump = prev_clump + col + 1;
+		bcol = col + 1;
+		do {
+		    Rk.mean = prev_in[bcol];
+		    if (row > 0 && !CISNULL(row - 1, bcol) &&
+		        globals->calculate_similarity(&Ri, &Rk, globals) <= hspec2ad) {
+			cur_clump[col] = *temp_clump;
+			if (OLD == 0) {
+			    OLD = *temp_clump;
+			}
+			else {
+			    NEW = *temp_clump;
+
+			    /* threshold > 0 and diagonal requires a bit of extra work
+			     * because of bridge cells:
+			     * A similar to B, B similar to C, but A not similar to C
+			     * -> B is bridge cell */
+			    if (NEW != OLD) {
+				CELL *temp_clump2;
+
+				/* conflict! preserve NEW clump ID and change OLD clump ID.
+				 * Must go back to the left in the current row and to the right
+				 * in the previous row to change all the clump values as well.
+				 */
+
+				/* left of the current row from 1 to col - 1 */
+				temp_clump2 = cur_clump;
+				n = col - 1;
+				while (n-- > 0) {
+				    temp_clump2++;	/* skip left edge */
+				    if (*temp_clump2 == OLD)
+					*temp_clump2 = NEW;
+				}
+
+				/* right of previous row from col - 1 to ncols */
+				temp_clump2 = prev_clump + col - 1;
+				n = ncols - col + 2;
+				while (n-- > 0) {
+				    if (*temp_clump2 == OLD)
+					*temp_clump2 = NEW;
+				    temp_clump2++;
+				}
+
+				/* modify the OLD index */
+				index[OLD] = NEW;
+
+				OLD = NEW;
+				NEW = 0;
+			    }
+			}
+		    }
+		    temp_clump--;
+		} while (bcol-- > col - 1);
+	    }
+	    else {
+		/* check above */
+		Rk.mean = prev_in[col];
+
+		if (row > 0 && !CISNULL(row - 1, col) &&
+		    globals->calculate_similarity(&Ri, &Rk, globals) <= hspec2ad) {
+		    temp_clump = prev_clump + col;
+		    cur_clump[col] = *temp_clump;
+		    if (OLD == 0) {
+			OLD = *temp_clump;
+			}
+		    else {
+			NEW = *temp_clump;
+			if (NEW != OLD) {
+
+			    /* conflict! preserve NEW clump ID and change OLD clump ID.
+			     * Must go back to the left in the current row and to the right
+			     * in the previous row to change all the clump values as well.
+			     */
+
+			    /* left of the current row from 1 to col - 1 */
+			    temp_clump = cur_clump;
+			    n = col - 1;
+			    while (n-- > 0) {
+				temp_clump++;	/* skip left edge */
+				if (*temp_clump == OLD)
+				    *temp_clump = NEW;
+			    }
+
+			    /* right of previous row from col + 1 to ncols */
+			    temp_clump = prev_clump + col;
+			    n = ncols - col;
+			    while (n-- > 0) {
+				temp_clump++;	/* skip col */
+				if (*temp_clump == OLD)
+				    *temp_clump = NEW;
+			    }
+
+			    /* modify the OLD index */
+			    index[OLD] = NEW;
+
+			    OLD = NEW;
+			    NEW = 0;
+			}
+		    }
+		}
+	    }
+
+	    if (NEW == 0 || OLD == NEW) {	/* ok */
+		if (OLD == 0) {
+		    /* start a new clump */
+		    if (label == cellmax)
+			G_fatal_error(_("Too many objects: integer overflow"));
+
+		    label++;
+		    cur_clump[col] = label;
+		    if (label >= nalloc) {
+			nalloc += INCR;
+			index =
+			    (CELL *) G_realloc(index,
+					       nalloc * sizeof(CELL));
+		    }
+		    index[label] = label;
+		}
+	    }
+	    /* else the relabelling above failed */
+	}
+
+	/* write initial clump IDs */
+	/* this works also with writing out cur_clump, but only 
+	 * prev_clump is complete and will not change any more */
+	if (row > 0) {
+	    if (write(cfd, prev_clump + 1, csize) != csize)
+		G_fatal_error(_("Unable to write to temp file"));
+	}
+
+	/* switch the buffers so that the current buffer becomes the previous */
+	temp_in = cur_in;
+	cur_in = prev_in;
+	prev_in = temp_in;
+
+	temp_clump = cur_clump;
+	cur_clump = prev_clump;
+	prev_clump = temp_clump;
+    }
+    /* write last row with initial clump IDs */
+    if (write(cfd, prev_clump + 1, csize) != csize)
+	G_fatal_error(_("Unable to write to temp file"));
+    G_percent(1, 1, 1);
+
+    G_free(prev_in[0]);
+    G_free(cur_in[0]);
+    G_free(prev_in);
+    G_free(cur_in);
+
+    /* generate a renumbering scheme */
+    G_message(_("Generating renumbering scheme..."));
+    G_debug(1, "%d initial labels", label);
+    /* allocate final clump ID */
+    renumber = (CELL *) G_malloc((label + 1) * sizeof(CELL));
+    renumber[0] = 0;
+    cat = 0;
+    G_percent(0, label, 1);
+    for (n = 1; n <= label; n++) {
+	G_percent(n, label, 1);
+	OLD = n;
+	NEW = index[n];
+	if (OLD != NEW) {
+	    renumber[n] = 0;
+	    /* find valid clump ID */
+	    while (OLD != NEW) {
+		OLD = NEW;
+		NEW = index[OLD];
+	    }
+	    index[n] = NEW;
+	}
+	else
+	    /* set final clump id */
+	    renumber[n] = ++cat;
+    }
+
+    if (cat > cellmax - globals->max_rid)
+	G_fatal_error(_("Too many objects: integer overflow"));
+
+    /* rewind temp file */
+    lseek(cfd, 0, SEEK_SET);
+
+    /****************************************************
+     *                      PASS 2                      *
+     * apply renumbering scheme to initial clump labels *
+     ****************************************************/
+
+    /* the input raster is no longer needed, 
+     * using instead the temp file with initial clump labels */
+
+    G_message(_("Assigning final region IDs..."));
+    for (row = 0; row < nrows; row++) {
+
+	G_percent(row, nrows, 2);
+    
+	if (read(cfd, cur_clump, csize) != csize)
+	    G_fatal_error(_("Unable to read from temp file"));
+
+	temp_clump = cur_clump;
+
+	for (col = 0; col < ncols; col++) {
+	    if (!(FLAG_GET(globals->null_flag, row, col))) {
+		out_cell = renumber[index[*temp_clump]] + globals->max_rid;
+
+		Segment_put(&globals->rid_seg, (void *)&out_cell, row, col);
+	    }
+	    temp_clump++;
+	}
+    }
+    G_percent(1, 1, 1);
+
+    close(cfd);
+    unlink(cname);
+
+    /* free */
+    G_free(prev_clump);
+    G_free(cur_clump);
+    G_free(index);
+    G_free(renumber);
+
+    G_message(_("Found %d clumps"), cat);
+    globals->max_rid += cat;
+
+    return globals->max_rid;
+}

imagery/i.segment/i.segment.html

@@ -1,6 +1,6 @@
 <h2>DESCRIPTION</h2>
 Image segmentation or object recognition is the process of grouping 
-similar pixels into unique segments, also referred to as objects. 
+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 
@@ -28,7 +28,7 @@ distance will be smaller than to any other neighbor), and</li>
 process is repeated until no merges are made during a complete pass.</li>
 </ol>
 
-<h3>Similarity and Threshold</h3>
+<h4>Similarity and Threshold</h4>
 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.
@@ -63,7 +63,7 @@ and <tt>smoothness = Perimeter Length / Bounding Box</tt>. The
 perimeter length is estimated as the number of pixel sides the segment
 has.
 
-<h3>Seeds</h3>
+<h4>Seeds</h4>
 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. If the seeds are the results of a previous segmentation 
@@ -71,25 +71,41 @@ with lower threshold, hierarchical segmentation can be performed. 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.
+
+<h4>Maximum number of segments</h4>
+
+The current limit with CELL storage used for segment IDs is 2
+billion starting segment IDs. Segment IDs are assigned whenever a yet 
+unprocessed pixel is merged with another segment. Integer overflow can 
+happen for computational regions with more than 2 billion cells and 
+very low threshold values, resulting in many segments. If integer 
+overflow occurs durin region growing, starting segments can be used 
+(created by initial classification or other methods).
+
+<h4>Goodness of Fit</h4>
+The <b>goodness</b> of fit for each pixel is calculated as 1 - distance 
+of the pixel to the object it belongs to. The distance is calculated 
+with the selected <b>similarity</b> method. A value of 1 means 
+identical values, perfect fit, and a value of 0 means maximum possible 
+distance, worst possible fit.
+
+<h3>Mean shift</h3>
+Mean shift image segmentation consists of 2 steps: anisotrophic 
+filtering and 2. clustering. For anisotrophic filtering new cell values 
+are calculated from all pixels not farther than <b>hs</b> pixels away 
+from the current pixel and with a spectral difference not larger than 
+<b>hr</b>. That means that pixels that are too different from the current 
+pixel are not considered in the calculation of new pixel values. 
+<b>hs</b> and <b>hr</b> are the spatial and spectral (range) bandwidths 
+for anisotrophic filtering. Cell values are iteratively recalculated 
+(shifted to the segment's mean) until the maximum number of iterations 
+is reached or until the largest shift is smaller than <b>threshold</b>.
+
 <p>
-It is expected that the <b>minsize</b> 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 will also 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:
-<ol>
-  <li>Use starting seed pixels. (Maximum 2 billion pixels can be 
-labeled with positive integers.)</li>
-  <li>Use starting seed segments. (By initial classification or other 
-methods.)</li>
-</ol>
+If input bands have been reprojected, they should not be reprojected 
+with bilinear resampling because that method causes smooth transitions 
+between objects. More appropriate methods are bicubic or lanczos 
+resampling.
 
 <h3>Boundary Constraints</h3>
 Boundary constraints limit the adjacency of pixels and segments. 
@@ -98,18 +114,13 @@ 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 <b>minsize</b> greater 
+To reduce the salt and pepper effect, a <b>minsize</b> 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.
-
-<h3>Goodness of Fit</h3>
-The <b>goodness</b> of fit for each pixel is calculated as 1 - distance 
-of the pixel to the object it belongs to. The distance is calculated 
-with the selected <b>similarity</b> method. A value of 1 means 
-identical values, perfect fit, and a value of 0 means maximum possible 
-distance, worst possible fit.
+most similar neighbor. A minimum segment size larger than 1 is 
+recommended when using adaptive bandwith selected with the <b>-a</b> 
+flag.
 
 <h2>EXAMPLES</h2>
 

imagery/i.segment/iseg.h

@@ -67,9 +67,8 @@ struct globals
     /* 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 */
+    char *gof;			/* indicator for segment heterogeneity / goodness of fit */
+    char *bsuf;			/* suffix to be appended to input bands */
 
     /* general segmentation */
     int method;			/* Segmentation method code */
@@ -88,6 +87,8 @@ struct globals
 
     /* mean shift */
     double hs, hr;		/* spectral and range bandwidth */
+    int ms_adaptive;		/* use adaptive bandwidth */
+    int ms_progressive;		/* use progressive bandwidth */
 
     /* region info */
     int nrows, ncols;
@@ -159,6 +160,9 @@ int region_growing(struct globals *);
 /* mean_shift.c */
 int mean_shift(struct globals *);
 
+/* cluster.c */
+CELL cluster_bands(struct globals *globals);
+
 /* watershed.c */
 int watershed(struct globals *);
 

imagery/i.segment/main.c

@@ -54,12 +54,12 @@ int main(int argc, char *argv[])
     if (write_ids(&globals) != TRUE)
 	G_fatal_error(_("Error in writing IDs"));
 
-    if (globals.method == ORM_RG && globals.out_band) {
+    if (globals.method == ORM_RG && globals.gof) {
 	if (write_gof_rg(&globals) != TRUE)
 	    G_fatal_error(_("Error in writing goodness of fit"));
     }
 
-    if (globals.method == ORM_MS) {
+    if (globals.method == ORM_MS && globals.bsuf) {
 	if (write_bands_ms(&globals) != TRUE)
 	    G_fatal_error(_("Error in writing new band values"));
     }

imagery/i.segment/mean_shift.c

@@ -13,6 +13,8 @@
 #include <grass/rbtree.h>	/* Red Black Tree library functions */
 #include "iseg.h"
 
+int remove_small_clumps(struct globals *globals);
+
 /* standard gauss funtion:
  * a * exp(-(x - m)^2 / (2 * stddev^2)
  * a is not needed because the sum of weights is calculated for each 
@@ -33,82 +35,740 @@ static double gauss_kernel(double diff2, double var)
 
 int mean_shift(struct globals *globals)
 {
-    int row, col, t;
-    int n_changes;
-    double alpha2;
-    struct ngbr_stats Rin, Rout;
-    double diff2;
-
-    G_fatal_error(_("Mean shift is not yet implemented"));
-    return FALSE;
-
+    int row, col, t, n;
+    int mwrow, mwrow1, mwrow2, mwnrows, mwcol, mwcol1, mwcol2, mwncols, radiusc;
+    double hspat, hspec, hspat2, hspec2, sigmaspat2, sigmaspec2;
+    double hspecad, hspecad2;
+    double ka2;
+    double w, wsum;
+    LARGEINT n_changes;
+    double alpha2, maxdiff2;
+    struct ngbr_stats Rin, Rout, Rn;
+    double diff, diff2;
+    SEGMENT *seg_tmp;
+    double mindiff, mindiffzero, mindiffavg, mindiffzeroavg;
+    double avgdiff, avgdiffavg;
+    LARGEINT nvalid, count;
+    int do_gauss = 0, do_adaptive, do_progressive;
 
     Rin.mean = G_malloc(globals->datasize);
     Rout.mean = G_malloc(globals->datasize);
-
-    /* TODO: need another segment structure holding output
-     * for mean shift, output of the previous iteration becomes input of the next iteration
-     * initially, input and output are original band values */
+    Rn.mean = G_malloc(globals->datasize);
 
     alpha2 = globals->alpha * globals->alpha;
+    do_adaptive = globals->ms_adaptive;
+    do_progressive = globals->ms_progressive;
     
+    globals->candidate_count = 0;
+    flag_clear_all(globals->candidate_flag);
+
+    /* Set candidate flag to true/1 for all non-NULL cells */
+    for (row = globals->row_min; row < globals->row_max; row++) {
+	for (col = globals->col_min; col < globals->col_max; col++) {
+	    if (!(FLAG_GET(globals->null_flag, row, col))) {
+		
+		FLAG_SET(globals->candidate_flag, row, col);
+		globals->candidate_count++;
+	    }
+	}
+    }
+
+    /* spatial bandwidth */
+    hspat = globals->hs;
+    if (hspat < 1) {
+	hspat = 1.5;
+	globals->hs = hspat;
+    }
+
+    hspat2 = hspat * hspat;
+    sigmaspat2 = hspat2 / 9.;
+    radiusc = hspat;	/* radius in cells truncated to integer */
+    mwnrows = mwncols = radiusc * 2 + 1;
+
+    /* estimate spectral bandwidth for given spatial bandwidth */
+    mindiffavg = mindiffzeroavg = 0;
+    avgdiffavg = 0;
+    nvalid = 0;
+
+    G_message(_("Estimating spectral bandwidth for spatial bandwidth %g..."), hspat);
+    G_percent_reset();
+    for (row = globals->row_min; row < globals->row_max; row++) {
+	G_percent(row - globals->row_min,
+		  globals->row_max - globals->row_min, 4);
+
+	mwrow1 = row - radiusc;
+	mwrow2 = mwrow1 + mwnrows;
+	if (mwrow1 < globals->row_min)
+	    mwrow1 = globals->row_min;
+	if (mwrow2 > globals->row_max)
+	    mwrow2 = globals->row_max;
+
+	for (col = globals->col_min; col < globals->col_max; col++) {
+	    if ((FLAG_GET(globals->null_flag, row, col)))
+		continue;
+
+	    /* get current band values */
+	    Segment_get(globals->bands_in, (void *)Rin.mean,
+			row, col);
+
+	    mwcol1 = col - radiusc;
+	    mwcol2 = mwcol1 + mwncols;
+	    if (mwcol1 < globals->col_min)
+		mwcol1 = globals->col_min;
+	    if (mwcol2 > globals->col_max)
+		mwcol2 = globals->col_max;
+	    
+	    /* get minimum spectral distance for this cell */
+	    count = 0;
+	    mindiff = globals->max_diff;
+	    mindiffzero = globals->max_diff;
+	    avgdiff = 0;
+	    for (mwrow = mwrow1; mwrow < mwrow2; mwrow++) {
+		for (mwcol = mwcol1; mwcol < mwcol2; mwcol++) {
+		    if ((FLAG_GET(globals->null_flag, mwrow, mwcol)))
+			continue;
+		    if (mwrow == row && mwcol == col)
+			continue;
+
+		    diff = mwrow - row;
+		    diff2 = diff * diff;
+		    diff = mwcol - col;
+		    diff2 += diff * diff;
+
+		    if (diff2 <= hspat2) {
+
+			Segment_get(globals->bands_in, (void *)Rn.mean,
+				    mwrow, mwcol);
+
+			/* get spectral distance */
+			diff2 = (globals->calculate_similarity)(&Rin, &Rn, globals);
+
+			if (mindiff > diff2)
+			    mindiff = diff2;
+			if (mindiffzero > diff2 && diff2 > 0)
+			    mindiffzero = diff2;
+			avgdiff += sqrt(diff2);
+			count++;
+		    }
+		}
+	    }
+	    if (count) {
+		nvalid++;
+		if (mindiff > 0)
+		    mindiffavg += sqrt(mindiff);
+		mindiffzeroavg += sqrt(mindiffzero);
+		if (avgdiff > 0)
+		    avgdiffavg += avgdiff / count;
+	    }
+	}
+    }
+    G_percent(1, 1, 1);
+    if (!nvalid) {
+	G_fatal_error(_("Empty moving windows"));
+    }
+
+    mindiffavg /= nvalid;
+    mindiffzeroavg /= nvalid;
+    avgdiffavg /= nvalid;
+    G_debug(1, "Average minimum difference to neighbours: %g", mindiffavg);
+    G_debug(1, "Average minimum difference excl zero to neighbours: %g", mindiffzeroavg);
+    G_debug(1, "Average average difference to neighbours: %g", avgdiffavg);
+
+    /* use avgdiffavg as hspec for adaptive bandwidth */
+
+    hspec = globals->hr;
+    if (hspec < 0 || hspec >= 1) {
+	hspec = sqrt(avgdiffavg / 10.0);
+	hspec = avgdiffavg;
+	hspec = mindiffzeroavg;
+	
+	if (do_progressive)
+	    G_message(_("Initial range bandwidth: %g"), hspec);
+	else
+	    G_message(_("Estimated range bandwidth: %g"), hspec);
+	globals->hr = hspec;
+    }
+    else {
+	G_message(_("Estimated range bandwidth: %g"), mindiffzeroavg);
+    }
+    if (do_adaptive) {
+	/* bandwidth is now standard deviation for adaptive bandwidth 
+	 * using a gaussian function with range bandwith used as 
+	 * bandwidth for the gaussian function
+	 * the aim is to produce similar but improved results with 
+	 * adaptive bandwidth
+	 * thus increase bandwidth */
+	hspec = sqrt(hspec);
+    }
+
+    hspec2 = hspec * hspec;
+    sigmaspec2 = hspec2 / 9.;
+
+    if (!do_progressive) {
+	G_message(_("Spatial bandwidth: %g"), hspat);
+	G_message(_("Range bandwidth: %g"), hspec);
+    }
+
+    G_debug(4, "Starting to process %ld candidate cells",
+	    globals->candidate_count);
+
     t = 0;
     n_changes = 1;
+    maxdiff2 = 0;
     while (t < globals->end_t && n_changes > 0) {
 
 	G_message(_("Processing pass %d..."), ++t);
 
-	n_changes = 0;
-	globals->candidate_count = 0;
-	flag_clear_all(globals->candidate_flag);
+	/* cells within an object should become more similar with each pass
+	 * therefore the spectral bandwidth could be decreased
+	 * and the spatial bandwidth could be increased */
 
-	/* Set candidate flag to true/1 for all non-NULL cells */
-	for (row = globals->row_min; row < globals->row_max; row++) {
-	    for (col = globals->col_min; col < globals->col_max; col++) {
-		if (!(FLAG_GET(globals->null_flag, row, col))) {
-		    
-		    FLAG_SET(globals->candidate_flag, row, col);
-		    globals->candidate_count++;
-		}
-	    }
+	/* spatial bandwidth: double the area covered by the moving window
+	 * area = M_PI * hspat * hspat
+	 * new hspat = sqrt(M_PI * hspat * hspat * 2 / M_PI)
+	 * no good, too large increases */
+
+	if (do_progressive) {
+	    if (t > 1)
+		hspat *= 1.1;
+	    hspat2 = hspat * hspat;
+	    sigmaspat2 = hspat2 / 9.;
+	    radiusc = hspat;	/* radius in cells truncated to integer */
+	    mwnrows = mwncols = radiusc * 2 + 1;
+
+	    /* spectral bandwidth: reduce by 0.7 */
+	    if (t > 1)
+		hspec *= 0.9;
+	    hspec2 = hspec * hspec;
+	    sigmaspec2 = hspec2 / 9.;
+
+	    G_verbose_message(_("Spatial bandwidth: %g"), hspat);
+	    G_verbose_message(_("Range bandwidth: %g"), hspec);
 	}
 
-	G_debug(4, "Starting to process %ld candidate cells",
-		globals->candidate_count);
+	n_changes = 0;
+	maxdiff2 = 0;
+
+	/* swap input and output */
+	seg_tmp = globals->bands_in;
+	globals->bands_in = globals->bands_out;
+	globals->bands_out = seg_tmp;
 
 	/*process candidate cells */
 	G_percent_reset();
 	for (row = globals->row_min; row < globals->row_max; row++) {
 	    G_percent(row - globals->row_min,
 	              globals->row_max - globals->row_min, 4);
+
+	    mwrow1 = row - radiusc;
+	    mwrow2 = mwrow1 + mwnrows;
+	    if (mwrow1 < globals->row_min)
+		mwrow1 = globals->row_min;
+	    if (mwrow2 > globals->row_max)
+		mwrow2 = globals->row_max;
+
 	    for (col = globals->col_min; col < globals->col_max; col++) {
-		if (!(FLAG_GET(globals->candidate_flag, row, col)))
+		if ((FLAG_GET(globals->null_flag, row, col)))
 		    continue;
 
 		/* get current band values */
-		Segment_get(&globals->bands_seg, (void *)Rin.mean,
+		Segment_get(globals->bands_in, (void *)Rin.mean,
 			    row, col);
+		
+		/* init output */
+		for (n = 0; n < globals->nbands; n++)
+		    Rout.mean[n] = 0.;
 
-		/* adapt initial spatial and range bandwiths */
+		mwcol1 = col - radiusc;
+		mwcol2 = mwcol1 + mwncols;
+		if (mwcol1 < globals->col_min)
+		    mwcol1 = globals->col_min;
+		if (mwcol2 > globals->col_max)
+		    mwcol2 = globals->col_max;
 
-		/* calculate new band values */
+		hspecad2 = hspec2;
 		
+		if (do_adaptive) {
+		    /* adapt initial range bandwith */
+		    
+		    ka2 = hspec2; 	/* OTB: conductance parameter */
+		    
+		    avgdiff = 0;
+		    count = 0;
+		    for (mwrow = mwrow1; mwrow < mwrow2; mwrow++) {
+			for (mwcol = mwcol1; mwcol < mwcol2; mwcol++) {
+			    if ((FLAG_GET(globals->null_flag, mwrow, mwcol)))
+				continue;
+			    if (mwrow == row && mwcol == col)
+				continue;
+
+			    diff = mwrow - row;
+			    diff2 = diff * diff;
+			    diff = mwcol - col;
+			    diff2 += diff * diff;
+
+			    if (diff2 <= hspat2) {
+
+				Segment_get(globals->bands_in, (void *)Rn.mean,
+					    mwrow, mwcol);
+
+				/* get spectral distance */
+				diff2 = (globals->calculate_similarity)(&Rin, &Rn, globals);
+
+				avgdiff += sqrt(diff2);
+				count++;
+			    }
+			}
+		    }
+		    hspecad2 = 0;
+		    if (avgdiff > 0) {
+			avgdiff /= count;
+			hspecad = hspec;
+			/* OTB-like, contrast enhancing */
+			hspecad = exp(-avgdiff * avgdiff / (2 * ka2)) * avgdiff;
+			/* preference for large regions, from Perona Malik 1990 
+			 * if the settings are right, it could be used to reduce noise */
+			/* hspecad = 1 / (1 + (avgdiff * avgdiff / (2 * hspec2))); */
+			hspecad2 = hspecad * hspecad;
+			G_debug(1, "avg spectral diff: %g", avgdiff);
+			G_debug(1, "initial hspec2: %g", hspec2);
+			G_debug(1, "adapted hspec2: %g", hspecad2);
+		    }
+		}
+		
+		/* actual mean shift */
+		wsum = 0;
+		for (mwrow = mwrow1; mwrow < mwrow2; mwrow++) {
+		    for (mwcol = mwcol1; mwcol < mwcol2; mwcol++) {
+			if ((FLAG_GET(globals->null_flag, mwrow, mwcol)))
+			    continue;
+			diff = mwrow - row;
+			diff2 = diff * diff;
+			diff = mwcol - col;
+			diff2 += diff * diff;
+
+			if (diff2 <= hspat2) {
+			    w = 1;
+			    if (do_gauss)
+				w = gauss_kernel(diff2, sigmaspat2);
+
+			    Segment_get(globals->bands_in, (void *)Rn.mean,
+					mwrow, mwcol);
+
+			    /* check spectral distance */
+			    diff2 = (globals->calculate_similarity)(&Rin, &Rn, globals);
+			    if (diff2 <= hspecad2) {
+				if (do_gauss)
+				    w *= gauss_kernel(diff2, sigmaspec2);
+				wsum += w;
+				for (n = 0; n < globals->nbands; n++)
+				    Rout.mean[n] += w * Rn.mean[n];
+			    }
+			}
+		    }
+		}
+		
+		if (wsum > 0) {
+		    for (n = 0; n < globals->nbands; n++)
+			Rout.mean[n] /= wsum;
+		}
+		else {
+		    for (n = 0; n < globals->nbands; n++)
+			Rout.mean[n] = Rin.mean[n];
+		}
+
+		/* put new band values */
+		Segment_put(globals->bands_out, (void *)Rout.mean,
+			    row, col);
+
 		/* if the squared difference between old and new band values 
 		 * is larger than alpha2, then increase n_changes */
 		
 		diff2 = (globals->calculate_similarity)(&Rin, &Rout, globals);
 		if (diff2 > alpha2)
 		    n_changes++;
+		if (maxdiff2 < diff2)
+		    maxdiff2 = diff2;
 	    }
 	}
+	G_percent(1, 1, 1);
+	G_message(_("Changes > threshold: %d, largest change: %g"), n_changes, sqrt(maxdiff2));
     }
     if (n_changes > 1)
 	G_message(_("Mean shift stopped at %d due to reaching max iteration limit, more changes may be possible"), t);
     else
 	G_message(_("Mean shift converged after %d iterations"), t);
-    
+
     /* identify connected components */
+    cluster_bands(globals);
+
+    /* remove small regions */
+    remove_small_clumps(globals);
     
     return TRUE;
 }
 
+static int cmp_rc(const void *first, const void *second)
+{
+    struct rc *a = (struct rc *)first, *b = (struct rc *)second;
+
+    if (a->row == b->row)
+	return (a->col - b->col);
+
+    return (a->row - b->row);
+}
+
+static int find_best_neighbour(struct globals *globals, int row, int col,
+                               int this_id, struct NB_TREE *nbtree,
+			       int *reg_size, struct ngbr_stats **Rbest,
+			       int *best_n_row, int *best_n_col)
+{
+    int rown, coln, n, count;
+    int neighbors[8][2];
+    struct rc next, ngbr_rc;
+    struct rclist rilist;
+    int no_check;
+    int ngbr_id;
+    struct RB_TREE *visited;
+    struct ngbr_stats Ri, Rk, *Rfound;
+    double sim, best_sim;
+    int best_n_id;
+    int have_Ri;
+
+    Ri.mean = G_malloc(globals->datasize);
+    Rk.mean = G_malloc(globals->datasize);
+
+    nbtree_clear(nbtree);
+    FLAG_UNSET(globals->candidate_flag, row, col);
+
+    visited = rbtree_create(cmp_rc, sizeof(struct rc));
+    ngbr_rc.row = row;
+    ngbr_rc.col = col;
+    rbtree_insert(visited, &ngbr_rc);
+
+    /* breadth-first search */
+    next.row = row;
+    next.col = col;
+    rclist_init(&rilist);
+    count = 1;
+    best_n_id = -1;
+    best_sim = 2;
+    
+    do {
+	have_Ri = 0;
+	globals->find_neighbors(next.row, next.col, neighbors);
+	n = globals->nn - 1;
+	do {
+	    rown = neighbors[n][0];
+	    coln = neighbors[n][1];
+	    no_check = 0;
+	    if (rown < globals->row_min || rown >= globals->row_max ||
+		coln < globals->col_min || coln >= globals->col_max)
+		no_check = 1;
+	    if (!no_check && (FLAG_GET(globals->null_flag, rown, coln)))
+		no_check = 1;
+
+	    ngbr_rc.row = rown;
+	    ngbr_rc.col = coln;
+
+	    if (!no_check && !rbtree_find(visited, &ngbr_rc)) {
+		rbtree_insert(visited, &ngbr_rc);
+
+		/* get neighbor ID */
+		Segment_get(&globals->rid_seg, (void *) &ngbr_id, rown, coln);
+		/* same neighbour */
+		if (ngbr_id == this_id) {
+		    count++;
+		    rclist_add(&rilist, rown, coln);
+		    FLAG_UNSET(globals->candidate_flag, rown, coln);
+		}
+		else { /* different neighbour */
+		    /* compare to this cell next.row, next.col */
+		    if (!have_Ri) {
+			Segment_get(globals->bands_out, (void *) Ri.mean, next.row, next.col);
+			have_Ri = 1;
+		    }
+		    Segment_get(globals->bands_out, (void *) Rk.mean, rown, coln);
+		    sim = globals->calculate_similarity(&Ri, &Rk, globals);
+		    if (best_sim > sim) {
+			best_sim = sim;
+			best_n_id = ngbr_id;
+			*best_n_row = rown;
+			*best_n_col = coln;
+		    }
+
+		    /* find in neighbor tree */
+		    Rk.id = ngbr_id;
+		    if ((Rfound = nbtree_find(nbtree, &Rk))) {
+			Rfound->count++;
+			if (*Rbest && (*Rbest)->count < Rfound->count)
+			    *Rbest = Rfound;
+		    }
+		    else {
+			Rk.count = 1;
+			Rk.row = rown;
+			Rk.col = coln;
+			nbtree_insert(nbtree, &Rk);
+			if (!(*Rbest))
+			    *Rbest = nbtree_find(nbtree, &Rk);
+		    }
+		}
+	    }
+	} while (n--);    /* end do loop - next neighbor */
+    } while (rclist_drop(&rilist, &next));   /* while there are cells to check */
+
+    rclist_destroy(&rilist);
+    rbtree_destroy(visited);
+    G_free(Ri.mean);
+    G_free(Rk.mean);
+    
+    *reg_size = count;
+
+    return best_n_id;
+}
+
+static int check_reg_size(struct globals *globals, int minsize, int row, int col)
+{
+    int rown, coln, n;
+    int neighbors[8][2];
+    int this_id;
+    int ngbr_id;
+    LARGEINT reg_size;
+    struct RB_TREE *visited;
+    struct rc next, ngbr_rc;
+    struct rclist rilist;
+    int no_check;
+    
+    if (!(FLAG_GET(globals->candidate_flag, row, col)))
+	return minsize;
+
+    FLAG_UNSET(globals->candidate_flag, row, col);
+
+    visited = rbtree_create(cmp_rc, sizeof(struct rc));
+    ngbr_rc.row = row;
+    ngbr_rc.col = col;
+    rbtree_insert(visited, &ngbr_rc);
+
+    /* get this ID */
+    Segment_get(&globals->rid_seg, (void *) &this_id, row, col);
+
+    /* breadth-first search */
+    next.row = row;
+    next.col = col;
+    rclist_init(&rilist);
+    reg_size = 1;
+
+    do {
+	globals->find_neighbors(next.row, next.col, neighbors);
+	n = globals->nn - 1;
+	do {
+	    rown = neighbors[n][0];
+	    coln = neighbors[n][1];
+	    no_check = 0;
+	    if (rown < globals->row_min || rown >= globals->row_max ||
+		coln < globals->col_min || coln >= globals->col_max)
+		no_check = 1;
+	    if (!no_check && (FLAG_GET(globals->null_flag, rown, coln)))
+		no_check = 1;
+
+	    ngbr_rc.row = rown;
+	    ngbr_rc.col = coln;
+
+	    if (!no_check && !rbtree_find(visited, &ngbr_rc)) {
+		rbtree_insert(visited, &ngbr_rc);
+
+		/* get neighbour ID */
+		Segment_get(&globals->rid_seg, (void *) &ngbr_id, rown, coln);
+		/* same neighbour */
+		if (ngbr_id == this_id) {
+		    reg_size++;
+		    rclist_add(&rilist, rown, coln);
+		    FLAG_UNSET(globals->candidate_flag, rown, coln);
+		}
+	    }
+	} while (n--);    /* end do loop - next neighbor */
+    } while (rclist_drop(&rilist, &next));   /* while there are cells to check */
+
+    rclist_destroy(&rilist);
+    rbtree_destroy(visited);
+
+    return reg_size;
+}
+
+static int update_rid(struct globals *globals, int row, int col, int new_id)
+{
+    int rown, coln, n;
+    int neighbors[8][2];
+    int this_id;
+    int ngbr_id;
+    struct rc next;
+    struct rclist rilist;
+    int no_check;
+
+    /* get this ID */
+    Segment_get(&globals->rid_seg, (void *) &this_id, row, col);
+    Segment_put(&globals->rid_seg, (void *) &new_id, row, col);
+
+    /* breadth-first search */
+    next.row = row;
+    next.col = col;
+    rclist_init(&rilist);
+
+    do {
+	globals->find_neighbors(next.row, next.col, neighbors);
+	n = globals->nn - 1;
+	do {
+	    rown = neighbors[n][0];
+	    coln = neighbors[n][1];
+	    no_check = 0;
+	    if (rown < globals->row_min || rown >= globals->row_max ||
+		coln < globals->col_min || coln >= globals->col_max)
+		no_check = 1;
+	    if (!no_check && (FLAG_GET(globals->null_flag, rown, coln)))
+		no_check = 1;
+
+	    if (!no_check) {
+		/* get neighbour ID */
+		Segment_get(&globals->rid_seg, (void *) &ngbr_id, rown, coln);
+		/* same neighbour */
+		if (ngbr_id == this_id) {
+		    rclist_add(&rilist, rown, coln);
+		    Segment_put(&globals->rid_seg, (void *) &new_id, rown, coln);
+		}
+	    }
+	} while (n--);    /* end do loop - next neighbor */
+    } while (rclist_drop(&rilist, &next));   /* while there are cells to check */
+
+    rclist_destroy(&rilist);
+
+    return 1;
+}
+
+int remove_small_clumps(struct globals *globals)
+{
+    int row, col, i;
+    struct NB_TREE *nbtree;
+    int this_id;
+    struct ngbr_stats *Rbest;
+    int best_n_id, best_n_row, best_n_col;
+    int reg_size;
+    CELL *renumber, n_regions, min_rid;
+
+    /* two possible modes
+     * best (most similar) neighbor
+     * neighbor with longest shared boundary 
+     */
+
+    if (globals->min_segment_size < 2)
+	return 0;
+
+    G_message(_("Merging segments smaller than %d cells..."), globals->min_segment_size);
+
+    /* init renumber */
+    renumber = (CELL *) G_malloc(sizeof(CELL) * (globals->max_rid + 1));
+    for (i = 0; i <= globals->max_rid; i++)
+	renumber[i] = 0;
+
+    /* clear candidate flag */
+    flag_clear_all(globals->candidate_flag);
+
+    min_rid = globals->max_rid;
+
+    /* Set candidate flag to true/1 for all non-NULL cells */
+    for (row = globals->row_min; row < globals->row_max; row++) {
+	for (col = globals->col_min; col < globals->col_max; col++) {
+	    if (!(FLAG_GET(globals->null_flag, row, col))) {
+		FLAG_SET(globals->candidate_flag, row, col);
+		Segment_get(&globals->rid_seg, (void *) &this_id, row, col);
+		/* renumber is region size */
+		if (renumber[this_id] <= globals->min_segment_size) {
+		    renumber[this_id]++;
+		    if (min_rid > this_id)
+			min_rid = this_id;
+		}
+	    }
+	}
+    }
+    min_rid--;
+
+    nbtree = nbtree_create(globals->nbands, globals->datasize);
+
+    /* go through all cells */
+    G_percent_reset();
+    for (row = globals->row_min; row < globals->row_max; row++) {
+	G_percent(row - globals->row_min,
+		  globals->row_max - globals->row_min, 2);
+	for (col = globals->col_min; col < globals->col_max; col++) {
+	    if ((FLAG_GET(globals->null_flag, row, col)))
+		continue;
+	    if (!(FLAG_GET(globals->candidate_flag, row, col)))
+		continue;
+
+	    /* get this ID */
+	    Segment_get(&globals->rid_seg, (void *) &this_id, row, col);
+
+	    reg_size = renumber[this_id];
+	    best_n_id = 1;
+
+	    while (reg_size < globals->min_segment_size && best_n_id > 0) {
+
+		Rbest = NULL;
+		reg_size = 1;
+		best_n_row = best_n_col = -1;
+
+		best_n_id = find_best_neighbour(globals, row, col, this_id,
+		                                nbtree, &reg_size, &Rbest,
+						&best_n_row, &best_n_col);
+
+		/* Rbest: pointer to most common neighbour
+		 * best_n_id, best_n_[row|col]: most similar neighbour */
+
+		if (reg_size < globals->min_segment_size && best_n_id > 0) {
+		    /* update rid */
+		    update_rid(globals, row, col, best_n_id);
+		    /* mark as merged */
+		    renumber[best_n_id] += renumber[this_id];
+		    reg_size = renumber[best_n_id];
+		    renumber[this_id] = 0;
+		    this_id = best_n_id;
+		}
+	    }
+	}
+    }
+    G_percent(1, 1, 1);
+
+    n_regions = 0;
+    /* renumber becomes new region ID */
+    for (i = 1; i <= globals->max_rid; i++) {
+	if (renumber[i] > 0) {
+	    renumber[i] = ++n_regions;
+	}
+    }
+
+    G_message(_("Renumbering remaining %d segments..."), n_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)))
+		continue;
+
+	    /* get this ID */
+	    Segment_get(&globals->rid_seg, (void *) &this_id, row, col);
+	    
+	    if (Rast_is_c_null_value(&this_id) || this_id < 1)
+		continue;
+
+	    this_id = renumber[this_id] + min_rid;
+	    Segment_put(&globals->rid_seg, (void *) &this_id, row, col);
+	}
+    }
+    G_percent(1, 1, 1);
+    
+    globals->max_rid = n_regions + min_rid;
+    G_free(renumber);
+    nbtree_clear(nbtree);
+
+    return 1;
+}

imagery/i.segment/parse_args.c

@@ -11,18 +11,24 @@ int parse_args(int argc, char *argv[], struct globals *globals)
 {
     struct Option *group, *seeds, *bounds, *output,
                   *method, *similarity, *threshold, *min_segment_size,
+		  *hs, *hr, *bsuf,
 #ifdef _OR_SHAPE_
 		  *shape_weight, *smooth_weight,
 #endif
 		   *mem;
-    struct Flag *diagonal, *weighted;
-    struct Option *outband, *endt;
+    struct Flag *diagonal, *weighted, *ms_a, *ms_p;
+    struct Option *gof, *endt;
 
     /* required parameters */
     group = G_define_standard_option(G_OPT_I_GROUP);
 
     output = G_define_standard_option(G_OPT_R_OUTPUT);
 
+    bsuf = G_define_standard_option(G_OPT_R_OUTPUT);
+    bsuf->key = "band_suffix";
+    bsuf->required = NO;
+    bsuf->label = _("Suffix for output bands with modified band values");
+
     threshold = G_define_option();
     threshold->key = "threshold";
     threshold->type = TYPE_DOUBLE;
@@ -32,6 +38,21 @@ int parse_args(int argc, char *argv[], struct globals *globals)
 
     /* optional parameters */
 
+    hs = G_define_option();
+    hs->key = "radius";
+    hs->type = TYPE_DOUBLE;
+    hs->required = NO;
+    hs->answer = "1.5";
+    hs->label = _("Spatial radius in number of cells");
+    hs->description = _("Must be >= 1, only cells within spatial bandwidth are considered for mean shift");
+
+    hr = G_define_option();
+    hr->key = "hr";
+    hr->type = TYPE_DOUBLE;
+    hr->required = NO;
+    hr->label = _("Range (spectral) bandwidth [0, 1]");
+    hr->description = _("Only cells within range (spectral) bandwidth are considered for mean shift. Range bandwidth is used as conductance parameter for adaptive bandwidth");
+
     method = G_define_option();
     method->key = "method";
     method->type = TYPE_STRING;
@@ -97,7 +118,6 @@ int parse_args(int argc, char *argv[], struct globals *globals)
     endt->key = "iterations";
     endt->type = TYPE_INTEGER;
     endt->required = NO;
-    endt->answer = "50";
     endt->description = _("Maximum number of iterations");
     endt->guisection = _("Settings");
 
@@ -107,6 +127,7 @@ int parse_args(int argc, char *argv[], struct globals *globals)
     seeds->key = "seeds";
     seeds->required = NO;
     seeds->description = _("Name for input raster map with starting seeds");
+    seeds->guisection = _("Settings");
 
     /* Polygon constraints. */
     bounds = G_define_standard_option(G_OPT_R_INPUT);
@@ -115,12 +136,14 @@ int parse_args(int argc, char *argv[], struct globals *globals)
     bounds->label = _("Name of input bounding/constraining raster map");
     bounds->description =
 	_("Must be integer values, each area will be segmented independent of the others");
+    bounds->guisection = _("Settings");
 
-    outband = G_define_standard_option(G_OPT_R_OUTPUT);
-    outband->key = "goodness";
-    outband->required = NO;
-    outband->description =
+    gof = G_define_standard_option(G_OPT_R_OUTPUT);
+    gof->key = "goodness";
+    gof->required = NO;
+    gof->description =
 	_("Name for output goodness of fit estimate map");
+    gof->guisection = _("Settings");
 
     diagonal = G_define_flag();
     diagonal->key = 'd';
@@ -134,6 +157,19 @@ int parse_args(int argc, char *argv[], struct globals *globals)
 	_("Weighted input, do not perform the default scaling of input raster maps");
     weighted->guisection = _("Settings");
 
+    ms_a = G_define_flag();
+    ms_a->key = 'a';
+    ms_a->label = _("Use adaptive bandwidth for mean shift");
+    ms_a->description = _("Range (spectral) bandwidth is adapted for each moving window");
+    ms_a->guisection = _("Settings");
+
+    ms_p = G_define_flag();
+    ms_p->key = 'p';
+    ms_p->label = _("Use progressive bandwidth for mean shift");
+    ms_p->description =
+	_("Spatial bandwidth is increased, range (spectral) bandwidth is decreased in each iteration");
+    ms_p->guisection = _("Settings");
+
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
@@ -146,10 +182,37 @@ int parse_args(int argc, char *argv[], struct globals *globals)
     else
 	G_fatal_error("Invalid output raster name");
 
+    globals->bsuf = bsuf->answer;
+
     globals->alpha = atof(threshold->answer);
     if (globals->alpha <= 0 || globals->alpha >= 1)
 	G_fatal_error(_("Threshold should be > 0 and < 1"));
 
+    globals->hs = -1;
+    if (hs->answer) {
+	globals->hs = atof(hs->answer);
+	if (globals->hs < 1) {
+	    G_fatal_error(_("Option '%s' must be >= 1"), hs->key);
+	}
+    }
+
+    globals->hr = -1;
+    if (hr->answer) {
+	globals->hr = atof(hr->answer);
+	if (globals->hr < 0) {
+	    G_warning(_("Negative value %s for option '%s': disabling"),
+	              hr->answer, hr->key);
+	    globals->hr = -1;
+	}
+	if (globals->hr >= 1) {
+	    G_warning(_("Value %s for option '%s' is >= 1: disabling"),
+	              hr->answer, hr->key);
+	    globals->hr = -1;
+	}
+    }
+    globals->ms_adaptive = ms_a->answer;
+    globals->ms_progressive = ms_p->answer;
+
     /* segmentation methods */
     if (strcmp(method->answer, "region_growing") == 0) {
 	globals->method = ORM_RG;
@@ -166,7 +229,7 @@ int parse_args(int argc, char *argv[], struct globals *globals)
     else
 	G_fatal_error(_("Unable to assign segmentation method"));
 
-    G_debug(1, "segmentation method: %s", method->answer);
+    G_debug(1, "segmentation method: %s (%d)", method->answer, globals->method);
 
     /* distance methods for similarity measurement */
     if (strcmp(similarity->answer, "euclidean") == 0)
@@ -257,26 +320,33 @@ int parse_args(int argc, char *argv[], struct globals *globals)
     }
 
     /* debug help */
-    if (outband->answer == NULL)
-	globals->out_band = NULL;
+    if (gof->answer == NULL)
+	globals->gof = NULL;
     else {
-	if (G_legal_filename(outband->answer) == TRUE)
-	    globals->out_band = outband->answer;
+	if (G_legal_filename(gof->answer) == TRUE)
+	    globals->gof = gof->answer;
 	else
 	    G_fatal_error(_("Invalid output raster name for goodness of fit"));
     }
 
-    if (endt->answer) {
+    if (!endt->answer) {
+	globals->end_t = 50;
+	if (globals->method == ORM_MS)
+	    globals->end_t = 10;
+	G_message(_("Maximum number of iterations set to %d"),
+		  globals->end_t);
+    }
+    else {
 	if (atoi(endt->answer) > 0)
 	    globals->end_t = atoi(endt->answer);
 	else {
 	    globals->end_t = 50;
+	    if (globals->method == ORM_MS)
+		globals->end_t = 10;
 	    G_warning(_("Invalid number of iterations, %d will be used"),
 	              globals->end_t);
 	}
     }
-    else
-	globals->end_t = 50;
 
     if (mem->answer && atoi(mem->answer) > 10)
 	globals->mb = atoi(mem->answer);

imagery/i.segment/region_growing.c

@@ -128,13 +128,6 @@ static int compare_rc(const void *first, const void *second)
     return 0;
 }
 
-static int compare_ints(const void *first, const void *second)
-{
-    int a = *(int *)first, b = *(int *)second;
-
-    return (a < b ? -1 : (a > b));
-}
-
 static int compare_double(double first, double second)
 {
     if (first < second)

imagery/i.segment/write_output.c

@@ -79,7 +79,7 @@ int write_gof_rg(struct globals *globals)
     struct History hist;
     DCELL *min, *max;
 
-    mean_fd = Rast_open_new(globals->out_band, FCELL_TYPE);
+    mean_fd = Rast_open_new(globals->gof, FCELL_TYPE);
     meanbuf = Rast_allocate_f_buf();
 
     /* goodness of fit for each cell: 1 = good fit, 0 = bad fit */
@@ -190,11 +190,11 @@ int write_gof_rg(struct globals *globals)
 
     Rast_init_colors(&colors);
     Rast_make_grey_scale_fp_colors(&colors, mingood, 1);
-    Rast_write_colors(globals->out_band, G_mapset(), &colors);
+    Rast_write_colors(globals->gof, G_mapset(), &colors);
 
-    Rast_short_history(globals->out_band, "raster", &hist);
+    Rast_short_history(globals->gof, "raster", &hist);
     Rast_command_history(&hist);
-    Rast_write_history(globals->out_band, &hist);
+    Rast_write_history(globals->gof, &hist);
 
     G_free(meanbuf);
 
@@ -227,12 +227,13 @@ int write_bands_ms(struct globals *globals)
     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);
+	G_asprintf(&name[n], "%s%s", globals->Ref.file[n].name, globals->bsuf);
 	out_fd[n] = Rast_open_new(name[n], DCELL_TYPE);
     }
     
     Rout.mean = G_malloc(globals->datasize);
 
+    G_message(_("Writing out shifted band values..."));
 
     for (row = 0; row < globals->nrows; row++) {