基於Python深度圖生成3D點雲
theme: channing-cyan highlight: agate
攜手創作，共同成長！這是我參與「掘金日新計劃 · 8 月更文挑戰」的第2天，點選檢視活動詳情 >>
前言
廢話不多說，直接開造。這裡的話我們有兩個目標，第一個是如何把一個2維圖片上的點對映到3維空間。第二就是如何生成3D點雲。當然實際上這是一個大問題，因為只要解決了第一個問題，第二個問題就是“送分”
二維RGB影象
在說到3D點雲之前我們就不得不說到RGB影象，也就是一些二維影象。
影象以像RGB三個通道的形式進行儲存。也就是這樣：
平時我們看到的就是左邊的2D影象，實際上是以右邊的形式儲存的。
那麼從我們的矩陣角度來看的話，大概是這個樣子的：
```python [ [[r,g,b],[r,g,b],[r,g,b],[r,g,b]], [[r,g,b],[r,g,b],[r,g,b],[r,g,b]], ...
] ```
這裡面存在了兩組資訊， 第一組就是我們最容易忽略的位置資訊，也就是畫素點的位置。 第二組就是我們的色彩資訊。
顯然如果我們想要生成3D點雲的話，我們的色彩資訊是必要的。
而位置資訊，在二維平面上，是在3為空間上的對映
那麼要想得到3D點雲，那麼我們就必須得到完整的位置座標，也就是x,y,z
成像原理
那麼說到這裡的話，我們就不得不說到，攝像機是如何把咱們的影像映射出來的了。
大概他是這樣成像的：
這裡咱們有三個東西需要注意一下。
首先是我們的攝像機位置 之後是相片的位置也就是中間那個 最後是我們物品實際位置
所以我們實際上的一個二維影象，就是一個投影，在實際上的話，我們可以理解為深度為1的空間。 但是這個投影
所以的話，用數學矩陣的形式表示為：
其中u,v為二維平面上的任意點座標，1為深度，Zc就是咱們的二維平面上的點。
其中R,T是外參矩陣，並且我們的世界座標原點就是相機的座標原點
所以，R，T可以取下面的矩陣，（具體原理的話比較複雜，我也不是很懂為什麼會這樣，我回頭再補補數學）
帶入得到：
最後轉換得到這個公式：
那麼接下來咱們就對這個公式進行處理，來生成我們的點雲
資料準備
前面通過公式我們可以瞭解到，一個二維影象，具備了RGB資訊，還有對應的那個畫素點在二維空間的“投影”。我們需要將這個點重新還原到三維空間。所以這裡需要使用到深度資訊，而這個是需要一些專業攝像機才有的，比如iphone的，或者機器人的：
有了這個深度圖片，我們就可以還原座標了。其實說到這兒，我們其實可以發現，RGB圖片的話其實只是提供了一下色彩如果要還原為3D點雲的話，因為深度圖和對應的rgb的矩陣的位置是一一對應的。
圖片載入
說到這裡，我們就可以開始我們的正式編碼了。
首先是咱們的圖片載入階段，這裡的話主要是我們先需要一些資訊，和讀取圖片。
```python
def __init__(self, rgb_file, depth_file, save_ply, camera_intrinsics=[784.0, 779.0, 649.0, 405.0]):
self.rgb_file = rgb_file
self.depth_file = depth_file
self.save_ply = save_ply
self.rgb = cv2.imread(rgb_file)
self.depth = cv2.imread(self.depth_file, -1)
print("your depth image shape is:",self.depth.shape)
self.width = self.rgb.shape[1]
self.height = self.rgb.shape[0]
self.camera_intrinsics = camera_intrinsics
self.depth_scale = 1000
``` 下面的這個東西，就是咱們公式當中的u0,v0,dx,dy
python
camera_intrinsics=[784.0, 779.0, 649.0, 405.0]
這個要注意的是，這個玩意呢，是不同的裝置相機有不同的引數，要根據自己的設定。
一般來說相機會給這樣的矩陣：
我們對應起來就可以了。
之後是要注意的是這個深度圖是uint16的，所以讀取的時候加個-1
這裡還有個東西： 深度圖與比例因子（scale_factor） 在程式碼中是：
python
self.depth_scale = 1000
深度圖對應的尺度因子是深度圖中儲存的值與真實深度（單位為m）的比例 通常情況下，深度值以毫米為單位儲存在16位無符號整數中(0~65535)，因此要獲得以米為單位的z值，深度圖的畫素需要除以比例因子1000。不過不同相機的的比例因子可能不同，我這裡的話是1000，這個看自己的實際的。
演算法實現
之後就是咱們的演算法實現了，這個主要是位置換算，得到真實的x,y,z
```python
depth = np.asarray(self.depth, dtype=np.uint16).T
# depth[depth==65535]=0
self.Z = depth / self.depth_scale
fx, fy, cx, cy = self.camera_intrinsics
X = np.zeros((self.width, self.height))
Y = np.zeros((self.width, self.height))
for i in range(self.width):
X[i, :] = np.full(X.shape[1], i)
self.X = ((X - cx / 2) * self.Z) / fx
for i in range(self.height):
Y[:, i] = np.full(Y.shape[0], i)
self.Y = ((Y - cy / 2) * self.Z) / fy
```
之後,x,y,z就算好了。
注意的的是我們計算完成後得到的x，y , x 的每一組向量都是 寬×高
生成點雲
現在咱們已經得到了具體的座標。那麼接下來是要生成點雲的，我們要把剩下的色彩資訊搞上去。 程式碼很簡單：
python
data_ply = np.zeros((6, self.width * self.height))
data_ply[0] = self.X.T.reshape(-1)
data_ply[1] = -self.Y.T.reshape(-1)
data_ply[2] = -self.Z.T.reshape(-1)
img = np.array(self.rgb, dtype=np.uint8)
data_ply[3] = img[:, :, 0:1].reshape(-1)
data_ply[4] = img[:, :, 1:2].reshape(-1)
data_ply[5] = img[:, :, 2:3].reshape(-1)
self.data_ply = data_ply
t2 = time.time()
print('calcualte 3d point cloud Done.', t2 - t1)
之後就是儲存檔案了。 目前點雲的主要儲存格式包括：pts、LAS、PCD、.xyz 和. pcap 等
例如：
pts 點雲檔案格式是最簡便的點雲格式，直接按 XYZ 順序儲存點雲資料， 可以是整型或者浮點型。
LAS 是鐳射雷達資料（LiDAR），儲存格式比 pts 複雜，旨在提供一種開放的格式標準，允許不同的硬體和軟體提供商輸出可互操作的統一格式。LAS 格式點雲截圖，其中 C：class(所屬類)，F：flight(航線號)，T：time(GPS 時間)，I：intensity(回波強度)，R：return(第幾次回波)，N：number of return(回波次數)，A：scan angle(掃描角)，RGB：red green blue(RGB 顏色值)。
等等，我們這裡的是ply 全名為多邊形檔案（Polygon File Format）或 史丹佛三角形檔案（Stanford Triangle Format）。. 該格式主要用以儲存立體掃描結果的三維數值，透過多邊形片面的集合描述三維物體，與其他格式相較之下這是較為簡單的方法。. 它可以儲存的資訊包含顏色、透明度、表面法向量、材質座標與資料可信度，並能對多邊形的正反兩面設定不同的屬性。
格式為
頭部
頂點列表
面片列表
其他元素列表
```python start = time.time() float_formatter = lambda x: "%.4f" % x points = [] for i in self.data_ply.T: points.append("{} {} {} {} {} {} 0\n".format (float_formatter(i[0]), float_formatter(i[1]), float_formatter(i[2]), int(i[3]), int(i[4]), int(i[5])))
file = open(self.save_ply, "w")
file.write('''ply
format ascii 1.0
element vertex %d
property float x
property float y
property float z
property uchar red
property uchar green
property uchar blue
property uchar alpha
end_header
%s
''' % (len(points), "".join(points)))
file.close()
end = time.time()
print("Write into .ply file Done.", end - start)
```
點雲顯示
這個的話可以使用軟體
也可以使用open3d模組
完整程式碼
```python import cv2 import numpy as np import open3d as o3d
import time
class point_cloud_generator():
def __init__(self, rgb_file, depth_file, save_ply, camera_intrinsics=[784.0, 779.0, 649.0, 405.0]):
self.rgb_file = rgb_file
self.depth_file = depth_file
self.save_ply = save_ply
self.rgb = cv2.imread(rgb_file)
self.depth = cv2.imread(self.depth_file, -1)
print("your depth image shape is:",self.depth.shape)
self.width = self.rgb.shape[1]
self.height = self.rgb.shape[0]
self.camera_intrinsics = camera_intrinsics
self.depth_scale = 1000
def compute(self):
t1 = time.time()
depth = np.asarray(self.depth, dtype=np.uint16).T
# depth[depth==65535]=0
self.Z = depth / self.depth_scale
fx, fy, cx, cy = self.camera_intrinsics
X = np.zeros((self.width, self.height))
Y = np.zeros((self.width, self.height))
for i in range(self.width):
X[i, :] = np.full(X.shape[1], i)
self.X = ((X - cx / 2) * self.Z) / fx
for i in range(self.height):
Y[:, i] = np.full(Y.shape[0], i)
self.Y = ((Y - cy / 2) * self.Z) / fy
data_ply = np.zeros((6, self.width * self.height))
data_ply[0] = self.X.T.reshape(-1)
data_ply[1] = -self.Y.T.reshape(-1)
data_ply[2] = -self.Z.T.reshape(-1)
img = np.array(self.rgb, dtype=np.uint8)
data_ply[3] = img[:, :, 0:1].reshape(-1)
data_ply[4] = img[:, :, 1:2].reshape(-1)
data_ply[5] = img[:, :, 2:3].reshape(-1)
self.data_ply = data_ply
t2 = time.time()
print('calcualte 3d point cloud Done.', t2 - t1)
def write_ply(self):
start = time.time()
float_formatter = lambda x: "%.4f" % x
points = []
for i in self.data_ply.T:
points.append("{} {} {} {} {} {} 0\n".format
(float_formatter(i[0]), float_formatter(i[1]), float_formatter(i[2]),
int(i[3]), int(i[4]), int(i[5])))
file = open(self.save_ply, "w")
file.write('''ply
format ascii 1.0
element vertex %d
property float x
property float y
property float z
property uchar red
property uchar green
property uchar blue
property uchar alpha
end_header
%s
''' % (len(points), "".join(points)))
file.close()
end = time.time()
print("Write into .ply file Done.", end - start)
def show_point_cloud(self):
pcd = o3d.io.read_point_cloud(self.save_ply)
o3d.visualization.draw([pcd])
if name == 'main': camera_intrinsics = [378.998657, 378.639862, 321.935120, 240.766663] rgb_file = "data/1.jpg" depth_file = "data/1.png" save_ply = "data.ply" a = point_cloud_generator(rgb_file=rgb_file, depth_file=depth_file, save_ply=save_ply, camera_intrinsics=camera_intrinsics ) a.compute() a.write_ply() a.show_point_cloud()
```
效果如下：
總結
這裡的話其實還是從一張深度圖到3d點雲，實際上有時候我們可能只需一組座標然後還原，那麼這部分的話，要實現的話也不難，不過要重新做一點轉換，公式還是那個公式。這裡的話就不能在多說了，而且這裡還有一點值得一提的是這個通過這種方式得到的3d點雲其實怎麼說呢，並不是特別準，所以這邊還是要更加牛批的人工智慧演算法的，目前有個整合玩意可以玩玩的是middlepipe。這方面的話Google還是得是他。
- 還在調API寫所謂的AI“女友”，嘮了嘮了，教你基於python咱們“new”一個（深度學習）
- Java前後端分離實戰Auto2.0使用者登入註冊--基本的使用者登入 郵箱驗證
- 卡爾曼濾波器（目標跟蹤一）（上）
- 手把手教你如何自制目標檢測框架（從理論到實現）
- 基於Python深度圖生成3D點雲
- Pandas基礎使用（機器學習基礎）
- CEC2017基礎函式說明Python版本
- 全國空氣質量爬取實戰
- 智慧演算法整合測試平臺V0.1實戰開發
- DDPG神經網路實戰（基於強化學習優化粒子群演算法）
- 關於強化學習優化粒子群演算法的論文解讀（全）
- 關於強化學習優化粒子群演算法的論文解讀（上）
- 基於多種群機制的PSO演算法（優化與探索三 *混合種群思想優化多種群與廣義PSO求解JSP）
- 基於多種群機制的PSO演算法Python實現（優化與探索二）
- 基於多種群機制的PSO演算法Python實現
- 520桌面手勢告白
- 嘿~瞎話一下，為啥要用Sigmoid？！