圆的散点拟合, 已知圆的采样点, 求圆的圆心和半径.

UTF8gbsn

本文介绍一个圆的拟合方法. 参照论文为 Coope, I. D. (1993). Circle fitting
by linear and nonlinear least squares. Journal of Optimization Theory
and Applications, 76(2), 381–388. https://doi.org/10.1007/BF00939613

问题描述

首先还是先来描述一下问题,
假如我们有一个采样数据 A = { a 1 , a 2 , ⋯   , a N } A=\{a_1,a_2, \cdots, a_N\} A={a1​,a2​,⋯,aN​}, 其中$a_i$,
是n维空间中点的坐标向量. 我们希望去拟合一个圆心$x$, 和一个半径$r$,
注意$x$是一个n维向量而$r$是一个标量. 注意重申一下我们的已知集合$A$,
未知$x,r$.

求解过程

首先, 我们需要写出我们的最优化化方程

$$\arg \underset{x,r}{\min }=\sum _{i=1}^N[\Vert x-a_i{\Vert }_2^2-r^2{]}^2$$

接下来我们来分析下中括号内部的这个东西怎么来写.

$$\begin{array}{rl}\Vert x-a_i{\Vert }_2^2-r^2& =(x-a_i{)}^T(x-a_i)-r^2\\ & =(x^T-a_i^T)(x-a_i)-r^2\\ & =x^{T}x-2a_i^{T}x+a_i^T{a}_i-r^2\end{array}$$

对于上式最后的一个等式, 我们可以构造一个成新的形式

$$x^{T}x-2x^T{a}_i+a_i^T{a}_i-r^2=a_i^T{a}_i-\left(\begin{array}{cc}{a}_i^T& 1\end{array}\right)\left(\begin{array}{c}2x\\ r^2-x^{T}x\end{array}\right)$$

接下来令 $$y=\left(\begin{array}{c}2x\\ r^2-x^{T}x\end{array}\right),b_i=\left(\begin{array}{c}{a}_i\\ 1\end{array}\right)$$ 可得
arg ⁡ min ⁡ x , r = ∑ i = 1 N [ ∥ x − a i ∥ 2 2 − r 2 ] 2 ⇒ arg ⁡ min ⁡ y ∑ i = 1 N { a i T a i − b i T y } 2 \arg\min_{x,r}=\sum_{i=1}^{N}[\|x-a_i\|_2^2-r^2]^2 \Rightarrow \arg\min_{y}\sum_{i=1}^{N}\{a_i^Ta_i-b_i^Ty\}^2 argx,rmin​=i=1∑N​[∥x−ai​∥22​−r2]2⇒argymin​i=1∑N​{aiT​ai​−biT​y}2
再次令
B T = { b 1 , b 2 , ⋯   , b N } , d T = { ∥ a 1 ∥ 2 2 , ∥ a 2 ∥ 2 2 , ⋯   , ∥ a N ∥ 2 2 } B^T=\{b_1,b_2, \cdots, b_N\}, d^T=\{\|a_1\|_2^2, \|a_2\|_2^2, \cdots, \|a_N\|_2^2\} BT={b1​,b2​,⋯,bN​},dT={∥a1​∥22​,∥a2​∥22​,⋯,∥aN​∥22​}
可得, 最终的优化方程为 $$\arg \underset{y}{\min }\Vert By-d{\Vert }_2^2$$

至此, 我们实际上就得到了一个非常简单的线性最小二乘法方程.
解得$y=(B^{T}B{)}^{+}{B}^{T}d$之后. 我们就可以很顺利方便的求出$x,r$.($C^{+}$
是$C$的伪逆)

x i = 1 2 y i , i ∈ { 1 , 2 , ⋯   , n } , r = y n + 1 + x T x x_i=\frac{1}{2}y_i, i\in\{1,2,\cdots, n\}, r=\sqrt{y_{n+1}+x^Tx} xi​=21​yi​,i∈{1,2,⋯,n},r=yn+1​+xTx ​

python 代码实现

我们在2D上实现这个算法. 也就是拟合一个平面圆圈.

• 首先, 产生试验数据, 生成一个圆并进行均匀采样,
  采样的时候随机加入噪声.

def genPoints(r, c, n):
    '''
    r: radius
    c: center
    n: sample's number
    '''
    if n <= 10:
        n = 10

    points = []
    for k in range(n):
        angle = 2*np.math.pi * k / n
        x = np.math.cos(angle) * r + np.random.rand()*0.01*r
        y = np.math.sin(angle) * r + np.random.rand()*0.01*r

        p = [c[0]+x, c[1]+y]
        points.append(p)

    return points

• 进行线程方程组求解.

def circleFit(points):

    # generates d
    d = []
    for p in points:
        d.append(np.linalg.norm(p)**2)
    d = np.array(d)

    # generates B
    B = []
    for p in points:
        v = p+[1]
        B.append(v)
    B = np.array(B)

    tB = B.T
    M = tB.dot(B)
    invM = np.linalg.inv(M)

    td = d.T
    sd = td.dot(B)
    y = invM.dot(sd)

    center = [y[0]/2, y[1]/2]
    radius = np.math.sqrt(y[2]+np.linalg.norm(center)**2)

    return center, radius

• 调用代码和测试代码如下

points = genPoints(10, [2,7], 2000)
c, r = circleFit(points)
print(c, r)
# [2.0485439006362554, 7.051224621760529] 10.000426704550165