CN106296749A - 基于l1范数约束的rgb‑d图像本征分解方法 - Google Patents

Abstract

本发明涉及计算机应用、图像处理技术，为实现易于求解，并能够展现出更多的细节，获取更好的图像处理结果，本发明采用的技术方案是，基于L1范数约束的RGB‑D图像本征分解方法，具体包括以下步骤：1)通过输入的彩色图及其对应的深度图计算得到图像每一点的色度、表面法向量；2)建立出优化方程的数据项部分E_data：3)建立整个优化方程的正则项E_reg；4)建立总的优化方程；5)由D、N、C三个部分得到亮度部分的log值矩阵S；6)得到反照率图a、亮度图s。本发明主要应用于图像处理场合。

Description

基于L1范数约束的RGB-D图像本征分解方法

技术领域

本发明属于计算机应用领域，具体讲，涉及基于L1范数约束的RGB-D图像本征分解方法。

背景技术

本征图像分解问题是一个在计算机视觉和图形图像领域长期存在的问题,旨在将输入的图像分解成几张不同成分图像，这些图像分别描述输入图像中的物体的本质的材质特性和环境光照，描述图像中真实的物理世界。最普遍的分解方法是将图像分解成一个albedo(反照率)图像和一个shading(亮度)图像。本征图像分解得到的结果将使很多应用变得可能，如图像重新光照、图像中物体表面重新上色、材质转换以及基于图像的材质编辑，这些应用当中很多都可以用于增强现实的场景当中。从一幅图片中估计物体的本征反照率是由Land和McCann提出的(Land E H,Mccann J J.Lightness and RetinexTheory.In Journal of the Optical Society of America,1971)，他们提出的Retinex模型奠定了本征图像分解的基础。Retinex模型提出了一个比较理想的蒙德里安世界——图像满足局部一致性，即一个平面画布的图像是由一块块的图像块形成的，在每个块内，反照率和亮度是恒定不变的，并且图像块所接收到的光照是缓慢变化的。那么由于光照是缓慢变化的，在这种图像的log域中，大的导数产生的地方可以被假设为相应的反照率边界。基于这个假设,Land和McCann提出了一个通过沿着两点间的路径进行积分来计算两点间相应的反照率的算法。这个算法被Horn从一维扩展到了两维(B.K.Hom.Determining lightnessfrom an image.In Computer graphics and image processing,1974)，之后Finlayson等人(G.D.Finlayson,S.D.Hordley,and M.S.Drew.Removing shadows from images[M].inComputer Vision—ECCV,2002)将同样的方法应用到移除图像阴影的算法中，通过对比色信息进行计算，然后把这些信息应用于寻找由阴影引起的导数，并且把这些由阴影引起的导数设置为0，剩下的导数进行分类处理，分类的时候将其分为亮度导数或反照率导数，再进行积分，就可以获得不含阴影的图像。尽管在这个问题上已经有很多的人做了大量的研究，想要高质量的分解结果依然是一个非常具有挑战性的任务。传统的算法解决每个像素上的shading和albedo,例如彩色Retinex理论,材质线索或者训练分类器。之后，人工辅助、图像序列法等方法也被引入，来提高分解质量。

RGB-D成像设备的商业化和普及，提供了一个机会来更进一步地研究本征图像分解问题，有可能在不需要人工辅助的情况下获得更高质量的分解结果。Qifeng Chen和Vladlen Koltun(Qifeng Chen and Vladlen Koltun.A simple model for intrinsicimage decomposition with depth cues.In ICCV,2013)通过对图像成像过程的更加细节性地分析，建立了一个效果很好的RGB-D图像本征分解模型。在他们的方法中，将图像首先分解为四个部分，直接光照部分，非直接光照部分，光照颜色部分，反照率部分，最后用前三个部分合成shading部分。并且实验结果表明，他们的方法要优于之前提出的其他的比较好的本征图像分解方法，平均的误差值比较小，大大提高了分解结果的质量。不过该方法中也存在着一些不足，比如在shading(亮度)图像中可能出现比较多地彩色信息、过度地估计光照颜色、亮度图中细 节可能有些模糊等。

发明内容

为克服现有技术的不足，实现易于求解，并能够展现出更多的细节，本发明采用的技术方案是，基于L1范数约束的RGB-D图像本征分解方法，具体包括以下步骤：

1)通过输入的彩色图及其对应的深度图计算得到图像每一点的色度、表面法向量，计算色度的公式为：其中,i_p表示输入的彩色图i在p点处的值，ch(i_p)表示p点处的色度值，r_p、g_p、b_p分别表示输入的彩色图在p点处红、绿、蓝三个通道上的值；

2)根据i_p＝a_pd_pn_pc_p，其中i_p、a_p、n_p、d_p、c_p分别表示输入的彩色图i、反照率图a、直接光照图d、非直接光照图n、光照颜色图c在p点处的值，对等式的两边求log，则得到：I_p＝A_p+D_p+N_p+C_p，其中I_p、A_p、D_p、N_p、C_p分别表示i_p、a_p、d_p、n_p、c_p的log值，由此建立出优化方程的数据项部分E_data：

$E_{data}=\underset{p}{\Σ}\left|\right|\sqrt{lum\left(i_p\right)}(I_p-A_p-\stackrel{\→}{1}{D}_p-\stackrel{\→}{1}{N}_p-C_p)|{|}_2^2---\left(1\right)$

其中，||.||₂表示L2范数,表示L2范数的平方，表示全1的向量，通过相乘，将D_p、N_p变为三维以便于计算。像素p的权重是通过像素p处的亮度lum(i_p)来计算，为了方便之后的求解，对数据项中的变量，分别用矩阵来定义，然后用矩阵形式来表示整个数据项，于是，数据项形式化为

$E_{data}=\left|\right|W(I-A-DK-NK-C)|{|}_2^2---\left(2\right)$

其中，

$W=diag(\sqrt{lum\left(i_1\right)},...,\sqrt{lum\left(i_n\right)})---\left(3\right)$

K=[1 1 1] (4)

diag(.)表示对角阵，所以W为一个n×n的对角矩阵，n表示图像的像素总个数，I、A、D、N、C分别表示输入的彩色图、反照率部分、直接关照部分、非直接光照部分、光照颜色组分的log值矩阵，也就是要求的未知数；

3)建立整个优化方程的正则项E_reg为：

E_reg＝∑_{j∈{A,D,N,N',C}}λ_jE_j (5)

E_j即E_A、E_D、E_N、E_N'、E_C，分别表示A、D、N、C的正则项，其中对于非直接关照的log值矩阵N，正则项由E_N、E_N'构成，λ_j即λ_A、λ_D、λ_N、λ_N'、λ_C分别表示正则项中A、D、N、C的正则项所占的权重；

4)建立总的优化方程为：

$\begin{array}{c}E=E_{data}+E_{reg}\\ =E_{data}+\underset{j\∈\{A,D,N,N^{\′},C\}}{\mathrm{\Σ}}{\mathrm{\λ}}_j{E}_j\\ =\left|\right|W(I-A-DK-NK-C)|{|}_2^2+{\mathrm{\λ}}_A|\left|Q_{1}A\right|{|}_1+{\mathrm{\λ}}_D\left|\right|Q_{2}D|{|}_1\\ +{\mathrm{\λ}}_N\left|\right|Q_{3}N|{|}_2^2+{\mathrm{\λ}}_{N^{\′}}|\left|Q_{3}N\right|{|}_2^2+{\mathrm{\λ}}_C\left|\right|Q_{5}C|{|}_2^2\end{array}---\left(16\right)$

通过最小化上式总能量E来求解各个分量成分的log值矩阵：A、D、N、C；

5)由D、N、C三个部分得到亮度部分的log值矩阵S：

S＝DK+NK+C (29)

其中，K＝[1 1 1]，亮度图s中每个像素上的log值，都是由D与K相乘、N与K相乘、C这三个矩阵对应位置的值相加得到的；

6)对反照率组分log值矩阵A和亮度部分的log值矩阵S的每一点处的值求自然指数： a_p、s_p、A_p、S_p分别表示反照率图A、亮度图S、反照率图log值矩阵A、亮度图log值矩阵S在p点处的值，e⁽ _. ⁾表示自然指数，由此得到反照率图a、亮度图s。

步骤3)中求解各正则项的具体步骤是：

3-1)建立反照率部分的正则项：其中，A_p、A_q分别表示反照率log值矩阵A在p、q点处的值，||.||₁表示L1范数，M_A是将每个像素简单地连接到图像中k个随机的点而得到的一个两两成对的像素点的集合，权重α_pq调节了正则项的力度，是由p、q两点之间色度的差值，以及亮度值的差异来构成的：

其中ch(i_p)表示输入彩色图上p点处的色度，lum(i_p)表示输入彩色图i_p处的亮度，max(.)表示求最大值，对于反照率部分的正则项，为了之后求解的方便，仍然要把它变成矩阵形式，对M_A中的每一对像素点，定义边e_pq，表示像素点p与像素点q之间的连接，则得到边的集合L_A：

L_A＝{e_p,q|(p,q)∈M_A} (6)

然后，定义一个矩阵Q₁，与L_A有相同的行数，列数为彩色图i中所有像素点的个数，矩阵Q₁中的每一行与L_A中的每一条边对应，每一列对应于图像I中的点，在Q₁的每一行，只有两个非零值，假设Q₁中第r^th行对应的边为e_pq,则第r行第p列值为第q列的值为反照率组分A的正则项则表示为：

E_A＝||Q₁A||₁ (7)

3-2)建立直接光照部分的正则项：D_p、D_q分别表示直接光照log值矩阵D在p、q点处的值，M_D是两两连接的“像素对”的集合：对于每个像素p，计算一个特征向量(x,y,z,n_x,n_y,n_z)；向量(x,y,z)是三维空间点的坐标，由p像素点在图像中 的坐标以及对应的深度信息来获得；向量(n_x,n_y,n_z)是在p点处表面的法线向量，由p点和临近的点的深度值来计算得到；由此，图像中所有的像素都被放置到一个六维的特征空间，为了使这些特征值规范化，在(x,y,z)三维做了白化变换，然后，对每个像素p,在这个特征空间里寻找k个最近点，对每个临近点q，添加像素对{p,q}到“像素对”集合M_D中，对于M_D中的每一对像素点，定义边e_pq，表示像素p与像素q之间的连接，则得到边的集合L_D：

L_D＝{e_pq|(p,q)∈M_D} (8)

类似于Q₁，定义一个矩阵矩阵Q₂中的每一行与L_D中的每一条边对应，每一列对应于图像i中的点，在Q₂的每一行，只有两个非零值，假设Q₂中第r行对应的边为e_pq,则第r行第p列值为1,第q列的值为-1，那么，直接光照组分D的正则项表示为：

E_D＝||Q₂D||₁ (9)

3-3)建立非直接光照部分的正则项E_N与E_N'，其中N_p、N_q分别表示非直接光照log值矩阵N在p、q点处的值，M_N是两两连接的“像素对”的集合：对N中的每个像素p,在矩阵N中寻找k个最近点，对每个临近点q，添加像素对{p,q}到“像素对”集合M_N中。将E_N写成矩阵形式的过程为：对于M_N中的像素点对{p,q}，定义边e_pq表示像素点p与像素点q之间的连接，则得到边的集合L_N：

L_N＝{e_pq|(p,q)∈M_N} (10)

根据边的集合L_N,构建矩阵Q₃,矩阵Q₃的每一行对应于L_N中的每一条边，每一列对应于图像中的每个像素点，每行只有两列值不等于0，然后，E_N就被表示为：

$E_N=\left|\right|Q_{3}N|{|}_2^2---\left(11\right)$

其中，||.||₂表示L2范数,表示L2范数的平方。

构建E_N'＝∑_pN_p ²，将E_N'表示为矩阵形式，首先构造矩阵Q₄:

Q₄＝diag(1,1,....,1) (12)

即Q₄是一个n×n的单位矩阵，E_N'就表示为：

$E_{N^{\′}}=\left|\right|Q_{4}N|{|}_2^2---\left(13\right)$

3-4)光照颜色部分的正则项为：C_p、C_q分别表示光照颜色log值矩阵在p、q点处的值，权重γ_p,q调整了约束项的力度，是根据p与q在三维空间内的位置的欧几里得距离来计算得到的: 表示该点的位置坐标，M_c是两两连接的“像素对”的集合，建立的方法为：连接每个像素点p到图像中k个随机的像素点；

对于M_C中的每一对像素点，定义边e_pq，表示像素点p与像素点q之间的连接，得到的边的集合L_C如下：

L_C＝{e_pq|(p,q)∈M_C} (14)

随后，定义矩阵Q₅，与L_C有相同的行数，列数为输入彩色图像i中所有像素点的个数，矩阵Q₅中的每一行与L_C中的每一条边对应，每一列对应于图像i中的点，在Q₅的每一行，只有两个非零值，假设Q₅中第r行对应的边为e_pq,则第r行第p列值为第q列的值为构建Q₅时使用是因为在求解的过程中用到的是Q₅ ^TQ₅，其中.^T表示矩阵的转置，那么，C的正则项则表示为：

$E_C=\left|\right|Q_{5}C|{|}_2^2---\left(15\right).$

步骤4)具体是利用增广拉格朗日方法进行最终求解，包括以下步骤：

4-1)把优化方程进行转换，令B₁＝Q₁A，B₂＝Q₂D，转换之后的能量最小化方程为：

$\begin{array}{c}{\mathrm{argmin}}_{X=(A,D,N,C)}\left|\right|W(I-A-DK-NK-C)|{|}_2^2\\ +{\mathrm{\λ}}_A\left|\right|B_1|{|}_1+{\mathrm{\λ}}_D|\left|B_2\right|{|}_1+{\mathrm{\λ}}_N\left|\right|Q_{3}N|{|}_2^2\\ +{\mathrm{\λ}}_{N^{\′}}\left|\right|Q_{4}N|{|}_2^2+{\mathrm{\λ}}_C|\left|Q_{5}C\right|{|}_2^2\end{array}---\left(17\right)$

s.t B₁＝Q₁A,B₂＝Q₂D

其中，X为A、D、N、C的一个组合，求解方法的目的是找到最优的X,使得能量函数的值最小；

4-2)用拉格朗日乘子法来对(4-1)中转换后的式子进行转换，列出相应的增广拉格朗日方程：

$\begin{array}{c}E(B_1,B_2,Y_1,Y_2,{\mathrm{\&mu;}}_1,{\mathrm{\&mu;}}_2)=\left|\right|W(I-A-D\&times;K-N\&times;K-C)|{|}_2^2\\ +{\mathrm{\&lambda;}}_A\left|\right|B_1|{|}_1+<Y_1,B_1-Q_{1}A>+\frac{{\mathrm{\&mu;}}_1}{2}||B_1-Q_{1}A|{|}_2^2\\ +{\mathrm{\&lambda;}}_D\left|\right|B_2|{|}_1+<Y_2,B_2-Q_{2}D>+\frac{{\mathrm{\&mu;}}_2}{2}||B_2-Q_{2}D|{|}_2^2\\ +{\mathrm{\&lambda;}}_N\left|\right|Q_{3}N|{|}_2^2+{\mathrm{\&lambda;}}_{N^{\&prime;}}|\left|Q_{4}N\right|{|}_2^2+{\mathrm{\&lambda;}}_C\left|\right|Q_{5}C|{|}_2^2\end{array}---\left(18\right)$

其中，μ₁、μ₂为两个常量，值为正数，Y₁、Y₂是拉格朗日乘子，<·,·>表示两个矩阵的内积,E是目标方程的总能量；

4-3)使用交替求解的迭代求解过程来求解，获得最优的A、D、N、C，第k+1次迭代过程如下:

${B_1}^{(k+1)}={\mathrm{argmin}}_{B_1}{\mathrm{\&lambda;}}_A\left|\right|B_1\left|\right|+<{Y_1}^{\left(k\right)},B_1-Q_1{A}^{\left(k\right)}>+\frac{{{\mathrm{\&mu;}}_1}^{\left(k\right)}}{2}\left|\right|B_1-Q_1{A}^{\left(k\right)}|{|}_2^2---\left(19\right)$

$\begin{array}{c}{A}^{(k+1)}={\mathrm{argmin}}_A\left|\right|W(I-A-D^{\left(k\right)}K-N^{\left(k\right)}K-C^{\left(k\right)}|{|}_2^2+<{Y_1}^{\left(k\right)},{B_1}^{(k+1)}>\\ +\frac{{{\mathrm{\&mu;}}_1}^{\left(k\right)}}{2}\left|\right|{B_1}^{(k+1)}-Q_{1}A|{|}_2^2\end{array}---\left(20\right)$

${B_2}^{(k+1)}={\mathrm{argmin}}_{B_2}{\mathrm{\&lambda;}}_D\left|\right|B_2\left|\right|+<{Y_2}^{\left(k\right)},B_2-Q_2{D}^{\left(k\right)}>+\frac{{{\mathrm{\&mu;}}_2}^{\left(k\right)}}{2}\left|\right|B_2-Q_2{D}^{\left(k\right)}|{|}_2^2---\left(21\right)$

$\begin{array}{c}{D}^{(k+1)}={\mathrm{argmin}}_D\left|\right|W(I-A^{\left(k\mathrm{+1}\right)}-DK-N^{\left(k\right)}K-C^{\left(k\right)}|{|}_2^2+<{Y_2}^{\left(k\right)},{B_2}^{(k+1)}>\\ +\frac{{{\mathrm{\&mu;}}_2}^{\left(k\right)}}{2}\left|\right|{B_2}^{(k+1)}-Q_{2}D|{|}_2^2\end{array}---\left(22\right)$

$\begin{array}{c}{N}^{(k+1)}={\mathrm{argmin}}_N\left|\right|W(I-A^{(k+1)}-D^{(k+1)}K-NK-C^{\left(k\right)}|{|}_2^2+{\mathrm{\&lambda;}}_N\left|\right|Q_{3}N|{|}_2^2\\ +{\mathrm{\&lambda;}}_{N^{\&prime;}}\left|\right|Q_{4}N|{|}_2^2\end{array}---\left(23\right)$

$C^{(k+1)}={\mathrm{argmin}}_C\left|\right|W(I-A^{(k+1)}-D^{(k+1)}K-N^{(k+1)}K-C|{|}_2^2+{\mathrm{\&lambda;}}_C\left|\right|Q_{4}C|{|}_2^2---\left(24\right)$

Y₁ ^(k+1)＝Y₁ ^(k)+(B₁ ^(k+1)-Q₁A^(k+1)) (25)

Y₂ ^(k+1)＝Y₂ ^(k)+(B₂ ^(k+1)-Q₂D^(k+1)) (26)

μ₁ ^(k+1)＝ρ₁μ₁ ^(k)ρ₁＞1 (27)

μ₂ ^(k+1)＝ρ₂μ₂ ^(k)ρ₂＞1 (28)

其中，.^(k)、.^(k+1)分别表示第k次、第k+1次迭代得到的结果，ρ₁、ρ₂为两个大于1常数系数，随着迭代次数的增加，使得μ₁、μ₂不断增大，对变量B₁、B₂、A、D、N、C这六个变量不断迭代最小化，不断地更新μ₁、μ₂、Y₁、Y₂直到总能量方程收敛，得到最优的A、D、N、C。

与已有技术相比，本发明的技术特点与效果：

本发明方法针对单张RGB-D图像进行本征分解，将一幅彩色图分解为反照率图、直接光照图、非直接光照图、光照颜色图，提出了使用L1范数表示反照率部分和直接光照部分的正则项，用tv-l₁模型进行建模，使用交替迭代求解，具有以下特点：

1、程序简单，易于实现。

2、由于0范数的非凸特性，使得求解变得非常困难，本发明采用L0范数的最优凸近似L1范数进行约束，L1范数最小化是凸优化问题，可以进行线性方程的求解。

3、分解得到的亮度图(shading)部分展现出更多的细节。

附图说明

本发明上述的和/或附加的方面和优点从下面结合附图对实施例的描述中将变得明显和容易理解:

图1为本发明方法的总框图；

图2为一张图片的分解效果图。

具体实施方式

本发明提出一种基于L1范数约束的RGB-D图像本征分解方法。具体包括以下步骤：

1)通过输入的彩色图及其对应的深度图计算得到图像每一点的色度、表面法向量，计算色度的公式为：其中,i_p表示输入的彩色图i在p点处的值，ch(i_p)表示p点处的色度值，r_p、g_p、b_p分别表示输入的彩色图在p点处红、绿、蓝三个通道上的值；

2)根据i_p＝a_pd_pn_pc_p，其中i_p、a_p、n_p、d_p、c_p分别表示输入的彩色图i、反照率图a、直接光照图d、非直接光照图n、光照颜色图c在p点处的值，对等式的两边求log，则得到：I_p＝A_p+D_p+N_p+C_p，其中I_p、A_p、D_p、N_p、C_p分别表示i_p、a_p、d_p、n_p、c_p的log值，由此建立出优化方程的数据项部分E_data：

$E_{data}=\underset{p}{\&Sigma;}\left|\right|\sqrt{lum\left(i_p\right)}(I_p-A_p-\stackrel{\&RightArrow;}{1}{D}_p-\stackrel{\&RightArrow;}{1}{N}_p-C_p)|{|}_2^2---\left(1\right)$

其中，||.||₂表示L2范数,表示L2范数的平方，表示全1的向量，通过相乘，将D_p、N_p变为三维以便于计算。像素p的权重是通过像素p处的亮度lum(i_p)来计算，为了方便之后的求解，对数据项中的变量，分别用矩阵来定义，然后用矩阵形式来表示整个数据项，于是，数据项形式化为

$E_{data}=\left|\right|W(I-A-DK-NK-C)|{|}_2^2---\left(2\right)$

其中，

$W=diag(\sqrt{lum\left(i_1\right)},...,\sqrt{lum\left(i_n\right)})---\left(3\right)$

K＝[1 1 1] (4)

diag(.)表示对角阵，所以W为一个n×n的对角矩阵，n表示图像的像素总个数，I、A、D、N、C分别表示输入的彩色图、反照率部分、直接关照部分、非直接光照部分、光照颜色组分的log值矩阵，也就是要求的未知数；

3)建立整个优化方程的正则项E_reg为：

E_reg＝∑_{j∈{A,D,N,N',C}}λ_jE_j (5)

E_j即E_A、E_D、E_N、E_N'、E_C，分别表示A、D、N、C的正则项，其中对于非直接关照N，正则项由E_N、E_N'构成，λ_j即λ_A、λ_D、λ_N、λ_N'、λ_C分别表示正则项中A、D、N、C的正则项所占的权重；

3-1)建立反照率部分的正则项：其中，A_p、A_q分别表示反照率log值矩阵A在p、q点处的值，||.||₁表示L1范数，M_A是将每个像素简单地连接到图像中k个随机的点而得到的一个两两成对的像素点的集合，权重α_pq调节了正则项的力度，是由p、q两点之间色度的差值，以及亮度值的差异来构成的：

其中ch(i_p)表示输入彩色图上p点处的色度，lum(i_p)表示输入彩色图i_p处的亮度，max(.)表示求最大值，对于反照率部分的正则 项，为了之后求解的方便，仍然要把它变成矩阵形式，对M_A中的每一对像素点，定义边e_pq，表示像素点p与像素点q之间的连接，则可以得到边的集合L_A：

L_A＝{e_pq|(p,q)∈M_A} (6)

然后，定义一个矩阵Q₁，与L_A有相同的行数，列数为彩色图i中所有像素点的个数，矩阵Q₁中的每一行与L_A中的每一条边对应，每一列对应于图像I中的点，在Q₁的每一行，只有两个非零值，假设Q₁中第r^th行对应的边为e_pq,则第r行第p列值为第q列的值为反照率组分A的正则项则表示为：

E_A＝||Q₁A||₁ (7)

3-2)建立直接光照部分的正则项：D_p、D_q分别表示直接光照log值矩阵D在p、q点处的值，M_D是两两连接的“像素对”的集合：对于每个像素p，计算一个特征向量(x,y,z,n_x,n_y,n_z)；向量(x,y,z)是三维空间点的坐标，由p像素点在图像中的坐标以及对应的深度信息来获得；向量(n_x,n_y,n_z)是在p点处表面的法线向量，由p点和临近的点的深度值来计算得到；由此，图像中所有的像素都被放置到一个六维的特征空间，为了使这些特征值规范化，在(x,y,z)三维做了白化变换，然后，对每个像素p,在这个特征空间里寻找k个最近点，对每个临近点q，添加像素对{p,q}到“像素对”集合M_D中，对于M_D中的每一对像素点，定义边e_pq，表示像素p与像素q之间的连接，则得到边的集合L_D：

L_D＝{e_pq|(p,q)∈M_D} (8)

类似于Q₁，定义一个矩阵矩阵Q₂中的每一行与L_D中的每一条边对应，每一列对应于图像i中的点，在Q₂的每一行，只有两个非零值，假设Q₂中第r行对应的边为e_pq,则第r行第p列值为1,第q列的值为-1，那么，直接光照组分D的正则项可以表示为：

E_D＝||Q₂D||₁ (9)

3-3)建立非直接光照部分的正则项E_N与E_N'，其中N_p、N_q分别表示非直接光照log值矩阵在p、q点处的值，M_N是两两连接的“像素对”的集合：对N中的每个像素p,在矩阵N中寻找k个最近点，对每个临近点q，添加像素对{p,q}到“像素对”集合M_N中。将E_N写成矩阵形式的过程为：对于M_N中的像素点对{p,q}，定义边e_pq表示像素点p与像素点q之间的连接，则得到边的集合L_N：

L_N＝{e_pq|(p,q)∈M_N} (10)

根据边的集合L_N,构建矩阵Q₃,矩阵Q₃的每一行对应于L_N中的每一条边，每一列对应于图像中的每个像素点，每行只有两列值不等于0，然后，E_N就被表示为：

$E_N=\left|\right|Q_{3}N|{|}_2^2---\left(11\right)$

其中，||.||₂表示L2范数,表示L2范数的平方。

构建E_N'＝∑_pN_p ²，将E_N'表示为矩阵形式，首先构造矩阵Q₄:

Q₄＝diag(1,1,....,1) (12)

即Q₄是一个n×n的单位矩阵，E_N'就表示为：

$E_{N^{\&prime;}}=\left|\right|Q_{4}N|{|}_2^2---\left(13\right)$

3-4)光照颜色部分的正则项为：C_p、C_q分别表示光照颜色log值矩阵在p、q点处的值，权重γ_p,q调整了约束项的力度，是根据p与q在三维空间内的位置的欧几里得距离来计算得到的: 表示该点的位置坐标，M_c是两两连接的“像素对”的集合，建立的方法为：连接每个像素点p到图像中k个随机的像素点；

对于M_C中的每一对像素点，定义边e_pq，表示像素点p与像素点q之间的连接，得到的边的集合L_C如下：

L_C＝{e_pq|(p,q)∈M_C} (14)

随后，定义矩阵Q₅，与L_C有相同的行数，列数为输入彩色图像i中所有像素点的个数，矩阵Q₅中的每一行与L_C中的每一条边对应，每一列对应于图像i中的点，在Q₅的每一行，只有两个非零值，假设Q₅中第r行对应的边为e_pq,则第r行第p列值为第q列的值为构建Q₅时使用是因为在求解的过程中用到的是Q₅ ^TQ₅，其中.^T表示矩阵的转置，那么，C的正则项则表示为：

$E_C=\left|\right|Q_{5}C|{|}_2^2---\left(15\right);$

4)建立总的优化方程为：

$\begin{array}{c}E=E_{data}+E_{reg}\\ =E_{data}+\underset{j\&Element;\{A,D,N,N^{\&prime;},C\}}{\mathrm{\&Sigma;}}{\mathrm{\&lambda;}}_j{E}_j\\ =\left|\right|W(I-A-DK-NK-C)|{|}_2^2+{\mathrm{\&lambda;}}_A|\left|Q_{1}A\right|{|}_1+{\mathrm{\&lambda;}}_D\left|\right|Q_{2}D|{|}_1\\ +{\mathrm{\&lambda;}}_N\left|\right|Q_{3}N|{|}_2^2+{\mathrm{\&lambda;}}_{N^{\&prime;}}|\left|Q_{3}N\right|{|}_2^2+{\mathrm{\&lambda;}}_C\left|\right|Q_{5}C|{|}_2^2\end{array}---\left(16\right)$

通过最小化上式总能量E来求解各个分量成分的log值矩阵：A、D、N、C。具体来说，利用增广拉格朗日方法进行最终求解，包括以下步骤：

4-1)把优化方程进行转换，令B₁＝Q₁A，B₂＝Q₂D，转换之后的能量最小化方程为：

$\begin{array}{c}{\mathrm{argmin}}_{X=(A,D,N,C)}\left|\right|W(I-A-DK-NK-C)|{|}_2^2\\ +{\mathrm{\&lambda;}}_A\left|\right|B_1|{|}_1+{\mathrm{\&lambda;}}_D|\left|B_2\right|{|}_1+{\mathrm{\&lambda;}}_N\left|\right|Q_{3}N|{|}_2^2\\ +{\mathrm{\&lambda;}}_{N^{\&prime;}}\left|\right|Q_{4}N|{|}_2^2+{\mathrm{\&lambda;}}_C|\left|Q_{5}C\right|{|}_2^2\end{array}---\left(17\right)$

s.t B₁＝Q₁A,B₂＝Q₂D

其中，X为A、D、N、C的一个组合，求解方法的目的是找到最优的X,使得能量函数的值最小。

4-2)用拉格朗日乘子法来对(4-1)中转换后的式子进行转换，列出相应的增广拉格朗日方程：

$\begin{array}{c}E(B_1,B_2,Y_1,Y_2,{\mathrm{\&mu;}}_1,{\mathrm{\&mu;}}_2)=\left|\right|W(I-A-D\&times;K-N\&times;K-C)|{|}_2^2\\ +{\mathrm{\&lambda;}}_A\left|\right|B_1|{|}_1+<Y_1,B_1-Q_{1}A>+\frac{{\mathrm{\&mu;}}_1}{2}||B_1-Q_{1}A|{|}_2^2\\ +{\mathrm{\&lambda;}}_D\left|\right|B_2|{|}_1+<Y_2,B_2-Q_{2}D>+\frac{{\mathrm{\&mu;}}_2}{2}||B_2-Q_{2}D|{|}_2^2\\ +{\mathrm{\&lambda;}}_N\left|\right|Q_{3}N|{|}_2^2+{\mathrm{\&lambda;}}_{N^{\&prime;}}|\left|Q_{4}N\right|{|}_2^2+{\mathrm{\&lambda;}}_C\left|\right|Q_{5}C|{|}_2^2\end{array}---\left(18\right)$

其中，μ₁、μ₂为两个常量，值为正数，Y₁、Y₂是拉格朗日乘子，<·,·>表示两个矩阵的内积,E是目标方程的总能量。

4-3)使用交替求解的迭代求解过程来求解，获得最优的A、D、N、C，第k+1次迭代过程如下:

${B_1}^{(k+1)}={\mathrm{argmin}}_{B_1}{\mathrm{\&lambda;}}_A\left|\right|B_1\left|\right|+<{Y_1}^{\left(k\right)},B_1-Q_1{A}^{\left(k\right)}>+\frac{{{\mathrm{\&mu;}}_1}^{\left(k\right)}}{2}\left|\right|B_1-Q_1{A}^{\left(k\right)}|{|}_2^2---\left(19\right)$

$\begin{array}{c}{A}^{(k+1)}={\mathrm{argmin}}_A\left|\right|W(I-A-D^{\left(k\right)}K-N^{\left(k\right)}K-C^{\left(k\right)}|{|}_2^2+<{Y_1}^{\left(k\right)},{B_1}^{(k+1)}>\\ +\frac{{{\mathrm{\&mu;}}_1}^{\left(k\right)}}{2}\left|\right|{B_1}^{(k+1)}-Q_{1}A|{|}_2^2\end{array}---\left(20\right)$

${B_2}^{(k+1)}={\mathrm{argmin}}_{B_2}{\mathrm{\&lambda;}}_D\left|\right|B_2\left|\right|+<{Y_2}^{\left(k\right)},B_2-Q_2{D}^{\left(k\right)}>+\frac{{{\mathrm{\&mu;}}_2}^{\left(k\right)}}{2}\left|\right|B_2-Q_2{D}^{\left(k\right)}|{|}_2^2---\left(21\right)$

$\begin{array}{c}{D}^{(k+1)}={\mathrm{argmin}}_D\left|\right|W(I-A^{\left(k\mathrm{+1}\right)}-DK-N^{\left(k\right)}K-C^{\left(k\right)}|{|}_2^2+<{Y_2}^{\left(k\right)},{B_2}^{(k+1)}>\\ +\frac{{{\mathrm{\&mu;}}_2}^{\left(k\right)}}{2}\left|\right|{B_2}^{(k+1)}-Q_{2}D|{|}_2^2\end{array}---\left(22\right)$

$\begin{array}{c}{N}^{(k+1)}={\mathrm{argmin}}_N\left|\right|W(I-A^{(k+1)}-D^{(k+1)}K-NK-C^{\left(k\right)}|{|}_2^2+{\mathrm{\&lambda;}}_N\left|\right|Q_{3}N|{|}_2^2\\ +{\mathrm{\&lambda;}}_{N^{\&prime;}}\left|\right|Q_{4}N|{|}_2^2\end{array}---\left(23\right)$

$C^{(k+1)}={\mathrm{argmin}}_C\left|\right|W(I-A^{(k+1)}-D^{(k+1)}K-N^{(k+1)}K-C|{|}_2^2+{\mathrm{\&lambda;}}_C\left|\right|Q_{4}C|{|}_2^2---\left(24\right)$

Y₁ ^(k+1)＝Y₁ ^(k)+(B₁ ^(k+1)-Q₁A^(k+1)) (25)

Y₂ ^(k+1)＝Y₂ ^(k)+(B₂ ^(k+1)-Q₂D^(k+1)) (26)

μ₁ ^(k+1)＝ρ₁μ₁ ^(k)ρ₁＞1 (27)

μ₂ ^(k+1)＝ρ₂μ₂ ^(k)ρ₂＞1 (28)

其中，.^(k)、.^(k+1)分别表示第k次、第k+1次迭代得到的结果，ρ₁、ρ₂为两个大于1常数系数，随着迭代次数的增加，使得μ₁、μ₂不断增大。对变量B₁、B₂、A、D、N、C这六个变量不断迭代最小化，不断地更新μ₁、μ₂、Y₁、Y₂直到总能量方程收敛，得到最优的A、D、N、C；

5)由D、N、C三个部分得到亮度部分的log值矩阵S：

S＝DK+NK+C (29)

其中，K＝[1 1 1]，亮度图s中每个像素上的log值，都是由DK、NK、C这三个矩阵对应位置的值相加得到的；

6)对反照率组分log值矩阵A和亮度部分的log值矩阵S的每一点处的值求自然指数： a_p、s_p、A_p、S_p分别表示反照率图A、亮度图S、反照率图log值矩阵A、亮度图log值矩阵S在p点处的值，e^(.)表示自然指数，由此得到反照率图a、亮度图s。

下表为对来自MPI-Sintel数据集中的五张图片的实验结果的定量评估数据，

	图像1	图像2	图像3	图像4	图像5	平均值
							DSSIM(反照率)	0.2808	0.1905	0.3785	0.2793	0.3302	0.2799
DSSIM(亮度)	0.2687	0.1822	0.2628	0.2787	0.2609	0.2609
							LMSE(反照率)	0.0152	0.0137	0.0211	0.0080	0.0165	0.0149
LMSE(亮度)	0.0287	0.0116	0.0115	0.0170	0.0400	0.0214

其中有两种定量评估标准：LMSE为局部均方误差，把图片分成若干个小的区域(在本实验中分为10个小的窗口)，然后在局部区域分别计算结果与标定数据对应的点之间的误差；DSSIM是由SSIM计算得到的：DSSIM＝(1-SSIM)/2，其中SSIM(structural similaritymeasurement)计算的是实验结果与标定数据之间的结构的相似度。

Claims (3)

1.一种基于L1范数约束的RGB-D图像本征分解方法，其特征是，包括以下步骤：

1)通过输入的彩色图及其对应的深度图计算得到图像每一点的色度、表面法向量，计算色度的公式为：其中,i_p表示输入的彩色图i在p点处的值，ch(i_p)表示p点处的色度值，r_p、g_p、b_p分别表示输入的彩色图在p点处红、绿、蓝三个通道上的值；

2)根据i_p＝a_pd_pn_pc_p，其中i_p、a_p、n_p、d_p、c_p分别表示输入的彩色图i、反照率图a、直接光照图d、非直接光照图n、光照颜色图c在p点处的值，对等式的两边求log，则得到：I_p＝A_p+D_p+N_p+C_p，其中I_p、A_p、D_p、N_p、C_p分别表示i_p、a_p、d_p、n_p、c_p的log值，由此建立出优化方程的数据项部分E_data：

其中，||.||₂表示L2范数,表示L2范数的平方，表示全1的向量，通过相乘，将D_p、N_p变为三维以便于计算，像素p的权重是通过像素p处的亮度lum(i_p)来计算，为了方便之后的求解，对数据项中的变量，分别用矩阵来定义，然后用矩阵形式来表示整个数据项，于是，数据项形式化为

其中，

K＝[1 1 1] (4)

diag(.)表示对角阵，所以W为一个n×n的对角矩阵，n表示图像的像素总个数，I、A、D、N、C分别表示输入的彩色图、反照率部分、直接关照部分、非直接光照部分、光照颜色组分的log值矩阵，也就是要求的未知数；

3)建立整个优化方程的正则项E_reg为：

E_reg＝∑_{j∈{A,D,N,N',C}}λ_jE_j (5)

E_j即E_A、E_D、E_N、E_N'、E_C，分别表示A、D、N、C的正则项，其中对于非直接关照的log值矩阵N，正则项由E_N、E_N'构成，λ_j即λ_A、λ_D、λ_N、λ_N'、λ_C分别表示正则项中A、D、N、C的正则项所占的权重；

4)建立总的优化方程为：

通过最小化上式总能量E来求解各个分量成分的log值矩阵：A、D、N、C；

5)由D、N、C三个部分得到亮度部分的log值矩阵S：

S＝DK+NK+C (29)

其中，K＝[1 1 1]，亮度图s中每个像素上的log值，都是由D与K相乘、N与K相乘、C这三个矩阵对应位置的值相加得到的；

6)对反照率组分log值矩阵A和亮度部分的log值矩阵S的每一点处的值求自然指数：a_p、s_p、A_p、S_p分别表示反照率图A、亮度图S、反照率图log值矩阵A、亮度图log值矩阵S在p点处的值，e^(.)表示自然指数，由此得到反照率图a、亮度图s。

2.如权利要求1所述的基于L1范数约束的图像本征分解方法，其特征是，步骤3)中求解各正则项的具体步骤是，

3-1)建立反照率部分的正则项：其中，A_p、A_q分别表示反照率log值矩阵A在p、q点处的值，||.||₁表示L1范数，M_A是将每个像素简单地连接到图像中k个随机的点而得到的一个两两成对的像素点的集合，权重α_pq调节了正则项的力度，是由p、q两点之间色度的差值，以及亮度值的差异来构成的：其中ch(i_p)表示输入彩色图上p点处的色度，lum(i_p)表示输入彩色图i_p处的亮度，max(.)表示求最大值，对于反照率部分的正则项，为了之后求解的方便，仍然要把它变成矩阵形式，对M_A中的每一对像素点，

定义边e_pq，表示像素点p与像素点q之间的连接，则得到边的集合L_A：

L_A＝{e_p,q|(p,q)∈M_A} (6)

然后，定义一个矩阵Q₁，与L_A有相同的行数，列数为彩色图i中所有像素点的个数，矩阵Q₁中的每一行与L_A中的每一条边对应，每一列对应于图像I中的点，在Q₁的每一行，只有两个非零值，假设Q₁中第r^th行对应的边为e_pq,则第r行第p列值为第q列的值为反照率组分A的正则项则表示为：

E_A＝||Q₁A||₁ (7)

3-2)建立直接光照部分的正则项：D_p、D_q分别表示直接光照log值矩阵D在p、q点处的值，M_D是两两连接的“像素对”的集合：对于每个像素p，计算一个特征向量(x,y,z,n_x,n_y,n_z)；向量(x,y,z)是三维空间点的坐标，由p像素点在图像中的坐标以及对应的深度信息来获得；向量(n_x,n_y,n_z)是在p点处表面的法线向量，由p点和临近的点的深度值来计算得到；由此，图像中所有的像素都被放置到一个六维的特征空间，为了使这些特征值规范化，在(x,y,z)三维做了白化变换，然后，对每个像素p,在这个特征空间里寻找k个最近点，对每个临近点q，添加像素对{p,q}到“像素对”集合M_D中，对于M_D中的每一对像素点，定义边e_pq，表示像素p与像素q之间的连接，则得到边的集合L_D：

L_D＝{e_pq|(p,q)∈M_D} (8)

类似于Q₁，定义一个矩阵矩阵Q₂中的每一行与L_D中的每一条边对应，每一列对应于图像i中的点，在Q₂的每一行，只有两个非零值，假设Q₂中第r行对应的边为e_pq,则第r行第p列值为1,第q列的值为-1，那么，直接光照组分D的正则项表示为：

E_D＝||Q₂D||₁ (9)

3-3)建立非直接光照部分的正则项E_N与E_N'，其中N_p、N_q分别表示非直接光照log值矩阵N在p、q点处的值，M_N是两两连接的“像素对”的集合：对N中的每个像素p,在矩阵N中寻找k个最近点，对每个临近点q，添加像素对{p,q}到“像素对”集合M_N中；将E_N写成矩阵形式的过程为：对于M_N中的像素点对{p,q}，定义边e_pq表示像素点p与像素点q之间的连接，则得到边的集合L_N：

L_N＝{e_pq|(p,q)∈M_N} (10)

根据边的集合L_N,构建矩阵Q₃,矩阵Q₃的每一行对应于L_N中的每一条边，每一列对应于图像中的每个像素点，每行只有两列值不等于0，然后，E_N就被表示为：

其中，||.||₂表示L2范数,表示L2范数的平方；

构建E_N'＝∑_pN_p ²，将E_N'表示为矩阵形式，首先构造矩阵Q₄:

Q₄＝diag(1,1,....,1) (12)

即Q₄是一个n×n的单位矩阵，E_N'就表示为：

3-4)光照颜色部分的正则项为：C_p、C_q分别表示光照颜色log值矩阵在p、q点处的值，权重γ_p，q调整了约束项的力度，是根据p与q在三维空间内的位置的欧几里得距离来计算得到的: 表示该点的位置坐标，M_c是两两连接的“像素对”的集合，建立的方法为：连接每个像素点p到图像中k个随机的像素点；

对于M_C中的每一对像素点，定义边e_pq，表示像素点p与像素点q之间的连接，得到的边的集合L_C如下：

L_C＝{e_pq|(p,q)∈M_C} (14)

随后，定义矩阵Q₅，与L_C有相同的行数，列数为输入彩色图像i中所有像素点的个数，矩阵Q₅中的每一行与L_C中的每一条边对应，每一列对应于图像i中的点，在Q₅的每一行，只有两个非零值，假设Q₅中第r行对应的边为e_pq,则第r行第p列值为第q列的值为构建Q₅时使用是因为在求解的过程中用到的是Q₅ ^TQ₅，其中.^T表示矩阵的转置，那么，C的正则项则表示为：

3.如权利要求1所述的基于L1范数约束的RGB-D图像本征分解方法，其特征是，步骤4)具体是利用增广拉格朗日方法进行最终求解，包括以下步骤：

4-1)把优化方程进行转换，令B₁＝Q₁A，B₂＝Q₂D，转换之后的能量最小化方程为：

其中，X为A、D、N、C的一个组合，求解方法的目的是找到最优的X,使得能量函数的值最小；

4-2)用拉格朗日乘子法来对(4-1)中转换后的式子进行转换，列出相应的增广拉格朗日方程：

其中，μ₁、μ₂为两个常量，值为正数，Y₁、Y₂是拉格朗日乘子，<·,·>表示两个矩阵的内积,E是目标方程的总能量；

4-3)使用交替求解的迭代求解过程来求解，获得最优的A、D、N、C，第k+1次迭代过程如下:

Y₁ ^(k+1)＝Y₁ ^(k)+(B₁ ^(k+1)-Q₁A^(k+1)) (25)

Y₂ ^(k+1)＝Y₂ ^(k)+(B₂ ^(k+1)-Q₂D^(k+1)) (26)

μ₁ ^(k+1)＝ρ₁μ₁ ^(k)ρ₁>1 (27)

μ₂ ^(k+1)＝ρ₂μ₂ ^(k)ρ₂>1 (28)

其中，.^(k)、.^(k+1)分别表示第k次、第k+1次迭代得到的结果，ρ₁、ρ₂为两个大于1常数系数，随着迭代次数的增加，使得μ₁、μ₂不断增大，对变量B₁、B₂、A、D、N、C这六个变量不断迭代最小化，不断地更新μ₁、μ₂、Y₁、Y₂直到总能量方程收敛，得到最优的A、D、N、C。