GRASS GIS 8 Programmer's Manual  8.5.0dev(2025)-b35369f6b7
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
xmedian.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/calc.h>
7 
8 /**********************************************************************
9 median(x1,x2,..,xn)
10  return median of arguments
11 **********************************************************************/
12 
13 #define SIZE_THRESHOLD 32
14 
15 static int icmp(const void *aa, const void *bb)
16 {
17  const CELL *a = aa;
18  const CELL *b = bb;
19 
20  return *a - *b;
21 }
22 
23 static int fcmp(const void *aa, const void *bb)
24 {
25  const FCELL *a = aa;
26  const FCELL *b = bb;
27 
28  if (*a < *b)
29  return -1;
30  if (*a > *b)
31  return 1;
32  return 0;
33 }
34 
35 static int dcmp(const void *aa, const void *bb)
36 {
37  const DCELL *a = aa;
38  const DCELL *b = bb;
39 
40  if (*a < *b)
41  return -1;
42  if (*a > *b)
43  return 1;
44  return 0;
45 }
46 
47 int f_median(int argc, const int *argt, void **args)
48 {
49  int size = argc * Rast_cell_size(argt[0]);
50  int i, j;
51  bool use_heap = false;
52 
53  if (argc < 1)
54  return E_ARG_LO;
55 
56  for (i = 1; i <= argc; i++)
57  if (argt[i] != argt[0])
58  return E_ARG_TYPE;
59 
60  switch (argt[0]) {
61  case CELL_TYPE: {
62  CELL stack_array[SIZE_THRESHOLD];
63  CELL *a = stack_array;
64 
65  if (argc > SIZE_THRESHOLD) {
66  a = G_malloc(size);
67  use_heap = true;
68  }
69 
70  CELL *res = args[0];
71  CELL **argv = (CELL **)&args[1];
72  CELL *a1 = &a[(argc - 1) / 2];
73  CELL *a2 = &a[argc / 2];
74 
75  for (i = 0; i < columns; i++) {
76  int nv = 0;
77 
78  for (j = 0; j < argc && !nv; j++) {
79  if (IS_NULL_C(&argv[j][i]))
80  nv = 1;
81  else
82  a[j] = argv[j][i];
83  }
84 
85  if (nv)
86  SET_NULL_C(&res[i]);
87  else {
88  qsort(a, argc, sizeof(CELL), icmp);
89  res[i] = (*a1 + *a2) / 2;
90  }
91  }
92 
93  if (use_heap) {
94  G_free(a);
95  }
96  return 0;
97  }
98  case FCELL_TYPE: {
99  FCELL stack_array[SIZE_THRESHOLD];
100  FCELL *a = stack_array;
101 
102  if (argc > SIZE_THRESHOLD) {
103  a = G_malloc(size);
104  use_heap = true;
105  }
106 
107  FCELL *res = args[0];
108  FCELL **argv = (FCELL **)&args[1];
109  FCELL *a1 = &a[(argc - 1) / 2];
110  FCELL *a2 = &a[argc / 2];
111 
112  for (i = 0; i < columns; i++) {
113  int nv = 0;
114 
115  for (j = 0; j < argc && !nv; j++) {
116  if (IS_NULL_F(&argv[j][i]))
117  nv = 1;
118  else
119  a[j] = argv[j][i];
120  }
121 
122  if (nv)
123  SET_NULL_F(&res[i]);
124  else {
125  qsort(a, argc, sizeof(FCELL), fcmp);
126  res[i] = (*a1 + *a2) / 2;
127  }
128  }
129 
130  if (use_heap) {
131  G_free(a);
132  }
133  return 0;
134  }
135  case DCELL_TYPE: {
136  DCELL stack_array[SIZE_THRESHOLD];
137  DCELL *a = stack_array;
138 
139  if (argc > SIZE_THRESHOLD) {
140  a = G_malloc(size);
141  use_heap = true;
142  }
143 
144  DCELL *res = args[0];
145  DCELL **argv = (DCELL **)&args[1];
146  DCELL *a1 = &a[(argc - 1) / 2];
147  DCELL *a2 = &a[argc / 2];
148 
149  for (i = 0; i < columns; i++) {
150  int nv = 0;
151 
152  for (j = 0; j < argc && !nv; j++) {
153  if (IS_NULL_D(&argv[j][i]))
154  nv = 1;
155  else
156  a[j] = argv[j][i];
157  }
158 
159  if (nv)
160  SET_NULL_D(&res[i]);
161  else {
162  qsort(a, argc, sizeof(DCELL), dcmp);
163  res[i] = (*a1 + *a2) / 2;
164  }
165  }
166 
167  if (use_heap) {
168  G_free(a);
169  }
170  return 0;
171  }
172  default:
173  return E_INV_TYPE;
174  }
175 }
@ 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_median(int argc, const int *argt, void **args)
Definition: xmedian.c:47
#define SIZE_THRESHOLD
Definition: xmedian.c:13