From 92fd99d5468ddad84a9593cd4146c9fc45d8175d Mon Sep 17 00:00:00 2001 From: Phil Rosenberg Date: Thu, 11 May 2023 21:27:04 +0100 Subject: [PATCH] Fix point in polygon --- src/plfill.c | 97 ++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 76 insertions(+), 21 deletions(-) mode change 100644 => 100755 src/plfill.c diff --git a/src/plfill.c b/src/plfill.c old mode 100644 new mode 100755 index eda932808..b286bb90a --- a/src/plfill.c +++ b/src/plfill.c @@ -58,13 +58,15 @@ // enum PL_CrossedStatus { + PL_CROSSED = 0x0, PL_NOT_CROSSED = 0x1, PL_NEAR_A1 = 0x2, PL_NEAR_A2 = 0x4, PL_NEAR_B1 = 0x8, PL_NEAR_B2 = 0x10, PL_NEAR_PARALLEL = 0x20, - PL_PARALLEL = 0x40 + PL_PARALLEL = 0x40, + PL_NEAR_LINE = 0x80 }; struct point @@ -116,6 +118,8 @@ static int notcrossed( PLINT *xintersect, PLINT *yintersect, PLINT xA1, PLINT yA1, PLINT xA2, PLINT yA2, PLINT xB1, PLINT yB1, PLINT xB2, PLINT yB2 ); +static int +notcrossedhorizontalray(PLINT x1, PLINT y1, PLINT x2, PLINT y2, PLINT xp, PLINT yp); //-------------------------------------------------------------------------- // void plfill() @@ -837,7 +841,9 @@ plP_plfclp( PLINT *x, PLINT *y, PLINT npts, // if ( iclp == 0 ) { - if ( inside_lb ) + //some of inside_xx may be 2 indicating we are unsure. Check at least + //one is definitely inside + if ( (inside_lb == 1 || inside_lu == 1 || inside_rb == 1 || inside_ru == 1) ) { xclp[0] = (short) xmin; yclp[0] = (short) ymin; xclp[1] = (short) xmax; yclp[1] = (short) ymin; @@ -1260,37 +1266,23 @@ notpointinpolygon( PLINT n, PLINT_VECTOR x, PLINT_VECTOR y, PLINT xp, PLINT yp ) #ifdef NEW_NOTPOINTINPOLYGON_CODE int i, im1, ifnotcrossed; int count_crossings = 0; - PLINT xmin, xout, yout, xintersect, yintersect; - - // Determine a point outside the polygon - - xmin = x[0]; - xout = x[0]; - yout = y[0]; - for ( i = 1; i < n; i++ ) - { - xout = MAX( xout, x[i] ); - xmin = MIN( xmin, x[i] ); - } - // + 10 to make sure completely outside. - xout = xout + ( xout - xmin ) + 10; - // Determine whether the line between (xout, yout) and (xp, yp) intersects - // one of the polygon segments. + // Determine whether the horizontal ray from (xp, yp) intersects + // each of the polygon segments. im1 = n - 1; for ( i = 0; i < n; i++ ) { if ( !( x[im1] == x[i] && y[im1] == y[i] ) ) { - ifnotcrossed = notcrossed( &xintersect, &yintersect, + ifnotcrossed = notcrossedhorizontalray( x[im1], y[im1], x[i], y[i], - xp, yp, xout, yout ); + xp, yp ); if ( !ifnotcrossed ) count_crossings++; - else if ( ifnotcrossed & ( PL_NEAR_A1 | PL_NEAR_A2 | PL_NEAR_B1 | PL_NEAR_B2 ) ) + else if ( ifnotcrossed ==2 ) return 1; } im1 = i; @@ -2019,6 +2011,69 @@ notcrossed( PLINT * xintersect, PLINT * yintersect, return status; } +//this is an optimised version of notcrossed, but where the second line is a semi-infinite +//horizontal ray in the positive direction from point xp,yp. +//returns 0 if the line crosses the ray, returns 1 if the line does not cross the ray +//returns 2 if we are not sure, due to floating point accuracy. +//If the line ends exactly on the ray, we assume the end of the line is actually just slightly +//higher than the ray. +//There is exactly one floating point division in the calculation, so we should be pretty +//accurate. +int +notcrossedhorizontalray(PLINT x1, PLINT y1, PLINT x2, PLINT y2, PLINT xp, PLINT yp) +{ + PLFLT dy1f, dy2f, xpf; + PLFLT xcrossupper, xcrosslower, dxbydy, epsilon; + if (xp > x1 && xp > x2) // the line is all to the left of the ray, return PL_NOT_CROSSED + return PL_NOT_CROSSED; + dy1f = (PLFLT)(y1 - yp); + dy2f = (PLFLT)(y2 - yp); + if (dy1f * dy2f > 0.0) //the line is above or below the ray, return PL_NOT_CROSSED + return PL_NOT_CROSSED; + if (xp < x1 && xp < x2) //the line is to the right of the ray, check yp isn't exacly on y1 or y2 + { + + //if we are exacly on y1 or y2, do further checks + if (y1 == yp) //assume y1 is slightly higher than it really is + { + if (y2 < yp) + return PL_CROSSED; + else + return PL_NOT_CROSSED; //this will include y1==y2, which should return 1 + } + if (y2 == yp) //assume y1 is slightly higher than it really is + { + if (y1 < yp) + return PL_CROSSED; + else + return PL_NOT_CROSSED; //this will include y1==y2, which should return 1 + } + return 0; //the line is well within the vertical limit of the line + } + if (y1 == y2) //horizontal lines never cross the ray - this is compatible with the points slightly higher. + return PL_NOT_CROSSED; + //the above logic will deal with almost all cases. We are left with only scenarios where the point is within the + //bounding box of the line. we just need to find the x value where the line crosses y=yp. + //In the equation below abs(dy1-dy2) must be at least 1, so it is numerically stable. + //also, we have mostly integer maths, so we are pretty darned accuarate +#ifdef PL_DOUBLE + epsilon = DBL_EPSILON; +#else + epsilon = FLT_EPSILON; +#endif + + xpf = (PLFLT)xp; + dxbydy = (PLFLT)(x1 - x2) / (dy1f - dy2f) * (-dy1f); + xcrossupper = (dxbydy * epsilon) * (-dy1f); + xcrosslower = (dxbydy / epsilon) * (-dy1f); + + if ((xcrossupper - xpf) * (xcrosslower - xpf) < 0.0) //the point is within the calcualtion accuaracy of the line + return PL_NEAR_LINE; + if (xpf < xcrosslower) + return PL_CROSSED; + return PL_NOT_CROSSED; +} + #ifdef USE_FILL_INTERSECTION_POLYGON // Decide if polygon has a positive orientation or not. // Note the simple algorithm given in -- 2.34.1