最小二乘法平面拟合

最小二乘法平面拟合在介绍平面拟合之前，我先给大家介绍一下有关平面的相关知识(相关介绍来自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.Zangle誉'NormaltoPkne3〕YAngle.^.AngleAngle口Flatness: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.平面相关知识：1槽武蓟讽

万理i：4r* idi*卫'^4+CJs-*■ ■Q1^j/K.FTJFAx-i呼亠Et+H+Grr0E Jt PT-ATJ4片⑪尸砒*耳b*冃q向■%*"血思册己知点为:p心一“和已烈亘取£方程为_一匚&=上二Z二匚空m n p苴中J*0」=Iffi以直城「乂主=・2二口£p tn tt陷、儿冷H均活知舞答如下："+H:卄童啟£为平面_-吋+ +z+^■:-dKij』〕点片过平屈苴平湖前向黄为片={*.M}工平面中的点则：d=Te1航十门■片十齣h算出了平间方理.直线号平面相交于职何*g将F跷化为善議方裡$f・ffs^=n'I+j>rI…优人平tfl方程稈：d1Rt+xa)*hb(n11+yo)+i+j=0算法如下:Ax+》+良+0=0Ax+Ax+》+良+0=0Ax+By~1■抒.+D=0日=Ax^By+左+£)使其残差最水兀护二乞伍诂砂+針前二迟⑷；+妙+莖+D)m詈二艺(百兀+砂+云+口)*丁二oBC 厂迟［Ay+®氏+D)勺二0N:Fro?^A丄十&+貶〉(4H1) .吃宓-2?2>=-电/-壬送£-加吃旷2>2>)'斷二吃砲-迟犬>冲1二号？-为送诃二吃砂■迟亡吃旷2>2>=_4£才迟兀&卜曲吃^_2>2>咕吃严2>2?恣=2?片迟吃>3皓吃y-2>J—汛二£坷+占直I-毛二眉遹十月岛Ax+5'.嗫+ “+D',r=04迟心+召迟戸+乞坪+刀乞，=0

贬> 亠乞"+门乞1=。贬:/+眨>+迟矽+口乞"0加(1)

运◎+迄b+乞坪+□》>=0加⑵

◎忑+吃y+Dny踽h(3)_5__叼勺…咛］&%—冬A2＞十迟卄ED

~irVB源代码：OptionExplicitPublicConstPI=3.1415926535897PublicTypetagPointxAsDoubleyAsDoublezAsDoubleEndTypePublicTypetagLine2DkAsDouble 'Slope,Kisthe〃K〃of〃y二kx+b〃bAsDouble 'intercept,Bisthe〃B〃of〃y二kx+b〃AngleAsDouble '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-axis

zAsDouble'Thedistancefromtheorigintothecentroid,asmeasuredalongthez-axis+CzAsDouble'Thedistancefromtheorigintothecentroid,asmeasuredalongthez-axis+C=0z's'coefficient'coefficient'Z+A*x+B*yaxAsDoublebyAsDoubleczAsDoubledAsDoubleAngleAsDoublexAnAsDoubleYAnasDoublezAnAsDoubleFlatAsDoublecoefficientisjust1ofXofY'ConstantCMinFlatAsDoubleEndType'*************************************************************' 函数名：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,zAsDoubleDimXYAsDouble,XZAsDouble,YZAsDoubleDimX2AsDouble,Y2AsDoubleDimaAsDouble,bAsDouble,cAsDouble,dAsDoubleDima1AsDouble,b1AsDouble,z1AsDoubleDima2AsDouble,b2AsDouble,z2AsDoubleDimn1AstagVector'{.ax,by,1}s1Dimn2AstagVector'{0,0,N}XYplanes2Dimn3AstagVector'lineprojectedplaneDimSLineAstagVectorDimxLineAstagVector'{1,0,0}DimyLineAstagVector'{1,0,0}DimzLineAstagVector'{1,0,0}DimVectorPlaneAstagVector

xLine.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+.x2Y2=Y2+.y2EndWithNextiz1=n*XZ-x*z'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入2a=(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.hereMiissetedtobe1and.zisjusttheaverageofz.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)AsLongrtv3.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)'MsgBoxtmptmpSqrl=Sqr(vl.a八2+vl.b八2+vl.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 -Long 0：点到平面距离(有正有负)，其它值：点到平面距离(绝对值)' 返回值：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'**************************************