GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-2653d7b383
diglib/poly.c
Go to the documentation of this file.
1 /*****************************************************************************
2  *
3  * MODULE: Vector library
4  *
5  * AUTHOR(S): Original author CERL, probably Dave Gerdes.
6  * Update to GRASS 5.7 Radim Blazek.
7  * Update to GRASS 7.0 Markus Metz
8  *
9  * PURPOSE: Lower level functions for reading/writing/manipulating vectors.
10  *
11  * COPYRIGHT: (C) 2009 by the GRASS Development Team
12  *
13  * This program is free software under the GNU General Public
14  * License (>=v2). Read the file COPYING that comes with GRASS
15  * for details.
16  *
17  *****************************************************************************/
18 
19 #include <math.h>
20 #include <grass/vector.h>
21 
22 #ifndef HUGE_VAL
23 #define HUGE_VAL 9999999999999.0
24 #endif
25 
26 /*
27  * fills BPoints (must be inited previously) by points from input
28  * array LPoints
29  *
30  * each input LPoints[i] must have at least 2 points
31  *
32  * returns number of points or -1 on error
33  */
34 
35 int dig_get_poly_points(int n_lines, struct line_pnts **LPoints,
36  int *direction, /* line direction: > 0 or < 0 */
37  struct line_pnts *BPoints)
38 {
39  register int i, j, point, start, end, inc;
40  struct line_pnts *Points;
41  int n_points;
42 
43  BPoints->n_points = 0;
44 
45  if (n_lines < 1) {
46  return 0;
47  }
48 
49  /* Calc required space */
50  n_points = 0;
51  for (i = 0; i < n_lines; i++) {
52  Points = LPoints[i];
53  n_points += Points->n_points - 1; /* each line from first to last - 1 */
54  }
55  n_points++; /* last point */
56 
57  if (0 > dig_alloc_points(BPoints, n_points))
58  return (-1);
59 
60  point = 0;
61  j = 0;
62  for (i = 0; i < n_lines; i++) {
63  Points = LPoints[i];
64  if (direction[i] > 0) {
65  start = 0;
66  end = Points->n_points - 1;
67  inc = 1;
68  }
69  else {
70  start = Points->n_points - 1;
71  end = 0;
72  inc = -1;
73  }
74 
75  for (j = start; j != end; j += inc) {
76  BPoints->x[point] = Points->x[j];
77  BPoints->y[point] = Points->y[j];
78  point++;
79  }
80  }
81  /* last point */
82  BPoints->x[point] = Points->x[j];
83  BPoints->y[point] = Points->y[j];
84 
85  BPoints->n_points = n_points;
86 
87  return (BPoints->n_points);
88 }
89 
90 /*
91  * calculate signed area size for polygon
92  *
93  * points must be closed polygon with first point = last point
94  *
95  * returns signed area, positive for clockwise, negative for
96  * counterclockwise, 0 for degenerate
97  */
98 int dig_find_area_poly(struct line_pnts *Points, double *totalarea)
99 {
100  int i;
101  double *x, *y;
102  double tot_area;
103 
104  x = Points->x;
105  y = Points->y;
106 
107  /* line integral: *Points do not need to be pruned */
108  /* surveyor's formula is more common, but more prone to
109  * fp precision limit errors, and *Points would need to be pruned */
110  tot_area = 0.0;
111  for (i = 1; i < Points->n_points; i++) {
112  tot_area += (x[i] - x[i - 1]) * (y[i] + y[i - 1]);
113  }
114  *totalarea = 0.5 * tot_area;
115 
116  return (0);
117 }
118 
119 /*
120  * find orientation of polygon (clockwise or counterclockwise)
121  * in theory faster than signed area for > 4 vertices, but is not robust
122  * against special cases
123  * use dig_find_area_poly instead
124  *
125  * points must be closed polygon with first point = last point
126  *
127  * this code uses bits and pieces from softSurfer and GEOS
128  * (C) 2000 softSurfer (www.softsurfer.com)
129  * (C) 2006 Refractions Research Inc.
130  *
131  * copes with partially collapsed boundaries and 8-shaped isles
132  * the code is long and not much faster than dig_find_area_poly
133  * it can be written much shorter, but that comes with speed penalty
134  *
135  * returns orientation, positive for CW, negative for CCW, 0 for degenerate
136  */
137 double dig_find_poly_orientation(struct line_pnts *Points)
138 {
139  unsigned int pnext, pprev, pcur = 0;
140  unsigned int lastpoint = Points->n_points - 1;
141  double *x, *y, orientation;
142 
143  x = Points->x;
144  y = Points->y;
145 
146  /* first find leftmost highest vertex of the polygon */
147  for (pnext = 1; pnext < lastpoint; pnext++) {
148  if (y[pnext] < y[pcur])
149  continue;
150  else if (y[pnext] == y[pcur]) { /* just as high */
151  if (x[pnext] > x[pcur]) /* but to the right */
152  continue;
153  if (x[pnext] ==
154  x[pcur]) { /* duplicate point, self-intersecting polygon ? */
155  pprev = (pcur == 0 ? lastpoint - 1 : pcur - 1);
156  if (y[pnext - 1] < y[pprev])
157  continue;
158  }
159  }
160  pcur = pnext; /* a new leftmost highest vertex */
161  }
162 
163  /* Points are not pruned, so ... */
164  pnext = pcur;
165  pprev = pcur;
166 
167  /* find next distinct point */
168  do {
169  if (pnext < lastpoint - 1)
170  pnext++;
171  else
172  pnext = 0;
173  } while (pnext != pcur && x[pcur] == x[pnext] && y[pcur] == y[pnext]);
174 
175  /* find previous distinct point */
176  do {
177  if (pprev > 0)
178  pprev--;
179  else
180  pprev = lastpoint - 1;
181  } while (pprev != pcur && x[pcur] == x[pprev] && y[pcur] == y[pprev]);
182 
183  /* orientation at vertex pcur == signed area for triangle pprev, pcur, pnext
184  * rather use robust determinant of Olivier Devillers? */
185  orientation = (x[pnext] - x[pprev]) * (y[pcur] - y[pprev]) -
186  (x[pcur] - x[pprev]) * (y[pnext] - y[pprev]);
187 
188  if (orientation)
189  return orientation;
190 
191  /* orientation is 0, can happen with dirty boundaries, next check */
192  /* find rightmost highest vertex of the polygon */
193  pcur = 0;
194  for (pnext = 1; pnext < lastpoint; pnext++) {
195  if (y[pnext] < y[pcur])
196  continue;
197  else if (y[pnext] == y[pcur]) { /* just as high */
198  if (x[pnext] < x[pcur]) /* but to the left */
199  continue;
200  if (x[pnext] ==
201  x[pcur]) { /* duplicate point, self-intersecting polygon ? */
202  pprev = (pcur == 0 ? lastpoint - 1 : pcur - 1);
203  if (y[pnext - 1] < y[pprev])
204  continue;
205  }
206  }
207  pcur = pnext; /* a new rightmost highest vertex */
208  }
209 
210  /* Points are not pruned, so ... */
211  pnext = pcur;
212  pprev = pcur;
213 
214  /* find next distinct point */
215  do {
216  if (pnext < lastpoint - 1)
217  pnext++;
218  else
219  pnext = 0;
220  } while (pnext != pcur && x[pcur] == x[pnext] && y[pcur] == y[pnext]);
221 
222  /* find previous distinct point */
223  do {
224  if (pprev > 0)
225  pprev--;
226  else
227  pprev = lastpoint - 1;
228  } while (pprev != pcur && x[pcur] == x[pprev] && y[pcur] == y[pprev]);
229 
230  /* orientation at vertex pcur == signed area for triangle pprev, pcur, pnext
231  * rather use robust determinant of Olivier Devillers? */
232  orientation = (x[pnext] - x[pprev]) * (y[pcur] - y[pprev]) -
233  (x[pcur] - x[pprev]) * (y[pnext] - y[pprev]);
234 
235  if (orientation)
236  return orientation;
237 
238  /* orientation is 0, next check */
239  /* find leftmost lowest vertex of the polygon */
240  pcur = 0;
241  for (pnext = 1; pnext < lastpoint; pnext++) {
242  if (y[pnext] > y[pcur])
243  continue;
244  else if (y[pnext] == y[pcur]) { /* just as low */
245  if (x[pnext] > x[pcur]) /* but to the right */
246  continue;
247  if (x[pnext] ==
248  x[pcur]) { /* duplicate point, self-intersecting polygon ? */
249  pprev = (pcur == 0 ? lastpoint - 1 : pcur - 1);
250  if (y[pnext - 1] > y[pprev])
251  continue;
252  }
253  }
254  pcur = pnext; /* a new leftmost lowest vertex */
255  }
256 
257  /* Points are not pruned, so ... */
258  pnext = pcur;
259  pprev = pcur;
260 
261  /* find next distinct point */
262  do {
263  if (pnext < lastpoint - 1)
264  pnext++;
265  else
266  pnext = 0;
267  } while (pnext != pcur && x[pcur] == x[pnext] && y[pcur] == y[pnext]);
268 
269  /* find previous distinct point */
270  do {
271  if (pprev > 0)
272  pprev--;
273  else
274  pprev = lastpoint - 1;
275  } while (pprev != pcur && x[pcur] == x[pprev] && y[pcur] == y[pprev]);
276 
277  /* orientation at vertex pcur == signed area for triangle pprev, pcur, pnext
278  * rather use robust determinant of Olivier Devillers? */
279  orientation = (x[pnext] - x[pprev]) * (y[pcur] - y[pprev]) -
280  (x[pcur] - x[pprev]) * (y[pnext] - y[pprev]);
281 
282  if (orientation)
283  return orientation;
284 
285  /* orientation is 0, last check */
286  /* find rightmost lowest vertex of the polygon */
287  pcur = 0;
288  for (pnext = 1; pnext < lastpoint; pnext++) {
289  if (y[pnext] > y[pcur])
290  continue;
291  else if (y[pnext] == y[pcur]) { /* just as low */
292  if (x[pnext] < x[pcur]) /* but to the left */
293  continue;
294  if (x[pnext] ==
295  x[pcur]) { /* duplicate point, self-intersecting polygon ? */
296  pprev = (pcur == 0 ? lastpoint - 1 : pcur - 1);
297  if (y[pnext - 1] > y[pprev])
298  continue;
299  }
300  }
301  pcur = pnext; /* a new rightmost lowest vertex */
302  }
303 
304  /* Points are not pruned, so ... */
305  pnext = pcur;
306  pprev = pcur;
307 
308  /* find next distinct point */
309  do {
310  if (pnext < lastpoint - 1)
311  pnext++;
312  else
313  pnext = 0;
314  } while (pnext != pcur && x[pcur] == x[pnext] && y[pcur] == y[pnext]);
315 
316  /* find previous distinct point */
317  do {
318  if (pprev > 0)
319  pprev--;
320  else
321  pprev = lastpoint - 1;
322  } while (pprev != pcur && x[pcur] == x[pprev] && y[pcur] == y[pprev]);
323 
324  /* orientation at vertex pcur == signed area for triangle pprev, pcur, pnext
325  * rather use robust determinant of Olivier Devillers? */
326  orientation = (x[pnext] - x[pprev]) * (y[pcur] - y[pprev]) -
327  (x[pcur] - x[pprev]) * (y[pnext] - y[pprev]);
328 
329  return orientation; /* 0 for degenerate */
330 }
int dig_alloc_points(struct line_pnts *, int)
allocate room for 'num' X and Y arrays in struct line_pnts
Definition: struct_alloc.c:336
int dig_get_poly_points(int n_lines, struct line_pnts **LPoints, int *direction, struct line_pnts *BPoints)
Definition: diglib/poly.c:35
int dig_find_area_poly(struct line_pnts *Points, double *totalarea)
Definition: diglib/poly.c:98
double dig_find_poly_orientation(struct line_pnts *Points)
Definition: diglib/poly.c:137
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