三维激光点云包围盒算法压缩

%求点云的平均密度 %以求出单位体积内所包含的平均数据点个数来表征点云的密度。 %(1)确定点云数据的最小长方体包围盒 clear all %清除所有记录 clc; in = importdata('Chair.mat');%加载点云数据 %% plot of the current point cloud figure(1); %图名 hold on %你在当前图的轴（坐标系）中画了一幅图，再画另一幅图时，原来的图还在，与新图共存，都看得到 axis equal %是将坐标轴的长度单位设成相等 title('Points Cloud','fontsize',14) %图名与字体尺寸 plot3(in(:,1),in(:,2),in(:,3),'g.') %g代表绿色 grid %添加网格 view(-37.5,30) %方位角和俯视角，方位角相当于球坐标中的经度，俯视角相当于球坐标中的纬度 xlabel('X');ylabel('Y');zlabel('Z'); %给XYZ坐标加标签 % [N,b]=size(in);e=0.8;K=1; %N为行数即点云总数，b为3 max_x=max(in(:,1));min_x=min(in(:,1)); max_y=max(in(:,2));min_y=min(in(:,2)); max_z=max(in(:,3));min_z=min(in(:,3)); Lx=max_x-min_x;%点云数据的长方体包围盒的边长 Ly=max_y-min_y; Lz=max_z-min_z; %(2)确定单位小立方体栅格中数据点的个数 V=Lx*Ly*Lz; %最大包围盒的大小 %则单位小立方体栅格中数据点的个数为： n=N/(V+eps); % 其中 N 为点云中数据点的总数。eps=eps(1)，是1的精度。 Ls=(e*K/(n+eps))^(1/3);%确定子立方体栅格的边长 Ls，λ是比例因子，用于调节子 %立方体栅格的边长，K 为邻近点的个数 mm=floor(Lx/(Ls+eps))+1; nn=floor(Ly/(Ls+eps))+1; ll=floor(Lz/(Ls+eps))+1;%将点云数据分割为 m × n×l个小立方体栅格 hh=floor((in(:,1)-min_x)/(Ls+eps))+1; %floor向下取整 jj=floor((in(:,2)-min_y)/(Ls+eps))+1; kk=floor((in(:,3)-min_z)/(Ls+eps))+1; %确定各个点云在在那个方格子中 %求包围盒中心坐标 center=zeros(mm*nn*ll,3); for i=1:mm center_x=min_x+Ls*(i-0.5); for j=1:nn center_y=min_y+Ls*(j-0.5); for k=1:ll center_z=min_z+Ls*(k-0.5); center((i-1)*nn*ll+(j-1)*ll+k,1)=center_x; center((i-1)*nn*ll+(j-1)*ll+k,2)=center_y; center((i-1)*nn*ll+(j-1)*ll+k,3)=center_z; end end end %求中心点与各个点云之间的距离 distance=zeros(N,4); for i=1:N distance(i,1)=hh(i); distance(i,2)=jj(i); distance(i,3)=kk(i); distance(i,4)=sqrt((center((hh(i)-1)*nn*ll+(jj(i)-1)*ll+kk(i),1)-in(i,1))^2+(center((hh(i)-1)*nn*ll+(jj(i)-1)*ll+kk(i),2)-in(i,2))^2+(center((hh(i)-1)*nn*ll+(jj(i)-1)*ll+kk(i),3)-in(i,3))^2); end %找出距离最小值作为包围盒的重心，只保留重心 %{ Yes=[]; k=1; for i=1:N for j=1:N if distance(i,1:3)==distance(j,1:3) Yes=[Yes;distance(j,:)];%可以找到相同的放在B end end strcut(k).Yes=Yes;k=k+1; Yes=[]; end Result=[]; for k=1:length(strcut) A= strcut(k).Yes(:,4); MX=max(A); Result=[Result;strcut(k).Yes(1,1:3),max(A)]; end [Y]=unique(Result,'rows');%删除相同的行，Y里面就是最后的结果 %} [Y mm nn]=unique(distance(:,1:3),'rows');%m体现b中元素在原向量（矩阵a）中的位置(最后出现的位置)；j体现原向量（矩阵a）在b中的位置（第一次出现的位置） num=length(mm); X=zeros(num,1); for i=1:num X(i)=max(distance(nn==i,4)); end Y=[Y X]; %% plot of the current point cloud figure(2); %图名 hold off %你在当前图的轴（坐标系）中画了一幅图，再画另一幅图时，原来的图还在，与新图共存，都看得到 axis equal %是将坐标轴的长度单位设成相等 title('Points Cloud','fontsize',14) %图名与字体尺寸 plot3(Y(:,1),Y(:,2),Y(:,3),'g.') %g代表绿色 grid %添加网格 view(-37.5,30) %方位角和俯视角，方位角相当于球坐标中的经度，俯视角相当于球坐标中的纬度 xlabel('X');ylabel('Y');zlabel('Z'); %给XYZ坐标加标签 % %% Run program p=Y(:,1:3); [t]=MyCrust(p); [w]=MyCrust(in); %% plot of the oyput triangulation figure(3) hold on title('Output Triangulation','fontsize',14) axis equal trisurf(w,in(:,1),in(:,2),in(:,3),'facecolor','c','edgecolor','b')%plot della superficie trattata grid view(-37.5,30) xlabel('X');ylabel('Y');zlabel('Z'); %给XYZ坐标加标签 % %% plot of the oyput triangulation figure(4) hold on title('Output Triangulation','fontsize',14) axis equal trisurf(t,p(:,1),p(:,2),p(:,3),'facecolor','c','edgecolor','b')%plot della superficie trattata grid view(-37.5,30) xlabel('X');ylabel('Y');zlabel('Z'); %给XYZ坐标加标签 %