In this reference we have an algorithm to determine the smallest circle containing a convex polygon. Follows a python code which uses this algorithm to find the smallest semi-circle container. The focused example is the same referenced in the OP in the cited paper. I hope the script is self explained. I am a early python programmer...
import math
import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt
from shapely.geometry import Polygon
data0 = [[2.30, 0.15],[0.63, 0.41],[0.37, 0.59],[0.79, 1.47],[2.32, 1.87],[3.6107, 0.72],[2.73, 0.14]]
def sub(p1,p2):
return list(map(lambda i, j: i-j,p1,p2))
def add(p1,p2):
return list(map(lambda i, j: i+j,p1,p2))
def cline(p,v,u):
v = [element * u for element in v]
return list(map(lambda i, j: i+j,p,v))
def max_secant(data):
n = len(data)
dmax = 0
for i in range(n):
for j in range(i):
d = LA.norm(sub(data[i],data[j]))
if d > dmax:
dmax = d
i0 = i
j0 = j
return (i0, j0)
def verify(data, feasible):
internal = True
error = 0.005
for i in range(len(data)):
dif = LA.norm(sub(data[i],feasible[0]))-feasible[1]
if dif > error:
internal = False
return internal
def polar_form(triangle):
(x1,y1) = triangle[0]
(x2,y2) = triangle[1]
(x3,y3) = triangle[2]
M = np.array([[2*(x2-x1),2*(y2-y1)],[2*(x2-x3),2*(y2-y3)]])
b = np.array([-(x1**2-x2**2+y1**2-y2**2),-(x3**2-x2**2+y3**2-y2**2)])
(x0, y0) = list(np.linalg.solve(M,b))
r = LA.norm([x1-x0,y1-y0])
return [[x0,y0], r]
def collect_triangles(data, i0, j0):
triangs = []
for i in range(len(data)):
if i not in [i0, j0]:
triangs.append([data[i0],data[i],data[j0]])
return triangs
def rotate(data):
data0 = []
n = len(data)
dummy = data[0]
for i in range(n-1):
data0.append(data[i+1])
data0.append(dummy)
return data0
def take_extremals(data):
breaks = []
sant = 1
v = sub(data[1],data[0])
n = len(data)
for i in range(1,n-2):
s = np.sign(np.dot(v,sub(data[i+1],data[i])))
if (sant != s):
breaks.append(i)
sant = s
if len(breaks) == 1:
breaks.append(n-1)
return breaks
def mirror(data, p, v):
reflected = []
vn = LA.norm(v)
n = len(data)
v = [v[0]/vn,v[1]/vn]
for i in range(n):
v0 = np.dot(sub(data[i],p),v)
v1 = [v[0]*v0,v[1]*v0]
v2 = add(p, v1)
v2 = [2*v2[0],2*v2[1]]
pr = sub(v2, data[i])
reflected.append(pr)
return reflected
def glue(data1, data2):
sdata = []
n1 = len(data1)
for i in range(n1):
sdata.append(data1[i])
n2 = len(data2)
for i in range(n2):
sdata.append(data2[n2-i-1])
return sdata
def select(data, k1, k2):
datas = []
for i in range(k1, k2+1):
datas.append(data[i])
return datas
def best_circle(data1):
(k1, k2) = take_extremals(data1)
p0b = data1[0]
vb = sub(data1[1],data1[0])
datas = select(data1, k1, k2)
datam = mirror(datas,p0b,vb)
dataf = glue(datam, datas)
(i0, j0) = max_secant(dataf)
v = sub(dataf[i0],dataf[j0])
r = 0.5*LA.norm(v)
p1 = add(dataf[i0],dataf[j0])
p0 = [element*0.5 for element in p1]
triangles = collect_triangles(dataf,i0,j0)
polar = []
polar.append([p0,r])
for i in range(len(triangles)):
polar.append(polar_form(triangles[i]))
feasible = []
for i in range(len(polar)):
if verify(dataf,polar[i]):
feasible.append(polar[i])
bestr = math.inf
for i in range(len(feasible)):
[pc, r] = feasible[i]
if r < bestr:
bestr = r
bestcirc = feasible[i]
return(bestcirc, p0b, vb)
########################
#### main program ####
########################
data1 = data0
circmin = math.inf
for i in range(len(data0)):
(circ, p0x, vx) = best_circle(data1)
if circ[1] < circmin:
circmin = circ[1]
bestcirc = circ
p0b = p0x
vb = vx
data1 = rotate(data1)
print(bestcirc)
#############################
#### plotting the result ####
#############################
(figure, axes) = plt.subplots()
(cx,cy) = bestcirc[0]
r = bestcirc[1]
poly = Polygon(data0)
(x, y) = poly.exterior.xy
xmin = cx - 1.1*r
xmax = cx + 1.1*r
ymin = cy - 1.1*r
ymax = cy + 1.1*r
axes.set_xlim((xmin,xmax))
axes.set_ylim((ymin,ymax))
uncolored_circle = plt.Circle( (cx,cy), r, fill = False)
axes.set_aspect( 1 )
axes.add_artist( uncolored_circle )
plt.plot(x,y)
v12 = vb
nv12 = LA.norm(v12)
v12 = [v12[0]/nv12,v12[1]/nv12]
s1x = cx - v12[0]*r
s1y = cy - v12[1]*r
s2x = cx + v12[0]*r
s2y = cy + v12[1]*r
plt.plot([s1x,s2x],[s1y,s2y])
plt.title( 'Result' )
plt.show()