Smallest enclosing circle for a convex polygon

Question

The smallest enclosing circle of a set of 2D points is known to be computable in linear time by the Megiddo or Welzl methods. (https://en.wikipedia.org/wiki/Smallest-circle_problem). But the first is not practical and the second is randomized.

Is there an ad-hoc solution that meets the $O(N)$ bound when the points are the vertices of a convex polygon ? (Note that as the convex hull of a simple polygon is found in linear time by Melkman's algorithm, the convex requirement can be replaced by simple.)

Answer 1

As a naive idea we begin with an assertion. Given a convex polygon, the two most distant vertices are contained in the smallest enclosing circle. With this idea we follow considering all the triangles that can be formed using the two more distant vertices and other point pertaining to the polygon. The smallest circle should contain all the triangles. Also we should consider the anomalous "triangle" formed by the two vertices most distant. Follows a MATHEMATICA script which computes the best circle. At the beginning in data there are a sequence of cases which can be tested.

data = {{1.7`, 0.27`}, {2.84`, 0.75`}, {3.96`, 0.51`}, {4.93`, -0.61`}, {3.62`, -0.87`}, {2.17`, -0.76`}};
data = {{3.33`, 0.54`}, {3.82`, 0.81`}, {5.06`, 0.83`}, {5.61`, 0.26`}, {5.67`, -0.66`}, {4.82`, -0.76`}, {3.17`, -0.23`}};
data = {{2.41`, 0.34`}, {2.75`, 0.76`}, {3.53`, 0.96`}, {4.61`, 0.84`}, {5.69`, 0.42`}, {5.94`, -0.7`}, {3.39`, -0.87`}, {2.63`, -0.36`}};
data = {{2.65`, 0.87`}, {4.24`, 1.39`}, {5.82`, -0.1`}, {5.22`, -1.92`}, {3.07`, -1.89`}, {2.23`,-0.52`}};
data = {{2.63`, 1.11`}, {3.37`, 1.44`}, {4.8`, 1.29`}, {5.9`, -1.45`}, {4.35`, -2.07`}, {3.04`, -2.12`}, {2.32`, -1.31`}, {2.01`, 0.`}, {2.18`, 0.61`}};
data = {{5.23`, 1.2`}, {4.36`, 1.59`}, {3.34`, 1.51`}, {2.51`,  0.93`}, {2.12`, -0.82`}, {2.54`, -1.69`}, {3.34`, -2.3`}, {4.53`,-2.2`}, {5.73`, -1.47`}, {5.89`, -0.32`}, {5.62`, 0.65`}};
data = {{3.38`, -2.37`}, {2.71`, -2.02`}, {2.27`, -1.21`}, {2.25`, -0.36`}, {2.43`, 0.27`}, {2.87`, 0.91`}, {3.47`, 1.14`}, {4.1`, 1.22`}, {4.75`, 1.18`}, {5.26`, 0.87`}, {5.76`, 0.4`}, {6.01`, -0.65`}, {5.91`, -1.15`}, {5.63`, -1.68`}, {4.91`, -2.15`}, {4.07`, -2.48`}};
error = 0.00001;

verify[data_List, solution_] := Module[{internal = True, i}, 
   For[i = 1, i <= Length[data], i++, 
      If[Norm[data[[i]] - solution[[1]]] - solution[[2]] > error, internal = False]
      ]; 
 Return[internal]
 ]

 cr[tri_List] := Module[{x1, x2, x3, y1, y2, y3, x0, y0, r},
    {x1, y1} = tri[[1]]; {x2, y2} = tri[[2]]; 
    {x3, y3} = tri[[3]]; 
    r = Sqrt[(x1 - x0)^2 + (y1 - y0)^2]; 
    Return[{{x0, y0}, r} /. Solve[{x1^2 - x2^2 + y1^2 - y2^2 - 2 x0 (x1 - x2) - 2 y0 (y1 - y2) == 0,
                                   x3^2 - x2^2 + y3^2 - y2^2 - 2 x0 (x3 - x2) - 2 y0 (y3 - y2) == 0}, {x0, y0}][[1]]]
  ]

n = Length[data];
dmax = 0;
For[i = 1, i <=  n, i++, 
  For[j = 1, j < i, j++, d = Norm[data[[i]] - data[[j]]]; 
    If[d > dmax, dmax = d; i0 = i; j0 = j]]]
    v = data[[i0]] - data[[j0]];
    p0 = (data[[i0]] + data[[j0]])/2;
    r = Norm[v]/2;
tri = {};
For[i = 1, i <= n, i++, 
  If[(i != i0) && (i != j0), 
    AppendTo[tri, {data[[i0]], data[[i]], data[[j0]]}]]];
circs = Table[cr[tri[[k]]], {k, 1, Length[tri]}]
nc = Length[circs]
If[verify[data, {p0, r}], feasible = {{p0, r}}, feasible = {}]
For[i = 1, i <= nc, i++,
  choosed = verify[data, circs[[i]]];
  If[choosed,
    AppendTo[feasible, circs[[i]]]
  ]
]
rmin = 1000;
For[i = 1, i <= Length[feasible], i++, 
   If[feasible[[i, 2]] < rmin, rmin = feasible[[i, 2]]; ibest = i]]
best = feasible[[ibest]]


xmin = Min[Transpose[data][[1]]] - 1;
xmax = Max[Transpose[data][[1]]] + 1;
ymin = Min[Transpose[data][[2]]] - 1;
ymax = Max[Transpose[data][[2]]] + 1;
aspectratio = (ymax - ymin)/(xmax - xmin);
pc = best[[1]]
r = best[[2]]
gr0 = ListLinePlot[Join[data, {data[[1]]}], PlotStyle -> {Thick, Blue}];
gr1 = ParametricPlot[u data[[i0]] + (1 - u) data[[j0]], {u, 0, 1}, PlotStyle -> Red];
gr3 = Graphics[Circle[pc, r]];
gr4 = Graphics[{Red, PointSize[0.02], Point[pc]}];
gr5 = Table[ListLinePlot[tri[[k]]], {k, 1, Length[tri]}];
Show[gr0, gr1, gr3, gr4, gr5, AspectRatio -> aspectratio,PlotRange -> {{xmin, xmax}, {ymin, ymax}}]

Follows a python version to the essential procedures excluding plotting

import numpy as np
from numpy import linalg as LA

data =[[3.38, -2.37],
   [2.71, -2.02],
   [2.27, -1.21],
   [2.25, -0.36],
   [2.43, 0.27],
   [2.87, 0.91],
   [3.47, 1.14],
   [4.1, 1.22],
   [4.75, 1.18],
   [5.26,0.87],
   [5.76, 0.4],
   [6.01, -0.65],
   [5.91, -1.15],
   [5.63, -1.68],
   [4.91, -2.15],
   [4.07, -2.48]]

data = [[3.33, 0.54],
    [3.82, 0.81],
    [5.06, 0.83],
    [5.61, 0.26],
    [5.67, -0.66],
    [4.82, -0.76],
    [3.17, -0.23]]

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 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.00001
    for i in range(len(data)):
        if LA.norm(sub(data[i],feasible[0]))-feasible[1] > 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

(i0, j0) = max_secant(data)
v = sub(data[i0],data[j0])
r = 0.5*LA.norm(v)
p1 = add(data[i0],data[j0])
p0 = [element*0.5 for element in p1]
triangles = collect_triangles(data,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(data,polar[i]):
        feasible.append(polar[i])
    
print("*** best circle ***")
    
bestr = 10000
for i in range(len(feasible)):
    [pc, r] = feasible[i]
    if r < bestr:
        bestr = r
        bestcirc = feasible[i]
    
print(bestcirc)