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