1 问题
平面上有两个点:起始点 \(\mathbf{p}_1\)和结束点 \(\mathbf{p}_2\)。现在要在这两点之间构建一条曲线方程，方程要求如下:
① 曲线必须过 \(\mathbf{p}_1\)、 \(\mathbf{p}_2\)。
② 曲线在 \(\mathbf{p}_1\) 位置处的方向为 \(\mathbf{n}_1\)。(这里 \(\mathbf{n}_1\cdot (\mathbf{p}_1 – \mathbf{p}_2) \neq 0\))
③ 曲线在 \(\mathbf{p}_2\) 位置处的方向为 \(\mathbf{n}_2\)。(这里 \(\mathbf{n}_2\cdot (\mathbf{p}_1 – \mathbf{p}_2) \neq 0\))
④ 曲线在 \(\mathbf{p}_1\)和 \(\mathbf{p}_2\) 间尽量接近直线。
如下图，想找一条红色曲线，使得曲线过两端点，且在两端点处斜率保持一阶连续，且曲线尽量接近直线形态。
2 求解
2.1 坐标变换
这里以\(\mathbf{p}_1\)为原点，以\(\mathbf{p}_1\)指向\(\mathbf{p}_2\)方向为x正轴，建立局部右手坐标系。因此可以得到变换矩阵
(1)
$$
T=\begin{bmatrix}R&-R\mathbf{p}_1 \end{bmatrix}
$$
公式里\(R\)为\(\mathbf{p}_1\)、\(\mathbf{p}_2\)构建的SO2变换矩阵。
2.2 计算斜率
这里设\(\mathbf{n}_1\)、\(\mathbf{n}_2\)在局部坐标系下对应的斜率为\(k_1、k_2\)，则可以按如下方式获得
(2)
$$
\begin{cases}
\bar{\mathbf{n}}_1&=R\mathbf{n}_1=\begin{bmatrix}q_x &q_y\end{bmatrix}
\\
\bar{\mathbf{n}}_2&=R\mathbf{n}_2=\begin{bmatrix}p_x &p_y\end{bmatrix}
\\
\\
k_1&=\dfrac{q_y}{q_x}
\\
\\
k_2&=\dfrac{p_y}{p_x}
\end{cases}
$$
2.3 计算位置
设\(\mathbf{p}_1、\mathbf{p}_2\)在新坐标系下的表示为 \(\bar{\mathbf{p}}_1、\bar{\mathbf{p}}_2\)，则可以得到
(3)
$$
\begin{cases}
\bar{\mathbf{p}}_1&=R(\mathbf{p}_1-\mathbf{p}_1)=\begin{bmatrix}0&0\end{bmatrix}
\\
\bar{\mathbf{p}}_2&=R(\mathbf{p}_2-\mathbf{p}_1)=\begin{bmatrix}s_x&s_y\end{bmatrix}
\end{cases}
$$
设临时变量\(k_3\)
(4)
$$
k_3=\dfrac{s_y}{s_x}
$$
2.4 多项式模型
在局部坐标系下，这里假设曲线方程为如下的4次多项式
(5)
$$
y=k_1x+a_2x^2+a_3x^3+a_4x^4
$$
这里的目的就是获取曲线方程系数 \(a_2,a_3,a_4\)。
很明显，上面的曲线方程必过 \(\bar{\mathbf{p}}_1\)点，且在这点的斜率必为\(k_1\)。由于方程必须过 \(\bar{\mathbf{p}}_2\)点且在其斜率为\(k_2\)，因此将这条件代入(5)可以得到
(6)
$$
\begin{cases}
k_1s_x+a_2s_x^2+a_3s_x^3+a_4s_x^4&=s_y
\\
k_1+2a_2s_x+3a_3s_x^2+4a_4s_x^3&=k_2
\end{cases}
$$
求解公式(6)可以得到
(7)
$$
\begin{cases}
a_2&=\dfrac{1}{s_x}(a_4s_x^3+3k_3-k_2-2k_1)
\\
\\
a_3&=\dfrac{1}{s_x^2}(k_1+k_2-2k_3-2a_4s_x^3)
\end{cases}
$$
问题中的强约束①②③已经使用，而④的约束还未使用，且还有未知数\(a_4\)待求解。由于约束④要求曲线方程在区间\((0,s_x)\)内尽量接近直线。一种很自然的思想是在这个区间均匀采样 \(n-1\)个点，然后使得这些点的趋近直线，因此可构建如下最小化问题
(8)
$$
\begin{aligned}
\mathbf{min}:e&=\sum_{i=1}^{n-1}(k_1(\dfrac{i}{n}s_x)+a_2(\dfrac{i}{n}s_x)^2+a_3(\dfrac{i}{n}s_x)^3+a_4(\dfrac{i}{n}s_x)^4-\dfrac{i}{n}s_y)^2
\\
&=\sum_{i=1}^{n-1}(k_1(\dfrac{i}{n}s_x)+\dfrac{1}{s_x}(a_4s_x^3+3k_3-k_2-2k_1)(\dfrac{i}{n}s_x)^2\\&+\dfrac{1}{s_x^2}(k_1+k_2-2k_3-2a_4s_x^3)(\dfrac{i}{n}s_x)^3+a_4(\dfrac{i}{n}s_x)^4-\dfrac{i}{n}s_y)^2
\\
&=\sum_{i=1}^n(a_4s_x^4((\dfrac{i}{n})^4-2(\dfrac{i}{n})^3+(\dfrac{i}{n})^2)+s_x((k_1-k_3)(\dfrac{i}{n})+(3k_3-k_2-2k_1)(\dfrac{i}{n})^2+(k_1+k_2-2k_3)(\dfrac{i}{n})^3))^2
\end{aligned}
$$
这里设置如下中间变量
(9)
$$
\begin{aligned}
b_i&=(\dfrac{i}{n})^4-2(\dfrac{i}{n})^3+(\dfrac{i}{n})^2
\\
c_i&=(k_1-k_3)(\dfrac{i}{n})+(3k_3-k_2-2k_1)(\dfrac{i}{n})^2+(k_1+k_2-2k_3)(\dfrac{i}{n})^3)
\\
\alpha&=\sum_{i=1}^{n-1}b_i^2=\dfrac{n^8+20n^2-21}{630n^7}
\\
\beta&=\sum_{i=1}^{n-1}b_ic_i=(k_1-k_3)\dfrac{n^4-1}{60n^3}+(3k_3-k_2-2k_1)\dfrac{2n^6-7n^2+5}{210n^5}+(k_1+k_2-2k_3)\dfrac{n^6-7n^2+6}{168n^5}
\end{aligned}
$$
将(9)代入(8)有
(10)
$$
\begin{aligned}
\mathbf{min}:e&=\sum_{i=1}^{n-1}(a_4b_is_x^4+c_is_x)^2
\\
&=a_4s_x^8\sum_{i=1}^{n-1}b_i^2+2a_4s_x^5\sum_{i=1}^{n-1}b_ic_i+s_x^2\sum_{i=1}^{n-1}c_i^2
\\
&=\alpha a_4^2s_x^8+2\beta a_4s_x^5+s_x^2\sum_{i=1}^{n-1}c_i^2
\\
&=\alpha s_x^8(a_4+\dfrac{\beta}{\alpha s_x^3})^2+s_x^2\sum_{i=1}^{n-1}c_i^2-\dfrac{\beta^2s_x^2}{\alpha}
\\
&\Downarrow
\\
a_4&=-\dfrac{\beta}{\alpha s_x^3}
\end{aligned}
$$
特别地，当 \(n\to \infty\)有
(11)
$$
\lim_{n\to \infty}a_4=\dfrac{9}{4s_x^3}(k_2-k_1)
$$
将(11)代入(7)可以得到
(12)
$$
\begin{cases}
a_2&=\dfrac{1}{4s_x}(12k_3+5k_2-17k_1)
\\
\\
a_3&=\dfrac{1}{2s_x^2}(11k_1-7k_2-4k_3)
\\
\\
a_4&=\dfrac{9}{4s_x^3}(k_2-k_1)
\end{cases}
$$
最终，通过(12)计算得到公式(5)曲线方程系数。
3 变种
上面(8)构建的目标是点的位置接近直线，这里考虑曲线中间线的斜率与直线斜率一致，因此可以构造如下最小二乘目标函数
(13)$$
\mathbf{min}:e=\sum_{i=1}^{n-1}(k_1+2a_2(\dfrac{i}{n}s_x)+3a_3(\dfrac{i}{n}s_x)^2+4a_4(\dfrac{i}{n}s_x)^3-k_3)^2
$$
按上面类似处理方式，可以最终得到如下结果
(14)
$$
\begin{cases}
a_2&=\dfrac{3}{4s_x}(4k_3+k_2-5k_1)
\\
\\
a_3&=\dfrac{1}{2s_x^2}(9k_1-5k_2-4k_3)
\\
\\
a_4&=\dfrac{7}{4s_x^3}(k_2-k_1)
\end{cases}
$$
4 后记
使用中，需要注意，模型是在局部坐标系下建立的方程。如果想整体位置更接近直线，可考虑使用(12)的系数，如果想整体斜率更接近直线，可考虑(14)系数。如果想计算一些信息，比如位置、方向，可以先通过曲线方程，获取局部坐标系下的表示方式，然后再回转回原坐标系。对于一些不变的信息，比如长度、曲率等，则可以直接使用局部坐标系下的结果。
还有，注意本函数与样条函数的区别。本函数除了要求要过对应点以及斜率连续外，还要求曲线内部尽量直，也就说在满足约束外，让创建的曲线曲率尽量小。这在一些特殊场景里有所应用。
目前 MathSword 已经将本功能集成到函数DataSmoothLineIn2Point。