根据三点估算导数/Bessel Tangent方法

在曲线拟合问题中，通常需要根据已知曲线上的离散点，估算出曲线在端点处的导数（严格来说是导矢），常用的一种导数估计方法称为“Bessel Tangent”方法。

已知点\(P_{i-1},P_i,P_{i+1}\)，求一种导数（矢）估算算法，计算\(P_{i-1}^\prime,P_i^\prime,P_{i+1}^\prime\)。Bessel tangent 算法的主要思路是，插值一条通过三点的二次曲线，然后通过计算曲线上三点的导数的方式给出三点的估计导数。

令 \(P_{i-1},P_{i},P_{i+1}\) 的参数分别为 \(t_{i-1},t_{i},t_{i+1}\)

插值的曲线的表达式如下：

$$P(t)=\frac{(t-t_i)(t-t_{i+1})}{(t_{i-1}-t_i)(t_{i-1}-t_{i+1})}(P_{i-1}-P_{i})+\frac{(t-t_{i-1})(t-t_{i})}{(t_{i+1}-t_{i-1})(t_{i+1}-t_{i})}(P_{i+1}-P_{i})+P_{i}\tag1$$

1) \(P_i\)导数

$$P(t_i)^\prime=\frac{t_i-t_{i+1}}{(t_{i-1}-t_i)(t_{i-1}-t_{i+1})}(P_{i-1}-P_{i})+\frac{t_i-t_{i-1}}{(t_{i+1}-t_{i-1})(t_{i+1}-t_{i})}(P_{i+1}-P_{i})$$

令 \(\triangle t_i = t_{i+1}-t_i;\triangle P_i=\frac{P_{i+1}-P_i}{\triangle t_i}\)

$$\begin{align}
P(t_i)^\prime&=\frac{t_i-t_{i+1}}{t_{i-1}-t_{i+1}}\triangle P_{i-1}+\frac{t_i-t_{i-1}}{t_{i+1}-t_{i-1}}\triangle P_{i}\\
&= \frac{t_i-t_{i+1}}{t_{i-1}-t_{i+1}}\triangle P_{i-1}+\frac{t_i-t_{i-1}}{t_{i+1}-t_{i-1}}\triangle P_{i} \\
&= \frac{\triangle t_{i}}{t_{i+1}-t_{i}+t_{i}-t_{i-1}}\triangle P_{i-1} + \frac{\triangle t_{i-1}}{t_{i+1}-t_{i}+t_{i}-t_{i-1}}\triangle P_{i}\\
&=\frac{\triangle t_{i}}{\triangle t_{i}+\triangle t_{i-1}}\triangle P_{i-1}+\frac{\triangle t_{i-1}}{\triangle t_{i}+\triangle t_{i-1}}\triangle P_{i} \tag2
\end{align}$$即\(P_i^\prime\)是\(\triangle P_{i-1}\)，\(\triangle P_{i}\)的加权平均。

2)\(P_{i-1},P_{i+1}\)导数

根据公式(1)

$$\begin{align}
P(t_{i-1})^\prime&=\frac{2t_{i-1}-2t_{i+1}+t_{i+1}-t_{i}}{t_{i-1}-t_{i+1}}\triangle P_{i-1}+\frac{t_{i-1}-t_{i}}{t_{i+1}-t_{i-1}}\triangle P_{i}\\
&=2\triangle P_{i-1}-(\frac{\triangle t_{i}}{\triangle t_{i}+\triangle t_{i-1}}\triangle P_{i-1}+\frac{\triangle t_{i-1}}{\triangle t_{i}+\triangle t_{i-1}}\triangle P_{i})\\
&=2\triangle P_{i-1}-P(t_i)^\prime \tag3
\end{align}$$

$$\begin{align}
P(t_{i+1})^\prime&=-(\frac{\triangle t_{i}}{\triangle t_{i}+\triangle t_{i-1}}\triangle P_{i-1}+\frac{\triangle t_{i-1}}{\triangle t_{i}+\triangle t_{i-1}}\triangle P_{i} )+2\triangle P_i \\
&=2\triangle P_i-P(t_{i})^\prime \tag4
\end{align}$$

参数化的选择

对于三点的参数可以按照累计弦长化确定参数$$ \left\{
\begin{aligned}
t_{i-1} & = 0 \\
t_{i} & = t_0+|P_{i-1}P_i| \\
t_{i+1} & = t_1+|P_iP_{i+1}|
\end{aligned}
\right.
$$

我对算法的C++实现是：

 /*!
  *\brief Bessel Tangent 方法估算导数
 *\ param const Point & p0,p1,p2 待求点
 *\ param Vec3d & p0deriv,p1deriv,p2deriv 三点上导数的估计值
 *\ Returns:   
 */
void BesselTanget(const Point& p0,const Point& p1,const Point& p2, 
				  Vec3d& p0deriv,Vec3d& p1deriv,Vec3d& p2deriv)
{
	double delta_t1 =PointDistance(p2,p1);
	double delta_t0 = PointDistance(p1,p0);
	Vec3d delta_p0 = (p1-p0)/delta_t0;
	Vec3d delta_p1 = (p2-p1)/delta_t1;
	double sum = delta_t0+delta_t1;
	p1deriv = delta_t1/sum * delta_p0 + delta_t0/sum * delta_p1;
	p0deriv = 2*delta_p0 - p1deriv;
	p2deriv = 2*delta_p1 - p1deriv;
}

参考资料：

The Essentials of CAGD : Chapter 11 Working with B Spline Curves 

于谦. 基于切向约束的二次B样条插值曲线的研究[D]. 山东大学, 2009.