GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
xnmedian.c
Go to the documentation of this file.
1 
2 #include <stdlib.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 static int icmp(const void *aa, const void *bb)
15 {
16  const CELL *a = aa;
17  const CELL *b = bb;
18 
19  return *a - *b;
20 }
21 
22 static int fcmp(const void *aa, const void *bb)
23 {
24  const FCELL *a = aa;
25  const FCELL *b = bb;
26 
27  if (*a < *b)
28  return -1;
29  if (*a > *b)
30  return 1;
31  return 0;
32 }
33 
34 static int dcmp(const void *aa, const void *bb)
35 {
36  const DCELL *a = aa;
37  const DCELL *b = bb;
38 
39  if (*a < *b)
40  return -1;
41  if (*a > *b)
42  return 1;
43  return 0;
44 }
45 
46 int f_nmedian(int argc, const int *argt, void **args)
47 {
48  static void *array;
49  static int alloc;
50  int size = argc * Rast_cell_size(argt[0]);
51  int i, j;
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  if (size > alloc) {
61  alloc = size;
62  array = G_realloc(array, size);
63  }
64 
65  switch (argt[0]) {
66  case CELL_TYPE:
67  {
68  CELL *res = args[0];
69  CELL **argv = (CELL **) &args[1];
70  CELL *a = array;
71  CELL a1;
72  CELL *resc;
73 
74  for (i = 0; i < columns; i++) {
75  int n = 0;
76 
77  for (j = 0; j < argc; j++) {
78  if (IS_NULL_C(&argv[j][i]))
79  continue;
80  a[n++] = argv[j][i];
81  }
82 
83  resc = &res[i];
84 
85  if (!n)
86  SET_NULL_C(resc);
87  else {
88  qsort(a, n, sizeof(CELL), icmp);
89  *resc = a[n / 2];
90  if ((n & 1) == 0) {
91  a1 = a[(n - 1) / 2];
92  if (*resc != a1)
93  *resc = (*resc + a1) / 2;
94  }
95  }
96  }
97 
98  return 0;
99  }
100  case FCELL_TYPE:
101  {
102  FCELL *res = args[0];
103  FCELL **argv = (FCELL **) &args[1];
104  FCELL *a = array;
105  FCELL a1;
106  FCELL *resc;
107 
108  for (i = 0; i < columns; i++) {
109  int n = 0;
110 
111  for (j = 0; j < argc; j++) {
112  if (IS_NULL_F(&argv[j][i]))
113  continue;
114  a[n++] = argv[j][i];
115  }
116 
117  resc = &res[i];
118 
119  if (!n)
120  SET_NULL_F(resc);
121  else {
122  qsort(a, n, sizeof(FCELL), fcmp);
123  *resc = a[n / 2];
124  if ((n & 1) == 0) {
125  a1 = a[(n - 1) / 2];
126  if (*resc != a1)
127  *resc = (*resc + a1) / 2;
128  }
129  }
130  }
131 
132  return 0;
133  }
134  case DCELL_TYPE:
135  {
136  DCELL *res = args[0];
137  DCELL **argv = (DCELL **) &args[1];
138  DCELL *a = array;
139  DCELL a1;
140  DCELL *resc;
141 
142  for (i = 0; i < columns; i++) {
143  int n = 0;
144 
145  for (j = 0; j < argc; j++) {
146  if (IS_NULL_D(&argv[j][i]))
147  continue;
148  a[n++] = argv[j][i];
149  }
150 
151  resc = &res[i];
152 
153  if (!n)
154  SET_NULL_D(resc);
155  else {
156  qsort(a, n, sizeof(DCELL), dcmp);
157  *resc = a[n / 2];
158  if ((n & 1) == 0) {
159  a1 = a[(n - 1) / 2];
160  if (*resc != a1)
161  *resc = (*resc + a1) / 2;
162  }
163  }
164  }
165 
166  return 0;
167  }
168  default:
169  return E_INV_TYPE;
170  }
171 }
#define CELL_TYPE
Definition: raster.h:11
#define SET_NULL_C(x)
Definition: calc.h:32
double DCELL
Definition: gis.h:614
#define IS_NULL_F(x)
Definition: calc.h:29
#define IS_NULL_C(x)
Definition: calc.h:28
int columns
Definition: calc.c:12
#define DCELL_TYPE
Definition: raster.h:13
double b
Definition: r_raster.c:39
size_t Rast_cell_size(RASTER_MAP_TYPE)
Returns size of a raster cell in bytes.
Definition: alloc_cell.c:39
int f_nmedian(int argc, const int *argt, void **args)
Definition: xnmedian.c:46
Definition: calc.h:12
float FCELL
Definition: gis.h:615
#define IS_NULL_D(x)
Definition: calc.h:30
int CELL
Definition: gis.h:613
#define G_realloc(p, n)
Definition: defs/gis.h:114
#define FCELL_TYPE
Definition: raster.h:12
#define SET_NULL_F(x)
Definition: calc.h:33
#define SET_NULL_D(x)
Definition: calc.h:34