GRASS GIS 8 Programmer's Manual  8.5.0dev(2025)-f25cea218f
xnmedian.c
Go to the documentation of this file.
1 #include <stdlib.h>
2 #include <stdbool.h>
3 
4 #include <grass/gis.h>
5 #include <grass/raster.h>
6 #include <grass/raster.h>
7 #include <grass/calc.h>
8 
9 /**********************************************************************
10 median(x1,x2,..,xn)
11  return median of arguments
12 **********************************************************************/
13 
14 #define SIZE_THRESHOLD 32
15 
16 static int icmp(const void *aa, const void *bb)
17 {
18  const CELL *a = aa;
19  const CELL *b = bb;
20 
21  return *a - *b;
22 }
23 
24 static int fcmp(const void *aa, const void *bb)
25 {
26  const FCELL *a = aa;
27  const FCELL *b = bb;
28 
29  if (*a < *b)
30  return -1;
31  if (*a > *b)
32  return 1;
33  return 0;
34 }
35 
36 static int dcmp(const void *aa, const void *bb)
37 {
38  const DCELL *a = aa;
39  const DCELL *b = bb;
40 
41  if (*a < *b)
42  return -1;
43  if (*a > *b)
44  return 1;
45  return 0;
46 }
47 
48 int f_nmedian(int argc, const int *argt, void **args)
49 {
50  int size = argc * Rast_cell_size(argt[0]);
51  int i, j;
52  bool use_heap = false;
53 
54  if (argc < 1)
55  return E_ARG_LO;
56 
57  for (i = 1; i <= argc; i++)
58  if (argt[i] != argt[0])
59  return E_ARG_TYPE;
60 
61  switch (argt[0]) {
62  case CELL_TYPE: {
63  CELL stack_array[SIZE_THRESHOLD];
64  CELL *a = stack_array;
65 
66  if (argc > SIZE_THRESHOLD) {
67  a = G_malloc(size);
68  use_heap = true;
69  }
70 
71  CELL *res = args[0];
72  CELL **argv = (CELL **)&args[1];
73  CELL a1;
74  CELL *resc;
75 
76  for (i = 0; i < columns; i++) {
77  int n = 0;
78 
79  for (j = 0; j < argc; j++) {
80  if (IS_NULL_C(&argv[j][i]))
81  continue;
82  a[n++] = argv[j][i];
83  }
84 
85  resc = &res[i];
86 
87  if (!n)
88  SET_NULL_C(resc);
89  else {
90  qsort(a, n, sizeof(CELL), icmp);
91  *resc = a[n / 2];
92  if ((n & 1) == 0) {
93  a1 = a[(n - 1) / 2];
94  if (*resc != a1)
95  *resc = (*resc + a1) / 2;
96  }
97  }
98  }
99 
100  if (use_heap) {
101  G_free(a);
102  }
103  return 0;
104  }
105 
106  case FCELL_TYPE: {
107  FCELL stack_array[SIZE_THRESHOLD];
108  FCELL *a = stack_array;
109 
110  if (argc > SIZE_THRESHOLD) {
111  a = G_malloc(size);
112  use_heap = true;
113  }
114 
115  FCELL *res = args[0];
116  FCELL **argv = (FCELL **)&args[1];
117  FCELL a1;
118  FCELL *resc;
119 
120  for (i = 0; i < columns; i++) {
121  int n = 0;
122 
123  for (j = 0; j < argc; j++) {
124  if (IS_NULL_F(&argv[j][i]))
125  continue;
126  a[n++] = argv[j][i];
127  }
128 
129  resc = &res[i];
130 
131  if (!n)
132  SET_NULL_F(resc);
133  else {
134  qsort(a, n, sizeof(FCELL), fcmp);
135  *resc = a[n / 2];
136  if ((n & 1) == 0) {
137  a1 = a[(n - 1) / 2];
138  if (*resc != a1)
139  *resc = (*resc + a1) / 2;
140  }
141  }
142  }
143 
144  if (use_heap) {
145  G_free(a);
146  }
147  return 0;
148  }
149 
150  case DCELL_TYPE: {
151  DCELL stack_array[SIZE_THRESHOLD];
152  DCELL *a = stack_array;
153 
154  if (argc > SIZE_THRESHOLD) {
155  a = G_malloc(size);
156  use_heap = true;
157  }
158 
159  DCELL *res = args[0];
160  DCELL **argv = (DCELL **)&args[1];
161  DCELL a1;
162  DCELL *resc;
163 
164  for (i = 0; i < columns; i++) {
165  int n = 0;
166 
167  for (j = 0; j < argc; j++) {
168  if (IS_NULL_D(&argv[j][i]))
169  continue;
170  a[n++] = argv[j][i];
171  }
172 
173  resc = &res[i];
174 
175  if (!n)
176  SET_NULL_D(resc);
177  else {
178  qsort(a, n, sizeof(DCELL), dcmp);
179  *resc = a[n / 2];
180  if ((n & 1) == 0) {
181  a1 = a[(n - 1) / 2];
182  if (*resc != a1)
183  *resc = (*resc + a1) / 2;
184  }
185  }
186  }
187 
188  if (use_heap) {
189  G_free(a);
190  }
191  return 0;
192  }
193 
194  default:
195  return E_INV_TYPE;
196  }
197 }
@ E_INV_TYPE
Definition: calc.h:15
@ E_ARG_TYPE
Definition: calc.h:13
@ E_ARG_LO
Definition: calc.h:11
#define IS_NULL_C(x)
Definition: calc.h:26
#define SET_NULL_D(x)
Definition: calc.h:32
int columns
Definition: calc.c:11
#define SET_NULL_C(x)
Definition: calc.h:30
#define IS_NULL_F(x)
Definition: calc.h:27
#define IS_NULL_D(x)
Definition: calc.h:28
#define SET_NULL_F(x)
Definition: calc.h:31
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
#define G_malloc(n)
Definition: defs/gis.h:94
size_t Rast_cell_size(RASTER_MAP_TYPE)
Returns size of a raster cell in bytes.
Definition: alloc_cell.c:38
float FCELL
Definition: gis.h:636
double DCELL
Definition: gis.h:635
int CELL
Definition: gis.h:634
double b
Definition: r_raster.c:39
#define FCELL_TYPE
Definition: raster.h:12
#define DCELL_TYPE
Definition: raster.h:13
#define CELL_TYPE
Definition: raster.h:11
int f_nmedian(int argc, const int *argt, void **args)
Definition: xnmedian.c:48
#define SIZE_THRESHOLD
Definition: xnmedian.c:14