Largest circle inside a non-convex polygon

Question

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.

Solution 1

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.

Solution 2

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

Solution 3

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.

Solution 4

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.

Solution 5

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

View more solutions

Comments

[Plow]

Seems like an algorithm that is fairly straight-forward to implement, which is exactly what I am looking for. According to the article there is no guarantee that the solution found is an absolute maximum however (for my particular case this may not be a problem).

[brainjam]

You want the Voronoi diagram of the edges, not the vertices. See, for example valis.cs.uiuc.edu/~sariel/research/CG/applets/medial_axis/…. The edge Voronoi diagram has curved segments, the vertex Voronoi diagram has only straight lines. Another name for what you want is "medial axis". Here's a site that discusses the difference: groups.csail.mit.edu/graphics/classes/6.838/S98/meetings/m25‌​/…

[Doc Brown]

I think this algorithm can be modified to find the absolute maximum for sure. The idea is to calcualte two values for each rectangle: a lower limit for the maximum distance from the polygon edge (the maximum distance of the 4 vertices of the rectangle), and an upper limit (by adding 0.5*sqrt(rect_size_x^2 + rect_size_y^2). Then, run a subdividing search which keeps all non-processed candidate rectangles in a priority queue (ordered descending by the upper limit) and throws away every rectangle with an upper limit below the biggest lower limit found so far.

[Blablaenzo]

To bad the link the broken ... another reference: arxiv.org/pdf/1212.3193.pdf

[Zaid Mirza]

is this available in Java?

[David Faivre]

I needed this in C#, so ported it: gist.github.com/dfaivre/acfef42cdbf411555956e9eba65dd30d

[StayOnTarget]

Related: stackoverflow.com/questions/1203135/…

[SeldomSeenSlim]

Great answer! This explanation allowed me to implement the solution in code in just a few minutes.

[J S]

Seems to require convexity of the polygon. For nonconvex polygons with or without holes, I could immediately construct examples where all intersections of any such set of half-planes will be empty, because there could be two edges that are "back to back".

[J S]

Ifs there a correctness proof or a quality estimate? This could clearly run into a local minimum if the points are not well chosen.

[StayOnTarget]

This is the same as the most popular answer for non-convex stackoverflow.com/a/40464906/3195477

[André Sousa]

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