GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
xmedian.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_median(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 = &a[(argc - 1) / 2];
72  CELL *a2 = &a[argc / 2];
73 
74  for (i = 0; i < columns; i++) {
75  int nv = 0;
76 
77  for (j = 0; j < argc && !nv; j++) {
78  if (IS_NULL_C(&argv[j][i]))
79  nv = 1;
80  else
81  a[j] = argv[j][i];
82  }
83 
84  if (nv)
85  SET_NULL_C(&res[i]);
86  else {
87  qsort(a, argc, sizeof(CELL), icmp);
88  res[i] = (*a1 + *a2) / 2;
89  }
90  }
91 
92  return 0;
93  }
94  case FCELL_TYPE:
95  {
96  FCELL *res = args[0];
97  FCELL **argv = (FCELL **) &args[1];
98  FCELL *a = array;
99  FCELL *a1 = &a[(argc - 1) / 2];
100  FCELL *a2 = &a[argc / 2];
101 
102  for (i = 0; i < columns; i++) {
103  int nv = 0;
104 
105  for (j = 0; j < argc && !nv; j++) {
106  if (IS_NULL_F(&argv[j][i]))
107  nv = 1;
108  else
109  a[j] = argv[j][i];
110  }
111 
112  if (nv)
113  SET_NULL_F(&res[i]);
114  else {
115  qsort(a, argc, sizeof(FCELL), fcmp);
116  res[i] = (*a1 + *a2) / 2;
117  }
118  }
119 
120  return 0;
121  }
122  case DCELL_TYPE:
123  {
124  DCELL *res = args[0];
125  DCELL **argv = (DCELL **) &args[1];
126  DCELL *a = array;
127  DCELL *a1 = &a[(argc - 1) / 2];
128  DCELL *a2 = &a[argc / 2];
129 
130  for (i = 0; i < columns; i++) {
131  int nv = 0;
132 
133  for (j = 0; j < argc && !nv; j++) {
134  if (IS_NULL_D(&argv[j][i]))
135  nv = 1;
136  else
137  a[j] = argv[j][i];
138  }
139 
140  if (nv)
141  SET_NULL_D(&res[i]);
142  else {
143  qsort(a, argc, sizeof(DCELL), dcmp);
144  res[i] = (*a1 + *a2) / 2;
145  }
146  }
147 
148  return 0;
149  }
150  default:
151  return E_INV_TYPE;
152  }
153 }
#define CELL_TYPE
Definition: raster.h:11
#define SET_NULL_C(x)
Definition: calc.h:32
int f_median(int argc, const int *argt, void **args)
Definition: xmedian.c:46
double DCELL
Definition: gis.h:603
#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
Definition: calc.h:12
float FCELL
Definition: gis.h:604
#define IS_NULL_D(x)
Definition: calc.h:30
int CELL
Definition: gis.h:602
#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