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

Question

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

Shenika

What if the set defined by the polygon is (highly) non-convex (en.wikipedia.org/wiki/Convex_set); is it allowed to have the center outside the polygon? – Megrim

Yes, but for the purpose of labeling, we would need to find a point inside. – Shenika

@Mikhil: to expand on @Pukku's comment, could you please post a "hard" aspect of this problem, i.e. a shape that would be difficult to label given "naive" answers such as center-of-mass? The ones I can think of easily are a giant U or the state of Florida (center of mass of these shapes are outside the boundary) – Acetic

Related: #22797020 Leaflet seems to have this capability built in – Affiant

See https://mcmap.net/q/303274/-largest-circle-inside-a-non-convex-polygon for an answer. – Kolyma

Answer 1

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.

Megrim

Also see #1110036 – Extend

I think these are your best bets by far. You could adapt the above by vertically stretching the polygon by a factor of 2 or 3, then searching for the largest circle contained in the stretched polygon. This will give you the largest ellipse contained within the polygon, which will give you the best spot to put your label. – Adige

The straight skeleton does not give you the center of the largest inscribed disk, but the medial axis (and Voronoi diagram) does: https://mcmap.net/q/303274/-largest-circle-inside-a-non-convex-polygon – Kolyma

Answer 2

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

Blithering

The straight skeleton does not give you the center of the largest inscribed disk, but the medial axis (and Voronoi diagram) does: https://mcmap.net/q/303274/-largest-circle-inside-a-non-convex-polygon – Kolyma

This answer really helped me! I needed this in Dart, so I ported it: pub.dev/packages/polylabel – Grous

Answer 3

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:

Vashtivashtia

Love it mate! Really clever! Now in terms of implementation, you and or anyone else solve? – Yclept

@MaraisRossouw I've posted an answer to a similar question to OP's that uses this method: https://mcmap.net/q/303275/-using-sql-spatial-data-c-to-find-the-quot-visual-quot-center-of-irregular-polygons – Coadjutrix

Sorry for "self-citing", but you can see my answer below (https://mcmap.net/q/301186/-what-is-the-fastest-way-to-find-the-quot-visual-quot-center-of-an-irregularly-shaped-polygon) for a python implementation (get_center_of_half_area_line) . – Diversity

Answer 4

Have you looked into using the centroid formula?

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

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

Mopup

Centroid == Center of Mass at Uniform Density – Liva

Answer 5

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.

Fraya

Good thinking, but doesn't always work, because the center of the bounding box can be far outside the polygon itself. !Center of bounding box outside polygon (img) – Acumen

Answer 6

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

Diversity

Do any of these happen to be the same algorithm as Polylabel? (mentioned above https://mcmap.net/q/301186/-what-is-the-fastest-way-to-find-the-quot-visual-quot-center-of-an-irregularly-shaped-polygon) – Affiant

Yes, it seems that get_furthest_point_from_edge implements the same function as Polylabel: "A fast algorithm for finding -- the most distant internal point from the polygon outline" (github.com/mapbox/polylabel) – Diversity

Thanks for the clarification. I've been using polylabel for a while mainly because it produced a point in a useful location, and seemed stable/reliable. But this comparison is more objective, overall, than the info here before. – Affiant

Answer 7

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.

Limemann

I think this suffers the same problem as my solution when it comes to highly non-convex shapes... – Liva

Yep, and without taking a weighted average it also over-emphasizes short edges, even if the polygon is convex. – Megrim

Answer 8

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.

Astrobiology

The straight skeleton does not give you the center of the largest inscribed disk, but the medial axis (and Voronoi diagram) does: https://mcmap.net/q/303274/-largest-circle-inside-a-non-convex-polygon – Kolyma

Answer 9

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.

Apuleius

Won't this be slow if a polygon has several triangulations? – Shenika

Answer 10

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.

Castiglione

It would select one of the sides. What is the desired behavior in that situation? – Castiglione

For a giant U, an acceptable solution is the middle of the lower thick section. – Shenika

If the lower thick section is the largest convex hull, then it would be selected. Is there some type of criteria for selected convex hull to be more of a square? – Castiglione

Wont the largest convex hull cover the entire U and be a rectangle? – Shenika

Oh, you would need to modify the algorithm to not include any interior vertices. – Castiglione

Answer 11

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

Acetic

GIS softwares like PostGIS have functions like ST_Buffer that do this. I don't know how, so quickly. – Shenika

Answer 12

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.

Galluses

Answer 13

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.

Maighdiln

Pathological test case #3: a symmetrical barbell-like shape with a thin rectangle and two big octagons on the ends. Centroid is within the polygon but the rectangle is a bad place to label since it may not fit. – Acetic

Answer 14

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"

Liva

No. See the figure #4 in the ESRI paper that the OP linked to. – Acetic

It seems that my assumption is what they used in #2; the only time it breaks down is under this condition: "However, this method provides an incorrect result if a polygon has holes" – Liva

No. Imagine a giant U. There are no holes, and the center of mass is not inside the boundary of the polygon. I think your answer is correct only for convex polygons. – Acetic

Thank you; it'd help if the asker gave us some boundary conditions to work with as well! – Liva

Answer 15

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

Blacksmith