最小二乘法平面擬合

最小二乘法平面擬合在介紹平面擬合之前，我先給大家介紹一下有關(guān)平面的相關(guān)知識(相關(guān)介紹來自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譽(yù)'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.平面相關(guān)知識：1槽武薊諷

萬理i：4r* idi*衛(wèi)'^4+CJs-*■ ■Q1^j/K.FTJFAx-i呼亠Et+H+Grr0E Jt PT-ATJ4片?尸砒*耳b*冃q向■%*"血思冊己知點(diǎn)為:p心一“和已烈亙?nèi)　攴匠虨開一匚&=上二Z二匚空m n p苴中J*0」=Iffi以直城「乂主=?2二口￡p tn tt陷、兒冷H均活知舞答如下："+H:卄童啟￡為平面_-吋+ +z+^■:-dKij』〕點(diǎn)片過平屈苴平湖前向黃為片={*.M}工平面中的點(diǎn)則：d=Te1航十門■片十齣h算出了平間方理.直線號平面相交于職何*g將F蹺化為善議方裡$f?ffs^=n'I+j>rI…優(yōu)人平tfl方程稈：d1Rt+xa)*hb(n11+yo)+i+j=0算法如下:Ax+》+良+0=0Ax+Ax+》+良+0=0Ax+By~1■抒.+D=0日=Ax^By+左+￡)使其殘差最水兀護(hù)二乞伍詁砂+針前二遲⑷；+妙+莖+D)m詈二藝(百兀+砂+云+口)*丁二oBC 廠遲［Ay+?氏+D)勺二0N:Fro?^A丄十&+貶〉(4H1) .吃宓-2?2>=-電/-壬送￡-加吃曠2>2>)'斷二吃砲-遲犬>沖1二號？-為送訶二吃砂■遲亡吃曠2>2>=_4￡才遲兀&卜曲吃^_2>2>咕吃嚴(yán)2>2?恣=2?片遲吃>3皓吃y-2>J—汛二￡坷+占直I-毛二眉遹十月島Ax+5'.囁+ “+D',r=04遲心+召遲戸+乞坪+刀乞，=0

貶> 亠乞"+門乞1=。貶:/+眨>+遲矽+口乞"0加(1)

運(yùn)◎+迄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'*************************************************************' 函數(shù)名：PlaneSet' 功能：求擬合平面（相關(guān)參數(shù)）'參數(shù)：dataRaw-tagPoint型 自定義點(diǎn)類型（x,y,z）,數(shù)組' PlaneSet-tagPlane 其值返回為平面圓相關(guān)參數(shù)'返回值：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'函數(shù)名：VectorMulti' 功能： 求兩向量的積，結(jié)果存放于參數(shù)RtV3中' 參數(shù)： 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'************************************************************'函數(shù)名：VectorN'功能： 求兩向量的數(shù)量積，結(jié)果存放于參數(shù)RtV3中'參數(shù)：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'*************************************************************'函數(shù)名：Intersect'功能：求兩向量之間的夾角'參數(shù)：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'*************************************************************'函數(shù)名：PointToPlane'功能：求點(diǎn)到平面的距離'參數(shù)： dataRaw-tagPoint型 自定義點(diǎn)類型(x,y,z)' Plane-tagPlaneDouble' RtnDistance-Double 其值為點(diǎn)到直線的距離.' AbsDist -Long 0：點(diǎn)到平面距離(有正有負(fù))，其它值：點(diǎn)到平面距離(絕對值)' 返回值：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'**************************************