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