Write an algorithm in Python to find the smallest circle that can enclose a polygon.
The polygon is defined as a set of coordinates that map the vertices inside a 2D plane.
Here is an answer using cvxpy and scipy.
By default, this finds a circle that encloses a square. There are options for generating regular polygons (specifying number of points and radius) or random points (specifying number of points, minimum and maximum coordinates). SciPy is used to find the convex hull of the polygon points.
import sys
import math
import random
import cvxpy as cp
from scipy.spatial import ConvexHull
def get_convex_hull(vertices_set):
'''
Compute the convex hull of the polygon.
The convex hull of a set of points is the smallest convex set that
contains all of the points. The convex hull of a concave polygon contains
fewer vertices than the original polygon. This should reduce cvxpy
constraints and run time for concave polygons with many vertices.
Uses scipy.ConvexHull.
:param vertices_set: Set of (x, y) tuples representing the polygon vertices
:return: hull_vertices Set of (x, y) tuples representing the convex hull
of the polygon vertices
'''
# Convert set to list for ConvexHull function.
vertices = list(vertices_set)
vertices_len = len(vertices)
# Compute the convex hull of the polygon.
hull = ConvexHull(vertices)
# Extract convex hull vertices.
hull_vertices = [vertices[hull_vertex] for hull_vertex in hull.vertices]
hull_vertices_set = set(hull_vertices)
hull_vertices_len = len(hull_vertices_set)
# print('Convex hull vertices:', hull_vertices_set)
# Warn about vertex removal.
if hull_vertices_len < vertices_len:
vertices_diff = vertices_len - hull_vertices_len
print(f'Warning - convex hull calculation removes {vertices_diff} vertices:')
print('t', vertices_set - hull_vertices_set)
return hull_vertices_set
def check_vertices(vertices):
'''
Check vertices for problems.
Must be at least 3 vertices.
Vertices must all be tuples.
Vertices must all be 2 dimensional.
:param vertices: Set of (x, y) tuples representing the polygon vertices
'''
# Must be at least 3 vertices.
if len(vertices) < 3:
print('Must be at least 3 vertices:', vertices)
sys.exit()
# Vertices must all be tuples.
if not all(isinstance(vertex, tuple) for vertex in vertices):
print('Vertices must all be tuples:', vertices)
sys.exit()
# Vertices must all be 2 dimensional.
vertex_sizes = [len(p) for p in vertices]
if not all(vertex_size == 2 for vertex_size in vertex_sizes):
vertex_size_min = min(vertex_sizes)
vertex_size_max = max(vertex_sizes)
print('Vertices must be 2 dimensional:', vertices)
print('Vertex sizes:', vertex_sizes)
print(f'Vertex size - min, max: {vertex_size_min}, {vertex_size_max}')
sys.exit()
def minimum_enclosing_circle(vertices):
'''
Find minimum enclosing circle for a polygon.
Uses cvxpy with Euclidean distance constraints.
Euclidean distance from each vertex to the circle center should be
less than or equal to circle radius.
The cvxpy objective is to minimise the enclosing circle radius.
:param vertices: Set of (x, y) tuples representing the polygon vertices
:return: Center (x, y) and radius of the minimal enclosing circle
'''
# Variables for the center of the circle (c_x, c_y) and the radius.
center = cp.Variable(2)
radius = cp.Variable(nonneg=True)
# Objective: Minimize the radius.
objective = cp.Minimize(radius)
# Constraints:
# Euclidean distance from each vertex to the center should be <= radius
# cp.hstack - Horizontal concatenation of arguments
# cp.norm - Euclidean distance
constraints = [cp.norm(cp.hstack(vertex - center), 2) <= radius
for vertex in vertices]
# Formulate the problem.
problem = cp.Problem(objective, constraints)
# Solve the problem.
problem.solve()
# Check problem status.
if problem.status != cp.OPTIMAL:
print('Solver failed to find a solution')
sys.exit()
# Return the center coordinates and the radius.
return center.value, radius.value
def describe_circle(center, radius):
'''
Print circle position and radius.
:param center: Center (x, y) coordinates of the minimal enclosing circle
:param radius: Radius of the minimal enclosing circle
'''
(c_x, c_y) = center
print(f'Circle center: ({c_x:.4f}, {c_y:.4f})')
print(f'Circle radius: {radius:.4f}')
def get_vertices(num_points=5, radius=None, min_coord=None, max_coord=None):
'''
Either generate regular polygon, generate random points are use hardcoded
polygon.
A polygon is defined as a set of coordinates that map the vertices
inside a 2D plane.
:param num_points: number of points
:param radius: radius
:param min_coord: minimum coordinate
:param max_coord: maximum coordinate
:return: Set of tuples containing 2D coordinates
'''
if radius is not None:
# Generate regular polygon.
vertices = generate_regular_polygon(num_points, radius)
elif min_coord is not None and max_coord is not None:
# Generate set of random 2D tuples.
vertices = generate_random_points(num_points, min_coord, max_coord)
else:
# Use hardcoded vertices.
# Triangle
# vertices = {(0, 0), (2, 0), (2, 2)}
# Square
vertices = {(0, 0), (0, 2), (2, 0), (2, 2)}
# Irregular pentagon
# vertices = {(0, 0), (0, 2), (2, 0), (2, 2), (1, 3)}
# Bow tie - concave
# vertices = {(0, 0), (1, 1), (2, 2), (2, 4), (1, 3), (0, 4)}
# 5-sided star
# vertices = {(0, 100), (59, 81), (95, 31), (36, -12), (59, -81),
# (-59, -81), (-36, -12), (-95, 31), (-59, 81)}
print('Hardcoded polygon:', vertices)
return vertices
def main():
'''
Standard main function.
Get hardcoded vertices or generate a polygon or random points on 2D plane.
Check for problems with points/vertices.
Calculate convex hull of points/vertices.
Find the minimal enclosing circle of the convex hull.
Describe the minimal enclosing circle.
'''
# Generate a polygon
# vertices = get_vertices(num_points=5, radius=5)
# Generate random points on 2D plane
# vertices = get_vertices(num_points=5, min_coord=0, max_coord=10)
# Get hardcoded vertices - square by default
vertices = get_vertices()
# Check for problems with points/vertices.
check_vertices(vertices)
# Calculate convex hull of the points/vertices with scipy.
convex_hull = get_convex_hull(vertices)
# Find the minimal enclosing circle of the convex hull with cvxpy.
cen, rad = minimum_enclosing_circle(convex_hull)
# Print description of minimal enclosing circle.
describe_circle(cen, rad)
def generate_regular_polygon(num_points, radius):
'''
Generates coordinates of a regular polygon.
:param num_points: number of points
:param radius: radius
:return: Set of tuples containing 2D coordinates of regular polygon
'''
# Sanity check parameters.
assert isinstance(radius, int)
assert isinstance(num_points, int)
assert radius > 0
assert num_points >= 3
# Generate regular polygon.
regular_coords = []
angle = 2 * math.pi / num_points
for i in range(num_points):
# coords or tuples containing x, y coordinates
coords = (radius * math.cos(i * angle), radius * math.sin(i * angle))
regular_coords.append(coords)
# Convert list to set.
regular_coords_set = set(regular_coords)
print('Generated regular polygon:', regular_coords_set)
return regular_coords_set
def generate_random_points(num_points, min_coord, max_coord):
'''
Generates at least 3 random 2D points for testing.
A set of n random points may not result in an n-sided polygon.
:param num_points: number of points
:param min_coord: minimum coordinate
:param max_coord: maximum coordinate
:return: Set of tuples containing random integer 2D coordinates
'''
# Sanity check parameters.
assert isinstance(min_coord, int)
assert isinstance(max_coord, int)
assert isinstance(num_points, int)
assert num_points >= 3
assert min_coord < max_coord
# Generate random 2D integer points.
random_coords = []
for _ in range(num_points):
# coords or tuples containing random integer coordinates
coords = (random.randint(min_coord, max_coord),
random.randint(min_coord, max_coord))
random_coords.append(coords)
# Convert list to set.
random_coords_set = set(random_coords)
print('Generated random 2D integer points:', random_coords_set)
return random_coords_set
if __name__ == "__main__":
# Call main function.
main()
For the default square, this returns:
Hardcoded polygon: {(2, 2), (0, 2), (2, 0), (0, 0)}
Circle center: (1.0000, 1.0000)
Circle radius: 1.4142
CVXPY may be overpowered for this problem.