曲率、曲率（对弧长）的导数以及曲率导数（对弧长）的导数的计算

笛卡尔坐标系与Frenet坐标系互转，可能需要曲率的导数信息。此出给出推导过程与计算式，方便以后写代码时查阅。

1. 曲线的表示形式

二维平面上的曲线有两种参数化形式，如下所示：

参数方程1
{xt=x(t)yt=y(t)\left\{\begin{matrix} x_t=x(t) \\ y_t=y(t) \end{matrix} \right. {xt=x(t)yt=y(t)
参数方程2
{xt=xtyt=y(xt)\left\{\begin{matrix} x_t=x_t \\ y_t=y(x_t) \end{matrix} \right. {xt=xtyt=y(xt)

以上两种参数方程都可以唯一确定一条二维平面内的曲线。因此，下文计算的曲率、曲率的导数以及曲率导数的导数的公式都有两种等价的形式。

2. 曲率计算公式及推导

先给出熟悉的曲率计算公式：
k=x′y′′−y′x′′(x′2+y′2)32（对应参数方程1）(1)k=\frac{x'y''-y'x''}{(x'^{2}+y'^{2})^{\frac{3}{2}}} （对应参数方程1）\tag{1}k=(x′2+y′2)23x′y′′−y′x′′（对应参数方程1）(1)以及：
k=y′′(1+y′2)32（对应参数方程2）(2)k=\frac{y''}{(1+y'^2)^{\frac{3}{2}}}（对应参数方程2）\tag{2}k=(1+y′2)23y′′（对应参数方程2）(2)

2.1 参数方程1曲率公式推导

假定点(xt,yt)(x_t, y_t)(xt,yt)处的切角为α\alphaα，则此点处曲线的斜率为tan⁡(α)\tan(\alpha)tan(α)。xt,ytx_t,y_txt,yt的变量都为ttt，假定ttt有一个小的增量Δt\Delta_tΔt，则，xtx_txt与yty_tyt相应的都有一个小的增量Δx\Delta_xΔx与Δy\Delta_yΔy。如下图所示，当Δt\Delta_tΔt很小时，ΔyΔx≈tan(α)\frac{\Delta_y}{\Delta_x}\approx tan(\alpha)ΔxΔy≈tan(α)。

当Δt→0\Delta_t\rightarrow 0Δt→0时，Δx\Delta_xΔx与Δy\Delta_yΔy为xtx_txt与yty_tyt在ttt处的微分，一般分别表示为dxdxdx与dydydy。此时有(其中x′x'x′与y′y'y′表示函数x(t),y(t)x(t),y(t)x(t),y(t)对ttt的导数):dydx=tan⁡(α)⇒dydtdxdt=tan⁡(α)⇒y′x′=tan⁡(α)\frac{dy}{dx}=\tan(\alpha)\Rightarrow \frac{\frac{dy}{dt}}{\frac{dx}{dt}}=\tan(\alpha)\Rightarrow \frac{y'}{x'}=\tan(\alpha)dxdy=tan(α)⇒dtdxdtdy=tan(α)⇒x′y′=tan(α)得：

y′x′=tan⁡(α)\frac{y'}{x'}=\tan(\alpha)x′y′=tan(α)

对上式两边分别求导得：
y′′x′−x′′y′x′2dt=(1+tan⁡2(α))dα=(x′2+y′2x′2)dα\frac{y''x'-x''y'}{x'^{2}}dt=(1+\tan^2(\alpha))d\alpha=(\frac{x'^2+y'^2}{x'^2})d\alphax′2y′′x′−x′′y′dt=(1+tan2(α))dα=(x′2x′2+y′2)dα
化简得：
y′′x′−x′′y′x′2+y′2dt=dα\frac{y''x'-x''y'}{x'^2+y'^2}dt=d\alphax′2+y′2y′′x′−x′′y′dt=dα

又ds=x′2+y′2dtds=\sqrt{x'^2+y'^2}dtds=x′2+y′2dt,代入上式得到式（1）所示的曲率计算公式（曲率为角度对弧度的导数）：

k=dαds=y′′x′−x′′y′(x′2+y′2)32k=\frac{d\alpha}{ds}=\frac{y''x'-x''y'}{(x'^2+y'^2)^{\frac{3}{2}}}k=dsdα=(x′2+y′2)23y′′x′−x′′y′

2.2 参数方程2曲率公式推导

此时曲线的自变量为xxx，曲线在点(xt,yt)(x_t, y_t)(xt,yt)处的导数为y′=f′(x)=tan⁡(α)y'=f'(x)=\tan(\alpha)y′=f′(x)=tan(α)得：
y′=tan⁡(α)y'=\tan(\alpha)y′=tan(α)

对上式两边分别求导得：
y′′dx=(1+tan⁡2(α))dα=(1+y′2)dαy''dx=(1+\tan^{2}(\alpha))d\alpha=(1+y'^2)d\alphay′′dx=(1+tan2(α))dα=(1+y′2)dα

化简得：
y′′1+y′2dx=dα\frac{y''}{1+y'^2}dx=d\alpha1+y′2y′′dx=dα

又ds=1+y′2dxds=\sqrt{1+y'2}dxds=1+y′2dx，代入上式得到式(2)所示的曲率计算公式：

k=y′′(1+y′2)32k=\frac{y''}{(1+y'^2)^{\frac{3}{2}}}k=(1+y′2)23y′′

2.3 小结

两种参数方程得到的曲率公式推导过程相似，最终公式形式也差不多。在表示曲线时，不同情况下用到的参数化方程不一样。为了简便，可以统一两种参数方程，令x(t)=tx(t)=tx(t)=t时，参数方程1就变成了参数方程2。此时，x′=1,x′′=0x'=1,x''=0x′=1,x′′=0，代入式(1)就得到式(2)。下文中，只求针对参数方程1的曲率导数k′k'k′以及曲率导数的导数k′′k''k′′。

3. 曲率的导数（或称为变化率）公式及推导

对曲率公式k=y′′x′−x′′y′(x′2+y′2)32k=\frac{y''x'-x''y'}{(x'^2+y'^2)^{\frac{3}{2}}}k=(x′2+y′2)23y′′x′−x′′y′两边分别求导：

k′=dkds=(y′′x′−x′′y′)′(x′2+y′2)32−((x′2+y′2)32)′(y′′x′−x′′y′)(x′2+y′2)3dtdsk'=\frac{dk}{ds}=\frac{\frac{(y''x'-x''y')'(x'^2+y'^2)^{\frac{3}{2}}-((x'^2+y'^2)^{\frac{3}{2}})'(y''x'-x''y')}{(x'^2+y'^2)^3}dt}{ds}k′=dsdk=ds(x′2+y′2)3(y′′x′−x′′y′)′(x′2+y′2)23−((x′2+y′2)23)′(y′′x′−x′′y′)dt

将ds=x′2+y′2dtds=\sqrt{x'^2+y'^2}dtds=x′2+y′2dt代入上式得：
k′=(y′′′x′+y′′x′′−x′′′y′−x′′y′′)(x′2+y′2)32−32(x′2+y′2)12(2x′x′′+2y′y′′)(y′′x′−x′′y′)(x′2+y′2)3(dtds)=(x′2+y′2)12[(y′′′x′−x′′′y′)(x′2+y′2)−3(x′x′′+y′y′′)(y′′x′−x′′y′)](x′2+y′2)31(x′2+y′2)12=(y′′′x′+y′′x′′−x′′′y′−x′′y′′)(x′2+y′2)32−32(x′2+y′2)12(2x′x′′+2y′y′′)(y′′x′−x′′y′)(x′2+y′2)3(dtds)=(y′′′x′−x′′′y′)(x′2+y′2)−3(x′x′′+y′y′′)(y′′x′−x′′y′)(x′2+y′2)3k'=\frac{(y'''x'+y''x''-x'''y'-x''y'')(x'^2+y'^2)^{\frac{3}{2}}-\frac{3}{2}(x'^2+y'^2)^{\frac{1}{2}}(2x'x''+2y'y'')(y''x'-x''y')}{(x'^2+y'^2)^3}(\frac{dt}{ds})\\=\frac{(x'^2+y'^2)^{\frac{1}{2}}[(y'''x'-x'''y')(x'^2+y'^2)-3(x'x''+y'y'')(y''x'-x''y')]}{(x'^2+y'^2)^3}\frac{1}{(x'^2+y'^2)^{\frac{1}{2}}}\\=\frac{(y'''x'+y''x''-x'''y'-x''y'')(x'^2+y'^2)^{\frac{3}{2}}-\frac{3}{2}(x'^2+y'^2)^{\frac{1}{2}}(2x'x''+2y'y'')(y''x'-x''y')}{(x'^2+y'^2)^3}(\frac{dt}{ds})\\=\frac{(y'''x'-x'''y')(x'^2+y'^2)-3(x'x''+y'y'')(y''x'-x''y')}{(x'^2+y'^2)^3}k′=(x′2+y′2)3(y′′′x′+y′′x′′−x′′′y′−x′′y′′)(x′2+y′2)23−23(x′2+y′2)21(2x′x′′+2y′y′′)(y′′x′−x′′y′)(dsdt)=(x′2+y′2)3(x′2+y′2)21[(y′′′x′−x′′′y′)(x′2+y′2)−3(x′x′′+y′y′′)(y′′x′−x′′y′)](x′2+y′2)211=(x′2+y′2)3(y′′′x′+y′′x′′−x′′′y′−x′′y′′)(x′2+y′2)23−23(x′2+y′2)21(2x′x′′+2y′y′′)(y′′x′−x′′y′)(dsdt)=(x′2+y′2)3(y′′′x′−x′′′y′)(x′2+y′2)−3(x′x′′+y′y′′)(y′′x′−x′′y′)

最后提取k′k'k′的计算公式如下：

k′=(y′′′x′−x′′′y′)(x′2+y′2)−3(x′x′′+y′y′′)(y′′x′−x′′y′)(x′2+y′2)3（对应参数方程1）(3)k'=\frac{(y'''x'-x'''y')(x'^2+y'^2)-3(x'x''+y'y'')(y''x'-x''y')}{(x'^2+y'^2)^3} （对应参数方程1） \tag{3}k′=(x′2+y′2)3(y′′′x′−x′′′y′)(x′2+y′2)−3(x′x′′+y′y′′)(y′′x′−x′′y′)（对应参数方程1）(3)

令x′=1,x′′=0,x′′′=0x'=1, x''=0, x'''=0x′=1,x′′=0,x′′′=0代入式(3)可得参数方程2的k′k'k′计算公式如下：

k′=y′′′+y′′′y′2−3y′y′′2(1+y′2)3(对应参数方程2)(4)k'=\frac{y'''+y'''y'^2-3y'y''^2}{(1+y'^2)^3} (对应参数方程2)\tag{4}k′=(1+y′2)3y′′′+y′′′y′2−3y′y′′2(对应参数方程2)(4)

4. 曲率导数的导数公式及推导

k′′=dk′ds=d(y′′′x′−x′′′y′)(x′2+y′2)−3(x′x′′+y′y′′)(y′′x′−x′′y′)(x′2+y′2)3dsk''=\frac{dk'}{ds}=\frac{d\frac{(y'''x'-x'''y')(x'^2+y'^2)-3(x'x''+y'y'')(y''x'-x''y')}{(x'^2+y'^2)^3}}{ds}k′′=dsdk′=dsd(x′2+y′2)3(y′′′x′−x′′′y′)(x′2+y′2)−3(x′x′′+y′y′′)(y′′x′−x′′y′)

算不下去了。。。。。

不过计算方式和曲率的导数一样。

5. C++与Python函数实现

计算曲率和曲率导数的c++函数如下：

// kappa = (dx * d2y - dy * d2x) / [(dx * dx + dy * dy)^(3/2)]
double CurveMath::ComputeCurvature(const double dx, const double d2x,const double dy, const double d2y) {const double a = dx * d2y - dy * d2x;auto norm_square = dx * dx + dy * dy;auto norm = std::sqrt(norm_square);const double b = norm * norm_square;return a / b;
}double CurveMath::ComputeCurvatureDerivative(const double dx, const double d2x,const double d3x, const double dy,const double d2y,const double d3y) {const double a = dx * d2y - dy * d2x;const double b = dx * d3y - dy * d3x;const double c = dx * d2x + dy * d2y;const double d = dx * dx + dy * dy;return (b * d - 3.0 * a * c) / (d * d * d);
}

计算曲率和曲率导数的python函数如下：

def ComputeCurvature(dx, ddx, dy, ddy):a = dx*ddy - dy*ddxnorm_square = dx*dx+dy*dynorm = sqrt(norm_square)b = norm*norm_squarereturn a/bdef ComputeCurvatureDerivative(dx, ddx, dddx, dy, ddy, dddy):a = dx*ddy-dy*ddxb = dx*dddy-dy*dddxc = dx*ddx+dy*ddyd = dx*dx+dy*dyreturn (b*d-3.0*a*c)/(d*d*d)def ComputeCurvature(dy, ddy):a = ddy norm_square = 1+dy*dynorm = sqrt(norm_square)b = norm*norm_squarereturn a/bdef ComputeCurvatureDerivative(dy, ddy, dddy):a = ddyb = dddyc = dy*ddyd = 1.0+dy*dyreturn (b*d-3.0*a*c)/(d*d*d)

END
by windSeS
2020.9.6