如何将圆拟合到半径受约束的一组点? [英] How to fit a circle to a set of points with a constrained radius?
问题描述
我有一组代表小圆弧的点。
I have a set of points that represent a small arc of a circle.
当前代码使用线性最小二乘法将圆与这些点拟合:
The current code fits a circle to these points using linear least-squares:
void fit_circle(const std::vector<cv::Point2d> &pnts,cv::Point2d ¢re, double &radius)
{
int cols = 3;
cv::Mat X( static_cast<int>(pnts.size()), cols, CV_64F );
cv::Mat Y( static_cast<int>(pnts.size()), 1, CV_64F );
cv::Mat C;
if (int(pnts.size()) >= 3 )
{
for (size_t i = 0; i < pnts.size(); i++)
{
X.at<double>(static_cast<int>(i),0) = 2 * pnts[i].x;
X.at<double>(static_cast<int>(i),1) = 2 * pnts[i].y;
X.at<double>(static_cast<int>(i),2) = -1.0;
Y.at<double>(static_cast<int>(i),0) = (pnts[i].x * pnts[i].x + pnts[i].y * pnts[i].y);
}
cv::solve(X,Y,C,cv::DECOMP_SVD);
}
std::vector<double> coefs;
C.col(0).copyTo(coefs);
centre.x = coefs[0];
centre.y = coefs[1];
radius = sqrt ( coefs[0] * coefs[0] + coefs[1] * coefs[1] - coefs[2] );
}
但是,我有一个限制,即半径必须在特定范围,
,即 minRadius< = radius< = maxRadius
。
因此,应选择半径在该范围内的最佳拟合圆,而不是整体最佳拟合。
So the best fitting circle with a radius within that range should be selected, instead of the overall best fit.
如何修改代码以添加此约束?
推荐答案
希望我能得到你的问题权。第一个想法可能是:使用具有约束的拟合算法。但是,我不太喜欢这些内容,因为它们非常耗时。我的方法是通过使用类型 r = r_min +(r_max-r_min)* 0.5 *(1 + tanh(x))
的函数替换半径来约束半径并进行拟合 x
。因此, x
可以在整个实数范围内变化,从而在给定限制内提供半径。所有这些仍然是连续函数的简单非线性拟合。
I hope I get your problem right. First idea would probably be: use a fit algorithm with constrains. I do not like those very much, though, as they are very time consuming. My approach would be to constrain the radius by replacing it with a function of type r=r_min+(r_max-r_min)*0.5*(1+tanh(x))
and fitting x
. So x
can vary in the whole range of real numbers giving you radii in the given limits. All remains a simple non-linear fit of continuous functions.
作为python示例,它看起来像:
As a python example it looks like:
import matplotlib
matplotlib.use('Qt4Agg')
from matplotlib import pyplot as plt
from random import random
from scipy import optimize
import numpy as np
def boxmuller(x0,sigma):
u1=random()
u2=random()
ll=np.sqrt(-2*np.log(u1))
z0=ll*np.cos(2*np.pi*u2)
z1=ll*np.cos(2*np.pi*u2)
return sigma*z0+x0, sigma*z1+x0
def random_arcpoint(x0,y0, r,f1,f2,s):
arc=f1+(f2-f1)*random()
rr,_=boxmuller(r,s)
return [x0+rr*np.cos(arc),y0+rr*np.sin(arc)]
def residuals(parameters,dataPoint):
xc,yc,Ri = parameters
distance = [np.sqrt( (x-xc)**2 + (y-yc)**2 ) for x,y in dataPoint]
res = [(Ri-dist)**2 for dist in distance]
return res
def residualsLim(parameters,dataPoint,rMin,rMax):
xc,yc,rP = parameters
Ri=rMin+(rMax-rMin)*0.5*(1+np.tanh(rP))
distance = [np.sqrt( (x-xc)**2 + (y-yc)**2 ) for x,y in dataPoint]
res = [(Ri-dist)**2 for dist in distance]
return res
def f0(phi,x0,y0,r):
return [x0+r*np.cos(phi),y0+r*np.sin(phi)]
def f_lim(phi,rMin,rMax,x0,y0,x):
rr=(rMin+(rMax-rMin)*0.5*(1+np.tanh(x)))
return [x0+rr*np.cos(phi), y0+rr*np.sin(phi)]
data=np.array([random_arcpoint(2.1,-1.2,11,1.0,3.144,.61) for s in range(200)])
estimate = [0, 0, 10]
bestFitValues_01, ier_01 = optimize.leastsq(residuals, estimate,args=(data))
print bestFitValues_01
bestFitValues_02, ier_02 = optimize.leastsq(residualsLim, estimate,args=(data,2,15))
print bestFitValues_02, 2+(15-2)*0.5*(1+np.tanh(bestFitValues_02[-1]))
bestFitValues_03, ier_03 = optimize.leastsq(residualsLim, estimate,args=(data,2,8))
print bestFitValues_03, 2+(8-2)*0.5*(1+np.tanh(bestFitValues_03[-1]))
bestFitValues_04, ier_04 = optimize.leastsq(residualsLim, estimate,args=(data,14,24))
print bestFitValues_04, 14+(24-14)*0.5*(1+np.tanh(bestFitValues_04[-1]))
pList=np.linspace(0,2*np.pi,35)
rList_01=np.array([f0(p,*bestFitValues_01) for p in pList])
rList_02=np.array([f_lim(p,2,15,*bestFitValues_02) for p in pList])
rList_03=np.array([f_lim(p,2,8,*bestFitValues_03) for p in pList])
rList_04=np.array([f_lim(p,14,24,*bestFitValues_04) for p in pList])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(data[:,0],data[:,1])
ax.plot(rList_01[:,0],rList_01[:,1])
ax.plot(rList_02[:,0],rList_02[:,1],linestyle="--")
ax.plot(rList_03[:,0],rList_03[:,1],linestyle="--")
ax.plot(rList_04[:,0],rList_04[:,1],linestyle="--")
plt.show()
>> [ 2.05070788 -1.12399476 11.02276442]
>> [ 2.05071109 -1.12399791 0.40958281] 11.0227695567
>> [ 3.32479577e-01 2.14732017e+00 6.60281574e+02] 8.0
>> [ 3.64934819 -4.14597378 -679.24155201] 14.0
如果半径在区间内,则结果与简单适合。否则 x
会变成一个很大的正数或负数,基本上意味着 r
是您指定的限制之一。很好的是,这是连续的。
If the radius is within the interval the result coincides with the simple fit. Otherwise x
becomes a very large positive or negative number, basically meaning that r
is one of your given limits. The nice thing is that this is continuous.
给定代码的结果看起来像
蓝点:随机数据,蓝色图形:简单拟合,黄色图形:在边界内适合 r
,绿色图形:与 r
大于上限,红色图形:适合 r
小于下限。
Results for the given code look like
Blue dots: random data, blue graph: simple fit, yellow graph: fit with r
inside boundaries, green graph: fit with r
larger than upper boundary, red graph: fit with r
smaller than lower boundary.
这篇关于如何将圆拟合到半径受约束的一组点?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!