6
$\begingroup$

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

$\endgroup$
7
  • 2
    $\begingroup$ Good question, but isn’t it better suited for the computer science stackexchange? $\endgroup$
    – Milten
    Apr 17 at 7:02
  • $\begingroup$ @Milten: Mathematics has the computational-geometry tag, and CG problems are not infrequent. $\endgroup$
    – user65203
    Apr 17 at 7:36
  • $\begingroup$ Fair enough then. $\endgroup$
    – Milten
    Apr 17 at 7:42
  • $\begingroup$ What would happen if you ran a non-randomized version of Welzl's algorithm on a convex polygon? $\endgroup$
    – Théophile
    Apr 22 at 13:59
  • 1
    $\begingroup$ @Théophile: interesting question. I don't know but I assume that if that lead to a linear-time solution for a convex polygon (or convex hull), that result would be well-known. $\endgroup$
    – user65203
    Apr 22 at 14:08
2
+100
$\begingroup$

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}}]

enter image description here enter image description here enter image description here

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)
$\endgroup$
6
  • $\begingroup$ Sorry, I am not after some code to solve the problem. I am asking if an O(N) solution exists. $\endgroup$
    – user65203
    Apr 21 at 21:43
  • $\begingroup$ This conclusion can be obtained after the code. Finding the base segment, the sequence of triangles and choosing the best. All this seems $O(N)$ $\endgroup$
    – Cesareo
    Apr 21 at 21:47
  • $\begingroup$ How do you find the two most distant vertices ? I see a double loop, which does not smell $O(N)$. $\endgroup$
    – user65203
    Apr 21 at 21:52
  • $\begingroup$ Yes. You are right. I will think if it is possible to reduce that. $\endgroup$
    – Cesareo
    Apr 21 at 22:03
  • 1
    $\begingroup$ @Yves: Here it is claimed that the diameter of a convex polygon can be computed in O(n): stackoverflow.com/a/321994/1187415 $\endgroup$
    – Martin R
    Apr 22 at 6:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy