diff --git a/imagery/i.cca/main.c b/imagery/i.cca/main.c index fbe2d11ae11..50deb3b628c 100644 --- a/imagery/i.cca/main.c +++ b/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); diff --git a/imagery/i.cca/stats.c b/imagery/i.cca/stats.c index 08f4e680751..96873bbf043 100644 --- a/imagery/i.cca/stats.c +++ b/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); diff --git a/imagery/i.cca/transform.c b/imagery/i.cca/transform.c index 0c3607215fe..a753c5c36cc 100644 --- a/imagery/i.cca/transform.c +++ b/imagery/i.cca/transform.c @@ -1,4 +1,5 @@ #include +#include #include #include @@ -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; }