80 LWDEBUGF(3,
"edge #%d (%g %g, %g %g)", i, p1->
x, p1->
y, p2->
x, p2->
y);
100 node->
radius = diameter / 2.0;
153 if ( u1 < u2 )
return -1;
154 if ( u1 > u2 )
return 1;
170 LWDEBUGF(4,
"calculating spherical center", dir);
194 double proportion = offset/
distance;
196 LWDEBUG(4,
"calculating cartesian center");
202 p1p2.
x = p2.
x - p1.
x;
203 p1p2.
y = p2.
y - p1.
y;
204 p1p2.
z = p2.
z - p1.
z;
207 p1p2.
x *= proportion;
208 p1p2.
y *= proportion;
209 p1p2.
z *= proportion;
212 pc.
x = p1.
x + p1p2.
x;
213 pc.
y = p1.
y + p1p2.
y;
214 pc.
z = p1.
z + p1p2.
z;
234 double offset1, dist, D, r1, ri;
235 int i, new_geom_type;
237 LWDEBUGF(3,
"called with %d nodes --", num_nodes);
244 new_center = c[0]->
center;
245 new_radius = c[0]->
radius;
249 for ( i = 1; i < num_nodes; i++ )
259 if ( ! new_geom_type )
267 if ( new_geom_type != c[i]->geom_type )
287 LWDEBUG(3,
" distance between centers is zero");
288 new_radius = r1 + 2*dist;
291 else if ( dist < fabs(r1 - ri) )
304 new_center = c[i]->
center;
310 LWDEBUG(3,
" calculating new center");
316 new_radius = D / 2.0;
319 offset1 = ri + (D - (2.0*r1 + 2.0*ri)) / 2.0;
320 LWDEBUGF(3,
" offset1 is %g", offset1);
332 LWDEBUGF(3,
" new center is (%g %g) new radius is %g", new_center.
lon, new_center.
lat, new_radius);
338 node->
center = new_center;
339 node->
radius = new_radius;
373 for ( i = 0; i < num_edges; i++ )
411 int num_children = num_nodes;
420 while( num_children > 1 )
422 for ( j = 0; j < num_children; j++ )
425 if ( inode_num == 0 )
428 inodes[inode_num] = nodes[j];
430 if ( inode_num == CIRC_NODE_SIZE-1 )
435 if ( inode_num == 0 )
438 nodes[num_parents++] = inodes[0];
447 num_children = num_parents;
503 LWDEBUGF(3,
"edge_distance_to_point=%g, node_radius=%g", d, node->
radius);
506 LWDEBUGF(3,
"entering this branch (%p)", node);
522 LWDEBUG(3,
" got stab line edge_intersection with this edge!");
528 LWDEBUG(3,
" rejecting stab line grazing by left-side edge");
533 LWDEBUG(3,
" accepting stab line intersection");
544 LWDEBUG(3,
"internal node calculation");
545 LWDEBUGF(3,
" calling circ_tree_contains_point on child %d!", i);
553 LWDEBUGF(3,
"skipping this branch (%p)", node);
581 double min_dist = FLT_MAX;
582 double max_dist = FLT_MAX;
587 double threshold_radians = 0.95 * threshold / spheroid->
radius;
592 if ( spheroid->
a == spheroid->
b )
616 if( *min_dist < threshold || *min_dist == 0.0 )
622 LWDEBUGF(4,
"pruning pair %p, %p", n1, n2);
629 if( max < *max_dist )
638 LWDEBUGF(4,
"n1 is polygon, testing if contains (%.5g,%.5g)", pt.
x, pt.
y);
654 LWDEBUGF(4,
"n2 is polygon, testing if contains (%.5g,%.5g)", pt.
x, pt.
y);
672 if ( n1->
p1 == n1->
p2 || n2->
p1 == n2->
p2 )
678 if ( n1->
p1 == n1->
p2 && n2->
p1 == n2->
p2 )
682 close1 = gp1; close2 = gp2;
686 else if ( n1->
p1 == n1->
p2 )
729 LWDEBUGF(4,
"edge_distance_to_edge returned %g", d);
796 printf(
"%*s[%d] C(%.5g %.5g) R(%.5g) ((%.5g %.5g),(%.5g,%.5g))",
797 3*depth + 6,
"NODE", node->
edge_num,
816 printf(
"%*s C(%.5g %.5g) R(%.5g)",
865 if ( lwpoly->
nrings == 1 )
873 for ( i = 0; i < lwpoly->
nrings; i++ )
908 for ( i = 0; i < lwcol->
ngeoms; i++ )
930 switch ( lwgeom->
type )
double sphere_distance(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Given two points on a unit sphere, calculate their distance apart in radians.
static CIRC_NODE * lwpoint_calculate_circ_tree(const LWPOINT *lwpoint)
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 ...
int lwtype_get_collectiontype(uint8_t type)
Given an lwtype number, what homogeneous collection can hold it?
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.
Two-point great circle segment from a to b.
#define PIR_B_TOUCH_RIGHT
void normalize(POINT3D *p)
Normalize to a unit vector.
int circ_tree_get_point(const CIRC_NODE *node, POINT2D *pt)
Returns a POINT2D that is a vertex of the input shape.
uint32_t lwgeom_get_type(const LWGEOM *geom)
Return LWTYPE number.
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...
void cart2geog(const POINT3D *p, GEOGRAPHIC_POINT *g)
Convert cartesion coordinates on unit sphere to spherical coordinates.
#define LWDEBUG(level, msg)
static int circ_center_spherical(const GEOGRAPHIC_POINT *c1, const GEOGRAPHIC_POINT *c2, double distance, double offset, GEOGRAPHIC_POINT *center)
Given the centers of two circles, and the offset distance we want to put the new center between them ...
static 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)
Point in spherical coordinates on the world.
CIRC_NODE * lwgeom_calculate_circ_tree(const LWGEOM *lwgeom)
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.
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)
double edge_distance_to_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *gp, GEOGRAPHIC_POINT *closest)
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)
void circ_tree_print(const CIRC_NODE *node, int depth)
int lwtype_is_collection(uint8_t type)
Determine whether a type number is a collection or not.
void lwpoly_pt_outside(const LWPOLY *poly, POINT2D *pt_outside)
static CIRC_NODE * circ_nodes_merge(CIRC_NODE **nodes, int num_nodes)
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)
void geog2cart(const GEOGRAPHIC_POINT *g, POINT3D *p)
Convert spherical coordinates to cartesion coordinates on unit sphere.
void vector_sum(const POINT3D *a, const POINT3D *b, POINT3D *n)
Calculate the sum of two vectors.
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.
static double circ_node_min_distance(const CIRC_NODE *n1, const CIRC_NODE *n2)
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
static CIRC_NODE * circ_node_leaf_point_new(const POINTARRAY *pa)
Return a point node (zero radius, referencing one point)
static CIRC_NODE * lwpoly_calculate_circ_tree(const LWPOLY *lwpoly)
static CIRC_NODE * lwcollection_calculate_circ_tree(const LWCOLLECTION *lwcol)
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.
void * lwalloc(size_t size)
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members) ...
#define LWDEBUGF(level, msg,...)
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.
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
struct circ_node ** nodes
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...
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.