1 问题
平面上有\(m+q\)个点(\(q>0\)),其分别为\(\mathbf{p}_1、\mathbf{p}_2、…、\mathbf{p}_m、\mathbf{s}_1、\mathbf{s}_2、…、\mathbf{s}_q\)。现在要构建一条满足如下条件的曲线方程:
①曲线在\(\mathbf{s}_1\)位置处的方向为\(\mathbf{k}\)。
②曲线必须通过 \(\mathbf{s}_1、\mathbf{s}_2、…、\mathbf{s}_q\)。
③ 曲线距离 \(\mathbf{p}_1、\mathbf{p}_2、…、\mathbf{p}_m\) 的距离尽量小。
2 求解
2.1 坐标变换
这里以 \(\mathbf{s}_1\) 为原点，然后在 \(m+q\) 个点里优先选择距离 \(\mathbf{s}_1\) 较远的一个点 \(\mathbf{p}\) 且 \(\mathbf{k}\cdot(\mathbf{p}-\mathbf{s}_1)\neq 0\),然后以\(\mathbf{s}_1\)指向\(\mathbf{p}\)的方向为\(x\)轴正向建立局部坐标系，因此可以得到变换矩阵
(1)
$$
T=\begin{bmatrix}R&-R\mathbf{s}_1\end{bmatrix}
$$
公式里\(R\)为\(\mathbf{s}_1、\mathbf{p}\)构建的SO2变换矩阵。
设\(\mathbf{p}_1、\mathbf{p}_2、…、\mathbf{p}_m\)在新坐标系下的表示为\((x_1,y_1)、(x_2,y_2)、…、(x_m,y_m)\),\(\mathbf{s}_1、\mathbf{s}_2、…、\mathbf{s}_q\)在新坐标系下的表示为\((0,0)、(\bar{x}_2,\bar{y}_2)、…、(\bar{x}_q,\bar{y}_q)\)。
设\(\mathbf{k}\)对应斜率为\(k\)，则
(2)
$$
\begin{cases}
\begin{bmatrix}
s_x
\\
s_y
\end{bmatrix}&=R(\mathbf{k}-\mathbf{s}_1)
\\
\\
k&=\dfrac{s_y}{s_x}
\end{cases}
$$
2.2 多项式模型
这里假设曲线方程为 \(n\) 次多项式（这里\(n>q\)）。即在局部坐标系下，曲线方程为
(3)
$$
y=kx+a_2x^2+a_3x^3+…+a_nx^n
$$
很明显，(3)严格满足第①点，并且 \(\mathbf{s}_1\)也在曲线上，根据第②点其它点在曲线上的条件可以得到
(4)
$$
\begin{cases}
\bar{y}_2&=k\bar{x}_2+a_2\bar{x}_2^2+a_3\bar{x}_2^3+…+a_n\bar{x}_2^n
\\
\bar{y}_3&=k\bar{x}_3+a_2\bar{x}_3^2+a_3\bar{x}_3^3+…+a_n\bar{x}_3^n
\\
…
\\
\bar{y}_q&=k\bar{x}_q+a_2\bar{x}_q^2+a_3\bar{x}_q^3+…+a_n\bar{x}_q^n
\end{cases}
$$
设置如下中间变量
(5)
$$
\begin{aligned}
\mathbf{a}_1&=\begin{bmatrix}
a_2
\\
a_3
\\
…
\\
a_{n-q}
\end{bmatrix},
\mathbf{a}_2=\begin{bmatrix}
a_{n-q+1}
\\
a_{n-q+2}
\\
…
\\
a_n
\end{bmatrix},
\mathbf{b}=\begin{bmatrix}
\bar{y}_2-k\bar{x}_2
\\
\bar{y}_3-k\bar{x}_3
\\
…
\\
\bar{y}_q-k\bar{x}_q
\end{bmatrix}
\\
\\
B_1&=\begin{bmatrix}
\bar{x}_2^2&\bar{x}_2^3&…&\bar{x}_2^{n-q}
\\
\bar{x}_3^2&\bar{x}_3^3&…&\bar{x}_3^{n-q}
\\
…
\\
\bar{x}_q^2&\bar{x}_q^3&…&\bar{x}_q^{n-q}
\end{bmatrix}
\\
\\
B_2&=\begin{bmatrix}
\bar{x}_2^{n-q+1}&\bar{x}_2^{n-q+2}&…&\bar{x}_2^n
\\
\bar{x}_3^{n-q+1}&\bar{x}_3^{n-q+2}&…&\bar{x}_3^n
\\
…
\\
\bar{x}_q^{n-q+1}&\bar{x}_q^{n-q+2}&…&\bar{x}_q^n
\end{bmatrix}
\\
\\
\mathbf{g}&=B_2^{-1}\mathbf{b}
\\
G&=B_2^{-1}B_1
\end{aligned}
$$
将(5)代入(4)可以得到
(6)
$$
\mathbf{a}_2=\mathbf{g}-G\mathbf{a}_1
$$
公式里的 \(\mathbf{g}、G\) 为关系矩阵，如果对求解效率有要求，实际处理时可以跳过矩阵求逆而直接按线性拓扑关系计算此矩阵。
接下来，为了满足第③个条件，一种很自然的思想是，构建一个最小化目标函数，使得所有点距离曲线距离的平方和最小。对于距离，比较直观的是欧式距离，由于计算点到多项式曲线的欧式距离算法复杂，因此这不是构建目标的首选方式。这里退而求其次，构建每个点在\(y\)方向距离最小(注意，如果所有点均在曲线上，这种方式构建的最小化目标函数和欧式距离构建的最小化目标函数是等价的)，也就是构建如下的最小化目标函数:
(7)
$$
\mathbf{min}:e=\sum_{i=1}^m(kx_i+a_2x_i^2+a_3x_i^3+…+a_nx_i^n-y_i)^2
$$
在导出公式前，先设置如下中间变量
(8)
$$
\begin{aligned}
\mathbf{c}&=\begin{bmatrix}
y_1-kx_1
\\
y_2-kx_2
\\
…
\\
y_m-kx_m
\end{bmatrix}
\\
\\
A_1&=\begin{bmatrix}
x_1^2&x_1^3&…&x_1^{n-q}
\\
x_2^2&x_2^3&…&x_2^{n-q}
\\
…
\\
x_m^2&x_m^3&…&x_m^{n-q}
\end{bmatrix}
\\
\\
A_2&=\begin{bmatrix}
x_1^{n-q+1}&x_1^{n-q+2}&…&x_1^n
\\
x_2^{n-q+1}&x_2^{n-q+2}&…&x_2^n
\\
…
\\
x_m^{n-q+1}&x_m^{n-q+2}&…&x_m^n
\end{bmatrix}
\\
\\
D&=A_1-A_2G
\\
\mathbf{d}&=\mathbf{c}-A_2\mathbf{g}
\end{aligned}
$$
对于公式(7)的最小化目标函数，等效于求解如下线性方程组的最小二乘解
(9)
$$
A_1\mathbf{a}_1+A_2\mathbf{a}_2=\mathbf{c}
$$
将(6)代入(9)可以得到
(10)
$$
D\mathbf{a}_1=\mathbf{d}
$$
求解(10)的最小二乘解，即是(7)的最佳解。因此通过(5)、(6)、(8)、(10)可以求得曲线方程所有系数。
3 带权求解
对于1中的第③个条件，假如 \(\mathbf{p}_1、\mathbf{p}_2、…、\mathbf{p}_m\) 影响到距离目标函数的权重分别为 \(w_1、w_2、…、w_m\)，则很容易得到，最小化目标函数的解将变为求解如下线性方程的最小二乘解
(11)
$$
WD\mathbf{a}_1=W\mathbf{d}
$$
公式中的\(W\)矩阵如下
(12)
$$
W=\begin{bmatrix}
w_1&0&…&0
\\
0&w_2&…&0
\\
…
\\
0&0&…&w_m
\end{bmatrix}
$$
因此通过(5)、(6)、(8)、(11)可以求得曲线方程所有系数。
4 后记
使用中，需要注意，模型是在局部坐标系下建立的方程。如果想计算一些信息，比如位置、方向，可以先通过曲线方程，获取局部坐标系下的表示方式，然后再回转回原坐标系。对于一些不变的信息，比如长度、曲率等，则可以直接使用局部坐标系下的结果。
还有，注意本函数与样条函数的区别。