Julia内核项目记录（二）

Julia内核项目

1、圆拟合函数🔥

1.1、线性最小二乘拟合原理🔥

这是最简单的圆拟合方法，圆拟合的方法熟知的还有非线性最小二乘，RANSAC拟合圆，样条函数拟合圆等等。还有基于图像的Hough圆拟合。线性最小二乘拟合其实在应用中足够普适了，为了所谓一点的精度采取非线性最小二乘，我觉得毫无必要，为了所谓的抗噪，完全可以采用RANSAC方法去弥补。

线性最小二乘拟合圆本质就是个数学问题。

圆的方程：$(x - a)^2 + (y - b)^2 = R^2$

拆开该方程，整理可以得到：

$x^2 + y^2 = R^2 +2ax + 2by +(a^2+b^2)$

设$z=x^2+y^2$，则：

$z = R^2 + 2ax + 2by + (a^2+b^2)$

换一种形式写：

$z = 2ax + 2by + C$

x和y就是样本点，用选定的样本点，计算$z=2ax+2by+C$，然后最小化平方误差，确定参数a，b，R

可以使用现成的包装库，但我选择手动实现。 这里需要用到了一个概念叫 正规方程 简单说就是专门解最小二乘的手段，一般表示如下： [ A^T A \mathbf{x} = A^T \mathbf{b} ]

如果 ( A^T A ) 可逆，则最优解为： [ \mathbf{x} = (A^T A)^{-1} A^T \mathbf{b} ]

在平面拟合圆中，点数不太多，一般小于10000点的时候，$A^T A$计算非常快，这会极大加速拟合的过程。

这里$x$就是$(a , b ,c)$ A为下式：

A的转置同理。

线性最小二乘拟合代码

原理清楚，代码就好写了

 @inbounds for i in 1:n
        xi = x[i]
        yi = y[i]
        xi2 = xi * xi
        yi2 = yi * yi
        Sx  += xi
        Sy  += yi
        Sxx += xi2
        Sxy += xi * yi
        Syy += yi2
        zi  = xi2 + yi2
        Sxz += xi * zi
        Syz += yi * zi
        Sz  += zi
    end

对于每一个样本点，依次计算需要用到的结果，累加到相应的变量中。把$(A^T A)^{-1} A^T$展开罢了~

# 使用 StaticArrays 构建正规方程
    AtA = zeros(MMatrix{3,3,T,9})
    Atz = zeros(MVector{3,T})
    
    AtA[1,1] = 4*Sxx
    AtA[1,2] = 4*Sxy
    AtA[1,3] = 2*Sx
    AtA[2,2] = 4*Syy
    AtA[2,3] = 2*Sy
    AtA[3,3] = n  
    
    Atz[1] = 2*Sxz
    Atz[2] = 2*Syz
    Atz[3] = Sz
    
    # 对称填充
    AtA[2,1] = AtA[1,2]
    AtA[3,1] = AtA[1,3]
    AtA[3,2] = AtA[2,3]

往里填就行，没啥好说的，就是别填错了。这里用到了StaticArrays，处理小数组更好，但我实验实际并不会提升特别的多性能，所以用普通数组也无所谓。

3、总结

优化过的代码，处理速度完全可以超过正常处理速度，实测提升了一半的性能。

正常写法

function circle_fit_3(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:Real}
    n = length(x);
    z = x.^2 +y.^2;
    A = [2x 2y fill(one(T), n)];
    coeff = A \ z ;  
    a, b, c = coeff[1], coeff[2], coeff[3];
    R = sqrt(c + a^2 + b^2);
    return a,b,R
end

但是必须认识到，在处理大量点云的时候，可能普通写法更稳妥，它采取了更稳定的解法。但实际应用中，似乎很难出现对密集点云拟合圆的场景，至少我没遇到😃