散乱数据点的三次多项式插值.doc

散乱数据点的三次多项式插值第10卷第j期1998年9月计算机辅助设计与图形学CAD&CGV0【.10,NO5Sep.,1998一散乱数据点的三次多项式插值/张彩明孙德汪嘉业(山东大学计算机科学系济南9.501001摘要对散乱毗驹I2分邢在每个;角邓上构造一个三班多项或曲面片,整体的C曲面由各三角形上的曲面片拼台而成.讨诧了整体C曲面需满足的条件组成的方程组的性质,并培出了求解方程蛆的方法.插值方法的多项式准确桌包括所有三次和小于三波的多项式关键词兰兰t芝,中图分类号TP391.22cADCUBICPoLYN0MIALINTERP0LAT10NoFSCATTEREDDATAP0INTSZHANGCaiMingSUNDe._FaWANGJia.Ye(DepartmoztComputerScie,ShandogUniversity,250100)AbstractThispaperpresentsamethodtOconstructcubicpolynomialinterpolationonthescattereddatapoints.Themethoddividesthegivenregionintotriangulation,andconstructsacubicpatchoneachtriangle.ThepropertieseftheequationsetcornposedoftheCconditionswhichthewholesurfacehastosatisfyisdiscussedandatechniqueforsolvingtheequationsetisgiven.Thepolynomialprecisionsetofthemethodincludesallthepolynomialsofdegreethreeorless.Keywordsscattereddatapoints,interpolation+polynomial,surfacepatch1引言构造插值曲面是数值计算,CG和CAD等领域中的一个重要问题.三角形阿格能较好地逼近不规则区域,所以常采用把给定区域三角化的方法构造插值曲面.多项式曲面由19961205收稿.19970324收到修改穑.奉文得到回国人员科研资助基金资助张彩明t1955年11月生,博士.副教授.主要研究方向为计算机图形学.计算机辅助设计,图象处理及模式识别.孙德法.硕士t主要研究方向为计算机圆形学和什算机辅助设计.汪嘉业,教授.主要研究方向为计算几何,计算机图形学和计算机辅助设计5期张彩明等:散乱数据点的三次多项式插值417于其构造简单,易于计算等特点而被广泛用于构造插值曲面.众所周知,在三角形网格上构造C连续的分片多项式曲面需要在各三角形上给出21个插值条件.文献1中方法需要的插值条件是三角形顶点处的函数值和一阶偏导数值,各三角形上构造的是4次多项式曲面片;文献2,3中方法除了文献Eli中的插值条件外,还需三角形边中点处的一阶法向导数值,各三角形上构造的是C连续的分片二次或三次多项式曲面片.由散乱分布的空间数据点构造插值曲面,所知道的只是三角形顶点处的函数值,虽然顶点处的偏导数及其它条件可由某些方法近似得到,但有时往往得不到理想的结果.用分片多项式曲面对散乱数据点插值已有许多研究成果,其中构造分片三次多项式曲面的关键是解一个线性方程组.由于方程组解的不唯一性而提出了多种解方程组的方法,例如文献E6,10,143,这些方法的缺点是构造的插值曲面的多项式精度较低.本文提出了一个构造具有较高精度的插值曲面的方法,该方法的多项式准确集包括所有三次或小于三次的多项式,对于用分片三次多项式构造插值曲面,这已是最高的多项式准确集.构造C连续的分片三次多项式曲面需解一个线性方程组,方程组的大小和给定插值点的个数成正比,当不是太大时,本文方法是可行的.在文献6,10,14中,各三角形上构造的是Beizer形式的三次曲面片,为计算方便,我们采用另一种形式的三次曲面片.文中首先讨论了相邻三角形上两曲面片在其公共边界上满足C连续的条件,然后讨论了由这些条件组成的方程组的性质,并给出了解方程组的方法,最后是一个例子.2相邻三角形上曲面片c连续的条件2.1三角形上三次曲面片FAx,y)的构造设丁.是三角形网格中具有顶点P一(却,Y)(p-i,J,),编号为t的三角形(图1)点P处给定的函数值为F(p=i,J,记点对边上的单位外法向为(M)口,P处的一阶偏导数为(F),(F)(p=i,J,设(,J,),(L),(LD)表示丁.的面积坐标(重心坐标)(,J)=(,J,).=-x.YJ+()+(z.一z.)y/(2S)(,J)=(厶(z,):一x,yk+(ykY.+.一xDy/(2S)(厶)一(厶(,)=-x.Y一+(一Y,)z+(,一0)y/(2S)(1)圉1其中S是丁.的面积.计算可知,丁.上对Fn(FD和()(p-i,)插值的三次曲面片(z.)是F,)一E(32(L)F+(z卸)(F)+(一)()(,J).I,J.+S(,J),(,J),(L)?(2)其中,(F)n(F)(p=i,)是待定参数.418计算机辅助设计与图形学1998芷2.2相邻三角形上曲面片C连续的条件_PI图2对图2中编号为1和2的两个三角形,设由式(2)定义的两个三次多项式曲面分别是Fl(,.).)和F2(z,.).).F1(z,.).)和Fz(,.).)在公共边的两个顶点处有相同的函数值和日一阶偏导数值,所以它们在公共边上是cc连续的,又因为F(,.).)和F:(z,)在公共边上的法向导数是二次多项式,所它们在公共边上c连续的条件是EaV(,)/a()I三一一aF(z,.).)/a()l一(3)其中(z,Y)一(,+)/Z,(十)/2),EaV/aN是F的法向导数,计算知(so,.).)/a()一(Ear/o(N)(F)+E/a(N)(F)/4+(3F+(以一墨)(F),+(挑一Y,)(F)Ea(L1)/O(N)/2+(Ea/a(N),(F)+E/a(N),(F)/4+(3F+(,一.72)(F)+(Y)(F)Ea(L1)/O(N)/2+AEa(L)/O(N)/由式(3)和对称性得EaV(so,)/a(1)+EaV2,y)/o(N2)l:=一(3F+(一)(F),+(一Y)(F)(Ea(L.),/a()+a(,J)/a(N)/Z+(3F+(一)(F)+(一)()(Ea(L)/a(N)+a(,J)/O(N)/2+SAEa(L)/O(N)/4+SAEa(Lz)/O(Nz),40整理得+=G,J,五,(4)其中G口,五,一4(3F,+(zk五)(F),+(一yi)(F)(Ea(L),/a()j+a(L2)/2)+(3F+(,一)(F)+(一)(F)(Ea(L)/a(N)+Ea(L)/o(Nz)/DJ_u-_-_J_u_一D=(墨)+(y).在每个三角形上都按式(2)定义了一张曲面片,且相邻三角形上的曲面片满足式(4),则由所有三角形上的曲面片拼台而成的整体曲面P,.).)是c连续的.由式(4)知,构造,(z,)需解一个方程组,该方程组记为E.称式(4)为P,)的一个约束条件,E中的未知量是待定参数.下一节将讨论方程组E所含约束条件的个数和待定参数的个数.5期张彩明等:散乱数据的三次多项式插值4193约束条件和待定参数的个数为简化描述引进两个符号:表示由个顶点构成的三角形网格;L表示AN的边界的顶点个数.定义1一个三角形网格称为的子网格,如果它是由删掉一些三角形形成的.定义2被称为具有A属性,如果只有一个三角形或有下列性质:(1)的边界上至少有_个不大于4的顶点P,在以P,为顶点的三角形的角中,任何两相邻角的角度之和不等于180.,如图4,1+2180.,2+34:180.;(2)删去P后形成的子三角形网格仍具有性质(1).图3和图4中的三角形网格具有A属性,图5中的三角形网格不具有A属性./图3图4图5从定义l可推得:引理1若,具有A属性,则至少有一个的子同格序列40=3,4,),使得4+可由A加上一顶点构成,且|+的三角形数目至多比A增加3.定理1设P(l,2,3,)是平面上个不重合顶点,由P(l,2,3,)构成的三角网格中三角形的个数为,为三角形网格边界的顶点个数,则M一2NL一2证明由欧拉公式得+一1,其中是中边的数目.由S=(3一,|)/2+,得M=2N-L一2.证毕.中相邻三角形公共边的条数是(3M+,|)/Zz,所以P(,)要满足的方程组E的约束条件的个数E一(3M4-L)/2一LN=3N一2L3(5)由式(2)知,每一三角形对应一个待定参数A,每个顶点P处有两个待定参数(F)和(F),因此,E中E个约束条件中待定参数的个数U=M+2N一4NL一2(6)E中待定参数和约束条件的个数之差,一EN4-L_-1(7)定理2若具有A属性,则有3一2,|一3个约束条件的方程组E是可解的.证明由引理1可知存在一子网序列(3,4,),我们用归纳法证明E.的系计算机辅助设计与图形学数矩阵的秩(简称E的秩)为3一2L一3(8)对于i一4(图6),若中有两个三角形,E的秩为1,若中有3个三角形棚由式(4)图6知,E4是A+A一GF,J,k,A+A.一GEl,k,f,力A.+A一GEl,显然E.的秩是3,所以对一4式(8)成立.设i=k一1时式(8)成立(5),即&一的秩是3(一1)一2厶一13,下证i=k时式(8)成立.根据引理1,可由j一加上一度不大于4的顶点得到,厶的三角形的个数等于一中三角形的个数加上(13).中的约束条件的个数比一中约束条件的个数多2一l,中待定参数的个数比一中待定参数的个数多+2,因为13,所以2一1+2.下证五的秩=五一的秩+2一1,只考虑一3的情况(一1和一2可类似证明).不失一般性,在图7中,设PP,尸属于j一,一加上顶点P.构成,)由一厅(,),F(,),Fz(,)和Fa(,)构成.由式(4)知,由一和下列方程组成A1+A.=GD,k,P_A+A一GEk,2,A.+A一GEl,3,m,(9)A+A一Gi,k,妇A2+AGi,k,f,m图7E中比B一中多的待定参数是A,A,A,(F)和(F)(图7),显然,若A,A,A,().和(F)可由五一中别的待定参数表示,则E的秩等于五一的秩加5.式(4)表明只有Gi,k,和Gi,k,m中有待定参数(F).和(F).,由假定具有A属性+18o.+岛18o.所以a(,J)./a()+a(,J.)./a()0a(,J)./a(N:)+厶)da(N)0固此,式(9)可写成A一GD,l,一AA2一G-k,2,f,一A5A3一GEl,3,m,一A65期张彩明等:散乱数据点的三次多项式插值421(xk一)(F).+(一)(F)一g(xz)(F)+(yz)(F).一g2(10)其中g一Dm(A+A2)/4一/(a(L)/O(N)+a(工)/O(N2)一3F.g=D(A+A,)/46r/(a()./a(Nz)十a(工.):/O(N)一3F一3F+(薯一4)(F)+(一)(F)(a(,J)/O(N-),+a(工)/O(N)一3F+(墨)(F),+(一)()口(厶)/a()+厶)/a()显然,式(10)的秩是5.注意到L厶一1(见图7),所以E的秩E一的秩+53(一1)一2Ll一3+53一2厶一3证毕4构造P,)在具有A属性的三角网格上构造尸(,),有+工+1个待定参数可自由取值,虽然这+三+1个参数不足以使尸,y)C连续,但可使尸,j,)更为光顺或者实现其它目标.下面给出一种构造尸(z,)的方法,方法描述如下:对于A()和(F)01,2,M,=1,2,),设B一At=1,2,且1=(F)(11)B+一(F)1,2,由式(4)一(6)知,尸,)应满足的约束方程是C1.】B1+C】.2B2+C1.EBE+C1.uBu:Cl,U+lC2.1B1+C2.2B2+C2.EB+C2.vBvC2,U+l(12)CEl1BI+CB2+C,B+CE.且一C,U+l设经高斯消去法计算,式(12)可写成B=C.+B+CI,+B.E斗+C.B0+c,uB0+C1l+B=c,+B+c,+B2+cuB+c,uB+c?(13)B=c.E+lB1+C.2B+2+c,uB+c.B+c.+其中BB】,B2,Bi=1,2,【,为简单,我们对式(13)增加等式BB(JE+1,E+2,U)并定义c一;E+1,E+2,u,JE+1.E+2,u,u+1这样式(13)可写成422计算机辅助设计与图形学1998正UB:一.B+c.i一1.2,U(14)j-+】由式(7)知,uE一+L+1,注意到式(11),所以式(】4)又可写成(F)一b+m一其中b一占r一1,2,+L+1at,U,.,C,i一1,2,U,一E+1,E+2,+1)把式(15)代入式(2)得十(v日1:bf+vN+Ubf+UNL+21(I5)dr.+d一L+2(L)(L)(L)l(16)l-】由上述讨论可知,若P(z,)由式(15)和(16)定义的曲面片构成,则对任意的6,(r-1,2.,-v+L+1),P(z,)是C连续的.下面讨论如何决定b(r一1,2,+十1).设7是序号为的三角形,的相邻顶点是PPP(图8),T上定义的曲面片是,),从逼近的角度讲,最好选择b使Fc(z,Y)一FI(p-t,z,f)尽可能小b可用最小二乘法确定,令图8R(6:,b2,.6+】)一(F(,Y)一F)=(32(L(z,)F)+(一L?2,3-N+L+】+.(舢+累加平方差R(6,b,bx+)(t-1,2,),得S(b】,b,bN-Lr+L)一R(6l,b一?,+】)(6l,bz,bu+L+)是变量b(r-1,2,N+L+1)的二次函数,b可由下列方程组决定M一+一A2一p+一U+6+F+_FL23.LF)6喙rL_l一+6+r一+5期张影明等:散乱数据点的三次多项式插值423(6b,b+一)/描一0r一1,2,+L+l(17)显然,若给定插值点是一三次曲面F(x,)上的点,则有一组(r=1,2,Nq-L,+1),使式(16)中的(,)满足F(,)=F,),且使s(,b,bu+)=0.因为式(17)有唯一解,所以得下面定理:定理3若各三角形上的曲面片由式(16)和(17)定义,则P(,)C可由各三角形上的曲面片拼和而成,P,)的多项式准确集是._y,zy,xy,Y.,Y,15例子本节的散乱数据点是由实际测量得到的某地区的一些地下水位点,散乱点在xy平面上划分成三角形网格如图9所示.在图10中,图(a)在每个三角形上用平面插值给定散乱点,图(b)是用本文方法构造的曲面插值给定散乱点.图9圈10计算机辅助设计与图形学1998盎6结论本文方法可用于对散乱分布的空间数据点插值,所构造的插值曲面具有三次多项式精度.对于用分片三次多项式曲面作插值函数,这已是最好的多项式插值精度.构造整体曲面需解一个线性方程组,当给定数据点个数不太大时,可用列主元素消除法获得精确解.参考文献1WangCY,ZhangCM.Polynomialofdegreefourinterpolariortontriangles.JournalofComputationalMathematics,1991,1912):1551612PowellMJD,binMA.Piecewisequadraticapproximations.ntanglesACMTransactJortOnMathematicalSoftware,1977,3(413l63253CiarletPG.SurlelementdecloughetTocher,RevueFrancaisedAotomatique,iMormatique.Rechercheoperationelle,1974,R一2,:19274AlfeldP.AdiscreteCinterpolantfortetrahedradata.RockyMouninJMath,l984,14:4165BarnhillRE,FarinG.Cquinticinterpolationovertriangles:Twoexplicitrepresentations.I眦JforNumMethinEng,l981,17:179317786DN,LevinD?RippaSSurhceinterpolationandsmoothingbythinplatesplines.In:ChuiCK.SchukerLL,WAJ13.5,ApproximatlonTheory1V,NewYork,AcademicPress,1983,4-ti497FarinG.SmoothInterpolationtOScattered3DDataBarnhillRE.BoehmW_eds,SurfaCeSinComputerAidedGeometricDesigrt,Aterdam,NorthHolland,l983,43