%% 非均匀网格法
clc
clear
close all
%% 载入数据
P = ascread('bun000.asc');
P = P{2};
figure;
x=P(1,1:1:end);
y=P(2,1:1:end);
z=P(3,1:1:end);
c=z+1;
scatter3(x,y,z,1.2,c,'filled');
colorbar
view(2)
title('原始数据')
X = P(1,:);
Y = P(2,:);
Z = P(3,:);
step=0.004;% 包围盒尺寸，控制精简率的唯一参数
std_th = 1.2;% 分割阈值，法向量标准差
tic;
nor = norcur(P,8);
% [idx00,dis00] = rangesearch(P',P',20);
% [~,nor] = My_ISS_ver1(P, 20, 0.4,0.1,idx00,dis00);
Xmin=min(X);Xmax=max(X);
Ymin=min(Y);Ymax=max(Y);
Zmin=min(Z);Zmax=max(Z);
jingjian = zeros(size(P));
ite = 0;
for ix=Xmin:step:Xmax
xindex=find((X>=ix&X<(ix+step)));
if isempty(xindex)==1
continue
end
Xselect=X(xindex);
Yselect=Y(xindex);
Zselect=Z(xindex);
nor_slt_dim1 = nor(:,xindex);
for iy=Ymin:step:Ymax
yindex=find(Yselect>=iy&Yselect<(iy+step));
if isempty(yindex)==1
continue
end
Xslt=Xselect(yindex);
Yslt=Yselect(yindex);
Zslt=Zselect(yindex);
nor_slt_dim2 = nor_slt_dim1(:,yindex);
for iz=Zmin:step:Zmax
zindex=find(Zslt>=iz&Zslt<(iz+step));
if isempty(zindex)==1
continue
end
Xsel=Xslt(zindex);
Ysel=Yslt(zindex);
Zsel=Zslt(zindex);
nor_slt_dim3 = nor_slt_dim2(:,zindex);
% 开始非均匀细分
if length(zindex)==1
ite = ite + 1;
jingjian(:,ite) = [Xsel;Ysel;Zsel];
continue
end
dev = angle_deviation(nor_slt_dim3);
if dev<std_th
ite = ite + 1;
jingjian(:,ite) = mean([Xsel;Ysel;Zsel],2);
% Distance=sqrt((Xsel-mean(Xsel)).^2+(Ysel-mean(Ysel)).^2+(Zsel-mean(Zsel)).^2);
% index=find(Distance==min(Distance));
% if(size(index,2)>1)
% index=index(1);
% end
% ite = ite + 1;
% jingjian(:,ite) = [Xsel(index);Ysel(index);Zsel(index)];
continue
end
if dev>=std_th
%输入：Xsel,Ysel,Zsel,nor_slt_dim3,ix,iy,iz,step,ite,jingjian
[jingjian,ite] = interation(Xsel,Ysel,Zsel,nor_slt_dim3,ix,iy,iz,step,ite,jingjian,std_th);
end
end
end
end
jingjian(:,ite+1:end) = [];% 删除矩阵jingjian中的空位置
figure;
x=jingjian(1,:);
y=jingjian(2,:);
z=jingjian(3,:);
c=z+1;
scatter3(x,y,z,1.2,c,'filled');
colorbar
view(2)
title('精简后的点云');
toc;
disp(['精简率为：',num2str((size(P,2)-size(jingjian,2))/size(P,2)*100),'%']);
function [jingjian,ite] = interation(Xsel,Ysel,Zsel,nor_slt_dim3,ix,iy,iz,step,ite,jingjian,std_th)
for iix = [ix,ix+step/2]
xidx = find((Xsel>=iix&Xsel<(iix+step/2)));
if isempty(xidx)
continue
end
x0 = Xsel(xidx);
y0 = Ysel(xidx);
z0 = Zsel(xidx);
nor0 = nor_slt_dim3(:,xidx);
for iiy = [iy,iy+step/2]
yidx = find((y0>=iiy&y0<(iiy+step/2)));
if isempty(yidx)
continue
end
x1 = x0(yidx);
y1 = y0(yidx);
z1 = z0(yidx);
nor1 = nor0(:,yidx);
for iiz = [iz,iz+step/2]
zidx = find((z1>=iiz&z1<(iiz+step/2)));
if isempty(zidx)
continue
end
x2 = x1(zidx);
y2 = y1(zidx);
z2 = z1(zidx);
nor2 = nor1(:,zidx);
if length(zidx)==1
ite = ite + 1;
jingjian(:,ite) = [x2;y2;z2];
continue
end
dev = angle_deviation(nor2);
if dev<std_th
ite = ite + 1;
jingjian(:,ite) = mean([x2;y2;z2],2);
else
[jingjian,ite] = interation(x2,y2,z2,nor2,iix,iiy,iiz,step/2,ite,jingjian,std_th);
end
end
end
end
end
function dev = angle_deviation(nor)
nor_bar = mean(nor,2);
nor_bar = nor_bar./sqrt(dot(nor_bar,nor_bar));
angle = zeros(1,size(nor,2));
for i = 1:size(nor,2)
angle(i) = acos(dot(nor(:,i),nor_bar));
end
dev = std(angle);
% if isnan(dev)
% error('!')
% end
end
评论0