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 秒。