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

程序员文章站 2022-03-31 20:53:12

...

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, \cdots, a_N\}$ , 其中 $a_i$ ,
是n维空间中点的坐标向量. 我们希望去拟合一个圆心 $x$ , 和一个半径 $r$ ,
注意 $x$ 是一个n维向量而 $r$ 是一个标量. 注意重申一下我们的已知集合 $A$ ,
未知 $x, r$ .

求解过程

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

$\arg\min_{x,r}=\sum_{i=1}^{N}[\|x-a_i\|_2^2-r^2]^2$

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

$\left. \begin{aligned} \|x-a_i\|_2^2-r^2&=(x-a_i)^T(x-a_i)-r^2\\ &=(x^T-a_i^T)(x-a_i)-r^2\\ &=x^Tx-2a_i^Tx+a_i^Ta_i-r^2\\ \end{aligned} \right.$

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

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

接下来令 $y=\left( \begin{array}{c} 2x \\ r^2-x^Tx \end{array} \right), b_i = \left( \begin{array}{c} a_i \\ 1 \end{array} \right)$ 可得
$\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$
再次令
$B^T=\{b_1,b_2, \cdots, b_N\}, d^T=\{\|a_1\|_2^2, \|a_2\|_2^2, \cdots, \|a_N\|_2^2\}$
可得, 最终的优化方程为 $\arg\min_{y}\|By-d\|^2_2$

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

$x_i=\frac{1}{2}y_i, i\in\{1,2,\cdots, n\}, r=\sqrt{y_{n+1}+x^Tx}$

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

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

问题描述

求解过程

python 代码实现

已知圆上两点坐标和半径求圆心

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