1 #include <stdio.h>
2 #include <math.h>
3 #include <stdlib.h>
4 #include <unistd.h>
5 #include <grass/gis.h>
6 #include <grass/raster.h>
7 #include <grass/glocale.h>
8 #include "method.h"
9
10 #define MEM 1024
11
12 /* function prototypes */
13 static int m_var(double *, int, double *);
14
15
o_var(const char * basemap,const char * covermap,const char * outputmap,int usecats,struct Categories * cats)16 int o_var(const char *basemap, const char *covermap, const char *outputmap,
17 int usecats, struct Categories *cats)
18 {
19 struct Popen stats_child, reclass_child;
20 FILE *stats, *reclass;
21 int first, mem, i, count;
22 long basecat, covercat, catb, catc;
23 double value, vari, x;
24 double *tab;
25
26
27 mem = MEM * sizeof(double);
28 tab = (double *)G_malloc(mem);
29
30 stats = run_stats(&stats_child, basemap, covermap, "-cn");
31 reclass = run_reclass(&reclass_child, basemap, outputmap);
32
33 first = 1;
34 while (read_stats(stats, &basecat, &covercat, &value)) {
35 if (first) {
36 first = 0;
37 catb = basecat;
38 catc = covercat;
39 i = 0;
40 count = 0;
41 }
42
43 if (basecat != catb) {
44 m_var(tab, count, &vari);
45 fprintf(reclass, "%ld = %ld %f\n", catb, catb, vari);
46 /*fprintf (stdout, "1. %ld = %ld %f\n", catb, catb, vari); */
47 catb = basecat;
48 catc = covercat;
49 count = 0;
50 }
51
52 if (usecats)
53 sscanf(Rast_get_c_cat((CELL *) &covercat, cats), "%lf", &x);
54 else
55 x = covercat;
56
57 for (i = 0; i < value; i++) {
58 if (count * sizeof(double) >= mem) {
59 mem += MEM * sizeof(double);
60 tab = (double *)G_realloc(tab, mem);
61 /* fprintf(stderr,"MALLOC: %d KB needed\n",(int)(mem/1024)); */
62 }
63 tab[count++] = x;
64 }
65
66 }
67
68 if (first) {
69 catb = catc = 0;
70 }
71
72 m_var(tab, count, &vari);
73 fprintf(reclass, "%ld = %ld %f\n", catb, catb, vari);
74 G_debug(5, "2. %ld = %ld %f", catb, catb, vari);
75
76 G_popen_close(&stats_child);
77 G_popen_close(&reclass_child);
78
79 return 0;
80 }
81
82
83 /***********************************************************************
84 *
85 * Given an array of data[1...n], this routine returns its variance.
86 *
87 ************************************************************************/
88
m_var(double * data,int n,double * vari)89 static int m_var(double *data, int n, double *vari)
90 {
91 double ave, ep, s;
92 int i;
93
94 if (n < 1) {
95 G_warning(_("o_var: No data in array"));
96 return (1);
97 }
98
99 *vari = 0.0;
100 ep = 0;
101 s = 0.0;
102
103
104 for (i = 0; i < n; i++) /* First pass to get the mean */
105 s += data[i];
106 ave = s / n;
107
108 for (i = 0; i < n; i++) {
109 s = data[i] - ave;
110 /*fprintf(stderr,"s: %lf data[i]: %lf ave: %lf n: %d\n",s,data[i],ave,n); */
111 *vari += s * s;
112 ep += s;
113 }
114
115 *vari = (*vari - ep * ep / n) / (n - 1);
116
117 return (0);
118 }
119