Fix issues displaying polygons with holes in Shapes #6654
Conversation
This is a partial fix for napari#5673. Current limitations: - only works when `triangle` is installed (`py-triangle` in conda-forge) - currently on macOS, this works fine running in IPython but it segfaults when called from plain Python. I don’t know whether it’s just an env thing. Summary: It turns out that what we actually need is a *constrained Delaunay triangulation*[1], that is, a Delaunay triangulation that preserves specific edges from the vertices. That’s what the 'p' option in `triangle.triangulate` does *but* it doesn’t work if the edges are actually not planar — that’s why the triangulation was failing for this polygon even if triangle was installed. This PR fixes that by finding edges that appear twice in the polygon and removing them. Additionally, certain triangles from the resulting triangulation need to be removed, and this is easy to do by running `skimage.measure.points_in_poly` on the triangle centers, which are fast to compute with NumPy indexing. [1]: http://en.wikipedia.org/wiki/Constrained_Delaunay_triangulation
|
Current output when it doesn't segfault 😅 |
|
Based on playing around with it in IPython, the segfault happens when calling |
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #6654 +/- ##
==========================================
- Coverage 93.01% 92.93% -0.09%
==========================================
Files 633 635 +2
Lines 59445 59860 +415
==========================================
+ Hits 55292 55628 +336
- Misses 4153 4232 +79 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
|
So as I understand, we need a better triangulate algorithm? |
Yes. I started to figure this out while writing out this question on the mpl discourse, which I thought would be a good place since mpl can draw these polygons no problem. Unfortunately, as you can see from the answer there, there's not really a lot of good options in Python, so I think we'll have to write our own or wrap our own. Not sure yet which is better. But at least I'm happy that I understand the problem space a bit better. |
|
btw @Czaki if you can pull this down and run it on your Linux and/or Windows machine to check whether the segfault happens there, that would be useful — maybe this problem can be fixed by improving the triangle binaries... (Though I am not really holding my breath...) |
|
@jni can you clarify, are you using arm64 binary of py-triangle from conda-forge? |
yup! |
(... Yes, I aliased my micromamba to µ, which you can type as alt-m on macOS. 😃) |
|
This pull request has been mentioned on Image.sc Forum. There might be relevant details there: https://forum.image.sc/t/some-edges-of-polygon-extend-to-infinity/90454/12 |
The issue was that the triangulation input cannot have repeated vertices. I fixed it by removing repeated vertices and translating all index coordinates to the unique vertices. This involves a bit more bookkeeping but it’s fine!
|
It's working! 😭 The issue was that the triangulation input cannot have repeated vertices. I fixed it by removing repeated vertices and translating all index coordinates to the unique vertices. |
| @njit(cache=True) | ||
| def reconstruct_polygons_from_edges( | ||
| vertices: CoordinateArray, edges: EdgeArray | ||
| ) -> list[CoordinateArray2D] | list[CoordinateArray3D]: |
| # In case of an open chain, exit. | ||
| break | ||
| polygons.append([vertices[x] for x in poly]) | ||
| return polygons |
There was a problem hiding this comment.
I'm having a hard time following this algorithm. I'm wondering if we can avoid the breaks inside nested if statements, which make it difficult to understand to what level of nesting we are breaking out. One small suggestion is changing:
for i in range(n_edges):
if not used[i]:
# huge loop bodyto
for i in range(n_edges):
if used[i]:
continue
# huge loop bodySame for the incident edges loop.
I'm also wondering about "found" vs "closed". It seems to me that open chains are just ignored? In what situations can this algorithm encounter open chains?
There was a problem hiding this comment.
Actually changing the nesting locally made it clear — open chains are just added anyway?
Again, maybe none of this would be surprising if I understood what this function is supposed to do... 😅
Co-authored-by: Tim Monko <timmonko@gmail.com>
for more information, see https://pre-commit.ci
Co-authored-by: Tim Monko <timmonko@gmail.com>
|
🚀 |
|
Haha thanks for dropping by @LucaMarconato! Unfortunately there is still a bug in some situations, see vispy/vispy#2666 😭 But we'll get there! 😅 (And the bug does not exist in PartSegCore-compiled-backend or bermuda, the latter of which will be available after #7747.) |
|
This issue has been mentioned on Image.sc Forum. There might be relevant details there: https://forum.image.sc/t/upcoming-napari-0-6-0-release/108954/10 |
|
This pull request has been mentioned on Image.sc Forum. There might be relevant details there: |
References
Partial fix forCloses #5673. (!!!)Closes #1653. 🥳
Closes #1058. 🥳
Closes #7799.
Depends on #7622
Description
(Updated 29/2/2024.)
(Updated 16/4/2025.)
This PR adds three preprocessing steps to polygon data during meshing in Shapes layers:
trianglebackend, which otherwise includes these triangles.)For the VisPy step, we go directly into the Triangulation class, because the PolygonData class completely ignores input edges when building a triangulation.
Current limitations:
only works whenEdit: nope, works with plain VisPy too!triangleis installed (py-trianglein conda-forge)currently on macOS, I am getting random segfaults when calling triangle.triangulate. I don't know what triggers them, but I suspect we're just going to have to implement our own triangulation, maybe with Numba.Edit: these were caused by repeated vertices, against which triangulation algorithms are generally not robust. Uniquifying the vertices fixed the problem!I still need to filter out the edges of the polygon for drawing, and maybe also for highlighting. Currently, the edges connecting the interior and exterior polygons are visible, which looks crap.Edit: this is now done!the polygons are processed twice, once for the interiors and once for the edges.Edit: fixed by @Czaki 🥳the processing is probably a fair bit of overhead when many datasets would get no benefit. It might be worthwhile adding a flag to turn these checks on/off.Edit: fixed by @Czaki 🥳Summary of the issue and strategies:
It turns out that what we actually need is a constrained Delaunay triangulation1, that is, a Delaunay triangulation that preserves specific edges from the vertices. That’s what the 'p' option in
triangle.triangulatedoes but it doesn’t work if the edges are actually not planar — that’s why the triangulation was failing for this polygon even if triangle was installed.This PR fixes that by finding edges that appear twice in the polygon and removing them.
Additionally, certain triangles from the resulting triangulation need to be removed, and this is easy to do by running
skimage.measure.points_in_polyon the triangle centers, which are fast to compute with NumPy indexing.