GRASS GIS 8 Programmer's Manual 8.2.0(2022)-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/**********************************************************************
10median(x1,x2,..,xn)
11 return median of arguments
12**********************************************************************/
13
14static 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
22static 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
34static 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
46int 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}
int columns
Definition: calc.c:12
double b
int f_nmedian(int argc, const int *argt, void **args)
Definition: xnmedian.c:46