用於雙目重建中的GPU編程：julia-cuda
作者：京東科技 李大沖
一、Julia是什麼
julia是2010年開始面世的語言，作為一個10後，Julia必然有前輩們沒有的特點。Julia被期望塑造成原生的有C++的運行速度、python的易交互性以及膠水性。最重要的是，julia也是nVidia-cuda官方選中的接口語言，這點讓cuda編程變得稍微容易一些。
二、項目背景
雙目立體視覺是比較常規獲取深度信息的算法。這是一種模仿人的雙眼的方法，根據對同一個點在不同的視角進行觀察，並在不同的觀察視角上的觀察差異推算出深度，這在圖1中更清晰一點。被動雙目的缺點是需要紋理特徵去判斷左右相機裏是否為同一個點。也就是，當沒有問題、紋理不夠明顯或者紋理重複時，這種方法難以工作。
被動雙目示意圖
觀測點在一個相機和在另外一個相機上的成像是有明顯可區分性的，二者的觀察差異（L1和L2）形成了視差，進而可形成深度信息。
被動雙目的缺點是某些場景下無效，主動雙目便出現了。如果我們投射一些人工設定編碼的圖案，根據實現約定通過解碼可以對視野內是所有點賦予唯一的標識符號，那被動雙目的問題不就解決了嗎。
事實上，人們很快發現可以使用空間上編碼和時間上編碼兩種方式。空間上編碼，就像人們觀察銀河命名星座一樣，去投射隨機生成的星星點點的光斑，沒有紋理的夜空照射進了完全可辨識的圖案。如同星座往往是一個空間範圍很大的概念，空間編碼的方式無法做到非常精細。所以，時間編碼的方式應運而生，時間編碼是投射一組圖案，讓被照射的區域，被照射的明暗變化情況是完全不一樣的。
像圖2中的二進制編碼（或行業內稱為格雷碼），再加上灰度變化更加豐富的編碼方式以及相機在豎方向上對齊之後是隻需要考慮橫軸方向上的差異這兩點，可以實現全視野的點完全可區分。
空間編碼和時間編碼
空間編碼如同銀河中的星星點點，通過組合星點變為星座，人們可以對深邃夜空不同位置依靠銀河區分開。時間編碼是投射多組圖案，讓每一個點是時間上都有唯一的變化情況。就像右圖，這16個方格都是完全可區分的。
三、效率問題
使用主動雙目的結構去重建，一個很大的影響運行效率的問題是需要對左右圖像求算視差（原理如圖1所示）。具體求算視差的操作是要對左相機裏的每一個點，去找右相機中同一行中與之距離最近的點的橫座標，然後繼續插值之後找到距離更近的插值座標，最後把這兩個座標的差異當作視差進行後續的計算。
使用時間編碼的主動雙目，按照實現約定的編碼進行解碼之後，其效果如圖3所示，這個環節就是在這對數據上進行處理。以此數據為例，有效區域大概是700*850，完全遍歷所需要的時間複雜度大概是O(700*850*850)，直接估算就很耗時。
時間編碼解碼後的圖像示意圖，左右相機拍攝的圖案編程瞭解碼值
接下來需要在這對數據上計算視差，而確定相同編碼值需要大量的遍歷，對相同編碼進一步的優化還要對每個點進行小範圍的插值處理，時間複雜度很大。
視差圖效果
1. 使用for訓練的方式
使用for循環，需要3層。經測試後，在TX2上，使用c++來完成大概需要7s多。
2. 使用numpy的矩陣操作來實現
-
代碼如下，使用numpy的廣播的原理來實現，然後通過一系列的過濾異常值後，再經過插值進一步的處理得到最終想要的結果。
-
這段代碼沒有進一步的去實現插值的過程，耗時情況為1.35s。
def disp(l_, r_, colstartR, colendR, colstartL, colendL, rowstart, rowend):
# 聚焦在有效值區域
l = l_[rowstart: rowend, colstartL: colendL]
r = r_[rowstart: rowend, colstartR: colendR]
# 算一下閾值。閾值幫助確定判斷太遠的點就不能當作是目標值
thres = (np.max(l) - np.min(l)) / (colendL - colstartL) * 10
h, w = l.shape
ll = l.T.reshape(w, h, 1)
# print(l.shape, r.shape)
# 加快運算，縮小圖片，4倍下采樣
quart_r = cv2.resize(r, (w // 4, h))
# 計算距離
distance = np.abs(ll - quart_r)
# 計算點位
idx_old = np.argmin(distance, axis=-1).T * 4
# 去除指向同一個點的重複值
_, d2 = np.unique(idx_old, return_index=True)
idx = np.zeros_like(idx_old).reshape(-1)
idx[d2] = idx_old.reshape(-1)[d2]
idx = idx.reshape(idx_old.shape)
# 計算最值
minV = np.min(distance, axis=-1).T
ori_idx = np.tile(np.arange(w).reshape(w, 1), h).T
disparity = ori_idx - idx + colstartL - colstartR
# 去除異常點
disparity[minV > thres] = 0# 大於閾值的
disparity[l < 1e-3] = 0 # 左圖點無效的
# 插值
# 方法：在最值附近的區域去-5：5的值，二次線性插值上採樣，採完之後再計算一遍上面的步驟，算完之後取小數進行計算
# 暫時沒有實現這一步
return disparity, l, r
3. 使用julia-cuda實現
這段主要是cuda核函數的部分，主要包括以下部分：
1）用到了512*700個cuda線程來實現，增加的並行計算抵消了單個點位計算耗時。
2）還用到了共享顯存，降低讀取數據的耗時，因為右圖中的每一行都會被左圖中的同行中的每一個元素加載一遍，對於元素for循環是840*840的時間空間複雜度。這兒加載一遍，就可以讓這些線程同時看到了。
3）對於上述過程中提到的插值一事，使用cuda的紋理內容輕鬆搞定，因為紋理內存的一個特性就是插值。即對於一個向量，可以取得到下標不是整數時的值，這些值就是自動插值出現的，且此過程的資源消耗非常低，因為這部分的初始設計是為了圖像渲染和可視化。
4）通過這部分代碼，實現相同的效果在TX2上的耗時是165ms，快了近10倍。
function bench_match_smem(cfg, phaL, phaR, w, h, winSize, pha_dif)
texarr2D = CuTextureArray(phaR)
tex2D = CuTexture(texarr2D; interpolation = CUDA.LinearInterpolation())
cp, minv, maxv = cfg.cpdiff, cfg.minv, cfg.maxv
colstart, colend = cfg.colstart, cfg.colend
rowstart, rowend = cfg.rowstart, cfg.rowend
mindis, maxdis = cfg.mindis, cfg.maxdis
col_ = cld((colend - colstart), 32) * 32
row_ = cld(rowend - rowstart, 32) * 32
stride = Int(cld((maxv - minv + 1), 512))
threadsPerBlock = round(Int32, cld((maxv - minv + 1) / stride, 32) * 32)
blocksPerGrid = row_
println("blocks = $blocksPerGrid threads = $threadsPerBlock left $(col_) right $(maxv - minv) h=$(row_)")
@cuda blocks = blocksPerGrid threads = threadsPerBlock shmem =
(threadsPerBlock * sizeof(Float32)) phaseMatch_smem!(cp, mindis, maxdis, minv, maxv, colstart, colend, rowstart, rowend, phaL, tex2D, threadsPerBlock,blocksPerGrid,stride, w, h, winSize, pha_dif)
CUDA.synchronize()
return
end
#進行立體相位匹配
function phaseMatch_smem!(cp, mindis, maxdis, minv, maxv, colstart, colend, rowstart, rowend, phaL, phaR, threadsPerBlock, blocksPerGrid, stride, w, h, winSize, pha_dif)
#---------------------------------------------------
# cp: diparity map
# mindis, maxdis 最近最遠視差
#minv, maxv 仿射變換計算得到的R圖中有效橫向範圍
# colstart, colend, rowstart, rowend仿射變換計算得到的左圖中有效橫向、豎向範圍
#phaL, phaR 左右圖像
#w, h圖像大小
#winSize, pha_dif 3*3的框； 閾值：約等於20個像素的平均相位距離和
# Set up shared memory cache for this current block.
#---------------------------------------------------
wh = fld(winSize, 2)
cache = @cuDynamicSharedMem(Float32, threadsPerBlock)
left_stride = 64
minv = max(1,minv)#必須是有效值，且是julia下的下標計數方式
colstart= max(1,wh*stride)
# 數據讀入共享內存
j = blockIdx().x + rowstart # 共用的行序號
i = threadIdx().x + colstart# 左圖的列序號
while(j <= min(rowend,h - wh))
#數據拷貝到共享內存中去，並將由threadsPerBlock共享
ri = (threadIdx().x - 1) * stride + minv
tid = threadIdx().x
while(tid <= threadsPerBlock && ri <= maxv)
cache[tid] = phaR[j, ri]
tid+=threadsPerBlock
ri+=threadsPerBlock
end
# synchronise threads
sync_threads()
maxv = min(fld(maxv - minv,stride) * stride + minv,maxv)
# 計算最小匹配項
while(i <= min(colend,w - (wh*stride)))
min_v = 10000
XR = -1
VV = phaL[j, i]
if(VV > 0.001f0)
kStart = max(minv, i - maxdis) + 1
kEnd = min(maxv, i - mindis) - 1
for k = kStart:kEnd #遍歷一整行
RK = cache[cld(k - minv + 1,stride)]#從0開始計數
if RK <= 0.001f0
continue
end
dif = abs(VV - RK)
if dif < pha_dif
sum = 0.0f0
sn = 1
for ki in 0:(winSize - 1)
R_local = cache[cld(k - minv + 1 - wh + ki,stride)]
(R_local < 1e-5) && continue
#phaL[j + kj - wh, i + ki - wh*stride]
VR = VV - R_local
sum = sum + abs(VR)# * VR
sn += 1
end
v = sum / sn
if v < min_v
min_v = v
XR = k
end
end
end
#需要作插值
#http://discourse.julialang.org/t/base-function-in-cuda-kernels/21866
if XR > 0
XR_new = bisection(VV, phaR, Float32(j), Float32(XR - 3), Float32(XR + 3))
# 注意，這裏直接做了視差處理了
state = (i - XR_new) > 0
@inbounds cp[j, i] = state ? (i - XR_new) : 0.0f0
end
end
i+=threadsPerBlock
end
j+=blocksPerGrid
end
sync_threads()
return
end
4. 使用金字塔原理再次加速
1）這部分和上述代碼的區別是，沒有讓少於列數目的線程進行多次運算，(上述代碼中有個自增操作，是讓線程進行了多次運算，原因是可分配線程總數不夠)
2）整體是金字塔模式，即相鄰者相似的原理。
3）時間消耗在TX2上降低到了72ms，相比於前一種方法又有了58%的降幅。
#進行立體相位匹配
function phaseMatch_smem!(cp, mindis, maxdis,
minv, maxv, colstart, colend,
rowstart, rowend, phaL, phaR,
threadsPerBlock, blocksPerGrid, stride,
w, h, winSize, pha_dif)
#---------------------------------------------------
# cp: diparity map
# mindis, maxdis 最近最遠視差
#minv, maxv 仿射變換計算得到的R圖中有效橫向範圍
# colstart, colend, rowstart, rowend仿射變換計算得到的左圖中有效橫向、豎向範圍
#phaL, phaR 左右圖像
#w, h圖像大小
#winSize, pha_dif 3*3的框； 閾值：約等於20個像素的平均相位距離和
# Set up shared memory cache for this current block.
#---------------------------------------------------
wh = fld(winSize, 2)
cache = @cuDynamicSharedMem(Float32, threadsPerBlock)
minv = max(1,minv)#必須是有效值，且是julia下的下標計數方式
colstart = max(1, colstart, wh*stride)
rowstart = max(1, rowstart, wh)#需要算3*3的矩陣，目前沒有計算，所以不用嚴格滿足>wh
# i 最大、最小區間，所需要的迭代次數，或者可以看成所需處理的步長
left_stride = Int(cld(colend - colstart, threadsPerBlock))
j = (blockIdx().x - 1) + rowstart # 共用的行序號
i = (threadIdx().x - 1) * left_stride + colstart# 左圖的列序號
while(j <= min(rowend,h - wh))
#數據拷貝到共享內存中去，並將由threadsPerBlock共享
ri = (threadIdx().x - 1) * stride + minv
tid = threadIdx().x
while(tid <= threadsPerBlock && ri <= maxv)
cache[tid] = phaR[j, ri]
tid+=threadsPerBlock
ri+=threadsPerBlock
end
# synchronise threads
sync_threads()
maxv = min(fld(maxv - minv,stride) * stride + minv,maxv)
# 計算最小匹配項
#end_i = i + left_stride
#while(i <= min(colend, w - (wh*stride))) #end_i - 1,
if(i <= min(colend, w - (wh*stride))) #end_i - 1,
min_v = 10000
XR = -1
VV = phaL[j, i]
if(VV > 0.001f0)
kStart = max(minv, i - maxdis)
kEnd = min(maxv, i - mindis)
for k = kStart:kEnd #遍歷一整行
RK = cache[cld(k - minv + 1,stride)]#從0開始計數
if RK <= 0.001f0
continue
end
dif = abs(VV - RK)
if dif < pha_dif
sum = 0.0f0
sn = 1
for ki in 0:(winSize - 1)
R_local = cache[cld(k - minv + 1 - wh + ki,stride)]
(R_local < 1e-5) && continue
#phaL[j + kj - wh, i + ki - wh*stride]
VR = VV - R_local
sum = sum + abs(VR)# * VR
sn += 1
end
#v = sqrt(sum)/ sn
v = sum / sn
if v < min_v
min_v = v
XR = k
end
end
end
#需要作插值
#http://discourse.julialang.org/t/base-function-in-cuda-kernels/21866
if XR > 0
for offset in 0:left_stride - 1
i += offset
VV = phaL[j, i]
XR_new = bisection(VV, phaR, Float32(j), Float32(XR - 3), Float32(XR + 3))
# 注意，這裏直接做了視差處理了
state = (i - XR_new) > 0
@inbounds cp[j, i] = state ? (i - XR_new) : 0.0f0
end
end
end
end
j+=blocksPerGrid
end
sync_threads()
return
end
四、總結
julia作為一個交互性強的語言，在上述目的的達成上，它至少是做到了。對於上述算是一個標準的cuda加速過程，用cuda-c進行編寫的話，也是類似的過程，典型操作。但是c++編寫的可視化、調試、最終的編譯要麻煩的多得多得多。julia達到高效率的目的的同時，讓編寫的過程沒有那麼痛苦，堪稱完美的一次體驗。
- 應用健康度隱患刨析解決系列之數據庫時區設置
- 對於Vue3和Ts的心得和思考
- 一文詳解擴散模型：DDPM
- zookeeper的Leader選舉源碼解析
- 一文帶你搞懂如何優化慢SQL
- 京東金融Android瘦身探索與實踐
- 微前端框架single-spa子應用加載解析
- cookie時效無限延長方案
- 聊聊前端性能指標那些事兒
- Spring竟然可以創建“重複”名稱的bean？—一次項目中存在多個bean名稱重複問題的排查
- 京東金融Android瘦身探索與實踐
- Spring源碼核心剖析
- 深入淺出RPC服務 | 不同層的網絡協議
- 安全測試之探索windows遊戲掃雷
- 關於數據庫分庫分表的一點想法
- 對於Vue3和Ts的心得和思考
- Bitmap、RoaringBitmap原理分析
- 京東小程序CI工具實踐
- 測試用例設計指南
- 當你對 redis 説你中意的女孩是 Mia