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

matlab优化的一个方向就是将对单个元素的操作转换成矩阵之间的运算，因为matlab的老本行就是矩阵运算，特别快
矩阵运算是matlab编程的核心，只要数据量上去了，什么运算都往矩阵上靠就是最快的
可谓—— matlab：面向矩阵编程

以下是我总结的几个矩阵运算
【matlab矩阵运算】01、把矩阵展开成向量
【matlab矩阵运算】02、矩阵A中减掉矩阵B所含行或者列
【matlab矩阵运算】03、矩阵中间隔取某一行或某一列
【matlab矩阵运算】04、两个矩阵对应元素操作：点乘.* 点左除./ 点右除.\ 点幂.^
【matlab矩阵运算】05、善用reshape之两个矩阵A和B你一行我一行插成新的矩阵D

下面两段代码，第一段都是对元素的操作
第二段痛定思痛改成矩阵之间的运算，速度提升很快

function [ divided_vertices ,divided_faces,divided_curv,divided_density] = subdivide_mesh( vertices,faces,curv,density  ,threshold )
%SUBDIVIDE_MESH 将曲面简单细分为1to4
%   输入：
%       vertices：细分前的点集n*3
%       faces：细分前的面集m*3
%       curv：细分前的曲率n*1
%       density：细分前密度n*1
%       threshold：网格面积小于这个就细分
%   输出：
%       divided_vertices：细分后的点集
%       divided_faces：细分后的面集
%       divided_curv：细分后的曲率
%       dicided_density：细分后密度
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,:);c1=curv(inx1);d1=density(inx1);
    v2=vertices(inx2,:);c2=curv(inx2);d2=density(inx2);
    v3=vertices(inx3,:);c3=curv(inx3);d3=density(inx3);
    s=my_area(v1,v2,v3);
    if s>threshold
        delete_flag=[delete_flag;m];
        %细分
        new1=(v1+v2)/2;new_c1=(c1+c2)/2;new_d1=(d1+d2)/2;
        new2=(v2+v3)/2;new_c2=(c2+c3)/2;new_d2=(d2+d3)/2;
        new3=(v3+v1)/2;new_c3=(c3+c1)/2;new_d3=(d3+d3)/2;
        [inx_new1,vertices,curv,density]=my_index(new1,vertices,new_c1,curv,new_d1,density);
        [inx_new2,vertices,curv,density]=my_index(new2,vertices,new_c2,curv,new_d2,density);
        [inx_new3,vertices,curv,density]=my_index(new3,vertices,new_c3,curv,new_d3,density);  
        
        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;
divided_curv=curv;
divided_density=density;

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,curv,density]=my_index(v,vertices,c,curv,d,density)
%输入一个点，返回该点在点集中的索引行
%如果该点存在，返回索引行，如果该点不存在，插入该点，返回索引行
mem = ismember(vertices,v,'rows');
inx=find(mem==1);
if isempty(inx)%如果inx为空
    vertices=[vertices;v];
    curv=[curv;c];
    density=[density;d];
    inx=size(vertices,1);
end
end

function [ divided_vertices ,divided_faces,divided_curv,divided_density] = subdivide( vertices,faces,curv,density  ,threshold )
%SUBDIVIDE_MESH 将曲面简单细分为1to4
%   输入：
%       vertices：细分前的点集n*3
%       faces：细分前的面集m*3
%       curv：细分前的曲率n*1
%       density：细分前密度n*1
%       threshold：网格面积小于这个就细分
%   输出：
%       divided_vertices：细分后的点集
%       divided_faces：细分后的面集
%       divided_curv：细分后的曲率
%       dicided_density：细分后密度
tic
%点、曲率、密度是一起变化的，所以把它们放在一起
data=[vertices,curv,density];
%1.取出超过阈值的面集，
%1.1先把面集的面积算出来
fa=faces(:,1);
fb=faces(:,2);
fc=faces(:,3);
%1.2取出三个点
p1=vertices(fa,:);
p2=vertices(fb,:);
p3=vertices(fc,:);
%1.3算出三个点之间的距离
dist12=power(sum((p1-p2).*(p1-p2),2),0.5);
dist23=power(sum((p2-p3).*(p2-p3),2),0.5);
dist31=power(sum((p3-p1).*(p3-p1),2),0.5);
%1.4算出这个面的面积
p=(dist12+dist23+dist31)/2;
area=power(p.*(p-dist12).*(p-dist23).*(p-dist31),0.5);
%1.5取出需要细分的面
need_sub_faces=faces(area>threshold,:);
num_face=size(need_sub_faces,1);

%2取出这些面里的边
n_s_f_a=need_sub_faces(:,1);
n_s_f_b=need_sub_faces(:,2);
n_s_f_c=need_sub_faces(:,3);
edge1=[n_s_f_a,n_s_f_b]';
edge2=[n_s_f_b,n_s_f_c]';
edge3=[n_s_f_c,n_s_f_a]';
edge=[edge1;edge2;edge3];
edge_1=reshape(edge,2,[])';%这样的话edge三条边为一组就在一个面上
%每条边有两个点
end_point1=edge_1(:,1);
end_point2=edge_1(:,2);
%每条边中间插一个点
new_data=(data(end_point1,:)+data(end_point2,:))/2;
data=[data;new_data];
%%接下来处理新面，这段又臭又长
%新点的索引
num_vertice=size(vertices,1);
inx=(1:(3*num_face))+num_vertice;
%完全由新点组成的面
new_face_1=reshape(inx,3,[])';
%半新半老的面，每个面有一个老点，两个新点
%先取出老点
old_faces_1=need_sub_faces(:,1);
old_faces_2=need_sub_faces(:,2);
old_faces_3=need_sub_faces(:,3);
%再取出新点
%新点按三个为一组，分别取一组中的第一个、第二个、第三个
inx_1=inx(1:3:end)';
inx_2=inx((1:3:end)+1)';
inx_3=inx((1:3:end)+2)';

%新点跟老点拼起来
new_face_2_1=[old_faces_1,inx_1,inx_3];
new_face_2_2=[old_faces_2,inx_1,inx_2];
new_face_2_3=[old_faces_3,inx_2,inx_3];
new_face_2=[new_face_2_1;new_face_2_2;new_face_2_3];
new_face=[new_face_1;new_face_2];

%去掉原始被分割的面
faces = setdiff(faces,need_sub_faces,'rows');
faces=[faces;new_face];

%接下来就是去重
%去重
[data,~,ic]=unique(data,'rows');
%new_face摊开
faces=faces';
faces=faces(:);
%摊开之后取ic里的索引
face_unique=ic(faces);
face_unique=reshape(face_unique,3,[]);
face_unique=face_unique';

%输出
divided_vertices =data(:,1:3);
divided_curv=data(:,4);
divided_density=data(:,5);
divided_faces=face_unique;

toc
end

>> [ divided_vertices_1 ,divided_faces_1,divided_curv_1,divided_density_1] = subdivide_mesh( v',f',curv,density  ,threshold );
时间已过 7.900510 秒。
>> [ divided_vertices_2 ,divided_faces_2,divided_curv_2,divided_density_2] = subdivide( v',f',curv,density  ,threshold );
时间已过 0.013031 秒。

数据量增大11倍的效果如下。可见矩阵运算对数据量增大不敏感，相反，元素运算对数据量增大相当敏感

>> [ divided_vertices_1 ,divided_faces_1,divided_curv_1,divided_density_1 ] = subdivide_mesh( v',f' ,curv,density,threshold  );
时间已过 682.153727 秒。
>> [ divided_vertices ,divided_faces,divided_curv,divided_density ] = subdivide( v',f' ,curv,density,threshold  );
时间已过 0.076227 秒。