方纯洁，赖小波，石 磊

浙江中医药大学医学技术学院 杭州 310053

基于Census变换的立体匹配算法优化及医学影像应用*

方纯洁，赖小波#，石 磊

浙江中医药大学医学技术学院 杭州 310053

#通信作者，男，1981年6月生，博士，副教授，研究方向：医学图像处理，E-mail：shopo@zcmu.edu.cn

立体匹配；Census非参数变换；快速实现；医学影像应用

目的：建立并优化一种基于Census非参数变换的立体匹配算法。方法：详细分析了基于Census非参数变换的立体匹配，运用移动窗口技术和高速缓存，采用SSE2指令进行数据流的并行处理从而优化算法结构，缩短取指时间，提高计算效率。用Middlebury标准图像对对该算法进行验证，并用于临床阴道镜子宫颈图像对的立体匹配。结果：快速实现了对Middlebury标准图像对以及子宫颈图像对的立体匹配。结论：该算法能够有效提高立体匹配的实时性，可应用于外科手术导航等计算机辅助外科诊疗系统。

双目立体视觉是计算机视觉的一个重要分支，由不同位置的两台摄像机或者一台摄像机经过移动和旋转对同一场景进行拍摄，再计算空间点在图像中对应点的视差，然后经过一系列投影逆变换，最后获得该空间点的三维坐标。与其他诸如透镜板三维成像、全息照相术等获取立体视觉的方法相比，双目立体视觉直接模拟人类双眼处理场景图像的方式，可靠简便，已经广泛应用于机器人导航和3D远程会议等领域。立体匹配是双目立体视觉中最核心、最重要的一个步骤，主要解决场景中同一空间点在参考图像上的投影点如何在匹配图像上匹配到其对应点的问题。目前研究人员已经提出了一些快速立体匹配算法[1]。其中一种方式是采用盒滤波法[2]，一般采用零均值归一化(ZNCC)或绝对差总和(SAD)等简单的互相关方法进行立体匹配[3-4]，匹配性能不会受匹配窗口尺寸的影响。文献[5]对基于互相关的立体匹配算法进行了改进，并在两种不同的硬件配置系统中进行了实验，能够实时得到移动机器人周围环境的视差图。Wang等[6]在动态规划(dynamic programming, DP)[7]算法架构中引入垂直方向的自适应匹配代价聚合步骤，可以在保证实时性的前提下获得高质量的视差图；但此算法存在的主要问题是不能有效融合水平与垂直方向上的连续性约束，导致视差图带有明显的条纹状瑕疵。Felzenszwalb等[8]采用分层匹配的策略有效降低了基于置信度传播(belief propagation, BP)[9]的立体匹配算法的复杂度；但是此算法完成一幅320×288像素的立体匹配用时需要1 s左右，难以达到实时性的要求。因此，如何在算法的执行时间与立体匹配质量之间进行权衡取舍是一项艰巨的任务。Census变换是一种非参数局部变换，由Zabih和Woodfill提出并最早应用于立体匹配[10]。作者提出了一种基于Census变换的立体匹配快速实现算法并进行了图像验证，报道如下。

1 材料与方法

1.1 材料 研究所使用的实际场景图像对来自Middlebury网站，该网站是Middlebury大学提供的一个权威的立体匹配算法测试平台，可以对算法的性能进行综合评估；标准图像对来自Middlebury图片库，用来检验算法，这些图片是在各种不同的环境下拍摄，外极线均已校正，并提供了标准视差图；子宫颈图像对来自中国科学院深圳先进技术研究院生物医学与健康工程研究所。全部计算在Intel Pentium Ⅳ 2.4 GHz、2.0 G内存的普通计算机上采用Visual Basic C++ 6.0编程实现。

1.2 算法结构 基于Census变换的立体匹配算法流程图如图1所示。图中第1个模块是为了在同一时间获取参数相同的图像对，通常情况下一个折射系统能够很容易满足此要求。第2个模块是为了校正由于镜头畸变而引起的图像失真，校正之后可以使用水平极线约束，使输入图像的每一行都与基线平行。在进行极线校正之前，左、右摄像机需进行离线标定。若使用基于Census变换的相似性测度函数，则第3个模块用于计算图像中每个像素的Census比特串。第4个模块用于获得视差空间图像，即匹配代价的计算。 第5个模块用于获得稳定的匹配结果，通常采用基于窗口的互相关法进行匹配代价的聚合，最佳匹配点通过在视差搜索范围内采用优胜全选(winner-take-all,WTA)策略获得。基于窗口的互相关匹配代价聚合假设窗口内所有像素点的视差连续，窗口大小可以根据实际情况调整。为了滤除视差图中由于误匹配以及遮挡引起的不可靠匹配点，引入第6个模块，进行左右一致性检测。第7个模块是视差图的后处理模块，根据需要可选，目的在于进一步提高视差图的质量，获得亚像素级别的视差图。

图1 Census变换立体匹配算法流程图

1.3 算法结构的改进 为了提高立体匹配的实时性，作者提出了一种基于Census非参数变换的立体匹配快速实现算法。假设在进行立体匹配之前，左、右图像均已经过预处理和极线校正。在匹配代价的聚合步骤中，采用基于局部矩形窗口的WTA策略对左右匹配窗口的Hamming距离分别进行运算。若直接进行Hamming距离的求和运算，由于在每一个步骤中，每个像素的灰度值都必须计算一次，因此具有较大的计算复杂度。移动窗口技术[11]可以消除重复计算，提高计算效率。如图2所示，通过移动窗口技术将匹配窗口内Hamming距离的求和分解为按列和按行的Hamming距离求和。具体步骤如下。步骤1：通过计算获得初次的行像素Hamming距离的和，并储存在计算机的一级缓存区中。步骤2：后一行像素Hamming距离的和通过已知的前一行像素Hamming距离的和加上后一个像素的Hamming距离同时减去前一行(列)第一个像素的Hamming距离得到。列像素Hamming距离的和用同样的方式获得。一旦行、列像素Hamming距离和的初始值算出，那么后一行、列像素Hamming距离和的计算复杂度是固定的，即两次加法运算和两次减法运算。因此，假设匹配窗口大小为n×n，那么计算复杂度从n降为次数固定的加法运算和减法运算，计算复杂度与窗口尺寸无关[12]。

图2 移动窗口技术

1.4 存储器组织的优化 针对立体匹配应用，通常有3种存储器组织方式可供选择。其中，x、y分别为像素点的横坐标和纵坐标，d为视差搜索范围。

在图3方式1中，像素点纵坐标y变化最快，其次为像素点横坐标x，最后才是视差搜索范围d。视差搜索范围d每改变一次，系统都需要重新读取左、右图像的数据。对于分辨率较高的图像，存取图像会消耗大量的时间。同样的，在方式2中，像素点纵坐标y变化最快，其次为视差搜索范围d，最后才是像素点横坐标x。在遍历不同视差值d时也存在与方式1同样的问题。在方式3中，视差搜索范围d变化最快，其次为像素点纵坐标y，最后才是像素点横坐标x。此时在遍历不同视差值d时，只需读取一次该像素的值，由于该像素值可以存放在缓存区，可大大减少访问系统内存的时间。

1.5 SSE2并行处理指令加速 Streaming SIMD Extensions 2(SSE2)指令集是Intel公司在SSE指令集的基础上发展起来的。相比于SSE，SSE2新增了144个指令，扩展了Multiple Media Extension(MMX)部分和SSE部分，这些指令提高了诸多应用程序的运行性能。随着整数指令从64位扩展到128 位，MMX技术部分引进的单指令多数据流(single instruction multiple data, SIMD)整数类型操作的有效执行率成倍提高。特别的，SSE2指令可让软件开发员极其灵活地实施算法，并在运行诸如MPEG-2、MP3、3D图形等之类的软件时能够增强性能。SSE2指令集由8个128位的通用寄存器构成，其容量是MMX寄存器的2倍。在指令处理速度保持不变的情况下，通过SSE2优化后的应用软件运行速度也将提高2倍。使用SSE2指令集可有效加速立体匹配的下列运算。①匹配代价的聚合：用来计算两匹配窗口内像素Hamming距离的和。经过SSE2指令集加速，一次可以同时处理8个16比特的像素。具体实现过程如下所示：

movdqa xmm3,[src]; //将src的值赋给xmm3

movdqa xmm1,[dst]; //将dst的值赋给xmm1

addps xmm1, xmm3; //xmm1的值加上xmm3

的值，和赋值给xmm1。

②提取匹配代价最小的视差值：比较每个待匹配点的匹配代价，把匹配代价最小的视差值作为当前像素点的实际视差值。其中，匹配代价最小值的计算可以采用下列指令：

pmin disp; //获取最小视差。

图3 3种存储器组织方式

1.6 实际应用 为了验证该算法的有效性，首先采用该方法对4组Middlebury图像库标准图像对Map、Teddy、Tsukuba和Cones分别进行匹配[1]，其中，图像大小为320×240像素，匹配窗口大小为11×11像素，变换窗口大小为7×7像素。其次，对来自临床阴道镜的子宫颈图像对采用该算法进行立体匹配。

2 结果

优化前后算法执行时间的对比情况见表1。从表1可以看出，该算法能够有效提高立体匹配的实时性。

表1 优化前后执行时间对比

子宫颈图像对视差图如图4所示。其中，视差搜索范围d分别为0～16和0～32，匹配窗口大小为11×11像素，变换窗口大小为7×7像素。图5为d=32时算法优化前后的帧频对比情况。从图4和图5可以看出，该算法即使对医学图像对进行立体匹配亦能获得较高质量的视差图，并具有实时性。

A:左图像;B:右图像;C:d=16;D:d=32。图4 子宫颈图像对视差图

图5 优化前后帧数对比

3 讨论

立体匹配大体上可以分为基于局部的立体匹配和基于全局的立体匹配两大类[13]。基于局部的立体匹配能够利用对应点自身以及局部邻域约束信息，计算效率高，在纹理丰富区域的匹配正确率较高，但是在遮挡区域或视差不连续区域易产生误匹配。与其相比，基于全局的立体匹配一般能获得精度较高的视差图，但相应的参数设置较难，算法复杂且费时。近年来，虽然广大的科研工作者提出了许多精度较高的立体匹配算法，但是，这些算法绝大多数是计算密集型的，若要从复杂的场景图像中获取可靠的三维信息，需要高端的硬件设备支持。此外，尽管立体匹配的精度是一个很重要的因素，但在实际应用中算法的实时性尤为重要。

Census变换以图像的信号强度为基础，围绕某一特定区域内中心像素的灰度值计算相似性测度。由于Census变换对图像的径向失真具有较强的鲁棒性，因此在立体匹配中应用较为广泛[14]。Census变换将变换窗口内各邻域像素灰度值相对于中心像素灰度值的差异转换为0或1的一维向量形式。若变换窗口大小为n×n，那么Census变换的计算结果为n-1位的比特串，在实际应用中这种属性能够有效消除强噪声或光照不均匀对立体匹配性能的影响。

对于一个实时的应用程序来说，存储器的组织方式非常重要。存储器组织方式可以分为一级缓存、二级缓存和系统内存；存储容量按顺序逐渐增大，存储速度却按顺序逐渐降低。在图像处理应用中，由于数据量十分庞大，仅仅依靠缓存难以实现存储大容量数据的需要，这时部分数据必须通过系统内存进行存储。根据文献[14]的描述，存储器的传输带宽限制了计算机的多媒体应用，因此只有有效利用缓存，从算法结构上减少访问系统内存的时间，才能够实现大容量数据的存储。方式1和方式2在当前视差搜索范围d搜索结束之前不会使用到额外载入的像素值，因此缓存命中率低；而方式3使用到了这些数据，从而提高了高速缓存的命中率。因此，在计算完Hamming距离的和之后，作者选用方式3在视差搜索范围d内进行Census变换立体匹配。

总之，为了克服大多数立体匹配算法计算复杂度偏高以及实时性不强的局限性，作者提出了一种基于Census变换的立体匹配快速实现算法。该算法运用移动窗口技术和存储器组织优化技术简化计算过程，提高计算效率，最后采用SSE2多媒体增强指令并行处理加速。实验结果表明，该优化算法能够有效改善立体匹配的实时性，即使对来自阴道镜的子宫颈图像对进行立体匹配也能够获得满意的效果，可应用于外科手术导航等对实时性有较高要求的计算机辅助外科诊疗系统。如何进一步将该算法应用于医学图像的三维重建将是课题组今后的研究方向。

[1]CHEN W,ZHANG MJ,XIONG ZH.Fast semi-global stereo matching via extracting disparity candidates from region boundaries[J].IET Computer Vision,2011,5(2):143

[2]PARK H,MITSUMINE H,FUJII M.Fast detection of robust features by reducing the number of box filtering in SURF[J].IEICE Trans Inf Syst,2011,E94D(3):725

[3]TIPPETTS BJ,LEE DJ,ARCHIBALD JK,et al.Dense disparity real-time stereo vision algorithm for resource-limited systems[J].IEEE Transacti Circuits Systems Video Technol,2011,21(10):1547

[4]CHIANG Hsiung,LIN Ting,HOU Lun.Development of a stereo vision measurement system for a 3D three-axial pneumatic parallel mechanism robot arm[J].Sensors(Basel),2011,11(2):2257

[5]OLIVER F, THIERRY V, ERIC T, et al. Real time correlation-based stereo: algorithm, implementations and applications[R].France:INRIA, 1993.

[6]WANG L,LIAO M,GONG M,et al.High-quality real-time stereo using adaptive cost aggregation and dynamic programming:Third International Symposium on 3D Data Processing,Visualization,and Transmission,Proceedings[C].Univ North Carolina,Chapel Hill,NC:IEEE,2006:798

[7]刘赫伟,汪增福.一种沿区域边界的动态规划立体匹配算法[J].模式识别与人工智能,2010,23(1):38

[8]FELZENSZWALB PF,HUTTENLOCHER DP.Efficient belief propagation for early vision[J].Int J Comput Vis,2006,70(1):41

[9]SUN J,ZHENG NN,SHUM HY.Stereo matching using belief propagation[J].IEEE Trans Pattern Anal Mach Intell,2003,25(7):787

[10]ZABIH R,WOODFILL J. Non-parametric local transforms for computing visual correspondence[C].Stockholm,Swedwn:Proceedings of the Third European Conference on Computer Vision,1994: 151

[11]MÜHLMANN K,MAIER D,HESSER J,et al.Calculating dense disparity maps from color stereo images, an efficient implementation[J].Int J Comput Vis,2002,47(1/3):79

[12]IBARRA-MANZANO MA,ALMANZA-OJEDA DL,DEVY M,et al.Stereo vision algorithm implementation in FPGA using census transform for effective resource optimization：Proceedings of The 2009 12th Euromicro Conference on Digital System Design,Architectures,Methods and Tools[C].Patras,Greece:IEEE,2009:799

[13]SCHARSTEIN D,SZELISKI R.A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[J].Int J Comput Vis,2002,47(1/3):7

[14]SEBOT J,DRACH-TEMAM N. Memory bandwidth: the true bottleneck of SIMD multimedia performance on a superscalar processor[Z].LNCS,2001,2150:439

(2016-08-04收稿 责任编辑徐春燕)

Research on a fast implementation algorithm based on Census transform stereo matching and its applications in medical imaging

FANGChunjie,LAIXiaobo,SHILei

CollegeofMedicalTechnology,ZhejiangChineseMedicalUniversity,Hangzhou310053

stereo matching;Census transform; fast implementation; applications in medical imaging

Aim: To propose a fast implementation algorithm based on Census transform stereo matching.Methods: Firstly, the Census transform stereo matching was investigated, and its implementation process was analyzed in detail. Secondly, in order to simplify the calculation process and improve the computational efficiency, the moving window technique and cache were employed to optimize the algorithm structure and shorten the fetch time of the instructions. Finally, the SSE2 instructions were adopted to realize the parallel processing for the data stream, and the proposed algorithm was tested with different image pairs.Results: Fast implementation for the stereo matching with the Middlebury standard image pairs and cervical image pair was achieved.Conclusion: The proposed algorithm can effectively improve the real time of the stereo matching based on Census transform, and it can be applied to the computer-aided surgical diagnosis and treatment system such as surgical navigation.

10.13705/j.issn.1671-6825.2017.02.010

*国家自然科学基金资助项目 61602419；浙江省自然科学基金资助项目 LY16F010008，LQ16F020003；浙江省教育厅科研项目

Y201431354；浙江中医药大学校级课题 2012ZY18

TP751.1