fix skyums algorithm

[ljwolf]

As I understood it when implementing it, the minimum bounding circle in Skyum's algorithm is found by removing boundary points on the convex hull. Specifically, we start with the largest circle formed by adjacent triples and remove the middle of the triple when it subtends an obtuse angle, recalculate angles/radii for the left and right part of the triple, and continue. A similar strategy focused on acute triples leads to my pointset paring algorithm in #155.

It has been a while since I read the original article, but given my comment about "confusing as hell defaults", I think #75 stems from what "lexicographic order" means in the algorithm:

I remember finding this confusing at the time---what is meant by "lexicographic" over the two keys? Do we sort first by angle, then radius? or radius, then angle?

Eleven years ago, I ordered them in the way they appear in the paper. However, as #75 notes, this can lead to cases where the radii of the MBC is too small when the termination criteria is hit on the angle.

Rereading the original paper, I think the proof for Lemma 2 implies a sort on radius, then angle:

the only way radius(a,b,c) is certainly maximal in this case is if we sort on radius first. When this is done, we're picking the largest radius triple that subtends an acute angle. In the case of tied radii/circle size, we pick the largest acute circle. This is guaranteed to contain the input points, so the MBC won't be sometimes too small anymore.

import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
import shapely
from pointpats.centrography import minimum_bounding_circle as skyum

points = np.array([(1, 2), (2, 3), (2, 2), (2, 1), (3, 4), (3, 3), (3, 1), (4, 4)]).astype(float)
geoms = gpd.GeoSeries.from_xy(*points.T)

# skyum outcome
(x,y), r = skyum(points) 

ax = geoms.plot(figsize=(12, 12))
gpd.GeoSeries([shapely.Point(x, y).buffer(r)]).boundary.plot(ax=ax, color='r', linewidth=1)
geoms.to_frame().dissolve().minimum_bounding_circle().boundary.plot(ax=ax, color='g', linewidth=1)
plt.show()

[martinfleis]

We can either fix our implementation in this way or pass everything onto shapely as discussed in #75. Might be worth merging this before doing the switch so there's at least in git history a reference implementation of Skyum that is correct.

[codecov]

Codecov Report

❌ Patch coverage is 70.00000% with 3 lines in your changes missing coverage. Please review.
✅ Project coverage is 73.5%. Comparing base (a798e9e) to head (194b654).
⚠️ Report is 5 commits behind head on main.

Files with missing lines	Patch %	Lines
pointpats/centrography.py	70.0%	3 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff           @@
##            main    #206     +/-   ##
=======================================
+ Coverage   72.9%   73.5%   +0.6%     
=======================================
  Files         12      12             
  Lines       2118    2111      -7     
=======================================
+ Hits        1543    1551      +8     
+ Misses       575     560     -15     

Files with missing lines	Coverage Δ
pointpats/centrography.py	86.0% <70.0%> (+2.0%)	⬆️

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[ljwolf]

Between minimum_bounding_circle and minimum_bounding_radius in shapely, we have everything we need at an acceptable performance. I will prepare that PR after this.

We could keep it as a separate skyum() function if helpful, but I am wary of "keeping things around" just for the sake of it. What do other maintainers think?

[jGaboardi]

but I am wary of "keeping things around" just for the sake of it. What do other maintainers think?

+1. If no longer used, let's remove. We can always put it back later if it turns out to be useful in future.