PostGIS  2.5.0beta2dev-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
2065  (ia->themeasure > ib->themeasure) ? 1 : ((ia->themeasure < ib->themeasure) ? -1 : 0);
2066 }
2067 
2071 int
2073 {
2074  const POINT2D *p1, *p2, *p3, *p4, *p01, *p02;
2075  int pnr1,pnr2,pnr3,pnr4, n1, n2, i, u, r, twist;
2076  double maxmeasure;
2077  n1= l1->npoints;
2078  n2 = l2->npoints;
2079 
2080  LWDEBUG(2, "lw_dist2d_pre_seg_seg is called");
2081 
2082  p1 = getPoint2d_cp(l1, list1[0].pnr);
2083  p3 = getPoint2d_cp(l2, list2[0].pnr);
2084  lw_dist2d_pt_pt(p1, p3, dl);
2085  maxmeasure = sqrt(dl->distance*dl->distance + (dl->distance*dl->distance*k*k));
2086  twist = dl->twisted; /*to keep the incoming order between iterations*/
2087  for (i =(n1-1); i>=0; --i)
2088  {
2089  /*we break this iteration when we have checked every
2090  point closer to our perpendicular "checkline" than
2091  our shortest found distance*/
2092  if (((list2[0].themeasure-list1[i].themeasure)) > maxmeasure) break;
2093  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*/
2094  {
2095  pnr1 = list1[i].pnr;
2096  p1 = getPoint2d_cp(l1, pnr1);
2097  if (pnr1+r<0)
2098  {
2099  p01 = getPoint2d_cp(l1, (n1-1));
2100  if (( p1->x == p01->x) && (p1->y == p01->y)) pnr2 = (n1-1);
2101  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*/
2102  }
2103 
2104  else if (pnr1+r>(n1-1))
2105  {
2106  p01 = getPoint2d_cp(l1, 0);
2107  if (( p1->x == p01->x) && (p1->y == p01->y)) pnr2 = 0;
2108  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*/
2109  }
2110  else pnr2 = pnr1+r;
2111 
2112 
2113  p2 = getPoint2d_cp(l1, pnr2);
2114  for (u=0; u<n2; ++u)
2115  {
2116  if (((list2[u].themeasure-list1[i].themeasure)) >= maxmeasure) break;
2117  pnr3 = list2[u].pnr;
2118  p3 = getPoint2d_cp(l2, pnr3);
2119  if (pnr3==0)
2120  {
2121  p02 = getPoint2d_cp(l2, (n2-1));
2122  if (( p3->x == p02->x) && (p3->y == p02->y)) pnr4 = (n2-1);
2123  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*/
2124  }
2125  else pnr4 = pnr3-1;
2126 
2127  p4 = getPoint2d_cp(l2, pnr4);
2128  dl->twisted=twist;
2129  if (!lw_dist2d_selected_seg_seg(p1, p2, p3, p4, dl)) return LW_FALSE;
2130 
2131  if (pnr3>=(n2-1))
2132  {
2133  p02 = getPoint2d_cp(l2, 0);
2134  if (( p3->x == p02->x) && (p3->y == p02->y)) pnr4 = 0;
2135  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*/
2136  }
2137 
2138  else pnr4 = pnr3+1;
2139 
2140  p4 = getPoint2d_cp(l2, pnr4);
2141  dl->twisted=twist; /*we reset the "twist" for each iteration*/
2142  if (!lw_dist2d_selected_seg_seg(p1, p2, p3, p4, dl)) return LW_FALSE;
2143 
2144  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*/
2145  }
2146  }
2147  }
2148 
2149  return LW_TRUE;
2150 }
2151 
2152 
2158 int
2159 lw_dist2d_selected_seg_seg(const POINT2D *A, const POINT2D *B, const POINT2D *C, const POINT2D *D, DISTPTS *dl)
2160 {
2161  LWDEBUGF(2, "lw_dist2d_selected_seg_seg [%g,%g]->[%g,%g] by [%g,%g]->[%g,%g]",
2162  A->x,A->y,B->x,B->y, C->x,C->y, D->x, D->y);
2163 
2164  /*A and B are the same point */
2165  if ( ( A->x == B->x) && (A->y == B->y) )
2166  {
2167  return lw_dist2d_pt_seg(A,C,D,dl);
2168  }
2169  /*U and V are the same point */
2170 
2171  if ( ( C->x == D->x) && (C->y == D->y) )
2172  {
2173  dl->twisted= ((dl->twisted) * (-1));
2174  return lw_dist2d_pt_seg(D,A,B,dl);
2175  }
2176 
2177  if ((lw_dist2d_pt_seg(A,C,D,dl)) && (lw_dist2d_pt_seg(B,C,D,dl)))
2178  {
2179  dl->twisted= ((dl->twisted) * (-1)); /*here we change the order of inputted geometrys and that we notice by changing sign on dl->twisted*/
2180  return ((lw_dist2d_pt_seg(C,A,B,dl)) && (lw_dist2d_pt_seg(D,A,B,dl))); /*if all is successful we return true*/
2181  }
2182  else
2183  {
2184  return LW_FALSE; /* if any of the calls to lw_dist2d_pt_seg goes wrong we return false*/
2185  }
2186 }
2187 
2188 /*------------------------------------------------------------------------------------------------------------
2189 End of New faster distance calculations
2190 --------------------------------------------------------------------------------------------------------------*/
2191 
2192 
2193 /*------------------------------------------------------------------------------------------------------------
2194 Functions in common for Brute force and new calculation
2195 --------------------------------------------------------------------------------------------------------------*/
2196 
2204 int
2205 lw_dist2d_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B, DISTPTS *dl)
2206 {
2207  POINT2D c;
2208  double r;
2209  /*if start==end, then use pt distance */
2210  if ( ( A->x == B->x) && (A->y == B->y) )
2211  {
2212  return lw_dist2d_pt_pt(p,A,dl);
2213  }
2214  /*
2215  * otherwise, we use comp.graphics.algorithms
2216  * Frequently Asked Questions method
2217  *
2218  * (1) AC dot AB
2219  * r = ---------
2220  * ||AB||^2
2221  * r has the following meaning:
2222  * r=0 P = A
2223  * r=1 P = B
2224  * r<0 P is on the backward extension of AB
2225  * r>1 P is on the forward extension of AB
2226  * 0<r<1 P is interior to AB
2227  */
2228 
2229  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) );
2230 
2231  /*This is for finding the maxdistance.
2232  the maxdistance have to be between two vertexes,
2233  compared to mindistance which can be between
2234  two vertexes vertex.*/
2235  if (dl->mode == DIST_MAX)
2236  {
2237  if (r>=0.5)
2238  {
2239  return lw_dist2d_pt_pt(p,A,dl);
2240  }
2241  if (r<0.5)
2242  {
2243  return lw_dist2d_pt_pt(p,B,dl);
2244  }
2245  }
2246 
2247  if (r<0) /*If p projected on the line is outside point A*/
2248  {
2249  return lw_dist2d_pt_pt(p,A,dl);
2250  }
2251  if (r>=1) /*If p projected on the line is outside point B or on point B*/
2252  {
2253  return lw_dist2d_pt_pt(p,B,dl);
2254  }
2255 
2256  /*If the point p is on the segment this is a more robust way to find out that*/
2257  if (( ((A->y-p->y)*(B->x-A->x)==(A->x-p->x)*(B->y-A->y) ) ) && (dl->mode == DIST_MIN))
2258  {
2259  dl->distance = 0.0;
2260  dl->p1 = *p;
2261  dl->p2 = *p;
2262  }
2263 
2264  /*If the projection of point p on the segment is between A and B
2265  then we find that "point on segment" and send it to lw_dist2d_pt_pt*/
2266  c.x=A->x + r * (B->x-A->x);
2267  c.y=A->y + r * (B->y-A->y);
2268 
2269  return lw_dist2d_pt_pt(p,&c,dl);
2270 }
2271 
2272 
2280 int
2281 lw_dist2d_pt_pt(const POINT2D *thep1, const POINT2D *thep2, DISTPTS *dl)
2282 {
2283  double hside = thep2->x - thep1->x;
2284  double vside = thep2->y - thep1->y;
2285  double dist = sqrt ( hside*hside + vside*vside );
2286 
2287  if (((dl->distance - dist)*(dl->mode))>0) /*multiplication with mode to handle mindistance (mode=1) and maxdistance (mode = (-1)*/
2288  {
2289  dl->distance = dist;
2290 
2291  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*/
2292  {
2293  dl->p1 = *thep1;
2294  dl->p2 = *thep2;
2295  }
2296  else
2297  {
2298  dl->p1 = *thep2;
2299  dl->p2 = *thep1;
2300  }
2301  }
2302  return LW_TRUE;
2303 }
2304 
2305 
2306 
2307 
2308 /*------------------------------------------------------------------------------------------------------------
2309 End of Functions in common for Brute force and new calculation
2310 --------------------------------------------------------------------------------------------------------------*/
2311 
2312 inline double
2313 distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2)
2314 {
2315  double hside = p2->x - p1->x;
2316  double vside = p2->y - p1->y;
2317 
2318  return hypot(hside, vside);
2319 }
2320 
2321 inline double
2322 distance2d_sqr_pt_pt(const POINT2D *p1, const POINT2D *p2)
2323 {
2324  double hside = p2->x - p1->x;
2325  double vside = p2->y - p1->y;
2326 
2327  return hside * hside + vside * vside;
2328 }
2329 
2330 double
2331 distance2d_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B)
2332 {
2333  double r,s;
2334 
2335  /*if start==end, then use pt distance */
2336  if ( ( A->x == B->x) && (A->y == B->y) )
2337  return distance2d_pt_pt(p,A);
2338 
2339  /*
2340  * otherwise, we use comp.graphics.algorithms
2341  * Frequently Asked Questions method
2342  *
2343  * (1) AC dot AB
2344  * r = ---------
2345  * ||AB||^2
2346  * r has the following meaning:
2347  * r=0 P = A
2348  * r=1 P = B
2349  * r<0 P is on the backward extension of AB
2350  * r>1 P is on the forward extension of AB
2351  * 0<r<1 P is interior to AB
2352  */
2353 
2354  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) );
2355 
2356  if (r<0) return distance2d_pt_pt(p,A);
2357  if (r>1) return distance2d_pt_pt(p,B);
2358 
2359 
2360  /*
2361  * (2)
2362  * (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
2363  * s = -----------------------------
2364  * L^2
2365  *
2366  * Then the distance from C to P = |s|*L.
2367  *
2368  */
2369 
2370  s = ( (A->y-p->y)*(B->x-A->x)- (A->x-p->x)*(B->y-A->y) ) /
2371  ( (B->x-A->x)*(B->x-A->x) +(B->y-A->y)*(B->y-A->y) );
2372 
2373  return FP_ABS(s) * sqrt(
2374  (B->x-A->x)*(B->x-A->x) + (B->y-A->y)*(B->y-A->y)
2375  );
2376 }
2377 
2378 /* return distance squared, useful to avoid sqrt calculations */
2379 double
2380 distance2d_sqr_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B)
2381 {
2382  double r,s;
2383 
2384  if ( ( A->x == B->x) && (A->y == B->y) )
2385  return distance2d_sqr_pt_pt(p,A);
2386 
2387  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) );
2388 
2389  if (r<0) return distance2d_sqr_pt_pt(p,A);
2390  if (r>1) return distance2d_sqr_pt_pt(p,B);
2391 
2392 
2393  /*
2394  * (2)
2395  * (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
2396  * s = -----------------------------
2397  * L^2
2398  *
2399  * Then the distance from C to P = |s|*L.
2400  *
2401  */
2402 
2403  s = ( (A->y-p->y)*(B->x-A->x)- (A->x-p->x)*(B->y-A->y) ) /
2404  ( (B->x-A->x)*(B->x-A->x) +(B->y-A->y)*(B->y-A->y) );
2405 
2406  return s * s * ( (B->x-A->x)*(B->x-A->x) + (B->y-A->y)*(B->y-A->y) );
2407 }
2408 
2409 
2410 
2415 int
2416 azimuth_pt_pt(const POINT2D *A, const POINT2D *B, double *d)
2417 {
2418  if ( A->x == B->x )
2419  {
2420  if ( A->y < B->y ) *d=0.0;
2421  else if ( A->y > B->y ) *d=M_PI;
2422  else return 0;
2423  return 1;
2424  }
2425 
2426  if ( A->y == B->y )
2427  {
2428  if ( A->x < B->x ) *d=M_PI/2;
2429  else if ( A->x > B->x ) *d=M_PI+(M_PI/2);
2430  else return 0;
2431  return 1;
2432  }
2433 
2434  if ( A->x < B->x )
2435  {
2436  if ( A->y < B->y )
2437  {
2438  *d=atan(fabs(A->x - B->x) / fabs(A->y - B->y) );
2439  }
2440  else /* ( A->y > B->y ) - equality case handled above */
2441  {
2442  *d=atan(fabs(A->y - B->y) / fabs(A->x - B->x) )
2443  + (M_PI/2);
2444  }
2445  }
2446 
2447  else /* ( A->x > B->x ) - equality case handled above */
2448  {
2449  if ( A->y > B->y )
2450  {
2451  *d=atan(fabs(A->x - B->x) / fabs(A->y - B->y) )
2452  + M_PI;
2453  }
2454  else /* ( A->y < B->y ) - equality case handled above */
2455  {
2456  *d=atan(fabs(A->y - B->y) / fabs(A->x - B->x) )
2457  + (M_PI+(M_PI/2));
2458  }
2459  }
2460 
2461  return 1;
2462 }
2463 
2464 
double distance2d_sqr_pt_pt(const POINT2D *p1, const POINT2D *p2)
Definition: measures.c:2322
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:228
#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)
Definition: measures.c:2331
double distance2d_sqr_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B)
Definition: measures.c:2380
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:2072
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:2416
#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:2205
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:2159
#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)
Definition: measures.c:2313
#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:2281
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