GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
iscatt_core.c
Go to the documentation of this file.
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  */
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  */
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 interaction exists
117  \return -1 if regions does not intersect
118  */
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 interaction exists
172  \return -1 if regions do not intersect
173  \return -2 resolution of regions is not same
174  */
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 (fabs(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 (fabs(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  */
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  */
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  */
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 != (row_size) / sizeof(unsigned char)) {
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) {
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  */
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  */
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  */
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  */
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  */
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 
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 
1001  G_free(perimeter.points);
1002  return 0;
1003 }
void Rast_get_null_value_row(int, char *, int)
Read or simulate null value row.
#define CELL_TYPE
Definition: raster.h:11
#define G_malloc(n)
Definition: defs/gis.h:112
int 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)
Compute scatter plots data.
Definition: iscatt_core.c:695
void Rast_set_null_value(void *, int, RASTER_MAP_TYPE)
To set one or more raster values to null.
Definition: null_val.c:98
int I_rasterize(double *polygon, int pol_n_pts, unsigned char val, struct Cell_head *rast_region, unsigned char *rast)
Wrapper for using of iclass perimeter rasterization by scatter plot. Warning: calls Rast_set_window...
Definition: iscatt_core.c:948
int make_perimeter(struct line_pnts *points, IClass_perimeter *perimeter, struct Cell_head *band_region)
Creates one perimeter from vector area.
2D/3D raster map header (used also for region)
Definition: gis.h:423
#define GRASS_EPSILON
Definition: gis.h:159
double west
Extent coordinates (west)
Definition: gis.h:475
CELL max
Definition: raster.h:228
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:149
#define NULL
Definition: ccmath.h:32
int * scatts_bands
Definition: imagery.h:177
int I_insert_patch_to_cat_rast(const char *patch_rast, struct Cell_head *cat_rast_region, const char *cat_rast)
Insert raster map patch into pgm file.
Definition: iscatt_core.c:236
int Rast_read_range(const char *, const char *, struct Range *)
Read raster range (CELL)
Definition: range.c:160
int I_apply_colormap(unsigned char *vals, unsigned char *vals_mask, unsigned nvals, unsigned char *colmap, unsigned char *col_vals)
Apply colromap to the raster.
Definition: iscatt_core.c:905
void Rast_set_window(struct Cell_head *)
Establishes &#39;window&#39; as the current working window.
CELL min
Definition: raster.h:227
void Rast_get_c_row(int, CELL *, int)
Get raster row (CELL type)
Feature geometry info - coordinates.
Definition: dig_structs.h:1675
double north
Extent coordinates (north)
Definition: gis.h:469
int n_a_cats
Definition: imagery.h:161
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
int n_a_scatts
Definition: imagery.h:175
int Rast_open_new(const char *, RASTER_MAP_TYPE)
Opens a new raster map.
Definition: raster/open.c:997
double south
Extent coordinates (south)
Definition: gis.h:471
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
Definition: line.c:149
struct scdScattData ** scatts_arr
Definition: imagery.h:183
int Rast_open_old(const char *, const char *)
Open an existing integer raster map (cell)
Definition: raster/open.c:112
int n_bands
Definition: imagery.h:158
int * cats_ids
Definition: imagery.h:162
int I_create_cat_rast(struct Cell_head *cat_rast_region, const char *cat_rast)
Create category raster conditions file. The file is used for holding selected areas from mapwindow...
Definition: iscatt_core.c:58
Definition: raster.h:225
int * cats_idxs
Definition: imagery.h:164
unsigned int * scatt_vals_arr
Definition: imagery.h:195
const char * G_find_raster(char *, const char *)
Find a raster map.
Definition: find_rast.c:55
double * y
Array of Y coordinates.
Definition: dig_structs.h:1684
int cols
Number of columns for 2D data.
Definition: gis.h:442
void Rast_get_cellhd(const char *, const char *, struct Cell_head *)
Read the raster header.
Definition: get_cellhd.c:41
void Rast_put_c_row(int, const CELL *)
Writes the next row for cell file (CELL version)
double ns_res
Resolution - north to south cell size for 2D data.
Definition: gis.h:463
int Rast_window_rows(void)
Number of rows in active window.
Definition: raster/window.c:85
void G_zero(void *, int)
Zero out a buffer, buf, of length i.
Definition: gis/zero.c:23
int Rast_window_cols(void)
Number of columns in active window.
void G_warning(const char *,...) __attribute__((format(printf
int CELL
Definition: gis.h:613
struct scScatts ** cats_arr
Definition: imagery.h:167
RASTER_MAP_TYPE Rast_get_map_type(int)
Determine raster type from descriptor.
Definition: raster/open.c:918
double east
Extent coordinates (east)
Definition: gis.h:473
#define _(str)
Definition: glocale.h:10
int RASTER_MAP_TYPE
Definition: raster.h:25
char * Rast_allocate_null_buf(void)
Allocates memory for a null buffer.
Definition: alloc_cell.c:121
const char * G_find_raster2(const char *, const char *)
Find a raster map (look but don&#39;t touch)
Definition: find_rast.c:76
double ew_res
Resolution - east to west cell size for 2D data.
Definition: gis.h:459
unsigned char * b_conds_arr
Definition: imagery.h:193
void Vect_destroy_line_struct(struct line_pnts *)
Frees all memory associated with a line_pnts structure, including the structure itself.
Definition: line.c:77
int rows
Number of rows for 2D data.
Definition: gis.h:438
int G_debug(int, const char *,...) __attribute__((format(printf
CELL * Rast_allocate_c_buf(void)
Allocate memory for a CELL type raster map.
Definition: alloc_cell.c:82
int I_merge_arrays(unsigned char *merged_arr, unsigned char *overlay_arr, unsigned rows, unsigned cols, double alpha)
Merge arrays according to opacity. Every pixel in array must be represented by 4 values (RGBA)...
Definition: iscatt_core.c:862
void Rast_close(int)
Close a raster map.
Definition: raster/close.c:99