From: Paul Ramsey <pramsey@cleverelephant.ca>
Date: Mon, 1 Oct 2012 22:10:33 +0000 (+0000)
Subject: Anal retentive code re-organization. Try and move the primitive computational geometr... 
X-Git-Tag: 2.1.0beta2~601
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=e92cf6071d6e9ca36a4c1c5f43a4c0a6fc80ee6d;p=postgis

Anal retentive code re-organization. Try and move the primitive computational geometry functions into lwalgorithm.c


git-svn-id: http://svn.osgeo.org/postgis/trunk@10356 b70326c6-7e19-0410-871a-916f4a2858ee
---

diff --git a/liblwgeom/cunit/cu_algorithm.c b/liblwgeom/cunit/cu_algorithm.c
index 774d369f1..c24b09e70 100644
--- a/liblwgeom/cunit/cu_algorithm.c
+++ b/liblwgeom/cunit/cu_algorithm.c
@@ -61,7 +61,7 @@ static int clean_cg_suite(void)
 */
 static void test_lw_segment_side(void)
 {
-	double rv = 0.0;
+	int rv = 0;
 	POINT2D p1, p2, q;
 
 	/* Vertical line at x=0 */
@@ -75,19 +75,19 @@ static void test_lw_segment_side(void)
 	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);
 
 }
 
diff --git a/liblwgeom/g_box.c b/liblwgeom/g_box.c
index dc32122fd..8c4d61ab3 100644
--- a/liblwgeom/g_box.c
+++ b/liblwgeom/g_box.c
@@ -321,85 +321,89 @@ size_t gbox_serialized_size(uint8_t flags)
 ** 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, &center);
-	
 	/* 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;
diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
index 20e5e864e..ca85d8345 100644
--- a/liblwgeom/liblwgeom.h.in
+++ b/liblwgeom/liblwgeom.h.in
@@ -1809,7 +1809,6 @@ extern char *lwmessage_truncate(char *str, int startpos, int endpos, int maxleng
  ******************************************************************************/
 
 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);
 
diff --git a/liblwgeom/liblwgeom_internal.h b/liblwgeom/liblwgeom_internal.h
index df90ead34..6203ad37c 100644
--- a/liblwgeom/liblwgeom_internal.h
+++ b/liblwgeom/liblwgeom_internal.h
@@ -121,6 +121,10 @@
 #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
@@ -203,17 +207,6 @@ enum CG_SEGMENT_INTERSECTION_TYPE {
 */
 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)
 */
@@ -354,8 +347,22 @@ int lwtin_is_closed(const LWTIN *tin);
 
 
 
+/*
+* 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.
diff --git a/liblwgeom/lwalgorithm.c b/liblwgeom/lwalgorithm.c
index 5b0dc465c..531d282b7 100644
--- a/liblwgeom/lwalgorithm.c
+++ b/liblwgeom/lwalgorithm.c
@@ -27,22 +27,265 @@ int signum(double n)
 /**
 * 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);
@@ -80,10 +323,10 @@ int lw_segment_envelope_intersects(const POINT2D *p1, const POINT2D *p2, const P
 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;
 	}
@@ -119,31 +362,31 @@ int lw_segment_intersects(const POINT2D *p1, const POINT2D *p2, const POINT2D *q
 	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;
diff --git a/liblwgeom/lwsegmentize.c b/liblwgeom/lwsegmentize.c
index 1019d5251..aa1c008ee 100644
--- a/liblwgeom/lwsegmentize.c
+++ b/liblwgeom/lwsegmentize.c
@@ -31,11 +31,6 @@ LWGEOM *lwmpolygon_desegmentize(LWMPOLY *mpoly);
 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.
@@ -74,59 +69,6 @@ lwgeom_has_arc(const LWGEOM *geom)
 	}
 }
 
-/*
- * 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;
-}
 
 
 /*******************************************************************************
@@ -171,7 +113,7 @@ lwcircle_segmentize(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32_t perQuad)
 	LWDEBUG(2, "lwcircle_calculate_gbox called.");
 
 	radius = lwcircle_center((POINT2D*)p1, (POINT2D*)p2, (POINT2D*)p3, &center);
-	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 )
@@ -528,22 +470,26 @@ lwgeom_segmentize(LWGEOM *geom, uint32_t perQuad)
 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, &center);
+	POINT2D *t1 = (POINT2D*)a1;
+	POINT2D *t2 = (POINT2D*)a2;
+	POINT2D *t3 = (POINT2D*)a3;
+	POINT2D *tb = (POINT2D*)b;
+	double radius = lwcircle_center(t1, t2, t3, &center);
 	double b_distance, diff;
 
 	/* Co-linear a1/a2/a3 */
 	if ( radius < 0.0 )
 		return LW_FALSE;
 
-	b_distance = distance2d_pt_pt((POINT2D*)b, &center);
+	b_distance = distance2d_pt_pt(tb, &center);
 	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. */
diff --git a/liblwgeom/lwtree.c b/liblwgeom/lwtree.c
index 30b51acb1..edb5d55bf 100644
--- a/liblwgeom/lwtree.c
+++ b/liblwgeom/lwtree.c
@@ -38,9 +38,9 @@ int rect_tree_contains_point(const RECT_NODE *node, const POINT2D *pt, int *on_b
 		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
 		{
diff --git a/liblwgeom/measures.c b/liblwgeom/measures.c
index 4bcd29fa5..053abef08 100644
--- a/liblwgeom/measures.c
+++ b/liblwgeom/measures.c
@@ -825,104 +825,6 @@ lw_dist2d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2,DISTPTS *dl)
 	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 
@@ -1687,61 +1589,6 @@ Functions in common for Brute force and new calculation
  *	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) */
-}
 
 
 /**
diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index e50270252..9b95c7d54 100644
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -714,7 +714,10 @@ ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt)
 		
 		/* 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);
@@ -726,6 +729,119 @@ ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt)
 			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);
 
 		/* 
@@ -758,6 +874,7 @@ ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt)
 		}
 		
 		seg1 = seg2;
+		seg2 = seg3;
 	}
 
 	/* Outside */