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.

[petrasovaa]

Claude's review:

Vect_find_poly_centroid_cog is declared non-static and has no public header entry. It's only used inside poly.c. Either mark it static (and drop the forward declaration at the top of the file) or, if the intent is to
expose it as a peer of Vect_find_poly_centroid, add it to the vector API header. Currently it's a globally visible symbol with no prototype in a header — a namespace-pollution hazard and a linker footgun if another
translation unit ever defines a same-named function.

[metzm]

Vect_find_poly_centroid_cog is now in the public header alongside Vect_find_poly_centroid.

There are arguments to not use at all Vect_find_poly_centroid in Vect_get_point_in_poly_isl and use only Vect_find_poly_centroid_cog followed by the hard way:

• Vect_find_poly_centroid does not consider isles inside an area
• Vect_find_poly_centroid calculates the center of gravity of the outer polygon line, not the area. For very thin areas, this can cause the centroid to be dangerously close to the polygon line, such that the distance to the polygon line is not safely representable due to floating point precision limits. Vect_find_poly_centroid_cog is in this regard more robust.

[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.