空间圆弧转三阶Bezier曲线的方法

在我们的实际应用中，经常遇到需要把空间圆弧转化为Bezier曲线的方法，下面是详细的matlab代码：

function [P1, P2] = arc_to_bezier_3d(P0, P3, C, direction)
% arc_to_bezier_3d: 将空间圆弧近似为三阶Bezier曲线
% 输入:
%   P0: 起点 [x0, y0, z0]
%   P3: 终点 [x3, y3, z3]
%   C:  圆心 [cx, cy, cz]
%   direction: 'CCW' 或 'CW' (从法向量方向看)
% 输出:
%   P1, P2: 两个控制点坐标

    % 向量从圆心到起点和终点
    v0 = P0 - C;
    v3 = P3 - C;

    % 计算半径（取平均）
    r0 = norm(v0);
    r3 = norm(v3);
    r = (r0 + r3) / 2;

    % 单位化
    v0 = v0 / r0;
    v3 = v3 / r3;

    % 计算圆心角 theta
    cos_theta = dot(v0, v3);
    cos_theta = max(-1, min(1, cos_theta)); % 防止浮点误差
    theta = acos(cos_theta);

    % 计算平面法向量（右手法则）
    n = cross(v0, v3);
    n_norm = norm(n);
    if n_norm < 1e-10
        error('P0, P3, C 三点共线，无法构成圆弧');
    end
    n = n / n_norm;

    % 根据方向调整法向量（CCW: 正方向；CW: 反方向）
    if strcmp(direction, 'CW')
        n = -n;
    end

    % 计算起点和终点的单位切向量（沿圆弧方向）
    % 切向量 = n × 半径向量（右手法则）
    t0 = cross(n, v0);  % 已单位化，因为 n 和 v0 正交且单位
    t3 = cross(n, v3);

    % 根据方向决定切向量方向
    if strcmp(direction, 'CW')
        t0 = -t0;
        t3 = -t3;
    end

    % 计算 alpha
    alpha = (4/3) * r * tan(theta / 4);

    % 构造控制点
    P1 = P0 + alpha * t0;
    P2 = P3 - alpha * t3;

end

函数调用可以用下面的代码：

% 示例：在XOY平面上的四分之一圆弧，逆时针
C = [0, 0, 0];
P0 = [1, 0, 0];
P3 = [0, 1, 0];
direction = 'CCW';

[P1, P2] = arc_to_bezier_3d(P0, P3, C, direction);

% 绘图验证
t = linspace(0, 1, 100)';
B = (1-t).^3 * P0 + 3*(1-t).^2.*t * P1 + 3*(1-t).*t.^2 * P2 + t.^3 * P3;

figure;
plot3(B(:,1), B(:,2), B(:,3), 'b-', 'LineWidth', 2); hold on;
% 绘制真实圆弧
theta_true = linspace(0, pi/2, 100);
X = cos(theta_true);
Y = sin(theta_true);
Z = zeros(size(theta_true));
plot3(X, Y, Z, 'r--', 'LineWidth', 1.5);
plot3(P0(1), P0(2), P0(3), 'go');
plot3(P3(1), P3(2), P3(3), 'go');
plot3([P0(1), P1(1)], [P0(2), P1(2)], [P0(3), P1(3)], 'k:');
plot3([P3(1), P2(1)], [P3(2), P2(2)], [P3(3), P2(3)], 'k:');
legend('Bezier Approximation', 'True Arc', 'Start/End', 'Control Lines');
title('3D Arc to Cubic Bezier Curve');
axis equal; grid on; hold off;