如何将圆拟合到具有约束半径的一组点？

Question

我有一组点代表一个小圆弧。

当前代码使用线性最小二乘法将圆拟合到这些点：

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。

因此应选择半径在该范围内的最佳拟合圆，而不是整体最佳拟合。

如何修改代码以添加此约束？

Answer 1

为什么不使用HoughCircles()？它允许您指定两个参数。示例取自链接文档，修改为添加最小和最大半径参数：

#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>
#include <math.h>

using namespace cv;

int main(int argc, char** argv)
{
    Mat img, gray;
    img = imread('ring.png')
    cvtColor(img, gray, COLOR_BGR2GRAY);
    // smooth it, otherwise a lot of false circles may be detected
    medianBlur(gray, gray, 7);
    vector<Vec3f> circles;
    int minRadius = 90
    int maxRadius = 100
    HoughCircles(gray, circles, HOUGH_GRADIENT,
                 1, 100, 20, 20, minRadius, maxRadius);
    for( size_t i = 0; i < circles.size(); i++ )
    {
         Point center(cvRound(circles[i][0]), cvRound(circles[i][1]));
         int radius = cvRound(circles[i][2]);
         // draw the circle center
         circle( img, center, 3, Scalar(0,255,0), -1, 8, 0 );
         // draw the circle outline
         circle( img, center, radius, Scalar(0,0,255), 3, 8, 0 );
    }
    namedWindow( "circles", 1 );
    imshow( "circles", img );
    return 0;
}

运行示例（注意，我是用 Python 做的，但参数相同）：

图片从a PhD student at CMU盗取。

Answer 2

使用来自

的坐标方程

https://de.wikipedia.org/wiki/Kreis

     G: a*x + b*y + c + y^2 + x^2 = 0
     with
     a: -2*xm
     b: -2*ym
     c: ym^2+xm^2-r^2

（英文维基百科似乎没有提到这个确切的等式，但提到了基本形式。）

要获得线性回归，请对 a/b/c 区分 G^2

然后你得到矩阵形式

           A        B         C      R        R
     G1:  +2*a*x^2 +2*b*x*y  +2*c*x +2*x^3    +2*x*y^2  = 0 
     G2:  +2*a*x*y +2*b*y^2  +2*c*y +2*y^3    +2*x^2*y  = 0
     G3:  +2*a*x   +2*b*y    +2*c   +2*y^2    +2*x^2    = 0

这个派生方程可以用你的点来求解 a/b/c。您可以从 a/b/c 重新替换以获取 xm/ym/r。要拒绝特定半径，只需检查其限制即可。

例子：

编辑：来源

#include <opencv2/opencv.hpp>

#define _USE_MATH_DEFINES
#include <math.h>

#include <vector>
#include <random>

void fit_circle(const std::vector<cv::Point2d> &pnts, cv::Point2d &centre,
                double &radius)
{
    /*
                   A        B         C      R        R
         G1:  +2*a*x^2 +2*b*x*y  +2*c*x +2*x^3    +2*x*y^2  = 0
         G2:  +2*a*x*y +2*b*y^2  +2*c*y +2*y^3    +2*x^2*y  = 0
         G3:  +2*a*x   +2*b*y    +2*c   +2*y^2    +2*x^2    = 0
    */

    static const int rows = 3;
    static const int cols = 3;
    cv::Mat LHS(rows, cols, CV_64F, 0.0);
    cv::Mat RHS(rows, 1, CV_64F, 0.0);
    cv::Mat solution(rows, 1, CV_64F, 0.0);

    if (pnts.size() < 3)
    {
        throw std::runtime_error("To less points");
    }

    for (int i = 0; i < static_cast<int>(pnts.size()); i++)
    {
        double x1 = pnts[i].x;
        double x2 = std::pow(pnts[i].x, 2);
        double x3 = std::pow(pnts[i].x, 3);
        double y1 = pnts[i].y;
        double y2 = std::pow(pnts[i].y, 2);
        double y3 = std::pow(pnts[i].y, 3);

        // col 0 = A / col 1 = B / col 2 = C
        // Row 0 = G1
        LHS.at<double>(0, 0) += 2 * x2;
        LHS.at<double>(0, 1) += 2 * x1 * y1;
        LHS.at<double>(0, 2) += 2 * x1;

        RHS.at<double>(0, 0) -= 2 * x3 + 2 * x1 * y2;

        // Row 1 = G2
        LHS.at<double>(1, 0) += 2 * x1 * y1;
        LHS.at<double>(1, 1) += 2 * y2;
        LHS.at<double>(1, 2) += 2 * y1;

        RHS.at<double>(1, 0) -= 2 * y3 + 2 * x2 * y1;

        // Row 2 = G3
        LHS.at<double>(2, 0) += 2 * x1;
        LHS.at<double>(2, 1) += 2 * y1;
        LHS.at<double>(2, 2) += 2;

        RHS.at<double>(2, 0) -= 2 * y2 + 2 * x2;
    }

    cv::solve(LHS, RHS, solution);

    std::vector<double> abc{solution.at<double>(0, 0),
                            solution.at<double>(1, 0),
                            solution.at<double>(2, 0)};

    centre.x = abc[0] / -2.0;
    centre.y = abc[1] / -2.0;
    radius = std::sqrt(
        std::abs(std::pow(centre.x, 2) + std::pow(centre.y, 2) - abc[2]));
}

int main(int argc, const char *argv[])
{
    try
    {
        std::random_device rd;
        std::mt19937 gen(rd());
        std::uniform_real_distribution<> dis(-0.5, 0.5);

        // Generate reandom points
        double radius_in = 25;
        double xm_in = 10;
        double ym_in = 20;

        std::vector<cv::Point2d> pnts;
        {
            for (double ang = 0; ang <= 90; ang += 10)
            {
                cv::Point2d p(
                    cos(ang / 180.0 * M_PI) * radius_in + xm_in + dis(gen),
                    sin(ang / 180.0 * M_PI) * radius_in + ym_in + dis(gen));
                pnts.push_back(p);
            }
        }

        cv::Point2d c_out;
        double radius_out = 0;
        fit_circle(pnts, c_out, radius_out);

        double dxm = c_out.x - xm_in;
        double dym = c_out.y - ym_in;
        double dradius = radius_out - radius_in;

        std::cout << "Deltas: " << dxm << " " << dym << " " << dradius
                  << std::endl;
        return EXIT_SUCCESS;
    }
    catch (const std::exception &ex)
    {
        std::cerr << ex.what() << std::endl;
    }

    return EXIT_FAILURE;
}

这输出 增量：0.180365 0.190016 -0.231563

---

编辑：添加源以生成有关算法的统计信息和统计信息

#include <opencv2/opencv.hpp>

#define _USE_MATH_DEFINES
#include <math.h>

#include <vector>
#include <random>

void fit_circle(const std::vector<cv::Point2d> &pnts, cv::Point2d &centre,
                double &radius)
{
    /*
                A        B         C      R        R
        G1:  +2*a*x^2 +2*b*x*y  +2*c*x +2*x^3    +2*x*y^2  = 0
        G2:  +2*a*x*y +2*b*y^2  +2*c*y +2*y^3    +2*x^2*y  = 0
        G3:  +2*a*x   +2*b*y    +2*c   +2*y^2    +2*x^2    = 0
    */

    static const int rows = 3;
    static const int cols = 3;
    cv::Mat LHS(rows, cols, CV_64F, 0.0);
    cv::Mat RHS(rows, 1, CV_64F, 0.0);
    cv::Mat solution(rows, 1, CV_64F, 0.0);

    if (pnts.size() < 3)
    {
        throw std::runtime_error("To less points");
    }

    for (int i = 0; i < static_cast<int>(pnts.size()); i++)
    {
        double x1 = pnts[i].x;
        double x2 = std::pow(pnts[i].x, 2);
        double x3 = std::pow(pnts[i].x, 3);
        double y1 = pnts[i].y;
        double y2 = std::pow(pnts[i].y, 2);
        double y3 = std::pow(pnts[i].y, 3);

        // col 0 = A / col 1 = B / col 2 = C
        // Row 0 = G1
        LHS.at<double>(0, 0) += 2 * x2;
        LHS.at<double>(0, 1) += 2 * x1 * y1;
        LHS.at<double>(0, 2) += 2 * x1;

        RHS.at<double>(0, 0) -= 2 * x3 + 2 * x1 * y2;

        // Row 1 = G2
        LHS.at<double>(1, 0) += 2 * x1 * y1;
        LHS.at<double>(1, 1) += 2 * y2;
        LHS.at<double>(1, 2) += 2 * y1;

        RHS.at<double>(1, 0) -= 2 * y3 + 2 * x2 * y1;

        // Row 2 = G3
        LHS.at<double>(2, 0) += 2 * x1;
        LHS.at<double>(2, 1) += 2 * y1;
        LHS.at<double>(2, 2) += 2;

        RHS.at<double>(2, 0) -= 2 * y2 + 2 * x2;
    }

    cv::solve(LHS, RHS, solution);

    std::vector<double> abc{solution.at<double>(0, 0),
                            solution.at<double>(1, 0),
                            solution.at<double>(2, 0)};

    centre.x = abc[0] / -2.0;
    centre.y = abc[1] / -2.0;
    radius = std::sqrt(
        std::abs(std::pow(centre.x, 2) + std::pow(centre.y, 2) - abc[2]));
}

int main(int argc, const char *argv[])
{
    int count = 100;
    if (argc == 2)
    {
        count = std::atoi(argv[1]);
    }

    try
    {
        std::random_device rd;
        std::mt19937 gen(rd());
        std::uniform_real_distribution<> jitter(-0.5, 0.5);
        std::uniform_real_distribution<> pos(-100, 100);

        std::uniform_real_distribution<> w_start(0, 360);
        std::uniform_real_distribution<> w_length(10, 180);
        std::uniform_real_distribution<> rad(10, 180);

        std::cout << "start\tlen\tdelta xm\tdelta ym\tdelta radius" << std::endl;

        for (int i = 0; i < count; ++i)
        {
            // Generate reandom points
            double radius_in = rad(gen);
            double xm_in = pos(gen);
            double ym_in = pos(gen);

            double start = w_start(gen);
            double len = w_length(gen);

            std::vector<cv::Point2d> pnts;
            {
                for (double ang = start; ang <= start + len; ang += 1)
                {
                    cv::Point2d p(
                        cos(ang / 180.0 * M_PI) * radius_in + xm_in + jitter(gen),
                        sin(ang / 180.0 * M_PI) * radius_in + ym_in + jitter(gen));
                    pnts.push_back(p);
                }
            }

            cv::Point2d c_out;
            double radius_out = 0;
            fit_circle(pnts, c_out, radius_out);

            double dxm = c_out.x - xm_in;
            double dym = c_out.y - ym_in;
            double dradius = radius_out - radius_in;

            std::cout << start << "\t" << len << "\t" << dxm << "\t" << dym << "\t" << dradius
                    << std::endl;
        }

        return EXIT_SUCCESS;
    }
    catch (const std::exception &ex)
    {
        std::cerr << ex.what() << std::endl;
    }

    return EXIT_FAILURE;
}

100000 个圈子

绘图已生成：

set multiplot layout 3,1 rowsfirst
set yrange [-50:50]
set grid 
plot 'stat.txt' using 2:3 title 'segment-length vs delta x' with points ps 0.1
set yrange [-50:50]
set grid
plot 'stat.txt' using 2:4 title 'segment-length vs delta y' with points ps 0.1
set yrange [-50:50]
set grid
plot 'stat.txt' using 2:5 title 'segment-length vs delta radius' with points ps 0.1

Answer 3

我希望我能解决您的问题。第一个想法可能是：使用带有约束的拟合算法。不过，我不太喜欢这些，因为它们非常耗时。我的方法是通过用r=r_min+(r_max-r_min)*0.5*(1+tanh(x)) 类型的函数替换它并拟合x 来限制半径。所以x 可以在整个实数范围内变化，从而为您提供给定范围内的半径。 All 仍然是连续函数的简单非线性拟合。

作为一个 python 示例，它看起来像：

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 是您给定的限制之一。好消息是这是连续的。

给定代码的结果如下所示 蓝点：随机数据，蓝色图表：简单拟合，黄色图表：与r 内部边界拟合，绿色图表：与r 拟合大于上边界，红色图表：与r 拟合小于下边界。