1 | /********************************************************************** |
---|
2 | * $Id: measures.c 2797 2008-05-31 09:56:44Z mcayland $ |
---|
3 | * |
---|
4 | * PostGIS - Spatial Types for PostgreSQL |
---|
5 | * http://postgis.refractions.net |
---|
6 | * Copyright 2001-2006 Refractions Research Inc. |
---|
7 | * |
---|
8 | * This is free software; you can redistribute and/or modify it under |
---|
9 | * the terms of the GNU General Public Licence. See the COPYING file. |
---|
10 | * |
---|
11 | **********************************************************************/ |
---|
12 | |
---|
13 | #include <math.h> |
---|
14 | #include <string.h> |
---|
15 | |
---|
16 | #include "liblwgeom.h" |
---|
17 | |
---|
18 | |
---|
19 | /* |
---|
20 | * pt_in_ring_2d(): crossing number test for a point in a polygon |
---|
21 | * input: p = a point, |
---|
22 | * pa = vertex points of a ring V[n+1] with V[n]=V[0] |
---|
23 | * returns: 0 = outside, 1 = inside |
---|
24 | * |
---|
25 | * Our polygons have first and last point the same, |
---|
26 | * |
---|
27 | */ |
---|
28 | int pt_in_ring_2d(POINT2D *p, POINTARRAY *ring) |
---|
29 | { |
---|
30 | int cn = 0; /* the crossing number counter */ |
---|
31 | int i; |
---|
32 | POINT2D v1, v2; |
---|
33 | |
---|
34 | #if INTEGRITY_CHECKS |
---|
35 | POINT2D first, last; |
---|
36 | |
---|
37 | getPoint2d_p(ring, 0, &first); |
---|
38 | getPoint2d_p(ring, ring->npoints-1, &last); |
---|
39 | if ( memcmp(&first, &last, sizeof(POINT2D)) ) |
---|
40 | { |
---|
41 | lwerror("pt_in_ring_2d: V[n] != V[0] (%g %g != %g %g)", |
---|
42 | first.x, first.y, last.x, last.y); |
---|
43 | |
---|
44 | } |
---|
45 | #endif |
---|
46 | |
---|
47 | LWDEBUGF(2, "pt_in_ring_2d called with point: %g %g", p->x, p->y); |
---|
48 | /* printPA(ring); */ |
---|
49 | |
---|
50 | /* loop through all edges of the polygon */ |
---|
51 | getPoint2d_p(ring, 0, &v1); |
---|
52 | for (i=0; i<ring->npoints-1; i++) |
---|
53 | { |
---|
54 | double vt; |
---|
55 | getPoint2d_p(ring, i+1, &v2); |
---|
56 | |
---|
57 | /* edge from vertex i to vertex i+1 */ |
---|
58 | if |
---|
59 | ( |
---|
60 | /* an upward crossing */ |
---|
61 | ((v1.y <= p->y) && (v2.y > p->y)) |
---|
62 | /* a downward crossing */ |
---|
63 | || ((v1.y > p->y) && (v2.y <= p->y)) |
---|
64 | ) |
---|
65 | { |
---|
66 | |
---|
67 | vt = (double)(p->y - v1.y) / (v2.y - v1.y); |
---|
68 | |
---|
69 | /* P.x <intersect */ |
---|
70 | if (p->x < v1.x + vt * (v2.x - v1.x)) |
---|
71 | { |
---|
72 | /* a valid crossing of y=p.y right of p.x */ |
---|
73 | ++cn; |
---|
74 | } |
---|
75 | } |
---|
76 | v1 = v2; |
---|
77 | } |
---|
78 | |
---|
79 | LWDEBUGF(3, "pt_in_ring_2d returning %d", cn&1); |
---|
80 | |
---|
81 | return (cn&1); /* 0 if even (out), and 1 if odd (in) */ |
---|
82 | } |
---|
83 | |
---|
84 | double distance2d_pt_pt(POINT2D *p1, POINT2D *p2) |
---|
85 | { |
---|
86 | double hside = p2->x - p1->x; |
---|
87 | double vside = p2->y - p1->y; |
---|
88 | |
---|
89 | return sqrt ( hside*hside + vside*vside ); |
---|
90 | |
---|
91 | /* the above is more readable |
---|
92 | return sqrt( |
---|
93 | (p2->x-p1->x) * (p2->x-p1->x) + (p2->y-p1->y) * (p2->y-p1->y) |
---|
94 | ); */ |
---|
95 | } |
---|
96 | |
---|
97 | /*distance2d from p to line A->B */ |
---|
98 | double distance2d_pt_seg(POINT2D *p, POINT2D *A, POINT2D *B) |
---|
99 | { |
---|
100 | double r,s; |
---|
101 | |
---|
102 | /*if start==end, then use pt distance */ |
---|
103 | if ( ( A->x == B->x) && (A->y == B->y) ) |
---|
104 | return distance2d_pt_pt(p,A); |
---|
105 | |
---|
106 | /* |
---|
107 | * otherwise, we use comp.graphics.algorithms |
---|
108 | * Frequently Asked Questions method |
---|
109 | * |
---|
110 | * (1) AC dot AB |
---|
111 | * r = --------- |
---|
112 | * ||AB||^2 |
---|
113 | * r has the following meaning: |
---|
114 | * r=0 P = A |
---|
115 | * r=1 P = B |
---|
116 | * r<0 P is on the backward extension of AB |
---|
117 | * r>1 P is on the forward extension of AB |
---|
118 | * 0<r<1 P is interior to AB |
---|
119 | */ |
---|
120 | |
---|
121 | r = ( (p->x-A->x) * (B->x-A->x) + (p->y-A->y) * (B->y-A->y) )/( (B->x-A->x)*(B->x-A->x) +(B->y-A->y)*(B->y-A->y) ); |
---|
122 | |
---|
123 | if (r<0) return distance2d_pt_pt(p,A); |
---|
124 | if (r>1) return distance2d_pt_pt(p,B); |
---|
125 | |
---|
126 | |
---|
127 | /* |
---|
128 | * (2) |
---|
129 | * (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay) |
---|
130 | * s = ----------------------------- |
---|
131 | * L^2 |
---|
132 | * |
---|
133 | * Then the distance from C to P = |s|*L. |
---|
134 | * |
---|
135 | */ |
---|
136 | |
---|
137 | s = ( (A->y-p->y)*(B->x-A->x)- (A->x-p->x)*(B->y-A->y) ) / |
---|
138 | ( (B->x-A->x)*(B->x-A->x) +(B->y-A->y)*(B->y-A->y) ); |
---|
139 | |
---|
140 | return LW_ABS(s) * sqrt( |
---|
141 | (B->x-A->x)*(B->x-A->x) + (B->y-A->y)*(B->y-A->y) |
---|
142 | ); |
---|
143 | } |
---|
144 | |
---|
145 | /* find the minimum 2d distance from AB to CD */ |
---|
146 | double distance2d_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D) |
---|
147 | { |
---|
148 | |
---|
149 | double s_top, s_bot,s; |
---|
150 | double r_top, r_bot,r; |
---|
151 | |
---|
152 | LWDEBUGF(2, "distance2d_seg_seg [%g,%g]->[%g,%g] by [%g,%g]->[%g,%g]", |
---|
153 | A->x,A->y,B->x,B->y, C->x,C->y, D->x, D->y); |
---|
154 | |
---|
155 | |
---|
156 | /*A and B are the same point */ |
---|
157 | if ( ( A->x == B->x) && (A->y == B->y) ) |
---|
158 | return distance2d_pt_seg(A,C,D); |
---|
159 | |
---|
160 | /*U and V are the same point */ |
---|
161 | |
---|
162 | if ( ( C->x == D->x) && (C->y == D->y) ) |
---|
163 | return distance2d_pt_seg(D,A,B); |
---|
164 | |
---|
165 | /* AB and CD are line segments */ |
---|
166 | /* from comp.graphics.algo |
---|
167 | |
---|
168 | Solving the above for r and s yields |
---|
169 | (Ay-Cy)(Dx-Cx)-(Ax-Cx)(Dy-Cy) |
---|
170 | r = ----------------------------- (eqn 1) |
---|
171 | (Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx) |
---|
172 | |
---|
173 | (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay) |
---|
174 | s = ----------------------------- (eqn 2) |
---|
175 | (Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx) |
---|
176 | Let P be the position vector of the intersection point, then |
---|
177 | P=A+r(B-A) or |
---|
178 | Px=Ax+r(Bx-Ax) |
---|
179 | Py=Ay+r(By-Ay) |
---|
180 | By examining the values of r & s, you can also determine some other limiting conditions: |
---|
181 | If 0<=r<=1 & 0<=s<=1, intersection exists |
---|
182 | r<0 or r>1 or s<0 or s>1 line segments do not intersect |
---|
183 | If the denominator in eqn 1 is zero, AB & CD are parallel |
---|
184 | If the numerator in eqn 1 is also zero, AB & CD are collinear. |
---|
185 | |
---|
186 | */ |
---|
187 | r_top = (A->y-C->y)*(D->x-C->x) - (A->x-C->x)*(D->y-C->y) ; |
---|
188 | r_bot = (B->x-A->x)*(D->y-C->y) - (B->y-A->y)*(D->x-C->x) ; |
---|
189 | |
---|
190 | s_top = (A->y-C->y)*(B->x-A->x) - (A->x-C->x)*(B->y-A->y); |
---|
191 | s_bot = (B->x-A->x)*(D->y-C->y) - (B->y-A->y)*(D->x-C->x); |
---|
192 | |
---|
193 | if ( (r_bot==0) || (s_bot == 0) ) |
---|
194 | { |
---|
195 | return ( |
---|
196 | LW_MIN(distance2d_pt_seg(A,C,D), |
---|
197 | LW_MIN(distance2d_pt_seg(B,C,D), |
---|
198 | LW_MIN(distance2d_pt_seg(C,A,B), |
---|
199 | distance2d_pt_seg(D,A,B)) |
---|
200 | ) |
---|
201 | ) |
---|
202 | ); |
---|
203 | } |
---|
204 | s = s_top/s_bot; |
---|
205 | r= r_top/r_bot; |
---|
206 | |
---|
207 | if ((r<0) || (r>1) || (s<0) || (s>1) ) |
---|
208 | { |
---|
209 | /*no intersection */ |
---|
210 | return ( |
---|
211 | LW_MIN(distance2d_pt_seg(A,C,D), |
---|
212 | LW_MIN(distance2d_pt_seg(B,C,D), |
---|
213 | LW_MIN(distance2d_pt_seg(C,A,B), |
---|
214 | distance2d_pt_seg(D,A,B)) |
---|
215 | ) |
---|
216 | ) |
---|
217 | ); |
---|
218 | |
---|
219 | } |
---|
220 | else |
---|
221 | return -0; /*intersection exists */ |
---|
222 | |
---|
223 | } |
---|
224 | |
---|
225 | /* |
---|
226 | * search all the segments of pointarray to see which one is closest to p1 |
---|
227 | * Returns minimum distance between point and pointarray |
---|
228 | */ |
---|
229 | double distance2d_pt_ptarray(POINT2D *p, POINTARRAY *pa) |
---|
230 | { |
---|
231 | double result = 0; |
---|
232 | int t; |
---|
233 | POINT2D start, end; |
---|
234 | |
---|
235 | getPoint2d_p(pa, 0, &start); |
---|
236 | |
---|
237 | for (t=1; t<pa->npoints; t++) |
---|
238 | { |
---|
239 | double dist; |
---|
240 | getPoint2d_p(pa, t, &end); |
---|
241 | dist = distance2d_pt_seg(p, &start, &end); |
---|
242 | if (t==1) result = dist; |
---|
243 | else result = LW_MIN(result, dist); |
---|
244 | |
---|
245 | if ( result == 0 ) return 0; |
---|
246 | |
---|
247 | start = end; |
---|
248 | } |
---|
249 | |
---|
250 | return result; |
---|
251 | } |
---|
252 | |
---|
253 | /* test each segment of l1 against each segment of l2. Return min */ |
---|
254 | double distance2d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2) |
---|
255 | { |
---|
256 | double result = 99999999999.9; |
---|
257 | char result_okay = 0; /*result is a valid min */ |
---|
258 | int t,u; |
---|
259 | POINT2D start, end; |
---|
260 | POINT2D start2, end2; |
---|
261 | |
---|
262 | LWDEBUGF(2, "distance2d_ptarray_ptarray called (points: %d-%d)", |
---|
263 | l1->npoints, l2->npoints); |
---|
264 | |
---|
265 | getPoint2d_p(l1, 0, &start); |
---|
266 | for (t=1; t<l1->npoints; t++) /*for each segment in L1 */ |
---|
267 | { |
---|
268 | getPoint2d_p(l1, t, &end); |
---|
269 | |
---|
270 | getPoint2d_p(l2, 0, &start2); |
---|
271 | for (u=1; u<l2->npoints; u++) /*for each segment in L2 */ |
---|
272 | { |
---|
273 | double dist; |
---|
274 | |
---|
275 | getPoint2d_p(l2, u, &end2); |
---|
276 | |
---|
277 | dist = distance2d_seg_seg(&start, &end, &start2, &end2); |
---|
278 | |
---|
279 | LWDEBUGF(4, "line_line; seg %i * seg %i, dist = %g\n",t,u,dist); |
---|
280 | |
---|
281 | if (result_okay) |
---|
282 | result = LW_MIN(result,dist); |
---|
283 | else |
---|
284 | { |
---|
285 | result_okay = 1; |
---|
286 | result = dist; |
---|
287 | } |
---|
288 | |
---|
289 | LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f", |
---|
290 | t, u, dist, result); |
---|
291 | |
---|
292 | if (result <= 0) return 0; /*intersection */ |
---|
293 | |
---|
294 | start2 = end2; |
---|
295 | } |
---|
296 | start = end; |
---|
297 | } |
---|
298 | |
---|
299 | return result; |
---|
300 | } |
---|
301 | |
---|
302 | /* true if point is in poly (and not in its holes) */ |
---|
303 | int pt_in_poly_2d(POINT2D *p, LWPOLY *poly) |
---|
304 | { |
---|
305 | int i; |
---|
306 | |
---|
307 | /* Not in outer ring */ |
---|
308 | if ( ! pt_in_ring_2d(p, poly->rings[0]) ) return 0; |
---|
309 | |
---|
310 | /* Check holes */ |
---|
311 | for (i=1; i<poly->nrings; i++) |
---|
312 | { |
---|
313 | /* Inside a hole */ |
---|
314 | if ( pt_in_ring_2d(p, poly->rings[i]) ) return 0; |
---|
315 | } |
---|
316 | |
---|
317 | return 1; /* In outer ring, not in holes */ |
---|
318 | } |
---|
319 | |
---|
320 | /* |
---|
321 | * Brute force. |
---|
322 | * Test line-ring distance against each ring. |
---|
323 | * If there's an intersection (distance==0) then return 0 (crosses boundary). |
---|
324 | * Otherwise, test to see if any point is inside outer rings of polygon, |
---|
325 | * but not in inner rings. |
---|
326 | * If so, return 0 (line inside polygon), |
---|
327 | * otherwise return min distance to a ring (could be outside |
---|
328 | * polygon or inside a hole) |
---|
329 | */ |
---|
330 | double distance2d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly) |
---|
331 | { |
---|
332 | POINT2D pt; |
---|
333 | int i; |
---|
334 | double mindist = 0; |
---|
335 | |
---|
336 | LWDEBUGF(2, "distance2d_ptarray_poly called (%d rings)", poly->nrings); |
---|
337 | |
---|
338 | for (i=0; i<poly->nrings; i++) |
---|
339 | { |
---|
340 | double dist = distance2d_ptarray_ptarray(pa, poly->rings[i]); |
---|
341 | if (i) mindist = LW_MIN(mindist, dist); |
---|
342 | else mindist = dist; |
---|
343 | |
---|
344 | LWDEBUGF(3, " distance from ring %d: %f, mindist: %f", |
---|
345 | i, dist, mindist); |
---|
346 | |
---|
347 | if ( mindist <= 0 ) return 0.0; /* intersection */ |
---|
348 | } |
---|
349 | |
---|
350 | /* |
---|
351 | * No intersection, have to check if a point is |
---|
352 | * inside polygon |
---|
353 | */ |
---|
354 | getPoint2d_p(pa, 0, &pt); |
---|
355 | |
---|
356 | /* |
---|
357 | * Outside outer ring, so min distance to a ring |
---|
358 | * is the actual min distance |
---|
359 | */ |
---|
360 | if ( ! pt_in_ring_2d(&pt, poly->rings[0]) ) return mindist; |
---|
361 | |
---|
362 | |
---|
363 | /* |
---|
364 | * Its in the outer ring. |
---|
365 | * Have to check if its inside a hole |
---|
366 | */ |
---|
367 | for (i=1; i<poly->nrings; i++) |
---|
368 | { |
---|
369 | if ( pt_in_ring_2d(&pt, poly->rings[i]) ) |
---|
370 | { |
---|
371 | /* |
---|
372 | * Its inside a hole, then the actual |
---|
373 | * distance is the min ring distance |
---|
374 | */ |
---|
375 | return mindist; |
---|
376 | } |
---|
377 | } |
---|
378 | |
---|
379 | return 0.0; /* Not in hole, so inside polygon */ |
---|
380 | } |
---|
381 | |
---|
382 | double distance2d_point_point(LWPOINT *point1, LWPOINT *point2) |
---|
383 | { |
---|
384 | POINT2D p1; |
---|
385 | POINT2D p2; |
---|
386 | |
---|
387 | getPoint2d_p(point1->point, 0, &p1); |
---|
388 | getPoint2d_p(point2->point, 0, &p2); |
---|
389 | |
---|
390 | return distance2d_pt_pt(&p1, &p2); |
---|
391 | } |
---|
392 | |
---|
393 | double distance2d_point_line(LWPOINT *point, LWLINE *line) |
---|
394 | { |
---|
395 | POINT2D p; |
---|
396 | POINTARRAY *pa = line->points; |
---|
397 | getPoint2d_p(point->point, 0, &p); |
---|
398 | return distance2d_pt_ptarray(&p, pa); |
---|
399 | } |
---|
400 | |
---|
401 | double distance2d_line_line(LWLINE *line1, LWLINE *line2) |
---|
402 | { |
---|
403 | POINTARRAY *pa1 = line1->points; |
---|
404 | POINTARRAY *pa2 = line2->points; |
---|
405 | return distance2d_ptarray_ptarray(pa1, pa2); |
---|
406 | } |
---|
407 | |
---|
408 | /* |
---|
409 | * 1. see if pt in outer boundary. if no, then treat the outer ring like a line |
---|
410 | * 2. if in the boundary, test to see if its in a hole. |
---|
411 | * if so, then return dist to hole, else return 0 (point in polygon) |
---|
412 | */ |
---|
413 | double distance2d_point_poly(LWPOINT *point, LWPOLY *poly) |
---|
414 | { |
---|
415 | POINT2D p; |
---|
416 | int i; |
---|
417 | |
---|
418 | getPoint2d_p(point->point, 0, &p); |
---|
419 | |
---|
420 | LWDEBUG(2, "distance2d_point_poly called"); |
---|
421 | |
---|
422 | /* Return distance to outer ring if not inside it */ |
---|
423 | if ( ! pt_in_ring_2d(&p, poly->rings[0]) ) |
---|
424 | { |
---|
425 | LWDEBUG(3, " not inside outer-ring"); |
---|
426 | |
---|
427 | return distance2d_pt_ptarray(&p, poly->rings[0]); |
---|
428 | } |
---|
429 | |
---|
430 | /* |
---|
431 | * Inside the outer ring. |
---|
432 | * Scan though each of the inner rings looking to |
---|
433 | * see if its inside. If not, distance==0. |
---|
434 | * Otherwise, distance = pt to ring distance |
---|
435 | */ |
---|
436 | for (i=1; i<poly->nrings; i++) |
---|
437 | { |
---|
438 | /* Inside a hole. Distance = pt -> ring */ |
---|
439 | if ( pt_in_ring_2d(&p, poly->rings[i]) ) |
---|
440 | { |
---|
441 | LWDEBUG(3, " inside an hole"); |
---|
442 | |
---|
443 | return distance2d_pt_ptarray(&p, poly->rings[i]); |
---|
444 | } |
---|
445 | } |
---|
446 | |
---|
447 | LWDEBUG(3, " inside the polygon"); |
---|
448 | |
---|
449 | return 0.0; /* Is inside the polygon */ |
---|
450 | } |
---|
451 | |
---|
452 | /* |
---|
453 | * Brute force. |
---|
454 | * Test to see if any rings intersect. |
---|
455 | * If yes, dist=0. |
---|
456 | * Test to see if one inside the other and if they are inside holes. |
---|
457 | * Find min distance ring-to-ring. |
---|
458 | */ |
---|
459 | double distance2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2) |
---|
460 | { |
---|
461 | POINT2D pt; |
---|
462 | double mindist = -1; |
---|
463 | int i; |
---|
464 | |
---|
465 | LWDEBUG(2, "distance2d_poly_poly called"); |
---|
466 | |
---|
467 | /* if poly1 inside poly2 return 0 */ |
---|
468 | getPoint2d_p(poly1->rings[0], 0, &pt); |
---|
469 | if ( pt_in_poly_2d(&pt, poly2) ) return 0.0; |
---|
470 | |
---|
471 | /* if poly2 inside poly1 return 0 */ |
---|
472 | getPoint2d_p(poly2->rings[0], 0, &pt); |
---|
473 | if ( pt_in_poly_2d(&pt, poly1) ) return 0.0; |
---|
474 | |
---|
475 | LWDEBUG(3, " polys not inside each other"); |
---|
476 | |
---|
477 | /* |
---|
478 | * foreach ring in Poly1 |
---|
479 | * foreach ring in Poly2 |
---|
480 | * if intersect, return 0 |
---|
481 | */ |
---|
482 | for (i=0; i<poly1->nrings; i++) |
---|
483 | { |
---|
484 | int j; |
---|
485 | for (j=0; j<poly2->nrings; j++) |
---|
486 | { |
---|
487 | double d = distance2d_ptarray_ptarray(poly1->rings[i], |
---|
488 | poly2->rings[j]); |
---|
489 | if ( d <= 0 ) return 0.0; |
---|
490 | |
---|
491 | /* mindist is -1 when not yet set */ |
---|
492 | if (mindist > -1) mindist = LW_MIN(mindist, d); |
---|
493 | else mindist = d; |
---|
494 | |
---|
495 | LWDEBUGF(3, " ring%i-%i dist: %f, mindist: %f", i, j, d, mindist); |
---|
496 | } |
---|
497 | |
---|
498 | } |
---|
499 | |
---|
500 | /* otherwise return closest approach of rings (no intersection) */ |
---|
501 | return mindist; |
---|
502 | |
---|
503 | } |
---|
504 | |
---|
505 | double distance2d_line_poly(LWLINE *line, LWPOLY *poly) |
---|
506 | { |
---|
507 | return distance2d_ptarray_poly(line->points, poly); |
---|
508 | } |
---|
509 | |
---|
510 | |
---|
511 | /*find the 2d length of the given POINTARRAY (even if it's 3d) */ |
---|
512 | double lwgeom_pointarray_length2d(POINTARRAY *pts) |
---|
513 | { |
---|
514 | double dist = 0.0; |
---|
515 | int i; |
---|
516 | POINT2D frm; |
---|
517 | POINT2D to; |
---|
518 | |
---|
519 | if ( pts->npoints < 2 ) return 0.0; |
---|
520 | for (i=0; i<pts->npoints-1;i++) |
---|
521 | { |
---|
522 | getPoint2d_p(pts, i, &frm); |
---|
523 | getPoint2d_p(pts, i+1, &to); |
---|
524 | dist += sqrt( ( (frm.x - to.x)*(frm.x - to.x) ) + |
---|
525 | ((frm.y - to.y)*(frm.y - to.y) ) ); |
---|
526 | } |
---|
527 | return dist; |
---|
528 | } |
---|
529 | |
---|
530 | /* |
---|
531 | * Find the 3d/2d length of the given POINTARRAY |
---|
532 | * (depending on its dimensions) |
---|
533 | */ |
---|
534 | double |
---|
535 | lwgeom_pointarray_length(POINTARRAY *pts) |
---|
536 | { |
---|
537 | double dist = 0.0; |
---|
538 | int i; |
---|
539 | POINT3DZ frm; |
---|
540 | POINT3DZ to; |
---|
541 | |
---|
542 | if ( pts->npoints < 2 ) return 0.0; |
---|
543 | |
---|
544 | /* compute 2d length if 3d is not available */ |
---|
545 | if ( ! TYPE_HASZ(pts->dims) ) return lwgeom_pointarray_length2d(pts); |
---|
546 | |
---|
547 | for (i=0; i<pts->npoints-1;i++) |
---|
548 | { |
---|
549 | getPoint3dz_p(pts, i, &frm); |
---|
550 | getPoint3dz_p(pts, i+1, &to); |
---|
551 | dist += sqrt( ( (frm.x - to.x)*(frm.x - to.x) ) + |
---|
552 | ((frm.y - to.y)*(frm.y - to.y) ) + |
---|
553 | ((frm.z - to.z)*(frm.z - to.z) ) ); |
---|
554 | } |
---|
555 | |
---|
556 | return dist; |
---|
557 | } |
---|
558 | |
---|
559 | /* |
---|
560 | * This should be rewritten to make use of the curve itself. |
---|
561 | */ |
---|
562 | double |
---|
563 | lwgeom_curvepolygon_area(LWCURVEPOLY *curvepoly) |
---|
564 | { |
---|
565 | LWPOLY *poly = (LWPOLY *)lwgeom_segmentize((LWGEOM *)curvepoly, 32); |
---|
566 | return lwgeom_polygon_area(poly); |
---|
567 | } |
---|
568 | |
---|
569 | /* |
---|
570 | * Find the area of the outer ring - sum (area of inner rings). |
---|
571 | * Could use a more numerically stable calculator... |
---|
572 | */ |
---|
573 | double |
---|
574 | lwgeom_polygon_area(LWPOLY *poly) |
---|
575 | { |
---|
576 | double poly_area=0.0; |
---|
577 | int i; |
---|
578 | POINT2D p1; |
---|
579 | POINT2D p2; |
---|
580 | |
---|
581 | LWDEBUGF(2, "in lwgeom_polygon_area (%d rings)", poly->nrings); |
---|
582 | |
---|
583 | for (i=0; i<poly->nrings; i++) |
---|
584 | { |
---|
585 | int j; |
---|
586 | POINTARRAY *ring = poly->rings[i]; |
---|
587 | double ringarea = 0.0; |
---|
588 | |
---|
589 | LWDEBUGF(4, " rings %d has %d points", i, ring->npoints); |
---|
590 | |
---|
591 | for (j=0; j<ring->npoints-1; j++) |
---|
592 | { |
---|
593 | getPoint2d_p(ring, j, &p1); |
---|
594 | getPoint2d_p(ring, j+1, &p2); |
---|
595 | ringarea += ( p1.x * p2.y ) - ( p1.y * p2.x ); |
---|
596 | } |
---|
597 | |
---|
598 | ringarea /= 2.0; |
---|
599 | |
---|
600 | LWDEBUGF(4, " ring 1 has area %lf",ringarea); |
---|
601 | |
---|
602 | ringarea = fabs(ringarea); |
---|
603 | if (i != 0) /*outer */ |
---|
604 | ringarea = -1.0*ringarea ; /* its a hole */ |
---|
605 | |
---|
606 | poly_area += ringarea; |
---|
607 | } |
---|
608 | |
---|
609 | return poly_area; |
---|
610 | } |
---|
611 | |
---|
612 | /* |
---|
613 | * Compute the sum of polygon rings length. |
---|
614 | * Could use a more numerically stable calculator... |
---|
615 | */ |
---|
616 | double lwgeom_polygon_perimeter(LWPOLY *poly) |
---|
617 | { |
---|
618 | double result=0.0; |
---|
619 | int i; |
---|
620 | |
---|
621 | LWDEBUGF(2, "in lwgeom_polygon_perimeter (%d rings)", poly->nrings); |
---|
622 | |
---|
623 | for (i=0; i<poly->nrings; i++) |
---|
624 | result += lwgeom_pointarray_length(poly->rings[i]); |
---|
625 | |
---|
626 | return result; |
---|
627 | } |
---|
628 | |
---|
629 | /* |
---|
630 | * Compute the sum of polygon rings length (forcing 2d computation). |
---|
631 | * Could use a more numerically stable calculator... |
---|
632 | */ |
---|
633 | double lwgeom_polygon_perimeter2d(LWPOLY *poly) |
---|
634 | { |
---|
635 | double result=0.0; |
---|
636 | int i; |
---|
637 | |
---|
638 | LWDEBUGF(2, "in lwgeom_polygon_perimeter (%d rings)", poly->nrings); |
---|
639 | |
---|
640 | for (i=0; i<poly->nrings; i++) |
---|
641 | result += lwgeom_pointarray_length2d(poly->rings[i]); |
---|
642 | |
---|
643 | return result; |
---|
644 | } |
---|
645 | |
---|
646 | double |
---|
647 | lwgeom_mindistance2d_recursive(uchar *lw1, uchar *lw2) |
---|
648 | { |
---|
649 | return lwgeom_mindistance2d_recursive_tolerance( lw1, lw2, 0.0 ); |
---|
650 | } |
---|
651 | |
---|
652 | double |
---|
653 | lwgeom_mindistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double tolerance) |
---|
654 | { |
---|
655 | LWGEOM_INSPECTED *in1, *in2; |
---|
656 | int i, j; |
---|
657 | double mindist = -1; |
---|
658 | |
---|
659 | in1 = lwgeom_inspect(lw1); |
---|
660 | in2 = lwgeom_inspect(lw2); |
---|
661 | |
---|
662 | for (i=0; i<in1->ngeometries; i++) |
---|
663 | { |
---|
664 | uchar *g1 = lwgeom_getsubgeometry_inspected(in1, i); |
---|
665 | int t1 = lwgeom_getType(g1[0]); |
---|
666 | double dist=tolerance; |
---|
667 | |
---|
668 | /* it's a multitype... recurse */ |
---|
669 | if ( t1 >= 4 ) |
---|
670 | { |
---|
671 | dist = lwgeom_mindistance2d_recursive_tolerance(g1, lw2, tolerance); |
---|
672 | if ( dist <= tolerance ) return tolerance; /* can't be closer */ |
---|
673 | if ( mindist == -1 ) mindist = dist; |
---|
674 | else mindist = LW_MIN(dist, mindist); |
---|
675 | continue; |
---|
676 | } |
---|
677 | |
---|
678 | for (j=0; j<in2->ngeometries; j++) |
---|
679 | { |
---|
680 | uchar *g2 = lwgeom_getsubgeometry_inspected(in2, j); |
---|
681 | int t2 = lwgeom_getType(g2[0]); |
---|
682 | |
---|
683 | if ( t1 == POINTTYPE ) |
---|
684 | { |
---|
685 | if ( t2 == POINTTYPE ) |
---|
686 | { |
---|
687 | dist = distance2d_point_point( |
---|
688 | lwpoint_deserialize(g1), |
---|
689 | lwpoint_deserialize(g2) |
---|
690 | ); |
---|
691 | } |
---|
692 | else if ( t2 == LINETYPE ) |
---|
693 | { |
---|
694 | dist = distance2d_point_line( |
---|
695 | lwpoint_deserialize(g1), |
---|
696 | lwline_deserialize(g2) |
---|
697 | ); |
---|
698 | } |
---|
699 | else if ( t2 == POLYGONTYPE ) |
---|
700 | { |
---|
701 | dist = distance2d_point_poly( |
---|
702 | lwpoint_deserialize(g1), |
---|
703 | lwpoly_deserialize(g2) |
---|
704 | ); |
---|
705 | } |
---|
706 | } |
---|
707 | else if ( t1 == LINETYPE ) |
---|
708 | { |
---|
709 | if ( t2 == POINTTYPE ) |
---|
710 | { |
---|
711 | dist = distance2d_point_line( |
---|
712 | lwpoint_deserialize(g2), |
---|
713 | lwline_deserialize(g1) |
---|
714 | ); |
---|
715 | } |
---|
716 | else if ( t2 == LINETYPE ) |
---|
717 | { |
---|
718 | dist = distance2d_line_line( |
---|
719 | lwline_deserialize(g1), |
---|
720 | lwline_deserialize(g2) |
---|
721 | ); |
---|
722 | } |
---|
723 | else if ( t2 == POLYGONTYPE ) |
---|
724 | { |
---|
725 | dist = distance2d_line_poly( |
---|
726 | lwline_deserialize(g1), |
---|
727 | lwpoly_deserialize(g2) |
---|
728 | ); |
---|
729 | } |
---|
730 | } |
---|
731 | else if ( t1 == POLYGONTYPE ) |
---|
732 | { |
---|
733 | if ( t2 == POLYGONTYPE ) |
---|
734 | { |
---|
735 | dist = distance2d_poly_poly( |
---|
736 | lwpoly_deserialize(g2), |
---|
737 | lwpoly_deserialize(g1) |
---|
738 | ); |
---|
739 | } |
---|
740 | else if ( t2 == POINTTYPE ) |
---|
741 | { |
---|
742 | dist = distance2d_point_poly( |
---|
743 | lwpoint_deserialize(g2), |
---|
744 | lwpoly_deserialize(g1) |
---|
745 | ); |
---|
746 | } |
---|
747 | else if ( t2 == LINETYPE ) |
---|
748 | { |
---|
749 | dist = distance2d_line_poly( |
---|
750 | lwline_deserialize(g2), |
---|
751 | lwpoly_deserialize(g1) |
---|
752 | ); |
---|
753 | } |
---|
754 | } |
---|
755 | else /* it's a multitype... recurse */ |
---|
756 | { |
---|
757 | dist = lwgeom_mindistance2d_recursive_tolerance(g1, g2, tolerance); |
---|
758 | } |
---|
759 | |
---|
760 | if (mindist == -1 ) mindist = dist; |
---|
761 | else mindist = LW_MIN(dist, mindist); |
---|
762 | |
---|
763 | LWDEBUGF(3, "dist %d-%d: %f - mindist: %f", |
---|
764 | i, j, dist, mindist); |
---|
765 | |
---|
766 | |
---|
767 | if (mindist <= tolerance) return tolerance; /* can't be closer */ |
---|
768 | |
---|
769 | } |
---|
770 | |
---|
771 | } |
---|
772 | |
---|
773 | if (mindist<0) mindist = 0; |
---|
774 | |
---|
775 | return mindist; |
---|
776 | } |
---|
777 | |
---|
778 | |
---|
779 | |
---|
780 | int |
---|
781 | lwgeom_pt_inside_circle(POINT2D *p, double cx, double cy, double rad) |
---|
782 | { |
---|
783 | POINT2D center; |
---|
784 | |
---|
785 | center.x = cx; |
---|
786 | center.y = cy; |
---|
787 | |
---|
788 | if ( distance2d_pt_pt(p, ¢er) < rad ) return 1; |
---|
789 | else return 0; |
---|
790 | |
---|
791 | } |
---|
792 | |
---|
793 | /* |
---|
794 | * Compute the azimuth of segment AB in radians. |
---|
795 | * Return 0 on exception (same point), 1 otherwise. |
---|
796 | */ |
---|
797 | int |
---|
798 | azimuth_pt_pt(POINT2D *A, POINT2D *B, double *d) |
---|
799 | { |
---|
800 | if ( A->x == B->x ) |
---|
801 | { |
---|
802 | if ( A->y < B->y ) *d=0.0; |
---|
803 | else if ( A->y > B->y ) *d=M_PI; |
---|
804 | else return 0; |
---|
805 | return 1; |
---|
806 | } |
---|
807 | |
---|
808 | if ( A->y == B->y ) |
---|
809 | { |
---|
810 | if ( A->x < B->x ) *d=M_PI/2; |
---|
811 | else if ( A->x > B->x ) *d=M_PI+(M_PI/2); |
---|
812 | else return 0; |
---|
813 | return 1; |
---|
814 | } |
---|
815 | |
---|
816 | if ( A->x < B->x ) |
---|
817 | { |
---|
818 | if ( A->y < B->y ) |
---|
819 | { |
---|
820 | *d=atan(fabs(A->x - B->x) / fabs(A->y - B->y) ); |
---|
821 | } |
---|
822 | else /* ( A->y > B->y ) - equality case handled above */ |
---|
823 | { |
---|
824 | *d=atan(fabs(A->y - B->y) / fabs(A->x - B->x) ) |
---|
825 | + (M_PI/2); |
---|
826 | } |
---|
827 | } |
---|
828 | |
---|
829 | else /* ( A->x > B->x ) - equality case handled above */ |
---|
830 | { |
---|
831 | if ( A->y > B->y ) |
---|
832 | { |
---|
833 | *d=atan(fabs(A->x - B->x) / fabs(A->y - B->y) ) |
---|
834 | + M_PI; |
---|
835 | } |
---|
836 | else /* ( A->y < B->y ) - equality case handled above */ |
---|
837 | { |
---|
838 | *d=atan(fabs(A->y - B->y) / fabs(A->x - B->x) ) |
---|
839 | + (M_PI+(M_PI/2)); |
---|
840 | } |
---|
841 | } |
---|
842 | |
---|
843 | return 1; |
---|
844 | } |
---|
845 | |
---|