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

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

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

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

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,

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

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

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

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?

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!

🙏

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.

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

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.

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

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.

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.

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.

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.