百度一圈没有找到合适的博客,通过外网找到了python实现,所以整理记录一下。
最小圆问题
The smallest-circle problem (also known as minimum covering circle problem, bounding circle problem, smallest enclosing circle problem)
有多种算法求解实现,具体参见*。
最容易理解,使用最多的是一些计算几何算法,大家可以参考这篇博文最小圆覆盖(经典算法【三点定圆)进行理论学习
代码实现
scipy.optimize.leastsq
转化将上述最小圆问题,转化为半径最小时的圆心定位问题,再通过scipy中的最小二乘算法进行求解即可。
展示一个例子
from numpy import *
# Coordinates of the 2D points
x = r_[ 9, 35, -13, 10, 23, 0]
y = r_[ 34, 10, 6, -14, 27, -10]
basename = 'circle'
x_m = mean(x)
y_m = mean(y)
# Decorator to count functions calls
import functools
def countcalls(fn):
"decorator function count function calls "
@functools.wraps(fn)
def wrapped(*args):
wrapped.ncalls +=1
return fn(*args)
wrapped.ncalls = 0
return wrapped
# == METHOD 2 ==
# Basic usage of optimize.leastsq
from scipy import optimize
method_2 = "leastsq"
def calc_R(xc, yc):
""" calculate the distance of each 2D points from the center (xc, yc) """
return sqrt((x-xc)**2 + (y-yc)**2)
@countcalls
def f_2(c):
""" calculate the algebraic distance between the 2D points and the mean circle centered at c=(xc, yc) """
Ri = calc_R(*c)
return Ri - Ri.mean()
center_estimate = x_m, y_m
center_2, ier = optimize.leastsq(f_2, center_estimate)
xc_2, yc_2 = center_2
Ri_2 = calc_R(xc_2, yc_2)
R_2 = Ri_2.mean()
residu_2 = sum((Ri_2 - R_2)**2)
residu2_2 = sum((Ri_2**2-R_2**2)**2)
ncalls_2 = f_2.ncalls
# == METHOD 2b ==
# Advanced usage, with jacobian
method_2b = "leastsq with jacobian"
def calc_R(xc, yc):
""" calculate the distance of each 2D points from the center c=(xc, yc) """
return sqrt((x-xc)**2 + (y-yc)**2)
@countcalls
def f_2b(c):
""" calculate the algebraic distance between the 2D points and the mean circle centered at c=(xc, yc) """
Ri = calc_R(*c)
return Ri - Ri.mean()
@countcalls
def Df_2b(c):
""" Jacobian of f_2b
The axis corresponding to derivatives must be coherent with the col_deriv option of leastsq"""
xc, yc = c
df2b_dc = empty((len(c), x.size))
Ri = calc_R(xc, yc)
df2b_dc[ 0] = (xc - x)/Ri # dR/dxc
df2b_dc[ 1] = (yc - y)/Ri # dR/dyc
df2b_dc = df2b_dc - df2b_dc.mean(axis=1)[:, newaxis]
return df2b_dc
center_estimate = x_m, y_m
center_2b, ier = optimize.leastsq(f_2b, center_estimate, Dfun=Df_2b, col_deriv=True)
xc_2b, yc_2b = center_2b
Ri_2b = calc_R(xc_2b, yc_2b)
R_2b = Ri_2b.mean()
residu_2b = sum((Ri_2b - R_2b)**2)
residu2_2b = sum((Ri_2b**2-R_2b**2)**2)
ncalls_2b = f_2b.ncalls
print ("""
Method 2b :
print "Functions calls : f_2b=%d Df_2b=%d""" % ( f_2b.ncalls, Df_2b.ncalls))
print (method_2 , xc_2 , yc_2 , R_2 , ncalls_2 , Ri_2.std() , residu_2 , residu2_2 )
参考
https://scipy-cookbook.readthedocs.io/items/Least_Squares_Circle.html
Project Nayuki
网页提供了不同语言的算法源程序,可以供大家参考学习
算法实现参考
https://www.nayuki.io/page/smallest-enclosing-circle