【论文解读】利用高光谱图像对场景反射率进行有效估计（Efficient Estimation of Reflectance Parameters from Imaging Spectropy）

---

前言

论文参考格式：

[1]Gu Lin,Robles-Kelly Antonio A,Zhou Jun. Efficient estimation of reflectance parameters from imaging spectroscopy.[J]. IEEE transactions on image processing : a publication of the IEEE Signal Processing Society,2013,22(9).

---

摘要

本文解决了从单组多光谱或高光谱图像中有效恢复反射率参数的问题。主要的技术点有三个 。
1.利用基于小波的估计器使用小波来恢复图像中的阴影系数；
2.通过约束优化方法恢复表面反射率和镜面系数；
3.最小二乘法更新光源功率谱。

Index Terms—Multispectral and hyperspectral imaging, scene analysis, inverse methods for imaging spectroscopy.

Ⅰ. 介绍

成像光谱数据可用于恢复场景形状和光度不变性，进而提供识别目的的手段。 此外，光度不变性可以与物体形状和光功率谱结合使用，以在数字媒体的可视化，创建和修改中进行高级应用。如下图所示 其中显示了在人像照片中替代皮肤反射率的结果，其流程为对高光谱数据进行皮肤识别后，改变皮肤区域反射率参数，而保持其他光度参数，如阴影参数、其他区域反射率、高光参数、光谱功率不变，最后对原光谱数据与修改后的光谱数据进行伪彩色合成。值得注意的是，后面油画上的皮肤并没有改变颜色，这正是因为高光谱数据所具有的特异性识别的优点。

尽管有前瞻性的应用，但从光谱数据中恢复光度不变量（阴影参数、表面反射率、镜面度、光谱功率）并不是一项简单的任务。这是因为每个场景都包含复杂的光源，对象不用弯曲程度和阴影情况会造成反射率和照明效果改变。对于三元色图像，即RGB图像，对场景的特异性辨识能力远远达不到高光谱图像。例如，带有RGB传感器的相机容易出现同色异谱现象，无法分析场景中与光照无关的材料的光度不变量。
在计算机视觉和图像处理中，文献大多都致力于对发亮或粗糙表面进行建模。通常依赖于双向反射率分布函数（Bidirectional Reflectance Distribution Function ,BRDF）。但是其缺少处理现实世界图像所需的通用性。大多数文献主要关注对粗糙表面对比度和细节表现力的提升，而对粗糙表面物体增亮的精确建模，仍是应用光学中极为困难但很活跃的领域。
现有的对于RGB图像建模方法对有几十个甚至上百个波段的光谱数据适用性差。 这主要是由于所涉及的高维图像数据的封闭形式解很复杂，而且与BRDF模型参数变得与波长相关。
在这篇文章中，我们提出一种基于小波的方法，该方法从图像辐射中恢复反射参数，即阴影参数，镜面度，光谱功率和表面反射率。

Ⅱ. 估计反射参数的方法

为了恢复反射参数，我们采用了双色模型。 长期以来，此模型已用于基于物理学的视觉中，用于表征非朗伯表面上的镜面反射。
该模型指出，图像中像素的辐射度可以表示如下：

$I\left({\lambda }_i,u\right)=g(u)L\left({\lambda }_i\right)S\left({\lambda }_i,u\right)+k(u)L\left({\lambda }_i\right)$ （1）

其中$L\left({\lambda }_i\right)$是光源波长功率谱，是第$i$个波段${\lambda }_i$的函数，$S\left({\lambda }_i,u\right)$是表面反射率，而$g(u)$和$k(u)$是阴影因子和镜面系数 在像素位置$u$ 它们与波长无关，并受照明光源和表面几何形状的控制。
在公式1中，右侧的第一项对应于漫反射，而第二项对应镜面反射。我们可以对模型进行如下直观的解释。对于镜面反射像素，即在阴影因子$g(u)$可以忽略不计的情况下，该表面起着反射镜的作用，其中辐射度与光源功率谱$L\left({\lambda }_i\right)$成正比。对于镜面反射像素，即$k(u)$接近于零，则表面将显示为阴影，并且反射率将成为光源上的乘法项。 公式1可变化为公式2如下：

$\frac{I\left({\lambda }_i,u\right)}{L\left({\lambda }_i\right)}=r\left({\lambda }_i,u\right)+k(u)$ （2）

其中$r\left({\lambda }_i,u\right)=g(u)S\left({\lambda }_i,u\right)$去掉了光源功率谱$L\left({\lambda }_i\right)$并与通过阴影因子$g(u)$的反射率成比例。
请注意，按照等式2中的模型，如果有其他参数，则可以很容易地重建镜面反射系数$k(u)$。 类似地，假设光源光谱和镜面反射系数$k(u)$可用，则可以将辐照度归一化并恢复$r\left({\lambda }_i,u\right)=g(u)S\left({\lambda }_i,u\right)$。 这揭示了阴影因子、镜面系数、表面反射率、光源功率谱的迭代解决方法。
现在，我们将假定反射模型参数的初始估计值。 使用这些初始值，我们可以通过三步迭代过程估算光度不变量参数。 首先，计算阴影因子，其次，使用约束优化方法计算表面反射率和镜面系数。 最终，可以从已经重建的上述参数恢复发光功率谱。 应用这三个交错的步骤，直到达到收敛为止。

A. 重建阴影因子

如上所述，我们采用小波来恢复等式2中的阴影因子。支持我们方法的思想是，可以将阴影因子$g(u)$看作是控制表面漫反射的因子。 结果，我们可以使用朗伯定律并设$g(u)=cos(\alpha ,u)$，其中α是光源方向和表面法线之间的夹角。 如果光源方向已知，则可以通过表面法线和光源方向的内积以简单的方式从表面法线计算出该角度。 在之后，将讨论如何消除与光源方向已知这个要求。 现在，我们假设已知光源方向，并专注于阴影因子的计算。
如果物体表面的3D表示已知，则阴影系数的计算就成为简单内积的问题。 不幸的是，恢复表面的3D表示是一项艰巨的任务。 小波可被用作基函数的完整正交集，因此尽管受到噪声破坏的影响，表面的3D表示仍可重建。
小波是在图像分析中使用的一组基函数。 由于关于小波分析的完整论述不在本文的讨论范围之内，因此我们想向感兴趣的读者推荐Mallat [2]的书

[2] S. G. Mallat, A Wavelet Tour of Signal Processing. San Diego, CA,USA: Academic Press, 1999.

或Refregier [3]的论文。

[3] A. Refregier, “Shapelets - I. A method for image analysis,” Monthly Notice R. Astronomical Society, vol. 338, no. 1, pp. 35–47, 2003.

3D表面的恢复是通过将表面法线场的倾斜与小波状体的倾斜相关性来实现的。具体可参考文献[4]。

[4] P. Kovesi, “Shapelets correlated with surface normals produce surfaces,”in Proc. Int. Conf. Comput. Vis., 2005, pp. 994–1001.

具体上，可以通过卷积来计算相关性，而可以使用图像平面上图像梯度的方向来重建表面法线的倾斜。 如[4]所示，3D表面重建$H$为小波 $tan\left(s_q\right)$的加权和，其中$s_q$是小波法线的倾斜度。 重建公式为：

$\begin{array}{rlr}H=& \sum _q|\tan (\alpha )|\ast \left|\tan \left(s_q\right)\right|\cdot \mid \cos \left({\tau }_f\right)\ast \cos \left({\tau }_q\right)& +\sin \left({\tau }_f\right)\ast \sin \left({\tau }_q\right)\mid \\ =& \sum _q|\tan (\alpha )|\ast \left|\tan \left(s_q\right)\right|\cdot C_{{\tau }_q}& \end{array}$ （3）

其中${\tau }_q$和${\tau }_f$分别对应于小波和目标对象的法线倾斜。 在公式3中，$\alpha$是表面法线与光源方向之间的角度，并且$\ast$表示卷积。
不幸的是，角度$\alpha$仍然是要估计的变量，因此，上述等式不能直接使用。 尽管如此，利用公式2，我们可以得出以下关系：
$\frac{1}{r\left({\lambda }_i,u\right)-r\left({\lambda }_j,u\right)}=\sec (\alpha ,u)\frac{1}{\left(S\left({\lambda }_i,u\right)-S\left({\lambda }_j,u\right)\right)}$ （4）
注意，$\tan (\alpha )$，$sec(\alpha )$和$sec(\alpha )-1$在$\left[0,\frac{\pi }{2}\right]$中具有类似的单调增加行为。 在下图中，我们展示了一个形状以及使用$\tan (\alpha )$和$sec(\alpha )-1$相应的3D重建。请注意，这两个函数产生的重建彼此之间非常相似，并且与实际结果一致。

结果，我们可以使用$sec(\alpha )-1$作为$\tan (\alpha )$的替代，并通过：

$H_{{\lambda }_i}\approx {\sum }_q{\mathbf{R}}_{{\lambda }_i}\ast \left|\tan \left(s_q\right)\right|\cdot C_{{\tau }_q}$ （5）

其中${\mathbf{R}}_{{\lambda }_i}$是一个矩阵，其索引为像素位置$u$，由
$\mid (S\left({\lambda }_i,u\right)-S\left({\lambda }_j,u\right)/\left(r\left({\lambda }_i,u\right)-r\left({\lambda }_j,u\right)\right)-1\mid$得出。
回想一下， $tan\left(s_q\right)$和$C_{{\tau }_q}$是已知的，可以从小波基函数和图像梯度轻松计算出。 从理论上讲，任意选择的参考频带${\lambda }_j$均可用于上述差异。 但是，由于测量噪声，$r\left({\lambda }_j,u\right)$和$S\left({\lambda }_j,u\right)$可能会偏离由公式2建模的理想情况。 在实践中，可以为图像中的每个波段进行表面重建，但会大大增加算法复杂度。
为解决表面重建中噪声的影响，我们假设测量误差是独立的并且为具有零均值的高斯噪声。可以证明，在这种假设下，估计值的方差可以通过平均大大减少。因此，我们使用由下式给出的$r(\cdot ,u)$和$S(\cdot ,u)$的平均值代替任意一个波段。

$r_u^{\ast }=\frac{1}{N}{\sum }_{j=1}^{N}r\left({\lambda }_j,u\right)S_u^{\ast }=\frac{1}{N}{\sum }_{j=1}^{N}S\left({\lambda }_j,u\right)$ （6）

其中$N$是图像中的波段总数。 对于公式5中的矩阵${\mathbf{R}}_{{\lambda }_i}$的计算，我们采用上述公式得出的结果$r_u^{\ast }$和$S_u^{\ast }$来替代对应于索引为J的任意波段$r\left({\lambda }_j,u\right)$和$S\left({\lambda }_j,u\right)$。
遵循每个波段的表面重建的类似原理,我们也可以将公式3中的$H$做一个替代，即取平均。利用下式中的表面法线场$H^{\ast }$，每一个像素的$cos(\alpha ,u)$，即阴影因子，可以计算得出。

$H^{\ast }=\frac{1}{N}{\sum }_{j=1}^N{H}_{{\lambda }_j}$ （7）

B. 表面反射率和镜面系数的计算

目前 可以利用光源功率谱和阴影因子重建反射率$S\left({\lambda }_j,u\right)$和镜面反射系数$k\left(u\right)$。 由于$cos(\alpha ,u)$是使用上一节中介绍的小波模型计算的，因此我们可以使用以下公式恢复表面反射率。

$S{\left({\lambda }_i,u\right)}^{\ast }=\frac{1}{\cos (\alpha ,u)}\left(\frac{I\left({\lambda }_i,u\right)}{L\left({\lambda }_i\right)}-k(u{)}^{\ast }\right)$ （8）

其中$k(u{)}^{\ast }$是重建的镜面反射系数，可用下式求得：

$k(u{)}^{\ast }={\min }_{k(u)\ge 0}\left\{{\sum }_{j=1}^N{\left|I\left({\lambda }_j,u\right)-L\left({\lambda }_j\right)\left(g(u)S{\left({\lambda }_j,u\right)}^{\ast }-k(u)\right)\right|}^2\right\}$

上面的表达式与前文介绍的双色模型以及镜面系数一致，都应为正实数值。 实际上，可以通过使用约束最小二乘法来实现镜面反射系数的恢复。在第下一节进行讨论。

C. 光源功率谱计算

为了估计照明情况，我们注意到公式1中的光源功率谱$L\left({\lambda }_i\right)$与像素$u$表面法线方向无关。 尽管相机感测到的辐射强度会受到阴影和其他变化（例如镜面系数）的影响，但光源的光谱功率分布不会改变。 这表明光源是整个图像中的全局光度变量。这个结论很重要，因为我们可以通过下式得到光源功率谱$L\left({\lambda }_i\right)$：

$L{\left({\lambda }_i\right)}^{\ast }=\arg {\min }_{L\left({\lambda }_i\right)}\{{\sum }_{u\in \mathcal{I}}\mid I\left({\lambda }_i,u\right)-{\left(\cos (\alpha ,u)S{\left({\lambda }_i,u\right)}^{\ast }+k(u{)}^{\ast }\right)L\left({\lambda }_i\right)|}^2\}$ （9）

$\mathcal{I}$是图像中所有像素的集合。 由于上述损失函数中唯一未知的是光谱$L\left({\lambda }_i\right)$，因此上述最小化变成了最小二乘平方，其解是唯一的，并且计算效率高。

\quad 1）模型简化：为了简化上述方法，注意到，在现实世界的图像中，许多像素通常显示出较大的镜面反射分量。 这在图像的镜面反射峰值处特别明显。 对于这些像素，实际上，镜面系数比阴影因子大得多，即 $k(u)\gg g(u)$。 这意味着对于那些镜面反射分量较大的像素，我们有$I(\lambda i,u)\approx k(u)L(\lambda i)$，这是通过将 $g(u)\approx 0$代入公式1得到的。
\quad 我们注意到，尽管根据双色反射模型中有漫反射项，但许多像素具有较大的镜面反射系数。 结果，对于包含$K$个镜面斑的集合$\mathcal{P}$的那些像素，我们可以忽略公式(1)中的漫反射项$g(u)L\left({\lambda }_i\right)$，并通过使下式来恢复光源功率谱，公式(10)是(9)的简化版本：

$L{\left({\lambda }_i\right)}^{\ast }=\arg {\min }_{L\left({\lambda }_i\right)}\left\{{\sum }_{u\in \mathcal{P}}{\left|I\left({\lambda }_i,u\right)-k(u{)}^{\ast }L\left({\lambda }_i\right)\right|}^2\right\}$ （10）

\quad 2）应用稳健统计：实际上，成像光谱有噪声破坏。 此外，对于我们的光源重建步骤，我们使用最小二乘最小化方案。 这种最小二乘法虽然有效，但对异常值很敏感。
\quad 为了解决这个问题，可以采用统计方法来代替最小二乘法。 统计方法通过引入对异常值具有鲁棒性的估计器，降低噪声影响。 可以使用许多可靠的估计器，例如maximum likelihood (M-estimators)，linear combination of order statistics (L-estimators)，rank transformation (R-estimators) 和 repeated M-estimators (RM-estimators)
为了将这些估计值应用于我们的光源恢复步骤，我们假设可以将这些像素位置$u\in \mathcal{P}$的图像辐射强度建模为真实镜面反射率$k(u)L\left({\lambda }_i\right)$ 和噪声分量$\eta (u)$。 图像辐射率变为$\mathcal{I}\left({\lambda }_i,u\right)=k(u)L\left({\lambda }_i\right)+\eta (u)$。 利用这个式子，我们可以设置误差补偿函数$g(\cdot )$，可替代公式10 如下所示：

L ( λ i ) ∗ = arg ⁡ min ⁡ L ( λ i ) { ∑ u ∈ P g ( I ( λ i , u ) − k ( u ) ∗ L ( λ ) } = arg ⁡ min ⁡ L ( λ i ) { g ( η ( u ) ) } \begin{aligned} L\left(\lambda_{i}\right)^{*} &=\arg \min _{L\left(\lambda_{i}\right)}\left\{\sum_{u \in \mathcal{P}} g\left(I\left(\lambda_{i}, u\right)-k(u)^{*} L(\lambda)\right\}\right.\\ &=\arg \min _{L\left(\lambda_{i}\right)}\{g(\eta(u))\} \end{aligned} L(λi​)∗​=argL(λi​)min​{u∈P∑​g(I(λi​,u)−k(u)∗L(λ)}=argL(λi​)min​{g(η(u))}​ （11）

D. 光源方向

在前面的部分中，我们假定光源方向$\overrightarrow{\mathbf{L}}$是已知的。 实际上，情况可能并非如此。 为了放宽限制，根据兰伯特定律(Lambert’s Law)，当光源方向为$\overrightarrow{\mathbf{L}}$时，索引为$u$的像素与$H^{\ast }$上的$p$点相对应，阴影系数$g(u)$为$\left({\overrightarrow{\mathbf{N}}}_p\cdot \overrightarrow{\mathbf{L}}\right)$，其中${\overrightarrow{\mathbf{N}}}_p$是$p$点在表面$H^{\ast }$上的的法线。
通过将$g(u)=\left({\overrightarrow{\mathbf{N}}}_p\cdot \overrightarrow{\mathbf{L}}\right)$代入公式1，我们得到$I\left({\lambda }_i,u\right)=\left({\overrightarrow{\mathbf{N}}}_p\cdot \overrightarrow{\mathbf{L}}\right)L\left({\lambda }_i\right)S\left({\lambda }_i,u\right)+k(u)L\left({\lambda }_i\right)$，其中我们使用的符号与第II-A节中引入的符号一致。
为了消除与已知光源方向相关的限制，可以相对于任意参考系$e={\left[e_x,e_y,e_z\right]}^T$测量场景中的光方向。如下图所示：

在图中，我们显示了参考轴，表面法线以及观察者和光源方向。这里我们让光源方向$\overrightarrow{\mathbf{L}}\equiv [0,0,1{]}^T$，使阴影系数满足关系$g(u{)}^{\ast }=\cos (\alpha ,u)=\left({\overrightarrow{\mathbf{N}}}_p\cdot \overrightarrow{\mathbf{L}}\right)$。
这与阴影法形成对比，后者的计算目标是H *给出的曲面。在此，$H^{\ast }$是通过来计算$g(u{)}^{\ast }=\cos (\alpha ,u)=\left({\overrightarrow{\mathbf{N}}}_p\cdot \overrightarrow{\mathbf{L}}\right)$的手段。 在从阴影法中，广泛使用参考轴是相机轴，即观察者方向对应于物体的深度。 如图右侧所示。相比之下，我们选择的参考轴$e$具有旋转表面的效果，不再相对于相机平面表示其深度。 选择光源轴为参考轴时，相对于观看者的表面$H^{\ast }$的深度将旋转。 但是，该旋转不会影响法线${\overrightarrow{\mathbf{N}}}_p$与光源方向$\overrightarrow{\mathbf{L}}$之间的角度$\alpha$。 由于我们的计算目标不是表面深度$H^{\ast }$而是角度$\alpha$，因此不需要光源先验的方向。

E. 将方法扩展到三色图像

我们的方法也可以应用于使用RGB相机捕获的三色图像。 在三色图像中，三个颜色通道中的每个通道都有其自己的光谱灵敏度函数，即，对于每个波段${\lambda }_i$，都有颜色匹配函数$Q_c(\cdot )$，其中 c = { R , G , B } c=\{R, G, B\} c={R,G,B}。可以用各波长的辐射度来表示某像素索引为$u$的$c$色响应，如下所示：

$I_c(u)={\sum }_{i=1}^{N}I\left({\lambda }_i,u\right)Q_c\left({\lambda }_i\right)$

在这里用$I_c(u)$来表示对应于具有颜色匹配函数$Q_c(\cdot )$的相机在像素$u$处的辐射$I\left({\lambda }_i,u\right)$总和。
当我们将上述方程式应用于双反射模型时，我们得到

$\begin{array}{rl}{I}_c(u)=& g(u)\left(\sum _{i=1}^{N}S\left({\lambda }_i,u\right)Q_c\left({\lambda }_i\right)L\left({\lambda }_i\right)\right)+k(u)\left(\sum _{i=1}^{N}L\left({\lambda }_i\right)Q_c\left({\lambda }_i\right)\right)\\ =& g(u)B_c(u)+k(u)L_c\end{array}$ （13）

其中

$B_c(u)={\sum }_{i=1}^{N}S\left({\lambda }_i,u\right)Q_c\left({\lambda }_i\right)L\left({\lambda }_i\right)L_c={\sum }_{i=1}^{N}L\left({\lambda }_i\right)Q_c\left({\lambda }_i\right)$

从公式13，我们可以得出结论，方法可以扩展到三元色图像，以重建上述变量$g(u)$，$B_c(u)$，$k(u)$和$L_c$。但是，这些变量与其波长索引的量之间存在两个重要区别。 首先，$L_c$不是真实的光源功率谱，而是它与颜色匹配函数$Q_c(\cdot )$的乘积和。 其次，变量$B_c(u)$是波长索引的反射率和颜色匹配函数加权的光源功率谱的乘积的和。 为了进一步恢复三色图像的反射率，即反照率${\rho }_c(u)$，对于像素为$u$的第$c$个颜色通道，我们使用以下公式

${\rho }_c(u)=\frac{B_c(u)}{L_c}$ (14)

Ⅲ. 实现方式

现在将注意力转向方法的实现。 为此，我们在算法1中给出了相应的伪代码。请注意，在图中，对于第二节中介绍的计算顺序，对计算顺序进行了一些更改。 输入光功率谱$L{\left({\lambda }_i\right)}^{\ast }$的初始估计值，图像 $\mathcal{I}$，用于计算小波的标度数$M$和阈值$ϵ$。 输出估计的光源功率谱$L_{new}^{\ast }$和数组$\mathbf{G}$，$\mathbf{K}$和$\mathbf{S}$中的光度参数，其索引u的值分别为为$g(u{)}^{\ast }$，$k(u{)}^{\ast }$和$S(\cdot ,u{)}^{\ast }$。算法伪代码如下:

在实现中，我们利用文献[5]中的过程选择一组“双反射”图像块。

[5] C. P. Huynh and A. Robles-Kelly, “A solution of the dichromatic model for multispectral photometric invariance,” Int. J. Comput. Vis., vol. 90, no. 1, pp. 1–27, 2010.

这是通过将图像细分成点阵来实现的，然后将点阵中每个矩形块中的像素拟合到二维超平面，即 一个双反射平面，由表面反射率和光源功率谱组成两个维度。 在这里，如果斑块包含一小部分像素，这些像素的辐射度与其在双反射平面上的投影紧密相关，则我们将其视为双反射性。 对于双反射平面计算的详细描述，我们请读者参考[6]。

[6] S. A. Shafer, “Using color to separate reflection components,” Color Res. Appl., vol. 10, no. 4, pp. 210–218, 1985.

有了这些图像块，再应用文献[7]中连续体去除，以便恢复光源的像素级估计。 光源的初始估计由所有连续去除的图像块的光谱的中位数给出。这遵循了中位数实际上是鲁棒估计量的观点。

[7] Z. Fu, A. Robles-Kelly, T. Caelli, and R. Tan, “On automatic absorption detection for imaging spectroscopy: A comparative study,” IEEE Trans. Geosci. Remote Sens., vol. 45, no. 11, pp. 3827–3844, Nov. 2007.

根据文献[3], 我们的伪代码通过计算小波的参数$s_k$和$C_{{\tau }_k}$开始。 在第3和第4行中，我们初始化$S_u^{\ast }$和$S\left({\lambda }_i,u\right)$，在第一次迭代中，$S\left({\lambda }_i,u\right)-S_u^{\ast }\equiv 1$。 注意到总是可以将$r\left({\lambda }_i,u\right)$限制在0到1之间。 可以对此假设进行物理解释。$r\left({\lambda }_i,u\right)$的有界性质反映了相机可以始终按比例缩放单位表面积的光的辐射功率，因此其最大值为1。在第5行中，我们初始化镜面反射系数$k(u)$，即：

在第5行中，我们初始化镜面反射系数$k(u)$。

$k(u{)}^{\ast }={\min }_i\left\{\frac{I\left({\lambda }_i,u\right)}{L\left({\lambda }_i\right)}\right\}$
这些参数的初始化允许在第8行中计算${\mathbf{R}}_{{\lambda }_i}(u)$。在第7行和第17行中，按照Section II-A中的方法计算$r_u^{\ast }$和$S_u^{\ast }$。在算法1中，我们利用$H^{\ast }$得到法线${\mathbf{N}}_p={\nabla H^{\ast }|}_{u\sim p}$ ，其中表达式$u\sim p$表示索引为$u\in I$的像素对应于该点$p\in H^{\ast }$ 。请注意：如果估算的光源功率谱的变化小于或等于阈值，则它将返回想求光度参数。 否则，它将转到第6行，以继续进行相应的迭代。 使用基于牛顿法的约束最小二乘优化产生更新后的镜面系数。

Ⅳ. 实验

本节中将说明本文方法对光度参数重建的有效性。 首先简要描述整节中使用的数据集。 然后，我们对成像光谱数据和三元色图像的光源功率谱和阴影因子的重建情况进行了定量分析。 我们显示了使用重建的反射率对皮肤识别的结果，对于三元色图像，则使用了朗伯反射率校正。
通过介绍光源和表面反射率替换的结果来完成本节。

A. 数据集

对于我们的实验，我们使用了四个高光谱和两个三元色图像数据集。 在下图中，我们显示了数据集的样本图像，其中高光谱图像已用伪彩色渲染。
前两个高光谱数据集是“人脸可见光”和“人脸NIR”图像。 这两个数据集均包含51个人类对象的图像，每个图像都是在10个具有不同方向和光谱功率的定向光源之一下捕获的。 光源有两种高度。 其中第一个放置在摄像头系统上方，第二个放置在与摄像头相同的高度。 调整灯光的主要方向，使其指向场景的中心。 图像是使用装有液晶可调滤波器的相机获取的，这些相机可以将多光谱图像分辨到10 nm。 人脸可见光数据库中的主题是在可见光（430–720 nm）波长范围内捕获的。 人脸NIR数据库中是在近红外（650–990 nm）波长范围内捕获的。
第三个高光谱数据集包括描绘风景的图像。 该景观数据集由14个场景组成，从城市到野外，其中包含天空，树木，植被和人造物体。 所有这些场景都被阳光完全照亮。 在每幅图像中，我们在采集时都将一个白色校准目标放在框架的一角，以便捕获照明功率谱的地面真实情况。 值得注意的是，为了避免在处理风景图像时我们的方法或替代方法有任何不公平的优势，在处理时已从图像中裁剪了此校准目标。
第四个高光谱数据集包括使用配备了声光可调谐滤波器的相机获得的图像。 该滤光片允许以1 nm的步长进行调谐。 数据集包括对应于五个类别的玩具和一组的图像。 通过围绕物体的垂直轴以10°的增量旋转对象，在十个视图中获取了每个玩具样本，而光谱数据仅在头和尾的两个不同视图中捕获。 在我们的数据集中，总共有62个玩具和32个硬币，这些硬币在多个视角上产生684个高光谱图像。 每个图像均包含550至1000 nm波长，间隔9nm，总共51个波段。 为了进行光度校准，我们还捕获了白色校准目标的图像，以获得整个场景中光源的功率谱。
其他两个数据集是三色的。第一个是蒙德里安图像数据集。其中包含743个对象的图像，这些对象在11种不同的光源下被照射，这些光源具有真实光源功率谱。 最后，我们将MSRA图像数据集用于我们的朗伯反射校正结果

B. 光源功率谱复原

首先将注意力转向光源功率谱的恢复。为此，我们使用了人脸可视化数据库，人脸近红外数据库，硬币和玩具数据库，景观数据库和蒙德里安数据集。 为了计算与光源功率谱真实值的的准确性，我们采用了两种相似性度量。 首先是光谱信息发散量 （Spectral Information Divergence，SID），这是一种信息理论量度，它考虑光谱之间的概率差异来评估相似性。SID被广泛用于处理高光谱图像数据。 本文还采用了欧几里德角 （Euclidean angle） ，其公式为：

$d_{\text{angle}}={\cos }^{-1}\left(\frac{{\sum }_{k=1}^{N}L\left({\lambda }_k\right)L_{gt}\left({\lambda }_k\right)}{\sqrt{\left({\sum }_{k=1}^{N}L{\left({\lambda }_k\right)}^2\right)}\sqrt{\left({\sum }_{k=1}^N{L}_{gt}{\left({\lambda }_k\right)}^2\right)}}\right)$ （15）

其中$L\left({\lambda }_i\right)$ 是估计的照度，而$L_{gt}\left({\lambda }_k\right)$是真实光源功率谱。
利用上一节所述具体实施方法，我们介绍了基于小波的估计器 （SBE） 得出的照度估计结果，其简化如II-C.1节所述。 为了使用数据的统计量，这里使用的正则化器是 Tukey双权函数 （Tukey bi-weight function） ，Andrews M估计器（the Andrews M-estimator）， 柯西正则化器（Cauchy regulariser） 和 Huber内核（Huber kernel）。 这些在公式11中用作惩罚函数$g(\cdot )$。在下图中，我们显示了实验中使用的惩罚函数的表达式。请注意，在下图的方程式中，我们将变量$\kappa$用作阈值参数。在所有实验中，我们将$\kappa$设置为相机动态范围的10％。

另外，光源的初始估计会影响我们的SBE方法产生的结果。 这是因为，如第三部分所述，我们的方法需要对光度参数进行初始估计。 结果，我们将双色块的中值使用与两种替代初始化方案产生的结果进行了比较。 其中第一个采用光源功率谱的统一初始化，这是$L\left({\lambda }_i\right)\equiv v$，其中$v$是图像中每个波段的最大辐射值集合的平均值。 第二种方案对每个双色块中的任一像素采用每个频带中的最大辐射度。 因此，在前一种情况下，我们在不使用均质斑块的情况下将每个波段的全局最大值设置为初始光源估计，而在后一种情况下，利用了光源的波段连续性。
在表I中，我们显示了在处理高光谱数据集和Mondrian三色数据集时使用前面提到的三种初始化方案的三种方法的光源恢复性能。 在表中，将每个数据集的最佳性能设置为粗体。 另外，提供了近红外（NIR）（即650 nm–990 nm）和可见（VIS）波段（即430nm-720nm）的测量误差。
人类受试者数据集的。 这将为使用NIR照明作为识别手段的NIR人脸识别系统提供一致的波长范围。 从表I中，我们观察到，一般来说，人类对象的NIR波段比可见光波段更准确。
此外，除NIR人脸图像之外的所有数据集，双色块中值方案都优于另外两种初始化方案。

对于我们所有的实验，我们将SBE方法的结果，在II-C.1节（SBE简化）中介绍的简化方法以及根据II-C.2节的统计数据与许多替代方法进行比较。 为了在考虑的方法之间提供公平的比较，采用均值，中位数，三边形，最佳25％和最差25％的平均值，以满足非高斯误差分布。 这一点特别重要，因为，正如我们稍后将看到的，跨数据集的误差测量值并不总是呈正态分布。
在表II–VI中，我们比较了我们的方法以及替代方法，通过Eculidean角和SID得出的性能度量。 为了清楚起见，在表中，对于每个误差测量，每个数据集最好Eculidean角性能使用了粗体字体。 请注意，我们的方法及其变体的性能优于人脸（VIS和NIR）数据集的替代方法。 当应用于景观数据集和蒙德里安数据集时，它们通常也与其他算法同等效果。
简化方法具有有效性：为了恢复光源功率谱，它使用了一个子集的像素，其图像辐射率与光源相关。 因此，通过采用不太容易受到噪声破坏的像素子集来估计光源。 这是因为大镜面系数的假设间接意味着，对于强度大的像素，光源几乎与图像辐射成正比，因此就噪声和入射到光源的光强度而言，SNR高。 尽管如此，我们的方法并未采用图像中最亮的像素，而是采用了 $I\left({\lambda }_i,u\right)\approx k(u)L\left({\lambda }_i\right)$ 的像素。
注意，这种简化不会影响反射率，阴影因子或镜面系数的恢复。 另外，统计方法通常会改善中位数和最佳平均值25％的结果。 此外，统计数据确实改善了人脸NIR和Mondrian数据集上总体均值的结果。 这意味着统计量使得误差分布偏斜较小，因此更符合高斯分布。

最后，在表VII中，我们显示了使用我们的方法与替代方法进行比较后，记录每个数据集中每个图像的平均处理时间（以秒为单位）。根据表中的时间，我们注意到，我们的简化方法比替代方法快约10.21倍。 这是由于我们更简单的优化方法以及小波估计器的使用有关。此外，对于我们的SBE（简化）方法，稳健估计器的应用在处理时间上略有增加。 这在高光谱数据集上尤为明显，增加了0.6秒，占计算总时间的1.01％。

C. 阴影因子恢复和朗伯反射校正

这一节展示对人脸数据集的重建阴影因子的定性结果。 我们还展示从MSRA数据集获取的三色图像上的朗伯反射校正的示例结果。 这意味着，对于多光谱数据，我们计算阴影因子$g(u)$，而对于三色图像，该任务涉及从发亮的表面去除镜面反射，并且要补偿粗糙物体上的部分变亮效果。由于阴影因子是表面法线和光源方向之间的角度的余弦值，因此多光谱图像上和三色数据上的处理方式相近。
在下图中，我们显示了阴影图，即与数据集中的对象相对应的图像细节中每个像素的$g(u)$值。 在该图中，左侧显示了输入多光谱图像的伪彩色图像，而其他三列分别显示了其他方法，我们的方法及其简化版本所产生的结果。

在图7中，我们显示了与我们从MSRA数据集获取的两个样本图像相对应的朗伯反射校正的结果。

此外，通过我们的方法恢复的阴影因子更加平滑，并能更好显示出阴影部分。 这在图6中对象的前额和图7中第一行的图像中尤为明显。还要注意，如另两幅图的第二列所示，由其他方法传递的阴影图像的背景有些噪点和扭曲。 此外，对于三元色图像，会受到MSRA数据集中纹理的影响。

D. 皮肤识别

现在，我们将注意力转向由我们的算法及其在识别任务中的应用所恢复的光谱反射率的不变性。 具体来说，我们专注于使用从图像中提取的表面反射率进行皮肤区域分割。 该任务可以看作是分类问题，其中皮肤和非皮肤光谱包含正和负类。 为了识别皮肤，我们人类数据集中每个对象的图像中选择一个皮肤区域和非皮肤区域。 我们使用在一个光源下拍摄的每个对象的单个图像，并将归一化的反射光谱存储为训练数据。随后，我们在训练集上训练了具有2级多项式内核的支持向量机（SVM）分类器。将所得的SVM模型应用于其余图像中的皮肤像素与非皮肤像素进行分类，这些像素是在各种照明光谱条件下获取的。 通过这种方式，我们可以断言我们的算法重建的光源功率谱与表面反射率具有鲁棒性。
在图8的左侧栏中，我们显示了用于训练的图像的样本，该样本带对应于数据集中的一个主题。 第二列显示了在不同照明条件下相同对象的测试图像。其他方法得到的测试图像中皮肤概率图显示在第三列中，而最右边两列显示了使用我们的方法及其简化版本计算的反射率产生的皮肤概率图。 在这些皮肤概率图中，较亮的像素更有可能是皮肤。图8提供了一个概念证明，即恢复的皮肤反射光谱对于光源变化和光源方向变化是不变的。 确实，我们方法的皮肤分割图显示了正确分类的皮肤像素的很大一部分。 此外，非皮肤的面部细节（例如眉毛和嘴巴）也可以与皮肤区分开。 相反，其他方法在皮肤区域产生了更多的错误，在场景其他材料中也产生了错误。

我们在表VIII中提供了更定量的分析。
在表中，我们通过5折交叉验证显示了上述皮肤分割方案的分类率（CR）和错误检测率（FDR）。 在表中，我们还显示了将原始辐射率用作通过我们的SBE方法或[15]中的替代方法获得的反射率的替代方法时的准确性。 分类率是正确分类的皮肤像素的总体百分比。 错误检测率是被错误分类为皮肤的非皮肤像素的百分比。 该表显示了在对由正面光源照射的图像进行测试时，数据集中所有对象的CR和FDR。 不出所料，通过本文的方法获得的反射率可实现最高的皮肤识别率。

E. 反射率和光源替换

最后，展示该方法对与图像编辑相关任务。 在本节中，我们说明如何将从光谱数据中恢复的光度参数用于替换场景中对象的反射率或更改光源功率谱。
为了替换物体的反射率或更改场景中的光源，我们首先使用我们的方法恢复光度参数。 一旦有了反射率，光源功率谱，阴影系数和镜面反射系数，我们就可以通过评估公式2中的模型来恢复新图像。 在这组实验中，我们使用了光谱仪来获取某些天然材料的反射率。 这些反射光谱被用作所考虑对象的替代光谱，以恢复编辑后的图像。 对于光源，我们使用了配备有积分球的光谱仪来获取现实世界中常见的光源的功率谱，例如日光，钨光源，氙气-克里普顿灯等。
我们已经在人类受试者数据集中对三个受试者的衣服进行了反射识别。 在这里，我们使用了通过我们的方法提取的反射率，并在考虑中的布料上应用了目标材料识别算法，以选择反射率将被替换的像素。 这些布料反射率替换操作的结果显示在下图的第二行中。请注意，对于第二个场景，使用替代反射率时折痕更为明显。 当然，所有结果没有进行进一步的编辑或润饰。 当使用替代反射率时，折痕显得更加明显的原因在于，这些现象确实在多光谱图像的某些波段中可见。 当由于T恤的深色（黑色）而产生伪彩色图像时，这些带的边缘作用很小。 尽管如此，对于具有米色的替代反射率，这些折痕会出现在伪彩色图像上。 这与以下观点一致：我们的方法可以为形状，反射率和其他光度学参数提供单独的输出，因此，反射率替换不会影响T恤衫中的折痕。

最后，我们在实验中使用的三个对象上显示了替换光源功率谱的结果。 在图9中，底部一行图像显示的是顶部图像经过光源功率谱替换的结果。 请注意，图像的色温会发生变化，趋向于由霓虹灯管引起的偏绿色调。

Ⅴ. CONCLUSION

本文提出了一种有效的基于小波的估计器，用于恢复双反射模型光度参数，以解决多光谱图像中的光度不变性问题。 重建过程基于小波，并使用三步过程，可以通过约束最小二乘优化来解决。 我们还介绍了该方法的简化版本，该方法采用镜面像素，以重建光源功率谱。 事实证明，这种简化在我们的实验中是有效的，并提供了不使用最小二乘优化的方法。 我们还运用了统计方法，并在光源功率谱重建，阴影因子计算和皮肤识别方面进行了许多实验。比较了本文方法和其他方法，并显示了简单的光源功率谱替换工作，这些示例说明了该方法在数字媒体生成，图像编辑和修饰中的潜力。