[GRASS-SVN] r54321 - grass/trunk/imagery/i.smap

svn_grass at osgeo.org svn_grass at osgeo.org
Mon Dec 17 05:33:42 PST 2012


Author: mmetz
Date: 2012-12-17 05:33:42 -0800 (Mon, 17 Dec 2012)
New Revision: 54321

Modified:
   grass/trunk/imagery/i.smap/bouman.h
   grass/trunk/imagery/i.smap/closefiles.c
   grass/trunk/imagery/i.smap/i.smap.html
   grass/trunk/imagery/i.smap/interp.c
   grass/trunk/imagery/i.smap/openfiles.c
   grass/trunk/imagery/i.smap/parse.c
   grass/trunk/imagery/i.smap/region.h
   grass/trunk/imagery/i.smap/segment.c
   grass/trunk/imagery/i.smap/write_img.c
Log:
i.smap: new goodness of fit output map

Modified: grass/trunk/imagery/i.smap/bouman.h
===================================================================
--- grass/trunk/imagery/i.smap/bouman.h	2012-12-17 00:46:37 UTC (rev 54320)
+++ grass/trunk/imagery/i.smap/bouman.h	2012-12-17 13:33:42 UTC (rev 54321)
@@ -7,6 +7,7 @@
 struct files
 {
     int output_fd;
+    int goodness_fd;
     struct Categories output_labels;
 
     int *band_fd;
@@ -20,6 +21,7 @@
 struct parms
 {
     char *output_map;
+    char *goodness_map;
     char *group;
     char *subgroup;
     char *sigfile;
@@ -66,7 +68,7 @@
 int create_output_labels(struct SigSet *, struct files *);
 
 /* write_img.c */
-int write_img(unsigned char **, int, int, struct SigSet *, struct parms *,
+int write_img(unsigned char **, float **, int, int, struct SigSet *, struct parms *,
 	      struct files *);
 #endif
 

Modified: grass/trunk/imagery/i.smap/closefiles.c
===================================================================
--- grass/trunk/imagery/i.smap/closefiles.c	2012-12-17 00:46:37 UTC (rev 54320)
+++ grass/trunk/imagery/i.smap/closefiles.c	2012-12-17 13:33:42 UTC (rev 54321)
@@ -21,5 +21,11 @@
     make_history(parms->output_map,
 		 parms->group, parms->subgroup, parms->sigfile);
 
+    if (files->goodness_fd >= 0) {
+	Rast_close(files->goodness_fd);
+	make_history(parms->goodness_map,
+		     parms->group, parms->subgroup, parms->sigfile);
+    }
+
     return 0;
 }

Modified: grass/trunk/imagery/i.smap/i.smap.html
===================================================================
--- grass/trunk/imagery/i.smap/i.smap.html	2012-12-17 00:46:37 UTC (rev 54320)
+++ grass/trunk/imagery/i.smap/i.smap.html	2012-12-17 13:33:42 UTC (rev 54321)
@@ -131,6 +131,11 @@
 excessively large regions are not formed.
 
 <p>
+The degree of misclassifications can be investigated with the goodness 
+of fit output map. Lower values indicate a better fit. The largest 5 to 
+15% of the goodness values may need some closer inspection.
+
+<p>
 The module <em>i.smap</em> does not support MASKed or NULL cells. Therefore 
 it might be necessary to create a copy of the classification results 
 using e.g. r.mapcalc:

Modified: grass/trunk/imagery/i.smap/interp.c
===================================================================
--- grass/trunk/imagery/i.smap/interp.c	2012-12-17 00:46:37 UTC (rev 54320)
+++ grass/trunk/imagery/i.smap/interp.c	2012-12-17 13:33:42 UTC (rev 54321)
@@ -9,12 +9,12 @@
 #define ML_PRECISION 1e-6
 
 static void seq_MAP_routine(unsigned char ***, struct Region *,
-			    LIKELIHOOD ****, int, double *);
+			    LIKELIHOOD ****, int, double *, float **);
 static double alpha_dec_max(double ***);
 static void print_N(double ***);
 static void print_alpha(double *);
 static void interp(unsigned char **, struct Region *, unsigned char **,
-		   LIKELIHOOD ***, int, double *, int, double ***, int);
+		   LIKELIHOOD ***, int, double *, int, double ***, int, float **);
 void MLE(unsigned char **, LIKELIHOOD ***, struct Region *, int);
 static int up_char(int, int, struct Region *, unsigned char **,
 		   unsigned char **);
@@ -23,8 +23,9 @@
 void seq_MAP(unsigned char ***sf_pym,	/* pyramid of segmentations */
 	     struct Region *region,	/* specifies image subregion */
 	     LIKELIHOOD **** ll_pym,	/* pyramid of class statistics */
-	     int M,		/* number of classes */
-	     double *alpha_dec	/* decimation parameters returned by seq_MAP */
+	     int M,		        /* number of classes */
+	     double *alpha_dec,	        /* decimation parameters returned by seq_MAP */
+	     float **goodness          /* goodness of fit */
     )
 {
     int repeat;
@@ -36,15 +37,16 @@
 	G_debug(1, "Pyramid constructed");
 
 	/* Perform sequential MAP segmentation using EM algorithm */
-	seq_MAP_routine(sf_pym, region, ll_pym, M, alpha_dec);
+	seq_MAP_routine(sf_pym, region, ll_pym, M, alpha_dec, goodness);
     }
 }
 
 static void seq_MAP_routine(unsigned char ***sf_pym,	/* pyramid of segmentations */
 			    struct Region *region,	/* specifies image subregion */
 			    LIKELIHOOD **** ll_pym,	/* pyramid of class statistics */
-			    int M,	/* number of classes */
-			    double *alpha_dec	/* decimation parameters returned by seq_MAP */
+			    int M,	                /* number of classes */
+			    double *alpha_dec,	        /* decimation parameters returned by seq_MAP */
+			    float **goodness           /* goodness of fit */
     )
 {
     int j, k;			/* loop index */
@@ -105,7 +107,7 @@
 	 * fixed number of iterations or until convergence.   */
 	do {
 	    interp(sf_pym[D], &(regionary[D]), sf_pym[D + 1], ll_pym[D], M,
-		   alpha, period[D], N, 1);
+		   alpha, period[D], N, 1, NULL);
 	    print_N(N);
 	    G_debug(4, "log likelihood = %f", log_like(N, alpha, M));
 	    for (j = 0; j < 3; j++)
@@ -119,8 +121,13 @@
 		diff1 += fabs(tmp[j] - alpha[j]);
 	    diff2 = log_like(N, alpha, M) - log_like(N, tmp, M);
 	} while ((diff1 > EM_PRECISION) && (diff2 > 0));
-	interp(sf_pym[D], &(regionary[D]), sf_pym[D + 1], ll_pym[D], M, alpha,
-	       1, N, 0);
+	/* get goodness of fit for D == 0 */
+	if (D == 0)
+	    interp(sf_pym[D], &(regionary[D]), sf_pym[D + 1], ll_pym[D], M, alpha,
+		   1, N, 0, goodness);
+	else
+	    interp(sf_pym[D], &(regionary[D]), sf_pym[D + 1], ll_pym[D], M, alpha,
+		   1, N, 0, NULL);
 	alpha_dec[D] = alpha_dec_max(N);
 
 	print_N(N);
@@ -193,7 +200,8 @@
 		      double *alpha,	/* transition probability parameters; alpha[3] */
 		      int period,	/* sampling period of interpolation */
 		      double ***N,	/* transition probability statistics; N[2][3][2] */
-		      int statflag	/* compute transition statistics if == 1 */
+		      int statflag,	/* compute transition statistics if == 1 */
+		      float **goodness /* cost of best class */
     )
 {
     int i, j;			/* pixel index */
@@ -255,6 +263,9 @@
 		}
 	    }
 	    sf1[i][j] = best;
+	    /* save cost as best fit indicator */
+	    if (goodness)
+		goodness[i][j] = mincost;
 
 	    /* if not on boundary, compute expectation of N */
 	    if ((!bflag) && (statflag)) {

Modified: grass/trunk/imagery/i.smap/openfiles.c
===================================================================
--- grass/trunk/imagery/i.smap/openfiles.c	2012-12-17 00:46:37 UTC (rev 54320)
+++ grass/trunk/imagery/i.smap/openfiles.c	2012-12-17 13:33:42 UTC (rev 54321)
@@ -39,5 +39,10 @@
     /* open output map */
     files->output_fd = open_cell_new(parms->output_map);
 
+    if (parms->goodness_map)
+	files->goodness_fd = Rast_open_new(parms->goodness_map, FCELL_TYPE);
+    else
+	files->goodness_fd = -1;
+
     return 0;
 }

Modified: grass/trunk/imagery/i.smap/parse.c
===================================================================
--- grass/trunk/imagery/i.smap/parse.c	2012-12-17 00:46:37 UTC (rev 54320)
+++ grass/trunk/imagery/i.smap/parse.c	2012-12-17 13:33:42 UTC (rev 54321)
@@ -7,7 +7,7 @@
 
 int parse(int argc, char *argv[], struct parms *parms)
 {
-    struct Option *group, *subgroup, *sigfile, *output;
+    struct Option *group, *subgroup, *sigfile, *output, *goodness;
     struct Option *blocksize;
     struct Flag *ml;
 
@@ -25,6 +25,11 @@
 
     output = G_define_standard_option(G_OPT_R_OUTPUT);
 
+    goodness = G_define_standard_option(G_OPT_R_OUTPUT);
+    goodness->key = "goodness";
+    goodness->required = NO;
+    goodness->label = _("Name of goodness of fit map (lower is better)");
+
     blocksize = G_define_option();
     blocksize->key = "blocksize";
     blocksize->description = _("Size of submatrix to process at one time");
@@ -46,6 +51,8 @@
     parms->group = group->answer;
     parms->subgroup = subgroup->answer;
     parms->sigfile = sigfile->answer;
+    
+    parms->goodness_map = goodness->answer;
 
     /* check all the inputs */
     if (!I_find_group(parms->group))

Modified: grass/trunk/imagery/i.smap/region.h
===================================================================
--- grass/trunk/imagery/i.smap/region.h	2012-12-17 00:46:37 UTC (rev 54320)
+++ grass/trunk/imagery/i.smap/region.h	2012-12-17 13:33:42 UTC (rev 54321)
@@ -21,7 +21,7 @@
 
 /* interp.c */
 void seq_MAP(unsigned char ***, struct Region *, LIKELIHOOD ****, int,
-	     double *);
+	     double *, float **);
 void MLE(unsigned char **, LIKELIHOOD ***, struct Region *, int);
 
 /* reg_util.c */

Modified: grass/trunk/imagery/i.smap/segment.c
===================================================================
--- grass/trunk/imagery/i.smap/segment.c	2012-12-17 00:46:37 UTC (rev 54320)
+++ grass/trunk/imagery/i.smap/segment.c	2012-12-17 13:33:42 UTC (rev 54321)
@@ -41,6 +41,7 @@
     LIKELIHOOD ****ll_pym;	/* pyramid of log likelihoods */
     unsigned char ***sf_pym;	/* pyramid of segmentations */
     int D;			/* number of levels in pyramid */
+    float **goodness;          /* goodness of fit */
     double *alpha_dec;		/* class transition probabilities */
     int i;
 
@@ -82,6 +83,16 @@
     /* allocate memory for segmentation pyramid */
     sf_pym = (unsigned char ***)get_pyramid(wd, ht, sizeof(char));
 
+    if (parms->goodness_map) {
+	goodness = G_malloc(ht * sizeof(float *));
+	goodness[0] = G_malloc(sizeof(float *) * ht * wd);
+	for (i = 1; i < ht; i++) {
+	    goodness[i] = goodness[i - 1] + wd;
+	}
+    }
+    else
+	goodness = NULL;
+
     /* tiled segmentation */
     init_reg(&region, wd, ht, block_size);
     extract_init(S);
@@ -103,13 +114,13 @@
 	else {
 	    for (i = 0; i < D; i++)
 		alpha_dec[i] = 1.0;
-	    seq_MAP(sf_pym, &region, ll_pym, nclasses, alpha_dec);
+	    seq_MAP(sf_pym, &region, ll_pym, nclasses, alpha_dec, goodness);
 	}
 
 
     } while (increment_reg(&region, wd, ht, block_size));
 
-    write_img(sf_pym[0], wd, ht, S, parms, files);
+    write_img(sf_pym[0], goodness, wd, ht, S, parms, files);
 
     return 0;
 }

Modified: grass/trunk/imagery/i.smap/write_img.c
===================================================================
--- grass/trunk/imagery/i.smap/write_img.c	2012-12-17 00:46:37 UTC (rev 54320)
+++ grass/trunk/imagery/i.smap/write_img.c	2012-12-17 13:33:42 UTC (rev 54321)
@@ -4,14 +4,20 @@
 
 #include "bouman.h"
 
-int write_img(unsigned char **img, int ncols, int nrows, struct SigSet *S,	/* class parameters */
+int write_img(unsigned char **img, float **goodness, int ncols, int nrows,
+              struct SigSet *S,	/* class parameters */
 	      struct parms *parms,	/* parms: command line parameters */
 	      struct files *files)
 {				/* files: contains file to output */
     int row, col;
+    FCELL *fcellbuf = NULL;
 
     G_message(_("Writing raster map <%s>..."), parms->output_map);
 
+    /* write goodness of fit */
+    if (parms->goodness_map)
+	fcellbuf = Rast_allocate_f_buf();
+
     for (row = 0; row < nrows; row++) {
 	G_percent(row, nrows, 2);
 	for (col = 0; col < ncols; col++) {
@@ -19,8 +25,13 @@
 		
 	    G_debug(3, "class: [%d] row/col: [%d][%d]", class, row, col);
 	    files->outbuf[col] = (CELL) S->ClassSig[class].classnum;
+	    
+	    if (parms->goodness_map)
+		fcellbuf[col] = goodness[row][col];
 	}
 	Rast_put_row(files->output_fd, files->outbuf, CELL_TYPE);
+	if (parms->goodness_map)
+	    Rast_put_row(files->goodness_fd, fcellbuf, FCELL_TYPE);
     }
     G_percent(nrows, nrows, 2);
 



More information about the grass-commit mailing list