多段圆弧拟合离散点实现切线连续
使用多段圆弧来拟合一个由离散点组成的曲线,并且保证切线连续。也就是说，生成的每一段圆弧之间在连接点处必须有一阶导数连续，也就是切线方向相同。
点集分割
确保每个段的终点是下一段的起点，相邻段共享连接点，避免连接点位于数据点之间。
int totalPoints = envelopePoints.size();
int m = (totalPoints - 1) / numArcs; // 每段间隔数，点数为m+1
for (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+1
gp_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 direction
if (i == 0) {
// First segment: use direction from first two points
gp_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 points
double 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 line
arc.radius = 0;
arc.endTangent = arc.startTangent;
arcs.push_back(arc);
currentTangent = arc.endTangent;
continue;
}
double s = numerator / denominator;
// Calculate circle parameters
double Ox = Sx - s * ty;
double Oy = Sy + s * tx;
arc.center = gp_Pnt(Ox, Oy);
arc.radius = std::abs(s);
// Calculate end tangent direction
gp_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;
}
关键点说明
- **点集分割：**通过均匀分割确保相邻段共享连接点，避免位置不连续。
- **圆弧约束拟合：**使用几何约束直接计算必须经过起点和终点的圆弧，保证位置连续。
- **切线连续性：**通过传递切线方向确保每段圆弧在连接点处切线方向一致。
- **特殊情况处理：**当分母接近零时，处理为直线段，避免计算错误。
此方案在保证切线连续的同时，简化了计算逻辑，适用于大多数离散点拟合场景。