GRASS GIS 8 Programmer's Manual 8.2.0(2022)-exported
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/**********************************************************************
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_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}
int columns
Definition: calc.c:12
double b
int f_median(int argc, const int *argt, void **args)
Definition: xmedian.c:46