python点云根据高程渲染颜色 点云提取高程

目标：

• 熟悉并掌握CloudCompare和matlab的基本命令。通过
  CloudCompare软件来显示点云数据，其中包括使用不同的分割法来显示点云、查看点云数据模型的不同统计量等，熟悉并掌握点云数据常用的存储格式.xyz，.ply和.las。
• 在此基础上，利用matlab平台的混合编程，通过不同的地物的特征，来提取目标点云数据。本次实习主要是根据道路的特征来提取点云数据中道路信息。

实验数据：

• 本实验是基于LAS点云数据进行道路地面特征点的提取。从本质上来说,LAS格式是一种二进制文件格式。其目的是提供一种开放的格式标准,允许不同的硬件和软件提供商输出可互操作的统一格式。现在LAS格式文件已成为LiDAR数据的工业标准格式。
• LAS文件按每条扫描线排列方式存放数据,包括激光点的三维坐标、多次回波信息、强度信息、扫描角度、分类信息、飞行航带信息、飞行姿态信息、项目信息、GPS信息、数据点颜色信息等。
• LAS格式定义中用到的数据类型遵循1999ANSI(AmericanNationalStandardsInstitute,美国国家标准化协会)C语言标准。一个符合LAS标准的LIDAR文件分为三个部分:公用文件头块(PUBLICHEADERBLOCK)、变量长度记(VARIABLELENGTHRECORDS)和点数据记录(POINTDATARECORD)。

数据参数介绍：

点云数据道路特征与原理：

• 在机载LiDAR点云数据中，道路点云与其他地物相比具有以下特征：①道路点属于地面点，高程低于 周围其他地物，如建筑物、植被等。利用该特征进行点 云滤波可区分出道路点与其他非地面点；②道路较周围其他地物平坦，可计算各点的法向量，设置合适的坡度阈值从地面点中提取候选道路点：③在几何形状上 道路呈长条带状分布且彼此连通。
• 因此，道路点集彼此相互连通，与离散点相比进行欧式聚类后，点数越多的类别属于道路的可能性越大。根据上述道路点特征，提出了一种如图l所示的道路提取方法。

提取方法

详细方法说明：

• 本实习主要利用matlab软件在研究区域的机载点云数据中提取出道路的部分的点云数据并在CloudCompare软件中显示出结果。LAS点云数据按每条扫描线排列方式存放数据,包括激光点的三维坐标、多次回波信息、强度信息、扫描角度、分类信息、飞行航带信息、飞行姿态信息、项目信息、GPS信息、数据点颜色信息等。
• 这里只用到激光点的三维坐标，以数据点为球心，创建一定半径的球体搜索范围内的点，求得搜索范围内的最小高程。最终得到整体研究区域的相对高程数据。
• 在CloudCompare软件中通过实验获得道路的高程阙值，得到道路的点云数据，剔除非地面点云。利用聚类的方法将之前获得的道路点云数据中干扰点排除，进一步优化结果数据。

下面是操作方法：

在CloudCompare中显示原始数据：

• 导入原始的las数据，加载数据时只保留三维坐标X、Y、Z，并保存为txt格式数据导出。

• 将刚刚保存到txt格式数据加载到CloudCompare软件中，在弹出的对话框中将Z数据改为Scala格式赋予颜色。

设置高程阈值

• 在MATLAB中，编写如下代码获得研究区域的相对高程数据，保存为txt格式。
  
• 在CloudCompare软件中按照2.2.1的操作显示相对高程的点云数据，并将相对高程值赋色。
  
• 手动调节相对高程的阈值。比较多次，确定阈值为多少时才能更好的使道路突显出来，此时得到阈值大概为3.6232275时效果最佳。
  

提取地面特征点数据：

• 在MATLAB中，高程阈值设置为4.1232275编写如下代码获得研究区域的地面特征点数据，保存为txt格式。
• clc; clear all; tic;
  [s,h,v]=lasread(‘F:\ARRACA_SanFranCoast_2010_10SEG5082.las’,‘xyzi’);
  ReferencePts(:,1)=s.X; ReferencePts(:,2)=s.Y; ReferencePts(:,3)=s.Z;
  XY=ReferencePts(:,1:2); XYZ=ReferencePts(:,1:3); Z=ReferencePts(:,3);
  road_d=1; f=fopen(‘F:\c.txt’,‘w’); TreeXY= kdtree_build(XY);
  TreeXYZ= KDTree(XYZ); num_points=length(ReferencePts);
  RH =zeros(num_points,1); for num=1:num_points
  pointxyz=ReferencePts(num,1:3);
  [IDX, DIST] = ball(TreeXYZ,pointxyz,2); Hb=ReferencePts(IDX,:);
  rect= XY(num,:);
  switch road_d;
  case road_d == 1;
  range1 = [rect(1)-(+inf),rect(1)+(+inf);rect(2),rect(2)];
  range_i1 = kdtree_range_query( TreeXY, range1 );
  find1 = size(range_i1,1);
  for Num_range1=1:find1;
  Zr_1(Num_range1,:) = Z(range_i1(Num_range1)😅
  end
  sort_Zr_1 = sort(Zr_1);
  Z_min = sort_Zr_1(1,:);
  RH(num,:) = Z(num,:) - Z_min;
  otherwise
  range2 = [rect(1),rect(1);rect(2)-(+inf),rect(2)+(+inf)];
  range_i2 = kdtree_range_query( TreeXY, range2 );
  find2 = size(range_i2,1);
  for Num_range2=1:find2;
  Zr_2(Num_range2,:) = Z(range_i2(Num_range2)😅;
  end
  sort_Zr_2 = sort(Zr_2);
  Zr_min = sort_Zr_2(1,:);
  RH(num,:) = Z(num,:) - Zr_min;
  end
  if (RH(num,:)<=4.1232275 &&RH(num,:)>=0)
  fprintf(f,’%.5f %.5f %.3f %f \r\n’,
  pointxyz(1),pointxyz(2),pointxyz(3),RH(num,:));
  end end toc;

• 提取的结果中包含一些非地面点的数据，可能是低矮的树等一些地物，采用聚类的方法去除这些碎步点干扰。聚类后提取道路的效果比之前的较好，减少了局部地物的干扰。