绘制飞机并进行后向面判定

5.绘制飞机并进行后向面判定

判断后向面：即最多能看到长方体的三个面，看不到的三个面用back-face-detection作消隐；另有定义是一个面的背面？

%define original point 3D
original_pos_x = 2;
original_pos_y = 3;
original_pos_z = 5;

%define aircraft(cuboid)
l = 60;  %length
w = 10;  %width
h = 10;  %height

%define coordinate matrix
[a,b,c]=meshgrid([0 1]);

%8 peak position
x_pos = l*a(:)+original_pos_x;
y_pos = w*b(:)+original_pos_y;
z_pos = h*c(:)+original_pos_z;

%visualization process
p=alphaShape(x_pos, y_pos, z_pos);
plot(p,'edgecolor','none','facecolor','y'), hold on;
xlabel('x');ylabel('y');zlabel('z');
axis([-10 70 -10 40 -10 40]), axis on;
camlight
%draw coordinate axis
plot3([-100, 100], [0, 0], [0, 0],'k','linewidth',0.5), hold on;
plot3([0, 0], [-100, 100], [0, 0],'k','linewidth',0.5), hold on;
plot3([0, 0], [0, 0], [-100, 100],'k','linewidth',0.5), hold on;
text(original_pos_x+l, 0, 0, 'X axis');
text(0, original_pos_y+w, 0, 'Y axis');
text(0, 0, original_pos_z+h, 'Z axis');

为方便观察，需要对每个顶点进行标记（该标记顺序按照$alphaShape$函数中坐标的顺序进行）

figure(2);
    p=alphaShape(x_pos, y_pos, z_pos);
    plot(p,'edgecolor','none','facecolor','c'), hold on;
    xlabel('x');ylabel('y');zlabel('z');
    axis([-10 70 -10 40 -10 40]), axis on; 
    title('观察顶点编号图');
    camlight
    %remark peaks selected
    text(original_pos_x, original_pos_y, original_pos_z, 'P1');
    text(original_pos_x, original_pos_y+w, original_pos_z, 'P2');
    text(original_pos_x, original_pos_y, original_pos_z+h, 'P5');
    text(original_pos_x+l, original_pos_y+w, original_pos_z, 'P4');
    text(original_pos_x+l, original_pos_y+w, original_pos_z+h, 'P8');
    text(original_pos_x+l, original_pos_y, original_pos_z, 'P3');
    text(original_pos_x, original_pos_y+w, original_pos_z+h, 'P6');
    text(original_pos_x+l, original_pos_y, original_pos_z+h, 'P7');

更清楚地显示：

接下来进行法向量的计算

%back-face detection
    %define view vector
    v_view = [1 1 1];
    %define matrix
    peak_mat = [x_pos y_pos z_pos];
    %define vectors
    v_12 = peak_mat(1, :)-peak_mat(2, :);
    v_15 = peak_mat(1, :)-peak_mat(5, :);
    v_13 = peak_mat(1, :)-peak_mat(3, :);
    v_86 = peak_mat(8, :)-peak_mat(6, :);
    v_87 = peak_mat(8, :)-peak_mat(7, :);
    v_84 = peak_mat(8, :)-peak_mat(4, :);
    %calculate surface normal vectors
    vn_s1265 = cross(v_15, v_12);
    vn_s1573 = cross(v_13, v_15);
    vn_s1243 = cross(v_12, v_13);
    vn_s8657 = cross(v_86, v_87);
    vn_s8624 = cross(v_84, v_86);
    vn_s8437 = cross(v_87, v_84);
    %vn times v_view
    v_view_1265 = dot(vn_s1265, v_view);
    v_view_1573 = dot(vn_s1573, v_view);
    v_view_1243 = dot(vn_s1243, v_view);
    v_view_8657 = dot(vn_s8657, v_view);
    v_view_8624 = dot(vn_s8624, v_view);
    v_view_8437 = dot(vn_s8437, v_view);

• 注意法向量应该是指向长方体外的，即右手定则逆时针（矢量的叉积的方向指向体外）
• 每个面按照4个顶点定义，如s1265

最后判断后向面

视线向量和平面法向量乘积为负，则为前向面，可以参考如下博客：

https://blog.csdn.net/JY_95/article/details/52516932?utm_medium=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromBaidu-1.control&depth_1-utm_source=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromBaidu-1.control

    %define view line
    k = 100;
    plot3([0, k*v_view(1, 1)], [0, k*v_view(1, 2)], [0, k*v_view(1, 3)], 'r', 'linewidth',1);

红色线为观察视线

法向量与视线向量点乘的结果为：

图片

可以看出，s1243, s1265, s1573都为负数，根据判定定理，此三面为前向面，剩下三个面为后向面

这一结果与验证的图像相符合

作旋转：

%rotation process
    %define parameters
    x_f = 0;
    y_f = 3;
    z_f = 2;
    theta = pi/2;
    pos = [x_pos y_pos z_pos ones(length(x_pos), 1)];
    %define trans matrix
    T1 = [1 0 0 0;
        0 1 0 0;
        0 0 1 0;
        -x_f -y_f -z_f 1;];
    T3 = [1 0 0 0;
        0 1 0 0;
        0 0 1 0;
        x_f y_f z_f 1;];
    T_rot_z = [cos(theta) sin(theta) 0 0;
               -sin(theta) cos(theta) 0 0;
               0 0 1 0;
               0 0 0 1;];
    T_rot_x = [1 0 0 0;
               0 cos(theta) sin(theta) 0;
               0 -sin(theta) cos(theta) 0
               0 0 0 1];
    T_rot_y = [cos(theta) 0 -sin(theta) 0;
               0 1 0 0;
               sin(theta) 0 cos(theta) 0;
               0 0 0 1];
    %transforming process
    pos_rot = pos*T1*T_rot_x*T3;
    
    p=alphaShape(x_pos, y_pos, z_pos);
    plot(p,'edgecolor','none','facecolor','y'), hold on;
    xlabel('x');ylabel('y');zlabel('z');
    axis([-50 50 -50 50 -50 50]), axis on;
    camlight
    title('飞机建模旋转图');
    %draw coordinate axis
    plot3([-100, 100], [0, 0], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [-100, 100], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [0, 0], [-100, 100],'k','linewidth',0.5), hold on;
    p1=alphaShape(pos_rot(:, 1), pos_rot(:, 2), pos_rot(:, 3));
    plot(p1,'edgecolor','none','facecolor','r');

代码汇总：

%define parameters
    %define original point 3D
    original_pos_x = 2;
    original_pos_y = 3;
    original_pos_z = 5;
    %define aircraft(cuboid)
    l = 60;  %length
    w = 10;  %width
    h = 10;  %height
    %define coordinate matrix
    [a,b,c]=meshgrid([0 1]);
    %8 peak position
    x_pos = l*a(:)+original_pos_x;
    y_pos = w*b(:)+original_pos_y;
    z_pos = h*c(:)+original_pos_z;

%back-face detection
    %define view vector
    v_view = [1 1 1];
    %define matrix
    peak_mat = [x_pos y_pos z_pos];
    %define vectors
    v_12 = peak_mat(1, :)-peak_mat(2, :);
    v_15 = peak_mat(1, :)-peak_mat(5, :);
    v_13 = peak_mat(1, :)-peak_mat(3, :);
    v_86 = peak_mat(8, :)-peak_mat(6, :);
    v_87 = peak_mat(8, :)-peak_mat(7, :);
    v_84 = peak_mat(8, :)-peak_mat(4, :);
    %calculate surface normal vectors
    vn_s1265 = cross(v_15, v_12);
    vn_s1573 = cross(v_13, v_15);
    vn_s1243 = cross(v_12, v_13);
    vn_s8657 = cross(v_86, v_87);
    vn_s8624 = cross(v_84, v_86);
    vn_s8437 = cross(v_87, v_84);
    %vn times v_view
    v_view_1265 = dot(vn_s1265, v_view);
    v_view_1573 = dot(vn_s1573, v_view);
    v_view_1243 = dot(vn_s1243, v_view);
    v_view_8657 = dot(vn_s8657, v_view);
    v_view_8624 = dot(vn_s8624, v_view);
    v_view_8437 = dot(vn_s8437, v_view);
    
%rotation process
    %define parameters
    x_f = 0;
    y_f = 3;
    z_f = 2;
    theta = pi/2;
    pos = [x_pos y_pos z_pos ones(length(x_pos), 1)];
    %define trans matrix
    T1 = [1 0 0 0;
        0 1 0 0;
        0 0 1 0;
        -x_f -y_f -z_f 1;];
    T3 = [1 0 0 0;
        0 1 0 0;
        0 0 1 0;
        x_f y_f z_f 1;];
    T_rot_z = [cos(theta) sin(theta) 0 0;
               -sin(theta) cos(theta) 0 0;
               0 0 1 0;
               0 0 0 1;];
    T_rot_x = [1 0 0 0;
               0 cos(theta) sin(theta) 0;
               0 -sin(theta) cos(theta) 0
               0 0 0 1];
    T_rot_y = [cos(theta) 0 -sin(theta) 0;
               0 1 0 0;
               sin(theta) 0 cos(theta) 0;
               0 0 0 1];
    %transforming process
    pos_rot = pos*T1*T_rot_x*T3;
    
%visualization process
figure(1);
    p=alphaShape(x_pos, y_pos, z_pos);
    plot(p,'edgecolor','none','facecolor','y'), hold on;
    xlabel('x');ylabel('y');zlabel('z');
    axis([-10 70 -10 40 -10 40]), axis on;
    camlight
    title('飞机初步建模图')
    %draw coordinate axis
    plot3([-100, 100], [0, 0], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [-100, 100], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [0, 0], [-100, 100],'k','linewidth',0.5), hold on;
    text(original_pos_x+l, 0, 0, 'X axis');
    text(0, original_pos_y+w, 0, 'Y axis');
    text(0, 0, original_pos_z+h, 'Z axis');
figure(2);
    p=alphaShape(x_pos, y_pos, z_pos);
    plot(p,'edgecolor','none','facecolor','c'), hold on;
    xlabel('x');ylabel('y');zlabel('z');
    axis([-10 70 -10 40 -10 40]), axis on; 
    title('观察顶点编号图');
    camlight
    %remark peaks selected
    text(original_pos_x, original_pos_y, original_pos_z, 'P1');
    text(original_pos_x, original_pos_y+w, original_pos_z, 'P2');
    text(original_pos_x, original_pos_y, original_pos_z+h, 'P5');
    text(original_pos_x+l, original_pos_y+w, original_pos_z, 'P4');
    text(original_pos_x+l, original_pos_y+w, original_pos_z+h, 'P8');
    text(original_pos_x+l, original_pos_y, original_pos_z, 'P3');
    text(original_pos_x, original_pos_y+w, original_pos_z+h, 'P6');
    text(original_pos_x+l, original_pos_y, original_pos_z+h, 'P7');
    plot3([-100, 100], [0, 0], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [-100, 100], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [0, 0], [-100, 100],'k','linewidth',0.5), hold on;
    text(original_pos_x+l, 0, 0, 'X axis');
    text(0, original_pos_y+2*w, 0, 'Y axis');
    text(0, 0, original_pos_z+2*h, 'Z axis');
figure(3);
    p=alphaShape(x_pos, y_pos, z_pos);
    plot(p,'edgecolor','none','facecolor','g'), hold on;
    xlabel('x');ylabel('y');zlabel('z');
    axis([-10 70 -10 40 -10 40]), axis on; 
    title('后向面验证图');
    camlight
    %remark peaks selected
    text(original_pos_x, original_pos_y, original_pos_z, 'P1');
    text(original_pos_x, original_pos_y+w, original_pos_z, 'P2');
    text(original_pos_x, original_pos_y, original_pos_z+h, 'P5');
    text(original_pos_x+l, original_pos_y+w, original_pos_z, 'P4');
    text(original_pos_x+l, original_pos_y+w, original_pos_z+h, 'P8');
    text(original_pos_x+l, original_pos_y, original_pos_z, 'P3');
    text(original_pos_x, original_pos_y+w, original_pos_z+h, 'P6');
    text(original_pos_x+l, original_pos_y, original_pos_z+h, 'P7');
    plot3([-100, 100], [0, 0], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [-100, 100], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [0, 0], [-100, 100],'k','linewidth',0.5), hold on;
    text(original_pos_x+l, 0, 0, 'X axis');
    text(0, original_pos_y+2*w, 0, 'Y axis');
    text(0, 0, original_pos_z+2*h, 'Z axis');
    %define view line
    k = 100;
    plot3([0, k*v_view(1, 1)], [0, k*v_view(1, 2)], [0, k*v_view(1, 3)], 'r', 'linewidth',1);
figure(4);
    p=alphaShape(x_pos, y_pos, z_pos);
    plot(p,'edgecolor','none','facecolor','y'), hold on;
    xlabel('x');ylabel('y');zlabel('z');
    axis([-50 50 -50 50 -50 50]), axis on;
    camlight
    title('飞机建模旋转图');
    %draw coordinate axis
    plot3([-100, 100], [0, 0], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [-100, 100], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [0, 0], [-100, 100],'k','linewidth',0.5), hold on;
    p1=alphaShape(pos_rot(:, 1), pos_rot(:, 2), pos_rot(:, 3));
    plot(p1,'edgecolor','none','facecolor','r');
    

%define parameters
    %define original point 3D
    original_pos_x = 2;
    original_pos_y = 3;
    original_pos_z = 5;
    %define aircraft(cuboid)
    l = 60;  %length
    w = 10;  %width
    h = 10;  %height
    %define coordinate matrix
    [a,b,c]=meshgrid([0 1]);
    %8 peak position
    x_pos = l*a(:)+original_pos_x;
    y_pos = w*b(:)+original_pos_y;
    z_pos = h*c(:)+original_pos_z;

%back-face detection
    %define view vector
    v_view = [1 1 1];
    %define matrix
    peak_mat = [x_pos y_pos z_pos];
    %define vectors
    v_12 = peak_mat(1, :)-peak_mat(2, :);
    v_15 = peak_mat(1, :)-peak_mat(5, :);
    v_13 = peak_mat(1, :)-peak_mat(3, :);
    v_86 = peak_mat(8, :)-peak_mat(6, :);
    v_87 = peak_mat(8, :)-peak_mat(7, :);
    v_84 = peak_mat(8, :)-peak_mat(4, :);
    %calculate surface normal vectors
    vn_s1265 = cross(v_15, v_12);
    vn_s1573 = cross(v_13, v_15);
    vn_s1243 = cross(v_12, v_13);
    vn_s8657 = cross(v_86, v_87);
    vn_s8624 = cross(v_84, v_86);
    vn_s8437 = cross(v_87, v_84);
    %vn times v_view
    v_view_1265 = dot(vn_s1265, v_view);
    v_view_1573 = dot(vn_s1573, v_view);
    v_view_1243 = dot(vn_s1243, v_view);
    v_view_8657 = dot(vn_s8657, v_view);
    v_view_8624 = dot(vn_s8624, v_view);
    v_view_8437 = dot(vn_s8437, v_view);
    
%rotation process
    %define parameters
    x_f = 0;
    y_f = 3;
    z_f = 2;
    theta = pi/2;
    pos = [x_pos y_pos z_pos ones(length(x_pos), 1)];
    %define trans matrix
    T1 = [1 0 0 0;
        0 1 0 0;
        0 0 1 0;
        -x_f -y_f -z_f 1;];
    T3 = [1 0 0 0;
        0 1 0 0;
        0 0 1 0;
        x_f y_f z_f 1;];
    T_rot_z = [cos(theta) sin(theta) 0 0;
               -sin(theta) cos(theta) 0 0;
               0 0 1 0;
               0 0 0 1;];
    T_rot_x = [1 0 0 0;
               0 cos(theta) sin(theta) 0;
               0 -sin(theta) cos(theta) 0
               0 0 0 1];
    T_rot_y = [cos(theta) 0 -sin(theta) 0;
               0 1 0 0;
               sin(theta) 0 cos(theta) 0;
               0 0 0 1];
    %transforming process
    pos_rot = pos*T1*T_rot_x*T3;
    
%visualization process
figure(1);
    p=alphaShape(x_pos, y_pos, z_pos);
    plot(p,'edgecolor','none','facecolor','y'), hold on;
    xlabel('x');ylabel('y');zlabel('z');
    axis([-10 70 -10 40 -10 40]), axis on;
    camlight
    title('飞机初步建模图')
    %draw coordinate axis
    plot3([-100, 100], [0, 0], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [-100, 100], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [0, 0], [-100, 100],'k','linewidth',0.5), hold on;
    text(original_pos_x+l, 0, 0, 'X axis');
    text(0, original_pos_y+w, 0, 'Y axis');
    text(0, 0, original_pos_z+h, 'Z axis');
figure(2);
    p=alphaShape(x_pos, y_pos, z_pos);
    plot(p,'edgecolor','none','facecolor','c'), hold on;
    xlabel('x');ylabel('y');zlabel('z');
    axis([-10 70 -10 40 -10 40]), axis on; 
    title('观察顶点编号图');
    camlight
    %remark peaks selected
    text(original_pos_x, original_pos_y, original_pos_z, 'P1');
    text(original_pos_x, original_pos_y+w, original_pos_z, 'P2');
    text(original_pos_x, original_pos_y, original_pos_z+h, 'P5');
    text(original_pos_x+l, original_pos_y+w, original_pos_z, 'P4');
    text(original_pos_x+l, original_pos_y+w, original_pos_z+h, 'P8');
    text(original_pos_x+l, original_pos_y, original_pos_z, 'P3');
    text(original_pos_x, original_pos_y+w, original_pos_z+h, 'P6');
    text(original_pos_x+l, original_pos_y, original_pos_z+h, 'P7');
    plot3([-100, 100], [0, 0], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [-100, 100], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [0, 0], [-100, 100],'k','linewidth',0.5), hold on;
    text(original_pos_x+l, 0, 0, 'X axis');
    text(0, original_pos_y+2*w, 0, 'Y axis');
    text(0, 0, original_pos_z+2*h, 'Z axis');
figure(3);
    p=alphaShape(x_pos, y_pos, z_pos);
    plot(p,'edgecolor','none','facecolor','g'), hold on;
    xlabel('x');ylabel('y');zlabel('z');
    axis([-10 70 -10 40 -10 40]), axis on; 
    title('后向面验证图');
    camlight
    %remark peaks selected
    text(original_pos_x, original_pos_y, original_pos_z, 'P1');
    text(original_pos_x, original_pos_y+w, original_pos_z, 'P2');
    text(original_pos_x, original_pos_y, original_pos_z+h, 'P5');
    text(original_pos_x+l, original_pos_y+w, original_pos_z, 'P4');
    text(original_pos_x+l, original_pos_y+w, original_pos_z+h, 'P8');
    text(original_pos_x+l, original_pos_y, original_pos_z, 'P3');
    text(original_pos_x, original_pos_y+w, original_pos_z+h, 'P6');
    text(original_pos_x+l, original_pos_y, original_pos_z+h, 'P7');
    plot3([-100, 100], [0, 0], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [-100, 100], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [0, 0], [-100, 100],'k','linewidth',0.5), hold on;
    text(original_pos_x+l, 0, 0, 'X axis');
    text(0, original_pos_y+2*w, 0, 'Y axis');
    text(0, 0, original_pos_z+2*h, 'Z axis');
    %define view line
    k = 100;
    plot3([0, k*v_view(1, 1)], [0, k*v_view(1, 2)], [0, k*v_view(1, 3)], 'r', 'linewidth',1);
figure(4);
    p=alphaShape(x_pos, y_pos, z_pos);
    plot(p,'edgecolor','none','facecolor','y'), hold on;
    xlabel('x');ylabel('y');zlabel('z');
    axis([-50 50 -50 50 -50 50]), axis on;
    camlight
    title('飞机建模旋转图');
    %draw coordinate axis
    plot3([-100, 100], [0, 0], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [-100, 100], [0, 0],'k','linewidth',0.5), hold on;
    plot3([0, 0], [0, 0], [-100, 100],'k','linewidth',0.5), hold on;
    p1=alphaShape(pos_rot(:, 1), pos_rot(:, 2), pos_rot(:, 3));
    plot(p1,'edgecolor','none','facecolor','r');