mapbox / polylabel

A fast algorithm for finding the pole of inaccessibility of a polygon (in JavaScript and C++)
Other
1.44k stars 151 forks source link

Here is a C# implementation #26

Closed scastria closed 7 years ago

scastria commented 7 years ago

Using the PriorityQueue implementation from BlackWasp:

  using System;
  using System.Collections.Generic;
  using PriorityQueue;

  namespace SkiaDemo1
  {
     public class PolyLabel
     {
        private const float EPSILON = 1E-8f;

        public static float[] GetPolyLabel(float[][][] polygon, float precision = 1f)
        {
           //Find the bounding box of the outer ring
           float minX = 0, minY = 0, maxX = 0, maxY = 0;
           for (int i = 0; i < polygon[0].Length; i++) {
              float[] p = polygon[0][i];
              if (i == 0 || p[0] < minX)
                 minX = p[0];
              if (i == 0 || p[1] < minY)
                 minY = p[1];
              if (i == 0 || p[0] > maxX)
                 maxX = p[0];
              if (i == 0 || p[1] > maxY)
                 maxY = p[1];
           }

           float width = maxX - minX;
           float height = maxY - minY;
           float cellSize = Math.Min(width, height);
           float h = cellSize / 2;

           //A priority queue of cells in order of their "potential" (max distance to polygon)
           PriorityQueue<float,Cell> cellQueue = new PriorityQueue<float, Cell>();

           if (FloatEquals(cellSize, 0))
              return new[] { minX, minY };

           //Cover polygon with initial cells
           for (float x = minX; x < maxX; x += cellSize) {
              for (float y = minY; y < maxY; y += cellSize) {
                 Cell cell = new Cell(x + h, y + h, h, polygon);
                 cellQueue.Enqueue(cell.Max, cell);
              }
           }

           //Take centroid as the first best guess
           Cell bestCell = GetCentroidCell(polygon);

           //Special case for rectangular polygons
           Cell bboxCell = new Cell(minX + width / 2, minY + height / 2, 0, polygon);
           if (bboxCell.D > bestCell.D)
              bestCell = bboxCell;

           int numProbes = cellQueue.Count;

           while (cellQueue.Count > 0) {
              //Pick the most promising cell from the queue
              Cell cell = cellQueue.Dequeue();

              //Update the best cell if we found a better one
              if (cell.D > bestCell.D) {
                 bestCell = cell;
              }

              //Do not drill down further if there's no chance of a better solution
              if (cell.Max - bestCell.D <= precision)
                 continue;

              //Split the cell into four cells
              h = cell.H / 2;
              Cell cell1 = new Cell(cell.X - h, cell.Y - h, h, polygon);
              cellQueue.Enqueue(cell1.Max, cell1);
              Cell cell2 = new Cell(cell.X + h, cell.Y - h, h, polygon);
              cellQueue.Enqueue(cell2.Max, cell2);
              Cell cell3 = new Cell(cell.X - h, cell.Y + h, h, polygon);
              cellQueue.Enqueue(cell3.Max, cell3);
              Cell cell4 = new Cell(cell.X + h, cell.Y + h, h, polygon);
              cellQueue.Enqueue(cell4.Max, cell4);
              numProbes += 4;
           }

           return (new[] { bestCell.X, bestCell.Y });
        }

        //Signed distance from point to polygon outline (negative if point is outside)
        private static float PointToPolygonDist(float x, float y, float[][][] polygon)
        {
           bool inside = false;
           float minDistSq = float.PositiveInfinity;

           for (int k = 0; k < polygon.Length; k++) {
              float[][] ring = polygon[k];

              for (int i = 0, len = ring.Length, j = len - 1; i < len; j = i++) {
                 float[] a = ring[i];
                 float[] b = ring[j];

                 if ((a[1] > y != b[1] > y) && (x < (b[0] - a[0]) * (y - a[1]) / (b[1] - a[1]) + a[0]))
                    inside = !inside;

                 minDistSq = Math.Min(minDistSq, GetSegDistSq(x, y, a, b));
              }
           }

           return ((inside ? 1 : -1) * (float)Math.Sqrt(minDistSq));
        }

        //Get squared distance from a point to a segment
        private static float GetSegDistSq(float px, float py, float[] a, float[] b)
        {
           float x = a[0];
           float y = a[1];
           float dx = b[0] - x;
           float dy = b[1] - y;

           if (!FloatEquals(dx, 0) || !FloatEquals(dy, 0)) {
              float t = ((px - x) * dx + (py - y) * dy) / (dx * dx + dy * dy);
              if (t > 1) {
                 x = b[0];
                 y = b[1];
              } else if (t > 0) {
                 x += dx * t;
                 y += dy * t;
              }
           }
           dx = px - x;
           dy = py - y;
           return (dx * dx + dy * dy);
        }

        //Get polygon centroid
        private static Cell GetCentroidCell(float[][][] polygon)
        {
           float area = 0;
           float x = 0;
           float y = 0;
           float[][] ring = polygon[0];

           for (int i = 0, len = ring.Length, j = len - 1; i < len; j = i++) {
              float[] a = ring[i];
              float[] b = ring[j];
              float f = a[0] * b[1] - b[0] * a[1];
              x += (a[0] + b[0]) * f;
              y += (a[1] + b[1]) * f;
              area += f * 3;
           }
           if (FloatEquals(area, 0))
              return (new Cell(ring[0][0], ring[0][1], 0, polygon));
           return (new Cell(x / area, y / area, 0, polygon));
        }

        private static bool FloatEquals(float a, float b)
        {
           return (Math.Abs(a - b) < EPSILON);
        }

        private class Cell
        {
           public float X { get; private set; }
           public float Y { get; private set; }
           public float H { get; private set; }
           public float D { get; private set; }
           public float Max { get; private set; }

           public Cell(float x, float y, float h, float[][][] polygon)
           {
              X = x;
              Y = y;
              H = h;
              D = PointToPolygonDist(X, Y, polygon);
              Max = D + H * (float)Math.Sqrt(2);
           }
        }
     }
  }
4sfaloth commented 7 years ago

First of all, thanks for the mapbox team for making their work available, and thanks to @scastria for implementing a C# version. I am trying to use this in my program, but am having some issues. Problem is, I believe I don't really understand what I should feed in as input in the polygon array. From my understanding polygon array should be a N x M x 2 array In the first dimension, I believe that the first entry (polygon[0]) corresponds to the outer contours. But what about the others? Is it supposed to be all child contours?

Thanks in advance

mourner commented 7 years ago

@4sfaloth exactly; the first is the outer ring and the others are holes. This corresponds to how polygons are encoded in GeoJSON.

@scastria thank you for providing the port! Feel free to publish a repository with it and the link to the original implementation. Closing since this is not an issue.

DarinMacRae commented 6 years ago

@scastria Just curious if you ended up putting this out on NuGet? I didn't see it directly but thought I'd ask.

scastria commented 6 years ago

I did not. sorry.

cdc-dpbrown commented 6 years ago

That is great! Thank you. Great timing. We are in the middle of porting polylabel over to be used with ESRI layers now. https://github.com/Epi-Info/Epi-Info-Community-Edition/commits/master

matheuscschenfeld commented 4 years ago

My polygon is composed by a List<Position> that contain vertices and in each one I have the coordinates latitude and longitude. My data base saves the polygon in Postgis, for example:

POLYGON((-53.26865374982731 -29.007876847119025,-53.26582133710758 -29.01012872728904,-53.26642215192692 -29.011348475222643,-53.26852500379459 -29.013619044560283,-53.27112138212101 -29.010729220379005))

How can I convert it to float[][][] polygon? Does anybody know a solution for that?

DaveInCaz commented 4 years ago

@matheuscschenfeld the code I use does something like this:

    float[][] ConvertPolygonToFloatArray()
    {
        // this is a float[][]
        var polygon = new float[number_of_points_in_polygon][];

        var pointCount = 0;
        foreach (var point in points_in_polygon)
        {
            polygon[pointCount] = new float[2];
            polygon[pointCount][0] = point.Data.X;
            polygon[pointCount][1] = point.Data.Y;
            pointCount++;
        }

        return polygon;
    }
S4NT14G0 commented 4 years ago

Here's an implementation for those using Unity3D. I wrote the heap implementation myself, but a priority queue will work as well.


using System;
using System.Collections.Generic;
using System.Linq;
using IsoFrame.Pathfinding;
using UnityEngine;

namespace PolyLabel
{
    public class PolyLabelNet
    {
        private const float EPSILON = 1E-8f;

        public static Vector2Int FindPoleOfIsolation(List<Vector2Int> edgePoints, List<Vector2Int> internalPoints, float precision=1)
        {
            int minX = int.MaxValue;
            int minY = int.MaxValue;
            int maxX = int.MinValue;
            int maxY = int.MinValue;

            foreach (Vector2Int edgePoint in edgePoints)
            {
                minX = minX > edgePoint.x ? edgePoint.x : minX;
                minY = minY > edgePoint.y ? edgePoint.y : minY;
                maxX = maxX < edgePoint.x ? edgePoint.x : maxX;
                maxY = maxY < edgePoint.y ? edgePoint.y : maxY;
            }

            var width = maxX - minX;
            var height = maxY - minY;
            var cellSize = Mathf.Min(width, height);
            var h = cellSize / 2;

            if (cellSize == 0) return new Vector2Int(minX, minY);

            Heap<Cell> cellQueue = new Heap<Cell>();

            for (float x = minX; x < maxX; x += cellSize)
            {
                for (float  y = minY; y < maxY; y += cellSize)
                {
                    cellQueue.Insert(new Cell(Mathf.RoundToInt(x + h), Mathf.RoundToInt(y + h), h, edgePoints, internalPoints));
                }
            }

            var bestCell = GetCentroidCell(edgePoints, internalPoints);

            var bboxCell = new Cell(minX + width / 2, minY + height / 2, 0, edgePoints, internalPoints);
            bestCell = bboxCell.D > bestCell.D ? bboxCell: bestCell;

            var numProbs = cellQueue.Size();

            while(cellQueue.Size() > 0)
            {
                var cell = cellQueue.Pop();

                if (cell.D > bestCell.D)
                {
                    bestCell = cell;
                }

                if (cell.Max - bestCell.D <= precision) continue;

                h = cell.H / 2;
                cellQueue.Insert(new Cell(cell.X - h, cell.Y - h, h, edgePoints, internalPoints));
                cellQueue.Insert(new Cell(cell.X + h, cell.Y - h, h, edgePoints, internalPoints));
                cellQueue.Insert(new Cell(cell.X - h, cell.Y + h, h, edgePoints, internalPoints));
                cellQueue.Insert(new Cell(cell.X + h, cell.Y + h, h, edgePoints, internalPoints));
                numProbs += 4;
            }

            return new Vector2Int(bestCell.X, bestCell.Y);
        }

        private static Cell GetCentroidCell(List<Vector2Int> edgePoints, List<Vector2Int> internalPoints)
        {
            float area = 0;
            float x = 0;
            float y = 0;

            for (int i = 0, len = edgePoints.Count, j = len - 1; i < len; j = i++)
            {
                var a = edgePoints[i];
                var b = edgePoints[j];
                var f = a.x * b.y - b.x * a.y;
                x += (a.x + b.x) * f;
                y += (a.y + b.y) * f;
                area += f * 3;
            }
            if (area == 0)
                return new Cell(edgePoints[0].x, edgePoints[0].y, 0, edgePoints, internalPoints);

            return new Cell(Mathf.RoundToInt(x / area), Mathf.RoundToInt(y / area), 0, edgePoints, internalPoints);
        }

        private static bool FloatEquals(float a, float b)
        {
            return (Math.Abs(a - b) < EPSILON);
        }

        private class Cell : IComparable<Cell>
        {
            public int X { get; private set; }
            public int Y { get; private set; }
            public int H { get; private set; }
            public float D { get; private set; }
            public float Max { get; private set; }

            public List<Vector2Int> InternalPoints { get; private set; }

            public Cell(int x, int y, int h, List<Vector2Int> edgePoints, List<Vector2Int> internalPoints)
            {
                X = x;
                Y = y;
                H = h;
                InternalPoints = internalPoints;
                D = PointToPolygonDist(X, Y, edgePoints);
                Max = D + H * (float)Math.Sqrt(2);
            }

            private float PointToPolygonDist(int x, int y, List<Vector2Int> edgePoints)
            {
                float minDist = float.PositiveInfinity;

                Vector2Int cellCenter = new Vector2Int(x,y);

                foreach (Vector2Int edgePoint in edgePoints)
                {
                    var distance = Vector2Int.Distance(cellCenter, edgePoint);
                    minDist = minDist < distance ? minDist : distance;                        
                }

                if (!InternalPoints.Contains(cellCenter))
                    minDist = minDist * -1;

                return minDist;
            }

            public int CompareTo(Cell other)
            {
                float comparisonValue = this.Max - other.Max;

                if (comparisonValue > 0)
                    return 1;
                else if (FloatEquals(comparisonValue, 0))
                    return 0;
                else
                    return -1;    
            }
        }

    }   
}
ilirvg commented 3 years ago

@S4NT14G0 I am not sure what internalPoints are supposed to be. In my case I have polygons that are constructed by only external vectors, is it safe to discard internalPoints?

MarcoEder commented 3 years ago

To anyone who is not familiar with the GeoJSON format like me

The method must be used like this (assuming that there is only a single contour without any holes):

float[ ][ ][ ] polygon = new float[1][ ][ ]; polygon[0] = ConvertPolygonToFloatArray();

@matheuscschenfeld the code I use does something like this:

    float[][] ConvertPolygonToFloatArray()
    {
        // this is a float[][]
        var polygon = new float[number_of_points_in_polygon][];

        var pointCount = 0;
        foreach (var point in points_in_polygon)
        {
            polygon[pointCount] = new float[2];
            polygon[pointCount][0] = point.Data.X;
            polygon[pointCount][1] = point.Data.Y;
            pointCount++;
        }

        return polygon;
    }
minhth1010 commented 3 years ago

Here's an implementation for those using Unity3D. I wrote the heap implementation myself, but a priority queue will work as well.

using System;
using System.Collections.Generic;
using System.Linq;
using IsoFrame.Pathfinding;
using UnityEngine;

namespace PolyLabel
{
    public class PolyLabelNet
    {
        private const float EPSILON = 1E-8f;

        public static Vector2Int FindPoleOfIsolation(List<Vector2Int> edgePoints, List<Vector2Int> internalPoints, float precision=1)
        {
            int minX = int.MaxValue;
            int minY = int.MaxValue;
            int maxX = int.MinValue;
            int maxY = int.MinValue;

            foreach (Vector2Int edgePoint in edgePoints)
            {
                minX = minX > edgePoint.x ? edgePoint.x : minX;
                minY = minY > edgePoint.y ? edgePoint.y : minY;
                maxX = maxX < edgePoint.x ? edgePoint.x : maxX;
                maxY = maxY < edgePoint.y ? edgePoint.y : maxY;
            }

            var width = maxX - minX;
            var height = maxY - minY;
            var cellSize = Mathf.Min(width, height);
            var h = cellSize / 2;

            if (cellSize == 0) return new Vector2Int(minX, minY);

            Heap<Cell> cellQueue = new Heap<Cell>();

            for (float x = minX; x < maxX; x += cellSize)
            {
                for (float  y = minY; y < maxY; y += cellSize)
                {
                    cellQueue.Insert(new Cell(Mathf.RoundToInt(x + h), Mathf.RoundToInt(y + h), h, edgePoints, internalPoints));
                }
            }

            var bestCell = GetCentroidCell(edgePoints, internalPoints);

            var bboxCell = new Cell(minX + width / 2, minY + height / 2, 0, edgePoints, internalPoints);
            bestCell = bboxCell.D > bestCell.D ? bboxCell: bestCell;

            var numProbs = cellQueue.Size();

            while(cellQueue.Size() > 0)
            {
                var cell = cellQueue.Pop();

                if (cell.D > bestCell.D)
                {
                    bestCell = cell;
                }

                if (cell.Max - bestCell.D <= precision) continue;

                h = cell.H / 2;
                cellQueue.Insert(new Cell(cell.X - h, cell.Y - h, h, edgePoints, internalPoints));
                cellQueue.Insert(new Cell(cell.X + h, cell.Y - h, h, edgePoints, internalPoints));
                cellQueue.Insert(new Cell(cell.X - h, cell.Y + h, h, edgePoints, internalPoints));
                cellQueue.Insert(new Cell(cell.X + h, cell.Y + h, h, edgePoints, internalPoints));
                numProbs += 4;
            }

            return new Vector2Int(bestCell.X, bestCell.Y);
        }

        private static Cell GetCentroidCell(List<Vector2Int> edgePoints, List<Vector2Int> internalPoints)
        {
            float area = 0;
            float x = 0;
            float y = 0;

            for (int i = 0, len = edgePoints.Count, j = len - 1; i < len; j = i++)
            {
                var a = edgePoints[i];
                var b = edgePoints[j];
                var f = a.x * b.y - b.x * a.y;
                x += (a.x + b.x) * f;
                y += (a.y + b.y) * f;
                area += f * 3;
            }
            if (area == 0)
                return new Cell(edgePoints[0].x, edgePoints[0].y, 0, edgePoints, internalPoints);

            return new Cell(Mathf.RoundToInt(x / area), Mathf.RoundToInt(y / area), 0, edgePoints, internalPoints);
        }

        private static bool FloatEquals(float a, float b)
        {
            return (Math.Abs(a - b) < EPSILON);
        }

        private class Cell : IComparable<Cell>
        {
            public int X { get; private set; }
            public int Y { get; private set; }
            public int H { get; private set; }
            public float D { get; private set; }
            public float Max { get; private set; }

            public List<Vector2Int> InternalPoints { get; private set; }

            public Cell(int x, int y, int h, List<Vector2Int> edgePoints, List<Vector2Int> internalPoints)
            {
                X = x;
                Y = y;
                H = h;
                InternalPoints = internalPoints;
                D = PointToPolygonDist(X, Y, edgePoints);
                Max = D + H * (float)Math.Sqrt(2);
            }

            private float PointToPolygonDist(int x, int y, List<Vector2Int> edgePoints)
            {
                float minDist = float.PositiveInfinity;

                Vector2Int cellCenter = new Vector2Int(x,y);

                foreach (Vector2Int edgePoint in edgePoints)
                {
                    var distance = Vector2Int.Distance(cellCenter, edgePoint);
                    minDist = minDist < distance ? minDist : distance;                        
                }

                if (!InternalPoints.Contains(cellCenter))
                    minDist = minDist * -1;

                return minDist;
            }

            public int CompareTo(Cell other)
            {
                float comparisonValue = this.Max - other.Max;

                if (comparisonValue > 0)
                    return 1;
                else if (FloatEquals(comparisonValue, 0))
                    return 0;
                else
                    return -1;    
            }
        }

    }   
}

Hi @S4NT14G0 , thanks for your code!

I can't find IsoFrame.Pathfinding namespace and how to create Polygon from (vertices, triangles)?

Have a good day!

NooBiTTo commented 11 months ago

Here's an implementation for those using Unity3D. I wrote the heap implementation myself, but a priority queue will work as well.

using System;
using System.Collections.Generic;
using System.Linq;
using IsoFrame.Pathfinding;
using UnityEngine;

namespace PolyLabel
{
    public class PolyLabelNet
    {
        private const float EPSILON = 1E-8f;

        public static Vector2Int FindPoleOfIsolation(List<Vector2Int> edgePoints, List<Vector2Int> internalPoints, float precision=1)
        {
            int minX = int.MaxValue;
            int minY = int.MaxValue;
            int maxX = int.MinValue;
            int maxY = int.MinValue;

            foreach (Vector2Int edgePoint in edgePoints)
            {
                minX = minX > edgePoint.x ? edgePoint.x : minX;
                minY = minY > edgePoint.y ? edgePoint.y : minY;
                maxX = maxX < edgePoint.x ? edgePoint.x : maxX;
                maxY = maxY < edgePoint.y ? edgePoint.y : maxY;
            }

            var width = maxX - minX;
            var height = maxY - minY;
            var cellSize = Mathf.Min(width, height);
            var h = cellSize / 2;

            if (cellSize == 0) return new Vector2Int(minX, minY);

            Heap<Cell> cellQueue = new Heap<Cell>();

            for (float x = minX; x < maxX; x += cellSize)
            {
                for (float  y = minY; y < maxY; y += cellSize)
                {
                    cellQueue.Insert(new Cell(Mathf.RoundToInt(x + h), Mathf.RoundToInt(y + h), h, edgePoints, internalPoints));
                }
            }

            var bestCell = GetCentroidCell(edgePoints, internalPoints);

            var bboxCell = new Cell(minX + width / 2, minY + height / 2, 0, edgePoints, internalPoints);
            bestCell = bboxCell.D > bestCell.D ? bboxCell: bestCell;

            var numProbs = cellQueue.Size();

            while(cellQueue.Size() > 0)
            {
                var cell = cellQueue.Pop();

                if (cell.D > bestCell.D)
                {
                    bestCell = cell;
                }

                if (cell.Max - bestCell.D <= precision) continue;

                h = cell.H / 2;
                cellQueue.Insert(new Cell(cell.X - h, cell.Y - h, h, edgePoints, internalPoints));
                cellQueue.Insert(new Cell(cell.X + h, cell.Y - h, h, edgePoints, internalPoints));
                cellQueue.Insert(new Cell(cell.X - h, cell.Y + h, h, edgePoints, internalPoints));
                cellQueue.Insert(new Cell(cell.X + h, cell.Y + h, h, edgePoints, internalPoints));
                numProbs += 4;
            }

            return new Vector2Int(bestCell.X, bestCell.Y);
        }

        private static Cell GetCentroidCell(List<Vector2Int> edgePoints, List<Vector2Int> internalPoints)
        {
            float area = 0;
            float x = 0;
            float y = 0;

            for (int i = 0, len = edgePoints.Count, j = len - 1; i < len; j = i++)
            {
                var a = edgePoints[i];
                var b = edgePoints[j];
                var f = a.x * b.y - b.x * a.y;
                x += (a.x + b.x) * f;
                y += (a.y + b.y) * f;
                area += f * 3;
            }
            if (area == 0)
                return new Cell(edgePoints[0].x, edgePoints[0].y, 0, edgePoints, internalPoints);

            return new Cell(Mathf.RoundToInt(x / area), Mathf.RoundToInt(y / area), 0, edgePoints, internalPoints);
        }

        private static bool FloatEquals(float a, float b)
        {
            return (Math.Abs(a - b) < EPSILON);
        }

        private class Cell : IComparable<Cell>
        {
            public int X { get; private set; }
            public int Y { get; private set; }
            public int H { get; private set; }
            public float D { get; private set; }
            public float Max { get; private set; }

            public List<Vector2Int> InternalPoints { get; private set; }

            public Cell(int x, int y, int h, List<Vector2Int> edgePoints, List<Vector2Int> internalPoints)
            {
                X = x;
                Y = y;
                H = h;
                InternalPoints = internalPoints;
                D = PointToPolygonDist(X, Y, edgePoints);
                Max = D + H * (float)Math.Sqrt(2);
            }

            private float PointToPolygonDist(int x, int y, List<Vector2Int> edgePoints)
            {
                float minDist = float.PositiveInfinity;

                Vector2Int cellCenter = new Vector2Int(x,y);

                foreach (Vector2Int edgePoint in edgePoints)
                {
                    var distance = Vector2Int.Distance(cellCenter, edgePoint);
                    minDist = minDist < distance ? minDist : distance;                        
                }

                if (!InternalPoints.Contains(cellCenter))
                    minDist = minDist * -1;

                return minDist;
            }

            public int CompareTo(Cell other)
            {
                float comparisonValue = this.Max - other.Max;

                if (comparisonValue > 0)
                    return 1;
                else if (FloatEquals(comparisonValue, 0))
                    return 0;
                else
                    return -1;    
            }
        }

    }   
}

Hi @S4NT14G0 , thanks for your code!

I can't find IsoFrame.Pathfinding namespace and how to create Polygon from (vertices, triangles)?

Have a good day!

Hi, I have the same issue. Did you find any solution?