PostGIS  2.5.1dev-r@@SVN_REVISION@@
lwtree.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) 2009-2012 Paul Ramsey <pramsey@cleverelephant.ca>
22  *
23  **********************************************************************/
24 
25 #include "liblwgeom_internal.h"
26 #include "lwgeom_log.h"
27 #include "lwtree.h"
28 #include "measures.h"
29 
30 static inline int
32 {
33  return node->type == RECT_NODE_LEAF_TYPE;
34 }
35 
36 /*
37 * Support qsort of nodes for collection/multi types so nodes
38 * are in "spatial adjacent" order prior to merging.
39 */
40 static int
41 rect_node_cmp(const void *pn1, const void *pn2)
42 {
43  GBOX b1, b2;
44  RECT_NODE *n1 = *((RECT_NODE**)pn1);
45  RECT_NODE *n2 = *((RECT_NODE**)pn2);
46  uint64_t h1, h2;
47  b1.flags = 0;
48  b1.xmin = n1->xmin;
49  b1.xmax = n1->xmax;
50  b1.ymin = n1->ymin;
51  b1.ymax = n1->ymax;
52 
53  b2.flags = 0;
54  b2.xmin = n2->xmin;
55  b2.xmax = n2->xmax;
56  b2.ymin = n2->ymin;
57  b2.ymax = n2->ymax;
58 
59  h1 = gbox_get_sortable_hash(&b1);
60  h2 = gbox_get_sortable_hash(&b2);
61  return h1 < h2 ? -1 : (h1 > h2 ? 1 : 0);
62 }
63 
68 void
70 {
71  int i;
72  if (!node) return;
73  if (!rect_node_is_leaf(node))
74  {
75  for (i = 0; i < node->i.num_nodes; i++)
76  {
77  rect_tree_free(node->i.nodes[i]);
78  node->i.nodes[i] = NULL;
79  }
80  }
81  lwfree(node);
82 }
83 
84 static int
86 {
87  const POINT2D *p1, *p2, *p3, *q1, *q2, *q3;
88  DISTPTS dl;
89  lw_dist2d_distpts_init(&dl, 1);
90  switch (n1->seg_type)
91  {
93  {
94  p1 = getPoint2d_cp(n1->pa, n1->seg_num);
95 
96  switch (n2->seg_type)
97  {
99  q1 = getPoint2d_cp(n2->pa, n2->seg_num);
100  lw_dist2d_pt_pt(q1, p1, &dl);
101  return dl.distance == 0.0;
102 
104  q1 = getPoint2d_cp(n2->pa, n2->seg_num);
105  q2 = getPoint2d_cp(n2->pa, n2->seg_num+1);
106  lw_dist2d_pt_seg(p1, q1, q2, &dl);
107  return dl.distance == 0.0;
108 
110  q1 = getPoint2d_cp(n2->pa, n2->seg_num*2);
111  q2 = getPoint2d_cp(n2->pa, n2->seg_num*2+1);
112  q3 = getPoint2d_cp(n2->pa, n2->seg_num*2+2);
113  lw_dist2d_pt_arc(p1, q1, q2, q3, &dl);
114  return dl.distance == 0.0;
115 
116  default:
117  lwerror("%s: unsupported segment type", __func__);
118  break;
119  }
120 
121  break;
122  }
123 
125  {
126  p1 = getPoint2d_cp(n1->pa, n1->seg_num);
127  p2 = getPoint2d_cp(n1->pa, n1->seg_num+1);
128 
129  switch (n2->seg_type)
130  {
131  case RECT_NODE_SEG_POINT:
132  q1 = getPoint2d_cp(n2->pa, n2->seg_num);
133  lw_dist2d_pt_seg(q1, p1, p2, &dl);
134  return dl.distance == 0.0;
135 
137  q1 = getPoint2d_cp(n2->pa, n2->seg_num);
138  q2 = getPoint2d_cp(n2->pa, n2->seg_num+1);
139  return lw_segment_intersects(p1, p2, q1, q2) > 0;
140 
142  q1 = getPoint2d_cp(n2->pa, n2->seg_num*2);
143  q2 = getPoint2d_cp(n2->pa, n2->seg_num*2+1);
144  q3 = getPoint2d_cp(n2->pa, n2->seg_num*2+2);
145  lw_dist2d_seg_arc(p1, p2, q1, q2, q3, &dl);
146  return dl.distance == 0.0;
147 
148  default:
149  lwerror("%s: unsupported segment type", __func__);
150  break;
151  }
152 
153  break;
154  }
156  {
157  p1 = getPoint2d_cp(n1->pa, n1->seg_num*2);
158  p2 = getPoint2d_cp(n1->pa, n1->seg_num*2+1);
159  p3 = getPoint2d_cp(n1->pa, n1->seg_num*2+2);
160 
161  switch (n2->seg_type)
162  {
163  case RECT_NODE_SEG_POINT:
164  q1 = getPoint2d_cp(n2->pa, n2->seg_num);
165  lw_dist2d_pt_arc(q1, p1, p2, p3, &dl);
166  return dl.distance == 0.0;
167 
169  q1 = getPoint2d_cp(n2->pa, n2->seg_num);
170  q2 = getPoint2d_cp(n2->pa, n2->seg_num+1);
171  lw_dist2d_seg_arc(q1, q2, p1, p2, p3, &dl);
172  return dl.distance == 0.0;
173 
175  q1 = getPoint2d_cp(n2->pa, n2->seg_num*2);
176  q2 = getPoint2d_cp(n2->pa, n2->seg_num*2+1);
177  q3 = getPoint2d_cp(n2->pa, n2->seg_num*2+2);
178  lw_dist2d_arc_arc(p1, p2, p3, q1, q2, q3, &dl);
179  return dl.distance == 0.0;
180 
181  default:
182  lwerror("%s: unsupported segment type", __func__);
183  break;
184  }
185 
186  break;
187  }
188  default:
189  return LW_FALSE;
190  }
191  return LW_FALSE;
192 }
193 
194 
195 /*
196 * Returns 1 if segment is to the right of point.
197 */
198 static inline int
199 rect_leaf_node_segment_side(RECT_NODE_LEAF *node, const POINT2D *q, int *on_boundary)
200 {
201  const POINT2D *p1, *p2, *p3;
202  switch (node->seg_type)
203  {
205  {
206  int side;
207  p1 = getPoint2d_cp(node->pa, node->seg_num);
208  p2 = getPoint2d_cp(node->pa, node->seg_num+1);
209 
210  side = lw_segment_side(p1, p2, q);
211 
212  /* Always note case where we're on boundary */
213  if (side == 0 && lw_pt_in_seg(q, p1, p2))
214  {
215  *on_boundary = LW_TRUE;
216  return 0;
217  }
218 
219  /* Segment points up and point is on left */
220  if (p1->y < p2->y && side == -1 && q->y != p2->y)
221  {
222  return 1;
223  }
224 
225  /* Segment points down and point is on right */
226  if (p1->y > p2->y && side == 1 && q->y != p2->y)
227  {
228  return 1;
229  }
230 
231  /* Segment is horizontal, do we cross first point? */
232  if (p1->y == p2->y && q->x < p1->x)
233  {
234  return 1;
235  }
236 
237  return 0;
238  }
240  {
241  int arc_side, seg_side;
242 
243  p1 = getPoint2d_cp(node->pa, node->seg_num*2);
244  p2 = getPoint2d_cp(node->pa, node->seg_num*2+1);
245  p3 = getPoint2d_cp(node->pa, node->seg_num*2+2);
246 
247  /* Always note case where we're on boundary */
248  arc_side = lw_arc_side(p1, p2, p3, q);
249  if (arc_side == 0)
250  {
251  *on_boundary = LW_TRUE;
252  return 0;
253  }
254 
255  seg_side = lw_segment_side(p1, p3, q);
256  if (seg_side == arc_side)
257  {
258  /* Segment points up and point is on left */
259  if (p1->y < p3->y && seg_side == -1 && q->y != p3->y)
260  {
261  return 1;
262  }
263 
264  /* Segment points down and point is on right */
265  if (p1->y > p3->y && seg_side == 1 && q->y != p3->y)
266  {
267  return 1;
268  }
269  }
270  else
271  {
272  /* Segment points up and point is on left */
273  if (p1->y < p3->y && seg_side == 1 && q->y != p3->y)
274  {
275  return 1;
276  }
277 
278  /* Segment points down and point is on right */
279  if (p1->y > p3->y && seg_side == -1 && q->y != p3->y)
280  {
281  return 1;
282  }
283 
284  /* Segment is horizontal, do we cross first point? */
285  if (p1->y == p3->y)
286  {
287  return 1;
288  }
289  }
290 
291  return 0;
292 
293  }
294  default:
295  {
296  lwerror("%s: unsupported seg_type - %d", __func__, node->seg_type);
297  return 0;
298  }
299  }
300 
301  return 0;
302 }
303 
304 /*
305 * Only pass the head of a ring. Either a LinearRing from a polygon
306 * or a CompoundCurve from a CurvePolygon.
307 * Takes a horizontal line through the ring, and adds up the
308 * crossing directions. An "inside" point will be on the same side
309 * ("both left" or "both right") of the edges the line crosses,
310 * so it will have a non-zero return sum. An "outside" point will
311 * be on both sides, and will have a zero return sum.
312 */
313 static int
314 rect_tree_ring_contains_point(RECT_NODE *node, const POINT2D *pt, int *on_boundary)
315 {
316  /* Only test nodes that straddle our stabline vertically */
317  /* and might be to the right horizontally */
318  if (node->ymin <= pt->y && pt->y <= node->ymax && pt->x <= node->xmax)
319  {
320  if (rect_node_is_leaf(node))
321  {
322  return rect_leaf_node_segment_side(&node->l, pt, on_boundary);
323  }
324  else
325  {
326  int i, r = 0;
327  for (i = 0; i < node->i.num_nodes; i++)
328  {
329  r += rect_tree_ring_contains_point(node->i.nodes[i], pt, on_boundary);
330  }
331  return r;
332  }
333  }
334  return 0;
335 }
336 
337 /*
338 * Only pass in the head of an "area" type. Polygon or CurvePolygon.
339 * Sums up containment of exterior (+1) and interior (-1) rings, so
340 * that zero is uncontained, +1 is contained and negative is an error
341 * (multiply contained by interior rings?)
342 */
343 static int
345 {
346  /* Can't do anything with a leaf node, makes no sense */
347  if (rect_node_is_leaf(node))
348  return 0;
349 
350  /* Iterate into area until we find ring heads */
351  if (node->i.ring_type == RECT_NODE_RING_NONE)
352  {
353  int i, sum = 0;
354  for (i = 0; i < node->i.num_nodes; i++)
355  sum += rect_tree_area_contains_point(node->i.nodes[i], pt);
356  return sum;
357  }
358  /* See if the ring encloses the point */
359  else
360  {
361  int on_boundary = 0;
362  int edge_crossing_count = rect_tree_ring_contains_point(node, pt, &on_boundary);
363  /* Odd number of crossings => contained */
364  int contained = (edge_crossing_count % 2 == 1);
365  /* External rings return positive containment, interior ones negative, */
366  /* so that a point-in-hole case nets out to zero (contained by both */
367  /* interior and exterior rings. */
368  if (node->i.ring_type == RECT_NODE_RING_INTERIOR)
369  {
370  return on_boundary ? 0 : -1 * contained;
371  }
372  else
373  {
374  return contained || on_boundary;
375  }
376 
377  }
378 }
379 
380 /*
381 * Simple containment test for node/point inputs
382 */
383 static int
385 {
386  if (pt->y < node->ymin || pt->y > node->ymax ||
387  pt->x < node->xmin || pt->x > node->xmax)
388  return LW_FALSE;
389  else
390  return LW_TRUE;
391 }
392 
393 /*
394 * Pass in arbitrary tree, get back true if point is contained or on boundary,
395 * and false otherwise.
396 */
397 int
399 {
400  int i, c;
401 
402  /* Object cannot contain point if bounds don't */
403  if (!rect_node_bounds_point(node, pt))
404  return 0;
405 
406  switch (node->geom_type)
407  {
408  case POLYGONTYPE:
409  case CURVEPOLYTYPE:
410  return rect_tree_area_contains_point(node, pt) > 0;
411 
412  case MULTIPOLYGONTYPE:
413  case MULTISURFACETYPE:
414  case COLLECTIONTYPE:
415  {
416  for (i = 0; i < node->i.num_nodes; i++)
417  {
418  c = rect_tree_contains_point(node->i.nodes[i], pt);
419  if (c) return LW_TRUE;
420  }
421  return LW_FALSE;
422  }
423 
424  default:
425  return LW_FALSE;
426  }
427 }
428 
429 /*
430 * For area types, doing intersects and distance, we will
431 * need to do a point-in-poly test first to find the full-contained
432 * case where an intersection exists without any edges actually
433 * intersecting.
434 */
435 static int
437 {
438  switch (node->geom_type)
439  {
440  case POLYGONTYPE:
441  case CURVEPOLYTYPE:
442  case MULTISURFACETYPE:
443  return LW_TRUE;
444 
445  case COLLECTIONTYPE:
446  {
447  if (rect_node_is_leaf(node))
448  return LW_FALSE;
449  else
450  {
451  int i;
452  for (i = 0; i < node->i.num_nodes; i++)
453  {
454  if (rect_tree_is_area(node->i.nodes[i]))
455  return LW_TRUE;
456  }
457  return LW_FALSE;
458  }
459  }
460 
461  default:
462  return LW_FALSE;
463  }
464 }
465 
467 {
468  RECT_NODE_SEG_UNKNOWN, // "Unknown"
469  RECT_NODE_SEG_POINT, // "Point"
470  RECT_NODE_SEG_LINEAR, // "LineString"
471  RECT_NODE_SEG_LINEAR, // "Polygon"
472  RECT_NODE_SEG_UNKNOWN, // "MultiPoint"
473  RECT_NODE_SEG_LINEAR, // "MultiLineString"
474  RECT_NODE_SEG_LINEAR, // "MultiPolygon"
475  RECT_NODE_SEG_UNKNOWN, // "GeometryCollection"
476  RECT_NODE_SEG_CIRCULAR, // "CircularString"
477  RECT_NODE_SEG_UNKNOWN, // "CompoundCurve"
478  RECT_NODE_SEG_UNKNOWN, // "CurvePolygon"
479  RECT_NODE_SEG_UNKNOWN, // "MultiCurve"
480  RECT_NODE_SEG_UNKNOWN, // "MultiSurface"
481  RECT_NODE_SEG_LINEAR, // "PolyhedralSurface"
482  RECT_NODE_SEG_LINEAR, // "Triangle"
483  RECT_NODE_SEG_LINEAR // "Tin"
484 };
485 
486 /*
487 * Create a new leaf node.
488 */
489 static RECT_NODE *
490 rect_node_leaf_new(const POINTARRAY *pa, int seg_num, int geom_type)
491 {
492  const POINT2D *p1, *p2, *p3;
493  RECT_NODE *node;
494  GBOX gbox;
495  RECT_NODE_SEG_TYPE seg_type = lwgeomTypeArc[geom_type];
496 
497  switch (seg_type)
498  {
499  case RECT_NODE_SEG_POINT:
500  {
501  p1 = getPoint2d_cp(pa, seg_num);
502  gbox.xmin = gbox.xmax = p1->x;
503  gbox.ymin = gbox.ymax = p1->y;
504  break;
505  }
506 
508  {
509  p1 = getPoint2d_cp(pa, seg_num);
510  p2 = getPoint2d_cp(pa, seg_num+1);
511  /* Zero length edge, doesn't get a node */
512  if ((p1->x == p2->x) && (p1->y == p2->y))
513  return NULL;
514  gbox.xmin = FP_MIN(p1->x, p2->x);
515  gbox.xmax = FP_MAX(p1->x, p2->x);
516  gbox.ymin = FP_MIN(p1->y, p2->y);
517  gbox.ymax = FP_MAX(p1->y, p2->y);
518  break;
519  }
520 
522  {
523  p1 = getPoint2d_cp(pa, 2*seg_num);
524  p2 = getPoint2d_cp(pa, 2*seg_num+1);
525  p3 = getPoint2d_cp(pa, 2*seg_num+2);
526  /* Zero length edge, doesn't get a node */
527  if ((p1->x == p2->x) && (p2->x == p3->x) &&
528  (p1->y == p2->y) && (p2->y == p3->y))
529  return NULL;
530  lw_arc_calculate_gbox_cartesian_2d(p1, p2, p3, &gbox);
531  break;
532  }
533 
534  default:
535  {
536  lwerror("%s: unsupported seg_type - %d", __func__, seg_type);
537  return NULL;
538  }
539  }
540 
541  node = lwalloc(sizeof(RECT_NODE));
542  node->type = RECT_NODE_LEAF_TYPE;
543  node->geom_type = geom_type;
544  node->xmin = gbox.xmin;
545  node->xmax = gbox.xmax;
546  node->ymin = gbox.ymin;
547  node->ymax = gbox.ymax;
548  node->l.seg_num = seg_num;
549  node->l.seg_type = seg_type;
550  node->l.pa = pa;
551  return node;
552 }
553 
554 
555 static void
557 {
558  if (rect_node_is_leaf(node))
559  lwerror("%s: call on leaf node", __func__);
560  node->xmin = FP_MIN(node->xmin, add->xmin);
561  node->xmax = FP_MAX(node->xmax, add->xmax);
562  node->ymin = FP_MIN(node->ymin, add->ymin);
563  node->ymax = FP_MAX(node->ymax, add->ymax);
564  node->i.nodes[node->i.num_nodes++] = add;
565  return;
566 }
567 
568 
569 static RECT_NODE *
571 {
572  RECT_NODE *node = lwalloc(sizeof(RECT_NODE));
573  node->xmin = seed->xmin;
574  node->xmax = seed->xmax;
575  node->ymin = seed->ymin;
576  node->ymax = seed->ymax;
577  node->geom_type = seed->geom_type;
579  node->i.num_nodes = 0;
581  node->i.sorted = 0;
582  return node;
583 }
584 
585 /*
586 * We expect the incoming nodes to be in a spatially coherent
587 * order. For incoming nodes derived from point arrays,
588 * the very fact that they are
589 * a vertex list implies a reasonable ordering: points nearby in
590 * the list will be nearby in space. For incoming nodes from higher
591 * level structures (collections, etc) the caller should sort the
592 * nodes using a z-order first, so that this merge step results in a
593 * spatially coherent structure.
594 */
595 static RECT_NODE *
596 rect_nodes_merge(RECT_NODE ** nodes, uint32_t num_nodes)
597 {
598  if (num_nodes < 1)
599  {
600  return NULL;
601  }
602 
603  while (num_nodes > 1)
604  {
605  uint32_t i, k = 0;
606  RECT_NODE *node = NULL;
607  for (i = 0; i < num_nodes; i++)
608  {
609  if (!node)
610  node = rect_node_internal_new(nodes[i]);
611 
612  rect_node_internal_add_node(node, nodes[i]);
613 
614  if (node->i.num_nodes == RECT_NODE_SIZE)
615  {
616  nodes[k++] = node;
617  node = NULL;
618  }
619  }
620  if (node)
621  nodes[k++] = node;
622  num_nodes = k;
623  }
624 
625  return nodes[0];
626 }
627 
628 /*
629 * Build a tree of nodes from a point array, one node per edge.
630 */
631 RECT_NODE *
632 rect_tree_from_ptarray(const POINTARRAY *pa, int geom_type)
633 {
634  int num_edges = 0, i = 0, j = 0;
635  RECT_NODE_SEG_TYPE seg_type = lwgeomTypeArc[geom_type];
636  RECT_NODE **nodes = NULL;
637  RECT_NODE *tree = NULL;
638 
639  /* No-op on empty ring/line/pt */
640  if ( pa->npoints < 1 )
641  return NULL;
642 
643  /* For arcs, 3 points per edge, for lines, 2 per edge */
644  switch(seg_type)
645  {
646  case RECT_NODE_SEG_POINT:
647  return rect_node_leaf_new(pa, 0, geom_type);
648  break;
650  num_edges = pa->npoints - 1;
651  break;
653  num_edges = (pa->npoints - 1)/2;
654  break;
655  default:
656  lwerror("%s: unsupported seg_type - %d", __func__, seg_type);
657  }
658 
659  /* First create a flat list of nodes, one per edge. */
660  nodes = lwalloc(sizeof(RECT_NODE*) * num_edges);
661  for (i = 0; i < num_edges; i++)
662  {
663  RECT_NODE *node = rect_node_leaf_new(pa, i, geom_type);
664  if (node) /* Not zero length? */
665  nodes[j++] = node;
666  }
667 
668  /* Merge the list into a tree */
669  tree = rect_nodes_merge(nodes, j);
670 
671  /* Free the old list structure, leaving the tree in place */
672  lwfree(nodes);
673 
674  /* Return top of tree */
675  return tree;
676 }
677 
678 LWGEOM *
680 {
681  LWGEOM *poly = (LWGEOM*)lwpoly_construct_envelope(0, node->xmin, node->ymin, node->xmax, node->ymax);
682  if (rect_node_is_leaf(node))
683  {
684  return poly;
685  }
686  else
687  {
688  int i;
690  lwcollection_add_lwgeom(col, poly);
691  for (i = 0; i < node->i.num_nodes; i++)
692  {
694  }
695  return (LWGEOM*)col;
696  }
697 }
698 
699 char *
701 {
702  LWGEOM *geom = rect_tree_to_lwgeom(node);
703  char *wkt = lwgeom_to_wkt(geom, WKT_ISO, 12, 0);
704  lwgeom_free(geom);
705  return wkt;
706 }
707 
708 void
709 rect_tree_printf(const RECT_NODE *node, int depth)
710 {
711  printf("%*s----\n", depth, "");
712  printf("%*stype: %d\n", depth, "", node->type);
713  printf("%*sgeom_type: %d\n", depth, "", node->geom_type);
714  printf("%*sbox: %g %g, %g %g\n", depth, "", node->xmin, node->ymin, node->xmax, node->ymax);
715  if (node->type == RECT_NODE_LEAF_TYPE)
716  {
717  printf("%*sseg_type: %d\n", depth, "", node->l.seg_type);
718  printf("%*sseg_num: %d\n", depth, "", node->l.seg_num);
719  }
720  else
721  {
722  int i;
723  for (i = 0; i < node->i.num_nodes; i++)
724  {
725  rect_tree_printf(node->i.nodes[i], depth+2);
726  }
727  }
728 }
729 
730 static RECT_NODE *
732 {
733  const LWPOINT *lwpt = (const LWPOINT*)lwgeom;
734  return rect_tree_from_ptarray(lwpt->point, lwgeom->type);
735 }
736 
737 static RECT_NODE *
739 {
740  const LWLINE *lwline = (const LWLINE*)lwgeom;
741  return rect_tree_from_ptarray(lwline->points, lwgeom->type);
742 }
743 
744 static RECT_NODE *
746 {
747  RECT_NODE **nodes;
748  RECT_NODE *tree;
749  uint32_t i, j = 0;
750  const LWPOLY *lwpoly = (const LWPOLY*)lwgeom;
751 
752  if (lwpoly->nrings < 1)
753  return NULL;
754 
755  nodes = lwalloc(sizeof(RECT_NODE*) * lwpoly->nrings);
756  for (i = 0; i < lwpoly->nrings; i++)
757  {
758  RECT_NODE *node = rect_tree_from_ptarray(lwpoly->rings[i], lwgeom->type);
759  if (node)
760  {
762  nodes[j++] = node;
763  }
764  }
765  tree = rect_nodes_merge(nodes, j);
766  tree->geom_type = lwgeom->type;
767  lwfree(nodes);
768  return tree;
769 }
770 
771 static RECT_NODE *
773 {
774  RECT_NODE **nodes;
775  RECT_NODE *tree;
776  uint32_t i, j = 0;
777  const LWCURVEPOLY *lwcol = (const LWCURVEPOLY*)lwgeom;
778 
779  if (lwcol->nrings < 1)
780  return NULL;
781 
782  nodes = lwalloc(sizeof(RECT_NODE*) * lwcol->nrings);
783  for (i = 0; i < lwcol->nrings; i++)
784  {
785  RECT_NODE *node = rect_tree_from_lwgeom(lwcol->rings[i]);
786  if (node)
787  {
788  /*
789  * In the case of arc circle, it's possible for a ring to consist
790  * of a single closed edge. That will arrive as a leaf node. We
791  * need to wrap that node in an internal node with an appropriate
792  * ring type so all the other code can try and make sense of it.
793  */
794  if (node->type == RECT_NODE_LEAF_TYPE)
795  {
796  RECT_NODE *internal = rect_node_internal_new(node);
797  rect_node_internal_add_node(internal, node);
798  node = internal;
799  }
800  /* Each subcomponent is a ring */
802  nodes[j++] = node;
803  }
804  }
805  /* Put the top nodes in a z-order curve for a spatially coherent */
806  /* tree after node merge */
807  qsort(nodes, j, sizeof(RECT_NODE*), rect_node_cmp);
808 
809  tree = rect_nodes_merge(nodes, j);
810 
811  tree->geom_type = lwgeom->type;
812  lwfree(nodes);
813  return tree;
814 
815 }
816 
817 static RECT_NODE *
819 {
820  RECT_NODE **nodes;
821  RECT_NODE *tree;
822  uint32_t i, j = 0;
823  const LWCOLLECTION *lwcol = (const LWCOLLECTION*)lwgeom;
824 
825  if (lwcol->ngeoms < 1)
826  return NULL;
827 
828  /* Build one tree for each sub-geometry, then below */
829  /* we merge the root notes of those trees to get a single */
830  /* top node for the collection */
831  nodes = lwalloc(sizeof(RECT_NODE*) * lwcol->ngeoms);
832  for (i = 0; i < lwcol->ngeoms; i++)
833  {
834  RECT_NODE *node = rect_tree_from_lwgeom(lwcol->geoms[i]);
835  if (node)
836  {
837  /* Curvepolygons are collections where the sub-geometries */
838  /* are the rings, and will need to doint point-in-poly */
839  /* tests in order to do intersects and distance calculations */
840  /* correctly */
841  if (lwgeom->type == CURVEPOLYTYPE)
843  nodes[j++] = node;
844  }
845  }
846  /* Sort the nodes using a z-order curve, so that merging the nodes */
847  /* gives a spatially coherent tree (near things are in near nodes) */
848  /* Note: CompoundCurve has edges already spatially organized, no */
849  /* sorting needed */
850  if (lwgeom->type != COMPOUNDTYPE)
851  qsort(nodes, j, sizeof(RECT_NODE*), rect_node_cmp);
852 
853  tree = rect_nodes_merge(nodes, j);
854 
855  tree->geom_type = lwgeom->type;
856  lwfree(nodes);
857  return tree;
858 }
859 
860 RECT_NODE *
862 {
863  switch(lwgeom->type)
864  {
865  case POINTTYPE:
866  return rect_tree_from_lwpoint(lwgeom);
867  case TRIANGLETYPE:
868  case CIRCSTRINGTYPE:
869  case LINETYPE:
870  return rect_tree_from_lwline(lwgeom);
871  case POLYGONTYPE:
872  return rect_tree_from_lwpoly(lwgeom);
873  case CURVEPOLYTYPE:
874  return rect_tree_from_lwcurvepoly(lwgeom);
875  case COMPOUNDTYPE:
876  case MULTICURVETYPE:
877  case MULTISURFACETYPE:
878  case MULTIPOINTTYPE:
879  case MULTILINETYPE:
880  case MULTIPOLYGONTYPE:
882  case TINTYPE:
883  case COLLECTIONTYPE:
884  return rect_tree_from_lwcollection(lwgeom);
885  default:
886  lwerror("%s: Unknown geometry type: %s", __func__, lwtype_name(lwgeom->type));
887  return NULL;
888  }
889  return NULL;
890 }
891 
892 /*
893 * Get an actual coordinate point from a tree to use
894 * for point-in-polygon testing.
895 */
896 static const POINT2D *
898 {
899  if (!node) return NULL;
900  if (rect_node_is_leaf(node))
901  return getPoint2d_cp(node->l.pa, 0);
902  else
903  return rect_tree_get_point(node->i.nodes[0]);
904 }
905 
906 static inline int
908 {
909  if (n1->xmin > n2->xmax || n2->xmin > n1->xmax ||
910  n1->ymin > n2->ymax || n2->ymin > n1->ymax)
911  {
912  return 0;
913  }
914  else
915  {
916  return 1;
917  }
918 }
919 
920 #if POSTGIS_DEBUG_LEVEL >= 4
921 static char *
922 rect_node_to_str(const RECT_NODE *n)
923 {
924  char *buf = lwalloc(256);
925  snprintf(buf, 256, "(%.9g %.9g,%.9g %.9g)",
926  n->xmin, n->ymin, n->xmax, n->ymax);
927  return buf;
928 }
929 #endif
930 
931 /*
932 * Work down to leaf nodes, until we find a pair of leaf nodes
933 * that intersect. Prune branches that do not intersect.
934 */
935 static int
937 {
938  int i, j;
939 #if POSTGIS_DEBUG_LEVEL >= 4
940  char *n1_str = rect_node_to_str(n1);
941  char *n2_str = rect_node_to_str(n2);
942  LWDEBUGF(4,"n1 %s n2 %s", n1, n2);
943  lwfree(n1_str);
944  lwfree(n2_str);
945 #endif
946  /* There can only be an edge intersection if the rectangles overlap */
947  if (rect_node_intersects(n1, n2))
948  {
949  LWDEBUG(4," interaction found");
950  /* We can only test for a true intersection if the nodes are both leaf nodes */
951  if (rect_node_is_leaf(n1) && rect_node_is_leaf(n2))
952  {
953  LWDEBUG(4," leaf node test");
954  /* Check for true intersection */
955  return rect_leaf_node_intersects(&n1->l, &n2->l);
956  }
957  else if (rect_node_is_leaf(n2) && !rect_node_is_leaf(n1))
958  {
959  for (i = 0; i < n1->i.num_nodes; i++)
960  {
962  return LW_TRUE;
963  }
964  }
965  else if (rect_node_is_leaf(n1) && !rect_node_is_leaf(n2))
966  {
967  for (i = 0; i < n2->i.num_nodes; i++)
968  {
970  return LW_TRUE;
971  }
972  }
973  else
974  {
975  for (j = 0; j < n1->i.num_nodes; j++)
976  {
977  for (i = 0; i < n2->i.num_nodes; i++)
978  {
979  if (rect_tree_intersects_tree_recursive(n2->i.nodes[i], n1->i.nodes[j]))
980  return LW_TRUE;
981  }
982  }
983  }
984  }
985  return LW_FALSE;
986 }
987 
988 
989 int
991 {
992  /*
993  * It is possible for an area to intersect another object
994  * without any edges intersecting, if the object is fully contained.
995  * If that is so, then any point in the object will be contained,
996  * so we do a quick point-in-poly test first for those cases
997  */
998  if (rect_tree_is_area(n1) &&
1000  {
1001  return LW_TRUE;
1002  }
1003 
1004  if (rect_tree_is_area(n2) &&
1006  {
1007  return LW_TRUE;
1008  }
1009 
1010  /*
1011  * Not contained, so intersection can only happen if
1012  * edges actually intersect.
1013  */
1014  return rect_tree_intersects_tree_recursive(n1, n2);
1015 }
1016 
1017 /*
1018 * For slightly more speed when doing distance comparisons,
1019 * where we only care if d1 < d2 and not what the actual value
1020 * of d1 or d2 is, we use the squared distance and avoid the
1021 * sqrt()
1022 */
1023 static inline double
1024 distance_sq(double x1, double y1, double x2, double y2)
1025 {
1026  double dx = x1-x2;
1027  double dy = y1-y2;
1028  return dx*dx + dy*dy;
1029 }
1030 
1031 static inline double
1032 distance(double x1, double y1, double x2, double y2)
1033 {
1034  return sqrt(distance_sq(x1, y1, x2, y2));
1035 }
1036 
1037 /*
1038 * The closest any two objects in two nodes can be is the smallest
1039 * distance between the nodes themselves.
1040 */
1041 static inline double
1043 {
1044  int left = n1->xmin > n2->xmax;
1045  int right = n1->xmax < n2->xmin;
1046  int bottom = n1->ymin > n2->ymax;
1047  int top = n1->ymax < n2->ymin;
1048 
1049  //lwnotice("rect_node_min_distance");
1050 
1051  if (top && left)
1052  return distance(n1->xmin, n1->ymax, n2->xmax, n2->ymin);
1053  else if (top && right)
1054  return distance(n1->xmax, n1->ymax, n2->xmin, n2->ymin);
1055  else if (bottom && left)
1056  return distance(n1->xmin, n1->ymin, n2->xmax, n2->ymax);
1057  else if (bottom && right)
1058  return distance(n1->xmax, n1->ymin, n2->xmin, n2->ymax);
1059  else if (left)
1060  return n1->xmin - n2->xmax;
1061  else if (right)
1062  return n2->xmin - n1->xmax;
1063  else if (bottom)
1064  return n1->ymin - n2->ymax;
1065  else if (top)
1066  return n2->ymin - n1->ymax;
1067  else
1068  return 0.0;
1069 
1070  return 0.0;
1071 }
1072 
1073 /*
1074 * The furthest apart the objects in two nodes can be is if they
1075 * are at opposite corners of the bbox that contains both nodes
1076 */
1077 static inline double
1079 {
1080  double xmin = FP_MIN(n1->xmin, n2->xmin);
1081  double ymin = FP_MIN(n1->ymin, n2->ymin);
1082  double xmax = FP_MAX(n1->xmax, n2->xmax);
1083  double ymax = FP_MAX(n1->ymax, n2->ymax);
1084  double dx = xmax - xmin;
1085  double dy = ymax - ymin;
1086  //lwnotice("rect_node_max_distance");
1087  return sqrt(dx*dx + dy*dy);
1088 }
1089 
1090 /*
1091 * Leaf nodes represent individual edges from the original shape.
1092 * As such, they can be either points (if original was a (multi)point)
1093 * two-point straight edges (for linear types), or
1094 * three-point circular arcs (for curvilinear types).
1095 * The distance calculations for each possible combination of primitive
1096 * edges is different, so there's a big case switch in here to match
1097 * up the right combination of inputs to the right distance calculation.
1098 */
1099 static double
1101 {
1102  const POINT2D *p1, *p2, *p3, *q1, *q2, *q3;
1103  DISTPTS dl;
1104 
1105  //lwnotice("rect_leaf_node_distance, %d<->%d", n1->seg_num, n2->seg_num);
1106 
1108 
1109  switch (n1->seg_type)
1110  {
1111  case RECT_NODE_SEG_POINT:
1112  {
1113  p1 = getPoint2d_cp(n1->pa, n1->seg_num);
1114 
1115  switch (n2->seg_type)
1116  {
1117  case RECT_NODE_SEG_POINT:
1118  q1 = getPoint2d_cp(n2->pa, n2->seg_num);
1119  lw_dist2d_pt_pt(q1, p1, &dl);
1120  break;
1121 
1122  case RECT_NODE_SEG_LINEAR:
1123  q1 = getPoint2d_cp(n2->pa, n2->seg_num);
1124  q2 = getPoint2d_cp(n2->pa, n2->seg_num+1);
1125  lw_dist2d_pt_seg(p1, q1, q2, &dl);
1126  break;
1127 
1129  q1 = getPoint2d_cp(n2->pa, n2->seg_num*2);
1130  q2 = getPoint2d_cp(n2->pa, n2->seg_num*2+1);
1131  q3 = getPoint2d_cp(n2->pa, n2->seg_num*2+2);
1132  lw_dist2d_pt_arc(p1, q1, q2, q3, &dl);
1133  break;
1134 
1135  default:
1136  lwerror("%s: unsupported segment type", __func__);
1137  }
1138  break;
1139  }
1140 
1141  case RECT_NODE_SEG_LINEAR:
1142  {
1143  p1 = getPoint2d_cp(n1->pa, n1->seg_num);
1144  p2 = getPoint2d_cp(n1->pa, n1->seg_num+1);
1145 
1146  switch (n2->seg_type)
1147  {
1148  case RECT_NODE_SEG_POINT:
1149  q1 = getPoint2d_cp(n2->pa, n2->seg_num);
1150  lw_dist2d_pt_seg(q1, p1, p2, &dl);
1151  break;
1152 
1153  case RECT_NODE_SEG_LINEAR:
1154  q1 = getPoint2d_cp(n2->pa, n2->seg_num);
1155  q2 = getPoint2d_cp(n2->pa, n2->seg_num+1);
1156  lw_dist2d_seg_seg(q1, q2, p1, p2, &dl);
1157  // lwnotice(
1158  // "%d\tLINESTRING(%g %g,%g %g)\t%d\tLINESTRING(%g %g,%g %g)\t%g\t%g\t%g",
1159  // n1->seg_num,
1160  // p1->x, p1->y, p2->x, p2->y,
1161  // n2->seg_num,
1162  // q1->x, q1->y, q2->x, q2->y,
1163  // dl.distance, state->min_dist, state->max_dist);
1164  break;
1165 
1167  q1 = getPoint2d_cp(n2->pa, n2->seg_num*2);
1168  q2 = getPoint2d_cp(n2->pa, n2->seg_num*2+1);
1169  q3 = getPoint2d_cp(n2->pa, n2->seg_num*2+2);
1170  lw_dist2d_seg_arc(p1, p2, q1, q2, q3, &dl);
1171  break;
1172 
1173  default:
1174  lwerror("%s: unsupported segment type", __func__);
1175  }
1176  break;
1177  }
1179  {
1180  p1 = getPoint2d_cp(n1->pa, n1->seg_num*2);
1181  p2 = getPoint2d_cp(n1->pa, n1->seg_num*2+1);
1182  p3 = getPoint2d_cp(n1->pa, n1->seg_num*2+2);
1183 
1184  switch (n2->seg_type)
1185  {
1186  case RECT_NODE_SEG_POINT:
1187  q1 = getPoint2d_cp(n2->pa, n2->seg_num);
1188  lw_dist2d_pt_arc(q1, p1, p2, p3, &dl);
1189  break;
1190 
1191  case RECT_NODE_SEG_LINEAR:
1192  q1 = getPoint2d_cp(n2->pa, n2->seg_num);
1193  q2 = getPoint2d_cp(n2->pa, n2->seg_num+1);
1194  lw_dist2d_seg_arc(q1, q2, p1, p2, p3, &dl);
1195  break;
1196 
1198  q1 = getPoint2d_cp(n2->pa, n2->seg_num*2);
1199  q2 = getPoint2d_cp(n2->pa, n2->seg_num*2+1);
1200  q3 = getPoint2d_cp(n2->pa, n2->seg_num*2+2);
1201  lw_dist2d_arc_arc(p1, p2, p3, q1, q2, q3, &dl);
1202  break;
1203 
1204  default:
1205  lwerror("%s: unsupported segment type", __func__);
1206  }
1207  break;
1208  }
1209  default:
1210  lwerror("%s: unsupported segment type", __func__);
1211  }
1212 
1213  /* If this is a new global minima, save it */
1214  if (dl.distance < state->min_dist)
1215  {
1216  state->min_dist = dl.distance;
1217  state->p1 = dl.p1;
1218  state->p2 = dl.p2;
1219  }
1220 
1221  return dl.distance;
1222 }
1223 
1224 
1225 static int
1226 rect_tree_node_sort_cmp(const void *a, const void *b)
1227 {
1228  RECT_NODE *n1 = *(RECT_NODE**)a;
1229  RECT_NODE *n2 = *(RECT_NODE**)b;
1230  if (n1->d < n2->d) return -1;
1231  else if (n1->d > n2->d) return 1;
1232  else return 0;
1233 }
1234 
1235 /* Calculate the center of a node */
1236 static inline void
1238 {
1239  pt->x = (n->xmin + n->xmax)/2;
1240  pt->y = (n->ymin + n->ymax)/2;
1241 }
1242 
1243 /*
1244 * (If necessary), sort the children of each node in
1245 * order of their distance from the enter of the other node.
1246 */
1247 static void
1249 {
1250  int i;
1251  RECT_NODE *n;
1252  POINT2D c1, c2, c;
1253  if (!rect_node_is_leaf(n1) && ! n1->i.sorted)
1254  {
1255  rect_node_to_p2d(n2, &c2);
1256  /* Distance of each child from center of other node */
1257  for (i = 0; i < n1->i.num_nodes; i++)
1258  {
1259  n = n1->i.nodes[i];
1260  rect_node_to_p2d(n, &c);
1261  n->d = distance2d_sqr_pt_pt(&c2, &c);
1262  }
1263  /* Sort the children by distance */
1264  n1->i.sorted = 1;
1265  qsort(n1->i.nodes,
1266  n1->i.num_nodes,
1267  sizeof(RECT_NODE*),
1269  }
1270  if (!rect_node_is_leaf(n2) && ! n2->i.sorted)
1271  {
1272  rect_node_to_p2d(n1, &c1);
1273  /* Distance of each child from center of other node */
1274  for (i = 0; i < n2->i.num_nodes; i++)
1275  {
1276  n = n2->i.nodes[i];
1277  rect_node_to_p2d(n, &c);
1278  n->d = distance2d_sqr_pt_pt(&c1, &c);
1279  }
1280  /* Sort the children by distance */
1281  n2->i.sorted = 1;
1282  qsort(n2->i.nodes,
1283  n2->i.num_nodes,
1284  sizeof(RECT_NODE*),
1286  }
1287 }
1288 
1289 static double
1291 {
1292  double min, max;
1293 
1294  /* Short circuit if we've already hit the minimum */
1295  if (state->min_dist < state->threshold || state->min_dist == 0.0)
1296  return state->min_dist;
1297 
1298  /* If your minimum is greater than anyone's maximum, you can't hold the winner */
1299  min = rect_node_min_distance(n1, n2);
1300  if (min > state->max_dist)
1301  {
1302  //lwnotice("pruning pair %p, %p", n1, n2);
1303  LWDEBUGF(4, "pruning pair %p, %p", n1, n2);
1304  return FLT_MAX;
1305  }
1306 
1307  /* If your maximum is a new low, we'll use that as our new global tolerance */
1308  max = rect_node_max_distance(n1, n2);
1309  if (max < state->max_dist)
1310  state->max_dist = max;
1311 
1312  /* Both leaf nodes, do a real distance calculation */
1313  if (rect_node_is_leaf(n1) && rect_node_is_leaf(n2))
1314  {
1315  return rect_leaf_node_distance(&n1->l, &n2->l, state);
1316  }
1317  /* Recurse into nodes */
1318  else
1319  {
1320  int i, j;
1321  double d_min = FLT_MAX;
1322  rect_tree_node_sort(n1, n2);
1323  if (rect_node_is_leaf(n1) && !rect_node_is_leaf(n2))
1324  {
1325  for (i = 0; i < n2->i.num_nodes; i++)
1326  {
1327  min = rect_tree_distance_tree_recursive(n1, n2->i.nodes[i], state);
1328  d_min = FP_MIN(d_min, min);
1329  }
1330  }
1331  else if (rect_node_is_leaf(n2) && !rect_node_is_leaf(n1))
1332  {
1333  for (i = 0; i < n1->i.num_nodes; i++)
1334  {
1335  min = rect_tree_distance_tree_recursive(n1->i.nodes[i], n2, state);
1336  d_min = FP_MIN(d_min, min);
1337  }
1338  }
1339  else
1340  {
1341  for (i = 0; i < n1->i.num_nodes; i++)
1342  {
1343  for (j = 0; j < n2->i.num_nodes; j++)
1344  {
1345  min = rect_tree_distance_tree_recursive(n1->i.nodes[i], n2->i.nodes[j], state);
1346  d_min = FP_MIN(d_min, min);
1347  }
1348  }
1349  }
1350  return d_min;
1351  }
1352 }
1353 
1354 double rect_tree_distance_tree(RECT_NODE *n1, RECT_NODE *n2, double threshold)
1355 {
1356  double distance;
1358 
1359  /*
1360  * It is possible for an area to intersect another object
1361  * without any edges intersecting, if the object is fully contained.
1362  * If that is so, then any point in the object will be contained,
1363  * so we do a quick point-in-poly test first for those cases
1364  */
1365  if (rect_tree_is_area(n1) &&
1367  {
1368  return 0.0;
1369  }
1370 
1371  if (rect_tree_is_area(n2) &&
1373  {
1374  return 0.0;
1375  }
1376 
1377  state.threshold = threshold;
1378  state.min_dist = FLT_MAX;
1379  state.max_dist = FLT_MAX;
1380  distance = rect_tree_distance_tree_recursive(n1, n2, &state);
1381  // *p1 = state.p1;
1382  // *p2 = state.p2;
1383  return distance;
1384 }
#define LINETYPE
Definition: liblwgeom.h:85
RECT_NODE_SEG_TYPE
Definition: lwtree.h:40
static void rect_tree_node_sort(RECT_NODE *n1, RECT_NODE *n2)
Definition: lwtree.c:1248
const POINTARRAY * pa
Definition: lwtree.h:50
#define MULTICURVETYPE
Definition: liblwgeom.h:94
static int rect_node_intersects(const RECT_NODE *n1, const RECT_NODE *n2)
Definition: lwtree.c:907
char * r
Definition: cu_in_wkt.c:24
char * lwgeom_to_wkt(const LWGEOM *geom, uint8_t variant, int precision, size_t *size_out)
WKT emitter function.
Definition: lwout_wkt.c:676
void lwfree(void *mem)
Definition: lwutil.c:244
int lw_arc_calculate_gbox_cartesian_2d(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, GBOX *gbox)
Definition: g_box.c:459
LWGEOM ** rings
Definition: liblwgeom.h:537
RECT_NODE_SEG_TYPE seg_type
Definition: lwtree.h:51
static double rect_tree_distance_tree_recursive(RECT_NODE *n1, RECT_NODE *n2, RECT_TREE_DISTANCE_STATE *state)
Definition: lwtree.c:1290
#define POLYGONTYPE
Definition: liblwgeom.h:86
double xmax
Definition: liblwgeom.h:295
#define CURVEPOLYTYPE
Definition: liblwgeom.h:93
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1144
#define COMPOUNDTYPE
Definition: liblwgeom.h:92
#define MULTIPOINTTYPE
Definition: liblwgeom.h:87
static double rect_node_max_distance(const RECT_NODE *n1, const RECT_NODE *n2)
Definition: lwtree.c:1078
static int rect_node_bounds_point(RECT_NODE *node, const POINT2D *pt)
Definition: lwtree.c:384
static void rect_node_internal_add_node(RECT_NODE *node, RECT_NODE *add)
Definition: lwtree.c:556
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
#define TRIANGLETYPE
Definition: liblwgeom.h:97
unsigned char geom_type
Definition: lwtree.h:68
#define POLYHEDRALSURFACETYPE
Definition: liblwgeom.h:96
static int rect_tree_is_area(const RECT_NODE *node)
Definition: lwtree.c:436
static RECT_NODE_SEG_TYPE lwgeomTypeArc[]
Definition: lwtree.c:466
static int rect_tree_node_sort_cmp(const void *a, const void *b)
Definition: lwtree.c:1226
static int rect_node_cmp(const void *pn1, const void *pn2)
Definition: lwtree.c:41
RECT_NODE_TYPE type
Definition: lwtree.h:67
static double rect_leaf_node_distance(const RECT_NODE_LEAF *n1, const RECT_NODE_LEAF *n2, RECT_TREE_DISTANCE_STATE *state)
Definition: lwtree.c:1100
#define FP_MIN(A, B)
uint32_t ngeoms
Definition: liblwgeom.h:509
POINTARRAY * point
Definition: liblwgeom.h:413
RECT_NODE_LEAF l
Definition: lwtree.h:76
int lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl)
Definition: measures.c:1497
static RECT_NODE * rect_tree_from_lwpoint(const LWGEOM *lwgeom)
Definition: lwtree.c:731
static int rect_tree_intersects_tree_recursive(RECT_NODE *n1, RECT_NODE *n2)
Definition: lwtree.c:936
uint32_t nrings
Definition: liblwgeom.h:457
static RECT_NODE * rect_tree_from_lwpoly(const LWGEOM *lwgeom)
Definition: lwtree.c:745
RECT_NODE_RING_TYPE ring_type
Definition: lwtree.h:60
int rect_tree_contains_point(RECT_NODE *node, const POINT2D *pt)
Definition: lwtree.c:398
int rect_tree_intersects_tree(RECT_NODE *n1, RECT_NODE *n2)
Test if two RECT_NODE trees intersect one another.
Definition: lwtree.c:990
double xmax
Definition: lwtree.h:70
POINT2D p1
Definition: measures.h:49
uint64_t gbox_get_sortable_hash(const GBOX *g)
Return a sortable key based on the center point of the GBOX.
Definition: g_serialized.c:258
static int rect_node_is_leaf(const RECT_NODE *node)
Definition: lwtree.c:31
unsigned int uint32_t
Definition: uthash.h:78
double x
Definition: liblwgeom.h:330
struct rect_node * nodes[RECT_NODE_SIZE]
Definition: lwtree.h:61
static RECT_NODE * rect_tree_from_lwcollection(const LWGEOM *lwgeom)
Definition: lwtree.c:818
LWPOLY * lwpoly_construct_envelope(int srid, double x1, double y1, double x2, double y2)
Definition: lwpoly.c:98
double ymin
Definition: liblwgeom.h:296
static int rect_tree_ring_contains_point(RECT_NODE *node, const POINT2D *pt, int *on_boundary)
Definition: lwtree.c:314
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:218
#define WKT_ISO
Definition: liblwgeom.h:2074
#define DIST_MIN
Definition: measures.h:41
double xmin
Definition: liblwgeom.h:294
#define LW_FALSE
Definition: liblwgeom.h:76
#define RECT_NODE_SIZE
Definition: lwtree.h:25
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
LWGEOM ** geoms
Definition: liblwgeom.h:511
int lw_dist2d_pt_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, DISTPTS *dl)
Definition: measures.c:1439
#define TINTYPE
Definition: liblwgeom.h:98
static double rect_node_min_distance(const RECT_NODE *n1, const RECT_NODE *n2)
Definition: lwtree.c:1042
LWGEOM * rect_tree_to_lwgeom(const RECT_NODE *node)
Definition: lwtree.c:679
static const POINT2D * rect_tree_get_point(const RECT_NODE *node)
Definition: lwtree.c:897
int lw_dist2d_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B, DISTPTS *dl)
lw_dist2d_comp from p to line A->B This one is now sending every occasion to lw_dist2d_pt_pt Before i...
Definition: measures.c:2205
double ymax
Definition: lwtree.h:72
double d
Definition: lwtree.h:73
POINTARRAY ** rings
Definition: liblwgeom.h:459
static int rect_tree_area_contains_point(RECT_NODE *node, const POINT2D *pt)
Definition: lwtree.c:344
int lw_dist2d_seg_seg(const POINT2D *A, const POINT2D *B, const POINT2D *C, const POINT2D *D, DISTPTS *dl)
Finds the shortest distance between two segments.
Definition: measures.c:1820
void rect_tree_printf(const RECT_NODE *node, int depth)
Definition: lwtree.c:709
POINT2D p2
Definition: measures.h:50
double ymax
Definition: liblwgeom.h:297
double y
Definition: liblwgeom.h:330
RECT_NODE_INTERNAL i
Definition: lwtree.h:75
uint8_t flags
Definition: liblwgeom.h:293
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:89
int seg_num
Definition: lwtree.h:52
char * rect_tree_to_wkt(const RECT_NODE *node)
Definition: lwtree.c:700
int lw_arc_side(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, const POINT2D *Q)
Definition: lwalgorithm.c:178
#define MULTISURFACETYPE
Definition: liblwgeom.h:95
static RECT_NODE * rect_tree_from_lwcurvepoly(const LWGEOM *lwgeom)
Definition: lwtree.c:772
double distance
Definition: measures.h:48
double xmin
Definition: lwtree.h:69
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
double ymin
Definition: lwtree.h:71
static double distance(double x1, double y1, double x2, double y2)
Definition: lwtree.c:1032
int lw_pt_in_seg(const POINT2D *P, const POINT2D *A1, const POINT2D *A2)
Returns true if P is between A1/A2.
Definition: lwalgorithm.c:95
static int rect_leaf_node_intersects(RECT_NODE_LEAF *n1, RECT_NODE_LEAF *n2)
Definition: lwtree.c:85
uint8_t type
Definition: liblwgeom.h:398
int lw_dist2d_seg_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl)
Calculate the shortest distance between an arc and an edge.
Definition: measures.c:1292
uint32_t nrings
Definition: liblwgeom.h:535
static RECT_NODE * rect_node_leaf_new(const POINTARRAY *pa, int seg_num, int geom_type)
Definition: lwtree.c:490
static double distance_sq(double x1, double y1, double x2, double y2)
Definition: lwtree.c:1024
double distance2d_sqr_pt_pt(const POINT2D *p1, const POINT2D *p2)
Definition: measures.c:2322
#define CIRCSTRINGTYPE
Definition: liblwgeom.h:91
static RECT_NODE * rect_tree_from_lwline(const LWGEOM *lwgeom)
Definition: lwtree.c:738
void rect_tree_free(RECT_NODE *node)
Recurse from top of node tree and free all children.
Definition: lwtree.c:69
double rect_tree_distance_tree(RECT_NODE *n1, RECT_NODE *n2, double threshold)
Return the distance between two RECT_NODE trees.
Definition: lwtree.c:1354
void * lwalloc(size_t size)
Definition: lwutil.c:229
int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
lw_segment_side()
Definition: lwalgorithm.c:64
Structure used in distance-calculations.
Definition: measures.h:46
void lw_dist2d_distpts_init(DISTPTS *dl, int mode)
Definition: measures.c:67
#define MULTILINETYPE
Definition: liblwgeom.h:88
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:94
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
RECT_NODE * rect_tree_from_lwgeom(const LWGEOM *lwgeom)
Create a tree index on top an LWGEOM.
Definition: lwtree.c:861
LWCOLLECTION * lwcollection_add_lwgeom(LWCOLLECTION *col, const LWGEOM *geom)
Appends geom to the collection managed by col.
Definition: lwcollection.c:187
static RECT_NODE * rect_node_internal_new(const RECT_NODE *seed)
Definition: lwtree.c:570
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
static RECT_NODE * rect_nodes_merge(RECT_NODE **nodes, uint32_t num_nodes)
Definition: lwtree.c:596
int lw_dist2d_pt_pt(const POINT2D *thep1, const POINT2D *thep2, DISTPTS *dl)
Compares incoming points and stores the points closest to each other or most far away from each other...
Definition: measures.c:2281
RECT_NODE * rect_tree_from_ptarray(const POINTARRAY *pa, int geom_type)
Definition: lwtree.c:632
static int rect_leaf_node_segment_side(RECT_NODE_LEAF *node, const POINT2D *q, int *on_boundary)
Definition: lwtree.c:199
static void rect_node_to_p2d(const RECT_NODE *n, POINT2D *pt)
Definition: lwtree.c:1237
#define COLLECTIONTYPE
Definition: liblwgeom.h:90
#define FP_MAX(A, B)
POINTARRAY * points
Definition: liblwgeom.h:424
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:364
int lw_segment_intersects(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2)
returns the kind of CG_SEGMENT_INTERSECTION_TYPE behavior of lineseg 1 (constructed from p1 and p2) a...
Definition: lwalgorithm.c:372
uint32_t npoints
Definition: liblwgeom.h:373