Re: [Plplot-devel] bug in notpointinpolygon

[Arjen M. <Arj...@de...>]

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...@gm...>
Sent: donderdag 11 mei 2023 22:42
To: PLplot development list <plp...@li...>
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 "Ser...@de...<mailto:Ser...@de...>" 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...@gm...<mailto:p.d...@gm...>> 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.