Hi Phil,

Thanks for analysing this. In the past I have seen awkward fills as well. The 
code is rather tricky and it seems you have uncovered yet another case where 
things may go wrong. The test suite should be able to capture if anything 
breaks with your change. So, if that is not the case, then go ahead merging it, 
I’d say.

Regards,

Arjen

From: Phil Rosenberg <p.d.rosenb...@gmail.com>
Sent: donderdag 11 mei 2023 22:42
To: PLplot development list <plplot-devel@lists.sourceforge.net>
Subject: Re: [Plplot-devel] bug in notpointinpolygon

Caution: This message was sent from outside of Deltares. Please do not click 
links or open attachments unless you recognize the source of this email and 
know the content is safe. Please report all suspicious emails to 
"servicedesk-...@deltares.nl<mailto:servicedesk-...@deltares.nl>" as an 
attachment.

I have managed to further isolate the issue and I've attached a minimum 
example. Internal integer pixel units are a factor, so it might not show up on 
all devices. This example shows up the error on my Windows machine with the 
wxWidgets device.

The example draws a single rectangle using plfill, but the rectangle itself is 
well outside the range of the axes. Despite the position of the rectangle, the 
result is that the whole plot area is filled.

The bug occurs when drawing polygons with a very high aspect ratio (tall and 
thin), outside the plot area in certain positions. It may even only happen for 
polygons one internal plot unit wide.

As part of the plP_plfclp code, we check if each of the bottom left corner of 
the plot area is within the polygon and if no edges of the polygon intersect 
the plot area. If both these conditions are true, the polygon must cover the 
whole plot area, so we fill it all. To check if the bottom left corner is 
within the polygon, we trace a ray from that point in a particular direction. 
If the ray crosses an even number of polygon edges the point is external and if 
it crosses an odd number of edges then the point is internal. The crossings are 
checked by a function called notcrossed and they are accumulated over a whole 
polygon by a function called notpointinpolygon.

However, there is some ambiguity due to floating point arithmetic when lines 
are close to parallel and intersections are close to vertices. Hence, 
notcrossed can return 0 (crossed) 1 (notcrossed) and a series of other values 
indicating that we are uncertain if we have crossed or not. One of those return 
values is PL_NEAR_PARALLEL(32), however, notpointinpolygon does not check for 
this and counts it as true, i.e. not crossed. In almost all scenarios, the 
return of PL_NEAR_PARALLEL for one edge, results in a PL_NEAR_POINT (where 
POINT is A1, A2, B1 or B2) return for another edge. This return value is picked 
up by notpointinpolygon and is reflected in its return value, so this bug 
almost never manifests.

However, it turns out that a polygon edge with size of 1 internal unit, if 
intersected, will always generate a PL_NEAR_PARALLEL return value, even if the 
lines are close to perpendicular. this is treated as not crossed. If the 
polygon has a large aspect ratio, the ray can intersect another edge 
significantly far from its ends that it does not trigger a PL_NEAR_POINT return 
value. notpointinpolygon therefore misses one edge intersection and incorrectly 
thinks the point is in the polygon.

Where this happens for the bottom left corner, and no edges cross the plot area 
we end up filling the whole box. Although it might seem like this is very 
unlikely, if you use plshades to do a contour plot  with very fine horizontal 
resolution and much coarser vertical resolution and zoom in on the y axis to 
crop some of the data out, then it turns out to be almost inevitable that at 
least one polygon will trigger the problem. In my case I'm plotting data from a 
meteorological lidar - high temporal resolution, coarser vertical resolution.

To fix this I've written an optimised version of notcrossed, called 
notcrossedhorizontalray. This uses a horizontal ray rather than a ray to a 
point near the polygon. It uses entirely boolean logic and integer maths except 
for one division, so it should be accurate. I updated notpointinpolygon to use 
the new function and recognise the unsure case. I added some extra checks in 
plP_plfclp. I attach a patch if anyone wants to test it. Clearly I'm not the 
first person to wrestle with this code as there is another commented out 
version. In my code USE_FILL_INTERSECTION_POLYGON is not defined, so there is a 
large chunk of code unused. I don't know if this gets used in some builds, but 
maybe it could be cleaned up if not.

If nobody objects I will commit the change to the repo

Thanks
Phil

On Wed, 10 May 2023 at 10:38, Phil Rosenberg 
<p.d.rosenb...@gmail.com<mailto:p.d.rosenb...@gmail.com>> wrote:
Hi again
I just found a bug in plP_plfclp.

I hit a scenario where during a plshades call, my whole window (including 
outside the x and y limits of the data) got filled with one of the colours. 
I've traced this to plP_plfclp, the drawable variable and notpointinpolygon.

What happens is that, while filling one polygon the bottom left corner gets 
incorrectly assigned as inside the polygon and drawable gets incorrectly 
assigned as false. This results in plP_plfclp thinking that the polygon 
encircles the whole plot window and the plot window gets filled with the colour.

I haven't had chance to check why drawable ends up as false, however, the error 
in notpointinpolygon comes from notcrossed returning PL_NEAR_PARALLEL. 
notpointinparallel doesn't deal with this special return value and instead 
treats it as returning true. In this particular case, it should be returning 
false and so everything goes wrong.

I suggest a couple of fixes.

in notpointinpolygon, we should check if the point is outside the polygon 
bounding box. This will catch many of the near parallel cases.

For the reminder we then need a better way to find a point to do the ray 
tracing or to catch the special codes and then choose another point. This is 
easier to do when we know the point is close to the polygon, because then 
moving our ray tracing point has a bigger effect. We can also afford to do some 
more complicated maths here because almost all cases will be picked up by the 
bounding box case.

I'll work on this a bit today and send round a patch for people to review if 
they wish before committing.

Thanks

Phil
DISCLAIMER: This message is intended exclusively for the addressee(s) and may 
contain confidential and privileged information. If you are not the intended 
recipient please notify the sender immediately and destroy this message. 
Unauthorized use, disclosure or copying of this message is strictly prohibited. 
The foundation 'Stichting Deltares', which has its seat at Delft, The 
Netherlands, Commercial Registration Number 41146461, is not liable in any way 
whatsoever for consequences and/or damages resulting from the improper, 
incomplete and untimely dispatch, receipt and/or content of this e-mail.
_______________________________________________
Plplot-devel mailing list
Plplot-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/plplot-devel

Reply via email to