source: trunk/liblwgeom/measures.c @ 2815

Last change on this file since 2815 was 2815, checked in by mcayland, 8 years ago

Split the basic geometry accessors into a separate static library liblwgeom.a; this potentially allows re-use of the liblwgeom functions from within PostGIS, or could be extended at a later date to include databases other than MySQL. This patch includes a change to the liblwgeom handler functions; instead of sprinkling init_pg_func()s around the source, I have changed the default liblwgeom handlers to make use of a callback to allow linked libraries to set their own handlers the first time any of them are called. I have also tidied up the parser API a little in liblwgeom.h, which means wktparse.h can be removed from all of the headers in the lwgeom/ directory, plus renamed wktunparse.c to lwgunparse.c to keep things similar to lwgparse.c. Finally, I renamed liblwgeom.c to lwutil.c to avoid confusion within the new interface. TODO: the liblwgeom Makefile has some gcc-specific options, but these can be fixed later - it seemed more important to make the warnings visible to developers.

File size: 18.3 KB
Line 
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 */
28int 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
84double 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 */
98double 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 */
146double 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 */
229double 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 */
254double 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) */
303int 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 */
330double 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
382double 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
393double 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
401double 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 */
413double 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 */
459double 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
505double 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) */
512double 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 */
534double
535lwgeom_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 */
562double
563lwgeom_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 */
573double
574lwgeom_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 */
616double 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 */
633double 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
646double
647lwgeom_mindistance2d_recursive(uchar *lw1, uchar *lw2)
648{
649  return lwgeom_mindistance2d_recursive_tolerance( lw1, lw2, 0.0 );
650}
651
652double
653lwgeom_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
780int
781lwgeom_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, &center) < 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 */
797int
798azimuth_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
Note: See TracBrowser for help on using the repository browser.