【matlab教程】20、简单网格细分

问题

matlab处理网格时，有时需要将网格细分
思路参考：https://hanspond.github.io/2018/11/27/%E7%AE%80%E5%8D%95%E7%BD%91%E6%A0%BC%E7%BB%86%E5%88%86%201to4%20Mesh%20Subdivision/index.html

将每一个面的每一条边取中点并连接，形成新的三角形，原来一个面细分成四个面
原博主的代码只适合封闭曲面的细分，如果曲面是不封闭的，则拓扑结构会乱，如：
代码：

%% created by H_P 20181119: subdivide mesh into 4 smaller meshes in a simple manner

function [ Divided_faces, Divided_vertices ] = sub_test( faces, vertices )
%SUB_TEST 此处显示有关此函数的摘要
%   此处显示详细说明

mesh_length_Large=length(faces);  
Divided_faces=1:mesh_length_Large*4*3;
Divided_faces=(reshape(Divided_faces,[3,mesh_length_Large*4]))';
Divided_vertices=repelem(vertices,4,1);

for n_for=1:mesh_length_Large
    v1=faces(n_for,1);
    v2=faces(n_for,2);
    v3=faces(n_for,3);
    new_1=(vertices(v1,:)+vertices(v2,:))/2;
    new_2=(vertices(v3,:)+vertices(v2,:))/2;
    new_3=(vertices(v1,:)+vertices(v3,:))/2;
    
    Divided_vertices((n_for-1)*12+2,:)=    new_1;
    Divided_vertices((n_for-1)*12+3,:)=    new_3;
    Divided_vertices((n_for-1)*12+4,:)=    new_1;
    Divided_vertices((n_for-1)*12+6,:)=    new_2;
    Divided_vertices((n_for-1)*12+7,:)=    new_2;
    Divided_vertices((n_for-1)*12+8,:)=    new_3;
    Divided_vertices((n_for-1)*12+10,:)=    new_1;
    Divided_vertices((n_for-1)*12+11,:)=    new_2;
    Divided_vertices((n_for-1)*12+12,:)=    new_3;
end
end

调用：

>> [v,f]=read_obj('test.obj');
>> [ Divided_faces, Divided_vertices ] = sub_test( f', v');
>> writeOBJ('test_sub.obj',Divided_vertices,Divided_faces);

结果：
细分前：

细分后：

解决

我将代码改进了一下，将分割条件改为网格面积大于阈值就分割，否则不分割
代码：

function [ divided_vertices ,divided_faces ] = subdivide_mesh01( vertices,faces ,threshold  )
%SUBDIVIDE_MESH 将曲面简单细分为1to4
%   输入：
%       vertices：细分前的点集n*3
%       faces：细分前的面集m*3
%       threshold：网格面积小于这个就细分
%   输出：
%       divided_vertices：细分后的点集
%       divided_faces：细分后的面集
tic%403秒
size_faces=size(faces,1);
delete_flag=[];
for m=1:size_faces
    inx1=faces(m,1);
    inx2=faces(m,2);
    inx3=faces(m,3);
    v1=vertices(inx1,:);
    v2=vertices(inx2,:);
    v3=vertices(inx3,:);
    s=my_area(v1,v2,v3);
    if s>threshold
        delete_flag=[delete_flag;m];
        %细分
        new1=(v1+v2)/2;
        new2=(v2+v3)/2;
        new3=(v3+v1)/2;
        [inx_new1,vertices]=my_index(new1,vertices);
        [inx_new2,vertices]=my_index(new2,vertices);
        [inx_new3,vertices]=my_index(new3,vertices);
        
        faces_temp=[inx3,inx_new2,inx_new3;
            inx_new1,inx_new2,inx_new3;
            inx_new1,inx_new2,inx2;
            inx_new3,inx_new1,inx1];
        faces=[faces;faces_temp];
    end

end
toc
faces(delete_flag,:)=[];%删除细分过后的原始面
divided_vertices=vertices;
divided_faces=faces;

end

function [dist]=my_dist(v1,v2)
%输入两个点，得到它们之间的距离
dist=power( (v1-v2)*((v1-v2)'),1/2);
end

function [area]=my_area(v1,v2,v3)
%输入三个点，得到由它们构成的三角形的面集
    a=my_dist(v1,v2);
    b=my_dist(v2,v3);
    c=my_dist(v1,v3);
    p=(a+b+c)/2;
    area=power(p*(p-a)*(p-b)*(p-c),1/2);
end

function [inx,vertices]=my_index(v,vertices)
%输入一个点，返回该点在点集中的索引行
%如果该点存在，返回索引行，如果该点不存在，插入该点，返回索引行
mem = ismember(vertices,v,'rows');
inx=find(mem==1);
if isempty(inx)%如果inx为空
    vertices=[vertices;v];
    inx=size(vertices,1);
end
end

调用：

>> [ divided_vertices ,divided_faces ] = subdivide_mesh01( v',f' ,0.06  );
时间已过 4.625456 秒。
>> writeOBJ('test_sub02.obj',divided_vertices,divided_faces);

结果：
细分前：

细分后：

优化

将元素运算换成矩阵运算，见文章
【matlab教程】21、matlab优化一：将对单个元素的操作转换成矩阵之间的运算