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