[GRASS-SVN] r69015 - sandbox/metz/i.segment.ms
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jul 21 06:03:39 PDT 2016
Author: mmetz
Date: 2016-07-21 06:03:39 -0700 (Thu, 21 Jul 2016)
New Revision: 69015
Added:
sandbox/metz/i.segment.ms/cluster.c
sandbox/metz/i.segment.ms/i.segment.ms.html
Removed:
sandbox/metz/i.segment.ms/i.segment.html
Modified:
sandbox/metz/i.segment.ms/Makefile
sandbox/metz/i.segment.ms/iseg.h
sandbox/metz/i.segment.ms/mean_shift.c
sandbox/metz/i.segment.ms/parse_args.c
sandbox/metz/i.segment.ms/region_growing.c
sandbox/metz/i.segment.ms/write_output.c
Log:
sandbox: i.segment + mean shift
Modified: sandbox/metz/i.segment.ms/Makefile
===================================================================
--- sandbox/metz/i.segment.ms/Makefile 2016-07-21 12:53:33 UTC (rev 69014)
+++ sandbox/metz/i.segment.ms/Makefile 2016-07-21 13:03:39 UTC (rev 69015)
@@ -1,6 +1,6 @@
MODULE_TOPDIR = ../..
-PGM = i.segment
+PGM = i.segment.ms
LIBES = $(IMAGERYLIB) $(RASTERLIB) $(SEGMENTLIB) $(GISLIB) $(BTREE2LIB)
DEPENDENCIES = $(IMAGERYDEP) $(RASTERDEP) $(SEGMENTDEP) $(GISDEP)
Added: sandbox/metz/i.segment.ms/cluster.c
===================================================================
--- sandbox/metz/i.segment.ms/cluster.c (rev 0)
+++ sandbox/metz/i.segment.ms/cluster.c 2016-07-21 13:03:39 UTC (rev 69015)
@@ -0,0 +1,456 @@
+
+/****************************************************************************
+ *
+ * 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(c) (((c) == 0 || (c) == ncols + 1 ? 1 : (FLAG_GET(globals->null_flag, row, (c - 1)))))
+#define CISNULLPREV(c) (row == 0 ? 1 : CISNULL(c))
+
+
+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(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(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 (!CISNULLPREV(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 (!CISNULLPREV(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;
+}
Property changes on: sandbox/metz/i.segment.ms/cluster.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Deleted: sandbox/metz/i.segment.ms/i.segment.html
===================================================================
--- sandbox/metz/i.segment.ms/i.segment.html 2016-07-21 12:53:33 UTC (rev 69014)
+++ sandbox/metz/i.segment.ms/i.segment.html 2016-07-21 13:03:39 UTC (rev 69015)
@@ -1,271 +0,0 @@
-<h2>DESCRIPTION</h2>
-Image segmentation or object recognition is the process of grouping
-similar pixels into unique segments, also referred 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 raster 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:
-
-<ol>
- <li>The pair is mutually most similar to each other (the similarity
-distance will be smaller than to any other neighbor), and</li>
- <li>The similarity must be lower than the input threshold. The
-process is repeated until no merges are made during a complete pass.</li>
-</ol>
-
-<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 given
-threshold value. During the final pass, however, if a minimum
-segment size of 2 or larger is given with the <b>minsize</b>
-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 <b>threshold</b> must be larger than 0.0 and smaller than 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. It is recommended to start with a low value, e.g. 0.01, and
-then perform hierarchical segmentation by using the output of the last
-run as <b>seeds</b> for the next run.
-
-<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: <tt>(1-radioweight) * {smoothness * smoothness
-weight + compactness * (1-smoothness weight)}</tt>,
-where <tt>compactness = Perimeter Length / sqrt( Area )</tt>
-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>
-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
-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.
-<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>
-
-<h3>Boundary Constraints</h3>
-Boundary constraints limit the adjacency of pixels and segments.
-Each unique value present in the <b>bounds</b> 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 <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.
-
-<h2>EXAMPLES</h2>
-
-<h3>Segmentation of RGB orthophoto</h3>
-
-This example uses the ortho photograph included in the NC Sample
-Dataset. Set up an imagery group:
-<div class="code"><pre>
-i.group group=ortho_group input=ortho_2001_t792_1m at PERMANENT
-</pre></div>
-
-<p>Set the region to a smaller test region (resolution taken from
-input ortho photograph).
-
-<div class="code"><pre>
-g.region -p raster=ortho_2001_t792_1m n=220446 s=220075 e=639151 w=638592
-</pre></div>
-
-Try out a low threshold and check the results.
-<div class="code"><pre>
-i.segment group=ortho_group output=ortho_segs_l1 threshold=0.02
-</pre></div>
-
-<center>
-<img src="i_segment_ortho_segs_l1.jpg">
-</center>
-
-<p>
-From a visual inspection, it seems this results in too many segments.
-Increasing the threshold, using the previous results as seeds,
-and setting a minimum size of 2:
-<div class="code"><pre>
-i.segment group=ortho_group output=ortho_segs_l2 threshold=0.05 seeds=ortho_segs_l1 min=2
-
-i.segment group=ortho_group output=ortho_segs_l3 threshold=0.1 seeds=ortho_segs_l2
-
-i.segment group=ortho_group output=ortho_segs_l4 threshold=0.2 seeds=ortho_segs_l3
-
-i.segment group=ortho_group output=ortho_segs_l5 threshold=0.3 seeds=ortho_segs_l4
-</pre></div>
-
-<center>
-<img src="i_segment_ortho_segs_l2_l5.jpg">
-</center>
-
-<p>
-The output <tt>ortho_segs_l4</tt> with <b>threshold</b>=0.2 still has
-too many segments, but the output with <b>threshold</b>=0.3 has too few
-segments. A threshold value of 0.25 seems to be a good choice. There
-is also some noise in the image, lets next force all segments smaller
-than 10 pixels to be merged into their most similar neighbor (even if
-they are less similar than required by our threshold):
-
-<p>Set the region to match the entire map(s) in the group.
-<div class="code"><pre>
-g.region -p raster=ortho_2001_t792_1m at PERMANENT
-</pre></div>
-
-<p>
-Run <em>i.segment</em> on the full map:
-
-<div class="code"><pre>
-i.segment group=ortho_group output=ortho_segs_final threshold=0.25 min=10
-</pre></div>
-
-<center>
-<img src="i_segment_ortho_segs_final.jpg">
-</center>
-
-<p>
-Processing the entire ortho image with nearly 10 million pixels took
-about 450 times more then for the final run.
-
-<h3>Segmentation of panchromatic channel</h3>
-
-This example uses the panchromatic channel of the Landsat7 scene included
-in the North Carolina sample dataset:
-
-<div class="code"><pre>
-# create group with single channel
-i.group group=singleband input=lsat7_2002_80
-
-# set computational region to Landsat7 PAN band
-g.region raster=lsat7_2002_80 -p
-
-# perform segmentation with minsize=5
-i.segment group=singleband threshold=0.05 minsize=5 \
- output=lsat7_2002_80_segmented_min5 goodness=lsat7_2002_80_goodness_min5
-
-# perform segmentation with minsize=100
-i.segment group=singleband threshold=0.05 minsize=100
- output=lsat7_2002_80_segmented_min100 goodness=lsat7_2002_80_goodness_min100
-</pre></div>
-
-<p>
-<center>
-<img src="i_segment_lsat7_pan.png"><br>
-Original panchromatic channel of the Landsat7 scene
-</center>
-
-<p>
-<center>
-<img src="i_segment_lsat7_seg_min5.png"><br>
-Segmented panchromatic channel, minsize=5
-</center>
-<p>
-<center>
-<img src="i_segment_lsat7_seg_min100.png"><br>
-Segmented panchromatic channel, minsize=100
-</center>
-
-<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 <em><a href="i.maxlik.html">i.maxlik</a></em>
-to ensure the segmentation output can be used as input for the
-existing classification functionality.</li>
-<li>Integration/workflow for <em>r.fuzzy</em> (Addon).</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://grasswiki.osgeo.org/wiki/GRASS_GSoC_2012_Image_Segmentation">project wiki</a>.
-<p>
-Information about
-<a href="http://grasswiki.osgeo.org/wiki/Image_classification">classification in GRASS</a>
-is at available on the wiki.
-
-<h2>SEE ALSO</h2>
-<em>
-<a href="g.gui.iclass.html">g.gui.iclass</a>,
-<a href="i.group.html">i.group</a>,
-<a href="i.maxlik.html">i.maxlik</a>,
-<a href="i.smap.html">i.smap</a>,
-<a href="r.kappa.html">r.kappa</a>
-</em>
-
-<h2>AUTHORS</h2>
-Eric Momsen - North Dakota State University<br>
-Markus Metz (GSoC Mentor)
-
-<p>
-<i>Last changed: $Date$</i>
-
Copied: sandbox/metz/i.segment.ms/i.segment.ms.html (from rev 69014, sandbox/metz/i.segment.ms/i.segment.html)
===================================================================
--- sandbox/metz/i.segment.ms/i.segment.ms.html (rev 0)
+++ sandbox/metz/i.segment.ms/i.segment.ms.html 2016-07-21 13:03:39 UTC (rev 69015)
@@ -0,0 +1,274 @@
+<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 raster 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:
+
+<ol>
+ <li>The pair is mutually most similar to each other (the similarity
+distance will be smaller than to any other neighbor), and</li>
+ <li>The similarity must be lower than the input threshold. The
+process is repeated until no merges are made during a complete pass.</li>
+</ol>
+
+<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.
+<p>
+During normal processing, merges are only allowed when the
+similarity between two segments is lower than the given
+threshold value. During the final pass, however, if a minimum
+segment size of 2 or larger is given with the <b>minsize</b>
+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 <b>threshold</b> must be larger than 0.0 and smaller than 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. It is recommended to start with a low value, e.g. 0.01, and
+then perform hierachical segmentation by using the output of the last
+run as <b>seeds</b> for the next run.
+
+<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: <tt>(1-radioweight) * {smoothness * smoothness
+weight + compactness * (1-smoothness weight)}</tt>,
+where <tt>compactness = Perimeter Length / sqrt( Area )</tt>
+and <tt>smoothness = Perimeter Length / Bounding Box</tt>. The
+perimeter length is estimated as the number of pixel sides the segment
+has.
+
+<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
+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>.
+
+<h3>Boundary Constraints</h3>
+Boundary constraints limit the adjacency of pixels and segments.
+Each unique value present in the <b>bounds</b> 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 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.
+
+<h2>EXAMPLES</h2>
+
+<h3>Segmentation of RGB orthophoto</h3>
+
+This example uses the ortho photograph included in the NC Sample
+Dataset. Set up an imagery group:
+<div class="code"><pre>
+i.group group=ortho_group input=ortho_2001_t792_1m at PERMANENT
+</pre></div>
+
+<p>Set the region to a smaller test region (resolution taken from
+input ortho photograph).
+
+<div class="code"><pre>
+g.region -p raster=ortho_2001_t792_1m n=220446 s=220075 e=639151 w=638592
+</pre></div>
+
+Try out a low threshold and check the results.
+<div class="code"><pre>
+i.segment group=ortho_group output=ortho_segs_l1 threshold=0.02
+</pre></div>
+
+<center>
+<img src="i_segment_ortho_segs_l1.jpg">
+</center>
+
+<p>
+From a visual inspection, it seems this results in too many segments.
+Increasing the threshold, using the previous results as seeds,
+and setting a minimum size of 2:
+<div class="code"><pre>
+i.segment group=ortho_group output=ortho_segs_l2 threshold=0.05 seeds=ortho_segs_l1 min=2
+
+i.segment group=ortho_group output=ortho_segs_l3 threshold=0.1 seeds=ortho_segs_l2
+
+i.segment group=ortho_group output=ortho_segs_l4 threshold=0.2 seeds=ortho_segs_l3
+
+i.segment group=ortho_group output=ortho_segs_l5 threshold=0.3 seeds=ortho_segs_l4
+</pre></div>
+
+<center>
+<img src="i_segment_ortho_segs_l2_l5.jpg">
+</center>
+
+<p>
+The output <tt>ortho_segs_l4</tt> with <b>threshold</b>=0.2 still has
+too many segments, but the output with <b>threshold</b>=0.3 has too few
+segments. A threshold value of 0.25 seems to be a good choice. There
+is also some noise in the image, lets next force all segments smaller
+than 10 pixels to be merged into their most similar neighbor (even if
+they are less similar than required by our threshold):
+
+<p>Set the region to match the entire map(s) in the group.
+<div class="code"><pre>
+g.region -p raster=ortho_2001_t792_1m at PERMANENT
+</pre></div>
+
+<p>
+Run <em>i.segment</em> on the full map:
+
+<div class="code"><pre>
+i.segment group=ortho_group output=ortho_segs_final threshold=0.25 min=10
+</pre></div>
+
+<center>
+<img src="i_segment_ortho_segs_final.jpg">
+</center>
+
+<p>
+Processing the entire ortho image with nearly 10 million pixels took
+about 450 times more then for the final run.
+
+<h3>Segmentation of panchromatic channel</h3>
+
+This example uses the panchromatic channel of the Landsat7 scene included
+in the North Carolina sample dataset:
+
+<div class="code"><pre>
+# create group with single channel
+i.group group=singleband input=lsat7_2002_80
+
+# set computational region to Landsat7 PAN band
+g.region raster=lsat7_2002_80 -p
+
+# perform segmentation with minsize=5
+i.segment group=singleband threshold=0.05 minsize=5 \
+ output=lsat7_2002_80_segmented_min5 goodness=lsat7_2002_80_goodness_min5
+
+# perform segmentation with minsize=100
+i.segment group=singleband threshold=0.05 minsize=100
+ output=lsat7_2002_80_segmented_min100 goodness=lsat7_2002_80_goodness_min100
+</pre></div>
+
+<p>
+<center>
+<img src="i_segment_lsat7_pan.png"><br>
+Original panchromatic channel of the Landsat7 scene
+</center>
+
+<p>
+<center>
+<img src="i_segment_lsat7_seg_min5.png"><br>
+Segmented panchromatic channel, minsize=5
+</center>
+<p>
+<center>
+<img src="i_segment_lsat7_seg_min100.png"><br>
+Segmented panchromatic channel, minsize=100
+</center>
+
+<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 <em><a href="i.maxlik.html">i.maxlik</a></em>
+to ensure the segmentation output can be used as input for the
+existing classification functionality.</li>
+<li>Integration/workflow for <em>r.fuzzy</em> (Addon).</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://grasswiki.osgeo.org/wiki/GRASS_GSoC_2012_Image_Segmentation">project wiki</a>.
+<p>
+Information about
+<a href="http://grasswiki.osgeo.org/wiki/Image_classification">classification in GRASS</a>
+is at available on the wiki.
+
+<h2>SEE ALSO</h2>
+<em>
+<a href="g.gui.iclass.html">g.gui.iclass</a>,
+<a href="i.group.html">i.group</a>,
+<a href="i.maxlik.html">i.maxlik</a>,
+<a href="i.smap.html">i.smap</a>,
+<a href="r.kappa.html">r.kappa</a>
+</em>
+
+<h2>AUTHORS</h2>
+Eric Momsen - North Dakota State University<br>
+Markus Metz (GSoC Mentor)
+
+<p>
+<i>Last changed: $Date$</i>
+
Modified: sandbox/metz/i.segment.ms/iseg.h
===================================================================
--- sandbox/metz/i.segment.ms/iseg.h 2016-07-21 12:53:33 UTC (rev 69014)
+++ sandbox/metz/i.segment.ms/iseg.h 2016-07-21 13:03:39 UTC (rev 69015)
@@ -88,6 +88,8 @@
/* 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 +161,9 @@
/* mean_shift.c */
int mean_shift(struct globals *);
+/* cluster.c */
+CELL cluster_bands(struct globals *globals);
+
/* watershed.c */
int watershed(struct globals *);
Modified: sandbox/metz/i.segment.ms/mean_shift.c
===================================================================
--- sandbox/metz/i.segment.ms/mean_shift.c 2016-07-21 12:53:33 UTC (rev 69014)
+++ sandbox/metz/i.segment.ms/mean_shift.c 2016-07-21 13:03:39 UTC (rev 69015)
@@ -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 @@
int mean_shift(struct globals *globals)
{
- int row, col, t;
- int n_changes;
- double alpha2;
- struct ngbr_stats Rin, Rout;
- double diff2;
+ 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;
- G_fatal_error(_("Mean shift is not yet implemented"));
- return FALSE;
-
-
Rin.mean = G_malloc(globals->datasize);
Rout.mean = G_malloc(globals->datasize);
+ Rn.mean = G_malloc(globals->datasize);
- /* TODO: need another segment structure holding output
- * for mean shift, output of the previous iteration becomes input of the next iteration
- * initially, input and output are original band values */
-
alpha2 = globals->alpha * globals->alpha;
+ 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.7;
+ 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, ®_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;
+}
Modified: sandbox/metz/i.segment.ms/parse_args.c
===================================================================
--- sandbox/metz/i.segment.ms/parse_args.c 2016-07-21 12:53:33 UTC (rev 69014)
+++ sandbox/metz/i.segment.ms/parse_args.c 2016-07-21 13:03:39 UTC (rev 69015)
@@ -11,11 +11,12 @@
{
struct Option *group, *seeds, *bounds, *output,
*method, *similarity, *threshold, *min_segment_size,
+ *hs, *hr, *ms_suf,
#ifdef _OR_SHAPE_
*shape_weight, *smooth_weight,
#endif
*mem;
- struct Flag *diagonal, *weighted;
+ struct Flag *diagonal, *weighted, *ms_a, *ms_p;
struct Option *outband, *endt;
/* required parameters */
@@ -23,6 +24,12 @@
output = G_define_standard_option(G_OPT_R_OUTPUT);
+ ms_suf = G_define_standard_option(G_OPT_R_OUTPUT);
+ ms_suf->key = "ms";
+ ms_suf->required = NO;
+ ms_suf->answer = "_ms";
+ ms_suf->label = _("Suffix for output bands (mean shift)");
+
threshold = G_define_option();
threshold->key = "threshold";
threshold->type = TYPE_DOUBLE;
@@ -32,6 +39,20 @@
/* optional parameters */
+ hs = G_define_option();
+ hs->key = "hs";
+ hs->type = TYPE_DOUBLE;
+ hs->required = NO;
+ hs->label = _("Spatial bandwidth 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 @@
endt->key = "iterations";
endt->type = TYPE_INTEGER;
endt->required = NO;
- endt->answer = "50";
endt->description = _("Maximum number of iterations");
endt->guisection = _("Settings");
@@ -134,6 +154,19 @@
_("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 +179,37 @@
else
G_fatal_error("Invalid output raster name");
+ globals->ms_suf = ms_suf->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 +226,7 @@
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)
@@ -266,17 +326,24 @@
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);
Modified: sandbox/metz/i.segment.ms/region_growing.c
===================================================================
--- sandbox/metz/i.segment.ms/region_growing.c 2016-07-21 12:53:33 UTC (rev 69014)
+++ sandbox/metz/i.segment.ms/region_growing.c 2016-07-21 13:03:39 UTC (rev 69015)
@@ -128,13 +128,6 @@
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)
Modified: sandbox/metz/i.segment.ms/write_output.c
===================================================================
--- sandbox/metz/i.segment.ms/write_output.c 2016-07-21 12:53:33 UTC (rev 69014)
+++ sandbox/metz/i.segment.ms/write_output.c 2016-07-21 13:03:39 UTC (rev 69015)
@@ -233,6 +233,7 @@
Rout.mean = G_malloc(globals->datasize);
+ G_message(_("Writing out shifted band values..."));
for (row = 0; row < globals->nrows; row++) {
More information about the grass-commit
mailing list