[GRASS-SVN] r51647 - grass/branches/develbranch_6/imagery/i.pca
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri May 18 04:16:08 EDT 2012
Author: mmetz
Date: 2012-05-18 01:16:07 -0700 (Fri, 18 May 2012)
New Revision: 51647
Modified:
grass/branches/develbranch_6/imagery/i.pca/description.html
grass/branches/develbranch_6/imagery/i.pca/local_proto.h
grass/branches/develbranch_6/imagery/i.pca/main.c
grass/branches/develbranch_6/imagery/i.pca/support.c
Log:
i.pca backport from trunk
Modified: grass/branches/develbranch_6/imagery/i.pca/description.html
===================================================================
--- grass/branches/develbranch_6/imagery/i.pca/description.html 2012-05-18 02:16:57 UTC (rev 51646)
+++ grass/branches/develbranch_6/imagery/i.pca/description.html 2012-05-18 08:16:07 UTC (rev 51647)
@@ -1,6 +1,5 @@
<h2>DESCRIPTION</h2>
-
<em>i.pca</em> is an image processing program based on the
algorithm provided by Vali (1990), that processes n
(n >= 2) input raster map layers and produces n output
@@ -33,12 +32,12 @@
<dd>If output is rescaled, the output raster will be of type CELL. If
the output is not rescaled, the output raster will be of type DCELL.
-<dt><B>-n</B> <EM>normalize input raster maps</EM>
-<dd>By default, the values of the input raster maps are centered for each
-map separately with <EM>x - mean</EM>. With <EM>-n</EM>, the input raster
-maps are normalized for each map separately with <EM>(x - mean) / stddev</EM>.
+<dt><b>-n</b> <em>normalize input raster maps</em>
+<dd>By default, the values of the input raster maps are centered for each
+map separately with <em>x - mean</em>. With <em>-n</em>, the input raster
+maps are normalized for each map separately with <em>(x - mean) / stddev</em>.
Normalizing is highly recommended when the input raster maps have different
-units, e.g. represent different environmental parameters.
+units, e.g. represent different environmental parameters.
</dl>
@@ -47,9 +46,8 @@
Richards (1986) gives a good example of the application of principal
components analysis (pca) to a time series of LANDSAT images of a burned
region in Australia.
-<p>
-Eigenvalue and eigenvector information is stored in the output maps'
-history files. View with <em>r.info -h</em>.
+<p>Eigenvalue and eigenvector information is stored in the output maps'
+history files. View with <em>r.info</em>.
<h2>EXAMPLE</h2>
@@ -74,14 +72,12 @@
Springer-Verlag, 1986.
<p>
-
Vali, Ali R.,
Personal communication,
Space Research Center,
University of Texas, Austin, 1990.
<p>
-
<em>
<a href="i.cca.html">i.cca</a><br>
<a href="i.class.html">i.class</a><br>
@@ -96,15 +92,12 @@
<h2>AUTHOR</h2>
David Satnik, GIS Laboratory
-<p>
-Major modifications for GRASS 4.1 were made by <br>
+<p>Major modifications for GRASS 4.1 were made by <br>
Olga Waupotitsch and
Michael Shapiro,
U.S.Army Construction Engineering
Research Laboratory
-<p>
-Rewritten for GRASS 6.x and major modifications by <br>
+<p>Rewritten for GRASS 6.x and major modifications by <br>
Brad Douglas
-<p>
-<i>Last changed: $Date$</i>
+<p><i>Last changed: $Date$</i>
Modified: grass/branches/develbranch_6/imagery/i.pca/local_proto.h
===================================================================
--- grass/branches/develbranch_6/imagery/i.pca/local_proto.h 2012-05-18 02:16:57 UTC (rev 51646)
+++ grass/branches/develbranch_6/imagery/i.pca/local_proto.h 2012-05-18 08:16:07 UTC (rev 51647)
@@ -2,6 +2,6 @@
#define __LOCAL_PROTO_H__
/* support.c */
-int write_support(int, char *, double **, double *);
+void write_support(int, char *, double **, double *);
#endif
Modified: grass/branches/develbranch_6/imagery/i.pca/main.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.pca/main.c 2012-05-18 02:16:57 UTC (rev 51646)
+++ grass/branches/develbranch_6/imagery/i.pca/main.c 2012-05-18 08:16:07 UTC (rev 51647)
@@ -5,10 +5,11 @@
*
* AUTHOR(S): Original author Center for Space Research (Uni. of TX)
* Rewritten by Brad Douglas <rez touchofmadness com>
+ * NULL value/MASK handling and speed up by Markus Metz
*
- * PURPOSE: Principal Component Analysis transform of satellite data.
+ * PURPOSE: Principal Component Analysis transform of raster data.
*
- * COPYRIGHT: (C) 2004-2008 by the GRASS Development Team
+ * COPYRIGHT: (C) 2004-2011 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
@@ -33,8 +34,7 @@
/* function prototypes */
static CELL round_c(double);
static int set_output_scale(struct Option *, int *, int *, int *);
-static int calc_mu(int *, double *, double *, int);
-static int calc_covariance(int *, double **, double *, double *, int);
+static int calc_mu_cov(int *, double **, double *, double *, int);
static int write_pca(double **, double *, double *, int *, char *, int,
int, int, int);
@@ -142,12 +142,9 @@
opt_in->answers[i]);
}
- G_verbose_message(_("Calculating covariance matrix..."));
- if (!calc_mu(inp_fd, mu, stddev, bands))
+ if (!calc_mu_cov(inp_fd, covar, mu, stddev, bands))
G_fatal_error(_("No non-null values"));
- calc_covariance(inp_fd, covar, mu, stddev, bands);
-
G_math_d_copy(covar[0], eigmat[0], bands*bands);
G_debug(1, "Calculating eigenvalues and eigenvectors...");
G_math_eigen(eigmat, eigval, bands);
@@ -237,86 +234,32 @@
}
-static int calc_mu(int *fds, double *mu, double *stddev, int bands)
+static int calc_mu_cov(int *fds, double **covar, double *mu,
+ double *stddev, int bands)
{
- int i;
+ int i, j;
int row, col;
int rows = G_window_rows();
int cols = G_window_cols();
off_t count = 0;
- double *sumsq;
DCELL **rowbuf = (DCELL **) G_malloc(bands * sizeof(DCELL *));
+ double **sum2 = (double **)G_calloc(bands, sizeof(double *));
+ double *sumsq, *sd, *sum;
- for (i = 0; i < bands; i++) {
- rowbuf[i] = G_allocate_d_raster_buf();
- }
-
if (stddev) {
- G_message(_("Computing means and standard deviations..."));
sumsq = (double *)G_calloc(bands, sizeof(double));
+ sd = (double *)G_calloc(bands, sizeof(double));
}
else {
- G_message(_("Computing means..."));
sumsq = NULL;
+ sd = NULL;
}
- for (row = 0; row < rows; row++) {
- G_percent(row, rows, 2);
- for (i = 0; i < bands; i++)
- G_get_d_raster_row(fds[i], rowbuf[i], row);
-
- for (col = 0; col < cols; col++) {
- /* ignore cells where any of the maps has null value */
- for (i = 0; i < bands; i++)
- if (G_is_d_null_value(&rowbuf[i][col]))
- break;
- if (i != bands)
- continue;
- count++;
- for (i = 0; i < bands; i++) {
- mu[i] += rowbuf[i][col];
- if (stddev)
- sumsq[i] += rowbuf[i][col] * rowbuf[i][col];
- }
- }
- }
- G_percent(1, 1, 1);
-
- if (count < 2)
- return 0;
-
for (i = 0; i < bands; i++) {
- if (stddev)
- /* mu[i] is here still the sum */
- stddev[i] = sqrt((sumsq[i] - mu[i] * mu[i] / count) /
- (count - 1));
-
- mu[i] = mu[i] / count;
- G_free(rowbuf[i]);
- }
-
- if (rowbuf)
- G_free(rowbuf);
- if (sumsq)
- G_free(sumsq);
-
- return 1;
-}
-
-
-static int calc_covariance(int *fds, double **covar, double *mu,
- double *stddev, int bands)
-{
- int i, j;
- int row, col;
- int rows = G_window_rows();
- int cols = G_window_cols();
- off_t count = 0;
- DCELL **rowbuf = (DCELL **) G_malloc(bands * sizeof(DCELL *));
-
- for (i = 0; i < bands; i++) {
rowbuf[i] = G_allocate_d_raster_buf();
+ sum2[i] = (double *)G_calloc(bands, sizeof(double));
}
+ sum = mu;
G_message(_("Computing covariance matrix..."));
@@ -334,38 +277,54 @@
continue;
count++;
for (i = 0; i < bands; i++) {
- DCELL dval1 = rowbuf[i][col];
- for (j = i; j < bands; j++) {
- DCELL dval2 = rowbuf[j][col];
+ sum[i] += rowbuf[i][col];
+ if (stddev)
+ sumsq[i] += rowbuf[i][col] * rowbuf[i][col];
- if (stddev) {
- covar[i][j] += (dval1 - mu[i]) * (dval2 - mu[j]) /
- (stddev[i] * stddev[j]);
- }
- else {
- covar[i][j] += (dval1 - mu[i]) * (dval2 - mu[j]);
- }
-
- }
+ for (j = 0; j <= i; j++)
+ sum2[i][j] += rowbuf[i][col] * rowbuf[j][col];
}
}
}
G_percent(1, 1, 1);
+ if (count < 2)
+ return 0;
+
for (i = 0; i < bands; i++) {
- for (j = i; j < bands; j++) {
- covar[i][j] = covar[i][j] / (count - 1);
+ if (stddev) {
+ sd[i] = sqrt(count * sumsq[i] - sum[i] * sum[i]);
+ stddev[i] = sqrt((sumsq[i] - sum[i] * sum[i] / count) /
+ (count - 1));
+ }
+ for (j = 0; j <= i; j++) {
+ if (stddev)
+ covar[i][j] = (count * sum2[i][j] - sum[i] * sum[j]) /
+ (sd[i] * sd[j]);
+ else
+ covar[i][j] = (sum2[i][j] - sum[i] * sum[j] / count) /
+ (count - 1);
G_debug(3, "covar[%d][%d] = %f", i, j, covar[i][j]);
if (j != i)
covar[j][i] = covar[i][j];
}
+ G_free(sum2[i]);
G_free(rowbuf[i]);
}
+ for (i = 0; i < bands; i++)
+ mu[i] = sum[i] / count;
+
G_free(rowbuf);
+
+ G_free(sum2);
+ if (sd)
+ G_free(sd);
+ if (sumsq)
+ G_free(sumsq);
- return 0;
+ return 1;
}
@@ -375,165 +334,147 @@
int scale, int scale_min, int scale_max)
{
int i, j;
- void *outbuf, *outptr;
- double min = 0.;
- double max = 0.;
- double old_range = 0.;
+ void **outbuf = (void **) G_malloc(bands * sizeof(void *));
+ void **outptr = (void **) G_malloc(bands * sizeof(void *));
+ double *min = (double *) G_malloc(bands * sizeof(double));
+ double *max = (double *) G_malloc(bands * sizeof(double));
+ double *old_range = (double *) G_calloc(bands, sizeof(double));
double new_range = 0.;
+ int pass;
int rows = G_window_rows();
int cols = G_window_cols();
/* why CELL_TYPE when scaling output ? */
- int outmap_type = (scale) ? CELL_TYPE : DCELL_TYPE;
+ int outmap_type = (scale) ? CELL_TYPE : FCELL_TYPE;
int outcell_mapsiz = G_raster_size(outmap_type);
- DCELL *d_buf;
+ int *out_fd = (int *) G_malloc(bands * sizeof(int));
+ DCELL **inbuf = (DCELL **) G_malloc(bands * sizeof(DCELL *));
/* 2 passes for rescale. 1 pass for no rescale */
int PASSES = (scale) ? 2 : 1;
- /* temporary row storage */
- d_buf = G_allocate_d_raster_buf();
-
- /* allocate memory for output row buffer */
- outbuf = G_allocate_raster_buf(outmap_type);
-
- if (!outbuf)
- G_fatal_error(_("Unable to allocate memory for raster row"));
-
+ /* allocate memory for row buffers */
for (i = 0; i < bands; i++) {
- char name[100];
- int out_fd;
- int pass;
+ char name[GNAME_MAX];
+ /* open output raster maps */
sprintf(name, "%s.%d", out_basename, i + 1);
+ out_fd[i] = G_open_raster_new(name, outmap_type);
- G_message(_("Transforming <%s>..."), name);
+ inbuf[i] = G_allocate_d_raster_buf();
+ outbuf[i] = G_allocate_raster_buf(outmap_type);
+ min[i] = max[i] = old_range[i] = 0;
+ }
- /* open a new file for output */
- out_fd = G_open_raster_new(name, outmap_type);
+ for (pass = 1; pass <= PASSES; pass++) {
+ int row, col;
+ int first = 1;
- for (pass = 1; pass <= PASSES; pass++) {
- void *rowbuf = NULL;
- int row, col;
+ if (scale && (pass == PASSES)) {
+ G_message(_("Rescaling to range %d,%d..."),
+ scale_min, scale_max);
- if (scale && (pass == PASSES)) {
- G_message(_("Rescaling <%s> to range %d,%d..."),
- name, scale_min, scale_max);
+ for (i = 0; i < bands; i++)
+ old_range[i] = max[i] - min[i];
+ new_range = (double)(scale_max - scale_min);
+ }
+ else {
+ G_message(_("Calculating principal components..."));
+ }
- old_range = max - min;
- new_range = (double)(scale_max - scale_min);
- }
+ for (row = 0; row < rows; row++) {
- for (row = 0; row < rows; row++) {
- void *rowptr;
+ G_percent(row, rows, 2);
- G_percent(row, rows, 2);
+ for (i = 0; i < bands; i++) {
+ G_get_d_raster_row(inp_fd[i], inbuf[i], row);
+ outptr[i] = outbuf[i];
+ }
+ for (col = 0; col < cols; col++) {
+ /* ignore cells where any of the maps has null value */
+ for (i = 0; i < bands; i++)
+ if (G_is_d_null_value(&inbuf[i][col]))
+ break;
+
+ if (i != bands) {
+ for (i = 0; i < bands; i++) {
+ G_set_null_value(outptr[i], 1, outmap_type);
+ outptr[i] =
+ G_incr_void_ptr(outptr[i], outcell_mapsiz);
+ }
+ continue;
+ }
- /* reset d_buf */
- for (col = 0; col < cols; col++)
- d_buf[col] = 0.;
+ for (i = 0; i < bands; i++) {
+ DCELL dval = 0.;
- for (j = 0; j < bands; j++) {
- RASTER_MAP_TYPE maptype =
- G_get_raster_map_type(inp_fd[j]);
-
- /* don't assume each image is of the same type */
- if (rowbuf)
- G_free(rowbuf);
- if (!(rowbuf = G_allocate_raster_buf(maptype)))
- G_fatal_error(_("Unable allocate memory for row buffer"));
-
- if (G_get_raster_row(inp_fd[j], rowbuf, row, maptype) < 0)
- G_fatal_error(_("Unable to read raster map row %d"), row);
-
- rowptr = rowbuf;
- outptr = outbuf;
-
- /* add into the output cell eigmat[i][j] * corresp cell
- * of j-th band for current j */
- for (col = 0; col < cols; col++) {
- DCELL dval;
-
- /* handle null cells */
- if (G_is_null_value(rowptr, maptype)) {
- G_set_null_value(outptr, 1, outmap_type);
- outptr =
- G_incr_void_ptr(outptr, outcell_mapsiz);
-
- rowptr =
- G_incr_void_ptr(rowptr,
- G_raster_size(maptype));
- continue;
- }
-
+ for (j = 0; j < bands; j++) {
/* corresp. cell of j-th band */
- dval = G_get_raster_value_d(rowptr, maptype);
if (stddev)
- d_buf[col] += eigmat[i][j] * ((dval - mu[j]) / stddev[j]);
+ dval += eigmat[i][j] *
+ ((inbuf[j][col] - mu[j]) / stddev[j]);
else
- d_buf[col] += eigmat[i][j] * (dval - mu[j]);
+ dval += eigmat[i][j] * (inbuf[j][col] - mu[j]);
+ }
- /* the cell entry is complete */
- if (j == (bands - 1)) {
- if (scale && (pass == 1)) {
- if ((row == 0) && (col == 0))
- min = max = d_buf[0];
- if (d_buf[col] < min)
- min = d_buf[col];
+ /* the cell entry is complete */
+ if (scale && (pass == 1)) {
+ if (first)
+ min[i] = max[i] = dval;
+ if (dval < min[i])
+ min[i] = dval;
- if (d_buf[col] > max)
- max = d_buf[col];
- }
- else if (scale) {
+ if (dval > max[i])
+ max[i] = dval;
+ }
+ else if (scale) {
- if (min == max) {
- G_set_raster_value_c(outptr, 1,
- CELL_TYPE);
- }
- else {
- /* map data to 0, (new_range-1) and then adding new_min */
- CELL tmpcell =
- round_c((new_range *
- (d_buf[col] -
- min) / old_range) +
- scale_min);
-
- G_set_raster_value_c(outptr, tmpcell,
- outmap_type);
- }
- }
- else { /* (!scale) */
-
- G_set_raster_value_d(outptr, d_buf[col],
- outmap_type);
- }
+ if (min[i] == max[i]) {
+ G_set_raster_value_c(outptr[i], 1,
+ CELL_TYPE);
}
+ else {
+ /* map data to 0, (new_range-1) and then adding new_min */
+ CELL tmpcell =
+ round_c((new_range * (dval - min[i]) /
+ old_range[i]) + scale_min);
- outptr =
- G_incr_void_ptr(outptr, outcell_mapsiz);
-
- rowptr =
- G_incr_void_ptr(rowptr, G_raster_size(maptype));
+ G_set_raster_value_c(outptr[i], tmpcell,
+ outmap_type);
+ }
}
- } /* for j = 0 to bands */
+ else { /* (!scale) */
- if (pass == PASSES) {
- G_put_raster_row(out_fd, outbuf, outmap_type);
+ G_set_raster_value_d(outptr[i], dval,
+ outmap_type);
+ }
+ outptr[i] = G_incr_void_ptr(outptr[i], outcell_mapsiz);
}
+ first = 0;
}
+ if (pass == PASSES) {
+ for (i = 0; i < bands; i++)
+ G_put_raster_row(out_fd[i], outbuf[i], outmap_type);
+ }
+ }
+ G_percent(1, 1, 1);
- G_percent(1, 1, 1);
-
- /* close output file */
- if (pass == PASSES)
- G_close_cell(out_fd);
+ /* close output file */
+ if (pass == PASSES) {
+ for (i = 0; i < bands; i++) {
+ G_close_cell(out_fd[i]);
+ G_free(inbuf[i]);
+ G_free(outbuf[i]);
+ }
}
}
- if (d_buf)
- G_free(d_buf);
- if (outbuf)
- G_free(outbuf);
+ G_free(inbuf);
+ G_free(outbuf);
+ G_free(outptr);
+ G_free(min);
+ G_free(max);
+ G_free(old_range);
return 0;
}
Modified: grass/branches/develbranch_6/imagery/i.pca/support.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.pca/support.c 2012-05-18 02:16:57 UTC (rev 51646)
+++ grass/branches/develbranch_6/imagery/i.pca/support.c 2012-05-18 08:16:07 UTC (rev 51647)
@@ -1,15 +1,15 @@
#include <string.h>
#include <grass/gis.h>
#include <grass/glocale.h>
+#include "local_proto.h"
-
/* function prototypes */
-static int write_history(int, char *, double **, double *);
+static void write_history(int, char *, double **, double *);
-int write_support(int bands, char *outname, double **eigmat, double *eigval)
+void write_support(int bands, char *outname, double **eigmat, double *eigval)
{
- char *mapset = G_mapset();
+ const char *mapset = G_mapset();
struct Colors colors;
struct FPRange range;
DCELL min, max;
@@ -23,14 +23,13 @@
if (G_raster_map_is_fp(outname, mapset))
G_mark_colors_as_fp(&colors);
- if (G_write_colors(outname, mapset, &colors) < 0)
- G_message(_("Unable to write color table for raster map <%s>"), outname);
+ G_write_colors(outname, mapset, &colors);
- return write_history(bands, outname, eigmat, eigval);
+ write_history(bands, outname, eigmat, eigval);
}
-static int write_history(int bands, char *outname, double **eigmat, double *eigval)
+static void write_history(int bands, char *outname, double **eigmat, double *eigval)
{
int i, j;
static int first_map = TRUE; /* write to stderr? */
@@ -52,29 +51,28 @@
sprintf(tmpeigen, "PC%d %9.2f (", i+1, eigval[i]);
for (j = 0; j < bands; j++) {
sprintf(tmpa, "%7.4f", eigmat[i][j]);
- G_strcat(tmpeigen, tmpa);
+ strcat(tmpeigen, tmpa);
if (j < (bands - 1) ){
sprintf(tmpa, ",");
- G_strcat(tmpeigen, tmpa);
+ strcat(tmpeigen, tmpa);
}
}
- G_strcat(tmpeigen, ") ");
-
+ strcat(tmpeigen, ") ");
+
sprintf(tmpa, "[%5.2f%%]", eigval[i] * 100/eigval_total);
- G_strcat(tmpeigen, tmpa);
+ strcat(tmpeigen, tmpa);
sprintf(hist.edhist[i + 1], tmpeigen);
/* write eigen values to screen */
if(first_map)
- G_message("%s", tmpeigen);
+ fprintf(stdout, "%s\n", tmpeigen);
}
- hist.edlinecnt = i + 1;
G_command_history(&hist);
/* only write to stderr the first time (this fn runs for every output map) */
first_map = FALSE;
- return G_write_history(outname, &hist);
+ G_write_history(outname, &hist);
}
More information about the grass-commit
mailing list