PointCloud是open3d中用于點(diǎn)云處理的類，封裝了包括幾何變換、數(shù)據(jù)濾波、聚類分割等一系列實(shí)用算法。如無特別說明，本例中所有例程均基于斯坦福兔子的點(diǎn)云模型，下載地址：斯坦福標(biāo)準(zhǔn)模型

# 此行代碼后面不再重復(fù)引入
import open3d as o3d
# 載入斯坦福兔子 rabbit.pcd文件需在當(dāng)前python工作的文件夾中
pcd = o3d.io.read_point_cloud("rabbit.pcd")	

讀取和清除點(diǎn)云

一般點(diǎn)云數(shù)據(jù)的讀取方法屬于open3d.io的內(nèi)容，但點(diǎn)云類也提供了一些生成點(diǎn)云的方法，最簡單的是創(chuàng)建一個點(diǎn)云對象后，為其點(diǎn)集進(jìn)行賦值，下例中xyz為Nx3的矩陣

import open3d as o3d
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(xyz)

如果想讓pcd恢復(fù)到?jīng)]有點(diǎn)云的狀態(tài)，可以通過pcd.clear()來清除。

此外，點(diǎn)云類還提供了兩個靜態(tài)方法，create_from_depth_image和create_from_rgbd_image，使得用戶可以通過深度圖像或rgbd圖像來生成點(diǎn)云。

由于手頭沒有深度圖像，所以暫且用常見的RGB圖像來演示一下這個功能，用下面這張圖像生成深度圖像，并繪制點(diǎn)云

import open3d as o3d

raw = o3d.io.read_image("test.jpg")
rgbd = o3d.geometry.RGBDImage.create_from_color_and_depth(raw, raw, convert_rgb_to_intensity=False)

inter = o3d.camera.PinholeCameraIntrinsic()
# 此為相機(jī)內(nèi)部參數(shù)，是我亂寫的 1920x1080為像素尺寸；(600,600)為x,y方向的焦距；(640,360)為相機(jī)位置
inter.set_intrinsics(1920, 1080, 600, 600, 640, 360)
pcd = o3d.geometry.PointCloud().create_from_rgbd_image(rgbd, inter)
o3d.visualization.draw_geometries([pcd])

繪圖結(jié)果為

感覺很有沖擊力。

點(diǎn)云的顏色可以通過paint_uniform_color來指定，后面會經(jīng)常用到這個方法。

點(diǎn)云屬性

點(diǎn)云類中封裝了一些方法，用以提取點(diǎn)云的特征

對象屬性

>>> pcd.dimension()
3           # 點(diǎn)云維度為3
>>> pcd.has_colors()
False       # 點(diǎn)云無顏色
>>> pcd.has_covariances()
False       # 對象中不含協(xié)方差矩陣
>>> pcd.has_normals()
False       # 不含法向量
>>> pcd.has_points()
True        # 點(diǎn)云中有點(diǎn)
>>> pcd.is_empty()
False       # 點(diǎn)云不為空
>>> pcd.get_geometry_type()
<Type.PointCloud: 1>    #幾何體的類型是點(diǎn)云

點(diǎn)云特征

通過下面三個函數(shù)，可以獲取點(diǎn)云對象$x,y,z$三個坐標(biāo)的最大值、最小值以及質(zhì)心。

>>> pcd.get_max_bound()
array([8.7769  , 9.210494, 4.985259])
>>> pcd.get_min_bound()
array([-6.793   , -6.222866, -7.082071])
>>> pcd.get_center()
array([-1.00036164e-06, -1.81378116e-08,  4.44904992e-07])

點(diǎn)云框線

點(diǎn)云邊框

open3d為點(diǎn)云對象提供兩種生成邊框的方法，分別是生成定向邊框的get_oriented_bounding_box和生成軸對齊邊框的get_axis_aligned_bounding_box。

obb = pcd.get_oriented_bounding_box()
obb.color = [1,0,0]
aabb = pcd.get_axis_aligned_bounding_box()
aabb.color = [0,1,0]

效果為

其中紅色為定向邊框，相當(dāng)于最高效地將斯坦福兔子圈在了框中；綠色邊框則橫平豎直，每條框線均與坐標(biāo)軸平行。

凸包

凸包的本質(zhì)也是一種點(diǎn)云框線，但其采取的框選方案要更加緊密，即選取點(diǎn)云中最外側(cè)的點(diǎn)連接成一個凸多面體。具體情況看下圖就能明白

hull, _ = pcd.compute_convex_hull()
hull_ls = o3d.geometry.LineSet.create_from_triangle_mesh(hull)
o3d.visualization.draw_geometries([pcd, hull_ls])

幾何變換

點(diǎn)云類提供了平移、縮放以及旋轉(zhuǎn)的空間變換方法，分別用到函數(shù)translate、scale以及rotate，示例如下：

import copy

# 此為原始點(diǎn)云
pcd0 = o3d.io.read_point_cloud("rabbit.pcd")

# 1號點(diǎn)云，沿x軸移動20米
pcd1 = copy.deepcopy(pcd0).translate((20, 0, 0))
pcd1.paint_uniform_color([1, 0, 0])  #設(shè)為紅色，100 表示rgb

# 2號點(diǎn)云，演示縮放
pcd2 = copy.deepcopy(pcd0).translate((40, 0, 0))
# 前一個參數(shù)為縮放比例；后一個參數(shù)為縮放后的位置，為必填項
pcd2.scale(0.5, center=pcd2.get_center())
pcd2.paint_uniform_color([0, 1, 1])

# 3號點(diǎn)云，通過歐拉角旋轉(zhuǎn)
pcd3 = copy.deepcopy(pcd0).translate((0, -20, 0))
# 采用Euler角的方法生成旋轉(zhuǎn)角，表示繞y軸旋轉(zhuǎn)90°
thEuler = pcd3.get_rotation_matrix_from_xyz((0,np.pi/2,0))
pcd3.rotate(thEuler)
pcd3.paint_uniform_color([0, 1, 0])

# 4號點(diǎn)云，通過軸角法旋轉(zhuǎn)
pcd4 = copy.deepcopy(pcd0).translate((20, -20, 0))
# 通過軸角表示法生成旋轉(zhuǎn)角，表示繞y軸旋轉(zhuǎn)60°
th = np.array([0, np.pi/3, 0]).T
thAxis = pcd4.get_rotation_matrix_from_axis_angle(th)
pcd4.rotate(thAxis)
pcd4.paint_uniform_color([0, 0, 1])

# 5號點(diǎn)云，通過四元數(shù)旋轉(zhuǎn)
pcd5 = copy.deepcopy(pcd0).translate((40, -20, 0))
# 通過四元數(shù)法生成轉(zhuǎn)角，表示繞x軸旋轉(zhuǎn)180°
quart = np.array([0,0,0,1]).T
thQuart = pcd5.get_rotation_matrix_from_quaternion(quart)
pcd5.rotate(thQuart)
pcd5.paint_uniform_color([1, 0, 1])

# 5繪圖函數(shù)可以輸入點(diǎn)云列表
pcds = [pcd0, pcd1, pcd2, pcd3, pcd4, pcd5]
o3d.visualization.draw_geometries(pcds)

效果如下

其中平移和縮放沒什么好說的，對于旋轉(zhuǎn)，需要知道歐拉角有三個分量，分別是橫滾、俯仰以及航向，代表繞x、y、z軸旋轉(zhuǎn)。

考慮一架正在飛行的飛機(jī)，以某一時刻前后為x軸、左右為y軸、上下為z軸。則航向角就是飛機(jī)前進(jìn)方向偏離的角度；俯仰角就是飛機(jī)頭尾姿態(tài)的俯仰；橫滾角描述的就是飛機(jī)翅膀的擺動。

歐拉角用于描述靜態(tài)角度是沒問題的，但用于表示旋轉(zhuǎn)時會導(dǎo)致萬向節(jié)死鎖，簡單來說就是飛機(jī)的航向角變化90°之后，其橫滾軸變成了俯仰軸，而四元數(shù)則沒有這個顧慮。

設(shè)歐拉角為$(\psi ,\theta ,?)$，則四元數(shù)可表示為

$$\mathbf{Q}=?\begin{array}{c}\cos (\psi /2)\\ 0\\ 0\\ \sin (\psi /2)\end{array}?\left[\begin{array}{cccc}\cos (\theta /2)& 0& \sin (\theta /2)& 0\end{array}\right]?\begin{array}{c}\cos (?/2)\\ \sin (?/2)\\ 0\\ 0\end{array}?$$

點(diǎn)云類中的生成旋轉(zhuǎn)矩陣的方法均為靜態(tài)方法，可在不建立對象的情況下調(diào)用。

在演示變換的過程中調(diào)用了三個，其前綴均為get_rotation_matrix_from_，結(jié)尾是axis_angle表示通過歐拉角生成旋轉(zhuǎn)矩陣；quaternion通過四元數(shù)；from_xyz則通過旋轉(zhuǎn)向量。

open3d支持通過不同順序的xyz數(shù)組創(chuàng)建旋轉(zhuǎn)矩陣，由于三個坐標(biāo)總計有6個組合，故而提供了6個靜態(tài)方法。例如，針對向量$(x,y,z)$，其創(chuàng)建旋轉(zhuǎn)矩陣的方法為get_rotation_matrix_from_xyz；對于$(y,x,z)$，只需將xyz換為yxz即可，非常方便記憶。

法線

對于不包含法線的點(diǎn)云，通過estimate_normals方法可以生成法線。

>>> pcd.estimate_normals()
>>> pcd.has_normals()
True    # 可見pcd已經(jīng)生成了法線
>>> o3d.visualization.draw_geometries([pcd], point_show_normal=True)

結(jié)果表明，把法線畫出來實(shí)在是太丑了，可謂san值狂掉，屬于精神污染了。

estimate_normals中有兩個輸入?yún)?shù)，其中search_param為k-d樹的索引參數(shù)，用以生成法線時索引鄰近點(diǎn)；fast_normal_computation為布爾型的參數(shù)，顧名思義，為True時將開啟加速，但可能會導(dǎo)致數(shù)值不穩(wěn)定。

點(diǎn)云類中對法線的操作主要是兩種，一是歸一化；二是定向。

normalize_normals可對法線進(jìn)行歸一化，并返回歸一化法線后的點(diǎn)云。

定向則有三個方法：

• orient_normals_consistent_tangent_plane：通過切平面定向
• orient_normals_to_align_with_direction：通過對其坐標(biāo)軸定向
• orient_normals_towards_camera_location：通過相機(jī)位置定向

特征計算

協(xié)方差矩陣

對于點(diǎn)云中的某點(diǎn)$r_i$，其協(xié)方差矩陣為

$$C_i=\frac{1}{N}\sum _{i=1}^N(r_i?\bar{r})?(r_i?\bar{r}{)}^T$$

點(diǎn)云對象通過estimate_covariances函數(shù)計算每個點(diǎn)的協(xié)方差矩陣。此外，compute_mean_and_covariance可計算整個點(diǎn)云及其協(xié)方差的均值。

>>> M_and_C = pcd.compute_mean_and_covariance()
>>> pprint(M_and_C)
(array([-1.00036164e-06, -1.81378116e-08,  4.44904992e-07]),
 array([[16.80007645, -5.77246083,  0.46407276],
       [-5.77246083, 17.24825378, -2.5936246 ],
       [ 0.46407276, -2.5936246 ,  7.9322634 ]]))

最鄰近距離和馬氏距離

點(diǎn)云類為單點(diǎn)提供了兩種距離屬性，分別是最鄰近距離和Mahalanobis距離，對應(yīng)的函數(shù)分別是compute_nearest_neighbor_distance和compute_mahalanobis_distance。

最鄰近距離即當(dāng)前點(diǎn)和離它最近的其他點(diǎn)的距離；Mahalanobis距離的定義則為

$$D_M(r_i)=\sqrt{(r_i?\bar{r}{)}^T{C}_i(r_i?\bar{r})}$$

其中$C_i$為該點(diǎn)對應(yīng)的協(xié)方差矩陣。

索引、采樣和濾波

點(diǎn)云索引

select_by_index(ids, invert)方法，可以通過點(diǎn)號對點(diǎn)云進(jìn)行選取。其中ids為點(diǎn)號列表，invert參數(shù)默認(rèn)為False，表示返回ids中的點(diǎn)；若將invert置為True，將返回ids之外的點(diǎn)。

import open3d as o3d
import numpy as np
pcd = o3d.io.read_point_cloud("rabbit.pcd")
idx = np.arange(10000)
# 索引對應(yīng)的點(diǎn)
pIn = pcd.select_by_index(idx)
pIn.paint_uniform_color([1, 0, 0])
# 索引外的點(diǎn)云
pOut = pcd.select_by_index(idx, invert=True)
pOut.paint_uniform_color([0, 1, 0])
o3d.visualization.draw_geometries([pIn, pOut])

效果為

無效點(diǎn)剔除

remove_non_finite_points方法用于剔除NaN或者Inf之類的無窮值。其兩個bool型的輸入?yún)?shù)remove_nan和remove_infinite分別用于指明這兩者。

統(tǒng)計濾波和鄰域濾波

這兩中濾波方法都是先得到符合要求的點(diǎn)索引，然后通過索引濾波，將這些點(diǎn)挑選出來，輸出輸出為濾波后的點(diǎn)云和點(diǎn)的索引號。

# 上接索引濾波的內(nèi)容

pcd1 = copy.deepcopy(pcd).translate((20, 0, 0))
pcd2 = copy.deepcopy(pcd).translate((40, 0, 0))

# 統(tǒng)計濾波，參數(shù)分別表示K鄰域點(diǎn)的個數(shù)和標(biāo)準(zhǔn)差乘數(shù)
sPcd, sInd = pcd1.remove_statistical_outlier(6, 2.0)

# 半徑濾波，輸入?yún)?shù)為鄰域球內(nèi)最少點(diǎn)數(shù)和鄰域半徑
rPcd, rInd = pcd2.remove_radius_outlier(9, 0.   )

o3d.visualization.draw_geometries([sPcd, rPcd])

效果如下

這兩種算法的邏輯是一樣的，對于某點(diǎn)$x$，選取距離$x$最近的一些點(diǎn)，如果這些點(diǎn)的標(biāo)準(zhǔn)差小于設(shè)定值，則符合統(tǒng)計濾波的標(biāo)準(zhǔn)；如果均小于鄰域半徑，則符合半徑濾波的標(biāo)準(zhǔn)。

下采樣

點(diǎn)云對象共有三種下采樣方案

其中隨機(jī)下采樣即根據(jù)采樣率，隨機(jī)選擇一些點(diǎn)；等序則根據(jù)輸入的采樣間隔，則取等間隔的點(diǎn)的序號。

體素采樣稍顯復(fù)雜，會構(gòu)建三維體素格網(wǎng)，然后輸出格網(wǎng)內(nèi)的點(diǎn)云質(zhì)心，而非原始數(shù)據(jù)；若存在法線或顏色，則通通取均值。

此外，還有voxel_down_sample_and_trace函數(shù)，在體素采樣的基礎(chǔ)上，可以規(guī)定體素的邊界范圍。

downpcd = pcd.voxel_down_sample(20)

聚類算法

DBSCAN聚類

DBSCAN，即Density-Based Spatial Clustering of Applications with Noise，基于密度的噪聲應(yīng)用空間聚類。

在DBSCAN算法中，將數(shù)據(jù)點(diǎn)分為三類：

1. 核心點(diǎn)：若樣本$x_i$的$\epsilon$鄰域內(nèi)至少包含了$M$個點(diǎn)，則為核心點(diǎn)
2. 邊界點(diǎn)：若樣本$x_i$的$\epsilon$鄰域內(nèi)包含的點(diǎn)數(shù)小于$M$，但在其他核心點(diǎn)的$\epsilon$鄰域內(nèi)，則為邊界點(diǎn)
3. 噪聲：既非核心點(diǎn)也非邊界點(diǎn)則為噪聲

可見，DBSCAN算法需要兩個參數(shù)，分別是鄰域半徑$\epsilon$和點(diǎn)數(shù)$M$。

在open3d中，提供了cluster_dbscan接口，示例如下

import open3d as o3d
import numpy as np
import matplotlib.pyplot as plt

pcd = o3d.io.read_point_cloud("rabbit.pcd")

eps = 0.5       # 同一聚類中最大點(diǎn)間距
M = 50          # 有效聚類的最小點(diǎn)數(shù)
Labels = np.array(pcd.cluster_dbscan(eps, M))
print(np.max(Labels))                   # 得到結(jié)果為3
cs = plt.get_cmap("jet")(Labels/3)      # 偽彩映射
cs[labels < 0] = 0  # labels = -1 的簇為噪聲，以黑色顯示
pcd.colors = o3d.utility.Vector3dVector(cs[:, :3])
o3d.visualization.draw_geometries([pcd])

結(jié)果如圖所示

由于斯坦福兔子的掃描點(diǎn)比較連續(xù)，所以分割效果不太還，為了更好地演示DBSCAN算法，可以用手骨模型，下載地址：斯坦福標(biāo)準(zhǔn)模型

RANSAC 平面分割

RANSAC，即RANdom SAmple Consensus，隨機(jī)抽樣一致算法。

以平面上的點(diǎn)集舉例，假設(shè)點(diǎn)集中有一條直線$L$，$L$外的點(diǎn)很少，均為噪聲。

那么第一步，隨機(jī)選取兩個點(diǎn)連成一條直線$\hat{L}$，那么這條直線有可能就是$L$，也有可能是噪聲連出來的莫名其妙的一條線。

接下來，隨機(jī)抽取點(diǎn)集中的一些點(diǎn)，如果隨機(jī)抽取的大部分點(diǎn)都落在$L$附近，那么就說明$\hat{L}$有很大的概率就是$L$；否則說明不太像是$L$。隨著抽取出的直線越來越多，最后可以得到最接近$L$的直線，從而完成了對點(diǎn)集的分割。

在Open3d中，提供了基于RANSAC算法的平面分割接口segment_plane

pcd = o3d.io.read_point_cloud("rabbit.pcd")

d = 0.2             # 內(nèi)點(diǎn)到平面模型的最大距離
n = 5               # 用于擬合平面的采樣點(diǎn)數(shù)
nIter = 50          # 最大迭代次數(shù)

# 返回模型系數(shù)plane和內(nèi)點(diǎn)索引ids，并賦值
plane, ids = pcd.segment_plane(d, n, nIter)

# 平面方程
[a, b, c, d] = plane

# 平面內(nèi)點(diǎn)點(diǎn)云
iCloud = pcd.select_by_index(ids)
iCloud.paint_uniform_color([0, 0, 1.0])

# 平面外點(diǎn)點(diǎn)云
oCloud = pcd.select_by_index(ids, invert=True)
oCloud.paint_uniform_color([1.0, 0, 0])

# 可視化平面分割結(jié)果
o3d.visualization.draw_geometries([iCloud, oCloud])

最后得到的結(jié)果為

本來以為平面會出現(xiàn)在兔子的底座上，沒想到會把兔子一分為二。文章來源地址http://www.zghlxwxcb.cn/news/detail-418330.html

到了這里，關(guān)于Open3d點(diǎn)云對象詳解的文章就介紹完了。如果您還想了解更多內(nèi)容，請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！