﻿ 三维离散曲线曲率挠率的微中心差分算法
«上一篇
 文章快速检索 高级检索

 智能系统学报  2019, Vol. 14 Issue (1): 194-206  DOI: 10.11992/tis.201802008 0

引用本文

MU Shengpeng, LI Hongjun, LI Shilin. An algorithm for estimating curvature and torsion of discrete curve in three-dimensional space based on microcentral difference[J]. CAAI Transactions on Intelligent Systems, 2019, 14(1): 194-206. DOI: 10.11992/tis.201802008.

文章历史

An algorithm for estimating curvature and torsion of discrete curve in three-dimensional space based on microcentral difference
MU Shengpeng , LI Hongjun , LI Shilin
College of Science, Beijing Forestry University, Beijing 100083, China
Abstract: The curvature and torsion of a 3D discrete curve reflect the degrees of its bending and distortion. To calculate these quantities accurately, following the definition of the derivative of the continuous curve, a microcentral difference algorithm, which is an extension to the one-side difference algorithm, is proposed based on the smoothing of the difference quotient. Compared with the one-side difference algorithm, the microcentral difference algorithm fails to prolong the running time but it remarkably improves the calculation accuracy. Several experiments are conducted by uniform sampling from six continuous curves, which are then compared with the five traditional algorithms of curvature and torsion. The experimental results are analyzed from three aspects: the influence of the sampling density on the accuracy of the algorithm, the efficiency of calculation, and the anti-noise performance. The experimental results show the good performance of the proposed microcentral difference algorithm.
Key words: curvature    torsion    algorithm comparison    discrete curve    microcentral difference algorithm    discrete geometry method    three-dimensional space    difference quotient    uniformly sampling

1 离散曲线的曲率挠率算法

 $\kappa (t) = \;\frac{{\left| {{ r}'(t) \times { r}''(t)} \right|}}{{{{\left| {{ r}'(t)} \right|}^3}}}\;$ (1)

 $\tau (t) = \frac{{({ r}'(t), { r}''(t), { r}'''(t))}}{{{\rm{|}}{ r}{\rm{'(}}t{\rm{)}} \times { r}''(t){{\rm{|}}^2}}}$ (2)

1.1 多项式曲线拟合法

2005年，Cazals等[16]提出了一种基于二次多项式曲线拟合(polynomial curve fitting，PCF)的方法来计算离散平面曲线在给定点 ${P_i}$ 的曲率、挠率。该方法对点 ${P_i}$ 及其邻域 $\{ {P_k}|i - n \leqslant k \leqslant i + n\}$ ( $n$ 为点 ${P_i}$ 一侧选取的邻近点个数)，在以 ${P_i}$ 为原点的局部坐标系中拟合一个二次多项式曲线。二次多项式曲线的数学模型为

 $f(x) = {A_0} + {A_1}x + {A_2}{x^2}$

1.2 B样条曲线拟合法

2017年，赵夫群等[17]用4次B样条曲线拟合(B-spline curve fitting，BSCF)方法将空间离散曲线拟合成光滑的空间曲线，从而以此来计算离散曲线的曲率和挠率。B样条曲线的定义为：给定 $n = m + k + 1$ 个顶点，则可以定义 $m + 1$ $k$ 次参数曲线第 $i$ 段B样条曲线函数可表示为

 $Q(t) = \displaystyle\sum\limits_{l = 0}^k {{P_{i + l}}{N_{l, k}}(t)}$

 ${N_{l, k}}(t) = \displaystyle\frac{1}{{k!}}\sum\limits_{j = 0}^{k - l} {{{( - 1)}^j}C_{k + 1}^j{{(t + k - l - j)}^k}}$

1.3 投影法

Mardia在1999年提出一种基于投影(projection)的算法[20]，对于三维空间中的离散曲线，将点投影到单位切线方向上。对于由 $n$ 个点构成的离散曲线，其点 ${P_i}$ 处的单位切向量定义为

 ${{ T}_i}{\rm{ = }}\frac{{{P_i}{P_{i + 1}}}}{{\left| {{P_i}{P_{i + 1}}} \right|}}, \quad i = 3, 4, \cdots , n - 2$

 ${Y_k} = \left( {\frac{{{{ T}_k} \cdot {{ N}_i}}}{{{{ T}_k} \cdot {{ T}_i}}}, \frac{{{{ T}_k} \cdot {{ B}_i}}}{{{{ T}_k} \cdot {{ T}_i}}}} \right),\quad k = i - 2, \cdot \cdot \cdot , i + 2$

 $\kappa ({P_i}) = \frac{{\left\| {{Y_i} - {Y_{i - 1}}} \right\|}}{{\left\| {{P_i} - {P_{i - 1}}} \right\|}}$
 $\tau ({P_i}) = \frac{{2a}}{{{{({b^2} + 1)}^{3/2}}}}\kappa ({P_i})$
1.4 单侧差分法

2016年Blankenburg等[21]通过利用一般参数方程的曲线导数的有限差商逼近，直接从离散点中采用单侧差分(one-sided difference，OSD)逼近计算出空间离散曲线在给定点 ${P_i}$ 的一阶差商、二阶差商和三阶差商，然后代入曲率和挠率的解析计算公式进行求解。在该方法中，定义差分为

 $\Delta {f_k}(t) = f(t + k\Delta t) - f(t + (k - 1)\Delta t)\;, \;k = 1, 2, 3$

 $f'(t) = \mathop {\lim }\limits_{\Delta t \to 0} \frac{{f(t{\rm{ + }}\Delta t) - f(t)}}{{\Delta t}}{\rm{ = }}\mathop {\lim }\limits_{\Delta t \to 0} \frac{{\Delta {f_1}(t)}}{{\Delta t}}$
 $\begin{gathered} f''(t) = \mathop {\lim }\limits_{\Delta t \to 0} \frac{{f(t{\rm{ + 2}}\Delta t) - 2f(t{\rm{ + }}\Delta t){\rm{ + }}f(t)}}{{\Delta {t^2}}} {\rm{ = }}\\ \;\;\;\;\;\;\;\;\mathop {\lim }\limits_{\Delta t \to 0} (\frac{{\Delta {f_2}(t)}}{{\Delta {t^2}}} - \frac{{\Delta {f_1}(t)}}{{\Delta {t^2}}}) \\ \end{gathered}$

 $f'''(t)\; = \;\;\mathop {\lim }\limits_{\Delta t \to 0} (\frac{{\Delta {f_3}(t)}}{{\Delta {t^3}}} - 2\frac{{\Delta {f_2}(t)}}{{\Delta {t^3}}} + \frac{{\Delta {f_1}(t)}}{{\Delta {t^3}}})$

1.5 离散几何法

An等[22]在2011年提出一种称为离散几何法(discrete geometry)的方法，该方法基于对称邻域的最小二乘的导数计算过程，可以归结为：利用离散导数的定义，采用约束最优化问题的求解方法来计算指定点在其对称邻域内的离散导数。

 ${y'_i} = \displaystyle\frac{{\displaystyle\sum\limits_{k = i - {n_1}}^{i + {n_1}} {({x_k} - {x_i})(} {y_k} - {y_i})}}{{\displaystyle\sum\limits_{k = i - {n_1}}^{i + {n_1}} {{{({x_k} - {x_i})}^2}} }}$

 ${y_i}^{\prime \prime } = \displaystyle\frac{{\displaystyle\sum\limits_{k = i - {n_2}}^{i + {n_2}} {({x_k} - {x_i})(} {{y'}_k} - {{y'}_i})}}{{\displaystyle\sum\limits_{k = i - {n_2}}^{i + {n_2}} {{{({x_k} - {x_i})}^2}} }}$

 $y_i^{(n)} = \frac{{\displaystyle\sum\limits_{k = i - n}^{i + n} {({x_k} - {x_i})(} y_k^{(n - 1)} - y_i^{(n - 1)})}}{{\displaystyle\sum\limits_{k = i - n}^{i + n} {{{({x_k} - {x_i})}^2}} }}$

 ${s_i} = \frac{{\displaystyle\sum\limits_{k = 1}^{i - 1} {\left| {{p_{k + 1}} - {p_k}} \right|} }}{{\displaystyle\sum\limits_{k = 1}^{n - 1} {\left| {{p_{k + 1}} - {p_k}} \right|} }},\quad 2 \leqslant i \leqslant n$

 ${\psi _d} = \{ {P_i}\, {\rm{\} = \{ }}\left( {x({s_i}), y({s_i}), z({s_i})} \right)\}$

 $\left\{ \begin{gathered} x'({s_i}) = \frac{{\displaystyle\sum\limits_{j = i - {n_1}}^{i + {n_1}} {({s_j} - {s_i})({x_j} - {x_i})} }}{{\displaystyle\sum\limits_{j = i - {n_1}}^{i + {n_1}} {{{({s_j} - {s_i})}^2}} }} \\ y'({s_i}) = \frac{{\displaystyle\sum\limits_{j = i - {n_1}}^{i + {n_1}} {({s_j} - {s_i})({y_j} - {y_i})} }}{{\displaystyle\sum\limits_{j = i - {n_1}}^{i + {n_1}} {{{({s_j} - {s_i})}^2}} }} \\ z'({s_i}) = \frac{{\displaystyle\sum\limits_{j = i - {n_1}}^{i + {n_1}} {({s_j} - {s_i})({z_j} - {z_i})} }}{{\displaystyle\sum\limits_{j = i - {n_1}}^{i + {n_1}} {{{({s_j} - {s_i})}^2}} }} \\ \end{gathered} \right.$

1.6 微中心差分法

 $f_{\rm q}'(t) = \frac{1}{{{k_1}}}\sum\limits_{i = 1}^{{k_1}} {\frac{{f(t) - f(t - i\Delta t)}}{{\Delta t}}}$

 $f_{\rm h}'(t) = \displaystyle\frac{1}{{{k_1}}}\sum\limits_{i = 1}^{{k_1}} {\frac{{f(t + i\Delta t) - f(t)}}{{\Delta t}}}$

 $\begin{gathered} \mathop f\nolimits' (t) = \frac{1}{2}\left( {f_{\rm q}'(t) + f_{\rm h}'(t)} \right) {\rm{ = }} \frac{1}{{2{k_1}}}\sum\limits_{i = 1}^{{k_1}} {\frac{{f(t + i\Delta t) - f(t - i\Delta t)}}{{\Delta t}}} \end{gathered}$ (3)

 $f''_{\rm q}(t) = \frac{1}{{{k_2}}}\sum\limits_{i = 1}^{{k_2}} {\frac{{\mathop f\nolimits' (t) - \mathop f\nolimits' (t - i\Delta t)}}{{\Delta t}}}$

 $f''_{\rm h} (t) = \frac{1}{{{k_2}}}\sum\limits_{i = 1}^{{k_2}} {\frac{{\mathop f\nolimits' (t + i\Delta t) - \mathop f\nolimits' (t)}}{{\Delta t}}}$

 $\begin{gathered} f'' (t) \!= \!\frac{1}{2}\left( {f''_{\rm q} (t) \!+\! f''_{\rm h} (t)} \right) \!= \! \frac{1}{{2{k_2}\Delta t}}\left( {\mathop f\nolimits' (t\! +\! {k_2}\Delta t) \!-\! \mathop f\nolimits' (t - {k_2}\Delta t)} \right) \\ \end{gathered}$ (4)

 $f'''_{\rm q}(t) = \frac{1}{{{k_3}}}\sum\limits_{i = 1}^{{k_3}} {\frac{{ f''(t) - f'' (t - i\Delta t)}}{{\Delta t}}}$

 $f'''_{\rm h}(t) = \frac{1}{{{k_3}}}\sum\limits_{i = 1}^{{k_3}} {\frac{{ f'' (t + i\Delta t) - f'' (t)}}{{\Delta t}}}$

 $\begin{gathered} f''' (t) \!\!=\!\! \frac{1}{2}\left( {f_{\rm q}'''(t) \!+\! f_{\rm h}'''(t)} \right) \!\!= \!\! \frac{1}{{2{k_3}\Delta t}}\left( {f'' (t \!+\! {k_2}\Delta t) - f'' (t \!-\!\!{k_2}\Delta t)} \right) \end{gathered}$ (5)

2 实验比较

2.1 实验数据

 Download: 图 1 用于实验的6条空间曲线(线颜色与挠率值对应) Fig. 1 Six 3D curves used in the experiment (The color is corresponding to torsion value)

1)三次挠曲线：

 $r(t) = (2t, {t^2}, {t^3}), \quad t \in [ - 2, 2]$

2) 多变式内接弹簧线：

 $r(t) = (2\cos (3t) + \cos (2t), 2\sin (3t) + \sin (2t), t), \quad t \in [0, 2{\text π} ]$

3) 循环曲线：

 $r(t) = (\cos (t), \sin (t), 2\sin (t)),\quad t \in [0, 2{\text π} ]$

4) Viviani曲线：

 $r(t) = ({\cos ^2}(t), \cos (t)\sin (t), \sin (t)),\quad t \in [0, 2{\text π} ]$

5) Clelia曲线：

 $r(t) = (\cos (2t)\cos (t), \cos (2t)\sin (t), \sin (2t)), \quad t \in [0, 2{\text π} ]$

6) 环面螺线：

 $r(t) = ((2 + \cos (t))\cos (2t), (2 + \cos (t))\sin (2t), t),\quad t \in [0, 2{\text π} ]$

2.2 误差分析及时间效率比较

 Download: 图 3 多变式内接弹簧线误差对比实验 Fig. 3 Error comparison experiment of changeable inner spring curve
 Download: 图 5 Viviani曲线误差对比实验 Fig. 5 Error comparison experiment of Viviani curve
 Download: 图 6 Clelia曲线误差对比实验 Fig. 6 Error comparison experiment of Clelia curve

 Download: 图 2 三次挠曲线误差对比实验 Fig. 2 Error comparison experiment of cubic deflection curve
 Download: 图 7 环面螺线误差对比实验 Fig. 7 Error comparison experiment of toroidal spiral curve

2.3 噪声数据实验

6种算法计算的曲率最大相对误差和曲率均方根误差分别如图8(a)(b)所示，从图中可以看出：噪声水平越高，算法的误差变化越大；各算法在噪声水平不大于0.01时，均能较准确地计算离散曲线的曲率值；噪声水平大于0.01时，其误差明显变大，尤其是单侧差分法和投影法，而本文的方法误差变化浮动相对较小。

 Download: 图 8 Clelia曲线的不同噪声水平下的MSE和RMSE( $n = 200$ ) Fig. 8 Max relative error and root mean square error of curvatures of Clelia curve under different noise levels ( $n = 200$ )

6种算法计算的挠率最大相对误差和均方根误差分别如图8(c)(d)所示，从图中可以看出：含噪声的离散曲线在计算挠率时变得更加敏感，当噪声水平大于0.001时，6种算法得到的挠率值与真实值就已有了明显的偏差，并且随着噪声水平变大，挠率的最大相对误差和均方根误差相比曲率也会更快地增大。与其他5种算法相比，本文的方法误差变化浮动相对较小。

3 结束语

1)根据2016年Blankenburg等提出的单侧差分方法和导数定义，结合差商平滑策略，设计了微中心差分法，该方法能够有效计算离散曲线的曲率和挠率。从整个计算过程可以发现，本文的方法不依赖于任何曲线模型，从离散曲线中直接计算所得，从而对各类离散曲线均具有较好的普适性。

2)概述了多项式曲线拟合法、B样条曲线拟合法、投影法、单侧差分法和离散几何法5种典型算法；详细阐述了本文提出的微中心差分法的基本思想和计算公式。采用6条离散曲线数据进行了不同算法的对比实验，实验表明，本文提出的微中心差分法的平均精度最高。

 [1] 李杰, 彭双春, 安宏雷, 等. 基于微分几何与李群的无人机编队会合方法[J]. 国防科技大学学报, 2013, 35(6): 157-164. LI Jie, PENG Shuangchun, AN Honglei, et al. UAVs formation rendezvous method based on differential geometry and Lie group[J]. Journal of national university of defense technology, 2013, 35(6): 157-164. DOI:10.3969/j.issn.1001-2486.2013.06.027 (0) [2] 单洪春. 基于空间特性的弹道特征抽取[J]. 飞行器测控学报, 2012, 31(2): 69-72. SHAN Hongchun. Trajectory characteristics extraction based on space property[J]. Journal of spacecraft TT&C technology, 2012, 31(2): 69-72. (0) [3] 葛婷, 符锌砂, 李海峰, 等. 公路三维线形设计及约束建模[J]. 华南理工大学学报 (自然科学版), 2016, 44(8): 91-97, 105. GE Ting, FU Xinsha, LI Haifeng, et al. Three-dimensional highway alignment design and constraints modeling[J]. Journal of South China university of technology (natural science edition), 2016, 44(8): 91-97, 105. (0) [4] 潘登, 郑应平. 路径约束条件下车辆行为的时空演化模型[J]. 物理学报, 2015, 64(7): 078902. PAN Deng, ZHENG Yingping. Spatiotemporal evolution model of vehicular movement behavior under path constraints[J]. Acta physica sinica, 2015, 64(7): 078902. (0) [5] PATRIKALAKIS N M, MAEKAWA T. Shape interrogation for computer aided design and manufacturing[M]. Berlin, Heidelberg: Springer, 2002. (0) [6] MOKHTARIAN F, MACKWORTH A K. A theory of multiscale, curvature-based shape representation for planar curves[J]. IEEE transactions on pattern analysis and machine intelligence, 1992, 14(8): 789-805. DOI:10.1109/34.149591 (0) [7] 袁修久, 王胜勇, 刘欣, 等. 推广的三维L-系统及在树木模拟中的应用[J]. 系统仿真学报, 2011, 23(11): 2308-2311. YUAN Xiujiu, WANG Shengyong, LIU Xin, et al. Improvement of three-dimensional L-system and its application to trees modeling[J]. Journal of system simulation, 2011, 23(11): 2308-2311. (0) [8] 张冠军, 朱翔, 李天匀. 基于级数变换法的椭圆柱壳受迫振动分析[J]. 哈尔滨工程大学学报, 2017, 38(4): 506-513. ZHANG Guanjun, ZHU Xiang, LI Tianyun. Forced vibration analysis of elliptic cylindrical shell based on the series transformation method[J]. Journal of Harbin engineering university, 2017, 38(4): 506-513. (0) [9] 孟道骥, 梁科. 微分几何[M]. 2版. 北京: 科学出版社, 2004: 6–28. MENG Daoji, LIANG Ke. Differential geometry[M]. 2nd ed. Beijing: Science Press, 2004: 6–28. (0) [10] DO CARMO M P. Differential geometry of curves and surfaces[M]. Englewood Cliffs, New Jersey: Prentice Hall, 1976: 175–178. (0) [11] BRAY H L, JAUREGUI J L. On curves with nonnegative torsion[J]. Archiv der mathematik, 2015, 104(6): 561-575. DOI:10.1007/s00013-015-0767-0 (0) [12] LIU Huili. Curves in three dimensional riemannian space forms[J]. Results in mathematics, 2014, 66(3/4): 469-480. (0) [13] MASSER D, ZANNIER U. Torsion points on families of products of elliptic curves[J]. Advances in mathematics, 2014, 259: 116-133. DOI:10.1016/j.aim.2014.03.016 (0) [14] 郑长波, 刘会立. 三维欧氏空间中的几类特殊球面曲线[J]. 浙江大学学报 (理学版), 2013, 40(2): 127-130. ZHENG Changbo, LIU Huili. Some special spherical curves in Euclidean 3-space[J]. Journal of Zhejiang university (science edition), 2013, 40(2): 127-130. (0) [15] 袁媛, 李静, 刘会立. 三维欧氏空间中的特殊曲线[J]. 东北大学学报 (自然科学版), 2017, 38(11): 1669-1672. YUAN Yuan, LI Jing, LIU Huili. Special curves in 3-dimensional euclidean space[J]. Journal of Northeastern university (natural science), 2017, 38(11): 1669-1672. (0) [16] CAZALS F, POUGET M. Estimating differential quantities using polynomial fitting of osculating jets[J]. Computer aided geometric design, 2005, 22(2): 121-146. DOI:10.1016/j.cagd.2004.09.004 (0) [17] 赵夫群, 周明全. 一种有效的秦俑碎块匹配算法[J]. 计算机系统应用, 2017, 26(3): 198-203. ZHAO Fuqun, ZHOU Mingquan. Effective blocks matching algorithm of terracotta warrior[J]. Computer systems and applications, 2017, 26(3): 198-203. (0) [18] 喻德生, 程程. 基于离散曲率的三次均匀B样条的局部光顺算法[J]. 浙江大学学报 (理学版), 2011, 38(5): 511-517. YU Desheng, CHENG Cheng. On local fairing algorithm for cubic B-spline with discrete curvature[J]. Journal of Zhejiang university (science edition), 2011, 38(5): 511-517. (0) [19] KEHTARNAVAZ N, DEFIGUEIREDO R J P. A 3-D contour segmentation scheme based on curvature and torsion[J]. IEEE transactions on pattern analysis and machine intelligence, 1988, 10(5): 707-713. DOI:10.1109/34.6780 (0) [20] MARDIA K V. Estimation of torsion[J]. Journal of applied statistics, 1999, 26(3): 373-381. DOI:10.1080/02664769922476 (0) [21] BLANKENBURG C, DAUL C, OHSER J. Parameter free torsion estimation of curves in 3D images[C]//2016 IEEE International Conference on Image Processing. Phoenix, AZ, USA, 2016: 1081–1085. (0) [22] AN Yi, SHAO Cheng, WANG Xiaoliang, et al. Geometric properties estimation from discrete curves using discrete derivatives[J]. Computers & graphics, 2011, 35(4): 916-930. (0) [23] 陈维桓. 微分几何[M]. 北京: 北京大学出版社, 2006: 23–47. (0) [24] MA Weiyin, KRUTH J P. Parameterization of randomly measured points for least squares fitting of B-spline curves and surfaces[J]. Computer-aided design, 1995, 27(9): 663-675. DOI:10.1016/0010-4485(94)00018-9 (0) [25] 周志华. 机器学习[M]. 北京: 清华大学出版社, 2016: 28−37. (0) [26] 刘伟, 王雅静, 陈文钢, 等. 动态光散射反演算法的评价指标[J]. 光学学报, 2015, 35(S1): 129001. LIU Wei, WANG Yajing, CHEN Wengang, et al. Evaluation criteria of inversion algorithm for dynamic light scattering[J]. Acta optica sinica, 2015, 35(S1): 129001. (0)