[SCM] PostGIS branch master updated. 3.6.0rc2-75-gdc860997c
git at osgeo.org
git at osgeo.org
Wed Oct 1 15:49:57 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, master has been updated
via dc860997c4d074ad8ca11e5f1c53aec7254d8ff0 (commit)
via 791660030a5071bf6c56d43066728200de50be5e (commit)
via 332bebf6c3eb13a5a0c1c3e3678d277919a20a7b (commit)
via 063388ceda8d264cd1c1ed468c722ebbaf76c031 (commit)
via 6e8b86bd9204e494e330b296156a0ddddf36cd50 (commit)
from 45af77d62d20b37bc4ed47c75896013888082941 (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 dc860997c4d074ad8ca11e5f1c53aec7254d8ff0
Merge: 791660030 45af77d62
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Wed Oct 1 15:32:36 2025 -0700
Merge branch 'master' into master-5989
commit 791660030a5071bf6c56d43066728200de50be5e
Merge: 332bebf6c 34ee8312c
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Wed Oct 1 14:17:21 2025 -0700
Merge branch 'master' into master-5989
commit 332bebf6c3eb13a5a0c1c3e3678d277919a20a7b
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Wed Oct 1 13:57:10 2025 -0700
Quiet uninit warning maybe?
diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index 818d20523..507e2283d 100644
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -962,7 +962,7 @@ ptarrayarc_raycast_intersections(const POINTARRAY *pa, const POINT2D *p, int *on
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;
+ POINT2D center = {0,0};
double radius, d;
GBOX gbox;
commit 063388ceda8d264cd1c1ed468c722ebbaf76c031
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Wed Oct 1 13:50:48 2025 -0700
Match up types in comparison
diff --git a/liblwgeom/measures.c b/liblwgeom/measures.c
index b99f0db21..62e2a3586 100644
--- a/liblwgeom/measures.c
+++ b/liblwgeom/measures.c
@@ -1622,7 +1622,7 @@ lw_dist2d_circle_intersections(
* @return int The number of intersection points found (1 or 2). Returns 0 if
* the centers are coincident or another error occurs.
*/
-static int
+static uint32_t
lw_dist2d_circle_circle_intersections(
const POINT2D *cA, double rA,
const POINT2D *cB, double rB,
@@ -1741,7 +1741,7 @@ lw_dist2d_arc_arc(
/*
* Find the two points the circles intersect at.
*/
- int npoints = lw_dist2d_circle_circle_intersections(
+ uint32_t npoints = lw_dist2d_circle_circle_intersections(
¢er_A, radius_A,
¢er_B, radius_B,
intersectionPts);
commit 6e8b86bd9204e494e330b296156a0ddddf36cd50
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Wed Oct 1 13:35:45 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/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..b99f0db21 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 int
+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.
+ */
+ int 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 a1e790dee..818d20523 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;
+ 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:
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 +++++++++++++++++++++------------
7 files changed, 612 insertions(+), 489 deletions(-)
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list