激光雷达点云滤波---根据TIN渐进加密算法编程

temp=load（‘C:\Users\李建柱\Documents\pix4dmapper\33\2_densification\point_cloud\33_densified_point_cloud.txt‘）;	
 A=temp（1:300001:3）;
 %以下为拟合曲面；
%{
 [qx qy]=size（A）;
 qmx=zeros（floor（qx/2）2）;
 qmx（1:floor（qx/2）1）=A（1:floor（qx/2）1）;
 qmx（1:floor（qx/2）2）=A（floor（qx/2）+1:floor（qx/2）*21）;
 qmy=zeros（floor（qx/2）2）;
 qmy（1:floor（qx/2）1）A（1:floor（qx/2）2）;
 qmy（1:floor（qx/2）2）=A（floor（qx/2）+1:floor（qx/2）*22）;
qmz=zeros（floor（qx/2）2）;
 qmz（1:floor（qx/2）1）=A（1:floor（qx/2）3）;
 qmz（1:floor（qx/2）2）=A（floor（qx/2）+1:floor（qx/2）*23）;
meshc（qmxqmyqmz）;
 %}

n=size（A1）;                     %确定点数n
min1=min（A（1:n1））; 				%确定x最小值
max1=max（A（1:n1））;				%确定x最大值
min2=min（A（1:n2））;				%确定y最小值
max2=max（A（1:n2））;					%确定y最大值
a=input（‘请输入研究区最大尺度‘）;                                  % 参数1
f=fix（（max1-min1）/a）+1;			%确定总共列数
g=fix（（max2-min2）/a）+1;
C=zeros（1n）;
D=zeros（1n）;
for e=1:n;
C（1e）=fix（（A（e1）-min1）/a）+1;			%确定每一个元素在哪一行
D（1e）=fix（（max2-A（e2））/a）+1;			%确定每一个元素在哪一列
end;
E=（D-1）*f+C;			%确定每一个点在第几个格网
clearvars -except A E n
F=zeros（nn）;
for b=1:n;
c=length（find（E==b））;
F（b1:c）=find（E==b）;					%对点进行分类，第一行非零元素为在第一个网格中的点在A中的行数，
end
W=all（F==02）;
F（W:）=[];
W=all（F==01）;
F（:W）=[];
clearvars -except A F n
F1=size（F1）;%确定F的行数
F2=size（F2）;%确定F列数
 G=zeros（F1F2）;
U=sum（F~=02）;%U为F每一行非0个数
for w=1:F1;
G（w1:U（w1））=A（F（w1:U（w1））3）;%G的第n行元素代表第n个网格中的高程
end
clearvars -except G A F n
G（G==0）=inf;  %G中为0的元素设为无穷大
[jJ]=min（G[]2）;           %J代表最小值在G中的下标j代表每行最小值
[K1k2]=size（F）;
K=zeros（K11）;
for i=1:K1;
K（i1）=F（iJ（i1））;               %由于FG对应元素相关，K为最小值在A中行数
end
L=zeros（K13）;
L（1:K11:3）=A（K1:3）;