GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-f8115df121
cell_stats.c
Go to the documentation of this file.
1 /*!
2  * \file lib/raster/cell_stats.c
3  *
4  * \brief Raster Library - Raster cell statistics
5  *
6  * (C) 2001-2009 GRASS Development Team
7  *
8  * This program is free software under the GNU General Public License
9  * (>=v2). Read the file COPYING that comes with GRASS for details.
10  *
11  * \author Original author CERL
12  */
13 
14 #include <stdlib.h>
15 
16 #include <grass/gis.h>
17 #include <grass/raster.h>
18 
19 #define INCR 10
20 #define SHIFT 6
21 
22 static const int NCATS = 1 << SHIFT;
23 
24 #define NODE struct Cell_stats_node
25 
26 static int next_node(struct Cell_stats *);
27 static void init_node(NODE *, int, int);
28 
29 /*!
30  * \brief Initialize cell stats
31  *
32  * This routine, which must be called first, initializes the
33  * Cell_stats structure.
34  *
35  * Set the count for NULL-values to zero.
36  *
37  * \param s pointer to Cell_stats structure
38  */
40 {
41  s->N = 0;
42  s->tlen = INCR;
43  s->node = (NODE *)G_malloc(s->tlen * sizeof(NODE));
44  s->null_data_count = 0;
45 }
46 
47 /*!
48  * \brief Add data to cell stats
49  *
50  * The <i>n</i> CELL values in the <i>data</i> array are inserted (and
51  * counted) in the Cell_stats structure.
52  *
53  * Look for NULLs and update the NULL-value count.
54  *
55  * \param cell raster values
56  * \param n number of values
57  * \param s pointer to Cell_stats structure which holds cell stats info
58  *
59  * \return 1 on failure
60  * \return 0 on success
61  */
62 int Rast_update_cell_stats(const CELL *cell, int n, struct Cell_stats *s)
63 {
64  CELL cat;
65  int p, q;
66  int idx, offset;
67  int N;
68  NODE *node, *pnode;
69  NODE *new_node;
70 
71  if (n <= 0)
72  return 1;
73 
74  node = s->node;
75 
76  /* first non-null node is special case */
77  if ((N = s->N) == 0) {
78  cat = *cell++;
79  while (Rast_is_c_null_value(&cat)) {
80  s->null_data_count++;
81  cat = *cell++;
82  n--;
83  }
84  if (n > 0) { /* if there are some non-null cells */
85  N = 1;
86  if (cat < 0) {
87  idx = -((-cat) >> SHIFT) - 1;
88  offset = cat + ((-idx) << SHIFT) - 1;
89  }
90  else {
91  idx = cat >> SHIFT;
92  offset = cat - (idx << SHIFT);
93  }
94  fflush(stderr);
95  init_node(&node[1], idx, offset);
96  node[1].right = 0;
97  n--;
98  }
99  }
100  while (n-- > 0) {
101  cat = *cell++;
102  if (Rast_is_c_null_value(&cat)) {
103  s->null_data_count++;
104  continue;
105  }
106  if (cat < 0) {
107  idx = -((-cat) >> SHIFT) - 1;
108  offset = cat + ((-idx) << SHIFT) - 1;
109  }
110  else {
111  idx = cat >> SHIFT;
112  offset = cat - (idx << SHIFT);
113  }
114 
115  q = 1;
116  while (q > 0) {
117  pnode = &node[p = q];
118  if (pnode->idx == idx) {
119  pnode->count[offset]++;
120  break;
121  }
122  if (pnode->idx > idx)
123  q = pnode->left; /* go left */
124  else
125  q = pnode->right; /* go right */
126  }
127  if (q > 0)
128  continue; /* found */
129 
130  /* new node */
131  N++;
132 
133  /* grow the tree? */
134  if (N >= s->tlen) {
135  node = (NODE *)G_realloc((char *)node,
136  sizeof(NODE) * (s->tlen += INCR));
137  pnode = &node[p]; /* realloc moves node, must reassign pnode */
138  }
139 
140  /* add node to tree */
141  init_node(new_node = &node[N], idx, offset);
142 
143  if (pnode->idx > idx) {
144  new_node->right = -p; /* create thread */
145  pnode->left = N; /* insert left */
146  }
147  else {
148  new_node->right = pnode->right; /* copy right link/thread */
149  pnode->right = N; /* add right */
150  }
151  } /* while n-- > 0 */
152  s->N = N;
153  s->node = node;
154 
155  return 0;
156 }
157 
158 static void init_node(NODE *node, int idx, int offset)
159 {
160  long *count;
161  int i;
162 
163  count = node->count = (long *)G_calloc(i = NCATS, sizeof(long));
164  while (i--)
165  *count++ = 0;
166  node->idx = idx;
167  node->count[offset] = 1;
168  node->left = 0;
169 }
170 
171 /*!
172  * \brief Random query of cell stats
173  *
174  * This routine allows a random query of the Cell_stats structure. The
175  * \p count associated with the raster value \p cat is
176  * set. The routine returns 1 if \p cat was found in the
177  * structure, 0 otherwise.
178  *
179  * Allow finding the count for the NULL-value.
180  *
181  * \param cat raster value
182  * \param[out] count count
183  * \param s pointer to Cell_stats structure which holds cell stats info
184  *
185  * \return 1 if found
186  * \return 0 if not found
187  */
188 int Rast_find_cell_stat(CELL cat, long *count, const struct Cell_stats *s)
189 {
190  int q;
191  int idx;
192  int offset;
193 
194  *count = 0;
195  if (Rast_is_c_null_value(&cat)) {
196  *count = s->null_data_count;
197  return (*count != 0);
198  }
199 
200  if (s->N <= 0)
201  return 0;
202 
203  /*
204  if (cat < 0)
205  {
206  idx = -(-cat/NCATS) - 1;
207  offset = cat - idx*NCATS - 1;
208  }
209  else
210  {
211  idx = cat/NCATS;
212  offset = cat - idx*NCATS;
213  }
214  */
215  if (cat < 0) {
216  idx = -((-cat) >> SHIFT) - 1;
217  offset = cat + ((-idx) << SHIFT) - 1;
218  }
219  else {
220  idx = cat >> SHIFT;
221  offset = cat - (idx << SHIFT);
222  }
223 
224  q = 1;
225  while (q > 0) {
226  if (s->node[q].idx == idx) {
227  *count = s->node[q].count[offset];
228  return (*count != 0);
229  }
230  if (s->node[q].idx > idx)
231  q = s->node[q].left; /* go left */
232  else
233  q = s->node[q].right; /* go right */
234  }
235  return 0;
236 }
237 
238 /*!
239  * \brief Reset/rewind cell stats
240  *
241  * The structure <i>s</i> is rewound (i.e., positioned at the first
242  * raster category) so that sorted sequential retrieval can begin.
243  *
244  * \param s pointer to Cell_stats structure which holds cell stats info
245  *
246  * \return 0
247  */
249 {
250  int q;
251 
252  if (s->N <= 0)
253  return 1;
254  /* start at root and go all the way to the left */
255  s->curp = 1;
256  while ((q = s->node[s->curp].left))
257  s->curp = q;
258  s->curoffset = -1;
259 
260  return 0;
261 }
262 
263 static int next_node(struct Cell_stats *s)
264 {
265  int q;
266 
267  /* go to the right */
268  s->curp = s->node[s->curp].right;
269 
270  if (s->curp == 0) /* no more */
271  return 0;
272 
273  if (s->curp < 0) { /* thread. stop here */
274  s->curp = -(s->curp);
275  return 1;
276  }
277 
278  while ((q = s->node[s->curp].left)) /* now go all the way left */
279  s->curp = q;
280 
281  return 1;
282 }
283 
284 /*!
285  * \brief Retrieve sorted cell stats
286  *
287  * Retrieves the next <i>cat, count</i> combination from the
288  * structure. Returns 0 if there are no more items, non-zero if there
289  * are more. For example:
290  *
291  \code
292  struct Cell_stats s;
293  CELL cat;
294  long count;
295 
296  // updating <b>s</b> occurs here
297 
298  Rast_rewind_cell_stats(&s);
299  while (Rast_next_cell_stat(&cat,&count,&s)
300  fprintf(stdout, "%ld %ld\n", (long) cat, count);
301  \endcode
302  *
303  * Do not return a record for the NULL-value
304  *
305  * \param cat raster value
306  * \param[out] count
307  * \param s pointer to Cell_stats structure which holds cell stats info
308  *
309  * \return 0 if there are no more items
310  * \return non-zero if there are more
311  */
312 int Rast_next_cell_stat(CELL *cat, long *count, struct Cell_stats *s)
313 {
314  int idx;
315 
316  /* first time report stats for null */
317  /* decided not to return stats for null in this function
318  static int null_reported = 0;
319  if(!null_reported && s->null_data_count > 0)
320  {
321  *count = s->null_data_count;
322  Rast_set_c_null_value(&cat,1);
323  null_reported = 1;
324  return 1;
325  }
326  */
327  if (s->N <= 0)
328  return 0;
329  for (;;) {
330  s->curoffset++;
331  if (s->curoffset >= NCATS) {
332  if (!next_node(s))
333  return 0;
334  s->curoffset = -1;
335  continue;
336  }
337  if ((*count = s->node[s->curp].count[s->curoffset])) {
338  idx = s->node[s->curp].idx;
339 
340  /*
341  if (idx < 0)
342  *cat = idx*NCATS + s->curoffset+1;
343  else
344  *cat = idx*NCATS + s->curoffset;
345  */
346  if (idx < 0)
347  *cat = -((-idx) << SHIFT) + s->curoffset + 1;
348  else
349  *cat = (idx << SHIFT) + s->curoffset;
350 
351  return 1;
352  }
353  }
354 }
355 
356 /*!
357  * \brief Get number of null values.
358  *
359  * Get a number of null values from stats structure.
360  *
361  * Note: when reporting values which appear in a map using
362  * Rast_next_cell_stats(), to get stats for null, call
363  * Rast_get_stats_for_null_value() first, since Rast_next_cell_stats() does
364  * not report stats for null.
365  *
366  * \param count count
367  * \param s pointer to Cell_stats structure which holds cell stats info
368  */
369 void Rast_get_stats_for_null_value(long *count, const struct Cell_stats *s)
370 {
371  *count = s->null_data_count;
372 }
373 
374 /*!
375  * \brief Free cell stats structure
376  *
377  * The memory associated with structure <i>s</i> is freed. This
378  * routine may be called any time after calling Rast_init_cell_stats().
379  *
380  * \param s pointer to Cell_stats structure
381  */
383 {
384  int i;
385 
386  for (i = 1; i <= s->N; i++)
387  G_free(s->node[i].count);
388  G_free(s->node);
389 }
int Rast_rewind_cell_stats(struct Cell_stats *s)
Reset/rewind cell stats.
Definition: cell_stats.c:248
int Rast_find_cell_stat(CELL cat, long *count, const struct Cell_stats *s)
Random query of cell stats.
Definition: cell_stats.c:188
int Rast_next_cell_stat(CELL *cat, long *count, struct Cell_stats *s)
Retrieve sorted cell stats.
Definition: cell_stats.c:312
int Rast_update_cell_stats(const CELL *cell, int n, struct Cell_stats *s)
Add data to cell stats.
Definition: cell_stats.c:62
#define NODE
Definition: cell_stats.c:24
void Rast_init_cell_stats(struct Cell_stats *s)
Initialize cell stats.
Definition: cell_stats.c:39
void Rast_get_stats_for_null_value(long *count, const struct Cell_stats *s)
Get number of null values.
Definition: cell_stats.c:369
void Rast_free_cell_stats(struct Cell_stats *s)
Free cell stats structure.
Definition: cell_stats.c:382
#define INCR
Definition: cell_stats.c:19
#define SHIFT
Definition: cell_stats.c:20
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
#define G_realloc(p, n)
Definition: defs/gis.h:96
#define G_calloc(m, n)
Definition: defs/gis.h:95
#define G_malloc(n)
Definition: defs/gis.h:94
#define Rast_is_c_null_value(cellVal)
Definition: defs/raster.h:401
#define N
Definition: e_intersect.c:926
int CELL
Definition: gis.h:625
int count
struct Cell_stats::Cell_stats_node * node
long null_data_count
Definition: raster.h:192
int N
Definition: raster.h:190
int curoffset
Definition: raster.h:193
int tlen
Definition: raster.h:189
int curp
Definition: raster.h:191