1 #include <grass/gis.h>
2 #include <grass/gmath.h>
3
4 #include "local_proto.h"
5
6 int
within(int samptot,int nclass,double * nsamp,double *** cov,double ** w,int bands)7 within(int samptot, int nclass, double *nsamp, double ***cov,
8 double **w, int bands)
9 {
10 int i, j, k;
11
12 /* Initialize within class covariance matrix */
13 for (i = 0; i < bands; i++)
14 for (j = 0; j < bands; j++)
15 w[i][j] = 0.0;
16
17 for (i = 0; i < nclass; i++)
18 for (j = 0; j < bands; j++)
19 for (k = 0; k < bands; k++)
20 w[j][k] += (nsamp[i] - 1) * cov[i][j][k];
21
22 for (i = 0; i < bands; i++)
23 for (j = 0; j < bands; j++)
24 w[i][j] = (1.0 / ((double)(samptot - nclass))) * w[i][j];
25
26 return 0;
27 }
28
29
30 int
between(int samptot,int nclass,double * nsamp,double ** mu,double ** p,int bands)31 between(int samptot, int nclass, double *nsamp, double **mu,
32 double **p, int bands)
33 {
34 int i, j, k;
35 double **tmp0, **tmp1, **tmp2;
36 double *newvec;
37
38 tmp0 = G_alloc_matrix(bands, bands);
39 tmp1 = G_alloc_matrix(bands, bands);
40 tmp2 = G_alloc_matrix(bands, bands);
41 newvec = G_alloc_vector(bands);
42
43 for (i = 0; i < nclass; i++)
44 for (j = 0; j < bands; j++)
45 newvec[j] += nsamp[i] * mu[i][j];
46 for (i = 0; i < bands; i++)
47 for (j = 0; j < bands; j++)
48 tmp1[i][j] = (newvec[i] * newvec[j]) / samptot;
49
50 for (k = 0; k < nclass; k++) {
51 product(mu[k], nsamp[k], tmp0, bands);
52 for (i = 0; i < bands; i++)
53 for (j = 0; j < bands; j++)
54 tmp2[i][j] += tmp0[i][j] - tmp1[i][j];
55 }
56
57 for (i = 0; i < bands; i++)
58 for (j = 0; j < bands; j++)
59 p[i][j] = tmp2[i][j] / (nclass - 1);
60
61 G_free_matrix(tmp0);
62 G_free_matrix(tmp1);
63 G_free_matrix(tmp2);
64 G_free_vector(newvec);
65
66 return 0;
67 }
68
69