OpenCASCADE 平面求交
eryar@163.com
OpenCASCADE提供了類IntAna_QuadQuadGeo用來計算兩個二次曲面quadric（球面、圓柱面、圓錐面及平面，平面是二次曲面的特例）之間的交線。他們之間可能的結果有：
l 一個點
l 一條或兩條直線
l 一個點和一條直線
l 圓
l 橢圓
l 拋物線
l 雙曲線
將源碼結合《高等數(shù)學》、《解析幾何》等書，可以來學習如何將理論付諸實踐。本文主要介紹這個類中兩個平面求交的源碼實現(xiàn)。從源碼中也可以看出OpenCASCADE官方開發(fā)人員的編碼習慣。
將源碼列出如下：
void IntAna_QuadQuadGeo::Perform (const gp_Pln& P1,
const gp_Pln& P2,
const Standard_Real TolAng,
const Standard_Real Tol)
{
Standard_Real A1, B1, C1, D1, A2, B2, C2, D2, dist1, dist2, aMVD;
//
done=Standard_False;
param2bis=0.;
//
P1.Coefficients(A1,B1,C1,D1);
P2.Coefficients(A2,B2,C2,D2);
//
gp_Vec aVN1(A1,B1,C1);
gp_Vec aVN2(A2,B2,C2);
gp_Vec vd(aVN1.Crossed(aVN2));
//
const gp_Pnt& aLocP1=P1.Location();
const gp_Pnt& aLocP2=P2.Location();
//
dist1=A2*aLocP1.X() + B2*aLocP1.Y() + C2*aLocP1.Z() + D2;
dist2=A1*aLocP2.X() + B1*aLocP2.Y() + C1*aLocP2.Z() + D1;
//
aMVD=vd.Magnitude();
if(aMVD <=TolAng) {
// normalles are collinear - planes are same or parallel
typeres = (Abs(dist1) <= Tol && Abs(dist2) <= Tol) ? IntAna_Same
: IntAna_Empty;
}
else {
Standard_Real denom, denom2, ddenom, par1, par2;
Standard_Real X1, Y1, Z1, X2, Y2, Z2, aEps;
//
aEps=1.e-16;
denom=A1*A2 + B1*B2 + C1*C2;
denom2 = denom*denom;
ddenom = 1. - denom2;
denom = ( Abs(ddenom) <= aEps ) ? aEps : ddenom;
par1 = dist1/denom;
par2 = -dist2/denom;
gp_Vec inter1(aVN1.Crossed(vd));
gp_Vec inter2(aVN2.Crossed(vd));
X1=aLocP1.X() + par1*inter1.X();
Y1=aLocP1.Y() + par1*inter1.Y();
Z1=aLocP1.Z() + par1*inter1.Z();
X2=aLocP2.X() + par2*inter2.X();
Y2=aLocP2.Y() + par2*inter2.Y();
Z2=aLocP2.Z() + par2*inter2.Z();
pt1=gp_Pnt((X1+X2)*0.5, (Y1+Y2)*0.5, (Z1+Z2)*0.5);
dir1 = gp_Dir(vd);
typeres = IntAna_Line;
nbint = 1;
//
//-------------------------------------------------------
// When the value of the angle between the planes is small
// the origin of intersection line is computed with error
// [ ~0.0001 ] that can not br considered as small one
// e.g.
// for {A~=2.e-6, dist1=4.2e-5, dist2==1.e-4} =>
// {denom=3.4e-12, par1=12550297.6, par2=32605552.9, etc}
// So,
// the origin should be refined if it is possible
//
Standard_Real aTreshAng, aTreshDist;
//
aTreshAng=2.e-6; // 1.e-4 deg
aTreshDist=1.e-12;
//
if (aMVD < aTreshAng) {
Standard_Real aDist1, aDist2;
//
aDist1=A1*pt1.X() + B1*pt1.Y() + C1*pt1.Z() + D1;
aDist2=A2*pt1.X() + B2*pt1.Y() + C2*pt1.Z() + D2;
//
if (fabs(aDist1)>aTreshDist || fabs(aDist2)>aTreshDist) {
Standard_Boolean bIsDone, bIsParallel;
IntAna_IntConicQuad aICQ;
//
// 1.
gp_Dir aDN1(aVN1);
gp_Lin aL1(pt1, aDN1);
//
aICQ.Perform(aL1, P1, TolAng, Tol);
bIsDone=aICQ.IsDone();
if (!bIsDone) {
return;
}
//
const gp_Pnt& aPnt1=aICQ.Point(1);
//----------------------------------
// 2.
gp_Dir aDL2(dir1.Crossed(aDN1));
gp_Lin aL2(aPnt1, aDL2);
//
aICQ.Perform(aL2, P2, TolAng, Tol);
bIsDone=aICQ.IsDone();
if (!bIsDone) {
return;
}
//
bIsParallel=aICQ.IsParallel();
if (bIsParallel) {
return;
}
//
const gp_Pnt& aPnt2=aICQ.Point(1);
//
pt1=aPnt2;
}
}
}
done=Standard_True;
}
要理解這個源碼，需要知道平面的一般方程：Ax+By+Cz+D=0，兩個平面之間的夾角等概念。通過源碼，可以看出計算兩個平面之間的交線的步驟如下：
l 獲取兩個平面的一般方程的系數(shù)：A、B、C、D，其中平面的法向量（A,B,C)為單位向量；
l 將兩個平面的法向量叉乘得到的向量vd為平面交線的方向；
l 分別計算一個平面上的點到另外一個平面的距離：dist1和dist2；
l 如果向量vd的大小小于指定的精度TolAng，則認為兩個平面平行沒有交線；如果兩個距離dist1和dist2小于指定的精度Tol，則認為兩個平面是相同的（重合）；
l 計算兩個平面的夾角denom；
l 根據(jù)兩個平面的夾角計算交線上的點；
l 后面是處理兩個平面夾角很小的情況；
l 最后得到交線上的點pt1和方向dir1
其實上面求交線上點的代碼不好理解，可以換成三個平面求交點的處理更好理解，如將交線的方向作為法向得到的一個平面與那兩個平面一起計算交點，這個交點就一定在交線上，相關代碼如下：
gp_Pln P3(vd.X(), vd.Y(), vd.Z(), 0.0);
IntAna_Int3Pln aTool(P1, P2, P3);
if (aTool.IsDone())
{
pt1 = aTool.Value();
}
因為三個平面求交點是用高斯消元法解三元一次方程組，性能沒有上面的代碼好。生活中到處都是選擇題，如何抉擇是個問題啊。
為了方便大家在移動端也能看到我的博文和討論交流，現(xiàn)已注冊微信公眾號，歡迎大家掃描下方二維碼關注。
