第14卷第1期 智能系统学报 Vol.14 No.I 2019年1月 CAAI Transactions on Intelligent Systems Jan.2019 D0:10.11992/tis.201802008 网络出版地址：http:/kns.cnki.net/kcms/detail/23.1538.TP.20180927.1124.006html 三维离散曲线曲率挠率的微中心差分算法 慕生鹏，李红军，李世林 (北京林业大学理学院，北京100083) 摘要：曲率和挠率是描述三维空间离散曲线的弯曲和扭曲程度的两个微分量。为了准确计算这两个微分量， 从连续曲线的导数定义出发，提出微中心差分算法进行三维空间离散曲线的曲率和挠率计算。该算法基于差 商平滑策略实现对单侧差分算法的一个有效扩展。与单侧差分算法相比，微中心差分算法不增加算法执行时 间，但在计算精度方面有显著提升。实验分析是通过6条曲线的均匀采样获取离散曲线数据，与5种常用的曲 率和挠率计算算法相比较，对这6种算法从采样密度对算法精度的影响、计算效率和抗噪声性能这3个方面进 行了对比分析。实验结果表明，微中心差分算法总体效果最好。 关键词：曲率；挠率；算法比较：离散曲线；微中心差分法；离散几何法；三维空间；差商；均匀采样 中图分类号：TP311文献标志码：A文章编号：1673-4785(2019)01-0194-13 中文引用格式：慕生鹏，李红军，李世林.三维离散曲线曲率挠率的微中心差分算法J.智能系统学报，2019,14(1)： 194-206. 英文引用格式：MU Shengpeng,LI Hongjun,LI Shilin.An algorithm for estimating curvature and torsion of discrete curve in three-dimensional space based on microcentral differencelJ].CAAI transactions on intelligent systems,2019,14(1):194-206. 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 differ- ence 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 tor- sion.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 per- formance of the proposed microcentral difference algorithm. Keywords:curvature;torsion;algorithm comparison;discrete curve;microcentral difference algorithm;discrete geo- metry method;three-dimensional space;difference quotient;uniformly sampling 在人工智能算法设计与自动控制的相关研究 视觉研究，以及植物生长模拟四、结构工程分析图 中，无人机路径规划四，弹道分析P,公路线性设计 等问题的解决都可以基于离散曲线的曲率和挠率 路径约束下的车辆行为研究，几何处理和机器 的分析。曲率和挠率是空间曲线在固有运动下的 收稿日期：2018-02-05.网络出版日期：2018-09-29. 不变量，直接描述了曲线在一点邻近的形状。 基金项目：国家自然科学基金项目(61372190)：中央高校基本 科研业务费专项资金项目(2015 ZCQ-LY-01). 其中，曲率揭示曲线在所在平面的弯曲程度，挠 通信作者：李红军.E-mail:lihongjun69@bjfu.edu.cn. 率则刻画曲线离开既定平面的扭曲程度。也就是
DOI: 10.11992/tis.201802008 网络出版地址: http://kns.cnki.net/kcms/detail/23.1538.TP.20180927.1124.006.html 三维离散曲线曲率挠率的微中心差分算法 慕生鹏，李红军，李世林 （北京林业大学 理学院，北京 100083） 摘 要：曲率和挠率是描述三维空间离散曲线的弯曲和扭曲程度的两个微分量。为了准确计算这两个微分量， 从连续曲线的导数定义出发，提出微中心差分算法进行三维空间离散曲线的曲率和挠率计算。该算法基于差 商平滑策略实现对单侧差分算法的一个有效扩展。与单侧差分算法相比，微中心差分算法不增加算法执行时 间，但在计算精度方面有显著提升。实验分析是通过 6 条曲线的均匀采样获取离散曲线数据，与 5 种常用的曲 率和挠率计算算法相比较，对这 6 种算法从采样密度对算法精度的影响、计算效率和抗噪声性能这 3 个方面进 行了对比分析。实验结果表明，微中心差分算法总体效果最好。 关键词：曲率；挠率；算法比较；离散曲线；微中心差分法；离散几何法；三维空间；差商；均匀采样 中图分类号：TP311 文献标志码：A 文章编号：1673−4785(2019)01−0194−13 中文引用格式：慕生鹏, 李红军, 李世林. 三维离散曲线曲率挠率的微中心差分算法[J]. 智能系统学报, 2019, 14(1): 194–206. 英文引用格式：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. 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. Keywords: curvature; torsion; algorithm comparison; discrete curve; microcentral difference algorithm; discrete geometry method; three-dimensional space; difference quotient; uniformly sampling 在人工智能算法设计与自动控制的相关研究 中，无人机路径规划[1] ，弹道分析[2] ，公路线性设计[3] ， 路径约束下的车辆行为研究[4] ，几何处理[5]和机器 视觉研究[6] ，以及植物生长模拟[7] 、结构工程分析[8] 等问题的解决都可以基于离散曲线的曲率和挠率 的分析。曲率和挠率是空间曲线在固有运动下的 不变量[9] ，直接描述了曲线在一点邻近的形状。 其中，曲率揭示曲线在所在平面的弯曲程度，挠 率则刻画曲线离开既定平面的扭曲程度。也就是 收稿日期：2018−02−05. 网络出版日期：2018−09−29. 基金项目：国家自然科学基金项目 (61372190)；中央高校基本 科研业务费专项资金项目 (2015ZCQ-LY-01). 通信作者：李红军. E-mail：lihongjun69@bjfu.edu.cn. 第 14 卷第 1 期 智 能 系 统 学 报 Vol.14 No.1 2019 年 1 月 CAAI Transactions on Intelligent Systems Jan. 2019
第1期 慕生鹏，等：三维离散曲线曲率挠率的微中心差分算法 ·195· 说，三维空间中的曲线是可以通过一根直线弯曲 ),()中要求3个函数x()、y()、()都是连续 (曲率)和扭曲（挠率）得到的。显然，曲率和挠率 函数，且都3阶可导。 能完全刻画曲线的局部行为。因此，有研究者 但是对于没有方程的三维离散曲线的曲率挠 根据曲线的曲率挠率特性，对空间中的曲线进行 率计算，研究人员提出的各种方法中有代表性的 分析1)、分类4或判断轨迹线的形状1。对于 方法有：多项式曲线拟合法、B样条曲线拟合法、 有解析或者参数方程的曲线，经典微分几何中有 投影法、单侧差分法、离散几何法5种。本文在 成熟方法计算曲线上任意一点的曲率和挠率。 分析这5种算法的基础上，提出了微中心差分法。 因为近些年应用日益广泛的离散空间曲线没 11多项式曲线拟合法 有解析方程，经典微分几何中的求解方法不能直 2005年，Cazals等16提出了一种基于二次多 接应用于离散曲线的情形，所以三维空间中离散 项式曲线拟合(polynomial curve fitting,PCF)的方 曲线的曲率和挠率计算成为近些年的热点研究问 法来计算离散平面曲线在给定点P,的曲率、挠 题。此外，在实际应用中，受人工、仪器、环境等 率。该方法对点P,及其邻域(Pi-n≤k≤i+n(n为 诸多因素影响，得到的离散曲线的数据（一个序 点P,一侧选取的邻近点个数)，在以P:为原点的局 列的三维坐标点)常常存在一定的偏差或者噪 部坐标系中拟合一个二次多项式曲线。二次多项 声，因此如何使算法更加准确和鲁棒，成为三维 式曲线的数学模型为 空间离散曲线的曲率和挠率计算的难点。 f(x)=Ao+Ax+Azx2 近些年提出的关于三维空间离散曲线的曲率 通过最小二乘拟合，可求解该二次多项式曲 和挠率计算方法大致可以分为2类。1)基于曲线 线及其系数A,(i=0,1,2),利用多项式系数A可计算 拟合的方法，这类算法主要是对给定点及其近邻 出点P的曲率k(P)。类似的可推广拟合三次或三 点拟合一段光滑曲线，然后根据这段解析曲线的 次以上的多项式曲线，从而利用式(1)、式(2)计 微分几何公式计算给定点的曲率和挠率。其本质 算出三维空间离散曲线的曲率、挠率。 就是用拟合曲线的曲率和挠率来估算离散曲线的 1.2B样条曲线拟合法 曲率和挠率。这类算法的典型算法有多项式拟合 2017年，赵夫群等1用4次B样条曲线拟合 法(polynomial curve fitting)Io和B样条曲线拟合 (B-spline curve fitting,BSCF)方法将空间离散曲线 法(B-spline curve fitting).19.2)基于离散结构的 拟合成光滑的空间曲线，从而以此来计算离散曲 方法，这类算法无需对离散曲线进行光滑曲线的 线的曲率和挠率。B样条曲线的定义为：给定 拟合，直接利用离散点之间的距离和内接多边形 n=m+k+1个顶点，则可以定义m+1段k次参数曲 等离散结构，运用差分形式或者离散逼近等策略 线第段B样条曲线函数可表示为 估算出离散曲线上各个点的曲率和挠率。基于离散 结构的方法中比较典型的有投影法(projection)P、 e0=∑PNad 单侧差分法(one-sided difference)2和离散几何法 =0 (discrete geometry)2。本文根据导数定义以及单 式中：i=0,1,…,n;{P40,P41,…,P4t}为定义第i段 侧差分法的设计思路，融合差商平滑策略，提出 曲线特征多边形的k+1个顶点；Nu(①为B样条基 了微中心差分(microcentral difference)算法。 底函数，表示为 t-l 1离散曲线的曲率挠率算法 N(0= 1.(-1yct+k-l- j=0 经典微分几何中，设曲线C的参数方程为 对于离散曲线上的每一个离散点P,选取点 r=r(0=(x(),y(0,z()(a≤t≤)，若一阶导函数 P,及其邻域{Pli-n≤k≤i+(n为点P,一侧选取的 (0在区间[α，B上恒不为零，即参数化曲线r()是 邻近点个数)内的点进行样条曲线拟合。类似的 正则的，则曲率k(0)的计算公式2到为 方法有：喻德生等8采用的3次均匀B样条进行 曲线拟合，Kehtarnavaz等I采用的5次B样条进 K)= r(t)xr"(t) (1) r()l' 行曲线拟合，从而利用经典微分几何的式(1)、式 挠率x()的计算公式为 (2)求解方法对离散曲线进行曲率、挠率的计算。 )=.r).r() 1.3投影法 Ir(t)xr"(t) (2) Mardia在1999年提出一种基于投影(projec 式(1)、式(2)中的符号“×”表示叉积；式(2)中 tion)的算法20，对于三维空间中的离散曲线，将 的分子表示3个微分向量的混合积。()=(x(), 点投影到单位切线方向上。对于由n个点构成的
说，三维空间中的曲线是可以通过一根直线弯曲 (曲率) 和扭曲 (挠率) 得到的。显然，曲率和挠率 能完全刻画曲线的局部行为[10]。因此，有研究者 根据曲线的曲率挠率特性，对空间中的曲线进行 分析[11-13] 、分类[14]或判断轨迹线的形状[15]。对于 有解析或者参数方程的曲线，经典微分几何中有 成熟方法计算曲线上任意一点的曲率和挠率。 因为近些年应用日益广泛的离散空间曲线没 有解析方程，经典微分几何中的求解方法不能直 接应用于离散曲线的情形，所以三维空间中离散 曲线的曲率和挠率计算成为近些年的热点研究问 题。此外，在实际应用中，受人工、仪器、环境等 诸多因素影响，得到的离散曲线的数据 (一个序 列的三维坐标点) 常常存在一定的偏差或者噪 声，因此如何使算法更加准确和鲁棒，成为三维 空间离散曲线的曲率和挠率计算的难点。 近些年提出的关于三维空间离散曲线的曲率 和挠率计算方法大致可以分为 2 类。1) 基于曲线 拟合的方法，这类算法主要是对给定点及其近邻 点拟合一段光滑曲线，然后根据这段解析曲线的 微分几何公式计算给定点的曲率和挠率。其本质 就是用拟合曲线的曲率和挠率来估算离散曲线的 曲率和挠率。这类算法的典型算法有多项式拟合 法 (polynomial curve fitting)[16]和 B 样条曲线拟合 法 (B-spline curve fitting)[17-19]。2) 基于离散结构的 方法，这类算法无需对离散曲线进行光滑曲线的 拟合，直接利用离散点之间的距离和内接多边形 等离散结构，运用差分形式或者离散逼近等策略 估算出离散曲线上各个点的曲率和挠率。基于离散 结构的方法中比较典型的有投影法 (projection)[20] 、 单侧差分法 (one-sided difference)[21]和离散几何法 (discrete geometry)[22]。本文根据导数定义以及单 侧差分法的设计思路，融合差商平滑策略，提出 了微中心差分 (microcentral difference) 算法。 1 离散曲线的曲率挠率算法 C r = r(t) = (x (t), y (t),z(t)) (α ⩽ t ⩽ β) r ′ (t) [ α, β] r(t) κ(t) 经典微分几何中，设曲线 的参数方程为 ，若一阶导函数 在区间 上恒不为零，即参数化曲线 是 正则的，则曲率 的计算公式[23]为 κ(t) = |r ′ (t)× r ′′(t)| |r ′ (t)| 3 (1) 挠率τ(t) 的计算公式[16]为 τ(t) = (r ′ (t),r ′′(t),r ′′′(t)) |r ′ (t)× r ′′(t)| 2 (2) 式 (1)、式 (2) 中的符号“×”表示叉积；式 (2) 中 的分子表示 3 个微分向量的混合积。r(t)=(x(t)， y(t)，z(t)) 中要求 3 个函数 x(t)、y(t)、z(t) 都是连续 函数，且都 3 阶可导。 但是对于没有方程的三维离散曲线的曲率挠 率计算，研究人员提出的各种方法中有代表性的 方法有：多项式曲线拟合法、B 样条曲线拟合法、 投影法、单侧差分法、离散几何法 5 种。本文在 分析这 5 种算法的基础上，提出了微中心差分法。 1.1 多项式曲线拟合法 Pi Pi {Pk |i−n ⩽ k ⩽ i+n} n Pi Pi 2005 年，Cazals 等 [16]提出了一种基于二次多 项式曲线拟合 (polynomial curve fitting，PCF) 的方 法来计算离散平面曲线在给定点 的曲率、挠 率。该方法对点 及其邻域 ( 为 点 一侧选取的邻近点个数)，在以 为原点的局 部坐标系中拟合一个二次多项式曲线。二次多项 式曲线的数学模型为 f(x) = A0 + A1 x+ A2 x 2 Ai(i = 0,1,2) Ai Pi κ(Pi) 通过最小二乘拟合，可求解该二次多项式曲 线及其系数 ，利用多项式系数 可计算 出点 的曲率 。类似的可推广拟合三次或三 次以上的多项式曲线，从而利用式 (1)、式 (2) 计 算出三维空间离散曲线的曲率、挠率。 1.2 B 样条曲线拟合法 n = m+k+1 m+1 k i 2017 年，赵夫群等[17]用 4 次 B 样条曲线拟合 (B-spline curve fitting，BSCF) 方法将空间离散曲线 拟合成光滑的空间曲线，从而以此来计算离散曲 线的曲率和挠率。B 样条曲线的定义为：给定 个顶点，则可以定义 段 次参数曲 线第 段 B 样条曲线函数可表示为 Q(t) = ∑k l=0 Pi+lNl,k(t) i = 0,1,··· ,n {Pi+0,Pi+1,··· ,Pi+k} i k+1 Nl,k(t) 式中： ； 为定义第 段 曲线特征多边形的 个顶点； 为 B 样条基 底函数，表示为 Nl,k(t) = 1 k! ∑k−l j=0 (−1)jC j k+1 (t+k−l− j) k Pi Pi {Pk |i−n ⩽ k ⩽ i+n} n Pi 对于离散曲线上的每一个离散点 ，选取点 及其邻域 ( 为点 一侧选取的 邻近点个数) 内的点进行样条曲线拟合。类似的 方法有：喻德生等[18]采用的 3 次均匀 B 样条进行 曲线拟合，Kehtarnavaz 等 [19]采用的 5 次 B 样条进 行曲线拟合，从而利用经典微分几何的式 (1)、式 (2) 求解方法对离散曲线进行曲率、挠率的计算。 1.3 投影法 n Mardia 在 1999 年提出一种基于投影 (projection) 的算法[20] ，对于三维空间中的离散曲线，将 点投影到单位切线方向上。对于由 个点构成的 第 1 期 慕生鹏，等：三维离散曲线曲率挠率的微中心差分算法 ·195·
·196· 智能系统学报 第14卷 离散曲线，其点P处的单位切向量定义为 上述离散导数也称为点P的一阶离散导数。 T-PF i i=3,4,…,n-2 类似的在点P处的二阶离散导数为 t挡 容易计算点P处的单位法向量N,利用T和 -xy-y) N,做叉积得到副法向量B,从而点P处的局部投影 =i- y"= 平面离散曲线定义为 (T&-N:T&.Ba ￥=7TT k=i-2.…,i+2 同理，在点P处的n阶离散导数为 然后以点Y,为中心，拟合平面二次曲线 ar2+bx,从而三维空间离散曲线点P,的曲率k(P)、 （一0- 挠率x(P)分别为 y=E K(P)= lY;-Y IIP:-P-il 2a 上述各阶离散导数是对连续可微函数在点 t(P)= Γ(b2+)P) P处的各阶导数的逼近。对于三维空间中的离散 1.4单侧差分法 曲线，为了计算其离散导数，首先要对离散曲线 2016年Blankenburg等2通过利用一般参数 进行参数化。可采用累积弦长参数化方法2来确 方程的曲线导数的有限差商逼近，直接从离散点 定离散点的参数值。在该方法中，离散曲线山在各 中采用单侧差分(one-sided difference,.OSD)逼近 点P,=(化，y,乙1≤i≤)处的累积弦长参数s,可确定为 计算出空间离散曲线在给定点P的一阶差商、二 阶差商和三阶差商，然后代入曲率和挠率的解析 ∑Ip1-pl 计算公式进行求解。在该方法中，定义差分为 S:= 2≤i≤n △f)=ft+k△)-ft+(k-1)△)，k=1,2,3 Pu-p 则一阶导数、二阶导数分别为 k=1 并且s1=0,从而离散曲线a为 f(t)=lim f(t+△)-f( △t 2=0 a={P}=(x(s(s,z(s)} f0=mf0+2a）-2f+a+f0_ 式中：，∈{s1,2,,Sn:坐标函数x=x(s)、y=(s)和 42 =(s)就是弦长s的离散函数。因此，平面上的离 ,△f()△fi() mAe 散导数的计算方法可以转化为用于计算三维空间 42 同理，三阶导数为 中的离散曲线上各点的导数。从而，采用上述方 广0=四60-20+9 法，对三维空间中的离散曲线山上的点P,其邻域 △3 △3 为{Pi-m≤k≤i+n1(n1为点P一侧选取的邻近点 从而根据曲线曲率、挠率的微分定义，将上 个数)，则点P处的离散导数(x(s),y(s),x(s)分别为 述各阶导直接代入曲率、挠率的计算式(1)、式 t四 (2)即可求得各点处的曲率、挠率。 (si-si)(xi-xi) 1.5离散几何法 x(S)= An等2在2011年提出一种称为离散几何法 (discrete geometry)的方法，该方法基于对称邻域 - 的最小二乘的导数计算过程，可以归结为：利用 ∑5-sw-w 离散导数的定义，采用约束最优化问题的求解方 y(s)= 法来计算指定点在其对称邻域内的离散导数。 对于平面内的离散曲线，点P:=(xy)及其邻 域{Pi-n,≤k≤i+n1(m1为点P,一侧选取的邻近点 个数)，则点P处的离散导数为 2s-s6-a z(S)= ,(x-x)0y-) k=i-n y= u- (-x)2 上述离散导数也称为点P的一阶离散导数。 类似的可计算点P,的二阶、三阶离散导数。从而
离散曲线，其点 Pi处的单位切向量定义为 Ti= PiPi+1 |PiPi+1| , i = 3,4,··· ,n−2 Pi Ni Ti Ni Bi Pi 容易计算点 处的单位法向量 ，利用 和 做叉积得到副法向量 ，从而点 处的局部投影 平面离散曲线定义为 Yk = ( Tk · Ni Tk ·Ti , Tk · Bi Tk ·Ti ) , k = i−2,· · ·,i+2 Yi ax2 +bx Pi κ(Pi) τ(Pi) 然后以点 为中心，拟合平面二次曲线 ，从而三维空间离散曲线点 的曲率 、 挠率 分别为 κ(Pi) = ∥Yi −Yi−1∥ ∥Pi − Pi−1∥ τ(Pi) = 2a (b 2 +1)3/2 κ(Pi) 1.4 单侧差分法 Pi 2016 年 Blankenburg 等 [21]通过利用一般参数 方程的曲线导数的有限差商逼近，直接从离散点 中采用单侧差分 (one-sided difference，OSD) 逼近 计算出空间离散曲线在给定点 的一阶差商、二 阶差商和三阶差商，然后代入曲率和挠率的解析 计算公式进行求解。在该方法中，定义差分为 ∆fk(t) = f(t+k∆t)− f(t+(k−1)∆t) , k = 1,2,3 则一阶导数、二阶导数分别为 f ′ (t) = lim ∆t→0 f(t+∆t)− f(t) ∆t = lim ∆t→0 ∆f1(t) ∆t f ′′(t) = lim ∆t→0 f(t+2∆t)−2 f(t+∆t)+f(t) ∆t 2 = lim ∆t→0 ( ∆f2(t) ∆t 2 − ∆f1(t) ∆t 2 ) 同理，三阶导数为 f ′′′(t) = lim ∆t→0 ( ∆f3(t) ∆t 3 −2 ∆f2(t) ∆t 3 + ∆f1(t) ∆t 3 ) 从而根据曲线曲率、挠率的微分定义，将上 述各阶导直接代入曲率、挠率的计算式 (1)、式 (2) 即可求得各点处的曲率、挠率。 1.5 离散几何法 An 等 [22]在 2011 年提出一种称为离散几何法 (discrete geometry) 的方法，该方法基于对称邻域 的最小二乘的导数计算过程，可以归结为：利用 离散导数的定义，采用约束最优化问题的求解方 法来计算指定点在其对称邻域内的离散导数。 Pi = (xi , yi) {Pk |i−n1 ⩽ k ⩽ i+n1} n1 Pi Pi 对于平面内的离散曲线，点 及其邻 域 ( 为点 一侧选取的邻近点 个数)，则点 处的离散导数为 y ′ i = ∑i+n1 k=i−n1 (xk − xi)(yk −yi) ∑i+n1 k=i−n1 (xk − xi) 2 Pi Pi 上述离散导数也称为点 的一阶离散导数。 类似的在点 处的二阶离散导数为 yi ′′ = ∑i+n2 k=i−n2 (xk − xi)(y ′ k −y ′ i ) ∑i+n2 k=i−n2 (xk − xi) 2 同理，在点 Pi处的n阶离散导数为 y (n) i = ∑i+n k=i−n (xk − xi)(y (n−1) k −y (n−1) i ) ∑i+n k=i−n (xk − xi) 2 Pi ψd Pi= (xi , yi ,zi) 1 ⩽ i ⩽ n si 上述各阶离散导数是对连续可微函数在点 处的各阶导数的逼近。对于三维空间中的离散 曲线，为了计算其离散导数，首先要对离散曲线 进行参数化。可采用累积弦长参数化方法[24]来确 定离散点的参数值。在该方法中，离散曲线 在各 点 ( ) 处的累积弦长参数 可确定为 si = ∑i−1 k=1 |pk+1 − pk | ∑n−1 k=1 |pk+1 − pk | , 2 ⩽ i ⩽ n 并且 s1 = 0 ，从而离散曲线 ψd为 ψd = {Pi} = {(x(si), y(si),z(si))} si ∈ {s1,s2,· · ·,sn} xi = x(si) yi=y(si) zi=z(si) si ψd Pi {Pk |i−n1 ⩽ k ⩽ i+n1} n1 Pi Pi (x ′ (si), y ′ (si),z ′ (si)) 式中： ；坐标函数 、 和 就是弦长 的离散函数。因此，平面上的离 散导数的计算方法可以转化为用于计算三维空间 中的离散曲线上各点的导数。从而，采用上述方 法，对三维空间中的离散曲线 上的点 ，其邻域 为 ( 为点 一侧选取的邻近点 个数)，则点 处的离散导数 分别为 x ′ (si) = ∑i+n1 j=i−n1 (sj − si)(xj − xi) ∑i+n1 j=i−n1 (sj − si) 2 y ′ (si) = ∑i+n1 j=i−n1 (sj − si)(yj −yi) ∑i+n1 j=i−n1 (sj − si) 2 z ′ (si) = ∑i+n1 j=i−n1 (sj − si)(zj −zi) ∑i+n1 j=i−n1 (sj − si) 2 Pi Pi 上述离散导数也称为点 的一阶离散导数。 类似的可计算点 的二阶、三阶离散导数。从而 ·196· 智 能 系 统 学 报 第 14 卷
第1期 慕生鹏，等：三维离散曲线曲率挠率的微中心差分算法 ·197· 利用式(1)、式(2)即可求得三维空间离散曲线的 式中k2为邻域半径。 曲率、挠率。 同理，其三阶前向差商定义为 1.6微中心差分法 单侧差分法2在进行离散曲线的曲率挠率计 g0=1f"0-f- △t 算时，会对给定的曲率挠率产生一定的偏移或者 三阶后向差商定义为 缩放，为了克服由单侧差分法引起的这种计算偏 差，结合文献[22]中的思想和导数定义，把其单侧 f"0=员 f"(t+iA)-f"() At 差分法修改为对目标离散点的双侧差商取平均的 从而其三阶差商为 方法来得到目标离散点的各阶导数，并参考文 献[25]中的性能度量算法思想和差分方式，将改 "0=2(g"0*"02k,A"+kA)-f"-kA》 进的方法命名为微中心差分法(microcentral differ- (5) 式中k3为邻域半径。 ence algorithm). 因此对于微中心差分法，为了得到离散曲线 首先对微中心差分法名称中的“微”和“中心” 上给定的目标离散点的曲率值，根据本文的计算 作简单的说明。“微”是指在局部范围的逐阶差商 方法，则需要知道目标离散点邻域内的离散点数 计算。比如，先用改进的离散导数定义计算一阶 N,=2(k1+k2)+1,即目标离散点两侧各需k1+k2)个 差商，如式(3)：再利用一阶差商值计算离散曲线 离散点；而对于目标离散点的挠率计算，则至少 的二阶差商，如式(4)：最后利用二阶差商值计算 需要知道目标离散点邻域内的离散点数N, 三阶差商，如式（⑤）。“中心”指在目标点邻域内， 进行差分计算时，以离散曲线的目标点为中心进 2k1+2+k3)+1,其中目标离散点两侧各(k1+k2+k) 个离散点，k1、k、k均为大于等于1的整数。根 行差分计算。另外，方法中涉及求目标点的各阶 差商，我们指定，每阶差商目标点单侧需要上一 据式(3)（⑤）得到其各阶差商后，则可根据曲线曲 阶差商值的个数，记为k(i=1,2,3),称为邻域半 率、挠率的微分定义，将上述各阶导数式(3)~ 径。根据上述说明，微中心差分算法的计算方法 (⑤)代入曲率、挠率的计算式(1)、式(2)，即可求得 具体如下。定义其一阶前向差商为 离散曲线上各点处的曲率、挠率。这里提出的微 1点f0-ft-i) 中心差分计算方法是对单侧差分方法的一个依据 △t 导数定义进行的一个扩展，理论上是符合连续函 数条件下可导的要求。 一阶后向差商为 总体来看，上述6种方法分为2种思路：1)一 1月f+i)-f@ f()= 对近邻点进行曲线拟合，它们之间的区别在于近 At 邻点数的确定以及拟合曲线模型的假设；2)对近 从而其一阶差商为 邻点进行离散求导，各方法之间的区别在于取点 f0=0+o)2台 1ft+i)-ft-i△ 的策略以及如何利用近邻点的信息得到各阶导函 △ 数的可靠逼近，从而得到不同的曲率和挠率的计 (3) 算方法。 式中k为邻域半径。 其二阶前向差商定义为 2实验比较 f0= 1f-fa-i) 本节通过实验对上述6种算法计算的曲率和 挠率的计算精度进行评估。为了便于比较，本文 二阶后向差商定义为 中实验方法的k1、k2、k均取值为1。算法评估指 1台f+i)-f@ 标参考文献[26]的评价指标，选取最大相对误差 0台 △t (max relative error,MRE)和均方根误差(root mean square error,.RMSE),即对于获得的离散曲线，在 从而其二阶差商定义为 每一个离散点P,处的偏差定义为。=9，-，其中 )-A@kA-F(-kA 9为实验方法的计算值，为对应的准确值。因 ) 此，对于给定的离散曲线，其最大相对误差为EE=
利用式 (1)、式 (2) 即可求得三维空间离散曲线的 曲率、挠率。 1.6 微中心差分法 单侧差分法[22]在进行离散曲线的曲率挠率计 算时，会对给定的曲率挠率产生一定的偏移或者 缩放，为了克服由单侧差分法引起的这种计算偏 差，结合文献[22]中的思想和导数定义，把其单侧 差分法修改为对目标离散点的双侧差商取平均的 方法来得到目标离散点的各阶导数，并参考文 献[25]中的性能度量算法思想和差分方式，将改 进的方法命名为微中心差分法 (microcentral difference algorithm)。 ki(i = 1,2,3) 首先对微中心差分法名称中的“微”和“中心” 作简单的说明。“微”是指在局部范围的逐阶差商 计算。比如，先用改进的离散导数定义计算一阶 差商，如式 (3)；再利用一阶差商值计算离散曲线 的二阶差商，如式 (4)；最后利用二阶差商值计算 三阶差商，如式 (5)。“中心”指在目标点邻域内， 进行差分计算时，以离散曲线的目标点为中心进 行差分计算。另外，方法中涉及求目标点的各阶 差商，我们指定，每阶差商目标点单侧需要上一 阶差商值的个数，记为 ，称为邻域半 径。根据上述说明，微中心差分算法的计算方法 具体如下。定义其一阶前向差商为 f ′ q (t) = 1 k1 ∑k1 i=1 f(t)− f(t−i∆t) ∆t 一阶后向差商为 f ′ h (t) = 1 k1 ∑k1 i=1 f(t+i∆t)− f(t) ∆t 从而其一阶差商为 f ′ (t) = 1 2 ( f ′ q (t)+ f ′ h (t) ) = 1 2k1 ∑k1 i=1 f(t+i∆t)− f(t−i∆t) ∆t (3) 式中 k1为邻域半径。 其二阶前向差商定义为 f ′′ q (t) = 1 k2 ∑k2 i=1 f ′ (t)− f ′ (t−i∆t) ∆t 二阶后向差商定义为 f ′′ h (t) = 1 k2 ∑k2 i=1 f ′ (t+i∆t)− f ′ (t) ∆t 从而其二阶差商定义为 f ′′(t)= 1 2 ( f ′′ q (t)+ f ′′ h (t) ) = 1 2k2∆t ( f ′ (t+k2∆t)− f ′ (t−k2∆t) ) (4) 式中 k2为邻域半径。 同理，其三阶前向差商定义为 f ′′′ q (t) = 1 k3 ∑k3 i=1 f ′′(t)− f ′′(t−i∆t) ∆t 三阶后向差商定义为 f ′′′ h (t) = 1 k3 ∑k3 i=1 f ′′(t+i∆t)− f ′′(t) ∆t 从而其三阶差商为 f ′′′(t)= 1 2 ( f ′′′ q (t)+ f ′′′ h (t) ) = 1 2k3∆t (f ′′(t+k2∆t)− f ′′(t−k2∆t)) (5) 式中 k3为邻域半径。 N1 = 2(k1 +k2)+1 (k1 +k2) (k1 +k2 +k3)+1 (k1 +k2 +k3) 因此对于微中心差分法，为了得到离散曲线 上给定的目标离散点的曲率值，根据本文的计算 方法，则需要知道目标离散点邻域内的离散点数 ，即目标离散点两侧各需 个 离散点；而对于目标离散点的挠率计算，则至少 需要知道目标离散点邻域内的离散点 数 N2= 2 ，其中目标离散点两侧各 个离散点，k1、k2、k3 均为大于等于 1 的整数。根 据式 (3)~(5) 得到其各阶差商后，则可根据曲线曲 率、挠率的微分定义，将上述各阶导数式 (3)~ (5) 代入曲率、挠率的计算式 (1)、式 (2)，即可求得 离散曲线上各点处的曲率、挠率。这里提出的微 中心差分计算方法是对单侧差分方法的一个依据 导数定义进行的一个扩展，理论上是符合连续函 数条件下可导的要求。 总体来看，上述 6 种方法分为 2 种思路：1) 一 对近邻点进行曲线拟合，它们之间的区别在于近 邻点数的确定以及拟合曲线模型的假设；2) 对近 邻点进行离散求导，各方法之间的区别在于取点 的策略以及如何利用近邻点的信息得到各阶导函 数的可靠逼近，从而得到不同的曲率和挠率的计 算方法。 2 实验比较 Pi e i φ=φi−φi φi φi EMRE = 本节通过实验对上述 6 种算法计算的曲率和 挠率的计算精度进行评估。为了便于比较，本文 中实验方法的 k1、k2、k3 均取值为 1。算法评估指 标参考文献[26]的评价指标，选取最大相对误差 (max relative error，MRE) 和均方根误差 (root mean square error，RMSE)，即对于获得的离散曲线，在 每一个离散点 处的偏差定义为 ，其中 为实验方法的计算值， 为对应的准确值。因 此，对于给定的离散曲线，其最大相对误差为 第 1 期 慕生鹏，等：三维离散曲线曲率挠率的微中心差分算法 ·197·
·198· 智能系统学报 第14卷 max(p,-)/pl)均方根误差为ERsE= 3)循环曲线： 最大相对误差体现算法的稳定性，均方根误差反 r(t)=(cos(t),sin(t),2sin(t),t∈[0,2π] 映算法的整体有效性。实验平台为个人笔记本电 4)Viviani曲线： r(t)=(cos2(t),cos(t)sin(t),sin(t)),tE[0,2] 脑，计算机配置是Inter(R)Core(TM)i5-7300HQ 5)Clelia曲线： CPU@2.50GHz处理器和8GB内存，程序运行环 r(t)=(cos(2t)cos(t),cos(2r)sin(t),sin(2r)),t[0,2] 境是MATLAB R2014b。 6)环面螺线： 2.1实验数据 r(0=(2+cos(t)cos(2),(2+cos(t)sin(2r),),t∈[0,2 实验用如图1所示的6条连续可微参数曲线 需要说明的是，解析曲线方程用于生成离散 进行分析。由于有曲线方程，可以用经典微分几何 曲线，并用于计算离散曲线上每个采样点的准确 的方法（式(1）、式(2)计算出曲线上任意一个点 的曲率值和挠率值。实验中对每条空间参数曲线 的准确曲率值和挠率值，有了准确值便于进行误 通过均匀采样获得相应的离散曲线，采样点数 差分析。实验分析用到的6条曲线参数方程如下： n在定义域内从50个点变化到500个点，用于测 1)三次挠曲线： 试算法的有效性。第1章介绍的5种常用的离散 r(t0=(2t2,f),te[-2,2] 曲线曲率和挠率计算方法以及本文提出的算法都 2)多变式内接弹簧线： 采用离散曲线上的点的坐标进行计算，都用不到 r(t)=(2cos(3t)+cos(2t),2sin(3t)+sin(21),t),tE [0,2n] 曲线方程。 10 14 1.2 10r 0.18 0 1.0 0.16 0.8 15 0.14 0.6 0.12 -10 4 0.4 0.10 0-4202 0.2 0.08 54-20 0.06 (a)三次挠曲线 (b)多变式内接弹簧线 2 1.0 0.7 0.8 0.6 0 0.5 8 0 1.0 0.4 -1 05 0.5 -0.5 0.2 1.0 0.2 -0.5 0.5 0.1 -10 0 0.50 (©)循环曲线 (d)Viviani曲线 ■0.7 0.6 0.5 1 0.20 0 0 0.4 -1 0.15 0.3 0.10 0.5 1.0 0.2 .0 0.1 -2 0.05 -1.0 0.50 (①环面螺线 (e)Clelia曲线 图1用于实验的6条空间曲线（线颜色与挠率值对应） Fig.1 Six 3D curves used in the experiment (The color is corresponding to torsion value) 2.2误差分析及时间效率比较 线拟合法(PCF)、投影法(Proj)、单侧差分法 通过上述三维空间6条离散曲线来定量地评 (OSD)、离散几何法(DG)这5种算法的曲率和挠 价本文提出的微中心差方法(MCD),并将其与第 率的计算结果进行比较。图2~7分别给出了这 1章介绍的B样条曲线拟合法(BSCF)、多项式曲 6条离散曲线随着采样密度增加，曲率和挠率的
max(|(φi−φi)/φi |) ERMSE= √ 1 n ∑n i=1 (φi−φi) 2 , 均方根误差为 。 最大相对误差体现算法的稳定性，均方根误差反 映算法的整体有效性。实验平台为个人笔记本电 脑，计算机配置是 Inter(R) Core(TM) i5-7300HQ CPU @2.50 GHz 处理器和 8 GB 内存，程序运行环 境是 MATLAB R2014b。 2.1 实验数据 实验用如图 1 所示的 6 条连续可微参数曲线 进行分析。由于有曲线方程，可以用经典微分几何 的方法 (式 (1)、式 (2)) 计算出曲线上任意一个点 的准确曲率值和挠率值，有了准确值便于进行误 差分析。实验分析用到的 6 条曲线参数方程如下： 1) 三次挠曲线： r(t) = (2t,t 2 ,t 3 ), t ∈ [−2,2] 2) 多变式内接弹簧线： r(t) = (2 cos(3t)+cos(2t),2 sin(3t)+sin(2t),t), t ∈ [0,2π] 3) 循环曲线： r(t) = (cos(t),sin(t),2 sin(t)), t ∈ [0,2π] 4) Viviani 曲线： r(t) = (cos2 (t), cos(t) sin(t),sin(t)), t ∈ [0,2π] 5) Clelia 曲线： r(t) = (cos(2t) cos(t), cos(2t) sin(t),sin(2t)), t ∈ [0,2π] 6) 环面螺线： r(t) = ((2+cos(t)) cos(2t),(2+cos(t)) sin(2t),t), t ∈ [0,2π] 需要说明的是，解析曲线方程用于生成离散 曲线，并用于计算离散曲线上每个采样点的准确 的曲率值和挠率值。实验中对每条空间参数曲线 通过均匀采样获得相应的离散曲线，采样点数 n 在定义域内从 50 个点变化到 500 个点，用于测 试算法的有效性。第 1 章介绍的 5 种常用的离散 曲线曲率和挠率计算方法以及本文提出的算法都 采用离散曲线上的点的坐标进行计算，都用不到 曲线方程。 2.2 误差分析及时间效率比较 通过上述三维空间 6 条离散曲线来定量地评 价本文提出的微中心差方法 (MCD)，并将其与第 1 章介绍的 B 样条曲线拟合法 (BSCF)、多项式曲 线拟合 法 (PCF)、投影 法 (Proj)、单侧差分 法 (OSD)、离散几何法 (DG) 这 5 种算法的曲率和挠 率的计算结果进行比较。图 2～7 分别给出了这 6 条离散曲线随着采样密度增加，曲率和挠率的 (a) 三次挠曲线 1.4 1.2 1.0 0.8 0.6 0.4 0.2 10 −10 −2 0 2 4 −4 z 0 4 2 0 y x (b) 多变式内接弹簧线 0.18 0.16 0.14 0.12 0.10 0.08 0.06 10 5 0 z y x 5 0 −5 −4 −2 0 2 4 (c) 循环曲线 1.0 0.8 0.6 0.4 0.2 0 z 2 1 0 1.0 1.0 0.5 0.5 0 −0.5 0 −0.5 −1.0 y x (d) Viviani 曲线 0.7 0.6 0.5 0.4 0.3 0.2 0.1 z 1 0 −1 0.5 0 −0.5 y 0 0.5 1.0 x (e) Clelia 曲线 0.7 0.6 0.5 0.4 0.3 0.2 0.1 z 1 −1 0 1 0 −1.0 −0.5 0 0.5 1.0 y x y x (f) 环面螺线 0.20 0.15 0.10 0.05 z 1 −1 0 4 4 2 2 0 −2 0 −2 −4 图 1 用于实验的 6 条空间曲线 (线颜色与挠率值对应) Fig. 1 Six 3D curves used in the experiment (The color is corresponding to torsion value) ·198· 智 能 系 统 学 报 第 14 卷