if ( lon < -1.0 * M_PI )
lon = 2.0 * M_PI + lon;
-
+
if ( lon == -2.0 * M_PI )
lon *= -1.0;
double zmin = FLT_MAX;
double zmax = -1 * FLT_MAX;
POINT3D pt;
-
+
/* Take a copy of the box corners so we can treat them as a list */
/* Elements are xmin, xmax, ymin, ymax, zmin, zmax */
memcpy(d, &(gbox->xmin), 6*sizeof(double));
-
+
/* Generate all 8 corner vectors of the box */
for ( i = 0; i < 8; i++ )
{
{
double angle, dotprod;
POINT3D pt_n;
-
+
pt_n.x = d[i / 2];
pt_n.y = d[2 + (i % 2)];
magnitude = sqrt(pt_n.x*pt_n.x + pt_n.y*pt_n.y);
}
}
}
-
+
/* Return the distance between the two furthest vectors */
return maxangle;
}
/* Take a copy of the box corners so we can treat them as a list */
/* Elements are xmin, xmax, ymin, ymax, zmin, zmax */
memcpy(d, &(gbox->xmin), 6*sizeof(double));
-
+
/* Zero out our return vector */
pt.x = pt.y = pt.z = 0.0;
for ( i = 0; i < 8; i++ )
{
POINT3D pt_n;
-
+
pt_n.x = d[i / 4];
pt_n.y = d[2 + ((i % 4) / 2)];
pt_n.z = d[4 + (i % 2)];
normalize(&pt_n);
-
+
pt.x += pt_n.x;
pt.y += pt_n.y;
- pt.z += pt_n.z;
+ pt.z += pt_n.z;
}
-
+
pt.x /= 8.0;
pt.y /= 8.0;
pt.z /= 8.0;
cart2geog(&pt, &g);
out->x = longitude_degrees_normalize(rad2deg(g.lon));
out->y = latitude_degrees_normalize(rad2deg(g.lat));
-
+
return LW_SUCCESS;
}
cross_product(v1, v2, &normal);
normalize(&normal);
cross_product(&normal, v1, &v3);
-
+
x = dot_product(v1, v2);
y = dot_product(v2, &v3);
-
+
angle = atan2(y, x);
return angle;
}
{
double p_dot = dot_product(P1, P2);
POINT3D P3;
-
+
/* If edge is really large, calculate a narrower equivalent angle A1/A3. */
if ( p_dot < 0 )
{
{
P3 = *P2;
}
-
+
/* Normals to the A-plane and B-plane */
cross_product(P1, &P3, normal);
normalize(normal);
double uxuy, uyuz, uxuz;
double ux2, uy2, uz2;
double rxx, rxy, rxz, ryx, ryy, ryz, rzx, rzy, rzz;
-
+
/* Need a unit vector normal to rotate around */
unit_normal(v1, v2, &u);
-
+
uxuy = u.x * u.y;
uxuz = u.x * u.z;
uyuz = u.y * u.z;
-
+
ux2 = u.x * u.x;
uy2 = u.y * u.y;
uz2 = u.z * u.z;
-
+
rxx = cos_a + ux2 * (1 - cos_a);
rxy = uxuy * (1 - cos_a) - u.z * sin_a;
rxz = uxuz * (1 - cos_a) + u.y * sin_a;
-
+
ryx = uxuy * (1 - cos_a) + u.z * sin_a;
ryy = cos_a + uy2 * (1 - cos_a);
ryz = uyuz * (1 - cos_a) - u.x * sin_a;
-
+
rzx = uxuz * (1 - cos_a) - u.y * sin_a;
rzy = uyuz * (1 - cos_a) + u.x * sin_a;
rzz = cos_a + uz2 * (1 - cos_a);
LWDEBUG(4, "point is on plane (dot product is zero)");
return 0;
}
-
+
if ( w < 0 )
return -1;
else
robust_cross_product(b, c, &normal2);
normalize(&normal1);
normalize(&normal2);
- return sphere_distance_cartesian(&normal1, &normal2);
+ return sphere_distance_cartesian(&normal1, &normal2);
}
/**
double area_radians = 0.0;
int side;
GEOGRAPHIC_EDGE e;
-
+
angle_a = sphere_angle(b,a,c);
angle_b = sphere_angle(a,b,c);
angle_c = sphere_angle(b,c,a);
-
+
area_radians = angle_a + angle_b + angle_c - M_PI;
/* What's the direction of the B/C edge? */
e.start = *a;
e.end = *b;
side = edge_point_side(&e, c);
-
+
/* Co-linear points implies no area */
if ( side == 0 )
return 0.0;
int side = edge_point_side(e, p);
if ( side == 0 )
return LW_TRUE;
-
+
return LW_FALSE;
}
{
double heading = 0.0;
double f;
-
+
/* Starting from the poles? Special case. */
if ( FP_IS_ZERO(cos(s->lat)) )
return (s->lat > 0.0) ? M_PI : 0.0;
{
lon2 = lon1 + atan2(sin(azimuth)*sin(d)*cos(lat1), cos(d)-sin(lat1)*sin(lat2));
}
-
+
if ( isnan(lat2) || isnan(lon2) )
return LW_FAILURE;
/* Initialize the box with the edge end points */
gbox_init_point3d(A1, gbox);
gbox_merge_point3d(A2, gbox);
-
+
/* Zero length edge, just return! */
if ( p3d_same(A1, A2) )
return LW_SUCCESS;
-
+
/* Error out on antipodal edge */
if ( FP_EQUALS(A1->x, -1*A2->x) && FP_EQUALS(A1->y, -1*A2->y) && FP_EQUALS(A1->z, -1*A2->z) )
{
lwerror("Antipodal (180 degrees long) edge detected!");
return LW_FAILURE;
}
-
+
/* Create A3, a vector in the plane of A1/A2, orthogonal to A1 */
unit_normal(A1, A2, &AN);
unit_normal(&AN, A1, &A3);
memset(X, 0, sizeof(POINT3D) * 6);
X[0].x = X[2].y = X[4].z = 1.0;
X[1].x = X[3].y = X[5].z = -1.0;
-
+
/* Initialize a 2-space origin point. */
O.x = O.y = 0.0;
/* What side of the line joining R1/R2 is O? */
o_side = lw_segment_side(&R1, &R2, &O);
-
+
/* Add any extrema! */
for ( i = 0; i < 6; i++ )
{
RX.x = dot_product(&(X[i]), A1);
RX.y = dot_product(&(X[i]), &A3);
normalize2d(&RX);
-
+
/* Any axis end on the side of R1/R2 opposite the origin */
/* is an extreme point in the arc, so we add the 3-space */
/* version of the point on R1/R2 to the gbox */
Xn.x = RX.x * A1->x + RX.y * A3.x;
Xn.y = RX.x * A1->y + RX.y * A3.y;
Xn.z = RX.x * A1->z + RX.y * A3.z;
-
+
gbox_merge_point3d(&Xn, gbox);
}
}
}
void lwpoly_pt_outside(const LWPOLY *poly, POINT2D *pt_outside)
-{
+{
/* Make sure we have boxes */
if ( poly->bbox )
{
}
+static int ptarray_segmentize_sphere_edge_recursive (
+ const POINT3D *p1, const POINT3D *p2, /* 3-space points we are interpolating beween */
+ const POINT4D *v1, const POINT4D *v2, /* real values and z/m values */
+ double d, double max_seg_length, /* current segment length and segment limit */
+ POINTARRAY *pa) /* write out results here */
+{
+ /* Reached the terminal leaf in recursion. Add */
+ /* the left-most point to the pointarray here */
+ /* We recurse down the left side first, so outputs should */
+ /* end up added to the array in order this way */
+ if (d <= max_seg_length)
+ {
+ GEOGRAPHIC_POINT g;
+ POINT4D p;
+ cart2geog(p1, &g);
+ p.x = rad2deg(g.lon);
+ p.y = rad2deg(g.lat);
+ p.z = v1->z;
+ p.m = v1->m;
+ return ptarray_append_point(pa, &p, LW_FALSE);
+ }
+ /* Find the mid-point and recurse on the left and then the right */
+ else
+ {
+ /* Calculate mid-point */
+ POINT3D mid;
+ mid.x = (p1->x + p2->x) / 2.0;
+ mid.y = (p1->y + p2->y) / 2.0;
+ mid.z = (p1->z + p2->z) / 2.0;
+ normalize(&mid);
+
+ /* Calculate z/m mid-values */
+ /* (ignore x/y, we get those from the 3-space calculations) */
+ POINT4D midv;
+ midv.z = (v1->z + v2->z) / 2.0;
+ midv.m = (v1->m + v2->m) / 2.0;
+ /* Recurse on the left first */
+ ptarray_segmentize_sphere_edge_recursive(p1, &mid, v1, &midv, d/2.0, max_seg_length, pa);
+ ptarray_segmentize_sphere_edge_recursive(&mid, p2, &midv, v2, d/2.0, max_seg_length, pa);
+ return LW_SUCCESS;
+ }
+}
+
/**
* Create a new point array with no segment longer than the input segment length (expressed in radians!)
* @param pa_in - input point array pointer
POINTARRAY *pa_out;
int hasz = ptarray_has_z(pa_in);
int hasm = ptarray_has_m(pa_in);
- int pa_in_offset = 0; /* input point offset */
- POINT4D p1, p2, p;
- GEOGRAPHIC_POINT g1, g2, g;
- double d;
-
+ POINT4D p1, p2;
+ POINT3D q1, q2;
+ GEOGRAPHIC_POINT g1, g2;
+ int i;
+
/* Just crap out on crazy input */
if ( ! pa_in )
- lwerror("ptarray_segmentize_sphere: null input pointarray");
- if ( max_seg_length <= 0.0 )
- lwerror("ptarray_segmentize_sphere: maximum segment length must be positive");
+ lwerror("%s: null input pointarray", __func__);
+ if ( max_seg_length <= 0.0 )
+ lwerror("%s: maximum segment length must be positive", __func__);
/* Empty starting array */
pa_out = ptarray_construct_empty(hasz, hasm, pa_in->npoints);
- /* Add first point */
- getPoint4d_p(pa_in, pa_in_offset, &p1);
- ptarray_append_point(pa_out, &p1, LW_FALSE);
- geographic_point_init(p1.x, p1.y, &g1);
- pa_in_offset++;
-
- while ( pa_in_offset < pa_in->npoints )
+ /* Simple loop per edge */
+ for (i = 1; i < pa_in->npoints; i++)
{
- getPoint4d_p(pa_in, pa_in_offset, &p2);
+ getPoint4d_p(pa_in, i-1, &p1);
+ getPoint4d_p(pa_in, i, &p2);
+ geographic_point_init(p1.x, p1.y, &g1);
geographic_point_init(p2.x, p2.y, &g2);
-
+
/* Skip duplicate points (except in case of 2-point lines!) */
- if ( (pa_in->npoints > 2) && p4d_same(&p1, &p2) )
- {
- /* Move one offset forward */
- p1 = p2;
- g1 = g2;
- pa_in_offset++;
+ if ((pa_in->npoints > 2) && p4d_same(&p1, &p2))
continue;
- }
/* How long is this edge? */
- d = sphere_distance(&g1, &g2);
-
- /* We need to segmentize this edge */
- if ( d > max_seg_length )
+ double d = sphere_distance(&g1, &g2);
+
+ if (d > max_seg_length)
{
- int nsegs = 1 + d / max_seg_length;
- int i;
- double dzz = 0, dmm = 0;
- double delta = d / nsegs;
-
- /* The independent Z/M values on the ptarray */
- if ( hasz ) dzz = (p2.z - p1.z) / nsegs;
- if ( hasm ) dmm = (p2.m - p1.m) / nsegs;
-
- g = g1;
- p = p1;
- for ( i = 0; i < nsegs - 1; i++ )
- {
- GEOGRAPHIC_POINT gn;
- double heading;
-
- /* Compute the current heading to the destination */
- heading = sphere_direction(&g, &g2, (nsegs-i) * delta);
- /* Move one increment forwards */
- sphere_project(&g, delta, heading, &gn);
- g = gn;
-
- p.x = rad2deg(g.lon);
- p.y = rad2deg(g.lat);
- if ( hasz )
- p.z += dzz;
- if ( hasm )
- p.m += dmm;
- ptarray_append_point(pa_out, &p, LW_FALSE);
- }
-
- ptarray_append_point(pa_out, &p2, LW_FALSE);
+ geog2cart(&g1, &q1);
+ geog2cart(&g2, &q2);
+ /* 3-d end points, XYZM end point, current edge size, min edge size */
+ ptarray_segmentize_sphere_edge_recursive(&q1, &q2, &p1, &p2, d, max_seg_length, pa_out);
}
- /* This edge is already short enough */
+ /* If we don't segmentize, we need to add first point manually */
else
{
- ptarray_append_point(pa_out, &p2, (pa_in->npoints==2)?LW_TRUE:LW_FALSE);
+ ptarray_append_point(pa_out, &p1, LW_TRUE);
}
-
- /* Move one offset forward */
- p1 = p2;
- g1 = g2;
- pa_in_offset++;
}
-
- return pa_out;
+ /* Always add the last point */
+ ptarray_append_point(pa_out, &p2, LW_TRUE);
+ return pa_out;
}
/**
LWPOLY *lwpoly_in, *lwpoly_out;
LWCOLLECTION *lwcol_in, *lwcol_out;
int i;
-
+
/* Reflect NULL */
if ( ! lwg_in )
return NULL;
-
+
/* Clone empty */
if ( lwgeom_is_empty(lwg_in) )
return lwgeom_clone(lwg_in);
-
+
switch (lwg_in->type)
{
case MULTIPOINTTYPE:
lwg_in->type, lwtype_name(lwg_in->type));
break;
}
-
+
lwerror("lwgeom_segmentize_sphere got to the end of the function, should not happen");
return NULL;
}
const POINT2D *p;
GEOGRAPHIC_POINT a, b, c;
double area = 0.0;
-
+
/* Return zero on nonsensical inputs */
if ( ! pa || pa->npoints < 4 )
return 0.0;
-
+
p = getPoint2d_cp(pa, 0);
geographic_point_init(p->x, p->y, &a);
p = getPoint2d_cp(pa, 1);
geographic_point_init(p->x, p->y, &b);
-
+
for ( i = 2; i < pa->npoints-1; i++ )
{
p = getPoint2d_cp(pa, i);
area += sphere_signed_area(&a, &b, &c);
b = c;
}
-
+
return fabs(area);
}
lwerror("Distance must be between 0 and %g", M_PI * spheroid->radius);
return NULL;
}
-
+
/* Convert to ta geodetic point */
x = lwpoint_get_x(r);
y = lwpoint_get_y(r);
geographic_point_init(x, y, &geo_source);
-
+
/* Try the projection */
if( spheroid_project(&geo_source, spheroid, distance, azimuth, &geo_dest) == LW_FAILURE )
{
lwerror("Unable to project from (%g %g) with azimuth %g and distance %g", x, y, azimuth, distance);
return NULL;
}
-
+
/* Build the output LWPOINT */
pa = ptarray_construct(0, 0, 1);
pt_dest.x = rad2deg(longitude_radians_normalize(geo_dest.lon));
x2 = lwpoint_get_x(s);
y2 = lwpoint_get_y(s);
geographic_point_init(x2, y2, &g2);
-
+
/* Same point, return NaN */
if ( FP_EQUALS(x1, x2) && FP_EQUALS(y1, y2) )
{
return NAN;
}
-
+
/* Do the direction calculation */
return spheroid_direction(&g1, &g2, spheroid);
}
gbox_init(&gbox1);
gbox_init(&gbox2);
-
+
assert(lwgeom1);
assert(lwgeom2);
-
+
LWDEBUGF(4, "entered function, tolerance %.8g", tolerance);
/* What's the distance to an empty geometry? We don't know.
{
return 0.0;
}
-
+
/* Not inside, so what's the actual distance? */
for ( i = 0; i < lwpoly->nrings; i++ )
{
int type1, type2;
GBOX gbox1, gbox2;
gbox1.flags = gbox2.flags = 0;
-
+
assert(lwgeom1);
assert(lwgeom2);
p = getPoint2d_cp(pa, 0);
ll2cart(p, &A1);
-
+
for ( i = 1; i < pa->npoints; i++ )
{
-
+
p = getPoint2d_cp(pa, i);
ll2cart(p, &A2);
-
+
edge_calculate_gbox(&A1, &A2, &edge_gbox);
/* Initialize the box */
{
gbox_merge(&edge_gbox, gbox);
}
-
+
A1 = A2;
}
{
if ( lwgeom_is_empty(geom) )
return LW_TRUE;
-
+
switch (geom->type)
{
case POINTTYPE:
/* Add in the vertical displacement if we're in 3D */
if ( hasz )
seglength = sqrt( (zb-za)*(zb-za) + seglength*seglength );
-
+
/* Add this segment length to the total */
length += seglength;
{
POINT3D AC; /* Center point of A1/A2 */
double min_similarity, similarity;
-
+
/* The normalized sum bisects the angle between start and end. */
vector_sum(A1, A2, &AC);
normalize(&AC);
-
+
/* The projection of start onto the center defines the minimum similarity */
min_similarity = dot_product(A1, &AC);
if ( FP_IS_ZERO(dp) )
return 0;
-
+
return dp < 0.0 ? -1 : 1;
}
double ab_dot;
int a1_side, a2_side, b1_side, b2_side;
int rv = PIR_NO_INTERACT;
-
+
/* Normals to the A-plane and B-plane */
unit_normal(A1, A2, &AN);
unit_normal(B1, B2, &BN);
-
+
/* Are A-plane and B-plane basically the same? */
ab_dot = dot_product(&AN, &BN);
if ( FP_EQUALS(fabs(ab_dot), 1.0) )
}
return rv;
}
-
+
/* What side of plane-A and plane-B do the end points */
/* of A and B fall? */
a1_side = dot_product_side(&BN, A1);
{
return PIR_INTERSECTS;
}
-
+
return PIR_NO_INTERACT;
}
/* Touches at B2, B1 is on what side? */
rv |= (b1_side < 0 ? PIR_B_TOUCH_RIGHT : PIR_B_TOUCH_LEFT);
}
-
+
return rv;
}
{
continue;
}
-
+
/* Our test point is on an edge end! Point is "in ring" by our definition */
if ( point3d_equals(&S1, &E1) )
{
return LW_TRUE;
}
-
+
/* Calculate relationship between stab line and edge */
inter = edge_intersects(&S1, &S2, &E1, &E2);
-
+
/* We have some kind of interaction... */
if ( inter & PIR_INTERSECTS )
{
{
return LW_TRUE;
}
-
+
/* It's a touching interaction, disregard all the left-side ones. */
/* It's a co-linear intersection, ignore those. */
if ( inter & PIR_B_TOUCH_RIGHT || inter & PIR_COLINEAR )
{
LWDEBUGF(4," edge (%d) did not cross", i);
}
-
+
/* Increment to next edge */
E1 = E2;
}
projects / postgis / commitdiff