Mapbox Visual Center Algorithm: arcpy (attempt #2)

Yesterday, I attempted to replicate the Mapbox Visual Center algorithm using ArcGIS arcpy (Python). Today, I actually looked at the source code (code, license) and attempted to fully port it over from JavaScript. I think it works properly, but you tell me.

def polylabel(polygon, centroid, precision, debug):
    cells = []
    precision = precision or 1.0
    extent = polygon.extent
    minX, minY, maxX, maxY = (extent.XMin, extent.YMin, extent.XMax, extent.YMax)
    width = extent.width
    height = extent.height
    cellSize = min([width, height])
    h = cellSize / 2
    cellQueue = []
    x = minX
    while x < maxX:
        y = minY
        while y < maxY:             cell = Cell(x+h, y+h, h, polygon)             cells.append(cell.geom)             cellQueue.append(cell)             y += cellSize         x += cellSize     bestCell = Cell(centroid[0], centroid[1], 0, polygon)     numProbes = len(cellQueue)     while len(cellQueue):         cellQueue.sort(key=lambda cell: cell.d)         cell = cellQueue.pop()         cells.append(cell.geom)         if cell.d > bestCell.d:
            bestCell = cell
        if cell.max - bestCell.d <= precision:
            continue
        h = cell.h/2
        cellQueue.append(Cell(cell.x-h, cell.y-h, h, polygon))
        cellQueue.append(Cell(cell.x+h, cell.y-h, h, polygon))
        cellQueue.append(Cell(cell.x-h, cell.y+h, h, polygon))
        cellQueue.append(Cell(cell.x+h, cell.y+h, h, polygon))
        numProbes += 4
    return [bestCell.x, bestCell.y, cells]
class Cell():
    def __init__(self, x, y, h, polygon):
        self.x = x
        self.y = y
        self.h = h
        self.d = pointToPolygonDist(x, y, polygon)
        self.max = self.d + self.h * math.sqrt(2)
        self.geom = arcpy.Polygon(arcpy.Array([arcpy.Point(x-h,y-h),arcpy.Point(x+h,y-h),arcpy.Point(x+h,y+h),arcpy.Point(x-h,y+h)]))
def pointToPolygonDist(x, y, polygon):
    point_geom = arcpy.PointGeometry(arcpy.Point(x, y),sr)
    polygon_outline = polygon.boundary()
    min_dist = polygon_outline.queryPointAndDistance(point_geom)
    sign_dist = min_dist[2] * ((min_dist[3]-0.5)*2)
    return sign_dist
fc = 'BEC_POLY selection'
sr = arcpy.Describe(fc).spatialReference
centroids = []
label_points = []
cells = []
with arcpy.da.SearchCursor(fc,['SHAPE@','SHAPE@XY']) as cursor:
    for row in cursor:
        best_point = polylabel(row[0],row[1],.75,'#')
        centroids.append(arcpy.PointGeometry(row[0].centroid,sr))
        label_points.append(arcpy.PointGeometry(arcpy.Point(best_point[0], best_point[1]),sr))
        for cell in best_point[2]:
            cells.append(cell)
arcpy.CopyFeatures_management(centroids,r'in_memory\centroids')
arcpy.CopyFeatures_management(label_points, r'in_memory\label_points')
arcpy.CopyFeatures_management(cells, r'in_memory\cells')

Advertisements