多段圆弧拟合离散点实现切线连续

使用多段圆弧来拟合一个由离散点组成的曲线,并且保证切线连续。也就是说，生成的每一段圆弧之间在连接点处必须有一阶导数连续，也就是切线方向相同。

点集分割

确保每个段的终点是下一段的起点，相邻段共享连接点，避免连接点位于数据点之间。

int totalPoints = envelopePoints.size();
int m = (totalPoints - 1) / numArcs; // 每段间隔数，点数为m+1for (int i = 0; i < numArcs; ++i) {int startIdx = i * m;int endIdx = (i == numArcs - 1) ? (totalPoints - 1) : (startIdx + m);std::vector<gp_Pnt> segmentPoints(envelopePoints.begin() + startIdx, envelopePoints.begin() + endIdx + 1);// 处理该段
}

圆弧拟合约束

问题描述

给定起点 $ S(S_x, S_y) $、终点 $E(E_x, E_y) $ 和起点处的切线方向 $ \mathbf{T}(t_x, t_y) $，求解经过 $ S $ 和 $ E $ 的圆弧，并保证在 $ S $ 处的切线方向为 $ \mathbf{T} $。

推导过程

1. 圆心位置约束

圆心 $ C(O_x, O_y) $ 必须位于 $ S $ 点处与切线方向 $ \mathbf{T} $ 垂直的直线上。设 $ \mathbf{N}(-t_y, t_x) $ 为 $ \mathbf{T} $ 的法线方向，则圆心可表示为：
$$C=S+s\cdot \mathbf{N}$$
即：
$$O_x=S_x-s\cdot t_y$$
$$O_y=S_y+s\cdot t_x$$
其中，$ s $ 是圆心到$ S $点沿法线方向的距离，正负号表示方向。

2. 圆心到终点的距离约束

圆心到终点 $ E $ 的距离必须等于到起点$ S $ 的距离：
$$\sqrt{(E_x-O_x{)}^2+(E_y-O_y{)}^2}=\sqrt{(S_x-O_x{)}^2+(S_y-O_y{)}^2}$$
平方后得到：
$$(E_x-O_x{)}^2+(E_y-O_y{)}^2=(S_x-O_x{)}^2+(S_y-O_y{)}^2$$

3. 代入圆心表达式

将 $ O_x $ 和 $ O_y $ 代入上式：
$$(E_x-(S_x-s\cdot t_y){)}^2+(E_y-(S_y+s\cdot t_x){)}^2=(S_x-(S_x-s\cdot t_y){)}^2+(S_y-(S_y+s\cdot t_x){)}^2$$
简化右边：
$$(s\cdot t_y{)}^2+(s\cdot t_x{)}^2=s^2(t_x^2+t_y^2)=s^2$$
左边展开：
$$(E_x-S_x+s\cdot t_y{)}^2+(E_y-S_y-s\cdot t_x{)}^2$$
$$=(E_x-S_x{)}^2+2s\cdot t_y(E_x-S_x)+s^2{t}_y^2+(E_y-S_y{)}^2-2s\cdot t_x(E_y-S_y)+s^2{t}_x^2$$
$$=(E_x-S_x{)}^2+(E_y-S_y{)}^2+s^2(t_x^2+t_y^2)+2s[t_y(E_x-S_x)-t_x(E_y-S_y)]$$
由于 $ t_x^2 + t_y^2 = 1 $，上式进一步简化为：
$$=(E_x-S_x{)}^2+(E_y-S_y{)}^2+s^2+2s[t_y(E_x-S_x)-t_x(E_y-S_y)]$$

4. 建立方程求解 $ s $

将左边和右边代入距离约束方程：
$$(E_x-S_x{)}^2+(E_y-S_y{)}^2+s^2+2s[t_y(E_x-S_x)-t_x(E_y-S_y)]=s^2$$
化简得到：
$$(E_x-S_x{)}^2+(E_y-S_y{)}^2+2s[t_y(E_x-S_x)-t_x(E_y-S_y)]=0$$
解得：
$$s=\frac{(E_x-S_x{)}^2+(E_y-S_y{)}^2}{-2[t_y(E_x-S_x)-t_x(E_y-S_y)]}$$

5. 计算圆心和半径

代入 $ s$ 的值，得到圆心坐标：
$$O_x=S_x-s\cdot t_y$$
$$O_y=S_y+s\cdot t_x$$
半径为：
$$R=|s|$$

6. 终点切线方向

圆弧在终点 $E $ 处的切线方向由圆心到$ E$ 的径向向量的垂直方向确定。径向向量为：
$$\mathbf{R}=(E_x-O_x,E_y-O_y)$$
切线方向为：
$${\mathbf{T}}_E=(R_y,-R_x)$$
归一化后得到终点处的切线方向。

7. 代码实现

// 计算必须经过起点和终点的圆弧参数
double Sx = arc.start.X();
double Sy = arc.start.Y();
double Ex = arc.end.X();
double Ey = arc.end.Y();
double tx = arc.startTangent.X();
double ty = arc.startTangent.Y();double numerator = (Ex - Sx) * (Ex - Sx) + (Ey - Sy) * (Ey - Sy); // SE_length squared
double denominator = 2 * (ty * (Sx - Ex) - tx * (Sy - Ey));if (std::abs(denominator) < 1e-6) {// 处理直线段情况arc.center = gp_Pnt(0, 0, 0);arc.radius = 0;arc.endTangent = arc.startTangent;arcs.push_back(arc);currentTangent = arc.endTangent;continue;
}double s = numerator / denominator;
double Ox = Sx - s * ty;
double Oy = Sy + s * tx;
arc.center = gp_Pnt(Ox, Oy);
arc.radius = std::abs(s);

切线连续性处理

计算每段圆弧终点处的切线方向，并传递给下一段作为起点切线方向。

// 计算终点切线方向
gp_Vec radialVec(arc.center, arc.end);
if (radialVec.Magnitude() < 1e-6) {arc.endTangent = currentTangent; // 避免除零
} else {radialVec.Normalize();arc.endTangent = gp_Dir(radialVec.Y(), -radialVec.X(), 0);
}
currentTangent = arc.endTangent;

完整代码实现

#include <vector>
#include <gp_Pnt.hxx>
#include <gp_Dir.hxx>
#include <gp_Vec.hxx>
#include <cmath>struct ArcSegment {gp_Pnt start;gp_Dir startTangent;gp_Pnt end;gp_Dir endTangent;gp_Pnt center;double radius;
};std::vector<ArcSegment> fitArcsWithTangentContinuity(const std::vector<gp_Pnt>& envelopePoints, int numArcs) {std::vector<ArcSegment> arcs;if (envelopePoints.size() < 2 || numArcs < 1) return arcs;int totalPoints = envelopePoints.size();int m = (totalPoints - 1) / numArcs; // 每段间隔数，点数为m+1gp_Dir currentTangent;for (int i = 0; i < numArcs; ++i) {int startIdx = i * m;int endIdx = (i == numArcs - 1) ? (totalPoints - 1) : (startIdx + m);if (startIdx >= envelopePoints.size() || endIdx >= envelopePoints.size() || startIdx > endIdx) {break; // Invalid segment, skip}std::vector<gp_Pnt> segmentPoints(envelopePoints.begin() + startIdx, envelopePoints.begin() + endIdx + 1);if (segmentPoints.size() < 2) {continue; // Not enough points to define a segment}ArcSegment arc;arc.start = segmentPoints.front();arc.end = segmentPoints.back();// Determine start tangent directionif (i == 0) {// First segment: use direction from first two pointsgp_Vec initialVec(segmentPoints[0], segmentPoints[1]);if (initialVec.Magnitude() < 1e-6) {currentTangent = gp_Dir(1, 0, 0); // Default direction} else {initialVec.Normalize();currentTangent = gp_Dir(initialVec);}arc.startTangent = currentTangent;} else {arc.startTangent = currentTangent;}double tx = arc.startTangent.X();double ty = arc.startTangent.Y();// Calculate parameters for constrained circle passing through start and end pointsdouble Sx = arc.start.X();double Sy = arc.start.Y();double Ex = arc.end.X();double Ey = arc.end.Y();double numerator = (Ex - Sx) * (Ex - Sx) + (Ey - Sy) * (Ey - Ey);double denominator = 2 * (ty * (Sx - Ex) - tx * (Sy - Ey));if (std::abs(denominator) < 1e-6) {// Handle straight line case (infinite radius)arc.center = gp_Pnt(0, 0, 0); // Not meaningful for straight linearc.radius = 0;arc.endTangent = arc.startTangent;arcs.push_back(arc);currentTangent = arc.endTangent;continue;}double s = numerator / denominator;// Calculate circle parametersdouble Ox = Sx - s * ty;double Oy = Sy + s * tx;arc.center = gp_Pnt(Ox, Oy);arc.radius = std::abs(s);// Calculate end tangent directiongp_Vec radialVec(arc.center, arc.end);if (radialVec.Magnitude() < 1e-6) {arc.endTangent = currentTangent; // Avoid division by zero} else {radialVec.Normalize();arc.endTangent = gp_Dir(radialVec.Y(), -radialVec.X(), 0);}currentTangent = arc.endTangent;arcs.push_back(arc);}return arcs;
}

关键点说明

1. **点集分割：**通过均匀分割确保相邻段共享连接点，避免位置不连续。
2. **圆弧约束拟合：**使用几何约束直接计算必须经过起点和终点的圆弧，保证位置连续。
3. **切线连续性：**通过传递切线方向确保每段圆弧在连接点处切线方向一致。
4. **特殊情况处理：**当分母接近零时，处理为直线段，避免计算错误。

此方案在保证切线连续的同时，简化了计算逻辑，适用于大多数离散点拟合场景。