GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-7d4cd0c030
iclass_perimeter.c
Go to the documentation of this file.
1 /*!
2  \file lib/imagery/iclass_perimeter.c
3
4  \brief Imagery library - functions for wx.iclass
5
6  Computation based on training areas for supervised classification.
7  Based on i.class module (GRASS 6).
8
9  Vector map with training areas is used to determine corresponding
10  cells by computing cells on area perimeter.
11
12  Copyright (C) 1999-2007, 2011 by the GRASS Development Team
13
14  This program is free software under the GNU General Public License
15  (>=v2). Read the file COPYING that comes with GRASS for details.
16
17  \author David Satnik, Central Washington University (original author)
18  \author Markus Neteler <neteler itc.it> (i.class module)
19  \author Bernhard Reiter <bernhard intevation.de> (i.class module)
21  \author Glynn Clements <glynn gclements.plus.com> (i.class module)
22  \author Hamish Bowman <hamish_b yahoo.com> (i.class module)
23  \author Jan-Oliver Wagner <jan intevation.de> (i.class module)
24  \author Anna Kratochvilova <kratochanna gmail.com> (rewriting for wx.iclass)
25  \author Vaclav Petras <wenzeslaus gmail.com> (rewriting for wx.iclass)
26  */
27
28 #include <stdlib.h>
29
30 #include <grass/vector.h>
31 #include <grass/raster.h>
32 #include <grass/glocale.h>
33
34 #include "iclass_local_proto.h"
35
36 #define extrema(x, y, z) (((x < y) && (z < y)) || ((x > y) && (z > y)))
37 #define non_extrema(x, y, z) (((x < y) && (y < z)) || ((x > y) && (y > z)))
38
39 /*!
40  \brief Creates perimeters from vector areas of given category.
41
42  \param Map vector map
43  \param layer_name layer name (within vector map)
44  \param category vector category (cat column value)
45  \param[out] perimeters list of perimeters
46  \param band_region region which determines perimeter cells
47
48  \return number of areas of given cat
49  \return -1 on error
50  */
51 int vector2perimeters(struct Map_info *Map, const char *layer_name,
52  int category, IClass_perimeter_list *perimeters,
54 {
55  struct line_pnts *points;
56
57  int nareas, nareas_cat, layer;
58
59  int i, cat, ret;
60
61  int j;
62
63  G_debug(3, "iclass_vector2perimeters():layer = %s, category = %d",
64  layer_name, category);
65
66  layer = Vect_get_field_number(Map, layer_name);
67  nareas = Vect_get_num_areas(Map);
68  if (nareas == 0)
69  return 0;
70
71  nareas_cat = 0;
72  /* find out, how many areas have given category */
73  for (i = 1; i <= nareas; i++) {
74  if (!Vect_area_alive(Map, i))
75  continue;
76  cat = Vect_get_area_cat(Map, i, layer);
77  if (cat < 0) {
78  /* no centroid, no category */
79  }
80  else if (cat == category) {
81  nareas_cat++;
82  }
83  }
84  if (nareas_cat == 0)
85  return 0;
86
87  perimeters->nperimeters = nareas_cat;
88  perimeters->perimeters =
89  (IClass_perimeter *)G_calloc(nareas_cat, sizeof(IClass_perimeter));
90
91  j = 0; /* area with cat */
92  for (i = 1; i <= nareas; i++) {
93  if (!Vect_area_alive(Map, i))
94  continue;
95  cat = Vect_get_area_cat(Map, i, layer);
96  if (cat < 0) {
97  /* no centroid, no category */
98  }
99  else if (cat == category) {
100  j++;
101
102  points = Vect_new_line_struct(); /* Vect_destroy_line_struct */
103  ret = Vect_get_area_points(Map, i, points);
104
105  if (ret <= 0) {
106  Vect_destroy_line_struct(points);
107  free_perimeters(perimeters);
108  G_warning(_("Get area %d failed"), i);
109  return -1;
110  }
111  if (make_perimeter(points, &perimeters->perimeters[j - 1],
112  band_region) <= 0) {
113  Vect_destroy_line_struct(points);
114  free_perimeters(perimeters);
115  G_warning(_("Perimeter computation failed"));
116  return -1;
117  }
118  Vect_destroy_line_struct(points);
119  }
120  }
121
122  /* Vect_close(&Map); */
123
124  return nareas_cat;
125 }
126
127 /*!
128  \brief Frees all perimeters in list of perimeters.
129
130  It also frees list of perimeters itself.
131
132  \param perimeters list of perimeters
133  */
134 void free_perimeters(IClass_perimeter_list *perimeters)
135 {
136  int i;
137
138  G_debug(5, "free_perimeters()");
139
140  for (i = 0; i < perimeters->nperimeters; i++) {
141  G_free(perimeters->perimeters[i].points);
142  }
143  G_free(perimeters->perimeters);
144 }
145
146 /*!
147  \brief Creates one perimeter from vector area.
148
149  \param points list of vertices represting area
150  \param[out] perimeter perimeter
151  \param band_region region which determines perimeter cells
152
153  \return 1 on success
154  \return 0 on error
155  */
156 int make_perimeter(struct line_pnts *points, IClass_perimeter *perimeter,
158 {
159  IClass_point *tmp_points;
160
161  IClass_point *vertex_points;
162
163  int i, first, prev, skip, next;
164
165  int count, vertex_count;
166
167  int np; /* perimeter estimate */
168
169  G_debug(5, "iclass_make_perimeter()");
170  count = points->n_points;
171
172  tmp_points =
173  (IClass_point *)G_calloc(count, sizeof(IClass_point)); /* TODO test */
174
175  for (i = 0; i < count; i++) {
176  G_debug(5, "iclass_make_perimeter(): points: x: %f y: %f", points->x[i],
177  points->y[i]);
178
179  /* This functions are no longer used because of the different behavior
180  of Rast_easting_to_col depending whether location is LL or not.
181  It makes problem in interactive scatter plot tool,
182  which defines its own coordinates systems for the plots and
183  therefore it requires the function to work always in same way
184  without hidden dependency on location type.
185
186  tmp_points[i].y = Rast_northing_to_row(points->y[i], band_region);
187  tmp_points[i].x = Rast_easting_to_col(points->x[i], band_region);
188  */
189
190  tmp_points[i].y =
191  (band_region->north - points->y[i]) / band_region->ns_res;
192  tmp_points[i].x =
193  (points->x[i] - band_region->west) / band_region->ew_res;
194  }
195
196  /* find first edge which is not horizontal */
197
198  first = -1;
199  prev = count - 1;
200  for (i = 0; i < count; prev = i++) {
201  /* non absurd polygon has vertexes with different y coordinate */
202  if (tmp_points[i].y != tmp_points[prev].y) {
203  first = i;
204  break;
205  }
206  }
207  if (first < 0) {
208  G_free(tmp_points);
209  G_warning(_("Invalid polygon"));
210  return 0;
211  }
212
213  /* copy tmp to vertex list collapsing adjacent horizontal edges */
214
215  /* vertex_count <= count, size of vertex_points is count */
216  vertex_points =
217  (IClass_point *)G_calloc(count, sizeof(IClass_point)); /* TODO test */
218  skip = 0;
219  vertex_count = 0;
220  i = first; /* stmt not necessary */
221
222  do {
223  if (!skip) {
224  vertex_points[vertex_count].x = tmp_points[i].x;
225  vertex_points[vertex_count].y = tmp_points[i].y;
226  vertex_count++;
227  }
228
229  prev = i++;
230  if (i >= count)
231  i = 0;
232  if ((next = i + 1) >= count)
233  next = 0;
234
235  skip = ((tmp_points[prev].y == tmp_points[i].y) &&
236  (tmp_points[next].y == tmp_points[i].y));
237  } while (i != first);
238
239  G_free(tmp_points);
240
241  /* count points on the perimeter */
242
243  np = 0;
244  prev = vertex_count - 1;
245  for (i = 0; i < vertex_count; prev = i++) {
246  np += abs(vertex_points[prev].y - vertex_points[i].y);
247  }
248
249  /* allocate perimeter list */
250
251  perimeter->points = (IClass_point *)G_calloc(np, sizeof(IClass_point));
252  if (!perimeter->points) {
253  G_free(vertex_points);
254  G_warning(_("Outlined area is too large."));
255  return 0;
256  }
257
258  /* store the perimeter points */
259
260  perimeter->npoints = 0;
261  prev = vertex_count - 1;
262  for (i = 0; i < vertex_count; prev = i++) {
263  edge2perimeter(perimeter, vertex_points[prev].x, vertex_points[prev].y,
264  vertex_points[i].x, vertex_points[i].y);
265  }
266
267  /*
268  * now decide which vertices should be included
269  * local extrema are excluded
270  * local non-extrema are included
271  * vertices of horizontal edges which are pseudo-extrema
272  * are excluded.
273  * one vertex of horizontal edges which are pseudo-non-extrema
274  * are included.
275  */
276
277  prev = vertex_count - 1;
278  i = 0;
279  do {
280  next = i + 1;
281  if (next >= vertex_count)
282  next = 0;
283
284  if (extrema(vertex_points[prev].y, vertex_points[i].y,
285  vertex_points[next].y))
286  skip = 1;
287  else if (non_extrema(vertex_points[prev].y, vertex_points[i].y,
288  vertex_points[next].y))
289  skip = 0;
290  else {
291  skip = 0;
292  if (++next >= vertex_count)
293  next = 0;
294  if (extrema(vertex_points[prev].y, vertex_points[i].y,
295  vertex_points[next].y))
296  skip = 1;
297  }
298
299  if (!skip)
301  vertex_points[i].y);
302
303  i = next;
304  prev = i - 1;
305  } while (i != 0);
306
307  G_free(vertex_points);
308
309  /* sort the edge points by row and then by col */
310  qsort(perimeter->points, (size_t)perimeter->npoints, sizeof(IClass_point),
311  edge_order);
312
313  return 1;
314 }
315
316 /*!
317  \brief Converts edge to cells.
318
319  It rasterizes edge given by two vertices.
320  Resterized points are added to perimeter.
321
322  \param perimeter perimeter
323  \param x0,y0 first edge point row and cell
324  \param x1,y1 second edge point row and cell
325
326  \return 1 on success
327  \return 0 on error
328  */
329 int edge2perimeter(IClass_perimeter *perimeter, int x0, int y0, int x1, int y1)
330 {
331  float m;
332
333  float x;
334
335  if (y0 == y1)
336  return 0;
337
338  x = x0;
339  m = (float)(x0 - x1) / (float)(y0 - y1);
340
341  if (y0 < y1) {
342  while (++y0 < y1) {
343  x0 = (x += m) + .5;
345  }
346  }
347  else {
348  while (--y0 > y1) {
349  x0 = (x -= m) + .5;
351  }
352  }
353
354  return 1;
355 }
356
357 /*!
358  \brief Adds point to perimeter.
359
360  \a perimeter has to have allocated space for \c points member.
361
362  \param perimeter perimeter
363  \param x,y point row and cell
364  */
365 void perimeter_add_point(IClass_perimeter *perimeter, int x, int y)
366 {
367  int n;
368
369  G_debug(5, "perimeter_add_point(): x: %d, y: %d", x, y);
370
371  n = perimeter->npoints++;
372  perimeter->points[n].x = x;
373  perimeter->points[n].y = y;
374 }
375
376 /*!
377  \brief Determines points order during sorting.
378
379  \param aa first IClass_point
380  \param bb second IClass_point
381  */
382 int edge_order(const void *aa, const void *bb)
383 {
384  const IClass_point *a = aa;
385
386  const IClass_point *b = bb;
387
388  if (a->y < b->y)
389  return -1;
390  if (a->y > b->y)
391  return 1;
392
393  if (a->x < b->x)
394  return -1;
395  if (a->x > b->x)
396  return 1;
397
398  return 0;
399 }
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
#define G_calloc(m, n)
Definition: defs/gis.h:95
void G_warning(const char *,...) __attribute__((format(printf
int G_debug(int, const char *,...) __attribute__((format(printf
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
plus_t Vect_get_num_areas(struct Map_info *)
Get number of areas in vector map.
Definition: level_two.c:87
int Vect_area_alive(struct Map_info *, int)
Check if area is alive or dead (topological level required)
int Vect_get_field_number(struct Map_info *, const char *)
Get field number of given field.
Definition: field.c:608
int Vect_get_area_points(struct Map_info *, int, struct line_pnts *)
Returns polygon array of points (outer ring) of given area.
int Vect_get_area_cat(struct Map_info *, int, int)
Find FIRST category of given field and area.
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
#define _(str)
Definition: glocale.h:10
void perimeter_add_point(IClass_perimeter *perimeter, int x, int y)
int vector2perimeters(struct Map_info *Map, const char *layer_name, int category, IClass_perimeter_list *perimeters, struct Cell_head *band_region)
Creates perimeters from vector areas of given category.
int edge2perimeter(IClass_perimeter *perimeter, int x0, int y0, int x1, int y1)
Converts edge to cells.
void free_perimeters(IClass_perimeter_list *perimeters)
Frees all perimeters in list of perimeters.
int make_perimeter(struct line_pnts *points, IClass_perimeter *perimeter, struct Cell_head *band_region)
Creates one perimeter from vector area.
#define extrema(x, y, z)
#define non_extrema(x, y, z)
int edge_order(const void *aa, const void *bb)
Determines points order during sorting.
int count
double b
Definition: r_raster.c:39
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 ns_res
Resolution - north to south cell size for 2D data.
Definition: gis.h:477
double west
Extent coordinates (west)
Definition: gis.h:489
Vector map info.
Definition: dig_structs.h:1243
Feature geometry info - coordinates.
Definition: dig_structs.h:1651
double * y
Array of Y coordinates.
Definition: dig_structs.h:1659
double * x
Array of X coordinates.
Definition: dig_structs.h:1655
int n_points
Number of points.
Definition: dig_structs.h:1667
#define x