[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(®ion, 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, ®ion, ll_pym, nclasses, alpha_dec);
+ seq_MAP(sf_pym, ®ion, ll_pym, nclasses, alpha_dec, goodness);
}
} while (increment_reg(®ion, 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