GRASS GIS 7 Programmer's Manual  7.5.svn(2018)-r72846
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
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,
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,
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)
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;
351  }
352  }
353  else {
354  while (--y0 > y1) {
355  x0 = (x -= m) + .5;
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 Vect_get_area_cat(const struct Map_info *Map, int area, int field)
Find FIRST category of given field and area.
plus_t Vect_get_num_areas(const struct Map_info *Map)
Get number of areas in vector map.
Definition: level_two.c:86
CELL cat
Definition: raster3d/cats.c:82
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:149
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:390
struct line_pnts * Vect_new_line_struct()
Creates and initializes a line_pnts structure.
Definition: line.c:45
double west
Extent coordinates (west)
Definition: gis.h:442
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
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.
void Vect_destroy_line_struct(struct line_pnts *p)
Frees all memory associated with a line_pnts structure, including the structure itself.
Definition: line.c:77
#define x
double * x
Array of X coordinates.
Definition: dig_structs.h:1680
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:436
double b
Definition: r_raster.c:39
#define non_extrema(x, y, z)
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
void perimeter_add_point(IClass_perimeter *perimeter, int x, int y)
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:430
#define extrema(x, y, z)
#define _(str)
Definition: glocale.h:13
int Vect_get_area_points(const struct Map_info *Map, int area, struct line_pnts *BPoints)
Returns polygon array of points (outer ring) of given area.
double ew_res
Resolution - east to west cell size for 2D data.
Definition: gis.h:426
int Vect_get_field_number(const struct Map_info *Map, const char *field)
Get field number of given field.
Definition: field.c:553
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204
int Vect_area_alive(const struct Map_info *Map, int area)
Check if area is alive or dead (topological level required)