OCCT尖点采样算法: 如何赋予 完美的“尖点” 一个半径

double xxx::CalculateRadiusFromThreePoints(const gp_Pnt& p1, const gp_Pnt& p2, const gp_Pnt& p3){// 检查三点是否重合或距离过近if (p1.IsEqual(p2, Precision::Confusion()) || p2.IsEqual(p3, Precision::Confusion()) || p1.IsEqual(p3, Precision::Confusion())){return 0.0;}gp_Vec v1(p1, p2);gp_Vec v2(p1, p3);// 检查三点是否共线if (v1.Magnitude() < Precision::Confusion() || v2.Magnitude() < Precision::Confusion() || v1.IsParallel(v2, Precision::Angular())){return 0.0;}double a = p1.Distance(p2);double b = p2.Distance(p3);double c = p3.Distance(p1);double s = (a + b + c) / 2.0;double area_squared = s * (s - a) * (s - b) * (s - c);if (area_squared <= Precision::SquareConfusion()) {return 0.0;}double area = sqrt(area_squared);return (a * b * c) / (4.0 * area);}bool xxx::AnalyzeSectionProfileBySpecification(const TopoDS_Shape& face,const double z_level,double& out_leadingEdgeRadius,double& out_trailingEdgeRadius,double& out_maxChordLength,double& out_twistAngle,double& out_maxThickness,TopoDS_Shape& out_sectionProfile){out_leadingEdgeRadius = 0.0;out_trailingEdgeRadius = 0.0;out_maxChordLength = 0.0;out_twistAngle = 0.0;out_maxThickness = 0.0;out_sectionProfile.Nullify();// 1. 输入验证if (face.IsNull() || face.ShapeType() != TopAbs_FACE) {std::cerr << \"错误: 输入的不是一个有效的面。\" << std::endl;return false;}try {// 2. 创建并执行截面操作gp_Pln sectionPlane(gp_Pnt(0, 0, z_level), gp_Dir(0, 0, 1));BRepAlgoAPI_Section sectionBuilder(face, sectionPlane);sectionBuilder.Build();if (!sectionBuilder.IsDone() || sectionBuilder.Shape().IsNull()) {std::cout << \"信息: 在z=\" << z_level << \"高度处截面为空或创建失败。\" << std::endl;return false;}// 3. 构建轮廓线BRepBuilderAPI_MakeWire wireMaker;TopExp_Explorer edgeExplorer(sectionBuilder.Shape(), TopAbs_EDGE);if (!edgeExplorer.More()) return false;for (; edgeExplorer.More(); edgeExplorer.Next()) {wireMaker.Add(TopoDS::Edge(edgeExplorer.Current()));}if (!wireMaker.IsDone()) return false;TopoDS_Wire sectionProfile = wireMaker.Wire();out_sectionProfile = sectionProfile;// 4. 离散化轮廓线获取采样点BRepAdaptor_CompCurve curveAdaptor(sectionProfile, Standard_True);const int SAMPLE_COUNT = 501;GCPnts_UniformAbscissa discretizer(curveAdaptor, SAMPLE_COUNT);if (!discretizer.IsDone()) return false;std::vector<gp_Pnt> profilePoints;profilePoints.reserve(SAMPLE_COUNT);for (int i = 1; i <= SAMPLE_COUNT; ++i) {profilePoints.push_back(curveAdaptor.Value(discretizer.Parameter(i)));}if (profilePoints.empty()) return false;// 5. 找到前缘点和后缘点的索引和坐标int leadingEdgeIndex = 0;int trailingEdgeIndex = 0;for (size_t i = 1; i < profilePoints.size(); ++i) {if (profilePoints[i].X() < profilePoints[leadingEdgeIndex].X()) leadingEdgeIndex = i;if (profilePoints[i].X() > profilePoints[trailingEdgeIndex].X()) trailingEdgeIndex = i;}const gp_Pnt& leadingEdgePoint = profilePoints[leadingEdgeIndex];const gp_Pnt& trailingEdgePoint = profilePoints[trailingEdgeIndex];// 6. 计算最大弦长out_maxChordLength = trailingEdgePoint.Distance(leadingEdgePoint); // 简化为前后缘点距离，可按需改回暴力法// 7. 计算弦线扭转角double dx = trailingEdgePoint.X() - leadingEdgePoint.X();double dy = trailingEdgePoint.Y() - leadingEdgePoint.Y();out_twistAngle = atan2(dy, dx) * 180.0 / M_PI;// 8. 计算半径 - 鲁棒混合算法try {auto calculateRadius = [&](const gp_Pnt& extremePoint, int extremePointIndex) -> double{double radius = 0.0;TopoDS_Edge bestEdge;double bestParam = 0.0;double minDistance = -1.0;TopExp_Explorer explorer(out_sectionProfile, TopAbs_EDGE);for (; explorer.More(); explorer.Next()) {TopoDS_Edge currentEdge = TopoDS::Edge(explorer.Current());TopLoc_Location loc; Standard_Real first, last;Handle(Geom_Curve) curve = BRep_Tool::Curve(currentEdge, loc, first, last);if (curve.IsNull()) continue;GeomAPI_ProjectPointOnCurve projector(extremePoint, curve, first, last);if (projector.NbPoints() > 0) {double currentDist = projector.LowerDistance();if (minDistance < 0 || currentDist < minDistance) {minDistance = currentDist;bestEdge = currentEdge;bestParam = projector.LowerDistanceParameter();}}}if (!bestEdge.IsNull() && minDistance < 1e-6) {TopLoc_Location loc; Standard_Real first, last;Handle(Geom_Curve) bestCurve = BRep_Tool::Curve(bestEdge, loc, first, last);if (!bestCurve.IsNull()) {gp_Pnt pt; gp_Vec d1, d2;bestCurve->D2(bestParam, pt, d1, d2);double d1_mag_sq = d1.Dot(d1);if (d1_mag_sq > Precision::SquareConfusion()) {double d1_mag = sqrt(d1_mag_sq);double curvature = d1.Crossed(d2).Magnitude() / (d1_mag_sq * d1_mag);if (curvature > Precision::Confusion()) {radius = 1.0 / curvature;}}}}if (radius < Precision::Confusion()) {std::cout << \"信息: 解析法计算半径失败，回退至三点近似法。\" << std::endl;if (profilePoints.size() >= 3) {// 寻找顺序排列的邻近点，并处理环绕情况const int step = 2; // 可调整的步长，避免点过于密集int prevIndex = (extremePointIndex - step + profilePoints.size()) % profilePoints.size();int nextIndex = (extremePointIndex + step) % profilePoints.size();radius = CalculateRadiusFromThreePoints(profilePoints[prevIndex], extremePoint, profilePoints[nextIndex]);}}return radius;};out_leadingEdgeRadius = calculateRadius(leadingEdgePoint, leadingEdgeIndex);out_trailingEdgeRadius = calculateRadius(trailingEdgePoint, trailingEdgeIndex);}catch (const Standard_Failure& e) {std::cerr << \"警告: 半径计算时发生OCCT异常。错误: \" << e.GetMessageString() << std::endl;}// 9. 计算最大厚度gp_Vec maxChordVec(leadingEdgePoint, trailingEdgePoint);gp_Dir maxChordDir = maxChordVec.Normalized();double maxDistUpper = 0.0;double maxDistLower = 0.0;gp_Dir perpDir = maxChordDir.Crossed(gp_Dir(0, 0, 1));for (const auto& pt : profilePoints) {gp_Vec vecFromChord(leadingEdgePoint, pt);double perpDistance = vecFromChord.Dot(perpDir);if (perpDistance > maxDistUpper) {maxDistUpper = perpDistance;}else if (perpDistance < maxDistLower) {maxDistLower = perpDistance;}}out_maxThickness = maxDistUpper - maxDistLower;return true;}catch (const Standard_Failure& e) {std::cerr << \"计算截面特征参数时发生OpenCASCADE错误: \" << e.GetMessageString() << std::endl;return false;}catch (...) {std::cerr << \"计算截面特征参数时发生未知错误。\" << std::endl;return false;}}

从数学上讲，一个完美的“尖点”确实没有半径（或者说曲率是无穷大，半径为0）。我们的算法之所以能算出一个值，并不是因为它“算出”了那个尖点的半径，而是因为它非常聪明地执行了一个**“B计划”**。

您可以把我们的算法想象成一个有两套方案的工程师：

A计划：精密测量法 (解析法)

这是首选方案。它就像一个高精度的激光测量仪。

1. 它会精确地找到截面线上X坐标最小的那个点（我们称之为前缘点）。
2. 然后，它会调用OpenCASCADE的数学库，去问这条曲线：“在你（曲线）经过这个前缘点时，你的精确曲率是多少？”

在尖点处，这个A计划会失败！

为什么呢？因为在一个尖点，曲线是“突然转折”的，它在那个点上并不平滑（数学上叫“不可导”或“二阶导数不存在”）。所以，当你用精密仪器去测量时，仪器会报错，因为它没法在那一点上算出一个有效的曲率。

程序上表现出来就是，计算出的曲率是0，所以半径也是0。

B计划：聪明近似法 (三点画圆的回退方案)

当我们的工程师发现A计划失败了（算出的半径是0），他立刻启动了B计划。这个计划的思路是：

“既然我不能直接测量这个尖点，那我就看看它附近的形状有多‘圆’。”

它是这样做的：

1. 找到三个点：

   • 点2: 就是那个尖尖的前缘点本身。
   • 点1: 在离散化的点集中，找到前缘点前面紧挨着的那个点。
   • 点3: 在离散化的点集中，找到前缘点后面紧挨着的那个点。

2. 用三个点画一个圆：
   现在我们有了三个顺序排列的点。只要这三个点不完全在一条直线上，那么有且仅有一个圆能同时穿过这三个点。

3. 测量这个圆的半径：
   CalculateRadiusFromThreePoints 这个函数的工作，就是计算出这个穿过三个点的虚拟圆的半径。

结论

所以，这个算法最终给出的“前缘半径”：

• 如果前缘是平滑的：它就是A计划算出的、那个点的真实、精确的曲率半径。
• 如果前缘是个尖点：它就是B计划算出的、那个由尖点和它前后两个邻近点构成的虚拟圆的半径。这个半径非常好地近似了那个尖锐区域的“圆润程度”。

简单来说：它算的不是尖点“本身”的半径，而是对尖点附近曲线弯曲程度的一个非常合理、非常稳定的近似值。 这就是为什么这个算法即使面对尖点也能给出有意义的结果。

浏阳新闻