GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
xmedian.c
Go to the documentation of this file.
1#include <stdlib.h>
2
3#include <grass/gis.h>
4#include <grass/raster.h>
5#include <grass/raster.h>
6#include <grass/calc.h>
7
8/**********************************************************************
9median(x1,x2,..,xn)
10 return median of arguments
11**********************************************************************/
12
13static int icmp(const void *aa, const void *bb)
14{
15 const CELL *a = aa;
16 const CELL *b = bb;
17
18 return *a - *b;
19}
20
21static int fcmp(const void *aa, const void *bb)
22{
23 const FCELL *a = aa;
24 const FCELL *b = bb;
25
26 if (*a < *b)
27 return -1;
28 if (*a > *b)
29 return 1;
30 return 0;
31}
32
33static int dcmp(const void *aa, const void *bb)
34{
35 const DCELL *a = aa;
36 const DCELL *b = bb;
37
38 if (*a < *b)
39 return -1;
40 if (*a > *b)
41 return 1;
42 return 0;
43}
44
45int f_median(int argc, const int *argt, void **args)
46{
47 static void *array;
48 static int alloc;
49 int size = argc * Rast_cell_size(argt[0]);
50 int i, j;
51
52 if (argc < 1)
53 return E_ARG_LO;
54
55 for (i = 1; i <= argc; i++)
56 if (argt[i] != argt[0])
57 return E_ARG_TYPE;
58
59 if (size > alloc) {
60 alloc = size;
61 array = G_realloc(array, size);
62 }
63
64 switch (argt[0]) {
65 case CELL_TYPE: {
66 CELL *res = args[0];
67 CELL **argv = (CELL **)&args[1];
68 CELL *a = array;
69 CELL *a1 = &a[(argc - 1) / 2];
70 CELL *a2 = &a[argc / 2];
71
72 for (i = 0; i < columns; i++) {
73 int nv = 0;
74
75 for (j = 0; j < argc && !nv; j++) {
76 if (IS_NULL_C(&argv[j][i]))
77 nv = 1;
78 else
79 a[j] = argv[j][i];
80 }
81
82 if (nv)
83 SET_NULL_C(&res[i]);
84 else {
85 qsort(a, argc, sizeof(CELL), icmp);
86 res[i] = (*a1 + *a2) / 2;
87 }
88 }
89
90 return 0;
91 }
92 case FCELL_TYPE: {
93 FCELL *res = args[0];
94 FCELL **argv = (FCELL **)&args[1];
95 FCELL *a = array;
96 FCELL *a1 = &a[(argc - 1) / 2];
97 FCELL *a2 = &a[argc / 2];
98
99 for (i = 0; i < columns; i++) {
100 int nv = 0;
101
102 for (j = 0; j < argc && !nv; j++) {
103 if (IS_NULL_F(&argv[j][i]))
104 nv = 1;
105 else
106 a[j] = argv[j][i];
107 }
108
109 if (nv)
110 SET_NULL_F(&res[i]);
111 else {
112 qsort(a, argc, sizeof(FCELL), fcmp);
113 res[i] = (*a1 + *a2) / 2;
114 }
115 }
116
117 return 0;
118 }
119 case DCELL_TYPE: {
120 DCELL *res = args[0];
121 DCELL **argv = (DCELL **)&args[1];
122 DCELL *a = array;
123 DCELL *a1 = &a[(argc - 1) / 2];
124 DCELL *a2 = &a[argc / 2];
125
126 for (i = 0; i < columns; i++) {
127 int nv = 0;
128
129 for (j = 0; j < argc && !nv; j++) {
130 if (IS_NULL_D(&argv[j][i]))
131 nv = 1;
132 else
133 a[j] = argv[j][i];
134 }
135
136 if (nv)
137 SET_NULL_D(&res[i]);
138 else {
139 qsort(a, argc, sizeof(DCELL), dcmp);
140 res[i] = (*a1 + *a2) / 2;
141 }
142 }
143
144 return 0;
145 }
146 default:
147 return E_INV_TYPE;
148 }
149}
int columns
Definition calc.c:11
double b
int f_median(int argc, const int *argt, void **args)
Definition xmedian.c:45