How to create MATLAB code for making bilinear coon's surface patch from cubic Bezier curves

Question

I have following data from which I need to create coon's patch. Matrix of polygon vertices given below:

P0u = \[0,0,0; 0.2357,0.2357,0.3333; 1.1785,0.2357,0.3333; 1.4142,0,0\];

P0w = \[1.4142,0,0; 1.1785,0.2357,0.3333; 1.1785,1.1785,0.3333; 1.4142,1.4142,0\];

P1u = \[1.4142,1.4142,0; 1.1785,1.1785,0.3333; 0.2357,1.1785,0.3333; 0,1.4142,0\];

P1w = \[0,1.4142,0; 0.2357,1.1785,0.3333; 0.2357,0.2357,0.3333; 0,0,0\];

What can I try next?

I have prepared partial MATLAB code but do not know how to proceed further. See the below code:

clc
clear
close all

%% Define input parameters
A1 = [0,0,0];
A2 = [0.2357,0.2357,0.3333];
A3 = [1.1785,0.2357,0.3333];
A4 = [1.4142,0,0];
B1 = [1.4142,0,0];
B2 = [1.1785,0.2357,0.3333];
B3 = [1.1785,1.1785,0.3333];
B4 = [1.4142,1.4142,0];
C1 = [1.4142,1.4142,0];
C2 = [1.1785,1.1785,0.3333];
C3 = [0.2357,1.1785,0.3333];
C4 = [0,1.4142,0];
D1 = [0,1.4142,0];
D2 = [0.2357,1.1785,0.3333];
D3 = [0.2357,0.2357,0.3333];
D4 = [0,0,0];

P0u = [A1;A2;A3;A4]; % Creating matrix for polygon vertices
P0w = [B1;B2;B3;B4]; % Creating matrix for polygon vertices
P1u = [C1;C2;C3;C4]; % Creating matrix for polygon vertices
P1w = [D1;D2;D3;D4]; % Creating matrix for polygon vertices


[r1,s1] = size(P0u); % getting size of matrix in terms of rows and columns
[r2,s2] = size(P0w); % getting size of matrix in terms of rows and columns
[r3,s3] = size(P1u); % getting size of matrix in terms of rows and columns
[r4,s4] = size(P1w); % getting size of matrix in terms of rows and columns

n1 = r1-1; % n+1 represents number of vertices of the polygon
n2 = r2-1; % n+1 represents number of vertices of the polygon
n3 = r3-1; % n+1 represents number of vertices of the polygon
n4 = r4-1; % n+1 represents number of vertices of the polygon

np = 20; % represents number of equi-distance points on the bezier curve
t = linspace(0,1,np);

%% Plot polygon

for k1 = 1:n1
    figure(1)
    plot([P0u(k1,1),P0u(k1+1,1)], [P0u(k1,2),P0u(k1+1,2)], 'r', 'LineWidth', 2)
hold all
grid on
view(45,45)
end

for k2 = 1:n2
    plot([P0w(k2,1),P0w(k2+1,1)], [P0w(k2,2),P0w(k2+1,2)], 'r', 'LineWidth', 2)
end

for k3 = 1:n3
    plot([P1u(k3,1),P1u(k3+1,1)], [P1u(k3,2),P1u(k3+1,2)], 'r', 'LineWidth', 2)
end

for k4 = 1:n4
    plot([P1w(k4,1),P1w(k4+1,1)], [P1w(k4,2),P1w(k4+1,2)], 'r', 'LineWidth', 2)
end

%% Generate the points on the bezier curve

for j1 = 1:np
    P1 = [0,0,0];
    for i1 = 0:n1
        M1(i1+1) = (factorial(n1)/(factorial(i1)*factorial(n1-i1)))*((t(j1))^i1)*((1-t(j1))^(n1-i1));
        P1 = P1 + P0u(i1+1,:)*M1(i1+1);
    end
    Q1(j1,:) = P1;
end

for j2 = 1:np
    P2 = [0,0,0];
    for i2 = 0:n2
        M2(i2+1) = (factorial(n2)/(factorial(i2)*factorial(n2-i2)))*((t(j2))^i2)*((1-t(j2))^(n2-i2));
        P2 = P2 + P0w(i2+1,:)*M2(i2+1);
    end
    Q2(j2,:) = P2;
end

for j3 = 1:np
    P3 = [0,0,0];
    for i3 = 0:n3
        M3(i3+1) = (factorial(n3)/(factorial(i3)*factorial(n3-i3)))*((t(j3))^i3)*((1-t(j3))^(n3-i3));
        P3 = P3 + P1u(i3+1,:)*M3(i3+1);
    end
    Q3(j3,:) = P3;
end

for j4 = 1:np
    P4 = [0,0,0];
    for i4 = 0:n4
        M4(i4+1) = (factorial(n4)/(factorial(i4)*factorial(n4-i4)))*((t(j4))^i4)*((1-t(j4))^(n4-i4));
        P4 = P4 + P1w(i4+1,:)*M4(i4+1);
    end
    Q4(j4,:) = P4;
end

%% Plot the bezier curve from the obtained points

for l1 = 1:np-1
    plot([Q1(l1,1),Q1(l1+1,1)],[Q1(l1,2),Q1(l1+1,2)], 'b', 'LineWidth', 2);
end

for l2 = 1:np-1
    plot([Q2(l2,1),Q2(l2+1,1)],[Q2(l2,2),Q2(l2+1,2)], 'b', 'LineWidth', 2);
end

for l3 = 1:np-1
    plot([Q3(l3,1),Q3(l3+1,1)],[Q3(l3,2),Q3(l3+1,2)], 'b', 'LineWidth', 2);
end

for l4 = 1:np-1
    plot([Q4(l4,1),Q4(l4+1,1)],[Q4(l4,2),Q4(l4+1,2)], 'b', 'LineWidth', 2);
end

Answer 1

%% Clearing workspace

clc
clear
close all

%% Defining input parameters
% Defining individual point on polygon

A1 = [0,0,0];
A2 = [0.2357,0.2357,0.3333];
A3 = [1.1785,0.2357,0.3333];
A4 = [1.4142,0,0];
B1 = [1.4142,0,0];
B2 = [1.1785,0.2357,0.3333];
B3 = [1.1785,1.1785,0.3333];
B4 = [1.4142,1.4142,0];
C1 = [1.4142,1.4142,0];
C2 = [1.1785,1.1785,0.3333];
C3 = [0.2357,1.1785,0.3333];
C4 = [0,1.4142,0];
D1 = [0,1.4142,0];
D2 = [0.2357,1.1785,0.3333];
D3 = [0.2357,0.2357,0.3333];
D4 = [0,0,0];

%% Plot polygon

P0u = [A1;A2;A3;A4]; % Creating matrix for polygon vertices
P0w = [B1;B2;B3;B4]; % Creating matrix for polygon vertices
P1u = [C1;C2;C3;C4]; % Creating matrix for polygon vertices
P1w = [D1;D2;D3;D4]; % Creating matrix for polygon vertices

[r,s] = size(P0u); % getting size of matrix in terms of rows and columns
[r,s] = size(P0w); % getting size of matrix in terms of rows and columns
[r,s] = size(P1u); % getting size of matrix in terms of rows and columns
[r,s] = size(P1w); % getting size of matrix in terms of rows and columns

n = r-1; % n+1 represents number of vertices of the polygon
np = 50; % represents number of equi-distance points on the bezier curve
t = linspace(0,1,np);

hold all
grid on
view(45,45)
figure(1)
xlabel('X-Axis')
ylabel('Y-Axis')
zlabel('Z-Axis')

for k = 1:n
    plot3([P0u(k,1),P0u(k+1,1)], [P0u(k,2),P0u(k+1,2)], [P0u(k,3),P0u(k+1,3)], 'r', 'LineWidth', 1);
    plot3([P0w(k,1),P0w(k+1,1)], [P0w(k,2),P0w(k+1,2)], [P0w(k,3),P0w(k+1,3)], 'r', 'LineWidth', 1);
    plot3([P1u(k,1),P1u(k+1,1)], [P1u(k,2),P1u(k+1,2)], [P1u(k,3),P1u(k+1,3)], 'r', 'LineWidth', 1);
    plot3([P1w(k,1),P1w(k+1,1)], [P1w(k,2),P1w(k+1,2)], [P1w(k,3),P1w(k+1,3)], 'r', 'LineWidth', 1);
end

%% Generate the points on the bezier curve

for j = 1:np
    P1 = [0,0,0];
    P2 = [0,0,0];
    P3 = [0,0,0];
    P4 = [0,0,0];
    for i = 0:n
        M(i+1) = (factorial(n))/(factorial(i)*factorial(n-i))*((t(j))^i)*((1-t(j))^(n-i));
        P1 = P1 + P0u(i+1,:).*M(i+1);
        P2 = P2 + P0w(i+1,:).*M(i+1);
        P3 = P3 + P1u(i+1,:).*M(i+1);
        P4 = P4 + P1w(i+1,:).*M(i+1);
    end
    Q1(j,:) = P1;
    Q2(j,:) = P2;
    Q3(j,:) = P3;
    Q4(j,:) = P4;
end
Q3 = flipud(Q3);
Q4 = flipud(Q4);
%% Plot the bezier curve from the obtained points
for l = 1:np-1
    plot3([Q1(l,1),Q1(l+1,1)],[Q1(l,2),Q1(l+1,2)], [Q1(l,3),Q1(l+1,3)], 'b', 'LineWidth', 1);
    plot3([Q2(l,1),Q2(l+1,1)],[Q2(l,2),Q2(l+1,2)], [Q2(l,3),Q2(l+1,3)], 'b', 'LineWidth', 1);
    plot3([Q3(l,1),Q3(l+1,1)],[Q3(l,2),Q3(l+1,2)], [Q3(l,3),Q3(l+1,3)], 'b', 'LineWidth', 1);
    plot3([Q4(l,1),Q4(l+1,1)],[Q4(l,2),Q4(l+1,2)], [Q4(l,3),Q4(l+1,3)], 'b', 'LineWidth', 1);
end

%% Bilinear surface

npu = 50; % Number of points in u-direction
npw = 50; % Number of points in w-direction
u = linspace(0,1,npu);
w = linspace(0,1,npw);
Px = zeros(npu,npw);
Py = Px; 
Pz = Px;

for i = 1:npu
    for j = 1:npw
        Px(i,j) = Q1(i,1)*(1-w(j)) + Q2(j,1)*(u(i)) + Q3(i,1)*w(j) + Q4(j,1)*(1-u(i)) - Q1(1,1)*(1-u(i))*(1-w(j)) - Q1(end,1)*(u(i))*(1-w(j)) - Q3(1,1)*(1-u(i))*(w(j)) -Q3 (end,1)*u(i)*w(j);
        Py(i,j) = Q1(i,2)*(1-w(j)) + Q2(j,2)*(u(i)) + Q3(i,2)*w(j) + Q4(j,2)*(1-u(i)) - Q1(1,2)*(1-u(i))*(1-w(j)) - Q1(end,2)*(u(i))*(1-w(j)) - Q3(1,2)*(1-u(i))*(w(j)) - Q3(end,2)*u(i)*w(j);
        Pz(i,j) = Q1(i,3)*(1-w(j)) + Q2(j,3)*(u(i)) + Q3(i,3)*w(j) + Q4(j,3)*(1-u(i)) - Q1(1,3)*(1-u(i))*(1-w(j)) - Q1(end,3)*(u(i))*(1-w(j)) - Q3(1,3)*(1-u(i))*(w(j)) - Q3(end,3)*u(i)*w(j); 
    end
end
surf(Px,Py,Pz)
xlabel('X-Axis')
ylabel('Y-Axis')
zlabel('Z-Axis')