GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-7413740dd8
xmedian.c
Go to the documentation of this file.
1 #include <stdlib.h>
2 
3 #include <grass/gis.h>
4 #include <grass/raster.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 static int icmp(const void *aa, const void *bb)
14 {
15  const CELL *a = aa;
16  const CELL *b = bb;
17 
18  return *a - *b;
19 }
20 
21 static int fcmp(const void *aa, const void *bb)
22 {
23  const FCELL *a = aa;
24  const FCELL *b = bb;
25 
26  if (*a < *b)
27  return -1;
28  if (*a > *b)
29  return 1;
30  return 0;
31 }
32 
33 static int dcmp(const void *aa, const void *bb)
34 {
35  const DCELL *a = aa;
36  const DCELL *b = bb;
37 
38  if (*a < *b)
39  return -1;
40  if (*a > *b)
41  return 1;
42  return 0;
43 }
44 
45 int f_median(int argc, const int *argt, void **args)
46 {
47  static void *array;
48  static int alloc;
49  int size = argc * Rast_cell_size(argt[0]);
50  int i, j;
51 
52  if (argc < 1)
53  return E_ARG_LO;
54 
55  for (i = 1; i <= argc; i++)
56  if (argt[i] != argt[0])
57  return E_ARG_TYPE;
58 
59  if (size > alloc) {
60  alloc = size;
61  array = G_realloc(array, size);
62  }
63 
64  switch (argt[0]) {
65  case CELL_TYPE: {
66  CELL *res = args[0];
67  CELL **argv = (CELL **)&args[1];
68  CELL *a = array;
69  CELL *a1 = &a[(argc - 1) / 2];
70  CELL *a2 = &a[argc / 2];
71 
72  for (i = 0; i < columns; i++) {
73  int nv = 0;
74 
75  for (j = 0; j < argc && !nv; j++) {
76  if (IS_NULL_C(&argv[j][i]))
77  nv = 1;
78  else
79  a[j] = argv[j][i];
80  }
81 
82  if (nv)
83  SET_NULL_C(&res[i]);
84  else {
85  qsort(a, argc, sizeof(CELL), icmp);
86  res[i] = (*a1 + *a2) / 2;
87  }
88  }
89 
90  return 0;
91  }
92  case FCELL_TYPE: {
93  FCELL *res = args[0];
94  FCELL **argv = (FCELL **)&args[1];
95  FCELL *a = array;
96  FCELL *a1 = &a[(argc - 1) / 2];
97  FCELL *a2 = &a[argc / 2];
98 
99  for (i = 0; i < columns; i++) {
100  int nv = 0;
101 
102  for (j = 0; j < argc && !nv; j++) {
103  if (IS_NULL_F(&argv[j][i]))
104  nv = 1;
105  else
106  a[j] = argv[j][i];
107  }
108 
109  if (nv)
110  SET_NULL_F(&res[i]);
111  else {
112  qsort(a, argc, sizeof(FCELL), fcmp);
113  res[i] = (*a1 + *a2) / 2;
114  }
115  }
116 
117  return 0;
118  }
119  case DCELL_TYPE: {
120  DCELL *res = args[0];
121  DCELL **argv = (DCELL **)&args[1];
122  DCELL *a = array;
123  DCELL *a1 = &a[(argc - 1) / 2];
124  DCELL *a2 = &a[argc / 2];
125 
126  for (i = 0; i < columns; i++) {
127  int nv = 0;
128 
129  for (j = 0; j < argc && !nv; j++) {
130  if (IS_NULL_D(&argv[j][i]))
131  nv = 1;
132  else
133  a[j] = argv[j][i];
134  }
135 
136  if (nv)
137  SET_NULL_D(&res[i]);
138  else {
139  qsort(a, argc, sizeof(DCELL), dcmp);
140  res[i] = (*a1 + *a2) / 2;
141  }
142  }
143 
144  return 0;
145  }
146  default:
147  return E_INV_TYPE;
148  }
149 }
@ 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
#define G_realloc(p, n)
Definition: defs/gis.h:96
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:627
double DCELL
Definition: gis.h:626
int CELL
Definition: gis.h:625
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:45