i.cca: fix indexing, buffer size issues and matrix computation

imagery/i.cca/main.c

@@ -156,10 +156,10 @@ int main(int argc, char *argv[])
         cov[i] = G_alloc_matrix(bands, bands);
     }
 
-    outbandmax = (CELL *)G_calloc(nclass, sizeof(CELL));
-    outbandmin = (CELL *)G_calloc(nclass, sizeof(CELL));
-    datafds = (int *)G_calloc(nclass, sizeof(int));
-    outfds = (int *)G_calloc(nclass, sizeof(int));
+    outbandmax = (CELL *)G_calloc(bands, sizeof(CELL));
+    outbandmin = (CELL *)G_calloc(bands, sizeof(CELL));
+    datafds = (int *)G_calloc(bands, sizeof(int));
+    outfds = (int *)G_calloc(bands, sizeof(int));
 
     /*
        Here is where the information regarding
@@ -169,13 +169,13 @@ int main(int argc, char *argv[])
      */
 
     samptot = 0;
-    for (i = 1; i <= nclass; i++) {
-        nsamp[i] = sigs.sig[i - 1].npoints;
+    for (i = 0; i < nclass; i++) {
+        nsamp[i] = sigs.sig[i].npoints;
         samptot += nsamp[i];
-        for (j = 1; j <= bands; j++) {
-            mu[i][j] = sigs.sig[i - 1].mean[j - 1];
-            for (k = 1; k <= j; k++)
-                cov[i][j][k] = cov[i][k][j] = sigs.sig[i - 1].var[j - 1][k - 1];
+        for (j = 0; j < bands; j++) {
+            mu[i][j] = sigs.sig[i].mean[j];
+            for (k = 0; k <= j; k++)
+                cov[i][j][k] = cov[i][k][j] = sigs.sig[i].var[j][k];
         }
     }
 
@@ -192,6 +192,18 @@ int main(int argc, char *argv[])
     G_math_egvorder(eigval, eigmat, bands);
     G_math_d_AB(eigmat, w, q, bands, bands, bands);
 
+    /**Normalize each row (eigenvector) of the transformation matrix q */
+    for (i = 0; i < bands; i++) {
+        double norm = 0.0;
+        for (j = 0; j < bands; j++)
+            norm += q[i][j] * q[i][j];
+        norm = sqrt(norm);
+        if (norm > 0.0) {
+            for (j = 0; j < bands; j++)
+                q[i][j] /= norm;
+        }
+    }
+
     for (i = 0; i < bands; i++) {
         G_verbose_message("%i. eigen value: %+6.5f", i, eigval[i]);
         G_verbose_message("eigen vector:");
@@ -200,15 +212,14 @@ int main(int argc, char *argv[])
     }
 
     /* open the cell maps */
-    for (i = 1; i <= bands; i++) {
+    for (i = 0; i < bands; i++) {
         outbandmax[i] = (CELL)0;
         outbandmin[i] = (CELL)0;
 
-        datafds[i] =
-            Rast_open_old(refs.file[i - 1].name, refs.file[i - 1].mapset);
+        datafds[i] = Rast_open_old(refs.file[i].name, refs.file[i].mapset);
 
-        snprintf(tempname, sizeof(tempname), "%s.%d", out_opt->answer, i);
-        outfds[i] = Rast_open_c_new(tempname);
+        snprintf(tempname, sizeof(tempname), "%s.%d", out_opt->answer, i + 1);
+        outfds[i] = Rast_open_new(tempname, DCELL_TYPE);
     }
 
     /* do the transform */
@@ -219,18 +230,18 @@ int main(int argc, char *argv[])
     Rast_init_colors(&color_tbl);
 
     /* close the cell maps */
-    for (i = 1; i <= bands; i++) {
+    for (i = 0; i < bands; i++) {
         Rast_close(datafds[i]);
         Rast_close(outfds[i]);
 
         if (outbandmin[i] < (CELL)0 || outbandmax[i] > (CELL)255) {
             G_warning(_("The output cell map <%s.%d> has values "
                         "outside the 0-255 range."),
-                      out_opt->answer, i);
+                      out_opt->answer, i + 1);
         }
 
         Rast_make_grey_scale_colors(&color_tbl, 0, outbandmax[i]);
-        snprintf(tempname, sizeof(tempname), "%s.%d", out_opt->answer, i);
+        snprintf(tempname, sizeof(tempname), "%s.%d", out_opt->answer, i + 1);
 
         /* write a color table */
         Rast_write_colors(tempname, G_mapset(), &color_tbl);

imagery/i.cca/stats.c

@@ -20,7 +20,7 @@ int within(int samptot, int nclass, double *nsamp, double ***cov, double **w,
 
     for (i = 0; i < bands; i++)
         for (j = 0; j < bands; j++)
-            w[i][j] = (1.0 / ((double)(samptot - nclass))) * w[i][j];
+            w[i][j] /= (samptot - nclass);
 
     return 0;
 }
@@ -37,20 +37,43 @@ int between(int samptot, int nclass, double *nsamp, double **mu, double **p,
     tmp2 = G_alloc_matrix(bands, bands);
     newvec = G_alloc_vector(bands);
 
+    /* Initialize to zero */
+    for (i = 0; i < bands; i++) {
+        newvec[i] = 0.0;
+        for (j = 0; j < bands; j++) {
+            tmp1[i][j] = 0.0;
+            tmp2[i][j] = 0.0;
+        }
+    }
+
+    /* Compute weighted sum for overall mean vector */
     for (i = 0; i < nclass; i++)
         for (j = 0; j < bands; j++)
             newvec[j] += nsamp[i] * mu[i][j];
-    for (i = 0; i < bands; i++)
-        for (j = 0; j < bands; j++)
-            tmp1[i][j] = (newvec[i] * newvec[j]) / samptot;
 
+    /* Compute N * mu_bar * mu_bar^T */
+    for (i = 0; i < bands; i++) {
+        double mu_bar_i = newvec[i] / samptot;
+        for (j = 0; j < bands; j++) {
+            double mu_bar_j = newvec[j] / samptot;
+            tmp1[i][j] = samptot * mu_bar_i * mu_bar_j;
+        }
+    }
+
+    /* Compute sum of n_i * mu_i * mu_i^T */
     for (k = 0; k < nclass; k++) {
         product(mu[k], nsamp[k], tmp0, bands);
         for (i = 0; i < bands; i++)
             for (j = 0; j < bands; j++)
-                tmp2[i][j] += tmp0[i][j] - tmp1[i][j];
+                tmp2[i][j] += tmp0[i][j];
     }
 
+    /* Subtract overall mean component */
+    for (i = 0; i < bands; i++)
+        for (j = 0; j < bands; j++)
+            tmp2[i][j] -= tmp1[i][j];
+
+    /* Normalize */
     for (i = 0; i < bands; i++)
         for (j = 0; j < bands; j++)
             p[i][j] = tmp2[i][j] / (nclass - 1);

imagery/i.cca/transform.c

@@ -1,4 +1,5 @@
 #include <stdlib.h>
+#include <limits.h>
 
 #include <grass/gis.h>
 #include <grass/gmath.h>
@@ -11,49 +12,60 @@ int transform(int *datafds, int *outfds, int rows, int cols, double **eigmat,
 {
     int i, j, k, l;
     double *sum;
-    CELL **rowbufs;
+    DCELL **rowbufs;
 
     sum = G_alloc_vector(bands);
-    rowbufs = (CELL **)G_calloc(bands, sizeof(CELL *));
+    rowbufs = (DCELL **)G_calloc(bands, sizeof(DCELL *));
 
     /* allocate row buffers for each band */
-    for (i = 0; i < bands; i++)
-        if ((rowbufs[i] = Rast_allocate_c_buf()) == NULL)
-            G_fatal_error(_("Unable to allocate cell buffers."));
+    for (i = 0; i < bands; i++) {
+        rowbufs[i] = Rast_allocate_d_buf();
+        if (!rowbufs[i])
+            G_fatal_error(_("Unable to allocate DCELL buffers."));
+        mins[i] = INT_MAX;
+        maxs[i] = INT_MIN;
+    }
 
     for (i = 0; i < rows; i++) {
-        /* get one row of data */
-        for (j = 0; j < bands; j++)
-            Rast_get_c_row(datafds[j], rowbufs[j], i);
+        /* read each input band as CELL, convert to DCELL */
+        for (j = 0; j < bands; j++) {
+            CELL *tmp = Rast_allocate_c_buf();
+            Rast_get_c_row(datafds[j], tmp, i);
+            for (l = 0; l < cols; l++)
+                rowbufs[j][l] = (DCELL)tmp[l];
+            G_free(tmp);
+        }
 
-        /* transform each cell in the row */
+        /* apply canonical transform */
         for (l = 0; l < cols; l++) {
             for (j = 0; j < bands; j++) {
                 sum[j] = 0.0;
                 for (k = 0; k < bands; k++) {
-                    sum[j] += eigmat[j][k] * (double)rowbufs[k][l];
+                    sum[j] += eigmat[j][k] * rowbufs[k][l];
                 }
             }
             for (j = 0; j < bands; j++) {
-                rowbufs[j][l] = (CELL)(sum[j] + 0.5);
-                if (rowbufs[j][l] > maxs[j])
-                    maxs[j] = rowbufs[j][l];
-                if (rowbufs[j][l] < mins[j])
-                    mins[j] = rowbufs[j][l];
+                rowbufs[j][l] = sum[j];
+                /*min/max are tracked as integer CELL */
+                if ((CELL)rowbufs[j][l] > maxs[j])
+                    maxs[j] = (CELL)rowbufs[j][l];
+                if ((CELL)rowbufs[j][l] < mins[j])
+                    mins[j] = (CELL)rowbufs[j][l];
             }
         }
 
-        /* output the row of data */
+        /* write each output band row as DCELL */
         for (j = 0; j < bands; j++)
-            Rast_put_row(outfds[j], rowbufs[j], CELL_TYPE);
+            Rast_put_d_row(outfds[j], rowbufs[j]);
     }
+
     for (i = 0; i < bands; i++)
         G_free(rowbufs[i]);
 
     G_free(rowbufs);
     G_free_vector(sum);
 
-    G_message(_("Transform completed.\n"));
+    G_message(_("Canonical transform completed."));
 
     return 0;
 }