Skip to main content
added 87 characters in body
Source Link
Cesareo
  • 111
  • 4

enter image description here

enter image description here

Source Link
Cesareo
  • 111
  • 4

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