准备写一个粗差和滤波的长文

以前的一个客户提出了修改要求,趁着五一假期看了看,深深觉得自己之前做的东西问题太多,有点想弥补,折腾了几下又觉得之前走得弯路又要重新走一遍,又想把这些理论的东西整理一下。

大概列了一下提纲,发现内容还是非常多,如果说写这个东西有什么好处的话,大概以下几点:

  1. 顺便宣传公司
  2. 培训客户或内部员工使用
  3. 写得开心

考虑到现在的工作性质,前两点可以忽略不计了,理由只能是为爱发电了。好在缺点也不多,一是费时,一是可能因为知识量不够而虎头蛇尾。先把提纲列好,然后一小节一小节地写吧。

  • 前言
    • 几何处理
    • 测量元素与测量点
    • 元素的拟合
  • 粗大误差
    • 粗差的定义与作用
    • 背景知识
      • 平均值
      • 方差与标准差
    • 基本算法
    • 拟合方法与粗大误差
    • 迭代
    • 过滤额外的点
    • 预滤波
    • 相关标准
  • 滤波
    • 滤波的定义与作用
    • 背景知识
      • 时域与频域
      • 傅里叶变换(分析)
      • 连续信号与离散信号
      • 采样
      • 插值
      • 卷积
      • 传递函数
    • 滤波的类型
      • 按通过频率
      • 按窗口函数
      • 周期与非周期
      • 有限信号与无限信号
    • 基本算法
      • 时域算法
      • 频域算法
    • 拟合方法与滤波
      • 相关标准
      • 其他领域近似概念
  • 计算机实现(C# and/or Python)
    • 相关类库
    • 数据点类
    • 平均值,方差,标准差
    • 过滤粗大误差
    • 泛型支持
    • 插值
    • 傅里叶变换与逆变换
    • 频域滤波算法
    • 时域滤波算法
    • 一个低通窗口的实例
    • 高通、低通、带通、带阻窗口的关系
    • 高斯窗口实现
    • 样条窗口实现
    • RC窗口实现
    • 线滤波实例
    • 圆滤波实例
    • 非整圆滤波的探讨

放弃使用MathNet的Spatial库

这个库也只是更新到0.4,而且更新的非常慢,感觉项目基本停滞。

这次主要的目的之一也是锻炼自己,那么参考代码重新实现也对我自己更有利。

里面的一些结构非常值得参考,一些函数可以直接拿来用,不过这个库还是像通常能见到的库的一样是纯几何的,而不是我想实现的测量上的,所以还是要保留自己的一些思路。

另外就是目前在补回归分析的课,还有解析几何的课(高数下)。最近要花一些时间在这些基础上,代码要少写一些。

仿射变换

实际上坐标测量里研究的变换只有两种,一种是旋转,上一篇文章那样的,另外一种就是平移。

基本坐标系(自然基)

\(O=\begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&1\end{bmatrix}\)

点P

\(\begin{bmatrix}4\\1\\0\\1\end{bmatrix}\)

代表坐标系的矩阵由\(3*3\)变成了\(4*4\),坐标的列向量也由3个元素加了一个1。

虽然介绍仿射变换方法的文章很多,但是原理性的介绍没找到,这种表达形式意味着空间由3维变成了4维,虽然第四维的坐标永远是1,如果有哪位朋友有相关文章提供阅读,感激不尽。

先将坐标系延基本坐标系平移\([0,-3,0]\),再延Z轴旋转45度,在新坐标系

\(A=\begin{bmatrix}\frac{\sqrt{2}}{2}&-\frac{\sqrt{2}}{2}&0&0\\\frac{\sqrt{2}}{2}&\frac{\sqrt{2}}{2}&0&-3\\0&0&1&0\\0&0&0&1\end{bmatrix}\)

下的坐标为点

\(P_{A}=\begin{bmatrix}4\sqrt{2}\\0\\0\\1\end{bmatrix}\)

满足

\[OP=AP_{A}\]

新坐标系的过渡矩阵(相对于自然基)和之前的相比,加了一列(基的个数+1),是新坐标系0点在基本坐标系的坐标值。

将基本坐标系平移\([-10,-10,-10]\)再延Y轴旋转90度。

\(B=\begin{bmatrix}0&0&1&-10\\0&1&0&-10\\-1&0&0&-10\\0&0&0&1\end{bmatrix}\)

在B坐标系下P点坐标

\(P_{B}=B^{-1}P=\begin{bmatrix}-10\\11\\14\\1\end{bmatrix}\)

按向量空间的定义,由A到B的过渡矩阵为\(A^{-1}B\)

\(\begin{bmatrix}0&\frac{\sqrt{2}}{2}&\frac{\sqrt{2}}{2}&-12.0208\\0&\frac{\sqrt{2}}{2}&-\frac{\sqrt{2}}{2}&2.12132\\-1&0&0&-10\\0&0&0&1\end{bmatrix}\)

 

B到A的过渡矩阵\(Q^{-1}\)为

\(\begin{bmatrix}0&0&-1&-10\\\frac{\sqrt{2}}{2}&\frac{\sqrt{2}}{2}&0&7\\\frac{\sqrt{2}}{2}&-\frac{\sqrt{2}}{2}&0&10\\0&0&0&1\end{bmatrix}\)

过渡矩阵的最后一行依旧是\(\begin{bmatrix}0,0,0,1\end{bmatrix}\)

这样点的变换在平移和旋转坐标系的数学关系就整理好了。

在向量空间里,所有的向量表示的都是向量(既有方向,又有长度-大小),这和空间上的点坐标是类似的,但是在测量或者建模时,点还有一个特性就是点的法向(矢量),这个参数只有方向,没有大小(总是单位长度1)。

当坐标系平移时,这个法向是不会随之平移的,只能当坐标系旋转时,法向才会随之旋转。

所以对于法向,坐标的最后一个值永远为0,这样在矩阵变换时,就不会受最后基的最后一个列向量影响了。

算是写完,消耗时间大约是2个晚上加周六周日的业余时间,当然周六周日我都有早上打篮球,周六还去理了发。之所以要把这些很简单的东西整理一遍,而不是把网上现成的材料拷贝保留下来,原因有两个。

1.自己整理一遍,把自以为理解的东西讲一遍,可以排查一下自己是否真正理解了,而且理解的比较全面了。

2.2013年我第一次尝试用C#解决实际的问题,在用代码实现的过程中,发现复杂度比我估计得要多(然后我并没有吸取教训,直到现在我预估依然很差)。而且,虽然矩阵的理论很完备,并不意味着算法是统一的,比如DirectX用行向量,OpenGL用列向量,再比如Math.Net里的DenseMatrix.OfArray构造的矩阵,将数组里的元素按列优先的顺序存入数组。这就导致我在没有足够理论储备的情况下只能用实际数据的测试的方法实现了需要的功能。等到2015年我因为其他事情要扩展这部分代码时,我发现我已经无法修改了。

所以这次借写新库的机会,将理论整理清楚,未来写代码时就不会现一会矩阵乘向量,一会向量乘矩阵,最后乱成一团的情况了。

using System;
using MathNet.Numerics.LinearAlgebra.Double;
using System.Diagnostics;

namespace CoordinateTest
{
    class Program
    {
        static void Main(string[] args)
        {
            dimension3();
            dimension4();
        }
        static void dimension3 ()
        {
            DenseMatrix O = DenseMatrix.OfArray(new double[,]
{
                {1,0,0 },
                {0,1,0 },
                {0,0,1 }
});
            double SqrtTwo = Math.Sqrt(2) / 2.0;
            DenseMatrix A = DenseMatrix.OfArray(new double[,]
            {
                {SqrtTwo,-SqrtTwo,0 },
                {SqrtTwo,SqrtTwo,0 },
                {0,0,1 }
            });
            DenseMatrix B = DenseMatrix.OfArray(new double[,]
            {
                {0,0,1 },
                {0,1,0 },
                {-1,0,0 }
            });

            DenseVector P = DenseVector.OfArray(new double[] { 1, 1, 0 });
            DenseVector PA = A.Inverse() * P as DenseVector;
            DenseVector PB = B.Inverse() * P as DenseVector;

            Debug.WriteLine(PA);
            Debug.WriteLine(PB);
        }
        static void dimension4()
        {
            DenseMatrix O = DenseMatrix.OfArray(new double[,]
            {
                {1,0,0,0 },
                {0,1,0,0 },
                {0,0,1,0 },
                {0,0,0,1 }
            });
            DenseVector P = DenseVector.OfArray(new double[] { 4, 1, 0, 1 });
            double SqrtTwo = Math.Sqrt(2) / 2.0;

            DenseMatrix A = DenseMatrix.OfArray(new double[,]
            {
                {SqrtTwo,-SqrtTwo,0,0 },
                {SqrtTwo,SqrtTwo,0,-3 },
                {0,0,1,0 },
                {0,0,0,1 }
            });
            DenseVector PA = DenseVector.OfArray(new double[] { 8* SqrtTwo, 0, 0, 1 });

            Debug.WriteLine(A * PA);
            DenseMatrix B = DenseMatrix.OfArray(new double[,]
            {
                {0,0,1,-10 },
                {0,1,0,-10 },
                {-1,0,0,-10 },
                {0,0,0,1 }
            });
            Debug.WriteLine(B.Inverse());

            DenseVector PB = B.Inverse() * P as DenseVector;
            Debug.WriteLine(PB);

            DenseMatrix Q = A.Inverse() * B as DenseMatrix;
            Debug.WriteLine(Q);

            Debug.WriteLine(Q * PB);
            Debug.WriteLine(A * Q);

            Debug.WriteLine(Q.Inverse());

            //for vector test
            DenseVector v1 = DenseVector.OfArray(new double[] { SqrtTwo, SqrtTwo, 0, 0 });

            Debug.WriteLine(A.Inverse() * v1);
        }
    }
}

 

输出为:

DenseVector 3-Double
1.41421
      0
      0

DenseVector 3-Double
0
1
1

DenseVector 4-Double
4
1
0
1

DenseMatrix 4x4-Double
0  0  -1  -10
0  1   0   10
1  0   0   10
0  0   0    1

DenseVector 4-Double
-10
 11
 14
  1

DenseMatrix 4x4-Double
 0  0.707107   0.707107  -12.0208
 0  0.707107  -0.707107   2.12132
-1         0          0       -10
 0         0          0         1

DenseVector 4-Double
    5.65685
4.44089E-16
          0
          1

DenseMatrix 4x4-Double
 0  0  1  -10
 0  1  0  -10
-1  0  0  -10
 0  0  0    1

DenseMatrix 4x4-Double
       0          0  -1  -10
0.707107   0.707107   0    7
0.707107  -0.707107   0   10
       0          0   0    1

DenseVector 4-Double
1
0
0
0

 

 

向量空间的基本算法验证

主要是基变换与坐标变换。

先不考虑坐标平移(仿射)的情况,只考虑坐标系旋转(0点不变的情况)。

笛卡尔坐标系是三个轴正交的坐标系,在向量空间里,其实就是维数为3的向量空间的一个基,且3个向量正交且都是单位向量。

假设存在一个基本坐标系(自然基),那么三个轴分别是

\(X_{O}=\begin{bmatrix}1\\0\\0\end{bmatrix}Y_{O}=\begin{bmatrix}0\\1\\0\end{bmatrix}Z_{O}=\begin{bmatrix}0\\0\\1\end{bmatrix}\)

对应的坐标系矩阵为

\(O=\begin{bmatrix}1&0&0\\0&1&0\\0&0&1\end{bmatrix}\)

将坐标系延Z轴旋转45度,对应的三个轴变为

\(X_{A}=\begin{bmatrix}\frac{\sqrt{2}}{2}\\\frac{\sqrt{2}}{2}\\0\end{bmatrix}Y_{A}=\begin{bmatrix}-\frac{\sqrt{2}}{2}\\\frac{\sqrt{2}}{2}\\0\end{bmatrix}Z_{A}=\begin{bmatrix}0\\0\\1\end{bmatrix}\)

得到坐标系

\(A=\begin{bmatrix}\frac{\sqrt{2}}{2}&-\frac{\sqrt{2}}{2}&0\\\frac{\sqrt{2}}{2}&\frac{\sqrt{2}}{2}&0\\0&0&1\end{bmatrix}\)

将坐标系延Y轴旋转90度,对应的三个轴变为

\(X_{B}=\begin{bmatrix}0\\0\\-1\end{bmatrix}Y_{B}=\begin{bmatrix}0\\1\\0\end{bmatrix}Z_{B}=\begin{bmatrix}1\\0\\0\end{bmatrix}\)

得到坐标系

\(B=\begin{bmatrix}0&0&1\\0&1&0\\-1&0&0\end{bmatrix}\)

在C#里用Math.Net实现

            DenseMatrix O = DenseMatrix.OfArray(new double[,]
            {
                {1,0,0 },
                {0,1,0 },
                {0,0,1 }
            });
            double SqrtTwo = Math.Sqrt(2) / 2.0;
            DenseMatrix A = DenseMatrix.OfArray(new double[,]
            {
                {SqrtTwo,-SqrtTwo,0 },
                {SqrtTwo,SqrtTwo,0 },
                {0,0,1 }
            });
            DenseMatrix B = DenseMatrix.OfArray(new double[,]
            {
                {0,0,-1 },
                {0,1,0 },
                {1,0,0 }
            });

定义一个点(向量),即整个向量空间集合里的一个元素,这个点P在基本坐标系下的坐标值为\(\left(1,1,0\right)\),即存在一组常数

\[\begin{cases}\lambda_{O1}=1\\ \lambda_{O2}=1\\ \lambda_{O3}=0\end{cases}\]

满足

\(\begin{align*}P&=\lambda_{O1}X_{O}+\lambda_{O2}Y_{O}+\lambda_{O3}Z_{O}\\&=1*\begin{bmatrix}1\\0\\0\end{bmatrix}+1*\begin{bmatrix}0\\1\\0\end{bmatrix}+0*\begin{bmatrix}0\\0\\1\end{bmatrix}\end{align*}\)

这几个变量满足关系

\[OP=AP_{A}=BP_{B}\]

那么P在A和B基上的坐标为:

\[P_{A}=A^{-1}P=\begin{bmatrix}\sqrt{2}\\0\\0\end{bmatrix}\]

\[P_{B}=B^{-1}P=\begin{bmatrix}0\\1\\1\end{bmatrix}\]

代码实现

            DenseVector P = DenseVector.OfArray(new double[] { 1, 1, 0 });
            DenseVector PA = A.Inverse() * P as DenseVector;
            DenseVector PB = B.Inverse() * P as DenseVector;

            Debug.WriteLine(PA);    DenseVector 3-Double 1.41421 0 0
            Debug.WriteLine(PB);    DenseVector 3-Double 0 1 1

过渡矩阵

\(B=AA^{-1}B\)

另\(Q=A^{-1}B\)

则\(B=AQ\)

过渡矩阵(基变换),实际是上两个坐标系之间的变化关系。

又因为

\[AP_{A}=BP_{B}\]

所以\[P_{A}=A^{-1}BP_{B}\]

\[P_{A}=QP_{B}\]

\[P_{B}=Q^{-1}P_{A}\]

值得注意的是过渡矩阵在坐标系变换和坐标变换起的作用是相反的。这也好理解,一个是坐标系转,一个是点在转,所以方向是相反的。

重温线性代数

之前决定动手写一个库,每当有空想动手写的时候却又卡住不能动,理论知识太欠缺,只能从基础补起,在京东上买了同济的一套数学教材,4本包括线性代数,高数上下,概率。其中线性代数是我最关注的,空间上的坐标变换在计算机实现都是靠矩阵。

花了两个晚上,草草看完,只挑对自己有用的部分做一下记录。

第一章是矩阵定义,性质,结合方程组的介绍,矩阵的变换。第二章是方阵的行列式(determinant),定义,性质,展开,矩阵求逆。

记得当初学的时候是先学的行列式,然后是矩阵,从数学史的发展,也是先有的行列式,然后才有矩阵。不过先介绍哪个都一样,教材都把行列式和矩阵结合线性方程组来介绍,完全不直观。

这两章的内容对我来说意义已经不大了,目标就是计算机实现,很多定义和算法都集成得很好了,由于实际问题都是有解的,所以很多为了验证有解而需要检查的步骤也不需要做了。

矩阵运算的代码,C# with Math.Net

            DenseMatrix A;
            DenseMatrix B;
            A = DenseMatrix.OfArray(
                new double[,]{
                { 3,-1,2},
                { 2,1,-2}
                }
                );
            B = DenseMatrix.OfArray(
                new double[,]{
                { 1,5,1},
                { -2,-1,0}
                }
                );
            Console.WriteLine((A + 2 * B).ToString());
            Console.WriteLine((3 * A - B).ToString());
            Console.WriteLine((A * B.Transpose()).ToString());
            Console.WriteLine((A.Transpose() * B));
            DenseMatrix A = DenseMatrix.OfArray(new double[,]
            {
                {1,2,3 },
                {2,3,1 },
                {3,1,2 }
            });
            Console.WriteLine(A.Inverse().ToString());

            A = DenseMatrix.OfArray(new double[,]
            {
                {1,-1,-1 },
                {-2,2,1 },
                {2,-1,3 }
            });
            Console.WriteLine(A.Solve(DenseVector.OfArray(new double[] { -1, 1, 1 })));

第三章和第五章是向量空间与线性空间,我实在不理解为什么国内几乎所有的教材都无时无刻不把线性方程组拿出来,看起来很容易对应上,但是并不能激发人的好奇啊,学习线性代数就为了给方程组求解?真不如直接跟解析几何对应上,直观,而且也能给学生一些期望值,知道学了以后能干什么。

这两章我并没有搞清区别,向量空间是线性空间的特殊子集,但是对我来说,几何三维空间用哪个并没有区别。

[table “1” not found /]

下一步的工作,用代码实现一些坐标和坐标系的实现。

三坐标转台算法的思考

明天要去一家客户解决转台问题,当其他同事搞了四个月还没有搞好,我也没有信心去一次就能搞定。所以我今天就在做一些准备,就想从原理上先解决一下转台的算法,这样到了现场就比较容易地分析问题出在哪里。

作为一个附加的设备,控制柜中并没有能直接读出坐标值的功能,读到的值依旧是XYZ三个轴上的读数头的示值。而转台唯一能提供的就是当前旋转的角度。

而校准转台时又可以得到转台相对于机器坐标系的轴线坐标\[\begin{cases}x=x_0+it\\y=y_0+jt\\z=z_0+kt\end{cases}(t\in R)\],那么可以得到当转台旋转角度为\(0^\circ\)时的坐标系仿射矩阵为(抱歉这里打成转置矩阵了,我用的算法库里行列就是反的╮(╯_╰)╭)\[A=\begin{bmatrix}\frac{k}{\sqrt{1-j^2}}&0&\frac{-i}{\sqrt{1-j^2}}&-x\\\frac{-i*j}{\sqrt{1-j^2}}&\sqrt{1-j^2}&\frac{-j*k}{\sqrt{1-j^2}}&-y\\i&j&k&-z\\0&0&0&1\end{bmatrix}^T\]。

当转台旋转了\(\theta^\circ\)时,坐标系为原坐标系乘以一个旋转矩阵\[A’=A*\begin{bmatrix}1&0&0&0\\0&\cos{\theta}&-\sin{\theta}&0\\0&\sin{\theta}&\cos{\theta}&0\\0&0&0&1\end{bmatrix}\]

因为三坐标取点是离散的,那么每个角度采的点只是\(\theta\)不同而已,那么就可以通过\(A’\)矩阵将每个点转换到转台的坐标系上。

使用Python模拟高斯滤波过程并生成视频

Python在这类工作上的方便程度真的是——爽。

这里模拟的是低通高斯滤波的过程。

代码如下:需要numpy 和matplotlib,如果想保存成视频文件还需要下载ffmpeg

# encoding: utf-8
'''
Created on 2014年5月1日

@author: zcxsun
'''
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
plt.rcParams['animation.ffmpeg_path'] = r'C:\sunxin\ffmpeg_x64\bin\ffmpeg.exe'

# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure()

speed = 1
size = 360
a = np.arange(size) *np.pi/180
b = np.zeros(size)
x = -0.1
for i in range(1,14,2):
    x = x*(i-2)/i
    b += x * np.sin(a*i)
b+= 0.02 * np.sin(a*50) + 0.01 * np.sin(a*150)
fix_signal = b

hamonic = 15

len = max(1, size/hamonic/2)
print len
scale = 0.3*len+0.8;
scale2x = -0.5/scale/scale
win =  np.exp( pow(np.arange(2*len+1)-len,2)*scale2x)
win/=sum(win)
move_gauss = win
print win.size
filtered_data = np.convolve(b, win, 'same')


#set the axis arrange
y_max = 0.4
ax = plt.axes(xlim=(0, 360), ylim=(-0.2, y_max))
line_gauss, = ax.plot([], [], lw=2)
line_filtered, = ax.plot([],[],lw=2)

plt.plot(fix_signal)

# initialization function: plot the background of each frame
def init():
    line_gauss.set_data([],[])
    line_filtered.set_data([],[])
    return line_gauss,line_filtered,
# animation function.  This is called sequentially
def animate(i):
    n = np.floor(i*speed)
    line_gauss.set_data(n+np.arange(0,win.size),win + (y_max-0.03-win.max()))
    line_filtered.set_data(np.arange(0,n),filtered_data[0:n])
    return line_gauss,line_filtered,
anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=int(np.floor(size/speed)), interval=10, blit=True)

FFwriter = animation.FFMpegWriter(fps = 15)
path = r"C:\Users\Public\Documents\filter_show_hamonic_" + str(hamonic) + ".mp4"
anim.save(path, writer = FFwriter)
#
# plt.show()
print 'finish'

 

根据凸轮升程值求轮廓线

通常此类工作由设计人员完成,根据需求计算出合适的轮廓。

但由于国内质量部门地位低下,以及CAD在设计环节上的缺失,经常拿不到CAD Model或者轮廓的数据,只能拿到升程值。
当然另外一个原因就是专门的凸轮检测仪是模拟测量,可以直接得到对应的升程值。

但是对于坐标检测,就需要将升程值转化为轮廓坐标,通过坐标检测出轮廓线,再转化回升程值。

已知基圆半径\(r_0\),滚子半径\(r_r\),对应角度\(n\)的升程值为\(s_n\),从动件与凸轮偏心量为\(e\),则滚子圆心轨迹坐标\(P_n\left(x_n,y_n\right)\)为\[x_n=(r_0+r_r+s_n)\sin(n)+e\cos(n)\]\[y_n=(r_0+r_r+s_n)\cos(n)+e\sin(n)\]

凸轮轮廓是滚子圆心轨迹的包络曲线,只要求得对应点的法向方向,沿法向方向减去滚子半径即可。如果知道升程的解析式,直接求导数即得,通常情况下,用相邻点近似求得该点切向。

对应点\(P_n\left(x_n,y_n\right)\)的切向矢量表示为

\[\begin{bmatrix}i_n\\j_n\end{bmatrix}=\vec{\begin{bmatrix}\Delta x\\\Delta y\end{bmatrix}}=\vec{\begin{bmatrix}x_{n+1}-x_n\\y_{n+1}-y_n\end{bmatrix}}\]

压力角\(\theta_n=\arctan (\Delta y/\Delta x)\)

凸轮轮廓线上的点\({P’}\left({x’}_n,{y’}_n\right)\)为\[{x’}_n=x_n\mp r_r\cos\theta_n\]\[{y’}_n=y_n\mp r_r\sin\theta_n\]

最小二乘法拟合一些基本空间元素的C#实现 —— (三)圆柱

同前一篇文章相同,算法主要根据这篇文章,修改了其中一些笔误并微调以适合C#代码。

首先是将圆柱参数化并简化该模型。

圆柱的轴线为一空间直线,由一点\(P\)和一空间向量\(V\)定义。

如果经过坐标变换,新坐标系零点为\(P\)且坐标z轴为\(V\),那么该轴线的未知量就由\(x_0,y_0,z_0,a,b,c\)减少为4个。

参数化

圆柱的轴线表示为经过点\(P\left(x_0,y_0,z_0\right)\)向量为\(V\left(a,b,c\right)\)的直线,半径为R

求初始值

我个人觉得这种求法偏繁琐,不过的确相当接近最后拟合值,迭代次数相当少。

\(a,b,c\)的初始值

对于圆柱面上任一点\(P_i\left(x_i,y_i,z_i\right)\)有\[A^2+B^2+C^2=R^2\]其中\[A=c\left(y_i-y_0\right)-b\left(z_i-z_0\right)\]\[B=a\left(z_i-z_0\right)-c\left(x_i-x_0\right)\]\[C=b\left(x_i-x_0\right)-a\left(y_i-y_0\right)\]

将上面的等式代入展开\[\left[c\left(y_i-y_0\right)-b\left(z_i-z_0\right)\right]^2+\left[a\left(z_i-z_0\right)-c\left(x_i-x_0\right)\right]^2+\left[b\left(x_i-x_0\right)-a\left(y_i-y_0\right)\right]^2=R^2\]这个等式很简单(不得不吐槽,当时看到这里我都疯了,现在手打这些公式更疯了),它等价于\[Ax_i^2+By_i^2+Cz_i^2+Dx_iy_i+Ex_iz_i+Fy_iz_i+Gx_i+Hy_i+Iz_i+J=0\]其中\begin{array}{lclp{12}}
A&=&\left(b^2+c^2\right)\\
B&=&\left(a^2+c^2\right)\\
C&=&\left(a^2+b^2\right)\\
D&=&-2ab\\
E&=&-2ac\\
F&=&-2bc\\
G&=&-2\left(b^2+c^2\right)x_0+2aby_0+2acz_0\\
H&=&2abx_0-2\left(a^2+c^2\right)y_0+2bcz_0\\
I&=&2acx_0+2bcy_0-2\left(a^2+b^2\right)z_0\\
J&=&\left(b^2+c^2\right)x_0^2+\left(a^2+c^2\right)y_0^2+\left(a^2+b^2\right)z_0^2-2bcy_0z_0-2acz_0x_0-2abx_0y_0-R^2
\end{array}

两边同除以\(A\)并将\(x\)二次项移到等式右边得
\[\frac{B}{A}y_i^2+\frac{C}{A}z_i^2+\frac{D}{A}x_iy_i+\frac{E}{A}x_iz_i+\frac{F}{A}y_iz_i+\frac{G}{A}x_i+\frac{G}{A}y_i+\frac{I}{A}z_i+\frac{J}{A}=-x^2\]

这是一个线性方程组
\[\begin{bmatrix}y_1^2 &z_1^2 & x_1y_1 & x_1z_1 & y_1z_1 & x_1 & y_1 & z_1 & 1 \\y_2^2 & z_2^2 & x_2y_2 & x_2z_2 & y_2z_2 & x_2 & y_2 & z_2 & 1 \\\vdots &\vdots &\vdots &\vdots &\vdots &\vdots &\vdots &\vdots &\vdots \\y_n^2 &z_n^2 & x_ny_n & x_nz_n & y_nz_n & x_n & y_n & z_n & 1 \end{bmatrix}\begin{bmatrix}B/A\\C/A\\D/A\\E/A\\F/A\\G/A\\H/A\\I/A\\J/A\end{bmatrix}=\begin{bmatrix}-x_1^2\\-x_2^2\\\vdots\\-x_n^2\end{bmatrix}\]

求解:

令\(B’=B/A,C’=C/A,D’=D/A,E’=E/A,F’=F/A,G’=G/A,H’=H/A,I’=I/A,J’=J/A\)
如果\(|D’|\),\(|E’|\)和\(|F’|\)近似等于0(意味着abc中有两个近似于0,向量与某坐标轴平行),如果\(B’\)近似于1向量为\(\left(0,0,1\right)\),如果\(C’\)近似于1向量为\(\left(0,1,0\right)\)

否则,令\[k=\frac{2}{1+B’+C’}\]则\[A=k,B=kB’,C=kC’,D=kD’,E=kE’,F=kF’,G=kG’,H=kH’,I=kI’,J=kJ’\]如果\(A,B\)近似于1,则\begin{array}{lcl}c’&=&\left(1-C\right)^{1/2}\\a’&=&E/-2c’\\b’&=&F/-2c’\end{array}
如果A近似1,B不近似1则\begin{array}{lcl}b’&=&\left(1-B\right)^{1/2}\\a’&=&D/-2b’\\c’&=&F/-2b’\end{array}
如果A不近似于1\begin{array}{lcl}a’&=&\left(1-A\right)^{1/2}\\b’&=&D/-2a’\\c’&=&E/-2a’\end{array}
将向量\(\left(a’,b’,c’\right)\)正规化得到单位向量\[\left(a,b,c\right)\]

注:近似于这个概念很模糊,由于我使用在机械测量中,通常设备采集原始点的精度为0.001左右,而点云本身的标准差在0.001~1,所以我在程序中将近似于定义为\(\Delta<0.001\)
上文两次提到近似于,第一次是\(DEF\),用来判断向量是否平行于某一轴。第二次是\(A,B\),实际上此时已经排除了向量平行于轴线的可能性,这里是用来判断向量更接近于哪个轴,从而避免近似0 的小数开平方造成的精度误差,所以在算法中我将这部分改为了比大小。
另外就是\(ABC\)三个值有可能略微大于1,用原文算法有可能出现负数开方造成无解。

\(x_0,y_0,z_0\)的初始值

根据\(x_0,y_0,z_0\)与\(G,H,I\)的关系,以及等式\(ax_0+by_0+cz_0=0\)(注:原文中使用原点到轴线的垂足作为起始点,因为个人需要,我将之替换为数据第一点投影到轴线的垂足\(ax_0+by_0+cz_0=ax_1+by_1+cz_1\)得到线性方程组

\[\begin{bmatrix}-2\left(b^2+c^2\right)&2ab&2ac\\2ab&-2\left(a^2+c^2\right)&2bc\\2ac&2bc&-2\left(a^2+b^2\right)\\a&b&c\end{bmatrix}\begin{bmatrix}x_0\\y_0\\z_0\end{bmatrix}=\begin{bmatrix}G\\H\\I\\ax_1+by_1+cz_1\end{bmatrix}\]

注:在实际使用中,这个初始值算法有一些问题,当圆柱的轴线非常接近于\((0,0,1)\)时(很常见,通常为了仪器测量方便都会先找正),该方程组变为了\[\begin{bmatrix}-2&0&0\\0&-2&0\\0&0&0\\0&0&1\end{bmatrix}\begin{bmatrix}x_0\\y_0\\z_0\end{bmatrix}=\begin{bmatrix}G\\H\\I\\z_1\end{bmatrix}\]这个一看就知道\(x_0=-\frac{G}{2},y_0=-\frac{H}{2},z_0=z_1\),但是由于第三行都是0的存在,使得矩阵的任何操作都有问题。

\(r\)的初始值

原文使用了前面\(J\)的等式。\[R^2=(b^2+c^2)x_0^2+(a^2+c^2)y_0^2+(a^2+b^2)z_0^2-2bcy_0z_0-2acx_0z_0-2abx_0y_0-J\]

注:依然是在我的实际使用中,发现J的解在特定情况下会是极大的正数(可能是由于圆柱的轴线长度小于直径的原因),于是又出现了负数开方无解的情况,例如这组数据,于是我用了所有测量点半径平均值来代替。

迭代算法

1.坐标变换

在\((x_0,y_0,z_0)\)处以\((a,b,c)\)为Z轴建立一个坐标系U并将所有实测点的一个副本变换到该坐标系下。

\[\begin{bmatrix}x_i\\y_i\\z_i\end{bmatrix}=\left(\begin{bmatrix}x_i\\y_i\\z_i\end{bmatrix}-\begin{bmatrix}x_0\\y_0\\z_0\end{bmatrix}\right)U\]

其中\[U=\begin{bmatrix}1&0&0\\0&c_1&s_1\\0&-s1&c1\end{bmatrix}\begin{bmatrix}c_2&0&s2\\0&1&0\\-s2&0&c2\end{bmatrix}\]

当\((a,b,c)\)为\((1,0,0)\)时\(s_1=0,c_1=1,s_2=-1,c_2=0\)否则\[\begin{array}{lcl}c_1&=&\frac{c}{\sqrt{b^2+c^2}}\\s_1&=&\frac{-b}{\sqrt{b^2+c^2}}\\c_2&=&\frac{cc_1-b_s1}{\sqrt{a^2+(cc_1-bs_1)^2}}\\s_2&=&\frac{-a}{\sqrt{a^2+(cc_1-bs_1)^2}}\end{array}\]

2.构造Jacobian矩阵

参考直线距离公式,并且将\(a=b=0,c=1,z=0\)代入得到\[\begin{bmatrix}-x_1/r_1&-y_1/r_1&-x_1z_1/r_1&-y_1z_1/r_1&-1\\-x_2/r_2&-y_2/r_2&-x_2z_2/r_2&-y_2z_2/r_2&-1\\\vdots&\vdots&\vdots&\vdots&\vdots\\-x_m/r_m&-y_m/r_m&-x_mz_m/r_m&-y_mz_m/r_m&-1\end{bmatrix}\begin{bmatrix}P_{x0}\\P_{y0}\\P_a\\P_b\\P_r\end{bmatrix}=-1\begin{bmatrix}r_1-r\\r_2-r\\\vdots\\r_m-r\end{bmatrix}\]

3.更新参数

\[\begin{bmatrix}x_0\\y_0\\z_0\end{bmatrix}=\begin{bmatrix}x_0\\y_0\\z_0\end{bmatrix}+\begin{bmatrix}P_{x_0}\\P_{y_0}\\-P_{x_0}P_a-P_{y_0}P_b\end{bmatrix}U^T\]

\[\begin{bmatrix}a\\b\\c\end{bmatrix}=\begin{bmatrix}P_a\\P_b\\1\end{bmatrix}U^T\]

\[r=r+P_r\]

4.确定迭代结束条件

重复1~3,直到满足终止条件。

终止条件3选1:

  • 新参数与旧参数对应相减,平方和小于退出值
  • 求所有点的偏差二乘和,新值与旧值相减,差小于退出值
  • 迭代次数超过退出值

包含直线、平面与圆柱算法的C#实例与说明

最小二乘法拟合一些基本空间元素的C#实现 —— (二)直线与平面

本文讨论并实现的算法主要(99%)基于这篇文章,更多参考请阅读这里

大部分能搜索到的关于直线最小二乘拟合的资料是基于2维空间的\[y=ax+ b\]

而从三维空间的角度上来看,这只是一种特例。本文主要讨论三维空间上的直线拟合。

参数化

三维直线的表示如下。\[\left(x,y,z\right)=\left(x_0,y_0,z_0\right) + t\left(a,b,c\right)\]

其中\(\left(x_0,y_0,z_0\right)\)为直线上任一已知点,\(\left(a,b,c\right)\)为向量,所以所有直线都可以使用一个点坐标及一个向量来表示。

对于测量点\(P_i\left(x_i,y_i,z_i\right)\)与直线的距离为\[d_i=\sqrt{{u_i}^2+{v_i}^2+{w_i}^2}\]

其中:\[u_i=c\left(y_i-y_0\right)-b\left(z_i-z_0\right)\]\[v_i=a\left(z_i-z_0\right)-c\left(x_i-x_0\right)\]\[w_i=b\left(x_i-x_0\right)-a\left(y_i-y_0\right)\]

算法

根据原文:拟合的直线过所有点的平均值,没有证明……

所以未知量其实是3个,构建一个矩阵\[A=\begin{bmatrix}x_1-x_0&y_1-y_0&z_1-z_0\\x_2-x_0&y_2-y_0&z_2-z_0\\\vdots & \vdots & \vdots \\x_i-x_0&y_i-y_0&z_i-z_0\end{bmatrix}\]

将矩阵A做奇异值分解(SVD),令:\[A=U_1SV^T\]

矩阵\(V\)的第一列(Column)即是直线向量\(\left(a,b,c\right)\)。

代码实现:

 public override void Estimate(List<Point3D> datas)
        {
            double sum_x = 0;
            double sum_y = 0;
            double sum_z = 0;
            foreach (Point3D temp in datas)
            {
                sum_x += temp.x;
                sum_y += temp.y;
                sum_z += temp.z;
            }
            sum_x /= datas.Count;
            sum_y /= datas.Count;
            sum_z /= datas.Count;

            DenseMatrix jacobian = new DenseMatrix(datas.Count, 3);
            foreach (Point3D temp in datas)
            {
                Vector<double> gradient = new DenseVector(3);
                gradient[0] = temp.x - sum_x;
                gradient[1] = temp.y - sum_y;
                gradient[2] = temp.z - sum_z;
                jacobian.SetRow(datas.IndexOf(temp), gradient);
            }
            Svd svd = jacobian.Svd(true);
            // get matrix of left singular vectors with first n columns of U
            Matrix<double> U1 = svd.U().SubMatrix(0, datas.Count, 0, 3);
            // get matrix of singular values
            Matrix<double> S = new DiagonalMatrix(3, 3, svd.S().ToArray());
            // get matrix of right singular vectors
            Matrix<double> V = svd.VT().Transpose();

            Vector<double> parameters = new DenseVector(3);
            parameters = V.Column(0);
            x = sum_x;
            y = sum_y;
            z = sum_z;
            i = parameters[0];
            j = parameters[1];
            k = parameters[2];
        }

平面的最小二乘拟合

放在一起的原因是这两种元素的算法是完全一样的,直线拟合最后一步得到的矩阵V第一列是直线的向量,而第三列就是平面的向量。

参数化

平面同样可以用一个点\(\left(x_0,y_0,z_0\right)\)和一个向量\(\left(a,b,c\right)\)来表示,只是垂直与该向量。

平面的公式为\[a\left(x-x_0\right)+b\left(y-y_0\right)+c\left(z-z_0\right)=0\)

对于测量点\(P_i\left(x_i,y_i,z_i\right)\)直接的距离为\[d_i=a\left(x-x_0\right)+b\left(y-y_0\right)+c\left(z-z_0\right)\]

算法

与直线的算法完全相同。

包含直线、平面与圆柱算法的C#实例与说明