困扰了数个星期----关于曲线拟合的问题（求助大家）

silitex

这个曲线拟合的问题困扰了数个星期之久！后来才猛然间想到这里有高人，望求解：
目前，我们的设备产品碰到这样的问题：
一个大屏的感应白板（非红外），报点时包含X及Y的坐标；使用感应笔进行手写绘画，但是所报上来的X及Y都可能有偏差，造成画圆形及波浪线时极为毛糙；领导要求采用驱动上曲线拟合进行校准，使之所画出来的曲线是平滑的。
目前一直没有找到合适的算法，包括国内国外的网站，如果大家有现成的源码，还望提供相应链接！
继续刚才的问题，因为没有找到合适的算法，所以只有从头开始研究，就是怎么拟合一个曲线的问题。找了几种现成的方法，已经把贝塞尔曲线排除了，因为贝塞尔曲线并不是经过理想的点；同时也把GdiPlus里面的DrawCurve给排除了，因为经过仔细测试，那个DrawCurve的算法感觉很简单，3点就定位了一条曲线，是一种分段样条曲线，后面的点对前面所画的点没有影响
另外再找最小二乘法曲线模拟，但一直没有找到C++的源码，也从来没有一个源码能够直接画曲线。
想请教大家有没有更好算法的曲线拟合？数学半吊子，研究了数个星期，闹到都麻木了额！期待有好的想法啊！

wayne

你这个需求叫做 circle fitting 吧，用这个关键词能搜出很多有用信息来，比如这个：
http://people.cas.uab.edu/~mosya/cl/CPPcircle.html

silitex

额，谢谢啊！
我这个需求，很多时候可能不一定画圆，比如画波浪线的时候；还有可能绘画的时候，往往画的是抛物线或者椭圆；所以我感觉应该是曲线拟合，英文叫做Curve fitting，我在你这个网页上找一下Curve fitting看看，能否找到点东西。

dianyancao

我刚开始做这个，提取图像的轮廓，然后用参数方程表示曲线，可用NURBS样条分段近似拟合曲线，设置合适的误差值来拟合
需要检测出曲线的角点，直线部分，其他部分可用样条分段近似，加油吧

dianyancao

又看了楼主的需求，应该是将用户提供的图形的图像，拟合成某个参数曲线，
比如拟合曲线$y_i=a \sin(k x_i+b)$，用户给出其图像，需要求解出该曲线的参数$a,k,b$
为此将该图像应用离散傅立叶变换，取出幅值系数最大的分量作为拟合结果可以作为一个解决方案

shikang999

造成画圆形及波浪线时极为毛糙

1、根据你说的用户画圆形你需要识别的问题，对于圆形或者椭圆或者二次抛物线等特殊曲线，确实属于拟合问题，其实你稍微变换一下相应公式，就是求解一个超定线性方程组，可采用极小最小二乘解求解，很简单，这比纯粹的非线性拟合快很多

2、根据你说的画波浪线，要使波浪线平滑，这个属于插值问题（不一定是拟合，当然拟合和插值我个人认为差不多），你可看下人家那些编程语言内置的画曲线的算法，这里你就自己搜索了。不过我建议可以试一下三次样条插值，毕竟几乎每本数值分析书讲到关键部分都会提到这样条插值！

希望能帮到你！

dianyancao

最小二乘法拟合直线
直线的一般方程可表示为\(ax+by+c=0\)，考虑直线上的一个点\(\textbf{z}_k\)，采用齐次坐标，
将点\(\textbf{z}_k\)表示为行向量\(\textbf{A}_k=(x_k,y_k,1)\),直线表示为列向量\(\textbf{z}=(a,b,c)^T\),则：
\[\textbf{A}_k\textbf{z}=0\]
并限定\(a^2+b^2+c^2=1\)，防止\(a,b,c\)同时为\(0\)，写成矩阵形式为：
\[\textbf{A}\textbf{z}=\textbf{0},||\textbf{A}||_2=1\]
当需要用上式拟合\(n\geq3\)个点时，可用最小二乘法拟合，即是让每个点到该直线的误差的平方和\(\epsilon(\textbf{z})\)最小化
\(\epsilon(\textbf{z})\)可表示为\(\epsilon(\textbf{z})=\sum_{k=1}^{n}(\textbf{A}_k\textbf{z})^2\)
对\(\epsilon(\textbf{z})\)求对于列向量\(\textbf{z}\)中每个分量\(\textbf{z}_i\)的偏导数，并令其等于\(0\)，可求得\(\epsilon(\textbf{z})\)的最小值，各偏导数为：
\[\frac{\partial \epsilon(\textbf{z})}{\partial \textbf{z}_{i}}=\sum_{k=1}^{n}2(\textbf{A}_k\textbf{z})(\textbf{A}_k)_i=0\]
将上式写成矩阵形式为：
\[\textbf{A}^T\textbf{A}\textbf{z}=\textbf{0}\]
可直接用高斯消元法求得该线性齐次方程组的非\(\textbf{0}\)解\(\textbf{z}\)，\(\textbf{z}\)在\(\textbf{A}^T\textbf{A}\)的秩为\(2\)时是唯一的，
求得\(\textbf{z}\)后将其归一化为\(\textbf{z}/||\textbf{z}||_2^2\)即解出直线的参数\((a,b,c)^T\)
为将该直线参数化，设\(\textbf{z}(t)=\textbf{M}(t,1,0)^T\)，其中\(\textbf{M}\)是\(3*3\)矩阵，\(\textbf{M}\)排除左上角\(2*2\)子矩阵外的部分为\(0\)，那么
\[\textbf{A}\textbf{z}(t)=\textbf{A}\textbf{M}(t,1,0)^T=\textbf{0}\]
上式的取值不依赖于\(t\)，因此
\[\textbf{A}\textbf{M}=\textbf{0}\]
限定\(\textbf{M}\)的左上角\(2*2\)子矩阵\(||\textbf{M}_{2*2}||_2=1\)，用高斯消元法求得\(\textbf{M}_{2*2}\)的一个非\(\textbf{0}\)解并将其归一化后为\(\textbf{L}=\textbf{M}_{2*2}/{||\textbf{M}_{2*2}||_2^2}\)，
即得到拟合出的直线为\(\textbf{z}(t)=\textbf{L}(t,1)\)，写成参数方程形式为：
\[x(t)=\textbf{L}_{11}t+\textbf{L}_{12} \\
y(t)=\textbf{L}_{21}t+\textbf{L}_{22}\]

下面考虑拟合二次曲线，二次曲线的一般方程可表示为
\[a_{00}+a_{01}y+a_{02}y^2+a_{10}x+a_{11}xy+a_{20}x^2=0\]
有相应的二次型矩阵和上式对应，为：
\[\textbf{z}^T\textbf{A}\textbf{z}=\textbf{0} \\
||\textbf{A}||_2=1\]
其中\(\textbf{z}\)的第\(k\)列\(\textbf{z}_k=(x_k,y_k,1)^T\)，\(\textbf{A}\)为对称矩阵，和二次曲线参数对应
用最小二乘法拟合该曲线，使得误差的平方和\(\epsilon(\textbf{A})\)最小化
\[\epsilon(\textbf{A})=\sum_{k=1}^{n}({\textbf{z}_k}^T\textbf{A}\textbf{z}_k)^2\]
令偏导数\(\frac{\partial \epsilon(\textbf{A})}{\partial \textbf{A}_{ij}}=0\)，其中\(i\leq j\)，得：
\[\frac{\partial \epsilon(\textbf{A})}{\partial \textbf{A}_{ij}}=\sum_{k=1}^{n}({\textbf{z}_k}^T\textbf{A}\textbf{z}_k)^2=\sum_{k=1}^{n}2({\textbf{z}_k}^T\textbf{A}\textbf{z}_k)(\textbf{z}_k)_i(\textbf{z}_k)_j=0\]
一共有\(6\)个偏导数构成线性齐次方程组，可求得其非\(\textbf{0}\)解\(\textbf{A}\)，那么拟合出的二次曲线一般方程为：
\[\textbf{z}^T\textbf{A}\textbf{z}=\textbf{0}\]
将对称矩阵\(\textbf{A}\)对角化，使得\(\textbf{A}=\textbf{Q}^T\varLambda \textbf{Q}\)，其中\(\textbf{Q}\)为正交矩阵，\(\varLambda\)为对角阵，\(\varLambda\)的对角元素是\(\textbf{A}\)的实特征值
令\(\textbf{A}=\textbf{Q}^T\varLambda^{\frac{1}{2}}\varLambda^{\frac{1}{2}}\textbf{Q}=\textbf{C}^T\textbf{C}\)，那么
\[\textbf{z}^T\textbf{A}\textbf{z}=\textbf{z}^T\textbf{C}^T\textbf{C}\textbf{z}=(\textbf{C}\textbf{z})^T\textbf{C}\textbf{z}=\textbf{0}\]
因此，\(||\textbf{C}\textbf{z}||_2=0\)，即
\[\textbf{C}\textbf{z}=\textbf{0}\]
为将该二次曲线参数化，设\(\textbf{z}(t)=\textbf{M}(t^2,t,1)^T\)，其中\(\textbf{M}\)是\(3*3\)仿射矩阵，那么有：
\(\textbf{C}\textbf{M}(t^2,t,1)^T=\textbf{0}\)，该式的取值和\(t\)无关，因此
\[\textbf{C}\textbf{M}=\textbf{0}\]
求得该方程关于\(\textbf{M}\)中仿射参数的非\(\textbf{0}\)解并归一化得到\(\textbf{L}_{2*3}=\textbf{M}_{2*3}/||\textbf{M}_{2*3}||_2^2\)
那么拟合出的二次曲线的参数方程为：
\[x(t)=\textbf{L}_{11}t^2+\textbf{L}_{12}t+\textbf{L}_{13} \\
y(t)=\textbf{L}_{21}t^2+\textbf{L}_{22}t+\textbf{L}_{23}\]

对于平面上的\(n\)次多项式曲线，其一般方程可表示为：
\[\sum_{p+q\leq n}c_{pq}x^py^q=0\]
提供了一组\((x_k,y_k,1)\)的值组成向量组\(\textbf{A}\)后，就给出了关于系数\(c_{pq}\)的线性方程组
可以通过高斯消元法求解系数\(c_{pq}\)，当\(\textbf{A}^T\textbf{A}\)的秩为\(n-1\)时系数\(c_{pq}\)的非\(\textbf{0}\)解是唯一的


我刚系统的了解了下最小二乘法，来自这本书《计算机视觉中的数学方法 吴福朝》[189]页

关于\(n\)次曲线的参数化，是不是总是可以通过关于\(t\)的\(\lceil \frac{(n+2)(n+1)}{4}\rceil=m\)个多项式基\(\{1,t,t^2,...,t^m\}\)来参数化呀？
对于\(n\geq3\)时应该如何参数化用系数\(c_{pq}\)表示的\(n\)次的曲线呀？

dianyancao

上面的参数化方法是错误的，考虑单位圆的上半部分
$$x(t)=t,y(t)=\sqrt{1-t^2}$$
好像不能通过多项式基函数来拟合？
通过增加基函数$\sqrt{1-t^2}$后能线性化求解\(\textbf{M}(t)\)，其中\(\textbf{M}\)是$t$的函数，
对于一般的多项式曲线，应该如何构造基函数来拟合呀？

silitex

前段时间正好有事，没有上线，非常感谢大家！我先消化一下！

silitex

shikang999 发表于 2015-1-10 16:29
造成画圆形及波浪线时极为毛糙

1、根据你说的用户画圆形你需要识别的问题，对于圆形或者椭圆或者二次抛 ...

三次样条插值好像不行，他们的X坐标好像必须递增，限制很多。目前我正尝试x与y都是2次的拟合算法。看能否搞得定。