图像处理基本算法-立体视觉

朴灿烈づ我的快乐病毒、 2022-04-04 08:28 263阅读 0赞

立体视觉是计算机视觉领域的一个重要课题，它的目的在于重构场景的三维几何信息。立体视觉的研究具有重要的应用价值，其应用包括移动机器人的自主导航系统，航空及遥感测量，工业自动化系统等。  
  
1. 引言  
立体视觉是计算机视觉领域的一个重要课题，它的目的在于重构场景的三维几何信息。立体视觉的研究具有重要的应用价值，其应用包括移动机器人的自主导航系统，航空及遥感测量，工业自动化系统等。  
一般而言，立体视觉的研究有如下三类方法:   
(1) 直接利用测距器（如激光测距仪）获得程距（range data）信息，建立三维描述的方法;   
(2) 仅利用一幅图象所提供的信息推断三维形状的方法;   
(3) 利用不同视点上的，也许是不同时间拍摄的，两幅或更多幅图象提供的信息重构三维结构的方法。   
第一类方法，也就是程距法 (range data method)，根据已知的深度图，用数值逼近的方法重建表面信息，根据模型建立场景中的物体描述，实现图象理解功能。这是一种主动方式的立体视觉方法，其深度图是由测距器(range finders)获得的，如结构光(structured light)、激光测距器(laser range finders) 等其他主动传感技术 (active sensing techniques)。这类方法适用于严格控制下的环境(tightly controlled domains)，如工业自动化的应用方面。  
第二类方法，依据光学成象的透视原理及统计假设，根据场景中灰度变化导出物体轮廓及表面，由影到形(shape from shading)，从而推断场景中的物体。线条图的理解就是这样的一个典型问题，曾经引起了普遍的重视而成为计算机视觉研究领域的一个焦点，由此产生了各种各样的线条标注法。这种方法的结果是定性的，不能确定位置等定量信息，该方法由于受到单一图象所能提供信息的局限性，存在难以克服的困难。  
第三类方法，利用多幅图象来恢复三维信息的方法，它是被动方式的。根据图象获取方式的区别又可以划分成普通立体视觉和通常所称的光流(optical flow)两大类。普通立体视觉研究的是由两摄像机同时拍摄下的两幅图象，而光流法中研究的是单个摄像机沿任一轨道运动时顺序拍下的两幅或更多幅图象。前者可以看作后者的一个特例，它们具有相同的几何构形，研究方法具有共同点。双目立体视觉是它的一个特例。  
立体视觉的研究由如下几部分组成:   
(1) 图象获取 (image acquisition),  
用作立体视觉研究的图象的获取方法是多种多样的，在时间、视点、方向上有很大的变动范围，直接受所应用领域的影响。立体视觉的研究主要集中在三个应用领域中，即自动测绘中的航空图片的解释，自主车的导引及避障，人类立体视觉的功能模拟。不同的应用领域涉及不同类的景物，就场景特征的区别来分，可以划分成两大类，一类是含有文明特征(cultural features)的景物，如建筑、道路等; 另一类是含有自然特征的景物和表面(natural objects and surfaces)， 如山、水、平原及树木等。不同类的景物的图象处理方法大不相同，各有其特殊性。  
总之，与图象获取相关的主要因素可归纳如下:   
(a) 场景领域 (scene domain)，   
(b) 计时 (timing)，  
(c) 时间（照明和阴影）(time of day (lighting and presence ofshadows))，   
(d) 成像形态（包括特殊的遮盖）(photometry (including special coverage)),  
(e) 分辨率 (resolution),  
(f) 视野 (field of view),  
(g) 摄像机的相对位置 (relative camera positioning).  
场景的复杂程度受如下因素的影响:   
(a) 遮掩 (occlusion),  
(b) 人工物体（直的边界,平的表面) (man-made objects (straight edge, flat surfaces)),  
(c) 均匀的纹理区域 (smoothly textured areas),  
(d) 含有重复结构的区域 (areas containing repetitive structure)。  
(2) 摄像机模型 (camera modeling),  
摄像机模型就是对立体摄像机组的重要的几何与物理特征的表示形式，它作为一个计算模型，根据对应点的视差信息，用于计算对应点所代表的空间点的位置。摄像机模型除了提供图象上对应点空间与实际场景空间之间的映射关系外，还可以用于约束寻找对应点时的搜索空间，从而降低匹配算法的复杂性，减小误匹配率。  
(3) 特征抽取 (feature acquisition),  
几乎是同一灰度的没有特征的区域是难以找到可靠匹配的，因而，绝大部分计算机视觉中的工作都包括某种形式的特征抽取过程，而且特征抽取的具体形式与匹配策略紧密相关。在立体视觉的研究中，特征抽取过程就是提取匹配基元的过程。  
(4) 图象匹配 (image matching),  
图象匹配是立体视觉系统的核心，是建立图象间的对应从而计算视差的过程，是极为重要的。  
(5) 深度计算 (distance(depth) determination),  
立体视觉的关键在于图象匹配，一旦精确的对应点建立起来，距离的计算相对而言只是一个简单的三角计算而已。然而，深度计算过程也遇到了显著的困难，尤其是当对应点具有某种程度的非精确性或不可靠性时。粗略地说，距离计算的误差与匹配的偏差成正比，而与摄像机组的基线长成反比。加大基线长可以减少误差，但是这又增大了视差范围和待匹配特征间的差别，从而使匹配问题复杂化了。为了解决这一问题出现了各种匹配策略，如由粗到精策略，松驰法等。   
在很多情况下，匹配精度通常是一个象素。但是，实际上区域相关法和特征匹配法都可以获得更好的精度。区域相关法要达到半个象素的精度需要对相关面进行内插。尽管有些特征抽取方法可以得到比一个象素精度更好的特征，但这直接依赖于所使用的算子类型，不存在普遍可用的方法。  
另一种提高精度的方法是采用一个象素精度的算法，但是利用多幅图象的匹配，通过多组匹配的统计平均结果获得较高精度的估计。每组匹配结果对于最后深度估计的贡献可以根据该匹配结果的可靠性或精度加权处理。  
总之，提高深度计算精度的途径有三条，各自涉及了一些附加的计算量:   
(a) 半象素精度估计 (subpixel estimation),  
(b) 加长基线长 (increased stereo baseline),  
(c) 几幅图的统计平均 (statistical averaging over several views)。  
(6) 内插 (interpolation).  
在立体视觉的应用领域中，一般都需要一个稠密的深度图。基于特征匹配的算法得到的仅是一个稀疏而且分布并不均匀的深度图。在这种意义下，基于区域相关匹配的算法更适合于获得稠密的深度图，但是该方法在那些几乎没有信息（灰度均匀）的区域上的匹配往往不可靠。因此，两类方法都离不开某种意义的内插过程。最为直接的将稀疏深度图内插成稠密的深度图的方法是将稀疏深度图看作为连续深度图的一个采样，用一般的内插方法（如样条逼近）来近似该连续深度图。当稀疏深度图足以反映深度的重要变化时，该方法可能是合适的。如起伏地貌的航空立体照片的处理中用这种方式的内插也许是比较合适的。但是这种方法在许多应用领域中，尤其是在有遮掩边界的图象的领域中，就不适用了。  
Grimson 指出可匹配特征的遗漏程度反映了待内插表面变化程度的相应限度，在这种基础上，他提出了一个内插过程\[2\]。换一角度来看，根据单幅图象的“由影到形”的技术，用已经匹配上的特征来建立轮廓条件和光滑的交接表面可以确保内插的有效性。这些方法结合起来，可以使内插过程达到合乎要求的目标。内插的另一种途径是在已有的几何模型与稀疏深度图之间建立映射关系，这是模型匹配过程。一般而言，要进行模型匹配，预先应将稀疏深度图进行聚类，形成若干子集，各自相应于一种特殊结构。然后找每一类的最佳对应模型，该模型为这种特殊结构(物体)提供参数和内插函数。如 Gennery用这种方法来发现立体对图片中的椭园结构，Moravec 用于为自主车探测地面。  
2. 双目立体视觉(Binocular Stereo Vision)  
2.1 双目立体视觉模型  
双目立体视觉理论建立在对人类视觉系统研究的基础上，通过双目立体图象的处理，获取场景的三维信息，其结果表现为深度图，再经过进一步处理就可得到三维空间中的景物，实现二维图象到三维空间的重构。Marr-Poggio-Grimson \[1\] 最早提出并实现了一种基于人类视觉系统的计算视觉模型及算法。双目立体视觉系统中，获取深度信息的方法比其它方式（如由影到形方法）较为直接，它是被动方式的，因而较主动方式（如程距法）适用面宽，这是它的突出特点。  
双目立体视觉系统中，深度信息的获得是分如下两步进行的:   
(1) 在双目立体图象间建立点点对应,  
(2) 根据对应点的视差计算出深度。  
第一部分，也就是对应点问题，是双目立体视觉的关键; 第二部分是摄像机模型问题。双目立体视觉模型中，双摄像机彼此参数一致，光轴平行且垂直于基线，构成一共极性 (epipolar) 结构，这样做是为了缩小对应的搜索空间，只有水平方向的视差，简化了对应过程，如下图所示。   
  
　   
  
如上图所示，设空间一点P(X,Y,Z)在两个平行放置的完全相同的摄象机中像点分别是(x1,y1).(x2,y2),则在知道基线长B和焦距f的情况下，可以计算出深度  
  
这是双目立体视觉的基本原理，即根据视差来恢复立体信息。  
2.2 匹配基元  
匹配基元是指匹配算法的最小匹配对象，它是由特征抽取算法产生的。在建立立体视觉系统时，必须根据环境的特点和应用的领域选择适当的匹配基元。匹配基元可以是：  
(1) 过零点 (zero-crossings),  
(2) 边界与线片段 (edge and line fragments),  
(3) 线性特征 (linear features),  
(4) 边缘轮廓 (object boundaries),  
(5) 兴趣算子抽取的特征点（如角点等）  
基元作为匹配算法处理的基本单位，是局部特征，应包含以下一些信息:   
(1) 维量（点、线、边界等） (dimensionality),  
(2) 尺度（空间频度，长短、大小、方向等）(size (spatial frequency)),  
(3) 亮度（对比度） (contrast),  
(4) 语义量 (semantic content),  
(5) 稠密度 (density of occurrence),  
(6) 简单可量度的分布特征 (easily measurable attributes),  
(7) 唯一性／突出性 (uniqueness/distinguishability)  
2.3 匹配算法  
匹配算法就是在两幅图象的匹配基元之间建立对应关系的过程，它是双目立体视觉系统的关键。实际上，任何计算机视觉系统中都包含一个作为其核心的匹配算法，因而对于匹配算法的研究是极为重要的。  
为了比较全面地考察匹配算法，这里不妨将双目立体视觉的匹配算法扩展到更一般的情况来分析：假设给定两幅同一环境的图象，这两幅图象可能由于摄取的时间、方位或方式的不同而有差别，如双目立体视觉系统所摄取的两幅图象、地图与遥感或航测图象等，如何找到彼此对应的部分? 对于这个问题，一般有两种考虑途径:   
(1) 灰度分布的相关性，  
(2) 特征分布的相似性。  
因而就有两类算法:   
(1) 基于灰度的算法 (intensity based)，  
(2) 基于特征的算法 (feature based)。  
如果按照控制策略分，有如下几种:   
(1) 粗到精多层次结构 (coarse-to-fine,hierarchical),  
(2) 引入约束条件的松驰法 (constraints, relaxation),  
(3) 多级表示的决策结构 (multilevel representation)。  
2.3.1 基于灰度的匹配算法  
基于灰度的算法是指图象处理中所称的区域相关方法 (area-correlation technique)，它是解决对应问题的一个最直观最简单的方法。在一幅图象中以一点为中心选定一区域（窗口），在另一幅图象中寻找与该区域相关系数最大的区域，把该找到的区域的中心认为是原来那区域中心的对应点。这里所说的图象包括经过某种特殊处理如Gauss滤波后的图象。  
这种算法计算量大，但可以得到整幅图象的视差图。该算法对噪音很敏感，考虑到计算量，窗口不宜开得过大，因而可能匹配的选择较大，误对应可能性大，不适于灰度分布均匀的图象，较适于灰度分布很复杂的图象，如自然景物等。采用该方法的关键在于排除或减轻噪音的影响。通常采用多层次相关对应及多幅图象的统计平均处理方式来实现。如 D. B. Gennery \[2\]采用九幅图象多级处理方式来实现对应求解。  
2.3.2 基于特征的匹配算法  
鉴于灰度区域相关方法的局限性，现在大部分研究集中在这方面。在许多环境（如有线条轮廓特征可寻的人工环境 (man-made structured world)）中，图象的特征是很有规律地分布的，反映了场景的核心，数量少，处理方便。基于特征的匹配算法特别适用于特殊的比较简单的环境如室内环境，具有速度快、精度高的特点，但对于自然环境，由于缺少显著的主导特征，该方法也遇到了很大困难。  
基于特征的双目立体视觉的对应算法，通过建立所选基元的对应关系，旨在获取一稀疏深度图，如果需要再经过内插等方法可以得到整幅深度图。这一类算法因各自采用的匹配基元不同而相异。概括而言，该类匹配算法都是建立在匹配基元之间的相似性度量基础上的。这种相似性度量被称为亲合性 (affinity)\[2\]， 它是以匹配基元的各项参数信息为依据的局部特征相似程度的度量。这种度量方法与摄像机模型相结合，可以大大减小匹配时的搜索空间。  
由于仅利用亲合性建立匹配是模糊的，可能匹配的空间仍旧很大（多对一的），因此有必要引入其它约束条件及控制策略来限制搜索空间，减小模糊程度。匹配算法中常引入的两种约束条件及控制策略是:   
(1) 共极性 (epipolar) （双目立体视觉模型特点）,  
(2) 连续性 (continuity),  
(3) 分层次的匹配策略（即由粗到精策略）(hierarchical (e.g.,coarse-fine) matching strategy)。  
这种引入约束的方法实际上是将有关环境模型的知识融于算法之中。  
这种算法的具体实现，可以采用概率度量、松驰法迭代或者聚类等模式识别算法来实现。作为最后结果的1-1 对应，可以利用启发式搜索方法从已经大大减小了的搜索空间中获得。这部分可望能利用现代 AI 研究的许多手段如专家系统等研究方法，作为承上启下，建立更高层次描述的先导。  
可以从以下几个角度来比较各种匹配算法，  
(1) 精度 (accuracy),  
(2) 可靠性（排除总体分类误差的程度）(reliability)，  
(3) 通用性（适于不同场景的能力）(available of performance models),  
(4) 预见性 (predictability),  
(5) 复杂性（设备及计算量的代价）(complexity (cost implementation,   
computational requirements))。  
立体视觉的匹配算法有:   
(1) Marr-Poggio-Grimson算法，以过零点为基元，采用由粗到精的控制策略，用精度较低层次的匹配来限定精度较高层次匹配的搜索空间，最后利用连续性约束通过迭代方式实现匹配过程。处理对象是自然景物的双目立体图象。  
(2) R. Nevatia-G.Medioni算法，以线片段 (segments) 为基元，以最小差别视差 (minimum differential disparity) 为基准，建立匹配过程。该基准实际上是连续性约束的一种表现形式，在对应线片段各自邻域内存在的对应线片段的视差与其视差相近。处理对象是人工环境的双目立体图象。  
(3) R. Y. Wong算法，旨在建立两类图象的对应关系，如航空照片、遥感图象与灰度图象之间的对应关系。以边界特征(edge feature)为依据采用顺序的 (sequential)、多层次结构 (hierarchical structure)的搜索策略实现匹配过程。  
(4) K. Price-R. Reddy算法，依据场景的线条特征模型，将自顶向下（人工智能）(top-down (artificial intelligence))与自底向上（模式识别）(bottom-up (pattern recognition)) 两种控制策略有效地结合起来，采用广义的相关方法进行匹配，旨在建立形态差别较大的两幅图象（一幅是参照图或参考模型，另一幅是待对应的图象）的对应关系。如机场模型与机场的航空照片之间的对应关系。  
(5) C. S. Clark-A. L. Luck-C. A. McNary算法，抽取线条轮廓特征建立模型，在模型间建立对应。适于存在较大差别的图象的匹配。  
(6) K. E. Price算法，用于在图象间建立区域对应。该算法利用区域间的相互关系，以松驰法为基本思想实现了多层次表示结构下的匹配过程。突出特点是匹配算法考虑了图象本身区域间的相互关系（如包含、子部分等）的匹配，具有类似于某种语义网络式的启发性。  
(7) R. Horaud-T. Skorads算法，以线条特征为匹配基元，每个线条特征不仅含有其本身的端点坐标及方向矢量信息，而且含有它同那些与其相邻的线条特征之间存在的相对位置及结构关系的信息。这些特征将每幅图象表示成为一个关系图，根据该关系图对于每个线条特征确定它在另一幅图象中的可能对应集合，以每组对应为一结点构造对应图，依据关系图的相容性通过利益函数（benefit function）确定最佳对应。它处理的对象是室内环境的双目立体图象。  
(8) W. Hoff-N. Ahuja算法，以过零点为最小特征，将特征匹配、轮廓检测以及表面内插这三个过程结合在一起，采用基于多层表示的由粗到精的控制策略，根据对于表面的光滑性约束重构三维表面。这是一种与传统方法大不相同的算法，适合于有纹理特征的环境如工作台上的物品，不适合于稀疏特征环境如室内环境。另外 S. I. Olsen提出的算法与此相似，它将表面的重构过程（reconstruction process）结合在对应匹配过程中，基于多重属性用松弛法进行匹配，逐步提高重构的视差表面与实际的视差数据的一致性。  
2.4 双目立体视觉系统  
双目立体视觉经过几十年的研究已经取得了显著了成果，出现了各种专门的硬件设计和视频速率（实时）的立体视觉系统，在理论和技术方面都比较成熟了。但是，从普遍的意义来讲，由于很难彻底地解决对应点问题，具体的立体视觉系统一般都是有针对性的、不是普遍适用的，还无法与人类的双目视觉系统相媲美。  
  
下图是SRI的集成在电路板上的双目立体视觉系统。CMU设计了Stereo Machine, 可以实时地获取深度信息。  
  
立体摄象机校准 Stereo Camera Calibration   
三维视觉   
Milan Sonka, 3D Vision  
集成在电路板上的立体摄象机对SRI Stereo Engine, Stereo head onboard   
立体几何模型 SRI Stereo Geometry  
双目立体视觉Introduction to Stereo Imaging -- Theory  
3. 结构光方法(Structured Light)  
将平面光束照射在物体上可以形成光带，光带的偏转数据反映了物体表面的三维形状信息，用这种方法可以精确地获取物体的三维信息。借助于一组平行的平面光，或将物体置于专门的旋转工作台上通过一束平面光，都可以利用偏转数据直接地计算出深度信息，称这种方法为结构光方法。结构光方法适合于限制条件下，局部范围内需要精确测量的情况，用于不规则表面的三维建模。  
结构光方法在工业上有重要的应用，例如从传送带上检测工件，工件的逆工程（Reverse engineering）；在图形建模方面也有重要的应用，如人体建模，包括头部等躯体模型，雕塑造型的数字化。实际上它是三维扫描仪的基本原理。  
如下图所示的装置，就是结构光方法的典型事例。  
  
详细可见：Our Active Stereo Vision System  
4. 激光雷达与程距数据(Range Data)处理  
激光雷达（Laser range finder）与结构光方法不同，它直接利用激光光速扫描物体，通过测量光束从发出到反射回来的时间差来计算深度信息。它提供的数据是深度图，称为程距数据（Range data）。激光雷达可以用于比较大范围的测量，如移动机器人可以用激光雷达信息来建立环境内模型，以实现自主导航、躲避障碍等功能。  
程距数据实际上就是深度图象，结构光方法和激光雷达得到的数据最后都是深度信息。程距数据处理主要是表面拟合，恢复物体的表面结构。  
5. 视觉临场感系统  
临场感（Telepresence）技术是新一代遥操作（Teleoperation）系统的重要组成部分。顾名思义，它的目的就是使人从远地遥控操作时具有在现场处实地操作式的身临其境的感觉。在理想情况下，这些感觉应该包括人的各种感官所能感受到的感觉，如视觉、听觉、触觉、味觉、体位觉、力感等。  
临场感系统因其面对的任务不同，所需的现场信息有所区别，其中，视觉通常是最重要的信息之一，其次才是听觉、触觉等。目前，临场感技术主要涉及视觉和听觉。  
临场感遥操作系统的主要优点是：将人与机器人有机地结合起来，能够恰到好处地发挥出各自的特长。机器代替人去危险或人不可能到达的区域去工作，而人的判断能力和决策水平又明显地提高了系统的整体智能水平。  
如下图所示，室外车辆上的立体摄象机将视频信号传回基地端，操作员通过立体眼睛观察环行屏幕，仿佛他亲自在车上一样能够具有身临其境的感觉。  
  
（参见：艾海舟、张朋飞、何克忠、江潍、张军宇，室外移动机器人的视觉临场感系统，机器人，22（1）：28-32，2000。）  
有关立体视觉的前沿工作请参见微软研究院张正友博士的网页，他是这方面的著名学者：[http://research.microsoft.com/~zhang/][http_research.microsoft.com_zhang]  
参考文献  
1.马松德、张正友，计算机视觉计算理论与算法基础，科学出版社，1998。  
2.艾海舟,关于双目立体视觉的研究,硕士论文,121页,1988.4.  
3.艾海舟, 关于移动机器人自主式系统的研究, 博士论文, 153页, 1991.3.

最近一直学习立体视觉，写了很多的代码，但是还没整理具体的算法。使用左右两张图片，计算深度图

左图如下：

![1358164779_2048.jpg][]

右图如下：

![1358164806_3812.jpg][]

//SAD算法

效果：

![1358164851_3510.jpg][]

//SSD：

![1358164921_6648.jpg][]

//ZSSD

#include<iostream>
    #include<cv.h>
    #include<highgui.h>
    
    using namespace std;
    int GetHammingWeight(unsigned int value);
    int main(){
        /*Half of the window size for the census transform*/
        int hWin = 11;
        int compareLength = (2*hWin+1)*(2*hWin+1);
     
        cout<<"hWin: "<<hWin<<";  "<<"compare length:  "<<compareLength<<endl;  
        cout<<"ZSSD test"<<endl;
        // char stopKey;
      /*  IplImage * leftImage = cvLoadImage("l2.jpg",0);
        IplImage * rightImage = cvLoadImage("r2.jpg",0);*/
    
        IplImage * leftImage = cvLoadImage("left.bmp",0);
        IplImage * rightImage = cvLoadImage("right.bmp",0);
    
        int imageWidth = leftImage->width;
        int imageHeight =leftImage->height;
    
        IplImage * SADImage = cvCreateImage(cvGetSize(leftImage),leftImage->depth,1);
        IplImage * MatchLevelImage = cvCreateImage(cvGetSize(leftImage),leftImage->depth,1);
    
        int minDBounds = 0;
        int maxDBounds = 31;
       
        cvNamedWindow("Left",1);
        cvNamedWindow("Right",1);
        cvNamedWindow("Census",1);
        cvNamedWindow("MatchLevel",1);
    
        cvShowImage("Left",leftImage);
        cvShowImage("Right",rightImage);
    
    
    
        /*Census Transform */
        int i,j ,m,n,k;
        unsigned char centerPixel = 0; 
        unsigned char neighborPixel = 0;
        int bitCount = 0;
        unsigned int bigger = 0;
    
    
      
        int sumLeft = 0;
        int sumRight = 0;
        int sum =0;
        
        int zSumLeft  = 0;
        int zSumRight = 0;
    
        unsigned int *matchLevel = new unsigned int[maxDBounds - minDBounds  + 1];
        int tempMin = 0;
        int tempIndex = 0;
       
        unsigned char* dst;
        unsigned char* leftSrc  = NULL;
        unsigned char* rightSrc = NULL;
    
        unsigned char leftPixel = 0;
        unsigned char rightPixel =0;
        unsigned char subPixel = 0;
        unsigned char meanLeftPixel  = 0;
        unsigned char meanRightPixel = 0;
     
    
        for(i = 0 ; i < leftImage->height;i++){
            for(j = 0; j< leftImage->width;j++){
    
                /*均值计算 */
                for (k = minDBounds;k <= maxDBounds;k++)
                {
                    sumLeft  = 0;
                    sumRight = 0;
                    for (m = i-hWin; m <= i + hWin;m++)
                    {
                        for (n = j - hWin; n <= j + hWin;n++)
                        {
                              if (m < 0 || m >= imageHeight || n <0 || n >= imageWidth )
                              {
                                  sumLeft += 0;
                              }else {
                                  leftSrc = (unsigned char*)leftImage->imageData 
                                      + m*leftImage->widthStep + n + k; 
                                  leftPixel = *leftSrc;
                                  sumLeft += leftPixel;
                              }
                              
                              
                              if (m < 0 || m >= imageHeight || n + k <0 || n +k >= imageWidth)
                              {
                                   sumRight += 0;
                              }else
                              { 
                                  rightSrc = (unsigned char*)rightImage->imageData 
                                               + m*rightImage->widthStep + n;
                                  rightPixel = *rightSrc;
                                  sumRight += rightPixel;   
                              }
    
                        }
                    }
    
                    meanLeftPixel  = sumLeft/compareLength;
                    meanRightPixel = sumRight/compareLength;
                     /*ZSSD*/
                   sum = 0;
                     for (m = i-hWin; m <= i + hWin;m++)
                    {
                        for (n = j - hWin; n <= j + hWin;n++)
                        {
                              if (m < 0 || m >= imageHeight || n <0 || n >= imageWidth )
                              {
                                  //zSumLeft += 0;
                                  leftPixel = 0;
                              }else {
                                  leftSrc = (unsigned char*)leftImage->imageData 
                                      + m*leftImage->widthStep + n + k; 
                                  leftPixel = *leftSrc;
                                  //zSumLeft += (leftPixel - meanLeftPixel)*(leftPixel -meanLeftPixel);
                              }
                              
                              
                              if (m < 0 || m >= imageHeight || n + k <0 || n +k >= imageWidth)
                              {
                                   //zSumRight += 0;
                                  rightPixel = 0;
                              }else
                              { 
                                  rightSrc = (unsigned char*)rightImage->imageData 
                                               + m*rightImage->widthStep + n;
                                  rightPixel = *rightSrc;
                                  // zSumRight += (rightPixel - meanRightPixel)*(rightPixel - meanRightPixel);   
                              }
    
                              sum += ((rightPixel - meanRightPixel)-(leftPixel -meanLeftPixel))
                                     *((rightPixel - meanRightPixel)-(leftPixel -meanLeftPixel));
                        }
                    }
    
    
                    matchLevel[k] = sum;
                    //cout<<sum<<endl;
                }
    
                /*寻找最佳匹配点*/
               // matchLevel[0] = 1000000;
    
                tempMin = 0;
                tempIndex = 0;
                for ( m = 1;m < maxDBounds - minDBounds + 1;m++)
                {
                    //cout<<matchLevel[m]<<endl;
                    if (matchLevel[m] < matchLevel[tempIndex])
                    {
                        tempMin = matchLevel[m];
                        tempIndex = m;
                    }
                }
                dst = (unsigned char *)SADImage->imageData + i*SADImage->widthStep + j;
                //cout<<"index: "<<tempIndex<<"  ";
    
                *dst = tempIndex*8;
                dst = (unsigned char *)MatchLevelImage->imageData + i*MatchLevelImage->widthStep + j;
                *dst = tempMin;
                //cout<<"min:  "<<tempMin<<"  ";
                //cout<< tempIndex<<"  " <<tempMin<<endl;
            }
           // cvWaitKey(0);
        }
       
        cvShowImage("Census",SADImage);
        cvShowImage("MatchLevel",MatchLevelImage);
        cvSaveImage("depth.jpg",SADImage);
        cvSaveImage("matchLevel.jpg",MatchLevelImage);
    
    
        cout<<endl<<"Over"<<endl;
        cvWaitKey(0);
        cvDestroyAllWindows();
        cvReleaseImage(&leftImage);
        cvReleaseImage(&rightImage);
        return 0;
    }

![1358164982_3041.jpg][]

//census

#include<iostream>
    #include<cv.h>
    #include<highgui.h>
    
    using namespace std;
    int GetHammingWeight(unsigned int value);
    int main(){
        /*Half of the window size for the census transform*/
        int hWin = 11;
        int bitlength = 0;
        if ((2*hWin+1)*(2*hWin+1)%32 == 0)
        {
            bitlength = (2*hWin+1)*(2*hWin+1)/32;
        }else {
            bitlength = (2*hWin+1)*(2*hWin+1)/32  + 1; 
        }  
        cout<<"hWin: "<<hWin<<";  "<<"bit length:  "<<bitlength<<endl;  
        cout<<"Census test"<<endl;
       // char stopKey;
        IplImage * leftImage = cvLoadImage("left.bmp",0);
        IplImage * rightImage = cvLoadImage("right.bmp",0);
        
        int imageWidth = leftImage->width;
        int imageHeight =leftImage->height;
    
        IplImage * CensusImage = cvCreateImage(cvGetSize(leftImage),leftImage->depth,1);
        IplImage * MatchLevelImage = cvCreateImage(cvGetSize(leftImage),leftImage->depth,1);
    
         int minDBounds = 0;
         int maxDBounds = 31;
        // int leftCensus[imageHeight][imageWidth][bitlength] = {0};  
         unsigned  int *leftCensus = new unsigned int[imageHeight*imageWidth*bitlength];
         unsigned  int *rightCensus = new unsigned int[imageHeight*imageWidth*bitlength]; 
         for (int i = 0;i < imageHeight*imageWidth*bitlength;i++)
         {
             leftCensus[i] = 0;
             rightCensus[i] = 0;
         }
         int pointCnt = 0;
       
      
    
    
        cvNamedWindow("Left",1);
        cvNamedWindow("Right",1);
        cvNamedWindow("Census",1);
        cvNamedWindow("MatchLevel",1);
    
        cvShowImage("Left",leftImage);
        cvShowImage("Right",rightImage);
      
     
    
        /*Census Transform */
       int i,j ,m,n,k,l;
       unsigned char centerPixel = 0; 
       unsigned char neighborPixel = 0;
       int bitCount = 0;
       unsigned int bigger = 0;
        for(i = 0 ; i < leftImage->height;i++){
            for(j = 0; j< leftImage->width;j++){
                 centerPixel = *((unsigned char *)leftImage->imageData + i*leftImage->widthStep + j); 
                 bitCount = 0;
                 
                 for (m = i - hWin; m <= i + hWin;m++)
                 {
                     for (n = j - hWin; n<= j+hWin;n++)
                     {
                         bitCount++;
                         if (m < 0 || m >= leftImage->height || n < 0 || n >= leftImage->width)
                         {
                             neighborPixel = 0;
                         }else{
                              neighborPixel = *((unsigned char *)leftImage->imageData + m*leftImage->widthStep + n); 
                         }
                         bigger = (neighborPixel > centerPixel)?1:0;
                         leftCensus[(i*imageWidth + j)*bitlength + bitCount/32] |= (bigger<<(bitCount%32));
                     }
                 }
    
            }
        }
    
        for(i = 0 ; i < rightImage->height;i++){
            for(j = 0; j< rightImage->width;j++){
                centerPixel = *((unsigned char *)rightImage->imageData + i*rightImage->widthStep + j); 
                bitCount = 0;
    
                for (m = i - hWin; m <= i + hWin;m++)
                {
                    for (n = j - hWin; n<= j+hWin;n++)
                    {
                        bitCount++;
                        if (m < 0 || m >= rightImage->height || n < 0 || n >= rightImage->width)
                        {
                            neighborPixel = 0;
                        }else{
                            neighborPixel = *((unsigned char *)rightImage->imageData + m*rightImage->widthStep + n); 
                        }
                        bigger = (neighborPixel > centerPixel)?1:0;
                        rightCensus[(i*imageWidth + j)*bitlength + bitCount/32] |= (bigger<<(bitCount%32));
                    }
                }
    
            }
        }
        int sum = 0;
        unsigned int *matchLevel = new unsigned int[maxDBounds - minDBounds  + 1];
        int tempMin = 0;
        int tempIndex = 0;
        unsigned  char *dst;
        unsigned char pixle = 0;
        for(i = 0 ; i < rightImage->height;i++){
            for(j = 0; j< rightImage->width;j++){
              
    
                for (k = minDBounds;k <= maxDBounds;k++)
                {
                    sum = 0;
                    for (l = 0;l< bitlength;l++)
                    {   
                        if (((i*imageWidth+j+k)*bitlength + l) < imageHeight*imageWidth*bitlength)
                        {
                            sum += GetHammingWeight(rightCensus[(i*imageWidth+j)*bitlength + l] 
                            ^ leftCensus[(i*imageWidth+j+k)*bitlength + l]);
                        }else {
                            //sum += 0;
                           // cout<<".";
                        }
                       
                    }
                    matchLevel[k] = sum;
                }
    
                /*寻找最佳匹配点*/
                tempMin = 0;
                tempIndex = 0;
                for ( m = 1;m < maxDBounds - minDBounds + 1;m++)
                {
                     if (matchLevel[m] < matchLevel[tempIndex])
                     {
                         tempMin = matchLevel[m];
                         tempIndex = m;
                     }
                }
               
                if (tempMin > (2*hWin+1)*(2*hWin+1)*0.2)
                {
                    tempMin = 0;
                    pointCnt++;
                }else{
                    tempMin = 255;
                }
                dst = (unsigned char *)CensusImage->imageData + i*CensusImage->widthStep + j;
                *dst = tempIndex*8;
    
                dst = (unsigned char *)MatchLevelImage->imageData + i*MatchLevelImage->widthStep + j;
                *dst = tempMin;
    
                //cout<< tempIndex<<"  " <<tempMin<<endl;;
    
            }
        }
        cout<<"pointCnt:  "<<pointCnt<<endl;
        cvShowImage("Census",CensusImage);
        cvShowImage("MatchLevel",MatchLevelImage);
        cvSaveImage("depth.jpg",CensusImage);
        cvSaveImage("matchLevel.jpg",MatchLevelImage);
    
        cvWaitKey(0);
        cvDestroyAllWindows();
        cvReleaseImage(&leftImage);
        cvReleaseImage(&rightImage);
        return 0;
    }
    
    
    int GetHammingWeight(unsigned int value)
    {
        if(value == 0) return 0;
    
        int a = value;
        int b = value -1;
        int c = 0;
    
        int count = 1;
        while(c = a & b)
        {
            count++;
            a = c;
            b = c-1;
        }
        return count;
    }

![1358165030_2256.jpg][]

//NCC

#include<iostream>
    #include<cv.h>
    #include<highgui.h>
    #include <cmath>
    using namespace std;
    
    int main(){
        /*Half of the window size for the census transform*/
        int hWin = 11;
        int compareLength = (2*hWin+1)*(2*hWin+1);
     
        cout<<"hWin: "<<hWin<<";  "<<"compare length:  "<<compareLength<<endl;  
        cout<<"NCC test"<<endl;
      /*  IplImage * leftImage = cvLoadImage("l2.jpg",0);
        IplImage * rightImage = cvLoadImage("r2.jpg",0);*/
    
        IplImage * leftImage = cvLoadImage("left.bmp",0);                           
        IplImage * rightImage = cvLoadImage("right.bmp",0);
    
        int imageWidth = leftImage->width;
        int imageHeight =leftImage->height;
    
        IplImage * NCCImage = cvCreateImage(cvGetSize(leftImage),leftImage->depth,1);
        IplImage * MatchLevelImage = cvCreateImage(cvGetSize(leftImage),leftImage->depth,1);
    
        int minDBounds = 0;
        int maxDBounds = 31;
       
        cvNamedWindow("Left",1);
        cvNamedWindow("Right",1);
        cvNamedWindow("Census",1);
        cvNamedWindow("MatchLevel",1);
    
        cvShowImage("Left",leftImage);
        cvShowImage("Right",rightImage);
    
    
    
        /*Census Transform */
        int i,j ,m,n,k;
        unsigned char centerPixel = 0; 
        unsigned char neighborPixel = 0;
        int bitCount = 0;
        unsigned int bigger = 0;
    
        unsigned int sum =0;
        unsigned int leftSquareSum = 0;
        unsigned int rightSquareSum = 0; 
       
        double *matchLevel = new double[maxDBounds - minDBounds  + 1];
        double tempMax = 0;
        int tempIndex = 0;
       
        unsigned char* dst;
        unsigned char* leftSrc  = NULL;
        unsigned char* rightSrc = NULL;
    
        unsigned char leftPixel = 0;
        unsigned char rightPixel =0;
        unsigned char subPixel = 0;
        unsigned char meanLeftPixel  = 0;
        unsigned char meanRightPixel = 0;
     
    
        for(i = 0 ; i < leftImage->height;i++){
            for(j = 0; j< leftImage->width;j++){
    
                /*均值计算 */
                for (k = minDBounds;k <= maxDBounds;k++)
                {
                    sum = 0;
                    leftSquareSum  = 0;
                    rightSquareSum = 0;
    
                     for (m = i-hWin; m <= i + hWin;m++)
                    {
                        for (n = j - hWin; n <= j + hWin;n++)
                        {
                              if (m < 0 || m >= imageHeight || n <0 || n >= imageWidth )
                              {
                                 
                                  leftPixel = 0;
                              }else {
                                  leftSrc = (unsigned char*)leftImage->imageData 
                                      + m*leftImage->widthStep + n + k; 
                                  leftPixel = *leftSrc;
                                  
                              }
                              
                              
                              if (m < 0 || m >= imageHeight || n + k <0 || n +k >= imageWidth)
                              {
                                  
                                  rightPixel = 0;
                              }else
                              { 
                                  rightSrc = (unsigned char*)rightImage->imageData 
                                               + m*rightImage->widthStep + n;
                                  rightPixel = *rightSrc;
                                   
                              }
    
                              sum +=  leftPixel*rightPixel;
                              leftSquareSum  += leftPixel*leftPixel;
                              rightSquareSum += rightPixel*rightPixel;
    
                        }
                    }
                    matchLevel[k] = (double)sum/(sqrt(double(leftSquareSum))*sqrt((double)rightSquareSum));
                
                }
    
             
                tempMax = 0;
                tempIndex = 0;
                for ( m = 1;m < maxDBounds - minDBounds + 1;m++)
                {
                  
                    if (matchLevel[m] > matchLevel[tempIndex])
                    {
                        tempMax = matchLevel[m];
                        tempIndex = m;
                    }
                }
                dst = (unsigned char *)NCCImage->imageData + i*NCCImage->widthStep + j;
               
                *dst = tempIndex*8;
                dst = (unsigned char *)MatchLevelImage->imageData + i*MatchLevelImage->widthStep + j;
                *dst = (unsigned char)(tempMax*255);
              
            }
    
        }
       
        cvShowImage("Census",NCCImage);
        cvShowImage("MatchLevel",MatchLevelImage);
        cvSaveImage("depth.jpg",NCCImage);
        cvSaveImage("matchLevel.jpg",MatchLevelImage);
    
    
        cout<<endl<<"Over"<<endl;
        cvWaitKey(0);
        cvDestroyAllWindows();
        cvReleaseImage(&leftImage);
        cvReleaseImage(&rightImage);
        return 0;
    }

![1358165069_6291.jpg][]

//DP

//例如，X=“ABCBDAB”，Y=“BCDB”是X的一个子序列
    #include <cstdio>
    #include <cstring>
    #include <iostream>
    #include<cv.h>
    #include<highgui.h>
    #include <cmath>
    
    using namespace std;
    const int Width =  1024;
    const int Height = 1024;
     int Ddynamic[Width][Width];
    
    int  main()
    {  
        /*Half of the window size for the census transform*/
        int hWin = 11;
        int compareLength = (2*hWin+1)*(2*hWin+1);
    
        cout<<"hWin: "<<hWin<<";  "<<"compare length:  "<<compareLength<<endl;  
        cout<<"belief propagation test"<<endl;
    
        IplImage * leftImage = cvLoadImage("l2.jpg",0);
        IplImage * rightImage = cvLoadImage("r2.jpg",0);
    
       // IplImage * leftImage = cvLoadImage("left.bmp",0);                           
       // IplImage * rightImage = cvLoadImage("right.bmp",0);
    
        int imageWidth = leftImage->width;
        int imageHeight =leftImage->height;
    
      
    
        IplImage * DPImage = cvCreateImage(cvGetSize(leftImage),leftImage->depth,1);
        //IplImage * MatchLevelImage = cvCreateImage(cvGetSize(leftImage),leftImage->depth,1);
    
        unsigned char * pPixel = NULL;
        unsigned char  pixel;
        for (int i = 0; i< imageHeight;i++)
        {
            for (int j =0; j < imageWidth;j++ )
            {
                pPixel = (unsigned char *)DPImage->imageData + i*DPImage->widthStep + j;
                *pPixel = 0;
            }
        }
    
        int minDBounds = 0;
        int maxDBounds = 31;
    
        cvNamedWindow("Left",1);
        cvNamedWindow("Right",1);
        cvNamedWindow("Depth",1);
        cvNamedWindow("MatchLevel",1);
    
        cvShowImage("Left",leftImage);
        cvShowImage("Right",rightImage);
    
        int minD = 0;
        int maxD = 31;
        //假设图像是经过矫正的，那么每次都只是需要搜搜同一行的内容
        int max12Diff = 10;
       
        for (int i = 0;i < imageWidth;i++)
        {
            Ddynamic[0][i] = 0;
            Ddynamic[i][0] = 0;
        }
    
        unsigned char * pLeftPixel  = NULL;
        unsigned char * pRightPixel = NULL;
        unsigned char leftPixel = 0;
        unsigned char rightPixel =0;
        int m,n,l;
    
        for (int i = 0 ; i < imageHeight;i++)
        {
            for (int j = 0; j<imageWidth;j++)
            {
                for (int k = j + minD; k <= j + maxD;k++)
                {
                    if (k <0 || k >= imageWidth)
                    {
    
                    }else {
                        pLeftPixel = (unsigned char*)leftImage->imageData + i*leftImage->widthStep + k;
                        pRightPixel= (unsigned char*)rightImage->imageData+i*rightImage->widthStep + j;
                        leftPixel  = *pLeftPixel;
                        rightPixel = *pRightPixel;
    
                        if (abs(leftPixel - rightPixel) <= max12Diff)
                        {
                            Ddynamic[j + 1][k + 1] = Ddynamic[j][k] +1; 
                        }else if (Ddynamic[j][k+1] > Ddynamic[j+1][k])
                        {
                            Ddynamic[j + 1][k + 1] = Ddynamic[j][k+1];
                        }else{
                            Ddynamic[j+1][k+1] = Ddynamic[j+1][k];
                        }
                        
                        //cout<<Ddynamic[j +1][k+1]<<"  ";
                    }
                   
                }
                 //cout<<"\n";
            }
            //逆向搜索，找出最佳路径
             m = imageWidth;
             n = imageWidth;
             l = Ddynamic[imageWidth][imageWidth];
            while( l>0 )
            {
                if (Ddynamic[m][n] == Ddynamic[m-1][n])  
                    m--;
                else if (Ddynamic[m][n] == Ddynamic[m][n-1])  
                    n--;
                else
                { 
                    //s[--l]=a[i-1];
                    pPixel = (unsigned char *)DPImage->imageData + i*DPImage->widthStep + m;
                    *pPixel = (n-m)*8;
                    l--;
                    m--; 
                    n--;
                }
            }
    
           //cvWaitKey(0);
    
        }
    
        cvShowImage("Depth",DPImage);
        cvSaveImage("depth.jpg",DPImage);
        cvWaitKey(0);
        return 0;
    }

![1358165116_1552.jpg][]

//贴一张清楚的

![1358169390_5203.png][]

//DP5

//引入概率公式
    //
    #include <cstdio>
    #include <cstring>
    #include <iostream>
    #include<cv.h>
    #include<highgui.h>
    #include <cmath>
    
    using namespace std;
    const int Width =  1024;
    const int Height = 1024;
     int Ddynamic[Width][Width];
    
     //使用钟形曲线作为匹配概率，差值越小则匹配的概率越大，最终的要求是使匹配的概率最大，概率曲线使用matlab生成
     int Probability[256] = {
        255, 255, 254, 252, 250, 247, 244, 240, 235, 230, 225, 219, 213, 206, 200, 192, 185, 178, 170, 162, 
        155, 147, 139, 132, 124, 117, 110, 103, 96, 89, 83, 77, 71, 65, 60, 55, 50, 46, 42, 38, 35, 31, 28, 
        25, 23, 20, 18, 16, 14, 13, 11, 10, 9, 8, 7, 6, 5, 4, 4, 3, 3, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
     };
    
    int  main()
    {  
      
        IplImage * leftImage = cvLoadImage("l2.jpg",0);
        IplImage * rightImage = cvLoadImage("r2.jpg",0);
    
        //IplImage * leftImage = cvLoadImage("left.bmp",0);                           
        //IplImage * rightImage = cvLoadImage("right.bmp",0);
    
        int imageWidth = leftImage->width;
        int imageHeight =leftImage->height;
    
        IplImage * DPImage = cvCreateImage(cvGetSize(leftImage),leftImage->depth,1);
        IplImage * effectiveImage = cvCreateImage(cvGetSize(leftImage),leftImage->depth,1);
        IplImage * FilterImage = cvCreateImage(cvGetSize(leftImage),leftImage->depth,1);
    
        unsigned char * pPixel = NULL;
        unsigned char  pixel;
        unsigned char * pPixel2 = NULL;
        unsigned char  pixel2;
        for (int i = 0; i< imageHeight;i++)
        {
            for (int j =0; j < imageWidth;j++ )
            {
                pPixel = (unsigned char *)DPImage->imageData + i*DPImage->widthStep + j;
                *pPixel = 0;
                pPixel = (unsigned char *)effectiveImage->imageData + i*effectiveImage->widthStep + j;
                *pPixel = 0;
            }
        }
    
     
    
        cvNamedWindow("Left",1);
        cvNamedWindow("Right",1);
        cvNamedWindow("Depth",1);
        cvNamedWindow("effectiveImage",1);
    
        cvShowImage("Left",leftImage);
        cvShowImage("Right",rightImage);
    
        int minD = 0;
        int maxD = 31;
        //假设图像是经过矫正的，那么每次都只是需要搜搜同一行的内容
        int max12Diff = 5;
       
        for (int i = 0;i < imageWidth;i++)
        {
            Ddynamic[0][i] = 0;
            Ddynamic[i][0] = 0;
        }
    
        unsigned char * pLeftPixel  = NULL;
        unsigned char * pRightPixel = NULL;
        unsigned char leftPixel = 0;
        unsigned char rightPixel =0;
        int m,n,l;
    
        int t1 = clock();
        for (int i = 0 ; i < imageHeight;i++)
        {
            for (int j = 0; j<imageWidth;j++)
            {
                for (int k = j + minD; k <= j + maxD;k++)
                {
                    if (k <0 || k >= imageWidth)
                    {
    
                    }else {
                        pLeftPixel = (unsigned char*)leftImage->imageData + i*leftImage->widthStep + k;
                        pRightPixel= (unsigned char*)rightImage->imageData+i*rightImage->widthStep + j;
                        leftPixel  = *pLeftPixel;
                        rightPixel = *pRightPixel;
                        //之前概率最大的点加上当前的概率
    
                        Ddynamic[j + 1][k + 1] = max(Ddynamic[j][k],max(Ddynamic[j][k+1],Ddynamic[j+1][k]))
                                                 + Probability[abs(leftPixel - rightPixel)];
                        /* if (abs(leftPixel - rightPixel) <= max12Diff)
                        {
                        Ddynamic[j + 1][k + 1] = Ddynamic[j][k] +1; 
                        }else if (Ddynamic[j][k+1] > Ddynamic[j+1][k])
                        {
                        Ddynamic[j + 1][k + 1] = Ddynamic[j][k+1];
                        }else{
                        Ddynamic[j+1][k+1] = Ddynamic[j+1][k];
                        }*/
                        
                        //cout<<Ddynamic[j +1][k+1]<<"  ";
                    }
                   
                }
                 //cout<<"\n";
            }
            //逆向搜索，找出最佳路径
             m = imageWidth;
             n = imageWidth;
             l = Ddynamic[imageWidth][imageWidth];
            while( m >= 1 && n >= 1)
            {
                pPixel = (unsigned char *)DPImage->imageData + i*DPImage->widthStep + m;
                *pPixel = (n-m)*8;
                //标记有效匹配点
                pPixel = (unsigned char *)effectiveImage->imageData + i*effectiveImage->widthStep + m;
                *pPixel = 255;
                if (Ddynamic[m-1][n] >= Ddynamic[m][n -1] && Ddynamic[m-1][n] >= Ddynamic[m-1][n -1])  
                    m--;
                else if (Ddynamic[m][n-1] >= Ddynamic[m-1][n] && Ddynamic[m][n -1] >= Ddynamic[m-1][n -1])  
                    n--;
                else
                { 
                    //s[--l]=a[i-1];       
                   // l -= Ddynamic[m][n];
                    m--; 
                    n--;
                }
            }
    
           //cvWaitKey(0);
        }
    
        //refine the depth image  7*7中值滤波
        //统计未能匹配点的个数
        int count = 0;
        for (int i = 0 ;i< imageHeight;i++)
        {
            for (int j= 0; j< imageWidth;j++)
            {
                pPixel = (unsigned char *)effectiveImage->imageData + i*effectiveImage->widthStep + j;
                pixel = *pPixel;
                if (pixel == 0)
                {
                    count++;
                }
            }
        }
        int t2 = clock();
        cout<<"dt: "<<t2-t1<<endl;
        cout<<"Count:  "<<count<<"  "<<(double)count/(imageWidth*imageHeight)<<endl;
        cvShowImage("Depth",DPImage);
        cvShowImage("effectiveImage",effectiveImage);
       // cvWaitKey(0);
    
    
        FilterImage = cvCloneImage(DPImage);
    
        //7*7中值滤波
        int halfMedianWindowSize = 3;
        int medianWindowSize = 2*halfMedianWindowSize + 1;
        int medianArray[100] = {0};
        count = 0;
        int temp = 0;
        int medianVal = 0;
    
        for (int i = halfMedianWindowSize + 1 ;i< imageHeight - halfMedianWindowSize;i++)
        {
            for (int j = halfMedianWindowSize; j< imageWidth - halfMedianWindowSize;j++)
            {
                pPixel = (unsigned char *)effectiveImage->imageData + i*effectiveImage->widthStep + j;
                pixel = *pPixel;
                if (pixel == 0)
                {
                    count = 0;
                    for (int m = i - halfMedianWindowSize ; m <= i + halfMedianWindowSize ;m++)
                    {
                        for (int n = j - halfMedianWindowSize; n <= j + halfMedianWindowSize ;n++)
                        {
                            pPixel2 = (unsigned char *)DPImage->imageData + m*DPImage->widthStep + n;
                            pixel2 = *pPixel2;
                            if (pixel2 != 0)
                            {
                                 medianArray[count] = pixel2;
                                 count++;
                            }
           
                        }
                        //排序
                        for (int k = 0; k< count;k++)
                        {
                            for (int l = k + 1; l< count;l++)
                            {
                                if (medianArray[l] < medianArray[l-1] )
                                {
                                    temp = medianArray[l];
                                    medianArray[l] = medianArray[l-1];
                                    medianArray[l-1] = temp;
                                }
                            }
                        }
                        medianVal = medianArray[count/2];
                        pPixel = (unsigned char *)FilterImage->imageData + i*DPImage->widthStep + j;
                        *pPixel = medianVal;
                    }
    
                }
            }
        }
        cvShowImage("Depth",DPImage);
        cvShowImage("effectiveImage",effectiveImage);
        cvShowImage("Filter",FilterImage);
        cvSaveImage("depth.jpg",DPImage);
        cvSaveImage("effective.jpg",effectiveImage);
    
        cvWaitKey(0);
        return 0;
    }

![1358165176_9196.jpg][]

[给我写信][Link 1]

[http_research.microsoft.com_zhang]: http://research.microsoft.com/~zhang/
[1358164779_2048.jpg]: https://img-my.csdn.net/uploads/201301/14/1358164779_2048.jpg
[1358164806_3812.jpg]: https://img-my.csdn.net/uploads/201301/14/1358164806_3812.jpg
[1358164851_3510.jpg]: https://img-my.csdn.net/uploads/201301/14/1358164851_3510.jpg
[1358164921_6648.jpg]: https://img-my.csdn.net/uploads/201301/14/1358164921_6648.jpg
[1358164982_3041.jpg]: https://img-my.csdn.net/uploads/201301/14/1358164982_3041.jpg
[1358165030_2256.jpg]: https://img-my.csdn.net/uploads/201301/14/1358165030_2256.jpg
[1358165069_6291.jpg]: https://img-my.csdn.net/uploads/201301/14/1358165069_6291.jpg
[1358165116_1552.jpg]: https://img-my.csdn.net/uploads/201301/14/1358165116_1552.jpg
[1358169390_5203.png]: https://img-my.csdn.net/uploads/201301/14/1358169390_5203.png
[1358165176_9196.jpg]: https://img-my.csdn.net/uploads/201301/14/1358165176_9196.jpg
[Link 1]: http://mail.qq.com/cgi-bin/qm_share?t=qm_mailme&email=ckVHREZGR0tAQDIDA1wRHR8