[SCM] PostGIS branch stable-3.6 updated. 3.6.0-11-gbe2a0ffda
git at osgeo.org
git at osgeo.org
Thu Oct 2 14:08:49 PDT 2025
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "PostGIS".
The branch, stable-3.6 has been updated
via be2a0ffda5b08dca0d94c1498d84b62b83fe76f9 (commit)
from d3f3a8bf9e23c32beb28a04c490930fd82cba95d (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit be2a0ffda5b08dca0d94c1498d84b62b83fe76f9
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Thu Oct 2 14:08:21 2025 -0700
Fix arc/arc distance to handle more cases such as contained,
circular, and other interesting arc/arc alignments that come up rarely
in ordinary GIS data.
Fix compoundcurve distance calculations, in particular the
containment test, which was failing for some cases.
Replaced containment test with a ray casting algorithm
that makes more sense for circularstrings that the
winding order algorithm did.
References #5989
diff --git a/NEWS b/NEWS
index f65470d6f..98e536d63 100644
--- a/NEWS
+++ b/NEWS
@@ -8,6 +8,7 @@ PostGIS 3.6.1
- #5987, ST_GeometryN fails for non-collections (Paul Ramsey)
- #5991, CircularString distance error (Paul Ramsey)
- #5994, Null pointer in ST_AsGeoJsonRow (Alexander Kukushkin)
+ - #5998, ST_Distance error on CurvePolygon (Paul Ramsey)
PostGIS 3.6.0
diff --git a/liblwgeom/cunit/cu_measures.c b/liblwgeom/cunit/cu_measures.c
index 610b2a593..58e7cb4ba 100644
--- a/liblwgeom/cunit/cu_measures.c
+++ b/liblwgeom/cunit/cu_measures.c
@@ -222,6 +222,11 @@ static void test_mindistance2d_tolerance(void)
DIST2DTEST(
"CURVEPOLYGON(CIRCULARSTRING(7874821 8715927,8907663 8715927,8844683 7750316,7937800 7750316,7874821 8715927))",
"POINT(5433865 8243495)", 2271704.2698450615, default_accepted_error);
+
+ /* Ticket 5989 */
+ DIST2DTEST(
+ "CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(0 0, -1 5, 0 10), (0 10, -10 10, -10 0, 0 0)))",
+ "POINT(-0.5 5)", 0.5, default_accepted_error);
}
static void
@@ -918,6 +923,25 @@ test_lw_dist2d_arc_arc(void)
B2.x = 0 ; B2.y = 1;
B3.x = 1 ; B3.y = 0;
+ /* Arc inside the semicircle */
+ lw_dist2d_distpts_init(&dl, DIST_MIN);
+ A1.x = 0; A1.y = .5;
+ A2.x = -0.3; A2.y = .5;
+ A3.x = 0; A3.y = .6;
+ rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
+ printf("%g\n", dl.distance);
+ CU_ASSERT_EQUAL( rv, LW_SUCCESS );
+ CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.271798, 0.000001);
+
+ /* Arc inside the semicircle */
+ lw_dist2d_distpts_init(&dl, DIST_MIN);
+ A1.x = -0.5; A1.y = .5;
+ A2.x = -0.4; A2.y = .2;
+ A3.x = 0; A3.y = 0;
+ rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
+ CU_ASSERT_EQUAL( rv, LW_SUCCESS );
+ CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.292893, 0.000001);
+
/* Arc above the unit semicircle */
lw_dist2d_distpts_init(&dl, DIST_MIN);
A1.x = -1; A1.y = 3;
@@ -954,6 +978,15 @@ test_lw_dist2d_arc_arc(void)
CU_ASSERT_EQUAL( rv, LW_SUCCESS );
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
+ /* Closed circle */
+ lw_dist2d_distpts_init(&dl, DIST_MIN);
+ A1.x = -2.0; A1.y = -0.1;
+ A2.x = 1.5; A2.y = -0.1;
+ A3.x = -2.0; A3.y = -0.1;
+ rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
+ CU_ASSERT_EQUAL( rv, LW_SUCCESS );
+ CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.480742, 0.0001);
+
/* Concentric: and fully parallel */
lw_dist2d_distpts_init(&dl, DIST_MIN);
A1.x = -2.0; A1.y = 0.0;
diff --git a/liblwgeom/cunit/cu_ptarray.c b/liblwgeom/cunit/cu_ptarray.c
index a28da5cce..99c689ac2 100644
--- a/liblwgeom/cunit/cu_ptarray.c
+++ b/liblwgeom/cunit/cu_ptarray.c
@@ -433,14 +433,14 @@ static void test_ptarrayarc_contains_point()
{
/* int ptarrayarc_contains_point(const POINTARRAY *pa, const POINT2D *pt) */
- LWLINE *lwline;
+ LWCIRCSTRING *lwcirc;
POINTARRAY *pa;
POINT2D pt;
int rv;
/*** Collection of semi-circles surrounding unit square ***/
- lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 -1, -2 0, -1 1, 0 2, 1 1, 2 0, 1 -1, 0 -2, -1 -1)"));
- pa = lwline->points;
+ lwcirc = lwgeom_as_lwcircstring(lwgeom_from_text("CIRCULARSTRING(-1 -1, -2 0, -1 1, 0 2, 1 1, 2 0, 1 -1, 0 -2, -1 -1)"));
+ pa = lwcirc->points;
/* Point in middle of square */
pt.x = 0;
@@ -473,9 +473,9 @@ static void test_ptarrayarc_contains_point()
CU_ASSERT_EQUAL(rv, LW_OUTSIDE);
/*** Two-edge ring made up of semi-circles (really, a circle) ***/
- lwline_free(lwline);
- lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0, 0 -1, -1 0)"));
- pa = lwline->points;
+ lwcircstring_free(lwcirc);
+ lwcirc = lwgeom_as_lwcircstring(lwgeom_from_text("CIRCULARSTRING(-1 0, 0 1, 1 0, 0 -1, -1 0)"));
+ pa = lwcirc->points;
/* Point outside */
pt.x = -1.5;
@@ -514,9 +514,9 @@ static void test_ptarrayarc_contains_point()
CU_ASSERT_EQUAL(rv, LW_BOUNDARY);
/*** Two-edge ring, closed ***/
- lwline_free(lwline);
- lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(1 6, 6 1, 9 7, 6 10, 1 6)"));
- pa = lwline->points;
+ lwcircstring_free(lwcirc);
+ lwcirc = lwgeom_as_lwcircstring(lwgeom_from_text("CIRCULARSTRING(1 6, 6 1, 9 7, 6 10, 1 6)"));
+ pa = lwcirc->points;
/* Point to left of ring */
pt.x = 20;
@@ -525,9 +525,9 @@ static void test_ptarrayarc_contains_point()
CU_ASSERT_EQUAL(rv, LW_OUTSIDE);
/*** One-edge ring, closed circle ***/
- lwline_free(lwline);
- lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 1 0, -1 0)"));
- pa = lwline->points;
+ lwcircstring_free(lwcirc);
+ lwcirc = lwgeom_as_lwcircstring(lwgeom_from_text("CIRCULARSTRING(-1 0, 1 0, -1 0)"));
+ pa = lwcirc->points;
/* Point inside */
pt.x = 0;
@@ -548,23 +548,23 @@ static void test_ptarrayarc_contains_point()
CU_ASSERT_EQUAL(rv, LW_BOUNDARY);
/*** Overshort ring ***/
- lwline_free(lwline);
- lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 1 0)"));
- pa = lwline->points;
+ lwcircstring_free(lwcirc);
+ lwcirc = lwgeom_as_lwcircstring(lwgeom_from_text("CIRCULARSTRING(-1 0, 1 0)"));
+ pa = lwcirc->points;
cu_error_msg_reset();
rv = ptarrayarc_contains_point(pa, &pt);
//printf("%s\n", cu_error_msg);
- ASSERT_STRING_EQUAL("ptarrayarc_contains_point called with even number of points", cu_error_msg);
+ ASSERT_STRING_EQUAL(cu_error_msg, "ptarrayarc_raycast_intersections called with even number of points");
/*** Unclosed ring ***/
- lwline_free(lwline);
- lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 1 0, 2 0)"));
- pa = lwline->points;
+ lwcircstring_free(lwcirc);
+ lwcirc = lwgeom_as_lwcircstring(lwgeom_from_text("CIRCULARSTRING(-1 0, 1 0, 2 0)"));
+ pa = lwcirc->points;
cu_error_msg_reset();
rv = ptarrayarc_contains_point(pa, &pt);
- ASSERT_STRING_EQUAL("ptarrayarc_contains_point called on unclosed ring", cu_error_msg);
+ ASSERT_STRING_EQUAL(cu_error_msg, "ptarrayarc_contains_point called on unclosed ring");
- lwline_free(lwline);
+ lwcircstring_free(lwcirc);
}
static void test_ptarray_scale()
diff --git a/liblwgeom/liblwgeom_internal.h b/liblwgeom/liblwgeom_internal.h
index e51062b00..04a3504e2 100644
--- a/liblwgeom/liblwgeom_internal.h
+++ b/liblwgeom/liblwgeom_internal.h
@@ -445,13 +445,15 @@ double lw_arc_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, PO
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);
+int lw_pt_on_segment(const POINT2D* p1, const POINT2D* p2, const POINT2D* p);
double lw_seg_length(const POINT2D *A1, const POINT2D *A2);
double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3);
int pt_in_ring_2d(const POINT2D *p, const POINTARRAY *ring);
int ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt);
int ptarrayarc_contains_point(const POINTARRAY *pa, const POINT2D *pt);
int ptarray_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number);
-int ptarrayarc_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number);
+int ptarray_raycast_intersections(const POINTARRAY *pa, const POINT2D *p, int *on_boundary);
+int ptarrayarc_raycast_intersections(const POINTARRAY *pa, const POINT2D *p, int *on_boundary);
int lwcompound_contains_point(const LWCOMPOUND *comp, const POINT2D *pt);
int lwgeom_contains_point(const LWGEOM *geom, const POINT2D *pt);
diff --git a/liblwgeom/lwalgorithm.c b/liblwgeom/lwalgorithm.c
index cb3727f15..a2cb2b037 100644
--- a/liblwgeom/lwalgorithm.c
+++ b/liblwgeom/lwalgorithm.c
@@ -90,7 +90,9 @@ lw_seg_length(const POINT2D *A1, const POINT2D *A2)
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);
+ int pt_side = lw_segment_side(A1, A3, P);
+ int arc_side = lw_segment_side(A1, A3, A2);
+ return pt_side == 0 || pt_side == arc_side;
}
/**
@@ -104,6 +106,28 @@ lw_pt_in_seg(const POINT2D *P, const POINT2D *A1, const POINT2D *A2)
((A1->y <= P->y && P->y < A2->y) || (A1->y >= P->y && P->y > A2->y));
}
+int
+lw_pt_on_segment(const POINT2D* p1, const POINT2D* p2, const POINT2D* p)
+{
+ // Check if the point is within the bounding box of the segment
+ if (p->x < FP_MIN(p1->x, p2->x) || p->x > FP_MAX(p1->x, p2->x) ||
+ p->y < FP_MIN(p1->y, p2->y) || p->y > FP_MAX(p1->y, p2->y))
+ {
+ return LW_FALSE;
+ }
+
+ // Check for collinearity using the 2D cross-product.
+ // (p.y - p1.y) * (p2.x - p1.x) - (p.x - p1.x) * (p2.y - p1.y)
+ double cross_product = (p->y - p1->y) * (p2->x - p1->x) - (p->x - p1->x) * (p2->y - p1->y);
+
+ if (fabs(cross_product) < DBL_EPSILON)
+ {
+ return LW_TRUE; // Point is collinear and within the bounding box
+ }
+
+ return LW_TRUE;
+}
+
/**
* Returns true if arc A is actually a point (all vertices are the same) .
*/
diff --git a/liblwgeom/lwcompound.c b/liblwgeom/lwcompound.c
index 90fee3560..1e0d5c2ee 100644
--- a/liblwgeom/lwcompound.c
+++ b/liblwgeom/lwcompound.c
@@ -184,63 +184,48 @@ int lwgeom_contains_point(const LWGEOM *geom, const POINT2D *pt)
return LW_FAILURE;
}
+/*
+ * Use a ray-casting count to determine if the point
+ * is inside or outside of the compound curve. Ray-casting
+ * is run against each component of the overall arc, and
+ * the even/odd test run against the total of all components.
+ * Returns LW_INSIDE / LW_BOUNDARY / LW_OUTSIDE
+ */
int
lwcompound_contains_point(const LWCOMPOUND *comp, const POINT2D *pt)
{
- uint32_t i;
- LWLINE *lwline;
- LWCIRCSTRING *lwcirc;
- int wn = 0;
- int winding_number = 0;
- int result;
+ int intersections = 0;
- for ( i = 0; i < comp->ngeoms; i++ )
+ if (lwgeom_is_empty(lwcompound_as_lwgeom(comp)))
+ return LW_OUTSIDE;
+
+ for (uint32_t j = 0; j < comp->ngeoms; j++)
{
- LWGEOM *lwgeom = comp->geoms[i];
- if ( lwgeom->type == LINETYPE )
+ int on_boundary = LW_FALSE;
+ const LWGEOM *sub = comp->geoms[j];
+ if (sub->type == LINETYPE)
{
- lwline = lwgeom_as_lwline(lwgeom);
- if ( comp->ngeoms == 1 )
- {
- return ptarray_contains_point(lwline->points, pt);
- }
- else
- {
- /* Don't check closure while doing p-i-p test */
- result = ptarray_contains_point_partial(lwline->points, pt, LW_FALSE, &winding_number);
- }
+ LWLINE *lwline = lwgeom_as_lwline(sub);
+ intersections += ptarray_raycast_intersections(lwline->points, pt, &on_boundary);
+ }
+ else if (sub->type == CIRCSTRINGTYPE)
+ {
+ LWCIRCSTRING *lwcirc = lwgeom_as_lwcircstring(sub);
+ intersections += ptarrayarc_raycast_intersections(lwcirc->points, pt, &on_boundary);
}
else
{
- lwcirc = lwgeom_as_lwcircstring(lwgeom);
- if ( ! lwcirc ) {
- lwerror("Unexpected component of type %s in compound curve", lwtype_name(lwgeom->type));
- return 0;
- }
- if ( comp->ngeoms == 1 )
- {
- return ptarrayarc_contains_point(lwcirc->points, pt);
- }
- else
- {
- /* Don't check closure while doing p-i-p test */
- result = ptarrayarc_contains_point_partial(lwcirc->points, pt, LW_FALSE, &winding_number);
- }
+ lwerror("%s: unsupported type %s", __func__, lwtype_name(sub->type));
}
-
- /* Propagate boundary condition */
- if ( result == LW_BOUNDARY )
+ if (on_boundary)
return LW_BOUNDARY;
-
- wn += winding_number;
}
- /* Outside */
- if (wn == 0)
- return LW_OUTSIDE;
-
- /* Inside */
- return LW_INSIDE;
+ /*
+ * Odd number of intersections means inside.
+ * Even means outside.
+ */
+ return (intersections % 2) ? LW_INSIDE : LW_OUTSIDE;
}
LWCOMPOUND *
diff --git a/liblwgeom/measures.c b/liblwgeom/measures.c
index bc240eb45..62e2a3586 100644
--- a/liblwgeom/measures.c
+++ b/liblwgeom/measures.c
@@ -1540,37 +1540,156 @@ lw_dist2d_pt_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const P
return LW_TRUE;
}
-/* Auxiliary function to calculate the distance between 2 concentric arcs*/
-int lw_dist2d_arc_arc_concentric(const POINT2D *A1,
- const POINT2D *A2,
- const POINT2D *A3,
- double radius_A,
- const POINT2D *B1,
- const POINT2D *B2,
- const POINT2D *B3,
- double radius_B,
- const POINT2D *CENTER,
- DISTPTS *dl);
+
+
+/**
+ * @brief Calculates the intersection points of a circle and line.
+ *
+ * This function assumes the center and test point are distinct.
+ * Finds the line between the circle center and test point, and
+ * the two intersection points on the circle defined by that line.
+ * If those points fall within the provided arc (A1,A2,A3) then
+ * the points are added to the provided array (I) and the array
+ * length counter is incremented.
+ *
+ * @param A1 Point of arc
+ * @param A2 Point of arc
+ * @param A3 Point of arc
+ * @param center_A Center of arc A circle
+ * @param radius_A Radius of arc A circle
+ * @param P Point to use in calculating intersection line
+ * @param I [out] Pointer to an return value array
+ * @param ni [out] Pointer to return array size counter
+ * @return int
+ */
+static int
+lw_dist2d_circle_intersections(
+ const POINT2D *A1, const POINT2D *A2, const POINT2D *A3,
+ const POINT2D *center_A,
+ double radius_A,
+ const POINT2D *P,
+ POINT2D *I,
+ uint32_t *ni)
+{
+ POINT2D R;
+
+ // If the test point is on the center of the other
+ // arc, some other point has to be closer, by definition.
+ if (p2d_same(center_A, P))
+ return 0;
+
+ // Calculate vector from the center to the pt
+ double dir_x = center_A->x - P->x;
+ double dir_y = center_A->y - P->y;
+ double dist = sqrt(dir_x * dir_x + dir_y * dir_y);
+
+ // Normalize the direction vector to get a unit vector (length = 1)
+ double unit_x = dir_x / dist;
+ double unit_y = dir_y / dist;
+
+ // Calculate the two intersection points on the circle
+ // Point 1: Move from center_A along the unit vector by distance radius_A
+ R.x = center_A->x + unit_x * radius_A;
+ R.y = center_A->y + unit_y * radius_A;
+ if (lw_pt_in_arc(&R, A1, A2, A3))
+ I[(*ni)++] = R;
+
+ // Point 2: Move from center_A against the unit vector by distance radius_A
+ R.x = center_A->x - unit_x * radius_A;
+ R.y = center_A->y - unit_y * radius_A;
+ if (lw_pt_in_arc(&R, A1, A2, A3))
+ I[(*ni)++] = R;
+
+ return 0;
+}
+
+
+/**
+ * @brief Calculates the intersection points of two overlapping circles.
+ *
+ * This function assumes the circles are known to intersect at one or two points.
+ * Specifically, the distance 'd' between their centers must satisfy:
+ * d < (rA + rB) and d > fabs(rA - rB)
+ * If these conditions are not met, the results are undefined.
+ *
+ * @param cA [in] Pointer to the center point of the first circle (A).
+ * @param rA [in] The radius of the first circle (A).
+ * @param cB [in] Pointer to the center point of the second circle (B).
+ * @param rB [in] The radius of the second circle (B).
+ * @param I [out] Pointer to an array of at least 2 POINT2D structs to store the results.
+ * I[0] will contain the first intersection point.
+ * If a second exists, it will be in I[1].
+ * @return int The number of intersection points found (1 or 2). Returns 0 if
+ * the centers are coincident or another error occurs.
+ */
+static uint32_t
+lw_dist2d_circle_circle_intersections(
+ const POINT2D *cA, double rA,
+ const POINT2D *cB, double rB,
+ POINT2D *I)
+{
+ // Vector from center A to center B
+ double dx = cB->x - cA->x;
+ double dy = cB->y - cA->y;
+
+ // Distance between the centers
+ double d = sqrt(dx * dx + dy * dy);
+
+ // 'a' is the distance from the center of circle A to the point P,
+ // which is the base of the right triangle formed by cA, P, and an intersection point.
+ double a = (rA * rA - rB * rB + d * d) / (2.0 * d);
+
+ // 'h' is the height of that right triangle.
+ double h_squared = rA * rA - a * a;
+
+ // Due to floating point errors, h_squared can be slightly negative.
+ // This happens when the circles are perfectly tangent. Clamp to 0.
+ if (h_squared < 0.0)
+ h_squared = 0.0;
+
+ double h = sqrt(h_squared);
+
+ // Find the coordinates of point P
+ double Px = cA->x + a * (dx / d);
+ double Py = cA->y + a * (dy / d);
+
+ // The two intersection points are found by moving from P by a distance 'h'
+ // in directions perpendicular to the line connecting the centers.
+ // The perpendicular vector to (dx, dy) is (-dy, dx).
+
+ // Intersection point 1
+ I[0].x = Px - h * (dy / d);
+ I[0].y = Py + h * (dx / d);
+
+ // If h is very close to 0, the circles are tangent and there's only one intersection point.
+ if (FP_IS_ZERO(h))
+ return 1;
+
+ // Intersection point 2
+ I[1].x = Px + h * (dy / d);
+ I[1].y = Py - h * (dx / d);
+
+ return 2;
+}
+
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)
+lw_dist2d_arc_arc(
+ const POINT2D *A1, const POINT2D *A2, const POINT2D *A3,
+ const POINT2D *B1, const POINT2D *B2, const POINT2D *B3,
+ DISTPTS *dl)
{
- POINT2D CA, CB; /* Center points of arcs A and B */
+ POINT2D pts_A[16]; /* Points on A that might be the nearest */
+ POINT2D pts_B[16]; /* Points on B that might be the nearest */
+ uint32_t ai = 0, bi = 0; /* Number of points in pts_A and pts_B */
+ POINT2D center_A, center_B; /* Center points of arcs A and B */
double radius_A, radius_B, d; /* Radii of arcs A and B */
- POINT2D D; /* Mid-point between the centers CA and CB */
- int pt_in_arc_A, pt_in_arc_B; /* Test whether potential intersection point is within the arc */
+ int is_disjoint, is_overlapping, is_contained, is_same_center;
+ POINT2D intersectionPts[2];
if (dl->mode != DIST_MIN)
lwerror("lw_dist2d_arc_arc only supports mindistance");
- /* TODO: Handle case where arc is closed circle (A1 = A3) */
-
/* What if one or both of our "arcs" is actually a point? */
if (lw_arc_is_pt(B1, B2, B3) && lw_arc_is_pt(A1, A2, A3))
return lw_dist2d_pt_pt(B1, A1, dl);
@@ -1580,8 +1699,8 @@ lw_dist2d_arc_arc(const POINT2D *A1,
return lw_dist2d_pt_arc(A1, B1, B2, B3, dl);
/* Calculate centers and radii of circles. */
- radius_A = lw_arc_center(A1, A2, A3, &CA);
- radius_B = lw_arc_center(B1, B2, B3, &CB);
+ radius_A = lw_arc_center(A1, A2, A3, ¢er_A);
+ radius_B = lw_arc_center(B1, B2, B3, ¢er_B);
/* Two co-linear arcs?!? That's two segments. */
if (radius_A < 0 && radius_B < 0)
@@ -1595,275 +1714,112 @@ lw_dist2d_arc_arc(const POINT2D *A1,
if (radius_B < 0)
return lw_dist2d_seg_arc(B1, B3, A1, A2, A3, dl);
- /* Center-center distance */
- d = distance2d_pt_pt(&CA, &CB);
+ /* Circle relationships */
+ d = distance2d_pt_pt(¢er_A, ¢er_B);
+ is_disjoint = (d > (radius_A + radius_B));
+ is_contained = (d < fabs(radius_A - radius_B));
+ is_same_center = p2d_same(¢er_A, ¢er_B);
+ is_overlapping = ! (is_disjoint || is_contained || is_same_center);
- /* Concentric arcs */
- if (FP_EQUALS(d, 0.0))
- return lw_dist2d_arc_arc_concentric(A1, A2, A3, radius_A, B1, B2, B3, radius_B, &CA, dl);
+ /*
+ * Prime the array of potential closest points with the
+ * arc end points, which frequently participate in closest
+ * points.
+ */
+ pts_A[ai++] = *A1;
+ pts_A[ai++] = *A3;
+ pts_B[bi++] = *B1;
+ pts_B[bi++] = *B3;
- /* Make sure that arc "A" has the bigger radius */
- if (radius_B > radius_A)
+ /*
+ * Overlapping circles might have a zero distance
+ * case if the circle intersection points are inside both
+ * arcs.
+ */
+ if (is_overlapping)
{
- const POINT2D *tmp;
- POINT2D TP; /* Temporary point P */
- double td;
- tmp = B1;
- B1 = A1;
- A1 = tmp;
- tmp = B2;
- B2 = A2;
- A2 = tmp;
- tmp = B3;
- B3 = A3;
- A3 = tmp;
- TP = CB;
- CB = CA;
- CA = TP;
- td = radius_B;
- radius_B = radius_A;
- radius_A = td;
- }
-
- /* Circles touch at a point. Is that point within the arcs? */
- if (d == (radius_A + radius_B))
- {
- D.x = CA.x + (CB.x - CA.x) * radius_A / d;
- D.y = CA.y + (CB.y - CA.y) * radius_A / d;
-
- pt_in_arc_A = lw_pt_in_arc(&D, A1, A2, A3);
- pt_in_arc_B = lw_pt_in_arc(&D, B1, B2, B3);
-
- /* Arcs do touch at D, return it */
- if (pt_in_arc_A && pt_in_arc_B)
+ /*
+ * Find the two points the circles intersect at.
+ */
+ uint32_t npoints = lw_dist2d_circle_circle_intersections(
+ ¢er_A, radius_A,
+ ¢er_B, radius_B,
+ intersectionPts);
+ for (uint32_t i = 0; i < npoints; i++)
{
- lw_dist2d_distpts_set(dl, 0.0, &D, &D);
- return LW_TRUE;
+ /*
+ * If an intersection point is contained in both
+ * arcs, that is a location of zero distance, so
+ * we are done calculating.
+ */
+ if (lw_pt_in_arc(&intersectionPts[i], A1, A2, A3) &&
+ lw_pt_in_arc(&intersectionPts[i], B1, B2, B3))
+ {
+ lw_dist2d_distpts_set(dl, 0.0, &intersectionPts[i], &intersectionPts[i]);
+ return LW_TRUE;
+ }
}
}
- /* Disjoint or contained circles don't intersect. Closest point may be on */
- /* the line joining CA to CB. */
- else if (d > (radius_A + radius_B) /* Disjoint */ || d < (radius_A - radius_B) /* Contained */)
+
+ /*
+ * Join the circle centers and find the places that
+ * line intersects the circles. Where those places
+ * are in the arcs, they are potential sites of the
+ * closest points.
+ */
+ if (is_disjoint || is_contained || is_overlapping)
{
- POINT2D XA, XB; /* Points where the line from CA to CB cross their circle bounds */
-
- /* Calculate hypothetical nearest points, the places on the */
- /* two circles where the center-center line crosses. If both */
- /* arcs contain their hypothetical points, that's the crossing distance */
- XA.x = CA.x + (CB.x - CA.x) * radius_A / d;
- XA.y = CA.y + (CB.y - CA.y) * radius_A / d;
- XB.x = CB.x + (CA.x - CB.x) * radius_B / d;
- XB.y = CB.y + (CA.y - CB.y) * radius_B / d;
-
- pt_in_arc_A = lw_pt_in_arc(&XA, A1, A2, A3);
- pt_in_arc_B = lw_pt_in_arc(&XB, B1, B2, B3);
-
- /* If the nearest points are both within the arcs, that's our answer */
- /* the shortest distance is at the nearest points */
- if (pt_in_arc_A && pt_in_arc_B)
+ if (!is_same_center)
{
- return lw_dist2d_pt_pt(&XA, &XB, dl);
- }
- }
- /* Circles cross at two points, are either of those points in both arcs? */
- /* http://paulbourke.net/geometry/2circle/ */
- else if (d < (radius_A + radius_B))
- {
- POINT2D E, F; /* Points where circle(A) and circle(B) cross */
- /* Distance from CA to D */
- double a = (radius_A * radius_A - radius_B * radius_B + d * d) / (2 * d);
- /* Distance from D to E or F */
- double h = sqrt(radius_A * radius_A - a * a);
+ /* Add points on A that intersect line from center_A to center_B */
+ lw_dist2d_circle_intersections(
+ A1, A2, A3,
+ ¢er_A, radius_A, ¢er_B,
+ pts_A, &ai);
- /* Location of D */
- D.x = CA.x + (CB.x - CA.x) * a / d;
- D.y = CA.y + (CB.y - CA.y) * a / d;
-
- /* Start from D and project h units perpendicular to CA-D to get E */
- E.x = D.x + (D.y - CA.y) * h / a;
- E.y = D.y + (D.x - CA.x) * h / a;
-
- /* Crossing point E contained in arcs? */
- pt_in_arc_A = lw_pt_in_arc(&E, A1, A2, A3);
- pt_in_arc_B = lw_pt_in_arc(&E, B1, B2, B3);
-
- if (pt_in_arc_A && pt_in_arc_B)
- {
- lw_dist2d_distpts_set(dl, 0.0, &E, &E);
- return LW_TRUE;
+ /* Add points on B that intersect line from center_B to center_A */
+ lw_dist2d_circle_intersections(
+ B1, B2, B3,
+ ¢er_B, radius_B, ¢er_A,
+ pts_B, &bi);
}
- /* Start from D and project h units perpendicular to CA-D to get F */
- F.x = D.x - (D.y - CA.y) * h / a;
- F.y = D.y - (D.x - CA.x) * h / a;
+ /* Add points on A that intersect line to B1 */
+ lw_dist2d_circle_intersections(
+ A1, A2, A3,
+ ¢er_A, radius_A, B1,
+ pts_A, &ai);
- /* Crossing point F contained in arcs? */
- pt_in_arc_A = lw_pt_in_arc(&F, A1, A2, A3);
- pt_in_arc_B = lw_pt_in_arc(&F, B1, B2, B3);
+ /* Add points on A that intersect line to B3 */
+ lw_dist2d_circle_intersections(
+ A1, A2, A3,
+ ¢er_A, radius_A, B3,
+ pts_A, &ai);
- if (pt_in_arc_A && pt_in_arc_B)
- {
- lw_dist2d_distpts_set(dl, 0.0, &F, &F);
- return LW_TRUE;
- }
- }
- else
- {
- lwerror("lw_dist2d_arc_arc: arcs neither touch, intersect nor are disjoint! INCONCEIVABLE!");
- return LW_FALSE;
+ /* Add points on B that intersect line to A1 */
+ lw_dist2d_circle_intersections(
+ B1, B2, B3,
+ ¢er_B, radius_B, A1,
+ pts_B, &bi);
+
+ /* Add points on B that intersect line to A3 */
+ lw_dist2d_circle_intersections(
+ B1, B2, B3,
+ ¢er_B, radius_B, A3,
+ pts_B, &bi);
}
- /* Closest point is in the arc A, but not in the arc B, so */
- /* one of the B end points must be the closest. */
- if (pt_in_arc_A && !pt_in_arc_B)
- {
- lw_dist2d_pt_arc(B1, A1, A2, A3, dl);
- lw_dist2d_pt_arc(B3, A1, A2, A3, dl);
- return LW_TRUE;
- }
- /* Closest point is in the arc B, but not in the arc A, so */
- /* one of the A end points must be the closest. */
- else if (pt_in_arc_B && !pt_in_arc_A)
- {
- lw_dist2d_pt_arc(A1, B1, B2, B3, dl);
- lw_dist2d_pt_arc(A3, B1, B2, B3, dl);
- return LW_TRUE;
- }
- /* Finally, one of the end-point to end-point combos is the closest. */
- else
- {
- lw_dist2d_pt_pt(A1, B1, dl);
- lw_dist2d_pt_pt(A1, B3, dl);
- lw_dist2d_pt_pt(A3, B1, dl);
- lw_dist2d_pt_pt(A3, B3, dl);
- return LW_TRUE;
- }
+ /*
+ * Now just brute force check all pairs of participating
+ * points, to find the pair that is closest together.
+ */
+ for (uint32_t i = 0; i < ai; i++)
+ for (uint32_t j = 0; j < bi; j++)
+ lw_dist2d_pt_pt(&pts_A[i], &pts_B[j], dl);
return LW_TRUE;
}
-int
-lw_dist2d_arc_arc_concentric(const POINT2D *A1,
- const POINT2D *A2,
- const POINT2D *A3,
- double radius_A,
- const POINT2D *B1,
- const POINT2D *B2,
- const POINT2D *B3,
- double radius_B,
- const POINT2D *CENTER,
- DISTPTS *dl)
-{
- int seg_size;
- double dist_sqr, shortest_sqr;
- const POINT2D *P1;
- const POINT2D *P2;
- POINT2D proj;
-
- if (radius_A == radius_B)
- {
- /* Check if B1 or B3 are in the same side as A2 in the A1-A3 arc */
- seg_size = lw_segment_side(A1, A3, A2);
- if (seg_size == lw_segment_side(A1, A3, B1))
- {
- lw_dist2d_distpts_set(dl, 0.0, B1, B1);
- return LW_TRUE;
- }
- if (seg_size == lw_segment_side(A1, A3, B3))
- {
- lw_dist2d_distpts_set(dl, 0.0, B3, B3);
- return LW_TRUE;
- }
- /* Check if A1 or A3 are in the same side as B2 in the B1-B3 arc */
- seg_size = lw_segment_side(B1, B3, B2);
- if (seg_size == lw_segment_side(B1, B3, A1))
- {
- lw_dist2d_distpts_set(dl, 0.0, A1, A1);
- return LW_TRUE;
- }
- if (seg_size == lw_segment_side(B1, B3, A3))
- {
- lw_dist2d_distpts_set(dl, 0.0, A3, A3);
- return LW_TRUE;
- }
- }
- else
- {
- /* Check if any projection of B ends are in A*/
- seg_size = lw_segment_side(A1, A3, A2);
-
- /* B1 */
- proj.x = CENTER->x + (B1->x - CENTER->x) * radius_A / radius_B;
- proj.y = CENTER->y + (B1->y - CENTER->y) * radius_A / radius_B;
-
- if (seg_size == lw_segment_side(A1, A3, &proj))
- {
- lw_dist2d_distpts_set(dl, fabs(radius_A - radius_B), &proj, B1);
- return LW_TRUE;
- }
- /* B3 */
- proj.x = CENTER->x + (B3->x - CENTER->x) * radius_A / radius_B;
- proj.y = CENTER->y + (B3->y - CENTER->y) * radius_A / radius_B;
- if (seg_size == lw_segment_side(A1, A3, &proj))
- {
- lw_dist2d_distpts_set(dl, fabs(radius_A - radius_B), &proj, B3);
- return LW_TRUE;
- }
-
- /* Now check projections of A in B */
- seg_size = lw_segment_side(B1, B3, B2);
-
- /* A1 */
- proj.x = CENTER->x + (A1->x - CENTER->x) * radius_B / radius_A;
- proj.y = CENTER->y + (A1->y - CENTER->y) * radius_B / radius_A;
- if (seg_size == lw_segment_side(B1, B3, &proj))
- {
- lw_dist2d_distpts_set(dl, fabs(radius_A - radius_B), &proj, A1);
- return LW_TRUE;
- }
-
- /* A3 */
- proj.x = CENTER->x + (A3->x - CENTER->x) * radius_B / radius_A;
- proj.y = CENTER->y + (A3->y - CENTER->y) * radius_B / radius_A;
- if (seg_size == lw_segment_side(B1, B3, &proj))
- {
- lw_dist2d_distpts_set(dl, fabs(radius_A - radius_B), &proj, A3);
- return LW_TRUE;
- }
- }
-
- /* Check the shortest between the distances of the 4 ends */
- shortest_sqr = dist_sqr = distance2d_sqr_pt_pt(A1, B1);
- P1 = A1;
- P2 = B1;
-
- dist_sqr = distance2d_sqr_pt_pt(A1, B3);
- if (dist_sqr < shortest_sqr)
- {
- shortest_sqr = dist_sqr;
- P1 = A1;
- P2 = B3;
- }
-
- dist_sqr = distance2d_sqr_pt_pt(A3, B1);
- if (dist_sqr < shortest_sqr)
- {
- shortest_sqr = dist_sqr;
- P1 = A3;
- P2 = B1;
- }
-
- dist_sqr = distance2d_sqr_pt_pt(A3, B3);
- if (dist_sqr < shortest_sqr)
- {
- shortest_sqr = dist_sqr;
- P1 = A3;
- P2 = B3;
- }
-
- lw_dist2d_distpts_set(dl, sqrt(shortest_sqr), P1, P2);
- return LW_TRUE;
-}
/**
Finds the shortest distance between two segments.
diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index 82a3835e4..2bc6126b5 100644
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -742,17 +742,285 @@ ptarray_is_closed_z(const POINTARRAY *in)
}
/**
-* Return LW_INSIDE if the point is inside the POINTARRAY,
-* LW_OUTSIDE if it is outside, and LW_BOUNDARY if it is on
-* the boundary.
-* LW_INSIDE == 1, LW_BOUNDARY == 0, LW_OUTSIDE == -1
-*/
+ * The following is based on the "Fast Winding Number Inclusion of a Point
+ * in a Polygon" algorithm by Dan Sunday.
+ * http://softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm#Winding%20Number
+ *
+ * Return:
+ * - LW_INSIDE (1) if the point is inside the POINTARRAY
+ * - LW_BOUNDARY (0) if it is on the boundary.
+ * - LW_OUTSIDE (-1) if it is outside
+ */
int
ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt)
{
- return ptarray_contains_point_partial(pa, pt, LW_TRUE, NULL);
+ const POINT2D *seg1, *seg2;
+ int wn = 0;
+
+ 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");
+
+ for (uint32_t i = 1; i < pa->npoints; i++)
+ {
+ double side, ymin, ymax;
+
+ seg2 = getPoint2d_cp(pa, i);
+
+ /* Zero length segments are ignored. */
+ if (p2d_same(seg1, seg2))
+ {
+ seg1 = seg2;
+ continue;
+ }
+
+ ymin = FP_MIN(seg1->y, seg2->y);
+ ymax = FP_MAX(seg1->y, seg2->y);
+
+ /* Only test segments in our vertical range */
+ if (pt->y > ymax || pt->y < ymin)
+ {
+ seg1 = seg2;
+ 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 LW_BOUNDARY;
+ }
+
+ /*
+ * 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 increment.
+ */
+ 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;
+ }
+
+ /* wn == 0 => Outside, wn != 0 => Inside */
+ return wn == 0 ? LW_OUTSIDE : LW_INSIDE;
}
+int
+ptarray_raycast_intersections(const POINTARRAY *pa, const POINT2D *p, int *on_boundary)
+{
+ // A valid linestring must have at least 2 point
+ if (pa->npoints < 2)
+ lwerror("%s called on invalid linestring", __func__);
+
+ int intersections = 0;
+ double px = p->x;
+ double py = p->y;
+
+ // Iterate through each edge of the polygon
+ for (uint32_t i = 0; i < pa->npoints-1; ++i)
+ {
+ const POINT2D* p1 = getPoint2d_cp(pa, i);
+ const POINT2D* p2 = getPoint2d_cp(pa, i+1);
+
+ /* Skip zero-length edges */
+ if (p2d_same(p1, p2))
+ continue;
+
+ /* --- Step 1: Check if the point is ON the boundary edge --- */
+ if (lw_pt_on_segment(p1, p2, p))
+ {
+ *on_boundary = LW_TRUE;
+ return 0;
+ }
+
+ /* --- Step 2: Perform the Ray Casting intersection test --- */
+
+ /*
+ * Check if the horizontal ray from p intersects the edge (p1, p2).
+ * This is the core condition for handling vertices correctly:
+ * - One vertex must be strictly above the ray (py < vertex.y)
+ * - The other must be on or below the ray (py >= vertex.y)
+ */
+ if (((p1->y <= py) && (py < p2->y)) || ((p2->y <= py) && (py < p1->y)))
+ {
+ /*
+ * Calculate the x-coordinate where the edge intersects the ray's horizontal line.
+ * Formula: x = x1 + (y - y1) * (x2 - x1) / (y2 - y1)
+ */
+ double x_intersection = p1->x + (py - p1->y) * (p2->x - p1->x) / (p2->y - p1->y);
+
+ /*
+ * If the intersection point is to the right of our test point,
+ * it's a valid "crossing".
+ */
+ if (x_intersection > px)
+ {
+ intersections++;
+ }
+ }
+ }
+
+ return intersections;
+}
+
+
+/**
+ * @brief Calculates the intersection points of a circle and a horizontal line.
+ *
+ * The equation of a circle is (x - cx)^2 + (y - cy)^2 = r^2.
+ * The equation of the horizontal line is y = y_line.
+ * Substituting y_line into the circle equation gives:
+ * (x - cx)^2 = r^2 - (y_line - cy)^2
+ * This function solves for x.
+ *
+ * @param center A pointer to the center point of the circle.
+ * @param radius The radius of the circle.
+ * @param ray The y-coordinate of the horizontal line.
+ * @param i0 A pointer to a POINT2D to store the first intersection point.
+ * @param i1 A pointer to a POINT2D to store the second intersection point.
+ *
+ * @return The number of intersection points found (0, 1, or 2).
+ */
+static int
+circle_raycast_intersections(const POINT2D *center, double radius, double ray, POINT2D *i0, POINT2D *i1)
+{
+ // Calculate the vertical distance from the circle's center to the horizontal line.
+ double dy = ray - center->y;
+
+ // If the absolute vertical distance is greater than the radius, there are no intersections.
+ if (fabs(dy) > radius)
+ return 0;
+
+ // Use the Pythagorean theorem to find the horizontal distance (dx) from the
+ // center's x-coordinate to the intersection points.
+ // dx^2 + dy^2 = radius^2 => dx^2 = radius^2 - dy^2
+ double dx_squared = radius * radius - dy * dy;
+
+ // Case 1: One intersection (tangent)
+ // This occurs when the line just touches the top or bottom of the circle.
+ // dx_squared will be zero. We check against a small epsilon for floating-point safety.
+ if (FP_EQUALS(dx_squared, 0.0))
+ {
+ i0->x = center->x;
+ i0->y = ray;
+ return 1;
+ }
+
+ // Case 2: Two intersections
+ // The line cuts through the circle.
+ double dx = sqrt(dx_squared);
+
+ // The first intersection point has the smaller x-value.
+ i0->x = center->x - dx;
+ i0->y = ray;
+
+ // The second intersection point has the larger x-value.
+ i1->x = center->x + dx;
+ i1->y = ray;
+
+ return 2;
+}
+
+
+int
+ptarrayarc_raycast_intersections(const POINTARRAY *pa, const POINT2D *p, int *on_boundary)
+{
+ int intersections = 0;
+ double px = p->x;
+ double py = p->y;
+
+ assert(on_boundary);
+
+ // A valid circular arc must have at least 3 vertices (circle).
+ if (pa->npoints < 3)
+ lwerror("%s called on invalid circularstring", __func__);
+
+ if (pa->npoints % 2 == 0)
+ lwerror("%s called with even number of points", __func__);
+
+ // Iterate through each arc of the circularstring
+ for (uint32_t i = 1; i < pa->npoints-1; i +=2)
+ {
+ const POINT2D* p0 = getPoint2d_cp(pa, i-1);
+ const POINT2D* p1 = getPoint2d_cp(pa, i);
+ const POINT2D* p2 = getPoint2d_cp(pa, i+1);
+ POINT2D center = {0,0};
+ double radius, d;
+ GBOX gbox;
+
+ // Skip zero-length arc
+ if (lw_arc_is_pt(p0, p1, p2))
+ continue;
+
+ // --- Step 1: Check if the point is ON the boundary edge ---
+ if (p2d_same(p0, p) || p2d_same(p1, p) || p2d_same(p2, p))
+ {
+ *on_boundary = LW_TRUE;
+ return 0;
+ }
+
+ // Calculate some important pieces
+ radius = lw_arc_center(p0, p1, p2, ¢er);
+
+ d = distance2d_pt_pt(p, ¢er);
+ if (FP_EQUALS(d, radius) && lw_pt_in_arc(p, p0, p1, p2))
+ {
+ *on_boundary = LW_TRUE;
+ return 0;
+ }
+
+ // --- Step 2: Perform the Ray Casting intersection test ---
+
+ // Only process arcs that our ray crosses
+ lw_arc_calculate_gbox_cartesian_2d(p0, p1, p2, &gbox);
+ if ((gbox.ymin <= py) && (py < gbox.ymax))
+ {
+ // Point of intersection on the circle that defines the arc
+ POINT2D i0, i1;
+
+ // How many points of intersection are there
+ int iCount = circle_raycast_intersections(¢er, radius, py, &i0, &i1);
+
+ // Nothing to see here
+ if (iCount == 0)
+ continue;
+
+ // Cannot think of a case where a grazing is not a
+ // no-op
+ if (iCount == 1)
+ continue;
+
+ // So we must have 2 intersections
+ // Only increment the counter for intersections to the right
+ // of the test point
+ if (i0.x > px && lw_pt_in_arc(&i0, p0, p1, p2))
+ intersections++;
+
+ if (i1.x > px && lw_pt_in_arc(&i1, p0, p1, p2))
+ intersections++;
+ }
+ }
+
+ return intersections;
+}
/*
* The following is based on the "Fast Winding Number Inclusion of a Point
@@ -765,8 +1033,7 @@ ptarray_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int chec
int wn = 0;
uint32_t i;
double side;
- const POINT2D *seg1;
- const POINT2D *seg2;
+ const POINT2D *seg1, *seg2;
double ymin, ymax;
seg1 = getPoint2d_cp(pa, 0);
@@ -855,160 +1122,16 @@ ptarray_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int chec
int
ptarrayarc_contains_point(const POINTARRAY *pa, const POINT2D *pt)
{
- return ptarrayarc_contains_point_partial(pa, pt, LW_TRUE /* Check closed*/, NULL);
-}
+ int on_boundary = LW_FALSE;
+ int intersections;
+ if (!ptarray_is_closed_2d(pa))
+ lwerror("%s called on unclosed ring", __func__);
-int
-ptarrayarc_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number)
-{
- int wn = 0;
- uint32_t i;
- int 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) == 0 )
- {
- lwerror("ptarrayarc_contains_point called with even number of points");
- return LW_OUTSIDE;
- }
-
- /* Check for not an arc ring (always have >= 3 points) */
- if ( pa->npoints < 3 )
- {
- lwerror("ptarrayarc_contains_point called too-short pointarray");
- return LW_OUTSIDE;
- }
-
- /* Check for unclosed case */
- seg1 = getPoint2d_cp(pa, 0);
- seg3 = getPoint2d_cp(pa, pa->npoints-1);
- if ( check_closed && ! p2d_same(seg1, seg3) )
- {
- lwerror("ptarrayarc_contains_point called on unclosed ring");
- return LW_OUTSIDE;
- }
- /* OK, it's closed. Is it just one circle? */
- else if ( p2d_same(seg1, seg3) && pa->npoints == 3 )
- {
- double radius, d;
- POINT2D c;
- seg2 = getPoint2d_cp(pa, 1);
-
- /* Wait, it's just a point, so it can't contain anything */
- if ( lw_arc_is_pt(seg1, seg2, seg3) )
- return LW_OUTSIDE;
-
- /* See if the point is within the circle radius */
- radius = lw_arc_center(seg1, seg2, seg3, &c);
- d = distance2d_pt_pt(pt, &c);
- if ( FP_EQUALS(d, radius) )
- return LW_BOUNDARY; /* Boundary of circle */
- else if ( d < radius )
- return LW_INSIDE; /* Inside circle */
- else
- return LW_OUTSIDE; /* Outside circle */
- }
- else if ( p2d_same(seg1, pt) || p2d_same(seg3, pt) )
- {
- return LW_BOUNDARY; /* Boundary case */
- }
-
- /* Start on the ring */
- seg1 = getPoint2d_cp(pa, 0);
- for ( i=1; i < pa->npoints; i += 2 )
- {
- seg2 = getPoint2d_cp(pa, i);
- seg3 = getPoint2d_cp(pa, i+1);
-
- /* Catch an easy boundary case */
- if( p2d_same(seg3, pt) )
- return LW_BOUNDARY;
-
- /* Skip arcs that have no size */
- if ( lw_arc_is_pt(seg1, seg2, seg3) )
- {
- seg1 = seg3;
- continue;
- }
-
- /* Only test segments in our vertical range */
- lw_arc_calculate_gbox_cartesian_2d(seg1, seg2, seg3, &gbox);
- if ( pt->y > gbox.ymax || pt->y < gbox.ymin )
- {
- seg1 = seg3;
- continue;
- }
-
- /* Outside of horizontal range, and not between end points we also skip */
- if ( (pt->x > gbox.xmax || pt->x < gbox.xmin) &&
- (pt->y > FP_MAX(seg1->y, seg3->y) || pt->y < FP_MIN(seg1->y, seg3->y)) )
- {
- seg1 = seg3;
- continue;
- }
-
- side = lw_arc_side(seg1, seg2, seg3, pt);
-
- /* On the boundary */
- if ( (side == 0) && lw_pt_in_arc(pt, seg1, seg2, seg3) )
- {
- return LW_BOUNDARY;
- }
-
- /* Going "up"! Point to left of arc. */
- if ( side < 0 && (seg1->y <= pt->y) && (pt->y < seg3->y) )
- {
- wn++;
- }
-
- /* Going "down"! */
- if ( side > 0 && (seg3->y <= pt->y) && (pt->y < seg1->y) )
- {
- wn--;
- }
-
- /* Inside the arc! */
- if ( pt->x <= gbox.xmax && pt->x >= gbox.xmin )
- {
- POINT2D C;
- double radius = lw_arc_center(seg1, seg2, seg3, &C);
- double d = distance2d_pt_pt(pt, &C);
-
- /* On the boundary! */
- if ( d == radius )
- return LW_BOUNDARY;
-
- /* Within the arc! */
- if ( d < radius )
- {
- /* Left side, increment winding number */
- if ( side < 0 )
- wn++;
- /* Right side, decrement winding number */
- if ( side > 0 )
- wn--;
- }
- }
-
- seg1 = seg3;
- }
-
- /* Sent out the winding number for calls that are building on this as a primitive */
- if ( winding_number )
- *winding_number = wn;
-
- /* Outside */
- if (wn == 0)
- {
- return LW_OUTSIDE;
- }
-
- /* Inside */
- return LW_INSIDE;
+ intersections = ptarrayarc_raycast_intersections(pa, pt, &on_boundary);
+ if (on_boundary)
+ return LW_BOUNDARY;
+ else
+ return (intersections % 2) ? LW_INSIDE : LW_OUTSIDE;
}
/**
-----------------------------------------------------------------------
Summary of changes:
NEWS | 1 +
liblwgeom/cunit/cu_measures.c | 33 +++
liblwgeom/cunit/cu_ptarray.c | 42 ++--
liblwgeom/liblwgeom_internal.h | 4 +-
liblwgeom/lwalgorithm.c | 26 ++-
liblwgeom/lwcompound.c | 73 +++---
liblwgeom/measures.c | 508 +++++++++++++++++++----------------------
liblwgeom/ptarray.c | 415 +++++++++++++++++++++------------
8 files changed, 613 insertions(+), 489 deletions(-)
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list