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