最小二乘法平面拟合.docx
《最小二乘法平面拟合.docx》由会员分享,可在线阅读,更多相关《最小二乘法平面拟合.docx(14页珍藏版)》请在冰豆网上搜索。
最小二乘法平面拟合
最小二乘法平面拟合
在介绍平面拟合之前,我先给大家介绍一下有关平面的相关知识(相关介绍来自QVPak3D,日本三丰)
DefinitionofthePlaneFeature
Aplanefeatureisreportedastheprojectionofthecentroidofthepointsusedtofittheplane,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’snormalvectorandthex-axisofthecurrentcoordinatesystem.(X-Angle=arccosinek).Thex-angleisapositivenumberbetween0and180degrees.
Y-angle:
Theanglebetweentheplane’snormalvectorandthey-axisofthecurrentcoordinatesystem.(Y-Angle=arccosinel).They-angleisapositivenumberbetween0and180degrees.
Z-angle:
Theanglebetweentheplane’snormalvectorandthez-axisofthecurrentcoordinatesystem.(Z-Angle=arccosinem).Thez-angleisapositivenumberbetween0and180degrees.
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.
平面相关知识:
算法如下:
VB源代码:
OptionExplicit
PublicConstPI=3.1415926535897
PublicTypetagPoint
xAsDouble
yAsDouble
zAsDouble
EndType
PublicTypetagLine2D
kAsDouble 'Slope,Kis the"K"of"y=kx+b"
bAsDouble 'intercept,Bis the"B"of"y=kx+b"
AngleAsDouble 'arctg(k) '0to180deg
StraightnessAsDouble
RSQAsDouble
EndType
PublicTypetagLine3D
'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}
xAsDouble
yAsDouble
zAsDouble
x0AsDouble'AddedonJul.1st,2009
y0AsDouble'AddedonJul.1st,2009
z0AsDouble
mAsDouble '(x-x0)/m=(y-y0)/n=(z-z0)/p,sameas:
z=a+bx,z=c+dy
nAsDouble
pAsDouble 'vectors={m,n,p}
AngleAsDouble
xAnAsDouble
yAnAsDouble
zAnAsDouble
StraightnessAsDouble
MinStAsDouble
MaxStAsDouble
EndType
PublicTypetagCircle
xAsDouble
yAsDouble
zAsDouble
rAsDouble
dAsDouble
CircularAsDouble
EndType
PublicTypetagVector
aAsDouble
bAsDouble
cAsDouble
EndType
TypetagPlane
xAsDouble 'Thedistancefromtheorigintothecentroid,asmeasuredalongthex-axis.
yAsDouble 'Thedistancefromtheorigintothecentroid,asmeasuredalongthey-axis
zAsDouble 'Thedistancefromtheorigintothecentroid,asmeasuredalongthez-axis
'Z+A*x+B*y+C=0 z'scoefficientisjust1
axAsDouble 'coefficient ofX
byAsDouble 'coefficient ofY
czAsDouble
dAsDouble 'ConstantC
AngleAsDouble
xAnAsDouble
YAnasDouble
zAnAsDouble
FlatAsDouble
MinFlatAsDouble
EndType
'*************************************************************
' 函数名:
PlaneSet
' 功能:
求拟合平面(相关参数)
' 参数:
dataRaw -tagPoint型 自定义点类型(x,y,z),数组
' PlaneSet–tagPlane 其值返回为平面圆相关参数
' 返回值:
Long型,失败为0,成功为-1
'*************************************************************
PublicFunctionPlaneSet(dataRaw()AstagPoint,PlaneAstagPlane)AsLong'z+Ax+BY+C=0
PlaneSet=0
DimlbAsLong,ubAsLong,MaxFAsDouble,MinFAsDouble,tmpAsDouble
lb=LBound(dataRaw):
ub=UBound(dataRaw)
Ifub-lb+1<4ThenExitFunction
DimiAsLong,nAsLong
n=ub-lb+1
DimxAsDouble,yAsDouble,zAsDouble
DimXYAsDouble,XZAsDouble,YZAsDouble
DimX2AsDouble,Y2AsDouble
DimaAsDouble,bAsDouble,cAsDouble,dAsDouble
Dima1AsDouble,b1AsDouble,z1AsDouble
Dima2AsDouble,b2AsDouble,z2AsDouble
Dimn1AstagVector'{.ax,by,1} s1
Dimn2AstagVector'{0,0,N}XYplane s2
Dimn3AstagVector'lineprojectedplane
DimSLineAstagVector
DimxLineAstagVector'{1,0,0}
DimyLineAstagVector'{1,0,0}
DimzLineAstagVector'{1,0,0}
DimVectorPlaneAstagVector
xLine.a=1:
xLine.b=0:
xLine.c=0
yLine.a=0:
yLine.b=1:
yLine.c=0
zLine.a=0:
zLine.b=0:
zLine.c=1
Fori=lbToub
WithdataRaw(i)
x=x+.x
y=y+.y
z=z+.z
XY=XY+.x*.y
XZ=XZ+.x*.z
YZ=YZ+.y*.z
X2=X2+.x^2
Y2=Y2+.y^2
EndWith
Nexti
z1=n*XZ-x*z 'e=z-Ax-By-C z=Ax+By+D
a1=n*X2-x^2
b1=n*XY-x*y
z2=n*YZ-y*z
a2=n*XY-x*y
b2=n*Y2-y^2
a=(z1*b2-z2*b1)/(a1*b2-a2*b1)
b=(a1*z2-a2*z1)/(a1*b2-a2*b1)
c=1
d=(z-a*x-b*y)/n
WithPlane
.x=x/(ub-lb+1)
.y=y/(ub-lb+1)
.z=z/(ub-lb+1)
'sum(Mi*Ri)/sum(Mi),Miismass. here Miissetedtobe1and.zisjusttheaverageofz
.Ax=-a
.By=-b
.Cz=1
.d=-d 'z=Ax+By+D-----Ax+By+Z+D=0
VectorPlane.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=1
SLine.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-.Angle
Fori=lbToub
PointToPlanedataRaw(i),Plane,tmp,0
Ifi=lbThen
MaxF=tmp:
MinF=tmp
Else
IfMaxF IfMinF>tmpThenMinF=tmp
EndIf
Nexti
.MaxFlat=MaxF
.MinFlat=MinF
.Flat=MaxF-MinF
EndWith
PlaneSet=-1
EndFunction
' 函数名:
VectorMulti
' 功能:
求两向量的积,结果存放于参数rtV3中
' 参数:
v1 -tagVector
' v2 -tagVector
' rtV3 -tagVector
'*************************************************************
PublicFunctionVectorMulti(v1AstagVector,v2AstagVector,rtv3AstagVector)AsLong
'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) 'k
rtv3.a=(v1.b*v2.c-v2.b*v1.c)
rtv3.b=-(