用于从对象的图像生成立体地形图的方法、系统和计算机程序产品

摘要

一种方法，用于从对象(107)的多个图像获得所述对象的表面的模型，其中所述对象由一个或多个照射源照射，该方法包括：·-获得图像的集合，每个图像包括像素阵列且表示对象的类似视图，其中所述类似视图在不同的照射条件下获得；·-依据下列未知的模型分量定义表达所述图像的已知像素强度值的模型函数：·-第一模型分量(A)，表示所述对象在每一个像素位置的反射率，并且是对于所述多个图像都相同的反射率值；·-第二模型分量(L)，表示每个图像的照射强度，其是对于每个图像的所有像素都相同的照射源强度值；·-第三模型分量(V)，表示特定的照射方向，对于每个图像是不同的，并且是对于每个图像的所有像素都相同的照射向量；·-第四模型分量(N)，表示所述对象表面在每一个像素位置的表面法向量方向，并且是对于所有图像都是相同的法线向量；·-执行一个或多个序列的最小化运算来最小化从所述图像的集合获得的像素值与使用所述模型函数计算得到的像素值之间的差函数，每个最小化运算通过允许所述模型分量(A、L、V、N)的其中之一变动而其他模型分量保持不变来执行；-输出所述第四模型分量(N)作为所述对象表面的所述模型。

说明书

技术领域

本发明涉及一种用于从对象的多个图像生成该对象表面的立体地形图(raisedreliefmap)的方法、系统和计算机程序产品。

背景技术

在图像处理领域，光度立体法是用于生成对象的表面的立体地形图的公知方法。光度立体法通过使用对象的多个图像进行操作，这些图像是在一个或多个照射源提供不同方向的照射时，以固定位置处的相机拍摄的。光度立体法的目标是通过对与拍摄的图像相关的数据应用特定的算法，来提供对象的表面的立体地形图。在生物分析的领域，例如，对象包括微生物。使用这样的方法得到的例如为微生物的对象的立体地形图允许鉴别该微生物的若干特点，例如其表面的纹理和地形(relief)的类型。因此，该算法产生该对象的表面的立体地形图或深度地图(deepnessmap)，也被称作深度图(depthmap)。

此外，在生物分析的特定领域中，执行生物样品的体外分析，以确定可能存在于生物样品中的微生物类型。特定环境与生物样品关联，以允许微生物的生长。特定环境例如可以位于皮式培养皿中。在包括生物样品的皮式培养皿培养以后，皮式培养皿的表面可以被成像系统照射，如在已公布的欧洲专利申请EP2520923中描述的那样。成像系统包括位于生物样品周围的若干照射源。每个照射源可以从不同方向照射生物样品。固定位置相机可以拍摄皮式培养皿的同一个表面的不同图像。与不同照射方向相关的多个照射源代表复杂的照射条件。

在现有技术中，如在参考文献[12]中公开的光度立体法包括这样的算法：当照射发生时，其考虑了与描述灰度图像的图像特性的三个分量有关的三个矩阵。第一矩阵A涉及表示入射光被对象的表面反射的比例的反射率值。当整个表面假定由相同材料制成时，矩阵A对于所有像素包括相同值。第二矩阵V表示对于被相机捕捉的每个图像而言照射源的位置。第三矩阵N表示在三维空间中每个图像的每个像素的法向量值(normalvalue)。然而，这样的算法仅在如果照射源为点光源并且照射源的方向已知时才起作用。在生物分析领域中，举例来说，照射源不是点光源并且照射源的位置或者未知或者不易于建模。因此，所述算法不能被应用到这样的图像。另外，生物对象表面的光学性质或取向不是已知的。因此，用于这样的对象的反射率不是已知的。此外，在该算法中没有考虑到与照射源的强度相关的参数。因此，在重构处理期间，与对象表面的特性相关的信息可能缺乏或不够，从而导致对对象表面的错误表示。结果，这样的算法并不总是适用于生物分析领域中拍摄的图像的光度立体重构，也不适用于由固定位置相机在复杂的照射条件或未知照射方向下拍摄的任何图像的光度立体重构。

参考文献[13]公开了一种用于光度立体法的改进算法。该算法考虑了与值L相关的第四矩阵，值L表示照射源的各种强度。然而，这样的算法仅在如果照射源的位置和/或图像的特定像素的法向量已知时才起作用。在一些情形中，这样的信息或者不可利用，或者难以获得。因此，近似值对应用该算法来说是所必需的。结果，如果如生物分析领域中通常的情形那样，照射条件是未知的，那么这样的算法不适于生成对象表面的精确图像。另外，近似值可能在图像重构处理结束时导致对对象表面图像的错误表示。

结果，存在这样的需求：改进标准光度立体法来允许生成对象表面的立体地形图，以在对象没有先前信息可利用且与多个具有不同照射方向的照射源关联的情况下实现该对象的表面的重构。

发明目的

本发明的目的是克服与现有技术关联的至少一些问题。

本发明的进一步目的是提供一种方法、系统和计算机程序产品，用于利用光度立体法产生对象的表面的立体地形图或模型，而不需要使用与照射对象的照射源的位置和/或方向相关的任何信息且不需要使用与对象的表面的法向量相关的任何信息。

发明内容

本发明提供一种用于产生对象的表面的立体地形图的方法、系统和计算机程序产品，如在所附权利要求中陈述的那样。

根据本发明的一方面，提供了一种方法，用于从对象的多个图像获得所述对象的表面的模型，其中所述对象由一个或多个照射源照射，该方法包括：

-获得图像的集合，每个图像包括像素阵列且表示所述对象的类似视图，其中所述类似视图在不同的照射条件下获得；

-依据下列未知的模型分量定义表达所述图像的已知像素强度值的模型函数：

-第一模型分量，表示所述对象在每一个像素位置的反射率，并且所述第一模型分量是对于所述多个图像都相同的反射率值；

-第二模型分量，表示用于每个图像的照射强度，所述第二模型分量是对于每个图像的所有像素都相同的照射源强度值；

-第三模型分量，表示特定的照射方向，所述第三模型分量是对于每个图像都不同但对于每个图像的所有像素都相同的照射向量；

-第四模型分量，表示所述对象表面在每一个像素位置的表面法向量方向，所述第四模型分量是对于所有图像都相同的法线向量；

-执行一个或多个序列的最小化运算来最小化从所述图像的集合获得的像素值与使用所述模型函数计算得到的像素值之间的差函数，每个最小化运算通过允许所述模型分量其中之一变动而其他模型分量保持不变来执行；

-输出所述第四模型分量作为所述对象表面的所述模型。

优选地，所述序列的最小化运算包括迭代循环，在所述迭代循环中每个所述模型分量在后续迭代之前轮流改变至少一次，所述迭代循环重复直到满足结束条件。

优选地，所述迭代循环包括一个或多个内部循环，其中所述模型分量的子集在所述模型分量的另一个子集变动之前变动多次。

优选地，一个内部循环包括改变所述第三模型分量和所述第四模型分量多次。

优选地，一个内部循环包括改变所述第一模型分量和所述第二模型分量。

优选地，改变所述第三模型分量和所述第四模型分量的最小化运算在改变所述第一模型分量和所述第二模型分量的最小化被执行前执行多次。

优选地，所述第一模型分量包括第一模型子分量和第一结构分量，其中所述第一结构分量调适所述第一模型子分量来复制所述对象的表面特性。

优选地，所述第二模型分量包括第二模型子分量和第二结构分量，其中所述第二结构分量调适所述第二模型子分量来复制所述照射源的照射条件。

优选地，所述第三模型分量包括第三模型子分量和第三结构分量，其中所述第三结构分量调适所述第三模型子分量来复制所述照射源的所述照射条件。

优选地，所述第一结构分量表示对于所述多个图像都相同的系数，该系数包括对于所述多个图像的每个像素都相同的反射率值。

优选地，所述第二结构分量表示对于来自所述多个图像中一个图像的所述像素都相同的强度系数。

优选地，所述第三结构分量表示与所述多个图像中一个图像关联的一个照射源的照射方向，这对于来自所述多个图像中一个图像的所有像素都是相同的。

优选地，所述一个或多个照射源包括红色通道、绿色通道和蓝色通道。

优选地，所述第一模型分量对于用于所述多个图像的每个像素和每个色彩通道都是相同的，并且所述第二模型分量对于用于每个图像的所有像素的每个色彩通道都是相同的。

优选地，所述第一结构分量包括对于用于多个图像的每个色彩通道的每个像素都相同的反射率。

优选地，所述第二结构分量表示对于用于每个图像的所有像素的每个色彩通道都相同的强度系数。

优选地，模型函数包括修改后的模型函数，其中与所述第一、第二、第三和第四模型分量相关的所述模型函数等同于与所述第一、第二和第三模型子分量和所述第四模型分量相关的所述修改后的模型函数。

根据本发明的第二方面，提供了一种计算机程序产品，包括指令，所述指令被执行时使得可编程数据处理设备来执行该方法的步骤。

根据本发明的第三方面，提供了一种成像系统，用于从对象的多个图像获得所述对象的模型，该成像系统包括：

-一个或多个照射源，在不同的照射方向上照射所述对象，其中每个图像与确定的照射方向关联；

-图像捕捉设备，用于检测所述多个图像中每个图像的每个像素的像素值；

-数据处理器，包括第一存储部件，用于存储下列未知的模型分量：

-第一模型分量，表示所述对象在每一个像素位置的反射率，所述第一模型分量是对于所述多个图像都相同的反射率值；

-第二模型分量，表示用于每个图像的照射强度，所述第二模型分量是对于每个图像的所有像素都相同的照射源强度值；

-第三模型分量，表示特定的照射方向，所述第三模型分量是对于每个图像都不同的但对于每个图像的所有像素都相同的照射向量；

-第四模型分量，表示所述对象表面在每一个像素位置的表面法向量方向，所述第四模型分量是法向向量，其是对于所有图像都相同的法线向量；

其中所述数据处理器执行一个或多个序列的最小化运算，以最小化在所述图像的集合上接收的像素值与使用定义的模型函数计算得到的像素值之间的差函数，所述定义的模型函数依据所述未知的模型分量表示图像的所述已知像素值；

并且

其中通过允许所述模型分量其中之一变动而其他模型分量保持不变时来执行每个最小化运算并且

用于输出所述法向量模型分量作为所述对象表面的所述模型。

优选地，所述数据处理器还包括第二存储部件，用于存储结构分量以调适一个或多个所述模型子分量来复制所述对象的表面特性和所述照射源的照射条件。

附图说明

现在将以举例的方式作出对附图的参考，其中：

-图1示出根据本发明一实施例的成像系统的简图，其示出用于对象的各种类型的照射方向；

-图2示出根据本发明一实施例的成像系统的简图，其示出成像系统的模块；

-图3示出根据本发明一实施例的模拟光度立体法的数学公式；

-图4示出根据本发明另一个实施例的模拟光度立体法的数学公式；

-图5示出根据本发明一实施例的方法的流程图；

-图6示出根据本发明一实施例的第一对象的初始图像的例子；

-图7示出根据本发明一实施例的第二对象的初始图像的例子；

-图8示出根据本发明一实施例的第三对象的初始图像的例子；

-图9示出根据本发明一实施例的第四对象的初始图像的例子；

-图10示出根据本发明一实施例的图6的第一对象的立体地形图；

-图11示出根据本发明一实施例的图7的第二对象的立体地形图；

-图12示出根据本发明一实施例的图8的第三对象的立体地形图；

-图13示出根据本发明一实施例的图9的第四对象的立体地形图；

-图14示出根据本发明一实施例的第五对象的初始图像的例子；

-图15示出根据本发明一实施例的、考虑到第五对象的三个初始图像时的图14的第五对象的立体地形图；

-图16示出根据本发明一实施例的、考虑到第五对象的五个初始图像时的图14的第五对象的立体地形图；

-图17示出根据本发明一实施例的、考虑到第五对象的十个初始图像时的图14的第五对象的立体地形图；

-图18、图19、图20和图21示出根据本发明一实施例的第六对象的四个初始图像，每个初始图像与不同的照射方向相关；

-图22示出根据本发明一实施例的考虑到图18、图19、图20和图21的初始图像时的第六对象的立体地形图的俯视图；

-图23示出根据本发明一实施例的、如图22示出的第六对象的立体地形图的侧视图；

-图24、图25、图26和图27分别示出根据本发明一实施例的各种对象的立体地形图的侧视图；

-图28示出一表，其用于比较在本发明的第一实施例和第二实施例中以本发明的方法获得的不同重构表面；

-图29a、图29b和图29c分别涉及初始图像、参考标准法线(groundtruthnormal)图像和深度图图像；

-图29d示出一表，其用于在考虑到不同的初始化运算时比较以本发明的方法获得的不同重构表面；

-图30示出一表，其用于在考虑到不同的结束条件时比较以本发明的方法获得的不同重构表面；

-图31、图32、图33和图34示出根据本预处理一实施例的同一对象的若干图像，每个图像包括在不同位置的明亮区域；

-图35、图36和图37分别在红色/绿色/蓝色空间示出根据本预处理一实施例的图31的表示；

-图38示出根据本预处理一实施例的一曲线图，表示一图像的所选择像素的像素值与对于特定色彩通道的其他图像中的对应像素的其他像素值的比值；

-图39示出根据本预处理一实施例的用于从一图像除去明亮区域的S型(Sigmoid)函数的表示；

-图40示出根据本预处理一实施例的用于从一图像同时除去明亮区域和暗区域的S型函数的表示；

-图41示出根据本预处理的一个方面的方法的流程图；

-图42、图43、图44和图45示出根据本预处理一实施例的图31、图32、图33和图34的同一对象的若干修正后图像，各明亮区域被修正；

-图46、图47和图48分别示出根据本预处理一实施例的、图42在红色/绿色/蓝色空间的表示；

-图49示出包括明亮区域的对象的表示；

-图50示出根据本预处理一实施例的在图49中示出的对象的细节的表示，其中表现了多个明亮区域；

-图51示出根据本预处理一实施例的、在图52中示出的细节的修正后图像的表示；

-图52、图53和图54示出根据本预处理一实施例的同一对象的三个不同图像，其中对象包括特定的明亮区域和暗区域；

-图55、图56和图57分别涉及根据本预处理一实施例的图52、图53和图54的修正后图像；

-图58示出根据本预处理一实施例的、与在图53中示出的图像的红色通道相关的初始图像I1的表示；

-图59示出根据本预处理一实施例的、与在图26中示出的图像的红色通道相关的参考图像I2的表示；

-图60示出根据本预处理一实施例的、包括实线和虚线的曲线图，实线表示与在图58中示出的初始图像I1的感兴趣像素的线有关的强度值或像素值，虚线表示与在图59中示出的参考图像I2的感兴趣像素的线相关的强度值或像素值；

-图61示出根据本预处理一实施例的、使用在图59中示出的参考图像I2修正在图58中示出的初始图像I1而得到的修正后图像I1c；

-图62示出根据本预处理一实施例的、包括实线和虚线的曲线图，实线表示平均值W，并且虚线表示用于在图58中示出的初始图像I1和在图59中示出的参考图像I2的修正后函数f(W)；

-图63示出根据本预处理一实施例的、包括实线和虚线的曲线图，实线表示用于在图58中示出的初始图像I1的感兴趣的像素线的强度值或像素值，虚线表示用于在图61中示出的修正图像I1c的感兴趣的像素线的强度值或像素值；

-图64示出根据本预处理一实施例的、包括实线和虚线的曲线图，其中实线表示作为在图58中示出的初始图像I1和在图59中示出的参考图像I2的比值而估算的平均值W，并且虚线表示作为如图61中示出的修正后图像I1c与如图59中示出的参考图像I2的比值而估算的修正后的平均值Wc；

-图65示出根据本预处理一实施例的、以不同类型的虚线示出的f(W)的三个表示，其中k＝0.75且α＝100而为固定值，并且其中阈值V＝1、1.12和1.5；

-图66示出根据本预处理一实施例的、以不同类型的虚线示出的f(W)的三个表示，其中v＝1.25且α＝100而为固定值，并且其中k＝1、0.75和0.5；

-图67示出根据本预处理一实施例的、以不同类型的虚线示出的f(W)的三个表示，其中V＝1.25且k＝0.75而为固定值，并且其中α＝10³、10²和10；

-图68表示根据本预处理一实施例的在不同照射方向捕捉的对象的人工表面；

-图69表示根据本预处理一实施例的、其中明亮区域需要修正的其中一个初始图像；

-图70表示根据本发明一实施例的、由与空间上均匀分布的七个照射源对应的七个图像的镜面分量的叠加形成的人工图像；

-图71a表示根据本发明一实施例的、使用从空间上均匀分布的七个照射源获取的七个参考图像修正的图69的修正后图像；

-图71b表示根据本发明一实施例的、沿着在图43a中示出的水平白线的用于散射分量和用于软修正结果的像素强度的分布；

-图72表示根据本发明一实施例的、由与空间上随机分布的七个照射源对应的七个图像的镜面分量的叠加形成的人工图像；

-图73a表示根据本发明一实施例的、使用从空间上随机分布的七个照射源获取的七个参考图像修正的图69的修正后图像；

-图73b表示根据本发明一实施例的、沿着在图73a中示出的水平白线的用于散射分量和用于软修正结果的像素强度的分布；

-图74表示根据本发明一实施例的、由与位置紧邻的七个照射源对应的七个图像的镜面分量的叠加形成的人工图像；

-图75a表示根据本发明一实施例的、使用从位置紧邻的七个照射源获取的七个参考图像修正的图41的修正后图像；

-图75b表示根据本发明一实施例的、沿着在图75a中示出的水平白线的用于散射分量和用于软修正结果的像素强度的分布；

-图76表示根据本发明一实施例的一个图像的镜面分量，其对应于与要被修正的镜面区域相对的照射源；

-图77a表示根据本发明一实施例的、使用从与要被修正的镜面区域相对的一个照射源获取的一个参考图像修正的图69的修正后图像；

-图77b表示根据本发明一实施例的、沿着在图77a中示出的水平白线的用于散射分量和用于软修正结果的像素强度的分布；

-图78表示根据本发明一实施例的、由与空间上随机分布的九十九个照射源对应的九十九个图像的镜面分量的叠加而形成的人工图像；

-图79a表示根据本发明一实施例的、使用从空间上随机分布的九十九个照射源获取的九十九个参考图像修正的图69的修正后图像；

-图79b表示根据本发明一实施例的、沿着在图79a中示出的水平白线的用于散射分量和用于软修正结果的像素强度的分布；

-图80a、图80b、图80c以及图80d表示根据本预处理一实施例的圆形对象的四个初始图像，包括处于不同位置而没有任何重叠区域的明亮区域；

-图81a、图81b、图81c和图81d表示根据本预处理一实施例的图80a、图80b、图80c和图80d的圆形对象的图像，其中明亮区域已经使用现有技术中公知的修复处理进行修正；

-图82a、图82b、图82c和图82d表示根据本预处理一实施例的经过应用与软修正相关的本发明方法之后的图80a、图80b、图80c和图80d的圆形对象的修正后图像；

-图83表示根据本预处理一实施例的经过不对明亮区域做任何修正的表面重构方法处理之后的图80a、图80b、图80c和图80d的圆形对象的表面的等高线的示意图；

-图84表示根据本预处理一实施例的经过以现有技术中公知的修复处理重构之后的图80a、图80b、图80c和图80d的圆形对象的表面的等高线的示意图；

-图85表示根据本预处理一实施例的使用现有技术的科尔曼(Coleman)和简(Jain)算法的图80a、图80b、图80c和图80d的圆形对象的表面的等高线的示意图；

-图86表示根据本预处理一实施例的、在经过根据本发明一实施例的重构方法处理之后的对象的等高线的示意图，该方法中带有V_高、k_高、k_低和α_高、α_低的特定设定；

-图87表示根据本预处理一实施例的四个图像的叠加的示意图，其中每个图像包括镜面区域，并且其中没有镜面区域与另一个镜面区域重叠；

-图88表示根据本预处理一实施例的四个图像的叠加的示意图，其中每个图像包括镜面区域，并且其中所有镜面区域至少部分地与另一个镜面区域重叠；

-图89表示四个图像的叠加的示意图，其中每个图像包括镜面区域，并且其中至少一个镜面区域完全与另一个镜面区域重叠。

具体实施方式

基于特定的例子，下面的描述以充分清楚且完整的方式公开了本发明。然而，描述不应该被理解为将保护范围限制为下面描述的特定实施例和例子。

在下面的描述中，术语"像素"涉及图像中的位置，并且术语“像素值”和“像素强度值”涉及在所述图像中该像素的对应捕捉的强度的检测值。

在下面的描述中，术语“初始图像”和“图像”涉及对象的图像。一对象可以与多个初始图像相关联。在本发明中，初始图像的数量将应该至少等于三。对象的初始图像对应于照射源的特定照射方向或对应于多个照射源的主照射方向。另外，各个初始图像由固定位置相机在确定位置拍摄。描述的实施例包括三个初始图像，即图像1、2和3，由具有固定位置的相机在不同的照射方向下拍摄。应该理解本发明也适用于处理多个初始图像。初始图像包括多个像素和处于RGB空间中的三个色彩通道，即红、绿和蓝。作为例子，下面的计算在每个初始图像中仅考虑一个像素，例如像素1。应该理解计算也可以应用于初始图像的多个像素。术语“ni”涉及对象的初始图像的数量。术语“np”涉及每个初始图像中的像素的数量。

图像可以定义为两个分量的组合。第一分量是散射分量，其与光线进行没有任何优选方向的反射的那部分对应。第二部件是镜面分量，其与光线在优选方向进行反射的那部分对应。本发明的目标为将镜面分量去除来重新得到(retrieve)散射分量。

另外，在下面的描述中，词语“对象”可以涉及在各个照射方向被照射的任何对象并且其中对象的初始图像由固定位置相机拍摄。对象可以对应于一个或多个例如为微生物的小尺寸的对象或例如为一瓶水、动物等更大尺寸的对象等。

术语“阴影”涉及与由于光线的位置和对象的表面而导致的对象上的光线分布对应的信息。对于本发明的图像，该信息可以用来以本发明的立体光度法来重构对象表面。

在下面的描述中，对象的表面以处于正交空间x、y、z的立体地形图表示，其中z是x和y的函数，例如z＝f(x,y)。

图1示出包括一个或多个基本单元102的成像系统100。每个基本单元102包括用于产生与RGB色彩通道相关的红色(R)、绿色(G)和蓝色(B)照射的光学器件和控制电路103。每个基本单元102包括相应的红色、绿色、蓝色照射源104、105、106，它们是可完全独立地被控制。照射源104、105、106例如可以产生环形照射、逆环形照射(inverseannularillumination)或横向照射。

成像系统100可以可选地包括用于保持对象107的保持单元101。照射源104、105、106位于对象上方。

考虑到照射源104、105、106的取向，照射源104、105、106的光轴被布置为使得它们在不垂直于对象107的方向上施加照射。

成像系统100也包括图像捕捉设备108，例如相机，其指向保持单元101。来自任何基本单元102的照射源104、105、106的任何组合的照射能指向该对象。然后图像捕捉设备108可以捕捉已经被照射的任何对象107的初始图像。

如在图1中示出的，环形光束109也可以照射对象107并且其由一个基本单元102生成。

照射源104、105、106可以为任何所选的类型，例如工作在红光(R)、绿光(G)和蓝光(B)频率的发光二极管(LED)；简单的白光源；紫外线(UV)源或任何其他适当的辐射源。在任何位置的光源的数量可以与这里所示的和描述的不同。每个RGB照射由工作在各自频率的三个LED构成。可替换地，一个LED可以生成RGB照射。对于不同的照射源，可能存在RGBLED的不同组合。可替换地，可以布置二极管线或二极管条来生成一个照射方向。

除了在各个单元102中控制它的功能的控制电路和光学器件103之外，也可能存在全面控制系统110，其至少包括计算机、显示单元、处理模块和图像增强算法以及任何其他的处理或技术。

控制系统110可以被用来控制用于特定应用的照射源的选择。另外，控制系统可以对不同的应用施加不同的图像增强技术和图像处理过程。图像增强技术是用于例如在生物领域增强初始图像的质量或用于制作为分析图像内容的专家可见的相关信息的方法和技术。例如，图像处理涉及施加到图像矩阵的数学运算，来提供来自初始图像的信息的提取，以例如提供决策支持或使能自动决策。

控制系统110可以被用来执行用于成像系统100的任何其他功能和/或控制操作。这些包括但不限于：

-控制亮度的水平；

-控制红、绿、蓝分量的平衡；

-控制曝光时间；

-控制照射组合；

-测试系统；

-校准系统；以及

-基于分析的使用和目标的任何其他适当的控制。

本发明的目标是提出改进的光度立体法，用于在不知道对象的表面的特性和照射条件时从对象的多个初始图像生成为对象的表面立体地形图的图像，立体地形图也被称作凸版地图(raisedmap)。

如图2中所示，图像捕捉设备108包括强度检测设备112，用于选择图像中的每个像素，并确定该像素的对应像素值或强度值I。强度检测设备112检测或观察确定数量的图像中每个像素对于每个不同的色彩通道R、G、B的强度，以生成表示用于所有初始图像的所有像素的检测到或观察到或捕捉到的强度值的矩阵。

在第一实施例中，强度检测设备112选择例如为图像1的初始图像。然后强度检测设备112选择特定的色彩通道，例如色彩通道R。强度检测设备112检测图像1中用于色彩通道R的每个像素的像素值。当考虑在图1中的像素1的特定例子时，检测到的像素值可以称为IR1(像素1)。

强度检测设备112以相同的方式为图像1继续处理色彩通道G和B。以类似的方式，检测到的像素值例如可以称为用于色彩通道G的IG1(像素1)和用于色彩通道B的IB1(像素1)。

然后强度检测设备112为图像2和3执行相同的步骤。强度检测设备112总是在所有其他初始图像中选择对应像素来确定在所有不同的初始图像中的同一像素的对应强度。应该理解术语“同一像素”意指该像素在所有初始图像中可以位于相同的空间位置。因此，图像1中的像素1与图像2和3中的像素1位于相同的位置。

检测到的像素值参看下列注释，例如当考虑像素1时∶

-对于图像2：

IR2(像素1)

IR2(像素1)

IB2(像素1)；

-对于图像3：

IR3(像素1)

IG3(像素1)

IB3(像素1)。

强度检测设备112可以将检测到的强度存储在对应的矩阵中，对应的矩阵例如为矩阵IR1，其表示用于图像1和关联的色彩通道R的所有像素值。

以类似方式，图像1中与色彩通道G关联的所有像素的像素值被组织成矩阵IG1，并且图像1与色彩通道B相关的像素值被排列成矩阵IB1。

以相同的方式，对应于图像1的像素的像素值被排列成矩阵IR1、IG1和IB1，对于所有图像的所有像素的检测到的像素值被排列成对应的矩阵，其例如称为：

IR2，用于色彩通道R和图像2

IR3，用于色彩通道R和图像3

IG2，用于色彩通道G和图像2

IG3，用于色彩通道G和图像3

IB2，用于色彩通道B和图像2

IB3，用于色彩通道B和图像3。

由于在对应矩阵中集聚的像素值的数量，强度检测设备112可以经过展开步骤来将该矩阵转换为包括数据线(lineofdata)的向量，并且因此便于基于像素值的进一步计算。因此，矩阵I具有数量与初始图像的数量相等的行，其中初始图像包括与三个色彩通道相关的三条线。矩阵I具有数量与对于每个初始图像都相同的像素的数量相等的列。

成像系统100也包括数据处理器114来基于像素值和确定的数学公式执行特定的计算。

成像系统100还包括用于存储模型分量的数据库的第一存储部件116，该数据库用于存储四个模型分量。

第一模型分量表示与第一矩阵A关联的初始图像的像素的反射率值。第一矩阵A包括像素的反射率值，每个反射率值与被对象的表面的每一个像素反射的入射光的比例相关。因此，该反射率值表示依赖于对象的表面的材料的光学性质。第一矩阵A有等于3×ni的数量的行，其中每行表示初始图像用于特定色彩通道的所有像素的反射率值。因此，三行表示三个色彩通道R、G、B。第一矩阵A包括数量与在初始图像中的像素的数量相等的行。

第二模型分量表示各个照射源的强度，并且与第二矩阵L关联。对于每次图像获取，仅使用一个照射源。因此，第二矩阵L是其中每行具有非空值的矩阵，非空值与对于每个初始图像都是独特的照射源的强度值相关。第二矩阵L包含数量与3×ni相等的行，来表示每个初始图像用于每个色彩通道的强度。第二矩阵L包括数量与数量np的值相等的列。

第三模型分量V表示每个照射源的照射方向。该第三模型分量包括第三矩阵V，第三矩阵V包括与标准正交空间(x、y、z)相关的三个分量Vx、Vy和Vz。因此，向量(Vx、Vy、Vz)表示在图像获取处理期间照射源的方向。第三矩阵V强制该条件：对于每个初始图像仅有一个用于所有RGB通道的特定或主要照射方向。举例来说，例如为LED的同一个照射源可以被用于三个色彩通道R、G、B，这三个色彩通道R、G、B对于同一个图像都有相同的照射方向。然后该照射源对于每个色彩通道发射特定波长。因此，V^T包括三行和数量等于3×ni的列。在本发明中，由于基于下面指出的数学函数f和g的计算约束，转置矩阵V^T将被使用。然后将第三模型分量V以转置矩阵V^T的形式应用在下面的函数f和g中。因此，第三矩阵V有数量等于3×ni的行，其中一行表示一色彩通道，并且由于初始图像被表示为具有为R、G和B的三个色彩通道，所以初始图像被表示为三行。第三矩阵V包括三列，其中每列涉及标准正交空间x、y、z的轴。

第四模型分量表示图像的像素的三维法向量或3D法向量的位置。3D法向量是垂直于像素的正切平面的向量。由于法向量或表面取向涉及表面的特定特性，所以这样的特性不依赖于照射条件。因此，当考虑同一像素时，在所有图像中该法线向量相同。第四模型分量与第四矩阵N关联。第四矩阵N包括数量等于3的行，这些行与标准正交空间x、y、z的每个轴相关。第四矩阵N包括列，这些列表示在标准正交空间(x、y、z)中像素的每个法向量的坐标。第四矩阵N包括数量与数量np相等的列。

如下面示出的第一数学公式(1)基于上述的四个模型分量A、L、V和N来模拟光度立体模型：

f(A,L,V,N)＝A⊙L⊙(V^TN)(1)

其中⊙表示以分素(entrywise)方式的乘法运算，其中矩阵的每个值与另一个矩阵的对应值进行运算。

A、L、V和N必须被以本发明的方法重新得到来获得f(A,L,V,N)＝I。

图3示意性地示出每个矩阵A、L、V和N的结构。

数学公式(1)可以包括可选矩阵C，可选矩阵C在考虑到大尺寸的图像具有不均匀的照射条件时被使用。然后数学公式(1)可以涉及可选数学公式(1a)：

f(A,L,C,V,N)＝A⊙L⊙C⊙(V^TN)(1a)

将会理解，根据本发明一实施例，初始图像涉及小尺寸的图像，其中光线的分布可以被认为是均匀的。因此，被使用的数学公式没有包括可选矩阵C并参见数学公式(1)。

由于大量的数据在用于四个模型分量的所有矩阵A、N、L和V^T之间被处理，所以上述数学公式(1)还可以包括特定的计算最优化来避免增加计算时间。

因此，数据处理器114还包括第二存储器118用于存储结构分量的数据库，该数据库用于存储对计算进行最优化的至少三个结构分量。结构分量W通过考虑并模拟初始图像的物理特性和约束来调适该数学公式(1)，这通过考虑照射条件和该对象的表面的特性来实现。

第一结构分量W_A表示用于多个初始图像的相同反射率系数并仅包括为0和1的值。由于多个初始图像的每个像素对于相同色彩通道R、G或B具有相同的反射率值，所以如下所示，矩阵W_A具有数量等于3×ni的行，例如：

$>W_A=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\\ 1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\\ 1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)$>其中ni＝3

第二结构分量W_L表示对于每个色彩通道多个初始图像的所有像素都相同或相等的强度系数。

因此，因为每个照射源为相等的强度，那么W_L为长度例如等于数量np的向量，举例来说:

W_L＝[111(...)111]

第三结构分量W_V表示对于每个像素和用于每个色彩通道，与来自多个初始图像的一个图像关联的一个照射源的照射方向。

当考虑到三个初始图像时W_V显示如下。

$>W_V=\left(\begin{array}{ccc}1& 0& 0\\ 1& 0& 0\\ 1& 0& 0\\ 0& 1& 0\\ 0& 1& 0\\ 0& 1& 0\\ 0& 0& 1\\ 0& 0& 1\\ 0& 0& 1\end{array}\right)$>

最优化过程被执行来避免用于每个矩阵A、L和V的相等行和列的重复，并通过将每个矩阵A、L和V分解为模型子分量和对应结构分量W的相乘来执行，这些模型子分量分别为矩阵A'、L'和V'，并且这些对应结构分量分别为W_A、W_V和W_L。

因此，这些模型分量如下(2)：

对于矩阵A：A＝W_A.A'

对于矩阵L：L＝L'.W_L

对于矩阵V：V＝Wv.V'^T

矩阵A'、V'和L'表示用于一个图像的与像素相关的唯一的数据集。矩阵W_A、W_V和W_L分别设置用于矩阵A'、V'和L'的系数，以在计算(2)中示出的乘法时重新得到矩阵A、V和L的内容。

因此，数据处理器114也运算对应的函数g，函数g现在包括如在下面数学公式(3)中示出的结构分量和模型子分量：

g(A',L',V',N)＝(W_A×A')⊙(L'×W_L)⊙((W_V×V'^T)⊙N)(3)

其中f(A,L,V,N)＝g(A',L',V',N)

本发明的目标是确定矩阵A'、L'、V'和N的值，来获取用于g(A',L',V',N')的结果值。最好的结果值被认为是尽可能接近I的值，其表示取自图像捕捉设备108的初始图像的强度值。因此，函数f或g具有已知的结果，其等于I的值。函数g(A',L',V',N)包括为W_A、W_L和W_V的已知结构分量。模型分量A'、L'、V'和N必须被赋值来获取尽可能接近观测值I的g的结果值。

图4示意性地示出每个模型分量A'、L'、V'、N和结构分量W_A、W_L和W_V的结构。

包括模型分量A'、L'、V'、N和对应的结构分量W_A、W_L和W_v的函数g(A',L',V',N)被用于描述根据如本发明一实施例那样的图5的方法。根据图5的该方法还可以使用函数f(A,L,V,N)。

因此，本发明提供一种方法，用于以特定方式确定A'、L'、V'或A、L、V和N的值来最小化观测得到的强度I和对应的计算强度g或f之间的误差。然后被观测对象的表面的立体地形图或表示可以通过使用与第四模型分量有关的获得的矩阵N来生成。

该方法包括基于如图5中所示的步骤来运行一算法。该算法涉及交替最优化(alternatingoptimization)，其属于现有试探算法(heuristicalgorithm)的一种，如参考文献[1]、2]和[3]中公开的那样，其包含通过连续设定除了一个变量之外的所有变量来对目标函数进行的最小化。交替最优化的试探原理被经常用于当在一定条件下(例如，非负性)一个针对观察矩阵因子分解为两个矩阵时提供盲源分离。

根据本发明，该算法包括四个分量矩阵表示，并且包括用于每个模型分量的附加约束。

该算法包括两个子最小化过程，如在下面描述的和在图5中详细示出的那样。

该算法包括与主循环变量i和结束条件关联的主循环。在本实施例中，用于主循环的结束条件包括主循环的迭代的数量n_it。数量n_it例如可以等于10。该主循环包括内部循环，内部循环为第一循环和第二循环。第一循环与第一循环变量j关联并且第二循环与第二循环变量k关联。

在步骤500，方法的模型分量A'、L'、N被初始化为等于对应的确定值A'_O、L'_O、N_O。

在步骤502，表示该方法的主循环变量的i的值被初始化为0。

在步骤504，将n_it的值与主循环变量i比较来确定n_it的值是否大于主循环变量i的值。

如果迭代循环的数量n_it的值大于主循环变量i的值，那么处理进程跟随步骤506。

在步骤506，第一循环变量j被初始化为0。

在步骤508，第一循环变量j的值与关联到第一循环的预定值p比较来确定p的值(例如p＝5)是否大于j的值。如果p的值大于j的值，那么该方法跟随步骤510。因此，j<p表示结束条件。

下面的步骤510和512表示第一子最小化过程。因此，j<p是用于第一最小化方法的结束条件。

在步骤510，应用例如为最小化运算的最优化过程。该最小化运算可以包括例如为二次损失(quadraticloss)函数的最小二乘最小化(leastsquareminimization)技术的应用。在本实施例中，应用二次损失函数来获得观测得到的强度I和使用函数g计算得到的强度之间的二次误差。在步骤510，该二次损失函数以模型分量V'的值的变动为基础。在步骤510，V'的特定值被赋值来提供观测得到的强度I和函数g之间的最小差值。模型分量A'、L'、N保持不变。V'的赋值被确定。

在步骤512，二次损失函数类似地以模型分量N的值的变动为基础。在步骤512，N的特定值被确定来提供观测得到的强度I和函数g之间的最小差值。模型子分量A'、L'、V'保持不变。值V'相当于来自步骤510的V'的赋值。N的赋值被确定。

在步骤514，第一循环变量j从j增加到j+1。

在步骤508，j的增加值与p的值相比较。第一循环运行直到j的值大于p的值。当j的值小于p的值时，步骤510再次以第一循环的新的迭代运行。V'的新的赋值由最小化运算确定来获得最优化的赋值V'，其给出I和g之间的最小差值，其中g考虑到在来自第一循环的上述迭代中的步骤512的N的赋值。

在第一循环的新的迭代中，在步骤512，N的新的赋值与步骤512的N的所有之前的赋值比较来确定来赋值N，其给出I和g之间的最小差值，其中g考虑到来自第一循环的上述迭代中的步骤510的V'的赋值。

第一循环提供用于V'和N的最优化赋值。

当j的值大于p的值时，该方法跟随在步骤516开始的第二循环。

在步骤516，第二循环变量k被初始化为0。在步骤518，第二循环变量k的值与关联到第二循环的预定值p比较来确定p的值(例如p＝5)是否大于k的值。如果p的值大于k的值，那么该方法跟随步骤520。

因此，k<p是结束条件。

下面的步骤520和522表示第二子最小化过程。因此，k<p表示用于第二最小化方法的结束条件。

在步骤520，二次损失函数应用于模型分量L'。在步骤520，L'的特定值被赋值来提供观测得到的强度I和函数g之间的最小差值。模型子分量A'、V'和模型分量N保持不变。V'和N的值相当于来自第一循环的V'的最优化赋值，这两者都给出了I和g之间的最小差值。L'的赋值被确定。

在步骤522，二次损失函数应用于模型分量A'。在步骤522，A'的特定值被赋值来提供观测得到的强度I和函数g之间的最小差值。模型子分量L'、V'和模型分量N保持不变。L'的值相当于来自步骤520的L'的赋值，并且V'和N的值相当于来自第一循环的V'和N的最优化赋值，这两者都给出了I和g之间的最小差值。A'的赋值被确定。

在步骤524，第二循环变量k增加到k+1。

第二循环运行，直到k的值大于p的值。当k的值小于p的值时，步骤520再次以第二循环的新的迭代运行。L'的新的赋值由最小化运算确定来获得最优化的赋值L'，其给出I和g之间的最小差值，其中g考虑到在第二循环的之前迭代中来自步骤520的A'的赋值。

当k的值大于p的值时，该方法跟随步骤526，其中主循环变量i增加到i+1。

然后该方法跟随步骤504来比较循环迭代的数量n_it与主循环变量i的值。该主循环运行，直到迭代的数量n_it小于主循环变量i的值。其中当主循环变量i的值大于重复的数量n_it的值时，然后该方法停止。

n_it表示用于本方法的特定结束条件。

可替换结束条件可以包括与在观测得到的强度I和计算得到的强度f(A,L,V,N)或g(A',L',V',N)之间获得的差值关联的阈值。因此，例如，当I和f或I和g之间的计算得到的差值等于或小于阈值或阈差值时，该方法停止。

另一个可替换结束条件可以以阈差值为基础。当使用阈差值时，在同一模型分量基于两个连续的迭代的两个值之间的差值等于或小于该阈值时，该方法停止。

该方法的目标是分开地和迭代地考虑用于每个模型分量A'、L'、V'和N的值，该值给出了在二次损失函数的计算中的最小差值。该方法基于主迭代循环、第一迭代循环和第二迭代循环来最优化模型分量A'、L'、V'和N的值的确定。主循环、第一循环和第二循环重复的目标是达到最好的收敛。因此，对于每个步骤510、512、520和522，数据处理器114在第一循环和第二循环的每个迭代处并为每个模型分量运行比较步骤。实际上，每次当增大特定的模型分量时，数据处理器114计算二次损失函数来获得差值，此时数据处理器114比较该获得的差值与之前获得的任意现有差值。因此，数据处理器114可以确定最小化差值并保有该最小化差值用于进一步的计算。

迭代p的数量的值是确定的数量，其可以对于j和k是相同的。第一循环涉及模型分量V和N，并且第二循环涉及用于改进该方法的最优化的模型分量A和L。

模型分量V、N涉及在对象的表面的方向和照射源104、105、106的方向之间的相互影响。这样的相互影响确定两个方向之间的角度。由于对应矩阵的大小，与这种角度的确定有关的计算运行相对较复杂地。

因此，第一循环提供与模型分量V和N或V'和N相关的计算来最优化与第二循环关联的计算和运算时间。第二循环提供与模型分量A和L相关的计算。A相当于与对象的表面关联的加权系数。L相当于与照射源关联的加权系数。因此，由于对应矩阵的大小，与这种模型分量有关的计算运行相对较容易。

模型分量A、L、V、N或A'、L'、V'、N在这样的第一循环和第二循环中的设置提供了该方法更好的收敛。因此，当与用于执行例如为最小化运算的最优化过程的模型分量A、L、V、N或A'、L'、V'、N的不同设置相比，计算时间减少了。

可替换设置可以为：V'、N、L'、A'或V'、N、A'、L'或V'、L'、A'、N或V'、L'、N、A'或V'、A'、N、L'或V'、A'、L'、N。

当主循环、第一循环和第二循环的所有迭代都已经进行过时，该方法停止。

数据处理器114重新得到由本发明的方法确定的N的值。然后通过采用在参考文献[4]、[5]、[6]和[7]中公开的算法中的一种来考虑用于每个像素的N的值，数据处理器114一个像素接一个像素地建构表面。

因此，数据处理器114提供被考虑的对象的表面的立体地形图。

图6、图7、图8和图9每个示出第一对象的初始图像。图10、图11、图12和图13示出根据本发明的图6、图7、图8和图9的对象的表面的立体地形图。在本实施例中，循环迭代的数量是n_it＝10，且第一循环变量和第二循环变量的值是p＝5。

本方法确保当考虑很多良好质量的图像时，提供模型分量A'、L'、V'和N的更好赋值。另外，在图像获取期间该照射源均匀的空间分布提供了更好的结果。图14表示为茶叶袋的对象的初始图像。图15示出考虑到在图14中示出的对象的三个初始图像时图14的对象的立体地形图。图16示出考虑到在图14中示出的对象的五个初始图像时图14的对象的立体地形图。图17示出考虑到在图14中示出的对象的十个初始图像时图14的对象的立体地形图。从目视检查来看，图17中的立体地形图包括对象的细节比在图15和图16中的立体地形图更多。当与初始图像比较时，立体地形图的表面的分析得到的可见结果显示了立体地形图的改进质量，且没有对一般表面方面产生显著影响。

图18、图19、图20和图21示出为微生物的第一对象的初始图像。图22示出图18、图19、图20和图21的对象的表面的立体地形图的俯视图。

图23表示图22的对象的侧视图。

图24、图25、图26和图27表示根据本发明的例如为微生物的各种对象的立体地形图。

上述描述的最优化涉及本发明方法的第一实施例，其提供了用于确定表面法向量和用于重构深度图的方案。

在上述第一实施例中，在该最小化步骤期间，没有考虑值约束和质量约束。每个模型分量可以利用最小二乘法进行赋值。在迭代过程的执行期间，在每个子最小化之后获得的值可以被投射到约束中。对于反射率分量A，全部小于零的值可看作零，并且大于一的值被设定为一。可以对光源强度执行类似的运算。L矩阵中小于零的值被设定为零或极小值。照射方向V和表面法向量N仅在第一最小化循环的赋值之后被归一化。紧跟在归一化之后，辅助的表面法向量可以通过平均值滤波器以3×3的滑动窗口滤波以实施质量平滑性约束。由于算法的迭代性质，这样的后赋值约束应用是可能的。如在下面的第二实施例中描述的并与约束的赋值有关的附加步骤可以被设置在该方法的最后以提供第一光度立体模型分量的近似。

在第二实施例中，本发明在最优化方法的执行期间，提供用于考虑光度立体模型的值约束和质量约束的方案。

为了集成用于光度立体分量的值和质量限制，用于发现反射率值、光源强度、光源方向和表面法向量矩阵的最小化函数考虑了辅助等式约束和不等式约束。用于每个模型分量的最小化运算描述如下。

在图5中，应使用用于V的所示公式来建立包括范数约束(normconstraint)的光源方向矩阵。

受约束于$>v_{k,x}^2+v_{k,y}^2+v_{k,z}^2=1$>

其中k＝1,..,ni，并且向量(v_k,x,v_k,y,v_k,z)相当于在第k个图像的获取期间使用的照射方向。

为了获得具有值范数约束和质量约束的表面法向量矩阵，用于在图5中的矩阵N的最小化过程应该以下面的条件使用：

$>\left(\begin{array}{c}{n_{x,j}}^2+{n_{y,j}}^2+{n_{z,j}}^2=1,\\ {n_{x,j}}^2+{n_{y,j}}^2=0,\\ \underset{p=-1}{\overset{1}{\Sigma }}\underset{l=-1}{\overset{1}{\Sigma }}{n}_j^T{n}_j^{(p,l)}=8,\\ withp\ne 0,l\ne 0\end{array}\right)$>

其中j＝1,...,np，nj相当于第j个像素的3D法向量，并且表示在初始非向量化图像中的第j个像素的邻居。

光源强度向量应该只包含那些大于0且小于图像最大可能值H的值。因此，如图5所示的矩阵L'应该以下面的条件使用∶

受约束于I_k>0且I_k<H，

其中k＝1,..,3ni，ni是输入彩色图像的数量，并且I_k相当于在第k个输入图像的获取期间使用的光源强度。

反射率矩阵应当只包含来自区间[0；1]的值。因此，用于在图5中的反射率近似值的最小化函数A'_i应该以下面的条件使用：

受约束于a_kj≥0且a_kj≤1，

其中k＝1,..,3ni，j＝1,..,np，ni是输入的彩色图像的数量并且np是在每个输入图像中的像素的数量。

具有不等式约束的最小化问题的解析，即反射率A'和光源强度L'可以通过分析用于包括拉格朗日乘数[9]的最优化误差函数的如参考文献[8]所示的Karush-Kuhn-Tucker(KKT)条件来获得。由于用于照射方向和表面法向量的等式约束的非线性性质，用于子最小化的最佳方案可以利用用于损失函数且包括拉格朗日乘数的Levenberg-Marquardt算法获得，如参考文献[10]所示。

光度立体模型解析提供了单元法向量，即在每个像素中的表面取向。为了从这些法向量获得深度图，必须采用一个现有的算法。深度图允许比法向量矩阵提供的结果表面的目视分析更容易的结果表面的目视分析。在试验期间，参考文献[4]、[5]、[6]和[7]中出现的算法被测试。结果参考文献[[4]中出现的综合算法是所有给出例子中最快的和最有效率的。本算法基于来自两个视角的表面偏导数表达式∶作为离散导数(邻近点之间的简单差值)和在表面法向量的基础上计算得到的导数。本算法因此被用于在交替最优化之后获得的法向量和可视的深度图之间的过渡。在所有的下列说明中，使用深度图的重构的质量被分析。

图28示出用于简单半球形对象和用于来自如参考文献[11]所示的MIT本征图像数据集的对象的重构结果的比较。对于当前的例子，来自图5的交错法被使用，且对于每次最小化具有10次全局迭代和5次内部迭代。获得的结果示出实施约束的必要性。对于例如为半球形的简单对象，没有可变值约束的交替最优化是足够的。然而，更复杂的表面要求使用约束。如图28所示，不同形式的约束实施方式可以比较不同的处理时间t：在每个没有约束的可变赋值之后，约束直接集成到每个最小化过程，并且投射到约束的空间。第二实施方式给出了接近于第一实施方式的结果，但是显著减少了赋值的处理时间。因此，在进一步的试验中，优选赋值的第二方式。在64位系统、英特尔(R)核心(TM)i7-2600，CPU3.40GHz，RAM为8G，使用MatlabR2012b的环境的情况下表示赋值的时间。

在图5中，步骤500与基于用于A'、L'和N的确定值的初始化步骤有关。

当没有初始信息可供使用时，交替最优化算法所必需的模型分量的初始化应该在输入图像的基础上进行。试验示出一些变量比其他变量的初始化的影响更大。最重要的是法向量矩阵N的初始化。所有需要初始值的例如为A'和L'的其他变量可以被选择为一个或任意值的矩阵。图29a涉及对象的初始图像。图29b涉及基于参考标准法向量(groundtruthnormal)的矩阵的图像。图29c涉及基于深度图的图像。图29d示出包括多个结果图像的表，多个结果图像被布置在表中，用于比较基于用于矩阵N的不同初始化标准的不同图像。每个初始化标准涉及表中的特定行，如下所述。

第一行No涉及到在交替过程的第一步骤使用的N的初始值。Nx相当于表示的色彩通道R，Ny相当于该色彩通道G，并且Nz相当于色彩通道B。

第二行“N”涉及由交替最优化算法重新得到时矩阵N。第二行涉及以R、G、B通道编码和以不同于行No的初始化得到的法向量的结果值。

第三行“SF”涉及对象的结果表面，其相当于基于结果法向量从行N获得的深度图。

第一列“Ones”涉及单元矩阵。

第二列“Randn”涉及正态分布的随机值的矩阵，例如遵循高斯分布的值的矩阵。

第三列“Rand”涉及均匀分布的随机值的矩阵。

第四列“对于Med的3D法向量”涉及包括用于法向量的不同类型初始化的矩阵。初始法向量等于每个像素中位强度值的3D法向量。看来结合了第四列的标准的算法提供了图像的最优质量。

图29d包含对于矩阵N的不同初始化的重构结果。因此，具有等于每个像素中位强度值的3D法向量的初始化法向量的重构表面显著地好于其他的初始化法向量。

在赋值例子中，反射率矩阵A'被选择为适当大小的[3×np]的单元矩阵。光源强度向量L'被作为所有图像通道的3ni个中值的向量。表面法向量的初始矩阵N'被定义为每一个像素强度中值的3D法向量。因此，矩阵N'为所需的[3×np]大小。该描述的初始化不是唯一可能的选择，而是其与提出的交错法结合，已经被证明对于获得合格表面是有效的。

在图5中示出的优选实施例中，主循环的结束条件是以迭代循环的定义次数为基础。可替换的结束条件可以包括初始数据和赋值变量之间的小的二次误差，以及从一个迭代到另一迭代时这些变量(或只有一个最重要的变量)的小的变化。图30包括示出三个不同行的表，这些行涉及不同的对象，并且其中第一列涉及参考标准表面，第二列涉及第一结束条件，第三列涉及第二结束条件，且第四列涉及第三结束条件和对应的处理时间t。结果表明有可能利用所有提到的标准获得满意的表面。取决于重构的表面复杂度、初始图像大小和质量，类似的重构结果从一个迭代到另一个迭代以等于10次算法迭代的it_max的值、分析模型以及法向量变量而得到。然后Err和Var与以经验确定的阈值进行比较。根据期待的结果和赋值时间约束，可以使用这些标准的一个。

可以被选择的最简单标准是迭代的数量。主循环中迭代的数量等于10，并且在每个子最小化中迭代的数量是5。对于这些值，获得了算法收敛，并且结果表面的质量对于进一步的分析是足够的。

用于改进初始图像的质量的附加特定图像预处理可以被与本发明一起实施，来在将本发明的方法应用到将要在本发明中使用的所述图像上之前，改进初始图像的质量。这样的附加图像预处理可以去除初始图像内的任何暗区域或明亮区域，这样的区域包括像素值，这些像素值可以导致深度地图的错误表示。

数据处理器114可以基于上述矩阵，计算与用于不同色彩通道的特定图像中的像素的像素值与在用于对应色彩通道的其他图像中的对应像素的像素值之间的比较有关的特定比值。为1的比值相当于有相同像素值的两个像素的像素值的比较。因此，当比值不同于1时，该比值确定在一图像中的第一像素值与在一个或多个其他图像中的同一像素的第二像素值之间的差异。实际上，如果该比值高于1，这意味着上述的像素相当于较明亮像素。例如，该比值可能为1.20，其表示等于1的比值的20％的增加。如果该比值小于1，这意味着上述的像素相当于较暗的像素。

对于图像1：

举例来说，比值包括IR1/IR2来在与用于相同色彩通道的图像2比较时确定在与色彩通道R关联的图像1中明亮区域的存在。

举例来说，该结果比值可以设置成如下的对应的比值矩阵，其中每个向量一个像素接一个像素地相关于图像1中的所有像素的比值，并且其中矩阵的每列对应于有关图像中的特定像素。

对于图像1：

IR1_比值＝(IR1/IR2；IR1/IR3)

IR2_比值＝(IR2/IR1；IR2/IR3)

IR3_比值＝(IR3/IR1；IR3/IR2)

对于图像2和图像3：

数据处理器114对于图像2和图像3中的每个像素迭代地执行类似的比值计算来获得用于图像2和3，以及分别用于对于每个图像1、2、3的其他色彩通道G和B的类似的比值矩阵。

当考虑N个图像的多样性时，数据处理器114为用于每个色彩通道的所有N个图像的所有像素迭代地执行这样的比值计算。

数据处理器114也可以基于上述获得的比值矩阵来计算为平均值的加权值W。

平均值相当于基于统计学测量的广义平均值。该广义平均聚集数据集的值来提供表征或表示数据集的结果数据值。该结果数据值与该数据集的任何数据值的类型一样。因此，该结果数据值可与该数据集的任何数据值比较。广义平均数可以涉及集中趋势测量，例如计算得到的比值的算术平均(即平均值(average))、加权算术平均、几何平均或中值。

在第一实施例中，计算用于每个图像的每个像素和用于每个色彩通道的加权值W，并且加权值W可以被定义为如下所示：

$>W=\underset{j=1}{\overset{p}{\Sigma }}{c}_j\frac{1}{I_{refj}}$>

在上述公式中，c_j是允许计算加权平均的系数。在算术平均的情形中，它等于1/p，其中p是参考图像的数量；

I_refj是多个的p个参考图像中的第j个图像；

I是要被修正的图像。

结果加权值W可以在如下的对应加权矩阵中设置：

对于色彩通道R：

WR1＝用于图像1中的所有像素的IR1_比值的平均

WR2＝用于图像2中的所有像素的IR2_比值的平均

WR3＝用于图像3中的所有像素的IR3_比值的平均

在上述例子中，WR1是集聚了用于图像1的红色通道的所有像素的值W的加权矩阵。

获得用于色彩通道G和B的类似加权矩阵。

当考虑N个图像的多样性时，数据处理器114为用于每个色彩通道的所有N个图像的所有像素迭代地执行这样的加权矩阵计算。

图38示出用于如图31所示的图像1的比值WR1的表示，其包括在图像1以它们的各自位置来识别的6000个像素。在该例子中，算术平均数被用来基于用作如图32、图33和图34中所示的参考图像的三个其他图像计算平均值WR1。因此，WR1的高值对应于位于像素1500到像素2000之间的明亮区域。

加权值W与等于1的阈值V关联来比较该加权值W与该阈值V，以确定在哪种情况中与加权值关联的图像包含明亮区域。

当加权值W等于阈值1时，这意味着在两个考虑的图像中同一像素的强度值是相同的。

当加权值W大于阈值1时，这意味着该像素可能涉及镜面像素。

当用于特定像素的加权值W小于阈值1时，这意味着该像素没有涉及镜面像素。

当加权值W接近值0时，这意味着该像素可能涉及暗区域。

加权值W可以与特定的边界值(marginvalue)Vm关联，例如与用于明亮区域的Vm_高和用于暗区域的Vm_低关联来限定两个对应的阈值区间。这样的阈值区间使不与镜面区域或阴影区域相关的特定的明亮区域和暗区域能被保持在图像中。实际上，举例来说，由于照射源与对象因为特定的照射方向或对象的特定性质而导致的相互影响，这样的特定明亮区域和暗区域可能出现。因此，当考虑阈值区间时，这样的特定区域不被认为是需要以本发明的方法修正的明亮区域或暗区域，并且在图像中保持不变。Vm_高和Vm_低的值是限定值，比如Vm_高等于0.2且Vm_低等于0.8。

第一阈值区间涉及W的较高值，其接近值1，且指示不存在明亮像素。基于边界值Vm_高的值，第一阈值区间涉及V从1到V_高的值，其中V_高相当于用于明亮像素的V的预定值。在本例子中，其中Vm_高等于0.2，第一阈值区间相当于V从1到1+Vm_高即1.2的值。因此，当加权值W不等于1并被包括在第一阈值区间中时，这意味着相关像素没有被认为包括对应于镜面像素的任何明亮像素。因此，具有由例如由于对象的特定性质的光线效果产生的高强度的像素值可以被考虑用于像素修正，但是在像素值修正以后保持几乎不变。

第二阈值区间涉及W的较小值，其接近值0，且指示存在暗像素。基于边界值Vm_低的值，第二阈值区间相当于V从V_低到1的值，其中V_低相当于用于暗像素的V的预定值。在本例子中，其中Vm_低等于0.8，第二阈值区间相当于V从1到1-Vm_低即0.2的值。因此，当加权值W不等于1并被包括在第二阈值区间中时，这意味着有关的像素没有被认为是包括对应于暗区域的任何暗像素。因此，具有由例如由于对象的特定性质或皮式培养皿材料的光线效果产生的低强度的像素值可以被考虑用于像素修正，但是在像素修正以后保持几乎不变。

阈值V表示来确定以本预处理的方法进行的像素值修正是否是需要的请求。

数据处理器114可以基于上述比值、加权值W和阈值V计算被确定的公式，以提供修正后像素值或修正后强度值Ic来替换考虑的像素的被确定的像素值I，该公式例如可以采用下列形式∶

$>Ic(im\_number,p\_number)=\frac{I(im\_number,p\_number)}{W_{im\_number,p\_number}^{f\left(W_{im\_number,p\_number}\right)}}$>

其中：

im_number是用于来自多个的N个图像的修正后图像的下标；

p_number是用于修正后的像素的下标；且

f(W)提供对于要被修正的明亮像素和暗像素的持续修正。S型函数可以使能这样的持续。

f(W)是例如为下面的S型函数的函数，并且可以表示大于0且小于1的值：

$>f\left(W_{im\_number,p\_number}\right)=\frac{k}{1+e^{-\alpha *(W_{im\_number,p\_number}-V)}}$>

其中∶

-k是与如下定义的修正的水平的调节相关的系数；

-V是与阈值区间在上面一起描述的阈值；

-α是定义S型函数的过渡的锐度的参数；

-W是上面描述的权重值。

在上面的公式中，k是用于调节修正的水平的预定系数，来调节由对象的修正后图像中的所讨论像素导致的失真，并且用于调节在图像内的其他像素上的像素修正的影响的水平；k的值大于0且小于1。k的高值提供了像素的显著修正。举例来说，如图66所示，对于明亮像素和暗像素，例如为k＝0.75的k的预定值可以分别称为k_高和k_低。

基于上面的公式，当权重值W相对高于阈值V时，函数f(W)接近k的值，其提供相关像素值的修正。

以类似的方式，当权重值W相对低于阈值V时，函数f(W)接近值0，其也提供相关像素值的严格修正，这没有显著地改变像素值。

在上面的公式中，α是定义S型函数(sigmoidfunction)的过渡的锐度的参数。这意味着α可以变动来调节该修正以平滑在如图67中示出的图像中对象的修正后区域与没有修正区域之间的可见过渡。过渡的平滑性改进了在图像上的修正的可见效果，并且也改进了用于处理的图像的质量，处理例如为用于对象表面重构的进一步的光度立体处理。用于明亮像素和暗像素的α的预定值可以分别被称作α_高和α_低。

基于对于参数α、k和V的上述区别，上述f(W)的通用公式可以以如下的详细方式写出：

图39、图65、图66和图67示出不同实施例，包括用于参数k、f(W)、V_高和V_低以及α的值的不同集合。

图65表示三条虚线，每条涉及f(W)的表示，其中k＝0.75且α＝100为固定值，并且其中V＝1、1.2和1.5，从而示出V对f(W)的影响。

图66表示三条虚线，每条涉及f(W)的表示，其中v＝1.25且α＝100为固定值，并且其中k＝1、0.75和0.5，从而示出k对f(W)的影响。

图39表示三条虚线，每条涉及f(W)的表示，其中V和k是定值，并且其中α＝10³、10²和10，从而示出α对f(W)的影响。

图67表示三条虚线，每条涉及f(W)的表示，其中V＝1.25且k＝0.75为固定值，并且其中α＝10³、10²和10，从而也示出α对f(W)的影响。

图40示出具有与明亮区域和暗区域关联的α和k的不同值的f(W)的表示，其中α_高＝100且α_低＝50以及k_高＝0.6且k_低＝0.8。

如图41所示，附加预处理的方法包括若干步骤1302、1304、1306以及1308，如下面所述。该方法包括步骤1302，用于图像捕捉设备108，来获取多个的N个图像。在所述实施例中，必须通过不同的照射源以不同的照射方向拍摄至少两个图像。

在步骤1302，强度检测设备112考虑每个不同图像1、2和3，并且检测用于每个色彩通道的每个像素的强度，来获取用于每个像素、用于每个图像并且用于每个色彩通道的像素值IR。

在步骤1304，数据处理器114为同一像素计算对于相同的色彩通道，在一图像中考虑的该像素的强度与在其他图像中的对应像素的强度之间的比值，来获得相应的比值。

在步骤1306，数据处理器114计算用于每个像素、用于每个图像并且用于每个色彩通道的每个计算得到的比值的相应平均值。

在步骤1308，数据处理器114基于上述公式计算用于每个图像中的每个像素且用于每个色彩通道的修正后的像素值Ic。

图42、图43、图44以及图45分别示出图31、图32、图33和图34的对象的修正后图像，其中每个图像中已经使用如用于比值估算的参考图像的三个其他图像来修正明亮区域。算术平均已经被用作平均值来用于计算W。图46、图47和图48分别在红色、绿色和蓝色通道中示出图42的修正后图像。

在图46、图47、图48示出的修正的例子中，上述参数k_高、k_低、V_高、V_低和α_高以及α_低可以如下示出的那样确定：

k_高＝k_低＝0.75

V_高＝1.2且V_低＝0.2

α_高＝100且α_低＝50

对于严格修正，参数V_高可以等于1.2。

对于软修正，参数V_高可以等于1.2并且参数α可以等于10，与k的不同值组合。当k＝1时，可见该修正为最大。

当k接近0时，没有提供修正。

参数k_高、k_低、V_高、V_低和α_高以及α_低的这些值基于计算得到的比值的算术平均而被经验性地确定。在所有本说明书的例子中，上面示出的参数集合被使用。

图49到图50以及图52到图54涉及存储在如在A.Hertzmann,S.M.Seitz,“Example-BasedPhotometricStereo:ShapeReconstructionwithGeneralVaryingBRDF”,IEEETransactionsonpatternanalysisandmachineintelligence,vol27,8,p.1254-1263,2005(A.Hertzmann，S.M.Seitz，“基于例子的光度立体：以一般性改变BRDF的形状重构，IEEE模式分析与机器智能学报，27卷，8，1254-1263页，2005)中提及的公共数据库中的图像。图49对应于在公共数据库中可得到的表示猫的对象的八个图像中的一个图像。图52、图53和图54对应于在公共数据库中可得到的表示鱼的对象的十四个图像中的三个图像。

图49示出具有多个明亮区域的表示猫的对象的图像。

图50示出位于猫的眼睛周围和其上的特定明亮区域。在以附加的预处理的方法处理在图49中示出的图像之后，得到在图51中所示的修正后图像。用于获得修正后图像的方法通过使用作为参考图像的该猫的其他七个图像来执行。算术平均被用于计算平均值W，并且参数k_高、k_低、V_高、V_低和α_高以及α_低的集合与上面提及的那些相同。视觉检测指示出现在图50中示出的图像的明亮区域没有出现在图51中所示的图像中。另外，作为应用根据本附加预处理的方法的结果，进一步的视觉检测指示图50的表示猫的眼睛尚未经修改的细节。

在图52、图53和图54示出的其他例子中，多个明亮区域和暗区域可以被观察到。如在图55、图56和图57中示出的，暗区域和明亮区域与已经通过应用根据本发明的方法去除。用于获取修正后图像的方法包括使用鱼的其他十三个图像作为参考图像，其中该十三个图像不同于要被修正的图像。算数平均被用于计算平均值W，并且参数k_高、k_高、V_高、V_高和α_高以及α_高的集合与上面提及的那些相同。另外，作为应用根据本附加预处理的方法的结果，视觉检测指示对象的细节没有被修改。

为了表示初始像素值、W值、f(W)值以及修正后像素值，基于图53和图54的图像的例子现在将被描述。例如为红色通道的一个色彩通道被考虑。图53和图54的图像的每个红色通道被分别表示在图58和图59中。图58和图59示出表示感兴趣的像素线的白线，用于示出W和f(W)的值。

图60示出表示用于在图58中示出的初始图像I1和在图59中示出的参考图像I2的像素线的像素值或强度值的曲线图。

图61示出使用图59上的参考图像来对图58中的初始图像I1进行修正的修正后图像I1c。使用在图59中示出的一个参考图像计算平均值W。因此，平均值等于在图58中示出的初始图像I1与在图59中示出的参考图像I2的比值。

图62示出表示平均值W和修正后的函数f(W)的曲线图。f(W)使用如上提及的相同参数k_高、k_低、V_高、V_高和α_高以及α_高进行计算。图62示出f(W)的值没有等于零，这指示要被修正的像素。这些值既对应于明亮区域，即，W的相对高值，也对应于暗区域，即，W的相对低值。

图63示出表示用于在图58中示出的初始图像I1和在图61中示出的修正后图像I1c的感兴趣像素线的像素值的曲线图。

图64示出表示平均值W和修正后的平均值Wc的曲线图，平均值W被估计为在图58中示出的初始图像I1和在图59中示出的参考图像I2的比值，修正后的平均值Wc被估计为在图61中示出的修正后图像I1c和在图59中示出的参考图像I2的比值。

如在图63中所示，在修正后图像I1c中的明亮区域的强度值低于在图60中示出的图像I1的明亮区域的强度值。以类似的方式，如在图63中示出的，在修正后图像I1c中的暗区域的强度值高于如在图60中示出的图像I1中的暗区域的强度值。另外，当与如在图60中示出的图像I1中的一些像素的强度值比较时，对于图63中的修正后图像I1c，表征对象纹理和表面的像素的强度值保持不变。

能注意到在初始图像I1和参考图像I2上都明亮或都暗的一些区域没有被检测到并且因此没有被修正。结果，这些区域的强度值分别异常地保持高或低。然而，使用更多带有不会彼此完全重叠的明亮区域或暗区域的参考图像优化了本发明的方法，其至少被引导为避免这样的漏检。

可替换地，代替于用于上述参数的特定值的确定，用于调适参数值的自动系统可以被实施。

每个图像上的镜面区域的位置对于改进以根据本发明的方法修正后的图像的质量非常重要。实际上，如已经提及的，每个镜面区域、即明亮区域或暗区域，不应该完全与另一个镜面区域重叠来提供改进质量的修正后图像。

因此，图87表示对象580的四个图像的叠加的示意图，其中图像1包括镜面区域590，图像2包括镜面区域591，图像3包括镜面区域592，并且图像4包括镜面区域593。没有镜面区域彼此重叠。根据本预设处理的该方法提供了如在图87中那样定位的镜面区域的最好的结果。

图88表示对象580的四个图像的叠加的示意图，其中图像1包括镜面区域601，图像2包括镜面区域602，图像3包括镜面区域603，并且图像4包括镜面区域604。每个镜面区域601、602、603和604至少部分地与另一个镜面区域重叠。

图89表示对象580的四个图像的叠加的示意图，其中图像1包括镜面区域611，图像2包括镜面区域612，图像3包括镜面区域613，并且图像4包括镜面区域614。至少一个镜面区域611完全与镜面区域612、613和614重叠。

图68表示在不同照射方向捕捉的对象的人工表面。

图69表示一个初始图像，其中明亮区域的修正是需要的。

图70表示对应于空间上均匀分布的七个照射源的七个图像的镜面分量的叠加形成的人工图像。

图71a表示使用从空间上均匀分布的七个照射源获取的七个参考图像修正的图69的修正后图像。

图71b表示沿着在图71a中的水平白线的像素强度的分布，对于散射分量以实线示出且对于基于软修正的结果分量以虚线示出。由于明亮区域均匀地定位，修正后图像的质量令人满意。这通过跟随图43b中的图形上的散射分量的结果分量示出。另外，不需要修正的区域没有特定影响。

图72表示由与空间上随机分布的七个照射源对应的七个图像的镜面分量的叠加形成的人工图像。

图73a表示使用从空间上随机分布的七个照射源获取的七个参考图像对图41修正的修正后图像。

图73b表示沿着在图73a中示出的的水平白线的像素强度的分布，对于散射分量以实线表示，对于基于软修正的结果分量以虚线表示。当照射源是随机定位时，修正后图像的质量也是令人满意的。

图74表示由与位置紧邻的七个照射源对应的七个图像的镜面分量的叠加形成的人工图像。

图75a表示使用从位置紧邻的七个照射源获取的七个参考图像修正图41的修正后图像；

图75b表示沿着在图75a中示出的水平白线的像素强度的分布，对于散射分量以实线表示，对于基于软修正的结果分量以虚线表示。由于在图74中示出的明亮区域重叠，所以修正处理可能没有完成。实际上，重叠的明亮区域可能作为具有错误强度的区域而未被检测到。

图76表示一个图像的镜面分量，其对应于与要被修正的镜面区域相对的照射源；

图77a表示使用从与要被修正的镜面区域相对的一个照射源获取的一个参考图像修正的图69的修正后图像。

图77b表示沿着在图77a中示出的水平白线的像素强度的分布，对于散射分量以实线表示，并且对于基于软修正的结果分量以虚线表示。尽管没有获取对象表面的光线的一般性分布，修正后图像仍然具有良好的质量。

图78表示由与空间上随机分布的九十九个照射源对应的九十九个图像的镜面分量的叠加而形成的人工图像。

图79a表示使用从空间上随机分布的九十九个照射源获取的九十九个参考图像修正的图69的修正后图像。

图79b表示沿着在图79a中示出的水平白线的像素强度的分布，对于散射分量以实线表示，并且对于基于软修正的结果分量以虚线表示。修正后图像具有良好的质量。

因此，当考虑到在空间上随机或均匀分布的相当大数量的图像时，获得用于图像修正的最好结果。

可见使用从相对的光源方向获得的两个初始图像产生明亮区域，明亮区域彼此远离。因此，之前限定的权重值W以最优化方式应用并且修正后图像具有改进的质量。

当使用在紧邻的光源方向获得的多个初始图像(例如一打图像)时，明亮区域的人工的和全部的或部分的重叠区域可以呈现。因此，如之前定义的权重值W可能不能正确地检测明亮区域或以有效的方式修正这样的明亮区域。

图80a、图80b、图80c和图80d表示圆形对象的四个初始图像，包括处于不同位置而没有任何重叠区域的明亮区域。

图81a、图81b、图81c和图81d表示圆形对象的图像，其中明亮区域已经使用现有技术中公知的修复处理进行修正。

图82a、图82b、图82c和图82d表示经过应用使用软修正的本发明方法之后的圆形对象的修正后图像。可见当与以公知修复技术获得的相应修正后图像81a、81b、81c和81d相比较，图82a、图82b、图82c和图82d的修正后图像具有改进的质量。

可以提供示意性比较来比较用于去除包括若干圆形线的对象的图像上的镜面区域的现有技术。

图83示出经过无需修正明亮区域的表面重构处理之后的圆形对象的表面的等高线的示意图。其示出一些圆形线有拉伸。

图84示出经过以现有技术中公知的修复处理重构之后的圆形对象的表面的等高线的示意图。可见这些圆形线仍然被拉伸。修正后图像具有很差的质量。

图85示出根据使用现有技术公知的科尔曼(Coleman)和简(Jain)算法的圆形对象的表面的等高线的示意图。圆形线大幅度地被拉伸和修改。修正后图像具有很差的质量。

图86示出在经过使用本发明方法获得的重构处理之后的对象的等高线的示意图，其中参数V_高＝1.2、k_高＝k_低＝0.5并且α_高＝α_低＝10。

如图86中所示，圆形线没有被拉伸并且每个圆形线之间的距离保持恒定。根据附加的预处理的用于去除镜面区域的方法提供了最好的结果。

第一实施例包括：

$>W=\underset{j=1}{\overset{p}{\Sigma }}{c}_j\frac{1}{I_{refj}}$>

以及

$>f\left(W_{im\_number,p\_number}\right)=\frac{k}{1+e^{-\alpha *(W_{im\_number,p\_number}-V)}}$>

算法内进行如下运算：

算法1：软镜面修正算法。

1：输入

I-要被修正的图像；

M-在图像I中的像素的数量；

{I_ref}-参考图像的序号；

P-参考图像的数量；

K、α、V-S型函数参数。

2：输出：

Ic-具有修正后的镜面区域的图像。

3：对于p_number＝1,…,M进行

4：对于l＝1，…，p进行

5：

6：结束

7：WR(p_number)←AggregRatio(IR(p_number))

8：

9：I_c(p_number)←I(p_number)/W(p_number)f(p_number)

10：结束

上述的第一实施例提供了在其中用户可以基于参考图像的集合为每个参数k、V和α定义调适后的设定的情况下具有提高质量的结果。

在第二实施例中，权重参数W被计算如下：

$>W_{p\_number}=\underset{l=1}{\overset{p}{\Sigma }}{c}_l(\frac{I(p\_number)}{I_{ref}^l(p\_number)}-LinApprox(I,I_{ref}^l)+1),$>

基于在第二实施例中计算的W，现在f(W)仅基于参数α的值涉及可调适的修正。

I涉及要被修正的图像的像素强度。在第二实施例中，图像的数量im_number＝1，并且因此可以不需考虑进一步的标注。当图像被展开、即向量化时，图像I中的每个像素被称作I(p_number)。另外，在第二实施例的本描述中，图像I和用于修正的图像被假定为灰阶的，即单通道图像。因此，不需要使用IR、IG和IB标注。图像I是在算法中使用的唯一可变表示。P个参考图像的集合也被用于修正从相同的固定相机位置和在不同的光源位置处拍摄的图像I、图像{I_ref}。例如，当图像的初始集合包括3个图像I₁、I₂、I₃时，图像I＝I₁的第一修正与参数图像的集合{I_ref}＝{I₂，I₃}关联。以类似的方式，图像I＝I₂的第二修正与参数图像的集合{I_ref}＝{I₁，I₃}关联。以类似的方式，图像I＝I₃的第三修正与参数图像的集合{I_ref}＝{I₁，I₂}关联。参考图像可以相当于初始序列除了要被修正的图像之外的所有其他图像。用于每个图像I的一般性光源分布或位置可以被视为低频分量并且大致接近强度的附加平面。该平面的取向和对应转换依赖于光源的位置和强度。如果图像I具有基本的阴影，那么每个图像的相关阴影可以被表示为比值的矩阵的线性近似，称为如下面详述的LinApprox，其中涉及点到点除法。因此，当从IR_ratio，的赋值减去近似平面时，一般性阴影被去除。当没有镜面区域时，该减法的值趋向于零。

通过使用线性或鲁棒的回归系数的平面，是比值的矩阵的近似。C₁表示用于广义平均值的图像的权重。线性近似的系数可以通过最小化下列的误差得到：

$>\underset{a_0,a_1,a_2}{\min }\left|\right|\frac{I(p\_number)}{I_{ref}^l(p\_number)}-(a_0+a_{1}x+a_{2}y)|{|}_2^2.$>

其中a₀、a₁、a₂是定义近似平面的等式的系数。

将比值和近似平面之间的残差(redisual)加一提供了以W作为简单比值进行运算的可能性。W(p_number)的值越接近1，I(p_number)和之间的差值越小。

在该第二实施例中，赋予可调适比值消除了使用阈值V的需要。获得的近似平面、，可以被视为接近每个像素(p_number)并用于每个参考图像I的阈值的矩阵。修正后图像像素I_c(p_number)可以被赋值为：

$>I_c(p\_number)=\frac{I(p\_number)}{W{(p\_number)}^{g(p\_number)}},$>

其中k是修正的程度，并且α定义修正后的区域和没有修正的区域之间的过渡的锐度。参数K和α对于软镜面去除算法具有相同的意思而不需要可调适的阈值。

第二实施例以相应的第二算法运算如下：

算法2：可调适软镜面修正算法。

1：输入

I-要被修正的图像；

M-在图像I中的像素的数量；

{I_ref}-参考图像的序号；

P-参考图像的数量；

K、α-S型函数参数

2：输出：

Ic-具有修正后的镜面区域的图像。

3：对于l＝1，…，p进行

4：对于p_number＝1,…,M进行

5：

6：结束

7：S^l←LinApprox(IR^l)

8：K^l←IR^l-S^l+1

9：结束

10：对于p_number＝1,…,M进行

11：WR(p_number)←AggregRatio(K(p_number))

12：

13：I_c(p_number)←I(p_number)/W(p_number)g(p_number)

14：结束

算法1和算法2，涉及分别用于提供对象修正后图像的第一实施例和第二实施例，具有相同的修正原理。然而，在算法1中，预定义的阈值V被预先包括到修正函数赋值(第8行，算法1)。对于算法2，该值不再作为比值IR的矩阵线性近似(第7行，算法2)以及从比值的初始矩阵减去(第8行，算法2)的运算结果而被减去(第12行，算法2)。这类似于对于每个图像像素具有合适的阈值而不是具有用于图像的所有像素的固定值。由于识别用于线性近似S(第7行，算法2)的比值IR的矩阵的全部值的必要性，图像修正被划分为两个循环(第3-9行，算法2以及第10-14行，算法2)。所有其他用于算法1的运算保持像素相关性(pixel-related)。

第二实施例提供了具有提高质量的结果，其中要被修正的图像涉及的照射条件显著不同于参考图像的照射条件。

因此，具有高强度的像素值为像素修正所考虑，并且比具有产生高于Vm_高的权重值W的像素值修正更少，高强度的像素值可以由例如因对象的特定性质所致的光线效应产生。

因此，具有低强度的像素值为像素修正所考虑，并且比具有产生高于Vm_高的权重值W的像素值修正更少，低强度的像素值可以由例如因对象或皮氏培养皿的特定性质所致的光线效应产生。

本发明提供了用于替换被处理的图像的所有像素值的方案。本发明提供了包括例如为镜面像素的明亮像素和例如为阴影像素的暗像素的特定区域或图像局部的像素值的同时修正。本发明提供了一种用于在这样的区域修正对应的像素值而不需要关于照射方向和/或照射源的位置的预定信息的方案。本发明提供了一种基于图像的多样性、通过循序考虑所有图像中的所有像素的方法。这意味着来自所有图像的像素信息被用来提供对于与镜面像素或阴影像素关联的任何像素的修正。本发明允许包括与图像关联的关于特定照射方向的特定信息的特定区域的像素值保持不变。实际上，这样的特定区域的每个像素值被处理并被替换为相同的像素值，因为不同的像素值对于该像素值来说是不需要的。特定区域的像素值可以涉及相关图像的阴影信息。因此，在本发明中，图像中的所有像素被考虑并且被合并来提供替换，还有像素值的修正。这方面允许对象的表面的取向在处理光度立体法时例如在该方法的进一步的步骤中被可靠地估计。由于本发明提供了基于初始图像的多样性的修正后图像的多样性，因此进一步的光度立体法的处理也得到改进。在特定的实施例中，如果不需要进一步的表面重构处理，最少只需要对象的两个图像来修正对象的仅有一个的图像的像素值。

另外，本预处理适于恢复图像的高光和阴影区域，即使当考虑到具有实质大小并包括有复杂的像素纹理的区域时。

为了最优化本预处理的方案，很多个图像可以被使用，并且很多个照射源可以均匀位于对象的周围。因此，与对象的像素值有关的增长的信息量可以被考虑到。结果，像素值的修正可以得到改进。此外，为了最优化像素值的修正，感兴趣区域必须不能与在所有其他图像中的另一个感兴趣区域完全重叠。

本预处理的方法以最优化方式应用到与镜面半球形对象相关的图像，例如微生物的图像。

参考文献

[1]J.C.BEZDEK,R.J.HATHAWAY,SomeNotesonAlternatingOptimization,Proceedingsofthe2002AFSSInternationalConferenceonFuzzySystems,2002,pp.288–300(对交替最优化的几点注意，2002年AFSS模糊系统国际会议论文集，2002年，288-300页).

[2]I.CSISZ′AR,G.TUSN′ADY,InformationGeometryandAlternatingMinimizationProcedures,Statistics&Decisions,SupplementIssue1(1984)205–237(信息几何和交替最小化过程，统计与决策，增刊第1期(1984)，205-237页).

[3]R.J.HATHAWAY,Y.HU,J.C.BEZDEK,LocalConvergenceofTri-levelAlternatingOptimization,NeuralParallel&ScientificComp.9(2001)19–28(三水平交替最优化的局部收敛，神经并行与科学计算，9(2001)19-28页).

[4]R.BASRI,D.W.JACOBS,I.KEMELMACHER,PhotometricStereowithGeneral,UnknownLighting,InternationalJournalofComputerVision72(3)(2007)239–257(具有通用的、未知的照明的光度立体，计算机视觉国际杂志72(3)(2007)239-257页).

[5]R.T.FRANKOT,R.CHELLAPPA,AMethodforEnforcingIntegrabilityinShapefromShadingAlgorithms,IEEETransactionsonPatternAnalysisandMachineIntelligence10(4)(1988)439–451(以阴影算法施行形状积分的方法，IEEE模式分析与机器智能学报10(4)(1988)439-451页).

[6]E.COLEMAN,R.JAIN,Obtaining3-DimensionalShapeofTexturedandSpecularSurfacesUsingFour-SourcePhotometry,ComputerGraphicsandImageProcessing18(1982)309–328(采用四源测光获得有质感和镜面的表面的3维形状，计算机图形和图像处理，18(1982)309-328页).

[7]P.KOVESI,ShapeletsCorrelatedwithSurfaceNormalsProduceSurfaces,ProceedingsoftheInternationalConferenceonComputerVision,2005,pp.994–1001(与表面法线生成表面相关的Shapelets，计算机视觉国际会议论文集，2005年，994-1001页).

[8]S.BOYD,L.VANDENBERGHE,ConvexOptimization,CambridgeUniversityPress,NewYork,NY,USA,2004(凸优化，剑桥大学出版社，纽约，NY，美国，2004年).

[9]R.T.ROCKAFELLAR,LagrangeMultipliersandOptimality,SIAMRev.35(2)(1993)183–238(拉格朗日乘子和最优化，SIAM评论，35(2)(1993)183-238页).

[10]K.LEVENBERG,AMethodfortheSolutionofCertainProblemsinLeastSquares,Quart.AppliedMath.2(1944)164–168(用于某些问题的最小二乘解决方案的方法Quart.应用数学，2(1944)164-168页).

[11]R.GROSSE,M.K.JOHNSON,E.H.ADELSON,W.T.FREEMAN,GroundtruthDatasetandBaselineEvaluationsforIntrinsicImageAlgorithms,InternationalConferenceonComputerVision,2009,pp.2335–2342(用于内在图像算法的参考标准数据集和基线评估，计算机视觉国际会议，2009年，2335–2342页).

[12]Woodham,R.J.(1980):“Photometricmethodfordeterminingsurfaceorientationfrommultipleimages”,Opticalengineering19,pp.513-531(“用于从多个图像确定表面取向的光度方法”，光学工程19，513-531页).

[13]Hayakawa,H.(1994):“Photometricstereounderalightsourcewitharbitrarymotion”,SocietyofAmerica11,pp.3079-3088(“在任意运动的光源下的光度立体”，美国协会11，3079-3088页).