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

Question

我有一系列的矩形，知道它们各自4个角的确切位置。我想绘制它们并通过它们中的每一个来创建类似于矩形横截面的3D管道。这些点也不应限于遵循直轴。它应该灵活地采取偏差。我还希望两个端点被修补关闭。

我在你的网站上看到了一个类似的问题：“如何用椭圆形放置在MATLAB或Python中创建一个3d空心管道？”。答案给我留下了深刻的印象，但对于省略号和圆圈来说却很遗憾。我试图让它与矩形一起工作，但无法弄清楚所需的逻辑。我也尝试将所有东西拼凑起来，但这导致产生锋利的边缘，这是我不想要的。我的代码看起来像这样：

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

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

---

如果你有一组具有良好密度的点，你可以使用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,:)];

Answer 2

通过数据文件中的点插值无法轻松完成。我建议设计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。 Youtube上还有许多教程和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样条曲线的控制点和顶点。

1. u中的值的数量决定了最终曲线上的点数。
2. 三阶B样条曲线就足够了。
3. 曲线必须通过V1，V2，...，V10。
4. 要通过V1和V10，B样条曲线必须是钳位类型。因此，C1 = V1和C10 = V10。
5. 要通过V2，C2必须在线路上通过V2和V3，C3必须等于V2。因此，C2 - V2 = d*(V2-V3)和C3 = V2。
6. 要通过V3，V4和V5，必须满足以下条件：C3 = V2; C4 = V3; C5 = V5; C6 = V6。
7. 要通过V6，C7必须在线路上通过V5和V6，同时C6等于V6。因此，C7-V6 = d*(V6-V5)。
8. 要通过V7，C8必须在线路上通过V7和V8，C9必须等于V7。因此，C8-V7 = d*(V7-V8); C9 = V7;。
9. 从C10 = V10和C9 = V7开始，曲线经过V8和V9。
10. 最后，结可以是均匀的，也可以通过控制点之间的弦长来估计。

---

Update: 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