一种面向光场图像的深度估计方法与流程

本发明涉及计算机视觉与数字图像处理技术领域，尤其是一种面向光场图像的深度估计方法。

背景技术：

20世纪上半叶，a.gershun等人已经提出了光场概念，用以描述光在三维空间中的辐射特性，然而，光场成像技术相对滞后于理论概念的发展。陆续出现的光场成像设备主要有相机阵列、光场传感器、微透镜阵列等。r.ng等人于2005年设计的手持式光场相机，成为微透镜阵列光场采集设备的代表，直接推动了光场相机(如lytro、raytrix等)快速进入民用市场，同时也在学术界掀起了光场成像的应用研究热潮。

基于图像分析方法获取场景深度信息(物体与成像系统的距离)是计算机视觉领域的一个重要问题，在机器人导航、三维重建、摄影测量等方面有着广泛应用。与传统相机相比，基于微透镜阵列的光场相机能够同时记录光线的位置和方向，由一次曝光采集的光场数据既可以计算得到多视角的子孔径图像又可以生成焦点堆栈图像，为深度估计提供了便利条件。近年来，基于光场成像的深度估计成为cvpr、iccv等计算机视觉顶级会议讨论的话题之一，不少学者开展了该问题的研究：尹晓艮等人直接使用halcon软件的聚焦深度函数进行对焦点堆栈图像进行深度获取，并应用于3d重建。lin等人对编码焦点堆栈进行了系统研究，并利用稀疏编码技术从焦点堆栈中恢复全聚焦图像和全分辨率深度数据。wanner等人提出全局一致性深度标记算法，在全变分框架下，利用epi图像梯度的结构张量求解纹理斜率，并使用最小化能量函数进行全局约束，实现深度信息的估计。suzuki等人就视差范围受限问题，提出在结构张量分析之前先对epi图像进行裁剪，对不同裁剪系数的结构张量分析结果综合生成最终视差图。

中国专利cn104899870a，公布日2015.09.09，申请号201510251234.8，基于光场数据分布的深度估计方法。利用焦点堆栈图像，分别提取同一个宏像素的强度范围，进而选出最小强度范围对应的重聚焦光场图像，以该重聚焦光场图像的焦距作为该宏像素的场景深度。该发明利用lytro等光场相机采集的场景纹理及空间信息，获得细节丰富，特征清晰，准确性、一致性高的场景深度估计。中国专利cn104966289a，公布日2015.10.07，申请号201510324490.5，一种基于4d光场的深度估计方法。该发明重点解决了初始深度值的可靠性问题、深度过平滑问题、深度值融合问题。中国专利cn105551050a，公布日2016.05.04，申请号201511019609.4，一种基于光场的图像深度估计方法。该发明检测显著特征点并将显著特征点匹配作为深度优化的强约束条件，避免了插值过程引起的视点混淆，提升了深度估计准确性。中国专利cn106651943a，公布日2017.05.10，申请号201611262452.2，一种基于遮挡几何互补模型的光场相机深度估计方法。该发明利用相对的两个视角中，所看到遮挡区域具有互补关系的性质，实现较好的遮挡结果和准确深度图。

当前，对光场相机这类新型成像设备的认识以及基于光场成像的深度估计方法的研究均处于新兴探索阶段。现有基于光场的深度估计方法较少关注深度估计对噪声和遮挡的鲁棒性，而这对复杂场景的深度估计是至关重要的。

技术实现要素：

本发明所要解决的技术问题在于，提供一种面向光场图像的深度估计方法，能够由光场数据得到准确的深度估计结果，并且对存在噪声和遮挡的复杂场景具有较高鲁棒性。

为解决上述技术问题，本发明提供一种面向光场图像的深度估计方法，包括如下步骤：

(1)根据光场相机参数信息，将光场原始图像解码为4d光场数据，并提取中心子孔径图像其中(i^c,j^c)表示中心子孔径图像的视角坐标；

(2)由步骤(1)得到的4d光场数据，以视差δd为步长，计算生成焦点堆栈图像序列其中，n表示深度估计设定的视差标签数，根据光场相机的视差范围参数[dmin,dmax]，求得

(3)对中心子孔径图像上的每个位置坐标(k^*,l^*)处的像素p，以视差(γ·δd，γ∈{1,2,…,n})为横坐标，该视差对应的重聚焦图像上(k^*,l^*)处像素的灰度值为纵坐标，绘制像素p对应的聚焦曲线φp(γ·δd)；

(4)对中心子孔径图像上的每个位置坐标(k^*,l^*)处的像素p，检测步骤(3)得到的聚焦曲线(φp(γ·δd))的局部对称性，取最优局部对称中心的横坐标值作为像素p对应的视差估计值(d1p)，针对每个像素执行此处理后，得到视差图d1；

(5)对中心子孔径图像上的每个位置坐标(k^*,l^*)处的像素p，由步骤(1)得到的4d光场数据，分别计算生成关于p的水平、垂直方向的epi图像

(6)对中心子孔径图像上的每个位置坐标(k^*,l^*)处的像素p，在步骤(5)得到的epi图像上分别检测对应直线的斜率sh、sv，由此得到两个视差值，求两者均值作为p的估计视差d2p；对每个像素进行估计视差处理后得到视差图d2；

(7)对中心子孔径图像上的每个位置(k^*,l^*)处的像素p，由和d1绘制聚焦曲线片段再由和d2绘制聚焦曲线片段其中d1p、d2p分别表示视差图d1、d2上p点对应位置处的视差，记号[(a):(b):(c)]表示在a到c范围内以b为步长生成的等差数列；

(8)对中心子孔径图像上的每个位置坐标(k^*,l^*)处的像素p，在集合上计算聚焦曲线片段与聚焦曲线φp(γ·δd)的匹配度同样，在集合上计算聚焦曲线片段与聚焦曲线φp(γ·δd)的匹配度若则p对应的最终深度dp赋值为d1p，否则，dp赋值为d2p；求得每个像素的视差后，得到最终视差图d。

优选的，步骤(1)中，4d光场记为l:(i,j,k,l)→l(i,j,k,l)，其中，(i,j)表示微透镜图像的像素索引坐标，(k,l)表示微透镜中心的索引坐标，l(i,j,k,l)表示通过像素和微透镜中心的光线的辐射强度；提取中心子孔径图像的方法是抽取每个微透镜图像的中心像素，按微透镜位置索引排列得到二维图像，即

优选的，步骤(2)中，生成焦点堆栈图像序列具体为：基于光场数字重聚焦技术，即把光场重新投影积分到新的像平面上，等效于所有子孔径图像的平移和叠加，计算公式如下：

其中，k'、l'分别表示重聚焦图像像素在水平、垂直方向的坐标；ni、nj分别表示微透镜所覆盖子图像在行列方向上的像素数；nk、nl分别表示微透镜阵列在水平、垂直方向上的微透镜个数；interpn()表示常见的插值函数，返回在指定位置的插值，第1个参数是一个二维矩阵，第2参数是插入位置的水平坐标，第3个参数是插入位置的垂直坐标；β是比例因子；当k'∈[1:1:nk]，l'∈[1:1:nl]取遍集合中每个整数值时，得到重聚焦图像当γ取遍集合{1,2,…,n}时，求得重聚焦图像序列

优选的，步骤(4)中，检测聚焦曲线的局部对称性时，局部范围是指取一个窗口w，窗口大小w设定为对聚焦曲线(φp(γ·δd))上横坐标对应的每个位置，在区间范围内计算关于γ·δd镜像对称位置的灰度差异，作为该处的镜像对称性度量mp(γ·δd)，mp(γ·δd)越小，局部对称性越好；具体计算方法如下：

这里，ρ()是常用距离函数，定义为：

其中σ是控制因子。

优选的，步骤(5)中，由4d光场数据生成水平方向和垂直方向epi图像的过程可以看作如下映射：即固定4d光场中的两个维度，变化另外两个维度所得到的二维切片图像；固定i、k维度时，令i＝i^*＝i^c，k＝k^*；固定j、l维度时，令j＝j^*＝j^c，l＝l^*。

优选的，步骤(6)中，在上检测直线斜率的方法是模板匹配法，具体包括如下步骤：

(61)在角度区间[1°,179°]上，等步长生成n个元素的角度序列其第m个元素记为θm，度数为

(62)在的图像坐标系下，生成经过点(l^*,j^c)，且斜率为tanθm的直线，选择该直线经过的所有像素，构成集合sp(θm)；

(63)计算步骤(62)所得直线的灰度一致度cp(θm)，计算方法是：

其中，|sp(θm)|表示集合sp(θm)的基数，即集合中元素的个数；

(64)对于重复步骤(62)和(63)，计算一致度cp(θ)；取使一致度达到最小的角度作为p所在直线的倾斜角度即直线的斜率为

采用与(61)～(64)类似的步骤，可以在epi图像上计算直线斜率sv；进一步得到两个估计视差值计算公式是其中b是微透镜间距，由相机参数提供或标定得到，p的估计视差为

优选的，步骤(7)中，由中心子孔径图像和视差图d1绘制聚焦曲线片段的方法具体包括如下步骤：

(71)对于重复步骤(72)～(73)，得到p在视差为d1p附近的聚焦曲线片段

(72)在中心子孔径图像上，在以像素p为中心(位置坐标记为(k^*,l^*))、边长为w的窗口区域内，选择满足下列条件的像素q(位置坐标记为(k^q,l^q))：这里ε是误差门限值，通常设定为将所有满足条件的像素q记入集合cv；

(73)计算集合cv中的所有像素的均值，作为聚焦曲线片段上横坐标为视差f时的纵坐标值。

优选的，步骤(8)中，计算与φp(γ·δd)匹配度是在共同的视差集合上进行的，即在聚焦曲线φp(γ·δd)上取一个片段(即γ·δd＝f)与进行一一比较，将集合简记为a，匹配度的计算方法是：

同理，可计算匹配度

优选的，以步骤(8)得到的视差图d为优化的基础，根据聚焦曲线的匹配度建立数据项约束edata(d)，以视差连续性为指导增加平滑项约束esmooth(d)，由此构建能量函数e(d)＝edata(d)+esmooth(d)，采用迭代优化技术求解使能量函数取得最小的视差图作为最终视差图。

优选的，数据项约束的定义为：

其中，λ1是控制因子；

平滑项约束的定义为：

其中，n表示中心子孔径图像或视差图上所有相邻位置构成的集合，{p,q}表示一对相邻位置，λ、λ2、σ表示控制因子。

本发明的有益效果为：(1)本发明针对光场图像提出了一种深度估计方法，该方法使聚焦和视差两种方式估计的视差结果相互补充，提高了深度估计的准确性；(2)本发明的聚焦方式深度估计基于光场的数字重聚焦，重聚焦过程是积分过程，降低了噪声对的影响；另外，视差方式深度估计以水平方向和垂直方向的视差估计的均值作为结果，提高了对噪声的鲁棒性；(3)本发明结合中心子孔径图像和视差图生成聚焦曲线片段，并用该片段与聚焦曲线的匹配度作为视差值取舍的依据，提高了深度估计对遮挡的鲁棒性。

附图说明

图1为本发明的方法流程示意图。

图2为本发明中的中心子孔径图像示例示意图。

图3为本发明中的聚焦于鹦鹉的重聚焦图像示例示意图。

图4为本发明中的聚焦于建筑物的重聚焦图像示例示意图。

图5为本发明中的聚焦曲线示例示意图。

图6为本发明中的epi图像示例示意图。

图7为本发明中放大的epi图像示例示意图。

图8为本发明中在epi图像上选择像素的示意图。

图9为本发明与现有方法的深度估计结果对比示例示意图。

具体实施方式

如图1所示，一种面向光场图像的深度估计方法，包括如下步骤：

步骤1，解码光场图像得到4d光场数据、中心子孔径图像。即：根据光场相机参数信息，将光场原始图像解码为4d光场数据，并提取中心子孔径图像

如步骤1所述，光场相机(如lytro等)所采集的原始图像通常是12位bayer格式图像，而在光场分析和处理中通常采用2pp(twoplaneparameterization)形式表示的4d光场，因此，需要根据光场相机提供的参数或进一步标定修正的参数，将原始图像解码为4d光场数据。4d光场通常记为l:(i,j,k,l)→l(i,j,k,l)，其中，(i,j)表示微透镜所覆盖子图像的像素索引坐标，(k,l)表示微透镜中心的索引坐标，l(i,j,k,l)表示通过像素和微透镜中心的光线的辐射强度。提取中心子孔径图像的方法是抽取每个微透镜图像的中心像素，按微透镜位置索引排列得到二维图像，即这里，(i^c,j^c)表示中心子孔径图像的视角坐标。例如每个微透镜覆盖了9×9的微透镜图像，微透镜阵列的行列数均为380，那么i^c＝j^c＝5，将每个微透镜图像的(5,5)处的像素抽取出来，按照其所在微透镜的索引坐标组合成分辨率为380×380的图像即是中心子孔径图像图2给出了lytro相机采集的光场数据得到中心子孔径图像示例。

步骤2，由4d光场数据生成焦点堆栈图像序列。具体地，由步骤1得到的4d光场数据，以视差δd为步长，计算生成焦点堆栈图像序列其中，n表示深度估计设定的视差标签数，根据光场相机的视差范围参数[dmin,dmax]，求得

步骤2所述的生成焦点堆栈图像序列的方法基于光场数字重聚焦技术，即把光场重新投影积分到新的像平面上，等效于所有子孔径图像的平移和叠加。计算公式如下：

其中，k'、l'分别表示重聚焦图像像素在水平、垂直方向的坐标；ni、nj分别表示微透镜覆盖的子图像(微透镜图像)在行列方向上的像素数；nk、nl分别表示微透镜阵列水平、垂直方向上的微透镜个数；interpn()是常见的插值函数，返回在指定位置的插值，第1个参数是一个二维矩阵，第2参数是插入位置的水平坐标，第3个参数是插入位置的垂直坐标；β是比例因子。

当k'∈[1:1:nk]，l'∈[1:1:nl]取遍集合中每个整数值时，得到重聚焦图像当γ取遍集合{1,2,…,n}时，求得重聚焦图像序列例如本实施例中，设定深度估计有256个视差标签，即令n＝256，则生成的重聚焦图像序列中有256幅图像，图3、图4是同一场景的两幅光场重聚焦图像示例，图3聚焦于鹦鹉，图4聚焦于建筑物。

步骤3，对中心子孔径图像上的每个像素，绘制聚焦曲线。具体地，对中心子孔径图像上的每个位置坐标(k^*,l^*)处的像素p，以视差(γ·δd，γ∈[1,n])为横坐标，该视差对应的重聚焦图像上(k^*,l^*)处像素的灰度值为纵坐标，绘制像素p对应的聚焦曲线φp(γ·δd)。图5给出了一幅聚焦曲线的示例。

步骤4，检测聚焦曲线的局部对称中心，得视差图d1。具体地，对中心子孔径图像上的每个位置坐标(k^*,l^*)处的像素p，检测步骤3得到的聚焦曲线(φp(γ·δd))的局部对称性，取最优局部对称中心的横坐标值作为像素p对应的视差估计值(d1p)，针对每个像素执行此处理后，得到视差图d1。

步骤4所述的检测聚焦曲线的局部对称性时，局部范围是指取一个窗口w，窗口大小|w|通常设定为对聚焦曲线(φp(γ·δd))上横坐标对应的每个位置，在区间范围内计算关于γ·δd镜像对称位置的灰度差异，作为该处的镜像对称性度量mp(γ·δd)，mp(γ·δd)越小，局部对称性越好；具体计算方法如下：

这里，ρ()是常用距离函数，定义为：

其中σ是控制因子。

步骤5，由4d光场数据生成每个像素对应的水平方向、垂直方向的epi图像。具体地，对中心子孔径图像上的每个位置坐标(k^*,l^*)处的像素p，由步骤1得到的4d光场数据，分别计算生成关于p的水平、垂直方向的epi图像

步骤5所述的由4d光场数据生成水平方向和垂直方向epi图像的过程可以看作如下映射：即固定4d光场中的两个维度，变化另外两个维度所得到的二维切片图像；固定i、k维度时，令i＝i^*＝i^c，k＝k^*；固定j、l维度时，令j＝j^*＝j^c，l＝l^*。图6给出了epi图像示例，该epi图像是图2中虚线位置对应的epi图像。

步骤6，检测epi图像上直线的斜率，得到视差图d2。具体地，对中心子孔径图像上的每个位置(k^*,l^*)处的像素p，在步骤5得到的epi图像上分别检测对应直线的斜率sh、sv，由此得到两个视差值，求两者均值作为p的估计视差d2p；对每个像素进行估计视差处理后得到视差图d2。

步骤6所述的在上检测直线斜率的方法是模板匹配法，具体步骤包括：

(601)在角度区间[1°,179°]上，等步长生成n个元素的角度序列其第m个元素记为θm，度数为

(602)在的图像坐标系下，生成经过点(l^*,j^c)，且斜率为tanθm的直线，选择该直线经过的所有像素，构成集合sp(θm)；

(603)计算步骤(602)所得直线的灰度一致度cp(θm)，计算方法是：

其中，|sp(θm)|表示集合sp(θm)的基数，即集合中元素的个数。

(604)对于重复步骤(602)和(603)，计算一致度cp(θ)；取使一致度达到最小的角度作为p所在直线的倾斜角度即直线的斜率为

采用与(601)～(604)类似的步骤，可以在epi图像上计算直线斜率sv。进一步得到两个估计视差值计算公式是其中b是微透镜间距，由相机参数提供或标定得到。p的估计视差为

如图7是一幅epi图像的放大示意图，每个格子表示一个像素；图8给出在epi图像上选择像素集合的示例，根据假定斜率的直线所通过的像素，选择像素集合。

步骤7，基于中心子孔径图像和两个视差图，对每个像素绘制两种聚焦曲线片段。具体地，对中心子孔径图像上的每个位置坐标(k^*,l^*)处的像素p，由和d1绘制聚焦曲线片段再由和d2绘制聚焦曲线片段其中d1p、d2p分别表示视差图d1、d2上p点对应位置处的视差，记号[(a):(b):(c)]表示在a到c范围内以b为步长生成的等差数列。

步骤7所述的由中心子孔径图像和视差图d1绘制聚焦曲线片段的方法包括以下步骤：

(701)对于重复步骤(702)～(703)，得到p在视差为d1p附近的聚焦曲线片段

(702)在中心子孔径图像上，在以像素p为中心(位置坐标记为(k^*,l^*))、边长为w的窗口区域内，选择满足下列条件的像素q(位置坐标记为(k^q,l^q))：这里ε是误差门限值，通常设定为将所有满足条件的像素q记入集合cv；

(703)计算集合cv中的所有像素的均值，作为聚焦曲线片段上横坐标为视差f时的纵坐标值。

由中心子孔径图像和视差图d2绘制聚焦曲线片段的方法与绘制的方法类似。

步骤8，计算两种聚焦曲线片段与聚焦曲线的匹配度，由匹配度确定最终视差图d。具体地，对中心子孔径图像上的每个位置(k^*,l^*)处的像素p，在集合上计算聚焦曲线片段与聚焦曲线φp(γ·δd)的匹配度同样，在集合上计算聚焦曲线片段与聚焦曲线φp(γ·δd)的匹配度若则p对应的最终深度dp赋值为d1p，否则，dp赋值为d2p；求得每个像素的视差后，得到最终视差图d。

步骤8所述的计算与φp(γ·δd)匹配度是在共同的视差集合上进行的，即在聚焦曲线φp(γ·δd)上取一个片段(即γ·δd＝f)与进行一一比较，将集合简记为a，匹配度的计算方法是：

同理，可计算匹配度

在上述的基于光场图像的深度估计方法中，优选地，进一步还包括对以上步骤得到的结果进行优化的步骤9。

步骤9，优选地，执行基于能量最小化的深度优化。具体地，以步骤8得到的视差图d为优化的基础，根据聚焦曲线的匹配度建立数据项约束edata(d)，以视差连续性为指导增加平滑项约束esmooth(d)，由此构建能量函数e(d)＝edata(d)+esmooth(d)，采用迭代优化技术求解使能量函数取得最小的视差图作为最终视差图。

步骤9所述的数据项约束的定义为：

其中，λ1是控制因子；

平滑项约束的定义为：

其中，n表示中心子孔径图像或视差图上所有相邻位置构成的集合，{p,q}表示一对相邻位置，λ、λ2、σ表示控制因子。

本实施例中，采用图割法优化该能量函数，得到优化的视差图。

本实施例中两种典型的光场图像深度估计方法与本发明进行对比分析，一种是tao方法：tao于2013年提出，发表于cvpr；一种是jeon方法：jeon于2015年提出，发表于cvpr。

本实施例使用mousnier等提供的lytro光场数据集测试发明对真实场景光场数据的性能。该数据集包括30组lytro摄像数据，其中，25组是各地的室内外场景，3组包含运动模糊，1组包含长时间曝光，1组平面摄影。图9给出了其中5组典型场景的深度估计结果，第一列是场景的中心子孔径图像，第二至四列分别是本发明的结果、tao方法的结果、jeon方法的结果；上面两行展示两组室内场景、中间一行是局部细节场景、下面两行是典型室外场景。

从本实施例的分析可以明显看出，本发明不论对室外噪声场景还是对复杂遮挡的细节场景，均能较好估计深度信息。