[GRASS-SVN] r68794 - sandbox/bo/i.segment.gsoc2016/i.segment
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Jun 29 10:32:55 PDT 2016
Author: hao2309
Date: 2016-06-29 10:32:55 -0700 (Wed, 29 Jun 2016)
New Revision: 68794
Modified:
sandbox/bo/i.segment.gsoc2016/i.segment/mean_shift.c
Log:
added clustering/clump function to mean_shift
Modified: sandbox/bo/i.segment.gsoc2016/i.segment/mean_shift.c
===================================================================
--- sandbox/bo/i.segment.gsoc2016/i.segment/mean_shift.c 2016-06-28 18:51:55 UTC (rev 68793)
+++ sandbox/bo/i.segment.gsoc2016/i.segment/mean_shift.c 2016-06-29 17:32:55 UTC (rev 68794)
@@ -5,6 +5,7 @@
#include <string.h>
#include <float.h>
#include <math.h>
+#include <fcntl.h>
#include <time.h>
#include <grass/gis.h>
#include <grass/glocale.h>
@@ -13,6 +14,8 @@
#include <grass/rbtree.h> /* Red Black Tree library functions */
#include "iseg.h"
+#define INCR 1024
+
/* standard gauss funtion:
* a * exp(-(x - m)^2 / (2 * stddev^2)
* a is not needed because the sum of weights is calculated for each
@@ -30,232 +33,380 @@
return exp(-diff2 / (2 * var));
}
-int clustering2(struct globals *globals){
- int class_val, px_count, first_change;/* the value of class */
- int row, col, nrows, ncols, stop_row, stop_col,edge_row, edge_col,diag;
- struct ngbr_stats Rin, Rstop;
- nrows = globals->row_max;
- ncols = globals->col_max;
- CELL temp_id;
- double cluster_threshold,cluster_threshold2;
- int **indicator_matrix;
- /* the indicator matrix to mark if the pixel has been clumped */
+static int cmp_cells(DCELL **a, int acol, DCELL **b, int bcol,
+ DCELL *rng, int n, double threshold2)
+{
+ int i;
+ double diff, diff2;
+
+ diff2 = 0;
+ for (i = 0; i < n; i++) {
+ if (Rast_is_d_null_value(&b[i][bcol]))
+ return 0;
+ diff = a[i][acol] - b[i][bcol];
+ /* normalize with the band's range */
+ if (rng[i])
+ diff /= rng[i];
+ diff2 += diff * diff;
+ }
+ /* normalize difference to the range [0, 1] */
+ diff2 /= n;
+
+ return (diff2 <= threshold2);
+}
+
+int clump(struct globals *globals)
+{
+ register int row, col,nin,diag;
+ register int i, n;
+ /* input */
+ DCELL **prev_in, **cur_in, **temp_in;
+ int bcol;
+ DCELL *rng, maxdiff;
+ double threshold,thresh2;
+ struct ngbr_stats Rin, Rabove, Rleft;
+ /* output */
+ CELL OLD, NEW;
+ CELL *temp_clump, *out_cell, *out_cellp;
+ CELL *prev_clump, *cur_clump;
+ CELL *index, *renumber;
+ CELL label;
+ int nrows, ncols;
+ int isnull;
+ int len;
+ int nalloc;
+ char *cname;
+ int cfd, csize;
+ CELL cat;
+ nin = globals->nbands;
+ threshold = globals->ms_range_bandwidth;
+ G_message(_("%d-band clumping with threshold %g"), nin, threshold);
+
Rin.mean = G_malloc(globals->datasize);
- Rstop.mean = G_malloc(globals->datasize);
- indicator_matrix = G_malloc(nrows * sizeof(int*));
- diag=0;
+ nrows = globals->row_max;
+ ncols = globals->col_max;
- for (row = globals->row_min; row < nrows; row++) {
- indicator_matrix[row] = G_malloc(ncols * sizeof(int));
- for (col = globals->col_min; col < ncols; col++) {
- indicator_matrix[row][col] = 0;
- }
- }
+ diag = globals->diagonal;
- G_message(_("Pass 1 of 2 for clustering"));
+ thresh2 = threshold * threshold;
+
+ /* allocate clump index */
+ nalloc = INCR;
+ index = (CELL *) G_malloc(nalloc * sizeof(CELL));
+ index[0] = 0;
+ renumber = NULL;
+
+ /* allocate DCELL buffers two columns larger than current window */
+ len = (ncols + 2) * sizeof(DCELL);
+ prev_in = (DCELL **) G_malloc(sizeof(DCELL *) * nin);
+ cur_in = (DCELL **) G_malloc(sizeof(DCELL *) * nin);
+ rng = G_malloc(sizeof(DCELL) * nin);
+
+ maxdiff = 0;
+ for (i = 0; i < nin; i++) {
+ DCELL min, max; /* min/max values of each input raster */
- px_count=0;
- class_val=0;
- stop_row=0;
- stop_col=0;
- cluster_threshold = 0.05;
- cluster_threshold2= pow(cluster_threshold ,2);
+ max = 0;
+ min =1;
- while (px_count<nrows*ncols){
- class_val+=1;
- Segment_get(globals->bands_out, (void *)Rstop.mean, stop_row, stop_col);
- first_change =0;/* mark the first time non-clump */
-
-
- /* get the first stop row when */
for (row = globals->row_min; row < nrows; row++) {
- G_percent(row - globals->row_min,
- globals->row_max - globals->row_min, 4);
for (col = globals->col_min; col < ncols; col++) {
if (!(FLAG_GET(globals->candidate_flag, row, col)))
continue;
-
Segment_get(globals->bands_out, (void *)Rin.mean, row, col);
- // G_message(_("row: %d, col: %d, temp_id is: %d"),row,col,temp_id);
-
- if ((globals->calculate_similarity)(&Rin, &Rstop, globals) <= cluster_threshold2 && indicator_matrix[row][col]==0){
- indicator_matrix[row][col]+=1;
-
- /* if the seg is not in the edg of the image */
- // if (row!= globals->row_min && col!=globals->col_min && row!=nrows-1 && col!=ncols-1){
-
- // /* If this seg is saperate from surrounded segments */
- // if (indicator_matrix[row+1][col]+indicator_matrix[row][col+1]+indicator_matrix[row-1][col]+indicator_matrix[row][col-1]==0){
- // class_val+=1;
-
- // if (diag){
- // if (indicator_matrix[row-1][col-1]+indicator_matrix[row+1][col-1]+indicator_matrix[row-1][col+1]+indicator_matrix[row+1][col+1]==0){
- // class_val+=1;
- // }
- // }
- // }
- // }
-
- Segment_put(&globals->rid_seg, &class_val, row, col);
-
- px_count+=1;
- }
- else if(first_change == 0 && indicator_matrix[row][col]==0){
- first_change = 1;
- stop_col = col;
- stop_row = row;
- }
- }
-
+ if (Rin.mean[i]<min)
+ min = Rin.mean[i];
+ if (Rin.mean[i]>max)
+ max = Rin.mean[i];
}
}
-
- return TRUE;
-}
+ rng[i] = max - min;
+ maxdiff += rng[i] * rng[i];
+ prev_in[i] = (DCELL *) G_malloc(len);
+ cur_in[i] = (DCELL *) G_malloc(len);
-int clustering(struct globals *globals){
- int row, col, n3, nin;
- int diag;/* if using diag pixel */
- int class_val;
- struct ngbr_stats Rin, Rleft, Rabove;
- CELL seg_id, left_new_val, above_new_val, *temp_clump;
- CELL old_val, new_val;
- double cluster_threshold, cluster_threshold2;
- /*******************************************************
- * Part 2 cluster similar values *
- * Refer clump.c to cluster super-pixels *
- *******************************************************/
-
- Rleft.mean = G_malloc(globals->datasize);
- Rin.mean = G_malloc(globals->datasize);
- Rabove.mean = G_malloc(globals->datasize);
-
+ /* fake a previous row which is all NULL */
+ Rast_set_d_null_value(prev_in[i], ncols + 2);
- G_percent_reset();
- nin = globals->nbands;
- diag=0;
- class_val=0;
- cluster_threshold = 0.05;
- cluster_threshold2= pow(cluster_threshold ,2);
- G_message(_("Pass 1 of 2 for clustering"));
- 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->candidate_flag, row, col)))
- continue;
- old_val=new_val=0;
-
- /* get the current pixel value */
+ /* set left and right edge to NULL */
+ Rast_set_d_null_value(&cur_in[i][0], 1);
+ Rast_set_d_null_value(&cur_in[i][ncols + 1], 1);
+ }
+ G_debug(1, "maximum possible difference: %g", maxdiff);
+
+ /* 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);
+ out_cell = (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);
+ label = 0;
+
+ /****************************************************
+ * PASS 1 *
+ * pass thru the input, create initial clump labels *
+ ****************************************************/
+
+ G_message(_("Pass 1 of 2..."));
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 2);
+ for (col = 0; col < ncols; col++){
Segment_get(globals->bands_out, (void *)Rin.mean, row, col);
-
- /* get left pixel value if not the first col*/
- if(col != globals->col_min){
- Segment_get(globals->bands_out, (void *)Rleft.mean, row, (col-1));
+ for (i = 0; i < nin; i++) {
+ cur_in[i][col+1] = Rin.mean[i];
}
-
- /* get above pixel value if not the first row*/
- if(row != globals->row_min){
- Segment_get(&globals->bands_seg2, (void *)Rabove.mean, row-1, col);
+ }
+
+ for (col = 1; col <= ncols; col++) {
+ isnull = 0;
+ for (i = 0; i < nin; i++) {
+ if (Rast_is_d_null_value(&cur_in[i][col])) { /* don't clump NULL data */
+ cur_clump[col] = 0;
+ isnull = 1;
+ break;
+ }
+ }
+ if (isnull)
+ continue;
+
+ /*
+ * 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;
+ /* same clump as to the left */
+ if (cmp_cells(cur_in, col, cur_in, col - 1, rng, nin, thresh2)) {
+ 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 {
+ if (cmp_cells(cur_in, col, prev_in, bcol, rng, nin, thresh2)) {
+ cur_clump[col] = *temp_clump;
+ if (OLD == 0) {
+ OLD = *temp_clump;
}
-
- if (col == globals->col_min && row == globals->row_min){
- class_val+=1;
- seg_id=class_val;
- }
- else if (row == globals->row_min){
- /* compare with left value and set if same cluster */
- Segment_get(&globals->bands_seg2, (void *)Rleft.mean, row, col-1);
- if ((globals->calculate_similarity)(&Rin, &Rleft, globals) <= cluster_threshold2){
- G_debug(1, "row: %d, col: %d, left_similarity is: %f",row,col,(globals->calculate_similarity)(&Rin, &Rleft, globals));
- Segment_get(&globals->rid_seg, (void *)&left_new_val, row, col-1);
- seg_id = left_new_val;
- }
- /* end compare left */
- else{
- class_val+=1;
- seg_id=class_val;
- }
- }
- else if(col == globals->col_min){
- /* compare with above value and set if same cluster */
- if ((globals->calculate_similarity)(&Rin, &Rabove, globals) <= cluster_threshold2){
- G_debug(1, "row: %d, col: %d, above_similarity is: %f",row,col,(globals->calculate_similarity)(&Rin, &Rabove, globals));
- Segment_get(&globals->rid_seg, (void *)&above_new_val, row-1, col);
- seg_id = above_new_val;
- }
- /* end compare above */
- else{
- class_val+=1;
- seg_id=class_val;
- }
- }
- else{
- if ((globals->calculate_similarity)(&Rin, &Rleft, globals) <= cluster_threshold2){
- Segment_get(&globals->rid_seg, (void *)&left_new_val, row, col-1);
- old_val=seg_id = left_new_val;
- // G_message(_("row: %d, col: %d, oldval is: %d, newval is: %d, compare to left"),row,col,old_val,new_val);
- }
-
- if(diag){
- // /* check above right, center, left, in that order */
- }
- else{
- /* compare with above value and set if same cluster */
- if ((globals->calculate_similarity)(&Rin, &Rabove, globals) <= cluster_threshold2){
-
- Segment_get(&globals->rid_seg, (void *)&above_new_val, row-1, col);
- temp_clump = &above_new_val;
- seg_id = *temp_clump;
- // G_message(_("row: %d, col: %d, oldval is: %d, newval is: %d, compare to above"),row,col,old_val,new_val);
- if (old_val==0){
- old_val = *temp_clump;
- }
- else{
- new_val = *temp_clump;
- if(new_val != old_val){
- /* 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.
- */
- // G_message(_("row: %d, col: %d, oldval is: %d, newval is: %d, conflict"),row,col,old_val,new_val);
- /* left of the current row from 1 to col - 1 */
- for(n3 = globals->col_min+1; n3< col; n3++){
- Segment_get(&globals->rid_seg, (void *)temp_clump, row, n3);
- if (*temp_clump==old_val){
- Segment_put(&globals->rid_seg, &new_val, row, n3);
- }
- }
- // /* right of previous row from col + 1 to ncols */
- for(n3 = col+1; n3< globals->col_max; n3++){
- Segment_get(&globals->rid_seg, (void *)temp_clump, row-1, n3);
- if (*temp_clump==old_val){
- Segment_put(&globals->rid_seg, &new_val, row-1, n3);
- }
- }
- }
+ 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;
}
+ /* (drive left, change all old to new)yb_note*/
+
+ /* 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;
+ }
}
- if (new_val = 0 || old_val == new_val){
- if (old_val ==0){
- class_val+=1;
- seg_id = class_val;
- // G_message(_("row: %d, col: %d, oldval is: %d, newval is: %d, class++"),row,col,old_val,new_val);
- }
+ }
+ temp_clump--;
+ } while (bcol-- > col - 1);
+ }
+ else {
+ /* check above */
+ if (cmp_cells(cur_in, col, prev_in, col, rng, nin, thresh2)) {
+ 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;
}
- Segment_put(&globals->rid_seg, &seg_id, row, col);
+ }
}
+ }
+
+ if (NEW == 0 || OLD == NEW) { /* ok */
+ if (OLD == 0) {
+ /* start a new clump */
+ 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 */
}
- return TRUE;
+
+ /* 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);
+
+ /* 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 = 1;
+ 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];
+ if (NEW == n)
+ G_fatal_error("Circular relabelling for %d", n);
+ }
+ index[n] = NEW;
+ }
+ else
+ /* set final clump id */
+ renumber[n] = cat++;
+ }
+
+ /* 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(_("Pass 2 of 2..."));
+ 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;
+ out_cellp = out_cell;
+
+ for (col = 0; col < ncols; col++) {
+ *out_cellp = renumber[index[*temp_clump]];
+ if (*out_cellp == 0)
+ Rast_set_c_null_value(out_cellp, 1);
+ temp_clump++;
+ out_cellp++;
+ }
+ for (col = 0; col < ncols; col++){
+ Segment_put(&globals->rid_seg, &out_cell[col], row, col);
+ }
+ }
+ G_percent(1, 1, 1);
+
+ close(cfd);
+ unlink(cname);
+
+
+ return TRUE;
}
-
int mean_shift(struct globals *globals)
{
int row, col, t, n1, n2;
@@ -417,7 +568,7 @@
G_message(_("Mean shift converged after %d iterations"), t);
- if (clustering(globals) != TRUE)
+ if (clump(globals) != TRUE)
G_fatal_error(_("Error in clustering segments"));
/*******************************************************
* Part 3 identify connected components *
More information about the grass-commit
mailing list