用于双目重建中的GPU编程：julia-cuda

一、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达到高效率的目的的同时，让编写的过程没有那么痛苦，堪称完美的一次体验。