[GRASS-SVN] r49080 - grass/trunk/imagery/i.pca

[svn_grass at osgeo.org]

Author: mmetz
Date: 2011-11-03 11:02:42 -0700 (Thu, 03 Nov 2011)
New Revision: 49080

Modified:
   grass/trunk/imagery/i.pca/main.c
Log:
NULL/MASK fix, code cleanup

Modified: grass/trunk/imagery/i.pca/main.c
===================================================================
--- grass/trunk/imagery/i.pca/main.c	2011-11-03 16:05:29 UTC (rev 49079)
+++ grass/trunk/imagery/i.pca/main.c	2011-11-03 18:02:42 UTC (rev 49080)
@@ -34,7 +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 *, int);
+static off_t calc_mu(int *, double *, int);
 static int calc_covariance(int *, double **, double *, int);
 static int write_pca(double **, int *, char *, int, int, int, int);
 
@@ -53,6 +53,7 @@
     double **eigmat;
     int *inp_fd;
     int scale, scale_max, scale_min;
+    off_t cells;
 
     struct GModule *module;
     struct Option *opt_in, *opt_out, *opt_scale;
@@ -128,7 +129,7 @@
     }
 
     G_verbose_message(_("Calculating covariance matrix..."));
-    calc_mu(inp_fd, mu, bands);
+    cells = calc_mu(inp_fd, mu, bands);
 
     calc_covariance(inp_fd, covar, mu, bands);
 
@@ -136,7 +137,7 @@
 	for (j = 0; j < bands; j++) {
 	    covar[i][j] =
 		covar[i][j] /
-		((double)((Rast_window_rows() * Rast_window_cols()) - 1));
+		(cells - 1);
 	    G_debug(3, "covar[%d][%d] = %f", i, j, covar[i][j]);
 	}
     }
@@ -230,17 +231,20 @@
 }
 
 
-static int calc_mu(int *fds, double *mu, int bands)
+static off_t calc_mu(int *fds, double *mu, int bands)
 {
     int i;
     int rows = Rast_window_rows();
     int cols = Rast_window_cols();
+    off_t cells;
     void *rowbuf = NULL;
 
     for (i = 0; i < bands; i++) {
 	RASTER_MAP_TYPE maptype;
 	int row, col;
 	double sum = 0.;
+	
+	cells = 0;
 
 	maptype = Rast_get_map_type(fds[i]);
 
@@ -266,17 +270,18 @@
 		}
 
 		sum += Rast_get_d_value(ptr, maptype);
+		cells++;
 		ptr = G_incr_void_ptr(ptr, Rast_cell_size(maptype));
 	    }
 	}
 
-	mu[i] = sum / (double)(rows * cols);
+	mu[i] = sum / cells;
     }
 
     if (rowbuf)
 	G_free(rowbuf);
 
-    return 0;
+    return cells;
 }
 
 
@@ -360,19 +365,19 @@
     double new_range = 0.;
     int rows = Rast_window_rows();
     int cols = Rast_window_cols();
-    int cell_mapsiz = Rast_cell_size(CELL_TYPE);
-    int dcell_mapsiz = Rast_cell_size(DCELL_TYPE);
+    /* why CELL_TYPE when scaling output ? */
+    int outmap_type = (scale) ? CELL_TYPE : DCELL_TYPE;
+    int outcell_mapsiz = Rast_cell_size(outmap_type);
     DCELL *d_buf;
 
     /* 2 passes for rescale.  1 pass for no rescale */
     int PASSES = (scale) ? 2 : 1;
 
     /* temporary row storage */
-    d_buf = (DCELL *) G_malloc(cols * sizeof(double));
+    d_buf = Rast_allocate_d_buf();
 
     /* allocate memory for output row buffer */
-    outbuf = (scale) ? Rast_allocate_buf(CELL_TYPE) :
-	Rast_allocate_buf(DCELL_TYPE);
+    outbuf = Rast_allocate_buf(outmap_type);
 
     if (!outbuf)
 	G_fatal_error(_("Unable to allocate memory for raster row"));
@@ -387,12 +392,7 @@
 	G_message(_("Transforming <%s>..."), name);
 
 	/* open a new file for output */
-	if (scale)
-	    out_fd = Rast_open_c_new(name);
-	else {
-	    out_fd = Rast_open_fp_new(name);
-	    Rast_set_fp_type(DCELL_TYPE);
-	}
+	out_fd = Rast_open_new(name, outmap_type);
 
 	for (pass = 1; pass <= PASSES; pass++) {
 	    void *rowbuf = NULL;
@@ -435,15 +435,9 @@
 		    for (col = 0; col < cols; col++) {
 			/* handle null cells */
 			if (Rast_is_null_value(rowptr, maptype)) {
-			    if (scale) {
-				Rast_set_null_value(outptr, 1, CELL_TYPE);
-				outptr = G_incr_void_ptr(outptr, cell_mapsiz);
-			    }
-			    else {
-				Rast_set_null_value(outptr, 1, DCELL_TYPE);
-				outptr =
-				    G_incr_void_ptr(outptr, dcell_mapsiz);
-			    }
+			    Rast_set_null_value(outptr, 1, outmap_type);
+			    outptr =
+				G_incr_void_ptr(outptr, outcell_mapsiz);
 
 			    rowptr =
 				G_incr_void_ptr(rowptr,
@@ -483,19 +477,18 @@
 						scale_min);
 
 				    Rast_set_c_value(outptr, tmpcell,
-							 CELL_TYPE);
+							 outmap_type);
 				}
 			    }
 			    else {	/* (!scale) */
 
 				Rast_set_d_value(outptr, d_buf[col],
-						     DCELL_TYPE);
+						     outmap_type);
 			    }
 			}
 
-			outptr = (scale) ?
-			    G_incr_void_ptr(outptr, cell_mapsiz) :
-			    G_incr_void_ptr(outptr, dcell_mapsiz);
+			outptr = 
+			    G_incr_void_ptr(outptr, outcell_mapsiz);
 
 			rowptr =
 			    G_incr_void_ptr(rowptr, Rast_cell_size(maptype));
@@ -503,14 +496,11 @@
 		}		/* for j = 0 to bands */
 
 		if (pass == PASSES) {
-		    if (scale)
-			Rast_put_row(out_fd, outbuf, CELL_TYPE);
-		    else
-			Rast_put_row(out_fd, outbuf, DCELL_TYPE);
+		    Rast_put_row(out_fd, outbuf, outmap_type);
 		}
 	    }
 
-	    G_percent(row, rows, 2);
+	    G_percent(1, 1, 1);
 
 	    /* close output file */
 	    if (pass == PASSES)