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