最小二乘法平面拟合

最小二乘法平面拟合在介绍平面拟合之前，我先给大家介绍一下有关平面的相关知识(相关介绍来自QVPak3D,日本三丰)DefinitionofthePlaneFeatureAplanefeatureisreportedastheprojectionofthecentroidofthepointsusedtofittheplane,whichisprojectedontotheplanefeature,ameasurementofthedirectionmeasuredasanangle,ameasurementoftheflatnessoftheplaneandameasurementoftheparallelismoftheplane.IfmeasuredinaCartesiancoordinatesystem,thecoordinatesoftheplane'scentroidarereportedasfollows:X:Thedistancefromtheorigintothecentroid,asmeasuredalongthex-axis.Y:Thedistancefromtheorigintothecentroid,asmeasuredalongthey-axis.Z:Thedistancefromtheorigintothecentroid,asmeasuredalongthez-axis.IfmeasuredinaCylindricalcoordinatesystem,thecoordinatesoftheplane'scentroidarereportedasfollows:R:Thedistancefromthez-axisofthecoordinatesystemtothecentroid,asmeasuredwithinaplanewhichcontainsthecentroidandisorthogonaltothez-axisofthecoordinatesystem.A:Thedirection,measuredasanangle,betweenareferenceradiusvectorandaradiusvectorthatcontainsthecentroidandisprojectedontothexy-plane.Thereferenceradiusvectormaybeconsideredtobethex-axis.Z:Theheightfromtheorigintothecentroidinthecylindricalcoordinatesystem,asmeasuredalongthez-axis.Theotherattributesoftheplanefeatureare:Angle:Theanglebetweentheprojectionoftheplane’snormalvectorontothexy-planeandthex-axisofthecurrentcoordinatesystem.X-angle:Theanglebetweentheplane’snormalvectorandthe-xaxisofthecurrentcoordinatesystem.(X-Angle=arccosinek).Thex-angleisapositivenumberbetween0and180degrees.Y-angle:Theanglebetweentheplane’snormalvectorandthey-axisofthecurrentcoordinatesystem.(Y-Angle=arccosinel).They-angleisapositivenumberbetween0and180degrees.Z-angle:Theanglebetweentheplane’snormalvectorandthe-zaxisofthecurrentcoordinatesystem.(Z-Angle=arccosinem).Thez-angleisapositivenumberbetween0and180degrees.21.ZangleZ-；NormaltoPkne部YAngleAngleAngleAFlatness:Flatnessisaconditionforwhichanelementofasurfaceisinaplane.Flatnessisreportedasthewidthofthezoneformedbytwoclosestparallelplanesthatfullycontainthepointsetusedtofittheplanefeature.Avalueofzeroindicatesperfectflatness.Flatness(minimum):Thedistancefromthefittedplanetothemeasuredpointfarthestbelowthefittedplaininthepointset.Aboveandbelowaredeterminedbythedirectionoftheplanevector.SeeExplanationofMax/Mindistanceindifferentcases.Flatness(maximum):Thedistancefromthefittedplanetothemeasuredpointfarthestabovethefittedplaininthepointset.Aboveandbelowaredeterminedbythedirectionoftheplanevector.SeeExplanationofMax/Mindistanceindifferentcases.Parallelism:Theconditionofafeature,projectedtoacertainplane,beingequidistantatallelementsfromadatum(reference).Quantitatively,parallelismisdefinedastheabsolutedistantdifferencebetweenthefarthestandclosestpointsfromthedatum.Parallelismisevaluatedrelativetoareferencelineorxy-plane.Whenevaluatingasetofpointswithareferenceline,parallelismusestheprojectionsofthepointsandreferencelineontothexy-planeinthecurrentcoordinatesystem,orz/refplanefeature,andisspecifiedasazonetolerance.Thez/refplanefeatureisaplaneincludingthereferencelineandparallelto(orincluding)thez-axis.Whenevaluatingasetofpointswithaxy-plane,parallelismiscalculatedinthree-dimensionalspace.Parallelism(minimum):Thedistancefromthereferencedlineorplanetothepointinthepointsetwiththeleastvalue(leastpositivevalueifallevaluatedpointsarepositive,ormostnegativevalueifevaluatedpointsincludenegativevalues).SeeExplanationofMax/Mindistanceindifferentcases.Parallelism(maximum):Thedistancefromthereferencelineorplanetothepointinthepointsetwiththegreatestvalue(mostpositivevalueifevaluatedpointsincludepositivevalues,orleastnegativevalueifallevaluatedpointsarenegative).SeeExplanationofMax/Mindistanceindifferentcases.平面相关知识：o°J^+Thr+ +JIk+(E+LW)『阳:留Si早限由YM…1十LT=rF&'莒«=kJ哉乌联畚蚌聒袖豆鲜{**"■咔j"辑土蚤牌里业目寤覃割华Bl业＞1审喜1N斗顼/4w»=F：响馨岬山凰出*』,{W昌「叩叫=雄尊回风崖池直1!击??#W宵我陶0・F+E+£乌十k_厘立峰■『骑目J|『K444-4I4 阵暮曲N 珅 .古=亏二7==?「希单物蛔1=小0『*走箕亏二=可与二盘上：*雷£乌丫争瓦峰会。一.£■ A—n Jl—JL〔京侦那蚌罕哮£海w耳职忻w-j址*责土］也%*册-Xfr-sW-Xfr-sW-X*2-a0土打+』十=jt¥。■t/TwK*1址十/V_.i* ,p-tf*A拍♦iy][jp♦勺i由*明咨里站革g回g'•闻蛀皂g理割攵00古

祖町¥船i算法如下:AxAx+导+/+/?=0Ax+By+声.+D=0E=Ax^By+宓+也使其残差最小工F=£（虫+肉+窄+呀函=z（血;+导+就+巾）咤=。告二工（虫+力+了+□尸况二。3c 「=已（出+为+1「）勾=0N:Aj,¥S快+%+□r=Q旦工节+E2>"+工好+口2>=。

应2>+召2>+2#+口二1=，住+吃八以（i）应工仲+3二尸+工炉+□?>=CW⑵吃二+吃】+乏>+必=。AW伽澄）西口球3）农»+查>+眨？N阳1）吃球-»2>'=-也？-£工为）一曲吃或广£公^邑=吃遂-2>眼占1=以-»&初=吃殆£心A&^-2>2>=_4S/ZA2>）-&Aiy_2>2>气=套眼农=£寸1^心陀=宪了一2>[一苴1=&昭+B*[-弓二占丹+B%如磕挈^1-蜂1g=_一气备戌也!一叱A2>十位>十登>2_-'—VB源代码：OptionExplicitPublicConstPI=3.97PublicTypetagPointxAsDoubleyAsDoublezAsDoubleEndTypePublicTypetagLine2DkAsDouble 'Slope,Kisthe〃K〃of〃y=kx+b〃bAsDouble 'intercept,Bisthe〃B〃of"y=kx+bAngleAsDouble'arctg(k),0to180degStraightnessAsDoubleRSQAsDoubleEndTypePublicTypetagLine3D,3Dline'sformulaisshowingasfollowing.'1)|Ax+By+Z+D=0' |A1x+B1y+z+D1=0'2)(x-x0)/m=(y-y0)/n=(z-z0)/p (x-x0)/m=(y-y0)/n=z/1'3)x=mt+x0,y=nt+y0,z=pt+z0,Onlypoint'scoordinateis(a+b*Z,c+d*Z,Z),sotheline'svectoris{b,d,1}xAsDoubleyAsDoublezAsDoublex0AsDouble'AddedonJul.1st,2009y0AsDouble'AddedonJul.1st,2009z0AsDoublemAsDouble'(x-x0)/m=(y-y0)/n=(z-z0)/p,sameas:z=a+bx,z=c+dynAsDoublepAsDouble'vectors={m,n,p}AngleAsDoublexAnAsDoubleyAnAsDoublezAnAsDoubleStraightnessAsDoubleMinStAsDoubleMaxStAsDoubleEndTypePublicTypetagCirclexAsDoubleyAsDoublezAsDoublerAsDoubledAsDoubleCircularAsDoubleEndTypePublicTypetagVectoraAsDoublebAsDoublecAsDoubleEndTypeTypetagPlanexAsDouble'Thedistancefromtheorigintothecentroid,asmeasuredalongthex-axis.yAsDouble'Thedistancefromtheorigintothecentroid,asmeasuredalongthey-axiszAsDoubleThedistancefromtheorigintothecentroid,asmeasuredalongthez-axis'Z+A*x+B*y+C=0z'scoefficientisjust1axAsDouble 'coefficient ofXbyAsDouble 'coefficient ofYczAsDoubledAsDouble'ConstantCAngleAsDoublexAnAsDoubleYAnasDoublezAnAsDoubleFlatAsDoubleMinFlatAsDoubleEndType'*************************************************************'函数名：PlaneSet'功能：求拟合平面(相关参数)'参数：dataRaw-tagPoint型 自定义点类型(x,y,z),数组' PlaneSet-tagPlane 其值返回为平面圆相关参数' 返回值：Long型，失败为0,成功为-1'*************************************************************PublicFunctionPlaneSet(dataRaw()AstagPoint,PlaneAstagPlane)AsLong'z+Ax+BY+C=0PlaneSet=0DimlbAsLong,ubAsLong,MaxFAsDouble,MinFAsDouble,tmpAsDoublelb=LBound(dataRaw):ub=UBound(dataRaw)Ifub—lb+1<4ThenExitFunctionDimiAsLong,nAsLongn=ub-lb+1DimxAsDouble,yAsDouble,zAsDoubleDimXY As Double, XZ As Double, YZ As DoubleDimX2 As Double, Y2 As DoubleDimaAsDouble,bAsDouble,cAsDouble,dAsDoubleDima1 As Double, b1 As Double, z1 As DoubleDima2 As Double, b2 As Double, z2 As DoubleDimn1AstagVector'{.ax,by,1}siDimn2AstagVector'{0,0,N}XYplanes2Dimn3AstagVector'lineprojectedplaneDimSLineAstagVectorDimxLineAstagVector'{1,0,0}DimyLineAstagVector'{1,0,0}DimzLineAstagVector'{1,0,0}DimVectorPlaneAstagVector

TOC\o"1-5"\h\zxLine.a= 1: xLine.b = 0: xLine.c = 0yLine.a= 0: yLine.b = 1: yLine.c = 0zLine.a= 0: zLine.b = 0: zLine.c = 1Fori=lbToubWithdataRaw(i)x=x+.xy=y+.yz=z+.zXY=XY+.x*.yXZ=XZ+.x*.zYZ=YZ+.y*.zX2=X2+.x"2Y2=Y2+.y"2EndWithNextiz1=n*XZz，e=z-Ax-By-Cz=Ax+By+Da1=n*X2-x"2b1=n*XY-x*yz2=n*YZ-y*za2=n*XY-x*yb2=n*Y2-y"2x*hereMiissetedtobe1and.zhereMiissetedtobe1and.za=(z1*b2-z2*b1)/(a1*b2-a2*b1)b=(a1*z2-a2*z1)/(a1*b2-a2*b1)c=1d=(z-a*x-b*y)/nWithPlane.x=x/(ub-lb+1).y=y/(ub-lb+1).z=z/(ub-lb+1)'sum(Mi*Ri)/sum(Mi),Miismass.isjusttheaverageofz.Ax=-a.By=-b.Cz=1.d =-d 'z=Ax+By+D Ax+By+Z+D=0VectorPlane.a=.Ax:VectorPlane.b=.By:VectorPlane.c=1.xAn=Intersect(VectorPlane,xLine).yAn=Intersect(VectorPlane,yLine).zAn=Intersect(VectorPlane,zLine)n1.a=.Ax:n1.b=.By:n1.c=1SLine.a=.Ax:SLine.b=.By:SLine.c=0.Angle=Intersect(xLine,SLine),(xLine.A*SLine.A+xLine.A*SLine.B+xLine.C*SLine.C)IfSLine.b<0Then.Angle=360-.AngleFori=lbToubPointToPlanedataRaw(i),Plane,tmp,0Ifi=lbThenMaxF=tmp:MinF=tmpElseIfMaxF<tmpThenMaxF=tmpIfMinF>tmpThenMinF=tmpEndIfNexti.MaxFlat=MaxF.MinFlat=MinF.Flat=MaxF-MinFEndWithPlaneSet=-1EndFunction'函数名：VectorMulti'功能： 求两向量的积，结果存放于参数rtV3中'参数：v1-tagVector' v2-tagVector' rtV3-tagVector'*************************************************************PublicFunctionVectorMulti(v1AstagVector,v2AstagVector,rtv3AstagVector)AsLongTOC\o"1-5"\h\z,rtV3=v1*v2=|i j k|'|v1.A v1.B v1.C|'|v2.A v2.B v2.C|'rtv3.A=(v1.B*v2.c-v2.B*v1.C)'i'rtV3.B=-(v1.A*v2.C-v2.A*v1.C)'j'rtV3.C=(v1.A*v2.B-v2.A*V1.B)'krtv3.a=(v1.b*v2.c-v2.b*v1.c)rtv3.b=-(v1.a*v2.c-v2.a*v1.c)rtv3.c=(v1.a*v2.b-v2.a*v1.b)EndFunction'************************************************************'函数名：VectorN'功能： 求两向量的数量积，结果存放于参数rtV3中' 参数： v1-tagVector' v2-tagVector' rtV3-tagVector'*************************************************************PublicFunctionVectorN(v1AstagVector,v2AstagVector,rtv3AstagVector)AsLongTOC\o"1-5"\h\zrtv3.a = v1.a * v2.artv3.b = v1.b * v2.brtv3.c = v1.c * v2.cEndFunction'*************************************************************'函数名：Intersect'功能：求两向量之间的夹角' 参数：v1-tagVector' v2-tagVector' LinePlane-long 0:表示两直线之间的夹角，其它值:表示如线与平面之间,平面与平面之间的夹角(0~90)' 返回值：Double型，单位:度.'*************************************************************PublicFunctionIntersect(v1AstagVector,v2AstagVector,OptionalLinePlaneAsLong=0)AsDouble,LinePlane0:line-line,1:line--PlaneDimtmpAsDouble,tmpSqr1AsDouble,tmpSqr2AsDoubletmp=(v1.a*v2.a+v1.b*v2.b+v1.c*v2.c),MsgBoxtmptmpSqr1=Sqr(v1.a"2+v1.b"2+v1.c"2)tmpSqr2=Sqr(v2.a"2+v2.b"2+v2.c"2)IftmpSqr1<>0ThenIftmpSqr2<>0Thentmp=tmp/tmpSqr1/tmpSqr2Elsetmp=tmp/tmpSqr1EndIfElseIftmpSqr2<>0Thentmp=tmp/tmpSqr2Elsetmp=0EndIfEndIfIfLinePlane<>0Thentmp=Abs(tmp)EndIfIf-tmp*tmp+1<>0Thentmp=Atn(-tmp/Sqr(-tmp*tmp+1))+2*Atn(1)tmp=tmp/PI*180Elsetmp=90EndIfIntersect=tmpEndFunction'*************************************************************'函数名：PointToPlane'功能：求点到平面的距离'参数：dataRaw-tagPoint型自定义点类型(x,y,z)' Plane-tagPlaneDouble' RtnDistance-Double 其值为点到直线的距离.' AbsDist -Long0：点到平面距离(有正有负)，其它值：点到平面距离(绝对值)' 返回值：Long型，失败为0,成功为-1'*************************************************************PublicFunctionPointToPlane(dataRawAstagPoint,PlaneAstagPlane,rtnDistanceAsDouble,OptionalAbsDistAsLong=0)AsLongDimiAsLong,lbAsLong,ubAsLong,tmpAsDoubleWithPlanetmp=(.ax*dataRaw.x+.by*dataRaw.y+.cz*dataRaw.z+.d)/Sqr(.ax"2+.by"2+.cz"2)IfAbsDist<>0Thentmp=Abs(tmp)EndIfEndWithrtnDistance=tmpEndFunction'*************************************************************'函数名：PointToLine'功能：求点到直线的距离