[GRASS-SVN] r34871 - grass/trunk/imagery/i.smap
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Dec 14 09:29:10 EST 2008
Author: martinl
Date: 2008-12-14 09:29:09 -0500 (Sun, 14 Dec 2008)
New Revision: 34871
Added:
grass/trunk/imagery/i.smap/alpha_max.c
grass/trunk/imagery/i.smap/bouman.h
grass/trunk/imagery/i.smap/closefiles.c
grass/trunk/imagery/i.smap/decimate.c
grass/trunk/imagery/i.smap/eigen.c
grass/trunk/imagery/i.smap/history.c
grass/trunk/imagery/i.smap/i.smap.html
grass/trunk/imagery/i.smap/interp.c
grass/trunk/imagery/i.smap/invert.c
grass/trunk/imagery/i.smap/labels.c
grass/trunk/imagery/i.smap/local_proto.h
grass/trunk/imagery/i.smap/main.c
grass/trunk/imagery/i.smap/model.c
grass/trunk/imagery/i.smap/multialloc.c
grass/trunk/imagery/i.smap/opencell.c
grass/trunk/imagery/i.smap/openfiles.c
grass/trunk/imagery/i.smap/parse.c
grass/trunk/imagery/i.smap/read_block.c
grass/trunk/imagery/i.smap/read_sig.c
grass/trunk/imagery/i.smap/reg_util.c
grass/trunk/imagery/i.smap/region.h
grass/trunk/imagery/i.smap/segment.c
grass/trunk/imagery/i.smap/solve.c
grass/trunk/imagery/i.smap/write_img.c
Removed:
grass/trunk/imagery/i.smap/bouman/
grass/trunk/imagery/i.smap/shapiro/
Modified:
grass/trunk/imagery/i.smap/Makefile
Log:
i.smap: remove redundant 'bouman' and 'shapiro' directories
Modified: grass/trunk/imagery/i.smap/Makefile
===================================================================
--- grass/trunk/imagery/i.smap/Makefile 2008-12-14 14:20:47 UTC (rev 34870)
+++ grass/trunk/imagery/i.smap/Makefile 2008-12-14 14:29:09 UTC (rev 34871)
@@ -1,14 +1,14 @@
-
MODULE_TOPDIR = ../..
-SUBDIRS = \
- bouman \
- shapiro
+PGM = i.smap
-include $(MODULE_TOPDIR)/include/Make/Dir.make
+LIBES = $(IMAGERYLIB) $(GISLIB) $(GMATHLIB)
+DEPENDENCIES= $(IMAGERYDEP) $(GISDEP) $(GMATHDEP)
-default: parsubdirs
+include $(MODULE_TOPDIR)/include/Make/Module.make
+default: cmd
-shapiro: bouman
+
+
Copied: grass/trunk/imagery/i.smap/alpha_max.c (from rev 34868, grass/trunk/imagery/i.smap/bouman/alpha_max.c)
===================================================================
--- grass/trunk/imagery/i.smap/alpha_max.c (rev 0)
+++ grass/trunk/imagery/i.smap/alpha_max.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,161 @@
+#include <stdio.h>
+#include <math.h>
+#include "bouman.h"
+
+/* Global Variables: These variables are required as inputs *
+ * to the function called by ``solve''. Since the function is *
+ * called by a pointer, side inputs must be passed as global *
+ * variables root finding */
+static double ***N_glb; /* N_glb[2][3][2] rate statistics */
+static double *b_glb; /* line search direction */
+static double eps_glb; /* precision */
+static int M_glb; /* number of classes */
+
+
+
+void alpha_max(double ***N, /* Transition probability statistics; N2[2][3][2] */
+ double *a, /* Transition probability parameters; a[3] */
+ int M, /* Number of class types */
+ double eps /* Precision of convergence */
+ )
+{
+ double b[3];
+
+ /* b[0] = 0.4;
+ b[1] = 0.6/2;
+ b[2] = 0.0; */
+
+ b[0] = 3.0;
+ b[1] = 2.0;
+ b[2] = 0.0;
+ line_search(N, a, M, b, eps);
+}
+
+
+void line_search(
+ /* line_search determines the maximum value of the log likelihood *
+ * along the direction b, subject to the constraint a[0]+2*a[1]+a[3]<1. */
+ double ***N, /* Transition probability statistics; N2[2][3][2] */
+ double *a, /* Transition probability parameters; a[3] */
+ int M, /* Number of class types */
+ double *b, /* direction of search */
+ double eps /* Precision of convergence */
+ )
+{
+ int code; /* error code for solve subroutine */
+ double x; /* distance along line */
+ double max; /* maximum value for x */
+
+ normalize(b);
+
+ a[0] = eps * b[0];
+ a[1] = eps * b[1];
+ a[2] = eps * b[2];
+
+ /* enforce condition [1,2,1][a[0],a[1],a[2]]^t <1-eps */
+ max = (1 - eps) / (b[0] + 2 * b[1] + b[2]);
+
+ /* set global variables for solve routine */
+ N_glb = N;
+ b_glb = b;
+ eps_glb = eps;
+ M_glb = M;
+
+ /* minimize on line. Avoid singular boundary */
+ x = solve(func, eps, max, eps, &code);
+
+ /* If derivative was positive on line, x=max. */
+ if (code == 1)
+ x = max;
+
+ /* If derivative was negative on line. */
+ if (code == -1)
+ x = 0.0;
+
+ /* compute a */
+ a[0] = x * b[0];
+ a[1] = x * b[1];
+ a[2] = x * b[2];
+}
+
+
+int normalize(
+ /* normalize the vector b[3]. Return 0 if null vector */
+ double b[3]
+ )
+{
+ double norm;
+
+ norm = sqrt(b[0] * b[0] + b[1] * b[1] + b[2] * b[2]);
+ if (norm == 0)
+ return (0);
+ b[0] = b[0] / norm;
+ b[1] = b[1] / norm;
+ b[2] = b[2] / norm;
+ return (1);
+}
+
+double func(double x)
+{
+ double tmp[3], grad[3];
+
+ tmp[0] = x * b_glb[0];
+ tmp[1] = x * b_glb[1];
+ tmp[2] = x * b_glb[2];
+
+ gradient(grad, N_glb, tmp, M_glb);
+ return (b_glb[0] * grad[0] + b_glb[1] * grad[1] + b_glb[2] * grad[2]);
+}
+
+double log_like(
+ /* compute log likelihood being maximized. This subroutine *
+ * is useful for debuging since the log likelihood must be *
+ * monotonically decreasing. */
+ double ***N, /* transition statistics */
+ double a[3], /* transition parameters */
+ int M /* number of classes */
+ )
+{
+ double iM, tmp, sum;
+ int n1, n2, n3;
+
+ iM = 1.0 / M;
+ sum = 0;
+ for (n1 = 0; n1 <= 1; n1++)
+ for (n2 = 0; n2 <= 2; n2++)
+ for (n3 = 0; n3 <= 1; n3++) {
+ tmp =
+ log(a[0] * (n1 - iM) + a[1] * (n2 - 2.0 * iM) +
+ a[2] * (n3 - iM) + iM);
+ sum += N[n1][n2][n3] * tmp;
+ }
+ return (sum);
+}
+
+
+void gradient(
+ /* computes the gradient of the log likelihood being maximized. */
+ double grad[3], /* gradient vector; output */
+ double ***N, /* transition statistics; input */
+ double a[3], /* transition parameters; input */
+ int M /* number of classes */
+ )
+{
+ double iM, tmp;
+ int n1, n2, n3;
+
+ iM = 1.0 / M;
+ grad[0] = grad[1] = grad[2] = 0;
+ for (n1 = 0; n1 <= 1; n1++)
+ for (n2 = 0; n2 <= 2; n2++)
+ for (n3 = 0; n3 <= 1; n3++) {
+ tmp =
+ a[0] * (n1 - iM) + a[1] * (n2 - 2.0 * iM) + a[2] * (n3 -
+ iM) +
+ iM;
+ tmp = 1 / tmp;
+ grad[0] += tmp * (n1 - iM) * N[n1][n2][n3];
+ grad[1] += tmp * (n2 - 2.0 * iM) * N[n1][n2][n3];
+ grad[2] += tmp * (n3 - iM) * N[n1][n2][n3];
+ }
+}
Property changes on: grass/trunk/imagery/i.smap/alpha_max.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/bouman.h (from rev 34868, grass/trunk/imagery/i.smap/bouman/bouman.h)
===================================================================
--- grass/trunk/imagery/i.smap/bouman.h (rev 0)
+++ grass/trunk/imagery/i.smap/bouman.h 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,77 @@
+#include <math.h>
+#include <string.h>
+#include <grass/gis.h>
+
+#define LIKELIHOOD float
+
+struct files
+{
+ int output_fd;
+ struct Categories output_labels;
+
+ int *band_fd;
+ int nbands;
+
+ DCELL *cellbuf;
+ CELL *outbuf;
+ char *isdata;
+};
+
+struct parms
+{
+ char *output_map;
+ char *group;
+ char *subgroup;
+ char *sigfile;
+ int blocksize;
+ int quiet;
+ int ml;
+};
+
+/* parse.c */
+int parse(int, char *[], struct parms *);
+
+/* closefiles.c */
+int closefiles(struct parms *, struct files *);
+
+/* openfiles.c */
+int openfiles(struct parms *, struct files *);
+
+/* Suboutines in alpha_max.c */
+void alpha_max(double ***, double *, int, double);
+void line_search(double ***, double *, int, double *, double);
+int normalize(double[3]);
+double func(double);
+double log_like(double ***, double[3], int);
+void gradient(double[3], double ***, double[3], int);
+
+/* Subroutines in multialloc.c */
+char *multialloc(size_t, int, ...);
+void multifree(char *, int);
+unsigned char **get_img(int, int, size_t);
+void free_img(unsigned char **);
+
+/* Subroutine in solve.c */
+double solve(double (*)(double), double, double, double, int *);
+
+/* Subroutine in eigen.c */
+int eigen(double **, double *, int);
+
+/* Subroutine in invert.c */
+int invert(double **, int);
+
+#ifdef GRASS_IMAGERY_H
+int segment(struct SigSet *, struct parms *, struct files *);
+
+/* read_sig.c */
+int read_signatures(struct parms *, struct SigSet *);
+
+/* labels.c */
+int create_output_labels(struct SigSet *, struct files *);
+
+/* write_img.c */
+int write_img(unsigned char **, int, int, struct SigSet *, struct parms *,
+ struct files *);
+#endif
+
+/* Look for prototypes that use the Region structure in region.h */
Property changes on: grass/trunk/imagery/i.smap/bouman.h
___________________________________________________________________
Name: svn:mime-type
+ text/x-chdr
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/closefiles.c (from rev 34870, grass/trunk/imagery/i.smap/shapiro/closefiles.c)
===================================================================
--- grass/trunk/imagery/i.smap/closefiles.c (rev 0)
+++ grass/trunk/imagery/i.smap/closefiles.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,24 @@
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/imagery.h>
+#include "bouman.h"
+#include "local_proto.h"
+
+
+int closefiles(struct parms *parms, struct files *files)
+{
+ int n;
+
+ G_debug(1, "Creating support files for <%s>...", parms->output_map);
+
+ for (n = 0; n < files->nbands; n++)
+ G_close_cell(files->band_fd[n]);
+
+ G_close_cell(files->output_fd);
+ G_write_cats(parms->output_map, &files->output_labels);
+ make_history(parms->output_map,
+ parms->group, parms->subgroup, parms->sigfile);
+
+ return 0;
+}
Property changes on: grass/trunk/imagery/i.smap/closefiles.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/decimate.c (from rev 34870, grass/trunk/imagery/i.smap/bouman/decimate.c)
===================================================================
--- grass/trunk/imagery/i.smap/decimate.c (rev 0)
+++ grass/trunk/imagery/i.smap/decimate.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,229 @@
+#include <grass/gis.h>
+#include "bouman.h"
+#include "region.h"
+
+
+static void decimate(LIKELIHOOD ***, struct Region *, int, LIKELIHOOD ***,
+ double);
+static void up_ll(LIKELIHOOD *, int, double, LIKELIHOOD *);
+
+
+void make_pyramid(LIKELIHOOD **** ll_pym, /* log likelihood pyramid, ll_pym[scale][i][j][class] */
+ struct Region *region, /* specifies image subregion */
+ int M, /* number of classes */
+ double *alpha, /* decimation parameters */
+ int vlevel /* verbose operation */
+ )
+{
+ int D;
+ int wd, ht;
+ struct Region region_buff;
+
+ /* save region stucture */
+ copy_reg(region, ®ion_buff);
+
+ D = 0;
+ reg_to_wdht(region, &wd, &ht);
+ while ((wd > 2) && (ht > 2)) {
+ if (vlevel >= 2)
+ G_debug(1, "D = %d alpha = %f; 1-alpha = %f", D, alpha[D],
+ 1 - alpha[D]);
+ decimate(ll_pym[D], region, M, ll_pym[D + 1], alpha[D]);
+ dec_reg(region);
+ reg_to_wdht(region, &wd, &ht);
+ D++;
+ }
+
+ /* restore region structure */
+ copy_reg(®ion_buff, region);
+}
+
+static void decimate(
+ /* decimate statistics ll1 to form ll2 */
+ LIKELIHOOD *** ll1, /* log likelihood ll1[i][j][class], at fine resolution */
+ struct Region *region1, /* specifies image subregion */
+ int M, /* number of classes */
+ LIKELIHOOD *** ll2, /* loglikelihood ll1[i][j][class], at coarse resolution */
+ double alpha /* transition parameter */
+ )
+{
+ struct Region *region2, reg_spc; /* coarse resolution region */
+ int wflag, hflag; /* flags indicate odd number of pixels */
+ int i, j, m;
+ LIKELIHOOD *node; /* coarse resolution point */
+ LIKELIHOOD *pt1, *pt2, *pt3, *pt4; /* fine resolution neighbors */
+
+ region2 = ®_spc;
+
+ copy_reg(region1, region2);
+ dec_reg(region2);
+
+ wflag = region1->xmax & 1;
+ hflag = region1->ymax & 1;
+
+ for (i = region2->ymin; i < region2->ymax; i++)
+ for (j = region2->xmin; j < region2->xmax; j++) {
+ pt1 = ll1[2 * i][2 * j];
+ pt2 = ll1[2 * i][2 * j + 1];
+ pt3 = ll1[2 * i + 1][2 * j];
+ pt4 = ll1[2 * i + 1][2 * j + 1];
+
+ node = ll2[i][j];
+ for (m = 0; m < M; m++)
+ node[m] = 0.0;
+
+ up_ll(pt1, M, alpha, node);
+ up_ll(pt2, M, alpha, node);
+ up_ll(pt3, M, alpha, node);
+ up_ll(pt4, M, alpha, node);
+ }
+
+ if (wflag) {
+ for (i = region2->ymin; i < region2->ymax; i++) {
+ node = ll2[i][region2->xmax - 1];
+ for (m = 0; m < M; m++)
+ node[m] = 0.0;
+
+ pt1 = ll1[2 * i][region1->xmax - 1];
+ pt2 = ll1[2 * i + 1][region1->xmax - 1];
+ up_ll(pt1, M, alpha, node);
+ up_ll(pt2, M, alpha, node);
+ }
+ }
+
+ if (hflag) {
+ for (j = region2->xmin; j < region2->xmax; j++) {
+ node = ll2[region2->ymax - 1][j];
+ for (m = 0; m < M; m++)
+ node[m] = 0.0;
+
+ pt1 = ll1[region1->ymax - 1][2 * j];
+ pt2 = ll1[region1->ymax - 1][2 * j + 1];
+ up_ll(pt1, M, alpha, node);
+ up_ll(pt2, M, alpha, node);
+ }
+ }
+
+ if (hflag && wflag) {
+ node = ll2[region2->ymax - 1][region2->xmax - 1];
+ for (m = 0; m < M; m++)
+ node[m] = 0.0;
+
+ pt1 = ll1[region1->ymax - 1][region1->xmax - 1];
+ up_ll(pt1, M, alpha, node);
+ }
+}
+
+
+static void up_ll(LIKELIHOOD * pt1, /* array of log likelihood values, pt1[class] */
+ int M, double alpha, LIKELIHOOD * pt2 /* array of log likelihood values, pt2[class] */
+ )
+{
+ static int m;
+ static double sum, max, cprob[256];
+
+ if (alpha != 1.0) {
+ max = pt1[0];
+ for (m = 1; m < M; m++)
+ if (pt1[m] > max)
+ max = pt1[m];
+
+ sum = 0;
+ for (m = 0; m < M; m++) {
+ cprob[m] = exp(pt1[m] - max);
+ sum += cprob[m];
+ }
+
+ sum = ((1 - alpha) * sum) / M;
+
+ for (m = 0; m < M; m++)
+ pt2[m] += log(sum + alpha * cprob[m]) + max;
+ }
+ else {
+ for (m = 0; m < M; m++)
+ pt2[m] += pt1[m];
+ }
+}
+
+char ***get_pyramid(int w0, int h0, size_t size)
+{
+ char ***pym;
+ int D, wn, hn;
+
+ D = levels(w0, h0);
+ pym = (char ***)G_malloc((D + 1) * sizeof(char *));
+
+ D = 0;
+ wn = w0;
+ hn = h0;
+ pym[D] = (char **)get_img(wn, hn, size);
+ while ((wn > 2) && (hn > 2)) {
+ D++;
+ wn = wn / 2;
+ hn = hn / 2;
+ pym[D] = (char **)get_img(wn, hn, size);
+ }
+
+ return (pym);
+}
+
+
+void free_pyramid(char *pym, int wd, int ht)
+{
+ unsigned char ***pt;
+ int i, D;
+
+ pt = (unsigned char ***)pym;
+ D = levels(wd, ht);
+ for (i = 0; i <= D; i++)
+ free_img(pt[i]);
+ G_free(pym);
+}
+
+char ****get_cubic_pyramid(int w0, int h0, int M, size_t size)
+{
+ char ****pym;
+ int D, wn, hn;
+
+ D = levels(w0, h0);
+ pym = (char ****)G_malloc((D + 1) * sizeof(char *));
+
+ D = 0;
+ wn = w0;
+ hn = h0;
+ pym[D] = (char ***)multialloc(size, 3, hn, wn, M);
+ while ((wn > 2) && (hn > 2)) {
+ D++;
+ wn = wn / 2;
+ hn = hn / 2;
+ pym[D] = (char ***)multialloc(size, 3, hn, wn, M);
+ }
+
+ return (pym);
+}
+
+
+void free_cubic_pyramid(char *pym, int wd, int ht, int M)
+{
+ int i, D;
+ char **pt;
+
+ pt = (char **)pym;
+ D = levels(wd, ht);
+ for (i = 0; i <= D; i++)
+ multifree(pt[i], M);
+ G_free(pym);
+}
+
+int levels(int wd, int ht)
+{
+ int D;
+
+ D = 0;
+ while ((wd > 2) && (ht > 2)) {
+ D++;
+ wd = wd / 2;
+ ht = ht / 2;
+ }
+ return (D);
+}
Property changes on: grass/trunk/imagery/i.smap/decimate.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/eigen.c (from rev 34868, grass/trunk/imagery/i.smap/bouman/eigen.c)
===================================================================
--- grass/trunk/imagery/i.smap/eigen.c (rev 0)
+++ grass/trunk/imagery/i.smap/eigen.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,33 @@
+#include <grass/gis.h>
+
+int eigen(
+ /* Computes eigenvalues (and eigen vectors if desired) for *
+ * symmetric matices. */
+ double **M, /* Input matrix */
+ double *lambda, /* Output eigenvalues */
+ int n /* Input matrix dimension */
+ )
+{
+ int i, j;
+ double **a, *e;
+
+ a = G_alloc_matrix(n, n);
+ e = G_alloc_vector(n);
+
+ for (i = 0; i < n; i++)
+ for (j = 0; j < n; j++)
+ a[i][j] = M[i][j];
+
+ G_tred2(a, n, lambda, e);
+ G_tqli(lambda, e, n, a);
+
+ /* Returns eigenvectors */
+ /* for(i=0; i<n; i++)
+ for(j=0; j<n; j++)
+ M[i][j] = a[i][j]; */
+
+ G_free_matrix(a);
+ G_free_vector(e);
+
+ return 0;
+}
Property changes on: grass/trunk/imagery/i.smap/eigen.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/history.c (from rev 34868, grass/trunk/imagery/i.smap/shapiro/history.c)
===================================================================
--- grass/trunk/imagery/i.smap/history.c (rev 0)
+++ grass/trunk/imagery/i.smap/history.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,12 @@
+#include <grass/gis.h>
+
+void make_history(const char *name, const char *group, const char *subgroup, const char *sigfile)
+{
+ struct History hist;
+
+ if (G_read_history(name, G_mapset(), &hist) >= 0) {
+ sprintf(hist.datsrc_1, "Group/subgroup: %s/%s", group, subgroup);
+ sprintf(hist.datsrc_2, "Sigset file: %s", sigfile);
+ G_write_history(name, &hist);
+ }
+}
Property changes on: grass/trunk/imagery/i.smap/history.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/i.smap.html (from rev 34868, grass/trunk/imagery/i.smap/shapiro/i.smap.html)
===================================================================
--- grass/trunk/imagery/i.smap/i.smap.html (rev 0)
+++ grass/trunk/imagery/i.smap/i.smap.html 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,198 @@
+<H2>DESCRIPTION</H2>
+
+
+The <EM>i.smap</EM> program is used to segment
+multispectral images using a spectral class model known as
+a Gaussian mixture distribution. Since Gaussian mixture
+distributions include conventional multivariate Gaussian
+distributions, this program may also be used to segment
+multispectral images based on simple spectral mean and
+covariance parameters.
+
+<P>
+
+<EM>i.smap</EM> has two modes of operation. The first mode
+is the sequential maximum a posteriori (SMAP) mode
+[<A HREF="#ref1">1</A>,<A HREF="#ref2">2</A>]. The SMAP
+segmentation algorithm attempts to improve segmentation
+accuracy by segmenting the image into regions rather than
+segmenting each pixel separately
+(see <A HREF="#notes">NOTES</A>).
+
+
+<P>
+
+The second mode is the more conventional maximum likelihood (ML)
+classification which classifies each pixel separately,
+but requires somewhat less computation. This mode is selected with
+the <B>-m</B> flag (see <A HREF="#mflag.html">below</A>).
+
+<H2>OPTIONS</H2>
+
+<H3>Flags:</H3>
+
+<DL>
+
+<DT><B>-m</B>
+
+<DD>Use maximum likelihood estimation (instead of smap).
+Normal operation is to use SMAP estimation (see
+<A HREF="#notes">NOTES</A>).
+
+<DT><B>-q</B>
+
+<DD>Run quietly, without printing messages about program
+progress. Without this flag, messages will be printed (to
+stderr) as the program progresses.
+
+</DL>
+
+
+<H3>Parameters:</H3>
+
+<DL>
+<DT><B>group=</B><EM>name</EM>
+
+<DD>imagery group<BR>
+The imagery group that defines the image to be classified.
+
+<DT><B>subgroup=</B><EM>name</EM>
+
+<DD>imagery subgroup<BR>
+The subgroup within the group specified that specifies the
+subset of the band files that are to be used as image data
+to be classified.
+
+<DT><B>signaturefile=</B><EM>name</EM>
+
+<DD>imagery signaturefile<BR>
+The signature file that contains the spectral signatures (i.e., the
+statistics) for the classes to be identified in the image.
+This signature file is produced by the program
+<EM><A HREF="i.gensigset.html">i.gensigset</A></EM>
+(see <A HREF="#notes">NOTES</A>).
+
+<DT><B>blocksize=</B><EM>value</EM>
+
+<DD>size of submatrix to process at one time<BR>
+default: 128<BR>
+This option specifies the size of the "window" to be used when
+reading the image data.
+
+<P>
+
+This program was written to be nice about memory usage
+without influencing the resultant classification. This
+option allows the user to control how much memory is used.
+More memory may mean faster (or slower) operation depending
+on how much real memory your machine has and how much
+virtual memory the program uses.
+
+<P>
+
+The size of the submatrix used in segmenting the image has
+a principle function of controlling memory usage; however,
+it also can have a subtle effect on the quality of the
+segmentation in the smap mode. The smoothing parameters
+for the smap segmentation are estimated separately for each
+submatrix. Therefore, if the image has regions with
+qualitatively different behavior, (e.g., natural woodlands
+and man-made agricultural fields) it may be useful to use a
+submatrix small enough so that different smoothing
+parameters may be used for each distinctive region of the
+image.
+
+<P>
+
+The submatrix size has no effect on the performance of the
+ML segmentation method.
+
+<DT><B>output=</B><EM>name</EM>
+
+<DD>output raster map.<BR>
+The name of a raster map that will contain the
+classification results. This new raster map layer will
+contain categories that can be related to landcover
+categories on the ground.
+
+</DL>
+
+<H2>INTERACTIVE MODE</H2>
+
+If none of the arguments are specified on the command line,
+<EM>i.smap</EM> will interactively prompt for the names of
+the maps and files.
+
+<A NAME="notes"></A><H2>NOTES</H2>
+
+The SMAP algorithm exploits the fact that nearby pixels in
+an image are likely to have the same class. It works by
+segmenting the image at various scales or resolutions and
+using the coarse scale segmentations to guide the finer
+scale segmentations. In addition to reducing the number of
+misclassifications, the SMAP algorithm generally produces
+segmentations with larger connected regions of a fixed
+class which may be useful in some applications.
+
+<P>
+
+The amount of smoothing that is performed in the
+segmentation is dependent of the behavior of the data in
+the image. If the data suggests that the nearby pixels
+often change class, then the algorithm will adaptively
+reduce the amount of smoothing. This ensures that
+excessively large regions are not formed.
+
+<P>
+
+The module i.smap 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.
+<p>
+r.mapcalc MASKed_map=classification results
+
+<H2>REFERENCES</H2>
+
+<OL>
+<LI>C. Bouman and M. Shapiro,
+"Multispectral Image Segmentation using a Multiscale Image Model",
+<EM>Proc. of IEEE Int'l Conf. on Acoust., Speech and Sig. Proc.,</EM>
+pp. III-565 - III-568, San Francisco, California, March
+23-26, 1992.
+
+<LI>C. Bouman and M. Shapiro 1994,
+"A Multiscale Random Field
+Model for Bayesian Image Segmentation",
+<EM>IEEE Trans. on Image Processing., 3(2), 162-177"</EM>
+
+<LI>McCauley, J.D. and B.A. Engel 1995,
+"Comparison of Scene Segmentations: SMAP, ECHO and Maximum Likelyhood",
+<EM>IEEE Trans. on Geoscience and Remote Sensing, 33(6): 1313-1316.</EM>
+</OL>
+
+<H2>SEE ALSO</H2>
+
+<EM><A HREF="i.group.html">i.group</A></EM>
+for creating groups and subgroups
+
+<P>
+
+<EM><A HREF="r.mapcalc.html">r.mapcalc</A></EM>
+to copy classification result in order to cut out MASKed subareas
+
+<P>
+
+<EM><A HREF="i.gensigset.html">i.gensigset</A></EM>
+to generate the signature file required by this program
+
+<H2>AUTHORS</H2>
+
+<a href="http://dynamo.ecn.purdue.edu/~bouman/software/segmentation/">Charles Bouman,
+School of Electrical Engineering, Purdue University</a>
+
+<BR>
+Michael Shapiro,
+U.S.Army Construction Engineering
+Research Laboratory
+
+<p><i>Last changed: $Date$</i>
Property changes on: grass/trunk/imagery/i.smap/i.smap.html
___________________________________________________________________
Name: svn:mime-type
+ text/html
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/interp.c (from rev 34870, grass/trunk/imagery/i.smap/bouman/interp.c)
===================================================================
--- grass/trunk/imagery/i.smap/interp.c (rev 0)
+++ grass/trunk/imagery/i.smap/interp.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,382 @@
+#include <stdlib.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "bouman.h"
+#include "region.h"
+
+#define EM_PRECISION 1e-4
+#define ML_PRECISION 1e-6
+
+static void seq_MAP_routine(unsigned char ***, struct Region *,
+ LIKELIHOOD ****, int, double *, int);
+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);
+void MLE(unsigned char **, LIKELIHOOD ***, struct Region *, int);
+static int up_char(int, int, struct Region *, unsigned char **,
+ unsigned char **);
+
+
+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 vlevel /* verbose output */
+ )
+{
+ int repeat;
+
+ /* Repeat segmentation to get values for alpha_dec */
+ for (repeat = 0; repeat < 2; repeat++) {
+ /* Construct image log likelihood pyramid */
+ make_pyramid(ll_pym, region, M, alpha_dec, vlevel);
+ if (vlevel >= 2)
+ G_message(_("pyramid constructed."));
+
+ /* Perform sequential MAP segmentation using EM algorithm */
+ seq_MAP_routine(sf_pym, region, ll_pym, M, alpha_dec, vlevel);
+ }
+}
+
+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 vlevel /* verbose output */
+ )
+{
+ int j, k; /* loop index */
+ int wd, ht; /* width and height at each resolution */
+ int *period; /* sampling period at each resolution */
+ int D; /* number of resolutions -1 */
+ double ***N; /* transition probability statistics; N[2][3][2] */
+ double alpha[3]; /* transition probability parameters */
+ double tmp[3]; /* temporary transition probability parameters */
+ double diff1; /* change in parameter estimates */
+ double diff2; /* change in log likelihood */
+ struct Region *regionary; /* array of region stuctures */
+
+ /* determine number of resolutions */
+ D = levels_reg(region);
+
+ /* allocate memory */
+ if ((N = (double ***)multialloc(sizeof(double), 3, 2, 3, 2)) == NULL)
+ G_fatal_error(_("Unable to allocate memory"));
+
+ regionary = (struct Region *)G_malloc((D + 1) * sizeof(struct Region));
+ period = (int *)G_malloc(D * sizeof(int));
+
+ /* Compute the image region at each resolution. */
+ k = 0;
+ copy_reg(region, &(regionary[k]));
+ reg_to_wdht(&(regionary[k]), &wd, &ht);
+ while ((wd > 2) && (ht > 2)) {
+ copy_reg(&(regionary[k]), &(regionary[k + 1]));
+ dec_reg(&(regionary[k + 1]));
+ reg_to_wdht(&(regionary[k + 1]), &wd, &ht);
+ k++;
+ }
+
+ /* Compute sampling period for EM algorithm at each resolution. */
+ for (k = 0; k < D; k++) {
+ period[k] = (int)pow(2.0, (D - k - 2) / 2.0);
+ if (period[k] < 1)
+ period[k] = 1;
+ }
+
+ /* Compute Maximum Likelihood estimate at coarsest resolution */
+ MLE(sf_pym[D], ll_pym[D], &(regionary[D]), M);
+
+ /* Initialize the transition parameters */
+ alpha[0] = 0.5 * (3.0 / 7.0);
+ alpha[1] = 0.5 * (2.0 / 7.0);
+ alpha[2] = 0.0;
+
+ /* Interpolate the classification at each resolution */
+ for (D--; D >= 0; D--) {
+ if (vlevel >= 2)
+ G_debug(1, "Resolution = %d; period = %d", D, period[D]);
+
+ for (j = 0; j < 3; j++)
+ alpha[j] *= (1 - EM_PRECISION * 10);
+ if (vlevel >= 4)
+ print_alpha(alpha);
+ /* Apply EM algorithm to estimate alpha. Continue for *
+ * 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);
+ if (vlevel >= 4)
+ print_N(N);
+ if (vlevel >= 4)
+ G_debug(1, "log likelihood = %f", log_like(N, alpha, M));
+ for (j = 0; j < 3; j++)
+ tmp[j] = alpha[j];
+
+ alpha_max(N, alpha, M, ML_PRECISION);
+ if (vlevel >= 2)
+ print_alpha(alpha);
+ if (vlevel >= 4)
+ G_debug(1, "log likelihood = %f", log_like(N, alpha, M));
+
+ for (diff1 = j = 0; j < 3; j++)
+ 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);
+ alpha_dec[D] = alpha_dec_max(N);
+
+ if (vlevel >= 4)
+ print_N(N);
+ if (vlevel >= 2) {
+ alpha_max(N, alpha, M, ML_PRECISION);
+ print_alpha(alpha);
+ }
+ }
+
+ /* free up N */
+ G_free((char *)regionary);
+ G_free((char *)period);
+ multifree((char *)N, 3);
+}
+
+static double alpha_dec_max(double ***N)
+{
+ int i, j, k;
+ double N_marg[2]; /* Marginal transition rate */
+ double N_sum; /* total number of transitions counted */
+
+ for (k = 0; k < 2; k++) {
+ N_marg[k] = 0;
+ for (i = 0; i < 3; i++)
+ for (j = 0; j < 2; j++) {
+ N_marg[k] += N[k][i][j];
+ }
+ }
+
+ N_sum = N_marg[0] + N_marg[1];
+
+ if (N_sum == 0)
+ return (0.0);
+
+ return (N_marg[1] / N_sum);
+}
+
+static void print_N(
+ /* prints out class transition statistics */
+ double ***N)
+{
+ int n0, n1, n2;
+
+ G_debug(2, "Class transition statistics");
+ for (n0 = 0; n0 < 2; n0++) {
+ for (n1 = 0; n1 < 3; n1++) {
+ for (n2 = 0; n2 < 2; n2++)
+ G_debug(3, " %f", N[n0][n1][n2]);
+ }
+ }
+}
+
+
+static void print_alpha(
+ /* prints out transition parameters. */
+ double *alpha)
+{
+ G_debug(1, "Transition probabilities: %f %f %f; %f",
+ alpha[0], alpha[1], alpha[2],
+ 1.0 - alpha[0] - 2 * alpha[1] - alpha[2]);
+}
+
+
+static void interp(
+ /* Estimates finer resolution segmentation from coarser resolution
+ segmentation and texture statistics. */
+ unsigned char **sf1, /* finer resolution segmentation */
+ struct Region *region, /* image region */
+ unsigned char **sf2, /* coarser resolution segmentation */
+ LIKELIHOOD *** ll, /* Log likelihood ll[i][j][class] */
+ int M, /* number of classes */
+ 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 i, j; /* pixel index */
+ int m; /* class index */
+ int bflag; /* boundary flag */
+ int nn0, nn1, nn2; /* transition counts */
+ int *n0, *n1, *n2; /* transition counts for each possible pixel class */
+ unsigned char *nbr[8]; /* pointers to neighbors at courser resolution */
+ double cost, mincost; /* cost of class selection; minimum cost */
+ int best = 0; /* class of minimum cost selection */
+ double Constant, tmp;
+ double *pdf; /* propability density function of class selections */
+ double Z; /* normalizing costant for pdf */
+ double alpha0, alpha1, alpha2; /* transition probabilities */
+ double log_tbl[2][3][2]; /* log of transition probability */
+
+ /* allocate memory for pdf */
+ pdf = (double *)G_malloc(M * sizeof(double));
+ n0 = (int *)G_malloc(M * sizeof(int));
+ n1 = (int *)G_malloc(M * sizeof(int));
+ n2 = (int *)G_malloc(M * sizeof(int));
+
+ /* set constants */
+ alpha0 = alpha[0];
+ alpha1 = alpha[1];
+ alpha2 = alpha[2];
+ Constant = (1 - alpha0 - 2 * alpha1 - alpha2) / M;
+ if (Constant < 0)
+ G_fatal_error(_("Invalid parameter values"));
+
+ /* precompute logs and zero static vector */
+ for (nn0 = 0; nn0 < 2; nn0++)
+ for (nn1 = 0; nn1 < 3; nn1++)
+ for (nn2 = 0; nn2 < 2; nn2++) {
+ tmp = (alpha0 * nn0 + alpha1 * nn1 + alpha2 * nn2) + Constant;
+ if (tmp == 0)
+ log_tbl[nn0][nn1][nn2] = HUGE_VAL;
+ else
+ log_tbl[nn0][nn1][nn2] = -log(tmp);
+ if (statflag)
+ N[nn0][nn1][nn2] = 0;
+ }
+
+ /* classify points and compute expectation of N */
+ for (i = region->ymin; i < region->ymax; i += period)
+ for (j = region->xmin; j < region->xmax; j += period) {
+ /* compute minimum cost class */
+ mincost = HUGE_VAL;
+ bflag = up_char(i, j, region, sf2, nbr);
+ for (m = 0; m < M; m++) {
+ nn0 = n0[m] = (m == (*nbr[0]));
+ nn1 = n1[m] = (m == (*nbr[1])) + (m == (*nbr[2]));
+ nn2 = n2[m] = (m == (*nbr[3]));
+
+ pdf[m] = cost = log_tbl[nn0][nn1][nn2] - ll[i][j][m];
+ if (cost < mincost) {
+ mincost = cost;
+ best = m;
+ }
+ }
+ sf1[i][j] = best;
+
+ /* if not on boundary, compute expectation of N */
+ if ((!bflag) && (statflag)) {
+ Z = 0.0;
+ for (m = 0; m < M; m++) {
+ if (pdf[m] == HUGE_VAL)
+ pdf[m] = 0;
+ else
+ pdf[m] = exp(mincost - pdf[m]);
+ Z += pdf[m];
+ }
+ for (m = 0; m < M; m++)
+ N[n0[m]][n1[m]][n2[m]] += pdf[m] / Z;
+ }
+ }
+
+ G_free((char *)pdf);
+ G_free((char *)n0);
+ G_free((char *)n1);
+ G_free((char *)n2);
+}
+
+void MLE( /* computes maximum likelihood classification */
+ unsigned char **sf, /* segmentation classes */
+ LIKELIHOOD *** ll, /* texture statistics */
+ struct Region *region, /* image region */
+ int M /* number of classes */
+ )
+{
+ int i, j, m, best;
+ double max;
+
+ for (i = region->ymin; i < region->ymax; i++)
+ for (j = region->xmin; j < region->xmax; j++) {
+ max = ll[i][j][0];
+ best = 0;
+ for (m = 1; m < M; m++) {
+ if (max < ll[i][j][m]) {
+ max = ll[i][j][m];
+ best = m;
+ }
+ }
+ sf[i][j] = best;
+ }
+}
+
+
+static int up_char(
+ /* Computes list of pointers to nieghbors at next coarser resolution. *
+ * Returns flag when on boundary. */
+ int i, int j, /* fine resolution pixel location */
+ struct Region *region, /* fine resolution image region */
+ unsigned char **img, /* course resolution image */
+ unsigned char **pt /* list of pointers */
+ )
+{
+ static int xmax, ymax;
+ static int bflag; /* =1 when on boundary */
+ static int i2, j2; /* base indices at coarser level */
+ static int di, dj; /* displacements at coarser level */
+
+ /* create new xmax and ymax */
+ xmax = region->xmax;
+ ymax = region->ymax;
+
+ /* check for images of odd length */
+ if (xmax & 1) {
+ xmax--;
+ if (j == xmax)
+ j--;
+ }
+ if (ymax & 1) {
+ ymax--;
+ if (i == ymax)
+ i--;
+ }
+
+ /* compute indices */
+ di = ((i & 1) << 1) - 1;
+ dj = ((j & 1) << 1) - 1;
+ i2 = i >> 1;
+ j2 = j >> 1;
+
+ /* enforce an absorptive boundary */
+ bflag = 0;
+ if ((i == region->ymin) && region->free.top) {
+ di = 0;
+ bflag = 1;
+ }
+ if ((i == ymax - 1) && region->free.bottom) {
+ di = 0;
+ bflag = 1;
+ }
+ /* mod shapiro */
+ /* if((j==region->ymax)&®ion->free.left) { dj=0; bflag=1; } */
+ if ((j == region->xmin) && region->free.left) {
+ dj = 0;
+ bflag = 1;
+ }
+ /* end mod shapiro */
+ if ((j == xmax - 1) && region->free.right) {
+ dj = 0;
+ bflag = 1;
+ }
+
+ /* compute pointers */
+ pt[0] = img[i2] + j2;
+ pt[1] = img[i2] + j2 + dj;
+ pt[2] = img[i2 + di] + j2;
+ pt[3] = img[i2 + di] + j2 + dj;
+
+ return (bflag);
+}
Property changes on: grass/trunk/imagery/i.smap/interp.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/invert.c (from rev 34868, grass/trunk/imagery/i.smap/bouman/invert.c)
===================================================================
--- grass/trunk/imagery/i.smap/invert.c (rev 0)
+++ grass/trunk/imagery/i.smap/invert.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,38 @@
+#include <math.h>
+#include <grass/gis.h>
+
+int invert(
+ /* inverts a matrix of arbitrary size input as a 2D array. */
+ double **a, /* input/output matrix */
+ int n /* dimension */
+ )
+{
+ int status;
+ int i, j, *indx;
+ double **y, *col, d;
+
+ indx = G_alloc_ivector(n);
+ y = G_alloc_matrix(n, n);
+ col = G_alloc_vector(n);
+
+ if ((status = G_ludcmp(a, n, indx, &d))) {
+ for (j = 0; j < n; j++) {
+ for (i = 0; i < n; i++)
+ col[i] = 0.0;
+ col[j] = 1.0;
+ G_lubksb(a, n, indx, col);
+ for (i = 0; i < n; i++)
+ y[i][j] = col[i];
+ }
+
+ for (i = 0; i < n; i++)
+ for (j = 0; j < n; j++)
+ a[i][j] = y[i][j];
+ }
+
+ G_free_ivector(indx);
+ G_free_matrix(y);
+ G_free_vector(col);
+
+ return (status);
+}
Property changes on: grass/trunk/imagery/i.smap/invert.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/labels.c (from rev 34868, grass/trunk/imagery/i.smap/shapiro/labels.c)
===================================================================
--- grass/trunk/imagery/i.smap/labels.c (rev 0)
+++ grass/trunk/imagery/i.smap/labels.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,16 @@
+#include <grass/imagery.h>
+#include "bouman.h"
+
+int create_output_labels(struct SigSet *S, struct files *files)
+{
+ int n;
+ struct ClassSig *C;
+
+ G_init_cats((CELL) 0, S->title, &files->output_labels);
+ for (n = 0; n < S->nclasses; n++) {
+ C = &S->ClassSig[n];
+ G_set_cat((CELL) C->classnum, C->title, &files->output_labels);
+ }
+
+ return 0;
+}
Property changes on: grass/trunk/imagery/i.smap/labels.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/local_proto.h (from rev 34868, grass/trunk/imagery/i.smap/shapiro/local_proto.h)
===================================================================
--- grass/trunk/imagery/i.smap/local_proto.h (rev 0)
+++ grass/trunk/imagery/i.smap/local_proto.h 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,6 @@
+/* history.c */
+void make_history(const char *, const char *, const char *, const char *);
+
+/* opencell.c */
+int open_cell_old(const char *, const char *);
+int open_cell_new(const char *);
Property changes on: grass/trunk/imagery/i.smap/local_proto.h
___________________________________________________________________
Name: svn:mime-type
+ text/x-chdr
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/main.c (from rev 34870, grass/trunk/imagery/i.smap/shapiro/main.c)
===================================================================
--- grass/trunk/imagery/i.smap/main.c (rev 0)
+++ grass/trunk/imagery/i.smap/main.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,55 @@
+
+/****************************************************************************
+ *
+ * MODULE: i.smap
+ * AUTHOR(S): Michael Shapiro (USACERL) (original contributor)
+ * Markus Neteler <neteler itc.it>,
+ * Roberto Flor <flor itc.it>,
+ * Bernhard Reiter <bernhard intevation.de>,
+ * Brad Douglas <rez touchofmadness.com>,
+ * Glynn Clements <glynn gclements.plus.com>,
+ * Jan-Oliver Wagner <jan intevation.de>
+ * PURPOSE: segment multispectral images using a spectral class model
+ * known as a Gaussian mixture distribution
+ * COPYRIGHT: (C) 1999-2008 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 <stdlib.h>
+#include <unistd.h>
+#include <grass/imagery.h>
+#include <grass/glocale.h>
+#include "bouman.h"
+
+
+int main(int argc, char *argv[])
+{
+ struct parms parms; /* command line parms */
+ struct files files; /* file descriptors, io, buffers */
+ struct SigSet S;
+ struct GModule *module;
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ module->keywords = _("imagery, classification, supervised, SMAP");
+ module->description =
+ _("Performs contextual image classification "
+ "using sequential maximum a posteriori (SMAP) estimation.");
+
+ parse(argc, argv, &parms);
+ openfiles(&parms, &files);
+ read_signatures(&parms, &S);
+ create_output_labels(&S, &files);
+
+ segment(&S, &parms, &files);
+
+ closefiles(&parms, &files);
+
+ G_done_msg(" ");
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass/trunk/imagery/i.smap/main.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/model.c (from rev 34870, grass/trunk/imagery/i.smap/bouman/model.c)
===================================================================
--- grass/trunk/imagery/i.smap/model.c (rev 0)
+++ grass/trunk/imagery/i.smap/model.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,161 @@
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/imagery.h>
+#include "bouman.h"
+#include "region.h"
+
+#define PI M_PI
+
+
+void extract_init(struct SigSet *S)
+{
+ int m;
+ int i;
+ int b1, b2;
+ int nbands;
+ double *lambda;
+ struct ClassSig *C;
+ struct SubSig *SubS;
+
+ nbands = S->nbands;
+ /* allocate scratch memory */
+ lambda = (double *)G_malloc(nbands * sizeof(double));
+
+ /* invert matrix and compute constant for each subclass */
+
+ /* for each class */
+ for (m = 0; m < S->nclasses; m++) {
+ C = &(S->ClassSig[m]);
+
+ /* for each subclass */
+ for (i = 0; i < C->nsubclasses; i++) {
+ SubS = &(C->SubSig[i]);
+
+ /* Test for symetric matrix */
+ for (b1 = 0; b1 < nbands; b1++)
+ for (b2 = 0; b2 < nbands; b2++) {
+ if (SubS->R[b1][b2] != SubS->R[b2][b1])
+ G_warning(_("Nonsymetric covariance for class %d subclass %d"),
+ m + 1, i + 1);
+
+ SubS->Rinv[b1][b2] = SubS->R[b1][b2];
+ }
+
+ /* Test for positive definite matrix */
+ eigen(SubS->Rinv, lambda, nbands);
+ for (b1 = 0; b1 < nbands; b1++) {
+ if (lambda[b1] <= 0.0)
+ G_warning(_("Nonpositive eigenvalues for class %d subclass %d"),
+ m + 1, i + 1);
+ }
+
+ /* Precomputes the cnst */
+ SubS->cnst = (-nbands / 2.0) * log(2 * PI);
+ for (b1 = 0; b1 < nbands; b1++) {
+ SubS->cnst += -0.5 * log(lambda[b1]);
+ }
+
+ /* Precomputes the inverse of tex->R */
+ invert(SubS->Rinv, nbands);
+ }
+ }
+ G_free((char *)lambda);
+}
+
+
+void extract(DCELL *** img, /* multispectral image, img[band][i][j] */
+ struct Region *region, /* region to extract */
+ LIKELIHOOD *** ll, /* log likelihood, ll[i][j][class] */
+ struct SigSet *S /* class signatures */
+ )
+{
+ int i, j; /* column and row indexes */
+ int m; /* class index */
+ int k; /* subclass index */
+ int b1, b2; /* spectral index */
+ int no_data; /* no data flag */
+ int max_nsubclasses; /* maximum number of subclasses */
+ int nbands; /* number of spectral bands */
+ double *subll; /* log likelihood of subclasses */
+ double *diff;
+ double maxlike = 0.0L;
+ double subsum;
+ struct ClassSig *C;
+ struct SubSig *SubS;
+
+ nbands = S->nbands;
+
+ /* determine the maximum number of subclasses */
+ max_nsubclasses = 0;
+ for (m = 0; m < S->nclasses; m++)
+ if (S->ClassSig[m].nsubclasses > max_nsubclasses)
+ max_nsubclasses = S->ClassSig[m].nsubclasses;
+
+ /* allocate memory */
+ diff = (double *)G_malloc(nbands * sizeof(double));
+ subll = (double *)G_malloc(max_nsubclasses * sizeof(double));
+
+ /* Compute log likelihood at each pixel and for every class. */
+
+ /* for each pixel */
+ for (i = region->ymin; i < region->ymax; i++)
+ for (j = region->xmin; j < region->xmax; j++) {
+
+ /* Check for no data condition */
+ no_data = 1;
+ for (b1 = 0; (b1 < nbands) && no_data; b1++)
+ no_data = no_data && (G_is_d_null_value(&img[b1][i][j]));
+
+ if (no_data) {
+ for (m = 0; m < S->nclasses; m++)
+ ll[i][j][m] = 0.0;
+ }
+ else {
+ /* for each class */
+ for (m = 0; m < S->nclasses; m++) {
+ C = &(S->ClassSig[m]);
+
+ /* compute log likelihood for each subclass */
+ for (k = 0; k < C->nsubclasses; k++) {
+ SubS = &(C->SubSig[k]);
+ subll[k] = SubS->cnst;
+ for (b1 = 0; b1 < nbands; b1++) {
+ diff[b1] = img[b1][i][j] - SubS->means[b1];
+ subll[k] -=
+ 0.5 * diff[b1] * diff[b1] *
+ SubS->Rinv[b1][b1];
+ }
+ for (b1 = 0; b1 < nbands; b1++)
+ for (b2 = b1 + 1; b2 < nbands; b2++)
+ subll[k] -=
+ diff[b1] * diff[b2] * SubS->Rinv[b1][b2];
+ }
+
+ /* shortcut for one subclass */
+ if (C->nsubclasses == 1) {
+ ll[i][j][m] = subll[0];
+ }
+ /* compute mixture likelihood */
+ else {
+ /* find the most likely subclass */
+ for (k = 0; k < C->nsubclasses; k++) {
+ if (k == 0)
+ maxlike = subll[k];
+ if (subll[k] > maxlike)
+ maxlike = subll[k];
+ }
+
+ /* Sum weighted subclass likelihoods */
+ subsum = 0;
+ for (k = 0; k < C->nsubclasses; k++)
+ subsum +=
+ exp(subll[k] - maxlike) * C->SubSig[k].pi;
+
+ ll[i][j][m] = log(subsum) + maxlike;
+ }
+ }
+ }
+ }
+ G_free((char *)diff);
+ G_free((char *)subll);
+}
Property changes on: grass/trunk/imagery/i.smap/model.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/multialloc.c (from rev 34870, grass/trunk/imagery/i.smap/bouman/multialloc.c)
===================================================================
--- grass/trunk/imagery/i.smap/multialloc.c (rev 0)
+++ grass/trunk/imagery/i.smap/multialloc.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,124 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdarg.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+
+/* modified from dynamem.c on 4/29/91 C. Bouman */
+
+/*
+ * multialloc( s, d, dn, dn ....) allocates a d dimensional array, whose
+ * dimensions are stored in a list starting at d1. Each array element is
+ * of size s.
+ */
+
+char *multialloc(size_t s, /* individual array element size */
+ int d, /* number of dimensions */
+ ...)
+{
+ va_list ap; /* varargs list traverser */
+ int max, /* size of array to be declared */
+ *q; /* pointer to dimension list */
+ char **r, /* pointer to beginning of the array of the
+ * pointers for a dimension */
+ **s1, *t, *tree; /* base pointer to beginning of first array */
+ int i, j; /* loop counters */
+ int *d1; /* dimension list */
+
+ va_start(ap, d);
+ d1 = (int *)G_malloc(d * sizeof(int));
+
+ for (i = 0; i < d; i++)
+ d1[i] = va_arg(ap, int);
+
+ r = &tree;
+ q = d1; /* first dimension */
+ max = 1;
+ for (i = 0; i < d - 1; i++, q++) { /* for each of the dimensions
+ * but the last */
+ max *= (*q);
+ r[0] = (char *)G_malloc(max * sizeof(char **));
+ r = (char **)r[0]; /* step through to beginning of next
+ * dimension array */
+ }
+ max *= s * (*q); /* grab actual array memory */
+ r[0] = (char *)G_malloc(max * sizeof(char));
+
+ /*
+ * r is now set to point to the beginning of each array so that we can
+ * use it to scan down each array rather than having to go across and
+ * then down
+ */
+ r = (char **)tree; /* back to the beginning of list of arrays */
+ q = d1; /* back to the first dimension */
+ max = 1;
+ for (i = 0; i < d - 2; i++, q++) { /* we deal with the last
+ * array of pointers later on */
+ max *= (*q); /* number of elements in this dimension */
+ for (j = 1, s1 = r + 1, t = r[0]; j < max; j++) { /* scans down array for
+ * first and subsequent
+ * elements */
+
+ /* modify each of the pointers so that it points to
+ * the correct position (sub-array) of the next
+ * dimension array. s1 is the current position in the
+ * current array. t is the current position in the
+ * next array. t is incremented before s1 is, but it
+ * starts off one behind. *(q+1) is the dimension of
+ * the next array. */
+
+ *s1 = (t += sizeof(char **) * *(q + 1));
+ s1++;
+ }
+ r = (char **)r[0]; /* step through to begining of next
+ * dimension array */
+ }
+ max *= (*q); /* max is total number of elements in the
+ * last pointer array */
+
+ /* same as previous loop, but different size factor */
+ for (j = 1, s1 = r + 1, t = r[0]; j < max; j++)
+ *s1++ = (t += s * *(q + 1));
+
+ va_end(ap);
+ G_free((char *)d1);
+
+ return tree; /* return base pointer */
+}
+
+
+/*
+ * multifree releases all memory that we have already declared analogous to
+ * G_free () when using G_malloc ()
+ */
+void multifree(char *r, int d)
+{
+ char **p;
+ char *next = NULL;
+ int i;
+
+ for (p = (char **)r, i = 0; i < d; p = (char **)next, i++) {
+ if (p != NULL) {
+ next = *p;
+ G_free((char *)p);
+ }
+ }
+}
+
+
+unsigned char **get_img(int wd, int ht, size_t size)
+{
+ char *pt;
+
+ if ((pt = multialloc(size, 2, ht, wd)) == NULL)
+ G_fatal_error(_("Out of memory"));
+
+ return ((unsigned char **)pt);
+}
+
+
+void free_img(unsigned char **pt)
+{
+ multifree((char *)pt, 2);
+}
Property changes on: grass/trunk/imagery/i.smap/multialloc.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/opencell.c (from rev 34868, grass/trunk/imagery/i.smap/shapiro/opencell.c)
===================================================================
--- grass/trunk/imagery/i.smap/opencell.c (rev 0)
+++ grass/trunk/imagery/i.smap/opencell.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,34 @@
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+
+int open_cell_old(const char *name, const char *mapset)
+{
+ int fd;
+
+ if (mapset == NULL)
+ mapset = G_find_cell2(name, "");
+ fd = G_open_cell_old(name, mapset);
+ if (fd >= 0)
+ return fd;
+
+ G_fatal_error(_("Unable to open raster map <%s>"), name);
+
+ /* should not get here */
+ return -1;
+}
+
+int open_cell_new(const char *name)
+{
+ int fd;
+
+ fd = G_open_cell_new(name);
+ if (fd >= 0)
+ return fd;
+
+ G_fatal_error(_("Unable to create raster map <%s>"), name);
+
+ /* should not get here */
+ return -1;
+}
Property changes on: grass/trunk/imagery/i.smap/opencell.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/openfiles.c (from rev 34870, grass/trunk/imagery/i.smap/shapiro/openfiles.c)
===================================================================
--- grass/trunk/imagery/i.smap/openfiles.c (rev 0)
+++ grass/trunk/imagery/i.smap/openfiles.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,39 @@
+#include <stdlib.h>
+#include <grass/imagery.h>
+#include <grass/glocale.h>
+#include "bouman.h"
+#include "local_proto.h"
+
+
+int openfiles(struct parms *parms, struct files *files)
+{
+ struct Ref Ref; /* subgroup reference list */
+ int n;
+
+
+ if (!I_get_subgroup_ref(parms->group, parms->subgroup, &Ref))
+ G_fatal_error(_("Unable to read REF file for subgroup <%s> in group <%s>"),
+ parms->subgroup, parms->group);
+
+ if (Ref.nfiles <= 0)
+ G_fatal_error(_("Subgroup <%s> in group <%s> contains no raster maps"),
+ parms->subgroup, parms->group);
+
+ /* allocate file descriptors, and io buffer */
+ files->cellbuf = G_allocate_d_raster_buf();
+ files->outbuf = G_allocate_c_raster_buf();
+ files->isdata = G_malloc(G_window_cols());
+
+ files->nbands = Ref.nfiles;
+ files->band_fd = (int *)G_calloc(Ref.nfiles, sizeof(int));
+
+ /* open all group maps for reading */
+ for (n = 0; n < Ref.nfiles; n++)
+ files->band_fd[n] =
+ open_cell_old(Ref.file[n].name, Ref.file[n].mapset);
+
+ /* open output map */
+ files->output_fd = open_cell_new(parms->output_map);
+
+ return 0;
+}
Property changes on: grass/trunk/imagery/i.smap/openfiles.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/parse.c (from rev 34870, grass/trunk/imagery/i.smap/shapiro/parse.c)
===================================================================
--- grass/trunk/imagery/i.smap/parse.c (rev 0)
+++ grass/trunk/imagery/i.smap/parse.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,68 @@
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/imagery.h>
+#include "bouman.h"
+
+
+int parse(int argc, char *argv[], struct parms *parms)
+{
+ struct Option *group, *subgroup, *sigfile, *output;
+ struct Option *blocksize;
+ struct Flag *quiet;
+ struct Flag *ml;
+
+ group = G_define_standard_option(G_OPT_I_GROUP);
+
+ subgroup = G_define_standard_option(G_OPT_I_SUBGROUP);
+
+ sigfile = G_define_option();
+ sigfile->key = "signaturefile";
+ sigfile->label = _("Name of file containing signatures");
+ sigfile->description = _("Generated by i.gensigset");
+ sigfile->key_desc = "name";
+ sigfile->required = YES;
+ sigfile->type = TYPE_STRING;
+
+ output = G_define_standard_option(G_OPT_R_OUTPUT);
+
+ blocksize = G_define_option();
+ blocksize->key = "blocksize";
+ blocksize->description = _("Size of submatrix to process at one time");
+ blocksize->required = NO;
+ blocksize->type = TYPE_INTEGER;
+ blocksize->answer = "128";
+
+ ml = G_define_flag();
+ ml->key = 'm';
+ ml->description =
+ _("Use maximum likelihood estimation (instead of smap)");
+
+ quiet = G_define_flag();
+ quiet->key = 'q';
+ quiet->description = _("Run quietly");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ parms->quiet = quiet->answer;
+ parms->ml = ml->answer;
+
+ parms->output_map = output->answer;
+ parms->group = group->answer;
+ parms->subgroup = subgroup->answer;
+ parms->sigfile = sigfile->answer;
+
+ /* check all the inputs */
+ if (!I_find_group(parms->group))
+ G_fatal_error(_("Group <%s> not found"), parms->group);
+
+ if (!I_find_subgroup(parms->group, parms->subgroup))
+ G_fatal_error(_("Subgroup <%s> not found"), parms->subgroup);
+
+ if (sscanf(blocksize->answer, "%d", &parms->blocksize) != 1
+ || parms->blocksize <= 8)
+ parms->blocksize = 8;
+
+ return 0;
+}
Property changes on: grass/trunk/imagery/i.smap/parse.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/read_block.c (from rev 34870, grass/trunk/imagery/i.smap/bouman/read_block.c)
===================================================================
--- grass/trunk/imagery/i.smap/read_block.c (rev 0)
+++ grass/trunk/imagery/i.smap/read_block.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,24 @@
+#include <stdlib.h>
+#include <grass/imagery.h>
+#include <grass/glocale.h>
+#include "bouman.h"
+#include "region.h"
+
+int read_block(DCELL *** img, /* img[band][row[col] */
+ struct Region *region, struct files *files)
+{
+ int band, row, col;
+
+ for (band = 0; band < files->nbands; band++) {
+ for (row = region->ymin; row < region->ymax; row++) {
+ if (G_get_d_raster_row(files->band_fd[band], files->cellbuf, row)
+ < 0)
+ G_fatal_error(_("Unable to read raster map row %d"),
+ row);
+ for (col = region->xmin; col < region->xmax; col++)
+ img[band][row][col] = files->cellbuf[col];
+ }
+ }
+
+ return 0;
+}
Property changes on: grass/trunk/imagery/i.smap/read_block.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/read_sig.c (from rev 34870, grass/trunk/imagery/i.smap/shapiro/read_sig.c)
===================================================================
--- grass/trunk/imagery/i.smap/read_sig.c (rev 0)
+++ grass/trunk/imagery/i.smap/read_sig.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,32 @@
+#include <stdlib.h>
+#include <grass/imagery.h>
+#include <grass/glocale.h>
+#include "bouman.h"
+
+
+int read_signatures(struct parms *parms, struct SigSet *S)
+{
+ FILE *fd;
+ struct Ref Ref;
+
+ if (!I_get_subgroup_ref(parms->group, parms->subgroup, &Ref))
+ G_fatal_error(_("Unable to read REF file for subgroup <%s> in group <%s>"),
+ parms->subgroup, parms->group);
+
+ if (Ref.nfiles <= 0)
+ G_fatal_error(_("Subgroup <%s> in group <%s> contains no raster maps"),
+ parms->subgroup, parms->group);
+
+ fd = I_fopen_sigset_file_old(parms->group, parms->subgroup,
+ parms->sigfile);
+ if (fd == NULL)
+ G_fatal_error(_("Unable to read signature file <%s>"),
+ parms->sigfile);
+
+ if (I_ReadSigSet(fd, S) < 0 || Ref.nfiles != S->nbands)
+ G_fatal_error(_("Signature file <%s> is invalid"), parms->sigfile);
+
+ fclose(fd);
+
+ return 0;
+}
Property changes on: grass/trunk/imagery/i.smap/read_sig.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/reg_util.c (from rev 34868, grass/trunk/imagery/i.smap/bouman/reg_util.c)
===================================================================
--- grass/trunk/imagery/i.smap/reg_util.c (rev 0)
+++ grass/trunk/imagery/i.smap/reg_util.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,52 @@
+#include "bouman.h"
+#include "region.h"
+
+int levels_reg(struct Region *region)
+{
+ int D;
+ int wd, ht;
+ struct Region region_buff;
+
+ /* save region stucture */
+ copy_reg(region, ®ion_buff);
+
+ D = 0;
+ reg_to_wdht(region, &wd, &ht);
+ while ((wd > 2) && (ht > 2)) {
+ D++;
+ dec_reg(region);
+ reg_to_wdht(region, &wd, &ht);
+ }
+
+ /* restore region structure */
+ copy_reg(®ion_buff, region);
+
+ return (D);
+}
+
+void dec_reg(struct Region *region)
+{
+ region->xmin = region->xmin / 2;
+ region->xmax = region->xmax / 2;
+ region->ymin = region->ymin / 2;
+ region->ymax = region->ymax / 2;
+}
+
+void copy_reg(struct Region *region1, struct Region *region2)
+{
+ region2->xmin = region1->xmin;
+ region2->xmax = region1->xmax;
+ region2->ymin = region1->ymin;
+ region2->ymax = region1->ymax;
+
+ region2->free.left = region1->free.left;
+ region2->free.right = region1->free.right;
+ region2->free.top = region1->free.top;
+ region2->free.bottom = region1->free.bottom;
+}
+
+void reg_to_wdht(struct Region *region, int *wd, int *ht)
+{
+ *wd = region->xmax - region->xmin;
+ *ht = region->ymax - region->ymin;
+}
Property changes on: grass/trunk/imagery/i.smap/reg_util.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/region.h (from rev 34868, grass/trunk/imagery/i.smap/bouman/region.h)
===================================================================
--- grass/trunk/imagery/i.smap/region.h (rev 0)
+++ grass/trunk/imagery/i.smap/region.h 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,41 @@
+struct Region
+{
+ int xmin, xmax;
+ int ymin, ymax;
+ struct Free
+ {
+ int left;
+ int right;
+ int top;
+ int bottom;
+ } free;
+};
+
+/* decimate.c */
+void make_pyramid(LIKELIHOOD ****, struct Region *, int, double *, int);
+char ***get_pyramid(int, int, size_t);
+void free_pyramid(char *, int, int);
+char ****get_cubic_pyramid(int, int, int, size_t);
+void free_cubic_pyramid(char *, int, int, int);
+int levels(int, int);
+
+/* interp.c */
+void seq_MAP(unsigned char ***, struct Region *, LIKELIHOOD ****, int,
+ double *, int);
+void MLE(unsigned char **, LIKELIHOOD ***, struct Region *, int);
+
+/* reg_util.c */
+int levels_reg(struct Region *);
+void dec_reg(struct Region *);
+void copy_reg(struct Region *, struct Region *);
+void reg_to_wdht(struct Region *, int *, int *);
+
+/* read_block.c */
+int read_block(DCELL ***, struct Region *, struct files *);
+
+/* segment.c */
+#ifdef GRASS_IMAGERY_H
+/* model.c */
+void extract_init(struct SigSet *);
+void extract(DCELL ***, struct Region *, LIKELIHOOD ***, struct SigSet *);
+#endif
Property changes on: grass/trunk/imagery/i.smap/region.h
___________________________________________________________________
Name: svn:mime-type
+ text/x-chdr
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/segment.c (from rev 34870, grass/trunk/imagery/i.smap/bouman/segment.c)
===================================================================
--- grass/trunk/imagery/i.smap/segment.c (rev 0)
+++ grass/trunk/imagery/i.smap/segment.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,310 @@
+
+/****************************************************************************
+ *
+ * MODULE: i.smap segmentation
+ * AUTHOR(S): Charles Bouman,
+ * School of Electrical Engineering, Purdue University
+ * PURPOSE: image classification and segmentation
+ * COPYRIGHT: (C) 1999-2007 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 <stdlib.h>
+#include <unistd.h>
+#include <grass/imagery.h>
+#include <grass/glocale.h>
+#include "bouman.h"
+#include "region.h"
+
+
+static void init_reg(struct Region *, int, int, int);
+static int increment_reg(struct Region *, int, int, int);
+static int shift_img(DCELL ***, int, struct Region *, int);
+static int shift_ll(LIKELIHOOD ****, struct Region *, int);
+
+
+int segment(struct SigSet *S, /* class parameters */
+ struct parms *parms, struct files *files)
+{
+ int block_size; /* size of subregion blocks */
+ int ml; /* max likelihood? */
+ int quiet; /* be quiet when running? */
+
+ DCELL ***img; /* multispectral image, img[band][i][j] */
+ int last_row;
+ int wd, ht; /* image width and height */
+ struct Region region; /* specifies image subregion */
+ int nbands; /* number of bands */
+ int nclasses; /* number of classes */
+ LIKELIHOOD ****ll_pym; /* pyramid of log likelihoods */
+ unsigned char ***sf_pym; /* pyramid of segmentations */
+ int D; /* number of levels in pyramid */
+ double *alpha_dec; /* class transition probabilities */
+ int vlevel; /* level of verbose output */
+ int i;
+
+ quiet = parms->quiet; /* run quietly? */
+ ml = parms->ml; /* use maxl? */
+ block_size = parms->blocksize;
+
+ vlevel = quiet ? 0 : 1;
+ wd = G_window_cols(); /* get width from GRASS */
+ ht = G_window_rows(); /* get height from GRASS */
+
+ /* make blocksize a power of 2 */
+ if (block_size < 8)
+ block_size = 8;
+ for (i = 0; (block_size >> i) > 1; i++) {
+ }
+ block_size = 1 << i;
+
+/**** this code may stay the same ******/
+ nbands = S->nbands;
+ nclasses = S->nclasses;
+
+ /* Check for too many classes */
+ if (nclasses > 256)
+ G_fatal_error(_("Number of classes must be < 256"));
+
+ /* allocate alpha_dec parameters */
+ D = levels(block_size, block_size);
+ alpha_dec = (double *)G_malloc(D * sizeof(double));
+
+ /* allocate image block */
+ img =
+ (DCELL ***) multialloc(sizeof(DCELL), 3, nbands, block_size,
+ block_size);
+
+ /* allocate memory for log likelihood pyramid */
+ ll_pym =
+ (LIKELIHOOD ****) get_cubic_pyramid(block_size, block_size, nclasses,
+ sizeof(LIKELIHOOD));
+
+ /* allocate memory for segmentation pyramid */
+ sf_pym = (unsigned char ***)get_pyramid(wd, ht, sizeof(char));
+
+ /* tiled segmentation */
+ init_reg(®ion, wd, ht, block_size);
+ extract_init(S);
+ last_row = -1;
+ do {
+ if (vlevel >= 1 && last_row != region.ymin)
+ G_message(_("Processing rows %d-%d (of %d)..."),
+ region.ymin + 1, region.ymax, ht);
+ last_row = region.ymin;
+ shift_img(img, nbands, ®ion, block_size);
+ /* this reads grass images into the block defined in region */
+ read_block(img, ®ion, files);
+
+ shift_ll(ll_pym, ®ion, block_size);
+ extract(img, ®ion, ll_pym[0], S);
+
+ if (ml)
+ MLE(sf_pym[0], ll_pym[0], ®ion, nclasses);
+ else {
+ for (i = 0; i < D; i++)
+ alpha_dec[i] = 1.0;
+ seq_MAP(sf_pym, ®ion, ll_pym, nclasses, alpha_dec, vlevel);
+ }
+
+
+ } while (increment_reg(®ion, wd, ht, block_size));
+
+ write_img(sf_pym[0], wd, ht, S, parms, files);
+
+ return 0;
+}
+
+static void init_reg(struct Region *region, int wd, int ht, int block_size)
+{
+ region->xmin = 0;
+ region->ymin = 0;
+
+ region->xmax = block_size;
+ if (region->xmax > wd)
+ region->xmax = wd;
+
+ region->ymax = block_size;
+ if (region->ymax > ht)
+ region->ymax = ht;
+
+ region->free.left = 1;
+ region->free.top = 1;
+ region->free.right = 1;
+ region->free.bottom = 1;
+}
+
+static int increment_reg(struct Region *region, int wd, int ht,
+ int block_size)
+{
+
+ /* shift one block to right */
+ if (region->xmax < wd) {
+ region->xmin = region->xmax;
+ region->xmax = region->xmin + block_size;
+ if (region->xmax > wd)
+ region->xmax = wd;
+ }
+ else {
+ /* shift one block down */
+ if (region->ymax < ht) {
+ region->xmin = 0;
+ region->xmax = region->xmin + block_size;
+ if (region->xmax > wd)
+ region->xmax = wd;
+
+ region->ymin = region->ymax;
+ region->ymax = region->ymax + block_size;
+ if (region->ymax > ht)
+ region->ymax = ht;
+ }
+ /* finished shifting */
+ else {
+ return (0);
+ }
+ }
+
+ if (region->xmin == 0)
+ region->free.left = 1;
+ else
+ region->free.left = 0;
+
+ if (region->ymin == 0)
+ region->free.top = 1;
+ else
+ region->free.top = 0;
+
+ region->free.right = 1;
+ region->free.bottom = 1;
+
+ return (1);
+}
+
+static int shift_img(DCELL *** img, int nbands,
+ struct Region *region, int block_size)
+{
+ static int xoffset = 0;
+ static int yoffset = 0;
+ int xdelta;
+ int ystart, ystop, ydelta;
+ int b, i;
+
+ xdelta = region->xmin - xoffset;
+ ydelta = region->ymin - yoffset;
+ xoffset = region->xmin;
+ yoffset = region->ymin;
+
+ ystart = region->ymin;
+ ystop = ystart + block_size;
+
+ for (b = 0; b < nbands; b++) {
+ img[b] -= ydelta;
+ for (i = ystart; i < ystop; i++)
+ img[b][i] -= xdelta;
+ }
+
+ return 0;
+}
+
+static int shift_ll(LIKELIHOOD **** ll_pym,
+ struct Region *region, int block_size)
+{
+ static int first = 1;
+ static int xoffset[20];
+ static int yoffset[20];
+
+ int xdelta;
+ int ystart, ystop, ydelta;
+ int D;
+ int k, i;
+ int block_size_k;
+ struct Region region_buff;
+
+ /* if first time, set offsets to 0 */
+ if (first) {
+ D = levels(block_size, block_size);
+ for (k = 0; k <= D; k++)
+ xoffset[k] = yoffset[k] = 0;
+ first = 0;
+ }
+
+ /* save region information */
+ copy_reg(region, ®ion_buff);
+
+ /* subtract pointer offsets for each scale */
+ D = levels(block_size, block_size);
+ block_size_k = block_size;
+ for (k = 0; k <= D; k++) {
+ xdelta = region->xmin - xoffset[k];
+ ydelta = region->ymin - yoffset[k];
+ xoffset[k] = region->xmin;
+ yoffset[k] = region->ymin;
+
+ ystart = region->ymin;
+ ystop = ystart + block_size_k;
+
+ ll_pym[k] -= ydelta;
+ for (i = ystart; i < ystop; i++)
+ ll_pym[k][i] -= xdelta;
+ dec_reg(region);
+ block_size_k = block_size_k / 2;
+ }
+
+ /* replace region information */
+ copy_reg(®ion_buff, region);
+
+ return 0;
+}
+
+
+#ifdef COMMENTED_OUT
+int read_block(DCELL *** img, int wd, int ht, int nbands,
+ struct Region *region, char *infn)
+{
+ static first = 1;
+ static unsigned char ***img_buf;
+
+ int b, i, j;
+ int xstop, ystop;
+ char *infn_num;
+
+ /* allocate image buffer first time called */
+ if (first) {
+ img_buf =
+ (unsigned char ***)multialloc(sizeof(unsigned char), 3, nbands,
+ ht, wd);
+
+ /* allocate memory for name extension */
+ infn_num = (char *)G_malloc((strlen(optarg) + 10) * sizeof(char));
+
+ /* read data */
+ if (nbands == 1) {
+ read_img(img_buf[0], wd, ht, infn);
+ }
+ else {
+ for (b = 1; b <= nbands; b++) {
+ sprintf(infn_num, "%s.%d", infn, b);
+ read_img(img_buf[b - 1], wd, ht, infn_num);
+ }
+ }
+ first = 0;
+ }
+
+ /* copy data from buffer */
+ xstop = region->xmax;
+ if (xstop > wd)
+ xstop = wd;
+ ystop = region->ymax;
+ if (ystop > ht)
+ ystop = ht;
+ for (b = 0; b < nbands; b++)
+ for (i = region->ymin; i < ystop; i++)
+ for (j = region->xmin; j < xstop; j++)
+ img[b][i][j] = img_buf[b][i][j];
+
+ return 0;
+}
+#endif
Property changes on: grass/trunk/imagery/i.smap/segment.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/solve.c (from rev 34868, grass/trunk/imagery/i.smap/bouman/solve.c)
===================================================================
--- grass/trunk/imagery/i.smap/solve.c (rev 0)
+++ grass/trunk/imagery/i.smap/solve.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,61 @@
+
+double solve(
+ /* Solves equation (*f)(x) = 0 on x in [a,b]. Uses half interval method. */
+ /* Requires that (*f)(a) and (*f)(b) have opposite signs. */
+ /* Returns code=0 if signs are opposite. */
+ /* Returns code=1 if signs are both positive. */
+ /* Returns code=1 if signs are both negative. */
+ double (*f) (double), /* pointer to function to be solved */
+ double a, /* minimum value of solution */
+ double b, /* maximum value of solution */
+ double err, /* accuarcy of solution */
+ int *code /* error code */
+ )
+{
+ int signa, signb, signc;
+ double fa, fb, fc, c, signaling_nan();
+ double dist;
+
+ fa = (*f) (a);
+ signa = fa > 0;
+ fb = (*f) (b);
+ signb = fb > 0;
+
+ /* check starting conditions */
+ if (signa == signb) {
+ if (signa == 1)
+ *code = 1;
+ else
+ *code = -1;
+ return (0.0);
+ }
+ else
+ *code = 0;
+
+ /* half interval search */
+ if ((dist = b - a) < 0)
+ dist = -dist;
+ while (dist > err) {
+ c = (b + a) / 2;
+ fc = (*f) (c);
+ signc = fc > 0;
+ if (signa == signc) {
+ a = c;
+ fa = fc;
+ }
+ else {
+ b = c;
+ fb = fc;
+ }
+ if ((dist = b - a) < 0)
+ dist = -dist;
+ }
+
+ /* linear interpolation */
+ if ((fb - fa) == 0)
+ return (a);
+ else {
+ c = (a * fb - b * fa) / (fb - fa);
+ return (c);
+ }
+}
Property changes on: grass/trunk/imagery/i.smap/solve.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/trunk/imagery/i.smap/write_img.c (from rev 34870, grass/trunk/imagery/i.smap/shapiro/write_img.c)
===================================================================
--- grass/trunk/imagery/i.smap/write_img.c (rev 0)
+++ grass/trunk/imagery/i.smap/write_img.c 2008-12-14 14:29:09 UTC (rev 34871)
@@ -0,0 +1,27 @@
+#include <grass/imagery.h>
+#include <grass/glocale.h>
+#include "bouman.h"
+
+
+int write_img(unsigned char **img, 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;
+
+ G_message(_("Writing raster map <%s>..."), parms->output_map);
+
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 2);
+ for (col = 0; col < ncols; col++) {
+ int class = (int)img[row][col];
+
+ G_debug(3, "class: [%d] row/col: [%d][%d]", class, row, col);
+ files->outbuf[col] = (CELL) S->ClassSig[class].classnum;
+ }
+ G_put_raster_row(files->output_fd, files->outbuf, CELL_TYPE);
+ }
+ G_percent(nrows, nrows, 2);
+
+ return 0;
+}
Property changes on: grass/trunk/imagery/i.smap/write_img.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date Id
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
More information about the grass-commit
mailing list