[GRASS-SVN] r52217 - grass-addons/grass7/imagery/i.segment

svn_grass at osgeo.org svn_grass at osgeo.org
Mon Jun 25 20:33:26 PDT 2012


Author: momsen
Date: 2012-06-25 20:33:25 -0700 (Mon, 25 Jun 2012)
New Revision: 52217

Added:
   grass-addons/grass7/imagery/i.segment/times.txt
Modified:
   grass-addons/grass7/imagery/i.segment/create_isegs.c
   grass-addons/grass7/imagery/i.segment/iseg.h
   grass-addons/grass7/imagery/i.segment/main.c
   grass-addons/grass7/imagery/i.segment/open_files.c
   grass-addons/grass7/imagery/i.segment/parse_args.c
   grass-addons/grass7/imagery/i.segment/write_output.c
Log:
added mask/null handling of input bands, updated segmentation algorithm (delay candidate flag check until after find mutually best neighbor), general clean up and TODO resolution.

Modified: grass-addons/grass7/imagery/i.segment/create_isegs.c
===================================================================
--- grass-addons/grass7/imagery/i.segment/create_isegs.c	2012-06-26 03:22:09 UTC (rev 52216)
+++ grass-addons/grass7/imagery/i.segment/create_isegs.c	2012-06-26 03:33:25 UTC (rev 52217)
@@ -373,16 +373,24 @@
 	/*process candidate pixels */
 	G_verbose_message("Row percent complete for pass number %d: ", t);
 	/*check each pixel, start the processing only if it is a candidate pixel */
+	
+	/* for validation, select one of the two... could make this IFDEF or input parameter */
+	/* reverse order 
+	for (row = files->nrows; row >= 0; row--) {
+	    for (col = files->ncols; col >= 0 ; col--) {
+		G_percent(files->nrows-row, files->nrows, 1);
+end reverse order */
+/* "normal" order */
 	for (row = 0; row < files->nrows; row++) {
 	    for (col = 0; col < files->ncols; col++) {
+		G_percent(row, files->nrows, 1); /*end normal order */	/* TODO, can a message be included with G_percent? */
 
-		G_percent(row, files->nrows, 1);	/* TODO, can a message be included with G_percent? */
 
 		G_debug(4,
 			"Next starting pixel from next row/col, not from Rk");
 
 		segment_get(&files->out_seg, (void *)files->out_val, row, col);	/*TODO small time savings - if candidate_count reaches zero, bail out of these loops too? */
-		if (files->out_val[1] == 1) {	/* out_val[1] is the candidate pixel flag, want to process the 1's */
+		if (files->out_val[1] == 1 && files->out_val != NULL) {	/* out_val[1] is the candidate pixel flag, want to process the 1's */ /*TODO MASK handling - checking if we have a segment ID already, in open_files() we will put nulls in [0] slot... better/faster way to do this?*/
 		    G_debug(4, "going to free memory on linked lists...");
 		    /*free memory for linked lists */
 		    my_dispose_list(files->token, &Ri_head);
@@ -472,6 +480,7 @@
 			    G_debug(4,
 				    "simularity = %g for neighbor : row: %d, col %d.",
 				    tempsim, current->row, current->col);
+				    
 			    if (tempsim < Ri_similarity) {
 				Ri_similarity = tempsim;
 				Ri_bestn = current;	/*TODO want to point to the current pixel...confirm  when current changes need this to stay put! */
@@ -482,12 +491,19 @@
 			    }
 			}
 
-			if (Ri_bestn != NULL)
+			if (Ri_bestn != NULL){
 			    G_debug(4,
 				    "Lowest Ri_similarity = %g, for neighbor pixel row: %d col: %d",
 				    Ri_similarity, Ri_bestn->row,
 				    Ri_bestn->col);
 
+			    segment_get(&files->out_seg, (void *)files->out_val, Ri_bestn->row, Ri_bestn->col);
+			    if (files->out_val[1] == 0)
+				/* this check is important:
+				 * best neighbor is not a valid candidate, was already merged earlier in this time step */
+				Ri_similarity = threshold + 1;
+}
+
 			if (Ri_bestn != NULL && Ri_similarity < threshold) {	/* small TODO: should this be < or <= for threshold? */
 			    /* we'll have the neighbor pixel to start with. */
 			    G_debug(4, "3a: Working with Rk");
@@ -526,6 +542,7 @@
 
 			    for (current = Rkn_head; current != NULL; current = current->next) {	/* for each of Rk's neighbors */
 				tempsim = functions->calculate_similarity(Rk_head, current, files, functions);	/*TODO: need an error trap here, if something goes wrong with calculating similarity? */
+				
 				if (tempsim < Rk_similarity) {
 				    Rk_similarity = tempsim;
 				    break;	/* exit for Rk's neighbors loop here, we know that Ri and Rk aren't mutually best neighbors */
@@ -559,9 +576,17 @@
 
 			    }
 			}	/*end else (from if mutually best neighbors) */
-			else
+			else {
+			    /* no valid best neighbor for this Ri
+			     * exclude this Ri from further comparisons 
+			     * because we checked already Ri for a mutually best neighbor with all valid candidates
+			     * thus Ri can not be the mutually best neighbor later on during this pass
+			     * unfortunately this does happen sometimes */
+			    set_candidate_flag(Ri_head, 0, files);	/* TODO: error trap? */
+			    files->candidate_count++; /*first pixel was already set*/
 			    G_debug(4,
-				    "3b Rk didn't didn't exist or similarity was > threshold");
+				    "3b Rk didn't didn't exist, was not valid candidate, or similarity was > threshold");
+				} /*end else - Ri's best neighbor was not a candidate */
 		    }		/* end else - Ri did have neighbors */
 		    //          }           /*end pathflag do loop */
 		}		/*end if pixel is candidate pixel */
@@ -572,7 +597,7 @@
 
 	G_debug(4, "Finished one pass, t was = %d", t);
 	t++;
-    } while (t < 90 && endflag == 0);
+    } while (t <= functions->end_t && endflag == 0);
     /*end t loop *//*TODO, should there be a max t that it can iterate for?  Include t in G_message? */
 
     /* free memory *//*TODO: anything ? */
@@ -689,7 +714,8 @@
 
 		segment_get(&files->out_seg, (void *)files->out_val, pixel_neighbors[n][0], pixel_neighbors[n][1]);	/*TODO : do I need a second "out_val" data structure? */
 
-		if (files->out_val[1] == 1) {	/* valid candidate pixel */
+		if (files->out_val[1] == 1 || files->out_val[1] == 0) {	/* all pixels, not just valid pixels */
+		/* TODO: use -1 for NULL/MASKED pixels? */
 
 		    G_debug(5, "\tfiles->out_val[0] = %d Ri_seg_ID = %d",
 			    files->out_val[0], Ri_seg_ID);

Modified: grass-addons/grass7/imagery/i.segment/iseg.h
===================================================================
--- grass-addons/grass7/imagery/i.segment/iseg.h	2012-06-26 03:22:09 UTC (rev 52216)
+++ grass-addons/grass7/imagery/i.segment/iseg.h	2012-06-26 03:33:25 UTC (rev 52217)
@@ -34,8 +34,9 @@
     int nrows, ncols;
 
     /* files *//* TODO, for all map names, is this better for any reason, saw it in manual example: char name[GNAME_MAX]; */
-    char *out_name;		/* name of output raster map */
+    char *out_name;		/* name of output raster map */  
     char *seeds, *bounds_map, *bounds_mapset;	/* optional segment seeds and polygon constraints/boundaries */
+	char *out_band;		/* for debug */
 
     /* file processing */
     int nbands;			/* number of rasters in the image group */
@@ -68,9 +69,7 @@
      */
 
     /* Todo: Additional data to be stored in this structure:
-     * vector constraints
      * seeds (Put directly into output map?)
-     * Need to consider if they should be put in RAM or SEGMENT?
      */
 
 };
@@ -85,6 +84,8 @@
     int (*find_pixel_neighbors) (int, int, int[8][2], struct files *);	/*parameters: row, col, pixel_neighbors */
     double (*calculate_similarity) (struct pixels *, struct pixels *, struct files *, struct functions *);	/*parameters: two points (row,col) to compare */
 
+	/* for debug */
+	int end_t;
 };
 
 
@@ -101,9 +102,9 @@
 int ll_test(struct files *, struct functions *);
 int test_pass_token(struct pixels **, struct files *);
 int region_growing(struct files *, struct functions *);
-int find_segment_neighbors(struct pixels **, struct pixels **, int *, struct files *, struct functions *);	/* TODO: need data structure for Ri, Rin */
+int find_segment_neighbors(struct pixels **, struct pixels **, int *, struct files *, struct functions *);
 int set_candidate_flag(struct pixels *, int, struct files *);
-int merge_values(struct pixels *, struct pixels *, int, int, struct files *);	/* I assume this is a weighted mean? */
+int merge_values(struct pixels *, struct pixels *, int, int, struct files *);
 int find_four_pixel_neighbors(int, int, int[][2], struct files *);
 int find_eight_pixel_neighbors(int, int, int[8][2], struct files *);
 double calculate_euclidean_similarity(struct pixels *, struct pixels *,

Modified: grass-addons/grass7/imagery/i.segment/main.c
===================================================================
--- grass-addons/grass7/imagery/i.segment/main.c	2012-06-26 03:22:09 UTC (rev 52216)
+++ grass-addons/grass7/imagery/i.segment/main.c	2012-06-26 03:33:25 UTC (rev 52217)
@@ -3,7 +3,7 @@
  *
  * MODULE:       i.segment
  * AUTHOR(S):    Eric Momsen <eric.momsen at gmail com>
- * PURPOSE:      Provide short description of module here...
+ * PURPOSE:      Segments an image group.
  * COPYRIGHT:    (C) 2012 by Eric Momsen, and the GRASS Development Team
  *
  *               This program is free software under the GNU General Public
@@ -20,7 +20,7 @@
 
 #include <stdlib.h>
 #include <grass/gis.h>
-#include <grass/glocale.h>	/* message translation */
+#include <grass/glocale.h>
 #include "iseg.h"
 
 int main(int argc, char *argv[])
@@ -35,7 +35,7 @@
     G_add_keyword(_("imagery"));
     G_add_keyword(_("segmentation"));
     module->description =
-	_("Outputs a single segmention map (raster) based on input values in an image group.");
+	_("Outputs a single segmented map (raster) based on input values in an image group.");
 
     if (parse_args(argc, argv, &files, &functions) != 0)
 	G_debug(1, "Error in parse_args()");	/* TODO: should these be debug or G_fatal_error() or nested if statement? want to clean up mem and temp files */

Modified: grass-addons/grass7/imagery/i.segment/open_files.c
===================================================================
--- grass-addons/grass7/imagery/i.segment/open_files.c	2012-06-26 03:22:09 UTC (rev 52216)
+++ grass-addons/grass7/imagery/i.segment/open_files.c	2012-06-26 03:33:25 UTC (rev 52217)
@@ -10,8 +10,8 @@
 int open_files(struct files *files)
 {
     struct Ref Ref;		/* group reference list */
-    int *in_fd, bounds_fd;
-    int n, row, col, srows, scols, inlen, outlen, nseg;
+    int *in_fd, bounds_fd, null_check;
+    int n, s, row, col, srows, scols, inlen, outlen, nseg;
     DCELL **inbuf;		/* buffer array, to store lines from each of the imagery group rasters */
     CELL *boundsbuf;
     struct FPRange *fp_range;	/* for getting min/max values on each input raster */
@@ -22,8 +22,7 @@
 
     /* ****** open the input rasters ******* */
 
-    /* TODO, this is from i.smap/openfiles.c  lines 17-23 checked if subgroup had maps, does API handles the checks? */
-
+	/* TODO: I confirmed, the API does not check this.  Should this be checked in parse_args / validation ?  Or best to do here where I want to use the REF data? */
     if (!I_get_group_ref(files->image_group, &Ref))
 	G_fatal_error(_("Unable to read REF file for group <%s>"),
 		      files->image_group);
@@ -52,10 +51,10 @@
 	for (n = 0; n < Ref.nfiles; n++) {
 	    if (Rast_read_fp_range
 		(Ref.file[n].name, Ref.file[n].mapset, &fp_range[n]) < 0)
-		G_fatal_error(_("Unable to read fp range for raster map <%s>"), Ref.file[n].name);	/* TODO, still should close files? */
+		G_fatal_error(_("Unable to read fp range for raster map <%s>"), Ref.file[n].name);	/* TODO, still should free memory and close files? */
 	    Rast_get_fp_range_min_max(&(fp_range[n]), &min[n], &max[n]);
 	}
-	G_verbose_message("scaling, for first layer, min: %f, max: %f",
+	G_debug(1, "scaling, for first layer, min: %f, max: %f",
 			  min[0], max[0]);
     }
 
@@ -106,46 +105,51 @@
     if (segment_open(&files->no_check, G_tempfile(), files->nrows, files->ncols, srows, scols, sizeof(int), nseg) != 1)	/* todo could make this smaller ? just need 0 or 1 */
 	G_fatal_error("Unable to create flag temporary files");
 
-    /* load input bands to segment structure */
+    /* load input bands to segment structure and initialize output segmentation file */
     G_debug(1, "Reading input rasters into segmentation data files...");
 
     files->bands_val = (double *)G_malloc(inlen);
     files->second_val = (double *)G_malloc(inlen);
-
+    files->out_val = (int *)G_malloc(2 * sizeof(int));
+	s = 1; /* initial segment ID */
+	
     for (row = 0; row < files->nrows; row++) {
 	for (n = 0; n < Ref.nfiles; n++) {
 	    Rast_get_d_row(in_fd[n], inbuf[n], row);
 	}
 	for (col = 0; col < files->ncols; col++) {
+		/*tempval = 0; Doesn't work, no "null" for doubles in c */ /* want a number, not null */
+		null_check = 1; /*Assume there is data*/
 	    for (n = 0; n < Ref.nfiles; n++) {
+		/*tempval += inbuf[n][col]; */ /* if mask/null, adding a null value should set tempval to NULL */
+		if(Rast_is_d_null_value(&inbuf[n][col]))
+			null_check=-1;
 		if (files->weighted == 1)
 		    files->bands_val[n] = inbuf[n][col];	/*unscaled */
 		else
 		    files->bands_val[n] = (inbuf[n][col] - min[n]) / (max[n] - min[n]);	/*scaled version */
 	    }
 	    segment_put(&files->bands_seg, (void *)files->bands_val, row,
-			col);
-	}
-    }
-
-    /* Initialize output segmentation file */
-
-    /* TODO: Need to carry MASK over to the output file as well.  Need an additional flag that is set, during input read?  */
-
-    G_debug(1, "G_malloc size: = <%d>", (int)(2 * sizeof(int)));
-    files->out_val = (int *)G_malloc(2 * sizeof(int));
-
-    n = 1;
-    for (row = 0; row < files->nrows; row++) {
-	for (col = 0; col < files->ncols; col++) {
-	    files->out_val[0] = n;	/*starting segment number TODO: for seeds this will be different */
+			col); /* store input bands */
+			
+		if (null_check != -1){ /*good pixel*/
+		files->out_val[0] = s;	/*starting segment number TODO: for seeds this will be different */
 	    files->out_val[1] = 0;	/*flag */
-	    segment_put(&files->out_seg, (void *)files->out_val, row, col);
-	    n++;		/* sequentially number all pixels with their own segment ID */
 	}
+	else /*don't use this pixel*/
+	{
+		files->out_val[0] = -1;	/*starting segment number*/
+	    files->out_val[1] = -1;	/*flag */
+	}
+	    segment_put(&files->out_seg, (void *)files->out_val, row, col); /* initialize input */
+	    s++;		/* sequentially number all pixels with their own segment ID */
+	}
     }
 
     /* bounds/constraints */
+    /* TODO: You should also handle NULL cells in the bounds
+	 * raster map, I would suggest to replace NULL with min(bounds) - 1 or
+	 +* max(bounds) + 1. */
     if (files->bounds_map != NULL) {
 	if (segment_open
 	    (&files->bounds_seg, G_tempfile(), files->nrows, files->ncols,

Modified: grass-addons/grass7/imagery/i.segment/parse_args.c
===================================================================
--- grass-addons/grass7/imagery/i.segment/parse_args.c	2012-06-26 03:22:09 UTC (rev 52216)
+++ grass-addons/grass7/imagery/i.segment/parse_args.c	2012-06-26 03:33:25 UTC (rev 52217)
@@ -13,7 +13,8 @@
     /* reference: http://grass.osgeo.org/programming7/gislib.html#Command_Line_Parsing */
 
     struct Option *group, *seeds, *bounds, *output, *method, *threshold;	/* Establish an Option pointer for each option */
-    struct Flag *diagonal, *weighted;	/* Establish a Flag pointer for each option */
+    struct Flag *diagonal, *weighted;										/* Establish a Flag pointer for each option */
+	struct Option *outband, *endt; 		/* debugging parameters... TODO: leave in code or remove?  hidden options? */
 
     /* required parameters */
     group = G_define_standard_option(G_OPT_I_GROUP);	/* TODO: OK to require the user to create a group?  Otherwise later add an either/or option to give just a single raster map... */
@@ -31,7 +32,7 @@
     method->type = TYPE_STRING;
     method->required = YES;
     method->answer = "region_growing";
-    method->options = "region_growing, io_debug, ll_test";	/*io_debug just writes row+col to each output pixel, ll_test for testing linked list data structure */
+    method->options = "region_growing, io_debug, ll_test";	/* TODO at end, remove these from list: io_debug just writes row+col to each output pixel, ll_test for testing linked list data structure */
     method->description = _("Segmentation method.");
 
     /* optional parameters */
@@ -62,41 +63,37 @@
 	_("Optional bounding/constraining raster map, must be integer values, each area will be segmented independent of the others.");
     /*    bounds->description = _("Optional vector map with polygons to bound (constrain) segmentation."); */
     /* TODO: if performing second segmentation, will already have raster map from this module
-     * If have preexisting boundaries (landuse, etc) will have vector map
+     * If have preexisting boundaries (landuse, etc) will have vector map?
      * Seems we need to have it in raster format for processing, is it OK to have the user run v.to.rast before doing the segmentation?
      * Or for hierarchical segmentation, will it be easier to have the polygons?
-     * ....will start with rasters, quickest to implement.... */
+     *  */
 
     /* TODO input for distance function */
 
+	/* debug parameters */
+    endt = G_define_option();
+    endt->key = "endt";
+    endt->type = TYPE_INTEGER;
+    endt->required = YES; /* TODO, could put as optional, and if not supplied put in something large. */
+    endt->description = _("Maximum number of time steps to complete.");
+
+    outband = G_define_standard_option(G_OPT_R_INPUT);
+    outband->key = "final_mean";
+    outband->type = TYPE_STRING;
+    outband->required = YES;
+    outband->description =
+	_("debug - save band mean, currently implemented for only 1 band.");
+
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
-
-    /* G_debug(1, "For the option <%s> you chose: <%s>",
-       group->description, group->answer);
-
-       G_debug(1, "For the option <%s> you chose: <%s>",
-       seeds->description, seeds->answer);
-
-       G_debug(1, "For the option <%s> you chose: <%s>",
-       output->description, output->answer);
-
-       G_debug(1, "For the option <%s> you chose: <%s>",
-       method->description, method->answer);
-
-       G_debug(1, "For the option <%s> you chose: <%s>",
-       threshold->description, threshold->answer);
-
-       G_debug(1, "The value of the diagonal flag is: %d", diagonal->answer); */
-
     /* Validation */
 
     /* TODO: use checker for any of the data validation steps!? */
 
     /* ToDo The most important things to check are if the
        input and output raster maps can be opened (non-negative file
-       descriptor).  Do this here or in open_files?  */
+       descriptor).  ??? Do this here or in open_files?  */
 
 
     /* Check and save parameters */
@@ -112,7 +109,7 @@
     /* reference r.cost line 313 
        if (sscanf(opt5->answer, "%d", &maxcost) != 1 || maxcost < 0)
        G_fatal_error(_("Inappropriate maximum cost: %d"), maxcost); */
-    sscanf(threshold->answer, "%f", &functions->threshold);	/* Note: this gets changed after we know more at beginning of create_isegs() */
+    sscanf(threshold->answer, "%f", &functions->threshold);	/* Note: this threshold is scaled after we know more at the beginning of create_isegs() */
 
     /* segmentation methods:  0 = debug, 1 = region growing */
     /* TODO, instead of string compare, does the Option structure have these already numbered? */
@@ -145,23 +142,30 @@
 
     files->seeds = seeds->answer;
 
-    /*todo, check error trapping here, what if nothing is entered?  what if nothing is found? what if in same mapset */
     if (bounds->answer == NULL) {	/*no polygon constraints */
 	files->bounds_map = NULL;
     }
     else {			/* polygon constraints given */
-
-	files->bounds_map = bounds->answer;	/*todo, this needs to set files->bounds = NULL if no answer was given to the parameter */
+	files->bounds_map = bounds->answer;	
 	if ((files->bounds_mapset = G_find_raster(files->bounds_map, "")) == NULL) {	/* TODO, warning here:  parse_args.c:149:27: warning: assignment discards ‘const’ qualifier from pointer target type [enabled by default] */
 	    G_fatal_error(_("Segmentation constraint/boundary map not found."));
 	}
+	if (Rast_map_type(files->bounds_map, files->bounds_mapset) != CELL_TYPE) {
+		G_fatal_error(_("Segmentation constraint map must be CELL type (integers)"));
+	}
     }
-    /* todo, check raster type, needs to be CELL (integer) */
-    /*todo print out what these were before and after */
 
     /* other data */
     files->nrows = Rast_window_rows();
     files->ncols = Rast_window_cols();
 
+	/* debug help */
+    if (G_legal_filename(outband->answer) == 1)
+	files->out_band = outband->answer;	/* name of current means */
+    else
+	G_fatal_error("Invalid output raster name for means.");
+
+    sscanf(endt->answer, "%d", &functions->end_t);	/* Note: this threshold is scaled after we know more at the beginning of create_isegs() */
+
     return 0;
 }

Added: grass-addons/grass7/imagery/i.segment/times.txt
===================================================================
--- grass-addons/grass7/imagery/i.segment/times.txt	                        (rev 0)
+++ grass-addons/grass7/imagery/i.segment/times.txt	2012-06-26 03:33:25 UTC (rev 52217)
@@ -0,0 +1,41 @@
+Using version from around June 20, SVN 52181, it does not remember fragment membership.
+
+50 x 50 = 2500 pixels
+
+real	0m5.767s
+user	0m5.704s
+sys	0m0.028s
+
+ 0.00228 s/pixel
+
+
+100 x 100 = 10,000 pixels 
+
+real	4m8.335s
+user	4m7.363s
+sys	0m0.048s
+
+0.0248 s/pixel
+
+
+140 x 140 = 19,600 pixels
+
+real	150m7.509s
+user	149m27.656s
+sys	0m9.025s
+
+0.45 s/pixel
+
+
+Rough target:
+
+10,000 x 10,000 = 100,000,000 pixels
+
+10 hours?
+
+
+removing squareroot
+17s --> 4s at 50x50 resolution
+23 --> 17s at 70x70
+
+

Modified: grass-addons/grass7/imagery/i.segment/write_output.c
===================================================================
--- grass-addons/grass7/imagery/i.segment/write_output.c	2012-06-26 03:22:09 UTC (rev 52216)
+++ grass-addons/grass7/imagery/i.segment/write_output.c	2012-06-26 03:33:25 UTC (rev 52217)
@@ -9,36 +9,62 @@
 
 int write_output(struct files *files)
 {
-    int out_fd, row, col;
+    int out_fd, mean_fd, row, col; /* mean_fd for validiating/debug of means, todo, could add array... */
     CELL *outbuf;
+    DCELL *meanbuf;
 
     outbuf = Rast_allocate_c_buf();	/*hold one row of data to put into raster */
+	meanbuf = Rast_allocate_d_buf();
 
     /* Todo: return codes are 1 for these, need to check and react to errors? programmer's manual didn't include it... */
 
     segment_flush(&files->out_seg);	/* force all data to disk */
+    segment_flush(&files->bands_seg); /* TODO use IFDEF or just delete for all these parts? for debug/validation output */
 
     G_debug(1, "preparing output raster");
     /* open output raster map */
     out_fd = Rast_open_new(files->out_name, CELL_TYPE);
-
+	mean_fd = Rast_open_new(files->out_band, DCELL_TYPE);
+	
     G_debug(1, "start data transfer from segmentation file to raster");
     /* transfer data from segmentation file to raster */
     /* output segmentation file: each element includes the segment ID then processing flag(s).  So just need the first part of it. */
 
     for (row = 0; row < files->nrows; row++) {
-	G_percent(row, files->nrows, 1);	/*TODO: Why isn't this getting printed? */
+	G_percent(row, files->nrows, 1);
+	Rast_set_c_null_value(outbuf, files->ncols); /*set buffer to NULLs, only write those that weren't originally masked */
+	Rast_set_d_null_value(meanbuf, files->ncols);
 	for (col = 0; col < files->ncols; col++) {
 	    segment_get(&files->out_seg, (void *)files->out_val, row, col);
 	    G_debug(5, "outval[0] = %i", files->out_val[0]);
-	    outbuf[col] = files->out_val[0];	/*just want segment assignment, not the other processing flag(s) */
+	    if(files->out_val[0] >= 0) /* only write positive segment ID's, using -1 as indicator of Null/Masked pixels.  TODO: OK to use -1 as flag for this? */
+	    {outbuf[col] = files->out_val[0];	/*just want segment assignment, not the other processing flag(s) */
+	//	meanbuf[col] = files->out_val[0];
 	}
+		
+		segment_get(&files->bands_seg, (void *)files->bands_val, row, col);
+	    if(files->out_val[0] >= 0) 
+	    meanbuf[col] = files->bands_val[0];
+	    
+	}
 	Rast_put_row(out_fd, outbuf, CELL_TYPE);
+	Rast_put_row(mean_fd, meanbuf, DCELL_TYPE);
     }
 
+/* TODO: I don't understand the header/history/etc for raster storage.  Can/should we save any information about the segmentation
+ * settings used to create the raster?  What the input image group was?  Anything else the statistics module would need?
+ */
+
+/* TODO after we have a count of create segments... if(total segments == 0) G_warning(_("No segments were created. Verify threshold and region settings.")); */
+
     /* close and save file */
     Rast_close(out_fd);
+    Rast_close(mean_fd);
 
+	/* free memory */
+	G_free(outbuf);
+	G_free(meanbuf);
+
     return 0;
 }
 



More information about the grass-commit mailing list