如何将圆拟合到半径受约束的一组点？ [英] 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 &centre, 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屋！