PolyLabel (pole of inaccessibility of a polygon) - C#/NTS

PolyLabel.cs

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;
        }
    }
}

C# (with NetTopologySuite) implementation of mapbox/polylabel (A fast algorithm for finding the pole of inaccessibility of a polygon).

Notes

• This uses a SortedList for the priority queue. I did this so that I could simulate a pop by just removing the first element. A SortedSet could also work (and would remove the need for the duplicate Key, 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).

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.