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