在曲线拟合问题中，通常需要根据已知曲线上的离散点，估算出曲线在端点处的导数（严格来说是导矢），常用的一种导数估计方法称为“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.