CN102597744B - 照明谱恢复 - Google Patents

Abstract

本公开内容涉及电子图像(例如超谱或多谱图像)的处理。图像的处理包括：通过以下方式来进行对象识别和图像增强：基于这些类型的图像提供的波长索引的谱数据而利用从该图像可用的信息。使用基于二色反射模型和约束项的代价函数来估计(图21)这样的图像的照明谱。可以对图像的子图像执行该方法。还公开了一种用于选择这些子图像(图23)的方法。还公开了一种给定估计的照明谱而确定图像的光度测定的参数的方法(图22)。

Description

照明谱恢复

相关申请的交叉引用

本申请要求于2009年9月3日提交的来自澳大利亚临时专利申请No.2009904230的优先权，其内容通过引用合并于此。

本申请要求于2010年2月1日提交的来自澳大利亚临时专利申请No.2010900383的优先权，其内容通过引用合并于此。

本申请涉及2009年6月19日提交的国际(PCT)申请PCT/AU2009/000793并且公布为WO2009/152583，其内容通过引用合并于此。

本申请涉及与本申请共同提交、并且要求于2009年9月3日提交的澳大利亚临时专利申请No.2009904219以及于2010年4月22日提交的澳大利亚临时专利申请No.2010901785的优先权的国际(PCT)申请，其内容通过引用合并于此。

本发明涉及与本发明共同提交、并且要求于2009年9月3日提交的澳大利亚临时专利申请No.2009904218的优先权的国际(PCT)申请，其内容通过引用合并于此。

技术领域

本发明涉及电子图像(例如超谱或多谱图像)的处理。图像的处理包括对象识别和图像增强。具体地说，但并不限制于，还公开了一种用于确定图像的照明谱，估计光度测定的参数并且选择图像的子图像以用于图像处理的方法、软件和计算机。

背景技术

在多谱成像中，光度测定的不变量在材料标识[14]和形状分析的领域中提出很大机会和挑战。这归因于在穿过宽谱的几千个带中多谱和超谱感测设备(其传递波长索引的谱数据，即反射率(reflectance)或辐射率(radiance)数据)所获取的表面辐射率的信息丰富表现。如在三色成像中那样，不仅由光源和取景器方向，而且还由处于研究中的表面的材料特性来确定多谱反射率分布。

在计算机视觉中，为了表面分析和图像理解的目的，表面反射率的建模是关键重要性的课题。例如，Nayar和Bolle[42]已经使用光度测定的不变量来识别具有不同反射率特性的对象。

由于照明谱和材料反射率的恢复是相互独立的，因此色彩恒定性(colourconstancy)也是个问题。色彩恒定性是用于分辨固有材料反射率与它们在变化照明条件下捕获的三色彩色图像的能力。关于色彩恒定性的研究分支为两个主要趋势，其中一个依赖于发光体和材料反射率的统计，而另一个利用了表面材料的局部阴影和镜面性的基于物理的分析。

本发明涉及基于物理的色彩恒定性，其为了发光体和表面反射率估计的目的而分析光与材料相互作用的物理过程。基于物理的色彩恒定性的两种方法是Land的视网膜理论[35，36]和二色反射模型[48]。Land的视网膜理论已经激发了人类色彩恒定性的若干计算模型[7]。另一方面，二色模型描述了反射光作为主体反射和反射(高光)的组合，其将照明估计问题处理为从闪耀的表面分析高光[32，33，39，53]。

US7443539(Olympus)确定图像的谱数据特征并且需要分光计(42)和基准板(51)来补偿发光体。

本说明书中已经包括的文献、法案、材料、设备、物品等的任何讨论目的仅仅是提供本公开内容的上下文。并不看作承认由于在本申请的每一权利要求的优先权日期之前存在任何或所有这些内容，因此其形成现有技术基础的部分或是与本公开内容有关的领域中的公知常识。

发明内容

在第一方面中，提供一种计算机实现的方法，用于估计包括波长索引的谱数据(wavelength indexed spectral data)的点的图像的照明谱，所述方法包括：

基于以下来优化代价函数：

(i)波长索引的谱数据与基于针对不同二色反射模型参数值的二色反射模型的估计之间的差异，其中，所述二色反射模型至少包括照明谱参数，以及

(ii)用于基于强健内核或曲率一致性而约束模型参数值的项。

由于所述项，可以通过封闭形式来求解代价函数的优化。以此方式，本方法能够实现对于优化问题的高效解，并且具有更高的精度。

所述项可以基于：最小化阴影的梯度大小的变化或曲率大小变化。也就是说，所述项可以基于作为平滑性误差的测度的阴影的梯度大小或曲率变化。

由于所述项，与仅基于差异的平方的代价函数相比，通过向噪声破坏和图像获取误差提供减少的偏差，所述代价函数的优化可以是相对轮廓强健的。

由于所述项，由于噪声引起的波长索引的谱数据中的误差相对于基于二色模型的估计可以具有减少的差异。

内核的偏差基于图像中的点的阴影的确定的梯度大小与图像中的该点的局部邻域的阴影的确定的梯度大小之间的差异。可以对于每一点确定该情况。

内核函数可以在图像上估计的表面法线上使得代价函数相对轮廓是强健的。所述内核函数可以基于Huber估计者、Tukey双加权函数或Lee自适应势函数。

所述项可以将图像中的点的曲率约束为与图像中的该点的局部邻域的曲率相似。可以对于每一点确定该情况。

所述项可以偏好，即将较大权重分配给属于图像中的相同或相似形状或辐射率的点。

所述项可以基于图像中的阴影或辐射率所引起的曲率，例如，其可以使用主曲率。

所述项可以基于限制邻近谱波长之间的一个或多个模型参数值之间的差异。

在第二方面中，提供一种计算机实现的方法，用于估计包括波长索引的谱数据的点的图像的照明谱，所述方法包括：

基于以下来优化代价函数：

(i)所述波长索引的谱数据与基于针对不同二色反射模型参数值的二色反射模型的估计之间的差异，其中，所述二色反射模型至少包括照明谱参数，通过将所述波长索引的谱数据拟合为其中所述超平面约束照明谱参数值的所述二色反射模型的超平面来确定所述差异，以及

(ii)用于约束模型参数值的项。

基本向量可以跨越超平面约束照明谱参数值和反射率参数值。

由于变量的改变，可以通过封闭形式来求解照明谱、材料反射率、阴影和镜面性。

在第三方面中，提供一种计算机实现的方法，用于估计包括波长索引的谱数据的点的图像的照明谱，所述方法包括：

基于以下来优化代价函数：

(i)波长索引的谱数据与基于针对不同二色反射模型参数值的二色反射模型的估计之间的差异，其中，所述二色反射模型至少包括照明谱参数，以及

(ii)用于约束模型参数的项，

其中，迭代地优化所述代价函数，并且以均匀照明谱或基于仅基于所述差异的另一代价函数的解来初始化所述照明谱参数值。

这种初始化的优点在于，减少了估计的照明谱朝向不期望的平凡解的偏差。

另一代价函数可以根据Finlayson & Schaefer提出的方法[21]。

与第一、第二和第三方法相关，提供以下可选特征。

图像可以是表示捕获的场景的超谱或多谱图像。替代地，所述图像可以是从超谱或多谱数据得到的三色图像。替代地，所述图像可以是其中波长索引基于彩色通道的三色图像。

可以通过将波长索引的谱数据拟合为表示所述二色反射模型的超平面来确定基于所述二色反射模型的估计，所述超平面至少具有照明谱作为基础。超平面也可以具有漫射辐射谱作为基础。

所述二色反射模型可以还包括以下中的一个或多个：

表面阴影参数，

表面反射率参数，以及

镜面性参数。

优化代价函数还包括：迭代地优化针对参数的第一子集然后参数的第二子集的代价函数，其中，参数的第一子集包括照明谱参数，参数的第二子集包括表面阴影参数、表面反射率参数和镜面性参数中的一个或多个。

由于约束(即正则化)项没有照明谱分量，因此相对于参数的第一子集而优化代价函数可以独立于约束项。

相对于参数的第二子集而优化所述差异可以依赖于约束项。

在图像是子图像的集合的情况下，独立于其它子图像的参数的第二子集的值而优化子图像的参数的第二子集的值。可以根据下述第五方面来标识子图像。因此，子图像可以具有平滑阴影和均匀表面反射率(即同质的)。

代价函数可以评估图像的每一点(例如位置或像素)。代价函数可以是线性组合，约束项是几何项。代价函数可以包括用于平衡所述差异和约束项的平衡因子。

所述方法可以还包括：将优化代价函数的照明谱参数值识别为图像的照明谱的估计，例如提供估计作为输出，或存储该估计。所述方法可以包括：当照明谱参数值改变不大于一个或多个迭代之间的预先确定的量时，完成优化。

所述方法可以还包括：在与图像关联的进一步处理中使用图像的估计的照明谱，所述处理为以下中的一个或多个：

图像编辑：

形状估计；

材料识别或分类；以及

硬件校准。

在第四方面中，提供一种计算机实现的方法，用于估计包括波长索引的谱数据的点的图像的表面反射率、表面阴影和镜面性分量，其中，图像的照明谱是已知的，所述方法包括：

(a)对于图像中的每一点，基于二色反射模型、照明谱和波长索引的谱数据来估计该点的表面反射率；以及

(b)对于所述图像中的每一点，基于波长索引的谱数据、二色反射模型、照明谱以及估计的表面反射率来估计表面阴影和镜面性。

第四方面的优点是，可以通过以下方式来估计光度测定的参数：使用估计的谱反射率来估计或重新估计表面阴影和镜面性。这利用了表面反射率一般比表面阴影和镜面性更稳定。

可以对于图像的点的子集执行所述方法。

以上直接步骤(b)可以包括：

对于图像中的每一点，基于图像中的该点的局部邻域的波长索引的谱数据来估计该点的反射率谱；以及

将图像的每一点的波长索引的谱数据拟合到照明谱和反射率谱跨越的子空间上。

步骤(b)可以还包括：最小化波长索引的谱数据与照明谱和反射率谱跨越的子空间之间的差异。

估计该点的表面反射率的步骤(a)可以基于图像中的该点的局部邻域的波长索引的谱数据，从而局部邻域显现为缓慢地改变表面反射率。

所述方法可以还包括初始步骤：根据第一、第二或第三方面中的一个或多个的方法来估计照明谱。

在第五方面中，提供一种计算机实现的方法，用于选择用于在估计图像的照明谱中使用的图像的子图像，其中图像包括波长索引的谱数据的点，所述方法包括：

(a)标识图像的多个候选子图像；

(b)对于每一候选子图像，将基于二色反射模型的超平面拟合为候选子图像的波长索引的谱数据；以及

(c)选择具有不偏离于拟合超平面大于确定的量的波长索引的谱数据的候选子图像。

有利的是，可以在任何任意图像上以此方式自动地选择均匀反射率的子图像。然后可以对子图像使用分析图像的其它技术。

步骤(c)可以还包括：基于遍及波长的均值波长索引的谱数据的方差而选择也具有高对比度的子图像。

步骤(c)的确定的量可以是具有确定百分比的不偏离于超平面的波长索引的谱数据的候选子图像。

步骤(b)可以包括：拟合波长索引的谱数据的辐射率向量表示。

在第六方面中，提供一种计算机实现的方法，用于估计图像的照明谱，所述方法包括：

(a)根据第五方面的方法而选择图像的子图像；以及

(b)根据第一方面、第二方面或第三方面的方法估计图像的照明谱，其中，在估计中仅使用所选择的子图像。

在第七方面中，可以组合地执行第一方面、第二方面、第三方面、第四方面和第五方面的方法中的两个或更多个。第一方面、第二方面、第三方面、第四方面或第五方面的方法的可选特征也是组合的方法的可选特征。

在第八方面中，提供了当安装在计算机上时使得计算机执行以上第一方面、第二方面、第三方面、第四方面或第五方面中的任何一个或多个的方法。第一方面、第二方面、第三方面、第四方面或第五方面的方法的可选特征也是软件的可选特征。

在第九方面中，提供一种计算机，其包括处理器，用于执行以上第一方面、第二方面、第三方面、第四方面或第五方面的方法中的一个或多个的方法。第一方面、第二方面、第三方面、第四方面或第五方面的方法的可选特征也是计算机的可选特征。计算机的示例包括单机个人计算机或仅处理器(例如相机中合并的处理器)。

对于所有以上方面，点可以是一个或多个像素。波长可以是单一波长或波长带。

在一些实施例中，图像的特性/分量/参数(例如照明谱或阴影因子)的估计可以足够精确，以被看作确定。

附图说明

现将参照附图描述示例，其中：

图1示出以下说明书中使用的注记。

图2示出用于从同质块片(patch)的集合估计二色变量的算法1的伪代码。

图3示出用于估计发光体谱已知的图像的阴影(shading)、镜面性(secularity)和反射率(reflectance)的算法2的伪代码。

图4示出与Finalyson和Schaefer的方法(F & S)以及Sunshine等人的方法相比的以度为单位在可见和近红外谱下捕获的多谱面部图像数据库上具有精度与发光体估计方法中使用的块片的数量的关系的比较的表1。

图5示出与Finlayson和Schaefer的方法(F & S)和Sunshine等人的方法相比的以度为单位在多谱面部图像数据库上具有精度与照明估计示例中使用的块片的数量的关系的比较的表2，示出了针对可见(左边)和近红外范围(右边)的结果；

图6图形地示出与两种其它已知的方法相比的精度与多谱面部图像数据库上的示例使用的块片的数量的关系；

图7示出与Finlayson和Schaefer的方法(F & S)和Sunshine等人的方法相比的以度为单位的在可见和近红外谱下捕获的多谱室外图像数据库上精度与我们的发光体估计方法使用的块片的数量的关系的比较的表3；

图8图形地示出与Finlayson和Schaefer的方法(F & S)和Sunshine等人的方法相比的以度为单位的精度与多谱室外图像数据库上对于发光体估计示例中的我们使用的块片的数量的关系，示出了针对可见(左边)和近红外范围(右边)的结果；

图9示出与Finlayson和Schaefer的方法(F & S)和Sunshine等人的方法相比的以度为单位在合成的RGB面部图像数据库上其中α＝10000的发光体估计示例的精度的比较的表4，我们示出关于从Stiles和Burch的色彩匹配函数(行2-4)、佳能10D的谱灵敏度响应函数(行5-7)以及尼康D70相机(行8-10)的最终谱面部显像合成的RGB面部图像的结果；

图10图形地示出与Finlayson和Schaefer的方法(F & S)和Sunshine等人的方法相比的以度为单位的从多谱面部显像合成的RGB面部图像上精度与其中α＝10000的我们的发光体估计方法的图像块片的数量的关系，从左到右：通过(a)Stiles和Burch的色彩匹配函数(b)佳能10D相机传感器(c)尼康D70相机传感器捕获的仿真的RGB图像的结果；

图11示出与Finlayson & Schaefer的方法(F & S)和Sunshine等人的方法相比的以度为单位在合成的RGB室外图像数据库上具有0：＝100的我们的发光体估计方法的精度的比较的表5，我们示出关于从Stiles和Burch的色彩匹配函数(行2-4)、佳能100的谱灵敏度响应函数(行5-7)以及尼康070相机(行8-10)的多谱室外显像合成的RGB室外图像的结果；

图12图形地示出与Finlayson & Schaefer的方法(F & S)以及Sunshine等人的方法相比的以度为单位的在从多谱图像数据库合成的RGB室外图像上α＝100的照明估计示例的精度与所使用的块片的数量的关系，从左到右是通过(a)Stiles和Burch的色彩匹配函数[15](b)佳能10D相机传感器和(c)尼康D70相机传感器所捕获的仿真的RGB图像的结果；

图13示出与Finlayson和Schaefer的方法(F & S)以及Sunshine等人的方法相比的以度为单位的关于[4]中报告的Mondrian和镜面数据集合的α＝1000的我们的发光体估计示例的精度的比较，精度以度为单位测量；

图14图形地示出对于如[4]报告的具有8比特和16比特动态范围的Mondrian和镜面数据集合的α＝1000的我们的发光体估计示例的精度，我们的方法与Finlayson和Schaefer(F & S)以及Sunshine等人的方法进行比较；

图15图形地示出估计的发光体谱的精度与高斯噪声的标准偏差的关系。垂直轴示出来自对应真实状况(ground truth)的估计的谱的角度偏差，而水平轴示出高斯噪声的标准偏差作为显像的最大亮度的百分比，在左手图像中示出可见谱的性能，而在右手面板上示出与红外谱对应的性能；

图16示出高光移除结果。第一行：在670nm处捕获的初始图像。第二行：通过我们的方法所产生的对应的阴影映射。第三行：通过我们的方法所产生的对应镜面性图像。第四行：通过Ragheb和Hancock的方法[46]所产生的阴影映射；第五行：通过Ragheb和Hancock的方法[46]所产生的镜面性映射。

图17示出我们的方法的高光移除结果与Mallick等人的方法[40]之间的比较。第一列：输入图像。第二列和第三列：由我们的方法得到的漫射和镜面分量。第四列和第五列：Mallick等人的方法所得到的漫射和镜面分量；

图18示出皮肤分段结果。顶行：用于皮肤识别的采样训练图像，具有标记的皮肤区域(具有红色边界)和非皮肤区域(在蓝色边界中)。第二行：在不同照明方向下捕获的对应主体的测试图像。第三行至最后一行是使用不同特征获得的皮肤概率映射。第三行：使用给定真实状况发光体谱估计的反射率获得的。第四行：在估计发光体谱之后使用估计的反射率产生的。第五行：使用通过以真实状况发光体谱对辐射率谱进行归一化所获得的反射率产生的。第六行：使用输入图像的原始辐射率产生的。第七行：使用原始辐射率谱的前20个LDA分量恢复的。

图19示出分别使用用于分类的基于不同输入反射率的特征的具有若干皮肤像素识别方法的精度的表7。

图20图形地示出材料聚类，其中，由不同的灰色阴影来标记每种材料。顶行：在670nm示出的输入图像的带。第二行：通过聚类特征是给定真实状况发光体的估计的反射率而获得的材料聚类映射。第三行：通过使用以估计的发光体谱作为输入所恢复的反射率特征而得到的材料聚类映射。第四行：使用通过以真实状况发光体谱对输入辐射率图像进行归一化所获得的反射率的材料聚类映射。

图21是用于估计示例1的图像的照明谱的方法的流程图；

图22是用于估计示例1的图像的表面反射率、表面阴影和镜面性分量的方法的流程图；

图23是用于选择在估计示例1的图像的照明谱中使用的图像的子图像的方法的流程图；

图24是用于估计示例1的图像的照明谱的方法的流程图；以及

图25是能够实现所描述的方法的示例计算机系统。

具体实施方式

简介

在该示例中，我们使用基于二色模型的优化方法来确定多谱成像中的光度测定的不变性。我们抛出这样的问题：在结构优化设置中恢复发光体L的谱、表面反射率s和阴影g以及镜面k因子。使用多谱成像和图像块片的结构所提供的附加信息，我们恢复场景的二色参数L、g、S和k。为此，我们公式化将研究下的表面的二色误差和平滑性先验(smoothnesspriors)进行组合的目标代价函数。通过以坐标下降(coordinate descent)方式最小化该代价函数来恢复二色参数L、g、S和k。所使用的方法具有以下优点：在本质上具有通用性，容许对平滑性约束的强制，并且以直接方式扩展到三色设置。此外，目标函数相对于在最小化策略的每一交替步骤中待优化的变量的子集是凸的(convex)。这在我们的方法中产生用于每一个迭代中的优化封闭形式解。为了发光体谱L恢复、高光(highlight)移除和皮肤识别的目的，我们示出我们的方法的有效性。我们还将我们的结果与多个替选方案进行比较。

1.示例1

将参照图21至图24的流程图描述该示例。该示例集成了在统一的优化框架中从单个多谱图像恢复发光体L、光度测定的不变量(即材料反射率k)、阴影g和镜面性k因子。该示例将色彩恒定性问题从三色扩展到多谱和超谱显像，并且还带来若干优点。

(i)通过优化与二色模型的数据接近性，所述方法通常可应用于展现漫射和镜面反射的表面。

(ii)对发光体和表面反射率谱的参变形式或先验知识不进行假设。这与其中对发光体和表面反射率的有限线性模型或普通光源的色度进行假设的其它方法形成对照。

(iii)与使用二色模型的其它方法相比，该示例甚至对小数量的不同材料反射率谱也良好地执行。

(iv)与用于三色显像的基于二色平面的方法不同，该示例不需要预先分段的图像作为输入。替代地，自动二色块片选择处理确定要用于发光体估计的均匀反照率块片。

(v)噪声扰乱分析示出我们的照明估计方法比替选方案更精确，并且相对于所使用的表面块片的数量是稳定的。

(vi)在可以将表面几何上的任何正则化矩阵(regulariser)合并到所述方法的意义上，在此提出的优化框架是灵活和通用的。该示例包括用于阴影场的平滑性的强健正则化矩阵的两个实例。

(vii)与US 7443539相比，该示例不需要第二传感器，不需要基准板，可以仅使用软件来完成自动白平衡，从而可以通过更低的代价来提供全部解。

在此，我们解决这样的问题：通过对来自单个多谱图像的照明功率谱L、阴影g和镜面性k的估计来恢复光度测定的不变量(即材料/表面反射率S)。我们提出的方法假设场景是均匀地受照射的。这种假设对于广泛的情况(例如在自然太阳光或远距光源照射场景的情况下)是普遍并且有效的。基于二色反射模型[48]，我们抛出恢复问题作为在结构优化设置中的优化问题。使用多谱成像提供的附加信息以及自动选择的图像块片的结构，我们恢复场景的二色参数L、g、s和k。由于目标函数相对于待优化的每一变量子集是凸的，因此我们可以恢复迭代优化的封闭形式解。我们采用二次表面平滑性误差作为正则化矩阵，并且示出当使用替选正则化矩阵时可以如何获得封闭形式解。

2.反射模型参数的恢复

我们在此提出该示例中采用的目标函数。我们进一步关于用于根据二色反射模型恢复表面几何g和漫射反射率S的在此适用的优化方法进行阐述。

该示例使用基于反射率S在表面块片包括的所有这些像素上是均匀的平滑表面块片的处理的结构方法。如上所述，恢复用于图像的光度测定的参数的处理基于关注于减少二色模型产生的估计与波长索引的谱数据I之间的差异的优化方法。我们通过提供Shafer提出的二色模型的概述[48]而开始。随后，我们公式化相对于[48]中的模型的目标最小化函数，并且基于从具有均匀反射率的平滑图像块片绘制的辐射率结构而公式化优化策略。我们还提出我们的用于选择算法所使用的块片的策略，并且详细描述坐标下降优化过程。该优化策略基于关注于恢复光谱L、表面阴影g和表面反射率s特性从而恢复最佳的二色反射率参数(二色参数)L，S，g和k的交替的步骤。

2.1二色反射模型

我们采用Shafer[48]介绍的二色模型，从而将光谱功率、表面反射率和表面辐射率进行相关。该模型假设穿过被观测场景的空间域的均匀照明。按照该模型，表面辐射率分解为漫射和镜面分量。让其谱是L(λ)的发光体照射具有在像素位置u和波长λ处的谱数据(例如谱辐射率或谱反射率)I(λ，u)的图像中表示的对象。通过该组分，二色模型于是变为：

I(λ，u)＝g(u)L(λ)S(λ，u)+k(u)L(λ) (1)

在公式1中，阴影因子g(u)控制从对象反射的漫射光的部分，并唯一取决于表面几何。注意，对于完全的lambertian表面，我们具有即，表面法线和光方向之间的角度的余弦。另一方面，因子k(u)模拟引起场景中的镜面性的表面不规则性。利用该模型，我们目的在于从图像的谱辐射率I(λ，u)恢复位置u和波长λ处的阴影因子g(u)，镜面系数k(u)，光谱L(λ)和谱反射率S(λ，u)。

2.2目标函数

利用上述二色模型，为了优化的目的，我们进行定义我们的目标函数。该示例将其像素值与索引为波长λ∈{λ₁，...λ_n}的谱辐射率I(λ，u)的测量相对应的多谱图像作为输入。如前所述，我们的目标是使被观测的数据I(λ，u)适合于二色模型，以恢复参数g(u)，k(u)和S(λ，u)。一般来说，在此我们将多谱图像的二色代价函数I看作其二色误差和约束项(这里称作每个图像位置的几何正则化项(regularisation term)R(u))的加权和。

在公式2中，α是用作右手侧的二色误差与正则化项R(u)之间的平衡因子的常数。波长独立的正则化项R(u)与场景几何有关，并且将稍后阐述。

至此，我们将关注聚焦在方程2的解空间上。注意，在没有进一步约束的情况下最小化代价F(I)是个欠确定的问题。这是因为这样的事实：对于具有包含m个像素的n个谱带的图像，我们将必须对2m+n+m×n个变量最小化，而在方程2的和中仅具有m×n项。然而，我们注意到，如果模型应用于以相同材料制成的平滑表面(即反照率在考虑中的块片或图像区域是规则的)，则可以进一步约束该问题。这提出了两个约束。首先，表面上的所有位置共享共同漫射反射率。因此，假设均匀反照率表面P具有对于每一像素相同的反射率u∈P，S(λ_i，u)＝S_P(λ_i)。注意，该约束明显将未知S(λ_i，u)的数量从m×n减少到N×n，其中，N是场景中的表面反照率的数量。此外，块片几何的平滑变化允许我们将公式2中的正则化项R(u)公式化为阴影因子的函数g(u)。简而言之，平滑、均匀的反照率表面块片自然地提供约束，从而在提供几何正则化项R(u)的似真公式化的同时明显地减少未知数的数量。

按照以上原理，我们继续对手边的最小化问题施加约束。对于平滑、规则的表面块片，我们考虑用于表面块片的代价函数，并且写出

$F\left(P\right)\stackrel{\mathrm{\Δ}}{=}\underset{u\∈P}{\mathrm{\Σ}}[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{[I({\mathrm{\λ}}_i,u)-L\left({\mathrm{\λ}}_i\right)(g\left(u\right)S_P\left({\mathrm{\λ}}_i\right)+k\left(u\right))]}^2+\mathrm{\αR}\left(u\right)]$

如上，对于所有u∈P，我们具有S(λ_i，u)＝S_P(λ_i)。此外，块片上的平滑性约束意味着阴影因子g(u)应遍及P中的像素而平滑地变化。可以通过最小化阴影映射(shading map)的梯度大小的变化来有效地公式化该约束。这有效地排除经由正则化项的P的阴影映射中的不连续性。

$R\left(u\right)\stackrel{\mathrm{\Δ}}{=}{\left[\frac{\∂g\left(u\right)}{\∂x\left(u\right)}\right]}^2+{\left[\frac{\∂g\left(u\right)}{\∂y\left(u\right)}\right]}^2---\left(3\right)$

其中，对于像素位置u，变量x(u)和y(u)分别是列坐标和行坐标。

因此，通过使用图像I中的规则块片的集合P，我们可以通过将目标函数最小化为对于以下公式的替选来恢复二色模型参数：

作为对于F(I)的替选。

2.3光谱恢复

2.3.1同质块片选择120

以上我们将二色模型参数的恢复公式化为遍及表面块片集合P的优化过程。在此我们描述用于自动选择用于方程4中的代价函数的最小化的规则表面块片的我们的方法。在此提出的自动块片选择方法允许将我们的方法应用于任意图像。这与仅可应用于预先分段的图像的另外在文献[21、22、54、55]中的其它方法对照。

如下执行我们的块片选择策略。我们首先以类似格状的方式将图像再划分112为相等大小的块片。对于每一块片，我们将二维超平面拟合为块片中的像素的反射率向量或辐射率的谱数据。接下来，我们注意到，在完美地二色块片中，每一像素的波长索引的辐射率向量完美地位于该超平面(即二色平面)中。为了允许噪声效应，我们将二色块片看做包含其辐射率向量与[54]中的单数值分解(SVD)给定的它们的投影偏离的至少t_p像素百分比的块片。我们通过对与二色平面的角度偏差设置上界t_a来完成该操作，其中，t_p和t_a是全局参数。

然而，为了照明谱恢复的目的，并非所有这些块片可以是有用的。这是因为这样的事实：完美地漫射表面不提供关于发光体谱的任何信息。原因在于，用于这种表面的谱辐射率向量是一维的，仅由波长索引的漫射辐射率向量跨越。另一方面，二色模型意味着镜面性具有与发光体相同的谱，其中，镜面系数可以看做仅取决于表面几何的比例因子。

因此，对于二色模型参数的恢复，我们仅通过选择116被看作具有均匀反照率的区域当中具有最高对比度的区域而使用高镜面块片。我们通过在谱域上计算均值辐射率的方差而恢复每一块片的对比度。这些高镜面块片提供用于恢复光谱L的手段。这是因为这样的事实：对于具有均匀反照率的高镜面表面块片，表面漫射辐射率向量和发光体向量跨越辐射率空间中的超平面。这是色彩恒定性的特性，其中，为了光功率谱恢复的目的，方法[24、33、37]已经采用子空间投影。

2.3.2优化过程122

我们现在提出在使用图1呈现的注记的我们的方法中采用的优化过程。在此，我们采用迭代方法，从而发现产生方程4中的代价函数的最小值的变量L、S_P、g_P和k_P 100。在每一迭代，我们在分离的步骤中关于L和三元组{g_P，k_P，S_P}最小化代价函数。

实际上，在此提出的过程是关注于使代价函数最小化的坐标下降方法[6]。在图2所示的算法1的伪代码中2总结了我们的最小化策略的步骤序列。坐标下降方法包括两个交织的最小化步骤。在每一迭代，我们将二色变量索引为迭代号t，并且在交织的步骤中关于变量{g_P，k_P，S_P}、{L}的两个子集而优化客观函数。一旦前者变量可获得，则我们可以获得后者变量的优化值。我们在这两个步骤之间进行迭代，直到达到收敛。

算法通过将未知的光谱L(λ)初始化为无偏的均匀照明谱(见图21，100(c))而开始，如算法1的行1中指示的那样。在两次连续迭代之间，一旦发光体谱就角度而言不改变超过预先设置的全局阈值t_L的量，则终止以标识优化该代价函数102的照明谱参数L的值。在以下两个子部分中，我们示出可以采用以上两个优化步骤，以获得封闭形式的二色参数的优化值。

块片集合表面几何的恢复

在第一步骤中，我们估计在迭代t-1恢复的给定光谱L^t-1的几何和优化表面反射率106。这对应于算法1中的行3-5。注意，在迭代t，我们可以对于每一表面块片P单独地求解未知量 和这是因为，这些变量出现在方程4中的分离项中。因此，减少该步骤，以最小化

$F\left(P\right){|}_{L^{t-1}}=\underset{u\∈P}{\mathrm{\Σ}}[{\left|\right|I\left(u\right)-g\left(u\right)D_P^{t-1}-k\left(u\right)L^{t-1}\left|\right|}^2+\mathrm{\αR}\left(u\right)]---\left(5\right)$

其中，漫射辐射率向量是发光体和表面反射率谱的分量式乘法，||·||表示变元向量的L₂范式。

注意，上述最小化包括2|P|+n个未知量，其中，|P|是块片P中的像素的数量。因此，当表面面积较大时，变得在计算上难以处理。实际上，选择的块片仅需要大得足以从谱数据收集有用的统计数据。此外，如上所述，我们可以通过注意在同一表面中的所有像素处的谱辐射率向量位于漫射辐射率向量和光向量L^t-1跨越的2维子空间来进一步减少未知量的自由度。这是在关于色彩恒定性[21、22、54、55]的著作中已经利用的二色模型的特性。

当所有像素辐射率向量I(u)可获得，我们可以经由单数值分解(SVD)获得子空间Q。将来自该SVD运算的两个基础向量表示为z₁和z₂，因此，子空间是Q＝span(z₁，z₂)。由于因此我们可以将参数化放大到 $D_P^{t-1}=v{z}_1+z_2.$

类似地，光向量L^t-1∈Q也可以分解为L^t-1＝w₁z₁+w₁z₂，其中，w₁和w₂的值是两个未知标量。此外，二色平面假设还意味着，给定光向量L^t-1和表面漫射辐射率向量我们可以将任何像素辐射率I(u)分解为前两个向量的线性组合。换句话说，

$I\left(u\right)=g\left(u\right)D_P^{t-1}+k\left(u\right)L^{t-1}$

$=(g\left(u\right)v+k\left(u\right)w_1)z_1+(g\left(u\right)+k\left(u\right)w_2)z_2---\left(6\right)$

在已经获得基础向量z₁、z₂的情况下，我们可以计算像素辐射率I(u)到子空间Q的映射。通过投影而相对于该基础来完成该操作，从而获得标量τ₁(u)，τ₂(u)(见图21，100(b))，从而

I(u)＝τ₁(u)z₁+τ₂(u)z₂ (7)

此外，通过使得公式6和7的右手边相等，我们获得

$g\left(u\right)=\frac{w_2{\mathrm{\τ}}_1\left(u\right)-w_1{\mathrm{\τ}}_2\left(u\right)}{w_{2}v-w_1}---\left(8\right)$

$k\left(u\right)=\frac{{\mathrm{\τ}}_2\left(u\right)v-{\mathrm{\τ}}_1\left(u\right)}{w_{2}v-w_1}---\left(9\right)$

从公式8和9，我们注意到，g(u)和k(u)是v的单变有理函数。

此外，是关于v的线性函数。我们还观察到，项R(u)仅取决于g(u)。因此，方程5中的客观函数可以简化为v的单变有理函数。因此，将公式8和9代入方程5的右手边的第一项和第二项，我们具有

${F\left(P\right)|}_{L^{t-1}}=\underset{u\∈P}{\mathrm{\Σ}}{\left|\right|I\left(u\right)-\frac{w_2{\mathrm{\τ}}_1\left(u\right)-w_1{\mathrm{\τ}}_2\left(u\right)}{w_{2}v-w_1}(v{z}_1+z_2)-\frac{{\mathrm{\τ}}_2\left(u\right)v-{\mathrm{\τ}}_1\left(u\right)}{w_{2}v-w_1}{L}^{t-1}\left|\right|}^2$

$+\underset{u\∈P}{\mathrm{\Σ}}\frac{\mathrm{\α}}{{(w_{2}v-w_1)}^2}[{\left(\frac{\∂m\left(u\right)}{\∂x\left(u\right)}\right)}^2+{\left(\frac{\∂m\left(u\right)}{\∂y\left(u\right)}\right)}^2]$

$=\underset{u\∈P}{\mathrm{\Σ}}\frac{1}{{(w_{2}v-w_1)}^2}\left|\right|I\left(u\right)w_2-(w_2{\mathrm{\τ}}_1\left(u\right)-w_1{\mathrm{\τ}}_2\left(u\right))z_1-{\mathrm{\τ}}_2\left(u\right)L^{t-1})v$

$-(I\left(u\right)w_1-(w_2{\mathrm{\τ}}_1\left(u\right)-w_1{\mathrm{\τ}}_2\left(u\right))z_2-{\mathrm{\τ}}_1\left(u\right)L^{t-1}{\left|\right|}^2$

$+\frac{\mathrm{\α}}{{(w_{2}v-w_1)}^2}\underset{u\∈P}{\mathrm{\Σ}}[{\left(\frac{\∂m\left(u\right)}{\∂x\left(u\right)}\right)}^2+{\left(\frac{\∂m\left(u\right)}{\∂y\left(u\right)}\right)}^2]$

$=\underset{u\∈P}{\mathrm{\Σ}}{\left|\right|\frac{p\left(u\right)v-q\left(u\right)}{w_{2}v-w_1}\left|\right|}^2+\frac{\mathrm{\αN}}{{(w_{2}v-w_1)}^2}$

$=\underset{u\∈P}{\mathrm{\Σ}}{\left|\right|\frac{p\left(u\right)}{w_2}+\frac{\frac{w_1}{w_2}p\left(u\right)-q\left(u\right)}{w_{2}v-w_1}\left|\right|}^2+\frac{\mathrm{\αN}}{{(w_{2}v-w_1)}^2}$

$=\underset{u\&Element;P}{\mathrm{\&Sigma;}}\frac{{\left|\right|p\left(u\right)\left|\right|}^2}{w_2^2}+\frac{2}{w_{2}v-w_1}{\left|\right|}^2\underset{u\&Element;P}{\mathrm{\&Sigma;}}<\frac{p\left(u\right)}{w_2},\frac{w_1}{w_2}p\left(u\right)-q\left(u\right)>$

$+\frac{1}{{(w_{2}v-w_1)}^2}(\underset{u\&Element;P}{\mathrm{\&Sigma;}}{\left|\right|\frac{w_1}{w_2}p\left(u\right)-q\left(u\right)\left|\right|}^2+\mathrm{\&alpha;N})---\left(10\right)$

其中，<.，.>表示两个向量的内积，并且

m(u)＝w₂τ₁(u)-w₁τ₂(u)

p(u)＝I(u)w₂-(w₂τ₁(u)-w₁τ₂(u))z₁-τ₂(u)L^t-1

q(u)＝I(u)w₁-(w₂τ₁(u)-w₁τ₂(u))z₂-τ₁(u)L^t-1

$N=\underset{u\&Element;P}{\mathrm{\&Sigma;}}[{\left(\frac{\&PartialD;m\left(u\right)}{\&PartialD;x\left(u\right)}\right)}^2+{\left(\frac{\&PartialD;m\left(u\right)}{\&PartialD;y\left(u\right)}\right)}^2]$

注意，给定L^t-1，p(u)、q(u)、w₁和w₁是已知的。通过变量的改变，我们可以将方程10的右手边写为r的二次函数，r的最小值根据以下公式而得到：

$r^{*}=-\frac{{\mathrm{\&Sigma;}}_{u\&Element;P}<\frac{p\left(u\right)}{w_2},\frac{w_1}{w_2}p\left(u\right)-q\left(u\right)>}{{\mathrm{\&Sigma;}}_{u\&Element;P}{\left|\right|\frac{w_1}{w_2}p\left(u\right)-q\left(u\right)\left|\right|}^2+\mathrm{\&alpha;N}}---\left(11\right)$

这给定对应的最小者因此，给定发光体谱L^t-1，我们可以通过将v的优化值代入公式8和9而恢复g_P、k_P。

漫射辐射率分量计算为在波长λ处的谱反射率由 $S_P^t\left(\mathrm{\&lambda;}\right)=\frac{D_P^t\left(\mathrm{\&lambda;}\right)}{L^{t-1}\left(\mathrm{\&lambda;}\right)}$ 给定。

在该示例中，由于每一材料的二色模型暗示了谱辐射率或RGB图像色彩可以表示为2维超平面，因此我们使用2维超平面。替代地，可以使用可以包括与附加维度对应的来回反射或环境光的项的更复杂的模型。

发光体谱L的恢复

在每一迭代t的第二步骤中，我们求解给定和的L^t和由于方程4中的第二项R(u)是波长独立的，因此算法1的行6中的优化问题可以简化为最小化

$F^{*}\left(I\right){|}_{g^t,k^t}=\underset{P\&Element;P}{\mathrm{\&Sigma;}}\underset{u\&Element;P}{\mathrm{\&Sigma;}}{\left|\right|I\left(u\right)-g^t\left(u\right)D_P-k^t\left(u\right)L\left|\right|}^2$

$=\underset{P\&Element;P}{\mathrm{\&Sigma;}}\underset{u\&Element;P}{\mathrm{\&Sigma;}}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(I({\mathrm{\&lambda;}}_i,u)-g^t\left(u\right)D_P\left({\mathrm{\&lambda;}}_i\right)-k^t\left(u\right)L\left({\mathrm{\&lambda;}}_i\right))}^2---\left(12\right)$

其中，D_P＝L·S_P

由于客观函数12是二次的，并且因此相对于L和D_P是凸的，因此可以通过使得的各个偏微分等于零而获得这些变量的优化值。这些偏微分由以下公式给出：

$\frac{{\&PartialD;F}^{*}\left(I\right){|}_{g^t,k^t}}{\&PartialD;L\left({\mathrm{\&lambda;}}_i\right)}=-2\underset{P\&Element;P}{\mathrm{\&Sigma;}}\underset{u\&Element;P}{\mathrm{\&Sigma;}}(I({\mathrm{\&lambda;}}_i,u)-g^t\left(u\right)D_P\left({\mathrm{\&lambda;}}_i\right)-k^t\left(u\right)L\left({\mathrm{\&lambda;}}_i\right))k^t\left(u\right)$

$\frac{{\&PartialD;F}^{*}\left(I\right){|}_{g^t,k^t}}{\&PartialD;D_P\left({\mathrm{\&lambda;}}_i\right)}=-2\underset{u\&Element;P}{\mathrm{\&Sigma;}}(I({\mathrm{\&lambda;}}_i,u)-g^t\left(u\right)D_P\left({\mathrm{\&lambda;}}_i\right)-k^t\left(u\right)L\left({\mathrm{\&lambda;}}_i\right))g^t\left(u\right)$

使得以上两个公式等于零，我们获得

$L\left({\mathrm{\&lambda;}}_i\right)=\frac{{\mathrm{\&Sigma;}}_{P\&Element;P}{\mathrm{\&Sigma;}}_{u\&Element;P}[k^t\left(u\right)I({\mathrm{\&lambda;}}_i,u)-g^t\left(u\right)k^t\left(u\right)D_P\left({\mathrm{\&lambda;}}_i\right)]}{{\mathrm{\&Sigma;}}_{P\&Element;P}{\mathrm{\&Sigma;}}_{u\&Element;P}{\left(k^t\left(u\right)\right)}^2}---\left(13\right)$

$D_P\left({\mathrm{\&lambda;}}_i\right)=\frac{{\mathrm{\&Sigma;}}_{u\&Element;P}[g^t\left(u\right)I({\mathrm{\&lambda;}}_i,u)-g^t\left(u\right)k^t\left(u\right)L\left({\mathrm{\&lambda;}}_i\right)]}{{\mathrm{\&Sigma;}}_{u\&Element;P}{\left(g^t\left(u\right)\right)}^2}---\left(14\right)$

从公式13和14，发光体谱可以以封闭形式求解为：

$L^{*}\left({\mathrm{\&lambda;}}_i\right)=\frac{{\mathrm{\&Sigma;}}_{P\&Element;P}{\mathrm{\&Sigma;}}_{u\&Element;P}{k}^t\left(u\right)I({\mathrm{\&lambda;}}_i,u)-{\mathrm{\&Sigma;}}_{P\&Element;P}\left[\frac{\left({\mathrm{\&Sigma;}}_{u\&Element;P}{g}^t\left(u\right)k^t\left(u\right)\right)\left({\mathrm{\&Sigma;}}_{u\&Element;P}{g}^t\left(u\right)I({\mathrm{\&lambda;}}_i,u)\right)}{{\mathrm{\&Sigma;}}_{u\&Element;P}{\left(g^t\left(u\right)\right)}^2}\right]}{{\mathrm{\&Sigma;}}_{P\&Element;P}{\mathrm{\&Sigma;}}_{u\&Element;P}{\left(k^t\left(u\right)\right)}^2-{\mathrm{\&Sigma;}}_{P\&Element;P}\left[\frac{\left({\mathrm{\&Sigma;}}_{u\&Element;P}{g}^t\left(u\right)k^t\left(u\right){)}^2\right)}{{\mathrm{\&Sigma;}}_{u\&Element;P}{\left(g^t\left(u\right)\right)}^2}\right]}---\left(15\right)$

2.4阴影、反射率和镜面性恢复

注意，在以上的优化方案中，我们恢复用于发光体谱的恢复的每一块片P∈P中的像素的几何变量。这意味着，尽管我们仅已经计算了对于像素点u∈P的变量g(u)、k(u)和S(.，u)，但我们也已经恢复了发光体谱L。由于L是场景中的全局光度测定的变量，因此我们可以通过直接方式使用L来恢复其余二色变量。其余二色变量包括在每一像素处的阴影g、反射率s和镜面性k。

为此，我们假设输入场景包括具有缓慢地改变的反射率的平滑表面。换句话说，每一像素的邻域可以被看作以相同材料制成的局部平滑块片，即，邻域中的所有像素共享相同的谱反射率。给定发光体谱L，我们可以通过应用与图2的算法1中的行4对应的过程而获得在关注的像素处的邻域的阴影g、镜面性k和反射率s。这对应于在以上部分中提出的优化方法中使用的两个步骤中的第一步骤的应用。

在图3所示的算法2中总结该方法的伪代码。注意，具有缓慢地变化的反射率的平滑表面的假设可应用于一大类表面具有低文本程度或穿过边缘或附近封闭(nearocclusion)的场景。按照该假设，在每一像素处的反射率恢复为其周围邻域的共享的反射率106。为了估计阴影和镜面性108，我们可以应用它们的闭合形式公式，如公式8和9所示。这些公式在理想条件下得到精确解，这要求所有像素辐射率向量位于发光体谱L和漫射辐射率向量D跨越的同一二色超平面中。

然而，实际上，对于多谱图像来说，包含破坏这种假设并且数值上不稳定地呈现以上商表达式的噪声是普遍的。因此，为了强制越过像素的阴影因子的平滑变化，在获得谱反射率106之后，我们重新计算108阴影g和镜面性k系数。这是因为反射率谱一般比其它两个变量(即阴影g和镜面性k因子)更稳定的观测。具体地说，我们可以计算阴影g和镜面k系数作为源自像素辐射率I(u)在发光体谱L和漫射辐射率谱向量D跨越的子空间上的投影的系数。

当违反二色超平面假设(即发光体谱L对于材料的漫射辐射率谱D是共线的)时，这种框架中止。这呈现块片像素的辐射率谱跨越的子空间，以收缩(coppapse)到1维空间。因此，这些辐射率谱的单数值分解(SVD)没有成功找到子空间的两个基础向量。由于漫射分量D是发光体功率谱和材料反射率的乘积，因此该条件意味着材料具有均匀谱反射率。换句话说，故障情况仅当输入场景包含具有均匀反射率的单个材料(即类似灰色的阴影的材料)时发生。

这种故障情况实际上很稀少。实际上，当场景包含多于一个的材料时，随着从场景采样更多均匀反照率块片，存在更多的机会以在发光体谱L与表面漫射辐射率D谱之间引入非共线性。简而言之，我们的方法确保恢复具有多于一个独特反照率的场景上的二色模型参数L、g、s和k。

3施加平滑性约束(见图21，100(a))

在此，我们示出可以如何对优化处理施加平滑性约束。

以上我们解决了使用方程2中的正则化项R(u)强制阴影场g＝{g(u)}_u∈I上的平滑性约束的需要。在方程3中，我们提出鼓励阴影场的缓慢空间变化的正则化矩阵。存在用于在先前部分中介绍的优化框架中使用该正则化矩阵的两个原因。首先，给定发光体谱，正则化矩阵产生用于表面几何和反射率的闭合形式解。其次，正则化矩阵是在来自阴影方法的形状上施加的平滑性约束的回忆，并且因此，提供了文献中的(例如[58]中的)其它方法与上述优化方法之间的关联。然而，我们需要强调，以上优化过程绝不意味着框架不可应用于替选正则化矩阵。实际上，在可以根据手边的表面采用其它正则化函数的意义上，我们的目标函数是灵活的。

我们在阴影场上引入对于噪声和轮廓强健的并且适用于表面几何变化的多个替选正则化矩阵。为此，我们通过引入强健的正则化矩阵而开始。我们于是基于表面曲率和形状索引而提出扩展。

为了量化阴影的平滑性，一个选项是将阴影场的梯度处理为平滑性误差。在方程3中，我们引入平滑性的二次误差函数。然而，在特定情况下，强制方程2中引入的二次正则化矩阵导致过度平滑表面的不期望效果，参见来自阴影的形状的场中的[11、31]。

作为替选方法，我们利用来源于强健统计的场的内核函数。正规来说，强健内核函数ρ_σ(η)量化与余差(residual)η及其影响函数关联的能量，即测量对于阴影场的改变的灵敏度。每一余差进而被分配影响函数Γ_σ(η)定义的权重。因此，能量与如的影响函数的第一力矩(first moment)有关。图4的表1示出用于Tukey双权重[26]、Li的自适应势函数[38]和Huber的M估计器[30]的公式。

3.1强健阴影平滑性约束

在已经引入上述强健估计器的情况下，我们进行采用它们作为用于目标函数的正则化矩阵。在此，存在若干可能性。其中之一是通过定义关于阴影梯度的强健正则化矩阵而直接最小化几何变化。在此情况下，正则化矩阵R(u)由以下公式给定：

$R\left(u\right)={\mathrm{\&rho;}}_{\mathrm{\&sigma;}}\left(\right|\frac{\&PartialD;g}{\&PartialD;x}\left|\right)+{\mathrm{\&rho;}}_{\mathrm{\&sigma;}}\left(\right|\frac{\&PartialD;g}{\&PartialD;y}\left|\right)---\left(16\right)$

无论是否有效，以上公式都仍然采用阴影场的梯度作为平滑性的测度。

3.2曲率一致性

在此我们探索将曲率用作一致性的度量。我们可以替代地考虑由其曲率给定的手边的表面的固有特性。具体地说，Ferrie和Lagarde[17]已经使用主曲率的全局一致性来细化在来自阴影的形状中的表面估计。此外，确保曲率方向的一致性不一定意味着用于方位和深度的不连续性的大处罚。因此，这种度量可以避免过度平滑，这是二次平滑性误差的缺点。

可以通过将曲率一致性处理为副本(manifold)而在阴影场上定义曲率一致性。为了开始，我们使用Hessian矩阵来定义阴影场的结构

$H=\left(\begin{array}{cc}\frac{{\&PartialD;}^{2}g}{\&PartialD;x^2}& \frac{{\&PartialD;}^{2}g}{\&PartialD;x\&PartialD;y}\\ \frac{{\&PartialD;}^{2}g}{\&PartialD;x\&PartialD;y}& \frac{{\&PartialD;}^{2}g}{\&PartialD;y^2}\end{array}\right)$

副本的主曲率因此被定义为Hessian矩阵的特征值。设这些特征值由k₁和k₂表示，其中，k₁≥k₂。此外，我们可以使用主曲率来描述使用如下定义的形状索引[34]的局部拓扑

$\mathrm{\&phi;}=\frac{2}{\mathrm{\&pi;}}\arctan \left(\frac{k_1+k_2}{k_1-k_2}\right)---\left(17\right)$

以上观测是重要的，因为其允许将阴影场的平滑处理变为加权均值处理，其中，分配给像素的权重由关于本地邻域的局部拓扑的相似性(即形状索引)确定。有效地，想法是支持属于与关注的像素相同或相似的形状类的邻域中的像素。这是对于方程3中定义的二次平滑性项的改进，因为其避免了穿过不连续性的阴影因子的不加区分的平均化。也就是说，其由清晰度边缘保存。

对于每一像素u，我们考虑u周围的本地邻域N，并且给邻域中的每一像素u^＊分配权重为 $w\left(u^{*}\right)=\exp (-\frac{{(\mathrm{\&phi;}\left(u^{*}\right)-{\mathrm{\&mu;}}_{\mathrm{\&phi;}}\left(N\right))}^2}{2{\mathrm{\&sigma;}}_{\mathrm{\&phi;}}^2\left(N\right)}),$ 其中，μ_φ(N)和σ_φ(N)是在邻域N上的形状索引的均值和标准偏差。使用该加权处理，我们获得自适应加权平均正则化矩阵如下：

$R\left(u\right)={(g\left(u\right)-\frac{{\mathrm{\&Sigma;}}_{u^{*}\&Element;N}w\left(u^{*}\right)g\left(u^{*}\right)}{{\mathrm{\&Sigma;}}_{u^{*}\&Element;N}w\left(u^{*}\right)})}^2---\left(18\right)$

该方法可以被看作方程16中提出的具有固定内核的强健正则化矩阵函数的扩展。为了调节应用于邻域的平滑化的级别，我们考虑形状索引统计[34]，从而自适应地改变强健内核的宽度。自适应内核宽度背后的原理是，具有形状索引的大变化的邻域需要比具有较平滑的变化的邻域更强的平滑化。除了如以下定义的像素型的内核宽度之外，正则化矩阵函数确切地与方程16相同

其中，N是像素u周围的邻域，是N的基数，K_φ是归一化项。

通过内核宽度的以上公式，可以观察到，邻域内的形状索引的明显变化与小内核宽度对应，使得强健的正则化矩阵产生严重平滑化。相对而言，当形状索引变化较小时，由于较宽的内核宽度，产生较低的平滑化等级。

注意，使用该部分中较早引入的强健正则化矩阵作为二次正则化矩阵的替选不排除在部分2.3.2中描述的优化框架的应用性。实际上，正则化矩阵的改变仅影响方程10中的目标函数的公式化，其中，阴影因子g(u)可以表示为如方程8中给定的单变函数。由于所有以上强健正则化矩阵仅取决于阴影因子，因此所得目标函数仍然是变量的函数。此外，通过强健正则化矩阵的线性化，我们可以仍然在数值上将正则化矩阵表示为r变量的二次函数。随后，先前提出的闭合形式解表示得如初始地描述的那样。

4对于三色显像的自适应-示例2

我们在此提供本质上是超谱的以上示例1与该示例的三色显像之间的联系。也就是说，我们示出如何利用以上优化方法来从三色图像恢复二色参数。

我们将用于多谱图像的二色模型变换为用于三色显像的模型。让我们通过C_c(λ)表示三色传感器c的谱灵敏度函数(其中，c∈{R，G，B})。传感器c对于到达位置u处的谱辐射照度(irradiance)的响应由I_c(u)＝∫_ΩE(λ，u)C_c(λ)dλ给定，其中，E(λ，u)是图像辐射照度，Ω是进入光(incoming light)的谱。此外，已知图像辐射照度与场景辐射率I(λ，u)成正比，即E(λ，u)＝K_opt cos⁴β(u)I(λ，u)，其中，β(u)是进入光纤维(rayon)透镜的入射角，K_opt是仅取决于透镜的光学的常数[28]。因此，我们具有

$I_c\left(u\right)=K_{\mathrm{opt}}{\cos }^4\mathrm{\&beta;}\left(u\right)\underset{\mathrm{\&Omega;}}{\&Integral;}I(\mathrm{\&lambda;},u)C_c\left(\mathrm{\&lambda;}\right)\mathrm{d\&lambda;}$

$=K_{\mathrm{opt}}{\cos }^4\mathrm{\&beta;}\left(u\right)\underset{\mathrm{\&Omega;}}{\&Integral;}(g\left(u\right)L\left(\mathrm{\&lambda;}\right)S(\mathrm{\&lambda;},u)+k\left(u\right)L\left(\mathrm{\&lambda;}\right))C_c\left(\mathrm{\&lambda;}\right)\mathrm{d\&lambda;}$

$=K_{\mathrm{opt}}{\cos }^4\mathrm{\&beta;}\left(u\right)\underset{\mathrm{\&Omega;}}{\&Integral;}L\left(\mathrm{\&lambda;}\right)S(\mathrm{\&lambda;},u)C_c\left(\mathrm{\&lambda;}\right)\mathrm{d\&lambda;}+K_{\mathrm{opt}}{\cos }^4\mathrm{\&beta;}\left(u\right)k\left(u\right)\underset{\mathrm{\&Omega;}}{\&Integral;}L\left(\mathrm{\&lambda;}\right)C_c\left(\mathrm{\&lambda;}\right)\mathrm{d\&lambda;}$

$=g^{*}\left(u\right)D_c\left(u\right)+k^{*}\left(u\right)L_c$

其中，g^*(u)＝K_opt cos⁴β(u)g(u)并且k^*(u)＝K_opt cos⁴β(u)k(u)。

在此，我们注意到，D_c(u)＝∫_ΩL(λ)S(λ，u)Cc(λ)dλ和L_c＝∫_ΩL(λ)C_c(λ)dλ分别是与位置u对应的表面漫射色彩以及发光体色彩的c分量。

用于场景的三色图像I的二色代价函数公式化为

其中，R(u)是方程2中描述的空间变化正则化项。

方程20中的代价函数是方程2的特殊情况，其中，n＝3。因此，如上所述的恢复二色参数的方法可以应用于该情况，以恢复三色漫射色彩D(u)＝[D_R(u)，D_G(u)，D_B(u)]^T和发光体色彩L＝[L_R，L_G，L_B]²以及阴影和镜面因子g(u)和k(u)直至乘法因子。

5实现

为了恢复发光体谱、去除高光和皮肤分段124的目的，我们在此示出方法的实用性。

可以在图25所示的计算机系统上实现。在该计算机系统包括传感器200和计算机202。

在该示例中，传感器200是超谱相机，其能够捕获场景204的图像，并且将数字图像提供给计算机202的输入端口206。相机可以具有多个带，从而平衡计算代价与精度。相机可以具有少至五个带而高至几百个带。处理器将接收到的图像存储在本地存储器208中。可以如WO2009/152583中描述的紧凑表示形式来在该方法中存储并使用图像。处理器210能够对接收到的图像执行图21至图24所示的方法。这可以包括对先前存储和分析的超谱图像的数据库212的查询。在该示例中，数据库212对于计算机212是本地的，但可替代地可以对于计算机200是远程的。在另一替选方法中，计算机202可以合并到相机200中。

图21至图24的方法中的任何一个或多个的输出是

(i)图像的照明谱

(ii)图像中的点的表面反射率、表面阴影和镜面性，

(iii)选择的子图像的标识；

(iv)图像的修改版本；或

(v)材料标识指示，

然后作为输出对于被提供给输出端口212。该输出端口212可以连接到输出设备，例如显示设备或打印机。可替代地，输出端口212可以连接到能够对接收到的输出执行其它方法的本地或远程的另一计算机系统，或简单地连接到存储设备。

在多个图像数据库上实现该方法，从而验证恢复的二色参数的精度。我们的数据集合包括室内和室外、具有均匀和散乱背景的多谱和RGB图像、在自然和人工光照条件下。为此，我们已经获取户内的可见和近红外谱中捕获的两个多谱图像数据库。它们包括在人工光源下拍摄的室内图像和在自然日光和天空光条件下的室外图像。从这两个数据库，对于佳能10D和尼康D70相机传感器的谱灵敏度函数以及CIE标准RGB色彩匹配函数合成两个三色图像数据库[15]。除了这些数据库，我们也已经对我们的算法的性能与[3]中的Barnard等人报告的基准数据库的替选进行了比较。

室内数据库包括在具有变化的方向和谱功率的10个方向光源中的一个下捕获的51个人类主体的图像。光源划分为两行。它们中的第一行放置在相机系统上，第二行在与相机相同的高度处。调整光的主方向，从而指向场景的中心。已经使用一对OKSI Turnkey超谱相机获取显像。这些相机装配有液晶可调谐滤波器，这些滤波器允许在可见(430-720nm)和近红外(650-990nm)波长范围中分辨多谱图像达到10nmin。为了获得每一图像的真实状况发光体谱，我们已经测量从白校准目标(即考虑中的光源所照射的LabSphere Spectralon)反射的平均辐射率。使用相同的相机系统和校准目标，我们已经从来自一天中的不同时间的七个不同观看角度的四个不同视点捕获围场(paddock)的室外图像。

5.1照明谱恢复

对于我们关于照明谱恢复的试验，我们对我们的方法得到的结果与Finlayson和Schaefer[21]提出的色彩恒定性方法交付的结果进行比较。注意，[21]中的方法可以被看作方程2(见图21，100(c))中的代价函数的最小二乘解，其中，基于二色模型来估计发光体色彩，无需关于发光体统计的先前假设。虽然对三色显像执行他们的试验，但该方法可以以直接的方式而适用于多谱数据。他们的方法依赖于二色平面假设。也就是说二色模型意味着具有同质反射率的块片中的像素的二维色彩空间。使用该构思，照明估计转化为优化问题，从而最大化所有二色平面上的光色彩向量的总投影长度。几何学上，该方法预测发光体色彩作为二色平面的交点，当二色平面之间的角度较小时，这可能导致数值上不稳定的解。

Finlayson和Schaefer的方法可以如下适用于多谱图像。首先，我们采用我们的自动块片选择方法以提供同质块片作为用于他们的色彩恒定性算法的输入。其次，我们求解二色平面上的投影矩阵的和的特征系统。光色彩向量是与最大特征值对应的特征向量。

在此使用的另外替选方法类似于Sunshine等人提出的谱解卷积方法[52]以将吸收带特性恢复到表面材料化学。该方法利用反射率谱的上界包络，也称为反射率谱的闭连集(continuum)，其可以被看作没有任何吸收特征的反射率谱。对于发光体恢复，我们将估计的发光体谱看做在所有像素处的辐射率谱的闭连集。[52]中的著作假设闭连集按照对数反射率的比例是波数(即波长的倒数)的线性函数。利用该假设，于是将该参变形式拟合为辐射谱的闭连集，以恢复发光体。注意，所得到的发光体不依赖于块片选择，并且因此独立于块片的数量。

我们通过提供关于超谱显像的结果而开始。然后我们将我们的关注转向三色显像中的光色彩恢复。我们通过提供噪声扰动分析来结束这部分。

5.1.1多谱光谱恢复

如上所述，我们首先将我们的关注聚焦于将我们的二色参数恢复算法用于超谱显像中的发光体谱估计。为此，我们已经使用1、5、10、20、30、40和50个常规反照率的自动选择的块片来进行试验。每一块片具有20x20像素的大小。光谱恢复的精度被测量为n维中的估计的谱与真实状况谱之间的欧几里得偏差角度，其中，n是采样的波长的数量。然后这些结果与在相同数量的块片上通过Finlayson和Schaefer[21]以及Sunshine等人[52]的方法获得的结果进行比较。

图5中的表2示出在可见和红外谱范围中、室内面部数据库中的所有图像与选择的块片的数量的以度为单位的角度误差的均值和标准偏差。图6中绘制相似的统计，由误差条的中点和长度表示的角度误差的均值和标准偏差。也就是说，图6示出与Finlayson和Schaefer的方法(F & S)以及Sunshine等人的方法相比、以度为单位的多谱面部图像数据库上的精度相对我们的照明估计方法中使用的块片的数量。示出对于可见(左边)和近红外范围(右边)的结果。再次，注意，Sunshine等人的方法[52]与选择的块片的数量无关。

通过分配给方程2中的正则化项的权重α＝100000来报告结果。在该实验中，正则化项定义为阴影变化的平滑性，如方程3所示。为了获得α的优化值，我们执行与交叉确认中采用的栅格搜索相似的过程。该过程包括：对于不同参数值在数据库的随机采样的部分上若干次应用我们的算法，然后选择产生最高总体性能的值。

如图5的表2所示，当对于两个谱范围使用不超过20个同质块片时，我们的算法实现比替选方法更高的精度。表2示出与Finlayson和Schaefer的方法(F & S)以及Sunshine等人的方法相比、以度为单位的精度相对在可见和红外谱下捕获的多谱面部图像数据库上用于我们的发光体估计方法的块片的数量。

可注意到，甚至对于单个块片，我们的方法仍然明显优于Finlayson和Schaefer的方法[21]。该观查结果与如下的已知事实一致：二色平面方法及其变化需要至少两个同质表面以计算二色平面之间的交集。此外，估计的发光体的角度误差随着块片的数量从1增加到5而减少。然而，随着块片的数量增大到超过5，我们的方法的角度误差趋向于增加。虽然我们的方法在可见范围中保持比Finlayson和Schaefer的方法更精确，但在近红外范围中，当多于20个块片的情况下，我们的方法的精度稍微低于Finlayson和Schaefer的方法。然而，我们的方法对于小数量的同质块片能够实现合理的估计。最后，我们可以下结论：Sunshine等人的方法[52]通常要差于其它两个。

该示例对于增加的块片数量提供了精度改进的余量。或者该示例在计算红外光谱中优于替选方法。同时，其性能在可见范围中可与Finlayson和Schaefer[21]比较，并且好于Sunshine等人的方法[52]。此外，甚至在使用单个同质块片用于光估计时，我们的方法比其它方法更强健。

在图7和图8的表3中，在我们的数据集合中我们示出照明室外场景的自然光的恢复的谱的精度。在此，在正则化权重α＝10000000的情况想应用我们的算法。在此情况下，当使用多于20个随机选择的块片时，我们的方法比Finlayson和Schaefer的[21]给予了略微的改进。在5-15块片之间使用时，两种算法跨越整个数据库而实现相似的精度和误差方差。如上，在单个均匀反照率块片的情况下，我们的算法明显优于替选方法。图8还示出我们的方法相对于选择的块片的数量的增加的稳定性。还可注意到，由于越过场景的谱辐射率的较宽的变化，用于室外图像数据库的所有算法的精度低于用于面部数据库的所有算法的精度。

5.1.2三色光恢复

接下来，为了来自三色图像的发光体色彩估计的目的，我们将我们的关注转向我们的参数恢复方法的实用性。为此，我们从上述多谱面部和室外数据库生成RGB显像。通过仿真多个三色谱响应(包括CIE-1932色彩匹配函数[15]和用于尼康D70和佳能10D的相机传感器)来合成它们。此外，我们将我们的方法以及替选方法应用于Barnard等人描述的Mondrian和镜面图像数据集合[4]。我们还对我们的算法的性能与相同作者的[3]中描述的若干色彩恒定性方法进行比较。

为了示出改变α的值的效果，我们在α＝1000以及α＝100的情况下进行试验。在图9和图10的表4中，我们示出在α＝10000以及20x20像素的块片大小的情况下RGB面部图像上的光估计精度的结果。当块片的数量(或场景中的可用固有色彩的数量)是5或更少时，我们的方法在估计精度和稳定性方面优于替选方法。另一普遍趋势是，随着选择的块片的数量增加，我们的方法和Finlayson和Schaefer的[21]方法提高它们的精度。然而，当我们使用20个或更多的块片时，这种改进对于我们的方法是微小的。同时，当Finlayson和Schaefer的方法[21]应用于来自对于佳能10D和尼康D70传感器合成的图像的明显大数量的块片时，Finlayson和Schaefer的方法[21]趋向于实现更接近于真实状况的估计。有趣的是，对于模仿人类视觉感知的色彩匹配函数的仿真的图像，我们的方法在更稳定(具有更低的角度误差的方差)的同时，跨越所有数量的选择的块片实现与Finlayson和Schaefer的[21]相似精度。在所有我们的试验中，Sunshine等人的方法[52]是给予最差性能的方法。

图11和图12的表5示出室外RGB图像数据库上的发光体色彩估计的精度，其中α＝100以及块片大小为20x20像素。这些统计的主要趋势是，我们的方法实现明显高于其它方法实现的精度。性能的差异在于角度误差的若干标准偏差的量级。虽然随着块片的数量增加(大约20)，我们的方法的性能略微降低，但估计的稳定性对于佳能10D和尼康D70图像保持恒定，并且甚至改进针对CIE1932标准色彩匹配函数仿真的图像。还重要的是，我们的方法对使用CIE1932标准合成的显像执行得更好。这与对于ROB面部图像报告的结果一致。总之，当应用于具有色彩和纹理的高度变化的散乱场景时，我们的方法显现为优于替选方法。

接下来，我们将我们的关注转向[4]中报告的Mondrian和镜面图像数据集合上的发光体估计精度。为了考虑一些图像中的纹理密度的级别，我们选取10x10像素的块片大小，该大小足够小使得每一块片上的反照率均匀的假设仍然成立。在该实验中，如果块片像素的75％或更多从它们在块片的二色平面上的投影偏离小于1度，则块片被看作是同质反射率。我们还执行排除选择包含多于一个固有色彩的高对比块片以及包含由于高光的饱和像素的块片的准则。具体地说，我们在图像中按它们的对比度等级对块片进行排序，并且选择最高对比的块片，排除Mondriant数据集合的每一图像中的前10％的块片。对于镜面数据集合，我们排除每一图像中的块片的前30％，以提供更高级别的色彩饱和度。

在图13和图14的表6中，我们示出当使用1、5、10、20、30、40和50个块片时的精度。

我们的结果与先前试验一致，其示出在具有8比特标准动态范围和16比特扩展的动态范围的两个数据集合上我们的方法优于Finlayson和Schaefer的方法。这不仅由我们的方法产生较低均值的角度误差而且还由其性能的较低标准偏差反映出来。此外，我们的方法提供对于5或更多个选择的块片几乎恒定的角度误差的方差。此外，当在16比特数据集合而不是8比特数据集合上运行时，我们的方法的性能对于更多选择的块片有略微改进。在某种程度上，我们的估计者显现为对于输入图像的动态范围不敏感。这示出我们的方法对于场景反射率的变化更稳定和强健。与Barnard等人报告的基准方法[3]比较，我们的方法排名第二，仅低于[23]中提出的全域映射(gamut mapping)方法。[23]中提出的方法提供用于16比特图像的5∶6-7∶1度与用于8比特图像的6∶3-8∶3度之间的精度，如[3]中的表II所示。然而，对于已经经历包括分段、缩放和修剪操作的若干处理步骤的显像而报告这些结果。另一方面，我们的方法不需要任何前处理，此外，该方法能够恢复可应用于以三色作为特殊情况的超谱显像的所有二色模型参数。

5.1.3噪声扰动分析

我们检查我们的算法对于添加的图像噪声的强健性。为此，我们利用各种级别的加性高斯噪声来扰动多谱面部图像数据库。噪声具有最大图像亮度的0∶5与2％之间的、具有0.5％的增量的增加的标准偏差。在图15中，我们在可见(左手边)和红外(右手边)谱中绘制我们的算法和替选方法在各个级别噪声上的性能。对于我们的算法和Finlayson和Schaefer的算法[21]，我们采用从图像恢复的所有同质块片。用于我们的方法的正则化权重是α＝100000。如图所示，在可见谱中，我们的方法实现比其它两种方法更低的均值偏差角度，具有比[21]中的方法提供的更稳定的估计。此外，在红外谱中，我们的方法比Sunshine等人的方法[52]在恢复的发光体谱中胜过多于角度误差的两个标准偏差。

在随着的增加性能降级方面，Sunshine等人的方法最稳定，因为其考虑图像中的所有辐射率谱的上界，更少受高斯噪声的级别的影响。然而，该方法是三个替选方法中最不精确的方法。我们的方法和Finlayson & Schaefer的方法显现为随着噪声的级别而线性地降级，但Finlayson & Schaefer的方法以比我们的方法更慢的速率降级。这可以通过以下事实而解释：Finlayson & Schaefer的方法依赖于特征值分解，这等同于我们的方法具有零正则化项。虽然获得发光体的更强健的解，但他们的方法没有考虑表面几何和高光。我们的方法不仅可以估计发光体谱，而且能够计算所有二色参数，同时使用正则化保持合理的估计。

5.2高光检测/移除

在已经估计发光体谱的情况下，如上所述，我们可以采用上述过程来分离来自多谱显像中的反射的镜面分量的漫射。这在谱反射率在场景的小空间邻域内缓慢地变化的情况下是可行的。注意，该假设对于很多现实世界表面是有效的假设。因此，每个本地邻域可以被看作平滑同质表面。结果，在块片P中的位置u处的与几何形状无关的漫射分量被估计为k(u)L给定镜面性分量。

在此，我们对先前提出的室内人类面部图像数据库执行高光移除。如在先前部分中所述，我们通过估计发光体谱而开始。我们考虑每个图像像素周围大小为11x11的邻域，并且假设该邻域具有普通反射率。如上所述，阴影因子的平滑变化的实际执行承担将像素辐射重新投射到发光体谱和漫射辐射谱向量跨越的子空间上。在该试验中，我们采用使像素辐射与该二维子空间之间的距离的L₂范数最小化的投影。

图16示出对于数据集合中的多个采样面部图像估计的所得阴影和镜面性系数。顶行示出从高倾斜光方向照射的输入图像。第二和第三行分别示出从我们的估计方法得到的阴影和镜面系数。最后两行示出与[46]中的替选高光移除方法产生的相同结果。替选方法使用基于出自漫射照明中的Lambertian反射率的统计的概率框架。注意，替选方法仅可应用于符合Lambertian模型的单带灰度图像。因此，为了比较两种方法，我们将替选方法应用于每一个输入图像中的最亮的带。此外，为了符合[46]中关于光和观看者方向的共线性的假设，我们已经将替选方法应用于其中相机和发光体方向是线性依赖的数据集合中的面部显像。

如从图16观测的那样，我们的方法在典型高光位置(例如鼻子、眼睑和前额)处检测并且将镜面从漫射分量分离方面是成功的。此外，我们的方法产生捕获面部几何的变化的平滑不光滑的漫射图像。注意，我们的方法不需要先验的照射方向。另一方面，Ragheb和Hancock的方法[46]假设照射和观看方向的一致。因此，不能处理对于我们的方法示出的应用设置，其中，光源被放置在相对于相机的高倾斜位置处。如期望的那样，替选方法往往错过重要高光点，并且生成错误的镜面尖峰(specular spike)。

由于我们的方法对于输入表面块片假设均匀的反照率，因此其往往在高纹理化区域中并且沿着材料边界生成高光。然而，注意，粗糙和高纹理化区域的精细尺度浮雕纹理可能产生仅由精致测量检测到的镜面性，如Wang和Dana的著作[57]中讨论的那样。由于图10中的多谱图像的背景可以被看作图像中的纹理化区域，因此其可以产生我们的方法检测到的高光，如图16的第三行所示。

现在，我们将我们的关注转向我们的方法对于三色图像上的镜面性检测和移除的应用。我们将我们的方法的性能与采用偏微分方程来侵蚀(erode)在每个像素处的镜面性的另一高光移除方法[40]进行比较。

在图17中，我们将我们的方法与Mallick等人的高光移除方法[40]进行比较。如上所述，我们的方法对平滑同质表面比对纹理化区域(例如第二和第三行中的区域)执行得更好。在平滑表面上，虽然沿着闭合边界和材料边界可能存在错误镜面性，但我们的方法检测到大多数镜面尖峰。另一方面，替选方法产生更平滑的漫射分量，该漫射分量更直观地正确。此外，替选方法检测平滑表面上的更多镜面像素。

如图17可见，我们的方法能够检测期望的镜面尖峰，如第三列中所示。注意，我们的方法可能沿着边缘引入不期望的镜面性。这可以在第三行的镜面性映射中观测到。这是因为，在这些位置处的块片没有均匀的反照率。注意，第二行中的镜面性映射示出梨的潜在纹理变化，其可以是镜面性在水果皮上散射的原因。在第二列中，我们示出通过我们的方法恢复的漫射分量，其中，在镜面点处的漫射色彩与邻近非镜面像素接近。

5.3皮肤分段

最后，我们将我们的关注转向我们的算法恢复的谱反射率的照射不变性及其对于识别任务的应用。具体地说，我们关注对于皮肤识别使用通过上述方法提取的谱图像反射率。该任务可以被看作分类问题，其中，皮肤和非皮肤谱分别包括正类和负类。以此方式，我们可以断言在训练时间由我们的算法所恢复的以及在测试阶段由用于皮肤识别的方法所产生的发光体谱和表面反射率的强健性和一致性。

在该实验中，我们对作为用于分类的特征的使用我们的方法恢复的反射率谱所产生的皮肤识别性能与通过使用多个替选特征的分类器所产生的那些结果进行比较。为此，我们提出部分2.4中描述的过程所估计的我们的恢复的反射率的两个变化的结果。对于第一变化，提供真实状况发光体谱作为输入。对于第二变化，我们使用部分5.1.1中的试验获得的估计的发光体谱。通过比较这两种情况下的性能，我们可以评估当估计的发光体谱用作与真实状况进行比较时恢复的反射率的强健性。此外，我们还对这些变形与多个替选进行比较。这些中的第一个是通过以测量的真实状况发光体对原始图像辐射率谱进行归一化获得的谱反射率。第二情况是分类器应用于原始辐射率谱。最后，我们使用源自初始辐射率谱上经由线性判别分析(Linear Discriminant Analysis，LOA)执行子空间投影的主要分量。

对在先前描述的可见范围中捕获的面部图像数据库执行该试验。为了获得训练数据集合，我们从放置在主体的前面的高倾斜位置的光源下捕获的图像来选择皮肤和非皮肤区域。平均起来，存在从每个图像中的若干区域选择的856个皮肤像素和7796个非皮肤像素作为训练数据。随后，上述特征中的每一个从每个训练集合提取，并且用作对于具有径向基函数(Radial Basis Function，REF)内核的支持向量机(Support Vector Machine，SVM)分类器[13]的输入。此外，在训练时间使用5倍交叉验证来选择参数。为了分类皮肤相对非皮肤像素，所得SVM模型应用于同一主体的测试图像。已经在其它发光体条件下获取各自具有340x400像素大小的测试图像。

在图18中，我们提出使用上述输入特征获得皮肤分段映射。顶行示出多个采样主体的训练图像，其中红色矩形边界中包围的皮肤训练区域和蓝色矩形中包围的非皮肤训练区域。第二行从上到下示出用于顶行中的主体的测试图像。

注意，两行中的发光体方向和功率谱不同，如可以在阴影和阴影中观测的那样。实际上，放置在主体前面的高倾斜位置的光源照射训练图像，而具有与观看方向相同方向的正面光源照射测试图像。下面五行是使用上述特征训练的SMV产生的皮肤概率映射。在图中，较亮像素分类为更可能是皮肤。第三和第四行分别与我们的恢复的反射率的变化对应，其中，提供真实状况发光体和估计的发光体谱作为输入。第五、第六和第七行与通过分别以真实状况发光体谱、原始辐射率以及图像辐射率谱的前20个LDA分量对图像辐射率进行归一化所获得的反射率对应。

从图18，我们可以得出以下结论：实际上我们的方法恢复的皮肤反射率谱对于发光体功率和方向是不变的。这源于这样的事实：我们的方法提供的反射率特征产生视觉上最精确的皮肤映射。在很多情况下，正确地将非皮肤面部细节(例如眉毛和嘴)与皮肤区分。此外，两个反射率变化的结果良好一致。这是由于估计的发光体与真实状况之间的低差异，其通常在1度与3度之间偏离。

另一方面，用于关于第五行的结果的反射率特征虽然是发光体不变的，但仍然产生错误分类的皮肤像素。发光体谱和表面几何引起的变化可以解释这些特征获得的不良分类结果。这在面部边界和高光位置附近的像素处是明显的。由于以发光体功率对辐射率进行归一化没有实现几何独立性并且忽略了二色模型固有的镜面分量，因此这在第五行中是期望的。相对而言，我们的方法实现了没有几何和镜面性伪像的反射率的恢复。因此，其能够识别在入射余角(grazing angle)和镜面尖峰处的皮肤像素。此外，对于第六和第七行上的分类采用的归一化的原始辐射率谱及其LDA分量不是发光体不变的。因此，它们不能处理训练图像与测试图像之间的照明的改变。如最后两行所示，与我们的方法产生的反射率特征相比，这导致皮肤区域中的更多的负误识以及其它材料中的正误识。

为了提供量化的分析，在图19的表7中，我们示出按照分类率(classificationrate，CR)、正确检测率(correct detection rate，CDR)和错误检测率(false detectionrate，FDR)的上述皮肤分段机制的性能。正确检测率是正确地分类的皮肤像素的百分比。错误检测率是不正确地分类的非皮肤像素的百分比。分类率是正确地分类的皮肤和非皮肤像素的总百分比。该表示出在正面光源照射的数据集合中的所有主体的所有可见面部图像上的分段精度测量。该表的行与上述不同皮肤分类特征对应。如期望的那样，我们的方法恢复的反射率实现最高皮肤识别率。这与以上量化结果一致。此外，基于我们的方法的两个反射率变化即当使用估计的和真实状况光谱时)之间的整体性能差异(小于4％。这展现了我们的反射率估计方法对于输入发光体谱中的误差的强健性。如前，通过以发光体功率对辐射率进行归一化获得的反射率比原始辐射率谱及其LDA分量执行得更好。此外，辐射率特征及其LDA分量产生最多正误识和负误识。

接下来，针对多谱图像上的未监督的材料聚类(material clustering)我们检查我们的算法恢复的谱反射率的实用性。该试验可以被看作皮肤分段应用的延伸。此外，它是对于使用三色特征的皮肤聚类算法的补充，已经在文献[44]中另外描述。还对估计的反射率的聚类精度与测量的(真实状况)反射率的聚类精度进行比较。在此，我们对先前试验中提到的三个反射率特征基于确定性退火方法[27]执行聚类算法。这些特征包括在给定真实状况发光体谱的情况下估计的估计反射率、从估计的发光体谱产生的估计的反射率以及通过以真实状况发光体谱对辐射率谱进行归一化而获得的反射率。以用于所有材料的单个聚类来初始化聚类算法。随着算法进行，通过对包含不同材料的现有聚类进行划分来引入新的聚类。因此，所得多个聚类是与数据无关的，并且无需被指定为输入。

在图20中，我们示出若干人类主体的图像的聚类映射，其中，每一行与反射率特征对应。以不同的灰色阴影标记所得的材料聚类。实际上，通过估计的发光体谱估计的反射率特征与通过提供作为输入的真实状况发光体谱估计的反射率特征所产生的聚类结果之间存在高相似性级别，如行2和3所示。这再次展现了我们的反射率估计方法对于输入发光体谱中的误差的强健性。在这些聚类映射中，所有材料良好地彼此分离。此外，每一聚类内存在非常少的误分类的像素。在面部上，皮肤像素与邻近区域清楚地区分。注意，显示打印的面部(例如第三列中的面部)的背景区域正确地聚类为纸张。该结果说明，材料反射率的谱变化比三色色彩对于分类是更好的特征。注意，使用三色显像，实质上本来不可能使具有相同外观色彩的材料分离(例如真实面部与打印的面部)。在最后一行中，我们使用真实状况(测量的)反射率作为用于聚类算法的特征。与我们估计的反射率相比，测量的反射率产生更多噪声的聚类映射，其中，一些邻域中的像素分配给各种材料。换句话说，我们的反射率恢复方法通过减少原始反射率谱中的测量噪声而改进聚类性能。

以上示例已经示出可以如何估计图像的照明谱和光度测定的不变量。可以使用这些估计的一个或多个对图像执行的预处理和后处理。

图像的处理包括，但不限于：

图像编辑

白平衡

重新照明，例如用于将图像从温暖改变为寒冷，或改变图像的照明谱以将照明从日光改变为阳光或从阳光改变为钨的色彩操控

重新阴影

重新色彩化，例如基于已知色彩集合的特性或将一个图像的反射特性应用于另一图像将黑白图像改变为彩色

光源重定位

材料修改

高光移除

表面渲染

材料修改

形状估计

图像中捕获的场景或对象的形状的恢复

材料识别或分类

材料、图案或对象识别和分类

硬件校准

改进例如捕获图像的相机的光度测定的校准

以此方式处理图像可以被包括作为在相机200上或安装有用于在捕获之后处理图像202的软件的分离计算机上可用的功能。

在此描述的方法的应用包括以下领域：

数字摄影，例如图像编辑

制造，例如质量和生产控制。

产品分析，例如确定车辆是否已经处于事故中，以及

监督，例如面部标识和跟踪

在使用超谱相机的情况下，这些应用将一般具有这样的优点：减少产生相同分析所需的相机的数量，减少设置成本。

谱图像可以转换为彩色带呈现(例如RGB)，在该意义上，在此描述的方法可以用在彩色图像上。

本领域技术人员应理解，在不脱离权利要求中限定的范围的情况下，可以对具体实施例进行大量变形和/或修改。

应理解，可以使用各种技术来实现本公开的技术。例如，可以通过驻留在合适的计算机可读介质上的一系列计算机可执行指令来实现在此描述的方法。合适的计算机可读介质可以包括易失性(例如RAM)和/或非易失性(例如ROM、光盘)存储器、载波和传输介质。示例性载波可以采取沿着本地网络或公众可访问网络(例如因特网)传送数字数据流的电子、电磁或光信号的形式。

还应理解，除非如从以下讨论明显的另外具体地声明，应理解，贯穿说明书，使用例如“估计”或“处理”或“计算”或“演算”、“优化”或“确定”或“显示”等术语的讨论指的是将表示为计算机系统的寄存器和存储器内的物理(电子)量的数据处理并且变换相似地表示为计算机系统存储器或寄存器或其它这样的信息存储、传输或显示设备内的物理量的其它数据的计算机系统或类似的电子计算设备的动作或处理。

因此，这些实施例被看作在所有方面是说明性而非限制性的。

REFERENCES

[3]BARNARD，K.，MARTIN，L.，COATH，A.，AND FUNT，B.Y.A comparison ofcomputational color constancy Algorithms-Part II：Experiments with imagedata.IEEE Transactions on Image Processing II，9(2002)，985-996.

[4]BARNARD，K.，MARTIN，L.，FUNT，B.，AND COATH，A.A Data Set for ColourResearch.Color Research and Application 27，3(2002)，147-151.

[6]BOYD，S.，AND VANDENBERGHE，L.Convex Optimization.CambridgeUniversity Press，2004.

[7]BRAINARD，D.，AND WANDELL，B.Analysis of the retinex theory of colorvision.Journal of Optical Society America A 3(1986)，1651-1661.

[11]BROOKS，M.，AND HORN，B.Shape and source from shading.In MIT AI Memo(1985).

[13]CHANG，C.-c.，AND LIN，cr.LlBSVM：a library for Support VectorMachines，2001.

[14]CHANG，J.Y，LEE，K.M.，AND LEE，S.U.Shape from shading using graphcuts.In Proc.of the Int.Conf on Image Processing(2003).

[15]CIE.Commission Internationale de l′Eclairage Proceedings，1931.Cambridge University Press，Cambridge，1932.

[17]FERRIE，F.，AND LAGARDE，J.Curvature consistency improves localshading analysis.CVGIP：Image Understanding 55，I(1992)，95-105.

[21]FINLAYSON，G.D.，AND SCHAEFER，G.Convex and Non-convex IlluminantConstraints for Dichromatic Colour Constancy.CVPR I(200I)，598-604.

[22]FINLAYSON，G.D.，AND SCHAEFER，G.Solving for Colour Constancy usinga Constrained Dichromatic Reflection Model.International Journal of ComputerVision 42，3(2001)，127-144.

[23]FORSYTH，D.A.A novel algorithm for color constancy.InternationalJournal of Computer Vision 5，I(1990)，5-36.

[24]HEALEY，G.Estimating spectral reflectance using highlights.Imageand Vision Computing 9，5(October 1991)，333-337.

[26]HOAGLIN，D.c.，MOSTELLER，F.，AND TUKEY，J.W.Understanding Robust andExploratory Data Analysis.Wiley-Interscience，2000.

[27]HOFMANN，T，AND BUHMANN，M.Pairwise data clustering by deterministicannealing.IEEE Tansactions on Pattern Analysis and Machine Intelligence 19，I(1997)，1-14.

[28]HORN，B.K.P.Robot Vision.MIT Press，Cambridge，Massachusetts，1986.

[29]HORN，B.K.P.，AND BROOKS，M.1.The variational approach to shape fromshading.CVGIP 33，2(1986)，174-208.

[30]HUBER，P.1.Robust Statistics.Wiley-Interscience，1981.

[31]IKEUCHI，K.，AND HORN，B.Numerical shape from shading and occludingboundaries.Artificial Intelligence 17，1-3(August 1981)，141-184.

[32]KLINKER，G.，SHAFER，S.，AND KANADE，T A Physical Approach to ColorImage Understanding.International Journal of Computer Vision 4，I(1990)，7-38.

[33]KLINKER，G.1.，SHAFER，S.A.，AND KANADE，T The Measurement ofHighlights in Color Images.International Journal of Computer Vision 2(1988)，7-32.

[34]KOENDERINK，1.1.，AND VAN DOORN，A.1.Surface shape and curvaturescales.Image Vision Computing 10，8(1992)，557-565.

[35]LAND，E.，AND MCCANN，1.Lightness and retinex theory.Journal ofOptical Society America 61，I(1971)，I-II.

[36]LAND，E.H.Recent advances in retinex theory.Vision Research 26，I(1986).

[37]LEE，H.-C.Method for computing the scene-illuminant chromaticityfrom specular highlights.Journal of the Optical Society of America A 3(1986)，1694-1699.

[38]LI，S.Z.Discontinuous mrf prior and robust statistics：acomparative study.Image Vision Computing 13，3(1995)，227-233.

[39]LIN，S.，AND SHUM，H.-Y.Separation of Diffuse and SpecularReflection in Color Images.Computer Vision and Pattern Recognition，IEEEComputer Society Conference on I(200I)，341.

[40]MALLICK，S.P.，ZICKLER，T，BELHUMEUR，P.N.，AND KRIEGMAN，D.1.Specularity Removalin Images and Videos：A POE Approach.In ECCV(I)(2006)，pp.550-563.

[42]NAYAR，S.，AND BOLLE，R.Reflectance based objectrecognition.International Journal of Computer Vision 17，3(1996)，219-240.

[44]PHUNG，S.L.，BOUZERDOUM，A.，AND CHAI，D.Skin segmentation using colorpixel classification：Analysis and comparison.IEEE Trans.PatternAnal.Mach.Intell.27，I(2005)，148-154.

[45]POGGIO，T，AND TORRE，Y.Ill-posed problems and regularizationanalysis in early vision，Tech.rep.，Cambridge，MA，USA，1984.

[46]RAGHEB，H.，AND HANCOCK，E.R.Highlight removal using shape-from-shading.In European Conf on Compo Vision(2002)，no.2351 in LNCS，pp.626-641.

[48]SHAFER，S.A.Using color to separate reflection components.ColorResearch and Applications 10，4(1985)，210-218.

[52]SUNSHINE，1.M.，PIETERS，C.M.，AND PRATT，S.F.Deconvolution of MineralAbsorption Bands：An Improved Approach.Journal of Geophysical Research 95，B5(1990)，6955-6966.

[53]TAN，R.T，NISHINO，K.，AND IKEUCHI，K.Separating reflection componentsbased on chromaticity and noise analysis.IEEE Trans.PatternAnal.Mach.Intell.26，10(2004)，1373-1379.

[54]TOMINAGA，S.，AND WANDELL，B.A.Standard surface-reflectance modeland illuminant estimation.Journal of the Optical Society of America A 6(April1989)，576-584.

[55]TOMINAGA，S.，AND WANDELL，B.A.Component estimation of surfacespectral reflectance.Journal of the Optical Society of America A 7，2(February1990)，312-317.

[57]WANG，1.，AND DANA，K.1.Relief texture from specularities.IEEETransactions on Pattern Analysis and Machine Intelligence 28，3(2006)，446-457.

[58]WORTHINGTON，P.L.，AND HANCOCK，E.R.New constraints on data-closeness and needle map consistency for shape-from-shading.IEEE Transactionson Pattern Analysis and Machine Intelligence 21，12(1999)，1250-1267.

Claims (29)

1.一种计算机实现的用于估计包括波长索引的谱数据的点的图像的照明谱的方法，所述方法包括：

基于以下来优化代价函数：

(i)所述波长索引的谱数据与基于针对不同二色反射模型参数值的二色反射模型的估计之间的差异，其中，所述二色反射模型至少包括照明谱参数，以及

(ii)用于基于强健内核或曲率一致性而约束所述不同二色反射模型参数值中的至少一个的项。

2.如权利要求1所述的方法，其中，由于所述项，以封闭形式求解所述代价函数的优化。

3.如权利要求1或2所述的方法，其中，所述项基于最小化曲率大小变化或阴影的梯度大小的变化。

4.如权利要求3所述的方法，其中，所述项基于作为平滑性误差的测度的曲率大小或所述阴影的梯度大小。

5.如权利要求1所述的方法，其中，由于所述项，与仅基于差异的平方的代价函数相比，通过向噪声破坏和图像获取误差提供减少的偏差，所述代价函数的优化相对轮廓是强健的。

6.如权利要求1所述的方法，其中，由于所述项，由于噪声引起的所述波长索引的谱数据中的误差相对于基于所述二色反射模型的估计具有减少的差异。

7.如权利要求5所述的方法，其中，对于每一点，所述强健内核的偏差基于所述图像中的点的阴影的确定的梯度大小与所述图像中的该点的局部邻域的阴影的确定的梯度大小之间的差异。

8.如权利要求1所述的方法，其中，所述强健内核在所述图像上估计的表面法线中使得所述代价函数相对轮廓呈现为强健的。

9.如权利要求1所述的方法，其中，所述强健内核基于Huber估计者、Tukey双加权函数或Lee自适应势函数。

10.如权利要求1所述的方法，其中，所述项将所述图像中的点的曲率约束为与所述图像中的该点的局部邻域的曲率相似。

11.如权利要求1所述的方法，其中，所述项给属于所述图像中的相同或相似形状或辐射率的点分配较大的权重。

12.如权利要求1所述的方法，其中，所述项基于所述图像中的阴影或辐射率引起的曲率。

13.如权利要求1所述的方法，其中，所述项基于限制邻近波长之间的不同二色反射模型参数值中的一个或多个之间的差异。

14.如权利要求1所述的方法，其中，由于变量的变化或使用中间变量来表示所述代价函数，因而所述代价函数是封闭形式。

15.如权利要求1所述的方法，其中，迭代地优化所述代价函数，并且以均匀照明谱或基于仅基于所述差异的另一代价函数的解来初始化所述照明谱参数。

16.如权利要求15所述的方法，其中，所述另一代价函数是根据Finlayson&Schaefer在文章“Convex and Non-convex Illuminant Constraints for Dichromatic ColourConstancy”中提出的方法。

17.如权利要求1所述的方法，其中，所述方法还包括：将优化所述代价函数的照明谱参数识别为所述图像的照明谱的估计。

18.如权利要求1所述的方法，其中，所述方法还包括：在与所述图像关联的进一步处理中使用所述图像的照明谱，所述处理是以下中的一个或多个：

图像编辑：

形状估计；

材料识别或分类；以及

硬件校准。

19.如权利要求1所述的方法，其中，所述图像的照明谱是已知的，所述方法包括：

(a)对于所述图像中的每一点，基于二色反射模型、所述照明谱和所述波长索引的谱数据来估计该点的表面反射率；以及

(b)对于所述图像中的每一点，基于所述波长索引的谱数据、所述二色反射模型、所述照明谱以及所述表面反射率来估计表面阴影和镜面性。

20.如权利要求19所述的方法，其中，步骤(b)包括：

对于所述图像中的每一点，基于所述图像中的该点的局部邻域的波长索引的谱数据来估计该点的反射率谱；以及

将所述图像的每一点的波长索引的谱数据拟合到所述照明谱和所述反射率谱跨越的子空间上。

21.如权利要求20所述的方法，其中，步骤(b)还包括：最小化所述波长索引的谱数据与所述照明谱和所述反射率谱跨越的所述子空间之间的差异。

22.如权利要求19所述的方法，其中，估计该点的表面反射率的步骤(a)是基于所述图像中的该点的局部邻域的波长索引的谱数据，从而所述局部邻域显现为缓慢地改变表面反射率。

23.如权利要求1所述的方法，还包括：选择用于在估计图像的照明谱中使用的图像的子图像，所述方法包括：

(a)标识所述图像的多个候选子图像；

(b)对于每一候选子图像，将基于二色反射模型的超平面拟合为该候选子图像的波长索引的谱数据；以及

(c)选择具有不偏离于拟合所述超平面大于确定的量的波长索引的谱数据的候选子图像。

24.如权利要求23所述的方法，其中，步骤(c)还包括：基于遍及所述波长的均值波长索引的谱数据的方差而选择也具有高对比度的子图像。

25.如权利要求23所述的方法，其中，步骤(c)的所述确定的量是具有确定百分比的不偏离于所述超平面的波长索引的谱数据的候选子图像。

26.如权利要求23所述的方法，其中，步骤(b)包括：拟合所述波长索引的谱数据的辐射率向量表示。

27.如权利要求1所述的方法，其中所述图像是三色、超谱或多谱图像。

28.一种计算机实现的用于估计图像的照明谱的方法，所述图像包括波长索引的谱数据的点，所述方法包括：

(a)通过以下操作来选择在估计所述图像的照明谱中使用的所述图像的子图像:

(i)标识所述图像的多个候选子图像；

(ii)对于每一候选子图像，将基于二色反射模型的超平面拟合为所述候选子图像的波长索引的谱数据；以及

(iii)选择具有不偏离于拟合所述超平面大于确定的量的波长索引的谱数据的候选子图像；以及

(b)仅使用所选择的子图像，通过基于所述波长索引的谱数据与基于针对不同二色反射模型参数值的二色反射模型的估计之间的差异来优化代价函数而估计所述照明谱，其中，所述二色反射模型至少包括照明谱参数，并且还基于以下中的任一个来优化所述代价函数：

(i)用于基于强健内核或曲率一致性而约束所述不同二色反射模型参数值中的至少一个的项；以及

(ii)所述差异是通过将所述波长索引的谱数据拟合为所述二色反射模型的超平面而确定的，其中所述超平面约束照明谱参数值，以及用于约束模型参数值的项；以及

(iii)用于约束模型参数的项，其中，迭代地优化所述代价函数，并且以均匀照明谱或基于仅基于所述差异的另一代价函数的解来初始化所述照明谱参数值。

29.一种包括可操作为执行如前述权利要求中的任一项所述的方法的处理器的计算机。