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