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

Question

I have a series of rectangles and know the exact locations of their respective 4 corners. I would like to plot them and loft through each of them to create something like a 3D pipe of rectangular cross-section. The points should also not be limited to follow a straight axis. It should be flexible in taking deviations. I would also like two the ends to be patched up closed.

I have seen a similar question about lofting on your website under "How to loft with ellipses to create a 3d hollow pipe in MATLAB or Python?". The answer impressed me but was sadly for ellipses and circles. I tried to make it work with rectangles but could not figure out the logic required. I also tried patching up everything together but that led to generation of sharp edges, which I don't want. My code looks something like this:

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    

The end result should be like a very curvy pipe with rectangular cross-sections. There should not be any sharp corners.

PS: The MATLAB file is importing a notepad .txt document, where the coordinates will be inputted as shown below-

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

Desired variation in direction Image 2

A little more detailed representation pof the desired variaion in direction Image 3

Answer 1

This answer was originally posted in my first answer, but the first answer gets too crowded and difficult to track. Therefore, I now separate this answer to a second answer. This is also an alternative and different approach to solve the question.

---

If you have a set of points with good density, you can use MATLAB's built-in interpolation function interp1.

First, you need to parameterise the points. Chord Length is good enough here:

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

Then generate a new sets of u values for calculating the fitted points:

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

And then interplate the points with 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;

To summarise the for loop is now:

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

However, since you only have 10 points for each edge, you will get results like this:

Taking one of the edge curves:

The huge curvature between V1 and V2 is clearly not desirable. However, since the variation between V1 and V2 only needs to be linear, you can artificially add points between them. For example, adding two points:

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,:)];

And then do the same calculation, you get:

If you add even more:

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

But it is not perfect since you will get visible oscillation at V2 which is unavoidable because of continuity:

However, it is controllable by removing or adding points close to V2:

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

Note that you also get such oscillation with the B-spline method posted in the other answer, but it is more controllable and predictable by moving the 2nd control point closer or further to V2.

A lofted body generated with the interp1 method is like this:

To summarise, the full code to reproduce the above image is given below:

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);

If you have more points to constrain the interpolation, you can remove the lines involving 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,:)];