贝塞尔曲线最小二乘法拟合（随意切向/切向方向统一）------路适用于绝大多数的最小二乘法拟合

标题@TOC

注意

本文出现大量公式，但是博主保证，只要专心看下去，公式都不难的。

简介

最近参加了一个拟合边缘的项目，本来是刚接手的时候只有直线拟合和2、3次线性拟合，因为实验品是丰富多样的，再加上客户比较挑，两端之间的交界处那边希望看起来圆滑点，
所以以上这三种显然是无法满足的，于是就想到了贝塞尔曲线，最终效果还是蛮好的，因为公式是自己纯手推的，所以想到特地来趁热记录下，免得后天又忘天上去了。先提醒一下，本篇博客公式是比较多，但是仔细看的话公式大部分都是积分，并不涉及其他比较深奥的公式，仔细看就很简单，相信我。
话不多说，先上效果图

上图红色部分为拟合后的效果，并非是一段三次贝塞尔，而是好多段拼合起来的
一段的效果如下

显然不可行。
下图是源数据

只看红色部分哦！

贝塞尔曲线

贝塞尔曲线我就不过多介绍了，不清楚的同学可以问一问度娘，不过熟悉ps的同学应该知道里面的钢笔工具就是贝塞尔，写到这瞬间感觉自己快实现自动抠图了，

飘一下下就行了，自动抠图可谈不上，本次实验全是基于已有源数据的情况下进行的，源数据的获取不在本文范畴。
下面是三次贝塞尔曲线公式【（这里的公式存在一点问题，有一位细心的朋友指正了一下，但是因为后面的公式都是用这个错误的公式推导的，博主也懒得改了，所以各位同学只要学习方法即可，自己动手尝试的时候不要被我的公式误导了，正确的公式是$$B(t)=P_0(1-t{)}^3+3\ast P_{1}t(1-t{)}^2+3\ast P_2{t}^2(1-t)+P_3{t}^3,t\in (0,1).$$）】
$$B(t)=P_0(1-t{)}^3+P_{1}t(1-t{)}^2+P_2{t}^2(1-t)+P_3{t}^3,t\in (0,1).$$
你需要知道的是：
三次贝塞尔曲线包括四个控制点，其中P0和P3分别代表曲线的起始位置和结束位置，P0到P1的方向表示曲线起始位置的切线方向，P2到P3的方向表示曲线结尾位置的切线方向。

最小二乘法拟合

曲线两端不带切向约束的公式推导

最小二乘法就是计算拟合后的位置和源数据的距离最小的时候的参数的值，假设Bix和Biy为拟合后第i个点对应的x，y坐标，Xi和Yi为源数据中第i个点的x，y坐标。N个点。
最小二乘法即是求公式：

$$\sum _{i=1}^N(B_{ix}-X_i{)}^2+(B_{iy}-Y_i{)}^2$$
上述公式最小时各个参数的值即为求得的贝塞尔曲线的系数。
那让我们来思考一下这种情况下贝塞尔曲线需要计算几个系数呢？
首先首尾两个点并非是变量，而是每次拟合前要确定的也就是P0和P3是已知的量（如果未知的话或许也可以，但是我没有那样做，所以不知道最后效果如何，感兴趣的同学可以尝试一下）。P0和P3已知的情况下，曲线公式中就只剩下P1和P2作为未知量，再加上不带切向约束，P1以及P2的x，y坐标不存在线性关系，因此需要求的是四个系数：P1x，P1y，P2x和P2y。
我们假设：

$$t_1=(1-t{)}^3,t_2=t(1-t{)}^2,t_3=t^2(1-t),t_4=t^3.$$
其中 t为已知量，所以上述四个数值皆为已知量。(假设有N个源数据，第i个数据对应的 t 为 (i-1)/(N-1) )
上述需要求解最小值的公式即为

$$\sum _{i=1}^N(P_{0x}{t}_1+P_{1x}{t}_2+P_{2x}{t}_3+P_{3x}{t}_4-X_i{)}^2+(P_{0y}{t}_1+P_{1y}{t}_2+P_{2y}{t}_3+P_{3y}{t}_4-Y_i{)}^2$$
因为上述公式只有P1x，P1y，P2x和P2y为未知量，所以分别对这四个变量求偏导数，四个偏导数结果同时为0时即为所求的结果，此时的贝塞尔曲线就是最接近所给源数据的贝塞尔曲线。
偏导数分别为：

• 求P1x的偏导数为0的公式为
  $$\sum _{i=1}^{N}2\ast (P_{0x}{t}_1+P_{1x}{t}_2+P_{2x}{t}_3+P_{3x}{t}_4-X_i)\ast t_2=0$$
• 求P2x的偏导数为0的公式为
  $$\sum _{i=1}^{N}2\ast (P_{0x}{t}_1+P_{1x}{t}_2+P_{2x}{t}_3+P_{3x}{t}_4-X_i)\ast t_3=0$$
• 求P1y的偏导数为0的公式为
  $$\sum _{i=1}^{N}2\ast (P_{0y}{t}_1+P_{1y}{t}_2+P_{2y}{t}_3+P_{3y}{t}_4-Y_i)\ast t_2=0$$
• 求P2y的偏导数为0的公式为
  $$\sum _{i=1}^{N}2\ast (P_{0y}{t}_1+P_{1y}{t}_2+P_{2y}{t}_3+P_{3y}{t}_4-Y_i)\ast t_3=0$$

上述四个公式整理之后为

$$\sum _{i=1}^N(P_{1x}{t}_2+P_{2x}{t}_3)\ast t_2=\sum _{i=1}^N(X_i-P_{0x}{t}_1+P_{3x}{t}_4)\ast t_2$$

$$\sum _{i=1}^N(P_{1x}{t}_2+P_{2x}{t}_3)\ast t_3=\sum _{i=1}^N(X_i-P_{0x}{t}_1+P_{3x}{t}_4)\ast t_3$$

$$\sum _{i=1}^N(P_{1y}{t}_2+P_{2y}{t}_3)\ast t_2=\sum _{i=1}^N(Y_i-P_{0y}{t}_1+P_{3y}{t}_4)\ast t_2$$

$$\sum _{i=1}^N(P_{1y}{t}_2+P_{2y}{t}_3)\ast t_3=\sum _{i=1}^N(Y_i-P_{0y}{t}_1+P_{3y}{t}_4)\ast t_3$$
继续整理，方便后续化为矩阵理解计算

$$\sum _{i=1}^N{t}_2^2{P}_{1x}+t_2{t}_3{P}_{2x}+0P_{1y}+0P_{2y}=\sum _{i=1}^N(X_i-P_{0x}{t}_1+P_{3x}{t}_4)\ast t_2$$

$$\sum _{i=1}^N{t}_2{t}_3{P}_{1x}+t_3^2{P}_{2x}+0P_{1y}+0P_{2y}=\sum _{i=1}^N(X_i-P_{0x}{t}_1+P_{3x}{t}_4)\ast t_3$$

$$\sum _{i=1}^{N}0P_{1x}+0P_{2x}+t_2^2{P}_{1y}+t_2{t}_3{P}_{2y}=\sum _{i=1}^N(Y_i-P_{0y}{t}_1+P_{3y}{t}_4)\ast t_2$$

$$\sum _{i=1}^{N}0P_{1x}+0P_{2x}+t_2{t}_3{P}_{1y}+t_3^2{P}_{2y}=\sum _{i=1}^N(Y_i-P_{0y}{t}_1+P_{3y}{t}_4)\ast t_3$$

细心的同学可能就发现了，上面四个公式完全可以转换成一个矩阵。

$$\left[\begin{array}{cccc}{\sum }_{i=1}^N{t}_2^2& {\sum }_{i=1}^N{t}_2{t}_3& 0& 0\\ & & & \\ {\sum }_{i=1}^N{t}_2{t}_3& {\sum }_{i=1}^N{t}_3^2& 0& 0\\ & & & \\ 0& 0& {\sum }_{i=1}^N{t}_2^2& {\sum }_{i=1}^N{t}_2{t}_3\\ & & & \\ 0& 0& {\sum }_{i=1}^N{t}_2{t}_3& {\sum }_{i=1}^N{t}_3^2\end{array}\right]\left[\begin{array}{c}{P}_{1x}\\ \\ P_{2x}\\ \\ P_{1y}\\ \\ P_{2y}\end{array}\right]=\left[\begin{array}{c}{\sum }_{i=1}^N(X_i-P_{0x}{t}_1+P_{3x}{t}_4)\ast t_2\\ \\ {\sum }_{i=1}^N(X_i-P_{0x}{t}_1+P_{3x}{t}_4)\ast t_3\\ \\ {\sum }_{i=1}^N(Y_i-P_{0y}{t}_1+P_{3y}{t}_4)\ast t_2\\ \\ {\sum }_{i=1}^N(Y_i-P_{0y}{t}_1+P_{3y}{t}_4)\ast t_3\end{array}\right]$$

类似于

$$AX=B$$

其中X为要求的四个变量组成的矩阵，A和B可以根据已有数据计算得出，然后可以利用opencv的solve函数估计得出X的结果，从而得出四个变量的值。
至此，两端不带切向约束的贝塞尔曲线拟合完成。

曲线一端带切向约束的公式推导

前端切向约束：
一端带切向约束的贝塞尔曲线拟合的意思是曲线某一端规定好切线方向，在此种情况下进行拟合求解，它包括前端约束和后端约束。这里只介绍前端约束即可，后端约束基本一致。

同样是求下面公式最小值时各非确定系数的值

$$\sum _{i=1}^N(B_{ix}-X_i{)}^2+(B_{iy}-Y_i{)}^2$$

一端带切向约束的情况下，有几个系数是非确定的呢？

答案是三个。

很简单，和不带约束的情况一样，头尾控制点是已知确定的，只需要求取中间两个控制点。同时又因为前端的切线方向已知，因此只需要一个系数k即可求得第二个控制点，第三个控制点的x和y不存在什么线性关系，因此算作两个变量；综上所述，需要计算三个变量即可。
假设为k，P2x和P2y;
其中前端切向已知，就可以求出单位切向量vec，这样第二个控制点可以表示为

$$P_1=P_0+k\ast vec$$

所以上述最小二乘式中的相关变量可以表示为

$$B_{ix}=P_{0x}{t}_1+(P_{0x}+k\ast vec.x)t_2+P_{2x}{t}_3+P_{3x}{t}_4$$

$$B_{iy}=P_{0y}{t}_1+(P_{0y}+k\ast vec.y)t_2+P_{2y}{t}_3+P_{3y}{t}_4$$

最小二乘式即可替换为

$$\sum _{i=1}^N(P_{0x}{t}_1+(P_{0x}+k\ast vec.x)t_2+P_{2x}{t}_3+P_{3x}{t}_4-X_i{)}^2+(P_{0y}{t}_1+(P_{0y}+k\ast vec.y)t_2+P_{2y}{t}_3+P_{3y}{t}_4-Y_i{)}^2$$

对三个未知量求偏导
偏导数分别为：

• k的偏导数公式为
  $$\sum _{i=1}^{N}2\ast (P_{0x}{t}_1+(P_{0x}+k\ast vec.x)t_2+P_{2x}{t}_3+P_{3x}{t}_4-X_i)\ast t_2\ast vec.x+2\ast (P_{0y}{t}_1+(P_{0y}+k\ast vec.y)t_2+P_{2y}{t}_3+P_{3y}{t}_4-Y_i)\ast t_2\ast vec.y=0$$

• 求P2x的偏导数为0的公式为
  $$\sum _{i=1}^{N}2\ast (P_{0x}{t}_1+(P_{0x}+k\ast vec.x)t_2+P_{2x}{t}_3+P_{3x}{t}_4-X_i)\ast t_3=0$$

• 求P2y的偏导数为0的公式为
  $$\sum _{i=1}^{N}2\ast (P_{0y}{t}_1+(P_{0y}+k\ast vec.y)t_2+P_{2y}{t}_3+P_{3y}{t}_4-Y_i)\ast t_3=0$$

上述三个公式整理后得到

$$\sum _{i=1}^N((k\ast vec.x)t_2+P_{2x}{t}_3)\ast t_2\ast vec.x+((k\ast vec.y)t_2+P_{2y}{t}_3)\ast t_2\ast vec.y=\sum _{i=1}^N(X_i-P_{0x}{t}_1-P_{0x}{t}_2-P_{3x}{t}_4)\ast t_2\ast vec.x+(Y_i-P_{0y}{t}_1-P_{0y}{t}_2-P_{3y}{t}_4)\ast t_2\ast vec.y$$

$$\sum _{i=1}^N((k\ast vec.x)t_2+P_{2x}{t}_3)\ast t_3=\sum _{i=1}^N(X_i-P_{0x}{t}_1-P_{0x}{t}_2-P_{3x}{t}_4)\ast t_3$$

$$\sum _{i=1}^N((k\ast vec.y)t_2+P_{2y}{t}_3)\ast t_3=\sum _{i=1}^N(Y_i-P_{0y}{t}_1-P_{0y}{t}_2-P_{3y}{t}_4)\ast t_3$$

继续整理

$$\sum _{i=1}^{N}k\ast t_2^2(vec.x^2+vec.y^2)+P_{2x}{t}_3\ast t_2\ast vec.x+P_{2y}{t}_3\ast t_2\ast vec.y=\sum _{i=1}^N(X_i-P_{0x}{t}_1-P_{0x}{t}_2-P_{3x}{t}_4)\ast t_2\ast vec.x+(Y_i-P_{0y}{t}_1-P_{0y}{t}_2-P_{3y}{t}_4)\ast t_2\ast vec.y$$

$$\sum _{i=1}^{N}k\ast vec.x\ast t_2\ast t_3+P_{2x}{t}_3^2+P_{2y}\ast 0=\sum _{i=1}^N(X_i-P_{0x}{t}_1-P_{0x}{t}_2-P_{3x}{t}_4)\ast t_3$$

$$\sum _{i=1}^{N}k\ast vec.y\ast t_2\ast t_3+P_{2x}\ast 0+P_{2y}{t}_3^2=\sum _{i=1}^N(Y_i-P_{0y}{t}_1-P_{0y}{t}_2-P_{3y}{t}_4)\ast t_3$$

同上可以得出一个矩阵

$$\left[\begin{array}{ccc}{\sum }_{i=1}^N{t}_2^2(vec.x^2+vec.y^2)& {\sum }_{i=1}^N{t}_3\ast t_2\ast vec.x& {\sum }_{i=1}^N{t}_3\ast t_2\ast vec.y\\ & & \\ {\sum }_{i=1}^{N}vec.x\ast t_2\ast t_3& {\sum }_{i=1}^N{t}_3^2& 0\\ & & \\ {\sum }_{i=1}^{N}vec.y\ast t_2\ast t_3& 0& {\sum }_{i=1}^N{t}_3^2\end{array}\right]\left[\begin{array}{c}k\\ \\ P_{2x}\\ \\ P_{2y}\end{array}\right]=\left[\begin{array}{c}{\sum }_{i=1}^N(X_i-P_{0x}{t}_1-P_{0x}{t}_2-P_{3x}{t}_4)\ast t_2\ast vec.x+(Y_i-P_{0y}{t}_1-P_{0y}{t}_2-P_{3y}{t}_4)\ast t_2\ast vec.y\\ \\ {\sum }_{i=1}^N(X_i-P_{0x}{t}_1-P_{0x}{t}_2-P_{3x}{t}_4)\ast t_3\\ \\ {\sum }_{i=1}^N(Y_i-P_{0y}{t}_1-P_{0y}{t}_2-P_{3y}{t}_4)\ast t_3\end{array}\right]$$

同上一种拟合，利用opencv里的矩阵估计函数进行计算

两端切向都约束

两端切向都约束的情况下，就变成了两个变量需要确定，即k1和k2;
读者可以试试自行推导。

本文纯属个人见解，欢迎指正错误。

转载请标明出处！