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