PostGIS  2.2.7dev-r@@SVN_REVISION@@
measures.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  * Copyright 2001-2006 Refractions Research Inc.
6  * Copyright 2010 Nicklas Avén
7  * Copyright 2012 Paul Ramsey
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 <string.h>
15 #include <stdlib.h>
16 
17 #include "measures.h"
18 #include "lwgeom_log.h"
19 
20 
21 
22 /*------------------------------------------------------------------------------------------------------------
23 Initializing functions
24 The functions starting the distance-calculation processses
25 --------------------------------------------------------------------------------------------------------------*/
26 
27 LWGEOM *
28 lwgeom_closest_line(const LWGEOM *lw1, const LWGEOM *lw2)
29 {
30  return lw_dist2d_distanceline(lw1, lw2, lw1->srid, DIST_MIN);
31 }
32 
33 LWGEOM *
34 lwgeom_furthest_line(const LWGEOM *lw1, const LWGEOM *lw2)
35 {
36  return lw_dist2d_distanceline(lw1, lw2, lw1->srid, DIST_MAX);
37 }
38 
39 LWGEOM *
40 lwgeom_closest_point(const LWGEOM *lw1, const LWGEOM *lw2)
41 {
42  return lw_dist2d_distancepoint(lw1, lw2, lw1->srid, DIST_MIN);
43 }
44 
45 LWGEOM *
46 lwgeom_furthest_point(const LWGEOM *lw1, const LWGEOM *lw2)
47 {
48  return lw_dist2d_distancepoint(lw1, lw2, lw1->srid, DIST_MAX);
49 }
50 
51 
52 void
54 {
55  dl->twisted = -1;
56  dl->p1.x = dl->p1.y = 0.0;
57  dl->p2.x = dl->p2.y = 0.0;
58  dl->mode = mode;
59  dl->tolerance = 0.0;
60  if ( mode == DIST_MIN )
61  dl->distance = FLT_MAX;
62  else
63  dl->distance = -1 * FLT_MAX;
64 }
65 
69 LWGEOM *
70 lw_dist2d_distanceline(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
71 {
72  double x1,x2,y1,y2;
73 
74  double initdistance = ( mode == DIST_MIN ? FLT_MAX : -1.0);
75  DISTPTS thedl;
76  LWPOINT *lwpoints[2];
77  LWGEOM *result;
78 
79  thedl.mode = mode;
80  thedl.distance = initdistance;
81  thedl.tolerance = 0.0;
82 
83  LWDEBUG(2, "lw_dist2d_distanceline is called");
84 
85  if (!lw_dist2d_comp( lw1,lw2,&thedl))
86  {
87  /*should never get here. all cases ought to be error handled earlier*/
88  lwerror("Some unspecified error.");
89  result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
90  }
91 
92  /*if thedl.distance is unchanged there where only empty geometries input*/
93  if (thedl.distance == initdistance)
94  {
95  LWDEBUG(3, "didn't find geometries to measure between, returning null");
96  result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
97  }
98  else
99  {
100  x1=thedl.p1.x;
101  y1=thedl.p1.y;
102  x2=thedl.p2.x;
103  y2=thedl.p2.y;
104 
105  lwpoints[0] = lwpoint_make2d(srid, x1, y1);
106  lwpoints[1] = lwpoint_make2d(srid, x2, y2);
107 
108  result = (LWGEOM *)lwline_from_ptarray(srid, 2, lwpoints);
109  }
110  return result;
111 }
112 
116 LWGEOM *
117 lw_dist2d_distancepoint(const LWGEOM *lw1, const LWGEOM *lw2,int srid,int mode)
118 {
119  double x,y;
120  DISTPTS thedl;
121  double initdistance = FLT_MAX;
122  LWGEOM *result;
123 
124  thedl.mode = mode;
125  thedl.distance= initdistance;
126  thedl.tolerance = 0;
127 
128  LWDEBUG(2, "lw_dist2d_distancepoint is called");
129 
130  if (!lw_dist2d_comp( lw1,lw2,&thedl))
131  {
132  /*should never get here. all cases ought to be error handled earlier*/
133  lwerror("Some unspecified error.");
134  result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
135  }
136  if (thedl.distance == initdistance)
137  {
138  LWDEBUG(3, "didn't find geometries to measure between, returning null");
139  result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
140  }
141  else
142  {
143  x=thedl.p1.x;
144  y=thedl.p1.y;
145  result = (LWGEOM *)lwpoint_make2d(srid, x, y);
146  }
147  return result;
148 }
149 
150 
154 double
155 lwgeom_maxdistance2d(const LWGEOM *lw1, const LWGEOM *lw2)
156 {
157  LWDEBUG(2, "lwgeom_maxdistance2d is called");
158 
159  return lwgeom_maxdistance2d_tolerance( lw1, lw2, 0.0 );
160 }
161 
166 double
167 lwgeom_maxdistance2d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
168 {
169  /*double thedist;*/
170  DISTPTS thedl;
171  LWDEBUG(2, "lwgeom_maxdistance2d_tolerance is called");
172  thedl.mode = DIST_MAX;
173  thedl.distance= -1;
174  thedl.tolerance = tolerance;
175  if (lw_dist2d_comp( lw1,lw2,&thedl))
176  {
177  return thedl.distance;
178  }
179  /*should never get here. all cases ought to be error handled earlier*/
180  lwerror("Some unspecified error.");
181  return -1;
182 }
183 
187 double
188 lwgeom_mindistance2d(const LWGEOM *lw1, const LWGEOM *lw2)
189 {
190  LWDEBUG(2, "lwgeom_mindistance2d is called");
191  return lwgeom_mindistance2d_tolerance( lw1, lw2, 0.0 );
192 }
193 
198 double
199 lwgeom_mindistance2d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
200 {
201  DISTPTS thedl;
202  LWDEBUG(2, "lwgeom_mindistance2d_tolerance is called");
203  thedl.mode = DIST_MIN;
204  thedl.distance= FLT_MAX;
205  thedl.tolerance = tolerance;
206  if (lw_dist2d_comp( lw1,lw2,&thedl))
207  {
208  return thedl.distance;
209  }
210  /*should never get here. all cases ought to be error handled earlier*/
211  lwerror("Some unspecified error.");
212  return FLT_MAX;
213 }
214 
215 
216 /*------------------------------------------------------------------------------------------------------------
217 End of Initializing functions
218 --------------------------------------------------------------------------------------------------------------*/
219 
220 /*------------------------------------------------------------------------------------------------------------
221 Preprocessing functions
222 Functions preparing geometries for distance-calculations
223 --------------------------------------------------------------------------------------------------------------*/
224 
230 int
231 lw_dist2d_comp(const LWGEOM *lw1,const LWGEOM *lw2, DISTPTS *dl)
232 {
233  LWDEBUG(2, "lw_dist2d_comp is called");
234 
235  return lw_dist2d_recursive(lw1, lw2, dl);
236 }
237 
238 static int
240 {
241 
242  switch (g->type)
243  {
244  case MULTIPOINTTYPE:
245  case MULTILINETYPE:
246  case MULTIPOLYGONTYPE:
247  case COLLECTIONTYPE:
248  case MULTICURVETYPE:
249  case MULTISURFACETYPE:
250  case COMPOUNDTYPE:
252  return LW_TRUE;
253  break;
254 
255  default:
256  return LW_FALSE;
257  }
258 }
259 
263 int lw_dist2d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS *dl)
264 {
265  int i, j;
266  int n1=1;
267  int n2=1;
268  LWGEOM *g1 = NULL;
269  LWGEOM *g2 = NULL;
270  LWCOLLECTION *c1 = NULL;
271  LWCOLLECTION *c2 = NULL;
272 
273  LWDEBUGF(2, "lw_dist2d_comp is called with type1=%d, type2=%d", lwg1->type, lwg2->type);
274 
275  if (lw_dist2d_is_collection(lwg1))
276  {
277  LWDEBUG(3, "First geometry is collection");
278  c1 = lwgeom_as_lwcollection(lwg1);
279  n1 = c1->ngeoms;
280  }
281  if (lw_dist2d_is_collection(lwg2))
282  {
283  LWDEBUG(3, "Second geometry is collection");
284  c2 = lwgeom_as_lwcollection(lwg2);
285  n2 = c2->ngeoms;
286  }
287 
288  for ( i = 0; i < n1; i++ )
289  {
290 
291  if (lw_dist2d_is_collection(lwg1))
292  {
293  g1 = c1->geoms[i];
294  }
295  else
296  {
297  g1 = (LWGEOM*)lwg1;
298  }
299 
300  if (lwgeom_is_empty(g1)) return LW_TRUE;
301 
302  if (lw_dist2d_is_collection(g1))
303  {
304  LWDEBUG(3, "Found collection inside first geometry collection, recursing");
305  if (!lw_dist2d_recursive(g1, lwg2, dl)) return LW_FALSE;
306  continue;
307  }
308  for ( j = 0; j < n2; j++ )
309  {
310  if (lw_dist2d_is_collection(lwg2))
311  {
312  g2 = c2->geoms[j];
313  }
314  else
315  {
316  g2 = (LWGEOM*)lwg2;
317  }
318  if (lw_dist2d_is_collection(g2))
319  {
320  LWDEBUG(3, "Found collection inside second geometry collection, recursing");
321  if (!lw_dist2d_recursive(g1, g2, dl)) return LW_FALSE;
322  continue;
323  }
324 
325  if ( ! g1->bbox )
326  {
327  lwgeom_add_bbox(g1);
328  }
329  if ( ! g2->bbox )
330  {
331  lwgeom_add_bbox(g2);
332  }
333 
334  /*If one of geometries is empty, return. True here only means continue searching. False would have stoped the process*/
335  if (lwgeom_is_empty(g1)||lwgeom_is_empty(g2)) return LW_TRUE;
336 
337  if ( (dl->mode != DIST_MAX) &&
338  (! lw_dist2d_check_overlap(g1, g2)) &&
339  (g1->type == LINETYPE || g1->type == POLYGONTYPE) &&
340  (g2->type == LINETYPE || g2->type == POLYGONTYPE) )
341  {
342  if (!lw_dist2d_distribute_fast(g1, g2, dl)) return LW_FALSE;
343  }
344  else
345  {
346  if (!lw_dist2d_distribute_bruteforce(g1, g2, dl)) return LW_FALSE;
347  if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
348  }
349  }
350  }
351  return LW_TRUE;
352 }
353 
354 
355 int
357 {
358 
359  int t1 = lwg1->type;
360  int t2 = lwg2->type;
361 
362  switch ( t1 )
363  {
364  case POINTTYPE:
365  {
366  dl->twisted = 1;
367  switch ( t2 )
368  {
369  case POINTTYPE:
370  return lw_dist2d_point_point((LWPOINT *)lwg1, (LWPOINT *)lwg2, dl);
371  case LINETYPE:
372  return lw_dist2d_point_line((LWPOINT *)lwg1, (LWLINE *)lwg2, dl);
373  case POLYGONTYPE:
374  return lw_dist2d_point_poly((LWPOINT *)lwg1, (LWPOLY *)lwg2, dl);
375  case CIRCSTRINGTYPE:
376  return lw_dist2d_point_circstring((LWPOINT *)lwg1, (LWCIRCSTRING *)lwg2, dl);
377  case CURVEPOLYTYPE:
378  return lw_dist2d_point_curvepoly((LWPOINT *)lwg1, (LWCURVEPOLY *)lwg2, dl);
379  default:
380  lwerror("Unsupported geometry type: %s", lwtype_name(t2));
381  }
382  }
383  case LINETYPE:
384  {
385  dl->twisted = 1;
386  switch ( t2 )
387  {
388  case POINTTYPE:
389  dl->twisted=(-1);
390  return lw_dist2d_point_line((LWPOINT *)lwg2, (LWLINE *)lwg1, dl);
391  case LINETYPE:
392  return lw_dist2d_line_line((LWLINE *)lwg1, (LWLINE *)lwg2, dl);
393  case POLYGONTYPE:
394  return lw_dist2d_line_poly((LWLINE *)lwg1, (LWPOLY *)lwg2, dl);
395  case CIRCSTRINGTYPE:
396  return lw_dist2d_line_circstring((LWLINE *)lwg1, (LWCIRCSTRING *)lwg2, dl);
397  case CURVEPOLYTYPE:
398  return lw_dist2d_line_curvepoly((LWLINE *)lwg1, (LWCURVEPOLY *)lwg2, dl);
399  default:
400  lwerror("Unsupported geometry type: %s", lwtype_name(t2));
401  }
402  }
403  case CIRCSTRINGTYPE:
404  {
405  dl->twisted = 1;
406  switch ( t2 )
407  {
408  case POINTTYPE:
409  dl->twisted = -1;
410  return lw_dist2d_point_circstring((LWPOINT *)lwg2, (LWCIRCSTRING *)lwg1, dl);
411  case LINETYPE:
412  dl->twisted = -1;
413  return lw_dist2d_line_circstring((LWLINE *)lwg2, (LWCIRCSTRING *)lwg1, dl);
414  case POLYGONTYPE:
415  return lw_dist2d_circstring_poly((LWCIRCSTRING *)lwg1, (LWPOLY *)lwg2, dl);
416  case CIRCSTRINGTYPE:
417  return lw_dist2d_circstring_circstring((LWCIRCSTRING *)lwg1, (LWCIRCSTRING *)lwg2, dl);
418  case CURVEPOLYTYPE:
419  return lw_dist2d_circstring_curvepoly((LWCIRCSTRING *)lwg1, (LWCURVEPOLY *)lwg2, dl);
420  default:
421  lwerror("Unsupported geometry type: %s", lwtype_name(t2));
422  }
423  }
424  case POLYGONTYPE:
425  {
426  dl->twisted = -1;
427  switch ( t2 )
428  {
429  case POINTTYPE:
430  return lw_dist2d_point_poly((LWPOINT *)lwg2, (LWPOLY *)lwg1, dl);
431  case LINETYPE:
432  return lw_dist2d_line_poly((LWLINE *)lwg2, (LWPOLY *)lwg1, dl);
433  case CIRCSTRINGTYPE:
434  return lw_dist2d_circstring_poly((LWCIRCSTRING *)lwg2, (LWPOLY *)lwg1, dl);
435  case POLYGONTYPE:
436  dl->twisted = 1;
437  return lw_dist2d_poly_poly((LWPOLY *)lwg1, (LWPOLY *)lwg2, dl);
438  case CURVEPOLYTYPE:
439  dl->twisted = 1;
440  return lw_dist2d_poly_curvepoly((LWPOLY *)lwg1, (LWCURVEPOLY *)lwg2, dl);
441  default:
442  lwerror("Unsupported geometry type: %s", lwtype_name(t2));
443  }
444  }
445  case CURVEPOLYTYPE:
446  {
447  dl->twisted = (-1);
448  switch ( t2 )
449  {
450  case POINTTYPE:
451  return lw_dist2d_point_curvepoly((LWPOINT *)lwg2, (LWCURVEPOLY *)lwg1, dl);
452  case LINETYPE:
453  return lw_dist2d_line_curvepoly((LWLINE *)lwg2, (LWCURVEPOLY *)lwg1, dl);
454  case POLYGONTYPE:
455  return lw_dist2d_poly_curvepoly((LWPOLY *)lwg2, (LWCURVEPOLY *)lwg1, dl);
456  case CIRCSTRINGTYPE:
457  return lw_dist2d_circstring_curvepoly((LWCIRCSTRING *)lwg2, (LWCURVEPOLY *)lwg1, dl);
458  case CURVEPOLYTYPE:
459  dl->twisted = 1;
460  return lw_dist2d_curvepoly_curvepoly((LWCURVEPOLY *)lwg1, (LWCURVEPOLY *)lwg2, dl);
461  default:
462  lwerror("Unsupported geometry type: %s", lwtype_name(t2));
463  }
464  }
465  default:
466  {
467  lwerror("Unsupported geometry type: %s", lwtype_name(t1));
468  }
469  }
470 
471  /*You shouldn't being able to get here*/
472  lwerror("unspecified error in function lw_dist2d_distribute_bruteforce");
473  return LW_FALSE;
474 }
475 
476 
477 
478 
483 int
485 {
486  LWDEBUG(2, "lw_dist2d_check_overlap is called");
487  if ( ! lwg1->bbox )
488  lwgeom_calculate_gbox(lwg1, lwg1->bbox);
489  if ( ! lwg2->bbox )
490  lwgeom_calculate_gbox(lwg2, lwg2->bbox);
491 
492  /*Check if the geometries intersect.
493  */
494  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))
495  {
496  LWDEBUG(3, "geometries bboxes did not overlap");
497  return LW_FALSE;
498  }
499  LWDEBUG(3, "geometries bboxes overlap");
500  return LW_TRUE;
501 }
502 
507 int
509 {
510  POINTARRAY *pa1, *pa2;
511  int type1 = lwg1->type;
512  int type2 = lwg2->type;
513 
514  LWDEBUGF(2, "lw_dist2d_distribute_fast is called with typ1=%d, type2=%d", lwg1->type, lwg2->type);
515 
516  switch (type1)
517  {
518  case LINETYPE:
519  pa1 = ((LWLINE *)lwg1)->points;
520  break;
521  case POLYGONTYPE:
522  pa1 = ((LWPOLY *)lwg1)->rings[0];
523  break;
524  default:
525  lwerror("Unsupported geometry1 type: %s", lwtype_name(type1));
526  return LW_FALSE;
527  }
528  switch (type2)
529  {
530  case LINETYPE:
531  pa2 = ((LWLINE *)lwg2)->points;
532  break;
533  case POLYGONTYPE:
534  pa2 = ((LWPOLY *)lwg2)->rings[0];
535  break;
536  default:
537  lwerror("Unsupported geometry2 type: %s", lwtype_name(type1));
538  return LW_FALSE;
539  }
540  dl->twisted=1;
541  return lw_dist2d_fast_ptarray_ptarray(pa1, pa2, dl, lwg1->bbox, lwg2->bbox);
542 }
543 
544 /*------------------------------------------------------------------------------------------------------------
545 End of Preprocessing functions
546 --------------------------------------------------------------------------------------------------------------*/
547 
548 
549 /*------------------------------------------------------------------------------------------------------------
550 Brute force functions
551 The old way of calculating distances, now used for:
552 1) distances to points (because there shouldn't be anything to gain by the new way of doing it)
553 2) distances when subgeometries geometries bboxes overlaps
554 --------------------------------------------------------------------------------------------------------------*/
555 
560 int
562 {
563  const POINT2D *p1, *p2;
564 
565  p1 = getPoint2d_cp(point1->point, 0);
566  p2 = getPoint2d_cp(point2->point, 0);
567 
568  return lw_dist2d_pt_pt(p1, p2, dl);
569 }
574 int
576 {
577  const POINT2D *p;
578  LWDEBUG(2, "lw_dist2d_point_line is called");
579  p = getPoint2d_cp(point->point, 0);
580  return lw_dist2d_pt_ptarray(p, line->points, dl);
581 }
582 
583 int
585 {
586  const POINT2D *p;
587  p = getPoint2d_cp(point->point, 0);
588  return lw_dist2d_pt_ptarrayarc(p, circ->points, dl);
589 }
590 
596 int
598 {
599  const POINT2D *p;
600  int i;
601 
602  LWDEBUG(2, "lw_dist2d_point_poly called");
603 
604  p = getPoint2d_cp(point->point, 0);
605 
606  if (dl->mode == DIST_MAX)
607  {
608  LWDEBUG(3, "looking for maxdistance");
609  return lw_dist2d_pt_ptarray(p, poly->rings[0], dl);
610  }
611  /* Return distance to outer ring if not inside it */
612  if ( ptarray_contains_point(poly->rings[0], p) == LW_OUTSIDE )
613  {
614  LWDEBUG(3, "first point not inside outer-ring");
615  return lw_dist2d_pt_ptarray(p, poly->rings[0], dl);
616  }
617 
618  /*
619  * Inside the outer ring.
620  * Scan though each of the inner rings looking to
621  * see if its inside. If not, distance==0.
622  * Otherwise, distance = pt to ring distance
623  */
624  for ( i = 1; i < poly->nrings; i++)
625  {
626  /* Inside a hole. Distance = pt -> ring */
627  if ( ptarray_contains_point(poly->rings[i], p) != LW_OUTSIDE )
628  {
629  LWDEBUG(3, " inside an hole");
630  return lw_dist2d_pt_ptarray(p, poly->rings[i], dl);
631  }
632  }
633 
634  LWDEBUG(3, " inside the polygon");
635  if (dl->mode == DIST_MIN)
636  {
637  dl->distance = 0.0;
638  dl->p1.x = dl->p2.x = p->x;
639  dl->p1.y = dl->p2.y = p->y;
640  }
641  return LW_TRUE; /* Is inside the polygon */
642 }
643 
644 int
646 {
647  const POINT2D *p;
648  int i;
649 
650  p = getPoint2d_cp(point->point, 0);
651 
652  if (dl->mode == DIST_MAX)
653  lwerror("lw_dist2d_point_curvepoly cannot calculate max distance");
654 
655  /* Return distance to outer ring if not inside it */
656  if ( lwgeom_contains_point(poly->rings[0], p) == LW_OUTSIDE )
657  {
658  return lw_dist2d_recursive((LWGEOM*)point, poly->rings[0], dl);
659  }
660 
661  /*
662  * Inside the outer ring.
663  * Scan though each of the inner rings looking to
664  * see if its inside. If not, distance==0.
665  * Otherwise, distance = pt to ring distance
666  */
667  for ( i = 1; i < poly->nrings; i++)
668  {
669  /* Inside a hole. Distance = pt -> ring */
670  if ( lwgeom_contains_point(poly->rings[i], p) != LW_OUTSIDE )
671  {
672  LWDEBUG(3, " inside a hole");
673  return lw_dist2d_recursive((LWGEOM*)point, poly->rings[i], dl);
674  }
675  }
676 
677  LWDEBUG(3, " inside the polygon");
678  if (dl->mode == DIST_MIN)
679  {
680  dl->distance = 0.0;
681  dl->p1.x = dl->p2.x = p->x;
682  dl->p1.y = dl->p2.y = p->y;
683  }
684 
685  return LW_TRUE; /* Is inside the polygon */
686 }
687 
692 int
694 {
695  POINTARRAY *pa1 = line1->points;
696  POINTARRAY *pa2 = line2->points;
697  LWDEBUG(2, "lw_dist2d_line_line is called");
698  return lw_dist2d_ptarray_ptarray(pa1, pa2, dl);
699 }
700 
701 int
703 {
704  return lw_dist2d_ptarray_ptarrayarc(line1->points, line2->points, dl);
705 }
706 
718 int
720 {
721  const POINT2D *pt;
722  int i;
723 
724  LWDEBUGF(2, "lw_dist2d_line_poly called (%d rings)", poly->nrings);
725 
726  pt = getPoint2d_cp(line->points, 0);
727  if ( ptarray_contains_point(poly->rings[0], pt) == LW_OUTSIDE )
728  {
729  return lw_dist2d_ptarray_ptarray(line->points, poly->rings[0], dl);
730  }
731 
732  for (i=1; i<poly->nrings; i++)
733  {
734  if (!lw_dist2d_ptarray_ptarray(line->points, poly->rings[i], dl)) return LW_FALSE;
735 
736  LWDEBUGF(3, " distance from ring %d: %f, mindist: %f",
737  i, dl->distance, dl->tolerance);
738  /* just a check if the answer is already given */
739  if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE;
740  }
741 
742  /*
743  * No intersection, have to check if a point is
744  * inside polygon
745  */
746  pt = getPoint2d_cp(line->points, 0);
747 
748  /*
749  * Outside outer ring, so min distance to a ring
750  * is the actual min distance
751 
752  if ( ! pt_in_ring_2d(&pt, poly->rings[0]) )
753  {
754  return ;
755  } */
756 
757  /*
758  * Its in the outer ring.
759  * Have to check if its inside a hole
760  */
761  for (i=1; i<poly->nrings; i++)
762  {
763  if ( ptarray_contains_point(poly->rings[i], pt) != LW_OUTSIDE )
764  {
765  /*
766  * Its inside a hole, then the actual
767  * distance is the min ring distance
768  */
769  return LW_TRUE;
770  }
771  }
772  if (dl->mode == DIST_MIN)
773  {
774  dl->distance = 0.0;
775  dl->p1.x = dl->p2.x = pt->x;
776  dl->p1.y = dl->p2.y = pt->y;
777  }
778  return LW_TRUE; /* Not in hole, so inside polygon */
779 }
780 
781 int
783 {
784  const POINT2D *pt = getPoint2d_cp(line->points, 0);
785  int i;
786 
787  if ( lwgeom_contains_point(poly->rings[0], pt) == LW_OUTSIDE )
788  {
789  return lw_dist2d_recursive((LWGEOM*)line, poly->rings[0], dl);
790  }
791 
792  for ( i = 1; i < poly->nrings; i++ )
793  {
794  if ( ! lw_dist2d_recursive((LWGEOM*)line, poly->rings[i], dl) )
795  return LW_FALSE;
796 
797  if ( dl->distance<=dl->tolerance && dl->mode == DIST_MIN )
798  return LW_TRUE;
799  }
800 
801  for ( i=1; i < poly->nrings; i++ )
802  {
803  if ( lwgeom_contains_point(poly->rings[i],pt) != LW_OUTSIDE )
804  {
805  /* Its inside a hole, then the actual */
806  return LW_TRUE;
807  }
808  }
809 
810  if (dl->mode == DIST_MIN)
811  {
812  dl->distance = 0.0;
813  dl->p1.x = dl->p2.x = pt->x;
814  dl->p1.y = dl->p2.y = pt->y;
815  }
816 
817  return LW_TRUE; /* Not in hole, so inside polygon */
818 }
819 
828 int
830 {
831 
832  const POINT2D *pt;
833  int i;
834 
835  LWDEBUG(2, "lw_dist2d_poly_poly called");
836 
837  /*1 if we are looking for maxdistance, just check the outer rings.*/
838  if (dl->mode == DIST_MAX)
839  {
840  return lw_dist2d_ptarray_ptarray(poly1->rings[0], poly2->rings[0], dl);
841  }
842 
843 
844  /* 2 check if poly1 has first point outside poly2 and vice versa, if so, just check outer rings
845  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.*/
846  pt = getPoint2d_cp(poly1->rings[0], 0);
847  if ( ptarray_contains_point(poly2->rings[0], pt) == LW_OUTSIDE )
848  {
849  pt = getPoint2d_cp(poly2->rings[0], 0);
850  if ( ptarray_contains_point(poly1->rings[0], pt) == LW_OUTSIDE )
851  {
852  return lw_dist2d_ptarray_ptarray(poly1->rings[0], poly2->rings[0], dl);
853  }
854  }
855 
856  /*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*/
857  pt = getPoint2d_cp(poly2->rings[0], 0);
858  for (i=1; i<poly1->nrings; i++)
859  {
860  /* Inside a hole */
861  if ( ptarray_contains_point(poly1->rings[i], pt) != LW_OUTSIDE )
862  {
863  return lw_dist2d_ptarray_ptarray(poly1->rings[i], poly2->rings[0], dl);
864  }
865  }
866 
867  /*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*/
868  pt = getPoint2d_cp(poly1->rings[0], 0);
869  for (i=1; i<poly2->nrings; i++)
870  {
871  /* Inside a hole */
872  if ( ptarray_contains_point(poly2->rings[i], pt) != LW_OUTSIDE )
873  {
874  return lw_dist2d_ptarray_ptarray(poly1->rings[0], poly2->rings[i], dl);
875  }
876  }
877 
878 
879  /*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.*/
880  pt = getPoint2d_cp(poly1->rings[0], 0);
881  if ( ptarray_contains_point(poly2->rings[0], pt) != LW_OUTSIDE )
882  {
883  dl->distance = 0.0;
884  dl->p1.x = dl->p2.x = pt->x;
885  dl->p1.y = dl->p2.y = pt->y;
886  return LW_TRUE;
887  }
888 
889  pt = getPoint2d_cp(poly2->rings[0], 0);
890  if ( ptarray_contains_point(poly1->rings[0], pt) != LW_OUTSIDE )
891  {
892  dl->distance = 0.0;
893  dl->p1.x = dl->p2.x = pt->x;
894  dl->p1.y = dl->p2.y = pt->y;
895  return LW_TRUE;
896  }
897 
898 
899  lwerror("Unspecified error in function lw_dist2d_poly_poly");
900  return LW_FALSE;
901 }
902 
903 int
905 {
906  LWCURVEPOLY *curvepoly1 = lwcurvepoly_construct_from_lwpoly(poly1);
907  int rv = lw_dist2d_curvepoly_curvepoly(curvepoly1, curvepoly2, dl);
908  lwgeom_free((LWGEOM*)curvepoly1);
909  return rv;
910 }
911 
912 int
914 {
916  int rv = lw_dist2d_line_curvepoly((LWLINE*)circ, curvepoly, dl);
917  lwgeom_free((LWGEOM*)curvepoly);
918  return rv;
919 }
920 
921 
922 int
924 {
925  return lw_dist2d_line_curvepoly((LWLINE*)circ, poly, dl);
926 }
927 
928 int
930 {
931  return lw_dist2d_ptarrayarc_ptarrayarc(line1->points, line2->points, dl);
932 }
933 
934 static const POINT2D *
936 {
937  switch( geom->type )
938  {
939  case LINETYPE:
940  return getPoint2d_cp(((LWLINE*)geom)->points, 0);
941  case CIRCSTRINGTYPE:
942  return getPoint2d_cp(((LWCIRCSTRING*)geom)->points, 0);
943  case COMPOUNDTYPE:
944  {
945  LWCOMPOUND *comp = (LWCOMPOUND*)geom;
946  LWLINE *line = (LWLINE*)(comp->geoms[0]);
947  return getPoint2d_cp(line->points, 0);
948  }
949  default:
950  lwerror("lw_curvering_getfirstpoint2d_cp: unknown type");
951  }
952  return NULL;
953 }
954 
955 int
957 {
958  const POINT2D *pt;
959  int i;
960 
961  LWDEBUG(2, "lw_dist2d_curvepoly_curvepoly called");
962 
963  /*1 if we are looking for maxdistance, just check the outer rings.*/
964  if (dl->mode == DIST_MAX)
965  {
966  return lw_dist2d_recursive(poly1->rings[0], poly2->rings[0], dl);
967  }
968 
969 
970  /* 2 check if poly1 has first point outside poly2 and vice versa, if so, just check outer rings
971  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.*/
972  pt = lw_curvering_getfirstpoint2d_cp(poly1->rings[0]);
973  if ( lwgeom_contains_point(poly2->rings[0], pt) == LW_OUTSIDE )
974  {
975  pt = lw_curvering_getfirstpoint2d_cp(poly2->rings[0]);
976  if ( lwgeom_contains_point(poly1->rings[0], pt) == LW_OUTSIDE )
977  {
978  return lw_dist2d_recursive(poly1->rings[0], poly2->rings[0], dl);
979  }
980  }
981 
982  /*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*/
983  pt = lw_curvering_getfirstpoint2d_cp(poly2->rings[0]);
984  for (i = 1; i < poly1->nrings; i++)
985  {
986  /* Inside a hole */
987  if ( lwgeom_contains_point(poly1->rings[i], pt) != LW_OUTSIDE )
988  {
989  return lw_dist2d_recursive(poly1->rings[i], poly2->rings[0], dl);
990  }
991  }
992 
993  /*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*/
994  pt = lw_curvering_getfirstpoint2d_cp(poly1->rings[0]);
995  for (i = 1; i < poly2->nrings; i++)
996  {
997  /* Inside a hole */
998  if ( lwgeom_contains_point(poly2->rings[i], pt) != LW_OUTSIDE )
999  {
1000  return lw_dist2d_recursive(poly1->rings[0], poly2->rings[i], dl);
1001  }
1002  }
1003 
1004 
1005  /*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.*/
1006  pt = lw_curvering_getfirstpoint2d_cp(poly1->rings[0]);
1007  if ( lwgeom_contains_point(poly2->rings[0], pt) != LW_OUTSIDE )
1008  {
1009  dl->distance = 0.0;
1010  dl->p1.x = dl->p2.x = pt->x;
1011  dl->p1.y = dl->p2.y = pt->y;
1012  return LW_TRUE;
1013  }
1014 
1015  pt = lw_curvering_getfirstpoint2d_cp(poly2->rings[0]);
1016  if ( lwgeom_contains_point(poly1->rings[0], pt) != LW_OUTSIDE )
1017  {
1018  dl->distance = 0.0;
1019  dl->p1.x = dl->p2.x = pt->x;
1020  dl->p1.y = dl->p2.y = pt->y;
1021  return LW_TRUE;
1022  }
1023 
1024  lwerror("Unspecified error in function lw_dist2d_curvepoly_curvepoly");
1025  return LW_FALSE;
1026 }
1027 
1028 
1029 
1034 int
1036 {
1037  int t;
1038  const POINT2D *start, *end;
1039  int twist = dl->twisted;
1040 
1041  LWDEBUG(2, "lw_dist2d_pt_ptarray is called");
1042 
1043  start = getPoint2d_cp(pa, 0);
1044 
1045  if ( !lw_dist2d_pt_pt(p, start, dl) ) return LW_FALSE;
1046 
1047  for (t=1; t<pa->npoints; t++)
1048  {
1049  dl->twisted=twist;
1050  end = getPoint2d_cp(pa, t);
1051  if (!lw_dist2d_pt_seg(p, start, end, dl)) return LW_FALSE;
1052 
1053  if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
1054  start = end;
1055  }
1056 
1057  return LW_TRUE;
1058 }
1059 
1064 int
1066 {
1067  int t;
1068  const POINT2D *A1;
1069  const POINT2D *A2;
1070  const POINT2D *A3;
1071  int twist = dl->twisted;
1072 
1073  LWDEBUG(2, "lw_dist2d_pt_ptarrayarc is called");
1074 
1075  if ( pa->npoints % 2 == 0 || pa->npoints < 3 )
1076  {
1077  lwerror("lw_dist2d_pt_ptarrayarc called with non-arc input");
1078  return LW_FALSE;
1079  }
1080 
1081  if (dl->mode == DIST_MAX)
1082  {
1083  lwerror("lw_dist2d_pt_ptarrayarc does not currently support DIST_MAX mode");
1084  return LW_FALSE;
1085  }
1086 
1087  A1 = getPoint2d_cp(pa, 0);
1088 
1089  if ( ! lw_dist2d_pt_pt(p, A1, dl) )
1090  return LW_FALSE;
1091 
1092  for ( t=1; t<pa->npoints; t += 2 )
1093  {
1094  dl->twisted = twist;
1095  A2 = getPoint2d_cp(pa, t);
1096  A3 = getPoint2d_cp(pa, t+1);
1097 
1098  if ( lw_dist2d_pt_arc(p, A1, A2, A3, dl) == LW_FALSE )
1099  return LW_FALSE;
1100 
1101  if ( dl->distance <= dl->tolerance && dl->mode == DIST_MIN )
1102  return LW_TRUE; /*just a check if the answer is already given*/
1103 
1104  A1 = A3;
1105  }
1106 
1107  return LW_TRUE;
1108 }
1109 
1110 
1111 
1112 
1116 int
1118 {
1119  int t,u;
1120  const POINT2D *start, *end;
1121  const POINT2D *start2, *end2;
1122  int twist = dl->twisted;
1123 
1124  LWDEBUGF(2, "lw_dist2d_ptarray_ptarray called (points: %d-%d)",l1->npoints, l2->npoints);
1125 
1126  if (dl->mode == DIST_MAX)/*If we are searching for maxdistance we go straight to point-point calculation since the maxdistance have to be between two vertexes*/
1127  {
1128  for (t=0; t<l1->npoints; t++) /*for each segment in L1 */
1129  {
1130  start = getPoint2d_cp(l1, t);
1131  for (u=0; u<l2->npoints; u++) /*for each segment in L2 */
1132  {
1133  start2 = getPoint2d_cp(l2, u);
1134  lw_dist2d_pt_pt(start, start2, dl);
1135  LWDEBUGF(4, "maxdist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);
1136  LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",
1137  t, u, dl->distance, dl->tolerance);
1138  }
1139  }
1140  }
1141  else
1142  {
1143  start = getPoint2d_cp(l1, 0);
1144  for (t=1; t<l1->npoints; t++) /*for each segment in L1 */
1145  {
1146  end = getPoint2d_cp(l1, t);
1147  start2 = getPoint2d_cp(l2, 0);
1148  for (u=1; u<l2->npoints; u++) /*for each segment in L2 */
1149  {
1150  end2 = getPoint2d_cp(l2, u);
1151  dl->twisted=twist;
1152  lw_dist2d_seg_seg(start, end, start2, end2, dl);
1153  LWDEBUGF(4, "mindist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);
1154  LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",
1155  t, u, dl->distance, dl->tolerance);
1156  if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
1157  start2 = end2;
1158  }
1159  start = end;
1160  }
1161  }
1162  return LW_TRUE;
1163 }
1164 
1168 int
1170 {
1171  int t, u;
1172  const POINT2D *A1;
1173  const POINT2D *A2;
1174  const POINT2D *B1;
1175  const POINT2D *B2;
1176  const POINT2D *B3;
1177  int twist = dl->twisted;
1178 
1179  LWDEBUGF(2, "lw_dist2d_ptarray_ptarrayarc called (points: %d-%d)",pa->npoints, pb->npoints);
1180 
1181  if ( pb->npoints % 2 == 0 || pb->npoints < 3 )
1182  {
1183  lwerror("lw_dist2d_ptarray_ptarrayarc called with non-arc input");
1184  return LW_FALSE;
1185  }
1186 
1187  if ( dl->mode == DIST_MAX )
1188  {
1189  lwerror("lw_dist2d_ptarray_ptarrayarc does not currently support DIST_MAX mode");
1190  return LW_FALSE;
1191  }
1192  else
1193  {
1194  A1 = getPoint2d_cp(pa, 0);
1195  for ( t=1; t < pa->npoints; t++ ) /* For each segment in pa */
1196  {
1197  A2 = getPoint2d_cp(pa, t);
1198  B1 = getPoint2d_cp(pb, 0);
1199  for ( u=1; u < pb->npoints; u += 2 ) /* For each arc in pb */
1200  {
1201  B2 = getPoint2d_cp(pb, u);
1202  B3 = getPoint2d_cp(pb, u+1);
1203  dl->twisted = twist;
1204 
1205  lw_dist2d_seg_arc(A1, A2, B1, B2, B3, dl);
1206 
1207  /* If we've found a distance within tolerance, we're done */
1208  if ( dl->distance <= dl->tolerance && dl->mode == DIST_MIN )
1209  return LW_TRUE;
1210 
1211  B1 = B3;
1212  }
1213  A1 = A2;
1214  }
1215  }
1216  return LW_TRUE;
1217 }
1218 
1222 int
1224 {
1225  int t, u;
1226  const POINT2D *A1;
1227  const POINT2D *A2;
1228  const POINT2D *A3;
1229  const POINT2D *B1;
1230  const POINT2D *B2;
1231  const POINT2D *B3;
1232  int twist = dl->twisted;
1233 
1234  LWDEBUGF(2, "lw_dist2d_ptarrayarc_ptarrayarc called (points: %d-%d)",pa->npoints, pb->npoints);
1235 
1236  if (dl->mode == DIST_MAX)
1237  {
1238  lwerror("lw_dist2d_ptarrayarc_ptarrayarc does not currently support DIST_MAX mode");
1239  return LW_FALSE;
1240  }
1241  else
1242  {
1243  A1 = getPoint2d_cp(pa, 0);
1244  for ( t=1; t < pa->npoints; t += 2 ) /* For each segment in pa */
1245  {
1246  A2 = getPoint2d_cp(pa, t);
1247  A3 = getPoint2d_cp(pa, t+1);
1248  B1 = getPoint2d_cp(pb, 0);
1249  for ( u=1; u < pb->npoints; u += 2 ) /* For each arc in pb */
1250  {
1251  B2 = getPoint2d_cp(pb, u);
1252  B3 = getPoint2d_cp(pb, u+1);
1253  dl->twisted = twist;
1254 
1255  lw_dist2d_arc_arc(A1, A2, A3, B1, B2, B3, dl);
1256 
1257  /* If we've found a distance within tolerance, we're done */
1258  if ( dl->distance <= dl->tolerance && dl->mode == DIST_MIN )
1259  return LW_TRUE;
1260 
1261  B1 = B3;
1262  }
1263  A1 = A3;
1264  }
1265  }
1266  return LW_TRUE;
1267 }
1268 
1273 int
1274 lw_dist2d_seg_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl)
1275 {
1276  POINT2D C; /* center of arc circle */
1277  double radius_C; /* radius of arc circle */
1278  POINT2D D; /* point on A closest to C */
1279  double dist_C_D; /* distance from C to D */
1280  int pt_in_arc, pt_in_seg;
1281  DISTPTS dltmp;
1282 
1283  /* Bail out on crazy modes */
1284  if ( dl->mode < 0 )
1285  lwerror("lw_dist2d_seg_arc does not support maxdistance mode");
1286 
1287  /* What if the "arc" is a point? */
1288  if ( lw_arc_is_pt(B1, B2, B3) )
1289  return lw_dist2d_pt_seg(B1, A1, A2, dl);
1290 
1291  /* Calculate center and radius of the circle. */
1292  radius_C = lw_arc_center(B1, B2, B3, &C);
1293 
1294  /* This "arc" is actually a line (B2 is colinear with B1,B3) */
1295  if ( radius_C < 0.0 )
1296  return lw_dist2d_seg_seg(A1, A2, B1, B3, dl);
1297 
1298  /* Calculate distance between the line and circle center */
1300  if ( lw_dist2d_pt_seg(&C, A1, A2, &dltmp) == LW_FALSE )
1301  lwerror("lw_dist2d_pt_seg failed in lw_dist2d_seg_arc");
1302 
1303  D = dltmp.p1;
1304  dist_C_D = dltmp.distance;
1305 
1306  /* Line intersects circle, maybe arc intersects edge? */
1307  /* If so, that's the closest point. */
1308  /* If not, the closest point is one of the end points of A */
1309  if ( dist_C_D < radius_C )
1310  {
1311  double length_A; /* length of the segment A */
1312  POINT2D E, F; /* points of interection of edge A and circle(B) */
1313  double dist_D_EF; /* distance from D to E or F (same distance both ways) */
1314 
1315  dist_D_EF = sqrt(radius_C*radius_C - dist_C_D*dist_C_D);
1316  length_A = sqrt((A2->x-A1->x)*(A2->x-A1->x)+(A2->y-A1->y)*(A2->y-A1->y));
1317 
1318  /* Point of intersection E */
1319  E.x = D.x - (A2->x-A1->x) * dist_D_EF / length_A;
1320  E.y = D.y - (A2->y-A1->y) * dist_D_EF / length_A;
1321  /* Point of intersection F */
1322  F.x = D.x + (A2->x-A1->x) * dist_D_EF / length_A;
1323  F.y = D.y + (A2->y-A1->y) * dist_D_EF / length_A;
1324 
1325 
1326  /* If E is within A and within B then it's an interesction point */
1327  pt_in_arc = lw_pt_in_arc(&E, B1, B2, B3);
1328  pt_in_seg = lw_pt_in_seg(&E, A1, A2);
1329 
1330  if ( pt_in_arc && pt_in_seg )
1331  {
1332  dl->distance = 0.0;
1333  dl->p1 = E;
1334  dl->p2 = E;
1335  return LW_TRUE;
1336  }
1337 
1338  /* If F is within A and within B then it's an interesction point */
1339  pt_in_arc = lw_pt_in_arc(&F, B1, B2, B3);
1340  pt_in_seg = lw_pt_in_seg(&F, A1, A2);
1341 
1342  if ( pt_in_arc && pt_in_seg )
1343  {
1344  dl->distance = 0.0;
1345  dl->p1 = F;
1346  dl->p2 = F;
1347  return LW_TRUE;
1348  }
1349  }
1350 
1351  /* Line grazes circle, maybe arc intersects edge? */
1352  /* If so, grazing point is the closest point. */
1353  /* If not, the closest point is one of the end points of A */
1354  else if ( dist_C_D == radius_C )
1355  {
1356  /* Closest point D is also the point of grazing */
1357  pt_in_arc = lw_pt_in_arc(&D, B1, B2, B3);
1358  pt_in_seg = lw_pt_in_seg(&D, A1, A2);
1359 
1360  /* Is D contained in both A and B? */
1361  if ( pt_in_arc && pt_in_seg )
1362  {
1363  dl->distance = 0.0;
1364  dl->p1 = D;
1365  dl->p2 = D;
1366  return LW_TRUE;
1367  }
1368  }
1369  /* Line misses circle. */
1370  /* If closest point to A on circle is within B, then that's the closest */
1371  /* Otherwise, the closest point will be an end point of A */
1372  else
1373  {
1374  POINT2D G; /* Point on circle closest to A */
1375  G.x = C.x + (D.x-C.x) * radius_C / dist_C_D;
1376  G.y = C.y + (D.y-C.y) * radius_C / dist_C_D;
1377 
1378  pt_in_arc = lw_pt_in_arc(&G, B1, B2, B3);
1379  pt_in_seg = lw_pt_in_seg(&D, A1, A2);
1380 
1381  /* Closest point is on the interior of A and B */
1382  if ( pt_in_arc && pt_in_seg )
1383  return lw_dist2d_pt_pt(&D, &G, dl);
1384 
1385  }
1386 
1387  /* Now we test the many combinations of end points with either */
1388  /* arcs or edges. Each previous check determined if the closest */
1389  /* potential point was within the arc/segment inscribed on the */
1390  /* line/circle holding the arc/segment. */
1391 
1392  /* Closest point is in the arc, but not in the segment, so */
1393  /* one of the segment end points must be the closest. */
1394  if ( pt_in_arc & ! pt_in_seg )
1395  {
1396  lw_dist2d_pt_arc(A1, B1, B2, B3, dl);
1397  lw_dist2d_pt_arc(A2, B1, B2, B3, dl);
1398  return LW_TRUE;
1399  }
1400  /* or, one of the arc end points is the closest */
1401  else if ( pt_in_seg && ! pt_in_arc )
1402  {
1403  lw_dist2d_pt_seg(B1, A1, A2, dl);
1404  lw_dist2d_pt_seg(B3, A1, A2, dl);
1405  return LW_TRUE;
1406  }
1407  /* Finally, one of the end-point to end-point combos is the closest. */
1408  else
1409  {
1410  lw_dist2d_pt_pt(A1, B1, dl);
1411  lw_dist2d_pt_pt(A1, B3, dl);
1412  lw_dist2d_pt_pt(A2, B1, dl);
1413  lw_dist2d_pt_pt(A2, B3, dl);
1414  return LW_TRUE;
1415  }
1416 
1417  return LW_FALSE;
1418 }
1419 
1420 int
1421 lw_dist2d_pt_arc(const POINT2D* P, const POINT2D* A1, const POINT2D* A2, const POINT2D* A3, DISTPTS* dl)
1422 {
1423  double radius_A, d;
1424  POINT2D C; /* center of circle defined by arc A */
1425  POINT2D X; /* point circle(A) where line from C to P crosses */
1426 
1427  if ( dl->mode < 0 )
1428  lwerror("lw_dist2d_pt_arc does not support maxdistance mode");
1429 
1430  /* What if the arc is a point? */
1431  if ( lw_arc_is_pt(A1, A2, A3) )
1432  return lw_dist2d_pt_pt(P, A1, dl);
1433 
1434  /* Calculate centers and radii of circles. */
1435  radius_A = lw_arc_center(A1, A2, A3, &C);
1436 
1437  /* This "arc" is actually a line (A2 is colinear with A1,A3) */
1438  if ( radius_A < 0.0 )
1439  return lw_dist2d_pt_seg(P, A1, A3, dl);
1440 
1441  /* Distance from point to center */
1442  d = distance2d_pt_pt(&C, P);
1443 
1444  /* P is the center of the circle */
1445  if ( FP_EQUALS(d, 0.0) )
1446  {
1447  dl->distance = radius_A;
1448  dl->p1 = *A1;
1449  dl->p2 = *P;
1450  return LW_TRUE;
1451  }
1452 
1453  /* X is the point on the circle where the line from P to C crosses */
1454  X.x = C.x + (P->x - C.x) * radius_A / d;
1455  X.y = C.y + (P->y - C.y) * radius_A / d;
1456 
1457  /* Is crossing point inside the arc? Or arc is actually circle? */
1458  if ( p2d_same(A1, A3) || lw_pt_in_arc(&X, A1, A2, A3) )
1459  {
1460  lw_dist2d_pt_pt(P, &X, dl);
1461  }
1462  else
1463  {
1464  /* Distance is the minimum of the distances to the arc end points */
1465  lw_dist2d_pt_pt(A1, P, dl);
1466  lw_dist2d_pt_pt(A3, P, dl);
1467  }
1468  return LW_TRUE;
1469 }
1470 
1471 /* Auxiliary function to calculate the distance between 2 concentric arcs*/
1472 int lw_dist2d_arc_arc_concentric( const POINT2D *A1, const POINT2D *A2,
1473  const POINT2D *A3, double radius_A,
1474  const POINT2D *B1, const POINT2D *B2,
1475  const POINT2D *B3, double radius_B,
1476  const POINT2D *CENTER, DISTPTS *dl);
1477 
1478 int
1479 lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3,
1480  const POINT2D *B1, const POINT2D *B2, const POINT2D *B3,
1481  DISTPTS *dl)
1482 {
1483  POINT2D CA, CB; /* Center points of arcs A and B */
1484  double radius_A, radius_B, d; /* Radii of arcs A and B */
1485  POINT2D P; /* Temporary point P */
1486  POINT2D D; /* Mid-point between the centers CA and CB */
1487  int pt_in_arc_A, pt_in_arc_B; /* Test whether potential intersection point is within the arc */
1488 
1489  if ( dl->mode != DIST_MIN )
1490  lwerror("lw_dist2d_arc_arc only supports mindistance");
1491 
1492  /* TODO: Handle case where arc is closed circle (A1 = A3) */
1493 
1494  /* What if one or both of our "arcs" is actually a point? */
1495  if ( lw_arc_is_pt(B1, B2, B3) && lw_arc_is_pt(A1, A2, A3) )
1496  return lw_dist2d_pt_pt(B1, A1, dl);
1497  else if ( lw_arc_is_pt(B1, B2, B3) )
1498  return lw_dist2d_pt_arc(B1, A1, A2, A3, dl);
1499  else if ( lw_arc_is_pt(A1, A2, A3) )
1500  return lw_dist2d_pt_arc(A1, B1, B2, B3, dl);
1501 
1502  /* Calculate centers and radii of circles. */
1503  radius_A = lw_arc_center(A1, A2, A3, &CA);
1504  radius_B = lw_arc_center(B1, B2, B3, &CB);
1505 
1506  /* Two co-linear arcs?!? That's two segments. */
1507  if ( radius_A < 0 && radius_B < 0 )
1508  return lw_dist2d_seg_seg(A1, A3, B1, B3, dl);
1509 
1510  /* A is co-linear, delegate to lw_dist_seg_arc here. */
1511  if ( radius_A < 0 )
1512  return lw_dist2d_seg_arc(A1, A3, B1, B2, B3, dl);
1513 
1514  /* B is co-linear, delegate to lw_dist_seg_arc here. */
1515  if ( radius_B < 0 )
1516  return lw_dist2d_seg_arc(B1, B3, A1, A2, A3, dl);
1517 
1518  /* Center-center distance */
1519  d = distance2d_pt_pt(&CA, &CB);
1520 
1521  /* Concentric arcs */
1522  if ( FP_EQUALS(d, 0.0) )
1523  return lw_dist2d_arc_arc_concentric(A1, A2, A3, radius_A,
1524  B1, B2, B3, radius_B,
1525  &CA, dl);
1526 
1527  /* Make sure that arc "A" has the bigger radius */
1528  if ( radius_B > radius_A )
1529  {
1530  const POINT2D *tmp;
1531  tmp = B1; B1 = A1; A1 = tmp;
1532  tmp = B2; B2 = A2; A2 = tmp;
1533  tmp = B3; B3 = A3; A3 = tmp;
1534  P = CB; CB = CA; CA = P;
1535  d = radius_B; radius_B = radius_A; radius_A = d;
1536  }
1537 
1538  /* Circles touch at a point. Is that point within the arcs? */
1539  if ( d == (radius_A + radius_B) )
1540  {
1541  D.x = CA.x + (CB.x - CA.x) * radius_A / d;
1542  D.y = CA.y + (CB.y - CA.y) * radius_A / d;
1543 
1544  pt_in_arc_A = lw_pt_in_arc(&D, A1, A2, A3);
1545  pt_in_arc_B = lw_pt_in_arc(&D, B1, B2, B3);
1546 
1547  /* Arcs do touch at D, return it */
1548  if ( pt_in_arc_A && pt_in_arc_B )
1549  {
1550  dl->distance = 0.0;
1551  dl->p1 = D;
1552  dl->p2 = D;
1553  return LW_TRUE;
1554  }
1555  }
1556  /* Disjoint or contained circles don't intersect. Closest point may be on */
1557  /* the line joining CA to CB. */
1558  else if ( d > (radius_A + radius_B) /* Disjoint */ || d < (radius_A - radius_B) /* Contained */ )
1559  {
1560  POINT2D XA, XB; /* Points where the line from CA to CB cross their circle bounds */
1561 
1562  /* Calculate hypothetical nearest points, the places on the */
1563  /* two circles where the center-center line crosses. If both */
1564  /* arcs contain their hypothetical points, that's the crossing distance */
1565  XA.x = CA.x + (CB.x - CA.x) * radius_A / d;
1566  XA.y = CA.y + (CB.y - CA.y) * radius_A / d;
1567  XB.x = CB.x + (CA.x - CB.x) * radius_B / d;
1568  XB.y = CB.y + (CA.y - CB.y) * radius_B / d;
1569 
1570  pt_in_arc_A = lw_pt_in_arc(&XA, A1, A2, A3);
1571  pt_in_arc_B = lw_pt_in_arc(&XB, B1, B2, B3);
1572 
1573  /* If the nearest points are both within the arcs, that's our answer */
1574  /* the shortest distance is at the nearest points */
1575  if ( pt_in_arc_A && pt_in_arc_B )
1576  {
1577  return lw_dist2d_pt_pt(&XA, &XB, dl);
1578  }
1579  }
1580  /* Circles cross at two points, are either of those points in both arcs? */
1581  /* http://paulbourke.net/geometry/2circle/ */
1582  else if ( d < (radius_A + radius_B) )
1583  {
1584  POINT2D E, F; /* Points where circle(A) and circle(B) cross */
1585  /* Distance from CA to D */
1586  double a = (radius_A*radius_A - radius_B*radius_B + d*d) / (2*d);
1587  /* Distance from D to E or F */
1588  double h = sqrt(radius_A*radius_A - a*a);
1589 
1590  /* Location of D */
1591  D.x = CA.x + (CB.x - CA.x) * a / d;
1592  D.y = CA.y + (CB.y - CA.y) * a / d;
1593 
1594  /* Start from D and project h units perpendicular to CA-D to get E */
1595  E.x = D.x + (D.y - CA.y) * h / a;
1596  E.y = D.y + (D.x - CA.x) * h / a;
1597 
1598  /* Crossing point E contained in arcs? */
1599  pt_in_arc_A = lw_pt_in_arc(&E, A1, A2, A3);
1600  pt_in_arc_B = lw_pt_in_arc(&E, B1, B2, B3);
1601 
1602  if ( pt_in_arc_A && pt_in_arc_B )
1603  {
1604  dl->p1 = dl->p2 = E;
1605  dl->distance = 0.0;
1606  return LW_TRUE;
1607  }
1608 
1609  /* Start from D and project h units perpendicular to CA-D to get F */
1610  F.x = D.x - (D.y - CA.y) * h / a;
1611  F.y = D.y - (D.x - CA.x) * h / a;
1612 
1613  /* Crossing point F contained in arcs? */
1614  pt_in_arc_A = lw_pt_in_arc(&F, A1, A2, A3);
1615  pt_in_arc_B = lw_pt_in_arc(&F, B1, B2, B3);
1616 
1617  if ( pt_in_arc_A && pt_in_arc_B )
1618  {
1619  dl->p1 = dl->p2 = F;
1620  dl->distance = 0.0;
1621  return LW_TRUE;
1622  }
1623  }
1624  else
1625  {
1626  lwerror("lw_dist2d_arc_arc: arcs neither touch, intersect nor are disjoint! INCONCEIVABLE!");
1627  return LW_FALSE;
1628  }
1629 
1630  /* Closest point is in the arc A, but not in the arc B, so */
1631  /* one of the B end points must be the closest. */
1632  if ( pt_in_arc_A & ! pt_in_arc_B )
1633  {
1634  lw_dist2d_pt_arc(B1, A1, A2, A3, dl);
1635  lw_dist2d_pt_arc(B3, A1, A2, A3, dl);
1636  return LW_TRUE;
1637  }
1638  /* Closest point is in the arc B, but not in the arc A, so */
1639  /* one of the A end points must be the closest. */
1640  else if ( pt_in_arc_B && ! pt_in_arc_A )
1641  {
1642  lw_dist2d_pt_arc(A1, B1, B2, B3, dl);
1643  lw_dist2d_pt_arc(A3, B1, B2, B3, dl);
1644  return LW_TRUE;
1645  }
1646  /* Finally, one of the end-point to end-point combos is the closest. */
1647  else
1648  {
1649  lw_dist2d_pt_pt(A1, B1, dl);
1650  lw_dist2d_pt_pt(A1, B3, dl);
1651  lw_dist2d_pt_pt(A2, B1, dl);
1652  lw_dist2d_pt_pt(A2, B3, dl);
1653  return LW_TRUE;
1654  }
1655 
1656  return LW_TRUE;
1657 }
1658 
1659 int
1661  const POINT2D *A3, double radius_A,
1662  const POINT2D *B1, const POINT2D *B2,
1663  const POINT2D *B3, double radius_B,
1664  const POINT2D *CENTER, DISTPTS *dl)
1665 {
1666  int seg_size;
1667  double dist_sqr, shortest_sqr;
1668  const POINT2D *P1;
1669  const POINT2D *P2;
1670  POINT2D proj;
1671 
1672  if (radius_A == radius_B)
1673  {
1674  /* Check if B1 or B3 are in the same side as A2 in the A1-A3 arc */
1675  seg_size = lw_segment_side(A1, A3, A2);
1676  if (seg_size == lw_segment_side(A1, A3, B1))
1677  {
1678  dl->p1 = *B1;
1679  dl->p2 = *B1;
1680  dl->distance = 0;
1681  return LW_TRUE;
1682  }
1683  if (seg_size == lw_segment_side(A1, A3, B3))
1684  {
1685  dl->p1 = *B3;
1686  dl->p2 = *B3;
1687  dl->distance = 0;
1688  return LW_TRUE;
1689  }
1690  /* Check if A1 or A3 are in the same side as B2 in the B1-B3 arc */
1691  seg_size = lw_segment_side(B1, B3, B2);
1692  if (seg_size == lw_segment_side(B1, B3, A1))
1693  {
1694  dl->p1 = *A1;
1695  dl->p2 = *A1;
1696  dl->distance = 0;
1697  return LW_TRUE;
1698  }
1699  if (seg_size == lw_segment_side(B1, B3, A3))
1700  {
1701  dl->p1 = *A3;
1702  dl->p2 = *A3;
1703  dl->distance = 0;
1704  return LW_TRUE;
1705  }
1706  }
1707  else
1708  {
1709  /* Check if any projection of B ends are in A*/
1710  seg_size = lw_segment_side(A1, A3, A2);
1711 
1712  /* B1 */
1713  proj.x = CENTER->x + (B1->x - CENTER->x) * radius_A / radius_B;
1714  proj.y = CENTER->y + (B1->y - CENTER->y) * radius_A / radius_B;
1715 
1716  if (seg_size == lw_segment_side(A1, A3, &proj))
1717  {
1718  dl->p1 = proj;
1719  dl->p2 = *B1;
1720  dl->distance = fabs(radius_A - radius_B);
1721  return LW_TRUE;
1722  }
1723  /* B3 */
1724  proj.x = CENTER->x + (B3->x - CENTER->x) * radius_A / radius_B;
1725  proj.y = CENTER->y + (B3->y - CENTER->y) * radius_A / radius_B;
1726  if (seg_size == lw_segment_side(A1, A3, &proj))
1727  {
1728  dl->p1 = proj;
1729  dl->p2 = *B3;
1730  dl->distance = fabs(radius_A - radius_B);
1731  return LW_TRUE;
1732  }
1733 
1734  /* Now check projections of A in B */
1735  seg_size = lw_segment_side(B1, B3, B2);
1736 
1737  /* A1 */
1738  proj.x = CENTER->x + (A1->x - CENTER->x) * radius_B / radius_A;
1739  proj.y = CENTER->y + (A1->y - CENTER->y) * radius_B / radius_A;
1740  if (seg_size == lw_segment_side(B1, B3, &proj))
1741  {
1742  dl->p1 = proj;
1743  dl->p2 = *A1;
1744  dl->distance = fabs(radius_A - radius_B);
1745  return LW_TRUE;
1746  }
1747 
1748  /* A3 */
1749  proj.x = CENTER->x + (A3->x - CENTER->x) * radius_B / radius_A;
1750  proj.y = CENTER->y + (A3->y - CENTER->y) * radius_B / radius_A;
1751  if (seg_size == lw_segment_side(B1, B3, &proj))
1752  {
1753  dl->p1 = proj;
1754  dl->p2 = *A3;
1755  dl->distance = fabs(radius_A - radius_B);
1756  return LW_TRUE;
1757  }
1758  }
1759 
1760  /* Check the shortest between the distances of the 4 ends */
1761  shortest_sqr = dist_sqr = distance2d_sqr_pt_pt(A1, B1);
1762  P1 = A1;
1763  P2 = B1;
1764 
1765  dist_sqr = distance2d_sqr_pt_pt(A1, B3);
1766  if (dist_sqr < shortest_sqr)
1767  {
1768  shortest_sqr = dist_sqr;
1769  P1 = A1;
1770  P2 = B3;
1771  }
1772 
1773  dist_sqr = distance2d_sqr_pt_pt(A3, B1);
1774  if (dist_sqr < shortest_sqr)
1775  {
1776  shortest_sqr = dist_sqr;
1777  P1 = A3;
1778  P2 = B1;
1779  }
1780 
1781  dist_sqr = distance2d_sqr_pt_pt(A3, B3);
1782  if (dist_sqr < shortest_sqr)
1783  {
1784  shortest_sqr = dist_sqr;
1785  P1 = A3;
1786  P2 = B3;
1787  }
1788 
1789  dl->p1 = *P1;
1790  dl->p2 = *P2;
1791  dl->distance = sqrt(shortest_sqr);
1792 
1793  return LW_TRUE;
1794 }
1795 
1801 int
1802 lw_dist2d_seg_seg(const POINT2D *A, const POINT2D *B, const POINT2D *C, const POINT2D *D, DISTPTS *dl)
1803 {
1804  double s_top, s_bot,s;
1805  double r_top, r_bot,r;
1806 
1807  LWDEBUGF(2, "lw_dist2d_seg_seg [%g,%g]->[%g,%g] by [%g,%g]->[%g,%g]",
1808  A->x,A->y,B->x,B->y, C->x,C->y, D->x, D->y);
1809 
1810  /*A and B are the same point */
1811  if ( ( A->x == B->x) && (A->y == B->y) )
1812  {
1813  return lw_dist2d_pt_seg(A,C,D,dl);
1814  }
1815  /*U and V are the same point */
1816 
1817  if ( ( C->x == D->x) && (C->y == D->y) )
1818  {
1819  dl->twisted= ((dl->twisted) * (-1));
1820  return lw_dist2d_pt_seg(D,A,B,dl);
1821  }
1822  /* AB and CD are line segments */
1823  /* from comp.graphics.algo
1824 
1825  Solving the above for r and s yields
1826  (Ay-Cy)(Dx-Cx)-(Ax-Cx)(Dy-Cy)
1827  r = ----------------------------- (eqn 1)
1828  (Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)
1829 
1830  (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
1831  s = ----------------------------- (eqn 2)
1832  (Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)
1833  Let P be the position vector of the intersection point, then
1834  P=A+r(B-A) or
1835  Px=Ax+r(Bx-Ax)
1836  Py=Ay+r(By-Ay)
1837  By examining the values of r & s, you can also determine some other limiting conditions:
1838  If 0<=r<=1 & 0<=s<=1, intersection exists
1839  r<0 or r>1 or s<0 or s>1 line segments do not intersect
1840  If the denominator in eqn 1 is zero, AB & CD are parallel
1841  If the numerator in eqn 1 is also zero, AB & CD are collinear.
1842 
1843  */
1844  r_top = (A->y-C->y)*(D->x-C->x) - (A->x-C->x)*(D->y-C->y);
1845  r_bot = (B->x-A->x)*(D->y-C->y) - (B->y-A->y)*(D->x-C->x);
1846 
1847  s_top = (A->y-C->y)*(B->x-A->x) - (A->x-C->x)*(B->y-A->y);
1848  s_bot = (B->x-A->x)*(D->y-C->y) - (B->y-A->y)*(D->x-C->x);
1849 
1850  if ( (r_bot==0) || (s_bot == 0) )
1851  {
1852  if ((lw_dist2d_pt_seg(A,C,D,dl)) && (lw_dist2d_pt_seg(B,C,D,dl)))
1853  {
1854  dl->twisted= ((dl->twisted) * (-1)); /*here we change the order of inputted geometrys and that we notice by changing sign on dl->twisted*/
1855  return ((lw_dist2d_pt_seg(C,A,B,dl)) && (lw_dist2d_pt_seg(D,A,B,dl))); /*if all is successful we return true*/
1856  }
1857  else
1858  {
1859  return LW_FALSE; /* if any of the calls to lw_dist2d_pt_seg goes wrong we return false*/
1860  }
1861  }
1862 
1863  s = s_top/s_bot;
1864  r= r_top/r_bot;
1865 
1866  if (((r<0) || (r>1) || (s<0) || (s>1)) || (dl->mode == DIST_MAX))
1867  {
1868  if ((lw_dist2d_pt_seg(A,C,D,dl)) && (lw_dist2d_pt_seg(B,C,D,dl)))
1869  {
1870  dl->twisted= ((dl->twisted) * (-1)); /*here we change the order of inputted geometrys and that we notice by changing sign on dl->twisted*/
1871  return ((lw_dist2d_pt_seg(C,A,B,dl)) && (lw_dist2d_pt_seg(D,A,B,dl))); /*if all is successful we return true*/
1872  }
1873  else
1874  {
1875  return LW_FALSE; /* if any of the calls to lw_dist2d_pt_seg goes wrong we return false*/
1876  }
1877  }
1878  else
1879  {
1880  if (dl->mode == DIST_MIN) /*If there is intersection we identify the intersection point and return it but only if we are looking for mindistance*/
1881  {
1882  POINT2D theP;
1883 
1884  if (((A->x==C->x)&&(A->y==C->y))||((A->x==D->x)&&(A->y==D->y)))
1885  {
1886  theP.x = A->x;
1887  theP.y = A->y;
1888  }
1889  else if (((B->x==C->x)&&(B->y==C->y))||((B->x==D->x)&&(B->y==D->y)))
1890  {
1891  theP.x = B->x;
1892  theP.y = B->y;
1893  }
1894  else
1895  {
1896  theP.x = A->x+r*(B->x-A->x);
1897  theP.y = A->y+r*(B->y-A->y);
1898  }
1899  dl->distance=0.0;
1900  dl->p1=theP;
1901  dl->p2=theP;
1902  }
1903  return LW_TRUE;
1904 
1905  }
1906  lwerror("unspecified error in function lw_dist2d_seg_seg");
1907  return LW_FALSE; /*If we have come here something is wrong*/
1908 }
1909 
1910 
1911 /*------------------------------------------------------------------------------------------------------------
1912 End of Brute force functions
1913 --------------------------------------------------------------------------------------------------------------*/
1914 
1915 
1916 /*------------------------------------------------------------------------------------------------------------
1917 New faster distance calculations
1918 --------------------------------------------------------------------------------------------------------------*/
1919 
1927 int
1929 {
1930  /*here we define two lists to hold our calculated "z"-values and the order number in the geometry*/
1931 
1932  double k, thevalue;
1933  float deltaX, deltaY, c1m, c2m;
1934  POINT2D c1, c2;
1935  const POINT2D *theP;
1936  float min1X, max1X, max1Y, min1Y,min2X, max2X, max2Y, min2Y;
1937  int t;
1938  int n1 = l1->npoints;
1939  int n2 = l2->npoints;
1940 
1941  LISTSTRUCT *list1, *list2;
1942  list1 = (LISTSTRUCT*)lwalloc(sizeof(LISTSTRUCT)*n1);
1943  list2 = (LISTSTRUCT*)lwalloc(sizeof(LISTSTRUCT)*n2);
1944 
1945  LWDEBUG(2, "lw_dist2d_fast_ptarray_ptarray is called");
1946 
1947  max1X = box1->xmax;
1948  min1X = box1->xmin;
1949  max1Y = box1->ymax;
1950  min1Y = box1->ymin;
1951  max2X = box2->xmax;
1952  min2X = box2->xmin;
1953  max2Y = box2->ymax;
1954  min2Y = box2->ymin;
1955  /*we want the center of the bboxes, and calculate the slope between the centerpoints*/
1956  c1.x = min1X + (max1X-min1X)/2;
1957  c1.y = min1Y + (max1Y-min1Y)/2;
1958  c2.x = min2X + (max2X-min2X)/2;
1959  c2.y = min2Y + (max2Y-min2Y)/2;
1960 
1961  deltaX=(c2.x-c1.x);
1962  deltaY=(c2.y-c1.y);
1963 
1964 
1965  /*Here we calculate where the line perpendicular to the center-center line crosses the axes for each vertex
1966  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 */
1967  if ((deltaX*deltaX)<(deltaY*deltaY)) /*North or South*/
1968  {
1969  k = -deltaX/deltaY;
1970  for (t=0; t<n1; t++) /*for each segment in L1 */
1971  {
1972  theP = getPoint2d_cp(l1, t);
1973  thevalue = theP->y - (k * theP->x);
1974  list1[t].themeasure=thevalue;
1975  list1[t].pnr=t;
1976 
1977  }
1978  for (t=0; t<n2; t++) /*for each segment in L2*/
1979  {
1980  theP = getPoint2d_cp(l2, t);
1981  thevalue = theP->y - (k * theP->x);
1982  list2[t].themeasure=thevalue;
1983  list2[t].pnr=t;
1984 
1985  }
1986  c1m = c1.y-(k*c1.x);
1987  c2m = c2.y-(k*c2.x);
1988  }
1989 
1990 
1991  /*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
1992  and we find it's crossing the X-axes with z = x-(1/k)y */
1993  else /*West or East*/
1994  {
1995  k = -deltaY/deltaX;
1996  for (t=0; t<n1; t++) /*for each segment in L1 */
1997  {
1998  theP = getPoint2d_cp(l1, t);
1999  thevalue = theP->x - (k * theP->y);
2000  list1[t].themeasure=thevalue;
2001  list1[t].pnr=t;
2002  /* lwnotice("l1 %d, measure=%f",t,thevalue ); */
2003  }
2004  for (t=0; t<n2; t++) /*for each segment in L2*/
2005  {
2006  theP = getPoint2d_cp(l2, t);
2007  thevalue = theP->x - (k * theP->y);
2008  list2[t].themeasure=thevalue;
2009  list2[t].pnr=t;
2010  /* lwnotice("l2 %d, measure=%f",t,thevalue ); */
2011  }
2012  c1m = c1.x-(k*c1.y);
2013  c2m = c2.x-(k*c2.y);
2014  }
2015 
2016  /*we sort our lists by the calculated values*/
2017  qsort(list1, n1, sizeof(LISTSTRUCT), struct_cmp_by_measure);
2018  qsort(list2, n2, sizeof(LISTSTRUCT), struct_cmp_by_measure);
2019 
2020  if (c1m < c2m)
2021  {
2022  if (!lw_dist2d_pre_seg_seg(l1,l2,list1,list2,k,dl))
2023  {
2024  lwfree(list1);
2025  lwfree(list2);
2026  return LW_FALSE;
2027  }
2028  }
2029  else
2030  {
2031  dl->twisted= ((dl->twisted) * (-1));
2032  if (!lw_dist2d_pre_seg_seg(l2,l1,list2,list1,k,dl))
2033  {
2034  lwfree(list1);
2035  lwfree(list2);
2036  return LW_FALSE;
2037  }
2038  }
2039  lwfree(list1);
2040  lwfree(list2);
2041  return LW_TRUE;
2042 }
2043 
2044 int
2045 struct_cmp_by_measure(const void *a, const void *b)
2046 {
2047  LISTSTRUCT *ia = (LISTSTRUCT*)a;
2048  LISTSTRUCT *ib = (LISTSTRUCT*)b;
2049  return ( ia->themeasure>ib->themeasure ) ? 1 : -1;
2050 }
2051 
2055 int
2057 {
2058  const POINT2D *p1, *p2, *p3, *p4, *p01, *p02;
2059  int pnr1,pnr2,pnr3,pnr4, n1, n2, i, u, r, twist;
2060  double maxmeasure;
2061  n1= l1->npoints;
2062  n2 = l2->npoints;
2063 
2064  LWDEBUG(2, "lw_dist2d_pre_seg_seg is called");
2065 
2066  p1 = getPoint2d_cp(l1, list1[0].pnr);
2067  p3 = getPoint2d_cp(l2, list2[0].pnr);
2068  lw_dist2d_pt_pt(p1, p3, dl);
2069  maxmeasure = sqrt(dl->distance*dl->distance + (dl->distance*dl->distance*k*k));
2070  twist = dl->twisted; /*to keep the incomming order between iterations*/
2071  for (i =(n1-1); i>=0; --i)
2072  {
2073  /*we break this iteration when we have checked every
2074  point closer to our perpendicular "checkline" than
2075  our shortest found distance*/
2076  if (((list2[0].themeasure-list1[i].themeasure)) > maxmeasure) break;
2077  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*/
2078  {
2079  pnr1 = list1[i].pnr;
2080  p1 = getPoint2d_cp(l1, pnr1);
2081  if (pnr1+r<0)
2082  {
2083  p01 = getPoint2d_cp(l1, (n1-1));
2084  if (( p1->x == p01->x) && (p1->y == p01->y)) pnr2 = (n1-1);
2085  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*/
2086  }
2087 
2088  else if (pnr1+r>(n1-1))
2089  {
2090  p01 = getPoint2d_cp(l1, 0);
2091  if (( p1->x == p01->x) && (p1->y == p01->y)) pnr2 = 0;
2092  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*/
2093  }
2094  else pnr2 = pnr1+r;
2095 
2096 
2097  p2 = getPoint2d_cp(l1, pnr2);
2098  for (u=0; u<n2; ++u)
2099  {
2100  if (((list2[u].themeasure-list1[i].themeasure)) >= maxmeasure) break;
2101  pnr3 = list2[u].pnr;
2102  p3 = getPoint2d_cp(l2, pnr3);
2103  if (pnr3==0)
2104  {
2105  p02 = getPoint2d_cp(l2, (n2-1));
2106  if (( p3->x == p02->x) && (p3->y == p02->y)) pnr4 = (n2-1);
2107  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*/
2108  }
2109  else pnr4 = pnr3-1;
2110 
2111  p4 = getPoint2d_cp(l2, pnr4);
2112  dl->twisted=twist;
2113  if (!lw_dist2d_selected_seg_seg(p1, p2, p3, p4, dl)) return LW_FALSE;
2114 
2115  if (pnr3>=(n2-1))
2116  {
2117  p02 = getPoint2d_cp(l2, 0);
2118  if (( p3->x == p02->x) && (p3->y == p02->y)) pnr4 = 0;
2119  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*/
2120  }
2121 
2122  else pnr4 = pnr3+1;
2123 
2124  p4 = getPoint2d_cp(l2, pnr4);
2125  dl->twisted=twist; /*we reset the "twist" for each iteration*/
2126  if (!lw_dist2d_selected_seg_seg(p1, p2, p3, p4, dl)) return LW_FALSE;
2127 
2128  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*/
2129  }
2130  }
2131  }
2132 
2133  return LW_TRUE;
2134 }
2135 
2136 
2142 int
2143 lw_dist2d_selected_seg_seg(const POINT2D *A, const POINT2D *B, const POINT2D *C, const POINT2D *D, DISTPTS *dl)
2144 {
2145  LWDEBUGF(2, "lw_dist2d_selected_seg_seg [%g,%g]->[%g,%g] by [%g,%g]->[%g,%g]",
2146  A->x,A->y,B->x,B->y, C->x,C->y, D->x, D->y);
2147 
2148  /*A and B are the same point */
2149  if ( ( A->x == B->x) && (A->y == B->y) )
2150  {
2151  return lw_dist2d_pt_seg(A,C,D,dl);
2152  }
2153  /*U and V are the same point */
2154 
2155  if ( ( C->x == D->x) && (C->y == D->y) )
2156  {
2157  dl->twisted= ((dl->twisted) * (-1));
2158  return lw_dist2d_pt_seg(D,A,B,dl);
2159  }
2160 
2161  if ((lw_dist2d_pt_seg(A,C,D,dl)) && (lw_dist2d_pt_seg(B,C,D,dl)))
2162  {
2163  dl->twisted= ((dl->twisted) * (-1)); /*here we change the order of inputted geometrys and that we notice by changing sign on dl->twisted*/
2164  return ((lw_dist2d_pt_seg(C,A,B,dl)) && (lw_dist2d_pt_seg(D,A,B,dl))); /*if all is successful we return true*/
2165  }
2166  else
2167  {
2168  return LW_FALSE; /* if any of the calls to lw_dist2d_pt_seg goes wrong we return false*/
2169  }
2170 }
2171 
2172 /*------------------------------------------------------------------------------------------------------------
2173 End of New faster distance calculations
2174 --------------------------------------------------------------------------------------------------------------*/
2175 
2176 
2177 /*------------------------------------------------------------------------------------------------------------
2178 Functions in common for Brute force and new calculation
2179 --------------------------------------------------------------------------------------------------------------*/
2180 
2188 int
2189 lw_dist2d_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B, DISTPTS *dl)
2190 {
2191  POINT2D c;
2192  double r;
2193  /*if start==end, then use pt distance */
2194  if ( ( A->x == B->x) && (A->y == B->y) )
2195  {
2196  return lw_dist2d_pt_pt(p,A,dl);
2197  }
2198  /*
2199  * otherwise, we use comp.graphics.algorithms
2200  * Frequently Asked Questions method
2201  *
2202  * (1) AC dot AB
2203  * r = ---------
2204  * ||AB||^2
2205  * r has the following meaning:
2206  * r=0 P = A
2207  * r=1 P = B
2208  * r<0 P is on the backward extension of AB
2209  * r>1 P is on the forward extension of AB
2210  * 0<r<1 P is interior to AB
2211  */
2212 
2213  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) );
2214 
2215  /*This is for finding the maxdistance.
2216  the maxdistance have to be between two vertexes,
2217  compared to mindistance which can be between
2218  tvo vertexes vertex.*/
2219  if (dl->mode == DIST_MAX)
2220  {
2221  if (r>=0.5)
2222  {
2223  return lw_dist2d_pt_pt(p,A,dl);
2224  }
2225  if (r<0.5)
2226  {
2227  return lw_dist2d_pt_pt(p,B,dl);
2228  }
2229  }
2230 
2231  if (r<0) /*If p projected on the line is outside point A*/
2232  {
2233  return lw_dist2d_pt_pt(p,A,dl);
2234  }
2235  if (r>=1) /*If p projected on the line is outside point B or on point B*/
2236  {
2237  return lw_dist2d_pt_pt(p,B,dl);
2238  }
2239 
2240  /*If the point p is on the segment this is a more robust way to find out that*/
2241  if (( ((A->y-p->y)*(B->x-A->x)==(A->x-p->x)*(B->y-A->y) ) ) && (dl->mode == DIST_MIN))
2242  {
2243  dl->distance = 0.0;
2244  dl->p1 = *p;
2245  dl->p2 = *p;
2246  }
2247 
2248  /*If the projection of point p on the segment is between A and B
2249  then we find that "point on segment" and send it to lw_dist2d_pt_pt*/
2250  c.x=A->x + r * (B->x-A->x);
2251  c.y=A->y + r * (B->y-A->y);
2252 
2253  return lw_dist2d_pt_pt(p,&c,dl);
2254 }
2255 
2256 
2264 int
2265 lw_dist2d_pt_pt(const POINT2D *thep1, const POINT2D *thep2, DISTPTS *dl)
2266 {
2267  double hside = thep2->x - thep1->x;
2268  double vside = thep2->y - thep1->y;
2269  double dist = sqrt ( hside*hside + vside*vside );
2270 
2271  if (((dl->distance - dist)*(dl->mode))>0) /*multiplication with mode to handle mindistance (mode=1) and maxdistance (mode = (-1)*/
2272  {
2273  dl->distance = dist;
2274 
2275  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*/
2276  {
2277  dl->p1 = *thep1;
2278  dl->p2 = *thep2;
2279  }
2280  else
2281  {
2282  dl->p1 = *thep2;
2283  dl->p2 = *thep1;
2284  }
2285  }
2286  return LW_TRUE;
2287 }
2288 
2289 
2290 
2291 
2292 /*------------------------------------------------------------------------------------------------------------
2293 End of Functions in common for Brute force and new calculation
2294 --------------------------------------------------------------------------------------------------------------*/
2295 
2296 
2300 double
2302 {
2303  double hside = p2->x - p1->x;
2304  double vside = p2->y - p1->y;
2305 
2306  return sqrt ( hside*hside + vside*vside );
2307 
2308 }
2309 
2310 double
2312 {
2313  double hside = p2->x - p1->x;
2314  double vside = p2->y - p1->y;
2315 
2316  return hside*hside + vside*vside;
2317 
2318 }
2319 
2320 
2325 double
2326 distance2d_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B)
2327 {
2328  double r,s;
2329 
2330  /*if start==end, then use pt distance */
2331  if ( ( A->x == B->x) && (A->y == B->y) )
2332  return distance2d_pt_pt(p,A);
2333 
2334  /*
2335  * otherwise, we use comp.graphics.algorithms
2336  * Frequently Asked Questions method
2337  *
2338  * (1) AC dot AB
2339  * r = ---------
2340  * ||AB||^2
2341  * r has the following meaning:
2342  * r=0 P = A
2343  * r=1 P = B
2344  * r<0 P is on the backward extension of AB
2345  * r>1 P is on the forward extension of AB
2346  * 0<r<1 P is interior to AB
2347  */
2348 
2349  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) );
2350 
2351  if (r<0) return distance2d_pt_pt(p,A);
2352  if (r>1) return distance2d_pt_pt(p,B);
2353 
2354 
2355  /*
2356  * (2)
2357  * (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
2358  * s = -----------------------------
2359  * L^2
2360  *
2361  * Then the distance from C to P = |s|*L.
2362  *
2363  */
2364 
2365  s = ( (A->y-p->y)*(B->x-A->x)- (A->x-p->x)*(B->y-A->y) ) /
2366  ( (B->x-A->x)*(B->x-A->x) +(B->y-A->y)*(B->y-A->y) );
2367 
2368  return FP_ABS(s) * sqrt(
2369  (B->x-A->x)*(B->x-A->x) + (B->y-A->y)*(B->y-A->y)
2370  );
2371 }
2372 
2373 /* return distance squared, useful to avoid sqrt calculations */
2374 double
2375 distance2d_sqr_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B)
2376 {
2377  double r,s;
2378 
2379  if ( ( A->x == B->x) && (A->y == B->y) )
2380  return distance2d_sqr_pt_pt(p,A);
2381 
2382  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) );
2383 
2384  if (r<0) return distance2d_sqr_pt_pt(p,A);
2385  if (r>1) return distance2d_sqr_pt_pt(p,B);
2386 
2387 
2388  /*
2389  * (2)
2390  * (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
2391  * s = -----------------------------
2392  * L^2
2393  *
2394  * Then the distance from C to P = |s|*L.
2395  *
2396  */
2397 
2398  s = ( (A->y-p->y)*(B->x-A->x)- (A->x-p->x)*(B->y-A->y) ) /
2399  ( (B->x-A->x)*(B->x-A->x) +(B->y-A->y)*(B->y-A->y) );
2400 
2401  return s * s * ( (B->x-A->x)*(B->x-A->x) + (B->y-A->y)*(B->y-A->y) );
2402 }
2403 
2404 
2405 
2410 int
2411 azimuth_pt_pt(const POINT2D *A, const POINT2D *B, double *d)
2412 {
2413  if ( A->x == B->x )
2414  {
2415  if ( A->y < B->y ) *d=0.0;
2416  else if ( A->y > B->y ) *d=M_PI;
2417  else return 0;
2418  return 1;
2419  }
2420 
2421  if ( A->y == B->y )
2422  {
2423  if ( A->x < B->x ) *d=M_PI/2;
2424  else if ( A->x > B->x ) *d=M_PI+(M_PI/2);
2425  else return 0;
2426  return 1;
2427  }
2428 
2429  if ( A->x < B->x )
2430  {
2431  if ( A->y < B->y )
2432  {
2433  *d=atan(fabs(A->x - B->x) / fabs(A->y - B->y) );
2434  }
2435  else /* ( A->y > B->y ) - equality case handled above */
2436  {
2437  *d=atan(fabs(A->y - B->y) / fabs(A->x - B->x) )
2438  + (M_PI/2);
2439  }
2440  }
2441 
2442  else /* ( A->x > B->x ) - equality case handled above */
2443  {
2444  if ( A->y > B->y )
2445  {
2446  *d=atan(fabs(A->x - B->x) / fabs(A->y - B->y) )
2447  + M_PI;
2448  }
2449  else /* ( A->y < B->y ) - equality case handled above */
2450  {
2451  *d=atan(fabs(A->y - B->y) / fabs(A->x - B->x) )
2452  + (M_PI+(M_PI/2));
2453  }
2454  }
2455 
2456  return 1;
2457 }
2458 
2459 
double distance2d_sqr_pt_pt(const POINT2D *p1, const POINT2D *p2)
Definition: measures.c:2311
double lw_arc_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result)
Determines the center of the circle defined by the three given points.
Definition: lwalgorithm.c:213
#define LINETYPE
Definition: liblwgeom.h:71
GBOX * bbox
Definition: liblwgeom.h:382
int lw_dist2d_line_curvepoly(LWLINE *line, LWCURVEPOLY *poly, DISTPTS *dl)
Definition: measures.c:782
int lw_arc_is_pt(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
Returns true if arc A is actually a point (all vertices are the same) .
Definition: lwalgorithm.c:91
int lw_dist2d_circstring_circstring(LWCIRCSTRING *line1, LWCIRCSTRING *line2, DISTPTS *dl)
Definition: measures.c:929
double lwgeom_maxdistance2d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
Function handling max distance calculations and dfyllywithin calculations.
Definition: measures.c:167
#define MULTICURVETYPE
Definition: liblwgeom.h:80
LWGEOM * lw_dist2d_distancepoint(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
Function initializing closestpoint calculations.
Definition: measures.c:117
int lw_dist2d_arc_arc_concentric(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, double radius_A, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, double radius_B, const POINT2D *CENTER, DISTPTS *dl)
Definition: measures.c:1660
int lw_dist2d_comp(const LWGEOM *lw1, const LWGEOM *lw2, DISTPTS *dl)
This function just deserializes geometries Bboxes is not checked here since it is the subgeometries b...
Definition: measures.c:231
int lw_dist2d_circstring_poly(LWCIRCSTRING *circ, LWPOLY *poly, DISTPTS *dl)
Definition: measures.c:913
int lw_dist2d_point_curvepoly(LWPOINT *point, LWCURVEPOLY *poly, DISTPTS *dl)
Definition: measures.c:645
char * r
Definition: cu_in_wkt.c:24
void lwfree(void *mem)
Definition: lwutil.c:214
LWGEOM ** rings
Definition: liblwgeom.h:519
double distance2d_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B)
The old function nessecary for ptarray_segmentize2d in ptarray.c.
Definition: measures.c:2326
int npoints
Definition: liblwgeom.h:355
double distance2d_sqr_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B)
Definition: measures.c:2375
int lw_dist2d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS *dl)
This is a recursive function delivering every possible combinatin of subgeometries.
Definition: measures.c:263
int struct_cmp_by_measure(const void *a, const void *b)
Definition: measures.c:2045
int lw_dist2d_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS *dl)
point to point calculation
Definition: measures.c:561
#define POLYGONTYPE
Definition: liblwgeom.h:72
LWPOINT * lwpoint_make2d(int srid, double x, double y)
Definition: lwpoint.c:132
int mode
Definition: measures.h:27
double xmax
Definition: liblwgeom.h:277
LWLINE * lwline_from_ptarray(int srid, uint32_t npoints, LWPOINT **points)
Definition: lwline.c:214
#define CURVEPOLYTYPE
Definition: liblwgeom.h:79
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1050
#define COMPOUNDTYPE
Definition: liblwgeom.h:78
#define MULTIPOINTTYPE
Definition: liblwgeom.h:73
int lw_dist2d_circstring_curvepoly(LWCIRCSTRING *circ, LWCURVEPOLY *poly, DISTPTS *dl)
Definition: measures.c:923
LWGEOM * lwgeom_closest_line(const LWGEOM *lw1, const LWGEOM *lw2)
Definition: measures.c:28
POINT2D * p1
Definition: lwtree.h:12
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:50
#define POLYHEDRALSURFACETYPE
Definition: liblwgeom.h:82
int lw_dist2d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS *dl)
Definition: measures.c:597
double lwgeom_mindistance2d(const LWGEOM *lw1, const LWGEOM *lw2)
Function initialazing min distance calculation.
Definition: measures.c:188
LWGEOM ** geoms
Definition: liblwgeom.h:506
LWGEOM * lwgeom_furthest_line(const LWGEOM *lw1, const LWGEOM *lw2)
Definition: measures.c:34
POINTARRAY * point
Definition: liblwgeom.h:395
int32_t srid
Definition: liblwgeom.h:383
int lw_dist2d_fast_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2, DISTPTS *dl, GBOX *box1, GBOX *box2)
The new faster calculation comparing pointarray to another pointarray the arrays can come from both p...
Definition: measures.c:1928
int lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl)
Definition: measures.c:1479
int lw_dist2d_curvepoly_curvepoly(LWCURVEPOLY *poly1, LWCURVEPOLY *poly2, DISTPTS *dl)
Definition: measures.c:956
POINT2D * p2
Definition: lwtree.h:13
int lw_dist2d_pre_seg_seg(POINTARRAY *l1, POINTARRAY *l2, LISTSTRUCT *list1, LISTSTRUCT *list2, double k, DISTPTS *dl)
preparation before lw_dist2d_seg_seg.
Definition: measures.c:2056
LWGEOM * lw_dist2d_distanceline(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
Function initializing shortestline and longestline calculations.
Definition: measures.c:70
POINT2D p1
Definition: measures.h:25
double themeasure
Definition: measures.h:34
double tolerance
Definition: measures.h:29
int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox)
Calculate bounding box of a geometry, automatically taking into account whether it is cartesian or ge...
Definition: lwgeom.c:652
double x
Definition: liblwgeom.h:312
int p2d_same(const POINT2D *p1, const POINT2D *p2)
Definition: lwalgorithm.c:35
int lw_dist2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS *dl)
Function handling polygon to polygon calculation 1 if we are looking for maxdistance, just check the outer rings.
Definition: measures.c:829
LWCURVEPOLY * lwcurvepoly_construct_from_lwpoly(LWPOLY *lwpoly)
Construct an equivalent curve polygon from a polygon.
Definition: lwcurvepoly.c:40
static int lw_dist2d_is_collection(const LWGEOM *g)
Definition: measures.c:239
double ymin
Definition: liblwgeom.h:278
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:188
#define DIST_MIN
Definition: measures.h:17
double xmin
Definition: liblwgeom.h:276
#define LW_FALSE
Definition: liblwgeom.h:62
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, int n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:472
LWGEOM * lwgeom_closest_point(const LWGEOM *lw1, const LWGEOM *lw2)
Definition: measures.c:40
int azimuth_pt_pt(const POINT2D *A, const POINT2D *B, double *d)
Compute the azimuth of segment AB in radians.
Definition: measures.c:2411
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:61
int lw_dist2d_ptarray_ptarrayarc(const POINTARRAY *pa, const POINTARRAY *pb, DISTPTS *dl)
Test each segment of pa against each arc of pb for distance.
Definition: measures.c:1169
LWGEOM ** geoms
Definition: liblwgeom.h:493
int lw_dist2d_pt_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, DISTPTS *dl)
Definition: measures.c:1421
int lw_dist2d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS *dl)
line to polygon calculation Brute force.
Definition: measures.c:719
int lw_dist2d_point_circstring(LWPOINT *point, LWCIRCSTRING *circ, DISTPTS *dl)
Definition: measures.c:584
int lw_dist2d_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B, DISTPTS *dl)
lw_dist2d_comp from p to line A->B This one is now sending every occation to lw_dist2d_pt_pt Before i...
Definition: measures.c:2189
POINTARRAY ** rings
Definition: liblwgeom.h:441
int lw_dist2d_point_line(LWPOINT *point, LWLINE *line, DISTPTS *dl)
point to line calculation
Definition: measures.c:575
int lw_dist2d_line_circstring(LWLINE *line1, LWCIRCSTRING *line2, DISTPTS *dl)
Definition: measures.c:702
int lw_dist2d_pt_ptarrayarc(const POINT2D *p, const POINTARRAY *pa, DISTPTS *dl)
Search all the arcs of pointarray to see which one is closest to p1 Returns minimum distance between ...
Definition: measures.c:1065
int lw_dist2d_seg_seg(const POINT2D *A, const POINT2D *B, const POINT2D *C, const POINT2D *D, DISTPTS *dl)
Finds the shortest distance between two segments.
Definition: measures.c:1802
int lw_dist2d_pt_ptarray(const POINT2D *p, POINTARRAY *pa, DISTPTS *dl)
search all the segments of pointarray to see which one is closest to p1 Returns minimum distance betw...
Definition: measures.c:1035
POINT2D p2
Definition: measures.h:26
int nrings
Definition: liblwgeom.h:439
double ymax
Definition: liblwgeom.h:279
char * s
Definition: cu_in_wkt.c:23
double y
Definition: liblwgeom.h:312
double lwgeom_maxdistance2d(const LWGEOM *lw1, const LWGEOM *lw2)
Function initialazing max distance calculation.
Definition: measures.c:155
int lw_dist2d_poly_curvepoly(LWPOLY *poly1, LWCURVEPOLY *curvepoly2, DISTPTS *dl)
Definition: measures.c:904
tuple x
Definition: pixval.py:53
int twisted
Definition: measures.h:28
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:75
int lw_dist2d_distribute_bruteforce(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS *dl)
Definition: measures.c:356
int lw_pt_in_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
Returns true if P is on the same side of the plane partition defined by A1/A3 as A2 is...
Definition: lwalgorithm.c:71
int lw_dist2d_selected_seg_seg(const POINT2D *A, const POINT2D *B, const POINT2D *C, const POINT2D *D, DISTPTS *dl)
This is the same function as lw_dist2d_seg_seg but without any calculations to determine intersection...
Definition: measures.c:2143
#define MULTISURFACETYPE
Definition: liblwgeom.h:81
double distance
Definition: measures.h:24
#define DIST_MAX
Definition: measures.h:16
int lw_dist2d_check_overlap(LWGEOM *lwg1, LWGEOM *lwg2)
We have to check for overlapping bboxes.
Definition: measures.c:484
int lw_dist2d_line_line(LWLINE *line1, LWLINE *line2, DISTPTS *dl)
line to line calculation
Definition: measures.c:693
LWCOLLECTION * lwgeom_as_lwcollection(const LWGEOM *lwgeom)
Definition: lwgeom.c:143
double distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2)
The old function nessecary for ptarray_segmentize2d in ptarray.c.
Definition: measures.c:2301
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:70
int lw_dist2d_ptarrayarc_ptarrayarc(const POINTARRAY *pa, const POINTARRAY *pb, DISTPTS *dl)
Test each arc of pa against each arc of pb for distance.
Definition: measures.c:1223
void lwgeom_add_bbox(LWGEOM *lwgeom)
Compute a bbox if not already computed.
Definition: lwgeom.c:599
LWGEOM * lwgeom_furthest_point(const LWGEOM *lw1, const LWGEOM *lw2)
Definition: measures.c:46
#define FP_ABS(a)
static const POINT2D * lw_curvering_getfirstpoint2d_cp(LWGEOM *geom)
Definition: measures.c:935
int lw_pt_in_seg(const POINT2D *P, const POINT2D *A1, const POINT2D *A2)
Returns true if P is between A1/A2.
Definition: lwalgorithm.c:81
uint8_t type
Definition: liblwgeom.h:380
int lw_dist2d_seg_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl)
Calculate the shortest distance between an arc and an edge.
Definition: measures.c:1274
#define FP_EQUALS(A, B)
int ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt)
Return 1 if the point is inside the POINTARRAY, -1 if it is outside, and 0 if it is on the boundary...
Definition: ptarray.c:733
POINTARRAY * points
Definition: liblwgeom.h:428
#define CIRCSTRINGTYPE
Definition: liblwgeom.h:77
double lwgeom_mindistance2d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
Function handling min distance calculations and dwithin calculations.
Definition: measures.c:199
#define LW_OUTSIDE
void * lwalloc(size_t size)
Definition: lwutil.c:199
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members) ...
Definition: lwgeom.c:1297
int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
lw_segment_side()
Definition: lwalgorithm.c:50
Structure used in distance-calculations.
Definition: measures.h:22
void lw_dist2d_distpts_init(DISTPTS *dl, int mode)
Definition: measures.c:53
int lwgeom_contains_point(const LWGEOM *geom, const POINT2D *pt)
Definition: lwcompound.c:116
#define MULTILINETYPE
Definition: liblwgeom.h:74
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:81
int lw_dist2d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2, DISTPTS *dl)
test each segment of l1 against each segment of l2.
Definition: measures.c:1117
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:55
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:74
tuple y
Definition: pixval.py:54
int lw_dist2d_pt_pt(const POINT2D *thep1, const POINT2D *thep2, DISTPTS *dl)
Compares incomming points and stores the points closest to each other or most far away from each othe...
Definition: measures.c:2265
int lw_dist2d_distribute_fast(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS *dl)
Here the geometries are distributed for the new faster distance-calculations.
Definition: measures.c:508
#define COLLECTIONTYPE
Definition: liblwgeom.h:76
POINTARRAY * points
Definition: liblwgeom.h:406
int pnr
Definition: measures.h:35