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