作者：京东科技 李大冲
一、Julia是什么
julia是2010年开始面世的语言，作为一个10后，Julia必然有前辈们没有的特点。Julia被冀望塑造成原生的有C++的运行速度、python的易交互性以及胶水性。最重要的是，julia也是nVidia-cuda官网选中的接口语言，这点让cuda编程变得略微容易一些。
二、我的项目背景
双目平面视觉是比拟惯例获取深度信息的算法。这是一种模拟人的双眼的办法，依据对同一个点在不同的视角进行察看，并在不同的察看视角上的察看差别推算出深度，这在图1中更清晰一点。被动双目的毛病是须要纹理特色去判断左右相机里是否为同一个点。也就是，当没有问题、纹理不够显著或者纹理反复时，这种办法难以工作。
图 1 被动双目示意图。 观测点在一个相机和在另外一个相机上的成像是有显著可辨别性的，二者的察看差别（L1和L2）造成了视差，进而可造成深度信息
被动双目的毛病是某些场景下有效，被动双目便呈现了。如果咱们投射一些人工设定编码的图案，依据实现约定通过解码能够对视线内是所有点赋予惟一的标识符号，那被动双目的问题不就解决了吗。事实上，人们很快发现能够应用空间上编码和工夫上编码两种形式。空间上编码，就像人们察看河汉命名星座一样，去投射随机生成的星星点点的光斑，没有纹理的夜空照耀进了齐全可辨识的图案。如同星座往往是一个空间范畴很大的概念，空间编码的形式无奈做到十分精密。所以，工夫编码的形式应运而生，工夫编码是投射一组图案，让被照耀的区域，被照耀的明暗变动状况是齐全不一样的。像图2中的二进制编码（或行业内称为格雷码），再加上灰度变动更加丰盛的编码方式以及相机在竖方向上对齐之后是只须要思考横轴方向上的差别这两点，能够实现全视线的点齐全可辨别。
![星塵星雲星河| 何漢權- 灼見名家]()
图 2 空间编码和工夫编码。空间编码如同河汉中的星星点点，通过组合星点变为星座，人们能够对高深夜空不同地位依附河汉辨别开。工夫编码是投射多组图案，让每一个点是工夫上都有惟一的变动状况。就像右图，这16个方格都是齐全可辨别的。
三、效率问题
应用被动双目的构造去重建，一个很大的影响运行效率的问题是须要对左右图像求算视差（原理如图1所示）。具体求算视差的操作是要对左相机里的每一个点，去找右相机中同一行中与之间隔最近的点的横坐标，而后持续插值之后找到间隔更近的插值坐标，最初把这两个坐标的差别当作视差进行后续的计算。
应用工夫编码的被动双目，依照实现约定的编码进行解码之后，其成果如图3所示，这个环节就是在这对数据上进行解决。以此数据为例，无效区域大略是700850，齐全遍历所须要的工夫复杂度大略是O(700850*850)，间接估算就很耗时。
图 3 工夫编码解码后的图像示意图。左右相机拍摄的图案编程了解码值。接下来须要在这对数据上计算视差，而确定雷同编码值须要大量的遍历，对雷同编码进一步的优化还要对每个点进行小范畴的插值解决，工夫复杂度很大。
图 4 视差图成果
1. 应用for训练的形式
应用for循环，须要3层。经测试后，在TX2上，应用c++来实现大略须要7s多。
2. 应用numpy的矩阵操作来实现
1) 代码如下，应用numpy的播送的原理来实现，而后通过一系列的过滤异样值后，再通过插值进一步的解决失去最终想要的后果。
2) 这段代码没有进一步的去实现插值的过程，耗时状况为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, r3. 应用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() returnend #进行平面相位匹配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 #须要作插值 #https://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() returnend4. 应用金字塔原理再次减速
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 #须要作插值 #https://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() returnend四、总结
julia作为一个交互性强的语言，在上述目标的达成上，它至多是做到了。对于上述算是一个规范的cuda减速过程，用cuda-c进行编写的话，也是相似的过程，典型操作。然而c++编写的可视化、调试、最终的编译要麻烦的多得多得多。julia达到高效率的目标的同时，让编写的过程没有那么苦楚，堪称完满的一次体验。