OpenCasCade点向曲面投影重新实现OCC中重新实现点向曲面投影计算原理迭代初值估计边界处理OCC中重新实现点向曲面投影最近使用OCC提供的GeomAPI_ProjectPointOnSurf计算点到曲面的投影发现速度比较慢。查阅了论坛和源码发现求极值的方法有Extrema_ExtAlgo_Tree和Extrema_ExtAlgo_Grad两种分别在用Extrema_GenExtPS计算极值前估算大致的投影参数。默认的Extrema_ExtAlgo_Grad较慢改为Extrema_ExtAlgo_Tree后速度有明显提升但在特殊条件下会计算不出投影点。为了补上这个坑就重新实现了这个功能在此做个记录。计算原理计算原理比较简单最小化曲面和点之间的距离F(u,v) |S(u,v) - P|^2。根据这个目标函数求一阶导和二阶导即可源码中的Extrema_FuncPSNorm.cxx文件包含了相似的实现下面为了方便直接贴出math_MultipleVarFunctionWithHessian的派生源码classPointOnSurfaceDistance:publicmath_MultipleVarFunctionWithHessian{public:PointOnSurfaceDistance(constHandle(Geom_Surface)surface,constgp_Pntpoint):mySurface(surface),myPoint(point){}virtualStandard_IntegerNbVariables()const{return2;}virtualStandard_BooleanValue(constmath_VectorX,Standard_RealF){gp_Pnt S;mySurface-D0(X(1),X(2),S);FmyPoint.SquareDistance(S);returnStandard_True;}virtualStandard_BooleanGradient(constmath_VectorX,math_VectorG){gp_Pnt P;gp_Vec dU,dV;mySurface-D1(X(1),X(2),P,dU,dV);gp_Vecdiff(myPoint,P);G(1)2.0*diff.Dot(dU);G(2)2.0*diff.Dot(dV);returnStandard_True;}virtualStandard_BooleanValues(constmath_VectorX,Standard_RealF,math_VectorG){Value(X,F);Gradient(X,G);returnStandard_True;}virtualStandard_BooleanValues(constmath_VectorX,Standard_RealF,math_VectorG,math_MatrixH){gp_Pnt P;gp_Vec dU,dV;gp_Vec dUU,dUV,dVV;mySurface-D2(X(1),X(2),P,dU,dV,dUU,dVV,dUV);gp_Vecdiff(myPoint,P);Fdiff.SquareMagnitude();G(1)2.0*diff.Dot(dU);G(2)2.0*diff.Dot(dV);H(1,1)2.0*(dU.Dot(dU)diff.Dot(dUU));H(1,2)2.0*(dU.Dot(dV)diff.Dot(dUV));H(2,1)H(1,2);H(2,2)2.0*(dV.Dot(dV)diff.Dot(dVV));returnStandard_True;}protected:Handle(Geom_Surface)mySurface;gp_Pnt myPoint;};迭代初值估计有了上面的math_Function就可以选择相应的求解器进行迭代计算了这里为了充分利用二阶导选了牛顿迭代math_NewtonMinimum。但是牛顿迭代对初值比较敏感尤其是我用来测试的组合曲面初值给的不好的情况投影结果依然差很多。Extrema_ExtAlgo_Tree用的是UB树可以直接搬运Extrema_GenExtPS的源码但是我测了曲面10x10采样添加到树中查询的效率不咋样后面干脆换成BVH树。关于BVH树的介绍可以参考eryar大佬的博客。BVH树的查询速度还是很给力的想要更快的话可以用NCollection_CellFilter进行实现。下面还是直接贴BVH树查询点投影大致参数代码#includeBVH_BoxSet.hxx#includeBVH_Box.hxx#includeBVH_LinearBuilder.hxx#includeBVH_Traverse.hxxstructCell{Standard_Real uMin,uMax,vMin,vMax;gp_Pnt p;};classShapeSelector:publicBVH_TraverseStandard_Real,3,BVH_BoxSetStandard_Real,3,Cell,Standard_Boolean{public:ShapeSelector():myMinDist(RealLast()){}voidSetBox(constgp_Pntp){myTargetp;}CellGetSelected()const{returnmySol;}public://! Defines the rules for node rejection by bounding boxvirtualStandard_BooleanRejectNode(constBVH_Vec3dtheCornerMin,constBVH_Vec3dtheCornerMax,Standard_BooleantheIsInside)constStandard_OVERRIDE{gp_Pntp(theCornerMin.x(),theCornerMin.y(),theCornerMin.z());returnp.Distance(myTarget)myMinDist;}virtualStandard_BooleanAcceptMetric(constStandard_BooleantheIsInside)constStandard_OVERRIDE{returntheIsInside;}virtualStandard_BooleanAccept(constStandard_Integer theIndex,constStandard_BooleantheIsInside)Standard_OVERRIDE{BVH_BoxStandard_Real,3boxmyBVHSet-Box(theIndex);gp_Pntp(box.CornerMin()[0],box.CornerMin()[1],box.CornerMin()[2]);Standard_Real disp.Distance(myTarget);if(dismyMinDist)returnStandard_False;myMinDistdis;mySolmyBVHSet-Element(theIndex);returnStandard_True;}protected:Standard_Real myMinDist;gp_Pnt myTarget;Cell mySol;};classSurfaceCellFilter{public:SurfaceCellFilter(){};voidInit(constHandle(Geom_Surface)surface,intuDiv,intvDiv){m_surfacesurface;buildCells(uDiv,vDiv);}gp_Pnt2dQueryNearestUV(constgp_Pntquery){ShapeSelector selector;selector.SetBox(query);selector.SetBVHSet(m_boxSet);intnumselector.Select();Cell cm_boxSet.Element(0);if(num0)cselector.GetSelected();returngp_Pnt2d(0.5*(c.uMinc.uMax),0.5*(c.vMinc.vMax));}private:Handle(Geom_Surface)m_surface;BVH_BoxSetStandard_Real,3,Cellm_boxSet;voidbuildCells(intuDiv,intvDiv){usingBox3DBVH_BoxStandard_Real,3;Standard_Real u1,u2,v1,v2;m_surface-Bounds(u1,u2,v1,v2);Standard_Real du(u2-u1)/uDiv;Standard_Real dv(v2-v1)/vDiv;BVH_LinearBuilderStandard_Real,3builder;m_boxSet.SetBuilder(builder);m_boxSet.SetSize(uDiv*vDiv);for(inti0;iuDiv;i){for(intj0;jvDiv;j){Cell c;c.uMinu1i*du;c.uMaxc.uMindu;c.vMinv1j*dv;c.vMaxc.vMindv;c.pm_surface-Value(0.5*(c.uMinc.uMax),0.5*(c.vMinc.vMax));Box3D b;b.CornerMin()[0]c.p.X();b.CornerMin()[1]c.p.Y();b.CornerMin()[2]c.p.Z();b.CornerMax()[0]c.p.X();b.CornerMax()[1]c.p.Y();b.CornerMax()[2]c.p.Z();m_boxSet.Add(c,b);}}m_boxSet.Build();}};边界处理由于我测试的曲面有一段是平的或许是目标函数非凸数学基础差轻喷在这个区域的Hessian矩阵有负的特征值这会导致搜索的曲面参数超出范围因此需要设置math_NewtonMinimum的左右限。测试的时候和Extrema_ExtAlgo_Grad的结果对比发现投影结果在曲面边缘处的计算经常出现错误将牛顿迭代的搜索范围放宽之后可以得到缓解但是不能根本上解决计算错误的问题。由于任务着急直接用GeomLib::ExtendSurfByLength把投影曲面拓宽计算了这里埋个坑后面有空再解决(●’◡’●)目前测试计算速度能比Extrema_ExtAlgo_Tree的方法快上30%左右精度也没有问题。希望能帮到和我有一样问题的OCC用户也期待官方能早点把这个问题给优化好毕竟吐槽点向曲面投影慢的帖子早就在官方论坛发出来了。