python-pclopen3d读取显示pcdbin等格式点云数据
Posted weixin_46124984
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了python-pclopen3d读取显示pcdbin等格式点云数据相关的知识，希望对你有一定的参考价值。
第二章 python-pcl、open3d读取、显示pcd、bin格式点云数据
文章目录
前言
点云数据实际上就是许多组点的集合，每个点由x,y,z组成。当然理论上的只包含有3D坐标。
实际激光雷达获取的点云数据还会包含强度、反射率等等。但我们一般只用提取x,y,z来处理即可。
点云数据相比于其他传感器数据的核心优势就是在于 精准的深度信息。可惜获取具体的坐标信息。
环境
因为涉及到深度学习的应用，所以运用的整体环境是基于 python 使用pycharm编辑。
主要涉及库 python-pcl open3D mayavi numpy 等等用来处理点云数据并显示。
实际上点云的核心处理库应该是PCL 类似于图像处理中的Opencv 。后面也会有扩展。
一、点云数据类型
基于激光雷达录制的点云文件具有多种格式，如pcd、npy、ply、bin、pcap、lvx等等。
而读取点云数据也有多个库，python-pcl open3D mayavi numpy。
这里主要使用 python-pcl 和 open3D 两个核心库。
1.基于python-pcl 读取显示pcd、bin格式文件
代码如下（示例）：
import numpy as np
import pcl.pcl_visualization
pt = pcl.load("D://code-python//Data//lidar//2094.799809520.pcd")
#转为数组 形如 24000 x 3
points = pt.to_array()
#但是pcl显示是要 N*4 所以要扩展一列 取一列插入数组使 N*3 变为 N*4
x = points[:,0]
points = np.insert(points,3,x,axis=1)
# 这里对第四列进行赋值，它代表颜色值，根据你自己的需要赋值即可；
points[:, 3] = 255
# PointCloud_PointXYZRGB 需要点云数据是N*4，分别表示x,y,z,RGB ,其中RGB 用一个整数表示颜色；
color_cloud = pcl.PointCloud_PointXYZRGB(points)
visual = pcl.pcl_visualization.CloudViewing()
#窗口名
visual.ShowColorCloud(color_cloud, b'sta')
flag = True
while flag:
flag != visual.WasStopped()
这里就可以成功显示 颜色是根据自己设置
下面基于 numpy 读取bin格式点云文件
import pcl.pcl_visualization
import numpy as np
# lidar_path 指定一个kitti 数据的点云bin文件就行了
#bin文件为2进制文件
lidar_path = r'D:/code-python/Data/lidar/000000.bin'
# reshape成 N*4
points = np.fromfile(lidar_path, dtype=np.float32).reshape(-1, 4)
# 在这里对第四列进行赋值，它代表颜色值，根据你自己的需要赋值即可；
points[:,3] = 255
# PointCloud_PointXYZRGB 需要点云数据是N*4，分别表示x,y,z,RGB ,其中RGB 用一个整数表示颜色；
color_cloud = pcl.PointCloud_PointXYZRGB(points)
visual = pcl.pcl_visualization.CloudViewing()
visual.ShowColorCloud(color_cloud, b'cloud')
flag = True
while flag:
flag != visual.WasStopped()
2.基于open3d 读取显示pcd格式文件
代码如下（示例）：
import open3d as o3d
import numpy as np
point = o3d.io.read_point_cloud("D:/code-python/Data/lidar/000000.pcd")
o3d.visualization.draw_geometries([point])
#open3d显示的时候有个bug 要把python pycharm设置为高性能 系统显示里 图形设置
效果相对于 python-pcl的显示还是差一点，不过也挺不错。同时是可以调的，它也是一个成熟的3D处理库。有许多内置的处理函数。
3.解析pcap格式点云文件并通过python-pcl显示
这里参考一位大佬的代码：
作者：lonlon ago
https://zhuanlan.zhihu.com/p/158621756
代码如下：
# -*- coding: UTF-8 -*-
import dpkt
import collections # 有序字典需要的模块
import time
import numpy as np
import struct
# 安装了pcl , 可以使用它来进行可视化
import pcl.pcl_visualization
viewer = pcl.pcl_visualization.PCLVisualizering()#初始化一个对象
viewer.SetBackgroundColor(0, 0, 0) #颜色
viewer.AddCoordinateSystem()
viewer.InitCameraParameters()
# vlp 16 的参数
'''
https://blog.csdn.net/qq_34911636/article/details/89946329#commentBox
激光雷达每一帧的数据长度固定为1248字节，其中分别为前42字节的前数据包标识、12组数据包、4字节时间戳和最后两字节雷达型号参数。
12组数据包中前两字节为数据包的开始标识（0xFFEE）、接下去两字节为的旋转角度（当前角度）值和连续32*（2字节的距离值+1字节的激光反射强度值）字节的距离信息，
其中32*3字节分别为雷达两次获取探测信息，每个数据包开头所携带的旋转角度是指当前数据包前16*3字节对应的角度，而后16*3字节对应的旋转角度激光雷达没有直接给出，
需要通过计算前后两次旋转角度然后求取平均值获得。
1248 = 42 + 12*(2 + 2 + 32*(2+1)) + 4 + 2 =1248
雷达扫描频率为10Hz，每秒数据包在480帧左右，即每次扫描会产生48个左右的数据包，需要将分散的数据包数据合并称为一次扫描的点云数据
75个udp包产生一圈数据 vlp格式解析那篇文章 下面评论没理解 为什么是75 75*384*10 = 288000
若按照这样理解 角分辨率按照0.1度 360/0.1 = 3600 一圈转3600次 一次16个点 = 57600 10hz 1s 10圈 即 57600*10 = 576000 反正
有点乱
'''
DISTANCE_RESOLUTION = 0.002 # 距离数值分辨率 2mm转换为单位米
udp_package_num = 1
line_per_udp = 12 # 每个UDP 有多少列
point_per_udp_line = 32 # 每个UDP 的每列包含有多少个点
point_num_per_udp = point_per_udp_line * line_per_udp # 32*12=384
thetas_lines = [-15, 1, -13, 3, -11, 5, -9, 7, -7, 9, -5, 11, -3, 13, -1, 15] #垂直角度w代表值
thetas_point = thetas_lines * 2 * line_per_udp * udp_package_num #感觉是嵌套列表 列表乘以一个数字 [[x],[x],[x]...]
thetas_point = np.radians(thetas_point) #角度从度转为弧度 thetas_point为输入的角度 它返回一个数组, 其中包含输入数组中给定度数的等效弧度角。
thetas_point_cos = np.cos(thetas_point) #cos弧度
thetas_point_sin = np.sin(thetas_point) #sin弧度
data_fmt = '<' + (('H' + 'H' + 'HB' * point_per_udp_line) * line_per_udp + 'IH') * udp_package_num
base_range = np.array(range(2, point_per_udp_line*2+1, 2)) # 32， 距离值的基础索引
#range (start,stop,step) range(2,65,2) 32个数
angle_base_range = np.array([1])
d_range = []
r_range = []
angle_range = []
k = 0
data_gap = 2 + 2*point_per_udp_line # 每一列的长度 其实是2字节标识 2字节旋转角度 32*(2字节距离,1字节反射强度)点 这里貌似只计算旋转角度+距离66
for i in range(udp_package_num):
for j in range(line_per_udp):
d_range.append(base_range + k * data_gap + i * 2) # 66 是 HH + HB*32 d.range 列表增加,多个列表嵌套
r_range.append(base_range + k * data_gap + i * 2 + 1)
angle_range.append(angle_base_range + k * data_gap + i * 2)
k += 1
d_range = np.hstack(d_range) # 多个array组成的列表
r_range = np.hstack(r_range)
angle_range = np.hstack(angle_range)
# 水平角度插值，如果有角度跳变会怎么样？ 针对从360跳变到20的这部分拟合的并不是很好，误差很大；已经改正
x_index = np.arange(point_num_per_udp)
xp_index = np.arange(0, point_num_per_udp, point_per_udp_line) # array([ 0, 32, 64, 96, 128, 160, 192, 224, 256, 288, 320, 352])
def unpack_udp(data):
data_tuple = struct.unpack(data_fmt, data) # 原始格式是元祖，要转array，元祖不能索引
data_unpack = np.array(data_tuple, dtype=np.int64) # np.array会多耗时15毫秒 todo 这里为什么会报错？？ 可能是时间戳的数值太大了
distances = data_unpack[d_range] # 115200
refs = data_unpack[r_range] / 255
angles = data_unpack[angle_range]
angles = np.radians(angles / 100).astype(np.float32) # 除以100再弧度值
# 第一种处理角度的方式
# angles = np.tile(angles, (32, 1)).flatten('F') # 因为angle只有1个，数据有32个，需要复制32次
# 第二种方式
angles_interp = np.interp(x_index, xp_index, angles).astype(np.float32)
if angles[0] > angles[-1]: # 出现了角度的转折点
# replace_angle = np.linspace(0,20,32) # 针对从360跳变到20 的角度替换
change_index = np.argmax(angles)
replace_index = change_index * 32 + 1
interp_num_2 = int(angles[change_index+1]*32/40) # 每个UDP数据包之间的角度间隔为 40，每个包有32条线；
interp_num_1 = 32 - interp_num_2
replace_angle_1 = np.linspace(angles[change_index], 35999, interp_num_1) # 针对从360跳变到 20 的角度替换
replace_angle_2 = np.linspace(0, angles[change_index+1], interp_num_2) # 针对从360跳变到 20 的角度替换
angles_interp[replace_index:(replace_index+interp_num_1)] = replace_angle_1
angles_interp[(replace_index+interp_num_1):(replace_index+32)] = replace_angle_2
distances = distances * DISTANCE_RESOLUTION
x = distances * thetas_point_cos * np.sin(angles_interp)
y = distances * thetas_point_cos * np.cos(angles_interp)
z = distances * thetas_point_sin
raw_points = np.stack((x, y, z), axis=1).astype(np.float32)
# raw_points = np.stack((distances, angles_interp, refs, ), axis=1) # 也可以只要原始数据
print(type(x))
print(len(x)) #384一组 每个udp384一组
#print(x)
return raw_points
def main(file_path):
# f = open(file_path) # 此写法为python2之下，
f = open(file_path, mode='rb') #python3
try:
pcap = dpkt.pcap.Reader(f) # 先按.pcap格式解析，若解析不了，则按pcapng格式解析
except:
print("it is not pcap ... format, pcapng format...")
pcap = dpkt.pcapng.Reader(f)
# 接下来就可以对pcap做进一步解析了，记住在使用结束后最好使用f.close()关掉打开的文件，虽然程序运行结束后，
# 系统会自己关掉，但是养成好习惯是必不可少的。当前变量pcap中是按照“间戳：单包”的格式存储着各个单包
# 将时间戳和包数据分开，一层一层解析，其中ts是时间戳，buf存放对应的包
all_pcap_data = collections.OrderedDict() # 有序字典
# all_pcap_data_hex = collections.OrderedDict() # 有序字典,存十六进制形式
cir_point = []
i = 1
for (ts, buf) in pcap:
try:
eth = dpkt.ethernet.Ethernet(buf) # 解包，物理层
if not isinstance(eth.data, dpkt.ip.IP): # 解包，网络层，判断网络层是否存在，
continue
ip = eth.data
# if not isinstance(ip.data, dpkt.tcp.TCP): # 解包，判断传输层协议是否是TCP，即当你只需要TCP时，可用来过滤
# continue
if not isinstance(ip.data, dpkt.udp.UDP):#解包，判断传输层协议是否是UDP
continue
transf_data = ip.data # 传输层负载数据，基本上分析流量的人都是分析这部分数据，即应用层负载流量
if not len(transf_data.data): # 如果应用层负载长度为0，即该包为单纯的tcp包，没有负载，则丢弃
continue
if len( transf_data.data) != 1206: # 长度过滤 todo 为什么会有512字节的数据
continue
all_pcap_data[ts] = transf_data.data # 将时间戳与应用层负载按字典形式有序放入字典中，方便后续分析.
points = unpack_udp(transf_data.data)
if i % 76 != 0: # vlp16 每75个UDP数据包形成一圈数据
cir_point.append(points)
else:
cir_udp = np.vstack(cir_point)
print(cir_udp.shape) # 最后需要的一圈完整的点云数据 只包含了28000余个xyz坐标
cloud_all = pcl.PointCloud(cir_udp[:, 0:3].astype(np.float32)) # 可视化
viewer.AddPointCloud(cloud_all)
viewer.SpinOnce(100)
viewer.RemoveAllPointClouds(0)
cir_point = []
i += 1
except Exception as err:
print( "[error] %s" % err)
f.close()
if __name__ == '__main__':
#file_path="D:xxxxxx.pcap"
file_path = "2020-10-21-10-13-57_Velodyne-VLP-16-Data.pcap"
main(file_path)
这里其实都是数据格式解析，要知道怎么样的数据格式，用什么方式能够更好的解析。
当然好的工具使用起来也是很舒服，感谢开发这些库的大佬们。
总结
本文仅仅简单介绍了基于python的不同类型点云数据读取、显示方法，这些都比较简单。
ubuntu下一般都是录制bag 文件，并通过rviz显示，后续可以扩展一下。
从 *.FBX 文件读取和显示动画
【中文标题】从 *.FBX 文件读取和显示动画【英文标题】：Read and display animations from *.FBX file 【发布时间】：2016-12-21 10:54:55 【问题描述】：我想从fbx 文件中读取 3d 模型并将其显示在 iPhone 上的 openGL es 2.0 引擎中（不使用 Unity），并且我还想显示读取 3d 对象的动画。
如何从fbx 文件中获取动画？
据我所知，目前我能够获取姿势名称列表，其中包含转换矩阵、图层、堆栈、图层和许多曲线中姿势的完整列表。
如何组合所有信息以显示正确的动画？
我也尝试解析TakeInfo中的一些信息，但结果对我来说有点奇怪，例如：
FbxTakeInfo *ltakeInfo = pScene->GetTakeInfo(lAnimStack->GetName());
FbxTime start = ltakeInfo->mLocalTimeSpan.GetStart();
FbxTime end = ltakeInfo->mLocalTimeSpan.GetStop();
self.startTime = start.GetSecondDouble();
self.endTime = end.GetSecondDouble();
这里我得到了每个解析层的start = 0 和end = 0.014，所以我猜这是不正确的（我要显示的fbx 文件包含1 个网格和简单的5 秒持续时间动画）。
更新
经过几个小时的调查，我得到了以下信息：
作为参考，这是我要显示的测试 obj 的结构：
在这里你可以看到很多骨头（更具体 - 19） 我能够（如上所述）在列表中获得 19 个动画 obj（如骨骼/obj 的名称）和 19 组曲线，每组 151 帧（帧速率为 30 正好 5 秒的动画。- 30*5=150 + 1 个最后一个单位矩阵）。
如果我尝试将每个曲线组一个一个地用于我的网格（我只能解析 1 个网格），我会看到应用于所有网格的网格不同部分的动画（例如垂直旋转或水平平移），所以我认为每组中的这条曲线应该完全应用于特定的骨骼，因此我将为我的网格获得动画。问题是我不知道如何只为选定的骨骼顶点部分设置动画。
现在的问题是如何将所有动画分成特定于骨骼的单独组到整个 obj（因为我只有一个网格）？
如何从包含所有曲线组的列表中为每个帧获取 1 个全局曲线列表？
更新2
感谢@codetiger 的建议，我遵循评论中提供的链接中的说明，使用该技术，我能够检索具有开始和结束时间以及所需变换的骨骼特定垫的列表，但这与我之前得到的几乎相同与曲线 - 与曲线的唯一区别是我应该从 9 条曲线创建垫子（平移/缩放/旋转 xyz）而不是使用完整矩阵，但问题仍然存在 - 我如何将它们组合到 1 个全局矩阵中？
我使用的代码（找到几个链接）：
FbxNode* modelNode = _fbxScene->GetRootNode();
FbxAMatrix geometryTransform = GetGeometryTransformation(modelNode);
for (unsigned int deformerIndex = 0; deformerIndex < numOfDeformers; ++deformerIndex)
FbxSkin* currSkin = reinterpret_cast<FbxSkin*>(mesh->GetDeformer(deformerIndex, FbxDeformer::eSkin));
if (!currSkin)
continue;
unsigned int numOfClusters = currSkin->GetClusterCount();
for (unsigned int clusterIndex = 0; clusterIndex < numOfClusters; ++clusterIndex)
FbxCluster* currCluster = currSkin->GetCluster(clusterIndex);
std::string currJointName = currCluster->GetLink()->GetName();
FbxAMatrix transformMatrix;
FbxAMatrix transformLinkMatrix;
FbxAMatrix globalBindposeInverseMatrix;
currCluster->GetTransformMatrix(transformMatrix);
currCluster->GetTransformLinkMatrix(transformLinkMatrix);
globalBindposeInverseMatrix = transformLinkMatrix.Inverse() * transformMatrix * geometryTransform;
FbxAnimStack* currAnimStack = _fbxScene->GetSrcObject<FbxAnimStack>(0);
FbxString animStackName = currAnimStack->GetName();
char *mAnimationName = animStackName.Buffer();
FbxTakeInfo* takeInfo = _fbxScene->GetTakeInfo(animStackName);
FbxTime start = takeInfo->mLocalTimeSpan.GetStart();
FbxTime end = takeInfo->mLocalTimeSpan.GetStop();
FbxLongLong mAnimationLength = end.GetFrameCount(FbxTime::eFrames24) - start.GetFrameCount(FbxTime::eFrames24) + 1;
for (FbxLongLong i = start.GetFrameCount(FbxTime::eFrames24); i <= end.GetFrameCount(FbxTime::eFrames24); ++i)
FbxTime currTime;
currTime.SetFrame(i, FbxTime::eFrames24);
FbxAMatrix currentTransformOffset = modelNode->EvaluateGlobalTransform(currTime) * geometryTransform;
FbxAMatrix mat = currentTransformOffset.Inverse() * currCluster->GetLink()->EvaluateGlobalTransform(currTime);
在这里我收到 121 矩阵而不是 151 - 但是某些矩阵变换的持续时间超过 1 帧绘制的持续时间，所以我认为这里的 q-ty 也是正确的
float duration = end.GetSecondDouble() - start.GetSecondDouble(); //return 5 as and expected
我猜Autodesk SDK 有一种方法可以让每个 drawCall 获得 1 个全局变换
有什么建议吗？这可能吗？
为任何能够描述如何使用 Autodesk SDK 在 openGLES 2.0 中在 iPhone 上显示动画的人添加奖励...（抱歉，错字而不是 Facebook）
我能得到什么
blender中的原始对象
如果我分别绘制骨骼（相同的 VBO 具有不同的变换，并且每个骨骼只有索引）
【问题讨论】：
请查看此链接gamedev.stackexchange.com/questions/59419/… 我看到了这个链接 - 尝试使用描述的方法获取垫子，但我收到的所有东西 - 151 个垫子（fps 30 意味着 5 秒 - 正确） - 是单位矩阵......但应该是我了解了具有特定开始和结束时间的每个网格的全局变换垫 我还尝试在来自 Autodesk 的 Blender 和 FBX Review 应用程序中播放文件 - 一切正常......堆栈和图层的使用怎么样 - 你能建议如何处理这些项目吗？跨度> 对不起，我的知识仅限于从FBX读取骨骼和其他信息，对动画数据没有太多经验。让我们等待有人回应。 你从文件中读取动画信息了吗？赏金是给你关于如何在 iOS 上渲染动画的信息？ 【参考方案1】：这是在 OpenGL ES 中渲染动画网格的方法。这将为您提供有关您需要从文件中读取哪些详细信息的票价信息。
方法 1：（GPU 蒙皮）
此方法仅适用于基于硬件功能的有限数量的骨骼。通常这取决于您可以发送到着色器的矩阵数量。
使用 BindAttributes 将网格信息绑定到 GPU 一次，并将矩阵发送到制服中的着色器。
第 1 步：读取所有骨骼矩阵并创建一个矩阵数组，并将此数据以制服的形式发送到着色器。
第2步：在顶点着色器中，像下面的着色器一样计算gl_position
attribute vec3 vertPosition;
attribute vec3 vertNormal;
attribute vec2 texCoord1;
attribute vec3 vertTangent;
attribute vec3 vertBitangent;
attribute vec4 optionalData1;
attribute vec4 optionalData2;
uniform mat4 mvp, jointTransforms[jointsSize];
void main()
mat4 finalMatrix;
int jointId = int(optionalData1.x);
if(jointId > 0)
finalMatrix = jointTransforms[jointId - 1] * optionalData2.x;
jointId = int(optionalData1.y);
if(jointId > 0)
finalMatrix = finalMatrix + (jointTransforms[jointId - 1] * optionalData2.y);
jointId = int( optionalData1.z);
if(jointId > 0)
finalMatrix = finalMatrix + (jointTransforms[jointId - 1] * optionalData2.z);
jointId = int( optionalData1.w);
if(jointId > 0)
finalMatrix = finalMatrix + (jointTransforms[jointId - 1] * optionalData2.w);
gl_Position = mvp * finalMatrix * vec4(vertPosition, 1.0);
在这个着色器中，我在 optionalData1 和 optionalData2 属性中发送了权重绘制信息。数据是这样打包的：(BoneId1, BoneId2, BoneId3, BoneId4) 和 (Weight4Bone1, Weight4Bone2, Weight4Bone3, Weight4Bone4)。因此，在这种情况下，您最多只能有 4 个骨骼影响每个属性。片段着色器部分是常用的。
方法 2：（CPU 剥皮）
如果您无法忍受 GPU 蒙皮的限制，那么唯一的方法就是在 CPU 端进行计算。
第 1 步：通过将属于影响顶点的骨骼的矩阵相乘来计算当前帧中的顶点位置。
第 2 步：收集当前帧的新位置，并在 Vertex Attributes 中将信息发送到 GPU。
第 3 步：重新计算每一帧中的新顶点位置并更新属性缓冲区。
这种方法将计算的负载带到了 CPU 上。
【讨论】：
我想我遵循第二种方法，但不是通过直接获取矩阵，而是通过从图层获取动画曲线并从该曲线创建（曲线就像 tx ty tz sx sy sz rx ry rz）矩阵进行变换。 .. 当我尝试使用来自GetTransformMatrix 和其他函数的集群中的度量时 - 在屏幕上我看到一些奇怪的东西......
我接受你的回答 - 因为这是一个很好的开始，但老实说，我希望任何人都面临与我相同的问题，并且可以更详细地建议如何使用从 fbx 中提取的骨骼制作动画文件....以上是关于python-pclopen3d读取显示pcdbin等格式点云数据的主要内容，如果未能解决你的问题，请参考以下文章