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

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

1. 曲线的表示形式

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

• 参数方程1
  $$\{\begin{array}{c}{x}_t=x(t)\\ y_t=y(t)\end{array}$$

• 参数方程2
  $$\{\begin{array}{c}{x}_t=x_t\\ y_t=y(x_t)\end{array}$$

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

2. 曲率计算公式及推导

先给出熟悉的曲率计算公式：
$$k=\frac{x^{\prime }{y}^{\prime \prime }-y^{\prime }{x}^{\prime \prime }}{(x^{\prime 2}+y^{\prime 2}{)}^{\frac{3}{2}}}（对应参数方程1）\text{\hspace{1em}(1)}$$以及：
$$k=\frac{y^{\prime \prime }}{(1+y^{\prime 2}{)}^{\frac{3}{2}}}（对应参数方程2）\text{\hspace{1em}(2)}$$

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

假定点$(x_t,y_t)$处的切角为$\alpha$，则此点处曲线的斜率为$\tan (\alpha )$。$x_t,y_t$的变量都为$t$，假定$t$有一个小的增量${\Delta }_t$，则，$x_t$与$y_t$相应的都有一个小的增量${\Delta }_x$与${\Delta }_y$。如下图所示，当${\Delta }_t$很小时，$\frac{{\Delta }_y}{{\Delta }_x}\approx tan(\alpha )$。

当${\Delta }_t\to 0$时，${\Delta }_x$与${\Delta }_y$为$x_t$与$y_t$在$t$处的微分，一般分别表示为$dx$与$dy$。此时有(其中$x^{\prime }$与$y^{\prime }$表示函数$x(t),y(t)$对$t$的导数):$\frac{dy}{dx}=\tan (\alpha )\Rightarrow \frac{\frac{dy}{dt}}{\frac{dx}{dt}}=\tan (\alpha )\Rightarrow \frac{y^{\prime }}{x^{\prime }}=\tan (\alpha )$得：

$$\frac{y^{\prime }}{x^{\prime }}=\tan (\alpha )$$

对上式两边分别求导得：
$$\frac{y^{\prime \prime }{x}^{\prime }-x^{\prime \prime }{y}^{\prime }}{x^{\prime 2}}dt=(1+{\tan }^2(\alpha ))d\alpha =(\frac{x^{\prime 2}+y^{\prime 2}}{x^{\prime 2}})d\alpha$$
化简得：
$$\frac{y^{\prime \prime }{x}^{\prime }-x^{\prime \prime }{y}^{\prime }}{x^{\prime 2}+y^{\prime 2}}dt=d\alpha$$

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

$$k=\frac{d\alpha }{ds}=\frac{y^{\prime \prime }{x}^{\prime }-x^{\prime \prime }{y}^{\prime }}{(x^{\prime 2}+y^{\prime 2}{)}^{\frac{3}{2}}}$$

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

此时曲线的自变量为$x$，曲线在点$(x_t,y_t)$处的导数为$y^{\prime }=f^{\prime }(x)=\tan (\alpha )$得：
$$y^{\prime }=\tan (\alpha )$$

对上式两边分别求导得：
$$y^{\prime \prime }dx=(1+{\tan }^2(\alpha ))d\alpha =(1+y^{\prime 2})d\alpha$$

化简得：
$$\frac{y^{\prime \prime }}{1+y^{\prime 2}}dx=d\alpha$$

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

$$k=\frac{y^{\prime \prime }}{(1+y^{\prime 2}{)}^{\frac{3}{2}}}$$

2.3 小结

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

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

对曲率公式$k=\frac{y^{\prime \prime }{x}^{\prime }-x^{\prime \prime }{y}^{\prime }}{(x^{\prime 2}+y^{\prime 2}{)}^{\frac{3}{2}}}$两边分别求导：

$$k^{\prime }=\frac{dk}{ds}=\frac{\frac{(y^{\prime \prime }{x}^{\prime }-x^{\prime \prime }{y}^{\prime }{)}^{\prime }(x^{\prime 2}+y^{\prime 2}{)}^{\frac{3}{2}}-((x^{\prime 2}+y^{\prime 2}{)}^{\frac{3}{2}}{)}^{\prime }(y^{\prime \prime }{x}^{\prime }-x^{\prime \prime }{y}^{\prime })}{(x^{\prime 2}+y^{\prime 2}{)}^3}dt}{ds}$$

将$ds=\sqrt{x^{\prime 2}+y^{\prime 2}}dt$代入上式得：
$$\begin{gathered} k^{\prime }=\frac{(y^{\prime \prime \prime }{x}^{\prime }+y^{\prime \prime }{x}^{\prime \prime }-x^{\prime \prime \prime }{y}^{\prime }-x^{\prime \prime }{y}^{\prime \prime })(x^{\prime 2}+y^{\prime 2}{)}^{\frac{3}{2}}-\frac{3}{2}(x^{\prime 2}+y^{\prime 2}{)}^{\frac{1}{2}}(2x^{\prime }{x}^{\prime \prime }+2y^{\prime }{y}^{\prime \prime })(y^{\prime \prime }{x}^{\prime }-x^{\prime \prime }{y}^{\prime })}{(x^{\prime 2}+y^{\prime 2}{)}^3}(\frac{dt}{ds}) \\ =\frac{(x^{\prime 2}+y^{\prime 2}{)}^{\frac{1}{2}}[(y^{\prime \prime \prime }{x}^{\prime }-x^{\prime \prime \prime }{y}^{\prime })(x^{\prime 2}+y^{\prime 2})-3(x^{\prime }{x}^{\prime \prime }+y^{\prime }{y}^{\prime \prime })(y^{\prime \prime }{x}^{\prime }-x^{\prime \prime }{y}^{\prime })]}{(x^{\prime 2}+y^{\prime 2}{)}^3}\frac{1}{(x^{\prime 2}+y^{\prime 2}{)}^{\frac{1}{2}}} \\ =\frac{(y^{\prime \prime \prime }{x}^{\prime }+y^{\prime \prime }{x}^{\prime \prime }-x^{\prime \prime \prime }{y}^{\prime }-x^{\prime \prime }{y}^{\prime \prime })(x^{\prime 2}+y^{\prime 2}{)}^{\frac{3}{2}}-\frac{3}{2}(x^{\prime 2}+y^{\prime 2}{)}^{\frac{1}{2}}(2x^{\prime }{x}^{\prime \prime }+2y^{\prime }{y}^{\prime \prime })(y^{\prime \prime }{x}^{\prime }-x^{\prime \prime }{y}^{\prime })}{(x^{\prime 2}+y^{\prime 2}{)}^3}(\frac{dt}{ds}) \\ =\frac{(y^{\prime \prime \prime }{x}^{\prime }-x^{\prime \prime \prime }{y}^{\prime })(x^{\prime 2}+y^{\prime 2})-3(x^{\prime }{x}^{\prime \prime }+y^{\prime }{y}^{\prime \prime })(y^{\prime \prime }{x}^{\prime }-x^{\prime \prime }{y}^{\prime })}{(x^{\prime 2}+y^{\prime 2}{)}^3} \end{gathered}$$

最后提取$k^{\prime }$的计算公式如下：

$$k^{\prime }=\frac{(y^{\prime \prime \prime }{x}^{\prime }-x^{\prime \prime \prime }{y}^{\prime })(x^{\prime 2}+y^{\prime 2})-3(x^{\prime }{x}^{\prime \prime }+y^{\prime }{y}^{\prime \prime })(y^{\prime \prime }{x}^{\prime }-x^{\prime \prime }{y}^{\prime })}{(x^{\prime 2}+y^{\prime 2}{)}^3}（对应参数方程1）\text{\hspace{1em}(3)}$$

令$x^{\prime }=1,x^{\prime \prime }=0,x^{\prime \prime \prime }=0$代入式(3)可得参数方程2的$k^{\prime }$计算公式如下：

$$k^{\prime }=\frac{y^{\prime \prime \prime }+y^{\prime \prime \prime }{y}^{\prime 2}-3y^{\prime }{y}^{\prime \prime 2}}{(1+y^{\prime 2}{)}^3}(对应参数方程2)\text{\hspace{1em}(4)}$$

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

$$k^{\prime \prime }=\frac{d{k}^{\prime }}{ds}=\frac{d\frac{(y^{\prime \prime \prime }{x}^{\prime }-x^{\prime \prime \prime }{y}^{\prime })(x^{\prime 2}+y^{\prime 2})-3(x^{\prime }{x}^{\prime \prime }+y^{\prime }{y}^{\prime \prime })(y^{\prime \prime }{x}^{\prime }-x^{\prime \prime }{y}^{\prime })}{(x^{\prime 2}+y^{\prime 2}{)}^3}}{ds}$$

算不下去了。。。。。

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

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*ddx
    norm_square = dx*dx+dy*dy
    norm = sqrt(norm_square)
    b = norm*norm_square
    return a/b

def ComputeCurvatureDerivative(dx, ddx, dddx, dy, ddy, dddy):
    a = dx*ddy-dy*ddx
    b = dx*dddy-dy*dddx
    c = dx*ddx+dy*ddy
    d = dx*dx+dy*dy
    return (b*d-3.0*a*c)/(d*d*d)


def ComputeCurvature(dy, ddy):
    a = ddy 
    norm_square = 1+dy*dy
    norm = sqrt(norm_square)
    b = norm*norm_square
    return a/b

def ComputeCurvatureDerivative(dy, ddy, dddy):
    a = ddy
    b = dddy
    c = dy*ddy
    d = 1.0+dy*dy
    return (b*d-3.0*a*c)/(d*d*d)

END
by windSeS
2020.9.6