CN110009676B - 一种双目图像的本征性质分解方法 - Google Patents

Abstract

本发明公开了一种双目图像的本征性质分解方法，该方法以双目图像中的明暗分量为待求的未知量，结合本发明提供的Retinex约束、反射分量一致性约束、局部约束、非局部约束以及绝对尺度约束等构造一个二次型的目标函数，从而双目图像的本征性质分解问题转化成为一个求解二次型代价函数能量最小解的问题。经过理论推导和实例验证，该技术方案具有分解效果好，算法效率高，容易实现等优点。

Description

一种双目图像的本征性质分解方法

技术领域

本发明属于计算机视觉和计算机图形学领域，具体涉及一种用于双目图像的本征性质分解方法。

背景技术

图像的本征性质分解是将一副原始的输入图像分解成对应的明暗分量(ShadingComponent)和反射分量(Reflectance Component)，在计算机视觉和计算机图形学中这是一项非常重要的底层技术；这项技术将有助于提高计算机视觉中许多算法的准确性，例如图像的分割、识别、检测等。同样的，图像的本征分解算法也有利于计算机图形学中的渲染技术、场景的几何形状建模等。

随着3D电影在影视界的成功应用，双目相机和立体显示设备变得越来越流行。近年来。双目图像和视频的处理方法和技术逐渐成为学术及和产业界的热门课题；例如，双目图像的扭曲变换(Stereoscopic warping)、克隆(Cloning)、全景(Panorama)、以及防抖(Stabilization)等应用。然而，针对双目图像的本征性质分解问题至今尚未出现非常有效的方法。

发明内容

针对现有技术中的上述不足，本发明提供的双目图像本征性质分解方法解决了现有的相关技术中，没有对能有效地对双目图像进行本征性质分解的问题。

为了达到上述发明目的，本发明采用的技术方案为：一种双目图像的本征性质分解方法，包括以下步骤：

S1、通过双目相机获得待进行本征性质分解的双目图像；

S2、确定双目图像本征性质分解时的约束项；

S3、根据约束项确定双目图像本征性质分解的目标函数；

S4、确定目标函数的最优解，进而获得双目图像的明暗分量

S5、根据获得的双目图像的明暗分量，获得双目图像左右视图的明暗分量和反射分量，实现双目图像的本征性质分解。

本发明的有益效果为：本发明提供了一种双目图像的本征性质分解方法以双目图像中的明暗分量为待求的未知量，结合本发明提供的Retinex约束、反射分量一致性约束、局部约束、非局部约束以及绝对尺度约束等构造一个二次型的目标函数，从而双目图像的本征性质分解问题转化成为一个求解二次型代价函数能量最小解的问题。经过理论推导和实例验证，该技术方案具有分解效果好，算法效率高，容易实现等优点。

附图说明

图1为本发明中双目图像Overlap Region示意图。

图2为本发明提供的双目图像本征性质分解方法流程图。

图3为本发明中局部约束条件对应的效果示意图。

图4为本发明中非局部约束条件对应的效果示意图。

图5为通过本发明方法进行双目图像本征性质分解的示例图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

图像的本征性质分解处理可以用如下公式进行描述：

I(x,y)＝R(x,y)·S(x,y) (1)

式中，I表示输入图像；

R表示本征性质分解方法所得的反射分量；

S表示本征性质分解方法所得的明暗分量；

(x,y)表示图像的像素点坐标；

符号·表示反射分量和明暗分量对应像素点值的乘法运算；

双目图像的本征性质分解方法则是分别对双目相机获得的左右两幅图像同时做本征性质分解，得到左右两幅图所对应的反射分量和明暗分量，即：

I_left(x,y)＝R_left(x,y)·S_left(x,y) (2)

I_right(x,y)＝R_right(x,y)·S_right(x,y) (3)

式中，下标left和right分别为双目图像的左右图像标识符；

因双目图像中的左右两张图像是有一定关联的。一般而言，双目图像的左右两图之间存在一定的左右视角重叠区域，即Overlap Region，如图1所示。

结合公式(2)，(3)，从双目图像本征性质的分解目标可以看出，对双目图像本征性质的分解是一个严重的欠约束问题(Under Constrained Problem)。原因在于，已知量有I_left(x,y)，I_right(x,y)，而待求解的未知量是R_left(x,y)，S_left(x,y)，R_right(x,y)和R_right(x,y)。未知量的个数是已知量的个数的两倍。对欠约束问题的求解，关键的因素在于如何寻找到合理的约束条件，使得上述欠约束的问题的解以最大的概率逼近最优解。

为了解决上述问题，本发明提供了一种双目图像的本征性质分解方法，如图2所示，包括以下步骤：

S1、通过双目相机获得待进行本征性质分解的双目图像；

S2、确定双目图像本征性质分解时的约束项；

S3、根据约束项确定双目图像本征性质分解的目标函数；

S4、确定目标函数的最优解，进而获得双目图像的明暗分量

S5、根据获得的双目图像的明暗分量，获得双目图像左右视图的明暗分量和反射分量，实现双目图像的本征性质分解。

上述步骤S2中双目图像本征性质分解时的约束项包括Retinex约束、一致性约束、局部约束、非局部约束和绝对尺度约束；

(1)Retinex约束

Retinex约束作为图像本征性质分解最重要的约束，它可以比较合理地区分引起图像像素值变化的原因，即图像像素值变化主要是由于场景反射分量的变化还是由于明暗分量的变化。Retinex是基于这样的假设：在一副图像的梯度图像中，当梯度的变化超过一定的阈值时，这种变化归结于反射分量的变化；而当梯度的变化较小时，则是由于明暗分量的变化所引起的。为此，本发明设置有如下的Retinex约束项：

式中，m，n分别表示双目图像中的两个任意像素点；

表示像素点m的一个预设领域，例如4；

ω_m,n表示相邻的像素点m，n的色度距离；

ΔS_m,n表示明暗分量差量，且ΔS_m,n＝S_m-S_n，S_m和S_n分别表示双目图像中像素点m和像素点n的明暗分量值；

ΔI_m,n表示像素差量，且ΔI_m,n＝I_m-I_n，I_m和I_n分别表示双目图像中像素点m和像素点n的像素值；

式中，ω_m,n为相邻像素点m，n的色度距离，当小于阈值τ时，ω_m,n的值为100，反之则取0；在本发明中阈值τ的优选值为0.0009；符号||·||₂表示2范数，<m,n>为像素点m和n的内积。

(2)一致性约束

一致性约束旨在约束双目图像中左右视角重叠区域(Overlap Region)匹配点的反射分量应一致。在本发明中，左右视角的匹配通过通用的光流法(Optical Flow)来计算，依此得到左视角相对于右视角的运动矢量(Motion Vector)。一致性约束的含义是左右视角同一场景点的反射分量应该严格一致，不随外界光照和观测视角的变化而不同。由光流法计算获得运动矢量后，可以对左右视角匹配点的反射分量增加如下约束条件：

式中，m₁与m′代表左右视角重叠区域匹配的像素点对；

和分别表示左右视角匹配像素点对所对应的像素值；

和分别表示左右视角匹配像素点对所对应的明暗分量值；

(3)局部约束

一般而言，在一副图像中，局部的像素块具有近似的反射属性。基于图像的局部特点，本发明加入局部约束条件，更进一步地降低图像本征性质分解的不确定度。局部约束是通过使得局部具有相似反射属性的像素块的反射分量尽量保持一致。在本发明中，采用超像素技术(Super-Pixel)来对具有近似反射属性的像素块进行划分。使用超像素技术对输入图像进行局部划分的结果如图3所示。图3中子图(a)是原图，(b)是子图(a)中的一个方块区域经过Superpixel技术分割后处于同一个Superpixel块的放大图，子图(c)是子图(a)中的另一个方块区域经过Superpixel技术分割后处于同一个Superpixel块的放大图，子图(d)是子图(a)经过Superpixel技术所得的分割后的Superpixel图。在子图(d)中，具有相同颜色的块即认为具有近似的反射属性，通过约束该块中的像素点的反射分量的值一致，可以进一步降低图像本征性质分解的不确定性。局部约束条件如下所示：

式中，表示任意像素块；

Γ_l表示输入图像经过超像素划分后所得到的像素块集合；

m₂和n₂表示某一像素块中任意的两个像素点；

和分别表示像素点m₂和像素点n₂的像素值；

和的分别表示像素点m₂和像素点n₂的明暗分量值。

(4)非局部约束

与局部约束类似，反射属性的相似性不仅具有局部性，很有可能在图像的全局中存在这种相似性，这是非局部约束最直接的出发点；非局部约束的示意图如图4所示。图4中，子图(a)是原图，子图(b)是子图(a)中的一个方块区域经过Superpixel技术分割后处于不同一个Superpixel块的放大图，子图(c)是子图(a)中的另一个方块区域经过Superpixel技术分割后处于不同一个Superpixel块的放大图，子图(d)是子图(a)经过Superpixel技术所得的分割后的Superpixel图。

所述非局部约束E_nl(S)为：

式中，Γ_nl表示左右两个视角中的像素块集合；

和是Γ_nl集合中任意的两个超像素块；

m₃和n₃表示两个不同超像素块中任意像素点所组成的左右视角匹配像素点对；

和分别表示像素点m₃和像素点n₃的像素值；

和分别表示像素点m₃和像素点n₃的明暗分量值。

与局部约束中的像素点对不同的是，局部约束中的像素点对属于同一个超像素块，而非局部约束中的像素点对分别属于两个不同的超像素块。

(5)绝对尺度约束

由于图像的本征性质分解问题的欠定性，导致最后分解所得的明暗分量或反射分量可能超出合理的取值范围。此时需要对明暗分量或反射分量的绝对尺度做限定，使得最后所得的分解更复合实际情况。

绝对尺度约束的出发点在于：图像中像素值越大的地方往往意味着明暗分量的值越高。因此在本发明中，通过对图像中最亮的像素点集合加上绝对尺度约束，使得最终分解后所得的明暗分量归一化到单位值。

所述绝对尺度约束E_a(S)为：

式中，P_a表示图像中最亮像素点集合；

m₄表示集合P_a中的任意一个最亮像素点；

表示像素点m₄的明暗分量值。

通过该约束，尽量使得图像中最亮的像素点的明暗分量为单位值1。

上述步骤S3中的目标函数为：

式中：λ_r、λ_d、λ_l、λ_nl和λ_a分别为Retinex约束、一致性约束、局部约束、非局部约束和绝对尺度约束的权重。

通过调整各种不同约束条件的权重，可以满足不同约束条件对最终分解的结果不同影响的效果。上述目标函数中的5个约束代价函数均为二次型的，经过线性加权后，所得的总的目标代价函数同样也是二次型的，因此总的代价函数可以得到一个具有闭合形式的最优解。因此目标函数的最优解可以通过求解一个大型的稀疏线性方程组来优化最后总的代价函数，从而完成双目图像的本征性质分解。

由总的代价函数(10)可知，待优化的未知变量为S，即双目图像的明暗分量为待求的未知量。前面提及的5种不同约束项分别对明暗分量和反射分量做出了约束。尽管总的代价函数中存在反射分量的约束，但由于明暗分量和反射分量还在公式(2)、(3)的约束关系，因此关于反射分量的约束可以转化成与之对应的明暗分量的约束。故我们可以看到最后总的代价函数(10)并没有直接出现关于反射分量的约束，而是将反射分量(Reflectance)转化成I/S。这一转化是因明暗分量和反射分量还存在约束关系(2)和(3)而得出的。总的代价函数E(S)中并没有出现I/S，而是以I_m-S_m或I_n-S_n的形式体现，原因在于我们将明暗分量和反射分量的约束关系转化到Log域中(Log域中的“+”和“-”分别对应于乘法和除法运算)。

因此，上述步骤S4中将确定目标函数的最优解的过程转化对大型稀疏线性方程组求解的过程；

其中，大型稀疏线性方程组为Ax＝b；

式中，A为大型稀疏线性方程组的系数矩阵；

x为待求向量，即双目图像的明暗分量；

b为大型系数线性矩阵的列向量。

对大型稀疏线性方程组求解的方法具体为：

A1、依次将双目图像本征性质分解时的约束项进行转换，得到对应最优解条件下的约束关系；

A2、将转换后的约束关系依次加入到大型稀疏线性方程组中，构建满足约束关系的大型稀疏线性方程组；

A3、通过右除操作对构建完成的大型稀疏线性方程组进行处理，得到双目图像的明暗分量，即目标函数的最优解。

在上述步骤A1中：

(1)Retinex约束E_r(S)对应的约束关系为：

得到该约束关系的方法为：

在Retinex约束项的代价函数公式(4)中存在关于明暗分量和反射分量转化后的约束。公式中的ω_m,n是一个门限函数，函数的具体形式如公式(5)所示。它所起的作用是：当相邻的两个像素点m，n的Chromaticity的差别小于某一个阈值时，ω_m,n＝100，否则ω_m,n＝0。而Chromaticity的差别则由来衡量。其中，这里的“<m,n>”表示两相邻像素点m，n三个颜色通道的内积。“||m||₂||n||₂”则表示两相邻像素点m，n的像素值的2范数之积。ω_m,n在这里体现了Retinex的作用：输入图像中大的变化认为是由反射分量的变化引起的，而微小的变化则是由明暗分量的变化引起的。通过衡量两相邻像素点之间的差别，来判断输入图像局部的变化，进而赋予ω_m,n的为100或0，从而使得最终将这种变化归结到明暗分量或者反射分量上。例如，如果两相邻像素点的Chromaticiy差别小于某一个阈值τ时，ω_m,n＝100，按照Retinex的假设，此时图像相邻像素点的值变化是由明暗分量所引起的，反射分量前的权重系数ω_m,n非常大，因此微小的反射分量变化将引起Retinex约束代价函数的剧烈变化，在优化过程中，相邻像素点的像素值的微小变化将归结到明暗分量中。反之，则将相邻像素点的像素值的较大变化归结到反射分量的变化中。

利用多元函数的最值理论，以明暗分量为未知量，对Retinex约束项的代价函数求导可得：

在上式中，ΔS_m,n＝S_m-S_n，ΔI_m,n＝I_m-I_n，其分别代表两相邻像素点的明暗分量的差值和输入图像像素点差值。

利用最值理论的结论，要使该代价函数最小化，只需令上式等于0即可。

故对于任意像素点m，其与邻域内的像素点n的明暗分量存在如下关系：

(2)一致性约束E_d(S)对应的约束关系为：

得到该约束关系的方法为：

在本发明中，我们采用光流法来计算左右视角的匹配，依此可以得到左视角相对于右视角的运动矢量；

对一致性约束项代价函数公式(6)而言，根据最优化理论，同样以明暗分量为未知变量，对一致性约束项代价函数求导可得：

同样利用最值理论，只要上式导数为0，该一致性约束项的代价函数将取得最小值，由此对于左视角的任意像素点m₁，其与右视角匹配的像素点m'存在如下关系：

(3)局部约束E_l(S)对应的约束关系为：

得到该约束关系的方法具体为：

在一副图像中，局部的像素块具有近似的反射属性。基于图像的局部特点，我们考虑加入局部约束条件，这可以更进一步地降低图像本征性质分解的不确定性。在本发明中，我们采用经典的超像素技术首先对左右视角的输入图像做粗略的划分，这样一来可以得到一张相似反射属性的像素块聚成一个超像素的图像，如图3所示。然后加入如公式(7)所示的局部约束，以此来增加反射分量的约束。

同样采用最值理论，以明暗分量为未知量，对局部约束项的代价函数(7)求导可得：

令上式等于0，对任一个超像素块内任意两个像素点m₂，n₂存在如下关系：

(4)非局部约束E_nl(S)对应的约束关系为：

得到该约束方关系的方法为：

反射属性的相似性不仅具有局部性，在图像的全局中同样可能存在这种相似性。这是非局部约束的基本出发点。非局部约束的示意图如附图4所示。首先也是先通过超像素块技术先把左右视角的输入图像做一个粗略的划分，这样具有相似的反射属性的局部像素块被划分成一个超像素块，在局部约束项中我们已经加入了局部约束项，使得局部像素块中的像素点的反射分量具有接近的值。在此基础上，我们再对划分后的超像素块进行近似最近邻处理，这样一来我们可以获得与该超像素块具有接近反射属性的全局超像素块。然后在这两个非局部的超像素块中加入非局部约束项如代价函数(8)所示。如此可进一步降低非局部的反射分量的不确定性。

同样采用最值理论，以明暗分量为未知量，对非局部约束项代价函数(8)求导，可得到如下的关系：

在上式中，分别表示不同超像素块内任意两个像素点m₃，n₃的明暗分量，分别表示不同超像素块中任意两个像素点m₃，n₃的像素点值。令上式等于0，我们可以得到非局部项代价函数取得最小能量的条件：对于不同超像素块中的任意两个像素点m₃和n₃的明暗分量存在如下约束：

(5)绝对尺度约束E_a(S)对应的约束关系为：

得到该约束关系的方法为：

绝对尺度约束的出发点在于：图像中像素值越大的地方一般明暗分量的值越高。在本发明中，通过对图像中最亮的像素点加上绝对尺度约束，使得最终分解后所得的明暗分量归一化到单位值。因此我们在总的代价函数中加入绝对尺度约束项代价函数(9)。显然，绝对尺度约束项的目的在于对于图像中最亮的像素点的明暗分量的值应为单位值1。因此对最亮的像素点的明暗分量有如下约束：

上述步骤A2中，将转换后的约束关系依次加入到大型稀疏线性方程组中时，

对于Retinex约束E_r(S)对应的约束关系加入大型稀疏线性方程组的方法为：

令双目图像中像素点m的坐标为(x,y)，在左视图中该像素点m列向量化后的明暗分量S_(x,y)为x的第y×M+x个分量，在右视图中该像素点m列向量化后的明暗分量S_(x,y)为x的第M×N+y×M+x个分量；

式中，M和N分别为双目图像的高度和宽度；

对于左视图，若x_y×M+x存在n个邻域像素点，则在系数矩阵A的第y×M+x行、第y×M+x列的元素的基础上加上n×λ_r，在系数矩阵A第y×M+x行，第y×M+x-1列、第y×M+x+1列、第(y-1)×M+x列和第_(y+1)×M+x列的元素的基础上分别加上-λ_r；

在列向量b的第y×M+x行的元素的基础上加上

对于右视图，若x_y×M+x存在n个邻域像素点，则在系数矩阵A的第M×N+y×M+x行、第M×N+y×M+x列的元素的基础上加上n×λ_r，在系数矩阵A第M×N+y×M+x行，第M×N+y×M+x-1列、第M×N+y×M+x+1列、第M×N+(y-1)×M+x列和第M×N+(y+1)×M+x列的元素的基础上分别加上-λ_r；

在列向量b的第M×N+y×M+x行的元素的基础上加上

式中，ω_{(x,y),(x-1,y)}表示坐标为(x,y)的像素点与其邻域的坐标为(x-1,y)的像素点的色度距离；

i_leftlog(x,y)和i_leftlog(x-1,y)分别表示左视图中坐标为(x,y)的像素点的log域值和坐标为(x-1,y)的像素点的log域值；

根据所述Retinex约束E_r(S)对应的约束关系的加入方法，依次对左右视图的每一个像素点往系数矩阵A和列向量b中添加Retinex约束E_r(S)；

对于一致性约束E_d(S)对应的约束关系加入大型稀疏线性方程组的方法为：

令左视图中像素点m₁的坐标为(x1,y1)，其经过光流放计算得到的相对于右视图的运动矢量为v(Δx1,Δy1)，与像素点m₁所匹配的右视图像素点m′的坐标为(x1+Δx1,y1+Δy1)；

对于左视图，在系数矩阵A的第行、第列的元素的基础上加上权重λ_d，在系数矩阵A的第行、第y1×M+x1列的元素的基础上加上-λ_d；

在列向量b的第行的元素的基础上加上

对于右视图，在系数矩阵A的第行、第列的元素的基础上加上权重λ_d，在系数矩阵A的第行、第M1×N+y1×M+x1列的元素的基础上加上-λ_d；

在列向量b的第行的元素的基础上加上

根据所述一致性约束E_d(S)对应的约束关系加入方法依次将左视图和右视图中所有匹配点对的约束关系加入到系数矩阵A和列向量b中；

对于局部约束E_l(S)对应的约束关系加入大型稀疏线性方程组的方法为：

令双目图像中同一像素块的像素点m₂和像素点n₂的坐标分别为(x2,y2)和(u,v)；

对于左视图，在系数矩阵A的第v×M+u行、第v×M+u列的元素的基础上加上权重λ_l，在系数矩阵A的第v×M+u行、第v×M+u列的元素的基础上加上-λ_l；

在列向量b的第y×M+x行的元素的基础上加上λ_l×(i_leftlog(x2,y2)-i_leftlog(u,v))；

对于右视图，在系数矩阵A的第M×N+y2×M+x2行、第M×N+y2×M+x2列的元素的基础上加上权重λ_l，在系数矩阵A的第M×N+y2×M+x2行，第M×N+v×M+u列的元素的基础上加上-λ_l；

在列向量b的第M×N+y2×M+x2行的元素的基础上加上λ_l×(i_leftlog(x2,y2)-i_leftlog(u,v))；

根据所述局部约束E_l(S)对应的约束关系加入方法依次将左视图和右视图中所有超像素块的约束关系加入到系数矩阵A和列向量b中；

对于非局部约束E_nl(S)对应的约束关系加入大型稀疏线性方程组的方法为：

令双目图像中非同一像素块的像素点m₃和像素点n₃的坐标分别为(x3,y3)和(u1,v1)；

对于左视图，在系数矩阵A的第v1×M+u1行、第v1×M+u1列的元素的基础上加上权重λ_l，在系数矩阵A的第v1×M+u1行、第v1×M+u1列的元素的基础上加上-λ_l；

在列向量b的第y3×M+x3行的元素的基础上加上λ_l×(i_leftlog(x3,y3)-i_leftlog(u1,v1))；

对于右视图，在系数矩阵A的第M×N+y3×M+x3行、第M×N+y3×M+x3列的元素的基础上加上权重λ_l，在系数矩阵A的第M×N+y3×M+x3行，第M×N+v1×M+u1列的元素的基础上加上-λ_l；

在列向量b的第M×N+y3×M+x3行的元素的基础上加上λ_l×(i_leftlog(x3,y3)-i_leftlog(u1,v1))；

根据所述非局部约束E_nl(S)对应的约束关系加入方法依次将左视图和右视图中所有超像素块的约束关系加入到系数矩阵A和列向量b中；

对于绝对尺度约束E_a(S)对应的约束关系加入大型稀疏线性方程组的规则为：

令双目图像中的最亮的像素点m₄的坐标为(x4,y4)；

对于左视图，在系数矩阵A的第y4×M+x4行、第y4×M+x4列的元素的基础上加上权重λ_a；

在列向量b的第y4×M+x4行的元素的基础上加上0；

对于右视图，在系数矩阵A的第M×N+y4×M+x4行、第M×N+y4×M+x4列的元素的基础上加上权重λ_a；

在列向量b的第M×N+y4×M+x4行的元素的基础上加上0。

上述步骤A3中，构建完成的大型系数线性方程组为一个2×M×N行，2×M×N列的对称方程，通过调用Matlab中的右除操作，使x＝A\b，进而得到双目图像的明暗分量x。

上述步骤S5具体为：

S51、从获得的明暗分量x中，取前M×N个分量，并按照列主序的方式将该M×N维的列向量转成一个M行，N列的左视图的明暗分量I_leftshading；

取明暗分量x的后M×N个分量，按照列主序的方式将该M×N维的列向量转成一个M行，N列的右视图的明暗分量I_rightshading；

S52、根据左右视图的明暗分量，通过图像成像约束，分别得到左视图的反射分量I_{leftreflectance}和右视图的反射分量I_{rightreflectance}；

S53、依次对左视图的明暗分量I_leftshading、右视图的明暗分量I_rightshading、左视图的反射分量I_{leftreflectance}和右视图的反射分量I_{rightreflectance}进行反归一化操作的，得到双目图像左右视图的明暗分量和反射分量，实现双目图像的本征性质分解。

在本发明的一个实施例中，提供了将上述五种不同的约束关系转换成构建一个大型稀疏线性方程组，通过解一个大型稀疏线性方程组来获得双目图像明暗分量和反射分量的具体过程：

第1步：对左右视角的双目图像I_leftinput和I_rightinput，对其分别进行归一化处理，分别得到归一化的左右视角图像I_{leftnormalized}和I_{rightnormalized}，

第2步：利用Chromaticity图像的计算公式：分别得到左右视角的Chromaticity图I_leftchroma和I_leftchroma。

在公式中，I_normalized、||I_normalized||₂和I_chroma分别表示归一化后的输入图像、归一化后的输入图像的二范数灰度图和归一化后的输入图像的Chromaticity图，这里的除法操作是逐像素点的除法操作；

即

第3步：将归一化后的左右视角图像I_{leftnormalized}和I_{rightnormalized}转化为与之相对应的灰度图I_leftgray和I_rightgray；

RGB三通道的彩色图转成灰度图的计算公式为：I_gray＝0.299*R+0.587*G+0.114*B。

根据这个计算公式，可以得到左右视角的灰度图I_leftgray和I_rightgray。

第4步：在左右视角的灰度图的基础上，再计算其Log域的值，分别得到左右视角Log域的图像i_leftlog和i_rightlog，分别为：i_leftlog＝log(I_leftgray)、i_rightlog＝log(I_rightgray)。

第5步：初始化一个2×M×N行，2×M×N列的大型稀疏矩阵A_{2×M×N,2×M×N},一个2×M×N行，1列的待求向量x_2×M×N，以及一个2×M×N行，1列的向量b_2×M×N。它们分别表示大型稀疏线性方程组的系数矩阵，待求的双目图像的左右明暗分量列向量化后所组成的列向量以及大型稀疏线性方程组的增广向量。

需要说明的是：这里假设输入的双目图像的高和宽分别为M和N，由于明暗分量为单色图，因此单一视角的明暗分量列向量化后为一个M×N行，1列的列向量，待求的向量x_2×M×N为左右视角明暗分量列向量化后所组成的一个2×M×N行，1列的列向量。

第6步：将权重系数λ_r、λ_d、λ_l、λ_nl、λ_a分别设为：1.0，1.0，1.0，1.0，1000.0，将Chromaticity差别的阈值τ设为0.0003。

第7步：根据Retinex约束项代价函数要取得最优解的结论，我们将这一约束关系填入到大型稀疏线性方程组Ax＝b中，这里主要的工作是构建矩阵A和列向量b。

在这一大型稀疏线性方程组中，将左右视角的明暗分量S作为待求的未知变量x。下面将详细讲解如何将这一约束关系转化到大型稀疏线性方程组Ax＝b中。

以左视角图像中的任意一个像素点m为例，假设像素点m的坐标为(x,y)，则它的4邻域像素点的坐标为：(x-1,y)，(x+1,y)，(x,y-1)，(x,y+1)。若像素点m为边界点，则它的4领域像素点只存在2个或3个，具体根据邻域像素点的坐标是否超越输入图像的边界为准。

取像素点m的其中一个邻域像素(x-1,y)为例，首先计算像素点(x,y)与其邻域像素点(x-1,y)的Chromaticity的差值。

点(x,y)与点(x-1,y)的Chromaticity差值为：

如果Δ的值小于设定的阈值τ，此时ω_{(x,y),(x-1,y)}＝100，反之，则ω_{(x,y),(x-1,y)}＝0。结合前面计算所得的图像Log域的值，可以将具体化为：

然后根据上述这一结论开始给系数矩阵A和列向量b在恰当的位置填入恰当的值。

具体方式如下：

像素点m的位置坐标为(x,y)，其任意一个像素点坐标为(x-1,y)，那么列向量化以后的明暗分量S_(x,y)为x的第y×M+x个分量，即S_(x,y)对应于x_y×M+x，同理明暗分量S_(x-1,y)对应于x_y×M+x-1。故将变换成如下形式：

类似的，可以得到像素点(x,y)与其它几个邻域像素点相似的约束关系如下所示：

将它们同时相加可得：

合并后得到一个方程，在此基础上，还需要对合并后的方程乘上加权系数，得到如下形式：

然后按照合并后的方程的系数根据如下规则填入到线性方程组Ax＝b中。

如果x_y×M+x存在n个邻域像素点，则在系数矩阵A的第_y×M+x行，第_y×M+x列上加上n×λ_r，同时，在系数矩阵A的第_y×M+x行，第_y×M+x-1、_y×M+x+1、_(y-1)×M+x、_(y+1)×M+x列中分别加上-λ_r，在列向量b的第_y×M+x行上加上值：

如果像素点(x,y)在图像的边界上，则上述方程组可能只有2个或3个，具体情况根据其邻域的像素点坐标是否超越了图像的边界而定。一般情况下，任意像素点均存在上述4个线性方程。因此在系数矩阵A和列向量b赋值时应根据像素点(x,y)的领域像素点个数情况考虑。

前面讨论的情况适用于左视图的像素点的Retinex约束。

对于右视图的Retinex约束而言，情况非常类似，唯一的区别在于：对于右视图的同一个位置坐标(x,y)像素点而言，其明暗分量对应于x_M×N+y×M+x,而非x_y×M+x，同理右视图像素点(x,y)的任意邻域(x-1,y)的明暗分量对应于x_{M×N+y×M+x-1}。可以发现右视图上的像素点的明暗分量与左视图相同位置的明暗分量在待求变量x的下标正好相差M×N。原因在于左视图的明暗分量正好是待求变量x的前M×N个分量。

按照上述方法，依次对左右视图的每一个像素点往系数矩阵A和列向量b添加Retinex约束；到此完成了Retinex约束项的转化。

第8步：根据反射分量一致性约束项所得出的结论，要往大型稀疏线性方程组Ax＝b添加的约束关系为

根据已有的光流算法，首先计算出左右视图的运动矢量，依此可以得到左右视图的匹配像素点对。

以左视图的某一个像素点m₁为例，假设像素点m₁的坐标为(x1,y1)，经过光流法计算得到的相对于右视图的运动矢量为v(Δx1,Δy1)，因此与像素点m所匹配的右视图像素点m′1的坐标为(x1+Δx1,y1+Δy1)。因此约束条件可以变换成

注意上式中的符号表示向上取整运算，例如

然后将该方程按照添加Retinex约束类似的方式添加到方程组Ax＝b中去。

例如，对于系数矩阵A而言，在A的第_y1×M+x1行，第_y1×M+x1列元素中加上λ_d，在A的第_y1×M+x1行，第列元素中加上-λ_d，且在列向量b的第y1×M+x1行元素中加上

在我们的实现中，为了使得分解算法的速度更快，往往将矩阵A兑成一个对称正定的方阵，因此将变换成如下形式：类似的，在系数矩阵A的第行，第列元素中加上λ_d，在系数矩阵A的第行，第y1×M+x1列元素中加上-λ_d，在列向量b的第行元素中加上

按照这种方式，依次对左右视图所有匹配点对的约束关系加入到系数矩阵A和列向量b的元素中。

第9步：根据局部约束项所得的结论，需往大型稀疏线性方程组Ax＝b添加一系列约束这里需要先利用超像素技术对左右视图进行超像素块划分，得到如附图3所示的超像素块图。在此基础上，依次对左右视图所有超像素块中的任意两个像素点添加约束：

以左视图中的其中一个超像素块中的两个像素点m₂和n₂为例，假设像素点m₂的坐标为(x2,y2)，像素点n₂的坐标为(u,v)。因此约束方程可以变换成如下形式：

λ_l×(x_y2×M+x2-x_v×M+u)＝λ_l×(i_leftlog(x2,y2)-i_leftlog(u,v)) (28)

仿照前面添加约束方程方法，在系数矩阵A的第y2×M+x2行，第y2×M+x2列元素中加上λ_l，在系数矩阵A的第y2×M+x2行，第v×M+u列元素中加上-λ_l，在列向量b的第y2×M+x2行元素中加上λ_l×(i_leftlog(x2,y2)-i_leftlog(u,v))。

出于同样的考虑，在系数矩阵A的第v×M+u行，第v×M+u列元素中加上λ_l，在系数矩阵A的第v×M+u行，第y2×M+x2列元素中加上-λ_l，在列向量b的第v×M+u行元素中加上λ_l×(i_leftlog(u,v)-i_leftlog(x2,y2))。

按照上述方法，依次对左右视图的所有超像素块中的局部约束方程添加到系数矩阵A和列向量b的元素中。

第10步：与局部约束项相似，非局部约束项的约束方程按照相似的方式添加到系数矩阵A和列向量b中。

非局部约束方程为：该约束方程与局部约束方程的区别在于：像素点m₃和n₃不在同一个超像素块中，而是处于两个近似的超像素块中。因此在添加该约束方程时需要注意像素点m₃和n₃的坐标，但是在往系数矩阵A和列向量b的元素中添加值的方式是类似的。

按照上述方法，依次对左右视图的非局部约束方程添加到系数矩阵A和列向量b的元素中。

第11步：最后添加绝对尺度约束方程。该约束方程为：由于将图像转换到Log域中，因此约束方程变成：再考虑绝对尺度权重后，约束方程变换成：在添加该约束方程时，需要首先通过计算寻找到左右视图中的最亮像素点m的坐标。

寻找最亮像素点的方法：依次遍历左右视图i_leftgray和i_rightgray，分别找到i_leftgray和i_rightgray图像中值最大的像素点m₄和m₄'，假设他们的坐标分别为(x4,y4)和(x4',y4')。

对于左视图最亮的像素点m而言，需要添加的约束方程为：变形后得到如下形式：

λ_a×x_y4×M+x4＝0 (29)

因此，在系数矩阵A的第y4×M+x4行，第y4×M+x4列元素中加上λ_a，在列向量b的第y4×M+x4行元素中加上0。

第12步：完成第7～11步后，所有的约束方程均加入到大型线性方程组Ax＝b后，此时的方程组Ax＝b成为一个2×M×N行，2×M×N列的对称方程。该方程可以通过调用Matlab的右除操作“\”，即x＝A\b。此后即可得到左右视角的明暗分量x。

第13步：从求解出来的变量x中，取出前M×N个分量，然后按照列主序的方式将该M×N维的列向量转成一个M行，N列的明暗图像I_leftshading，取变量x后的后M×N个分量，同样按照列主序的方式将该M×N维的列向量转成一个M行，N列的明暗图像I_rightshading。

第14步：根据图像成像约束：I＝R×S，分别得到左右视图的反射分量：I_{leftreflectance}＝I_{leftnormalized}/I_leftshading，I_{rightreflectance}＝I_{rightnormalized}/I_rightshading。

第15步：最后分别对左右视图的明暗分量I_leftshading和I_rightshading，反射分量、I_{leftreflectance}和I_{rightreflectance}进行反归一化操作，得到双目图像的明暗分量和反射分量，实现双目图像本征性质的分解。

图5给出了通过上述本征性质分解方法得到5组示例图，中5从第1列到第6列依次为左视图、右视图、左明暗分量图、右明暗分量图、左反射分量图、右反射分量图。

本发明的有益效果为：本发明提供的双目图像的本征性质分解方法以双目图像中的明暗分量为待求的未知量，结合本发明提供的Retinex约束、反射分量一致性约束、局部约束、非局部约束以及绝对尺度约束等构造一个二次型的目标函数，从而双目图像的本征性质分解问题转化成为一个求解二次型代价函数能量最小解的问题。经过理论推导和实例验证，该技术方案具有分解效果好，算法效率高，容易实现等优点。

Claims (8)

1.一种双目图像的本征性质分解方法，其特征在于，包括以下步骤：

S1、通过双目相机获得待进行本征性质分解的双目图像；

S2、确定双目图像本征性质分解时的约束项；

S3、根据约束项确定双目图像本征性质分解的目标函数；

S4、确定目标函数的最优解，进而获得双目图像的明暗分量；

S5、根据获得的双目图像的明暗分量，获得双目图像左右视图的明暗分量和反射分量，实现双目图像的本征性质分解；

所述步骤S2中双目图像本征性质分解时的约束项包括Retinex约束、一致性约束、局部约束、非局部约束和绝对尺度约束；

所述Retinex约束E_r(S)为：

式中，m，n分别表示双目图像中的两个任意像素点；

表示像素点m的一个预设领域；

ω_m,n表示相邻的像素点m，n的色度距离；

ΔS_m,n表示明暗分量差量，且ΔS_m,n＝S_m-S_n，S_m和S_n分别表示双目图像中像素点m和像素点n的明暗分量值；

ΔI_m,n表示像素差量，且ΔI_m,n＝I_m-I_n，I_m和I_n分别表示双目图像中像素点m和像素点n的像素值；

所述一致性约束E_d(S)为：

式中，m₁与m′代表左右视角重叠区域匹配的像素点对；

和分别表示左右视角匹配像素点对所对应的像素值；

和分别表示左右视角匹配像素点对所对应的明暗分量值；

所述局部约束E_l(S)为：

式中，P_l ⁱ表示任意像素块；

Γ_l表示输入图像经过超像素划分后所得到的像素块集合；

m₂和n₂表示某一像素块中任意的两个像素点；

和分别表示像素点m₂和像素点n₂的像素值；

和的分别表示像素点m₂和像素点n₂的明暗分量值；

所述非局部约束E_nl(S)为：

式中，Γ_nl表示左右两个视角中的像素块集合；

和是Γ_nl集合中任意的两个超像素块；

m₃和n₃表示两个不同超像素块中任意像素点所组成的左右视角匹配像素点对；

和分别表示像素点m₃和像素点n₃的像素值；

和分别表示像素点m₃和像素点n₃的明暗分量值；

所述绝对尺度约束E_a(S)为：

式中，P_a表示图像中最亮像素点集合；

m₄表示集合P_a中的任意一个最亮像素点；

表示像素点m₄的明暗分量值。

2.根据权利要求1所述的双目图像本征性质分解方法，其特征在于，所述步骤S3中的目标函数为：

式中：λ_r、λ_d、λ_l、λ_nl和λ_a分别为Retinex约束、一致性约束、局部约束、非局部约束和绝对尺度约束的权重。

3.根据权利要求2所述的双目图像的本征性质分解方法，其特征在于，所述步骤S4中，将确定目标函数的最优解的过程转化对大型稀疏线性方程组求解的过程；

其中，大型稀疏线性方程组为Ax＝b；

式中，A为大型稀疏线性方程组的系数矩阵；

x为待求向量，即双目图像的明暗分量；

b为大型系数线性矩阵的列向量。

4.根据权利要求3所述的双目图像的本征性质分解方法，其特征在于，对大型稀疏线性方程组求解的方法具体为：

A1、依次将双目图像本征性质分解时的约束项进行转换，得到对应最优解条件下的约束关系；

A2、将转换后的约束关系依次加入到大型稀疏线性方程组中，构建满足约束关系的大型稀疏线性方程组；

A3、通过右除操作对构建完成的大型稀疏线性方程组进行处理，得到双目图像的明暗分量，即目标函数的最优解。

5.根据权利要求4所述的双目图像的本征性质分解方法，其特征在于，所述步骤A1中，所述Retinex约束E_r(S)对应的约束关系为：

所述一致性约束E_d(S)对应的约束关系为：

所述局部约束E_l(S)对应的约束关系为：

所述非局部约束E_nl(S)对应的约束关系为：

所述绝对尺度约束E_a(S)对应的约束关系为：

6.根据权利要求5所述的双目图像的本征性质分解方法，其特征在于，所述步骤A2中，将转换后的约束关系依次加入到大型稀疏线性方程组中时，

对于Retinex约束E_r(S)对应的约束关系加入大型稀疏线性方程组的方法为：

令双目图像中像素点m的坐标为(x,y)，在左视图中该像素点m列向量化后的明暗分量S_(x,y)为x的第y×M+x个分量，在右视图中该像素点m列向量化后的明暗分量S_(x,y)为x的第M×N+y×M+x个分量；

式中，M和N分别为双目图像的高度和宽度；

对于左视图，若x_y×M+x存在n个邻域像素点，则在系数矩阵A的第y×M+x行、第y×M+x列的元素的基础上加上n×λ_r，在系数矩阵A第y×M+x行，第y×M+x-1列、第y×M+x+1列、第(y-1)×M+x列和第(y+1)×M+x列的元素的基础上分别加上-λ_r；

在列向量b的第y×M+x行的元素的基础上加上

对于右视图，若x_y×M+x存在n个邻域像素点，则在系数矩阵A的第M×N+y×M+x行、第M×N+y×M+x列的元素的基础上加上n×λ_r，在系数矩阵A第M×N+y×M+x行，第M×N+y×M+x-1列、第M×N+y×M+x+1列、第M×N+(y-1)×M+x列和第M×N+(y+1)×M+x列的元素的基础上分别加上-λ_r；

在列向量b的第M×N+y×M+x行的元素的基础上加上

式中，ω_{(x,y),(x-1,y)}表示坐标为(x,y)的像素点与其邻域的坐标为(x-1,y)的像素点的色度距离；

i_leftlog(x,y)和i_leftlog(x-1,y)分别表示左视图中坐标为(x,y)的像素点的log域值和坐标为(x-1,y)的像素点的log域值；

根据所述Retinex约束E_r(S)对应的约束关系的加入方法，依次对左右视图的每一个像素点往系数矩阵A和列向量b中添加Retinex约束E_r(S)；

对于一致性约束E_d(S)对应的约束关系加入大型稀疏线性方程组的方法为：

令左视图中像素点m₁的坐标为(x1,y1)，其经过光流放计算得到的相对于右视图的运动矢量为v(Δx1,Δy1)，与像素点m₁所匹配的右视图像素点m′的坐标为(x1+Δx1,y1+Δy1)；

对于左视图，在系数矩阵A的第行、第列的元素的基础上加上权重λ_d，在系数矩阵A的第行、第y1×M+x1列的元素的基础上加上-λ_d；

在列向量b的第行的元素的基础上加上

对于右视图，在系数矩阵A的第行、第列的元素的基础上加上权重λ_d，在系数矩阵A的第行、第M1×N+y1×M+x1列的元素的基础上加上-λ_d；

在列向量b的第行的元素的基础上加上

根据所述一致性约束E_d(S)对应的约束关系加入方法依次将左视图和右视图中所有匹配点对的约束关系加入到系数矩阵A和列向量b中；

对于局部约束E_l(S)对应的约束关系加入大型稀疏线性方程组的方法为：

令双目图像中同一像素块的像素点m₂和像素点n₂的坐标分别为(x2,y2)和(u,v)；

对于左视图，在系数矩阵A的第v×M+u行、第v×M+u列的元素的基础上加上权重λ_l，在系数矩阵A的第v×M+u行、第v×M+u列的元素的基础上加上-λ_l；

在列向量b的第y×M+x行的元素的基础上加上λ_l×(i_leftlog(x2,y2)-i_leftlog(u,v))；

对于右视图，在系数矩阵A的第M×N+y2×M+x2行、第M×N+y2×M+x2列的元素的基础上加上权重λ_l，在系数矩阵A的第M×N+y2×M+x2行，第M×N+v×M+u列的元素的基础上加上-λ_l；

在列向量b的第M×N+y2×M+x2行的元素的基础上加上λ_l×(i_leftlog(x2,y2)-i_leftlog(u,v))；

根据所述局部约束E_l(S)对应的约束关系加入方法依次将左视图和右视图中所有超像素块的约束关系加入到系数矩阵A和列向量b中；

对于非局部约束E_nl(S)对应的约束关系加入大型稀疏线性方程组的方法为：

令双目图像中非同一像素块的像素点m₃和像素点n₃的坐标分别为(x3,y3)和(u1,v1)；

对于左视图，在系数矩阵A的第v1×M+u1行、第v1×M+u1列的元素的基础上加上权重λ_l，在系数矩阵A的第v1×M+u1行、第v1×M+u1列的元素的基础上加上-λ_l；

在列向量b的第y3×M+x3行的元素的基础上加上λ_l×(i_leftlog(x3,y3)-i_leftlog(u1,v1))；

对于右视图，在系数矩阵A的第M×N+y3×M+x3行、第M×N+y3×M+x3列的元素的基础上加上权重λ_l，在系数矩阵A的第M×N+y3×M+x3行，第M×N+v1×M+u1列的元素的基础上加上-λ_l；

在列向量b的第M×N+y3×M+x3行的元素的基础上加上λ_l×(i_leftlog(x3,y3)-i_leftlog(u1,v1))；

根据所述非局部约束E_nl(S)对应的约束关系加入方法依次将左视图和右视图中所有超像素块的约束关系加入到系数矩阵A和列向量b中；

对于绝对尺度约束E_a(S)对应的约束关系加入大型稀疏线性方程组的规则为：

令双目图像中的最亮的像素点m₄的坐标为(x4,y4)；

对于左视图，在系数矩阵A的第y4×M+x4行、第y4×M+x4列的元素的基础上加上权重λ_a；

在列向量b的第y4×M+x4行的元素的基础上加上0；

对于右视图，在系数矩阵A的第M×N+y4×M+x4行、第M×N+y4×M+x4列的元素的基础上加上权重λ_a；

在列向量b的第M×N+y4×M+x4行的元素的基础上加上0。

7.根据权利要求6所述的双目图像的本征性质分解方法，其特征在于，所述步骤A3中，构建完成的大型系数线性方程组为一个2×M×N行，2×M×N列的对称方程，通过调用Matlab中的右除操作，使x＝A\b，进而得到双目图像的明暗分量x。

8.根据权利要求7所述的双目图像的本征性质分解方法，其特征在于，所述步骤S5具体为：

S51、从获得的明暗分量x中，取前M×N个分量，并按照列主序的方式将该M×N维的列向量转成一个M行，N列的左视图的明暗分量I_leftshading；

取明暗分量x的后M×N个分量，按照列主序的方式将该M×N维的列向量转成一个M行，N列的右视图的明暗分量I_rightshading；

S52、根据左右视图的明暗分量，通过图像成像约束，分别得到左视图的反射分量I_{leftreflectance}和右视图的反射分量I_{rightreflectance}；

S53、依次对左视图的明暗分量I_leftshading、右视图的明暗分量I_rightshading、左视图的反射分量I_{leftreflectance}和右视图的反射分量I_{rightreflectance}进行反归一化操作的，得到双目图像左右视图的明暗分量和反射分量，实现双目图像的本征性质分解。