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