*/
static void test_lw_segment_side(void)
{
- double rv = 0.0;
+ int rv = 0;
POINT2D p1, p2, q;
/* Vertical line at x=0 */
q.y = 1.5;
rv = lw_segment_side(&p1, &p2, &q);
//printf("left %g\n",rv);
- CU_ASSERT(rv < 0.0);
+ CU_ASSERT(rv < 0);
/* On the right */
q.x = 2.0;
rv = lw_segment_side(&p1, &p2, &q);
//printf("right %g\n",rv);
- CU_ASSERT(rv > 0.0);
+ CU_ASSERT(rv > 0);
/* On the line */
q.x = 0.0;
rv = lw_segment_side(&p1, &p2, &q);
//printf("on line %g\n",rv);
- CU_ASSERT_EQUAL(rv, 0.0);
+ CU_ASSERT_EQUAL(rv, 0);
}
** Compute cartesian bounding GBOX boxes from LWGEOM.
*/
-static int lwcircle_calculate_gbox_cartesian(const POINT4D *p1, const POINT4D *p2, const POINT4D *p3, GBOX *gbox)
+int lwcircle_calculate_gbox_cartesian_2d(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, GBOX *gbox)
{
POINT2D xmin, ymin, xmax, ymax;
- POINT2D center;
- int p2_side;
- double radius;
+ POINT2D C;
+ int A2_side;
+ double radius_A;
- LWDEBUG(2, "lwcircle_calculate_gbox called.");
+ LWDEBUG(2, "lwcircle_calculate_gbox_cartesian_2d called.");
+
+ radius_A = lwcircle_center(A1, A2, A3, &C);
- radius = lwcircle_center((POINT2D*)p1, (POINT2D*)p2, (POINT2D*)p3, ¢er);
-
/* Negative radius signals straight line, p1/p2/p3 are colinear */
- if (radius < 0.0)
+ if (radius_A < 0.0)
{
- gbox->xmin = FP_MIN(p1->x, p3->x);
- gbox->ymin = FP_MIN(p1->y, p3->y);
- gbox->zmin = FP_MIN(p1->z, p3->z);
- gbox->xmax = FP_MAX(p1->x, p3->x);
- gbox->ymax = FP_MAX(p1->y, p3->y);
- gbox->zmax = FP_MAX(p1->z, p3->z);
+ gbox->xmin = FP_MIN(A1->x, A3->x);
+ gbox->ymin = FP_MIN(A1->y, A3->y);
+ gbox->xmax = FP_MAX(A1->x, A3->x);
+ gbox->ymax = FP_MAX(A1->y, A3->y);
return LW_SUCCESS;
}
-
+
/* Matched start/end points imply circle */
- if ( p1->x == p3->x && p1->y == p3->y )
+ if ( A1->x == A3->x && A1->y == A3->y )
{
- gbox->xmin = center.x - radius;
- gbox->ymin = center.y - radius;
- gbox->zmin = FP_MIN(p1->z,p2->z);
- gbox->mmin = FP_MIN(p1->m,p2->m);
- gbox->xmax = center.x + radius;
- gbox->ymax = center.y + radius;
- gbox->zmax = FP_MAX(p1->z,p2->z);
- gbox->mmax = FP_MAX(p1->m,p2->m);
+ gbox->xmin = C.x - radius_A;
+ gbox->ymin = C.y - radius_A;
+ gbox->xmax = C.x + radius_A;
+ gbox->ymax = C.y + radius_A;
return LW_SUCCESS;
}
/* First approximation, bounds of start/end points */
- gbox->xmin = FP_MIN(p1->x, p3->x);
- gbox->ymin = FP_MIN(p1->y, p3->y);
- gbox->zmin = FP_MIN(p1->z, p3->z);
- gbox->mmin = FP_MIN(p1->m, p3->m);
- gbox->xmax = FP_MAX(p1->x, p3->x);
- gbox->ymax = FP_MAX(p1->y, p3->y);
- gbox->zmax = FP_MAX(p1->z, p3->z);
- gbox->mmax = FP_MAX(p1->m, p3->m);
+ gbox->xmin = FP_MIN(A1->x, A3->x);
+ gbox->ymin = FP_MIN(A1->y, A3->y);
+ gbox->xmax = FP_MAX(A1->x, A3->x);
+ gbox->ymax = FP_MAX(A1->y, A3->y);
/* Create points for the possible extrema */
- xmin.x = center.x - radius;
- xmin.y = center.y;
- ymin.x = center.x;
- ymin.y = center.y - radius;
- xmax.x = center.x + radius;
- xmax.y = center.y;
- ymax.x = center.x;
- ymax.y = center.y + radius;
-
-
+ xmin.x = C.x - radius_A;
+ xmin.y = C.y;
+ ymin.x = C.x;
+ ymin.y = C.y - radius_A;
+ xmax.x = C.x + radius_A;
+ xmax.y = C.y;
+ ymax.x = C.x;
+ ymax.y = C.y + radius_A;
+
/* Divide the circle into two parts, one on each side of a line
joining p1 and p3. The circle extrema on the same side of that line
as p2 is on, are also the extrema of the bbox. */
-
- p2_side = signum(lw_segment_side((POINT2D*)p1, (POINT2D*)p3, (POINT2D*)p2));
- if ( p2_side == signum(lw_segment_side((POINT2D*)p1, (POINT2D*)p3, &xmin)) )
+ A2_side = lw_segment_side(A1, A3, A2);
+
+ if ( A2_side == lw_segment_side(A1, A3, &xmin) )
gbox->xmin = xmin.x;
- if ( p2_side == signum(lw_segment_side((POINT2D*)p1, (POINT2D*)p3, &ymin)) )
+ if ( A2_side == lw_segment_side(A1, A3, &ymin) )
gbox->ymin = ymin.y;
- if ( p2_side == signum(lw_segment_side((POINT2D*)p1, (POINT2D*)p3, &xmax)) )
+ if ( A2_side == lw_segment_side(A1, A3, &xmax) )
gbox->xmax = xmax.x;
- if ( p2_side == signum(lw_segment_side((POINT2D*)p1, (POINT2D*)p3, &ymax)) )
+ if ( A2_side == lw_segment_side(A1, A3, &ymax) )
gbox->ymax = ymax.y;
return LW_SUCCESS;
}
+
+static int lwcircle_calculate_gbox_cartesian(const POINT4D *p1, const POINT4D *p2, const POINT4D *p3, GBOX *gbox)
+{
+ int rv;
+
+ LWDEBUG(2, "lwcircle_calculate_gbox_cartesian called.");
+
+ rv = lwcircle_calculate_gbox_cartesian_2d((POINT2D*)p1, (POINT2D*)p2, (POINT2D*)p3, gbox);
+ gbox->zmin = FP_MIN(p1->z, p3->z);
+ gbox->mmin = FP_MIN(p1->m, p3->m);
+ gbox->zmax = FP_MAX(p1->z, p3->z);
+ gbox->mmax = FP_MAX(p1->m, p3->m);
+ return rv;
+}
+
int ptarray_calculate_gbox_cartesian(const POINTARRAY *pa, GBOX *gbox )
{
int i;
******************************************************************************/
int lwgeom_has_arc(const LWGEOM *geom);
-double lwcircle_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result);
LWGEOM *lwgeom_segmentize(LWGEOM *geom, uint32_t perQuad);
LWGEOM *lwgeom_desegmentize(LWGEOM *geom);
#define SIZE_GET(varsize) (((varsize) >> 2) & 0x3FFFFFFF)
#define SIZE_SET(varsize, size) (((varsize) & 0x00000003)|(((size) & 0x3FFFFFFF) << 2 ))
+/*
+* Tolerance used to determine equality.
+*/
+#define EPSILON_SQLMM 1e-8
/*
* Internal prototypes
*/
int lw_segment_intersects(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2);
-/*
-* What side of the line formed by p1 and p2 does q fall?
-* Returns < 0 for left and > 0 for right and 0 for co-linearity
-*/
-double lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q);
-
-/*
-* Do the envelopes of the the segments intersect?
-*/
-int lw_segment_envelope_intersects(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2);
-
/*
* Get/Set an enumeratoed ordinate. (x,y,z,m)
*/
+/*
+* What side of the line formed by p1 and p2 does q fall?
+* Returns -1 for left and 1 for right and 0 for co-linearity
+*/
+int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q);
+int lwcircle_calculate_gbox_cartesian_2d(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, GBOX *gbox);
+double lwcircle_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result);
int lw_pt_in_seg(const POINT2D *P, const POINT2D *A1, const POINT2D *A2);
+int lw_pt_in_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT2D *A3);
+int lw_arc_is_pt(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3);
+double lw_seg_length(const POINT2D *A1, const POINT2D *A2);
+double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3);
+double lw_arc_side(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, const POINT2D *Q);
+int pt_in_ring_2d(const POINT2D *p, const POINTARRAY *ring);
int ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt);
+int ptarray_arc_contains_point(const POINTARRAY *pa, const POINT2D *pt);
/**
* Split a line by a point and push components to the provided multiline.
/**
* lw_segment_side()
*
-* Return < 0.0 if point Q is left of segment P
-* Return > 0.0 if point Q is right of segment P
-* Return = 0.0 if point Q in on segment P
+* Return -1 if point Q is left of segment P
+* Return 1 if point Q is right of segment P
+* Return 0 if point Q in on segment P
*/
-double lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
+int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
{
double side = ( (q->x - p1->x) * (p2->y - p1->y) - (p2->x - p1->x) * (q->y - p1->y) );
if ( FP_IS_ZERO(side) )
+ return 0;
+ else
+ return signum(side);
+}
+
+/**
+* Returns the length of a linear segment
+*/
+double
+lw_seg_length(const POINT2D *A1, const POINT2D *A2)
+{
+ return sqrt((A1->x-A2->x)*(A1->x-A2->x)+(A1->y-A2->y)*(A1->y-A2->y));
+}
+
+/**
+* Returns true if P is on the same side of the plane partition
+* defined by A1/A3 as A2 is. Only makes sense if P has already been
+* determined to be on the circle defined by A1/A2/A3.
+*/
+int
+lw_pt_in_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
+{
+ return lw_segment_side(A1, A3, A2) == lw_segment_side(A1, A3, P);
+}
+
+/**
+* Returns true if P is between A1/A2. Only makes sense if P has already been
+* deterined to be on the line defined by A1/A2.
+*/
+int
+lw_pt_in_seg(const POINT2D *P, const POINT2D *A1, const POINT2D *A2)
+{
+ return ((A1->x <= P->x && P->x < A2->x) || (A1->x >= P->x && P->x > A2->x)) ||
+ ((A1->y <= P->y && P->y < A2->y) || (A1->y >= P->y && P->y > A2->y));
+}
+
+/**
+* Returns true if arc A is actually a point (all vertices are the same) .
+*/
+int
+lw_arc_is_pt(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
+{
+ if ( A1->x == A2->x && A2->x == A3->x &&
+ A1->y == A2->y && A2->y == A3->y )
+ return LW_TRUE;
+ else
+ return LW_FALSE;
+}
+
+/**
+* Returns the length of a circular arc segment
+*/
+double
+lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
+{
+ POINT2D C;
+ double radius_A, circumference_A;
+ int a2_side, clockwise;
+ double a1, a2, a3;
+ double angle;
+
+ if ( lw_arc_is_pt(A1, A2, A3) )
return 0.0;
+
+ radius_A = lwcircle_center(A1, A2, A3, &C);
+
+ /* Co-linear! Return linear distance! */
+ if ( radius_A < 0 )
+ {
+ return sqrt((A1->x-A3->x)*(A1->x-A3->x) + (A1->y-A3->y)*(A1->y-A3->y));
+ }
+
+ /* Closed circle! Return the circumference! */
+ circumference_A = M_PI * 2 * radius_A;
+ if ( p2d_same(A1, A3) )
+ return circumference_A;
+
+ /* Determine the orientation of the arc */
+ a2_side = lw_segment_side(A1, A3, A2);
+
+ /* The side of the A1/A3 line that A2 falls on dictates the sweep
+ direction from A1 to A3. */
+ if ( a2_side == -1 )
+ clockwise = LW_TRUE;
+ else
+ clockwise = LW_FALSE;
+
+ /* Angles of each point that defines the arc section */
+ a1 = atan2(A1->y - C.y, A1->x - C.x);
+ a2 = atan2(A2->y - C.y, A2->x - C.x);
+ a3 = atan2(A3->y - C.y, A3->x - C.x);
+
+ /* What's the sweep from A1 to A3? */
+ if ( clockwise )
+ {
+ if ( a1 > a3 )
+ angle = a1 - a3;
+ else
+ angle = 2*M_PI + a1 - a3;
+ }
else
- return side;
+ {
+ if ( a3 > a1 )
+ angle = a3 - a1;
+ else
+ angle = 2*M_PI + a3 - a1;
+ }
+
+ /* Length as proportion of circumference */
+ return circumference_A* (angle / (2*M_PI));
}
+double lw_arc_side(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, const POINT2D *Q)
+{
+ POINT2D C;
+ double radius_A;
+ double side_Q, side_A2;
+ double d;
+
+ side_Q = lw_segment_side(A1, A3, Q);
+ radius_A = lwcircle_center(A1, A2, A3, &C);
+ side_A2 = lw_segment_side(A1, A3, A2);
+
+ /* Linear case */
+ if ( radius_A < 0 )
+ return side_Q;
+
+ d = distance2d_pt_pt(Q, &C);
+
+ if ( d < radius_A && side_Q == side_A2 )
+ {
+ side_Q *= -1;
+ }
+
+ return side_Q;
+}
+
+/**
+* Determines the center of the circle defined by the three given points.
+* In the event the circle is complete, the midpoint of the segment defined
+* by the first and second points is returned. If the points are colinear,
+* as determined by equal slopes, then NULL is returned. If the interior
+* point is coincident with either end point, they are taken as colinear.
+*/
+double
+lwcircle_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result)
+{
+ POINT2D c;
+ double cx, cy, cr;
+ double temp, bc, cd, det;
+
+ c.x = c.y = 0.0;
+
+ LWDEBUGF(2, "lwcircle_center called (%.16f,%.16f), (%.16f,%.16f), (%.16f,%.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y);
+
+ /* Closed circle */
+ if (fabs(p1->x - p3->x) < EPSILON_SQLMM &&
+ fabs(p1->y - p3->y) < EPSILON_SQLMM)
+ {
+ cx = p1->x + (p2->x - p1->x) / 2.0;
+ cy = p1->y + (p2->y - p1->y) / 2.0;
+ c.x = cx;
+ c.y = cy;
+ *result = c;
+ cr = sqrt(pow(cx - p1->x, 2.0) + pow(cy - p1->y, 2.0));
+ return cr;
+ }
+
+ temp = p2->x*p2->x + p2->y*p2->y;
+ bc = (p1->x*p1->x + p1->y*p1->y - temp) / 2.0;
+ cd = (temp - p3->x*p3->x - p3->y*p3->y) / 2.0;
+ det = (p1->x - p2->x)*(p2->y - p3->y)-(p2->x - p3->x)*(p1->y - p2->y);
+
+ /* Check colinearity */
+ if (fabs(det) < EPSILON_SQLMM)
+ return -1.0;
+
+
+ det = 1.0 / det;
+ cx = (bc*(p2->y - p3->y)-cd*(p1->y - p2->y))*det;
+ cy = ((p1->x - p2->x)*cd-(p2->x - p3->x)*bc)*det;
+ c.x = cx;
+ c.y = cy;
+ *result = c;
+ cr = sqrt((cx-p1->x)*(cx-p1->x)+(cy-p1->y)*(cy-p1->y));
+
+ LWDEBUGF(2, "lwcircle_center center is (%.16f,%.16f)", result->x, result->y);
+
+ return cr;
+}
+
+
+int
+pt_in_ring_2d(const POINT2D *p, const POINTARRAY *ring)
+{
+ int cn = 0; /* the crossing number counter */
+ int i;
+ POINT2D v1, v2;
+
+ POINT2D first, last;
+
+ getPoint2d_p(ring, 0, &first);
+ getPoint2d_p(ring, ring->npoints-1, &last);
+ if ( memcmp(&first, &last, sizeof(POINT2D)) )
+ {
+ lwerror("pt_in_ring_2d: V[n] != V[0] (%g %g != %g %g)",
+ first.x, first.y, last.x, last.y);
+ return LW_FALSE;
+
+ }
+
+ LWDEBUGF(2, "pt_in_ring_2d called with point: %g %g", p->x, p->y);
+ /* printPA(ring); */
+
+ /* loop through all edges of the polygon */
+ getPoint2d_p(ring, 0, &v1);
+ for (i=0; i<ring->npoints-1; i++)
+ {
+ double vt;
+ getPoint2d_p(ring, i+1, &v2);
+
+ /* edge from vertex i to vertex i+1 */
+ if
+ (
+ /* an upward crossing */
+ ((v1.y <= p->y) && (v2.y > p->y))
+ /* a downward crossing */
+ || ((v1.y > p->y) && (v2.y <= p->y))
+ )
+ {
+
+ vt = (double)(p->y - v1.y) / (v2.y - v1.y);
+
+ /* P.x <intersect */
+ if (p->x < v1.x + vt * (v2.x - v1.x))
+ {
+ /* a valid crossing of y=p.y right of p.x */
+ ++cn;
+ }
+ }
+ v1 = v2;
+ }
+
+ LWDEBUGF(3, "pt_in_ring_2d returning %d", cn&1);
+
+ return (cn&1); /* 0 if even (out), and 1 if odd (in) */
+}
-int lw_segment_envelope_intersects(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2)
+static int
+lw_seg_interact(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2)
{
double minq=FP_MIN(q1->x,q2->x);
double maxq=FP_MAX(q1->x,q2->x);
int lw_segment_intersects(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2)
{
- double pq1, pq2, qp1, qp2;
+ int pq1, pq2, qp1, qp2;
/* No envelope interaction => we are done. */
- if (!lw_segment_envelope_intersects(p1, p2, q1, p2))
+ if (!lw_seg_interact(p1, p2, q1, p2))
{
return SEG_NO_INTERSECTION;
}
LWDEBUGF(4, "qp1=%.15g qp2=%.15g", qp1, qp2);
/* Second point of p or q touches, it's not a crossing. */
- if ( pq2 == 0.0 || qp2 == 0.0 )
+ if ( pq2 == 0 || qp2 == 0 )
{
return SEG_NO_INTERSECTION;
}
/* First point of p touches, it's a "crossing". */
- if ( pq1 == 0.0 )
+ if ( pq1 == 0 )
{
- if ( FP_GT(pq2,0.0) )
+ if ( pq2 > 0 )
return SEG_CROSS_RIGHT;
else
return SEG_CROSS_LEFT;
}
/* First point of q touches, it's a crossing. */
- if ( qp1 == 0.0 )
+ if ( qp1 == 0 )
{
- if ( FP_LT(pq1,pq2) )
+ if ( pq1 < pq2 )
return SEG_CROSS_RIGHT;
else
return SEG_CROSS_LEFT;
}
/* The segments cross, what direction is the crossing? */
- if ( FP_LT(pq1,pq2) )
+ if ( pq1 < pq2 )
return SEG_CROSS_RIGHT;
else
return SEG_CROSS_LEFT;
LWGEOM *lwgeom_desegmentize(LWGEOM *geom);
-/*
- * Tolerance used to determine equality.
- */
-#define EPSILON_SQLMM 1e-8
-
/*
* Determines (recursively in the case of collections) whether the geometry
* contains at least on arc geometry or segment.
}
}
-/*
- * Determines the center of the circle defined by the three given points.
- * In the event the circle is complete, the midpoint of the segment defined
- * by the first and second points is returned. If the points are colinear,
- * as determined by equal slopes, then NULL is returned. If the interior
- * point is coincident with either end point, they are taken as colinear.
- */
-double
-lwcircle_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result)
-{
- POINT2D c;
- double cx, cy, cr;
- double temp, bc, cd, det;
-
- c.x = c.y = 0.0;
-
- LWDEBUGF(2, "lwcircle_center called (%.16f,%.16f), (%.16f,%.16f), (%.16f,%.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y);
-
- /* Closed circle */
- if (fabs(p1->x - p3->x) < EPSILON_SQLMM &&
- fabs(p1->y - p3->y) < EPSILON_SQLMM)
- {
- cx = p1->x + (p2->x - p1->x) / 2.0;
- cy = p1->y + (p2->y - p1->y) / 2.0;
- c.x = cx;
- c.y = cy;
- *result = c;
- cr = sqrt(pow(cx - p1->x, 2.0) + pow(cy - p1->y, 2.0));
- return cr;
- }
-
- temp = p2->x*p2->x + p2->y*p2->y;
- bc = (p1->x*p1->x + p1->y*p1->y - temp) / 2.0;
- cd = (temp - p3->x*p3->x - p3->y*p3->y) / 2.0;
- det = (p1->x - p2->x)*(p2->y - p3->y)-(p2->x - p3->x)*(p1->y - p2->y);
-
- /* Check colinearity */
- if (fabs(det) < EPSILON_SQLMM)
- return -1.0;
-
-
- det = 1.0 / det;
- cx = (bc*(p2->y - p3->y)-cd*(p1->y - p2->y))*det;
- cy = ((p1->x - p2->x)*cd-(p2->x - p3->x)*bc)*det;
- c.x = cx;
- c.y = cy;
- *result = c;
- cr = sqrt((cx-p1->x)*(cx-p1->x)+(cy-p1->y)*(cy-p1->y));
-
- LWDEBUGF(2, "lwcircle_center center is (%.16f,%.16f)", result->x, result->y);
-
- return cr;
-}
/*******************************************************************************
LWDEBUG(2, "lwcircle_calculate_gbox called.");
radius = lwcircle_center((POINT2D*)p1, (POINT2D*)p2, (POINT2D*)p3, ¢er);
- p2_side = signum(lw_segment_side((POINT2D*)p1, (POINT2D*)p3, (POINT2D*)p2));
+ p2_side = lw_segment_side((POINT2D*)p1, (POINT2D*)p3, (POINT2D*)p2);
/* Matched start/end points imply circle */
if ( p1->x == p3->x && p1->y == p3->y )
static int pt_continues_arc(const POINT4D *a1, const POINT4D *a2, const POINT4D *a3, const POINT4D *b)
{
POINT2D center;
- double radius = lwcircle_center((POINT2D*)a1, (POINT2D*)a2, (POINT2D*)a3, ¢er);
+ POINT2D *t1 = (POINT2D*)a1;
+ POINT2D *t2 = (POINT2D*)a2;
+ POINT2D *t3 = (POINT2D*)a3;
+ POINT2D *tb = (POINT2D*)b;
+ double radius = lwcircle_center(t1, t2, t3, ¢er);
double b_distance, diff;
/* Co-linear a1/a2/a3 */
if ( radius < 0.0 )
return LW_FALSE;
- b_distance = distance2d_pt_pt((POINT2D*)b, ¢er);
+ b_distance = distance2d_pt_pt(tb, ¢er);
diff = fabs(radius - b_distance);
LWDEBUGF(4, "circle_radius=%g, b_distance=%g, diff=%g, percentage=%g", radius, b_distance, diff, diff/radius);
/* Is the point b on the circle? */
if ( diff < EPSILON_SQLMM )
{
- int a2_side = signum(lw_segment_side((POINT2D*)a1, (POINT2D*)a3, (POINT2D*)a2));
- int b_side = signum(lw_segment_side((POINT2D*)a1, (POINT2D*)a3, (POINT2D*)b));
+ int a2_side = lw_segment_side(t1, t3, t2);
+ int b_side = lw_segment_side(t1, t3, tb);
/* Is the point b on the same side of a1/a3 as the mid-point a2 is? */
/* If not, it's in the unbounded part of the circle, so it continues the arc, return true. */
if ( rect_node_is_leaf(node) )
{
double side = lw_segment_side(node->p1, node->p2, pt);
- if ( side == 0.0 )
+ if ( side == 0 )
*on_boundary = LW_TRUE;
- return (side < 0.0 ? -1 : 1 );
+ return (side < 0 ? -1 : 1 );
}
else
{
return LW_TRUE;
}
-/**
-* Returns true if P is on the same side of the plane partition
-* defined by A1/A3 as A2 is. Only makes sense if P has already been
-* determined to be on the circle defined by A1/A2/A3.
-*/
-static int
-lw_pt_in_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
-{
- return signum(lw_segment_side(A1, A3, A2)) == signum(lw_segment_side(A1, A3, P));
-}
-
-/**
-* Returns true if P is between A1/A2. Only makes sense if P has already been
-* deterined to be on the line defined by A1/A2.
-*/
-int
-lw_pt_in_seg(const POINT2D *P, const POINT2D *A1, const POINT2D *A2)
-{
- return ((A1->x <= P->x && P->x < A2->x) || (A1->x >= P->x && P->x > A2->x)) ||
- ((A1->y <= P->y && P->y < A2->y) || (A1->y >= P->y && P->y > A2->y));
-}
-
-/**
-* Returns true if arc A is actually a point (all vertices are the same) .
-*/
-static int
-lw_arc_is_pt(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
-{
- if ( A1->x == A2->x && A2->x == A3->x &&
- A1->y == A2->y && A2->y == A3->y )
- return LW_TRUE;
- else
- return LW_FALSE;
-}
-
-/**
-* Returns the length of a circular arc segment
-*/
-double
-lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
-{
- POINT2D C;
- double radius_A, circumference_A;
- int a2_side, clockwise;
- double a1, a2, a3;
- double angle;
-
- if ( lw_arc_is_pt(A1, A2, A3) )
- return 0.0;
-
- radius_A = lwcircle_center(A1, A2, A3, &C);
-
- /* Co-linear! Return linear distance! */
- if ( radius_A < 0 )
- {
- return sqrt((A1->x-A3->x)*(A1->x-A3->x) + (A1->y-A3->y)*(A1->y-A3->y));
- }
-
- /* Closed circle! Return the circumference! */
- circumference_A = M_PI * 2 * radius_A;
- if ( p2d_same(A1, A3) )
- return circumference_A;
-
- /* Determine the orientation of the arc */
- a2_side = signum(lw_segment_side(A1, A3, A2));
-
- /* The side of the A1/A3 line that A2 falls on dictates the sweep
- direction from A1 to A3. */
- if ( a2_side == -1 )
- clockwise = LW_TRUE;
- else
- clockwise = LW_FALSE;
-
- /* Angles of each point that defines the arc section */
- a1 = atan2(A1->y - C.y, A1->x - C.x);
- a2 = atan2(A2->y - C.y, A2->x - C.x);
- a3 = atan2(A3->y - C.y, A3->x - C.x);
-
- /* What's the sweep from A1 to A3? */
- if ( clockwise )
- {
- if ( a1 > a3 )
- angle = a1 - a3;
- else
- angle = 2*M_PI + a1 - a3;
- }
- else
- {
- if ( a3 > a1 )
- angle = a3 - a1;
- else
- angle = 2*M_PI + a3 - a1;
- }
-
- /* Length as proportion of circumference */
- return circumference_A* (angle / (2*M_PI));
-}
-
/**
* Calculate the shortest distance between an arc and an edge.
* Line/circle approach from http://stackoverflow.com/questions/1073336/circle-line-collision-detection
* Our polygons have first and last point the same,
*
*/
-int
-pt_in_ring_2d(const POINT2D *p, const POINTARRAY *ring)
-{
- int cn = 0; /* the crossing number counter */
- int i;
- POINT2D v1, v2;
-
- POINT2D first, last;
-
- getPoint2d_p(ring, 0, &first);
- getPoint2d_p(ring, ring->npoints-1, &last);
- if ( memcmp(&first, &last, sizeof(POINT2D)) )
- {
- lwerror("pt_in_ring_2d: V[n] != V[0] (%g %g != %g %g)",
- first.x, first.y, last.x, last.y);
- return LW_FALSE;
-
- }
-
- LWDEBUGF(2, "pt_in_ring_2d called with point: %g %g", p->x, p->y);
- /* printPA(ring); */
-
- /* loop through all edges of the polygon */
- getPoint2d_p(ring, 0, &v1);
- for (i=0; i<ring->npoints-1; i++)
- {
- double vt;
- getPoint2d_p(ring, i+1, &v2);
-
- /* edge from vertex i to vertex i+1 */
- if
- (
- /* an upward crossing */
- ((v1.y <= p->y) && (v2.y > p->y))
- /* a downward crossing */
- || ((v1.y > p->y) && (v2.y <= p->y))
- )
- {
-
- vt = (double)(p->y - v1.y) / (v2.y - v1.y);
-
- /* P.x <intersect */
- if (p->x < v1.x + vt * (v2.x - v1.x))
- {
- /* a valid crossing of y=p.y right of p.x */
- ++cn;
- }
- }
- v1 = v2;
- }
-
- LWDEBUGF(3, "pt_in_ring_2d returning %d", cn&1);
-
- return (cn&1); /* 0 if even (out), and 1 if odd (in) */
-}
/**
/* Zero length segments are ignored. */
if ( seg1->x == seg2->x && seg1->y == seg2->y )
+ {
+ seg1 = seg2;
continue;
+ }
ymin = FP_MIN(seg1->y, seg2->y);
ymax = FP_MAX(seg1->y, seg2->y);
continue;
}
+ side = lw_segment_side(seg1, seg2, pt);
+
+ /*
+ * A point on the boundary of a ring is not contained.
+ * WAS: if (fabs(side) < 1e-12), see #852
+ */
+ if ( (side == 0) && lw_pt_in_seg(pt, seg1, seg2) )
+ {
+ return 0;
+ }
+
+ /*
+ * If the point is to the left of the line, and it's rising,
+ * then the line is to the right of the point and
+ * circling counter-clockwise, so incremement.
+ */
+ if ( (side < 0) && (seg1->y <= pt->y) && (pt->y < seg2->y) )
+ {
+ wn++;
+ }
+
+ /*
+ * If the point is to the right of the line, and it's falling,
+ * then the line is to the right of the point and circling
+ * clockwise, so decrement.
+ */
+ else if ( (side > 0) && (seg2->y <= pt->y) && (pt->y < seg1->y) )
+ {
+ wn--;
+ }
+
+ seg1 = seg2;
+ }
+
+ /* Outside */
+ if (wn == 0)
+ {
+ return -1;
+ }
+
+ /* Inside */
+ return 1;
+}
+
+/**
+* Return -1 if the point is inside the POINTARRAY, 1 if it is outside,
+* and 0 if it is on the boundary.
+*/
+int
+ptarray_arc_contains_point(const POINTARRAY *pa, const POINT2D *pt)
+{
+ int wn = 0;
+ int i;
+ double side;
+ const POINT2D *seg1;
+ const POINT2D *seg2;
+ const POINT2D *seg3;
+ GBOX gbox;
+
+ /* Check for not an arc ring (always have odd # of points) */
+ if ( pa->npoints % 2 )
+ lwerror("ptarray_arc_contains_point called with even number of points");
+
+ /* Check for unclosed */
+ seg1 = getPoint2d_cp(pa, 0);
+ seg2 = getPoint2d_cp(pa, pa->npoints-1);
+
+ if ( ! p2d_same(seg1, seg2) )
+ lwerror("ptarray_contains_point called on unclosed ring");
+
+ if ( pa->npoints == 3 && p2d_same(seg1, seg2) )
+ {
+ /* TODO handle circlular case */
+ }
+
+ /* Start on the ring */
+ seg2 = getPoint2d_cp(pa, 1);
+ for ( i=2; i < pa->npoints; i++ )
+ {
+ seg3 = getPoint2d_cp(pa, i);
+
+ if ( lw_arc_is_pt(seg1, seg2, seg3) )
+ {
+ seg1 = seg2;
+ seg2 = seg3;
+ continue;
+ }
+
+ lwcircle_calculate_gbox_cartesian_2d(seg1, seg2, seg3, &gbox);
+
+ /* Only test segments in our vertical range */
+ if ( pt->y > gbox.ymax || pt->y < gbox.ymin )
+ {
+ seg1 = seg2;
+ seg2 = seg3;
+ continue;
+ }
+
+ side = lw_segment_side(seg1, seg3, pt);
+
+
+ /* Going "up"! */
+ if ( (seg1->y <= pt->y) && (pt->y < seg3->y) )
+ {
+ }
+
+ /* Going "down"! */
+ if ( (seg2->y <= pt->y) && (pt->y < seg1->y) )
+ {
+ }
+
+
+
side = lw_segment_side(seg1, seg2, pt);
/*
}
seg1 = seg2;
+ seg2 = seg3;
}
/* Outside */
projects / postgis / commitdiff