关于c++:用于双目重建中的GPU编程juliacuda

作者：京东科技 李大冲

一、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, 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 
                #须要作插值
                #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()
    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

                #须要作插值
                #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()
    return
end

四、总结

julia 作为一个交互性强的语言，在上述目标的达成上，它至多是做到了。对于上述算是一个规范的 cuda 减速过程，用 cuda- c 进行编写的话，也是相似的过程，典型操作。然而 c ++ 编写的可视化、调试、最终的编译要麻烦的多得多得多。julia 达到高效率的目标的同时，让编写的过程没有那么苦楚，堪称完满的一次体验。