PostGIS  2.5.7dev-r@@SVN_REVISION@@
lwgeodetic_tree.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 (C) 2012-2015 Paul Ramsey <pramsey@cleverelephant.ca>
22  * Copyright (C) 2012-2015 Sandro Santilli <strk@kbt.io>
23  *
24  **********************************************************************/
25 
26 #include "liblwgeom_internal.h"
27 #include "lwgeodetic_tree.h"
28 #include "lwgeom_log.h"
29 
30 
31 /* Internal prototype */
32 static CIRC_NODE* circ_nodes_merge(CIRC_NODE** nodes, int num_nodes);
33 static double circ_tree_distance_tree_internal(const CIRC_NODE* n1, const CIRC_NODE* n2, double threshold, double* min_dist, double* max_dist, GEOGRAPHIC_POINT* closest1, GEOGRAPHIC_POINT* closest2);
34 
35 
39 static inline int
41 {
42  return (node->num_nodes == 0);
43 }
44 
49 void
51 {
52  uint32_t i;
53  if ( ! node ) return;
54 
55  if (node->nodes)
56  {
57  for (i = 0; i < node->num_nodes; i++)
58  circ_tree_free(node->nodes[i]);
59  lwfree(node->nodes);
60  }
61  lwfree(node);
62 }
63 
64 
68 static CIRC_NODE*
69 circ_node_leaf_new(const POINTARRAY* pa, int i)
70 {
71  POINT2D *p1, *p2;
72  POINT3D q1, q2, c;
73  GEOGRAPHIC_POINT g1, g2, gc;
74  CIRC_NODE *node;
75  double diameter;
76 
77  p1 = (POINT2D*)getPoint_internal(pa, i);
78  p2 = (POINT2D*)getPoint_internal(pa, i+1);
79  geographic_point_init(p1->x, p1->y, &g1);
80  geographic_point_init(p2->x, p2->y, &g2);
81 
82  LWDEBUGF(3,"edge #%d (%g %g, %g %g)", i, p1->x, p1->y, p2->x, p2->y);
83 
84  diameter = sphere_distance(&g1, &g2);
85 
86  /* Zero length edge, doesn't get a node */
87  if ( FP_EQUALS(diameter, 0.0) )
88  return NULL;
89 
90  /* Allocate */
91  node = lwalloc(sizeof(CIRC_NODE));
92  node->p1 = p1;
93  node->p2 = p2;
94 
95  /* Convert ends to X/Y/Z, sum, and normalize to get mid-point */
96  geog2cart(&g1, &q1);
97  geog2cart(&g2, &q2);
98  vector_sum(&q1, &q2, &c);
99  normalize(&c);
100  cart2geog(&c, &gc);
101  node->center = gc;
102  node->radius = diameter / 2.0;
103 
104  LWDEBUGF(3,"edge #%d CENTER(%g %g) RADIUS=%g", i, gc.lon, gc.lat, node->radius);
105 
106  /* Leaf has no children */
107  node->num_nodes = 0;
108  node->nodes = NULL;
109  node->edge_num = i;
110 
111  /* Zero out metadata */
112  node->pt_outside.x = 0.0;
113  node->pt_outside.y = 0.0;
114  node->geom_type = 0;
115 
116  return node;
117 }
118 
122 static CIRC_NODE*
124 {
125  CIRC_NODE* tree = lwalloc(sizeof(CIRC_NODE));
126  tree->p1 = tree->p2 = (POINT2D*)getPoint_internal(pa, 0);
127  geographic_point_init(tree->p1->x, tree->p1->y, &(tree->center));
128  tree->radius = 0.0;
129  tree->nodes = NULL;
130  tree->num_nodes = 0;
131  tree->edge_num = 0;
132  tree->geom_type = POINTTYPE;
133  tree->pt_outside.x = 0.0;
134  tree->pt_outside.y = 0.0;
135  return tree;
136 }
137 
142 static int
143 circ_node_compare(const void* v1, const void* v2)
144 {
145  POINT2D p1, p2;
146  unsigned int u1, u2;
147  CIRC_NODE *c1 = *((CIRC_NODE**)v1);
148  CIRC_NODE *c2 = *((CIRC_NODE**)v2);
149  p1.x = rad2deg((c1->center).lon);
150  p1.y = rad2deg((c1->center).lat);
151  p2.x = rad2deg((c2->center).lon);
152  p2.y = rad2deg((c2->center).lat);
153  u1 = geohash_point_as_int(&p1);
154  u2 = geohash_point_as_int(&p2);
155  if ( u1 < u2 ) return -1;
156  if ( u1 > u2 ) return 1;
157  return 0;
158 }
159 
166 static int
167 circ_center_spherical(const GEOGRAPHIC_POINT* c1, const GEOGRAPHIC_POINT* c2, double distance, double offset, GEOGRAPHIC_POINT* center)
168 {
169  /* Direction from c1 to c2 */
170  double dir = sphere_direction(c1, c2, distance);
171 
172  LWDEBUGF(4,"calculating spherical center", dir);
173 
174  LWDEBUGF(4,"dir is %g", dir);
175 
176  /* Catch sphere_direction when it barfs */
177  if ( isnan(dir) )
178  return LW_FAILURE;
179 
180  /* Center of new circle is projection from start point, using offset distance*/
181  return sphere_project(c1, offset, dir, center);
182 }
183 
191 static int
192 circ_center_cartesian(const GEOGRAPHIC_POINT* c1, const GEOGRAPHIC_POINT* c2, double distance, double offset, GEOGRAPHIC_POINT* center)
193 {
194  POINT3D p1, p2;
195  POINT3D p1p2, pc;
196  double proportion = offset/distance;
197 
198  LWDEBUG(4,"calculating cartesian center");
199 
200  geog2cart(c1, &p1);
201  geog2cart(c2, &p2);
202 
203  /* Difference between p2 and p1 */
204  p1p2.x = p2.x - p1.x;
205  p1p2.y = p2.y - p1.y;
206  p1p2.z = p2.z - p1.z;
207 
208  /* Scale difference to proportion */
209  p1p2.x *= proportion;
210  p1p2.y *= proportion;
211  p1p2.z *= proportion;
212 
213  /* Add difference to p1 to get approximate center point */
214  pc.x = p1.x + p1p2.x;
215  pc.y = p1.y + p1p2.y;
216  pc.z = p1.z + p1p2.z;
217  normalize(&pc);
218 
219  /* Convert center point to geographics */
220  cart2geog(&pc, center);
221 
222  return LW_SUCCESS;
223 }
224 
225 
230 static CIRC_NODE*
232 {
233  CIRC_NODE *node = NULL;
234  GEOGRAPHIC_POINT new_center, c1;
235  double new_radius;
236  double offset1, dist, D, r1, ri;
237  uint32_t i, new_geom_type;
238 
239  LWDEBUGF(3, "called with %d nodes --", num_nodes);
240 
241  /* Can't do anything w/ empty input */
242  if ( num_nodes < 1 )
243  return node;
244 
245  /* Initialize calculation with values of the first circle */
246  new_center = c[0]->center;
247  new_radius = c[0]->radius;
248  new_geom_type = c[0]->geom_type;
249 
250  /* Merge each remaining circle into the new circle */
251  for ( i = 1; i < num_nodes; i++ )
252  {
253  c1 = new_center;
254  r1 = new_radius;
255 
256  dist = sphere_distance(&c1, &(c[i]->center));
257  ri = c[i]->radius;
258 
259  /* Promote geometry types up the tree, getting more and more collected */
260  /* Go until we find a value */
261  if ( ! new_geom_type )
262  {
263  new_geom_type = c[i]->geom_type;
264  }
265  /* Promote singleton to a multi-type */
266  else if ( ! lwtype_is_collection(new_geom_type) )
267  {
268  /* Anonymous collection if types differ */
269  if ( new_geom_type != c[i]->geom_type )
270  {
271  new_geom_type = COLLECTIONTYPE;
272  }
273  else
274  {
275  new_geom_type = lwtype_get_collectiontype(new_geom_type);
276  }
277  }
278  /* If we can't add next feature to this collection cleanly, promote again to anonymous collection */
279  else if ( new_geom_type != lwtype_get_collectiontype(c[i]->geom_type) )
280  {
281  new_geom_type = COLLECTIONTYPE;
282  }
283 
284 
285  LWDEBUGF(3, "distance between new (%g %g) and %i (%g %g) is %g", c1.lon, c1.lat, i, c[i]->center.lon, c[i]->center.lat, dist);
286 
287  if ( FP_EQUALS(dist, 0) )
288  {
289  LWDEBUG(3, " distance between centers is zero");
290  new_radius = r1 + 2*dist;
291  new_center = c1;
292  }
293  else if ( dist < fabs(r1 - ri) )
294  {
295  /* new contains next */
296  if ( r1 > ri )
297  {
298  LWDEBUG(3, " c1 contains ci");
299  new_center = c1;
300  new_radius = r1;
301  }
302  /* next contains new */
303  else
304  {
305  LWDEBUG(3, " ci contains c1");
306  new_center = c[i]->center;
307  new_radius = ri;
308  }
309  }
310  else
311  {
312  LWDEBUG(3, " calculating new center");
313  /* New circle diameter */
314  D = dist + r1 + ri;
315  LWDEBUGF(3," D is %g", D);
316 
317  /* New radius */
318  new_radius = D / 2.0;
319 
320  /* Distance from cn1 center to the new center */
321  offset1 = ri + (D - (2.0*r1 + 2.0*ri)) / 2.0;
322  LWDEBUGF(3," offset1 is %g", offset1);
323 
324  /* Sometimes the sphere_direction function fails... this causes the center calculation */
325  /* to fail too. In that case, we're going to fall back to a cartesian calculation, which */
326  /* is less exact, so we also have to pad the radius by (hack alert) an arbitrary amount */
327  /* which is hopefully always big enough to contain the input edges */
328  if ( circ_center_spherical(&c1, &(c[i]->center), dist, offset1, &new_center) == LW_FAILURE )
329  {
330  circ_center_cartesian(&c1, &(c[i]->center), dist, offset1, &new_center);
331  new_radius *= 1.1;
332  }
333  }
334  LWDEBUGF(3, " new center is (%g %g) new radius is %g", new_center.lon, new_center.lat, new_radius);
335  }
336 
337  node = lwalloc(sizeof(CIRC_NODE));
338  node->p1 = NULL;
339  node->p2 = NULL;
340  node->center = new_center;
341  node->radius = new_radius;
342  node->num_nodes = num_nodes;
343  node->nodes = c;
344  node->edge_num = -1;
345  node->geom_type = new_geom_type;
346  node->pt_outside.x = 0.0;
347  node->pt_outside.y = 0.0;
348  return node;
349 }
350 
354 CIRC_NODE*
356 {
357  int num_edges;
358  int i, j;
359  CIRC_NODE **nodes;
360  CIRC_NODE *node;
361  CIRC_NODE *tree;
362 
363  /* Can't do anything with no points */
364  if ( pa->npoints < 1 )
365  return NULL;
366 
367  /* Special handling for a single point */
368  if ( pa->npoints == 1 )
369  return circ_node_leaf_point_new(pa);
370 
371  /* First create a flat list of nodes, one per edge. */
372  num_edges = pa->npoints - 1;
373  nodes = lwalloc(sizeof(CIRC_NODE*) * pa->npoints);
374  j = 0;
375  for ( i = 0; i < num_edges; i++ )
376  {
377  node = circ_node_leaf_new(pa, i);
378  if ( node ) /* Not zero length? */
379  nodes[j++] = node;
380  }
381 
382  /* Special case: only zero-length edges. Make a point node. */
383  if ( j == 0 ) {
384  lwfree(nodes);
385  return circ_node_leaf_point_new(pa);
386  }
387 
388  /* Merge the node list pairwise up into a tree */
389  tree = circ_nodes_merge(nodes, j);
390 
391  /* Free the old list structure, leaving the tree in place */
392  lwfree(nodes);
393 
394  return tree;
395 }
396 
402 static void
403 circ_nodes_sort(CIRC_NODE** nodes, int num_nodes)
404 {
405  qsort(nodes, num_nodes, sizeof(CIRC_NODE*), circ_node_compare);
406 }
407 
408 
409 static CIRC_NODE*
410 circ_nodes_merge(CIRC_NODE** nodes, int num_nodes)
411 {
412  CIRC_NODE **inodes = NULL;
413  int num_children = num_nodes;
414  int inode_num = 0;
415  int num_parents = 0;
416  int j;
417 
418  /* TODO, roll geom_type *up* as tree is built, changing to collection types as simple types are merged
419  * TODO, change the distance algorithm to drive down to simple types first, test pip on poly/other cases, then test edges
420  */
421 
422  while( num_children > 1 )
423  {
424  for ( j = 0; j < num_children; j++ )
425  {
426  inode_num = (j % CIRC_NODE_SIZE);
427  if ( inode_num == 0 )
428  inodes = lwalloc(sizeof(CIRC_NODE*)*CIRC_NODE_SIZE);
429 
430  inodes[inode_num] = nodes[j];
431 
432  if ( inode_num == CIRC_NODE_SIZE-1 )
433  nodes[num_parents++] = circ_node_internal_new(inodes, CIRC_NODE_SIZE);
434  }
435 
436  /* Clean up any remaining nodes... */
437  if ( inode_num == 0 )
438  {
439  /* Promote solo nodes without merging */
440  nodes[num_parents++] = inodes[0];
441  lwfree(inodes);
442  }
443  else if ( inode_num < CIRC_NODE_SIZE-1 )
444  {
445  /* Merge spare nodes */
446  nodes[num_parents++] = circ_node_internal_new(inodes, inode_num+1);
447  }
448 
449  num_children = num_parents;
450  num_parents = 0;
451  }
452 
453  /* Return a reference to the head of the tree */
454  return nodes[0];
455 }
456 
457 
461 int circ_tree_get_point(const CIRC_NODE* node, POINT2D* pt)
462 {
463  if ( circ_node_is_leaf(node) )
464  {
465  pt->x = node->p1->x;
466  pt->y = node->p1->y;
467  return LW_SUCCESS;
468  }
469  else
470  {
471  return circ_tree_get_point(node->nodes[0], pt);
472  }
473 }
474 
476 {
477  POINT3D center3d;
479  // if (node->radius >= M_PI) return LW_FAILURE;
480  geog2cart(&(node->center), &center3d);
481  vector_scale(&center3d, -1.0);
482  cart2geog(&center3d, &g);
483  pt->x = rad2deg(g.lon);
484  pt->y = rad2deg(g.lat);
485  return LW_SUCCESS;
486 }
487 
488 
495 int circ_tree_contains_point(const CIRC_NODE* node, const POINT2D* pt, const POINT2D* pt_outside, int* on_boundary)
496 {
497  GEOGRAPHIC_POINT closest;
498  GEOGRAPHIC_EDGE stab_edge, edge;
499  POINT3D S1, S2, E1, E2;
500  double d;
501  uint32_t i, c;
502 
503  /* Construct a stabline edge from our "inside" to our known outside point */
504  geographic_point_init(pt->x, pt->y, &(stab_edge.start));
505  geographic_point_init(pt_outside->x, pt_outside->y, &(stab_edge.end));
506  geog2cart(&(stab_edge.start), &S1);
507  geog2cart(&(stab_edge.end), &S2);
508 
509  LWDEBUG(3, "entered");
510 
511  /*
512  * If the stabline doesn't cross within the radius of a node, there's no
513  * way it can cross.
514  */
515 
516  LWDEBUGF(3, "working on node %p, edge_num %d, radius %g, center POINT(%g %g)", node, node->edge_num, node->radius, rad2deg(node->center.lon), rad2deg(node->center.lat));
517  d = edge_distance_to_point(&stab_edge, &(node->center), &closest);
518  LWDEBUGF(3, "edge_distance_to_point=%g, node_radius=%g", d, node->radius);
519  if ( FP_LTEQ(d, node->radius) )
520  {
521  LWDEBUGF(3,"entering this branch (%p)", node);
522 
523  /* Return the crossing number of this leaf */
524  if ( circ_node_is_leaf(node) )
525  {
526  int inter;
527  LWDEBUGF(3, "leaf node calculation (edge %d)", node->edge_num);
528  geographic_point_init(node->p1->x, node->p1->y, &(edge.start));
529  geographic_point_init(node->p2->x, node->p2->y, &(edge.end));
530  geog2cart(&(edge.start), &E1);
531  geog2cart(&(edge.end), &E2);
532 
533  inter = edge_intersects(&S1, &S2, &E1, &E2);
534 
535  if ( inter & PIR_INTERSECTS )
536  {
537  LWDEBUG(3," got stab line edge_intersection with this edge!");
538  /* To avoid double counting crossings-at-a-vertex, */
539  /* always ignore crossings at "lower" ends of edges*/
540 
541  if ( inter & PIR_B_TOUCH_RIGHT || inter & PIR_COLINEAR )
542  {
543  LWDEBUG(3," rejecting stab line grazing by left-side edge");
544  return 0;
545  }
546  else
547  {
548  LWDEBUG(3," accepting stab line intersection");
549  return 1;
550  }
551  }
552  }
553  /* Or, add up the crossing numbers of all children of this node. */
554  else
555  {
556  c = 0;
557  for ( i = 0; i < node->num_nodes; i++ )
558  {
559  LWDEBUG(3,"internal node calculation");
560  LWDEBUGF(3," calling circ_tree_contains_point on child %d!", i);
561  c += circ_tree_contains_point(node->nodes[i], pt, pt_outside, on_boundary);
562  }
563  return c % 2;
564  }
565  }
566  else
567  {
568  LWDEBUGF(3,"skipping this branch (%p)", node);
569  }
570 
571  return 0;
572 }
573 
574 static double
576 {
577  double d = sphere_distance(&(n1->center), &(n2->center));
578  double r1 = n1->radius;
579  double r2 = n2->radius;
580 
581  if ( d < r1 + r2 )
582  return 0.0;
583 
584  return d - r1 - r2;
585 }
586 
587 static double
589 {
590  return sphere_distance(&(n1->center), &(n2->center)) + n1->radius + n2->radius;
591 }
592 
593 double
594 circ_tree_distance_tree(const CIRC_NODE* n1, const CIRC_NODE* n2, const SPHEROID* spheroid, double threshold)
595 {
596  double min_dist = FLT_MAX;
597  double max_dist = FLT_MAX;
598  GEOGRAPHIC_POINT closest1, closest2;
599  /* Quietly decrease the threshold just a little to avoid cases where */
600  /* the actual spheroid distance is larger than the sphere distance */
601  /* causing the return value to be larger than the threshold value */
602  double threshold_radians = 0.95 * threshold / spheroid->radius;
603 
604  circ_tree_distance_tree_internal(n1, n2, threshold_radians, &min_dist, &max_dist, &closest1, &closest2);
605 
606  /* Spherical case */
607  if ( spheroid->a == spheroid->b )
608  {
609  return spheroid->radius * sphere_distance(&closest1, &closest2);
610  }
611  else
612  {
613  return spheroid_distance(&closest1, &closest2, spheroid);
614  }
615 }
616 
617 
618 /***********************************************************************
619 * Internal node sorting routine to make distance calculations faster?
620 */
621 
622 struct sort_node {
624  double d;
625 };
626 
627 static int
628 circ_nodes_sort_cmp(const void *a, const void *b)
629 {
630  struct sort_node *node_a = (struct sort_node *)(a);
631  struct sort_node *node_b = (struct sort_node *)(b);
632  if (node_a->d < node_b->d) return -1;
633  else if (node_a->d > node_b->d) return 1;
634  else return 0;
635 }
636 
637 static void
638 circ_internal_nodes_sort(CIRC_NODE **nodes, uint32_t num_nodes, const CIRC_NODE *target_node)
639 {
640  uint32_t i;
641  struct sort_node sort_nodes[CIRC_NODE_SIZE];
642 
643  /* Copy incoming nodes into sorting array and calculate */
644  /* distance to the target node */
645  for (i = 0; i < num_nodes; i++)
646  {
647  sort_nodes[i].node = nodes[i];
648  sort_nodes[i].d = sphere_distance(&(nodes[i]->center), &(target_node->center));
649  }
650 
651  /* Sort the nodes and copy the result back into the input array */
652  qsort(sort_nodes, num_nodes, sizeof(struct sort_node), circ_nodes_sort_cmp);
653  for (i = 0; i < num_nodes; i++)
654  {
655  nodes[i] = sort_nodes[i].node;
656  }
657  return;
658 }
659 
660 /***********************************************************************/
661 
662 static double
663 circ_tree_distance_tree_internal(const CIRC_NODE* n1, const CIRC_NODE* n2, double threshold, double* min_dist, double* max_dist, GEOGRAPHIC_POINT* closest1, GEOGRAPHIC_POINT* closest2)
664 {
665  double max;
666  double d, d_min;
667  uint32_t i;
668 
669  LWDEBUGF(4, "entered, min_dist=%.8g max_dist=%.8g, type1=%d, type2=%d", *min_dist, *max_dist, n1->geom_type, n2->geom_type);
670 /*
671  circ_tree_print(n1, 0);
672  circ_tree_print(n2, 0);
673 */
674 
675  /* Short circuit if we've already hit the minimum */
676  if( *min_dist < threshold || *min_dist == 0.0 )
677  return *min_dist;
678 
679  /* If your minimum is greater than anyone's maximum, you can't hold the winner */
680  if( circ_node_min_distance(n1, n2) > *max_dist )
681  {
682  LWDEBUGF(4, "pruning pair %p, %p", n1, n2);
683  return FLT_MAX;
684  }
685 
686  /* If your maximum is a new low, we'll use that as our new global tolerance */
687  max = circ_node_max_distance(n1, n2);
688  LWDEBUGF(5, "max %.8g", max);
689  if( max < *max_dist )
690  *max_dist = max;
691 
692  /* Polygon on one side, primitive type on the other. Check for point-in-polygon */
693  /* short circuit. */
694  if ( n1->geom_type == POLYGONTYPE && n2->geom_type && ! lwtype_is_collection(n2->geom_type) )
695  {
696  POINT2D pt;
697  circ_tree_get_point(n2, &pt);
698  LWDEBUGF(4, "n1 is polygon, testing if contains (%.5g,%.5g)", pt.x, pt.y);
699  if ( circ_tree_contains_point(n1, &pt, &(n1->pt_outside), NULL) )
700  {
701  LWDEBUG(4, "it does");
702  *min_dist = 0.0;
703  geographic_point_init(pt.x, pt.y, closest1);
704  geographic_point_init(pt.x, pt.y, closest2);
705  return *min_dist;
706  }
707  }
708  /* Polygon on one side, primitive type on the other. Check for point-in-polygon */
709  /* short circuit. */
710  if ( n2->geom_type == POLYGONTYPE && n1->geom_type && ! lwtype_is_collection(n1->geom_type) )
711  {
712  POINT2D pt;
713  circ_tree_get_point(n1, &pt);
714  LWDEBUGF(4, "n2 is polygon, testing if contains (%.5g,%.5g)", pt.x, pt.y);
715  if ( circ_tree_contains_point(n2, &pt, &(n2->pt_outside), NULL) )
716  {
717  LWDEBUG(4, "it does");
718  geographic_point_init(pt.x, pt.y, closest1);
719  geographic_point_init(pt.x, pt.y, closest2);
720  *min_dist = 0.0;
721  return *min_dist;
722  }
723  }
724 
725  /* Both leaf nodes, do a real distance calculation */
726  if( circ_node_is_leaf(n1) && circ_node_is_leaf(n2) )
727  {
728  double d;
729  GEOGRAPHIC_POINT close1, close2;
730  LWDEBUGF(4, "testing leaf pair [%d], [%d]", n1->edge_num, n2->edge_num);
731  /* One of the nodes is a point */
732  if ( n1->p1 == n1->p2 || n2->p1 == n2->p2 )
733  {
734  GEOGRAPHIC_EDGE e;
735  GEOGRAPHIC_POINT gp1, gp2;
736 
737  /* Both nodes are points! */
738  if ( n1->p1 == n1->p2 && n2->p1 == n2->p2 )
739  {
740  geographic_point_init(n1->p1->x, n1->p1->y, &gp1);
741  geographic_point_init(n2->p1->x, n2->p1->y, &gp2);
742  close1 = gp1; close2 = gp2;
743  d = sphere_distance(&gp1, &gp2);
744  }
745  /* Node 1 is a point */
746  else if ( n1->p1 == n1->p2 )
747  {
748  geographic_point_init(n1->p1->x, n1->p1->y, &gp1);
749  geographic_point_init(n2->p1->x, n2->p1->y, &(e.start));
750  geographic_point_init(n2->p2->x, n2->p2->y, &(e.end));
751  close1 = gp1;
752  d = edge_distance_to_point(&e, &gp1, &close2);
753  }
754  /* Node 2 is a point */
755  else
756  {
757  geographic_point_init(n2->p1->x, n2->p1->y, &gp1);
758  geographic_point_init(n1->p1->x, n1->p1->y, &(e.start));
759  geographic_point_init(n1->p2->x, n1->p2->y, &(e.end));
760  close1 = gp1;
761  d = edge_distance_to_point(&e, &gp1, &close2);
762  }
763  LWDEBUGF(4, " got distance %g", d);
764  }
765  /* Both nodes are edges */
766  else
767  {
768  GEOGRAPHIC_EDGE e1, e2;
770  POINT3D A1, A2, B1, B2;
771  geographic_point_init(n1->p1->x, n1->p1->y, &(e1.start));
772  geographic_point_init(n1->p2->x, n1->p2->y, &(e1.end));
773  geographic_point_init(n2->p1->x, n2->p1->y, &(e2.start));
774  geographic_point_init(n2->p2->x, n2->p2->y, &(e2.end));
775  geog2cart(&(e1.start), &A1);
776  geog2cart(&(e1.end), &A2);
777  geog2cart(&(e2.start), &B1);
778  geog2cart(&(e2.end), &B2);
779  if ( edge_intersects(&A1, &A2, &B1, &B2) )
780  {
781  d = 0.0;
782  edge_intersection(&e1, &e2, &g);
783  close1 = close2 = g;
784  }
785  else
786  {
787  d = edge_distance_to_edge(&e1, &e2, &close1, &close2);
788  }
789  LWDEBUGF(4, "edge_distance_to_edge returned %g", d);
790  }
791  if ( d < *min_dist )
792  {
793  *min_dist = d;
794  *closest1 = close1;
795  *closest2 = close2;
796  }
797  return d;
798  }
799  else
800  {
801  d_min = FLT_MAX;
802  /* Drive the recursion into the COLLECTION types first so we end up with */
803  /* pairings of primitive geometries that can be forced into the point-in-polygon */
804  /* tests above. */
805  if ( n1->geom_type && lwtype_is_collection(n1->geom_type) )
806  {
808  for ( i = 0; i < n1->num_nodes; i++ )
809  {
810  d = circ_tree_distance_tree_internal(n1->nodes[i], n2, threshold, min_dist, max_dist, closest1, closest2);
811  d_min = FP_MIN(d_min, d);
812  }
813  }
814  else if ( n2->geom_type && lwtype_is_collection(n2->geom_type) )
815  {
817  for ( i = 0; i < n2->num_nodes; i++ )
818  {
819  d = circ_tree_distance_tree_internal(n1, n2->nodes[i], threshold, min_dist, max_dist, closest1, closest2);
820  d_min = FP_MIN(d_min, d);
821  }
822  }
823  else if ( ! circ_node_is_leaf(n1) )
824  {
826  for ( i = 0; i < n1->num_nodes; i++ )
827  {
828  d = circ_tree_distance_tree_internal(n1->nodes[i], n2, threshold, min_dist, max_dist, closest1, closest2);
829  d_min = FP_MIN(d_min, d);
830  }
831  }
832  else if ( ! circ_node_is_leaf(n2) )
833  {
835  for ( i = 0; i < n2->num_nodes; i++ )
836  {
837  d = circ_tree_distance_tree_internal(n1, n2->nodes[i], threshold, min_dist, max_dist, closest1, closest2);
838  d_min = FP_MIN(d_min, d);
839  }
840  }
841  else
842  {
843  /* Never get here */
844  }
845 
846  return d_min;
847  }
848 }
849 
850 
851 
852 
853 
854 void circ_tree_print(const CIRC_NODE* node, int depth)
855 {
856  uint32_t i;
857 
858  if (circ_node_is_leaf(node))
859  {
860  printf("%*s[%d] C(%.5g %.5g) R(%.5g) ((%.5g %.5g),(%.5g,%.5g))",
861  3*depth + 6, "NODE", node->edge_num,
863  node->radius,
864  node->p1->x, node->p1->y,
865  node->p2->x, node->p2->y
866  );
867  if ( node->geom_type )
868  {
869  printf(" %s", lwtype_name(node->geom_type));
870  }
871  if ( node->geom_type == POLYGONTYPE )
872  {
873  printf(" O(%.5g %.5g)", node->pt_outside.x, node->pt_outside.y);
874  }
875  printf("\n");
876 
877  }
878  else
879  {
880  printf("%*s C(%.5g %.5g) R(%.5g)",
881  3*depth + 6, "NODE",
883  node->radius
884  );
885  if ( node->geom_type )
886  {
887  printf(" %s", lwtype_name(node->geom_type));
888  }
889  if ( node->geom_type == POLYGONTYPE )
890  {
891  printf(" O(%.5g %.5g)", node->pt_outside.x, node->pt_outside.y);
892  }
893  printf("\n");
894  }
895  for ( i = 0; i < node->num_nodes; i++ )
896  {
897  circ_tree_print(node->nodes[i], depth + 1);
898  }
899  return;
900 }
901 
902 
903 static CIRC_NODE*
905 {
906  CIRC_NODE* node;
907  node = circ_tree_new(lwpoint->point);
908  node->geom_type = lwgeom_get_type((LWGEOM*)lwpoint);;
909  return node;
910 }
911 
912 static CIRC_NODE*
914 {
915  CIRC_NODE* node;
916  node = circ_tree_new(lwline->points);
917  node->geom_type = lwgeom_get_type((LWGEOM*)lwline);
918  return node;
919 }
920 
921 static CIRC_NODE*
923 {
924  uint32_t i = 0, j = 0;
925  CIRC_NODE** nodes;
926  CIRC_NODE* node;
927 
928  /* One ring? Handle it like a line. */
929  if ( lwpoly->nrings == 1 )
930  {
931  node = circ_tree_new(lwpoly->rings[0]);
932  }
933  else
934  {
935  /* Calculate a tree for each non-trivial ring of the polygon */
936  nodes = lwalloc(lwpoly->nrings * sizeof(CIRC_NODE*));
937  for ( i = 0; i < lwpoly->nrings; i++ )
938  {
939  node = circ_tree_new(lwpoly->rings[i]);
940  if ( node )
941  nodes[j++] = node;
942  }
943  /* Put the trees into a spatially correlated order */
944  circ_nodes_sort(nodes, j);
945  /* Merge the trees pairwise up to a parent node and return */
946  node = circ_nodes_merge(nodes, j);
947  /* Don't need the working list any more */
948  lwfree(nodes);
949  }
950 
951  /* Metadata about polygons, we need this to apply P-i-P tests */
952  /* selectively when doing distance calculations */
953  node->geom_type = lwgeom_get_type((LWGEOM*)lwpoly);
954  lwpoly_pt_outside(lwpoly, &(node->pt_outside));
955 
956  return node;
957 }
958 
959 static CIRC_NODE*
961 {
962  uint32_t i = 0, j = 0;
963  CIRC_NODE** nodes;
964  CIRC_NODE* node;
965 
966  /* One geometry? Done! */
967  if ( lwcol->ngeoms == 1 )
968  return lwgeom_calculate_circ_tree(lwcol->geoms[0]);
969 
970  /* Calculate a tree for each sub-geometry*/
971  nodes = lwalloc(lwcol->ngeoms * sizeof(CIRC_NODE*));
972  for ( i = 0; i < lwcol->ngeoms; i++ )
973  {
975  if ( node )
976  nodes[j++] = node;
977  }
978  /* Put the trees into a spatially correlated order */
979  circ_nodes_sort(nodes, j);
980  /* Merge the trees pairwise up to a parent node and return */
981  node = circ_nodes_merge(nodes, j);
982  /* Don't need the working list any more */
983  lwfree(nodes);
984  node->geom_type = lwgeom_get_type((LWGEOM*)lwcol);
985  return node;
986 }
987 
988 CIRC_NODE*
990 {
991  if ( lwgeom_is_empty(lwgeom) )
992  return NULL;
993 
994  switch ( lwgeom->type )
995  {
996  case POINTTYPE:
997  return lwpoint_calculate_circ_tree((LWPOINT*)lwgeom);
998  case LINETYPE:
999  return lwline_calculate_circ_tree((LWLINE*)lwgeom);
1000  case POLYGONTYPE:
1001  return lwpoly_calculate_circ_tree((LWPOLY*)lwgeom);
1002  case MULTIPOINTTYPE:
1003  case MULTILINETYPE:
1004  case MULTIPOLYGONTYPE:
1005  case COLLECTIONTYPE:
1007  default:
1008  lwerror("Unable to calculate spherical index tree for type %s", lwtype_name(lwgeom->type));
1009  return NULL;
1010  }
1011 
1012 }
uint32_t lwtype_get_collectiontype(uint8_t type)
Given an lwtype number, what homogeneous collection can hold it?
Definition: lwgeom.c:1120
#define COLLECTIONTYPE
Definition: liblwgeom.h:91
#define LW_FAILURE
Definition: liblwgeom.h:79
#define MULTILINETYPE
Definition: liblwgeom.h:89
#define LINETYPE
Definition: liblwgeom.h:86
uint8_t * getPoint_internal(const POINTARRAY *pa, uint32_t n)
Definition: ptarray.c:1750
#define LW_SUCCESS
Definition: liblwgeom.h:80
unsigned int geohash_point_as_int(POINT2D *pt)
Definition: lwalgorithm.c:657
#define MULTIPOINTTYPE
Definition: liblwgeom.h:88
uint32_t lwgeom_get_type(const LWGEOM *geom)
Return LWTYPE number.
Definition: lwgeom.c:923
int lwtype_is_collection(uint8_t type)
Determine whether a type number is a collection or not.
Definition: lwgeom.c:1093
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:85
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:90
void lwfree(void *mem)
Definition: lwutil.c:244
#define POLYGONTYPE
Definition: liblwgeom.h:87
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:218
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members)
Definition: lwgeom.c:1393
void * lwalloc(size_t size)
Definition: lwutil.c:229
#define FP_MIN(A, B)
#define FP_EQUALS(A, B)
#define FP_LTEQ(A, B)
void normalize(POINT3D *p)
Normalize to a unit vector.
Definition: lwgeodetic.c:615
void vector_scale(POINT3D *n, double scale)
Scale a vector out by a factor.
Definition: lwgeodetic.c:487
void cart2geog(const POINT3D *p, GEOGRAPHIC_POINT *g)
Convert cartesian coordinates on unit sphere to spherical coordinates.
Definition: lwgeodetic.c:414
int lwpoly_pt_outside(const LWPOLY *poly, POINT2D *pt_outside)
Definition: lwgeodetic.c:1527
int sphere_project(const GEOGRAPHIC_POINT *r, double distance, double azimuth, GEOGRAPHIC_POINT *n)
Given a starting location r, a distance and an azimuth to the new point, compute the location of the ...
Definition: lwgeodetic.c:1314
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition: lwgeodetic.c:180
double sphere_distance(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Given two points on a unit sphere, calculate their distance apart in radians.
Definition: lwgeodetic.c:948
double sphere_direction(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e, double d)
Given two points on a unit sphere, calculate the direction from s to e.
Definition: lwgeodetic.c:975
int edge_intersection(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e2, GEOGRAPHIC_POINT *g)
Returns true if an intersection can be calculated, and places it in *g.
Definition: lwgeodetic.c:1127
void vector_sum(const POINT3D *a, const POINT3D *b, POINT3D *n)
Calculate the sum of two vectors.
Definition: lwgeodetic.c:465
uint32_t edge_intersects(const POINT3D *A1, const POINT3D *A2, const POINT3D *B1, const POINT3D *B2)
Returns non-zero if edges A and B interact.
Definition: lwgeodetic.c:3540
double edge_distance_to_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *gp, GEOGRAPHIC_POINT *closest)
Definition: lwgeodetic.c:1218
void geog2cart(const GEOGRAPHIC_POINT *g, POINT3D *p)
Convert spherical coordinates to cartesian coordinates on unit sphere.
Definition: lwgeodetic.c:404
double edge_distance_to_edge(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e2, GEOGRAPHIC_POINT *closest1, GEOGRAPHIC_POINT *closest2)
Calculate the distance between two edges.
Definition: lwgeodetic.c:1269
#define rad2deg(r)
Definition: lwgeodetic.h:80
double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
Computes the shortest distance along the surface of the spheroid between two points.
Definition: lwspheroid.c:191
#define PIR_COLINEAR
Definition: lwgeodetic.h:88
#define PIR_INTERSECTS
Definition: lwgeodetic.h:87
#define PIR_B_TOUCH_RIGHT
Definition: lwgeodetic.h:91
double circ_tree_distance_tree(const CIRC_NODE *n1, const CIRC_NODE *n2, const SPHEROID *spheroid, double threshold)
static CIRC_NODE * lwpoly_calculate_circ_tree(const LWPOLY *lwpoly)
CIRC_NODE * circ_tree_new(const POINTARRAY *pa)
Build a tree of nodes from a point array, one node per edge.
int circ_tree_get_point(const CIRC_NODE *node, POINT2D *pt)
Returns a POINT2D that is a vertex of the input shape.
static CIRC_NODE * lwcollection_calculate_circ_tree(const LWCOLLECTION *lwcol)
static double circ_node_max_distance(const CIRC_NODE *n1, const CIRC_NODE *n2)
static int circ_nodes_sort_cmp(const void *a, const void *b)
static double circ_tree_distance_tree_internal(const CIRC_NODE *n1, const CIRC_NODE *n2, double threshold, double *min_dist, double *max_dist, GEOGRAPHIC_POINT *closest1, GEOGRAPHIC_POINT *closest2)
static CIRC_NODE * lwpoint_calculate_circ_tree(const LWPOINT *lwpoint)
static CIRC_NODE * lwline_calculate_circ_tree(const LWLINE *lwline)
void circ_tree_print(const CIRC_NODE *node, int depth)
static CIRC_NODE * circ_node_leaf_new(const POINTARRAY *pa, int i)
Create a new leaf node, storing pointers back to the end points for later.
static int circ_center_spherical(const GEOGRAPHIC_POINT *c1, const GEOGRAPHIC_POINT *c2, double distance, double offset, GEOGRAPHIC_POINT *center)
Given the centers of two circles, and the offset distance we want to put the new center between them ...
static int circ_node_is_leaf(const CIRC_NODE *node)
Internal nodes have their point references set to NULL.
void circ_tree_free(CIRC_NODE *node)
Recurse from top of node tree and free all children.
static void circ_internal_nodes_sort(CIRC_NODE **nodes, uint32_t num_nodes, const CIRC_NODE *target_node)
int circ_tree_contains_point(const CIRC_NODE *node, const POINT2D *pt, const POINT2D *pt_outside, int *on_boundary)
Walk the tree and count intersections between the stab line and the edges.
static int circ_center_cartesian(const GEOGRAPHIC_POINT *c1, const GEOGRAPHIC_POINT *c2, double distance, double offset, GEOGRAPHIC_POINT *center)
Where the circ_center_spherical() function fails, we need a fall-back.
CIRC_NODE * lwgeom_calculate_circ_tree(const LWGEOM *lwgeom)
int circ_tree_get_point_outside(const CIRC_NODE *node, POINT2D *pt)
static CIRC_NODE * circ_nodes_merge(CIRC_NODE **nodes, int num_nodes)
static double circ_node_min_distance(const CIRC_NODE *n1, const CIRC_NODE *n2)
static CIRC_NODE * circ_node_leaf_point_new(const POINTARRAY *pa)
Return a point node (zero radius, referencing one point)
static int circ_node_compare(const void *v1, const void *v2)
Comparing on geohash ensures that nearby nodes will be close to each other in the list.
static void circ_nodes_sort(CIRC_NODE **nodes, int num_nodes)
Given a list of nodes, sort them into a spatially consistent order, then pairwise merge them up into ...
static CIRC_NODE * circ_node_internal_new(CIRC_NODE **c, uint32_t num_nodes)
Create a new internal node, calculating the new measure range for the node, and storing pointers to t...
#define CIRC_NODE_SIZE
Datum distance(PG_FUNCTION_ARGS)
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
GEOGRAPHIC_POINT start
Definition: lwgeodetic.h:63
GEOGRAPHIC_POINT end
Definition: lwgeodetic.h:64
Two-point great circle segment from a to b.
Definition: lwgeodetic.h:62
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:53
uint32_t ngeoms
Definition: liblwgeom.h:510
LWGEOM ** geoms
Definition: liblwgeom.h:512
uint8_t type
Definition: liblwgeom.h:399
POINTARRAY * points
Definition: liblwgeom.h:425
POINTARRAY * point
Definition: liblwgeom.h:414
uint8_t type
Definition: liblwgeom.h:410
POINTARRAY ** rings
Definition: liblwgeom.h:460
uint32_t nrings
Definition: liblwgeom.h:458
double y
Definition: liblwgeom.h:331
double x
Definition: liblwgeom.h:331
double z
Definition: liblwgeom.h:343
double x
Definition: liblwgeom.h:343
double y
Definition: liblwgeom.h:343
uint32_t npoints
Definition: liblwgeom.h:374
double radius
Definition: liblwgeom.h:321
double a
Definition: liblwgeom.h:316
double b
Definition: liblwgeom.h:317
uint32_t num_nodes
POINT2D * p2
struct circ_node ** nodes
double radius
POINT2D * p1
POINT2D pt_outside
GEOGRAPHIC_POINT center
uint32_t geom_type
Note that p1 and p2 are pointers into an independent POINTARRAY, do not free them.
CIRC_NODE * node
unsigned int uint32_t
Definition: uthash.h:78