1 #include <stdlib.h>
2 #include "global.h"
3
4 /* hash definitions (these should be prime numbers) ************* */
5 #define HASHSIZE 7307
6 #define HASHMOD 89
7
8 static CELL *values;
9 static struct Node *node_pool;
10 static int node_pool_count;
11 static CELL *value_pool;
12 static int value_pool_count;
13
14 #define NODE_INCR 32
15 #define VALUE_INCR 32
16
17 static struct Node **sorted_list;
18
19 struct Node
20 {
21 CELL *values;
22 struct Node *left;
23 struct Node *right;
24 struct Node *list;
25 long count;
26 double area;
27 };
28
29 static struct Node **hashtable;
30 static struct Node *node_list = NULL;
31 static int node_count = 0, total_count = 0;
32
initialize_cell_stats(int n)33 int initialize_cell_stats(int n)
34 {
35 int i;
36
37 /* record nilfes first */
38 nfiles = n;
39
40 /* allocate a pool of value arrays */
41 value_pool_count = 0;
42 allocate_values();
43
44 /* set Node pool to empty */
45 node_pool_count = 0;
46
47
48 /* empty the has table */
49 hashtable = (struct Node **)G_malloc(HASHSIZE * sizeof(struct Node *));
50 for (i = 0; i < HASHSIZE; i++)
51 hashtable[i] = NULL;
52
53 return 0;
54 }
55
allocate_values(void)56 int allocate_values(void)
57 {
58 value_pool_count = VALUE_INCR;
59 value_pool = (CELL *) G_calloc(nfiles * value_pool_count, sizeof(CELL));
60 values = value_pool;
61
62 return 0;
63 }
64
NewNode(double area)65 struct Node *NewNode(double area)
66 {
67 struct Node *node;
68
69 if (node_pool_count <= 0)
70 node_pool = (struct Node *)G_calloc(node_pool_count =
71 NODE_INCR, sizeof(struct Node));
72 node = &node_pool[--node_pool_count];
73 node->count = 1;
74 node->area = area;
75 node->values = values;
76
77 if (--value_pool_count <= 0)
78 allocate_values();
79 else
80 values += nfiles;
81
82 node->left = node->right = NULL;
83 node->list = node_list;
84 node_list = node;
85 node_count++;
86
87 return node;
88 }
89
90
91 /* Essentially, Rast_quant_add_rule() treats the ranges as half-open,
92 * i.e. the values range from low (inclusive) to high (exclusive).
93 * While half-open ranges are a common concept (e.g. floor() behaves
94 * the same way), the range of a GRASS raster is closed, i.e. both the
95 * low and high values are inclusive.
96 * Therefore the quantized max FP cell gets put in the nsteps+1'th bin
97 * and we need to manually place it back in the previous bin. */
fix_max_fp_val(CELL * cell,int ncols)98 void fix_max_fp_val(CELL *cell, int ncols)
99 {
100 while (ncols-- > 0) {
101 if (cell[ncols] > nsteps)
102 cell[ncols] = (CELL)nsteps;
103 /* { G_debug(5, ". resetting %d to %d", cell[ncols], nsteps); } */
104 }
105 return;
106 }
107
108
109 /* we can't compute hash on null values, so we change all
110 * nulls to max+1, set NULL_CELL to max+1, and later compare
111 * with NULL_CELL to chack for nulls */
reset_null_vals(CELL * cell,int ncols)112 void reset_null_vals(CELL *cell, int ncols)
113 {
114 while (ncols-- > 0) {
115 if (Rast_is_c_null_value(&cell[ncols]))
116 cell[ncols] = NULL_CELL;
117 }
118 return;
119 }
120
121
update_cell_stats(CELL ** cell,int ncols,double area)122 int update_cell_stats(CELL ** cell, int ncols, double area)
123 {
124 register int i;
125 register int hash;
126 register struct Node *q, *p = NULL;
127 register int dir = 0;
128
129 while (ncols-- > 0) {
130 /* copy this cell to an array, compute hash */
131
132 hash = values[0] = cell[0][ncols];
133
134 for (i = 1; i < nfiles; i++)
135 hash = hash * HASHMOD + (values[i] = cell[i][ncols]);
136
137 if (hash < 0)
138 hash = -hash;
139
140 hash %= HASHSIZE;
141
142 /* look it up and update/insert */
143 if ((q = hashtable[hash]) == NULL) {
144 hashtable[hash] = NewNode(area);
145 }
146 else {
147 while (1) {
148 for (i = 0; i < nfiles; i++) {
149 if (values[i] < q->values[i]) {
150 dir = -1;
151 p = q->left;
152 break;
153 }
154 if (values[i] > q->values[i]) {
155 dir = 1;
156 p = q->right;
157 break;
158 }
159 }
160
161 if (i == nfiles) { /* match */
162 q->count++;
163 q->area += area;
164 total_count++;
165 break;
166 }
167 else if (p == NULL) {
168 if (dir < 0)
169 q->left = NewNode(area);
170 else
171 q->right = NewNode(area);
172 break;
173 }
174 else
175 q = p;
176 }
177 }
178 }
179
180 return 0;
181 }
182
node_compare(const void * pp,const void * qq)183 static int node_compare(const void *pp, const void *qq)
184 {
185 struct Node *const *p = pp, *const *q = qq;
186 register int i;
187 register const CELL *a, *b;
188
189 a = (*p)->values;
190 b = (*q)->values;
191 for (i = nfiles; --i >= 0;) {
192 if (*a < *b)
193 return -1;
194 else if (*a > *b)
195 return 1;
196 a++, b++;
197 }
198
199 return 0;
200 }
201
node_compare_count_asc(const void * pp,const void * qq)202 static int node_compare_count_asc(const void *pp, const void *qq)
203 {
204 struct Node *const *p = pp, *const *q = qq;
205 long a, b;
206
207 a = (*p)->count;
208 b = (*q)->count;
209
210 if (a < b)
211 return -1;
212 return (a > b);
213 }
214
node_compare_count_desc(const void * pp,const void * qq)215 static int node_compare_count_desc(const void *pp, const void *qq)
216 {
217 struct Node *const *p = pp, *const *q = qq;
218 long a, b;
219
220 a = (*p)->count;
221 b = (*q)->count;
222
223 if (a > b)
224 return -1;
225 return (a < b);
226 }
227
sort_cell_stats(int do_sort)228 int sort_cell_stats(int do_sort)
229 {
230 struct Node **q, *p;
231
232 if (node_count <= 0)
233 return 0;
234
235 G_free(hashtable); /* make a bit more room */
236 sorted_list = (struct Node **)G_calloc(node_count, sizeof(struct Node *));
237 for (q = sorted_list, p = node_list; p; p = p->list)
238 *q++ = p;
239
240 if (do_sort == SORT_DEFAULT)
241 qsort(sorted_list, node_count, sizeof(struct Node *), node_compare);
242 else if (do_sort == SORT_ASC)
243 qsort(sorted_list, node_count, sizeof(struct Node *), node_compare_count_asc);
244 else if (do_sort == SORT_DESC)
245 qsort(sorted_list, node_count, sizeof(struct Node *), node_compare_count_desc);
246
247 return 0;
248 }
249
print_node_count(void)250 int print_node_count(void)
251 {
252 fprintf(stdout, "%d nodes\n", node_count);
253
254 return 0;
255 }
256
print_cell_stats(char * fmt,int with_percents,int with_counts,int with_areas,int with_labels,char * fs)257 int print_cell_stats(char *fmt, int with_percents, int with_counts,
258 int with_areas, int with_labels, char *fs)
259 {
260 int i, n, nulls_found;
261 struct Node *node;
262 CELL tmp_cell, null_cell;
263 DCELL dLow, dHigh;
264 char str1[50], str2[50];
265
266 if (no_nulls)
267 total_count -= sorted_list[node_count - 1]->count;
268
269 Rast_set_c_null_value(&null_cell, 1);
270 if (node_count <= 0) {
271 fprintf(stdout, "0");
272 for (i = 1; i < nfiles; i++)
273 fprintf(stdout, "%s%s", fs, no_data_str);
274 if (with_areas)
275 fprintf(stdout, "%s0.0", fs);
276 if (with_counts)
277 fprintf(stdout, "%s0", fs);
278 if (with_percents)
279 fprintf(stdout, "%s0.00%%", fs);
280 if (with_labels)
281 fprintf(stdout, "%s%s", fs, Rast_get_c_cat(&null_cell, &labels[i]));
282 fprintf(stdout, "\n");
283 }
284 else {
285 for (n = 0; n < node_count; n++) {
286 node = sorted_list[n];
287
288 if (no_nulls || no_nulls_all) {
289 nulls_found = 0;
290 for (i = 0; i < nfiles; i++)
291 /*
292 if (node->values[i] || (!raw_output && is_fp[i]))
293 break;
294 */
295 if (node->values[i] == NULL_CELL)
296 nulls_found++;
297
298 if (nulls_found == nfiles)
299 continue;
300
301 if (no_nulls && nulls_found)
302 continue;
303 }
304
305 for (i = 0; i < nfiles; i++) {
306 if (node->values[i] == NULL_CELL) {
307 fprintf(stdout, "%s%s", i ? fs : "", no_data_str);
308 if (with_labels && !(raw_output && is_fp[i]))
309 fprintf(stdout, "%s%s", fs,
310 Rast_get_c_cat(&null_cell, &labels[i]));
311 }
312 else if (raw_output || !is_fp[i] || as_int) {
313 fprintf(stdout, "%s%ld", i ? fs : "",
314 (long)node->values[i]);
315 if (with_labels && !is_fp[i])
316 fprintf(stdout, "%s%s", fs,
317 Rast_get_c_cat((CELL*) &(node->values[i]),
318 &labels[i]));
319 }
320 else { /* find out which floating point range to print */
321
322 if (cat_ranges)
323 Rast_quant_get_ith_rule(&labels[i].q, node->values[i],
324 &dLow, &dHigh, &tmp_cell,
325 &tmp_cell);
326 else {
327 dLow = (DMAX[i] - DMIN[i]) / nsteps *
328 (double)(node->values[i] - 1) + DMIN[i];
329 dHigh = (DMAX[i] - DMIN[i]) / nsteps *
330 (double)node->values[i] + DMIN[i];
331 }
332 if (averaged) {
333 /* print averaged values */
334 sprintf(str1, "%10f", (dLow + dHigh) / 2.0);
335 G_trim_decimal(str1);
336 G_strip(str1);
337 fprintf(stdout, "%s%s", i ? fs : "", str1);
338 }
339 else {
340 /* print intervals */
341 sprintf(str1, "%10f", dLow);
342 sprintf(str2, "%10f", dHigh);
343 G_trim_decimal(str1);
344 G_trim_decimal(str2);
345 G_strip(str1);
346 G_strip(str2);
347 fprintf(stdout, "%s%s-%s", i ? fs : "", str1, str2);
348 }
349 if (with_labels) {
350 if (cat_ranges)
351 fprintf(stdout, "%s%s", fs,
352 labels[i].labels[node->values[i]]);
353 else
354 fprintf(stdout, "%sfrom %s to %s", fs,
355 Rast_get_d_cat(&dLow, &labels[i]),
356 Rast_get_d_cat(&dHigh, &labels[i]));
357 }
358 }
359
360 }
361 if (with_areas) {
362 fprintf(stdout, "%s", fs);
363 fprintf(stdout, fmt, node->area);
364 }
365 if (with_counts)
366 fprintf(stdout, "%s%ld", fs, (long)node->count);
367 if (with_percents)
368 fprintf(stdout, "%s%.2f%%", fs,
369 (double)100 * node->count / total_count);
370 fprintf(stdout, "\n");
371 }
372 }
373
374 return 0;
375 }
376