宣城网站推广全国疫情又严重了
使用多段圆弧来拟合一个由离散点组成的曲线,并且保证切线连续。也就是说，生成的每一段圆弧之间在连接点处必须有一阶导数连续，也就是切线方向相同。
点集分割
确保每个段的终点是下一段的起点，相邻段共享连接点，避免连接点位于数据点之间。
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 ⋅ N C = S + s \cdot \mathbf{N} C=S+s⋅N
即：
O x = S x − s ⋅ t y O_x = S_x - s \cdot t_y Ox=Sx−s⋅ty
O y = S y + s ⋅ t x O_y = S_y + s \cdot t_x Oy=Sy+s⋅tx
其中，$ s $ 是圆心到$ S $点沿法线方向的距离，正负号表示方向。
2. 圆心到终点的距离约束
圆心到终点 $ E $ 的距离必须等于到起点$ S $ 的距离：
( E x − O x ) 2 + ( E y − O y ) 2 = ( S x − O x ) 2 + ( S y − O y ) 2 \sqrt{(E_x - O_x)^2 + (E_y - O_y)^2} = \sqrt{(S_x - O_x)^2 + (S_y - O_y)^2} (Ex−Ox)2+(Ey−Oy)2=(Sx−Ox)2+(Sy−Oy)2
平方后得到：
( E x − O x ) 2 + ( E y − O y ) 2 = ( 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 (Ex−Ox)2+(Ey−Oy)2=(Sx−Ox)2+(Sy−Oy)2
3. 代入圆心表达式
将 $ O_x $ 和 $ O_y $ 代入上式：
( E x − ( S x − s ⋅ t y ) ) 2 + ( E y − ( S y + s ⋅ t x ) ) 2 = ( S x − ( S x − s ⋅ t y ) ) 2 + ( S y − ( S y + s ⋅ t x ) ) 2 (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 (Ex−(Sx−s⋅ty))2+(Ey−(Sy+s⋅tx))2=(Sx−(Sx−s⋅ty))2+(Sy−(Sy+s⋅tx))2
简化右边：
( s ⋅ t y ) 2 + ( s ⋅ t x ) 2 = s 2 ( t x 2 + t y 2 ) = s 2 (s \cdot t_y)^2 + (s \cdot t_x)^2 = s^2 (t_x^2 + t_y^2) = s^2 (s⋅ty)2+(s⋅tx)2=s2(tx2+ty2)=s2
左边展开：
( E x − S x + s ⋅ t y ) 2 + ( E y − S y − s ⋅ t x ) 2 (E_x - S_x + s \cdot t_y)^2 + (E_y - S_y - s \cdot t_x)^2 (Ex−Sx+s⋅ty)2+(Ey−Sy−s⋅tx)2
= ( E x − S x ) 2 + 2 s ⋅ t y ( E x − S x ) + s 2 t y 2 + ( E y − S y ) 2 − 2 s ⋅ t x ( E y − S y ) + s 2 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 =(Ex−Sx)2+2s⋅ty(Ex−Sx)+s2ty2+(Ey−Sy)2−2s⋅tx(Ey−Sy)+s2tx2
= ( E x − S x ) 2 + ( E y − S y ) 2 + s 2 ( t x 2 + t y 2 ) + 2 s [ t y ( E x − S x ) − t x ( E y − S y ) ] = (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)] =(Ex−Sx)2+(Ey−Sy)2+s2(tx2+ty2)+2s[ty(Ex−Sx)−tx(Ey−Sy)]
由于 $ t_x^2 + t_y^2 = 1 $，上式进一步简化为：
= ( E x − S x ) 2 + ( E y − S y ) 2 + s 2 + 2 s [ t y ( E x − S x ) − t x ( E y − S y ) ] = (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)] =(Ex−Sx)2+(Ey−Sy)2+s2+2s[ty(Ex−Sx)−tx(Ey−Sy)]
4. 建立方程求解 $ s $
将左边和右边代入距离约束方程：
( E x − S x ) 2 + ( E y − S y ) 2 + s 2 + 2 s [ t y ( E x − S x ) − t x ( E y − S y ) ] = s 2 (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 (Ex−Sx)2+(Ey−Sy)2+s2+2s[ty(Ex−Sx)−tx(Ey−Sy)]=s2
化简得到：
( E x − S x ) 2 + ( E y − S y ) 2 + 2 s [ t y ( E x − S x ) − t x ( E y − S y ) ] = 0 (E_x - S_x)^2 + (E_y - S_y)^2 + 2s [t_y (E_x - S_x) - t_x (E_y - S_y)] = 0 (Ex−Sx)2+(Ey−Sy)2+2s[ty(Ex−Sx)−tx(Ey−Sy)]=0
解得：
s = ( E x − S x ) 2 + ( E y − S y ) 2 − 2 [ t y ( E x − S x ) − t x ( E y − S y ) ] 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)]} s=−2[ty(Ex−Sx)−tx(Ey−Sy)](Ex−Sx)2+(Ey−Sy)2
5. 计算圆心和半径
代入 $ s$ 的值，得到圆心坐标：
O x = S x − s ⋅ t y O_x = S_x - s \cdot t_y Ox=Sx−s⋅ty
O y = S y + s ⋅ t x O_y = S_y + s \cdot t_x Oy=Sy+s⋅tx
半径为：
R = ∣ s ∣ R = |s| R=∣s∣
6. 终点切线方向
圆弧在终点 $E $ 处的切线方向由圆心到$ E$ 的径向向量的垂直方向确定。径向向量为：
R = ( E x − O x , E y − O y ) \mathbf{R} = (E_x - O_x, E_y - O_y) R=(Ex−Ox,Ey−Oy)
切线方向为：
T E = ( R y , − R x ) \mathbf{T}_E = (R_y, -R_x) TE=(Ry,−Rx)
归一化后得到终点处的切线方向。
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;
}
关键点说明
- **点集分割：**通过均匀分割确保相邻段共享连接点，避免位置不连续。
- **圆弧约束拟合：**使用几何约束直接计算必须经过起点和终点的圆弧，保证位置连续。
- **切线连续性：**通过传递切线方向确保每段圆弧在连接点处切线方向一致。
- **特殊情况处理：**当分母接近零时，处理为直线段，避免计算错误。
此方案在保证切线连续的同时，简化了计算逻辑，适用于大多数离散点拟合场景。