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)
20  \author Brad Douglas <rez touchofmadness.com>(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,
53  struct Cell_head *band_region)
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,
157  struct Cell_head *band_region)
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)
300  perimeter_add_point(perimeter, vertex_points[i].x,
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;
344  perimeter_add_point(perimeter, x0, y0);
345  }
346  }
347  else {
348  while (--y0 > y1) {
349  x0 = (x -= m) + .5;
350  perimeter_add_point(perimeter, x0, y0);
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)
Adds point to perimeter.
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