What is the fastest way to find the "visual" center of an irregularly shaped polygon?

I need to find a point that is a visual center of an irregularly shaped polygon. By visual center, I mean a point that appears to be in the center of a large area of the polygon visually. The application is to put a label inside the polygon.

Here is a solution that uses inside buffering:

https://web.archive.org/web/20150708063910/http://proceedings.esri.com/library/userconf/proc01/professional/papers/pap388/p388.htm

If this is to be used, what is an effective and fast way to find the buffer? If any other way is to be used, which is that way?

A good example of really tough polygons is a giant thick U (written in Arial Black or Impact or some such font).

I have found a very good solution to this from MapBox called Polylabel. The full source is available on their Github too.

Essentially it tries to find the visual centre of the polygon as T Austin said.

Certain details suggest this may be a practical solution: > Unfortunately, calculating [the ideal solution ] is both complex > and slow. The published solutions to the problem require either > Constrained Delaunay Triangulation or computing a straight skeleton as > preprocessing steps — both of which are slow and error-prone. > > For our use case, we don’t need an exact solution — we’re willing to > trade some precision to get more speed. When we’re placing a label on > a map, it’s more important for it to be computed in milliseconds than > to be mathematically perfect.

A quick note about usage though. The source code works great for Javascript out of the box however if you intend on using this with a "normal" polygon then you should wrap it in an empty array as the functions here take GeoJSONPolygons rather than normal polygons i.e.

var myPolygon = [[x1, y1], [x2, y2], [x3, y3]];
var center = polylabel([myPolygon]);

If you can convert the polygon into a binary image, then you can use the foundation that exists in the field of image processing, e.g.: A Fast Skeleton Algorithm on Block Represented Binary Images.

But this is not really reasonable in the general case, because of discretization errors and extra work.

However, maybe you find these useful:

• Straight skeleton of a simple polygon
• Determining the Skeleton of a Simple Polygon in (Almost) Linear Time

EDIT: Maybe you want to look for the point that is the center of the largest circle contained in the polygon. It is not necessarily always in the observed centre, but most of the time would probably give the expected result, and only in slightly pathological cases something that is totally off.

How about:

If the centroid of the polygon is inside the polygon then use it, else:

1. Extend a line from the centroid through the polygon dividing the polygon into two halves of equal area

2. The "visual center" is the point half way between the nearest point where the line touches the perimeter and the next point cutting the perimeter in the direction going away from the centroid

Here are a couple of pictures to illustrate it:

Have you looked into using the centroid formula?

http://en.wikipedia.org/wiki/Centroid

http://en.wikipedia.org/wiki/K-means_algorithm

Here are six different approaches I have tried.

1. cv2 based center of mass (get_center_of_mass)
2. shapely based representative point (get_representative_point)
3. cv2 + skimage.skeleton based center of mass of the skeletonized shape (get_skeleton_center_of_mass)
4. scipy based furthest distance to border (get_furthest_point_from_edge)
5. cv2 based version of the previous furthest distance to border -algorithm (get_furthest_point_from_edge_cv2)
6. the "center point of half-area line" algorithm proposed in this thread by @T.Austin (get_center_of_half_area_line)

Let's begin with imports and some helper functions

import numpy as np
import cv2
from shapely.geometry import Polygon, LineString, MultiLineString, Point, MultiPoint, GeometryCollection
from skimage.morphology import skeletonize, medial_axis
from scipy.optimize import minimize_scalar
from scipy.ndimage.morphology import distance_transform_edt
import matplotlib.pyplot as plt

H, W = 300, 300

def get_random_contour():
    xs = np.random.randint(0, W, 4)
    ys = np.random.randint(0, H, 4)
    cnt = np.array([[x,y] for x,y in zip(xs,ys)])
    mask = draw_contour_on_mask((H,W), cnt)
    cnt, _ = cv2.findContours(mask, 1, 2)
    cnt = cnt[0]
    return cnt

def draw_contour_on_mask(size, cnt, color:int = 255):
    mask = np.zeros(size, dtype='uint8')
    mask = cv2.drawContours(mask, [cnt], -1, color, -1)
    return mask

def get_center_of_mass(cnt):
    M = cv2.moments(cnt)
    cx = int(M['m10']/M['m00'])
    cy = int(M['m01']/M['m00'])
    return cx, cy

def split_mask_by_line(mask, centroid:tuple, theta_degrees:float, eps:float = 1e-4):
    h, w = mask.shape[:2]
    mask = mask[..., None]
    cx, cy = centroid
    # convert theta first to radians and then to line slope(s)
    theta_degrees = np.atleast_1d(theta_degrees)
    theta_degrees = np.clip(theta_degrees, -90+eps, 90-eps)
    theta_rads = np.radians(theta_degrees)
    slopes = np.tan(theta_rads)[:, None]
    # define the line(s)
    x = np.arange(w, dtype="int32")
    y = np.int32(slopes * (x - cx) + cy)
    _y = np.arange(h, dtype="int32")
    # split the input mask into two halves by line(s)
    m = (y[..., None] <= _y).T
    m1 = (m * mask).sum((0,1))
    m2 = ((1 - m) * mask).sum((0,1))
    m2 = m2 + eps if m2==0 else m2
    # calculate the resultant masks ratio
    ratio = m1/m2
    return (x.squeeze(), y.squeeze()), ratio

def get_half_area_line(mask, centroid: tuple, eps: float = 1e-4):
    # find the line that splits the input mask into two equal area halves
    minimize_fun = lambda theta: abs(1. - split_mask_by_line(mask, centroid, theta, eps=eps)[1].item())
    bounds = np.clip((-90, 90), -90 + eps, 90 - eps)
    res = minimize_scalar(minimize_fun, bounds=bounds, method='bounded')
    theta_min = res.x
    line, _ = split_mask_by_line(mask, centroid, theta_min)
    return line

Now let's define the functions for finding the visual center

def get_representative_point(cnt):
    poly = Polygon(cnt.squeeze())
    cx = poly.representative_point().x
    cy = poly.representative_point().y
    return cx, cy

def get_skeleton_center_of_mass(cnt):
    mask = draw_contour_on_mask((H,W), cnt)
    skel = medial_axis(mask//255).astype(np.uint8) #<- medial_axis wants binary masks with value 0 and 1
    skel_cnt,_ = cv2.findContours(skel,1,2)
    skel_cnt = skel_cnt[0]
    M = cv2.moments(skel_cnt) 
    if(M["m00"]==0): # this is a line
        cx = int(np.mean(skel_cnt[...,0]))
        cy = int(np.mean(skel_cnt[...,1]))
    else:
        cx = int(M['m10']/M['m00'])
        cy = int(M['m01']/M['m00'])
    return cx, cy

def get_furthest_point_from_edge(cnt):
    mask = draw_contour_on_mask((H,W), cnt)
    d = distance_transform_edt(mask)
    cy, cx = np.unravel_index(d.argmax(), d.shape)
    return cx, cy

def get_furthest_point_from_edge_cv2(cnt):
    mask = draw_contour_on_mask((H,W), cnt)
    dist_img = cv2.distanceTransform(mask, distanceType=cv2.DIST_L2, maskSize=5).astype(np.float32)
    cy, cx = np.where(dist_img==dist_img.max())
    cx, cy = cx.mean(), cy.mean()  # there are sometimes cases where there are multiple values returned for the visual center
    return cx, cy

def get_center_of_half_area_line(cnt):
    mask = draw_contour_on_mask((H,W), cnt, color=1)
    # get half-area line that passes through centroid
    cx, cy = get_center_of_mass(mask)
    line = get_half_area_line(mask, centroid=(cx, cy))
    line = LineString(np.array(list(zip(line))).T.reshape(-1, 2))
    # find the visual center
    contours, _ = cv2.findContours(mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
    contours = [c for c in contours if cv2.contourArea(c) > 5]
    polys = [Polygon(c.squeeze(1)) for c in contours if len(c) >= 3]  # `Polygon` must have at least 3 points
    cpoint = Point(cx, cy)
    points = []
    for poly in polys:
        isect = poly.intersection(line)
        if isect.is_empty:
            # skip when intersection is empty: this can happen for masks that consist of multiple disconnected parts
            continue
        if isinstance(isect, (MultiLineString, GeometryCollection)):
            # take the line segment intersecting with `poly` that is closest to the centroid point
            isect = isect.geoms[np.argmin([g.distance(cpoint) for g in isect.geoms])]
        if isinstance(isect, Point):
            # sometimes the intersection can be a singleton point
            points.append(isect)
            continue
        isect = isect.boundary
        if poly.intersects(cpoint):
            points = [isect]
            break
        else:
            points.append(isect)

    if len(points) == 0:
        # multiple reasons for this one:
        # - if len(polys)==0
        # - if len(polys)==1, but for some reason the line does not intersect with polygon
        # - if the above search does not match with any points

        return cx, cy

    points = points[np.argmin([p.distance(cpoint) for p in points])]
    if isinstance(points, Point):
        return np.array(points.xy)
    
    points = [np.array(p.xy).tolist() for p in points.geoms]
    visual_center = np.average(points, (0, 2))
    return visual_center

Here's my analysis on the topic:

• get_center_of_mass is the fastest but, as mentioned in this thread, the center of mass can be located outside the shape for non-convex shapes.
• get_representative_point is also fast but the identified point, although always guaranteed to stay inside the shape (or with minor edits even multiple disconnected shapes!), does not have much if anything to do with the center of the object
• get_skeleton_center_of_mass returns a perceptually nice center point, but is slow and requires logic for disconnected shapes
• get_furthest_point_from_edge is relatively fast, generalizes easily to disconnected shapes and the center point is visually pleasing
• get_furthest_point_from_edge_cv performs otherwise similarly as get_furthest_point_from_edge but is an order of magnitude faster
• get_center_of_half_area_line performs neatly: the result is usually closest to where I myself would annotate the visual center. Unfortunately, at least my implementation is quite slow.

rows = 4
cols = 4
markers = ['x', '+', "*", "o", '^', "v"]
colors = ['r','b','g','orange', 'purple', 'lime']
functions = [
    get_center_of_mass, 
    get_representative_point, 
    get_skeleton_center_of_mass, 
    get_furthest_point_from_edge,
    get_furthest_point_from_edge_cv2,
    get_center_of_half_area_line
]

plt.figure(figsize=(2*cols, 2*rows, ))
for i in range(rows*cols): 
    cnt = get_random_contour()
    mask = draw_contour_on_mask((H,W), cnt)
    
    plt.subplot(cols,rows, i+1)
    plt.imshow(mask, cmap='gray')
    for c, m, f in zip(colors, markers, functions):
        l = f.__name__
        cx, cy = f(cnt)
        plt.scatter(cx, cy, c=c, s=100, label=l, marker=m, alpha=0.7)

plt.tight_layout()    
plt.legend(loc=3)
plt.show()

Here's how the algorithms, run on 100 random examples, compare in speed:

N_EXAMPLES = 100
cnts = [get_random_contour() for _ in range(N_EXAMPLES)]

for fn in functions:
    print(fn.__name__+":")
    %time _ = [fn(cnt) for cnt in cnts]
    print("~ "*40)

get_center_of_mass:
CPU times: user 2.35 ms, sys: 777 µs, total: 3.13 ms
Wall time: 1.91 ms
~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 
get_representative_point:
CPU times: user 15.7 ms, sys: 0 ns, total: 15.7 ms
Wall time: 14.8 ms
~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 
get_skeleton_center_of_mass:
CPU times: user 6.52 s, sys: 104 ms, total: 6.62 s
Wall time: 6.62 s
~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 
get_furthest_point_from_edge:
CPU times: user 413 ms, sys: 63 µs, total: 413 ms
Wall time: 413 ms
~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 
get_furthest_point_from_edge_cv2:
CPU times: user 47.8 ms, sys: 0 ns, total: 47.8 ms
Wall time: 47.8 ms
~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 
get_center_of_half_area_line:
CPU times: user 1.66 s, sys: 0 ns, total: 1.66 s
Wall time: 1.66 s
~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 

Compute the centre position (x,y) of each edge of the polygon. You can do this by finding the difference between the positions of the ends of each edge. Take the average of each centre in each dimension. This will be the centre of the polygon.

Centroid method has already been suggested multiple times. I think this is an excellent resource that describes the process (and many other useful tricks with polygons) very intuitively:

http://paulbourke.net/geometry/polygonmesh/centroid.pdf

Also, for placing a simple UI label, it might be sufficient to just calculate the bounding box of the polygon (a rectangle defined by the lowest and highest x and y coordinates of any vertex in the polygon), and getting its center at:

{
    x = min_x + (max_x - min_x)/2,
    y = min_y + (max_y - min_y)/2
}

This is a bit faster than calculating the centroid, which might be significant for a real-time or embedded application.

Also notice, that if your polygons are static (they don't change form), you could optimize by saving the result of the BB center / center of mass calculation (relative to e.g. the first vertex of the polygon) to the data structure of the polygon.

I am not saying that this is the fastest, but it will give you a point inside the polygon. Calculate the Straight Skeleton. The point you are looking for is on this skeleton. You could pick the one with the shortest normal distance to the center of the bounding box for example.

How about finding the "incircle" of the polygon (the largest circle that fits inside it), and then centering the label at the center of that? Here are a couple of links to get you started:

http://www.mathopenref.com/polygonincircle.html
https://nrich.maths.org/discus/messages/145082/144373.html?1219439473

This won't work perfectly on every polygon, most likely; a polygon that looked like a C would have the label at a somewhat unpredictable spot. But the advantage would be that the label would always overlap a solid part of the polygon.

I think if you broke the polygon back into it's vertices, and then applied a function to find the largest convex hull , and then find the center off that convex hull, it would match closely with the "apparent" center.

Finding the largest convex hull given the vertices: Look under the Simple Polygon paragraph.

Average the vertices of the convex hull to find the center.

If I understand the point of the paper you linked to (quite an interesting problem, btw), this "inside buffering" technique is somewhat analogous to modeling the shape in question out of a piece of sugar that's being dissolved by acid from the edges in. (e.g. as the buffer distance increases, less of the original shape remains) The last bit remaining is the ideal spot to place a label.

How to accomplish this in an algorithm unfortunately isn't very clear to me....

Could you place the label at the naive center (of the bounding box, perhaps), and then move it based on the intersections of the local polygon edges and the label's BB? Move along the intersecting edges' normals, and if multiple edges intersect, sum their normals for movement?

Just guessing here; in this sort of problem I would probably try to solve iteratively as long as performance isn't too much of a concern.

Not much time to elaborate or test this right now, but I'll try to do more when I get a chance.

Use centroids as your primary method. Test to see if the centroid is within the polygon; if not, draw a line through the nearest point and on to the other side of the polygon. At the midpoint of the section of that line that is within the polygon, place your label.

Because the point that is nearest to the centroid is likely to bound a fairly large area, I think this might give results similar to Kyralessa's incircles. Of course, this might go berserk if you had a polygon with holes. In that case, the incircles would probably fair much better. On the other hand, it defaults to the (quick?) centroid method for the typical cases.

This problem would probably be analogous to finding the "center of mass" assuming a uniform density.

EDIT: This method won't work if the polygon has "holes"

you can use Center of Mass (or Center of Gravity) method which is used in civil engineering, here is a useful link from wikipedia:

http://en.wikipedia.org/wiki/Center_of_mass