mapbox / polylabel Goto Github PK

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

License: Other

JavaScript 44.12% Makefile 4.24% C++ 51.64%

polylabel's Introduction

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.

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).
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).
Calculate the distance from the centroid of the polygon and pick it as the first "best so far".
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.
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;
}

Ports to other languages

andrewharvey/geojson-polygon-labels (CLI)
Twista/python-polylabel (Python)
Shapely (Python)
polylabelr (R)
polylabel-rs (Rust)
polylabel-java (Java)
php-polylabel (PHP)
dart_polylabel (Dart)
ruby-polylabel (Ruby)
Polylabel.jl (Julia)

polylabel's People

Contributors

Stargazers

Watchers

Forkers

johndpope cuulee asw456 hoestag templeblock developmentseed maturion annewallace ksobon petabytekr andrewharvey diamond0cat giohappy liusongwei alantgeo testbigorg gijs zhihuagui wsf1214 jalarrea rubythonode rejdasuki arturpasevic mapycz lixuanxian exlimit stefda veltman redbmk chongzhang-vis tanmancan mapoor andriimz pedro antmeehan fxfactorial nusken twista auserj qodbtn41 lp249839965 warmleaf jolars jwkvam zbennett10 kyuhwas sbrk buckaroo-pm ran125 sada-sam cartoecho projectcheshire rossrogers opencai ychuan1115 anbinhou ying-tch rubybrowncoat donilz frogcat chiehrosswang barrycug divisionof kcimit skymg frederic-schwarz andersjensen60 habx jutaz robinpowered az09 nerik tronyel panjarda yansj2016 changzheng1017 parkingboss anujsinghwd paulweeeel elocke codeyu001 ohjack stonetip datavis-tech felixyr lai3d fzbrill phlicess kalabasa lyrillind alexswan michaelowendyer erisonliang hollasch linecode kanyu tunchz tuanna59 hoorang xzw168

polylabel's Issues

Not getting the rectangle center when passing a rectangle as polygon

I'm using this to add labels to random shape polygons, mainly because I can sometimes get "L" shaped polygons and using the barycenter would give me points outside of the polygon. It works pretty good on those cases, but on very simple cases, I'm getting bad results.

On rectangular shapes, I most of the time get a point close to one of the 4 sides. I would expect a point in the dead center in those cases.

Coords wrong with complex L shape

Seems the coords generated from polylabel are wrong for this L shape

https://stackblitz.com/edit/mapbox-simple-map-q8vjbb?file=index.ts

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?

New update breaks Webpack build

Hello! While updating the dependencies of one of our apps we just noticed that polylabel got an update. Sadly this update breaks the module for our frontends, because with tinyqueues export syntax the imported module is no longer the class constructor, which instead lies in require('tinyqueue').default. This is the error message we are getting:

Uncaught TypeError: Queue is not a constructor

We've fixed it for now by pinning tinyqueue to version 1.2.2, but this isn't a great solution. Any tips on how we could fix this? Or is there potentially any way for you to fix this on your side? (e. g. by checking if a default property exists on the import)

Thank you!

Next release for node module

Hi there.

I'm curious when you're planning the next release for the node module? Currently the module won't work due to a problem with tinyqueue. It looks like it's been fixed in the most recent commit on master at L29 of polylabel.js.

If you aren't planning a release, can you create a stable branch that can be used.

Cheers!

[NaN, NaN] as a result of passing "Spain" geometry from TopoJSON

My task is to calculate center of each given continent/country/state based on topo.json coordinates I have for leaflet.
When I'm calculating this for Canada, US, GB, Sweden or even Italy - everything works great, but when I'm passing "Spain's" coordinates, polylabel returns value of [NaN, NaN].

In which cases polylabel can return this value and what is causing such behavior?

Only thing I've spotted here, is that first and last point is the same. Removing it prior to calculation didn't gave any results.

UPD: This also happens for a few other countries. Looking simply at coordinates is not making any sense for me, as I'm not familiar with polylabel's algorithm and principles yet.

Here is the data I'm trying to calculate center coordinates for:

[[-1.7945079745609398,43.38543254241482],[-1.7485686907125277,43.291246520623524],[-1.6394628915725526,43.25514187893686],[-1.559069144837835,43.284967452504105],[-1.4040240618494515,43.24258374269802],[-1.435607319495233,43.03537449475717],[-1.2748198260257944,43.04636286396616],[-0.9819568914921781,42.95688614326443],[-0.760874087971704,42.945897774055446],[-0.7292908303259189,42.876828024741826],[-0.5627609263754323,42.77322340077141],[-0.3072236599686491,42.83130478087604],[-0.16366339794236495,42.77793270186097],[-0.04020157259976287,42.68531644709953],[0.16939640995861183,42.726130389875756],[0.27563100385806294,42.668049009771124],[0.6431452746453452,42.683746680069675],[0.6575013008479758,42.83758384899546],[0.8125463838363594,42.83130478087604],[0.9216521829763344,42.78421176998039],[1.0680836502431426,42.77636293483111],[1.1341213707752367,42.716711787696624],[1.3178785061688778,42.71357225363692],[1.4298555105493769,42.5958397263978],[1.4068858686251708,42.48595603430796],[1.4470827419925278,42.43415372232275],[1.639453493107748,42.4671188299497],[1.7054912136398421,42.50322347163636],[1.9294452224008403,42.43729325638245],[1.9638996852871493,42.366653740038984],[2.0529070477434423,42.35252583677029],[2.1218159735160604,42.41060721687492],[2.2768610565044476,42.42787465420333],[2.4634893971386127,42.33996770053145],[2.641504122051206,42.3729328081584],[2.934367056584822,42.47025836400941],[3.1812907072700263,42.431014188263035],[3.1726770915484543,42.36351420597927],[3.3191085588152625,42.322700263203046],[3.2616844540047474,42.23950261062074],[3.1525786548647723,42.256770047949146],[3.1209953972189943,42.14060728773988],[3.2128739649158113,42.06839800436656],[3.2358436068400174,41.96636314742599],[3.1927755282321293,41.88787479593324],[3.069313702889527,41.79996784226137],[2.779321973596435,41.646130673335584],[2.4491333709359857,41.53467721421589],[2.2969594931881296,41.46560746490227],[2.1304295892376395,41.30392146082722],[2.0615206634650214,41.274095887259975],[1.538961309689352,41.187758700617955],[1.2173863227504746,41.10456104803564],[1.1858030651046967,41.06217733822956],[0.9905611087489525,41.04020059981159],[0.7149254056584873,40.815723914542346],[0.8728416938873984,40.70740998948236],[0.7838343314311054,40.66188674561656],[0.6000771960374607,40.613223967691056],[0.5483955017080007,40.57241002491483],[0.350282340111729,40.285142658451385],[0.20672207808544485,40.172119432301834],[0.14355556279388182,40.07950317754039],[0.06603302129968824,40.03084039961489],[-0.19524665558814647,39.715317226614054],[-0.3215796861712761,39.515956813822484],[-0.3014812494875976,39.321305702120476],[-0.21534509227182497,39.18159643646339],[-0.209602681790777,39.07799181249297],[-0.12059531933448042,38.950840683074716],[-0.008618314953977801,38.869212797522266],[0.1521691785154573,38.83153838880575],[0.2354341304907024,38.757759338402565],[-0.014360725435029309,38.6243291408649],[-0.0890120616886989,38.53485242016317],[-0.33019330189285157,38.47049197193912],[-0.41345825386809665,38.365317580938836],[-0.5196928477675478,38.299387365684936],[-0.5053368215649172,38.20834087795335],[-0.6087002102238408,38.17851530438611],[-0.6546394940722529,37.99014326080352],[-0.7178060093638159,37.94305024990787],[-0.8527526556685245,37.72328286572819],[-0.7034499831611889,37.61967824175776],[-0.9245327866816666,37.55531779353371],[-1.0680930487079472,37.58200383304124],[-1.1657140268858193,37.54432942432473],[-1.3236303151147304,37.563166628682986],[-1.470061782381542,37.483108510160385],[-1.487289013824693,37.43287596520503],[-1.6710461492183377,37.363806215891415],[-1.8117352060040943,37.20682951290592],[-1.8978713632198634,36.94938772000972],[-2.1275677824619166,36.73589940394945],[-2.20509032395611,36.73746917097931],[-2.3630066121850213,36.83636449386017],[-2.57547579998392,36.811248221382485],[-2.658740751959165,36.69508546117322],[-2.756361730137037,36.68095755790453],[-2.9315052498091028,36.751597074248],[-2.9860581493790903,36.73746917097931],[-3.227239389583243,36.75473660830771],[-3.4138677302174116,36.69508546117322],[-3.6004960708515803,36.74217847206887],[-3.729700306675234,36.724911034740465],[-3.838806105815209,36.751597074248],[-3.9622679311578146,36.72648080177032],[-4.057017704095159,36.746887773158434],[-4.226418813286173,36.71235289850163],[-4.4331455906040205,36.71078313147177],[-4.525024158300845,36.580492467993814],[-4.723137319897113,36.492585514321945],[-4.90402325005023,36.50671341759063],[-4.99590181774705,36.45648087263528],[-5.185401363621747,36.41095762876949],[-5.334704036129079,36.18805071053009],[-5.337575241369606,36.14095769963444],[-5.357673678053285,36.14095769963444],[-5.438067424788002,36.161364671022554],[-5.446681040509581,36.057760047052135],[-5.610339739219544,36.00752750209678],[-5.796968079853709,36.08601585358952],[-5.917558699955787,36.18491117647038],[-6.038149320057865,36.18962047755994],[-6.158739940159943,36.304213470739356],[-6.2190352502109825,36.40310879362021],[-6.2333912764136095,36.57578316690425],[-6.394178769883048,36.637004081068596],[-6.431504438009881,36.751597074248],[-6.359724306996739,36.78613194890481],[-6.500413363782499,36.9603760892187],[-6.744465809227179,37.101655121905644],[-7.02871512803922,37.20525974587606],[-7.413456630269657,37.192701609637226],[-7.4766231455612235,37.494096879369366],[-7.5283048398906836,37.56630616274269],[-7.42494145123176,37.74996890523572],[-7.335934088775463,37.81275958642991],[-7.269896368243376,37.94932931802729],[-7.123464900976565,38.04037580575888],[-7.022972717558169,38.021538601400614],[-6.948321381304499,38.197352508744366],[-7.126336106217092,38.19107344062495],[-7.166532979584449,38.2727013261774],[-7.358903730699669,38.44537569946144],[-7.269896368243376,38.73735236701445],[-7.05455597520395,38.855084894253565],[-7.066040796166053,38.90374767217907],[-6.974162228469233,39.01363136426892],[-7.037328743760796,39.10938715309006],[-7.166532979584449,39.11409645417963],[-7.244055521078643,39.19572433973208],[-7.249797931559694,39.27107315716512],[-7.327320473053888,39.34171267350859],[-7.312964446851261,39.45787543371785],[-7.39910060406703,39.493980075404515],[-7.5168200189285805,39.59287539828537],[-7.557016892295941,39.67921258492739],[-7.341676499256515,39.671363749778116],[-7.034457538520272,39.69020095413638],[-7.002874280874487,39.80322418028593],[-6.925351739380297,39.87229392959955],[-6.8794124555318845,40.008863661196926],[-7.043071154241847,40.181538034480965],[-6.896639686975039,40.25531708488414],[-6.796147503556639,40.355782174794854],[-6.8564428136076785,40.44211936143687],[-6.813374734999794,40.64775884234787],[-6.836344376924,40.84084018702002],[-6.942578970823451,41.01665409436377],[-6.8191171454808455,41.054328503080285],[-6.646844831049307,41.267816819140556],[-6.563579879074062,41.27095635320026],[-6.319527433629382,41.420084221036475],[-6.2046792240083555,41.57078185590255],[-6.365466717477791,41.66339811066399],[-6.554966263352487,41.67438647987297],[-6.577935905276689,41.73874692809702],[-6.523383005706702,41.86746782454513],[-6.655458446770883,41.93339803979903],[-7.051684769963426,41.941246874948305],[-7.143563337660247,41.988339885843956],[-7.249797931559694,41.86432829048542],[-7.588600149941723,41.82508411473904],[-7.723546796246428,41.898863165142224],[-8.047992988425829,41.81723527958977],[-8.18006842949001,41.81095621147035],[-8.23175012381947,41.88630502890339],[-8.108288298476868,42.01031662426192],[-8.191553250452113,42.06211893624713],[-8.220265302857367,42.15316542397872],[-8.346598333440497,42.10136311199351],[-8.625105241771486,42.05113056703815],[-8.751438272354616,41.96950268148569],[-8.877771302937745,41.9067120002915],[-8.897869739621424,42.10293287902336],[-8.843316840051436,42.1641537931877],[-8.757180682835667,42.40746768281521],[-8.866286481975642,42.415316517964484],[-8.823218403367758,42.56915368689027],[-8.903612150102475,42.62566529996504],[-8.969649870634566,42.5785722890694],[-9.070142054052965,42.5958397263978],[-9.029945180685605,42.7072931855175],[-9.11321013266085,42.76537456562213],[-9.141922185066104,42.92235126860762],[-9.225187137041349,42.93961870593603],[-9.251027984206079,43.016537290398915],[-9.216573521319773,43.15467678902615],[-9.05291482260981,43.19706049883223],[-8.969649870634566,43.286537219533955],[-8.843316840051436,43.33833953151917],[-8.6940141675441,43.29281628765338],[-8.507385826909935,43.33049069636989],[-8.444219311618369,43.37601394023569],[-8.214522892376316,43.42467671816119],[-8.33224230723787,43.47176972905684],[-8.309272665313664,43.54554877946002],[-8.085318656552662,43.657002238579715],[-8.059477809387932,43.70880455056493],[-7.775228490575891,43.743339425221734],[-7.542660866093314,43.7307812889829],[-7.3675173464212484,43.68054874402754],[-7.229699494876016,43.56909528490784],[-6.770306656391909,43.56909528490784],[-6.649716036289831,43.58008365411683],[-6.509026979504075,43.54711854648987],[-6.22190645545151,43.59578132441537],[-6.069732577703647,43.56909528490784],[-5.906073878993684,43.587932489266095],[-5.854392184664224,43.66328130669913],[-5.707960717397416,43.56909528490784],[-5.423711398585375,43.56124644975856],[-5.179658953140695,43.48903716638524],[-4.837985529518139,43.44979299063887],[-4.734622140859216,43.41682788301191],[-4.4331455906040205,43.40269997974322],[-4.18909314515934,43.40583951380293],[-3.827321284853106,43.49688600153451],[-3.718215485713131,43.48275809826582],[-3.5545567870031682,43.51729297292263],[-3.3908980882932056,43.42467671816119],[-3.11239117996222,43.35403720181772],[-2.9286340445685752,43.43723485440003],[-2.756361730137037,43.44979299063887],[-2.4835972322870994,43.371304639146125],[-2.368749022666073,43.31479302607134],[-2.098855730056659,43.30380465686236],[-1.7945079745609398,43.38543254241482]]

Not working with odd shaped polygon

The point generated falls outside of this polygon:

{
    "type": "Polygon",
    "coordinates": [
        [
            [-93.02957807480158, 44.861837817995564],
            [-93.02959954453688, 44.86892493940014],
            [-93.0280974715503, 44.86877290567057],
            [-93.02829060933638, 44.869153052906924],
            [-93.0298784714326, 44.86912268451289],
            [-93.02992141090319, 44.86551834530052],
            [-93.03996358763281, 44.86551834530052],
            [-93.04249557909264, 44.86548797498883],
            [-93.04513491922891, 44.86503171802141],
            [-93.04513491922891, 44.86354124914875],
            [-93.04599316965135, 44.86284162150566],
            [-93.0443409881808, 44.86136631401641],
            [-93.0433968588172, 44.859586789464856],
            [-93.0427531260846, 44.857852803290726],
            [-93.0422595916675, 44.856818474236],
            [-93.0414441908841, 44.856240456857805],
            [-93.04043565231463, 44.85555596011782],
            [-93.0401567254189, 44.85552558455064],
            [-93.0401567254189, 44.85596669734617],
            [-93.04183037662474, 44.85695538396923],
            [-93.04234538077714, 44.857837584275984],
            [-93.04294617403916, 44.859678037043125],
            [-93.04382598402843, 44.86150321293619],
            [-93.04513491922891, 44.862735161220904],
            [-93.0445126562316, 44.86326746067597],
            [-93.04457697560595, 44.86483395885438],
            [-93.0424312597183, 44.86522941283937],
            [-93.03916965658472, 44.865275064010014],
            [-93.03000720001282, 44.865244629900275],
            [-93.02989994116788, 44.86189862612768],
            [-93.02957807480158, 44.861837817995564]
        ]
    ]
}

The generated point appears to be these coordinates:

{"x":-93.03840446234251,"y":44.863239663088805,"spatialReference":{"wkid":4326}}

The dark solid gray point is the generated point:

import polylabel from 'polylabel';
import {arcgisToGeoJSON, geojsonToArcGIS} from '@esri/arcgis-to-geojson-utils';
import {project} from 'esri/geometry/support/webMercatorUtils';
import {WGS84, WebMercator} from 'esri/geometry/SpatialReference';
import Point from 'esri/geometry/Point';

export default function(graphic, event, vm) {
    let polygon = project(graphic.geometry, WGS84);
    polygon = arcgisToGeoJSON(polygon.toJSON());
    const point = geojsonToArcGIS({
        type: 'Point',
        coordinates: polylabel(polygon.coordinates),
    });
    return new Point(project(point, WebMercator));
}

return best distance

it would be super useful in many cases to know what the best distance was, for example to decide whether or not to display a label, or to use in other computations.

Non-breaking change 1:
var r = [bestCell.x, bestCell.y]; r.distance = bestCell.d; return r;

Not-too-much-breaking change 2:
return [bestCell.x, bestCell.y, bestCell.d];

1 is probably nicer but some developers might find the array/object mix awkward to work with.

Result point is out of polygon

I'm trying to find the label point for Lago di Livigno lake in Italy. The result point is out of polygon bounds. See the source polygon and the picture. The coordinates of the result point are:

10.179204834796153,
46.593747060598524

livigno.geojson.txt

Support for MultiPolygons

I'm not sure if this is related to #25 or not, but I couldn't get a working solution from that issue.

From what I can tell this will always return [NaN, NaN] for a MultiPolygon. I have a few disconnected polygons, and for each one I can do polylabel(feature.geometry.coordinates) to get a point for each. However, if I want the pole of inaccessibility for all the features combined, the following both just give me [NaN, NaN]:

polylabel(turf.union(features).geometry.coordinates)
polylabel(features.map(feature => feature.geometry.coordinates))

Is this the intended behavior? Currently when I run into this situation I use polylabel for each one, then figure out which one is closest to the centroid of the unioned MultiPolygon. Maybe other option would be to get the pole of inaccessibility of the largest Polygon.

Python wrapper bug

Opening an issue to inform that polylabel() from the linked Python implementation of this library fails when polygons are sufficiently small. (This JavaScript implementation handles this edge case already). I have patched the issue locally and would like to contribute the code, but since the library is from a third-party and this official repository references it, I'm mentioning the issue here.

regions of inaccessibility

I'm sure this is outside the scope of this module, but I was wondering your thoughts if you wanted to weigh in.

I am interested in finding regions of inaccessibility, which I would define as all points that have a distance to the polygon close to that of the pole of inaccessibility's distance to the polygon. I think I could adapt this algorithm so that once the pole of inaccessibility is found, I could visit the cells again and prune cells which are guaranteed to not have any regions of inaccessibility, keep cells that are wholly regions of inaccessible and split cells that are in-between. Something akin to alpha beta pruning, IIRC.

For my application, points that are almost poles can be reasonable solutions. For example, like putting text on a polygon :)

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

Not an issue - I ported polylabel to C# to work with ESRI Polygons.

Not an issue - I ported polylabel to C# to work with ESRI Polygons and was wondering if someone could lay an eye on my check-in (Epi-Info/Epi-Info-Community-Edition@d7ef1a0) and clue me in where I may have gone wrong by inspection before I have to go into a major debug. Thank you. David

...best distance: -0.00777946996556627
found best -0.0478 after 10 probes
found best -0.0393 after 18 probes
num probes: 18
best distance: -0.0393367048984066
found best -0.1129 after 2 probes
found best -0.0921 after 6 probes
found best -0.0107 after 10 probes
num probes: 10
best distance: -0.0107072800214816
found best -0.1414 after 14 probes
found best -0.1232 after 18 probes
found best -0.0968 after 22 probes
found best -0.0015 after 22 probes
num probes: 22
best distance: -0.00151480641263599
found best -0.0529 after 18 probes
num probes: 18
best distance: -0.0528912840084305

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.

Javascript Usage

Hi there,

I might be too stupid to use it in javascript. Just wondering how I should use it in my javascript/html so that I can avoid the error information of 'request is not defined' or 'Queue is not a constructor'?

Regards,

polylabel as a command line utility

@mourner I wanted a command line utility to wrap up this library so I put together one at https://github.com/andrewharvey/geojson-polygon-labels

Not sure if you're interested in including a utility within the polylabel repo or a link in the README, or if not that's fine. Anyway just thought I would let you know.

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.

How to import using ES6 syntax?

Having issues wrapping my head around importing this module using a ES6 syntax.

Can you elaborate on how to import this module?

Which one of these would be ES6 valid and still functional?

const polylabel = require('polylabel')
import polylabel from 'polylabel'
import { polylabel } from 'polylabel'
import * as polylabel from 'polylabel'

Units for precision

What units are the precision in? The same as the input polygon? e.g. if I use a polygon in decimal degrees with precision 1.0 then I will get the pole within 1 decimal degree?

Algorithm returns external point

The point returned for the polygon
var polygon = [[[12.8569046,43.7153127],[12.85684,43.7153125],[12.8568395,43.7152068],[12.8568758,43.7152072],[12.8568766,43.7151763],[12.8570085,43.7151777],[12.8570083,43.7152309],[12.8569055,43.7152302],[12.8569046,43.7153127]]];
is not a Point inside the polygon itself.
var point = polylabel(polygon,1,true);
debug output:
found best 0 after 2 probes num probes: 2 best distance: -0.0000023558598186856006
sample image with the returned point outside the given polygon:

algorithm fails on polygons around 180 longitude

i'm surprised that through all this years nobody found this.
e.g. something like (lat, lon)
66.835993985609136,179.53481017154348,
66.570003009644324,-179.76022250492355,
66.751380624020030,-179.47256342463140,
67.039144593094719,179.70497469791349,
will give you wrong bounding box (like polygon goes along all earth surface)
66.570003009644324,-179.76022250492355,
67.039144593094719,179.70497469791349,
and of course wrong center
66.937436163870117,179.45104221247121,

Use ES-modules instead of CommonJS or AMD

Compiling with Angular CLI version greater than 10 results in a warning by the Angular compiler. This is caused by the usage of CommonJS or AMD dependencies within this package or its dependencies.

Can you provide a version of your module only using ES-modules? I'm currently using version 1.1.0 of polylabel.

Move queue instantiation after return condition

I suggest moving this line after the "if( cellSize === 0 ) return;". Won't make a big difference in most cases but I think it's a good practice not to create an instance that could never be used.

polylabel/polylabel.js

Line 27 in 51f09d0

var cellQueue = new Queue(null, compareMax);

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

Two small optimizations for getCentroidCell()

There are two opportunities for improving the performance of the centroid calculation, saving one multiply and one addition per polygon vertex.

First, when summing the vertex coordinates, we scale the sum of coordinates of the prior and the current vertex. But since we're just looping over all vertices, that just means that we'll be adding in a coordinate twice (once as the prior, and once as the current). Thus, we can just add each coordinate once, and scale by two outside the loop, saving an addition for every vertex.

Second, when accumulating the area, we scale f by three before adding. Instead of multiplying N times inside the loop, we can just scale the sum of fs by three after the loop is done, saving a multiply for every vertex.

Here's my proposed iteration:

function getCentroidCell(polygon) {
    var doubleArea = 0;
    var xSum = 0;
    var ySum = 0;
    var points = polygon[0];

    for (var i = 0, len = points.length, j = len - 1;  i < len;  j = i++) {
        var a = points[j]; // Swapped indices. Might as well traverse in native order.
        var b = points[i];
        var f = a[0] * b[1] - a[1] * b[0];

        xSum += f * a[0];
        ySum += f * a[1];
        doubleArea += f;
    }

    if (doubleArea === 0) return new Cell(points[0][0], points[0][1], 0, polygon);

    var sumScale = 2 / (3 * doubleArea);
    return new Cell(xSum * sumScale, ySum * sumScale, 0, polygon);
}

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

Using lat/lng coordinates in lieu of GeoJSON format

Hi,
Your algorithm looks great :-), but was wondering if I could use it using traditional google maps lat/lng?
Perhaps converting lat/lng to points.. would that work?
https://developers.google.com/maps/documentation/javascript/reference/image-overlay#Projection.fromLatLngToPoint

Thanks

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.

Label positioning is wrong on rooms that aren't perfectly square?

@mourner if you remember I previously posted having issues with L shaped rooms (#55) where you told me to decrease the tolerance to 0.00001 . But now I noticed another issue with rooms that aren't perfectly square:

https://stackblitz.com/edit/mapbox-simple-map-cbnty5?file=index.ts

not sure if this helps but i also posted this when i zoomed into the points:

MultiPolygon

(enhancement)

As the GeoJSON format is concerned, I figure it out that "type":"MultiPolygon" is not supported. I suppose this would take an other kind of algorithm (idk?). Anyway, thanks for sharing !

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

[enhancement] weight x, y differently?

Labels are usually nowhere near square; they are often an elongated bit of text. So the desired weighting (in a loose sense) of x-distances should be different from y-distances.

I'd suggest being able to pass in an aspect ratio. The algorithm modifications may be as simple as just weighting distance measurements.

polylabel hanging with degenerate case

A user added a polygon to my system that is a degenerate case. This causes function polylabel to spin forever. The polygon has almost no height at all, causing a minuscule height on this line:

var cellSize = Math.min(width, height);

height is 4.97e-14 and width is 0.0017.

This causes the following for loop to spin for a very long time:

for (var x = minX; x < maxX; x += cellSize) {
[...]

Precision is set to 1.0. Shouldn't cellSize be bounded in minimum size by the precision somehow?

To get top N 'Best point' in the polygon and have a min-distance between each other.

Hi
i have this question when i trying to place more label in the poly.how to To get top N 'Best point' in the polygon and have a min-distance between each other.have any suggestion?

regards

Handling winded polygons with modified poin-in-polygon test

In my case of applying polylabel code, I've found some strange spaces on floor plan imported from some external CAD software, which occurred to wind themselves around its boundary. I detected this by splitting duplicated vertices apart, as below.

Original algorithm regarded all points inside visible contour as laying "outside", and "interior" point has been calculated as somewhere on duplicated boundary edge.
To cover such case, I had to modify test if point lies inside polygon by checking its winding number algorithm.
I don't want to cover it in formal pull request, so below is the part of modified code with core taken from well-known geometry algorithms:

// signed distance from point to polygon outline (negative if point is outside)
// redesigned using winding number algorithm for point inside polygon test
// http://geomalgorithms.com/a03-_inclusion.html#wn_PnPoly()

function pointToPolygonDist(x, y, polygon) {
    var wn = 0;
    var minDistSq = Infinity;

    for (var k = 0; k < 1; k++) {
        var ring = polygon[k];

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

            var xi = a[0], yi = a[1];
            var xj = b[0], yj = b[1];
            if (yj <= y) {
                if (yi > y) {
                    if (isLeft([xj, yj], [xi, yi], [x,y]) > 0) {
                        wn++;
                    }
                }
            } else {
                if (yi <= y) {
                    if (isLeft([xj, yj], [xi, yi], [x, y]) < 0) {
                        wn--;
                    }
                }
            }
            minDistSq = Math.min(minDistSq, getSegDistSq(x, y, a, b));
        }
    }

    return minDistSq === 0 ? 0 : (wn != 0 ? 1 : -1) * Math.sqrt(minDistSq);
}

function isLeft(P0, P1, P2)
{
    var res = ( (P1[0] - P0[0]) * (P2[1] - P0[1])
            - (P2[0] -  P0[0]) * (P1[1] - P0[1]) );
    return res;
}

Handling multi-ring polygons

I don't think this handles - or intends to handle - polygons with multiple rings. Typically, polygon rings with counter-clockwise winding are 'solid' and clockwise windings represent 'holes'.

Any thoughts on the level of effort to extend this library to account for polygons with one or more holes?

Typings not working.

I have tried to work with the typings you provided but its not working.

So I tried to put some modificationsand I ended up with the following code:

declare namespace polylabel {
/**
 * Polylabel returns the pole of inaccessibility coordinate in [x, y] format.
 * 
 * @name polylabel
 * @function
 * @param {Array<number>} polygon - Given polygon coordinates in GeoJSON-like format
 * @param {number} precision - Precision (1.0 by default)
 * @param {boolean} debug - Debugging for Console 
 * @return {Array<number>}
 * @example
 * var p = polylabel(polygon, 1.0);
 */
function polylabel(polygon: number[][][], precision?: number, debug?: boolean): number[];
}

declare module "polylabel" {
	export = polylabel.polylabel;
}

Still I'm getting the following error:

cannot get a expression whose type lacks a call signature

how to use in browser?

With the source using require I'm guessing this is supposed to be an npm component? Can it be run in the browser? I went and found tinyqueue and added it to my test html page (along with a few edits to remove request and module, etc), then ran this:

var polygon = [[30, 10], [40, 40], [20, 40], [10, 20], [30, 10]];
var p = polylabel(polygon, 1.0);

But the result is [NaN, NaN], with no errors at all. Note I also tried a non-closed polygon and got the same result.

Any suggestions?

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

Accuracy depends on binary heap implementation

Testing a non-degenerate Polygon:

[
    (0.0, 0.0),
    (4.0, 0.0),
    (4.0, 1.0),
    (1.0, 1.0),
    (1.0, 4.0),
    (0.0, 4.0),
    (0.0, 0.0)
];

Gives different values when tested against this implementation ((0.5, 0.5)), than against the implementation in Shapely 1.6b2 ((0.5, 1.5)), and my Rust implementation ((0.5, 1.5)). This is due to the fact that the binary heap implementations return elements with equal values in a different order.

This version, using tinyqueue:

Popped cell with max distance: 1.8284271247461903, 2, 2
Splitting
Popped cell with max distance: 1.4142135623730951, 1, 1
Found greater distance: 0, -0.3571428571428572
New best coords: (1, 1)
Splitting
Popped cell with max distance: 1.4142135623730951, 3, 1
Splitting
Popped cell with max distance: 1.4142135623730951, 1, 3
Splitting
Popped cell with max distance: 1.2071067811865475, 0.5, 0.5
Found greater distance: 0.5, 0
New best coords: (0.5, 0.5)
Popped cell with max distance: 1.2071067811865475, 3.5, 0.5
Popped cell with max distance: 1.2071067811865475, 1.5, 0.5
Popped cell with max distance: 1.2071067811865475, 0.5, 2.5
Popped cell with max distance: 1.2071067811865475, 0.5, 1.5
Popped cell with max distance: 1.2071067811865475, 0.5, 3.5
Popped cell with max distance: 1.2071067811865475, 2.5, 0.5
Popped cell with max distance: 0.20710678118654757, 2.5, 1.5
Popped cell with max distance: 0.20710678118654757, 1.5, 2.5
Popped cell with max distance: 0.20710678118654757, 1.5, 1.5
Popped cell with max distance: 0.20710678118654757, 3.5, 1.5
Popped cell with max distance: 0.20710678118654757, 1.5, 3.5
Popped cell with max distance: -0.5857864376269049, 3, 3

The Rust version:

Popped cell with max distance: (1.8284271247461903, 2, 2)
Splitting
Popped cell with max distance: (1.4142135623730951, 1, 1)
Found greater distance: (-0, -0.3571428571428572)
New best coords: (1, 1)
Splitting
Popped cell with max distance: (1.4142135623730951, 1, 3)
Splitting
Popped cell with max distance: (1.4142135623730951, 3, 1)
Splitting
Popped cell with max distance: (1.2071067811865475, 0.5, 1.5)
Found greater distance: (0.5, -0)
New best coords: (0.5, 1.5)
Popped cell with max distance: (1.2071067811865475, 0.5, 2.5)
Popped cell with max distance: (1.2071067811865475, 1.5, 0.5)
Popped cell with max distance: (1.2071067811865475, 2.5, 0.5)
Popped cell with max distance: (1.2071067811865475, 3.5, 0.5)
Popped cell with max distance: (1.2071067811865475, 0.5, 0.5)
Popped cell with max distance: (1.2071067811865475, 0.5, 3.5)
Popped cell with max distance: (0.20710678118654757, 1.5, 3.5)
Popped cell with max distance: (0.20710678118654757, 2.5, 1.5)
Popped cell with max distance: (0.20710678118654757, 3.5, 1.5)
Popped cell with max distance: (0.20710678118654757, 1.5, 1.5)
Popped cell with max distance: (0.20710678118654757, 1.5, 2.5)
Popped cell with max distance: (-0.5857864376269049, 3, 3)

I can't test the cpp version here, but it looks like we need a way to break ties between equal max distance values, because different binary heap implementations will internally order them differently…

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,

redundant calculation

In https://github.com/mapbox/polylabel/blob/master/index.js#L95, we iterate each segments to calculate distance from point to segment.

Support we have a polygon with 4 nodes as fellows. We should see that point 0 and 4 are identical. So I think the calculation of point to segment 4-0 is unnecessary.

However, for compatible reason, it is appropriate to keep this redundant calculation, in case of non-closed polygons. Am I correct?

Consider also distance from centroid when calculating label position

In cases where there are several equally good solutions, algorithm's output is not always the most elegant one. The simplest case is a rectangle: the set of solutions lies on a line and algorithm returns a point that coincides with one of the line's endpoints. Current implementation covers such case when calculating initial best cell, however these lines can occur also in more complex polygons containing parallel edges.

I propose using a cost function when estimating cell quality, which would consider both distance from polygon and distance from centroid in such a way it would prefer points closer to centroid, e.g.:
0.5 * (distFromPoly + distFromPoly / (distFromCtr/distFromPoly+1); distFromPoly >= 0
distFromPoly; distFromPoly < 0

So for example, if pia point is very far away from centroid and its distance from polygon is not much greater than centroid's distance, we would get centroid as optimal point. Weights on both values can be manipulated by cost function. What are optimal set of weights is of course subjective.

Cost function would be then used when comparing two cells in while loop and max attribute of a cell can be then calculated in a following way (java code):
// for distance from polygon take max possible distance
// for distance from mass centre take minimal possible distance
double maxDist = dist + h*SQRT2;
if (maxDist < 0) {
return maxDist;
}
double distFromCtr = Math.max(distance(this, massCentreCell) - h*SQRT2, 0);
return costFunction.compute(distFromCtr, maxDist);

I've implemented this in Java and I'm very satisfied with the results.

Special case for rectangular polygons comment is misleading

After computing the polygon centroid to see an initial guess, polylabel computes an additional special case, with the comment “special case for rectangular polygons”. The comment puzzled me, because for rectangular polygons, this would simply compute the centroid through a different method, yielding a redundant candidate point.

I then realized that this computes the bounding-box center point as a legitimate candidate for any polygon. There are indeed cases where this would yield a better solution than a centroid. I'd suggest replacing this comment with something more accurate, like “test the bounding-box center as a possible better candidate”.

suboptimal centroid

I'm getting suboptimal centroids for non-skewed rectangular polygons (horizontal and vertical).

Coordinates:
32.71997,-117.19310
32.71997,-117.21065
32.72408,-117.21065
32.72408,-117.19310

Posts about polylabel in Chinese

Hi, @mourner I have post two articles to describe your polylabel. Should I append it in README?
http://jingsam.github.io/2016/09/24/polylabel.html
http://jingsam.github.io/2016/09/26/polylabel2.html

Multipolygon/MultiLineString/MultiPoint support?

Any interest in supporting MultiX geojson, or maybe a workaround, like should I be doing it for each child and then averaging it all together?

IOS SDK Version

Hey- this is fantastic! Thanks so much for developing it. Any chance it's going to be integrated (or already is) with the Mapbox IOS SDK?

Make the algorithm faster with a spatial index

Currently the algorithm is relatively slow because checking each probe is expensive — you have to loop through all polygon points to do two things:

Determine if the probe is inside or outside the polygon.
Find the distance from the probe to the polygon.

In theory, we could make this much faster by indexing polygon segments with RBush (or any fast RTree index in other impementations), and then do the above with:

A zero-height bbox query that loops through all intersecting edges (very fast for all practical polygons) to determine if inside or outside.
A kNN-like tree search to find the closest edge for point-to-polygon calculation. Can be borrowed from Concaveman.

The current algorithm is O(NM), where M is the number of probes.

The new one would be O(N log N) for indexing + O(M log N) for checking probes, or O((M + N) log N in average, which should be MUCH faster.

cc @urschrei @chau-intl