1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <math.h>
5 #include <unistd.h>
6 #include <grass/gis.h>
7 #include <grass/raster.h>
8 #include <grass/glocale.h>
9
10 #include "local_proto.h"
11
12
13 /*--------------------------------------------------------
14 HISTOGRAM ANALYSIS
15 Define un factor de escala = hist_n/100 con objeto
16 de dividir el entero 1 por 100/hist_n partes y
17 aumentar la precision.
18
19 Afecta al almacenamiento en el histograma pero
20 modifica el calculo de quantiles y momentos.
21 --------------------------------------------------------*/
22
23 /* Global variable
24 allow use as parameter in the command line */
25 int hist_n = 100; /* interval of real data 100/hist_n */
26
hist_put(double t,int hist[])27 void hist_put(double t, int hist[])
28 {
29 int i;
30
31 /* scale factor */
32 i = (int)(t * ((double)hist_n / 100.));
33
34 if (i < 1)
35 i = 1;
36 if (i > hist_n)
37 i = hist_n;
38
39 hist[i - 1] += 1;
40 }
41
42 /* histogram moment */
moment(int n,int hist[],int k)43 double moment(int n, int hist[], int k)
44 {
45 int i, total;
46 double value, mean;
47
48 k = 0;
49
50 total = 0;
51 mean = 0.;
52 for (i = 0; i < hist_n; i++) {
53 total += hist[i];
54 mean += (double)(i * hist[i]);
55 }
56 mean /= ((double)total); /* histogram mean */
57
58 value = 0.;
59 for (i = 0; i < hist_n; i++) {
60 value += (pow((i - mean), n) * ((double)hist[i]));
61 }
62 value /= (double)(total - k);
63
64 return (value / pow((double)hist_n / 100., n) );
65 }
66
67 /* Real data quantile */
quantile(double q,int hist[])68 double quantile(double q, int hist[])
69 {
70 int i, total;
71 double value, qmax, qmin;
72
73 total = 0;
74 for (i = 0; i < hist_n; i++) {
75 total += hist[i];
76 }
77
78 value = 0;
79 qmax = 1.;
80 for (i = hist_n - 1; i >= 0; i--) {
81 qmin = qmax - (double)hist[i] / (double)total;
82 if (q >= qmin) {
83 value = (q - qmin) / (qmax - qmin) + (i - 1);
84 break;
85 }
86 qmax = qmin;
87 }
88
89 /* remove scale factor */
90 return (value / ((double)hist_n / 100.));
91 }
92
93 /*--------------------------------------------------------
94 FILTER HOLES OF CLOUDS
95 This a >=50% filter of 3x3
96 if >= 50% vecinos cloud then pixel set to cloud
97 --------------------------------------------------------*/
98
pval(void * rast,int i)99 int pval(void *rast, int i)
100 {
101 void *ptr = (void *)((CELL *) rast + i);
102
103 if (Rast_is_c_null_value(ptr))
104 return 0;
105 else
106 return (int)((CELL *) rast)[i];
107 }
108
filter_holes(Gfile * out)109 void filter_holes(Gfile * out)
110 {
111 int row, col, nrows, ncols;
112
113 void *arast, *brast, *crast;
114 int i, pixel[9], cold, warm, shadow, nulo, lim;
115
116 Gfile tmp;
117
118 nrows = Rast_window_rows();
119 ncols = Rast_window_cols();
120
121 if (nrows < 3 || ncols < 3)
122 return;
123
124 /* Open to read */
125 if ((out->fd = Rast_open_old(out->name, "")) < 0)
126 G_fatal_error(_("Unable to open raster map <%s>"), out->name);
127
128 arast = Rast_allocate_buf(CELL_TYPE);
129 brast = Rast_allocate_buf(CELL_TYPE);
130 crast = Rast_allocate_buf(CELL_TYPE);
131
132 /* Open to write */
133 sprintf(tmp.name, "_%d.BBB", getpid());
134 tmp.rast = Rast_allocate_buf(CELL_TYPE);
135 if ((tmp.fd = Rast_open_new(tmp.name, CELL_TYPE)) < 0)
136 G_fatal_error(_("Unable to create raster map <%s>"), tmp.name);
137
138 G_important_message(_("Filling small holes in clouds..."));
139
140 /* Se puede acelerar creandolos nulos y luego arast = brast
141 brast = crast y cargando crast solamente
142 G_set_f_null_value(cell[2], ncols);
143 */
144
145 for (row = 0; row < nrows; row++) {
146 G_percent(row, nrows, 2);
147 /* Read row values */
148 if (row != 0) {
149 Rast_get_c_row(out->fd, arast, row - 1);
150 }
151 Rast_get_c_row(out->fd, brast, row);
152 if (row != (nrows - 1)) {
153 Rast_get_c_row(out->fd, crast, row + 1);
154 }
155 /* Analysis of all pixels */
156 for (col = 0; col < ncols; col++) {
157 pixel[0] = pval(brast, col);
158 if (pixel[0] == 0) {
159 if (row == 0) {
160 pixel[1] = -1;
161 pixel[2] = -1;
162 pixel[3] = -1;
163 if (col == 0) {
164 pixel[4] = -1;
165 pixel[5] = pval(brast, col + 1);
166 pixel[6] = -1;
167 pixel[7] = pval(crast, col);
168 pixel[8] = pval(crast, col + 1);
169 }
170 else if (col != (ncols - 1)) {
171 pixel[4] = pval(brast, col - 1);
172 pixel[5] = pval(brast, col + 1);
173 pixel[6] = pval(crast, col - 1);
174 pixel[7] = pval(crast, col);
175 pixel[8] = pval(crast, col + 1);
176 }
177 else {
178 pixel[4] = pval(brast, col - 1);
179 pixel[5] = -1;
180 pixel[6] = pval(crast, col - 1);
181 pixel[7] = pval(crast, col);
182 pixel[8] = -1;
183 }
184 }
185 else if (row != (nrows - 1)) {
186 if (col == 0) {
187 pixel[1] = -1;
188 pixel[2] = pval(arast, col);
189 pixel[3] = pval(arast, col + 1);
190 pixel[4] = -1;
191 pixel[5] = pval(brast, col + 1);
192 pixel[6] = -1;
193 pixel[7] = pval(crast, col);
194 pixel[8] = pval(crast, col + 1);
195 }
196 else if (col != (ncols - 1)) {
197 pixel[1] = pval(arast, col - 1);
198 pixel[2] = pval(arast, col);
199 pixel[3] = pval(arast, col + 1);
200 pixel[4] = pval(brast, col - 1);
201 pixel[5] = pval(brast, col + 1);
202 pixel[6] = pval(crast, col - 1);
203 pixel[7] = pval(crast, col);
204 pixel[8] = pval(crast, col + 1);
205 }
206 else {
207 pixel[1] = pval(arast, col - 1);
208 pixel[2] = pval(arast, col);
209 pixel[3] = -1;
210 pixel[4] = pval(brast, col - 1);
211 pixel[5] = -1;
212 pixel[6] = pval(crast, col - 1);
213 pixel[7] = pval(crast, col);
214 pixel[8] = -1;
215 }
216 }
217 else {
218 pixel[6] = -1;
219 pixel[7] = -1;
220 pixel[8] = -1;
221 if (col == 0) {
222 pixel[1] = -1;
223 pixel[2] = pval(arast, col);
224 pixel[3] = pval(arast, col + 1);
225 pixel[4] = -1;
226 pixel[5] = pval(brast, col + 1);
227 }
228 else if (col != (ncols - 1)) {
229 pixel[1] = pval(arast, col - 1);
230 pixel[2] = pval(arast, col);
231 pixel[3] = pval(arast, col + 1);
232 pixel[4] = pval(brast, col - 1);
233 pixel[5] = pval(brast, col + 1);
234 }
235 else {
236 pixel[1] = pval(arast, col - 1);
237 pixel[2] = pval(arast, col);
238 pixel[3] = -1;
239 pixel[4] = pval(brast, col - 1);
240 pixel[5] = -1;
241 }
242 }
243
244 cold = warm = shadow = nulo = 0;
245 for (i = 1; i < 9; i++) {
246 switch (pixel[i]) {
247 case IS_COLD_CLOUD:
248 cold++;
249 break;
250 case IS_WARM_CLOUD:
251 warm++;
252 break;
253 case IS_SHADOW:
254 shadow++;
255 break;
256 default:
257 nulo++;
258 break;
259 }
260 }
261 lim = (int)(cold + warm + shadow + nulo) / 2;
262
263 /* Entra pixel[0] = 0 */
264 if (nulo < lim) {
265 if (shadow >= (cold + warm))
266 pixel[0] = IS_SHADOW;
267 else
268 pixel[0] =
269 (warm > cold) ? IS_WARM_CLOUD : IS_COLD_CLOUD;
270 }
271 }
272 if (pixel[0] != 0) {
273 ((CELL *) tmp.rast)[col] = pixel[0];
274 }
275 else {
276 Rast_set_c_null_value((CELL *) tmp.rast + col, 1);
277 }
278 }
279 Rast_put_row(tmp.fd, tmp.rast, CELL_TYPE);
280 }
281 G_percent(1, 1, 1);
282
283 G_free(arast);
284 G_free(brast);
285 G_free(crast);
286 Rast_close(out->fd);
287
288 G_free(tmp.rast);
289 Rast_close(tmp.fd);
290
291 G_remove("cats", out->name);
292 G_remove("cell", out->name);
293 G_remove("cellhd", out->name);
294 G_remove("cell_misc", out->name);
295 G_remove("hist", out->name);
296
297 G_rename("cats", tmp.name, out->name);
298 G_rename("cell", tmp.name, out->name);
299 G_rename("cellhd", tmp.name, out->name);
300 G_rename("cell_misc", tmp.name, out->name);
301 G_rename("hist", tmp.name, out->name);
302
303 return;
304 }
305