双曲线定位算法

2024-05-24

双曲线定位算法(精选7篇)

双曲线定位算法 篇1

0 引言

无线传感网络(wireless sensor network,WSN)中实现其应用的关键是如何能够准确定位节点的位置信息。无线传感网络中的节点定位包括锚节点和未知节点定位。其中,锚节点是指少量带有GPS定位装置的节点,由于加装GPS设备能量消耗高,因此无法广泛使用,未知节点需要通过锚节点来进行自身定位[1,2,3]。根据定位过程中测量节点间的距离,定位算法分为基于测距(range-based)的定位算法以及基于无距离(range-free)的定位算法[4]。基于测距的算法需要测量相邻节点之间的距离,利用节点之间的实际距离来计算未知节点的位置,这种算法大多对硬件要求比较高,不太适合消耗功率和成本低的应用领域,目前常用的测距技术主要有RSSI、TOA、TDOA、AOA[5,6];基于无距离的算法无需测量节点间的实际距离,而是通过节点之间的估计距离来计算节点位置,可以降低节点的硬件要求,但是节点之间定位误差会增加,常用算法主要有质心算法、DV-HOP算法[7]。其中,DV-HOP算法主要通过距离向量路由协议来获得各个锚节点的信息,但是伴随着锚节点与未知节点之间跳数距离的增加,所测算的距离也会逐渐积累,从而对接点的定位精度造成了一定的影响。

本文首先分析目前DV-HOP算法3个步骤中存在的问题,然后将加权质心算法和二维双曲线概念引入到算法中,通过加权质心算法计算校正值,将接受信号RSS作为参考标准,有效减少误差,同时针对距离估算问题采用改进的二维双曲线算法,在二维双曲线基础上引入权值概念,使得估算距离更加精确。仿真实验表明,本文算法在校正值定位误差、最大估算误差精度都有一定程度降低,相比于文献算法有明显提高。

1 基本知识

1.1 DV-Hop算法

DV-HOP算法的主要思想是将未知节点与锚节点之间的距离通过网络平均跳数和跳数的乘积来进行决定,然后通过似然估计得到节点的位置。其算法分为3个步骤:

(1)信息广播与校正值计算。在无线传感网络中,每一个锚节点都会将自己的位置信息通过数据分组的形式广播出去,设定初始化跳数为0,当未知节点记录到每一个锚节点的最小跳数时,与已经存在跳数进行比较,如果大于已存在的跳数,则忽略不计,然后将跳数值加1,继续转发其余未知节点,这样每一个未知节点就可以得到锚节点的最小跳数。

(2)校正值广播。每一个锚节点将自身校正值以广播的形式播放出去,没有定位的未知节点得到各个锚节点的校正值。

(3)未知节点定位。通过计算锚节点的校正值和自身到锚节点的跳数计算锚节点的距离。然后再通过锚节点的坐标信息和到锚节点的距离信息,计算出未知节点的坐标信息。

1.2 DV-HOP算法的不足

DV-HOP算法的3个阶段都存在一定误差,具体如下:

(1)校正值方面的误差。在无线传感网络中,所有未知节点全部使用跳数与校正值的乘积来表示距离,计算出来的估计距离与真实距离存在着很大的偏差。假设待定节点与锚节点之间的跳数为1时,采用DV-HOP算法来计算两点的距离,如图1所示。

假设节点N为未知节点,A,B,C三个节点都为锚节点,3个锚节点之间的相互之间的距离如图1所示。其中,AB之间的跳数为6,BC之间的跳数为5,AC之间的跳数为5。各个锚节点能够计算出自己的校正值。其中,节点A的校正值为:(30+110)/(6+5)=12.72,节点B的校正值为:(30+60)/(6+5)=8.18,节点C的校正值为:(60+110)/(5+5)=17。根据最近未知节点到锚节点获取校正值的原则,节点N从节点C获得校正值是17,因此节点N到3个锚节点的距离为:NA为3*17=51,NB为3*17=51,NC为17*2=34。从节点N到节点A的直线距离为42,节点N到节点B的距离为45,节点N到节点C的距离为50,因此误差分别是60-51=9,60-45=15,60-34=26,误差率分别为:9/60=15%,15/60=25%,26/60=43.3%。由于在真实的环节中,无线传感节点的分配更加复杂,通过这种定位方法,容易造成未知节点N的定位位置偏离真实位置。

(2)未知节点距离测算不准。在无线传感网络中,未知节点的位置使用校正值与锚节点之间的跳数乘积来进行代替,这种方式显然忽视了复杂情况下的节点定位,存在非常大的误差。以图1为例,NC之间通过上述方面测出的距离结果为34,而实际之间直线的距离为45,误差率为(45-34)/34=32.4%,可以发现,NC之间的实际距离是无法改变的,NC之间的测算距离决定着误差率,很显然NC之间的测算距离越接近实际距离,误差率越小,显然采用目前这种测算方法,无法尽可能减少误差。

(3)最大估算问题。目前,在DV-HOP算法中广泛采用估算法来衡量节点位置,主要是因为算法简单,容易实现,但由于未知节点需要3个以上的锚节点来确定位置,未知节点与锚节点存在误差,这样容易造成一定范围内的误差累计。

2 基于改进的加权质心和二维双曲线的DV-HOP定位算法

2.1 改进的加权质心算法

针对DV-Hop算法中未知节点定位误差导致校正值方面不准确的问题,本文将接受信号强度(RSS)作为一个权重值,用作针对未知节点的定位进行改进。以图2为例,在一个未知节点N周围存在着不同的锚节点{A,B,C,D,E,…..K},设定锚节点的坐标为A(xA,yA),B(xB,yB),C(xC,yC)…..K(xK,yK),因此未知节点N采用表1的方式来记录未知节点接收周围锚节点发出的RSSI信号的强度,未知节点与锚节点之间的路径距离以及锚节点坐标信息等。

根据信道传输模型的特点可以得出:信号传输距离越远,其RSSI的值越弱,因此锚节点对未知节点坐标求解的影响比较大,由于未知节点处于所在区域的位置不能随意更改,所以未知节点与锚节点的距离和RSSI值作为参考就显得非常重要。本文从以上两个因素考虑,将RSSI强度值的大小对比为具体的数值,然后通过RSS数值与路径距离的比值作为参考系数,通过与锚节点的坐标的乘积来确定未知节点的坐标。计算如下:

其中,RSSIN→R表示未知节点到锚节点R之间的信号强度,SN→R表示未知节点到锚节点之间的路径距离。

2.2 改进的二维双曲线定位算法

通过前述估算算法可以发现,由于DV-HOP算法通常采用似然估计法计算节点位置,这种通过3个方程组相减的方法对坐标信息具有一定程度的损失,为了避免这种误差,采用改进的双曲线定位算法,为了更好地提高节点的定位精度,引入一个误差值ε。

假设未知节点坐标为(x,y),锚节点的坐标为K(xK,yK),则未知节点与锚节点之间的距离为:

步骤1:在距离测量的过程中获得一个近似距离length′N→i

步骤2:综合公式(2)、(3)得到:

步骤3:展开公式(4)得到:

步骤4:设定S=x2+y2,T=x12+y12代入公式(5)中,得到公式(6)

步骤5:将式(6)按照矩阵进行展开,得到公式(7):

步骤6:利用误差ε最大值设定锚节点权值。从收到的锚节点中找出距离未知节点最远的锚节点,记作lengthmax,由于距离最远导致产生的误差可能比较大,设定权值为1,然后计算未知节点与锚节点之间的距离,按照不同的路径距离设定权值集。各个锚节点的权值计算按照公式(8)计算:

步骤7:将公式(8)中各个锚节点的权值按照对角矩形形式表示为:

步骤8:通过矩阵来控制锚节点与未知节点之间不同距离对于节点定位的影响。构造目标函数F:

步骤9:使用最小二乘对公式(10)进行求解,得到K=(A-TA′)-1(A-TB′),其中,

2.3 算法流程

本文算法如图3所示。

3 算法仿真与分析

为进一步验证算法的有效性和实用性,本文在Windows平台上采用Matlab2010进行仿真实验。仿真环境选择100*100的区域,选择500个节点,其中100个为锚节点,400个为未知节点,所有节点都进行随机分布。

3.1 校正值定位误差

选取本文算法与文献[8]算法在校正值定位误差上的比较,进行100次仿真实验,如图4所示

通过图4可以发现,本文算法的误差范围在8%左右,而文献[10]算法误差近20%,在校正值计算上由于采用了改进的加权质心算法,使得误差范围缩小。

3.2 最大估算定位误差

选取本文算法与文献[8]算法对最大估算定位误差进行比较,如图5所示。

通过图5可以发现,本文算法误差范围在12%左右,而文献[8]算法误差近21%,在最大估算定位选择上由于采用了改进的双曲线算法,使得最大估算误差范围缩小。

3.3 改进后的整体算法的结果

图4、图5分别从校正值误差定位和最大估算定位对两种算法进行了比较。本文算法与文献[9]、文献[10]算法从整体上进行比较,如图6所示。

4 结语

无线传感网络中未知节点的定位一直是WSN中研究的重点,特别是未知节点定位的好坏直接影响到无线传感网络效率的高低。本文针对DV-HOP算法中校正值的误差,未知节点距离测算不准以及最大估算误差进行分析,针对校正值方面存在的问题采用了改进的加权质心算法,在最大估算算法上采用了改进二维双曲线的方法,改进后的算法在定位精度上有了明显提高。

参考文献

[1]AKYILDIZ I,SU W,SANKARASUBRAMANIAM Y,et,al.A survey on snesor network[J].IEEE Communicaitons Magazine,2002,40(8):102-144.

[2]NICLESCU D,AMERIC N L.Communication paradigms for sensor networks[J].IEEE Communications Magazine,2005,43(3):116-122.

[3]BOUKERCHE A,OLIVERIRA H A,NAKAMURA E F.Localization systens for wireless sensor network[J].IEEE Wirless Communications,2007,14(6):6-12.

[4]谭志,张卉.无线传感网络RSSI定位算法的研究与改进[J].北京邮电大学学报,2013,36(3):88-91.

[5]刘运杰,金明录,崔承毅.基于RSSI的无线传感器网络修正加权质心定位算法[J].传感技术学报,2010,23(5):717-721.

[6]何艳丽.无线传感网络质心定位算法研究[J].计算机仿真,2011,28(5):163-166.

[7]NICULESCU D,NATH B.DV-based positioning in adhoc network[J].Telecommunication Systems,2003,22(1):267-280.

[8]夏少波,连丽君,王鲁娜.基于DV-Hop定位算法的改进[J].计算机应用,2014,35(5):1247-1250.

[9]李道全,刘月月,孙付龙.基于DV-HOP的WSN网络节点定位算法[J].计算机仿真,2014,31(4):303-306.

[10]金纯,叶诚,韩志斌.无线传感器网络中Dv-Hop定位算法的改进[J].计算机工程与设计,2013,34(2):18-25.

[11]苏兵,薛伟杰,王洪元.一种WSN非测距定位DV-Hop算法的误差改进方法[J].计算机测量与控制,2013,21(5):1357-1359.

[12]周天绮,姜凤茹.基于MPSO-DV-Hop的无线传感器节点定位[J].计算机工程与应用,2013,49(23):52-55.

双曲线异型拱柱定位测量技术 篇2

关键词:双曲线,异型拱柱,定位测量

为获取设计图纸上双曲线拱柱准确的三维坐标位置, 需运用软件 (犀牛软件) 对图纸进行深化, 绘制出双曲线拱柱三维模型, 然后运用Auto CAD对双曲线三维模型定位, 获取准确的三维坐标, 再运用全站仪极坐标放样方法, 对双曲线拱柱进行现场实地测量放线工作。

1 工程概况

呼和浩特汽车客运东枢纽站工程 (图1) , 规划面积83 773 m2。由综合站务楼和生产服务大楼组成, 工程结构采用钢筋混凝土框架-剪力墙结构体系, 总建筑面积46 471 m2。生产服务大楼建筑面积19 559 m2, 下设一层地下室, 地上12层, 建筑高度49.45m;综合站务楼建筑面积26 912 m2, 下设一层地下室, 建筑高度22.4 m。综合站务楼为地下一层, 地上两层, 屋顶为钢结构, 站房南北两侧共计有对称弧形拱柱36根。

2 异型拱柱设计分析

(1) 根据设计图纸, 异型拱柱是在一个椭圆上沿着两个对称切线放线横切下来的两个相交椭圆形, 且拱柱本身为六边形, 不仅有一个本身的弧度还有一个整体倾斜的弧度, 综合站务楼又分地下室和非地下室, 两部分标高不同, 故拱柱投影在地下室和非地下室部分的投影也不同。设计图纸中给出的拱柱位置又是投影位置, 须通过计算和画图得出拱柱的实际位置。测量放线工作难度大、工作量大。

(2) 本工程拱柱是双曲线形状, 地下部分与地上部分为钢筋混凝土结构, 并与上面钢结构对接, 要求测量精度高, 拱柱模板采用GRC新型模板。无论模板还是钢筋均须测量对其定位, 钢筋是直的, 拱柱是曲线形状, 所以钢筋的走向须进行实际定位, GRC模板是预制模板, 分块进行拼装, 在现场测量须对每一块GRC模板进行定位, 避免拼装错误。

(3) 异型拱柱不仅是工程施工中的重点, 更是施工中的一大难点, 因为异型拱柱钢筋混凝土部分最后要与钢结构对接, 所以要求精度达到钢结构的测量施工精度, 且拱柱裸露在外面, 为达到设计效果, 要求所有拱柱须保持形状一致, 无偏差。

3 测量工艺流程

了解设计图纸画Auto CAD轴线图→运用犀牛软件和Auto CAD软件绘制三维模型得到定位图和拱柱坐标→控制网建立→拱柱实地测量定位放线。

4 异型拱柱三维建模及坐标计算方法

采用Auto CAD和犀牛软件配合建立三维模型。首先在Auto CAD里建立三维主要轴线, 导入犀牛软件中通过扫掠, 建立直观的三维模型, 然后把扫掠之后的三维模型导入Auto CAD里, 得到异型拱柱在平面上的具体位置, 在Auto CAD里把拱柱具体位置、尺寸标注出来, 最后在坐标图中把拱柱的具体位置坐标标注出来。

(1) 根据设计图纸上的轴线位置, 在Auto CAD里建立一个三维的轴线图, 画出拱柱中心线, 再画一个拱柱横切面的六边形并放在拱柱中心线上进行角度旋转, 所有的角度根据拱柱剖面图上给图的角度和距离计算得出 (图2) 。

(2) 把Auto CAD里画好的三维图导入犀牛软件中, 进行扫掠和填充得出拱柱三维模型 (图3) , 把拱柱的实际位置测设在拱柱的承台上, 在承台标高处把拱柱横切, 得到拱柱在空间中的实际位置。

(3) 把在犀牛软件中横切出来的拱柱实际位置导入Auto CAD中, 在Auto CAD中定位出拱柱的空间位置并标注出来。。

(4) 通过以上准备, 得出拱柱的定位图, 也知道了尺寸、距离, 根据角度和弧度和已知的控制点坐标, 计算拱柱的定位坐标。

5 施测过程

5.1 准备工作

施测人员通过对总平面图和设计说明的学习, 了解工程总体布局、工程特点、周围环境、建筑物的位置及坐标, 了解现场测量坐标与建筑物的关系, 水准点的位置和高程以及首层正负零的绝对标高。认真学习建筑施工图, 及时校对建筑物的平面、立面、剖面的尺寸、构造。熟悉拱柱三维建模图, 了解每块GRC模板的实际形状、拼接位置。施工图是整个工程放线的依据, 在熟悉图纸时, 着重掌握轴线的尺寸、层高。对比基础、楼层、建筑、结构几者之间轴线的尺寸, 查看其相关轴线及标高是否吻合, 有无矛盾存在。

5.2 测量仪器的选用

测量中所用的仪器和钢尺等器具 (表1) , 根据有关规定, 送至有校验资质的检测厂家进行校验, 检验合格后方可投入使用。

5.3 现场施测过程

(1) 先根据测绘院给的已知点在整个现场布控控制网。根据现场的实际情况在施工现场选择了4个点并按照规范埋设控制点。从已知坐标的2个控制点通过坐标反算得出两点的方位角, 并采用闭合导线的方法布设控制网, 得出各个控制点的坐标和方位角。

(2) 根据现场实际情况选择控制点进行拱柱放样并作为测站点。

(3) 拱柱的测量定位采用全站仪坐标法。本工程采用PENTAXR-322型全站仪, 仪器精度在误差允许范围之内。全站仪可通过坐标反算出要测设的点的距离和角度。所以在测站点架设仪器, 后视另一个点即可根据计算得出的拱柱坐标进行测设。在测设过程中, 为保证异型拱柱位置的准确性, 须进行多次后视, 取得限差最小值, 然后进行具体测设。在测设点达到精度要求以后进行做点, 再进行复核, 保证测设点的精度达到要求, 之后再进行下一个测设点的放样。

6 现场实际施测中的精度控制

在施工中, 从布设控制网开始即严格要求精度。布网时严格遵循先整体、后局部、高精度控制低精度的原则, 控制点要选在拘束度不大、安全、易保护的位置, 通视条件良好, 分布均匀。根据测绘院给的控制点在场区内引测4个控制点, 要求埋深1.5m, 用混凝土浇筑并以钢柱做标记, 测定高程作为工程定位放线依据。精度限差要求见表2。

在具体施工放线中, 一般的轴线限差在±5mm, 测角中误差在±12"。但由于设计要求限差必须小于±5mm且不大于±3mm, 所以在具体测设中, 先架设仪器须强制对中, 在对后视时在允许限差范围内多测几次后视, 得出限差最小值, 之后再进行拱柱的放线工作。放线时来控制角度偏差在±2"内进行测距。为保证拱柱位置偏差不超过±3mm须对测设点进行多次测设, 在±3mm以内做点, 做点之后再复测一遍, 复测不超过限差再进行下一个点的测设。

在测设点放样完以后, 为保证点不被破坏, 立即进行下一步的工序, 也为保证测设出来的拱柱形状的正确性, 利用钢筋做一个与拱柱平面形状一样的模型, 用来检验所测设的拱柱是否正确。

这样经过反复测设和复核之后的点, 精度大大提高, 保证了拱柱位置的准确性, 为下一步的施工提供最大限度的保证。

7 结束语

通过运用犀牛软件和Auto CAD软件的结合, 加上使用全站仪进行测量定位, 工作效率提高很多, 预计需要3个月的完成的工作量在1个月内便基本完成, 节省了很大一部分人力、财力、物力。且测量定位之后的异型拱柱也全部验收合格, 并且陆续进行下面的工序, 加快了整个工程的施工进度。

参考文献

[1]张正禄.工程测量学[M].武汉:武汉大学出版社, 2009.

曲线驱动的网格变形算法与应用 篇3

关键词:鞋楦,定制CAD,拉普拉斯网格变形,特征线

0 引言

鞋楦是制作鞋子的模具,其设计的好坏直接决定鞋子的合体性与穿着的舒适性。传统的鞋楦设计需要由设计师对物理样楦进行反复修改,设计周期长,对操作人员个人能力的依赖性强。随着三维激光扫描仪的推广普及和逆向工程技术的进步,快速获得标准样楦的数字模型成为可能,这样我们就可以基于数字样楦进行鞋楦的再设计。为此必须提出相应的楦型设计方法并开发配套的鞋楦定制CAD系统[1,2]。

文献[1]基于距离映射构建楦型编辑系统,并提出了一种鞋楦局部修改方法,但鞋楦变形不够自然。文献[2]从脚型数据出发,在抽取脚型参数的基础上,通过经验公式获得对应的鞋楦参数,从而构建个性化的鞋楦模型,但该方法难于通过界面交互直接对楦型进行编辑修改。

本文基于笔者自主研制开发的三维楦型/脚型激光扫描系统,以曲线驱动三角网格曲面发生变形为主体框架,构建楦型快速定制CAD系统。系统通过编辑特征线、截面轮廓曲线和侧影轮廓线等实现鞋楦再设计,不需点选许多约束变形点,操作方便,实时直观;直接通过交互编辑B-Spline曲线来修改轮廓外形,不需要计算源曲线与目标曲线间的对应点;可以利用成熟的B-Spline曲线编辑技术,在曲线形状调整满意后,再实施变形物体的自由变形,减少计算量,提高设计效率和显示速度。相对于自由变形(free form deformation,FFD)技术,笔者采用的基于微分坐标的变形方法,不但可以保留细节,还可以交互指定任意变形区域,控制楦型设计的影响区域,使设计更加灵活。

1 离散拉普拉斯网格变形

针对三维网格模型的变形与形状修改,Olga等[3]提出了基于离散拉普拉斯坐标的形状编辑算法,该算法将三维网格顶点的拉普拉斯坐标定义为顶点与其一阶邻域顶点集的质心之间的偏差矢量,在编辑过程中努力保持多边形的拉普拉斯坐标,从而保持原多边形的视觉特征。该方法具有许多优点,如保留细节变形、线性求解,可以很容易地添加各种线性约束条件,变形时尽可能保持原多边形的视觉特征,使得变形满足用户的各种需要等。

1.1离散拉普拉斯算子

给定一个具有n个顶点的三角网格模型M=(V,E,F),V为顶点集,E为边集,F为三角面集合。对于网格中任意一点,其相对坐标(δ坐标)可以线性表示为其绝对坐标与其直接邻域内顶点质心坐标的差[4,5],即

式中,wijvj点的权值;δi为点vi的拉普拉斯坐标(δ坐标);L为网格M的拉普拉斯算子。

为了便于表示和计算,网格模型的绝对坐标和相对坐标之间的转换采用矩阵形式表示。设网格模型M的邻接矩阵为A,I为单位矩阵,D为对角矩阵,其元素Dii =di,divi 1-ring邻域内的顶点个数。网格M的拉普拉斯算子矩阵L可以表示为

L=I-D-1A (2)

网格M上各个顶点的拉普拉斯坐标矩阵Δ

Δ=LV

1.2基于δ坐标的网格变形

我们的目标在于给定一初始网格模型,将其中一些顶点移动到目标位置,在尽可能保持其各顶点相对坐标和初始网格模型视觉特征的基础上,求解其变形后各顶点的绝对坐标。因为网格的拉普拉斯算子矩阵L的秩小于n,是奇异矩阵,也就是说δ坐标具有平移变换的不变性,矩阵L是不可逆的,不能仅由L唯一确定网格模型各顶点的绝对坐标。要唯一确定网格模型各顶点的绝对坐标,必须求解一个满秩的线性方程系统,为此必须增加相应的约束条件。

Q为初始网格模型,变形后的网格为Q′。Q中的部分顶点vjCb作为约束条件,其变形后的绝对坐标已知,vj=uj,jC,C={1,2,…,m}。这些点我们称之为“锚点”,这时的线性方程组扩展为

其中,V′为n×3阶矩阵,表示网格模型变形后各个顶点的绝对坐标;Gm×n阶矩阵,其中每行只有一个非零元素,对应于每个约束点,其值为

gkj={ωj=ikC01km,1jn

式中,ω为约束点的权值。

Fm×3阶矩阵,每列代表一个坐标分量,这样以此即可求出3个坐标分量Fk:

Fk=ω uj

约束条件的添加使系统变为超静定,但扩展后的拉普拉斯矩阵L˜是满秩的,因而可获得最小二乘意义上的唯一解:

V˜=argminV(LV-Δ2+jCω2vj-uj2) (4)

V=(L˜ΤL˜)-1L˜ΤB

矩阵L˜TL˜是稀疏且正定的,可以通过Cholesky分解快速求解该线性系统,并改善系统求解的稳定性。

2 曲线的获取

对楦型的修改主要是通过对鞋楦特征线(如鞋楦底板线、统口线等)和一些截面轮廓线(如跖趾围线、后跟弧线等)的编辑实现的,下面主要针对三角网格模型,说明其上的截面轮廓曲线和特征线的生成方法。

2.1截面轮廓线

首先,我们用点、面、半边和边的数据结构来表示三角网格模型,并建立其拓扑邻接关系。在进行三角网格模型拓扑重构后,即可利用模型中各个三角网格之间的拓扑信息计算平面与三角网格模型的交线。对于一个平面Γ,首先计算第一个与该平面相交的三角片T,得到交点坐标,然后根据局部邻接信息找到相邻的三角片,并求出交点,依次追踪,直到回到T,并得到一条有向封闭的轮廓环;重复上述过程,直到所有轮廓环计算完毕,并最终得到该层完整的截面轮廓[6]。

设已知平面Γ的方程为

A(x-x0)+B(y-y0)+C(z-z0)=0

其中,平面的法矢量为(A,B,C),且经过点(x0,y0,z0)。建立三维点数组Intersection来存储交点信息,其中,交点的ID记为交点所在边的ID,则求平面Γ与三角网格模型交点的具体步骤如下。

第一步,在完成模型拓扑重建工作后,遍历模型的实体边表,按序取出边表中的边,逐一判断,求出第一条与Γ相交的边,具体过程为:

(1)取第r条边,获取当前边的两端点Vs和Ve。

(2)分别将Vs和Ve的坐标值代入平面函数F(x0,y0,z0)=A(Vx-x0)+B(Vy-y0)+C(Vz-z0),可得对应的函数值分别为value_svalue_e,VxVyVz为点的坐标分量。

(3)对value_svalue_e进行判断:若value_s=0且value_e≠0,说明平面与边的交点即是端点Vs,将Vs加入数组Intersection;若value_e=0且value_s≠0,说明平面与边的交点即是端点Ve,将Ve加入数组Intersection;若value_e=0且value_s=0,说明整条边都位于平面上,此时认为边的中点为交点,将其加入数组Intersection;若value_evalue_s异号,说明该边与平面有交点,求交,并将求得的交点加入数组Intersection;若value_evalue_s同号,说明该边与平面没有交点,则从边表读入下一条边,进行求交判断。

(4)当边表遍历结束还未找到与平面相交的边时,这说明平面与模型不相交。

第二步,根据截线的封闭性,可记录第一条与平面相交的边ID为FirstEdge,将其作为求交循环结束的标志;取该边的两条半边中的任意一条作为求交循环的初始半边He_init。具体过程如下:

(1)由初始半边开始,将其作为当前半边。

(2)取当前半边的上一条半边He_p与下一条半边He_n进行判断:若He_p所在的边与平面有交点,根据平面最多与当前三角面片的两条边相交的原则,无需再判断He_n所在边是否与平面相交,只需直接求出He_p所在的边与平面的交点,并将交点加入Intersection即可;若He_p所在的边与平面没有交点,则取He_n所在边与平面求交,将求得交点加入Intersection

(3)判断交点所在边是否为初始边FirstEdge,若不是,则令求得交点的半边的伙伴半边He_a为当前半边,转第(2)步继续求交;若是,则求交结束。

以上求交方法的优点在于:①利用拓扑关系,使截面得到的交点集合是有序的,无须再进行排序,可直接获得首尾相连的有向封闭截线;②在三角面片与平面求交时,对某个三角面片只需计算一个边的交点,由面的邻接关系,可继承邻接面片的一个交点。

在得到截面轮廓上的有序交点集合后,采用三次非均匀B-Spline曲线对其拟合,便可利用成熟的B-Spline曲线编辑技术对截面曲线进行修改了。

2.2楦型特征线

这里的楦型特征线,主要是指楦底板线和鞋楦统口线,它们的形状和尺寸决定楦底和楦口的形状和尺寸,因而需要首先把这两条特征线提取出来,通过对它们的交互修改实现楦型的再设计。特征线的提取主要分为两步,即特征点的提取与特征线生成。

2.2.1 特征点提取

特征线的生成以提取出模型的特征顶点为前提,这里我们采用主曲率极值法来判断特征顶点[6,7]。首先计算点的法矢,对三角网格曲面上一顶点P,由其邻接三角形法矢的加权叠加计算该顶点的法矢;然后由P点及其邻接顶点在一局部坐标系内拟合出一个二次曲面S(u,v)=(u,v,h(u,v)),其中h(u,v)=au2+buv+cv2,如图1所示。曲面在P点的法曲率k中的极小值k1和极大值k2称为主曲率。k1、k2对应的法截线在P点的切线分别为m1、m2,这两条切线的方向称为主方向,两切线总是互相垂直的。

在计算出各顶点的主曲率之后,可通过对主曲率的极值进行判断得出特征点集,具体方法如图2所示。设m1及其反向延长线与三角形的交点分别为AB,当Pi点主曲率k1的绝对值大于AB两点在切线m1的法曲率k的绝对值时,则Pi点就是在m1方向上的曲率极值点,标示为特征点。AB点的法曲率可由VjVj+1在曲线m1的k1线性组合求得。该方法同样适用于k2,可判断Pi在曲线m2是否为特征点。经过上述判断得出的特征点还需经过进一步过滤,得出最后可应用于生成特征线的特征点集。

图3所示为对鞋楦三角网格模型进行特征点提取的示意,网格顶点数初始为18 064,最后提取的特征点数为3824。

(a)鞋楦三维模型 (b)提取的特征点

2.2.2 特征线生成

在得出鞋楦上的特征点集后,我们根据需要将其分离为统口线点集和楦底板线点集两部分。由于特征点集的分布类似于带噪声的测量点云,我们首先采用移动最小二乘法(moving least-squares,MLS)对特征点集进行细化和光顺;然后对细化后的特征点集进行排序和稀疏处理,使之有序化,便于采用B-Spline构造特征曲线。

设特征点集R经过MLS细化后所得点集为R′={Pi=(xi,yi)|i=1,2,…,N},从中随机选出一点Pj,根据需要,指定一个整数K,将其作为点PjK-邻近顶点个数,该值的大小决定了特征点集稀疏的程度。如果点集R′足够小,点Pj及其邻域内的点近似在一条直线上,搜寻PjK-邻域中距离Pj最远的点作为下一个搜寻点。点集排序与稀疏的算法流程如下:

(1)从点集中随机选出一点Pj作为初始点,计算其K-邻域Al,计算Al内各点到Pj的距离,存入数组dj,并按照从大到小的顺序对dj排序。把Pj存入排序后的点列数组NewP

(2)设Al中距离Pj最远的点为Pj+1,向量F=Pj+1-Pj表示点PjPj+1的方向(Pj+1、Pj分别为点Pj+1和Pj的坐标)。把Pj+1存入排序后的点列数组NewP

(3)计算点Pj+1的K-邻域NP,同样计算NP内各点到Pj+1的距离,并按照从大到小的顺序排序。

(4)从已排序的NP中循环找出一点Pl(坐标为Pl),计算向量FN=Pl-Pj+1与向量F的内积e=FN·F,判断e值的大小,如果e>0,则将Pl存入排序后的点列数组NewP,以Pl作为新的Pj+1转入(3)进行迭代;如果e≤0,则继续在NP中循环。如果点Pl位于初始点PjK-邻域A中,说明点集构成一个封闭曲线,迭代终止。

对点云进行有序化稀疏后,即可利用比较成熟的有序点曲线重构技术对其进行曲线重构,本文采用3次非均匀B样条作为重构曲线类型。首先采用积累弦长参数化方法对有序化的点列进行参数化,从而确定曲线的节点矢量,由节点矢量反算三次B样条插值曲线的控制顶点,求解控制顶点时采用抛物线边界条件。在计算出曲线的控制顶点后,即可用德布尔算法对其进行正算,得出曲线上的点,其具体计算细节请参照文献[8]。图4所示为由楦底板特征点生成的特征曲线。

3 曲线驱动的网格变形

3.1曲线与网格模型间对应点的计算

要使曲线的变形能够驱动鞋楦模型做相应的变形修改,必须使模型与曲线相关联,使模型上的一些点与曲线上的采样点建立一一的映射关系,这样在曲线修改变形的时候就能带动网格模型做相应的变形。对于截面轮廓曲线,在构建曲线过程中,已自然建立了模型与曲线的关联,而对于特征线,由于对特征点进行了稀疏,曲线不完全经过特征点集,必须重新计算曲线上采样点与模型上顶点之间的映射关系,使特征线与楦型相关联。

首先以模型的平均边长为参考值,在特征线上按照等参变化取一系列的采样点c(ti),在模型上找其对应的最近点Vi,如果c(ti)与Vi的距离小于某一指定阈值,认为Vi为模型上的曲线关联点之一,记录Vi点的ID和其在曲线上的对应参数值ti

这里有两个问题需要注意,一是对于曲线上的两个甚至三个相邻采样点,它们在模型上的最近点有可能是同一点,这时只记录距离最小的那一组对应点;二是对于曲线上的采样点c(ti)而言,其在模型上的最近点为Vi,但反过来,Vi在曲线上的投影点并不一定是c(ti),须进一步计算其在曲线上的f oot point点,修正其在曲线上的对应位置坐标。

如图5所示,设曲线上采样点c(ti)在模型上的最近点为Vi,通过牛顿迭代计算点Vi在曲线上的对应参数值t*。要计算t*,可通过最小化点Vi到曲线的距离获得,即使下式最小

g(t)=‖c(t)-Vi‖2 (5)

其中,c(t)为曲线上任一点的坐标;Vi为点Vi的坐标,采用牛顿迭代计算式(5)使g(t)最小以获得t*。设t(i)为当前迭代值,下一次迭代值为[9]

迭代的初始值取Vi在曲线上的初始对应点c(ti)的参数值ti

建立了曲线与模型间的关联之后,这些与曲线关联的模型顶点将作为位移约束加入到求解网格变形的线性系统中。

3.2变形区域指定

在对楦型进行编辑时,根据设计意图,必须在鞋楦模型上指定一个变形区域,也称为兴趣区域(region of interest,ROI)。一方面这是鞋楦编辑的需要,即需要设计者按照设计意图制定所编辑区域的影响范围,以保证鞋楦的其他部分不被修改;另一方面,由于网格模型编辑的计算量与顶点数有关,指定区域后,可以减小计算量。

变形区域的指定按照鞋楦编辑的功能模块指定,重点在于如何构建其邻接矩阵和处理其边界。鞋楦编辑中,楦底编辑、后跟弧线编辑、顶弧编辑以及底弧编辑等的变形区域相对比较规范,其指定方法在后文叙述。要对楦型进行局部编辑,需要交互指定一个局部变形区域,指定过程如下:

(1)建立每个三角面包含的三个半边结构,记录每个半边所对应的三角形的ID[10],如图6所示。

(2)设网格曲面的法矢方向指向模型外部,通过人机交互,用鼠标按逆时针方向点选模型上的一系列顶点,其中最后一点需和第一个点重合。

(3)由相邻的两点生成dijkstra近似测地线,得到一个由网格三角形的边组成的封闭环,构成所选区域的边界,标记边界点为已访问。

(4)从边界点中选择两相邻顶点ViVi+1,则ViVi+1构成一条半边,选择该半边所在三角形的第三个顶点P作为选择区域搜寻的种子点,如图7所示。

(5)从种子点P开始,依托模型的拓扑邻接信息,按照广度优先算法搜索边界内的所有顶点及其对应的三角形,直到边界内所有点都标记为已访问。

在获得选择区域内的顶点后,就可以根据其拓扑邻接信息建立该选择区域的邻接矩阵,并进一步建立其拉普拉斯算子矩阵L。其中需要特别处理的是所选区域的边界点,因为一方面在构建邻接矩阵时,边界上的点需要特殊处理,只计那些处于选择区域内的邻接点;另一方面,需要对这些边界点作出标记,以便在做楦型修改时将其作为边界约束加入到求解网格变形的线性系统中。

3.3曲线约束的网格编辑

我们的目标在于通过交互修改曲线来驱动鞋楦模型做相应的变形,达到按设计意图进行楦型定制的目的。在获得指定区域的拉普拉斯算子矩阵L后,我们将把边界约束和曲线修改引起的位移约束加入求解系统,构建式(3)。

所选区域边界上的顶点将作为“锚点”约束加入求解系统,即这些边界顶点在变形过程中保持位置不变,以保证变形区与非变形区域的平滑过渡。可以通过调整边界点的权值来调整过渡区的变化。

对于模型中与曲线直接关联的顶点,将作为位移约束加入求解系统。设模型中顶点v与曲线上的点c(ti)相关联,曲线修改后,点c(ti)变为q(坐标为q),则认为v点修改后的位置为

v′=v+(q-c(ti))

在边界约束和位移约束条件下求解式(4),即可求的修改后的鞋楦模型的顶点位置。

4 鞋楦CAD

基于前述方法,我们采用VC++6.0和OpenGL技术在微机平台上实现了基于样楦的鞋楦定制CAD系统。系统主要分为文件管理模块、楦型重构模块、楦型编辑设计模块、数据库管理模块和视图管理模块等。其中楦型编辑设计模块是本系统的主体模块,分为楦底板编辑、后跟编辑、底弧编辑、顶弧编辑、楦头编辑和局部区域编辑以及楦型放样等子模块。其他功能模块不是本文的主要内容,这里不做详细介绍。

各个模块的功能和操作面板基本相似,以楦型底板编辑模块为例,首先针对楦底编辑的要求,选定一变形区域作为设计变更的影响区域,如图8所示。

而对于特征线的生成,可以采用灵活的方法,可以按照前述方法由系统自动生成,也可以通过手工交互输入,甚至可以通过外部导入。曲线的修改方法如下:主要通过交互拖动曲线上的型值点不断调整曲线的形状,在曲线形状修改满意后,再预览整个鞋楦的变形情况。这样把曲线修改和网格变形分步处理,可以减少计算量,提高设计效率和实时显示速度。如果对预览模型的编辑结果不满意,还可以继续修改曲线,直到获得满意的设计效果,最后确认设计结果。图9所示为通过调整楦底板线进行楦型再设计的结果。

相较于自由变形技术,本文方法可以指定变形区域,不会发生变形干涉现象,在变形的同时可保持细节特征。

5 结论

(1)采样曲线驱动三维网格变形方法实现楦型的编辑修改,不需点选许多约束点,易操作,交互性好,符合设计师设计习惯,容易表示设计意图。

(2)可以任意指定设计影响区域,设计灵活;采样先修改曲线,满意后再驱动三维模型作相应的更新,计算效率高,实时显示。

(3)采用离散拉普拉斯网格变形方法,容易增加各种线性约束条件,使得变形满足用户的各种需要。

参考文献

[1]Leng J,Du R.A CAD Approach for Designing Cus-tomized Shoe Last[J].Computer-aided Design&Applications,2006,3(1/4):377-384.

[2]徐从富,刘勇,蒋云良,等.个性化鞋楦CAD系统的设计与实现[J].计算机辅助设计与图形学学报,2004,16(10):1437-1441.

[3]Olga S,Yaron L,Daniel C O,et al.Laplacian Surface Editing[C]//Proc.of Eurographics Symposium on Geometry Processing.Nice,France,2004:179-188.

[4]Marc A.Differential Coordinates for Local Mesh Morphing and Deformation[J].The Visual Comput-er,2003,19(2):105-114.

[5]Nealen A,Sorkine O,Alexa M,et al.A Sketch-based Interface for Detail-preserving Mesh Editing[J].ACM Transactions on Graphics,2005,24(3):1142-1147.

[6]上官宁.基于三维网格模型的鞋楦CAD关键技术研究[D].泉州:华侨大学,2008.

[7]刘胜兰,周儒荣,张丽艳.三角网格模型的特征线提取[J].计算机辅助设计与图形学学报,2003,15(4):445-448.

[8]施法中.计算机辅助几何设与非均匀有理B样条[M].北京:高等教育出版社,2001.

[9]Hu S M,Wallner J.A Second Order Algorithm for Orthogonal Projection onto Curves and Surfaces[J].Computer Aided Geometric Design,2005,22(3):251-260.

基于曲线插值算法的巷道自动建模 篇4

巷道系统是矿山三维虚拟场景的重要组成部分, 是构建数字矿山的基础。井下巷道纵横交错、错综复杂, 工作地点及资源条件不断变化, 如何立体、直观、准确地表现并反映井下巷道及其空间关系, 是矿山测绘科技工作者的重要研究课题, 也是煤矿安全、高效、合理开发的重要保障[1]。

由于每个矿井的实际巷道不同, 对矿井巷道用三维软件建模得出的模型不宜推广。本文基于对煤矿安全培训系统的研究, 编写出根据具体巷道数据信息生成不同矿井巷道的具体模拟系统。

1 巷道的网格模型

在三维图形程序设计中, 网格模型占有非常重要的地位, Direct3D支持的模型文件格式为*.X文件格式, *.X文件模型通常是由大量的三角形拼接成的网格 (mesh) 模型。

1.1 解析.X文件

三维模型建模软件很多, 而三维模型可以转换成Direct3D能够识别的数据文件格式 (后缀名为.X) , 这样在Diect3D中载入和渲染文件模型就很方便了。

对坐标系中的物体进行点采样 (Point Sample) , 这些采样点按一定顺序连接成为一系列的小平面 (三角形或共面的四边形, 五边形等) , 这些小平面称为图元 (Primitive) , 3D引擎会处理每一个图元, 称为一个独立的渲染单位。这样取样后的物体看起来像是由许许多多的三角形, 四边形或五边形组成的, 就像网一样, 称为一个网格 (Mesh) 。 纪录这些顶点数据和连线情况到一个文件中, 3D引擎读取这些数据, 依次渲染每一个图元, 就能在显示屏幕上再现物体。取样的点越多, 再现的物体也会越逼真, 要处理的数据量也越大。在D3D中, 纪录这些顶点数据和连线情况的文件称为X文件 (X File) 。它是以X作为文件名后缀的。.X文件主要是用来存储网格数据的, 它通常存储了三维模型的顶点坐标、颜色、法向量、纹理坐标以及动画帧等信息。

.X文件是由模板驱动 (Template-Driven) 的, 模板定义了如何存储一个数据对象, 这使得这种文件格式具有结构自由, 内容丰富, 易应用, 可移植性高等优点[2]。

1.2 编写.X文件自动生成模型

3ds Max或Maya等三维建模软件所制作的网格模型可导出为Direct3D所支持的.X文件模型, 利用三维建模工具制作的三维模型通常比较复杂, 即多边形数量很多, 而多边形数量越多, 图形渲染速度越慢, 所以在自动建模时, 在不明显影响视觉效果的前提下尽量减少多边形的数量。

一般来说巷道体的断面表现为梯形或直壁拱形, 其规律和断面点信息比较容易找出, 并且在尽量减少多边形数量而不影响视觉效果的前提下, 可以找出规律并按照矿井实际数据, 编写.X文件, 制作相对应的巷道网络结构, 提高了加载、渲染速度[3]。

2 巷道的三维建模及算法

巷道断面是巷道几何建模中的重要参数, 主要有拱形、矩形、梯形、斜梯形等形态, 本文是利用拱形断面来制作巷道模型的, 根据巷道底面中线上点的位置信息来求取巷道断面上点的位置信息, 从而应用Direct3D扩展库提供的函数来编程实现巷道自动建模。

2.1 拱形巷道断面

Direct3D中的坐标系是左手坐标系, 图1是初始拱形断面在Direct3D坐标系中的位置。

一个断面上设置12个点, 并给出这12个点的初始信息, 如图1所示, 其中undefined。

利用D3DX扩展函数库d3dx9.lib提供的D3DXMatrixRotationY ( ) 、D3DXMatrixScaling ( ) 、D3DXMatrixTranslation ( ) 等函数, 旋转、缩放、平移断面初始信息, 得到所需断面信息, 编程实现断面点信息的连接, 写好.X文件, 即可完成巷道的自动生成。

由图2可知, 在中线上加载两个断面可以生成一个直巷道, 但实际情况, 一个煤矿有许多不同的巷道。对于含分岔口的巷道, 作加载半断面处理, 即对于图1所示的拱形断面, 分半保存点信息, 点0至点5组成半断面, 点6至点11组成半断面。在分岔口处加载半断面即可构成巷道体, 如图3所示, 在A2-A, A-A1, B-B1, B-B2, C2-C, C-C1, D2-D, D-D1, E2-E, E-E1线上加载半断面, 所加载的半断面要由初始半断面经过平移、缩放、旋转得到, 加载上断面后, 写好.X文件, 巷道模型即可生成。为了点B1、B2处的拐弯过度的圆滑一些, 可以在A-A1与D-D1之间, C-C1与D-D2之间内插出几个半断面。内插的半断面可以根据改进的Bezier插值算法求出。

2.2 曲线插值算法

曲线Q (u) 有三个控制点P0, P1, P2, 控制点的坐标发生改变, 曲线的形状随之改变。

曲线以参数u (0≤u≤1) 来定义, 曲线上的每一个点都是通过一个三次多项式对每一个控制点进行比例调节来确定的, 曲线Q (u) 是由B样条曲线改进所得, undefined

插值曲线的矩阵形式:

undefinedundefined

利用曲线插值算法对巷道内插断面, 如图4所示, 内插断面的点U1, B, U2必须在巷道的中线上, 来保证巷道无缝连接[4]。

2.3 采用插值算法的效果

直接利用断面的加载生成的巷道模型的拐弯处很尖, 不符合实际。为了使模型更逼真一些, 采用曲线插值算法, 内插一些断面, 使巷道模型拐弯处过度的更圆滑一些。比较图5和图6可以看出采用插值算法后的效果。

3 结束语

Direct3D是一套非常优秀的高性能三维图形程序可编程接口, 它对游戏、三维图形程序开发提供了全方位的支持。本文所采用的编程写.X文件自动建模方法可以快速、逼真地构造出巷道模型系统。文中提到的曲线插值算法实现起来比较简单, 运行速率也比较高。实验证明, 这种巷道自动建模方法, 大大简化了模型的设计和程序的编制。

摘要:通过分析Direct3D所能识别的数据文件 (.X文件) 格式, 编程实现组织数据得到.X文件, 快速自动生成三维巷道模型。并针对巷道拐弯处的圆滑问题, 采用了改进的B样条曲线插值算法。实验表明, 基于曲线插值算法的巷道自动建模, 比建模软件建模易推广且运行速率高。

关键词:三维巷道,自动建模,插值算法

参考文献

[1]王德才, 杨关胜, 孙玉萍.精通DirectX3D图形与动画程序设计[M].人民邮电出版社, 2007.

[2]汪云甲, 伏永明.矿井巷道三维自动建模方法研究[J].武汉大学学报, 2006, 31 (12) .

[3][英]Alan Watt.3D计算机图形学[M].包宏, 译.北京:机械工业出版社, 2005.

双曲线定位算法 篇5

关键词:Bezier曲线,控制顶点,矢量

曲线是计算机图形学的重要内容, 它是描述物体的外型, 建立所画对象的数学模型的有力工具[1]。Bezier曲线因具有良好的几何性质, 能简洁、完美地描述和表达自由曲线曲面, 是最常用的曲线之一。

在利用Bezier曲线造型时, 如果阶次太高, 固然能表示复杂的造型, 但易造成计算复杂度增加, 控制顶点增多, 形状不易控制。而二次Bezier曲线表示能力有限, 又是平面曲线, 所以最常用的就是三次Bezier曲线。本文讨论了三次Bezier曲线传统递推算法实现, 又提出了一种新的递推算法实现, 使其绘制速度和绘制效果大大提高。

1 Bezier曲线的描述

定义:给定空间n+1个点的位置矢量Pi (i=0, 1, 2, …, n) , 则Bezier参数曲线上各点坐标的插值公式是:

其中, Pi构成该Bezier曲线的特征多边形, 是n次Bernstein基函数:

2 Bezier曲线的分割定理

Bezier曲线P (t) , 经中间某点分割而得到的两条曲线分别为Q和R, 这两段曲线同样是Bezier曲线, 它们的控制定点由de Casteljau算法产生, 即:P=Q+R

3 三次Bezier曲线的生成算法

曲线Q曲线P和曲线P (t) 相比, 显然其控制顶点离曲线更近, 那么, 再次把曲线Q分割, 分割后得到的曲线的控制定点和曲线将距离更近。依此类推, 如果不停的分割, 将会得到一系列分割后的Bezier曲线和其控制顶点, 当控制定点和曲线的距离小于指定的很小正数ε, 就可以用这些控制顶点的连线所组成的折线近似代替曲线。为方便计算, 设分割的比例系数t=1/2。另外, 还应注意一种特殊情况, 如图2所示, 当控制顶点所组成的控制多边形中相邻两边的夹角θ为锐角时, 控制多边形与Bezier曲线外形有较大的差距, 此时控制顶点的连线不宜作为Bezier曲线的近似曲线, 因而需直接进行下一步的曲线分割。

因此有以下两种情况:

(1) 当 时, 据下述3.1或3.2的方法进行处理。

(2) 当 时, 直接进行曲线分割分割。

3.1 传统递推算法实现。

如图3所示, d1'和d2'表示控制顶点到曲线的距离, 当max (d1', d2') 小于指定的很小正数ε, 即可以用这些控制顶点的连线所组成的折线近似代替曲线, 而不再对曲线进行再次分割。直接计算控制顶点P1P2到Bézier曲线P (t) 的距离很麻烦, 但是由Bézier曲线的凸包性可知,

d1和d2表示点Pi到线段P0P3的距离, 计算点到直线的距离要相对容易, 故在一定的误差范围内, 可用右端代替左端。

当d1和d2中最大值小于ε时, 就可以用控制顶点的连线代替曲线段。

3.2 新递推算法实现。

4所示, 图中阴影部分面积为A1、A2, 当max (A1、A2) 小于指定的很小正数ε, 即可以用这些控制顶点的连线所组成的折线近似代替曲线, 而不再对曲线进行再次分割。

阴影部分面积A1, A2如下:

如图5所示, 当满足条件曲线分割时A1, i、A2, i分别为

结束语

综上所述, 绘制三次Bezier曲线如果用传统递推算法实现需要两次用到矢量叉乘计算点到直线的距离, 而如果采用新的递推算法实现, 除了第一次计算A1、A2的面积需用矢量叉乘之外, A1, i+1只需用上一次的计算结果A1, i即可求得, 计算量降低了一半, 算法效率提高显著。经实验验证, 绘制出的曲线效果良好。

参考文献

[1]杜华.Bezier曲线算法研究与实现[J].考试周刊, 2010, 12:159-160.

[2]倪明田, 吴良芝.计算机图形学) [M].北京:北京大学出版社, 1999:218-221.

[3]孙家广.计算机图形学 (第三版) [M].北京:清华大学出版社, 1998:306-307.

[4]罗振东, 廖光裕.计算机图形学[M].上海:复旦大学出版社, 1993:114-128.

水库库容曲线算法的研究及应用 篇6

关键词:库容曲线,最小二乘法,均值插入法,三次样条采样法,二次曲线拟合法

1 概述

在水库的设计施工过程中,库容和水位是一组非常重要的参数[1],它们的准确与否将会影响到水库的正常运转。在水库及大坝设计、施工过程中,需要根据计算或测量得到一组离散的数据,推导并建立水位和库容、水位和流量的数学函数模型,并以此为依据计算水库的库容。传统的水库库容曲线确定方法主要有以下几种:一根据在不同高程下水库的不同面积求出对应的库容,并得出相应的库容曲线;二是在应用时是通过查表或直线内插;三是通过Execl进行曲线拟合。但是上述方法要么不方便,要么不准确。实现一种方便高效的水位库容计算方法正是该文的研究重点。

2 算法原理及实现

通过观察多个水库的水位和库容、水位和流量数据并结合实际工程施工的经验,我们发现水库的水位和库容、流量之间有着特定的函数关系。所以,我们通过以下三种曲线拟合方式来进行水库库容曲线的拟合。

2.1 均值插入法

在一个给定的区间内进行等距离插值,即在给定的区间[a,b]上等距离的插入n个节点值xi=x0+i*h(i=0,1,2,…,n-2,n-1;h为插值的步长值),x和y之间的函数关系为yi=f(xi)。应用拉格朗日插值定理[2],就可以得出在插值点t处的函数近似值z=f(t)。

曲线拟合效果如图1所示。

2.2 三次样条采样法

数学上将具有一定光滑性的分段多项式称为样条函数。具体地说,给定区间[a,b]的一个分划∇。

如果函数s(x)满足:

1)在每个小区间[ai,bi](i=0,1,…,n)上s(x)是k次多项式;

2)s(x)在[a,b]上具有k-1阶连续导数。

则称s(x)为关于分划∇的k次样条函数,其图形称为k次样条曲线。

给定平面上(p0,p1,p2,p3,…,pn-1)n个固定点,要求样条曲线依次通过这n个点,就可以确定p1,pn-1之间的样条函数公式。根据这个公式和一定的步长stepVal,就可以通过插值得到P1,P2之间的一系列点q1,q2,...,qn,然后将这些点依次用线段连接,就可以得到一条拟合出来的样条曲线。步长值越小,拟合出的点就越多,得到的样条曲线越逼近真实曲线的值。但是,数据量就越多,计算量也就更大。曲线拟合效果如图2所示

2.3 二次曲线拟合法

根据经验,水库的库容和水位的关系是一个二次函数曲线关系。因此,我们可以应用最小二乘法[3]理论来拟合水库库容和水位之间的函数关系,即拟合函数值与实际值之差的平方和最小的函数为所求的函数。应用最小二乘法,采用2次曲线拟合的方式拟合已有的采样点,根据这些采样点所描述的2次曲线的变化,求解二次曲线方程y=ax2+bx+c。此方程就是我们所要求解的水位库容曲线经验方程。

在曲线拟合过程中,根据给定的采样值(xi,yi)(i=0,1,…m),求出x和y的函数关系y=f(x)=ax2+bx+c,其中a,b,c为方程系数。在实际工程中,采样数据是有限的,而且会有一定的误差。所以在实际求解的过程中f(x)=ax2+bx+c不可能全部通过每一个实际的点(xi,yi)(i=0,1,…m),只能是在工程允许的误差范围内求得最优解,即:在给定点(xi,yi)上求满足δi=s(xi)-yi(i=0,1,…,m)的平方和取最小值时,yi对应的最优值。得到一组满足方程的(x,y)解之后,再根据水位库容关系方程求出其系数a,b,c即可。

曲线拟合基本步骤如下:

1)给出采样点值,选择合适的曲线类型

该文给出的一组实测值为:(4,7),(5,9),(6,11),(7,22)。根据相关项目的工程经验,在中小型水库中水位和库容的关系满足一个线性的二次方程关系,即:f(x)=ax2+bx+c。

2)按最小二乘法原理解方程,求出方程系数

根据上述最小二乘法原理可知,只需要求解满足误差最小的方程即可得出我们所求的库容曲线方程,即是求解

当方程组有解的时候,必须满足如下的极值必要条件:

根据内积定义,可以得出一个关于方程系数的一个线性方程组,用矩阵表示为:

根据最小二乘法原理,只有求出方程组的唯一解,就可以得出整个方程的系数,记为:ak=ak*,k=0,1,...,n,根据这组系数就可以得到最小二乘拟合曲线。该文解方程得到a,b,c的值分别为:a=2.25,b=20.05,c=51.65。

3)根据第二步得出的方程系数a,b,c确定库容曲线方程。该文拟合得出方程为:f(x)=2.25x2-20.05x+51.65与使用Execl拟合得到的方程完全一致。

曲线拟合效果如图3所示。

对比上述三种算法的模拟效果,均值插入法计算最简单快捷,但是得到的曲线平滑性最差;三次样条采样法和二次曲线拟合法得到的曲线都较为平滑,但是二次曲线拟合法得到的结果更加贴近实际值,而且拟合的结果和Execl拟合结果一致。

3 总结

本研究的三种曲线拟合算法,对计算机硬件的要求不高,易于掌握,与水情监控系统的集成也十分方便,可广泛应用于中小型水库工程的设计实践中。

但是考虑如下两点原因:一是库容实测曲线本来就是对面积实测曲线的近似积分求得;二是曲线拟合得到的水库水位库容曲线只是在理论上非常逼近实测库容曲线,而不是实测的库容实测曲线。因此,在实际利用时必须验证其误差,以及讨论工程对误差的允许量。

参考文献

[1]陈昌杰,姜和平.应用Excel软件拟合库容曲线[J].小水电,2002(6):19-20.

[2]周长发.C#数值计算算法编程[M].电子工业出版社,2007.

双曲线定位算法 篇7

1 二维Bezier曲线的参数表示

在曲线运动过程中求交点, 无疑需要借助有关算法在动作脚本中以代码实现。而要用代码求两Bezier曲线的交点, 首先需要将曲线表示成函数。但对Bezier曲线而言, 要用普通函数来表示是困难的, 幸可用参数方程。当给定平面上n+1个点di (i=0, 1……n) , 则n次Bezier曲线可表示为其中Cni是组合数, 即, 而P (t) 是参数为t时曲线的位置, 是一个二维矢量。容易看出方程右边实际上是一个n次多项式, 可整理成, , 这里系数ai是一个二维矢量。

鉴于AS3.0中无法重载运算符, 因而无法定义矢量运算类, 需将矢量方程化成标量方程, 形如:其中ai、bi为标量系数, 可从给定控制点di求得1。对于三次Bezier曲线, 两参数方程的系数, 其中 (xi, yi) 为控制点di的横、纵坐标。MB为三次Bezier曲线的4×4矩阵。下面给出根据控制点求参数方程系数的AS3.0代码如下:

2 二维Bezier曲线求交算法及其优化

2.1 扫描法 (步进法)

2.1.1 算法与代码

扫描法是最先想到的求交方法, 基本思路是用一条与横轴垂直的直线 (扫描线) 从两曲线相交区域的左端向右端扫描, 扫描线每到一个位置x, 以便计算两曲线纵坐标之差的绝对值dy=|y1-y2| (如图1) , 当dy小于给定的精度ε (一个小正数) , 便将点 (x, (y1+y2) /2) 近似为两曲线的一个交点。代码如下:

为了突出算法重点, 代码中没有给出两曲线相交区域左、右端的计算过程。这在ActionScript3.0中不难解决, 可用MovieClip类的getBounds () 、getRect () 函数获得两曲线的包围盒, 再用Rectangle类的intersection () 方法求出两包围盒的相交区域 (不妨记为rect) 。于是相交区域的左、右端即为:

2.1.2 存在问题

该算法虽然简单明了, 但实现时问题不少。

(1) 根据x值求y值的getYByX (x) 函数的实现问题

由Bezier曲线的参数方程可知, 要根据x值求y值, 需分两步:

(1) 先从⑴式中由给定的x值求出参数t值;

(2) 再从⑵式中根据t值求出y值。

对于三次Bezier曲线, 第 (1) 步要解关于t的一元三次方程, 并不简单, 所幸的是ActionScript3.0的BezierSegment类提供了根据方程系数求三次方根的函数getCubicRoota (a, b, c, d) 。对于更高阶的Bezier曲线, 就没有现存方法可用, 只有用代数方程的数值解法求得近似解。

(2) 无值、多值问题

上面第 (1) 步从方程⑴中求t值时, 理论上有n个根。剔除[0, 1]以外的值后, t值的个数在0与n之间, 事先无法确定。相应地第 (2) 步得到的y值的个数与t值个数相同, 也在0与n之间。因此, 在计算dy=|y1-y2|时, 要根据y值的个数作出相应的处理。如果y1或y2的个数有一个为0, 则回到循环开始进行下一循环;否则, 需对每个y1值与每个y2值两两求差取绝对值, 再行判断, 需用到双循环。

(3) 步长、精确度与交点个数

当步长step较大时, 运行速度较快, 但会发生交点漏失, 即明明相交的地方, 因步长较大被跨过, 而未被检测到, 使求得的交点个数小于实际交点数。

当步长step很小时, 虽可能避免交点漏失的情况, 但会带来两个问题, 一是运行速度减慢, 二会使交点个数虚增。即在一个交点附近, 多次被判断为相交, 使求得的交点个数大于实际交点数。而若减小精确度, 又可能导致交点漏失。

2.1.3 改进思路

对于问题⑴, 曾经想到将“对x值扫描”改为“对参数t扫描”, 这样根据t值求相应的y值时, 无须解n次方程, 只需将t值代入方程⑵即可求得, 且也避免了多值问题。但这仅仅得到了一条曲线的纵坐标值。由于t值相同时, 两条曲线的横坐标值未必相同, 因此另一条曲线的纵坐标值不能简单地用相同的参数t值来求, 而要根据一条曲线的x值求出另一条曲线的参数值 (为了区别第一条曲线的参数, 用s表示) , 再由s求出另一条曲线的y值。

这里, 根据一条曲线的x值求另一条曲线的参数值s还得解n次方程, 还会有多个y值。改进效果不明显, 放弃。

对于问题⑵, 既对参数t扫描仍然无法彻底解决多值问题, 又尝试了对纵坐标y扫描。一般来说, 给定x求y值时有多个值时, 反过来给定y值求x值, 多值问题会有所改善, 但仍然无法避免。比如, 当曲线为类似圆时, 无论前者还是后者, 都会有两个值。因此, 考虑到兼容性, 针对多值问题的这个双循环是免不了的了。

对于问题⑶, 想到的是在扫描过程中, 改变步长。即在不太可能相交时, 步长大一点, 在可能相交时, 步长小一点。前者可以提高速度, 后者可以避免交点漏失。想法虽好, 但问题是, 如何判断某一时刻是不太可能相交, 还是将要相交, 不太可能相交是, 步长取多大, 快要相交时步长又该取多大。一个简单的办法时, 根据dy (在横向扫描时) 的大小来决定步长。在代码2中计算出dy后, 添加以下代码:

这一办法, 虽然对求交速度和准确性有所改善, 但改善不大。为此又尝试了用两曲线的斜率、二阶导数、曲率作为决定下一步步长的依据, 但都无功而返。最后采取的办法是, 求扫描线处两曲线切线的交点, 比较切线交点与扫描线的有向距离来判断两曲线的相交趋势, 进而决定步长的大小。如图2中曲线A的切线LA与曲线B的切线LB相交于点P, 点P与扫描线X的距离d即是判断相交趋势并决定步长大小的依据。图2中是即将相交的情形, 可取步长为d乘以一个小于1的正数。如果P在扫描线的左边, 则说明暂时不会相交, 步长可适当取大一些。当然在具体实现中, 还有很多特殊情况需要考虑, 在此, 就不一一细说了。

这样一来, 虽然增加了求曲线切线及其切线交点的运算, 但可大大减少扫描的次数, 实验证明, 对提高求交速度和交点数量的准确性还是有一定效果的。

2.2 两分 (折半) 查找法

2.2.1 算法与代码

两分查找法是一种十分常用的算法, 通常用于在一个排序数组中寻找某个元素。笔者将其用于曲线求交, 得归功于《数学手册》中用两分查找法解代数方程的启示。若设两曲线相交区域的左、右端为L、R (见1.1.1) , 则两曲线交点的求解过程是:

(1) 在区间[L, R]中插入中点M, 将原区间分成左、右两个子区间[L, M]和[M, R] (如图3) 。

(2) 分别计算x=L、M、R时两曲线y值之差 (不是绝对值) , 记为ldy、mdy、rdy, 则左子区间两端两曲线y值之差为ldy, mdy;右子区间两端两曲线y值之差为mdy, rdy。

(3) 若某子区间两端y值之差同号, 则直接无视, 舍去这一子区间 (图3中[M, R]) 。若某子区间两端y值之差异号 (图3中[L, M]) , 在这一子区间中再插入中点M, 将此子区间分成两段, 重复 (2) , (3) , 直到子区间的长度小于给定精度ε, 或者|ldy|、|mdy|、|rdy|中有一个小于ε为止。若是|ldy|<ε, 则将L作为交点的横坐标, x=L时两y值的平均值作为交点的纵坐标;若是|mdy|或|rdy|<ε, 则将M或R作为交点的横坐标, 相应两y值的平均值作为交点的纵坐标。

2.2.2 存在问题

与前面的扫描法相比, 效率明显提高, 但缺陷仍然不可避免:

(1) 局限于只有一个交点的情形。若有多个交点, 仅依靠上述代码无法求全。如图4 (a) , 左、右两子区间, 两端y值之差均异号, 代码3中最后两个判断句均为真, 但实际执行时, 当执行了第一个判断, 取了左子区间后, 遇到continue语句, 就开始对左子区间继续进行两分查找, 而第二个判断则执行不到, 从而把右子区间给舍去了。因而, 只能求出左段包含的交点。而如果将continue语句去掉, 则进入下一循环时, L=R=M, 不满足循环条件, 结束循环, 于是一个交点也求不到。

(2) 在根据x值求y值时仍然面临着无值、多值问题。如图4 (b) , 当x=L时, 求曲线A的Y值时有两个, 求曲线B的Y值时, 无值, 从而无法求得ldy;当x=R时也有类似情况而无法求得rdy。

(3) 上下相切时 (如图4 (c) ) , ldy, mdy, rdy始终同号, 切点将无法求得。

2.2.3 改进方法

对于问题⑴, 改进办法是将两分查找法与分治法相结合。先将两曲线相交区域的x值域[L, R]分成若干区间, 对每一区间运用两分查找法。但随之而来的问题是对区间[L, R]如何划分。这里采用的方法是, 将相交区域的左、右端及两曲线所有位于相交区域内的控制点的横坐标值放入一个数组, 将数组升序排列。这样数组中每相邻两元素便构成一个区间, 再对每一区间运用两分查找法求出每一段中包含的交点。

对于问题 (2) , 类似于1.1.3中对相似问题的改进, 在横向分割存在无值、多值问题时, 改用纵向分割。即对相交区域的y值域进行分割, 这样就变成以x值求y值为根据y值求x值, 原来的问题在大部分情况下, 可以得到解决, 但仍有例外。

问题⑶在两分查找法中无法解决。

2.3 牛顿法

2.3.1 算法与代码

这也是受《数学手册》中牛顿法解代数方程的启示。本来的思路是:从曲线的一个端点出发画切线, 与x轴交于x1, 再画曲线在x1处的切线交x轴交于x2, ……直到切点与x轴的距离小于给定的精度ε为止。

该算法应用于曲线求交时, 步骤如下:

(1) 过两曲线的起点分别画切线LA、LB (在Bezier曲线中就是前两个控制点的连线) ;如图5所示。

(2) 求出两切线的交点 (记为P (x, y) ) 。

(3) 过P画垂线交曲线A、B于PA、PB。

(4) 计算PA、PB的距离dy, 如果dy<ε, 则将点 (x, (PA.y+PB.y) /2) 作为交点的近似。否则, 过PA作曲线A的切线, 过PB作曲线B的切线, 重复 (2) 、 (3) 、 (4) 。

这里事实上隐含了数学中常用的化曲为直思想, 将求曲线交点转化为求直线交点。AS3.0代码如下:

上述代码中根据曲线上一点求过该点切线的函数getTan (P) , 需先求出曲线在该点的斜率, 再根据给定点、给定的切线长度及斜率求出切线段另一端点的坐标, 即可得到切线。

2.3.2 存在问题

对于正常情况, 循环二、三步即可得到一个交点, 效率比两分查找法又有提高, 但2.2.2节中的前两个问题还是存在。

2.3.3 改进方法

要解决求多个交点的问题, 总体思路还是分而治之的分治法, 但分的方法与两分查找法不同。这里用的方法是, 选阶数较低曲线两端点的连线作为LA, 另一曲线的控制点两两相连所成直线为LB1, LB2……, 再分别对LA与LB1, LA与LB2, ……使用牛顿法求交 (即重复本算法描述中第 (2) 、 (3) 、 (4) 步。假定曲线ca、cb的次数为na、nb, 控制点数组为pa、pb, 且na

代码中函数by PP (p1, p2) 为根据两端点重新定义直线段。

无值、多值问题的解决。鉴于y值的个数取决于t值的个数, 因此把由x值求y值转变为由x值求t值, 得到t值组成的数组:tsx=crv.getTByX (P.x) 。为了避免过P点的垂线与曲线不相交, 从而tsx数组为空, 再根据该点的纵坐标求一下t值, 得到另一个t值组成的数组tsy=crv.getTByY (P.y) , 两个数组元素的个数在0-n之间, 既可能无值, 也可能多值。而只需要一个值, 到底取哪一个呢?约定一个原则:所取t值对应的曲线上点与给定点P的距离是所有t值中最小的。为此, 要做以下3步:

(1) 在数组tsx中找出与给定点P距离最近的t值tx及最小距离dx。

若数组tsx非空, 对tsx中的每一个t值, 求出相应的y值, 及与给定点P的纵坐标的距离, 取距离最小时的t值为tx最小距离记为dx。

若数组tsx为空, 用给定点P与曲线两端点的距离大小进行取舍, 若P与起点较近, 取t=0, dx为P与起点的距离;否则若P与终点较近, 则取t=1, dx为P与终点的距离。

(2) 在数组tsy中找出与给定点P距离最近的t值ty及最小距离dy。

若数组tsy非空, 对tsy中的每一个t值, 求出相应的x值, 及与给定点P的横坐标的距离, 取距离最小时的t值为ty最小距离记为dy。

若数组tsy为空, 用与tsx为空同样的方法取t=0或t=1。

(3) 比较dx与dy, 若dx≤dy, 取tx, 否则取ty。

据此, 可得到代码:

这样便从根本上解决了无值、多值问题, 对于la与lb的交点P, 总可以得到曲线ca、cb的唯一一个参数值ta、tb, 从而可以求出曲线上与参数对应的点pa, pb和斜率ka、kb, 也可求出pa、pb两点的距离, 过pa、pb的曲线的切线, 从而可以把代码4改写成

2.4 离散法

2.4.1 算法描述

这是王国瑾等给出的Bezier曲线求交方法, 其基本思想还是化曲为直, 将曲线求交转化成直线求交。但转化的方法与牛顿法不同, 这里采用的是分割化直的办法, 即将曲线分割成若干等次的Bezier曲线, 当分割得很小时, 用子曲线段两端点的连线代替子曲线段, 然后求一曲线的所有子曲线替代直线段与另一曲线的所有子曲线替代直线段的交点, 求得的交点即为原曲线交点的近似解。其堆栈算法描述如下:

(1) 初始化:

将两曲线的控制点数组、当前分割次数 (=0) 压入堆栈。

(2) 循环 (直到堆栈为空止)

上述算法已在王国瑾给出的算法上作了些许改进, 改进之处为包围盒相交性检查的位置。原算法在曲线分割后, 不管子曲线与另一曲线 (或其子曲线) 的包围盒是否相交, 不作判断, 直接压入堆栈, 在循环开始从堆栈弹出后再行检查。这样, 无疑增加了大量的堆栈压入、弹出操作, 浪费了机时;同时也会使堆栈的空间占用大量增加。据此, 我们将包围盒相交性检查放在曲线分割后, 压入堆栈之前, 对于包围盒不相交的子曲线段直接无视, 不再压入堆栈, 仅把包围盒相交的子曲线段压入堆栈。虽然在代码上增加了不少, 但实际检查次数并未增加, 且在减少堆栈空间占用的同时, 减少了大量无效的堆栈压入、弹出操作。实验证明, 速度明显加快。

2.4.2 算法实现

上述算法, 在实现时还有一些具体工作要做:

(1) 先验分割次数的计算

王国瑾给出了先验分割次数 (先验离散层数) 的计算公式:

其中:n为Bezier曲线的次数。

ε为给定的精度 (一个小正数) , 使分割后子曲线上任一点到子曲线段两端点连线的距离小于此值。

实际上, 为了简便, 可视曲线的弯曲程度直接给定一个先验分割次数。

(2) 曲线包围盒相交性判断

如果运用自定义的继承于MovieClip类的McCurve类, 可用MovieClip类的getBounds () , getRect () 方法获得曲线的包围盒。但鉴于本算法用到的仅是曲线的控制点, 而McCurve类包含了大量与本算法无关的属性, 建立众多该类的实例, 难免会增加空间占用。因此在该算法实现上, 我们没有用McCurve类。于是便不能直接应用上述获取包围盒的方法, 需另行定义获取方法, 当然这并不困难, 无非是求出所有控制点的纵、横坐标的最小值、最大值, 据此构造一个Rectangle类实例。进而可用Rectangle类的intersects方法判断两矩形是否相交。

(3) Bezier曲线的分割

这在不少资料上有介绍, 可在t1=0.5处执行de Casteljan递推算法。

Pi (i=0……n) 是n阶Bezier曲线的n+1控制点, t∈[0, 1]。所得的两条子Bezier曲线的控制点分别为:

P0i (i=0……n) , t∈[0, t1]和Pjn-j (j=0, ……n) , t∈[t1, 1]。

2.4.3 存在问题与改进办法

该算法求交速度快, 适用性好, 而且避免了getYByX () 或getXByY () 可能出现的无值、多值问题。但在两曲线相切或接近相切时, 也会出现交点个数虚增的情况。对此的改进办法是, 在算法情况1中, 求得一个非空交点后, 将堆栈中所有分割次数大于曲线次数-1的元素全部弹出, 不作处理, 这样既节约了时间, 又避免了交点个数的虚增。

另外一个需要注意的问题是, 在求连线交点时, 所得交点可能在线段之外, 这可能导致所求交点不在两曲线之上。改进办法是判断所求交点是否位于两连线段的包围盒相交区域内, 若是则返回该交点, 否则返回null。在AS3.0中可以借助Rectangle类的intersection方法求得两连线段包围盒的相交区域, 再用其contain Point (point) 判断点是否在相交区域内。

2.5 代数法 (解非线性方程组法)

若曲线由显函数表示, 如曲线A:y=f (x) , 曲线B:y=g (x) , 则求曲线A、B的交点只需解方程组即可。但Bezier曲线通常只能用参数方程表示, 现记曲线A的参数方程为:, 曲线B的参数方程为, s, t∈[0, 1]分别为曲线A、B的参数, na, nb分别为曲线A、B的次数, ai、bi、ai’、bi’为参数方程的系数。于是要求曲线A、B的交点, 则需解以下两元Max (na, nb) 次方程组:

当曲线的阶数≥2时, 传统的消元法, 矩阵法便无能为力, 因而常到此却步不前。偶而在《数学手册》中发现“牛顿法解非线性方程组”一节, 顿时感觉眼睛一亮, 真是踏破铁鞋无觅处, 得来全不费功夫。

2.5.1 算法描述

从一组近似解P0 (s0, t0) 开始用迭代公式:

直到据此求得的两曲线上的点PAn+1与PBn+1之间的距离小于给定的精度ε为止。

2.5.2 算法实现

作为迭代出发点的近似解的确定, 可以用牛顿法中的方法, 取次数较低曲线两端点的连线作为LA, 取另一曲线的控制点两两相连所成直线为LB1、LB2……, 再求LA与LB1、LB2……的交点, 根据每一个交点用1.3.3中的getTByP (crv, p) 求出两曲线的参数s0, t0, 作为方程组的近似解开始迭代。

各偏导数的计算。由于u, v均为多项式函数, 故求偏导数并不困难, 可用以下公式直接求得:

行列式的计算, 由于迭代过程中用到的都是两阶行列式, 因此只要对角相乘再相减即可。

实验表明, 该算法效率最高, 准确性最高, 但也会有交点数虚增的现象, 但造成虚增的原因不在于算法的核心部分, 而是由于分治法, 在不该分治时分治所致。

3 各种求交算法比较

对各种求交算法针对7组不同的曲线分别求交40次, 计算每一次求交的运行时间, 实验结果如下 (实验环境:Intel Pentium 4, 2.93GHz CPU, 736M内存, Microsoft Windows XP SP2, Adobe Flash Player10) :

可见, 诸方法中以解非线性方程组的方法效率最高, 也最稳定。平均求交时间仅0.654毫秒, 标准差仅0.236毫秒, 极差也只有0.75毫秒, 说明面对各种不同情况求交时间差异不会太大。而扫描法是效率最低的方法, 即便用了改变步长的方法, 其平均求交时间仍然高达12.8毫秒, 而且稳定性差, 求交时间的极差高达23.575毫秒, 标准差也达8.556毫秒。

参考文献

[1]《数学手册》编写组.数学手册[M].北京:高等教育出版社, 1979, 5.

[2]曹锋.基于Bezier曲线求交的曲线裁剪算法[J].计算机应用, 1998, 8:20-22.

[3]王国瑾, 等.计算机辅助几何设计[M].高等教育出版社、施普林格出版社, 2001.

上一篇:虚拟教学实验下一篇:文献结构