[BUG] Invalid polygons/closed linestrings when smoothing is enabled

[netthier]

When smoothing is enabled, ContourBuilder::contours and probably ContourBuilder::isobands can end up generating invalid (self-intersecting) polygons.
E.g.:
Without smoothing:

With smoothing:

Overlaid:

I think that this is invalid according to geo-types, and many other geometry standards don't allow this either.
Specifically,

The interior of the polygon must be a connected point-set. That is, any two distinct points in the interior must admit a curve between these two that lies in the interior.

And,

The exterior and interior rings must be valid LinearRings (see LineString).

Further, a closed LineString must not self-intersect. Note that its validity is not enforced, and operations and predicates are undefined on invalid LineStrings.

@mthh You can find the areas the screenshots are from in the dataset I sent you recently at about 557500, 5932680 in its CRS, contours were generated at elevation -5.

[mthh]

I'm a but surprised, are you getting this result with version 0.13?
I can't find this kind of degeneracy when generating contour polygons with your dataset and the smooth option (and I can't find the one you pointed out, but perhaps I generated the contour using a different threshold from yours).

The spikes that are created without smooth don't normally overflow to make self intersecting polygons but instead give the following result (yellow with smooth option and brown without smooth option):

I used to check results of contour-isobands with some validity tests (https://github.com/mthh/geo-validity-check), maybe the same should be done for this crate.

EDIT: I just saw your edit with the threshold your are using and I can reproduce the bug.

[mthh]

Any idea what the problem is here?

[mthh]

I've clipped your raster onto the problem area and I tried to generate the same contours with d3-contour and with the previous versions of this crate (using the threshold value of -5) and it looks like this is a bug that was introduced in 0.13.0.

See the following screenshots:

d3-contour / contour-rs 0.12.1, no smoothing

contour-rs 0.13.0, no smoothing

d3-contour / contour-rs 0.12.1, with smoothing

contour-rs 0.13.0, with smoothing

I'm currently investigating the exact cause.

[mthh]

The culprit change seems to be b32d16b.

[netthier]

Oof yeah that seems like a mistake with the as usize casts. I must've thought those were some other int type at the time. Ah lol, just noticed you proposed those changes. I should've looked more closely 😅
It's a bit worrying that seemingly no test caught this though.
I forgot why I changed it though 😔 The method isn't precise enough with f32, should probably go the same route as area then. i.e. always use f64 for intermediate results 🤔

[mthh]

I've tinkered a bit with the code of this function, notably using f64 for intermediate values.

With something like (point.x as f64 * 2.0 + point.y as f64 * ((self.dx + 1) * 4) as f64) as usize it doesn't seem to trigger the bug of #12 neither the bug discussed here. But we could go with another index computation (I've made other positive attempts in this direction, but we might as well stay close to the code it's ported to).

Anyway, feel free to propose a fix you're confident with ;)

It's a bit worrying that seemingly no test caught this though.

Well, now we have a test case to avoid this kind of regression in the future (if you allow me to take a small piece of your raster to include in the tests? otherwise I'll manually create values to reproduce it).
But yes the test case are still very simple and are mostly pulled from d3-contour. I have two small datasets in the fixtures folder that could also be used into test cases to catch future regressions.

[netthier]

That solution makes sense, maybe #17 could be an opportunity to explore alternative ways later 🤔
Might also make sense to yank 0.13.0. Feel free to take like a 10x10 array of plain, non-georeferenced values for the tests.

[mthh]

Alright, I will push a fix and publish a new version.

Feel free to take like a 10x10 array of plain, non-georeferenced values for the tests.

Thanks!

[mthh]

Thanks for the data, i took exactly 10x10 values on the problematic area, non-georeferenced, and I modified a bit the values before including them in a test case.