[WIP] Filter holes from alpha shape returns #457

ljwolf · 2022-02-01T17:12:06Z

Noted by @cornhundred in #456, we've got some surprising behavior for alpha shapes with holes. This attempts to fix this in the alpha_shape() function by (1) checking if there are any "holes" in the shapely polygons and, if there are, (2) we only keep polygons that contain an input point.

This solves the problem OK, but I think there's a bit of an issue:

Note that this results in two separate polygons-with-holes on the left, but on the right, we retain one of the holes because a point falls within it... I think this happens because the "correct" polygon for that hole would have a zero-area "cut in" that only includes that point. When we use shapely.ops.polygonize, we lose that cut in, and thus omit the vertex from the alpha shape.

Consider this box with a "cut in" from [1,1] to [.5,.5]:

lines = [
         LineString([[0, 0], [0, 1]]), # 1
         LineString([[0, 1], [1, 1]]), # 2
         LineString([[0, 0], [1, 0]]), # 6
         LineString([[1, 0], [1, 1]]), # 3
         LineString([[1, 1], [0.5, 0.5]]), # 4 
         LineString([[0.5, 0.5], [1, 1]]), # 5
        ]

upon polygonization, we lose that cut-in:

from shapely.ops import polygonize
geom = polygonize(lines)
geom.wkt

geom is now 'POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))' with no mention of the cut-in.
This is good, since a polygon with this kind of cut-in is an invalid polygon, but it complicates our testing here... not all holes are empty.

The other way to do this would be to do an all-pairs within check for the polygons using geoms.sindex.query_bulk(). I was hoping the point-in-poly woudl be faster, but... but you never know. I'll check that & re-submit.

codecov · 2022-02-02T15:25:26Z

Codecov Report

Merging #457 (ac17fc5) into master (8aa08b1) will increase coverage by 0.1%.
The diff coverage is 100.0%.

@@           Coverage Diff            @@
##           master    #457     +/-   ##
========================================
+ Coverage    78.7%   78.8%   +0.1%     
========================================
  Files         122     122             
  Lines       12884   12910     +26     
========================================
+ Hits        10142   10177     +35     
+ Misses       2742    2733      -9

Impacted Files	Coverage Δ
libpysal/cg/alpha_shapes.py	`61.6% <100.0%> (+2.1%)`	⬆️
libpysal/cg/tests/test_ashapes.py	`97.2% <100.0%> (+0.7%)`	⬆️
libpysal/_version.py	`40.7% <0.0%> (+2.7%)`	⬆️

ljwolf · 2022-02-02T15:28:28Z

OK, I've switched to a zorder search method:

Extract the exterior ring of all shapes geoms. Call this shells.
Compute shells.sindex.query_bulk(geoms, predicate='within') to determine which shells each geometry falls within.
Construct the "within" graph adjacency matrix from these edges.
Compute the row sum of this adjacency matrix. This represents the number of shells in which a geometry falls, self-inclusive, also known as the "z-order".
Select all geometries with an odd z-order, because all holes have an even z-order and all non-holes are odd. Reprising the code comment:

   # A z-order of 1 means the polygon is only inside if its own exterior. This means it's not a hole.
   # A z-order of 2 means the polygon is inside of exactly one other exterior. Because
   #   the hull generation method is restricted to be planar, this means the polygon is a hole. 
   # In general, an even z-order means that the polygon is always exactly matched to one exterior, 
   #   plus some number of intermediate exterior-hole pairs. Therefore, the polygon is a hole.
   # In general, an odd z-order means that there is an uneven number of exteriors. 
   #   This means the polygon is not a hole.

ljwolf · 2022-02-02T15:38:44Z

@cornhundred, this should be good to go, and solves your reference case. Is there any other known pattern you'd want to test on?

cornhundred · 2022-02-02T16:04:10Z

Thanks @ljwolf! The random matrix with holes should cover our use case. Does it return a MultiPolygon or a list of Polygons?

Also, I started to scope out how to allow alternate alpha shape metrics (radius, distance, etc) https://github.com/cornhundred/libpysal/tree/alt_alpha_metrics I'll make a pull request when I've cleaned it up more and meet the pull request requirements (testing etc). For our case we would like to set the threshold based on max triangle length.

ljwolf · 2022-02-02T16:20:11Z

It keeps the return type as a GeoSeries. Would you prefer a single (multi)polygon?

The alpha_shape_auto returns the single shapely polygon, so keeping the return type consistent might be a good idea. However, that'd involve changing the return type of alpha_shape... which would break things.

cornhundred · 2022-02-02T16:24:14Z

Ok, keeping it the way it is should be fine. It is easy enough to calculate a MultiPolygon from it. We also want to calculate properties of the individual polygons (area, number of holes, etc) so that would be easier if it returns a list of discontiguous polygons.

darribas · 2022-02-04T18:28:35Z

Hello,

Thanks for this @ljwolf for the PR and @cornhundred for bringing it up and chipping in!

I've looked over and it seems all fine, very clever trick!

A couple of thoughts:

Would it be worth to add a test to test both options in alpha_shape?
Right now the holes are discarded, which is probably fine in most applications. Would it be worth considering using the holes to generate "holed" shapes? I'm not sure that's the name, but a MultiPolygon where inside the hole is not part of the polygon? That could also be an option (filter_holes=False/'remove'/'include', or something like that).
It'd be great to add documentation for the filter_holes parameter in alpha_shape.

cornhundred · 2022-02-05T14:13:58Z

Hey @darribas, those are great ideas, I'm happy to try and help address them. We definitely care about the holes but I'm not sure if it makes sense to pass them as separate polygons rather than derive them later from polygons.

For some background, we're working with single-cell spatial transcriptomics data - see mouse liver data release and example image below

The holes in our case can be blood vessels (or artifacts from sectioning) and alpha shapes seem like a great way to quantify them. Also here's a Colab notebook from a brain dataset where we interactively visualize our data.

martinfleis · 2022-02-07T08:57:54Z

I'm not sure that's the name, but a MultiPolygon where inside the hole is not part of the polygon?

That can still be a Polygon as long as it has a single part (with N holes in it but not another part within that hole).

ljwolf · 2022-02-07T11:51:36Z

holed shapes

The shapes we return are already holed. We currently return polygons with holes and also the holes themselves as separate extra polygons... Below shows our current behavior:

Currently, we return all of these shapes mixed together in a single GeoSeries. We don't differentiate which is the "hull" (with holes) and which are extra hole polygons. Further, you can't immediately tell which polygons represent the "hull" and which are just extra "holes," hence the z-order search.

False/'ignore'/'include'

Upon reflection, I actually we should remove the option and always filter. Two reasons:

80% use case: People looking for the alpha_shape() are probably expecting one (or more) polygons that represent the alpha hull of the input point set. In this case, returning duplicate hole polygons (with or without labelling) is giving the user "extra" information they didn't ask for, aren't interested in, and (unless labelled) can't tell apart from what they are interested in: the alpha hull.
Promote "good" patterns: If they do want them, then we should push users towards standard geopandas functionality (e.g. result.interiors.explode()) to extract holes from polygons. The information they want is already there in the holey polygons in the filtered geoseries.

ljwolf · 2022-02-07T11:54:32Z

I've added a test_holes() case to the test_ashapes.py now. If we are OK w/ filtering (no option) then I will add a note/discussion in the docstring about holes and commit. Otherwise, I have the docstring for the option ready to go, and can modify the test easily to check the options.

cornhundred · 2022-02-07T12:52:37Z

Hi @ljwolf, I agree with the always filtering logic. Since the holes are properties of specific polygons within the alpha shape we would want users to obtain holes with metadata about their 'parent' polygon, but that should be outside the scope of the alpha shape method.

cornhundred · 2022-03-02T14:58:40Z

Hi, just checking on the status of this pull request. I can try to make an additional pull request to have the option for different alpha shape distance thresholds this week so maybe it makes sense to wait on that before making a potential new release.

ljwolf · 2022-03-02T15:03:13Z

A determination about the proper behavior of the option from a reviewer (e.g. @darribas @martinfleis @sjsrey) would get this moving again. Happy to raise this on our monthly developer meeting friday if there is no assent by then!

cornhundred · 2022-03-02T15:07:19Z

Great thanks @ljwolf appreciate it.

martinfleis

I agree with @ljwolf. I believe that we should return the polygons that represent the alpha shape and not those that represent holes. As said, those are already included as interiors so we would be just duplicating the same information.

Apart from that one note in the docstring, this seems to be ready to be merged.

martinfleis · 2022-03-02T15:37:46Z

libpysal/cg/alpha_shapes.py

+        no CRS included in the object. Also note that holes are included in addition
+        to exteriors. 


This is probably not the case with _filter_holes anymore.

Ah, so, this is slightly confusing.
There are two functions:

alpha_shape(xys, alpha), which returns a GeoSeries for a given alpha and input points xys.

alpha_geoms(alpha, triangles, radii, xys), which does the same, but using the delaunay triangulation and its corresponding circumcircle radii.

I have only filtered the alpha_shape() output. alpha_geoms() still can return holes.

I didn't want to add the z-order search into alpha_geoms() because it's used repeatedly in a performance-critical loop by alpha_shape_auto(). That function constructs a single output Polygon that covers the input point set. For its particular case, the holes returned by alpha_geoms() won't matter because of the search procedure used in alpha_shape_auto(). So, if the z-order search were introduced in alpha_geoms, it would introduce an unnecessary (but slight) performance penalty.

Thinking about this, the "right" solution would be to make alpha_geoms() private to alpha_shape_auto(), but we don't know if anyone's using that function.

I propose to solve this:

make an unfiltered private version of alpha_geoms() to preserve alpha_shape_auto() speed

add a filter to alpha_geoms(), too.

deprecate the public alpha_geoms() method.

Thoughts?

It seems reasonable to me to make alpha_geoms() private since users will be directed to use alpha_shape().

alpha_geoms is not a part of public API as far as I can tell. It is not mentioned in the API documentation and it is not listed in __all__ in this file. So I'd say that we do not have to deprecate it and can change as we wish.

Awesome. Then I'll add an underscore to the name to be extra careful, and then merge this.

Expect a point release with these changes tomorrow UK time @cornhundred

ljwolf · 2022-03-03T12:49:46Z

tests are now failing in another part of the code. Looks like a networkx issue in Python 3.7. will investigate & merge if I can figure it out.

ljwolf · 2022-03-03T14:22:27Z

OK, well, networkx dropped support for python 3.7; we'll need to drop that too, then, or skip the test.

ljwolf · 2022-03-03T14:26:25Z

Or, I suppose we could pin to networkx 2.6 for Python 3.7. Will investigate.

cornhundred · 2022-03-03T15:48:04Z

Thanks!

ljwolf · 2022-03-03T16:04:38Z

No problem! 4.6.1 should be live on pypi!

ljwolf added 3 commits February 1, 2022 16:50

add filtering if alpha shape is nonempty

029b2bb

switch to zorder search

934fd3d

remove whitespace

b2651ce

revise comment

a0e91e8

ljwolf requested review from sjsrey and darribas February 2, 2022 15:37

ljwolf added 2 commits February 7, 2022 11:50

remove filter option

5024a32

add test to check holes are processed correctly

2e2bec8

fix aliases and imports

e2ae58d

ljwolf force-pushed the alpha-multipolygon branch from 171f86f to e2ae58d Compare February 7, 2022 12:02

ljwolf mentioned this pull request Feb 7, 2022

Question on Ripleys stats pysal/pointpats#66

Open

ljwolf requested a review from martinfleis March 2, 2022 15:03

martinfleis approved these changes Mar 2, 2022

View reviewed changes

update the docstrings of alpha shape functions

e19691c

pin to networkx 2.6 on python 3.7

ac17fc5

ljwolf merged commit 0a4b079 into pysal:master Mar 3, 2022

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[WIP] Filter holes from alpha shape returns #457

[WIP] Filter holes from alpha shape returns #457

ljwolf commented Feb 1, 2022

codecov bot commented Feb 2, 2022 •

edited

ljwolf commented Feb 2, 2022 •

edited

ljwolf commented Feb 2, 2022

cornhundred commented Feb 2, 2022

ljwolf commented Feb 2, 2022

cornhundred commented Feb 2, 2022

darribas commented Feb 4, 2022

cornhundred commented Feb 5, 2022 •

edited

martinfleis commented Feb 7, 2022

ljwolf commented Feb 7, 2022 •

edited

ljwolf commented Feb 7, 2022

cornhundred commented Feb 7, 2022

cornhundred commented Mar 2, 2022

ljwolf commented Mar 2, 2022

cornhundred commented Mar 2, 2022

martinfleis left a comment

martinfleis Mar 2, 2022

ljwolf Mar 2, 2022 •

edited

cornhundred Mar 2, 2022

martinfleis Mar 2, 2022

ljwolf Mar 2, 2022

cornhundred Mar 2, 2022

ljwolf commented Mar 3, 2022

ljwolf commented Mar 3, 2022 •

edited

ljwolf commented Mar 3, 2022

cornhundred commented Mar 3, 2022

ljwolf commented Mar 3, 2022

		no CRS included in the object. Also note that holes are included in addition
		to exteriors.

[WIP] Filter holes from alpha shape returns #457

[WIP] Filter holes from alpha shape returns #457

Conversation

ljwolf commented Feb 1, 2022

codecov bot commented Feb 2, 2022 • edited

Codecov Report

ljwolf commented Feb 2, 2022 • edited

ljwolf commented Feb 2, 2022

cornhundred commented Feb 2, 2022

ljwolf commented Feb 2, 2022

cornhundred commented Feb 2, 2022

darribas commented Feb 4, 2022

cornhundred commented Feb 5, 2022 • edited

martinfleis commented Feb 7, 2022

ljwolf commented Feb 7, 2022 • edited

ljwolf commented Feb 7, 2022

cornhundred commented Feb 7, 2022

cornhundred commented Mar 2, 2022

ljwolf commented Mar 2, 2022

cornhundred commented Mar 2, 2022

martinfleis left a comment

Choose a reason for hiding this comment

martinfleis Mar 2, 2022

Choose a reason for hiding this comment

ljwolf Mar 2, 2022 • edited

Choose a reason for hiding this comment

cornhundred Mar 2, 2022

Choose a reason for hiding this comment

martinfleis Mar 2, 2022

Choose a reason for hiding this comment

ljwolf Mar 2, 2022

Choose a reason for hiding this comment

cornhundred Mar 2, 2022

Choose a reason for hiding this comment

ljwolf commented Mar 3, 2022

ljwolf commented Mar 3, 2022 • edited

ljwolf commented Mar 3, 2022

cornhundred commented Mar 3, 2022

ljwolf commented Mar 3, 2022

codecov bot commented Feb 2, 2022 •

edited

ljwolf commented Feb 2, 2022 •

edited

cornhundred commented Feb 5, 2022 •

edited

ljwolf commented Feb 7, 2022 •

edited

ljwolf Mar 2, 2022 •

edited

ljwolf commented Mar 3, 2022 •

edited