lib/vector: improve numerical stability for point-in-polygon calculations

[metzm]

What this PR does:

• fix and synchronise intersection calculations in dig_x_intersect(), Vect__intersect_x_line_with_poly() and Vect__intersect_y_line_with_poly()
• use Vect_point_in_poly() whenever possible for consistency and to reduce code maintenance burden
• avoid floating point precision errors with changed calculations

The motivation for this PR arises from two issues:

• Vect_get_point_in_poly_isl() fails for very small or very thin areas. This PR improves this fn by successfully calculating centroid's coordinates for more very small or very thin areas. However, because of fp precision limits, it is technically not possible to calculate centroid's coordinates for too small or too thin areas. In these cases it can not be guaranteed that the centroid is really inside the area.
• Even if tests in Vect_get_point_in_poly_isl() confirm that the calculated centroid's coordinates are indeed inside the given area, Vect_find_area() might in extreme cases later on assign this centroid to a different area, causing problems with both duplicate centroids and missing centroids.

I am not sure if this PR should be regarded as a bug fix or as an enhancement. Suggestions welcome!

Please excuse the many changes in this single PR, but they really belong together.

The group of PRs #7333, #7338, #7366, #7370 belong together.

[petrasovaa]

Could you add a test, e.g. based on the branch name, you have a v.overlay error this is fixing?

[metzm]

Here are test data for nc_spm:
testdata.zip

Unpack the two vectors in GRASS with

v.unpack in=ainput.pack out=ainput
v.unpack in=binput.pack out=binput

Test command is

v.overlay ainput=ainput binput=binput output=test_topoerror snap=1e-8 operator=or --v

With main, this gives the following topology building output of the intermediate vector:

[...]
Number of boundaries: 91768
Number of centroids: 0
Number of areas: 23827
Number of isles: 349
WARNING: Vect_get_point_in_poly_isl(): collapsed area
WARNING: Cannot calculate area centroid
WARNING: Vect_get_point_in_poly_isl(): collapsed area
WARNING: Cannot calculate area centroid
WARNING: Vect_get_point_in_poly_isl(): collapsed area
WARNING: Cannot calculate area centroid
WARNING: Vect_get_point_in_poly_isl(): collapsed area
WARNING: Cannot calculate area centroid
WARNING: Vect_get_point_in_poly_isl(): collapsed area
WARNING: Cannot calculate area centroid
[...]
WARNING: Number of duplicate centroids: 4

With this PR, the output is

[...]
Number of boundaries: 91768
Number of centroids: 0
Number of areas: 23827
Number of isles: 349
WARNING: Vect_get_point_in_poly_isl(): collapsed area
WARNING: Cannot calculate area centroid
[...]
WARNING: Number of duplicate centroids: 1

As I wrote in the description, for very small or very thin areas it is technically not possible to calculate a centroid that is guaranteed to be inside the area. For v.overlay a solution could be to add an option to remove small areas with Vect_remove_small_areas(), but that should go into a separate PR vor v.overlay.

[metzm]

The PR #7338 should be applied first before coming up with any tests for v.overlay regarding this PR, because PR #7338 fixes the v.overlay bug of writing out invalid centroids.

[echoix]

I don't feel confident enough to be a good reviewer here. Were the previous rounds of review enough ? If so, I could still rely on them to unlock you and have the PR approved

[metzm]

I don't feel confident enough to be a good reviewer here. Were the previous rounds of review enough ? If so, I could still rely on them to unlock you and have the PR approved

This is the conversation that still needs to be resolved #7333 (review) Maybe give @petrasovaa a few more days to resolve the conversation herself.

[metzm]

I don't feel confident enough to be a good reviewer here. Were the previous rounds of review enough ? If so, I could still rely on them to unlock you and have the PR approved

This is the conversation that still needs to be resolved #7333 (review) Maybe give @petrasovaa a few more days to resolve the conversation herself.

@echoix I have addressed the issue of Claude's review posted by @petrasovaa. The new function Vect_find_poly_centroid_cog() is now in the public header, avoiding a a namespace-pollution hazard.
However, @petrasovaa is not answering and that is the only change request blocking this PR. I need this PR to move forward. Maybe you feel confident enough to approve since the last outstanding change request is solved? ;-)

[nilason]

Code-wise this looks good to me, I'm not in a position to judge the algorithms.

Just a few last notes:

Using bool pre-C23 requires #include <stdbool.h>, this is probably done in one of the other includes, so using this in .c file not .h file, shouldn't necessarily be a cause for another round of CI runs.

Alternative way of printing the current function name (C99):

G_debug(3, "Vect_get_point_in_poly_isl(): the hard way");

G_debug(3, "%s: the hard way", __func__);

[metzm]

Using bool pre-C23 requires #include <stdbool.h>, this is probably done in one of the other includes, so using this in .c file not .h file, shouldn't necessarily be a cause for another round of CI runs.

I trust the bool usage because "All checks have passed".

Alternative way of printing the current function name (C99):

G_debug(3, "Vect_get_point_in_poly_isl(): the hard way");

G_debug(3, "%s: the hard way", __func__);

That would be a massive PR, replacing the actual function name with __func__ in all GRASS debug messages.

[nilason]

Using bool pre-C23 requires #include <stdbool.h>, this is probably done in one of the other includes, so using this in .c file not .h file, shouldn't necessarily be a cause for another round of CI runs.

I trust the bool usage because "All checks have passed".

Me too, that's why I just mentioned it.

Alternative way of printing the current function name (C99):

G_debug(3, "Vect_get_point_in_poly_isl(): the hard way");

G_debug(3, "%s: the hard way", __func__);

That would be a massive PR, replacing the actual function name with __func__ in all GRASS debug messages.

I obviously didn't mean to change all the cases in repo, but it would have been possible to use in this new code. But, nothing urgent, so no need to re-run the CI's for that. Next time if you like, just wanted to mention the possibility.

[petrasovaa]

v.overlay currently doesn't have a test. @metzm do you think this is covered by any other test? I think it is reasonable to ask to have at least some test coverage, especially now when tests can be easily generated by AI.

[metzm]

v.overlay currently doesn't have a test. @metzm do you think this is covered by any other test? I think it is reasonable to ask to have at least some test coverage, especially now when tests can be easily generated by AI.

I don't think this is covered by another test. These warnings WARNING: Vect_get_point_in_poly_isl(): collapsed area are typically ignored. Writing a test would be easy, you just count the number of centroids in the output. The difficult part is creating a test vector with these edge cases (really tiny polygons). I have to ask my colleagues how they created this beast ;-)

Such a test could come with another open pull request for v.overlay #7370 that also deals with cleaning up small areas. The new test would then count areas and centroids with default settings and with the new option to remove small areas enabled. This would require to merge this PR first to get correct results.