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计划算出的、那个由尖点和它前后两个邻近点构成的虚拟圆的半径。这个半径非常好地近似了那个尖锐区域的“圆润程度”。

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