CN107316324A - 基于cuda实现的实时立体匹配及优化的方法 - Google Patents

Abstract

本发明涉及计算机视觉的技术领域，更具体地，涉及基于CUDA实现的实时立体匹配及优化的方法。本发明是一种利用CUDA并行化处理来对左右输入图像进行密集的立体匹配，并且生成实时的视差图的方法。包括：对左右两图做census转换，生成一个字符序列，用汉明码距得到一个初始的cost，动态规划每个像素点8条路径，取一个最短的路径和来得到最终的cost，得到最初的一个密集的视差图；对左图利用k‑means算法进行超像素分割，得到一个个超像素平面块，利用超像素平面拟合来优化初始视差进行优化。本发明还涉及到多任务GPU并行加速，具体涉及在NVIDA的CUDA架构并行实现多个任务，属于GPGPU计算领域。通过GPU多线程处理优化极大的缩小了计算时间，得到实时的视差图。

Description

基于CUDA实现的实时立体匹配及优化的方法

技术领域

本发明涉及计算机视觉的技术领域，更具体地，涉及基于CUDA实现的实时立体匹配及优化的方法。

背景技术

立体视觉是计算机视觉中的一个重要课题，它可以使得人们获取到物体和场景中的深度信息，是后期3D重建和内容分析的基础。立体匹配是立体视觉的关键技术之一，通过标定好的双目相机拍摄的左右图片，经过矫正对齐之后，立体匹配是直接获取视差信息的关键性步骤。

立体匹配的直接目的是获取精确的视差信息，而匹配的精确度和时效性是衡量匹配方法的标准，如何平衡精准度和时效性是一个重要的课题。目前主流的立体匹配方法分为局部匹配，半全局匹配和全局匹配三种，精准度依次增加，同时消耗的时间也依次增加。

半全局匹配在精准度和时效性平衡性较好，因此我们采用半全局匹配，并着力提高精准度，降低运行时间。

在前景和后景的边缘部分，精准度难以把握，如何在前景膨胀和噪点中取舍得到好的效果是提高精确性的重要因素。

通过对图像进行超像素平面分割，是前后景边缘分割的重要手段，而分割方法也是计算机视觉领域的研究热点之一，在图像处理，对象识别中也有大量应用。通过快速而准确的超像素平面块分割，拟合，显著提高了视差的精确度。

GPGPU是GPU进行大规模计算的技术，而CUDA是NVIDA公司提供的GPGPU架构，GPU具有远高于CPU的浮点运算能力和内存带宽，有高度的并行性，不同于CPU程序只能一个一个的顺序执行，GPU的硬件设计支持多任务共享，程序可以通过多个线程同时运行。

因此将程序通过CUDA运行在GPU上，能极大缩短运行时间，提高时效性。

目前尚未发现在具备高精确度并且能输出实时的视差图的专利或者文献。

发明内容

本发明为克服上述现有技术所述的至少一种缺陷，提供基于CUDA实现的实时立体匹配及优化的方法，是兼顾精准和实时的立体匹配方法，加入了超像素分割用以拟合视差平面提高精确性，降低误匹配几率，利用了CUDA多线程并行化设计，将原本的算法程序计算时间由秒级别降低到毫米级别，能够达到实时的密集匹配。

本发明的技术方案：基于CUDA实现的实时立体匹配及优化的方法，其中，包括以下步骤：

S1.选定一个框大小，对左右图像除边缘部分的每个像素都做census转换为字符序列；

S2.在给定的视差范围内，选定左图像素点，遍历右图范围内各个字符序列，计算汉明码距，得到初始的cost；

S3.路径聚合，通过动态规划找到每一个点的路径聚合值最小，得到能量函数；

S4.用WTA选取能量函数最小所对应的视差值；

S5.左右一致性检查进行后期校验，得到初始的密集视差图；

S6.将左图每个像素由RGB颜色空间转换为CIElab颜色空间；

S7.通过k-means聚合算法划分出超像素平面块，迭代几次直至收敛，并融合超像素平面块；

S8.对每个超像素平面多次采样，计算平面参数；

S9.用得到的平面参数计算得到新的视差值；

S10.对新的视差做插值运算，减小黑色块，平滑视差图。

基于CUDA的编程算法都运用到以上步骤中。

与现有技术相比，有益效果是：在本发明中，多次使用CUDA流，共享内存块，线程同步，多线程共享并行化设计，并且优化了数据结构，都极大提高了运算效率。在对左右图像进行census转换，计算汉明码距得出cost，用动态规划法优化路径聚合的能量值，得到初始的密集视差图，基于对超像素平面的理解，用k-means聚类收敛得到各个超像素块，将每个超像素块看成是一个平面，我们计算出平面参数用于拟合新的视差值。整个一套流程因为有超像素分割和平面拟合的加入，使得视差图结果精确性得以提高，而选择使用CUDA并行化编程也使得实时性得到了保证，因此本发明在精确性和时效性方面做到了很好的平衡，这在目前来看并没有别的文献和专利做到过。

附图说明

图1表示整套算法的流程图。

图2表示路径聚合使用8条CUDA流并行时线程格的分配和计算方式。

图3表示k-means聚合得到的分割图。

图4表示最终生成的视差图渲染图。

具体实施方式

附图仅用于示例性说明，不能理解为对本专利的限制；为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。附图中描述位置关系仅用于示例性说明，不能理解为对本专利的限制。

CUDA使用线程格管理线程，核函数的每个副本都可以通过内置变量blockIdx来确定线程块的索引，gridDim得到线程块的数量,blockDim得到线程的数量，使用不同的gridDim和blockDim是获得不同并行化程度的关键。当需要多个线程块并且每个线程块又包含多个线程时，索引计算方法为：int tid＝threadIdx.x+blockIdx.x*blockDim.x；得到正确的索引tid才能定位到我们需要处理的线程，建立起映射关系，具体见图2.

GPU的内存结构中可读写的有：寄存器、共享内存和全局内存。CUDA上下文切换机制非常高效，几乎是零开销，访问寄存器的速度是非常快的，所以应尽量优先使用寄存器。CUDA共享内存缓冲区驻留在物理GPU上，所以访问速度也很快，共享内存的存在就可使运行线程块中的多个线程之间相互通信。全局内存因为GPU与CPU都可以对它进行写操作，速度较慢。因此在CUDA编程中，合理使用各种内存方式能够有效的提高效率和运行速度。

CUDA可以使用线程同步机制__syncthreads()；对线程块中的线程进行同步。该函数可以确保同一线程块中的所有线程在执行完该语句之前的所有语句之后才会执行其下一条语句，一般会和共享内存搭配出现。

我们也采用了CUDA流的方式，CUDA流表示一个GPU操作队列，并且该队列中的操作以添加到队列的先后顺序执行。使用CUDA流可以实现任务级的并行，内核函数之间可以并行处理，当GPU在执行核函数的同时，还可以在主机和设备之间交换数据。

步骤一中，我们选取的是n*m的块大小进行census转换。在CUDA程序中，选取了16*16的线程块，以方便后面用共享内存提高效率，图片的宽度是width，高度是height，选取的kernel大小分别是((width+15)/16，(height+15)/16)。共享内存的大小是(n+16)*(m+16)，先将数据从全局内存复制到共享内存，然后进行census转换，转换后同步线程，将字符串结果保存在gpu全局内存中。Census转换的思想是将选中的像素点，假设为(x，y)点的灰度值与它邻域块大小的各个像素点的灰度值进行比较，比(x,y)大的设为0，小的为1。

步骤二中，我们选取的视差范围是d，在cuda中选取32*8大小的线程块，使用(32+d)*8大小的共享内存，先将步骤一中得到每个像素点的字符序列从gpu全局内存复制到共享内存中，然后从右图中找d个视差范围内对应的字符序列，通过异或操作计算两对字符序列的汉明码距。汉明码距表示两个相同长度字对应位不同的数量，我们以d(x,y)表示两个字x,y之间的汉明距离。对两个字符串进行异或运算，并统计结果为1的个数，那么这个数就是汉明码距。遍历d个视差，选择汉明码距最小的视差存储到全局内存中，得到我们需要的初始的cost值。

步骤三中，路径聚合是为了优化能量函数从而查找到最小化的能量函数以确定视差值。

Lr(p,d)＝C(p,d)+min(Lr(p-r,d),

Lr(p-r,d-1)+P1,Lr(p-r,d+1)+P1,

在上述公式中，C(p,d)就是上个步骤我们得到的在p点视差为d时对用的cost值。而Lr(p,d)则表示在盘点在r路径上视差为d的路径cost和值，表述为每一个点的代价聚合值是“当前代价+min(路径相邻点的当前视差代价聚合值，路径相邻点的视差差值为1的代价聚合值+P1，路径相邻点的视差插值大于1的最小代价聚合值+P2)-路径相邻点的视差插值大于1的最小代价聚合值”，最后那一项是为了防止数字过大而设置。其中P1，P2分别是相邻像素点视差值相差为1，和大于1的惩罚系数，加入这两个正则化项是为了保证视差图平滑同时保持边缘。P1是对物体不平整表面的适应，P2是对灰度不连续(可能是多个物体)的适应。

想在2D的图像找到使得cost聚合值最小的视差图是个NP-hard的问题，那么单独在某方向r上(比如从左到右，从右到左等)选择匹配点的路径使得cost聚合值最小(即最终L最小)，该问题就变成多项式时间内可以解决的。我们采用具有全局性质的动态规划算法来优化计算出最小值。若要使p的代价想要最小，那么前提必须是邻域内的点q的代价最小，q想要代价最小，那么必须保证q的邻域点m的代价最小，如此传递下去分而治之得出最小值。

如果我们只沿着每一行求解，那么行间的约束完全考虑不到，q是p的邻域的点其实这个时候被弱化到了q是p的左侧点或者右侧点，这样的求优效果差并且会出现长尾效应。于是，我们选择采用8个方向路径进行优化。也就是r＝8，最后结果是8条路径的聚合值，也就是我们要求的能量函数。

在我们实现的程序中，使用了8个CUDA流，在kernel函数执行的同时，还可以传输数据，并且各个路径能同时进行运算，以便提高效率，分别使用了4个32*8大小的block，(height/8,1)大小的grid来计算水平路径和垂直路径，使用了4个32*8大小block，((height+width)/8,1)大小的grid来计算斜向路径，使用了32*8大小的共享内存来存储像素点的路径聚合cost值，在求取最小值时采用并行规约和syncthreads线程同步的方法，极大减小了线程束的分化，并且在循环中多次用pragma unroll指令使得循环充分展开并行操作，提高编译器效率，图2展示使用8个CUDA流，内核之间并行运行以及内存空间分配的关系图。

步骤四中，WTA(Winner-take-all)是一种竞争学习规则，用于无导师学习，用在这里是指只有最小的能量值会激活使用，即选取能量值最小的视差值。在本程序里依然选用32*8大小的线程块，128*8大小的共享内存块，先将上一步骤中得到的能量函数从全局内存复制到共享内存中，然后依次比较，将最小的结果进行并行规约以减少分支，再将视差值d保存在新的全局内存中。

步骤五中，在视差图中经常会采用到后续的处理措施：左右一致性检测。遮挡是只出现在一幅图像，而在另一幅图中看不到的那些点。遮挡点通常是一块连续的区域，LRCcheck的作用是实现遮挡检测得到左图对应的遮挡图像。对于左图中的一个点p，求得的视差值是d1，那么p在右图里的对应点应该是(p-d1)，(p-d1)的视差值记作d2。若|d1-d2|>threshold，p标记为遮挡点，在本程序中采用16*16的block和(width/16,height/16)大小的grid，对左右图像遍历，每个thread对应一个像素点，若|d1-d2|>1，则将视差设置为0保存在全局内存.

步骤六中，已经进入了超像素分割平面的阶段。在本发明中采用SLIC方法，即简单的线性迭代聚类。分簇的依据是像素之间的颜色相似性与邻近性，将彩色图像转化为CIELAB颜色空间和XY坐标下的5维特征向量，然后对5维特征向量构造距离度量标准，对图像像素进行局部聚类的过程。因此首先要将图像从RGB色彩空间转换到CIElab色彩空间。首先RGB到XYZ的转换：

[X,Y,Z]＝[M]*[R,G,B]，其中M为一3x3矩阵：

XYZ到Lab的转换：

L＝116*f(Y1)-16

a＝500*(f(X1)-f(Y1))

b＝200*(f(Y1)-f(Z1))

其中f是一个类似Gamma函数的校正函数：

当x>0.008856时，f(x)＝x^(1/3)；

当x<＝0.008856时，f(x)＝(7.787*x)+(16/116)

X1、Y1、Z1分别是线性归一化之后的XYZ值。

在CUDA程序中，使用16*16大小的线程块，使得每一个像素都对应一个thread线程，每个thread都进行上述函数变化，计算效率基本上达到理论峰值，之后再将得到的值保存在gpu的lab全局变量中。

步骤七中，k-means算法是一种简单高效并且应用广泛的簇聚类方法其主要是来计算数据聚集的算法，主要通过不断地取离种子点最近均值的算法。将图片初始化分割n块超像素，并将n个簇心点lab信息和xy坐标保存在gpu全局内存中，使用w*h大小的block，每个thread对应一个像素点，都查找2w*2h区域类的簇心点信息，计算lab色彩信息和两个点位置之间的欧氏距离，每个thread将最小的欧氏距离的簇心点更新复制到全局内存，迭代多次直到收敛；

D_lab＝(P_lab.x-C_lab.x)²+(P_lab.y-C_lab.y)²+(P_lab.z-C_lab.z)²

D_xy＝(P_xy.x-C_xy.x)²+(P_xy.y-C_xy.y)²

D＝D_lab+w*D_xy

以上公式D_lab表示在lab色彩空间中P点与2w*2h区域内C点即簇心点的欧式距离平方和，D_xy表示在图片中P点与簇心点C的坐标的欧氏距离平方和，w是一个权重参数，每个像素更新后的簇心点的信息即为平面标签。设置一个阈值，小于这个阈值的超像素平面会跟周围的平面计算平均Lab色彩距离，距离最小的进行平面融合。图3是进行了k-means并且融合后的超像素分割图。

步骤八中，在步骤七中已经对图像中每个像素都打了标签，表明每个像素都有对应的超像素块。已知三个点可以表示一个平面，所以在每个超像素块随机采样3个视差值不为0且不同的点.由采样到的三个点确定平面方程：

A＝(x1×z2-x2×z1)×(y2×z3-y3×z2)-(x2×z3-x3×z2)×(y1×z2-y2×z1)

B＝y1×z2-y2×z1；

P0＝((z2×d1-z1×d2)×(y2×z3-y3×z2)-(z3×d2-z2×d3)×(y1×z2-y2×z1))/A；

P1＝(z3×d2-z2×d3-P0×(x2×z3-x3×z2))/B；

P2＝(d1-P0×x_i-P1×y_i)/z_i；

以上公式表示的是计算平面参数的方程，若A<0，则平面参数P0＝0,P1＝0,P2＝-1.若B<0，则

P1＝(z3×d2-z2×d3-P0×(x2×z3-x3×z2))/B；

P2参数取决于Z_i的大小是否小于0，在初始化第一次求平面参数时，Z均为1，(x,y,d)分别表示像素的X坐标，Y坐标和视差值。

在得到平面参数后，即可求得平面方程：S＝P0×x+P1×y+P2；对剩下的平面点用方程进行拟合，找出一个平面中拟合点的个数，循环多次，选取拟合点最多的平面参数，保存到结果中。对所有的拟合点进行累加，得出新的参数，进行新一轮的平面参数计算，并将最终结果保存在数据结构中。

步骤九中，在步骤八中我们已经得到了图像中每个像素点的标签和平面参数，因为超像素平面块分割让我们可以将一个平面块认为视差值是变化很小的，因此我们用对应的平面参数方程来代替初始的视差值，从而达到优化的效果。新的视差方程就是：d＝P0×x+P1×y+P2；如果d在[0,255]区间内的话，就用这个新的d值代替原来初始的视差值。注意在超像素平面块的边缘线上，并不进行替换，依然采用初始视差值。

步骤十中，对视差图进行插值。为了使视差图更加平滑，我们对视差值为0的像素进行插值处理。首先找到视差值为0的区间，假设为[start,end]，然后比较start-1和end+1点的视差值，将较小的那个视差值***到[start,end]区间。经过插值运算，在视差为0的地方过渡更加平滑，在前后景分割的边缘因为前景膨胀而造成的黑边也能比较好的平滑，提高了视插图的准确性。

在本发明中，我们多次使用CUDA流，共享内存块，线程同步，多线程共享并行化设计，并且优化了数据结构，都极大提高了运算效率。在对左右图像进行census转换，计算汉明码距得出cost，用动态规划法优化路径聚合的能量值，得到初始的密集视差图，基于对超像素平面的理解，用k-means聚类收敛得到各个超像素块，将每个超像素块看成是一个平面，我们计算出平面参数用于拟合新的视差值。整个一套流程因为有超像素分割和平面拟合的加入，使得视差图结果精确性得以提高，而选择使用CUDA并行化编程也使得实时性得到了保证，因此本发明在精确性和时效性方面做到了很好的平衡，这在目前来看并没有别的文献和专利做到过。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (2)

1.基于CUDA实现的实时立体匹配及优化的方法，其特征在于，包括以下步骤：

S1. 选定一个框大小，对左右图像除边缘部分的每个像素都做census转换为字符序列；

S2. 在给定的视差范围内，选定左图像素点，遍历右图范围内各个字符序列，计算汉明码距，得到初始的cost；

S3. 路径聚合，通过动态规划找到每一个点的路径聚合值最小，得到能量函数；

S4. 用WTA选取能量函数最小所对应的视差值；

S5. 左右一致性检查进行后期校验，得到初始的密集视差图；

S6. 将左图每个像素由RGB颜色空间转换为CIElab颜色空间；

S7. 通过k-means聚合算法划分出超像素平面块，迭代几次直至收敛，并融合超像素平面块；

S8. 对每个超像素平面多次采样，计算平面参数；

S9. 用得到的平面参数计算得到新的视差值；

S10. 对新的视差做插值运算，减小黑色块，平滑视差图。

2.根据权利要求1所述的基于CUDA实现的实时立体匹配及优化的方法，其特征在于：所述的步骤中，基于CUDA的编程算法都运用到以上步骤中。