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

### Related tags

Geospatial Library polylabel

## polylabel

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.

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

• #### Here is a C# implementation

Using the PriorityQueue implementation from BlackWasp:

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

{
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

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

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()`:

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

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:

weyer1.json.txt

opened by stonetip 6

### 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

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?

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

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?

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

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

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

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

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

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

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

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)
###### 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.
###### 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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1 Jan 24, 2022
5 Jul 5, 2022
###### Point in Polygon with 'CGAL'

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

3 Nov 12, 2021