GRASS GIS 7 Programmer's Manual  7.2.2(2017)-exported
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 *a2;
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  if (!n)
84  SET_NULL_C(&res[i]);
85  else {
86  qsort(a, n, sizeof(CELL), icmp);
87  a1 = &a[(n - 1) / 2];
88  a2 = &a[n / 2];
89  res[i] = (*a1 + *a2) / 2;
90  }
91  }
92 
93  return 0;
94  }
95  case FCELL_TYPE:
96  {
97  FCELL *res = args[0];
98  FCELL **argv = (FCELL **) &args[1];
99  FCELL *a = array;
100  FCELL *a1;
101  FCELL *a2;
102 
103  for (i = 0; i < columns; i++) {
104  int n = 0;
105 
106  for (j = 0; j < argc; j++) {
107  if (IS_NULL_F(&argv[j][i]))
108  continue;
109  a[n++] = argv[j][i];
110  }
111 
112  if (!n)
113  SET_NULL_F(&res[i]);
114  else {
115  qsort(a, n, sizeof(FCELL), fcmp);
116  a1 = &a[(n - 1) / 2];
117  a2 = &a[n / 2];
118  res[i] = (*a1 + *a2) / 2;
119  }
120  }
121 
122  return 0;
123  }
124  case DCELL_TYPE:
125  {
126  DCELL *res = args[0];
127  DCELL **argv = (DCELL **) &args[1];
128  DCELL *a = array;
129  DCELL *a1;
130  DCELL *a2;
131 
132  for (i = 0; i < columns; i++) {
133  int n = 0;
134 
135  for (j = 0; j < argc; j++) {
136  if (IS_NULL_D(&argv[j][i]))
137  continue;
138  a[n++] = argv[j][i];
139  }
140 
141  if (!n)
142  SET_NULL_D(&res[i]);
143  else {
144  qsort(a, n, sizeof(DCELL), dcmp);
145  a1 = &a[(n - 1) / 2];
146  a2 = &a[n / 2];
147  res[i] = (*a1 + *a2) / 2;
148  }
149  }
150 
151  return 0;
152  }
153  default:
154  return E_INV_TYPE;
155  }
156 }
int columns
Definition: calc.c:12
double b
int f_nmedian(int argc, const int *argt, void **args)
Definition: xnmedian.c:46