GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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 
37 #define extrema(x,y,z) (((x<y)&&(z<y))||((x>y)&&(z>y)))
38 #define non_extrema(x,y,z) (((x<y)&&(y<z))||((x>y)&&(y>z)))
39 
40 /*!
41  \brief Creates perimeters from vector areas of given category.
42 
43  \param Map vector map
44  \param layer_name layer name (within vector map)
45  \param category vector category (cat column value)
46  \param[out] perimeters list of perimeters
47  \param band_region region which determines perimeter cells
48 
49  \return number of areas of given cat
50  \return -1 on error
51  */
52 int vector2perimeters(struct Map_info *Map, const char *layer_name,
53  int category, IClass_perimeter_list * perimeters,
54  struct Cell_head *band_region)
55 {
56  struct line_pnts *points;
57 
58  int nareas, nareas_cat, layer;
59 
60  int i, cat, ret;
61 
62  int j;
63 
64  G_debug(3, "iclass_vector2perimeters():layer = %s, category = %d",
65  layer_name, category);
66 
67  layer = Vect_get_field_number(Map, layer_name);
68  nareas = Vect_get_num_areas(Map);
69  if (nareas == 0)
70  return 0;
71 
72  nareas_cat = 0;
73  /* find out, how many areas have given category */
74  for (i = 1; i <= nareas; i++) {
75  if (!Vect_area_alive(Map, i))
76  continue;
77  cat = Vect_get_area_cat(Map, i, layer);
78  if (cat < 0) {
79  /* no centroid, no category */
80  }
81  else if (cat == category) {
82  nareas_cat++;
83  }
84  }
85  if (nareas_cat == 0)
86  return 0;
87 
88  perimeters->nperimeters = nareas_cat;
89  perimeters->perimeters =
90  (IClass_perimeter *) G_calloc(nareas_cat, sizeof(IClass_perimeter));
91 
92  j = 0; /* area with cat */
93  for (i = 1; i <= nareas; i++) {
94  if (!Vect_area_alive(Map, i))
95  continue;
96  cat = Vect_get_area_cat(Map, i, layer);
97  if (cat < 0) {
98  /* no centroid, no category */
99  }
100  else if (cat == category) {
101  j++;
102 
103  points = Vect_new_line_struct(); /* Vect_destroy_line_struct */
104  ret = Vect_get_area_points(Map, i, points);
105 
106  if (ret <= 0) {
107  Vect_destroy_line_struct(points);
108  free_perimeters(perimeters);
109  G_warning(_("Get area %d failed"), i);
110  return -1;
111  }
112  if (make_perimeter
113  (points, &perimeters->perimeters[j - 1], band_region) <= 0) {
114  Vect_destroy_line_struct(points);
115  free_perimeters(perimeters);
116  G_warning(_("Perimeter computation failed"));
117  return -1;
118  }
119  Vect_destroy_line_struct(points);
120  }
121 
122  }
123 
124  /* Vect_close(&Map); */
125 
126  return nareas_cat;
127 }
128 
129 /*!
130  \brief Frees all perimeters in list of perimeters.
131 
132  It also frees list of perimeters itself.
133 
134  \param perimeters list of perimeters
135  */
136 void free_perimeters(IClass_perimeter_list * perimeters)
137 {
138  int i;
139 
140  G_debug(5, "free_perimeters()");
141 
142  for (i = 0; i < perimeters->nperimeters; i++) {
143  G_free(perimeters->perimeters[i].points);
144  }
145  G_free(perimeters->perimeters);
146 }
147 
148 /*!
149  \brief Creates one perimeter from vector area.
150 
151  \param points list of vertices represting area
152  \param[out] perimeter perimeter
153  \param band_region region which determines perimeter cells
154 
155  \return 1 on success
156  \return 0 on error
157  */
158 int make_perimeter(struct line_pnts *points, IClass_perimeter * perimeter,
159  struct Cell_head *band_region)
160 {
161  IClass_point *tmp_points;
162 
163  IClass_point *vertex_points;
164 
165  int i, first, prev, skip, next;
166 
167  int count, vertex_count;
168 
169  int np; /* perimeter estimate */
170 
171  G_debug(5, "iclass_make_perimeter()");
172  count = points->n_points;
173 
174  tmp_points = (IClass_point *) G_calloc(count, sizeof(IClass_point)); /* TODO test */
175 
176  for (i = 0; i < count; i++) {
177  G_debug(5, "iclass_make_perimeter(): points: x: %f y: %f",
178  points->x[i], points->y[i]);
179 
180  /* This functions are no longer used because of the different behavior
181  of Rast_easting_to_col depending whether location is LL or not.
182  It makes problem in interactive scatter plot tool,
183  which defines its own coordinates systems for the plots and
184  therefore it requires the function to work always in same way
185  without hidden dependency on location type.
186 
187  tmp_points[i].y = Rast_northing_to_row(points->y[i], band_region);
188  tmp_points[i].x = Rast_easting_to_col(points->x[i], band_region);
189  */
190 
191  tmp_points[i].y = (band_region->north - points->y[i]) / band_region->ns_res;
192  tmp_points[i].x = (points->x[i] - band_region->west) / band_region->ew_res;
193  }
194 
195  /* find first edge which is not horizontal */
196 
197  first = -1;
198  prev = count - 1;
199  for (i = 0; i < count; prev = i++) {
200  /* non absurd polygon has vertexes with different y coordinate */
201  if (tmp_points[i].y != tmp_points[prev].y) {
202  first = i;
203  break;
204  }
205  }
206  if (first < 0) {
207  G_free(tmp_points);
208  G_warning(_("Invalid polygon"));
209  return 0;
210  }
211 
212  /* copy tmp to vertex list collapsing adjacent horizontal edges */
213 
214  /* vertex_count <= count, size of vertex_points is count */
215  vertex_points = (IClass_point *) G_calloc(count, sizeof(IClass_point)); /* TODO test */
216  skip = 0;
217  vertex_count = 0;
218  i = first; /* stmt not necssary */
219 
220  do {
221  if (!skip) {
222  vertex_points[vertex_count].x = tmp_points[i].x;
223  vertex_points[vertex_count].y = tmp_points[i].y;
224  vertex_count++;
225  }
226 
227  prev = i++;
228  if (i >= count)
229  i = 0;
230  if ((next = i + 1) >= count)
231  next = 0;
232 
233  skip = ((tmp_points[prev].y == tmp_points[i].y) &&
234  (tmp_points[next].y == tmp_points[i].y));
235  }
236  while (i != first);
237 
238  G_free(tmp_points);
239 
240  /* count points on the perimeter */
241 
242  np = 0;
243  prev = vertex_count - 1;
244  for (i = 0; i < vertex_count; prev = i++) {
245  np += abs(vertex_points[prev].y - vertex_points[i].y);
246  }
247 
248  /* allocate perimeter list */
249 
250  perimeter->points = (IClass_point *) G_calloc(np, sizeof(IClass_point));
251  if (!perimeter->points) {
252  G_free(vertex_points);
253  G_warning(_("Outlined area is too large."));
254  return 0;
255  }
256 
257  /* store the perimeter points */
258 
259  perimeter->npoints = 0;
260  prev = vertex_count - 1;
261  for (i = 0; i < vertex_count; prev = i++) {
262  edge2perimeter(perimeter, vertex_points[prev].x,
263  vertex_points[prev].y, vertex_points[i].x,
264  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
285  (vertex_points[prev].y, vertex_points[i].y,
286  vertex_points[next].y))
287  skip = 1;
288  else if (non_extrema
289  (vertex_points[prev].y, vertex_points[i].y,
290  vertex_points[next].y))
291  skip = 0;
292  else {
293  skip = 0;
294  if (++next >= vertex_count)
295  next = 0;
296  if (extrema
297  (vertex_points[prev].y, vertex_points[i].y,
298  vertex_points[next].y))
299  skip = 1;
300  }
301 
302  if (!skip)
303  perimeter_add_point(perimeter, vertex_points[i].x,
304  vertex_points[i].y);
305 
306  i = next;
307  prev = i - 1;
308  }
309  while (i != 0);
310 
311  G_free(vertex_points);
312 
313  /* sort the edge points by row and then by col */
314  qsort(perimeter->points, (size_t) perimeter->npoints,
315  sizeof(IClass_point), edge_order);
316 
317  return 1;
318 
319 }
320 
321 /*!
322  \brief Converts edge to cells.
323 
324  It rasterizes edge given by two vertices.
325  Resterized points are added to perimeter.
326 
327  \param perimeter perimeter
328  \param x0,y0 first edge point row and cell
329  \param x1,y1 second edge point row and cell
330 
331  \return 1 on success
332  \return 0 on error
333  */
334 int edge2perimeter(IClass_perimeter * perimeter, int x0, int y0, int x1,
335  int y1)
336 {
337  float m;
338 
339  float x;
340 
341  if (y0 == y1)
342  return 0;
343 
344  x = x0;
345  m = (float)(x0 - x1) / (float)(y0 - y1);
346 
347  if (y0 < y1) {
348  while (++y0 < y1) {
349  x0 = (x += m) + .5;
350  perimeter_add_point(perimeter, x0, y0);
351  }
352  }
353  else {
354  while (--y0 > y1) {
355  x0 = (x -= m) + .5;
356  perimeter_add_point(perimeter, x0, y0);
357  }
358  }
359 
360  return 1;
361 }
362 
363 /*!
364  \brief Adds point to perimeter.
365 
366  \a perimeter has to have allocated space for \c points member.
367 
368  \param perimeter perimeter
369  \param x,y point row and cell
370  */
371 void perimeter_add_point(IClass_perimeter * perimeter, int x, int y)
372 {
373  int n;
374 
375  G_debug(5, "perimeter_add_point(): x: %d, y: %d", x, y);
376 
377  n = perimeter->npoints++;
378  perimeter->points[n].x = x;
379  perimeter->points[n].y = y;
380 }
381 
382 /*!
383  \brief Determines points order during sorting.
384 
385  \param aa first IClass_point
386  \param bb second IClass_point
387  */
388 int edge_order(const void *aa, const void *bb)
389 {
390  const IClass_point *a = aa;
391 
392  const IClass_point *b = bb;
393 
394  if (a->y < b->y)
395  return -1;
396  if (a->y > b->y)
397  return 1;
398 
399  if (a->x < b->x)
400  return -1;
401  if (a->x > b->x)
402  return 1;
403 
404  return 0;
405 }
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
double west
Extent coordinates (west)
Definition: gis.h:475
int edge_order(const void *aa, const void *bb)
Determines points order during sorting.
int n_points
Number of points.
Definition: dig_structs.h:1692
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:149
int count
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.
#define x
#define G_calloc(m, n)
Definition: defs/gis.h:113
double * x
Array of X coordinates.
Definition: dig_structs.h:1680
int Vect_get_field_number(const struct Map_info *, const char *)
Get field number of given field.
Definition: field.c:599
Feature geometry info - coordinates.
Definition: dig_structs.h:1675
void free_perimeters(IClass_perimeter_list *perimeters)
Frees all perimeters in list of perimeters.
double north
Extent coordinates (north)
Definition: gis.h:469
double b
Definition: r_raster.c:39
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
#define non_extrema(x, y, z)
void perimeter_add_point(IClass_perimeter *perimeter, int x, int y)
Adds point to perimeter.
plus_t Vect_get_num_areas(const struct Map_info *)
Get number of areas in vector map.
Definition: level_two.c:86
int Vect_area_alive(const struct Map_info *, int)
Check if area is alive or dead (topological level required)
int edge2perimeter(IClass_perimeter *perimeter, int x0, int y0, int x1, int y1)
Converts edge to cells.
Vector map info.
Definition: dig_structs.h:1259
double * y
Array of Y coordinates.
Definition: dig_structs.h:1684
double ns_res
Resolution - north to south cell size for 2D data.
Definition: gis.h:463
#define extrema(x, y, z)
void G_warning(const char *,...) __attribute__((format(printf
#define _(str)
Definition: glocale.h:10
double ew_res
Resolution - east to west cell size for 2D data.
Definition: gis.h:459
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 Vect_get_area_points(const struct Map_info *, int, struct line_pnts *)
Returns polygon array of points (outer ring) of given area.
int G_debug(int, const char *,...) __attribute__((format(printf
int Vect_get_area_cat(const struct Map_info *, int, int)
Find FIRST category of given field and area.