# source:trunk/liblwgeom/measures.c@5807

Last change on this file since 5807 was 5807, checked in by nicklas, 6 years ago

pushing *uchar use out of measures.c as part of #308

• Property svn:keywords set to `Author Date Id Revision`
File size: 42.7 KB
Line
1
2/**********************************************************************
3 * \$Id: measures.c 5807 2010-08-11 09:24:58Z nicklas \$
4 *
5 * PostGIS - Spatial Types for PostgreSQL
6 * http://postgis.refractions.net
7 * Copyright 2001-2006 Refractions Research Inc.
8 *
9 * This is free software; you can redistribute and/or modify it under
10 * the terms of the GNU General Public Licence. See the COPYING file.
11 *
12 **********************************************************************/
13
14#include <math.h>
15#include <string.h>
16#include <stdlib.h>
17
18#include "measures.h"
19
20
21/*------------------------------------------------------------------------------------------------------------
22Initializing functions
23The functions starting the distance-calculation processses
24--------------------------------------------------------------------------------------------------------------*/
25
26/**
27Function initializing shortestline and longestline calculations.
28*/
29LWGEOM *
30lw_dist2d_distanceline(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode)
31{
32        double x1,x2,y1,y2;
33
34        double initdistance = ( mode == DIST2D_MIN ? MAXFLOAT : -1.0);
35        DISTPTS thedl;
36        LWPOINT *lwpoints[2];
37        LWGEOM *result;
38
39        thedl.mode = mode;
40        thedl.distance = initdistance;
41        thedl.tolerance = 0.0;
42
43        LWDEBUG(2, "lw_dist2d_distanceline is called");
44
45        if (!lw_dist2d_comp( lw1,lw2,&thedl))
46        {
47                /*should never get here. all cases ought to be error handled earlier*/
48                lwerror("Some unspecified error.");
49                result =(LWGEOM *)lwcollection_construct_empty(srid, 0, 0);
50        }
51
52        /*if thedl.distance is unchanged there where only empty geometries input*/
53        if (thedl.distance == initdistance)
54        {
55                LWDEBUG(3, "didn't find geometries to measure between, returning null");
56                result =(LWGEOM *)lwcollection_construct_empty(srid, 0, 0);
57        }
58        else
59        {
60                x1=thedl.p1.x;
61                y1=thedl.p1.y;
62                x2=thedl.p2.x;
63                y2=thedl.p2.y;
64
65                lwpoints[0] = make_lwpoint2d(srid, x1, y1);
66                lwpoints[1] = make_lwpoint2d(srid, x2, y2);
67
68                result = (LWGEOM *)lwline_from_lwpointarray(srid, 2, lwpoints);
69        }
70        return result;
71}
72
73/**
74Function initializing closestpoint calculations.
75*/
76LWGEOM *
77lw_dist2d_distancepoint(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode)
78{
79        double x,y;
80        DISTPTS thedl;
81        double initdistance = MAXFLOAT;
82        LWGEOM *result;
83
84        thedl.mode = mode;
85        thedl.distance= initdistance;
86        thedl.tolerance = 0;
87
88        LWDEBUG(2, "lw_dist2d_distancepoint is called");
89
90        if (!lw_dist2d_comp( lw1,lw2,&thedl))
91        {
92                /*should never get here. all cases ought to be error handled earlier*/
93                lwerror("Some unspecified error.");
94                result =(LWGEOM *)lwcollection_construct_empty(srid, 0, 0);
95        }
96        if (thedl.distance == initdistance)
97        {
98                LWDEBUG(3, "didn't find geometries to measure between, returning null");
99                result =(LWGEOM *)lwcollection_construct_empty(srid, 0, 0);
100        }
101        else
102        {
103                x=thedl.p1.x;
104                y=thedl.p1.y;
105                result = (LWGEOM *)make_lwpoint2d(srid, x, y);
106        }
107        return result;
108}
109
110
111/**
112Function initialazing max distance calculation
113*/
114double
115lwgeom_maxdistance2d(LWGEOM *lw1, LWGEOM *lw2)
116{
117        LWDEBUG(2, "lwgeom_maxdistance2d is called");
118
119        return lwgeom_maxdistance2d_tolerance( lw1, lw2, 0.0 );
120}
121
122/**
123Function handling max distance calculations and dfyllywithin calculations.
124The difference is just the tolerance.
125*/
126double
127lwgeom_maxdistance2d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance)
128{
129        /*double thedist;*/
130        DISTPTS thedl;
131        LWDEBUG(2, "lwgeom_maxdistance2d_tolerance is called");
132        thedl.mode = DIST2D_MAX;
133        thedl.distance= -1;
134        thedl.tolerance = tolerance;
135        if (lw_dist2d_comp( lw1,lw2,&thedl))
136        {
137                return thedl.distance;
138        }
139        /*should never get here. all cases ought to be error handled earlier*/
140        lwerror("Some unspecified error.");
141        return -1;
142}
143
144/**
145        Function initialazing min distance calculation
146*/
147double
148lwgeom_mindistance2d(LWGEOM *lw1, LWGEOM *lw2)
149{
150        LWDEBUG(2, "lwgeom_mindistance2d is called");
151        return lwgeom_mindistance2d_tolerance( lw1, lw2, 0.0 );
152}
153
154/**
155        Function handling min distance calculations and dwithin calculations.
156        The difference is just the tolerance.
157*/
158double
159lwgeom_mindistance2d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance)
160{
161        DISTPTS thedl;
162        LWDEBUG(2, "lwgeom_mindistance2d_tolerance is called");
163        thedl.mode = DIST2D_MIN;
164        thedl.distance= MAXFLOAT;
165        thedl.tolerance = tolerance;
166        if (lw_dist2d_comp( lw1,lw2,&thedl))
167        {
168                return thedl.distance;
169        }
170        /*should never get here. all cases ought to be error handled earlier*/
171        lwerror("Some unspecified error.");
172        return MAXFLOAT;
173}
174
175
176/*------------------------------------------------------------------------------------------------------------
177End of Initializing functions
178--------------------------------------------------------------------------------------------------------------*/
179
180/*------------------------------------------------------------------------------------------------------------
181Preprocessing functions
182Functions preparing geometries for distance-calculations
183--------------------------------------------------------------------------------------------------------------*/
184
185/**
186        This function just deserializes geometries
187        Bboxes is not checked here since it is the subgeometries
188        bboxes we will use anyway.
189*/
190int
191lw_dist2d_comp(LWGEOM *lw1, LWGEOM *lw2, DISTPTS *dl)
192{
193        LWDEBUG(2, "lw_dist2d_comp is called");
194
195        return lw_dist2d_recursive      ((LWCOLLECTION*) lw1,(LWCOLLECTION*) lw2, dl);
196}
197
198/**
199This is a recursive function delivering every possible combinatin of subgeometries
200*/
201int lw_dist2d_recursive(const LWCOLLECTION * lwg1,const LWCOLLECTION * lwg2,DISTPTS *dl)
202{
203        int i, j;
204        int n1=1;
205        int n2=1;
206        LWGEOM *g1, *g2;
207
208        LWDEBUGF(2, "lw_dist2d_comp is called with type1=%d, type2=%d", TYPE_GETTYPE(lwg1->type), TYPE_GETTYPE(lwg2->type));
209
210        if (lwgeom_is_collection(TYPE_GETTYPE(lwg1->type)))
211        {
212                LWDEBUG(3, "First geometry is collection");
213                n1=lwg1->ngeoms;
214        }
215        if (lwgeom_is_collection(TYPE_GETTYPE(lwg2->type)))
216        {
217                LWDEBUG(3, "Second geometry is collection");
218                n2=lwg2->ngeoms;
219        }
220
221        for ( i = 0; i < n1; i++ )
222        {
223
224                if (lwgeom_is_collection(TYPE_GETTYPE(lwg1->type)))
225                {
226                        g1=lwg1->geoms[i];
227                }
228                else
229                {
230                        g1=(LWGEOM*)lwg1;
231                }
232
233                if (lwgeom_is_empty(g1)) return LW_TRUE;
234
235                if (lwgeom_is_collection(TYPE_GETTYPE(g1->type)))
236                {
237                        LWDEBUG(3, "Found collection inside first geometry collection, recursing");
238                        if (!lw_dist2d_recursive((LWCOLLECTION*)g1, (LWCOLLECTION*)lwg2, dl)) return LW_FALSE;
239                        continue;
240                }
241                for ( j = 0; j < n2; j++ )
242                {
243                        if (lwgeom_is_collection(TYPE_GETTYPE(lwg2->type)))
244                        {
245                                g2=lwg2->geoms[j];
246                        }
247                        else
248                        {
249                                g2=(LWGEOM*)lwg2;
250                        }
251                        if (lwgeom_is_collection(TYPE_GETTYPE(g2->type)))
252                        {
253                                LWDEBUG(3, "Found collection inside second geometry collection, recursing");
254                                if (!lw_dist2d_recursive((LWCOLLECTION*) g1, (LWCOLLECTION*)g2, dl)) return LW_FALSE;
255                                continue;
256                        }
257
258                        if ( ! g1->bbox )
259                        {
260                                g1->bbox = lwgeom_compute_box2d(g1);
261                        }
262                        if ( ! g2->bbox )
263                        {
264                                g2->bbox = lwgeom_compute_box2d(g2);
265                        }
266                        /*If one of geometries is empty, return. True here only means continue searching. False would have stoped the process*/
267                        if (lwgeom_is_empty(g1)||lwgeom_is_empty(g2)) return LW_TRUE;
268
269                        if ((dl->mode==DIST2D_MAX)||(TYPE_GETTYPE(g1->type)==POINTTYPE) ||(TYPE_GETTYPE(g2->type)==POINTTYPE)||lw_dist2d_check_overlap(g1, g2))
270                        {
271                                if (!lw_dist2d_distribute_bruteforce(g1, g2, dl)) return LW_FALSE;
272                                if (dl->distance<=dl->tolerance && dl->mode == DIST2D_MIN) return LW_TRUE; /*just a check if  the answer is already given*/
273                        }
274                        else
275                        {
276                                if (!lw_dist2d_distribute_fast(g1, g2, dl)) return LW_FALSE;
277                        }
278                }
279        }
280        return LW_TRUE;
281}
282
283
284
285/**
286
287This function distributes the "old-type" brut-force tasks depending on type
288*/
289int
290lw_dist2d_distribute_bruteforce(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS *dl)
291{
292
293        int     t1 = TYPE_GETTYPE(lwg1->type);
294        int     t2 = TYPE_GETTYPE(lwg2->type);
295
296        LWDEBUGF(2, "lw_dist2d_distribute_bruteforce is called with typ1=%d, type2=%d", TYPE_GETTYPE(lwg1->type), TYPE_GETTYPE(lwg2->type));
297
298        if  ( t1 == POINTTYPE )
299        {
300                if  ( t2 == POINTTYPE )
301                {
302                        dl->twisted=1;
303                        return lw_dist2d_point_point((LWPOINT *)lwg1, (LWPOINT *)lwg2, dl);
304                }
305                else if  ( t2 == LINETYPE )
306                {
307                        dl->twisted=1;
308                        return lw_dist2d_point_line((LWPOINT *)lwg1, (LWLINE *)lwg2, dl);
309                }
310                else if  ( t2 == POLYGONTYPE )
311                {
312                        dl->twisted=1;
313                        return lw_dist2d_point_poly((LWPOINT *)lwg1, (LWPOLY *)lwg2,dl);
314                }
315                else
316                {
317                        lwerror("Unsupported geometry type: %s", lwtype_name(t2));
318                        return LW_FALSE;
319                }
320        }
321        else if ( t1 == LINETYPE )
322        {
323                if ( t2 == POINTTYPE )
324                {
325                        dl->twisted=(-1);
326                        return lw_dist2d_point_line((LWPOINT *)lwg2,(LWLINE *)lwg1,dl);
327                }
328                else if ( t2 == LINETYPE )
329                {
330                        dl->twisted=1;
331                        /*lwnotice("start line");*/
332                        return lw_dist2d_line_line((LWLINE *)lwg1,(LWLINE *)lwg2,dl);
333                }
334                else if ( t2 == POLYGONTYPE )
335                {
336                        dl->twisted=1;
337                        return lw_dist2d_line_poly((LWLINE *)lwg1,(LWPOLY *)lwg2,dl);
338                }
339                else
340                {
341                        lwerror("Unsupported geometry type: %s", lwtype_name(t2));
342                        return LW_FALSE;
343                }
344        }
345        else if ( t1 == POLYGONTYPE )
346        {
347                if ( t2 == POLYGONTYPE )
348                {
349                        dl->twisted=(1);
350                        return lw_dist2d_poly_poly((LWPOLY *)lwg1,(LWPOLY *)lwg2,dl);
351                }
352                else if ( t2 == POINTTYPE )
353                {
354                        dl->twisted=(-1);
355                        return lw_dist2d_point_poly((LWPOINT *)lwg2,  (LWPOLY *)lwg1, dl);
356                }
357                else if ( t2 == LINETYPE )
358                {
359                        dl->twisted=(-1);
360                        return lw_dist2d_line_poly((LWLINE *)lwg2,(LWPOLY *)lwg1,dl);
361                }
362                else
363                {
364                        lwerror("Unsupported geometry type: %s", lwtype_name(t2));
365                        return LW_FALSE;
366                }
367        }
368        else
369        {
370                lwerror("Unsupported geometry type: %s", lwtype_name(t1));
371                return LW_FALSE;
372        }
373        /*You shouldn't being able to get here*/
374        lwerror("unspecified error in function lw_dist2d_distribute_bruteforce");
375        return LW_FALSE;
376}
377
378/**
379
380We have to check for overlapping bboxes
381*/
382int
383lw_dist2d_check_overlap(LWGEOM *lwg1,LWGEOM *lwg2)
384{
385        LWDEBUG(2, "lw_dist2d_check_overlap is called");
386        if ( ! lwg1->bbox )
387                lwg1->bbox = lwgeom_compute_box2d(lwg1);
388        if ( ! lwg2->bbox )
389                lwg2->bbox = lwgeom_compute_box2d(lwg2);
390
391        /*Check if the geometries intersect.
392        */
393        if ((lwg1->bbox->xmax<lwg2->bbox->xmin||lwg1->bbox->xmin>lwg2->bbox->xmax||lwg1->bbox->ymax<lwg2->bbox->ymin||lwg1->bbox->ymin>lwg2->bbox->ymax))
394        {
395                LWDEBUG(3, "geometries bboxes did not overlap");
396                return LW_FALSE;
397        }
398        LWDEBUG(3, "geometries bboxes overlap");
399        return LW_TRUE;
400}
401
402/**
403
404Here the geometries are distributed for the new faster distance-calculations
405*/
406int
407lw_dist2d_distribute_fast(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS *dl)
408{
409        POINTARRAY *pa1, *pa2;
410        int     type1 = TYPE_GETTYPE(lwg1->type);
411        int     type2 = TYPE_GETTYPE(lwg2->type);
412
413        LWDEBUGF(2, "lw_dist2d_distribute_fast is called with typ1=%d, type2=%d", TYPE_GETTYPE(lwg1->type), TYPE_GETTYPE(lwg2->type));
414
415        switch (type1)
416        {
417        case LINETYPE:
418                pa1 = ((LWLINE *)lwg1)->points;
419                break;
420        case POLYGONTYPE:
421                pa1 = ((LWPOLY *)lwg1)->rings[0];
422                break;
423        default:
424                lwerror("Unsupported geometry1 type: %s", lwtype_name(type1));
425                return LW_FALSE;
426        }
427        switch (type2)
428        {
429        case LINETYPE:
430                pa2 = ((LWLINE *)lwg2)->points;
431                break;
432        case POLYGONTYPE:
433                pa2 = ((LWPOLY *)lwg2)->rings[0];
434                break;
435        default:
436                lwerror("Unsupported geometry2 type: %s", lwtype_name(type1));
437                return LW_FALSE;
438        }
439        dl->twisted=1;
440        return lw_dist2d_fast_ptarray_ptarray(pa1, pa2, dl, lwg1->bbox, lwg2->bbox);
441}
442
443/*------------------------------------------------------------------------------------------------------------
444End of Preprocessing functions
445--------------------------------------------------------------------------------------------------------------*/
446
447
448/*------------------------------------------------------------------------------------------------------------
449Brute force functions
450The old way of calculating distances, now used for:
4511)      distances to points (because there shouldn't be anything to gain by the new way of doing it)
4522)      distances when subgeometries geometries bboxes overlaps
453--------------------------------------------------------------------------------------------------------------*/
454
455/**
456
457point to point calculation
458*/
459int
460lw_dist2d_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS *dl)
461{
462        POINT2D p1;
463        POINT2D p2;
464
465        getPoint2d_p(point1->point, 0, &p1);
466        getPoint2d_p(point2->point, 0, &p2);
467
468        return lw_dist2d_pt_pt(&p1, &p2,dl);
469}
470/**
471
472point to line calculation
473*/
474int
475lw_dist2d_point_line(LWPOINT *point, LWLINE *line, DISTPTS *dl)
476{
477        POINT2D p;
478        POINTARRAY *pa = line->points;
479        LWDEBUG(2, "lw_dist2d_point_line is called");
480        getPoint2d_p(point->point, 0, &p);
481        return lw_dist2d_pt_ptarray(&p, pa, dl);
482}
483/**
484 * 1. see if pt in outer boundary. if no, then treat the outer ring like a line
485 * 2. if in the boundary, test to see if its in a hole.
486 *    if so, then return dist to hole, else return 0 (point in polygon)
487 */
488int
489lw_dist2d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS *dl)
490{
491        POINT2D p;
492        int i;
493
494        LWDEBUG(2, "lw_dist2d_point_poly called");
495
496        getPoint2d_p(point->point, 0, &p);
497
498        if (dl->mode == DIST2D_MAX)
499        {
500                LWDEBUG(3, "looking for maxdistance");
501                return lw_dist2d_pt_ptarray(&p, poly->rings[0], dl);
502        }
503        /* Return distance to outer ring if not inside it */
504        if ( ! pt_in_ring_2d(&p, poly->rings[0]) )
505        {
506                LWDEBUG(3, "first point not inside outer-ring");
507                return lw_dist2d_pt_ptarray(&p, poly->rings[0], dl);
508        }
509
510        /*
511         * Inside the outer ring.
512         * Scan though each of the inner rings looking to
513         * see if its inside.  If not, distance==0.
514         * Otherwise, distance = pt to ring distance
515         */
516        for (i=1; i<poly->nrings; i++)
517        {
518                /* Inside a hole. Distance = pt -> ring */
519                if ( pt_in_ring_2d(&p, poly->rings[i]) )
520                {
521                        LWDEBUG(3, " inside an hole");
522                        return lw_dist2d_pt_ptarray(&p, poly->rings[i], dl);
523                }
524        }
525
526        LWDEBUG(3, " inside the polygon");
527        if (dl->mode == DIST2D_MIN)
528        {
529                dl->distance=0.0;
530                dl->p1.x=p.x;
531                dl->p1.y=p.y;
532                dl->p2.x=p.x;
533                dl->p2.y=p.y;
534        }
535        return LW_TRUE; /* Is inside the polygon */
536}
537/**
538
539line to line calculation
540*/
541int
542lw_dist2d_line_line(LWLINE *line1, LWLINE *line2, DISTPTS *dl)
543{
544        POINTARRAY *pa1 = line1->points;
545        POINTARRAY *pa2 = line2->points;
546        LWDEBUG(2, "lw_dist2d_line_line is called");
547        return lw_dist2d_ptarray_ptarray(pa1, pa2, dl);
548}
549/**
550
551line to polygon calculation
552*/
553int
554lw_dist2d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS *dl)
555{
556        LWDEBUG(2, "lw_dist2d_line_poly is called");
557        return lw_dist2d_ptarray_poly(line->points, poly, dl);
558}
559
560/**
561
562Function handling polygon to polygon calculation
5631       if we are looking for maxdistance, just check the outer rings.
5642       check if poly1 has first point outside poly2 and vice versa, if so, just check outer rings
5653       check if first point of poly2 is in a hole of poly1. If so check outer ring of poly2 against that hole of poly1
5664       check if first point of poly1 is in a hole of poly2. If so check outer ring of poly1 against that hole of poly2
5675       If we have come all the way here we know that the first point of one of them is inside the other ones outer ring and not in holes so we check wich one is inside.
568 */
569int
570lw_dist2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS *dl)
571{
572
573        POINT2D pt;
574        int i;
575
576        LWDEBUG(2, "lw_dist2d_poly_poly called");
577
578        /*1     if we are looking for maxdistance, just check the outer rings.*/
579        if (dl->mode == DIST2D_MAX)
580        {
581                return lw_dist2d_ptarray_ptarray(poly1->rings[0],       poly2->rings[0],dl);
582        }
583
584
585        /* 2    check if poly1 has first point outside poly2 and vice versa, if so, just check outer rings
586        here it would be possible to handle the information about wich one is inside wich one and only search for the smaller ones in the bigger ones holes.*/
587        getPoint2d_p(poly1->rings[0], 0, &pt);
588        if ( !pt_in_ring_2d(&pt, poly2->rings[0]))
589        {
590                getPoint2d_p(poly2->rings[0], 0, &pt);
591                if (!pt_in_ring_2d(&pt, poly1->rings[0]))
592                {
593                        return lw_dist2d_ptarray_ptarray(poly1->rings[0],       poly2->rings[0],dl);
594                }
595        }
596
597        /*3     check if first point of poly2 is in a hole of poly1. If so check outer ring of poly2 against that hole of poly1*/
598        getPoint2d_p(poly2->rings[0], 0, &pt);
599        for (i=1; i<poly1->nrings; i++)
600        {
601                /* Inside a hole */
602                if ( pt_in_ring_2d(&pt, poly1->rings[i]) )
603                {
604                        return lw_dist2d_ptarray_ptarray(poly1->rings[i],       poly2->rings[0],dl);
605                }
606        }
607
608        /*4     check if first point of poly1 is in a hole of poly2. If so check outer ring of poly1 against that hole of poly2*/
609        getPoint2d_p(poly1->rings[0], 0, &pt);
610        for (i=1; i<poly2->nrings; i++)
611        {
612                /* Inside a hole */
613                if ( pt_in_ring_2d(&pt, poly2->rings[i]) )
614                {
615                        return lw_dist2d_ptarray_ptarray(poly1->rings[0],       poly2->rings[i],dl);
616                }
617        }
618
619
620        /*5     If we have come all the way here we know that the first point of one of them is inside the other ones outer ring and not in holes so we check wich one is inside.*/
621        getPoint2d_p(poly1->rings[0], 0, &pt);
622        if ( pt_in_ring_2d(&pt, poly2->rings[0]))
623        {
624                dl->distance=0.0;
625                dl->p1.x=pt.x;
626                dl->p1.y=pt.y;
627                dl->p2.x=pt.x;
628                dl->p2.y=pt.y;
629                return LW_TRUE;
630        }
631
632        getPoint2d_p(poly2->rings[0], 0, &pt);
633        if (pt_in_ring_2d(&pt, poly1->rings[0]))
634        {
635                dl->distance=0.0;
636                dl->p1.x=pt.x;
637                dl->p1.y=pt.y;
638                dl->p2.x=pt.x;
639                dl->p2.y=pt.y;
640                return LW_TRUE;
641        }
642
643
644        lwerror("Unspecified error in function lw_dist2d_poly_poly");
645        return LW_FALSE;
646}
647
648/**
649
650 * search all the segments of pointarray to see which one is closest to p1
651 * Returns minimum distance between point and pointarray
652 */
653int
654lw_dist2d_pt_ptarray(POINT2D *p, POINTARRAY *pa,DISTPTS *dl)
655{
656        int t;
657        POINT2D start, end;
658        int twist = dl->twisted;
659
660        LWDEBUG(2, "lw_dist2d_pt_ptarray is called");
661
662        getPoint2d_p(pa, 0, &start);
663
664        for (t=1; t<pa->npoints; t++)
665        {
666                dl->twisted=twist;
667                getPoint2d_p(pa, t, &end);
668                if (!lw_dist2d_pt_seg(p, &start, &end,dl)) return LW_FALSE;
669
670                if (dl->distance<=dl->tolerance && dl->mode == DIST2D_MIN) return LW_TRUE; /*just a check if  the answer is already given*/
671                start = end;
672        }
673
674        return LW_TRUE;
675}
676
677/**
678 * Brute force.
679 * Test line-ring distance against each ring.
680 * If there's an intersection (distance==0) then return 0 (crosses boundary).
681 * Otherwise, test to see if any point is inside outer rings of polygon,
682 * but not in inner rings.
683 * If so, return 0  (line inside polygon),
684 * otherwise return min distance to a ring (could be outside
685 * polygon or inside a hole)
686 */
687
688int
689lw_dist2d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, DISTPTS *dl)
690{
691        POINT2D pt;
692        int i;
693
694        LWDEBUGF(2, "lw_dist2d_ptarray_poly called (%d rings)", poly->nrings);
695
696        getPoint2d_p(pa, 0, &pt);
697        if ( !pt_in_ring_2d(&pt, poly->rings[0]))
698        {
699                return lw_dist2d_ptarray_ptarray(pa,poly->rings[0],dl);
700        }
701
702        for (i=1; i<poly->nrings; i++)
703        {
704                if (!lw_dist2d_ptarray_ptarray(pa, poly->rings[i], dl)) return LW_FALSE;
705
706                LWDEBUGF(3, " distance from ring %d: %f, mindist: %f",
707                         i, dl->distance, dl->tolerance);
708                if (dl->distance<=dl->tolerance && dl->mode == DIST2D_MIN) return LW_TRUE; /*just a check if  the answer is already given*/
709        }
710
711        /*
712         * No intersection, have to check if a point is
713         * inside polygon
714         */
715        getPoint2d_p(pa, 0, &pt);
716
717        /*
718         * Outside outer ring, so min distance to a ring
719         * is the actual min distance
720
721        if ( ! pt_in_ring_2d(&pt, poly->rings[0]) )
722        {
723                return ;
724        } */
725
726        /*
727         * Its in the outer ring.
728         * Have to check if its inside a hole
729         */
730        for (i=1; i<poly->nrings; i++)
731        {
732                if ( pt_in_ring_2d(&pt, poly->rings[i]) )
733                {
734                        /*
735                         * Its inside a hole, then the actual
736                         * distance is the min ring distance
737                         */
738                        return LW_TRUE;
739                }
740        }
741        if (dl->mode == DIST2D_MIN)
742        {
743                dl->distance=0.0;
744                dl->p1.x=pt.x;
745                dl->p1.y=pt.y;
746                dl->p2.x=pt.x;
747                dl->p2.y=pt.y;
748        }
749        return LW_TRUE; /* Not in hole, so inside polygon */
750}
751
752
753
754/**
755
756test each segment of l1 against each segment of l2.
757*/
758int
759lw_dist2d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2,DISTPTS *dl)
760{
761        int t,u;
762        POINT2D start, end;
763        POINT2D start2, end2;
764        int twist = dl->twisted;
765
766        LWDEBUGF(2, "lw_dist2d_ptarray_ptarray called (points: %d-%d)",l1->npoints, l2->npoints);
767
768        if (dl->mode == DIST2D_MAX)/*If we are searching for maxdistance we go straight to point-point calculation since the maxdistance have to be between two vertexes*/
769        {
770                for (t=0; t<l1->npoints; t++) /*for each segment in L1 */
771                {
772                        getPoint2d_p(l1, t, &start);
773                        for (u=0; u<l2->npoints; u++) /*for each segment in L2 */
774                        {
775                                getPoint2d_p(l2, u, &start2);
776                                lw_dist2d_pt_pt(&start,&start2,dl);
777                                LWDEBUGF(4, "maxdist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);
778                                LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",
779                                         t, u, dl->distance, dl->tolerance);
780                        }
781                }
782        }
783        else
784        {
785                getPoint2d_p(l1, 0, &start);
786                for (t=1; t<l1->npoints; t++) /*for each segment in L1 */
787                {
788                        getPoint2d_p(l1, t, &end);
789                        getPoint2d_p(l2, 0, &start2);
790                        for (u=1; u<l2->npoints; u++) /*for each segment in L2 */
791                        {
792                                getPoint2d_p(l2, u, &end2);
793                                dl->twisted=twist;
794                                lw_dist2d_seg_seg(&start, &end, &start2, &end2,dl);
795                                LWDEBUGF(4, "mindist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);
796                                LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",
797                                         t, u, dl->distance, dl->tolerance);
798                                if (dl->distance<=dl->tolerance && dl->mode == DIST2D_MIN) return LW_TRUE; /*just a check if  the answer is already given*/
799                                start2 = end2;
800                        }
801                        start = end;
802                }
803        }
804        return LW_TRUE;
805}
806
807/**
808
809Finds the shortest distance between two segments.
810This function is changed so it is not doing any comparasion of distance
811but just sending every possible combination further to lw_dist2d_pt_seg
812*/
813int
814lw_dist2d_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D, DISTPTS *dl)
815{
816        double  s_top, s_bot,s;
817        double  r_top, r_bot,r;
818
819        LWDEBUGF(2, "lw_dist2d_seg_seg [%g,%g]->[%g,%g] by [%g,%g]->[%g,%g]",
820                 A->x,A->y,B->x,B->y, C->x,C->y, D->x, D->y);
821
822        /*A and B are the same point */
823        if (  ( A->x == B->x) && (A->y == B->y) )
824        {
825                return lw_dist2d_pt_seg(A,C,D,dl);
826        }
827        /*U and V are the same point */
828
829        if (  ( C->x == D->x) && (C->y == D->y) )
830        {
831                dl->twisted= ((dl->twisted) * (-1));
832                return lw_dist2d_pt_seg(D,A,B,dl);
833        }
834        /* AB and CD are line segments */
835        /* from comp.graphics.algo
836
837        Solving the above for r and s yields
838                                (Ay-Cy)(Dx-Cx)-(Ax-Cx)(Dy-Cy)
839                   r = ----------------------------- (eqn 1)
840                                (Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)
841
842                        (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
843                s = ----------------------------- (eqn 2)
844                        (Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)
845        Let P be the position vector of the intersection point, then
846                P=A+r(B-A) or
847                Px=Ax+r(Bx-Ax)
848                Py=Ay+r(By-Ay)
849        By examining the values of r & s, you can also determine some other limiting conditions:
850                If 0<=r<=1 & 0<=s<=1, intersection exists
851                r<0 or r>1 or s<0 or s>1 line segments do not intersect
852                If the denominator in eqn 1 is zero, AB & CD are parallel
853                If the numerator in eqn 1 is also zero, AB & CD are collinear.
854
855        */
856        r_top = (A->y-C->y)*(D->x-C->x) - (A->x-C->x)*(D->y-C->y) ;
857        r_bot = (B->x-A->x)*(D->y-C->y) - (B->y-A->y)*(D->x-C->x) ;
858
859        s_top = (A->y-C->y)*(B->x-A->x) - (A->x-C->x)*(B->y-A->y);
860        s_bot = (B->x-A->x)*(D->y-C->y) - (B->y-A->y)*(D->x-C->x);
861
862        if  ( (r_bot==0) || (s_bot == 0) )
863        {
864                if ((lw_dist2d_pt_seg(A,C,D,dl)) && (lw_dist2d_pt_seg(B,C,D,dl)))
865                {
866                        dl->twisted= ((dl->twisted) * (-1));  /*here we change the order of inputted geometrys and that we  notice by changing sign on dl->twisted*/
867                        return ((lw_dist2d_pt_seg(C,A,B,dl)) && (lw_dist2d_pt_seg(D,A,B,dl))); /*if all is successful we return true*/
868                }
869                else
870                {
871                        return LW_FALSE; /* if any of the calls to lw_dist2d_pt_seg goes wrong we return false*/
872                }
873        }
874
875        s = s_top/s_bot;
876        r=  r_top/r_bot;
877
878        if (((r<0) || (r>1) || (s<0) || (s>1)) || (dl->mode == DIST2D_MAX))
879        {
880                if ((lw_dist2d_pt_seg(A,C,D,dl)) && (lw_dist2d_pt_seg(B,C,D,dl)))
881                {
882                        dl->twisted= ((dl->twisted) * (-1));  /*here we change the order of inputted geometrys and that we  notice by changing sign on dl->twisted*/
883                        return ((lw_dist2d_pt_seg(C,A,B,dl)) && (lw_dist2d_pt_seg(D,A,B,dl))); /*if all is successful we return true*/
884                }
885                else
886                {
887                        return LW_FALSE; /* if any of the calls to lw_dist2d_pt_seg goes wrong we return false*/
888                }
889        }
890        else
891        {
892                if (dl->mode == DIST2D_MIN)     /*If there is intersection we identify the intersection point and return it but only if we are looking for mindistance*/
893                {
894                        POINT2D theP;
895
896                        if (((A->x==C->x)&&(A->y==C->y))||((A->x==D->x)&&(A->y==D->y)))
897                        {
898                                theP.x = A->x;
899                                theP.y = A->y;
900                        }
901                        else if (((B->x==C->x)&&(B->y==C->y))||((B->x==D->x)&&(B->y==D->y)))
902                        {
903                                theP.x = B->x;
904                                theP.y = B->y;
905                        }
906                        else
907                        {
908                                theP.x = A->x+r*(B->x-A->x);
909                                theP.y = A->y+r*(B->y-A->y);
910                        }
911                        dl->distance=0.0;
912                        dl->p1=theP;
913                        dl->p2=theP;
914                }
915                return LW_TRUE;
916
917        }
918        lwerror("unspecified error in function lw_dist2d_seg_seg");
919        return LW_FALSE; /*If we have come here something is wrong*/
920}
921
922
923/*------------------------------------------------------------------------------------------------------------
924End of Brute force functions
925--------------------------------------------------------------------------------------------------------------*/
926
927
928/*------------------------------------------------------------------------------------------------------------
929New faster distance calculations
930--------------------------------------------------------------------------------------------------------------*/
931
932/**
933
934The new faster calculation comparing pointarray to another pointarray
935the arrays can come from both polygons and linestrings.
936The naming is not good but comes from that it compares a
937chosen selection of the points not all of them
938*/
939int
940lw_dist2d_fast_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2,DISTPTS *dl, BOX2DFLOAT4 *box1, BOX2DFLOAT4 *box2)
941{
942        /*here we define two lists to hold our calculated "z"-values and the order number in the geometry*/
943
944        double k, thevalue;
945        float   deltaX, deltaY, c1m, c2m;
946        POINT2D theP,c1, c2;
947        float min1X, max1X, max1Y, min1Y,min2X, max2X, max2Y, min2Y;
948        int t;
949        int n1 = l1->npoints;
950        int n2 = l2->npoints;
951        LISTSTRUCT list1[n1]; /* This could be a problem, var-size array not legal in C90 */
952        LISTSTRUCT list2[n2]; /* This could be a problem, var-size array not legal in C90 */
953
954        LWDEBUG(2, "lw_dist2d_fast_ptarray_ptarray is called");
955
956        max1X = box1->xmax;
957        min1X = box1->xmin;
958        max1Y = box1->ymax;
959        min1Y = box1->ymin;
960        max2X = box2->xmax;
961        min2X = box2->xmin;
962        max2Y = box2->ymax;
963        min2Y = box2->ymin;
964        /*we want the center of the bboxes, and calculate the slope between the centerpoints*/
965        c1.x = min1X + (max1X-min1X)/2;
966        c1.y = min1Y + (max1Y-min1Y)/2;
967        c2.x = min2X + (max2X-min2X)/2;
968        c2.y = min2Y + (max2Y-min2Y)/2;
969
970        deltaX=(c2.x-c1.x);
971        deltaY=(c2.y-c1.y);
972
973
974        /*Here we calculate where the line perpendicular to the center-center line crosses the axes for each vertex
975        if the center-center line is vertical the perpendicular line will be horizontal and we find it's crossing the Y-axes with z = y-kx */
976        if ((deltaX*deltaX)<(deltaY*deltaY))        /*North or South*/
977        {
978                k = -deltaX/deltaY;
979                for (t=0; t<n1; t++) /*for each segment in L1 */
980                {
981                        getPoint2d_p(l1, t, &theP);
982                        thevalue = theP.y-(k*theP.x);
983                        list1[t].themeasure=thevalue;
984                        list1[t].pnr=t;
985
986                }
987                for (t=0; t<n2; t++) /*for each segment in L2*/
988                {
989                        getPoint2d_p(l2, t, &theP);
990                        thevalue = theP.y-(k*theP.x);
991                        list2[t].themeasure=thevalue;
992                        list2[t].pnr=t;
993
994                }
995                c1m = c1.y-(k*c1.x);
996                c2m = c2.y-(k*c2.x);
997        }
998
999
1000        /*if the center-center line is horizontal the perpendicular line will be vertical. To eliminate problems with deviding by zero we are here mirroring the coordinate-system
1001         and we find it's crossing the X-axes with z = x-(1/k)y */
1002        else        /*West or East*/
1003        {
1004                k = -deltaY/deltaX;
1005                for (t=0; t<n1; t++) /*for each segment in L1 */
1006                {
1007                        getPoint2d_p(l1, t, &theP);
1008                        thevalue = theP.x-(k*theP.y);
1009                        list1[t].themeasure=thevalue;
1010                        list1[t].pnr=t;
1011                        /* lwnotice("l1 %d, measure=%f",t,thevalue ); */
1012                }
1013                for (t=0; t<n2; t++) /*for each segment in L2*/
1014                {
1015                        getPoint2d_p(l2, t, &theP);
1016
1017                        thevalue = theP.x-(k*theP.y);
1018                        list2[t].themeasure=thevalue;
1019                        list2[t].pnr=t;
1020                        /* lwnotice("l2 %d, measure=%f",t,thevalue ); */
1021                }
1022                c1m = c1.x-(k*c1.y);
1023                c2m = c2.x-(k*c2.y);
1024        }
1025
1026        /*we sort our lists by the calculated values*/
1027        qsort(list1, n1, sizeof(LISTSTRUCT), struct_cmp_by_measure);
1028        qsort(list2, n2, sizeof(LISTSTRUCT), struct_cmp_by_measure);
1029
1030        if (c1m < c2m)
1031        {
1032                if (!lw_dist2d_pre_seg_seg(l1,l2,list1,list2,k,dl)) return LW_FALSE;
1033        }
1034        else
1035        {
1036                dl->twisted= ((dl->twisted) * (-1));
1037                if (!lw_dist2d_pre_seg_seg(l2,l1,list2,list1,k,dl)) return LW_FALSE;
1038        }
1039        return LW_TRUE;
1040}
1041
1042int
1043struct_cmp_by_measure(const void *a, const void *b)
1044{
1045        LISTSTRUCT *ia = (LISTSTRUCT*)a;
1046        LISTSTRUCT *ib = (LISTSTRUCT*)b;
1047        return ( ia->themeasure>ib->themeasure ) ? 1 : -1;
1048}
1049
1050/**
1051        preparation before lw_dist2d_seg_seg.
1052*/
1053int
1054lw_dist2d_pre_seg_seg(POINTARRAY *l1, POINTARRAY *l2,LISTSTRUCT *list1, LISTSTRUCT *list2,double k, DISTPTS *dl)
1055{
1056        POINT2D p1, p2, p3, p4, p01, p02;
1057        int pnr1,pnr2,pnr3,pnr4, n1, n2, i, u, r, twist;
1058        double maxmeasure;
1059        n1=     l1->npoints;
1060        n2 = l2->npoints;
1061
1062        LWDEBUG(2, "lw_dist2d_pre_seg_seg is called");
1063
1064        getPoint2d_p(l1, list1[0].pnr, &p1);
1065        getPoint2d_p(l2, list2[0].pnr, &p3);
1066        lw_dist2d_pt_pt(&p1, &p3,dl);
1067        maxmeasure = sqrt(dl->distance*dl->distance + (dl->distance*dl->distance*k*k));
1068        twist = dl->twisted; /*to keep the incomming order between iterations*/
1069        for (i =(n1-1); i>=0; --i)
1070        {
1071                /*we break this iteration when we have checked every
1072                point closer to our perpendicular "checkline" than
1073                our shortest found distance*/
1074                if (((list2[0].themeasure-list1[i].themeasure)) > maxmeasure) break;
1075                for (r=-1; r<=1; r +=2) /*because we are not iterating in the original pointorder we have to check the segment before and after every point*/
1076                {
1077                        pnr1 = list1[i].pnr;
1078                        getPoint2d_p(l1, pnr1, &p1);
1079                        if (pnr1+r<0)
1080                        {
1081                                getPoint2d_p(l1, (n1-1), &p01);
1082                                if (( p1.x == p01.x) && (p1.y == p01.y)) pnr2 = (n1-1);
1083                                else pnr2 = pnr1; /* if it is a line and the last and first point is not the same we avoid the edge between start and end this way*/
1084                        }
1085
1086                        else if (pnr1+r>(n1-1))
1087                        {
1088                                getPoint2d_p(l1, 0, &p01);
1089                                if (( p1.x == p01.x) && (p1.y == p01.y)) pnr2 = 0;
1090                                else pnr2 = pnr1; /* if it is a line and the last and first point is not the same we avoid the edge between start and end this way*/
1091                        }
1092                        else pnr2 = pnr1+r;
1093
1094
1095                        getPoint2d_p(l1, pnr2, &p2);
1096                        for (u=0; u<n2; ++u)
1097                        {
1098                                if (((list2[u].themeasure-list1[i].themeasure)) >= maxmeasure) break;
1099                                pnr3 = list2[u].pnr;
1100                                getPoint2d_p(l2, pnr3, &p3);
1101                                if (pnr3==0)
1102                                {
1103                                        getPoint2d_p(l2, (n2-1), &p02);
1104                                        if (( p3.x == p02.x) && (p3.y == p02.y)) pnr4 = (n2-1);
1105                                        else pnr4 = pnr3; /* if it is a line and the last and first point is not the same we avoid the edge between start and end this way*/
1106                                }
1107                                else pnr4 = pnr3-1;
1108
1109                                getPoint2d_p(l2, pnr4, &p4);
1110                                dl->twisted=twist;
1111                                if (!lw_dist2d_selected_seg_seg(&p1, &p2, &p3, &p4, dl)) return LW_FALSE;
1112
1113                                if (pnr3>=(n2-1))
1114                                {
1115                                        getPoint2d_p(l2, 0, &p02);
1116                                        if (( p3.x == p02.x) && (p3.y == p02.y)) pnr4 = 0;
1117                                        else pnr4 = pnr3; /* if it is a line and the last and first point is not the same we avoid the edge between start and end this way*/
1118                                }
1119
1120                                else pnr4 = pnr3+1;
1121
1122                                getPoint2d_p(l2, pnr4, &p4);
1123                                dl->twisted=twist; /*we reset the "twist" for each iteration*/
1124                                if (!lw_dist2d_selected_seg_seg(&p1, &p2, &p3, &p4, dl)) return LW_FALSE;
1125
1126                                maxmeasure = sqrt(dl->distance*dl->distance + (dl->distance*dl->distance*k*k));/*here we "translate" the found mindistance so it can be compared to our "z"-values*/
1127                        }
1128                }
1129        }
1130
1131        return LW_TRUE;
1132}
1133
1134
1135/**
1136        This is the same function as lw_dist2d_seg_seg but
1137        without any calculations to determine intersection since we
1138        already know they do not intersect
1139*/
1140int
1141lw_dist2d_selected_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D, DISTPTS *dl)
1142{
1143        LWDEBUGF(2, "lw_dist2d_selected_seg_seg [%g,%g]->[%g,%g] by [%g,%g]->[%g,%g]",
1144                 A->x,A->y,B->x,B->y, C->x,C->y, D->x, D->y);
1145
1146        /*A and B are the same point */
1147        if (  ( A->x == B->x) && (A->y == B->y) )
1148        {
1149                return lw_dist2d_pt_seg(A,C,D,dl);
1150        }
1151        /*U and V are the same point */
1152
1153        if (  ( C->x == D->x) && (C->y == D->y) )
1154        {
1155                dl->twisted= ((dl->twisted) * (-1));
1156                return lw_dist2d_pt_seg(D,A,B,dl);
1157        }
1158
1159        if ((lw_dist2d_pt_seg(A,C,D,dl)) && (lw_dist2d_pt_seg(B,C,D,dl)))
1160        {
1161                dl->twisted= ((dl->twisted) * (-1));  /*here we change the order of inputted geometrys and that we  notice by changing sign on dl->twisted*/
1162                return ((lw_dist2d_pt_seg(C,A,B,dl)) && (lw_dist2d_pt_seg(D,A,B,dl))); /*if all is successful we return true*/
1163        }
1164        else
1165        {
1166                return LW_FALSE; /* if any of the calls to lw_dist2d_pt_seg goes wrong we return false*/
1167        }
1168}
1169
1170/*------------------------------------------------------------------------------------------------------------
1171End of New faster distance calculations
1172--------------------------------------------------------------------------------------------------------------*/
1173
1174
1175/*------------------------------------------------------------------------------------------------------------
1176Functions in common for Brute force and new calculation
1177--------------------------------------------------------------------------------------------------------------*/
1178
1179/**
1180 * pt_in_ring_2d(): crossing number test for a point in a polygon
1181 *      input:   p = a point,
1182 *               pa = vertex points of a ring V[n+1] with V[n]=V[0]
1183 *      returns: 0 = outside, 1 = inside
1184 *
1185 *      Our polygons have first and last point the same,
1186 *
1187 */
1188int
1189pt_in_ring_2d(const POINT2D *p, const POINTARRAY *ring)
1190{
1191        int cn = 0;    /* the crossing number counter */
1192        int i;
1193        POINT2D v1, v2;
1194
1195#if INTEGRITY_CHECKS
1196        POINT2D first, last;
1197
1198        getPoint2d_p(ring, 0, &first);
1199        getPoint2d_p(ring, ring->npoints-1, &last);
1200        if ( memcmp(&first, &last, sizeof(POINT2D)) )
1201        {
1202                lwerror("pt_in_ring_2d: V[n] != V[0] (%g %g != %g %g)",
1203                        first.x, first.y, last.x, last.y);
1204                return LW_FALSE;
1205
1206        }
1207#endif
1208
1209        LWDEBUGF(2, "pt_in_ring_2d called with point: %g %g", p->x, p->y);
1210        /* printPA(ring); */
1211
1212        /* loop through all edges of the polygon */
1213        getPoint2d_p(ring, 0, &v1);
1214        for (i=0; i<ring->npoints-1; i++)
1215        {
1216                double vt;
1217                getPoint2d_p(ring, i+1, &v2);
1218
1219                /* edge from vertex i to vertex i+1 */
1220                if
1221                (
1222                    /* an upward crossing */
1223                    ((v1.y <= p->y) && (v2.y > p->y))
1224                    /* a downward crossing */
1225                    || ((v1.y > p->y) && (v2.y <= p->y))
1226                )
1227                {
1228
1229                        vt = (double)(p->y - v1.y) / (v2.y - v1.y);
1230
1231                        /* P.x <intersect */
1232                        if (p->x < v1.x + vt * (v2.x - v1.x))
1233                        {
1234                                /* a valid crossing of y=p.y right of p.x */
1235                                ++cn;
1236                        }
1237                }
1238                v1 = v2;
1239        }
1240
1241        LWDEBUGF(3, "pt_in_ring_2d returning %d", cn&1);
1242
1243        return (cn&1);    /* 0 if even (out), and 1 if odd (in) */
1244}
1245
1246
1247/**
1248
1249lw_dist2d_comp from p to line A->B
1250This one is now sending every occation to lw_dist2d_pt_pt
1251Before it was handling occations where r was between 0 and 1 internally
1252and just returning the distance without identifying the points.
1253To get this points it was nessecary to change and it also showed to be about 10%faster.
1254*/
1255int
1256lw_dist2d_pt_seg(POINT2D *p, POINT2D *A, POINT2D *B, DISTPTS *dl)
1257{
1258        POINT2D c;
1259        double  r;
1260        /*if start==end, then use pt distance */
1261        if (  ( A->x == B->x) && (A->y == B->y) )
1262        {
1263                return lw_dist2d_pt_pt(p,A,dl);
1264        }
1265        /*
1266         * otherwise, we use comp.graphics.algorithms
1267         * Frequently Asked Questions method
1268         *
1269         *  (1)               AC dot AB
1270             *         r = ---------
1271             *               ||AB||^2
1272         *      r has the following meaning:
1273         *      r=0 P = A
1274         *      r=1 P = B
1275         *      r<0 P is on the backward extension of AB
1276         *      r>1 P is on the forward extension of AB
1277         *      0<r<1 P is interior to AB
1278         */
1279
1280        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) );
1281
1282        /*This is for finding the maxdistance.
1283        the maxdistance have to be between two vertexes,
1284        compared to mindistance which can be between
1285        tvo vertexes vertex.*/
1286        if (dl->mode == DIST2D_MAX)
1287        {
1288                if (r>=0.5)
1289                {
1290                        return lw_dist2d_pt_pt(p,A,dl);
1291                }
1292                if (r<0.5)
1293                {
1294                        return lw_dist2d_pt_pt(p,B,dl);
1295                }
1296        }
1297
1298        if (r<0)        /*If the first vertex A is closest to the point p*/
1299        {
1300                return lw_dist2d_pt_pt(p,A,dl);
1301        }
1302        if (r>1)        /*If the second vertex B is closest to the point p*/
1303        {
1304                return lw_dist2d_pt_pt(p,B,dl);
1305        }
1306
1307        /*else if the point p is closer to some point between a and b
1308        then we find that point and send it to lw_dist2d_pt_pt*/
1309        c.x=A->x + r * (B->x-A->x);
1310        c.y=A->y + r * (B->y-A->y);
1311
1312        return lw_dist2d_pt_pt(p,&c,dl);
1313}
1314
1315
1316/**
1317
1318Compares incomming points and
1319stores the points closest to each other
1320or most far away from each other
1321depending on dl->mode (max or min)
1322*/
1323int
1324lw_dist2d_pt_pt(POINT2D *thep1, POINT2D *thep2,DISTPTS *dl)
1325{
1326        double hside = thep2->x - thep1->x;
1327        double vside = thep2->y - thep1->y;
1328        double dist = sqrt ( hside*hside + vside*vside );
1329
1330        if (((dl->distance - dist)*(dl->mode))>0) /*multiplication with mode to handle mindistance (mode=1)  and maxdistance (mode = (-1)*/
1331        {
1332                dl->distance = dist;
1333
1334                if (dl->twisted>0)      /*To get the points in right order. twisted is updated between 1 and (-1) every time the order is changed earlier in the chain*/
1335                {
1336                        dl->p1 = *thep1;
1337                        dl->p2 = *thep2;
1338                }
1339                else
1340                {
1341                        dl->p1 = *thep2;
1342                        dl->p2 = *thep1;
1343                }
1344        }
1345        return LW_TRUE;
1346}
1347
1348
1349
1350
1351/*------------------------------------------------------------------------------------------------------------
1352End of Functions in common for Brute force and new calculation
1353--------------------------------------------------------------------------------------------------------------*/
1354
1355
1356/*Mixed functions*/
1357
1358
1359/**
1360
1361 true if point is in poly (and not in its holes)
1362 It's not used by postgis but since I don't know what else
1363 can be affectes in the world I don't dare removing it.
1364 */
1365int
1366pt_in_poly_2d(const POINT2D *p, const LWPOLY *poly)
1367{
1368        int i;
1369
1370        /* Not in outer ring */
1371        if ( ! pt_in_ring_2d(p, poly->rings[0]) ) return 0;
1372
1373        /* Check holes */
1374        for (i=1; i<poly->nrings; i++)
1375        {
1376                /* Inside a hole */
1377                if ( pt_in_ring_2d(p, poly->rings[i]) ) return 0;
1378        }
1379
1380        return 1; /* In outer ring, not in holes */
1381}
1382
1383
1384/**
1385The old function nessecary for ptarray_segmentize2d in ptarray.c
1386*/
1387double
1388distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2)
1389{
1390        double hside = p2->x - p1->x;
1391        double vside = p2->y - p1->y;
1392
1393        return sqrt ( hside*hside + vside*vside );
1394
1395        /* the above is more readable
1396           return sqrt(
1397                (p2->x-p1->x) * (p2->x-p1->x) + (p2->y-p1->y) * (p2->y-p1->y)
1398                );  */
1399}
1400
1401
1402
1403/**
1404
1405The old function nessecary for ptarray_segmentize2d in ptarray.c
1406*/
1407double
1408distance2d_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B)
1409{
1410        double  r,s;
1411
1412        /*if start==end, then use pt distance */
1413        if (  ( A->x == B->x) && (A->y == B->y) )
1414                return distance2d_pt_pt(p,A);
1415
1416        /*
1417         * otherwise, we use comp.graphics.algorithms
1418         * Frequently Asked Questions method
1419         *
1420         *  (1)               AC dot AB
1421                *         r = ---------
1422                *               ||AB||^2
1423         *      r has the following meaning:
1424         *      r=0 P = A
1425         *      r=1 P = B
1426         *      r<0 P is on the backward extension of AB
1427         *      r>1 P is on the forward extension of AB
1428         *      0<r<1 P is interior to AB
1429         */
1430
1431        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) );
1432
1433        if (r<0) return distance2d_pt_pt(p,A);
1434        if (r>1) return distance2d_pt_pt(p,B);
1435
1436
1437        /*
1438         * (2)
1439         *           (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
1440         *      s = -----------------------------
1441         *                      L^2
1442         *
1443         *      Then the distance from C to P = |s|*L.
1444         *
1445         */
1446
1447        s = ( (A->y-p->y)*(B->x-A->x)- (A->x-p->x)*(B->y-A->y) ) /
1448            ( (B->x-A->x)*(B->x-A->x) +(B->y-A->y)*(B->y-A->y) );
1449
1450        return LW_ABS(s) * sqrt(
1451                   (B->x-A->x)*(B->x-A->x) + (B->y-A->y)*(B->y-A->y)
1452               );
1453}
1454
1455
1456
1457
1458
1459
1460
1461/**
1462find the 2d length of the given POINTARRAY (even if it's 3d)
1463*/
1464double
1465lwgeom_pointarray_length2d(const POINTARRAY *pts)
1466{
1467        double dist = 0.0;
1468        int i;
1469        POINT2D frm;
1470        POINT2D to;
1471
1472        if ( pts->npoints < 2 ) return 0.0;
1473        for (i=0; i<pts->npoints-1; i++)
1474        {
1475                getPoint2d_p(pts, i, &frm);
1476                getPoint2d_p(pts, i+1, &to);
1477                dist += sqrt( ( (frm.x - to.x)*(frm.x - to.x) )  +
1478                              ((frm.y - to.y)*(frm.y - to.y) ) );
1479        }
1480        return dist;
1481}
1482
1483/**
1484 * Find the 3d/2d length of the given POINTARRAY
1485 * (depending on its dimensions)
1486 */
1487double
1488lwgeom_pointarray_length(const POINTARRAY *pts)
1489{
1490        double dist = 0.0;
1491        int i;
1492        POINT3DZ frm;
1493        POINT3DZ to;
1494
1495        if ( pts->npoints < 2 ) return 0.0;
1496
1497        /* compute 2d length if 3d is not available */
1498        if ( ! TYPE_HASZ(pts->dims) ) return lwgeom_pointarray_length2d(pts);
1499
1500        for (i=0; i<pts->npoints-1; i++)
1501        {
1502                getPoint3dz_p(pts, i, &frm);
1503                getPoint3dz_p(pts, i+1, &to);
1504                dist += sqrt( ( (frm.x - to.x)*(frm.x - to.x) )  +
1505                              ((frm.y - to.y)*(frm.y - to.y) ) +
1506                              ((frm.z - to.z)*(frm.z - to.z) ) );
1507        }
1508        return dist;
1509}
1510
1511/**
1512 * This should be rewritten to make use of the curve itself.
1513 */
1514double
1515lwgeom_curvepolygon_area(LWCURVEPOLY *curvepoly)
1516{
1517        LWPOLY *poly = (LWPOLY *)lwgeom_segmentize((LWGEOM *)curvepoly, 32);
1518        return lwgeom_polygon_area(poly);
1519}
1520
1521/**
1522 * Find the area of the outer ring - sum (area of inner rings).
1523 * Could use a more numerically stable calculator...
1524 */
1525double
1526lwgeom_polygon_area(const LWPOLY *poly)
1527{
1528        double poly_area=0.0;
1529        int i;
1530        POINT2D p1;
1531        POINT2D p2;
1532
1533        LWDEBUGF(2, "in lwgeom_polygon_area (%d rings)", poly->nrings);
1534
1535        for (i=0; i<poly->nrings; i++)
1536        {
1537                int j;
1538                POINTARRAY *ring = poly->rings[i];
1539                double ringarea = 0.0;
1540
1541                LWDEBUGF(4, " rings %d has %d points", i, ring->npoints);
1542
1543                if ( ! ring->npoints ) continue; /* empty ring */
1544                for (j=0; j<ring->npoints-1; j++)
1545                {
1546                        getPoint2d_p(ring, j, &p1);
1547                        getPoint2d_p(ring, j+1, &p2);
1548                        ringarea += ( p1.x * p2.y ) - ( p1.y * p2.x );
1549                }
1550
1551                ringarea  /= 2.0;
1552
1553                LWDEBUGF(4, " ring 1 has area %lf",ringarea);
1554
1555                ringarea  = fabs(ringarea);
1556                if (i != 0)     /*outer */
1557                        ringarea  = -1.0*ringarea ; /* its a hole */
1558
1559                poly_area += ringarea;
1560        }
1561
1562        return poly_area;
1563}
1564
1565/**
1566 * Compute the sum of polygon rings length.
1567 * Could use a more numerically stable calculator...
1568 */
1569double
1570lwgeom_polygon_perimeter(const LWPOLY *poly)
1571{
1572        double result=0.0;
1573        int i;
1574
1575        LWDEBUGF(2, "in lwgeom_polygon_perimeter (%d rings)", poly->nrings);
1576
1577        for (i=0; i<poly->nrings; i++)
1578                result += lwgeom_pointarray_length(poly->rings[i]);
1579
1580        return result;
1581}
1582
1583/**
1584 * Compute the sum of polygon rings length (forcing 2d computation).
1585 * Could use a more numerically stable calculator...
1586 */
1587double
1588lwgeom_polygon_perimeter2d(const LWPOLY *poly)
1589{
1590        double result=0.0;
1591        int i;
1592
1593        LWDEBUGF(2, "in lwgeom_polygon_perimeter (%d rings)", poly->nrings);
1594
1595        for (i=0; i<poly->nrings; i++)
1596                result += lwgeom_pointarray_length2d(poly->rings[i]);
1597
1598        return result;
1599}
1600
1601
1602int
1603lwgeom_pt_inside_circle(POINT2D *p, double cx, double cy, double rad)
1604{
1605        POINT2D center;
1606
1607        center.x = cx;
1608        center.y = cy;
1609
1610        if ( distance2d_pt_pt(p, &center) < rad ) return 1;
1611        else return 0;
1612
1613}
1614
1615/**
1616 * Compute the azimuth of segment AB in radians.
1617 * Return 0 on exception (same point), 1 otherwise.
1618 */
1619int
1620azimuth_pt_pt(const POINT2D *A, const POINT2D *B, double *d)
1621{
1622        if ( A->x == B->x )
1623        {
1624                if ( A->y < B->y ) *d=0.0;
1625                else if ( A->y > B->y ) *d=M_PI;
1626                else return 0;
1627                return 1;
1628        }
1629
1630        if ( A->y == B->y )
1631        {
1632                if ( A->x < B->x ) *d=M_PI/2;
1633                else if ( A->x > B->x ) *d=M_PI+(M_PI/2);
1634                else return 0;
1635                return 1;
1636        }
1637
1638        if ( A->x < B->x )
1639        {
1640                if ( A->y < B->y )
1641                {
1642                        *d=atan(fabs(A->x - B->x) / fabs(A->y - B->y) );
1643                }
1644                else /* ( A->y > B->y )  - equality case handled above */
1645                {
1646                        *d=atan(fabs(A->y - B->y) / fabs(A->x - B->x) )
1647                           + (M_PI/2);
1648                }
1649        }
1650
1651        else /* ( A->x > B->x ) - equality case handled above */
1652        {
1653                if ( A->y > B->y )
1654                {
1655                        *d=atan(fabs(A->x - B->x) / fabs(A->y - B->y) )
1656                           + M_PI;
1657                }
1658                else /* ( A->y < B->y )  - equality case handled above */
1659                {
1660                        *d=atan(fabs(A->y - B->y) / fabs(A->x - B->x) )
1661                           + (M_PI+(M_PI/2));
1662                }
1663        }
1664
1665        return 1;
1666}
1667
1668
Note: See TracBrowser for help on using the repository browser.