Last active
February 27, 2024 16:18
-
-
Save dfaivre/acfef42cdbf411555956e9eba65dd30d to your computer and use it in GitHub Desktop.
PolyLabel (pole of inaccessibility of a polygon) - C#/NTS
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
using System; | |
using System.Collections.Generic; | |
using System.Linq; | |
using GeoAPI.Geometries; | |
using NetTopologySuite.Geometries; | |
namespace Foo | |
{ | |
/// <summary> | |
/// Finds the largest area inside a polygon (ie, for a label) | |
/// | |
/// See: https://github.com/mapbox/polylabel/blob/master/polylabel.js | |
/// </summary> | |
public static class PolyLabel | |
{ | |
public static SearchCell FindMax( | |
IPolygon polygon, double precision = 1) | |
{ | |
var boundingBox = polygon.EnvelopeInternal; | |
var cellSize = System.Math.Min(boundingBox.Width, boundingBox.Height); | |
var searchSize = cellSize / 2; | |
// I'm not sure why the source has a search size of zero for the centroid cell. Oh well? | |
// DMF 2018-05-16 | |
var bestCell = new SearchCell(polygon.Centroid, searchSize: 0, polygon); | |
var cellCandidates = new SortedList<SearchCell, SearchCell>(); | |
for (var x = boundingBox.MinX; x < boundingBox.MaxX; x += cellSize) | |
{ | |
for (var y = boundingBox.MinY; y < boundingBox.MaxY; y += cellSize) | |
{ | |
var cell = new SearchCell(new Point(x + searchSize, y + searchSize), searchSize, polygon); | |
cellCandidates.Add(cell, cell); | |
} | |
} | |
while (cellCandidates.Any()) | |
{ | |
// pick the most promising cell from the queue | |
var cell = cellCandidates.Keys[0]; | |
cellCandidates.RemoveAt(0); | |
// update the best cell if we found a better one | |
if (cell.CentroidDistance > bestCell.CentroidDistance) | |
{ | |
bestCell = cell; | |
} | |
// do not drill down further if there's no chance of a better solution | |
if (cell.MaxPotentialDistance - bestCell.CentroidDistance <= precision) continue; | |
// else: split the cell and serach | |
var splitCells = cell.Split(); | |
foreach (var splitCell in splitCelrrrls) | |
{ | |
cellCandidates.Add(splitCell, splitCell); | |
} | |
} | |
return bestCell; | |
} | |
public class SearchCell : IComparable<SearchCell> | |
{ | |
private readonly IPolygon _searchPolygon; | |
private static readonly double Sqrt2 = System.Math.Sqrt(2); | |
public SearchCell(IPoint point, double searchSize, double centroidDistance, double maxPotentialDistance) | |
{ | |
Point = point; | |
SearchSize = searchSize; | |
CentroidDistance = centroidDistance; | |
MaxPotentialDistance = maxPotentialDistance; | |
} | |
public SearchCell(IPoint point, double searchSize, IPolygon searchPolygon) | |
{ | |
_searchPolygon = searchPolygon; | |
Point = point; | |
SearchSize = searchSize; | |
CentroidDistance = PointToPolygonDistanceSigned(point, searchPolygon); | |
MaxPotentialDistance = CentroidDistance + searchSize * Sqrt2; | |
} | |
public IPoint Point { get; } | |
public double SearchSize { get; } | |
public double CentroidDistance { get; } | |
public double MaxPotentialDistance { get; } | |
public SearchCell[] Split() | |
{ | |
var splitSearchSize = SearchSize / 2; | |
SearchCell CreateSplitCell(double x, double y) | |
{ | |
return new SearchCell( | |
new Point(x, y), splitSearchSize, _searchPolygon); | |
} | |
return new[] | |
{ | |
CreateSplitCell(Point.X - splitSearchSize, Point.Y - splitSearchSize), | |
CreateSplitCell(Point.X + splitSearchSize, Point.Y - splitSearchSize), | |
CreateSplitCell(Point.X + splitSearchSize, Point.Y + splitSearchSize), | |
CreateSplitCell(Point.X - splitSearchSize, Point.Y + splitSearchSize), | |
}; | |
} | |
public int CompareTo(SearchCell other) => MaxPotentialDistance > other.MaxPotentialDistance ? -1 : 1; | |
} | |
private static double PointToPolygonDistanceSigned( | |
IPoint point, IPolygon polygon) | |
{ | |
var inside = polygon.Intersects(point); | |
var minDistance = | |
polygon.InteriorRings | |
.Select(r => r.Distance(point)) | |
.Concat(new[] {polygon.ExteriorRing.Distance(point)}) | |
.Min(); | |
return (inside ? 1 : -1) * minDistance; | |
} | |
} | |
} |
Just FYI there is another C# implementation: mapbox/polylabel#26
Oh, and another? https://github.com/Epi-Info/Epi-Info-Community-Edition/blob/1965e5291980858a63d577f7224def90f61a7101/EpiDashboard/Mapping/PolyLabel.cs
Hi @dfaivre
Thanks for sharing this. Any idea why this is not working for the following polygon?
I am showing the marker at the center point of the Search Cell. I would like the marker to be in the visual center of the polygon.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
C# (with NetTopologySuite) implementation of mapbox/polylabel (A fast algorithm for finding the pole of inaccessibility of a polygon).
Notes
SortedList
for the priority queue. I did this so that I could simulate apop
by just removing the first element. ASortedSet
could also work (and would remove the need for the duplicateKey, Value
), but would seem to need to remove the element by reference instead of position.SearchCell
implements a custom comparer so that no two can be considered equal (SortedList
does not allow duplicates, so we trick it).