如何使用矩形放样以在 MATLAB 中创建灵活的 3D 封闭管道？

Question

我有一系列矩形并且知道它们各自 4 个角的确切位置。我想绘制它们并放样通过它们中的每一个以创建类似于矩形横截面的 3D 管道的东西。这些点也不应限于遵循直轴。它应该灵活地采取偏差。我也想把两端缝起来。

我在你的网站 "How to loft with ellipses to create a 3d hollow pipe in MATLAB or Python?" 下看到过类似的关于放样的问题。答案给我留下了深刻的印象，但对于椭圆和圆圈来说却令人遗憾。我试图让它与矩形一起工作，但无法弄清楚所需的逻辑。我也尝试过将所有东西拼凑在一起，但这会导致产生尖锐的边缘，这是我不想要的。我的代码看起来像这样：

A = importdata(filename);

[size_A, ~] = size(A.data);
axis vis3d;

for i=1:12:size_A-12

    X1 = A.data(i+1); X2 = A.data(i+4); X3 = A.data (i+7); X4= A.data (i+10);
    Y1 = A.data(i+2); Y2 = A.data(i+5); Y3 = A.data(i+8); Y4 = A.data(i+11);
    Z1 = A.data(i+3); Z2 = A.data(i+6); Z3 = A.data(i+9); Z4 = A.data(i+12);

    X= [X1;X2;X3;X4]; 
    Y= [Y1;Y2;Y3;Y4]; 
    Z= [Z1;Z2;Z3;Z4]; 


    plot3(X, Y, Z)
    patch(X, Y, Z, 'g'); %% for the particular planes


    if(i>1) %% for the patching between two planes

        A1= [ X1 X1 X2 X4; a1 X4 a2 X3; a2 a4 a3 a3; X2 a1 X3 a4];
        B1= [ Y1 Y1 Y2 Y4; b1 Y4 b2 Y3; b2 b4 b3 b3; Y2 b1 Y3 b4];
        C1= [ Z1 Z1 Z2 Z4; c1 Z4 c2 Z3; c2 c4 c3 c3; Z2 c1 Z3 c4];

        plot3(A1, B1, C1)
        patch(A1, B1, C1, 'g');

    end

    a1=X1; a2=X2; a3=X3; a4=X4;
    b1=Y1; b2=Y2; b3=Y3; b4=Y4; 
    c1=Z1; c2=Z2; c3=Z3; c4=Z4;

    figure(1)
    grid on
    axis equal
    hold on
    xlabel('x');
    ylabel('y');
    zlabel('z');

end

最终结果应该是一个非常弯曲的矩形横截面管道。不应有任何尖角。

PS: MATLAB文件正在导入记事本.txt文件，输入坐标如下图-

Number_of_sections= 7

X_coordinate= 60.00
Y_coordinate= 13.00
Z_coordinate= -7.50
X_coordinate= 60.00
Y_coordinate= -13.00
Z_coordinate= -7.50
X_coordinate= 60.00
Y_coordinate= -13.00
Z_coordinate= -12.50
X_coordinate= 60.00
Y_coordinate= 13.00
Z_coordinate= -12.50

X_coordinate= 95.00
Y_coordinate= 13.00
Z_coordinate= -7.50
X_coordinate= 95.00
Y_coordinate= -13.00
Z_coordinate= -7.50
X_coordinate= 95.00
Y_coordinate= -13.00
Z_coordinate= -12.50
X_coordinate= 95.00
Y_coordinate= 13.00
Z_coordinate= -12.50

X_coordinate= 95.50
Y_coordinate= 13.50
Z_coordinate= -7.50
X_coordinate= 95.50
Y_coordinate= -13.50
Z_coordinate= -7.50
X_coordinate= 95.50
Y_coordinate= -13.50
Z_coordinate= -12.50
X_coordinate= 95.50
Y_coordinate= 13.50
Z_coordinate= -12.50

X_coordinate= 96.00
Y_coordinate= 14.00
Z_coordinate= -7.50
X_coordinate= 96.00
Y_coordinate= -14.00
Z_coordinate= -7.50
X_coordinate= 96.00
Y_coordinate= -14.00
Z_coordinate= -12.50
X_coordinate= 96.00
Y_coordinate= 14.00
Z_coordinate= -12.50

X_coordinate= 96.50
Y_coordinate= 14.50
Z_coordinate= -7.50
X_coordinate= 96.50
Y_coordinate= -14.50
Z_coordinate= -7.50
X_coordinate= 96.50
Y_coordinate= -14.50
Z_coordinate= -12.50
X_coordinate= 96.50
Y_coordinate= 14.50
Z_coordinate= -12.50

X_coordinate= 97.00
Y_coordinate= 15.00
Z_coordinate= -7.50
X_coordinate= 97.00
Y_coordinate= -15.00
Z_coordinate= -7.50
X_coordinate= 97.00
Y_coordinate= -15.00
Z_coordinate= -12.50
X_coordinate= 97.00
Y_coordinate= 15.00
Z_coordinate= -12.50

X_coordinate= 99.00
Y_coordinate= 15.00
Z_coordinate= -7.50
X_coordinate= 99.00
Y_coordinate= -15.00
Z_coordinate= -7.50
X_coordinate= 99.00
Y_coordinate= -15.00
Z_coordinate= -12.50
X_coordinate= 99.00
Y_coordinate= 15.00
Z_coordinate= -12.50

X_coordinate= 99.250
Y_coordinate= 14.7500
Z_coordinate= -7.50
X_coordinate= 99.250
Y_coordinate= -15.2500
Z_coordinate= -7.50
X_coordinate= 99.250
Y_coordinate= -15.2500
Z_coordinate= -12.50
X_coordinate= 99.250
Y_coordinate= 14.7500
Z_coordinate= -12.50

X_coordinate= 99.50
Y_coordinate= 14.50
Z_coordinate= -7.50
X_coordinate= 99.50
Y_coordinate= -15.500
Z_coordinate= -7.50
X_coordinate= 99.50
Y_coordinate= -15.500
Z_coordinate= -12.50
X_coordinate= 99.50
Y_coordinate= 14.50
Z_coordinate= -12.50

X_coordinate= 100.0
Y_coordinate= 14.00
Z_coordinate= -7.50
X_coordinate= 100.0
Y_coordinate= -15.750
Z_coordinate= -7.50
X_coordinate= 100.0
Y_coordinate= -15.750
Z_coordinate= -12.50
X_coordinate= 100.0
Y_coordinate= 14.00
Z_coordinate= -12.50

X_coordinate= 110.0
Y_coordinate= 14.00
Z_coordinate= -7.50
X_coordinate= 110.0
Y_coordinate= -15.750
Z_coordinate= -7.50
X_coordinate= 110.0
Y_coordinate= -15.750
Z_coordinate= -12.50
X_coordinate= 110.0
Y_coordinate= 14.00
Z_coordinate= -12.50



X_coordinate= 110.0
Y_coordinate= 14.00
Z_coordinate= -7.50
X_coordinate= 110.0
Y_coordinate= -15.750
Z_coordinate= -7.50
X_coordinate= 120.0
Y_coordinate= -15.750
Z_coordinate= -12.50
X_coordinate= 120.0
Y_coordinate= 14.00
Z_coordinate= -12.50

X_coordinate= 110.0
Y_coordinate= 14.00
Z_coordinate= -5.50
X_coordinate= 110.0
Y_coordinate= -15.750
Z_coordinate= -5.50
X_coordinate= 120.0
Y_coordinate= -15.750
Z_coordinate= -7.50
X_coordinate= 120.0
Y_coordinate= 14.00
Z_coordinate= -7.50

X_coordinate= 80.0
Y_coordinate= 14.00
Z_coordinate= -5.50
X_coordinate= 80.0
Y_coordinate= -15.750
Z_coordinate= -5.50
X_coordinate= 80.0
Y_coordinate= -15.750
Z_coordinate= -7.50
X_coordinate= 80.0
Y_coordinate= 14.00
Z_coordinate= -7.50

X_coordinate= 70.0
Y_coordinate= 14.00
Z_coordinate= -5.50
X_coordinate= 70.0
Y_coordinate= -15.750
Z_coordinate= -5.50
X_coordinate= 80.0
Y_coordinate= -15.750
Z_coordinate= -5.50
X_coordinate= 80.0
Y_coordinate= 14.00
Z_coordinate= -5.50

X_coordinate= 70.0
Y_coordinate= 14.00
Z_coordinate= -2.50
X_coordinate= 70.0
Y_coordinate= -15.750
Z_coordinate= -2.50
X_coordinate= 80.0
Y_coordinate= -15.750
Z_coordinate= -2.50
X_coordinate= 80.0
Y_coordinate= 14.00
Z_coordinate= -2.50

方向图 2 的所需变化

方向图 3 中所需变化的更详细的表示

Answer 1

通过数据文件中的点进行插值并不容易。我建议使用 design B 样条曲线在平面 1 和平面 2 以及平面 6 和平面 7 之间插入点。平面 2 和平面 6、平面 7 和平面 10 之间的空间看起来非常线性：

clear
filename = 'datafile.txt';
A = importdata(filename);

vertices = A.data(2:end);
vertices = reshape(vertices, 12, [])';

vx = vertices(:, 1:3:10);
vy = vertices(:,2:3:11);
vz = vertices(:,3:3:12);
figure; patch(vx',vy',vz',1); axis equal;

没有简单的方法可以进行此类插值，因为您要确保沿曲线至少有 C1 连续性以避免出现任何尖锐边缘。 B 样条曲线在这里可能很有用，但如果您不熟悉它，您将很难对其进行编程。幸运的是，我在做一个需要曲面和曲线拟合的项目，并且手头有代码：

function [x,y,z] = bspline(u, ctrlp, k, knots)    
U = bspbasis(u, numel(ctrlp(:,1)), k, knots);

x=U*ctrlp(:,1);
y=U*ctrlp(:,2);
z=U*ctrlp(:,3);
end

function U = bspbasis(u, nctrlp, K, knots)
nu = numel(u);
umax = max(u);
index = 1:nctrlp;

% preallocating variables
U = zeros(nu,nctrlp);
N = zeros(nctrlp+1,K);

% Calculate the denominators for basis functions (k>2). -may be useful when
% the size of point data is substantial, so the calculation is not repeated.
d1 = zeros(nctrlp,K);
d2 = d1;
for m=2:K
    d1(:,m) = knots(index+m-1) - knots(index);
    d2(:,m) = knots(index+m) - knots(index+1);
end

knots = knots(:); 
knots1 = knots(1:nctrlp+1); 
knots2 = knots(2:nctrlp+2);
knotSize = size(knots1);
knotc1 = knots(nctrlp+1); 
knotc2 = knots(nctrlp+2);

for ui = 1:nu
    % k = 1
    u_ = u(ui)*ones(knotSize);
    NA = u_>=knots1 & u_<knots2;
    if u(ui) == umax && knotc1 == umax && knotc2 == umax
        NA(1:nctrlp+1) = 0;
        NA(end-1,1) = 1;
    end
    N(:,1) = NA;
    % k > 2
    for k = 2:K
        p1 = (u(ui) - knots(index)) ./ d1(:,k) .* N(index,k-1);
        p1(isnan(p1)) = 0;
        p2 = (knots(index+k) - u(ui)) ./ d2(:,k) .* N(index+1,k-1);
        p2(isnan(p2)) = 0;
        N(index,k) = p1 + p2;
    end
    U(ui,:)=N(1:end-1,k);
end
end

如果想看懂上面的代码，可以通读B-spline Wikipedia. There are also lots of tutorials on Youtube and interactive tools such as this one。

下面的代码将三阶 B 样条拟合到四组角顶点中的每一组。

u = linspace(0,1,500);
k = 3;

for i = 1:4
    ishift = (i-1) * 3 + 1;
    p = vertices(:,ishift:ishift+2);
        
    ctrlp = [p(1,:); [0 0 0]; p(2:3,:); p(5:6,:); zeros(2,3); p(7,:); p(8,:)];
    ctrlp(2,:) = 5*ctrlp(3,:) - 4*ctrlp(4,:);
    ctrlp(7,:) = 2*ctrlp(6,:) - ctrlp(5,:); 
    ctrlp(8,:) = 4*ctrlp(9,:) - 3*ctrlp(10,:);
    knots = getknots(ctrlp, k);
    [x_,y_,z_] = bspline(u, ctrlp, k, knots);
    x(:,i) = [p(:,1); x_];
    y(:,i) = [p(:,2); y_];
    z(:,i) = [p(:,3); z_];
    
    [~,I] = sort(x(:,i));
    x(:,i) = x(I,i);
    y(:,i) = y(I,i);
    z(:,i) = z(I,i);
end
c = repmat(1:numel(x)/4,4,1)';
xx=[x;flipud(x(:,[2,3,4,1]))];
yy=[y;flipud(y(:,[2,3,4,1]))];
zz=[z;flipud(z(:,[2,3,4,1]))];
cc=[c;flipud(c)];
figure; patch(xx,yy,zz,cc);

结果是一个平滑的放样曲面：

我想向您解释代码是如何工作的，但是在我苦苦思索了一个小时之后，我放弃了...相反，我在下面提供了一些关键点的摘要。

在总结中，使用了以下符号：C1、C2、...、C10 是 B 样条曲线的控制点，V1、V2、...、V10 是用于计算B 样条曲线。下图显示了用于计算第一条 B 样条曲线的控制点和顶点。

u 中值的数量决定了最终曲线上的点数。
三阶 B 样条曲线就足够了。
曲线必须通过 V1, V2, ..., V10。
要通过V1和V10，B样条曲线必须是clamped类型。因此，C1 = V1 和 C10 = V10.
要通过 V2，C2 必须在通过 V2 和 V3 的线上，并且 C3 必须等于 V2。因此，C2 - V2 = d*(V2-V3) 和 C3 = V2.
要通过V3、V4、V5，必须满足以下条件：C3 = V2; C4 = V3; C5 = V5; C6 = V6.
要通过V6，C7必须在通过V5和V6的线上，并且C6等于V6。因此，C7-V6 = d*(V6-V5).
要通过V7，C8必须在通过V7和V8的线上，并且C9必须等于V7。因此，C8-V7 = d*(V7-V8); C9 = V7;.
由于C10 = V10和C9 = V7，曲线通过V8和V9。
最后，节点可以是统一的，也可以通过控制点之间的弦长来估计。

更新：`getknots`

getknots函数：

function knots = getknots(ctrlp, k)
d = sqrt(sum(diff(ctrlp).^2, 2));
ds = cumsum(d)./sum(d);
knots = [zeros(k,1); ds(k-1:end); ones(k,1)];
end

bspline 方法的全部代码：

clear
filename = 'datafile.txt';
A = importdata(filename);

vertices = A.data(2:end);
vertices = reshape(vertices, 12, [])';

u = linspace(0,1,500);
k = 3;

for i = 1:4
    ishift = (i-1) * 3 + 1;
    p = vertices(:,ishift:ishift+2);

    ctrlp = [p(1,:); [0 0 0]; p(2:3,:); p(5:6,:); zeros(2,3); p(7,:); p(8,:)];
    ctrlp(2,:) = 5*ctrlp(3,:) - 4*ctrlp(4,:);
    ctrlp(7,:) = 2*ctrlp(6,:) - ctrlp(5,:); 
    ctrlp(8,:) = 4*ctrlp(9,:) - 3*ctrlp(10,:);
    knots = getknots(ctrlp, k);
    [x_,y_,z_] = bspline(u, ctrlp, k, knots);
    x(:,i) = [p(:,1); x_];
    y(:,i) = [p(:,2); y_];
    z(:,i) = [p(:,3); z_];

    [~,I] = sort(x(:,i));
    x(:,i) = x(I,i);
    y(:,i) = y(I,i);
    z(:,i) = z(I,i);
end
c = repmat(1:numel(x)/4,4,1)';
xx=[x;flipud(x(:,[2,3,4,1]))];
yy=[y;flipud(y(:,[2,3,4,1]))];
zz=[z;flipud(z(:,[2,3,4,1]))];
cc=[c;flipud(c)];
figure; patch(xx,yy,zz,cc);

function [x,y,z] = bspline(u, ctrlp, k, knots)    
U = bspbasis(u, numel(ctrlp(:,1)), k, knots);

x=U*ctrlp(:,1);
y=U*ctrlp(:,2);
z=U*ctrlp(:,3);
end

function U = bspbasis(u, nctrlp, K, knots)
nu = numel(u);
umax = max(u);
index = 1:nctrlp;

% preallocating variables
U = zeros(nu,nctrlp);
N = zeros(nctrlp+1,K);

% Calculate the denominators for basis functions (k>2). -may be useful when
% the size of point data is substantial, so the calculation is not repeated.
d1 = zeros(nctrlp,K);
d2 = d1;
for m=2:K
    d1(:,m) = knots(index+m-1) - knots(index);
    d2(:,m) = knots(index+m) - knots(index+1);
end

knots = knots(:); 
knots1 = knots(1:nctrlp+1); 
knots2 = knots(2:nctrlp+2);
knotSize = size(knots1);
knotc1 = knots(nctrlp+1); 
knotc2 = knots(nctrlp+2);

for ui = 1:nu
    % k = 1
    u_ = u(ui)*ones(knotSize);
    NA = u_>=knots1 & u_<knots2;
    if u(ui) == umax && knotc1 == umax && knotc2 == umax
        NA(1:nctrlp+1) = 0;
        NA(end-1,1) = 1;
    end
    N(:,1) = NA;
    % k > 2
    for k = 2:K
        p1 = (u(ui) - knots(index)) ./ d1(:,k) .* N(index,k-1);
        p1(isnan(p1)) = 0;
        p2 = (knots(index+k) - u(ui)) ./ d2(:,k) .* N(index+1,k-1);
        p2(isnan(p2)) = 0;
        N(index,k) = p1 + p2;
    end
    U(ui,:)=N(1:end-1,k);
end
end

function knots = getknots(ctrlp, k)
d = sqrt(sum(diff(ctrlp).^2, 2));
ds = cumsum(d)./sum(d);
knots = [zeros(k,1); ds(k-1:end); ones(k,1)];
end

Answer 2

这个答案最初发布在我的第一个答案中，但是第一个答案太拥挤而难以追踪。因此，我现在将这个答案分成第二个答案。这也是解决问题的另一种不同方法。

如果你有一组密度较好的点，可以使用MATLAB内置的插值函数interp1。

首先，您需要对点进行参数化。这里的弦长足够好：

clear
filename = 'datafile.txt';
A = importdata(filename);

vertices = A.data(2:end);
vertices = reshape(vertices, 12, [])';

for i = 1:4
    % take a set of vertices that will form an edge of the lofted body.
    ishift = (i-1) * 3 + 1;
    p = vertices(:,ishift:ishift+2);

    %calculate distance between two neighbouring points.
    dp = sqrt(sum(diff(p).^2,2)); 
    u = [0; cumsum(dp)/sum(dp)];

    % do more later
end

然后生成一组新的u值用于计算拟合点：

unew = unique([u; linspace(0,1,500)']);

然后用interp1插入点：

pnew = interp1(u, p, unew, 'spline');
figure; hold on;
plot3(pnew(:,1), pnew(:,2), pnew(:,3));
plot3(p(:,1), p(:,2), p(:,3),'o');
axis equal;

总结 for 循环现在是：

for i = 1:4
    ishift = (i-1) * 3 + 1;
    p = vertices(:,ishift:ishift+2);

    dp = sqrt(sum(diff(p).^2,2));
    u = [0; cumsum(dp)/sum(dp)];
    unew = unique([u; linspace(0,1,500)']);
    pnew = interp1(u, p, unew, 'spline');
    x(:,i) = pnew(:,1);
    y(:,i) = pnew(:,2);
    z(:,i) = pnew(:,3);
end

但是，由于每条边只有 10 个点，您将得到如下结果：

取其中一条边曲线：

V1和V2之间的巨大曲率显然是不可取的。但是，由于 V1 和 V2 之间的变化只需要是线性的，您可以在它们之间人为地添加点。比如加两点：

p = vertices(:,ishift:ishift+2);
insertedP = (p(2,:) - p(1,:)).*(0.33; 0.66);
insertedP = insertedP + repmat(p(1,:),size(insertedP,1),1);
p = [p(1,:); insertedP; p(2:end,:)];

再做同样的计算，得到：

如果再添加：

insertedP = (p(2,:) - p(1,:)).*(0.02:0.02:0.98)';

但这并不完美，因为你会在 V2 处看到可见的振荡，这是由于连续性而不可避免的：

但是，可以通过移除或添加靠近 V2 的点来控制：

insertedP = (p(2,:) - p(1,:)).*(0.02:0.02:0.68)';

请注意，您也可以使用其他答案中发布的 B 样条方法获得这种振荡，但通过将第二个控制点移近或移远至 V2，它更可控和可预测。

用interp1方法生成的放样体是这样的：

总而言之，下面给出了重现上图的完整代码：

clear
filename = 'datafile.txt';
A = importdata(filename);

vertices = A.data(2:end);
vertices = reshape(vertices, 12, [])';

for i = 1:4
    ishift = (i-1) * 3 + 1;
    p = vertices(:,ishift:ishift+2);
    insertedP = (p(2,:) - p(1,:)).*(0.02:0.02:0.98)';
    insertedP = insertedP + repmat(p(1,:),size(insertedP,1),1);
    p = [p(1,:); insertedP; p(2:end,:)];

    dp = sqrt(sum(diff(p).^2,2));
    u = [0; cumsum(dp)/sum(dp)];
    unew = unique([u; linspace(0,1,500)']);
    pnew = interp1(u, p, unew, 'spline');
    x(:,i) = pnew(:,1);
    y(:,i) = pnew(:,2);
    z(:,i) = pnew(:,3);
end
c = repmat(1:numel(x)/4,4,1)';
xx=[x;flipud(x(:,[2,3,4,1]))];
yy=[y;flipud(y(:,[2,3,4,1]))];
zz=[z;flipud(z(:,[2,3,4,1]))];
cc=[c;flipud(c)];
figure; patch(xx,yy,zz,cc);

如果你有更多的点来约束插值，你可以删除涉及insertP:

的行

insertedP = (p(2,:) - p(1,:)).*(0.02:0.02:0.98)';
insertedP = insertedP + repmat(p(1,:),size(insertedP,1),1);
p = [p(1,:); insertedP; p(2:end,:)];

如何使用矩形放样以在 MATLAB 中创建灵活的 3D 封闭管道？

How to loft with rectangles to create a flexible 3D closed pipe in MATLAB?

matlab

surface

rectangles

更新：`getknots`

如何使用矩形放样以在 MATLAB 中创建灵活的 3D 封闭管道？

How to loft with rectangles to create a flexible 3D closed pipe in MATLAB?

matlab

surface

rectangles

更新：getknots

更新：`getknots`