A fast algorithm for finding the pole of inaccessibility of a polygon (in JavaScript and C++)

Overview

polylabel Build Status

A fast algorithm for finding polygon pole of inaccessibility, the most distant internal point from the polygon outline (not to be confused with centroid), implemented as a JavaScript library. Useful for optimal placement of a text label on a polygon.

It's an iterative grid algorithm, inspired by paper by Garcia-Castellanos & Lombardo, 2007. Unlike the one in the paper, this algorithm:

  • guarantees finding global optimum within the given precision
  • is many times faster (10-40x)

How the algorithm works

This is an iterative grid-based algorithm, which starts by covering the polygon with big square cells and then iteratively splitting them in the order of the most promising ones, while aggressively pruning uninteresting cells.

  1. Generate initial square cells that fully cover the polygon (with cell size equal to either width or height, whichever is lower). Calculate distance from the center of each cell to the outer polygon, using negative value if the point is outside the polygon (detected by ray-casting).
  2. Put the cells into a priority queue sorted by the maximum potential distance from a point inside a cell, defined as a sum of the distance from the center and the cell radius (equal to cell_size * sqrt(2) / 2).
  3. Calculate the distance from the centroid of the polygon and pick it as the first "best so far".
  4. Pull out cells from the priority queue one by one. If a cell's distance is better than the current best, save it as such. Then, if the cell potentially contains a better solution that the current best (cell_max - best_dist > precision), split it into 4 children cells and put them in the queue.
  5. Stop the algorithm when we have exhausted the queue and return the best cell's center as the pole of inaccessibility. It will be guaranteed to be a global optimum within the given precision.

image

JavaScript Usage

Given polygon coordinates in GeoJSON-like format and precision (1.0 by default), Polylabel returns the pole of inaccessibility coordinate in [x, y] format.

var p = polylabel(polygon, 1.0);

TypeScript

TypeScript type definitions are available via npm install --save @types/polylabel.

C++ Usage

It is recommended to install polylabel via mason. You will also need to install its dependencies: geometry.hpp and variant.

#include <mapbox/polylabel.hpp>

int main() {
    mapbox::geometry::polygon<double> polygon = readPolygon(); // Get polygon data from somewhere.
    mapbox::geometry::point<double> p = mapbox::polylabel(polygon, 1.0);
    return 0;
}

Ports to other languages

Comments
  • Here is a C# implementation

    Here is a C# implementation

    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);
               }
            }
         }
      }
    
    opened by scastria 11
  • PolyLabel still places center point outside some irregular geometry

    PolyLabel still places center point outside some irregular geometry

    This little bit of code is quite useful. Thanks for writing it. I found a few cases where the irregular geometries of feature still confuse the code and place the "center point" outside of the geometry. I'm tracking down the reason as best I can and will submit whatever I find. If one of you guys has a quicker route to fixing this than myself I'd gladly take that as well.

    I'll attach the GEOJSON that can be used to reproduce the error.

    {"name":"dbo.ICOOR","type":"FeatureCollection" ,"features":[ { "type":"Feature","geometry":{"type":"Polygon","coordinates":[[[-123.273625,49.3610850000001],[-123.274258,49.3612660000001],[-123.274735,49.3615070000001],[-123.275153,49.3618200000001],[-123.275551,49.3621840000001],[-123.275737,49.3624590000001],[-123.275855,49.3627650000001],[-123.275861,49.363161],[-123.275747,49.36373],[-123.275616,49.3642110000001],[-123.275515,49.3645940000001],[-123.275492,49.3652970000001],[-123.275571,49.365664],[-123.275692,49.365957],[-123.275861,49.3662830000001],[-123.275912,49.3665710000001],[-123.275873,49.366651],[-123.275788,49.366861],[-123.275402,49.3678840000001],[-123.274969,49.36884],[-123.274615,49.3694620000001],[-123.274245,49.369894],[-123.273793,49.3703180000001],[-123.273244,49.370706],[-123.272172,49.3712310000001],[-123.271774,49.3714650000001],[-123.270416,49.37244],[-123.269881,49.372918],[-123.269704,49.373143],[-123.269568,49.373548],[-123.269857,49.374645],[-123.269254,49.3769780000001],[-123.268865,49.37735],[-123.268584,49.3775890000001],[-123.268146,49.3780580000001],[-123.266609,49.3792130000001],[-123.266404,49.3795280000001],[-123.266175,49.38041],[-123.265683,49.381131],[-123.265409,49.3813930000001],[-123.265189,49.381537],[-123.264846,49.3817000000001],[-123.26402,49.381945],[-123.263704,49.3821070000001],[-123.263567,49.3822060000001],[-123.263004,49.382594],[-123.262497,49.3830090000001],[-123.260943,49.3836880000001],[-123.260792,49.383787],[-123.260148,49.38449],[-123.258927,49.3853740000001],[-123.256158,49.3881960000001],[-123.255899,49.388673],[-123.255725,49.3896810000001],[-123.255442,49.3909140000001],[-123.254923,49.3918590000001],[-123.25454,49.392391],[-123.254239,49.392724],[-123.253676,49.393175],[-123.251794,49.3944300000001],[-123.25141,49.3947370000001],[-123.250274,49.3958620000001],[-123.249319,49.3967850000001],[-123.248542,49.3975500000001],[-123.247974,49.398107],[-123.247079,49.399123],[-123.246589,49.3994760000001],[-123.246152,49.399776],[-123.245835,49.4000160000001],[-123.245701,49.4001950000001],[-123.244317,49.40112],[-123.243809,49.401526],[-123.243411,49.4019500000001],[-123.242975,49.4024890000001],[-123.242682,49.4028960000001],[-123.242316,49.403486],[-123.242146,49.403828],[-123.242058,49.404096],[-123.242061,49.4044870000001],[-123.242172,49.405069],[-123.242238,49.405353],[-123.242264,49.4056130000001],[-123.242234,49.4058910000001],[-123.241764,49.4068940000001],[-123.241679,49.407319],[-123.241601,49.407744],[-123.241536,49.408215],[-123.244232,49.408173],[-123.270444999,49.404594],[-123.288034999,49.3848850000001],[-123.29908,49.375297],[-123.292535,49.359316],[-123.278626,49.3433360000001],[-123.277399,49.332682],[-123.262672,49.325225],[-123.23281,49.33428],[-123.200085,49.3308180000001],[-123.168373999,49.3238880000001],[-123.168541,49.3252180000001],[-123.168093,49.327971],[-123.168086,49.3288900000001],[-123.168071,49.3299030000001],[-123.168059999,49.3306890000001],[-123.168061999,49.331625],[-123.168042999,49.3325980000001],[-123.168028,49.333602],[-123.168017,49.3341960000001],[-123.168012,49.3346320000001],[-123.168019999,49.3355250000001],[-123.167994,49.3364250000001],[-123.167985,49.337278],[-123.168019999,49.3382270000001],[-123.168000999,49.338779],[-123.168003999,49.3391280000001],[-123.167997,49.340024],[-123.167974,49.3409340010001],[-123.167980999,49.341038],[-123.16799,49.3418370000001],[-123.167951,49.3419180000001],[-123.167921999,49.3420180000001],[-123.167897,49.3421890000001],[-123.167872,49.342403],[-123.167857,49.34272],[-123.167846999,49.343146],[-123.167832,49.3436790000001],[-123.167825999,49.3439060000001],[-123.167823,49.344071],[-123.169233,49.3439450000001],[-123.170898,49.343805],[-123.172093,49.3437750000001],[-123.173445,49.343779],[-123.174346,49.343826],[-123.174979,49.3438770000001],[-123.175864,49.343965],[-123.176912,49.3441060000001],[-123.177676,49.344212],[-123.178424,49.3442930000001],[-123.178881,49.344322],[-123.180088,49.3443990000001],[-123.180892,49.344454],[-123.183843,49.3446250000001],[-123.184365,49.344663],[-123.184884,49.3447350000001],[-123.185302,49.3448270000001],[-123.18575,49.344954],[-123.186236,49.3451400000001],[-123.186726,49.34536],[-123.186884,49.3454390000001],[-123.186984,49.3455000000001],[-123.187256,49.3456620000001],[-123.188453,49.3463530000001],[-123.18909,49.346664],[-123.189654,49.346906],[-123.190154,49.3470590000001],[-123.190794,49.3471740000001],[-123.191209,49.3472110000001],[-123.191936,49.3472640000001],[-123.19271,49.347311],[-123.193393,49.3473320000001],[-123.194114,49.3473450000001],[-123.194774,49.3473530000001],[-123.195678,49.3473320000001],[-123.196919,49.347296],[-123.198029,49.347223],[-123.201019,49.3469940000001],[-123.208887,49.3462920000001],[-123.210389,49.3461810000001],[-123.211813,49.3461300000001],[-123.212968,49.3461380000001],[-123.214052,49.346168],[-123.215005,49.346215],[-123.216324,49.346303],[-123.218799,49.346456],[-123.219965,49.3465630000001],[-123.222048,49.3468930000001],[-123.226067,49.347905],[-123.227447,49.348235],[-123.228339,49.348419],[-123.23052,49.348756],[-123.231702,49.348895],[-123.232933,49.3490100000001],[-123.233671,49.3490610000001],[-123.234529,49.349071],[-123.235306,49.349052],[-123.236338,49.349022],[-123.237242,49.3489860000001],[-123.238815,49.3489500000001],[-123.240605,49.3490200000001],[-123.241229,49.349116],[-123.242498,49.3494590000001],[-123.243132,49.349697],[-123.24466,49.3503480000001],[-123.245512,49.3506840000001],[-123.246005,49.35085],[-123.246716,49.351028],[-123.247141,49.351096],[-123.248264,49.3511810000001],[-123.249486,49.351251],[-123.250037,49.351275],[-123.251496,49.351575],[-123.252554,49.352064],[-123.253178,49.352445],[-123.253504,49.352742],[-123.25377,49.3530620000001],[-123.253922,49.3534050000001],[-123.254148,49.354139],[-123.253527,49.357174],[-123.253524,49.357559],[-123.25372,49.3581780000001],[-123.253804,49.3583550000001],[-123.25433,49.359093],[-123.254732,49.3594780000001],[-123.255032,49.359722],[-123.255319,49.359922],[-123.255492,49.3600200000001],[-123.256348,49.360347],[-123.257288,49.3606980000001],[-123.258604,49.361204],[-123.259263,49.3614260000001],[-123.260125,49.3616170000001],[-123.260657,49.361722],[-123.262496,49.3619770000001],[-123.263399,49.3620230000001],[-123.263785,49.362019],[-123.265522,49.361878],[-123.267278,49.3616260000001],[-123.268055,49.361477],[-123.270194,49.361094],[-123.271137,49.360977],[-123.27196,49.360958],[-123.273625,49.3610850000001]]]},"properties":{"IAREA_ID":8472}} ]}

    bug 
    opened by codingwithoutbugs 8
  • Algorithm fails for some polygons

    Algorithm fails for some polygons

    Algorithm does not take into account segment endpoints. Look at this polygon: http://www.openstreetmap.org/way/300726657 As you see there are lines that are far from center but are considered very close to it, for instance one that is formed by http://www.openstreetmap.org/node/3048346767 and http://www.openstreetmap.org/node/3048346768

    question 
    opened by andreynovikov 7
  • Unexpected results

    Unexpected results

    Hi folks!

    Not sure I'm doing something wrong, but getting a somewhat unexpected pole of inaccessibility for a U-like polygon. The red marker is what I got from polylabel():

    screen shot 2018-05-10 at 4 01 03 pm

    I tried changing the precision but didn't notice any difference. On version 1.0.2 of the js library.

    Polygon in GeoJSON below:

    {
      "type":"Polygon",
      "coordinates":[
        [
          [
            -122.4018194,
            37.7909722
          ],
          [
            -122.400845,
            37.7910899
          ],
          [
            -122.4007716,
            37.7907104
          ],
          [
            -122.4009996,
            37.7906828
          ],
          [
            -122.4010399,
            37.7908913
          ],
          [
            -122.4011183,
            37.7908818
          ],
          [
            -122.4011134,
            37.7908568
          ],
          [
            -122.4011346,
            37.7908122
          ],
          [
            -122.4011967,
            37.7908047
          ],
          [
            -122.4012382,
            37.7908417
          ],
          [
            -122.401243,
            37.7908668
          ],
          [
            -122.4015202,
            37.7908333
          ],
          [
            -122.4014781,
            37.7906159
          ],
          [
            -122.4017442,
            37.7905837
          ],
          [
            -122.4018194,
            37.7909722
          ]
        ]
      ]
    }
    

    Any thoughts/feedback appreciated!

    Thanks,

    question 
    opened by pedro 6
  • Counterintuitive Pole Placement

    Counterintuitive Pole Placement

    I may misunderstand the way this algorithm works. And from reading the description, it is an impressive bit of work. However, this is just one of several examples where the output is unexpected:

    image

    weyer1.json.txt

    opened by stonetip 6
  • Add default module export to support ES6 Typescript

    Add default module export to support ES6 Typescript

    Typescript

    This enables a more ES6 import syntax used in Typescript.

    current

    import * as polylabel from 'polylabel'
    

    with default module

    import polylabel from 'polylabel'
    

    Babel

    Babel does compile the default module without this addition, however it does some hacks in the background.

    sample code

    import polylabel from 'polylabel';
    
    const polygon = [[[3116,3071],[3118,3068],[3108,3102],[3751,927]]]
    polylabel(polygon)
    

    babel converts into

    var _polylabel = require('polylabel');
    
    var _polylabel2 = _interopRequireDefault(_polylabel);
    
    function _interopRequireDefault(obj) { return obj && obj.__esModule ? obj : { default: obj }; }
    var polygon = [[[3116, 3071], [3118, 3068], [3108, 3102], [3751, 927]]]; 
    (0, _polylabel2.default)(polygon);
    
    opened by DenisCarriere 6
  • Incorrect decision of best solution

    Incorrect decision of best solution

    image

    Building at (53.89824, 27.56131)

    Upper cell got center point very close to polygon line.

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

    .. and filtered out by precision

    opened by serbod 5
  • What kind of file format does polylabel accept?

    What kind of file format does polylabel accept?

    I tried with geojson and json. But it seems that content of the file must be specific. Can you please add in documentation what kind of file format does it accept and also, how should the content format be?

    opened by gabimoncha 5
  • Return Best Distance

    Return Best Distance

    This PR is a proposed solution for #2.

    A fresh take on #14.

    Summary of changes:

    • Adds a distance property to returned values.
    • Includes changes to tests.
    • Guarantee that the distance value is 0, not -0, for degenerate cases.
    • Guarantee that there is always a value included for distance (never gets undefined).
    • Stick to ES5 (use var, and Object.assign in tests instead of rest/spread).

    Kindly requesting review @mourner @fil . Feedback welcome!

    🙏

    opened by curran 4
  • Ported to R?

    Ported to R?

    Do you know if this algorithm been ported to R, and if not, if there are plans to do so?

    Since the code is C++, it should be relatively easy to do – although I admit having no clear idea of how to do it.

    opened by briatte 4
  • Infinite loop while trying to find most promising cell from the queue

    Infinite loop while trying to find most promising cell from the queue

    In my case I've used us states map and have been trying to find the center of the states, but faced an infinite loop for some states because getCentroidCell used to return structure like { x: NaN, y: NaN, d: NaN . . . }

    so when it tried to compare the numbers later, it had always the negative result

    bug 
    opened by demkalkov 4
  • The algorithm is not guaranteed to find a point that lies inside the polygon

    The algorithm is not guaranteed to find a point that lies inside the polygon

    If all initial cells' centers lie outside the polygon and the precision is large enough for cell.max - bestCell.d <= precision to return true, the algorithm will return a point that lies outside the polygon, as it will not drill down into any more cells. It may be desirable to always return a point that lies inside the polygon, regardless of the precision parameter, as the dimensions of the polygons may not be known in advance to provide an adequate precision value. This can easily be fixed by preventing drill down only if the best cell is inside the polygon:

    // do not drill down further if there's no chance of a better solution
    // but only if our best cell's center is inside the polygon
    if ((cell.max - bestCell.d <= precision) && (bestCell.d > 0))
    

    A more efficient approach would be to monitor every potential cell's .d value and set a flag to true if one that is positive is found. This would reduce the amount of cells with low potential in the queue but to implement it we need to introduce an extra variable in the algorithm.

    enhancement 
    opened by alex-petkov 0
  • endless loop when vertices lie on the edge joining two other vertices

    endless loop when vertices lie on the edge joining two other vertices

    Please try computing the center of this polygon:

    [[-300,-2700],[-300,-2300],[500,-2300],[500,-2872.0465053408525],[1101.5552288470312,-2029.8691849550091],[2700,-2500],[-631187325490387600,-5259894379086588000],[2806.312641353628,64.56245288212403],[1500,500],[1500,700],[1300,700],[-1500,0],[-800,-1400],[-800,-2000],[-1000,-2000],[-1000,-3000]]
    
    bug 
    opened by ecoinomist 0
  • Trivially reject subcells before adding them to the priority queue

    Trivially reject subcells before adding them to the priority queue

    Frequently a newly-generated quadrant cell cannot possibly beat the current best cell. In this case, it's wasted work to add it to the priority queue just so we can discard it later, particularly since the priority queue needs to find the correct location to insert the new cell.

    In my test cases, a pure square polygon yielded zero savings. All other test cases were able to immediately skip adding between 10% and 72% of the newly-generated cells. The maximium queue size shrunk down to 87%-4.5% of the original queue size.

    To accomplish this in the JavaScript version, you'd change

    cellQueue.push(new Cell(cell.x - h, cell.y - h, h, polygon));
    cellQueue.push(new Cell(cell.x + h, cell.y - h, h, polygon));
    cellQueue.push(new Cell(cell.x - h, cell.y + h, h, polygon));
    cellQueue.push(new Cell(cell.x + h, cell.y + h, h, polygon));
    

    to something like this:

    newCell0 = new Cell(cell.x - h, cell.y - h, h, polygon));
    newCell1 = new Cell(cell.x + h, cell.y - h, h, polygon));
    newCell2 = new Cell(cell.x - h, cell.y + h, h, polygon));
    newCell3 = new Cell(cell.x + h, cell.y + h, h, polygon));
    
    if (newcell0.max > bestCell.d) cellQueue.push(newCell0);
    if (newcell1.max > bestCell.d) cellQueue.push(newCell1);
    if (newcell2.max > bestCell.d) cellQueue.push(newCell2);
    if (newcell3.max > bestCell.d) cellQueue.push(newCell3);
    

    Verified with my test cases.

    enhancement 
    opened by hollasch 1
  • Use single covering cell instead of initial tiling to simplify setup

    Use single covering cell instead of initial tiling to simplify setup

    The current algorithm computes the bounding box of the polygon, selects the shortest side as the initial cell size, and then tiles the polygon with 1×1, 1×N or N×1 of these initial cells. Following that, it seeds the current best interior distance with the better of the bounds center or the polygon centroid.

    Since the remainder of the algorithm has all the machinery to test and subdivide cells, we can rely on this and simplify the initial setup.

    Instead of choosing the smallest side of the bounding box, choose the largest side, to compute the square cell that fully and tightly contains the polygon. Initialize the priority queue with this single cell, eliminating the tiling code from the original algorithm. Further, because this cell contains the bounds center (which is identical to the rectangular bounds center), you do not need to compute and choose between two initial best values — just use the polygon centroid distance.

    Tested in my implementation, and works great.

    enhancement 
    opened by hollasch 1
  • Possible optimization for getSegDistSq

    Possible optimization for getSegDistSq

    I believe there are some promising optimizations of the getSegDistSq() function. In the examples below, let's consider the current polylabel getSegDistSq() function as solution A.

    Solution B

    Instead of storing a polygon as an array of vertices, you can store it as an array of segments, where each segment has the following data:

    struct PolygonSegment {
        Point2D  Q;
        Vector2D D;
        double   invNormDSquared; // = 1/dot(D,D)
    }
    
    double getSegDistSq(const Point2D& P, const PolygonSegment& segment) {
        Vector2D V = P - segment.Q;
        double t = dot(V,segment.D) * segment.invNormDSquared;
    
        if (t >= 1)
            V -= segment.D;
        else if (t > 0)
            V -= t*segment.D;
    
        return dot(V,V);
    }
    

    Alternatively without the geometry types:

    struct PolygonSegment {
        double Qx, Qy;
        double Dx, Dy;
        double invNormDSquared;
    }
    
    double getSegDistSq(double Px, double Py, const PolygonSegment& segment) {
        auto Vx = Px - segment.Qx;
        auto Vy = Py - segment.Qy;
        auto t = ((Vx * segment.Dx) + (Vy * segment.Dy)) * segment.invNormDSquared;
    
        if (t >= 1) {
            Vx -= segment.Dx;
            Vy -= segment.Dy;
        } else if (t > 0) {
            Vx -= t * segment.Dx;
            Vy -= t * segment.Dy;
        }
    
        return Vx*Vx + Vy*Vy;
    }
    

    (I can show the derivation better if you're interested.)

    This saves up to 5 additions, 1 multiply and 1 division over the original solution (solution A).

    Solution C

    Note that you can shave off one additional multiply for all cases by increasing segment data by two doubles each:

    struct PolygonSegment {
        double Qx, Qy; // Starting point of segment
        double Ex, Ey; // End point of segment
        double Ux, Uy; // Unit vector along segment direction
        double length; // Segment length
    }
    
    double getSegDistSq(double Px, double Py, const PolygonSegment& segment) {
        auto Vx = Px - segment.Qx;
        auto Vy = Py - segment.Qy;
        auto t = (Vx * segment.Ux) + (Vy * segment.Uy);
    
        if (t > segment.length) {
            Vx = Px - segment.Ex;
            Vy = Px - segment.Ey;
        } else if (t > 0) {
            Vx -= t * segment.Ux;
            Vy -= t * segment.Uy;
        }
    
        return Vx*Vx + Vy*Vy;
    }
    

    This saves an additional multiply compared to solution B.

    Solution D

    But wait! We're trying to find the shortest distance between the point and the polygon. So I can ignore all segments that can't yield a closer intersection than the one I already have. We can achieve this by considering a circle that encloses the line segment, to determine us the closest possible distance between the point and the segment, like so:

    struct PolygonSegment {
        double Qx, Qy; // Starting point of segment
        double Ex, Ey; // End point of segment
        double Ux, Uy; // Unit vector along segment direction
        double length; // Segment length
    }
    
    double getSegDistSq(double Px, double Py, const PolygonSegment& segment, double& closestSoFarSquared) {
        auto Vx = Px - segment.Qx;
        auto Vy = Py - segment.Qy;
    
        double minDistSquared = Vx*Vx + Vy*Vy - segment.lengthSquared;
        if (minDistSquared > closesSoFarSquared)
            return;
    
        auto t = (Vx * segment.Ux) + (Vy * segment.Uy);
    
        if (t > segment.length) {
            Vx = Px - segment.Ex;
            Vy = Px - segment.Ey;
        } else if (t > 0) {
            Vx -= t * segment.Ux;
            Vy -= t * segment.Uy;
        }
    
        auto dSquared = Vx*Vx + Vy*Vy;
        if (dSquared < closesSoFarSquared)
            closesSoFarSquared = dSquared;
    
        return dSquared;
    }
    

    This leaves us with a tiny savings over solution A in the near case (1 division, 1 or 0 additions). However, the common case of being too far to hope for a closer solution leaves us with a significant savings: 5 additions instead of 11, 2 multiplies instead of 6, and no divisions compared to the orginal 1. For complex polygons, there's definite promise, and for simple polygons, we're still a bit ahead -- the worst we did was trade an addition for a division. Still pretty good!

    Results

    Here's the final tally of operations for each approach:

    | Version | Condition | Adds | Mults | Divs |
    |---------+-----------+------+-------+------|
    |   A     |   t < 0   |  11  |   6   |  1   |
    |   B     |   t < 0   |   6  |   5   |  0   |
    |   C     |   t < 0   |   6  |   4   |  0   |
    |   D     |   t < 0   |  10  |   6   |  0   |
    |---------+-----------+------+-------+------|
    |   A     | 0 < t < 1 |  13  |   8   |  1   |
    |   B     | 0 < t < 1 |   8  |   7   |  0   |
    |   C     | 0 < t < 1 |   8  |   6   |  0   |
    |   D     | 0 < t < 1 |  12  |   8   |  0   |
    |---------+-----------+------+-------+------|
    |   A     |   1 < t   |  10  |   6   |  1   |
    |   B     |   1 < t   |   7  |   5   |  0   |
    |   C     |   1 < t   |   7  |   4   |  0   |
    |   D     |   1 < t   |  11  |   6   |  0   |
    |---------+-----------+------+-------+------|
    |   D     |    Far    |   5  |   2   |  0   |
    +---------+-----------+------+-------+------+
    

    Note that I haven't tested this, either for correctness or for performance. If anyone has the requisite setup to try this out in the project code, I'm hopeful that the results would be good.

    question 
    opened by hollasch 6
  • Added polylabel article, updated README

    Added polylabel article, updated README

    Article copied from the original post on the Mapbox blog: https://blog.mapbox.com/a-new-algorithm-for-finding-a-visual-center-of-a-polygon-7c77e6492fbc This includes a new directory, images/, to hold the article figures. Some very light editing has also been done.

    Also updated the project README to reference the article. Cleaned up the formatting a bit to be useable in plaintext form.

    opened by hollasch 7
Releases(v1.1.0)
Owner
Mapbox
Mapbox is the location data platform for mobile and web applications. We're changing the way people move around cities and explore our world.
Mapbox
a language for fast, portable data-parallel computation

Halide Halide is a programming language designed to make it easier to write high-performance image and array processing code on modern machines. Halid

Halide 5.2k Jan 5, 2023
Interactive, thoroughly customizable maps in native Android, iOS, macOS, Node.js, and Qt applications, powered by vector tiles and OpenGL

Mapbox GL Native A C++ library that powers customizable vector maps in native applications on multiple platforms by taking stylesheets that conform to

Mapbox 4.2k Jan 6, 2023
C++ implementation of R*-tree, an MVR-tree and a TPR-tree with C API

libspatialindex Author: Marios Hadjieleftheriou Contact: [email protected] Revision: 1.9.3 Date: 10/23/2019 See http://libspatialindex.org for full doc

null 633 Dec 28, 2022
2D and 3D map renderer using OpenGL ES

Tangram ES Tangram ES is a C++ library for rendering 2D and 3D maps from vector data using OpenGL ES. It is a counterpart to Tangram. This repository

Tangram 750 Jan 1, 2023
Terrain Analysis Using Digital Elevation Models (TauDEM) software for hydrologic terrain analysis and channel network extraction.

TauDEM (Terrain Analysis Using Digital Elevation Models) is a suite of Digital Elevation Model (DEM) tools for the extraction and analysis of hydrolog

David Tarboton 191 Dec 28, 2022
A command line toolkit to generate maps, point clouds, 3D models and DEMs from drone, balloon or kite images. 📷

An open source command line toolkit for processing aerial drone imagery. ODM turns simple 2D images into: Classified Point Clouds 3D Textured Models G

OpenDroneMap 3.9k Jan 6, 2023
Computational geometry and spatial indexing on the sphere

S2 Geometry Library Overview This is a package for manipulating geometric shapes. Unlike many geometry libraries, S2 is primarily designed to work wit

Google 1.9k Dec 31, 2022
A C++17 image representation, processing and I/O library.

Selene Selene is a C++17 image representation, processing, and I/O library, focusing on ease of use and a clean, modern, type-safe API. Overview: Brie

Michael Hofmann 286 Oct 26, 2022
A lean, efficient, accurate geohash encoder and decoder library implemented in C

Geohash encoder/decoder in C A lean, efficient, accurate geohash encoder and decoder library implemented in C. It does not depend on the C standard li

Christopher Wellons 20 Nov 20, 2022
A library of distance and occlusion generation routines

Distance/Occlusion Library + Tool From left to right: original, signed distance with zero at 0.5, red/green SDF, delta vectors to closest boundary poi

Andrew Willmott 105 Nov 2, 2022
✔️The smallest header-only GUI library(4 KLOC) for all platforms

Welcome to GUI-lite The smallest header-only GUI library (4 KLOC) for all platforms. 中文 Lightweight ✂️ Small: 4,000+ lines of C++ code, zero dependenc

null 6.6k Jan 8, 2023
Pole-like Objects Mapping and Long-Term Robot Localization in Dynamic Urban Scenarios

Long Term Localization Pole-like Objects Mapping and Long-Term Robot Localization is an algorithm that makes robot or UAV locate itself in Dynamic Urb

Networked Robotics and Sytems Lab 90 Nov 30, 2022
The goal of insidesp is to do fast point in polygon classification, the sp way.

insidesp The goal of insidesp is to do fast point in polygon classification, the sp way. We are comparing a few ways of implementing this, essentially

diminutive 2 Nov 12, 2021
A C++ implementation of the graph algorithm of finding all strongly connected components with DFS

SinkSCC A C++ implementation of the graph algorithm of finding all strongly connected components with DFS. Details Each SCC of a graph can be treated

Schrodinger ZHU Yifan 2 Nov 2, 2021
Implementation of Dijkstra's algorithm for finding the shortest paths between nodes in a graph using Priority Queue data structure in C

Implementation of Dijkstra's algorithm for finding the shortest paths between nodes in a graph using Priority Queue data structure in C. File "airport

Artem Kolpakov 1 Jan 24, 2022
Shamir’s Secret Sharing Algorithm: Shamir’s Secret Sharing is an algorithm in cryptography created by Adi Shamir. The main aim of this algorithm is to divide secret that needs to be encrypted into various unique parts.

Shamir-s-Secret-Sharing-Algorithm-Cryptography Shamir’s Secret Sharing Algorithm: Shamir’s Secret Sharing is an algorithm in cryptography created by A

Pavan Ananth Sharma 5 Jul 5, 2022
Point in Polygon with 'CGAL'

The goal of insidecgal is to provide fast point in polygon lookup with cgal.

diminutive 3 Nov 12, 2021
GeneralPolygonClipper - General Polygon Clipper, aka GPC, version 2.33

GeneralPolygonClipper General Polygon Clipper, aka GPC, version 2.33 See GPC-README.pdf for more official information about this library. GPC was disc

Rick Brewster 43 Dec 25, 2022
Fast Binary Encoding is ultra fast and universal serialization solution for C++, C#, Go, Java, JavaScript, Kotlin, Python, Ruby, Swift

Fast Binary Encoding (FBE) Fast Binary Encoding allows to describe any domain models, business objects, complex data structures, client/server request

Ivan Shynkarenka 654 Jan 2, 2023
Lee Thomason 299 Dec 10, 2022