1 /*!
2 \file lib/imagery/iscatt_core.c
3
4 \brief Imagery library - wx.iscatt (wx Interactive Scatter Plot Tool) backend.
5
6 Copyright (C) 2013 by the GRASS Development Team
7
8 This program is free software under the GNU General Public License
9 (>=v2). Read the file COPYING that comes with GRASS for details.
10
11 \author Stepan Turek <stepan.turek@seznam.cz> (GSoC 2013, Mentor: Martin Landa)
12 */
13
14 #include <stdio.h>
15 #include <string.h>
16 #include <math.h>
17
18 #include <grass/gis.h>
19 #include <grass/vector.h>
20 #include <grass/raster.h>
21 #include <grass/imagery.h>
22 #include <grass/glocale.h>
23
24 #include "iclass_local_proto.h"
25
26 struct rast_row
27 {
28 CELL *row;
29 char *null_row;
30 struct Range rast_range; /*Range of whole raster. */
31 };
32
33 /*!
34 \brief Create pgm header.
35
36 Scatter plot internally generates pgm files.
37 These pgms have header in format created by this function.
38
39 \param region region to be pgm header generated for
40 \param [out] header header of pgm file
41 */
get_cat_rast_header(struct Cell_head * region,char * header)42 static int get_cat_rast_header(struct Cell_head *region, char *header)
43 {
44 return sprintf(header, "P5\n%d\n%d\n1\n", region->cols, region->rows);
45 }
46
47 /*!
48 \brief Create category raster conditions file.
49 The file is used for holding selected areas from mapwindow.
50 The reason why the standard GRASS raster is not used is that for every
51 modification (new area is selected) we would need to create whole new raster.
52 Instead of this scatter plot only modifies affected region in
53 the internal pgm file.
54
55 \param cat_rast_region region to be file generated for
56 \param[out] cat_rast path where the pgm raster file will be generated
57 */
I_create_cat_rast(struct Cell_head * cat_rast_region,const char * cat_rast)58 int I_create_cat_rast(struct Cell_head *cat_rast_region, const char *cat_rast)
59 {
60 FILE *f_cat_rast;
61 char cat_rast_header[1024]; /* TODO magic number */
62 int i_row, i_col;
63 int head_nchars;
64
65 unsigned char *row_data;
66
67 f_cat_rast = fopen(cat_rast, "wb");
68 if (!f_cat_rast) {
69 G_warning("Unable to create category raster condition file <%s>.",
70 cat_rast);
71 return -1;
72 }
73
74 head_nchars = get_cat_rast_header(cat_rast_region, cat_rast_header);
75
76 fwrite(cat_rast_header, sizeof(char), head_nchars / sizeof(char),
77 f_cat_rast);
78 if (ferror(f_cat_rast)) {
79 fclose(f_cat_rast);
80 G_warning(_
81 ("Unable to write header into category raster condition file <%s>."),
82 cat_rast);
83 return -1;
84 }
85
86 row_data =
87 (unsigned char *)G_malloc(cat_rast_region->cols *
88 sizeof(unsigned char));
89 for (i_col = 0; i_col < cat_rast_region->cols; i_col++)
90 row_data[i_col] = 0 & 255;
91
92 for (i_row = 0; i_row < cat_rast_region->rows; i_row++) {
93 fwrite(row_data, sizeof(unsigned char),
94 (cat_rast_region->cols) / sizeof(unsigned char), f_cat_rast);
95 if (ferror(f_cat_rast)) {
96 fclose(f_cat_rast);
97 G_warning(_
98 ("Unable to write into category raster condition file <%s>."),
99 cat_rast);
100 return -1;
101 }
102 }
103
104 fclose(f_cat_rast);
105 return 0;
106 }
107
108 /*!
109 \brief Find intersection region of two regions.
110
111 \param A pointer to intersected region
112 \param B pointer to intersected region
113 \param [out] intersec pointer to intersection region of regions A B
114 (relevant params of the region are: south, north, east, west)
115
116 \return 0 if interection exists
117 \return -1 if regions does not intersect
118 */
regions_intersecion(struct Cell_head * A,struct Cell_head * B,struct Cell_head * intersec)119 static int regions_intersecion(struct Cell_head *A, struct Cell_head *B,
120 struct Cell_head *intersec)
121 {
122
123 if (B->north < A->south)
124 return -1;
125 else if (B->north > A->north)
126 intersec->north = A->north;
127 else
128 intersec->north = B->north;
129
130 if (B->south > A->north)
131 return -1;
132 else if (B->south < A->south)
133 intersec->south = A->south;
134 else
135 intersec->south = B->south;
136
137 if (B->east < A->west)
138 return -1;
139 else if (B->east > A->east)
140 intersec->east = A->east;
141 else
142 intersec->east = B->east;
143
144 if (B->west > A->east)
145 return -1;
146 else if (B->west < A->west)
147 intersec->west = A->west;
148 else
149 intersec->west = B->west;
150
151 if (intersec->north == intersec->south)
152 return -1;
153
154 if (intersec->east == intersec->west)
155 return -1;
156
157 return 0;
158
159 }
160
161 /*!
162 \brief Get rows and cols numbers, which defines intersection of the regions.
163
164 \param A pointer to intersected region
165 \param B pointer to intersected region (A and B must have same resolution)
166 \param [out] A_bounds rows and cols numbers of A stored in
167 south, north, east, west, which defines intersection of A and B
168 \param [out] B_bounds rows and cols numbers of B stored in
169 south, north, east, west, which defines intersection of A and B
170
171 \return 0 if interection exists
172 \return -1 if regions do not intersect
173 \return -2 resolution of regions is not same
174 */
get_rows_and_cols_bounds(struct Cell_head * A,struct Cell_head * B,struct Cell_head * A_bounds,struct Cell_head * B_bounds)175 static int get_rows_and_cols_bounds(struct Cell_head *A, struct Cell_head *B,
176 struct Cell_head *A_bounds,
177 struct Cell_head *B_bounds)
178 {
179 float ns_res, ew_res;
180
181 struct Cell_head intersec;
182
183 /* TODO is it right check? */
184 if (abs(A->ns_res - B->ns_res) > GRASS_EPSILON) {
185 G_warning(
186 "'get_rows_and_cols_bounds' ns_res does not fit, A->ns_res: %f B->ns_res: %f",
187 A->ns_res, B->ns_res);
188 return -2;
189 }
190
191 if (abs(A->ew_res - B->ew_res) > GRASS_EPSILON) {
192 G_warning(
193 "'get_rows_and_cols_bounds' ew_res does not fit, A->ew_res: %f B->ew_res: %f",
194 A->ew_res, B->ew_res);
195 return -2;
196 }
197
198 ns_res = A->ns_res;
199 ew_res = A->ew_res;
200
201 if (regions_intersecion(A, B, &intersec) == -1)
202 return -1;
203
204 A_bounds->north =
205 ceil((A->north - intersec.north - ns_res * 0.5) / ns_res);
206 A_bounds->south =
207 ceil((A->north - intersec.south - ns_res * 0.5) / ns_res);
208
209 A_bounds->east = ceil((intersec.east - A->west - ew_res * 0.5) / ew_res);
210 A_bounds->west = ceil((intersec.west - A->west - ew_res * 0.5) / ew_res);
211
212 B_bounds->north =
213 ceil((B->north - intersec.north - ns_res * 0.5) / ns_res);
214 B_bounds->south =
215 ceil((B->north - intersec.south - ns_res * 0.5) / ns_res);
216
217 B_bounds->east = ceil((intersec.east - B->west - ew_res * 0.5) / ew_res);
218 B_bounds->west = ceil((intersec.west - B->west - ew_res * 0.5) / ew_res);
219
220 return 0;
221 }
222
223 /*!
224 \brief Insert raster map patch into pgm file.
225 \see I_create_cat_rast
226
227 Warning: calls Rast_set_window
228
229 \param patch_rast name of raster map
230 \param cat_rast_region region of category raster file
231 \param cat_rast path to category raster file
232
233 \return 0 on success
234 \return -1 on failure
235 */
I_insert_patch_to_cat_rast(const char * patch_rast,struct Cell_head * cat_rast_region,const char * cat_rast)236 int I_insert_patch_to_cat_rast(const char *patch_rast,
237 struct Cell_head *cat_rast_region,
238 const char *cat_rast)
239 {
240
241 FILE *f_cat_rast;
242 struct Cell_head patch_region, patch_bounds, cat_rast_bounds;
243 char cat_rast_header[1024];
244 int i_row, i_col, ncols, nrows, patch_col;
245 int head_nchars, ret;
246 int fd_patch_rast, init_shift, step_shift;
247 unsigned char *patch_data;
248
249 char *null_chunk_row;
250
251 const char *mapset;
252
253 f_cat_rast = fopen(cat_rast, "rb+");
254 if (!f_cat_rast) {
255 G_warning(_("Unable to open category raster conditions file <%s>."),
256 cat_rast);
257 return -1;
258 }
259
260 head_nchars = get_cat_rast_header(cat_rast_region, cat_rast_header);
261 if ((mapset = G_find_raster((char *)patch_rast, "")) == NULL) {
262 fclose(f_cat_rast);
263 G_warning(_("Unable to find patch raster <%s>."), patch_rast);
264 return -1;
265 }
266
267 Rast_get_cellhd(patch_rast, mapset, &patch_region);
268 Rast_set_window(&patch_region);
269
270 if ((fd_patch_rast = Rast_open_old(patch_rast, mapset)) < 0) {
271 fclose(f_cat_rast);
272 return -1;
273 }
274
275 ret =
276 get_rows_and_cols_bounds(cat_rast_region, &patch_region,
277 &cat_rast_bounds, &patch_bounds);
278 if (ret == -2) {
279 G_warning(_
280 ("Resolutions of patch <%s> and patched file <%s> are not same."),
281 patch_rast, cat_rast);
282
283 Rast_close(fd_patch_rast);
284 fclose(f_cat_rast);
285
286 return -1;
287 }
288 else if (ret == -1) {
289
290 Rast_close(fd_patch_rast);
291 fclose(f_cat_rast);
292
293 return 0;
294 }
295
296 ncols = cat_rast_bounds.east - cat_rast_bounds.west;
297 nrows = cat_rast_bounds.south - cat_rast_bounds.north;
298
299 patch_data = (unsigned char *)G_malloc(ncols * sizeof(unsigned char));
300
301 init_shift =
302 head_nchars + cat_rast_region->cols * cat_rast_bounds.north +
303 cat_rast_bounds.west;
304
305 if (fseek(f_cat_rast, init_shift, SEEK_SET) != 0) {
306 G_warning(_
307 ("Corrupted category raster conditions file <%s> (fseek failed)"),
308 cat_rast);
309
310 Rast_close(fd_patch_rast);
311 fclose(f_cat_rast);
312
313 return -1;
314 }
315
316 step_shift = cat_rast_region->cols - ncols;
317
318 null_chunk_row = Rast_allocate_null_buf();
319
320 for (i_row = 0; i_row < nrows; i_row++) {
321 Rast_get_null_value_row(fd_patch_rast, null_chunk_row,
322 i_row + patch_bounds.north);
323
324 for (i_col = 0; i_col < ncols; i_col++) {
325 patch_col = patch_bounds.west + i_col;
326
327 if (null_chunk_row[patch_col] != 1)
328 patch_data[i_col] = 1 & 255;
329 else {
330 patch_data[i_col] = 0 & 255;
331 }
332 }
333
334 fwrite(patch_data, sizeof(unsigned char),
335 (ncols) / sizeof(unsigned char), f_cat_rast);
336 if (ferror(f_cat_rast)) {
337 G_warning(_
338 ("Unable to write into category raster conditions file <%s>"),
339 cat_rast);
340
341 Rast_close(fd_patch_rast);
342 G_free(null_chunk_row);
343 fclose(f_cat_rast);
344
345 return -1;
346 }
347 if (fseek(f_cat_rast, step_shift, SEEK_CUR) != 0) {
348 G_warning(_
349 ("Corrupted category raster conditions file <%s> (fseek failed)"),
350 cat_rast);
351
352 Rast_close(fd_patch_rast);
353 G_free(null_chunk_row);
354 fclose(f_cat_rast);
355
356 return -1;
357 }
358 }
359
360 Rast_close(fd_patch_rast);
361 G_free(null_chunk_row);
362 fclose(f_cat_rast);
363 return 0;
364 }
365
366 /*!
367 \brief Updates scatter plots data in category by pixels which meets category conditions.
368
369 \param bands_rows data represents data describig one row from raster band
370 \param belongs_pix array which defines which pixels belongs to category
371 (1 value) and which not (0 value)
372 \param [out] scatts pointer to scScatts struct of type SC_SCATT_DATA,
373 which are modified according to values in belongs_pix
374 (represents scatter plot category)
375 */
update_cat_scatt_plts(struct rast_row * bands_rows,unsigned short * belongs_pix,struct scScatts * scatts)376 static void update_cat_scatt_plts(struct rast_row *bands_rows,
377 unsigned short *belongs_pix,
378 struct scScatts *scatts)
379 {
380 int i_scatt, array_idx, i_chunk_rows_pix, max_arr_idx;
381
382 CELL *b_1_row;
383 CELL *b_2_row;
384 char *b_1_null_row, *b_2_null_row;
385 struct rast_row b_1_rast_row, b_2_rast_row;
386
387 struct Range b_1_range, b_2_range;
388 int b_1_range_size;
389
390 int row_size = Rast_window_cols();
391
392 int *scatts_bands = scatts->scatts_bands;
393
394 for (i_scatt = 0; i_scatt < scatts->n_a_scatts; i_scatt++) {
395 b_1_rast_row = bands_rows[scatts_bands[i_scatt * 2]];
396 b_2_rast_row = bands_rows[scatts_bands[i_scatt * 2 + 1]];
397
398 b_1_row = b_1_rast_row.row;
399 b_2_row = b_2_rast_row.row;
400
401 b_1_null_row = b_1_rast_row.null_row;
402 b_2_null_row = b_2_rast_row.null_row;
403
404 b_1_range = b_1_rast_row.rast_range;
405 b_2_range = b_2_rast_row.rast_range;
406
407 b_1_range_size = b_1_range.max - b_1_range.min + 1;
408 max_arr_idx =
409 (b_1_range.max - b_1_range.min + 1) * (b_2_range.max -
410 b_2_range.min + 1);
411
412 for (i_chunk_rows_pix = 0; i_chunk_rows_pix < row_size;
413 i_chunk_rows_pix++) {
414 /* pixel does not belongs to scatter plot or has null value in one of the bands */
415 if (!belongs_pix[i_chunk_rows_pix] ||
416 b_1_null_row[i_chunk_rows_pix] == 1 ||
417 b_2_null_row[i_chunk_rows_pix] == 1)
418 continue;
419
420 /* index in scatter plot array */
421 array_idx =
422 b_1_row[i_chunk_rows_pix] - b_1_range.min +
423 (b_2_row[i_chunk_rows_pix] - b_2_range.min) * b_1_range_size;
424
425 if (array_idx < 0 || array_idx >= max_arr_idx) {
426 G_warning
427 ("Inconsistent data. Value computed for scatter plot is out of initialized range.");
428 continue;
429 }
430
431 /* increment scatter plot value */
432 ++scatts->scatts_arr[i_scatt]->scatt_vals_arr[array_idx];
433 }
434 }
435 }
436
437 /*!
438 \brief Computes scatter plots data from bands_rows.
439
440 \param scatt_conds pointer to scScatts struct of type SC_SCATT_CONDITIONS,
441 where are selected areas (conditions) stored
442 \param f_cats_rasts_conds file which stores selected areas (conditions) from
443 mapwindow see I_create_cat_rast and I_insert_patch_to_cat_rast
444 \param bands_rows data arrays of raster rows from analyzed raster bands
445 (all data in bands_rows and belongs_pix arrays represents same region (row))
446 \param[out] scatts pointer to scScatts struct of type SC_SCATT_DATA,
447 where are computed scatter plots stored
448 \param[out] fd_cats_rasts array of opened raster maps,
449 which every represents all selected pixels for category
450
451 \return 0 on success
452 \return -1 on failure
453 */
compute_scatts_from_chunk_row(struct scCats * scatt_conds,FILE ** f_cats_rasts_conds,struct rast_row * bands_rows,struct scCats * scatts,int * fd_cats_rasts)454 static int compute_scatts_from_chunk_row(struct scCats *scatt_conds,
455 FILE ** f_cats_rasts_conds,
456 struct rast_row *bands_rows,
457 struct scCats *scatts,
458 int *fd_cats_rasts)
459 {
460
461 int i_rows_pix, i_cat, i_scatt, n_pixs;
462 int cat_id, scatt_plts_cat_idx, array_idx, max_arr_idx;
463 char *b_1_null_row, *b_2_null_row;
464 struct rast_row b_1_rast_row, b_2_rast_row;
465 CELL *cat_rast_row;
466
467 struct scScatts *scatts_conds;
468 struct scScatts *scatts_scatt_plts;
469
470 struct Range b_1_range, b_2_range;
471 int b_1_range_size;
472
473 int *scatts_bands;
474
475 CELL *b_1_row;
476 CELL *b_2_row;
477 unsigned char *i_scatt_conds;
478
479 int row_size = Rast_window_cols();
480
481 unsigned short *belongs_pix =
482 (unsigned short *)G_malloc(row_size * sizeof(unsigned short));
483 unsigned char *rast_pixs =
484 (unsigned char *)G_malloc(row_size * sizeof(unsigned char));
485 cat_rast_row = Rast_allocate_c_buf();
486
487
488 for (i_cat = 0; i_cat < scatt_conds->n_a_cats; i_cat++) {
489 scatts_conds = scatt_conds->cats_arr[i_cat];
490
491 cat_id = scatt_conds->cats_ids[i_cat];
492
493 scatt_plts_cat_idx = scatts->cats_idxs[cat_id];
494 if (scatt_plts_cat_idx < 0)
495 continue;
496 scatts_scatt_plts = scatts->cats_arr[scatt_plts_cat_idx];
497
498 G_zero(belongs_pix, row_size * sizeof(unsigned short));
499
500 /* if category has no conditions defined, scatter plots without
501 any constraint are computed (default scatter plots) */
502 if (!scatts_conds->n_a_scatts && !f_cats_rasts_conds[i_cat]) {
503 for (i_scatt = 0; i_scatt < scatts_scatt_plts->n_a_scatts;
504 i_scatt++) {
505 /* all pixels belongs */
506 for (i_rows_pix = 0; i_rows_pix < row_size; i_rows_pix++)
507 belongs_pix[i_rows_pix] = 1;
508 }
509 }
510 /* compute belonging pixels for defined conditions */
511 else {
512 scatts_bands = scatts_conds->scatts_bands;
513
514 /* check conditions from category raster conditions file
515 (see I_create_cat_rast) */
516 if (f_cats_rasts_conds[i_cat]) {
517 n_pixs =
518 fread(rast_pixs, sizeof(unsigned char),
519 (row_size) / sizeof(unsigned char),
520 f_cats_rasts_conds[i_cat]);
521
522 if (ferror(f_cats_rasts_conds[i_cat])) {
523 G_free(rast_pixs);
524 G_free(belongs_pix);
525 G_warning(_
526 ("Unable to read from category raster condition file."));
527 return -1;
528 }
529 if (n_pixs != n_pixs) {
530 G_free(rast_pixs);
531 G_free(belongs_pix);
532 G_warning(_
533 ("Invalid size of category raster conditions file."));
534 return -1;
535
536 }
537
538 for (i_rows_pix = 0; i_rows_pix < row_size; i_rows_pix++) {
539 if (rast_pixs[i_rows_pix] != (0 & 255))
540 belongs_pix[i_rows_pix] = 1;
541 }
542 }
543
544 /* check conditions defined in scatter plots */
545 for (i_scatt = 0; i_scatt < scatts_conds->n_a_scatts; i_scatt++) {
546 b_1_rast_row = bands_rows[scatts_bands[i_scatt * 2]];
547 b_2_rast_row = bands_rows[scatts_bands[i_scatt * 2 + 1]];
548
549 b_1_row = b_1_rast_row.row;
550 b_2_row = b_2_rast_row.row;
551
552 b_1_null_row = b_1_rast_row.null_row;
553 b_2_null_row = b_2_rast_row.null_row;
554
555 b_1_range = b_1_rast_row.rast_range;
556 b_2_range = b_2_rast_row.rast_range;
557
558 b_1_range_size = b_1_range.max - b_1_range.min + 1;
559 max_arr_idx =
560 (b_1_range.max - b_1_range.min + 1) * (b_2_range.max -
561 b_2_range.min + 1);
562
563 i_scatt_conds =
564 scatts_conds->scatts_arr[i_scatt]->b_conds_arr;
565
566 for (i_rows_pix = 0; i_rows_pix < row_size; i_rows_pix++) {
567 /* pixels already belongs to category from category raster conditions
568 file or contains null value in one of the bands */
569 if (belongs_pix[i_rows_pix] ||
570 b_1_null_row[i_rows_pix] == 1 ||
571 b_2_null_row[i_rows_pix] == 1)
572 continue;
573
574 array_idx =
575 b_1_row[i_rows_pix] - b_1_range.min +
576 (b_2_row[i_rows_pix] -
577 b_2_range.min) * b_1_range_size;
578 if (array_idx < 0 || array_idx >= max_arr_idx) {
579 G_warning(_("Data inconsistent. "
580 "Value computed for scatter plot is out of initialized range."));
581 continue;
582 }
583 /* pixels meets condtion defined in scatter plot ->
584 belongs to scatter plot category */
585 if (i_scatt_conds[array_idx])
586 belongs_pix[i_rows_pix] = 1;
587 }
588 }
589 }
590
591 /* update category raster with belonging pixels */
592 if (fd_cats_rasts[i_cat] >= 0) {
593 Rast_set_null_value(cat_rast_row, Rast_window_cols(), CELL_TYPE);
594
595 for (i_rows_pix = 0; i_rows_pix < row_size; i_rows_pix++)
596 if (belongs_pix[i_rows_pix])
597 cat_rast_row[i_rows_pix] = belongs_pix[i_rows_pix];
598
599 Rast_put_c_row(fd_cats_rasts[i_cat], cat_rast_row);
600 }
601
602 /* update scatter plots with belonging pixels */
603 update_cat_scatt_plts(bands_rows, belongs_pix, scatts_scatt_plts);
604 }
605
606 G_free(cat_rast_row);
607 G_free(rast_pixs);
608 G_free(belongs_pix);
609
610 return 0;
611 }
612
613 /*!
614 \brief Get list of bands needed to be opened for analysis from scCats struct.
615 */
get_needed_bands(struct scCats * cats,int * b_needed_bands)616 static void get_needed_bands(struct scCats *cats, int *b_needed_bands)
617 {
618 /* results in b_needed_bands - array of bools - if item has value 1,
619 band (defined by item index) is needed to be opened */
620 int i_cat, i_scatt;
621
622 for (i_cat = 0; i_cat < cats->n_a_cats; i_cat++) {
623 for (i_scatt = 0; i_scatt < cats->cats_arr[i_cat]->n_a_scatts;
624 i_scatt++) {
625 G_debug(3, "Active scatt %d in catt %d", i_scatt, i_cat);
626
627 b_needed_bands[cats->cats_arr[i_cat]->scatts_bands[i_scatt * 2]] =
628 1;
629 b_needed_bands[cats->cats_arr[i_cat]->
630 scatts_bands[i_scatt * 2 + 1]] = 1;
631 }
632 }
633 return;
634 }
635
636 /*!
637 \brief Helper function for clean up.
638 */
free_compute_scatts_data(int * fd_bands,struct rast_row * bands_rows,int n_a_bands,int * bands_ids,int * fd_cats_rasts,FILE ** f_cats_rasts_conds,int n_a_cats)639 static void free_compute_scatts_data(int *fd_bands,
640 struct rast_row *bands_rows,
641 int n_a_bands, int *bands_ids,
642 int *fd_cats_rasts,
643 FILE ** f_cats_rasts_conds, int n_a_cats)
644 {
645 int i, band_id;
646
647 for (i = 0; i < n_a_bands; i++) {
648 band_id = bands_ids[i];
649 if (band_id >= 0) {
650 Rast_close(fd_bands[i]);
651 G_free(bands_rows[band_id].row);
652 G_free(bands_rows[band_id].null_row);
653 }
654 }
655
656 if (f_cats_rasts_conds)
657 for (i = 0; i < n_a_cats; i++)
658 if (f_cats_rasts_conds[i])
659 fclose(f_cats_rasts_conds[i]);
660
661 if (fd_cats_rasts)
662 for (i = 0; i < n_a_cats; i++)
663 if (fd_cats_rasts[i] >= 0)
664 Rast_close(fd_cats_rasts[i]);
665
666 }
667
668 /*!
669 \brief Compute scatter plots data.
670
671 If category has not defined category raster condition file and no scatter plot
672 exists with condition, default/full scatter plot is computed.
673 Warning: calls Rast_set_window
674
675 \param region analysis region, beaware that all input data must be prepared for this region
676 (bands (their ranges), cats_rasts_conds rasters...)
677 \param region function calls Rast_set_window for this region
678 \param scatt_conds pointer to scScatts struct of type SC_SCATT_CONDITIONS,
679 where are stored selected areas (conditions) in scatter plots
680 \param cats_rasts_conds paths to category raster conditions files representing
681 selected areas from mapwindow (conditions) in rasters for every category
682 \param cats_rasts_conds index in array represents corresponding category id
683 \param cats_rasts_conds for manipulation with category raster conditions file
684 see also I_id_scatt_to_bands and I_insert_patch_to_cat_rast
685 \param bands names of analyzed bands, order of bands is defined by their id
686 \param n_bands number of bands
687 \param[out] scatts pointer to scScatts struct of type SC_SCATT_DATA,
688 where are computed scatter plots stored
689 \param[out] cats_rasts array of raster maps names for every category
690 where will be stored all selected pixels
691
692 \return 0 on success
693 \return -1 on failure
694 */
I_compute_scatts(struct Cell_head * region,struct scCats * scatt_conds,const char ** cats_rasts_conds,const char ** bands,int n_bands,struct scCats * scatts,const char ** cats_rasts)695 int I_compute_scatts(struct Cell_head *region, struct scCats *scatt_conds,
696 const char **cats_rasts_conds, const char **bands,
697 int n_bands, struct scCats *scatts,
698 const char **cats_rasts)
699 {
700 const char *mapset;
701 char header[1024];
702
703 int fd_cats_rasts[scatt_conds->n_a_cats];
704 FILE *f_cats_rasts_conds[scatt_conds->n_a_cats];
705
706 struct rast_row bands_rows[n_bands];
707
708 RASTER_MAP_TYPE data_type;
709
710 int nrows, i_band, n_a_bands, band_id;
711 int i_row, head_nchars, i_cat, id_cat;
712
713 int fd_bands[n_bands];
714 int bands_ids[n_bands];
715 int b_needed_bands[n_bands];
716
717 Rast_set_window(region);
718
719 for (i_band = 0; i_band < n_bands; i_band++)
720 fd_bands[i_band] = -1;
721
722 for (i_band = 0; i_band < n_bands; i_band++)
723 bands_ids[i_band] = -1;
724
725 if (n_bands != scatts->n_bands || n_bands != scatt_conds->n_bands)
726 return -1;
727
728 G_zero(b_needed_bands, (size_t) n_bands * sizeof(int));
729
730 get_needed_bands(scatt_conds, &b_needed_bands[0]);
731 get_needed_bands(scatts, &b_needed_bands[0]);
732
733 n_a_bands = 0;
734
735 /* open band rasters, which are needed for computation */
736 for (band_id = 0; band_id < n_bands; band_id++) {
737 if (b_needed_bands[band_id]) {
738 G_debug(3, "Opening raster no. %d with name: %s", band_id,
739 bands[band_id]);
740
741 if ((mapset = G_find_raster2(bands[band_id], "")) == NULL) {
742 free_compute_scatts_data(fd_bands, bands_rows, n_a_bands,
743 bands_ids, NULL, NULL,
744 scatt_conds->n_a_cats);
745 G_warning(_("Unable to find raster <%s>"),
746 bands[band_id]);
747 return -1;
748 }
749
750 if ((fd_bands[n_a_bands] =
751 Rast_open_old(bands[band_id], mapset)) < 0) {
752 free_compute_scatts_data(fd_bands, bands_rows, n_a_bands,
753 bands_ids, NULL, NULL,
754 scatt_conds->n_a_cats);
755 G_warning(_("Unable to open raster <%s>"), bands[band_id]);
756 return -1;
757 }
758
759 data_type = Rast_get_map_type(fd_bands[n_a_bands]);
760 if (data_type != CELL_TYPE) {
761 G_warning(_("Raster <%s> type is not <%s>"), bands[band_id],
762 "CELL");
763 return -1;
764 }
765
766 bands_rows[band_id].row = Rast_allocate_c_buf();
767 bands_rows[band_id].null_row = Rast_allocate_null_buf();
768
769 if (Rast_read_range
770 (bands[band_id], mapset,
771 &bands_rows[band_id].rast_range) != 1) {
772 free_compute_scatts_data(fd_bands, bands_rows, n_a_bands,
773 bands_ids, NULL, NULL,
774 scatt_conds->n_a_cats);
775 G_warning(_("Unable to read range of raster <%s>"),
776 bands[band_id]);
777 return -1;
778 }
779
780 bands_ids[n_a_bands] = band_id;
781 ++n_a_bands;
782 }
783 }
784
785 /* open category rasters condition files and category rasters */
786 for (i_cat = 0; i_cat < scatts->n_a_cats; i_cat++) {
787 id_cat = scatts->cats_ids[i_cat];
788 if (cats_rasts[id_cat]) {
789 fd_cats_rasts[i_cat] =
790 Rast_open_new(cats_rasts[id_cat], CELL_TYPE);
791 }
792 else
793 fd_cats_rasts[i_cat] = -1;
794
795 if (cats_rasts_conds[id_cat]) {
796 f_cats_rasts_conds[i_cat] = fopen(cats_rasts_conds[id_cat], "r");
797 if (!f_cats_rasts_conds[i_cat]) {
798 free_compute_scatts_data(fd_bands, bands_rows, n_a_bands,
799 bands_ids, fd_cats_rasts,
800 f_cats_rasts_conds,
801 scatt_conds->n_a_cats);
802 G_warning(_
803 ("Unable to open category raster condition file <%s>"),
804 bands[band_id]);
805 return -1;
806 }
807 }
808 else
809 f_cats_rasts_conds[i_cat] = NULL;
810 }
811
812 head_nchars = get_cat_rast_header(region, header);
813 for (i_cat = 0; i_cat < scatt_conds->n_a_cats; i_cat++)
814 if (f_cats_rasts_conds[i_cat])
815 if (fseek(f_cats_rasts_conds[i_cat], head_nchars, SEEK_SET) != 0) {
816 G_warning(_
817 ("Corrupted category raster conditions file (fseek failed)"));
818 return -1;
819 }
820
821 nrows = Rast_window_rows();
822
823 /* analyze bands by rows */
824 for (i_row = 0; i_row < nrows; i_row++) {
825 for (i_band = 0; i_band < n_a_bands; i_band++) {
826 band_id = bands_ids[i_band];
827 Rast_get_c_row(fd_bands[i_band], bands_rows[band_id].row, i_row);
828 Rast_get_null_value_row(fd_bands[i_band],
829 bands_rows[band_id].null_row, i_row);
830 }
831 if (compute_scatts_from_chunk_row
832 (scatt_conds, f_cats_rasts_conds, bands_rows, scatts,
833 fd_cats_rasts) == -1) {
834 free_compute_scatts_data(fd_bands, bands_rows, n_a_bands,
835 bands_ids, fd_cats_rasts,
836 f_cats_rasts_conds,
837 scatt_conds->n_a_cats);
838 return -1;
839 }
840
841 }
842 free_compute_scatts_data(fd_bands, bands_rows, n_a_bands, bands_ids,
843 fd_cats_rasts, f_cats_rasts_conds,
844 scatt_conds->n_a_cats);
845 return 0;
846 }
847
848 /*!
849 \brief Merge arrays according to opacity.
850 Every pixel in array must be represented by 4 values (RGBA).
851
852 Implementd for speeding up of scatter plots rendering.
853
854 \param merged_arr array which will be overlayd with overlay_arr
855 \param overlay_arr array to be merged_arr overlaid with
856 \param rows number of rows for the both arrays
857 \param cols number of columns for the both arrays
858 \param alpha transparency (0-1) of the overlay array for merging
859
860 \return 0
861 */
I_merge_arrays(unsigned char * merged_arr,unsigned char * overlay_arr,unsigned rows,unsigned cols,double alpha)862 int I_merge_arrays(unsigned char *merged_arr, unsigned char *overlay_arr,
863 unsigned rows, unsigned cols, double alpha)
864 {
865 unsigned int i_row, i_col, i_b;
866 unsigned int row_idx, col_idx, idx;
867 unsigned int c_a_i, c_a;
868
869 for (i_row = 0; i_row < rows; i_row++) {
870 row_idx = i_row * cols;
871 for (i_col = 0; i_col < cols; i_col++) {
872 col_idx = 4 * (row_idx + i_col);
873 idx = col_idx + 3;
874
875 c_a = overlay_arr[idx] * alpha;
876 c_a_i = 255 - c_a;
877
878 merged_arr[idx] =
879 (c_a_i * (int)merged_arr[idx] + c_a * 255) / 255;
880
881 for (i_b = 0; i_b < 3; i_b++) {
882 idx = col_idx + i_b;
883 merged_arr[idx] =
884 (c_a_i * (int)merged_arr[idx] +
885 c_a * (int)overlay_arr[idx]) / 255;
886 }
887 }
888 }
889 return 0;
890 }
891
892 /*!
893 \brief Apply colromap to the raster.
894
895 Implementd for speeding up of scatter plots rendering.
896
897 \param vals array of values for applying the colormap
898 \param vals_mask maks of vals array
899 \param nvals number of items of vals_mask and vals array
900 \param colmap colour map to be applied
901 \param[out] col_vals output raster with applied color map (length is 4 * nvals (RGBA))
902
903 \return 0
904 */
I_apply_colormap(unsigned char * vals,unsigned char * vals_mask,unsigned nvals,unsigned char * colmap,unsigned char * col_vals)905 int I_apply_colormap(unsigned char *vals, unsigned char *vals_mask,
906 unsigned nvals, unsigned char *colmap,
907 unsigned char *col_vals)
908 {
909 unsigned int i_val;
910 int v, i, i_cm;
911
912 for (i_val = 0; i_val < nvals; i_val++) {
913 i_cm = 4 * i_val;
914
915 v = vals[i_val];
916
917 if (vals_mask && vals_mask[i_val])
918 for (i = 0; i < 4; i++)
919 col_vals[i_cm + i] = colmap[258 * 4 + i];
920 else if (v > 255)
921 for (i = 0; i < 4; i++)
922 col_vals[i_cm + i] = colmap[257 * 4 + i];
923 else if (v < 0)
924 for (i = 0; i < 4; i++)
925 col_vals[i_cm + i] = colmap[256 * 4 + i];
926 else
927 for (i = 0; i < 4; i++) {
928 col_vals[i_cm + i] = colmap[v * 4 + i];
929 }
930 }
931 return 0;
932 }
933
934 /*!
935 \brief Wrapper for using of iclass perimeter rasterization by scatter plot.
936 Warning: calls Rast_set_window
937
938 \param polygon array of polygon coordinates [x, y, x, y...]
939 \param pol_n_pts number of points in the polygon array
940 \param val value to be assigned to cells, which belong to plygon
941 \param rast_region region of raster
942 \param[out] rast raster to be pologyn rasterized in
943
944 \return 0 on success
945 \return 1 on failure
946 */
947
I_rasterize(double * polygon,int pol_n_pts,unsigned char val,struct Cell_head * rast_region,unsigned char * rast)948 int I_rasterize(double *polygon, int pol_n_pts, unsigned char val,
949 struct Cell_head *rast_region, unsigned char *rast)
950 {
951 int i;
952 int x0, x1, y;
953 int row, row_idx, i_col;
954
955 IClass_perimeter perimeter;
956
957 struct line_pnts *pol;
958
959 pol = Vect_new_line_struct();
960
961 for (i = 0; i < pol_n_pts; i++) {
962 Vect_append_point(pol, polygon[i * 2], polygon[i * 2 + 1], 0.0);
963 }
964
965 /* Rast_set_window(rast_region); */
966
967 make_perimeter(pol, &perimeter, rast_region);
968 for (i = 1; i < perimeter.npoints; i += 2) {
969 y = perimeter.points[i].y;
970 if (y != perimeter.points[i - 1].y) {
971 G_warning(_
972 ("prepare_signature: scan line %d has odd number of points."),
973 (i + 1) / 2);
974 return 1;
975 }
976
977 x0 = perimeter.points[i - 1].x;
978 x1 = perimeter.points[i].x;
979
980 if (x0 > x1) {
981 G_warning(_("signature: perimeter points out of order."));
982 return 1;
983 }
984
985 row = (rast_region->rows - y);
986 if (row < 0 || row >= rast_region->rows) {
987 continue;
988 }
989
990 row_idx = rast_region->cols * row;
991
992 for (i_col = x0; i_col <= x1; i_col++) {
993 if (i_col < 0 || i_col >= rast_region->cols) {
994 continue;
995 }
996 rast[row_idx + i_col] = val;
997 }
998 }
999
1000 Vect_destroy_line_struct(pol);
1001 G_free(perimeter.points);
1002 return 0;
1003 }
1004