本文共 1907 字，大约阅读时间需要 6 分钟。

转自：http://www.cnblogs.com/dotLive/archive/2007/04/06/524633.html

最小二乘法(least squares analysis)是一种 数学 优化 技术，它通过 最小化 误差 的平方和找到一组数据的最佳 函数 匹配。 最小二乘法是用最简的方法求得一些绝对不可知的真值，而令误差平方之和为最小。 最小二乘法通常用于 曲线拟合 (least squares fitting) 。这里有 拟合圆曲线 的公式推导过程 和 vc实现。








VC实现的代码：

void

 CViewActionImageTool::LeastSquaresFitting()

{

    if (m_nNum<3)
    {

        return;
    }

    int i=0;

    double X1=0;
    double Y1=0;
    double X2=0;
    double Y2=0;
    double X3=0;
    double Y3=0;
    double X1Y1=0;
    double X1Y2=0;
    double X2Y1=0;

    for (i=0;i<m_nNum;i++)
    {

        X1 = X1 + m_points[i].x;
        Y1 = Y1 + m_points[i].y;
        X2 = X2 + m_points[i].x*m_points[i].x;
        Y2 = Y2 + m_points[i].y*m_points[i].y;
        X3 = X3 + m_points[i].x*m_points[i].x*m_points[i].x;
        Y3 = Y3 + m_points[i].y*m_points[i].y*m_points[i].y;
        X1Y1 = X1Y1 + m_points[i].x*m_points[i].y;
        X1Y2 = X1Y2 + m_points[i].x*m_points[i].y*m_points[i].y;
        X2Y1 = X2Y1 + m_points[i].x*m_points[i].x*m_points[i].y;
    }

    double C,D,E,G,H,N;
    double a,b,c;
    N = m_nNum;
    C = N*X2 - X1*X1;
    D = N*X1Y1 - X1*Y1;
    E = N*X3 + N*X1Y2 - (X2+Y2)*X1;
    G = N*Y2 - Y1*Y1;
    H = N*X2Y1 + N*Y3 - (X2+Y2)*Y1;
    a = (H*D-E*G)/(C*G-D*D);
    b = (H*C-E*D)/(D*D-G*C);
    c = -(a*X1 + b*Y1 + X2 + Y2)/N;

    double A,B,R;
    A = a/(-2);
    B = b/(-2);
    R = sqrt(a*a+b*b-4*c)/2;

    m_fCenterX = A;
    m_fCenterY = B;
    m_fRadius = R;

    return;
}

编译运行后随便打开一个图片，当然最好是全白的图片，然后就点吧，大于三个点后就会开始拟合。红线画的圆为拟合的圆，深蓝的点为鼠标点击设置的样本点。单击鼠标右键清空样本集。

下面是matlab代码： （http://wenku.baidu.com/link?url=mPD8NOhWJUjQwl2rcXms8Mm2dhI5EY402uXztZeCjwfRKx44v6sg10XDJnlRBJU-mIYJbywZ_4lHL0HBiHVoouURf-DWiKxXgPwfM951lXK）

%%此程序用于对二值图像获得的圆形边缘进行圆拟合，计算出圆心坐标及半径 %最小二乘法进行曲线拟合 function [xc,yc,R] = cirfit(x,y)%x、y为坐标点，都是一组向量。如x=[x1,x2,x3,...,xn] n=length(x); xx=x.*x; yy=y.*y; xy=x.*y; A=[sum(x) sum(y) n;sum(xy) sum(yy)  sum(y);sum(xx) sum(xy) sum(x)]; B=[-sum(xx+yy);  -sum(xx.*y+yy.*y);  -sum(xx.*x+xy.*y)]; a=A\B; xc = -0.5*a(1); yc = -0.5*a(2); R = sqrt((a(1)^2+a(2)^2)/4-a(3)); end

%结束后返回的值便是圆心坐标及半径