MATLAB激光密集点云基础化操作（一）

前言

参考网站(MATLAB文档帮助中心)：https://ww2.mathworks.cn/help/index.html?s_tid=CRUX_lftnav

---

一、利用pcshow读取并显示点云------pcshow

pcshow（ptCloud） 使用位置显示点 以及存储在点云对象中的颜色。要可视化大型点云（大于一百万点），请使用 pcviewer obect。

pcshow（xyzPoints） 显示器 由矩阵指定的点。xyzPoints

pcshow（xyzPoints，color） 显示矩阵中包含的点，其中由 指定的颜色。xyzPointscolor

pcshow（xyzPoints，colorMap） 显示矩阵中包含的点，其中 由 指定的颜色。xyzPointscolorMap电脑显示（文件名）显示 存储在 指定的文件中的点云。filename

pcshow(___,Name=Value)指定选项 使用一个或多个名称-值参数以及 来自先前语法的参数。例如，将点云可视化的平面设置为 XY 平面。ViewPlane="XY"

ax = pcshow(___)返回 绘图轴。

clc
clear all
format long g;
%---------------读取显示密集点云并绘图------------------
data=pcread("test_clound.ply");
figure;
pcshow_map=pcshow(data);
title('密集点云pcshow展示')
xlabel('x(m)','Color',"R");
ylabel('y(m)','Color','R');
zlabel('z(m)','Color','R');
%---------------点云采样------------------------------
sift_data=pcdownsample(data, 'gridAverage', 0.01);
figure;
pcshow(sift_data);

二、利用scatter3函数显示点云------scatter3

scatter3(X,Y,Z) 在 X、Y 和 Z 指定的位置显示圆圈。要绘制一组坐标，请将 X、Y 和 Z 指定为等长向量。要在同一组坐标轴上绘制多组坐标，请将 X、Y 或 Z 中的至少一个指定为矩阵。

scatter3(X,Y,Z,S) 指定圆圈大小。要更改圆圈的大小，请将 S 指定为向量。要跨多组坐标指定不同的圆圈大小，请指定一个矩阵。

scatter3(X,Y,Z,S,C) 指定圆颜色。您可以为所有圆指定一种颜色，也可以更改颜色。例如，您可以通过将 C 指定为 "red" 来绘制所有红色圆。

scatter3(___,'filled') 使用前面的语法中的任何输入参数组合填充这些圆。

scatter3(___,markertype) 指定标记类型。

xyz_data=data.Location;
color=data.Color/255;
pc_x=xyz_data(:,1);
pc_y=xyz_data(:,2);
pc_z=xyz_data(:,3);
figure;
scatter3_map=scatter3(pc_x,pc_y,pc_z,0.1,color,'.');
%axis off
title('密集点云scatter3展示',"FontSize",10)

三、将点云按高程赋值色彩------mapshow/pointCloud

ptCloud = pointCloud（xyzPoints） 返回一个点云对象，其坐标由 指定。xyzPoints

ptCloud = pointCloud(___,Name=Value)除了前面的语法之外，还使用一个或多个名称-值参数设置选项。 例如，将点云的颜色设置为 红。Color=[1 0 0]

%--------------点云按高程赋值色---------------------
z=sift_data.Location(:,3);
z_min=sift_data.ZLimits(1);   %高程最小值
z_max=sift_data.ZLimits(2);   %高程最大值
jet_color=jet(256);
hsv_color=hsv(256);
len_z=length(z);
z_jet_color=zeros(len_z,3);
z_hsv_color=zeros(len_z,3);
for i=1:len_z  
    scale=floor(((z(i)-z_min)/(z_max-z_min))*255)+1;
    z_jet_color(i,:)=jet_color(scale,:);
    z_hsv_color(i,:)=hsv_color(scale,:);
end
jet_scale_point=pointCloud(sift_data.Location,'color',z_jet_color);
hsv_scale_point=pointCloud(sift_data.Location,'color',z_hsv_color);
subplot(1,2,1);
pcshow(jet_scale_point);
title('the jet color elevation map ');
subplot(1,2,2);
pcshow(hsv_scale_point);
title('the hsv color elevation map ');

四、按照x-y轴人机交互式绘制点云刨面图------input/find/scatter

x = input(prompt) 显示 prompt 中的文本并等待用户输入值后按 Return 键。用户可以输入 pi/4 或 rand(3) 之类的表达式，并可以使用工作区中的变量。

如果用户不输入任何内容直接按下 Return 键，则 input 会返回空矩阵。

如果用户在提示下输入无效的表达式，则 MATLAB 会显示相关的错误消息，然后重新显示提示。

txt = input(prompt,"s") 返回输入的文本，而不会将输入作为表达式来计算。

%------------------绘制刨面图--------------------------
%------沿着Y轴方向-------
x_text=input("请输入要切割的X值(601530.25 <x 602151.0.6):");
x_set=sift_data.Location(:,1)-x_text;
x_near_indices=find(x_set<1 & x_set>-1);
xyz_dat1=sift_data.Location(x_near_indices,:);
color_dat1=sift_data.Color(x_near_indices,:);
figure;
subplot(2,2,1);
pcshow(xyz_dat1,color_dat1);
title('the Cut block along Y axis clound map');
subplot(2,2,2);
scatter(xyz_dat1(:,2),xyz_dat1(:,3),0.01,'r','.');
%------沿着X轴方向-------
y_text=input("请输入要切割的Y值(3410054.3 <x< 3410676.8):");
y_set=sift_data.Location(:,2)-y_text;
y_near_indices=find(y_set<1 & y_set>-1);
xyz_dat2=sift_data.Location(y_near_indices,:);
color_dat2=sift_data.Color(y_near_indices,:);
subplot(2,2,3);
pcshow(xyz_dat2,color_dat2);
title('the Cut block along X axis clound map');
subplot(2,2,4);
scatter(xyz_dat2(:,1),xyz_dat2(:,3),0.01,'r','.');

五、按俯视将点云投影到平面------pcfitplane/pointCloud

model = pcfitplane（ptCloudIn，maxDistance） 将平面拟合到点云，该点云具有最大允许距离 

model = pcfitplane（ptCloudIn，maxDistance，referenceVector） 将平面拟合到具有附加方向约束的点云 由 1×3 输入指定。referenceVector

model = pcfitplane（ptCloudIn，maxDistance，referenceVector，maxAngularDistance）将平面拟合到具有指定最大角度的点云 距离。

[model,inlierIndices,outlierIndices] = pcfitplane(___)    将线性指数返回到 点云输入。

[___,meanError] = pcfitplane(___)此外，将内在点的距离的平均误差返回给 模型，使用上述任何语法。

[___] = pcfitplane(___,Name=Value)使用一个或多个名称-值参数以及任何 来自先前语法的参数的组合。

%---------------------点云投影平面----------------------
%  曲面拟合表达式：Ax+By+Cz+D=0
model=pcfitplane(data,0.2,[0,0,1],5);   %利用pcfitplane求解曲面拟合参数，返回值为4个
A=model.Parameters(:,1);
B=model.Parameters(:,2);
C=model.Parameters(:,3);
D=model.Parameters(:,4);
model_x=data.Location(:,1);
model_y=data.Location(:,2);
model_z=data.Location(:,3);
t=-(A.*model_x+B.*model_y+C.*model_z+D)/(A^2+B^2+C^2);
x_projection=A.*t+model_x;
y_projection=B.*t+model_y;
z_projection=C.*t+model_z;
pc_projection=pointCloud([x_projection,y_projection,z_projection],'Color',data.Color);
figure;
pcshow(pc_projection);
title('the Point cloud projection plane map');

六、点云差异可视化显示------pcshowpair

pcshowpair（ptCloudA，ptCloudB） create 描述两个输入点之间差异的可视化 云。使用洋红色混合显示差异 对于点云 A，绿色表示点云 B。

pcshowpair（ptCloudA，ptCloudB，Name=Value） 除了任何 来自先前语法的参数的组合。

将图形背景色设置为黄色。pchowpair(ptCloudA,ptCloudB,BackgroundColor="yellow")

ax = pcshowpair(___)返回 绘图轴到差异的可视化

%---------------------点云差异可视化显示----------------------
shrink_z=data.Location(:,3).*0.5;
pc_shrink=pointCloud([data.Location(:,1),data.Location(:,2),shrink_z(:)],"Color",data.Color);
figure;
pcshowpair(data,pc_shrink,"MarkerSize",1);
title('the difference in two-point cloud visualization map');
hold off;

七、查询点云点半径中的邻近区域点------findNeighborsInRadius

[index，dists] = findNeighborsInRadius（ptCloud，point，radius） 返回查询点半径内的邻居 输入点云。 
可以是无组织的点，也可以是有组织的点 云。查询点半径内的邻居使用 Kd 树计算 基于搜索算法。indicesptCloud

[index，dists] = findNeighborsInRadius（ptCloud，point，radius，camMatrix） 返回输入点云中查询点半径内的邻居。
输入 点云是由深度相机生成的有组织的点云。邻居 查询点的半径使用快速近似邻居搜索算法确定。该函数使用相机投影矩阵
来了解 相邻点之间的关系，从而加快了搜索速度。然而， 与基于 Kd 树的方法相比，结果的准确性较低

[indices,dists] = findNeighborsInRadius(___,Name,Value)除了输入之外，还使用一个或多个名称-值对参数指定选项 上述语法中的参数

%---------------------查询点云点半径中的邻近区域点----------------------
n=randi(size(data.Location,1));
centre_point=data.Location(n,:);  %随机抽取点云中的点
r=5;
[indices,dists]=findNeighborsInRadius(data,centre_point,r); %输入搜索半径
neighbor_xyz=data.Location(indices,:);
dd=[1,0,0];
point_color=repmat(dd,length(indices),1);
neighbor_pc=pointCloud([neighbor_xyz(:,1),neighbor_xyz(:,2),neighbor_xyz(:,3)],"Color",point_color);
figure;
pcshow(neighbor_pc);
hold on;
centre_pc=pointCloud([centre_point(:,1),centre_point(:,2),centre_point(:,3)],'Color',[1,1,0]);
pcshow(centre_pc,"MarkerSize",200);
%---------------------
n = [1, 1, 2];
phi = linspace(0,2*pi,360);
theta = atan(-(n(1)*cos(phi)+n(2)*sin(phi))/n(3));
x_r = (centre_point(:,1)+r*cos(theta).*cos(phi))';
y_r = (centre_point(:,2)+r*cos(theta).*sin(phi))';
z_r = (centre_point(:,3)+r*sin(theta))';
range_xyz=[x_r,y_r,z_r];
dd2=[1,0,1];
range_color=repmat(dd2,length(phi),1);
range_pc=pointCloud([range_xyz(:,1),range_xyz(:,2),range_xyz(:,3)],"Color",range_color);
pcshow(range_pc);
title('search for point clouds by radius map');
hold off;

八、根据欧氏距离将点云分割成集群------pcsegdist

labels = pcsegdist（ptCloud，minDistance）将点云分割成簇，不同聚类的点之间的最小欧氏距离。 将整数聚类标签分配给 点云，并返回所有点。minDistancepcsegdistlabels

[labels，numClusters] = pcsegdist（ptCloud，minDistance） 也返回簇的数量。

[___] = pcsegdist(___,Name=Value)使用名称-值参数设置属性。

将每个聚类中的最小和最大点数设置为 。labels = pcsegdist(ptCloud,minDistance,NumClusterPoints=[1,Inf])[1,Inf]

%---------------------根据欧氏距离将点云分割成集群----------------
maxDistance = 40;
referenceVector = [0 0 1];
[~,inliers,outliers] = pcfitplane(data,maxDistance,referenceVector);
ptCloudWithoutGround = select(data,outliers);
minDistance = 20;
[labels,numClusters] = pcsegdist(ptCloudWithoutGround,minDistance);
idxValidPoints = find(labels);
labelColorIndex = labels(idxValidPoints);
segmentedPtCloud = select(ptCloudWithoutGround,idxValidPoints);
figure;
colormap(hsv(numClusters));
pcshow(segmentedPtCloud.Location,labelColorIndex);
title('Point Cloud Clusters');

九、绘制指定半径的球型点云并保存为.ply点云文件------sphere/pcwrite

[X,Y,Z] = sphere(100);
R=100;
loc = [X(:),Y(:),Z(:)].*R;
ptCloud = pointCloud(loc);
pcshow(ptCloud);
title('Point Cloud');
%------------------保存为ply格式--------------
pcwrite(ptCloud,'cloud_copy.ply');