Largest circle inside a non-convex polygon

How can I find the largest circle that can fit inside a concave polygon?

A brute force algorithm is OK as long as it can handle polygons with ~50 vertices in real-time.

The key to solving this problem is first making an observation: the center of the largest circle that will fit inside an arbitrary polygon is the point that is:

1. Inside the polygon; and
2. Furthest from any point on the edges of the polygon.

Why? Because every point on the edge of a circle is equidistant from that center. By definition, the largest circle will have the largest radius and will touch the polygon on at least two points so if you find the point inside furthest from the polygon you've found the center of the circle.

This problem appears in geography and is solved iteratively to any arbitrary precision. Its called the Poles of Inaccessibility problem. See Poles of Inaccessibility: A Calculation Algorithm for the Remotest Places on Earth.

The basic algorithm works like this:

1. Define R as a rectilinear region from (x_min,y_min) to (x_max,y_max);
2. Divide R into an arbitrary number of points. The paper uses 21 as a heuristic (meaning divide the height and width by 20);
3. Clip any points that are outside the polygon;
4. For the remainder find the point that is furthest from any point on the edge;
5. From that point define a new R with smaller intervals and bounds and repeat from step 2 to get to any arbitrary precision answer. The paper reduces R by a factor of the square root of 2.

One note, How to test if a point is inside the polygon or not: The simplest solution to this part of the problem is to cast a ray to the right of the point. If it crosses an odd number of edges, it's within the polygon. If it's an even number, it's outside.

Also, as far as testing the distance to any edge there are two cases you need to consider:

1. The point is perpendicular to a point on that edge (within the bounds of the two vertices); or
2. It isn't.

(2) is easy. The distance to the edge is the minimum of the distances to the two vertices. For (1), the closest point on that edge will be the point that intersects the edge at a 90 degree angle starting from the point you're testing. See Distance of a Point to a Ray or Segment.

In case anyone is looking for a practical implementation, I designed a faster algorithm that solves this problem for a given precision and made it a JavaScript library. It's similar to the iterative grid algorithm described by @cletus, but it's guaranteed to obtain global optimum, and is also 20-40 times faster in practice.

Check it out: https://github.com/mapbox/polylabel

An O(n log(n)) algorithm:

1. Construct the Voronoi Diagram of the edges in P. This can be done with, for example, Fortunes algorithm.
2. For Voronoi nodes (points equidistant to three or more edges) inside P;
3. Find the node with the maximum distance to edges in P. This node is the centre of the maximum inscribed circle.

Summary: In theory, this can be done in O(n) time. In practice you can do it in O(n log n) time.

Generalized Voronoi diagrams.

If you consider the vertices and edges of the polygon as a set of sites and tessellate the interior into the "nearest neighbor cells" then you get the so-called (generalized) Voronoi diagram. The Voronoi diagram consists of nodes and edges connecting them. The clearance of a node is the distance to its defining polygon faces.

(Here the polygon even has holes; the principle still works.)

The key observation now is that the center of the maximum inscribed circle touches three faces (vertices or edges) of the polygon, and no other face can be closer. So the center has to lie on a Voronoi node, i.e, the node with the largest clearance.

In the example above the node that marks the center of the maximum inscribed circle touches two edges and a vertex of the polygon.

The medial axis, by the way, is the Voronoi diagram with those Voronoi edges removed that emanate from reflex vertices. Hence, the center of the maximum inscribed circle also lies on the medial axis.

Source: A blog article of mine that deals with generalizations of maximum inscribed circles at some point. There you can find more on Voronoi diagrams and their relation to maximum inscribed circles.

Algorithms & implementations.

You could actually compute the Voronoi diagram. A worst-case O(n log n) algorithm for points and segments is given by Fortune, A sweepline algorithm for Voronoi diagrams, SoCG'86. Held published the software package Vroni with an expected O(n log n) time complexity, which actually computes the maximum inscribed circle, too. And there seems to be an implementation in boost, too.

For simple polygons (i.e., without holes) a time-optimal algorithm that runs in O(n) time is due to Chin et al., Finding the Medial Axis of a Simple Polygon in Linear Time, 1999.

Brute force.

However, as you stated that you are fine with a brute-force algorithm: What about simply trying out all triplets of sites (vertices and edges). For each triplet you find candidate Voronoi nodes, i.e., equidistant loci to the three sites and check whether any other site would intersect the candidate maximum inscribed circle. If there is an intersection you dismiss the candidate. Take the greatest you can find over all triplets.

See chapter 3 in my Master thesis about more details on computing equidistant loci for three sites.

I implemented a piece of python code based on cv2 to get the maximum/largest inscribed circle inside mask/polygon/contours. It supports non-convex/hollow shape.

import cv2
import numpy as np
def get_test_mask():
    # Create an image
    r = 100
    mask = np.zeros((4 * r, 4 * r), dtype=np.uint8)

    # Create a sequence of points to make a contour
    vert = [None] * 6
    vert[0] = (3 * r // 2, int(1.34 * r))
    vert[1] = (1 * r, 2 * r)
    vert[2] = (3 * r // 2, int(2.866 * r))
    vert[3] = (5 * r // 2, int(2.866 * r))
    vert[4] = (3 * r, 2 * r)
    vert[5] = (5 * r // 2, int(1.34 * r))
    # Draw it in mask
    for i in range(6):
        cv2.line(mask, vert[i], vert[(i + 1) % 6], (255), 63)
    return mask


mask = get_test_mask()

"""
Get the maximum/largest inscribed circle inside mask/polygon/contours.
Support non-convex/hollow shape
"""
dist_map = cv2.distanceTransform(mask, cv2.DIST_L2, cv2.DIST_MASK_PRECISE)
_, radius, _, center = cv2.minMaxLoc(dist_map)

result = cv2.cvtColor(mask, cv2.COLOR_GRAY2BGR)
cv2.circle(result, tuple(center), int(radius), (0, 0, 255), 2, cv2.LINE_8, 0)

# minEnclosingCircle directly by cv2
contours, _ = cv2.findContours(mask, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)[-2:]
center2, radius2 = cv2.minEnclosingCircle(np.concatenate(contours, 0))
cv2.circle(result, (int(center2[0]), int(center2[1])), int(radius2), (0, 255, 0,), 2)

cv2.imshow("mask", mask)
cv2.imshow("result", result)
cv2.waitKey(0)

Red circle is max inscribed circle

Source: https://gist.github.com/DIYer22/f82dc329b27c2766b21bec4a563703cc

I used Straight Skeletons to place an image inside a polygon with three steps:

1. Find the straight skeleton using the Straight Skeleton algorithm (pic 1)
2. Base on the straight skeleton, find the largest circle (pic 2)
3. Draw the image inside that circle (pic 3)

Try it at: https://smartdiagram.com/simple-infographics-3d-charts-2/

An O(n log X) algorithm, where X depends on the precision you want.

Binary search for largest radius R for a circle:

At each iteration, for a given radius r, push each edge E, "inward" by R, to get E'. For each edge E', define half-plane H as the set of all points "inside" the the polygon (using E' as the boundary). Now, compute the intersection of all these half-planes E', which could be done in O(n) time. If the intersection is non-empty, then if you draw a circle with radius r using any point in the intersection as the center, it will be inside the given polygon.