文章目录
一、问题描述
二、推导步骤
三、MATLAB代码

一、问题描述
给定空间不共线的三个点 A , B , C ，推导空间有向圆弧路径 A B C关于路径标量 s ( s ∈ [ 0 , 1 ] ) 的参数方程。

三、MATLAB代码

%{
Function: solve_spatial_arc_params
Description: 求三点空间圆弧路径参数
Input: 空间三个点startPos, midPos, endPos
Output: 空间圆弧圆心arcCenter, 半径radius, 圆心角arcCentralAngle, 
        旋转变换矩阵rotateMatrix, 求解状态status(1表示成功，0表示失败)
Author: Marc Pony(marc_pony@163.com)
%}
function [arcCenter, radius, arcCentralAngle, rotateMatrix, status] = solve_spatial_arc_params(startPos, midPos, endPos)
x1 = startPos(1);
y1 = startPos(2);
z1 = startPos(3);
x2 = midPos(1);
y2 = midPos(2);
z2 = midPos(3);
x3 = endPos(1);
y3 = endPos(2);
z3 = endPos(3);

arcCenter = zeros(3, 1);
radius = 0.0;
arcCentralAngle = 0.0;
rotateMatrix = zeros(3, 3);
status = 1;

if abs((y1 - y2) * (z2 - z3) - (y2 - y3) * (z1 - z2)) < 1.0e-8 ...
   && abs((x2 - x3) * (z1 - z2) - (x1 - x2) * (z2 - z3)) < 1.0e-8 ...
   && abs((x1 - x2) * (y2 - y3) - (x2 - x3) * (y1 - y2)) < 1.0e-8
    status = 0;
    return;
else
    x4 = 0.5 * (x1 + x2);
    y4 = 0.5 * (y1 + y2);
    z4 = 0.5 * (z1 + z2);
    
    x5 = 0.5 * (x2 + x3);
    y5 = 0.5 * (y2 + y3);
    z5 = 0.5 * (z2 + z3);
    
    a11 = x2 - x1;
    a12 = y2 - y1;
    a13 = z2 - z1;
    b1 = x4 * a11 + y4 * a12 + z4 * a13;
    
    a21 = x3 - x2;
    a22 = y3 - y2;
    a23 = z3 - z2;
    b2 = x5 * a21 + y5 * a22 + z5 * a23;
    
    a31 = (y1 - y2) * (z2 - z3) - (y2 - y3) * (z1 - z2);
    a32 = (x2 - x3) * (z1 - z2) - (x1 - x2) * (z2 - z3);
    a33 = (x1 - x2) * (y2 - y3) - (x2 - x3) * (y1 - y2);
    b3 = x3 * a31 + y3 * a32 + z3 * a33;
    
    temp = a11 * (a22 * a33 - a23 * a32) + a12 * (a23 * a31 - a21 * a33) + a13 * (a21 * a32 - a22 * a31);
    x0 = ((a12 * a23 - a13 * a22) * b3 + (a13 * a32 - a12 * a33) * b2 + (a22 * a33 - a23 * a32) * b1) / temp;
    y0 = -((a11 * a23 - a13 * a21) * b3 + (a13 * a31 - a11 * a33) * b2 + (a21 * a33 - a23 * a31) * b1) / temp;
    z0 = ((a11 * a22 - a12 * a21) * b3 + (a12 * a31 - a11 * a32) * b2 + (a21 * a32 - a22 * a31) * b1) / temp;
    radius = sqrt((x0 - x1)^2 + (y0 - y1)^2 + (z0 - z1)^2);
    
    r11 = (x1 - x0) / radius;
    r21 = (y1 - y0) / radius;
    r31 = (z1 - z0) / radius;
    
    temp1 = (y1 - y2) * (z2 - z3) - (y2 - y3) * (z1 - z2);
    temp2 = (x2 - x3) * (z1 - z2) - (x1 - x2) * (z2 - z3);
    temp3 = (x1 - x2) * (y2 - y3) - (x2 - x3) * (y1 - y2);
    vectorZLength = sqrt(temp1 * temp1 + temp2 * temp2 + temp3 * temp3);
    r13 = temp1 / vectorZLength;
    r23 = temp2 / vectorZLength;
    r33 = temp3 / vectorZLength;
    
    r12 = r23 * r31 - r21 * r33;
    r22 = r11 * r33 - r13 * r31;
    r32 = r13 * r21 - r11 * r23;
    
    arcCentralAngle = atan2((x3 - x0) * r12 + (y3 - y0) * r22 + (z3 - z0) * r32, ...
                      (x3 - x0) * r11 + (y3 - y0) * r21 + (z3 - z0) * r31);
    if  arcCentralAngle < 0.0
        arcCentralAngle = arcCentralAngle + 2.0 * pi;
    end
    arcCenter = [x0; y0; z0];
    rotateMatrix = [r11, r12, r13; r21, r22, r23; r31, r32, r33];
end
end

clc
clear
close all

startPos = [12, 7, 20]'; %圆弧起点
midPos = [0, 12, 10]';   %圆弧中间点
endPos = [8, 0, 10]';    %圆弧终点

[arcCenter, radius, arcCentralAngle, rotateMatrix, status] = solve_spatial_arc_params(startPos, midPos, endPos);

count = 1000;
s = linspace(0, 1, count)';
xs = radius * cos(s * arcCentralAngle);
ys = radius * sin(s * arcCentralAngle);
x = zeros(count, 1);
y = zeros(count, 1);
z = zeros(count, 1);
for i = 1 : count
    x(i) = rotateMatrix(1, 1) * xs(i) + rotateMatrix(1, 2) * ys(i) + arcCenter(1);
    y(i) = rotateMatrix(2, 1) * xs(i) + rotateMatrix(2, 2) * ys(i) + arcCenter(2);
    z(i) = rotateMatrix(3, 1) * xs(i) + rotateMatrix(3, 2) * ys(i) + arcCenter(3);
end

figure('color', 'w')
plot3([startPos(1), midPos(1), endPos(1)], [startPos(2), midPos(2), endPos(2)], [startPos(3), midPos(3), endPos(3)], 'o')
hold on
plot3(x, y, z, 'r')
xlabel('x')
ylabel('y')
zlabel('z')
axis equal