CN102607455A - 基于光学显微镜和变化光照的微观形貌三维测量方法 - Google Patents

Abstract

本发明公开了一种基于光学显微镜和变换光照的微观形貌三维测量方法。该方法根据光学显微镜在固定视点不同光照方向下拍摄得到的多张微观图像重建微观结构表面的三维形貌。其过程如下：首先采用UPS方法和表面可积性约束获得具有GBR歧义的表面反射率和法线方向；然后用基于最小熵的GBR消歧法得到无歧义的表面反射率和法线方向；为了减少噪声的影响，在马尔可夫随机场模型下，利用图割法进一步优化表面法线方向；最后，根据表面法线方向，采用积分的方法重建微观结构表面的三维形貌。

Description

基于光学显微镜和变化光照的微观形貌三维测量方法

技术领域

本发明属于精密测量领域，具体涉及一种基于光学显微镜和变化光照的微观形貌三维测量方法。本发明将所提方法称为MPS方法。

背景技术

微观形貌三维测量是一个非常重要的研究领域，目前国内外的研究机构及企业已开发出多种原理的形貌形貌三维测量设备，使用比较广泛的有聚焦显微镜(Confocal microscopy)，白光干涉仪(white-light interferometry)，显微光栅投影(microscopic fringe projection)及显微立体视觉(lightmicroscope stereo vision)等。其中，聚焦显微镜和白光干涉仪可以测量纳米尺度的三维数据；而显微光栅投影和光学显微镜能够快速、稳定地测量范围在1cm²-几个cm²内的微观物体。上述设备大多结构复杂且价格昂贵，不适合用于很多对精度要求不高、仅需要微观三维形貌的领域(如：缺陷检测、皮肤医学、金属腐蚀学等)。因此，开发一种易于实施且成本低廉的微观形貌测量方法具有非常明确的现实应用需求。

光度立体视觉技术(Photometric Stereo)是宏观领域里常用的一种光学三维测量方法，其基本原理是在不同的光照条件下，从同一视点拍摄多张图像，根据每个像素的强度变化计算物体表面的法线方向，最后根据法线方向计算出每个像素对应的三维坐标，从而得到被测物体表面完整密集的三维形貌。由于此方法非常简单且易于实施，目前已广泛应用于宏观尺度上的形状测量、缺陷检测以及物体分类上。但是，目前这一方法还没有被用于微观领域。本发明将光度立体法引入微观领域，提出一种新的基于光学显微镜和变化光照的微观形貌三维测量方法：Microscopic PhotometricStereo(MPS)。

现有的光度立体视觉方法可分为两类：已标定的光度立体视觉方法(calibrated photometric stereo，CPS)和未标定的光度立体视觉方法(uncalibrated photometric stereo，UPS)。在CPS中，应首先使用金属球和漫反射球标定光源方向和强度。使用这些额外的标定对象进行标定的过程较为复杂，且由于光学显微镜的视场很小，额外的标定对象会使有效测量区域变小，从而影响微观形貌测量。因此，不需要额外标定对象的UPS方法更适用于MPS。在UPS中，由于光照条件是未知的，会使得计算得到的法线方向存在歧义，因此，必须采用适当的约束来消除这种UPS中固有的歧义问题。在UPS中，表面法线和反射率可以利用奇异值分解(singular valuedecomposition，SVD)的方法恢复为一个有9个参数的线性变换，然后根据表面可积的约束性可使这个9个参数的线性变换减少为3个参数的广义浅浮雕(generalized bas-relief，GBR)模糊问题[文献：P.N.Belhumeur，D.J.Kriegman，A.L.Yuille.The Bas-Relief Ambiguity.International Journal ofComputer Vision，1999，35(1)：33-44.]。GBR模糊问题是UPS中的重要问题，产生这一问题的原因在于不同的表面形貌和光源方向可能产生同样的灰度图像，因此当光源方向未知时，必须引入其它的约束来解决这一问题。比如：当光源强度和/或相对表面反射率已知时，可以消除GBR模糊；也可通过非朗伯反射模型，如托伦斯-斯帕罗模型或朗伯加镜面反散模型等，来解决。

现有的方法并不能直接适用于MPS，因为在MPS中，当光源任意移动时显微镜拍摄得到的多张显微图像的光强变化非常显著，且图像噪声非常严重。前一现象使得光源强度一致的假设变得不再可用；后一现象使得采用复杂的非朗伯反射模型时会产生严重的数值误差。因此，现有的依靠光源强度约束或非朗伯模型的GBR消歧法不适用于MPS。

“Fast approximate energy minimization via graph cuts”(Boykov Y，Veksler O and Zabih R.，IEEE Trans.Pattern Anal.Mach.Intell.2011，23：1222-1239)中公开了一种图割法，图割法是一种能量最小化的优化方法，被广泛应用于一阶马尔可夫随机场模型的计算和求解。

发明内容

本发明的目的在于提供一种基于光学显微镜和变换光照的微观形貌三维测量方法，此方法能够将普通的光学显微镜转化为一种微观形貌三维测量工具。

本发明提供的一种基于光学显微镜和变换光照的微观形貌三维测量方法，其特征在于，该方法包括下述步骤：

第1步在固定视点不同光照方向下拍摄多幅被测微观物体的微观图像；

第2步采用UPS方法和表面可积性约束获得具有GBR歧义的表面反射率和法线方向；

第3步用基于最小熵的GBR消歧法得到无歧义的表面反射率和法线方向；

第4步在马尔可夫随机场模型下，利用图割法优化表面法线方向；

第5步根据表面法线方向，采用积分的方法重建微观结构表面的三维形貌。

作为上述技术方案的改进，第3步优选的实现过程为：：

将m个像素作为采样点，从m个采样点中生成一个h-bin直方图{a_i}，i＝1，2，...，h，h表示直方图中的采样数目，h为小于等于m的正整数；a_i表示数值分布在第i个区间内的采样点数目，a_i≤m；采用最大释然估计算子，得到最小化目标函数

$\underset{u,v,\mathrm{\λ}}{\min }O(u,v,\mathrm{\λ})=-\underset{i=1}{\overset{h}{\mathrm{\Σ}}}\frac{a_i}{m}\log \frac{a_i}{m}$

采用模拟退火法进行最小化优化，确定参数u，v，λ；

为伪法线方向，

为伪光源强度与方向，根据

计算出无歧义的表面反射率表面法线方向

和光源强度与方向

其中，

$G=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ u& v& \mathrm{\λ}\end{array}\right),$ $G^{-1}=\frac{1}{\mathrm{\λ}}\left(\begin{array}{ccc}\mathrm{\λ}& 0& 0\\ 0& \mathrm{\λ}& 0\\ -u& -v& 1\end{array}\right).$

作为上述技术方案的进一步改进，第4步优选的实现过程为：

设(p，q)为一对相邻像素，X＝N，其中N＝[n₁，n₂，...，n_m]为被测对象各像素的法线方向；Y＝I，其中I＝[I₁，I₂，...，I_f]为拍摄的一组图像，则光度立体视觉技术的马尔可夫随机场模型表示为：

式I

其中，

为图像中像素p的匹配代价函数，为邻域内法线方向的兼容性函数，最佳概率根据上述式I的最大后验概率进行计算：

$E\left(N\right)=\underset{p}{\mathrm{\Σ}}D(n_p,I_p)+\underset{(p,q)}{\mathrm{\Σ}}V(n_p,n_q)=E_{\mathrm{data}}\left(N\right)+E_{\mathrm{smoothness}}\left(N\right)$ 式II

其中，函数D和V为计划使用图割法进行最小化的能量函数，分别代表图割法中的数据单元和平滑单元，数据单元D对应式I中的匹配代价函数

平滑单元V对应式I中的兼容性函数；

其中， $E_{\mathrm{data}}\left(N\right)=\underset{p\∈P}{\mathrm{\Σ}}\left|\right|{\stackrel{\→}{b}}_{{\mathrm{\α}}_p}-{\hat{b}}_p\left|\right|,$ $E_{\mathrm{smoothness}}\left(N\right)=t\underset{(p,q)\∈P}{\mathrm{\Σ}}\log (1+\frac{\left|\right|{\stackrel{\→}{b}}_{{\mathrm{\α}}_p}-{\stackrel{\→}{b}}_{{\mathrm{\α}}_q}\left|\right|}{{2\mathrm{\σ}}^2})$

其中，

为图割法优化过程中得到的像素p的法向，α_p为像素p处的α扩张系数，为图割法优化过程中得到的像素q的法向，α_q为像素q处的α扩张系数，t和σ均为0到1之间的小数，对能量函数进行最小化运算，得到所需的表面法线方向n。

本发明根据光学显微镜在固定视点不同光照方向下拍摄得到的多张微观图像重建微观结构表面的三维形貌。其过程如下：首先采用UPS方法和表面可积性约束获得具有GBR歧义的表面反射率和法线方向；然后用基于最小熵的GBR消歧法优化表面反射率和法线方向；为了减少噪声的影响，并进一步优化获得的表面法线方向，在马尔可夫随机场(Markov randomfield，MRF)模型下，利用图割法优化表面法线方向；最后，根据表面法线方向，采用积分的方法重建微观结构表面的三维形貌。

附图说明

图1是MPS方法原理图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

下面对本发明做进一步详细的说明。

(1)使用光学显微镜在固定视点不同光照方向下拍摄多幅被测微观物体的微观图像。

(2)利用未标定的光度立体视觉方法(UPS)获取具有GBR歧义的表面反射率与表面法线方向。

假设在同一视点、不同的光照条件下获取f幅图像，每幅图像具有m个像素。根据朗伯体反射模型，在第j幅图像的第i个像素处的图像灰度为：

I_ij＝ρ_in_i ^Ts_j      (1)

其中ρ_i和n_i为第i个像素处的表面反射率和法线方向，s_j包含了拍摄第j幅图像时光源的强度和方向。为了简化分析，将公式(1)表达为矩阵形式，设I为一个m×f的矩阵，其中m表示每幅图像的像素数，f表示图像的数目(要求f＞3)。因此，I中每一列代表一幅图像。

I＝BS        (2)

其中I(m×f)存储所有图像像素的灰度值；B(m×3)包含了每个像素处物体表面的反射率ρ_i和法线方向n_i，第i行为第i个像素处物体表面反射率与表面法线方向的乘积，即：b_i＝ρ_in_i；S(3×f)包含了所有图像拍摄时的光源强度和方向，第j列为拍摄第j幅图像时的光源s_j。

如果光源强度和方向已知，可以利用已标定的光度立体视觉方法，直接使用最小二乘法由公式(2)计算物体表面的反射率和法线方向B。然而，当光源强度和方向未知时，通过SVD分解(即奇异值分解)方法，只能获取具有3×3线性变换的B和S：

${\mathrm{UDV}}^T=I,\tilde{B}={\mathrm{EP}}_3,\tilde{S}=Q_{3}F.---\left(3\right)$

其中，U和V代表二个相互正交矩阵，D代表一对角矩阵，

为伪法线方向(pseudo-normals)，为伪光源强度与方向(pseudo-lights)；P₃和Q₃均为3×3矩阵，他们满足P₃Q₃＝D₃，D₃为3×3的对角矩阵，其对角元素为I^T的最大三个特征值的平方根；E是V的前三列，F是U的前三列。P₃和Q₃存在歧义：

其中A为任意未知的可逆矩阵。这种歧义是公式(2)所固有的，必须进行消除。现有的研究已证明，借助表面可积性约束，可以将这种九个参数的线性模糊(或线性歧义)A转变为三个参数的广义浅浮雕(generalized bas relief，GBR)变换。

GBR变换可表示为如下形式：

$G=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ u& v& \mathrm{\λ}\end{array}\right),$ $G^{-1}=\frac{1}{\mathrm{\λ}}\left(\begin{array}{ccc}\mathrm{\λ}& 0& 0\\ 0& \mathrm{\λ}& 0\\ -u& -v& 1\end{array}\right).---\left(4\right)$

u，v，λ为GBR变换的三个参数(这三个参数没有特定含义，用于表明G是一个下三角矩阵)。假设一个朗伯表面的高度为z＝f(x，y)(其中x，y为某像素的二维坐标)，表面反射率和法线方向分别为ρ(x，y)和n(x，y)，GBR变换对他们具有如下影响：

$\tilde{z}=\mathrm{\λz}+\mathrm{ux}+\mathrm{vy},$

$\widetilde{\mathrm{\ρ}}=\mathrm{\ρ}\left|\right|n^T{G}^{-1}\left|\right|,---\left(5\right)$

$\tilde{n}=n^T{G}^{-1}/\left|\right|n^T{G}^{-1}\left|\right|.$

其中，

和

为变换后的高度、表面反射率和表面法线方向。很明显，这种GBR变换会影响后续三维重建的精度，因此在进行三维重建之前，必须对这种GBR变换进行消除，即：GBR消歧。

(3)GBR消歧

本发明所提的GBR消歧法基于如下观察：绝大部分人造或自然物体，其表面反射率是由少量的占用主导地位的反射率组成的，GBR变换将使得这类物体的表面反射率的概率密度函数(Probability Density Function，PDF)变得更加平滑，从而增加表面反射率的熵。因此，本发明将GBR消歧问题转化为表面反射率的熵的最小化问题，通过获取能够具有最小熵的u，v，λ，来实现GBR消歧。

设

为使用UPS方法得到的GBR变换后的表面反射率与表面法线方向的乘积，

为待求的真实的表面反射率。很明显，ρ_i是GBR参数u，v，λ的函数。本发明采用直方图近似计算待测物体表面反射率的概率密度函数，假设从m个采样点(每个像素即为一个采样点)中生成一个h-bin直方图{a_i}，i＝1，2，...，h，h表示直方图中的采样数目，h为小于等于m的正整数；a_i表示数值分布在第i个区间内的采用点数目，a_i≤m。采用最大释然估计算子，能够得到如下最小化目标函数：

$\underset{u,v,\mathrm{\λ}}{\min }O(u,v,\mathrm{\λ})=-\underset{i=1}{\overset{h}{\mathrm{\Σ}}}\frac{a_i}{m}\log \frac{a_i}{m}---\left(6\right)$

这是一个典型的多参数优化问题，本发明采用模拟退火法进行最小化优化。一旦u，v，λ即G被确定，即可根据

计算出无歧义的表面反射率表面法线方向

和光源强度与方向

但是在实际测量过程中，由于拍摄得到的微观形貌的图像中存在严重的图像噪声，采用上述方法很难得到非常准确的GBR参数，计算的表面法线方向存在误差，会影响后续三维重建的精度。因此，需要进一步优化表面法线方向。

(4)表面法线方向优化

本发明在推导光度立体视觉技术的马尔可夫随机场模型(MarkovRandom Field，MRF)的基础上，采用图割法对表面法线方向进行优化。

马尔可夫网络(Markov network)包括两种节点：X和Y，X为隐藏变量集，Y为观测变量集。其后验概率分布P(X|Y)可表示为：

其中代表局部信息，φ(x_p，x_q)表示兼容性函数，(p，q)为一对相邻像素。设：X＝N，其中N＝[n₁，n₂，...，n_m]为被测对象各像素的法线方向；Y＝I，其中I＝[I₁，I₂，...，I_f]为拍摄的一组图像。则光度立体视觉技术的马尔可夫随机场模型可表示为：

其中，为图像中像素p的匹配代价函数(matching cost function)，为邻域内法线方向的兼容性函数(compatibility function)。最佳概率根据上述公式的最大后验概率进行计算：

其中，函数D和V为计划使用图割法进行最小化的能量函数，分别代表图割法中的数据单元和平滑单元，其中数据单元D对应公式(8)中的匹配代价函数

平滑单元V对应公式(8)中的兼容性函数。E(N)是表示根据最大后验概率计算得到的最佳概率。本发明中，log均是指以e为底的对数。

本发明构建的数据单元为：

$E_{\mathrm{data}}\left(N\right)=\underset{p\∈P}{\mathrm{\Σ}}\left|\right|I_p-{\widehat{\mathrm{\ρ}}}_p{\hat{n}}_p{\hat{s}}_p\left|\right|---\left(10\right)$

如果保持光源强度和方向

不变，则图像中某像素p的灰度主要取决于该像素的表面发射率

和法线方向

的乘积(即非单位化的法向

因此，能够使用非单位化的法向

与图割法优化过程中得到的法向

之间的差值来替代公式(10)，即：数据单元可表示为：

$E_{\mathrm{data}}\left(N\right)=\underset{p\∈P}{\mathrm{\Σ}}\left|\right|{\stackrel{\→}{b}}_{{\mathrm{\α}}_p}-{\hat{b}}_p\left|\right|---\left(11\right)$

其中为图割法优化过程中得到的像素p的法向，α_p为像素p处的α扩张系数，可根据图割法原理进行确定。

本发明构建的平滑单元为：

$E_{\mathrm{smoothness}}\left(N\right)=t\underset{(p,q)\∈P}{\mathrm{\Σ}}{\log }_e(1+\frac{\left|\right|{\stackrel{\→}{b}}_{{\mathrm{\α}}_p}-{\stackrel{\→}{b}}_{{\mathrm{\α}}_q}\left|\right|}{{2\mathrm{\σ}}^2}),---\left(12\right)$

其中

为图割法优化过程中得到的像素q的法向，α_q为像素q处的α扩张系数，t是正则化系数，σ是一个用来控制洛仑兹函数形状的参数。t和σ的数值根据被测物体的不连续性进行设置，t和σ均为0到1之间的小数，被测物体越平滑t和σ的取值越小。

构建出上述能量函数后，即可使用Kolmogorov开发的程序进行最小化计算，从而得到更优的表面法线方向n。

(5)表面三维重建

使用上述过程计算得到物体表面法线方向n后，即可使用现有的多种方法重建出微观形貌表面的三维数据。

设某点的高度值为z(x，y)，其法线方向为n(x，y)＝(r，t，-1)^T，其中r＝z_x，t＝z_y，z_x，z_y为z在x方向和y方向的偏微分，可表示为：

z_x＝z(x+1，y)-z(x，y)＝-n_x/n_z，

                                                   (13)

z_y＝z(x，y+1)-z(x，y)＝-n_y/n_z.

其中(n_x，n_y，n_z)为计算得到的法线方向。上述方程为过约束的线性方程，可以使用最小二乘法进行求解，得到被测物体的高度值z(x，y)。

上述具体实施方式的内容只是本发明的一个具体实例，本发明的内容并不局限于上述实例的内容。本领域一般技术人员根据本发明公开的内容，可以采用其它多种具体实施方式实施本发明，因此，凡是采用本发明的设计结构和思路，做一些简单的变化或更改的设计，都落入本发明保护的范围。

Claims (3)

1.一种基于光学显微镜和变换光照的微观形貌三维测量方法，其特征在于，该方法包括下述步骤：

第1步在固定视点不同光照方向下拍摄多幅被测微观物体的微观图像；

第2步采用UPS方法和表面可积性约束获得具有GBR歧义的表面反射率和法线方向；

第3步用基于最小熵的GBR消歧法得到无歧义的表面反射率和法线方向；

第4步在马尔可夫随机场模型下，利用图割法优化表面法线方向；

第5步根据表面法线方向，采用积分的方法重建微观结构表面的三维形貌。

2.根据权利要求1所述的微观形貌三维测量方法，其特征在于，第3步具体包括下述步骤：

将m个像素作为采样点，从m个采样点中生成一个h-bin直方图{a_i}，i＝1，2，...，h，h表示直方图中的采样数目，h为小于等于m的正整数；a_i表示数值分布在第i个区间内的采样点数目，a_i≤m；采用最大释然估计算子，得到最小化目标函数

$\underset{u,v,\mathrm{\λ}}{\min }O(u,v,\mathrm{\λ})=-\underset{i=1}{\overset{h}{\mathrm{\Σ}}}\frac{a_i}{m}\log \frac{a_i}{m}$

采用模拟退火法进行最小化优化，确定参数u，v，λ；

为伪法线方向，为伪光源强度与方向，根据

计算出无歧义的表面反射率

表面法线方向和光源强度与方向

其中，

$G=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ u& v& \mathrm{\λ}\end{array}\right),$ $G^{-1}=\frac{1}{\mathrm{\λ}}\left(\begin{array}{ccc}\mathrm{\λ}& 0& 0\\ 0& \mathrm{\λ}& 0\\ -u& -v& 1\end{array}\right)$

。

3.根据权利要求1所述的微观形貌三维测量方法，其特征在于，第4步具体包括下述步骤：

设(p，q)为一对相邻像素，X＝N，其中N＝[n₁，n₂，...，n_m]为被测对象各像素的法线方向；Y＝I，其中I＝[I₁，I₂，...，I_f]为拍摄的一组图像，则光度立体视觉技术的马尔可夫随机场模型表示为：

式I

其中，

为图像中像素p的匹配代价函数，为邻域内法线方向的兼容性函数，最佳概率根据上述式I的最大后验概率进行计算：

$E\left(N\right)=\underset{p}{\mathrm{\Σ}}D(n_p,I_p)+\underset{(p,q)}{\mathrm{\Σ}}V(n_p,n_q)=E_{\mathrm{data}}\left(N\right)+E_{\mathrm{smoothness}}\left(N\right)$ 式II

其中，函数D和V为计划使用图割法进行最小化的能量函数，分别代表图割法中的数据单元和平滑单元，数据单元D对应式I中的匹配代价函数

平滑单元V对应式I中的兼容性函数；

其中， $E_{\mathrm{data}}\left(N\right)=\underset{p\∈P}{\mathrm{\Σ}}\left|\right|{\stackrel{\→}{b}}_{{\mathrm{\α}}_p}-{\hat{b}}_p\left|\right|,E_{\mathrm{smoothness}}\left(N\right)=t\underset{(p,q)\∈P}{\mathrm{\Σ}}\log (1+\frac{\left|\right|{\stackrel{\→}{b}}_{{\mathrm{\α}}_p}-{\stackrel{\→}{b}}_{{\mathrm{\α}}_q}\left|\right|}{{2\mathrm{\σ}}^2})$

其中，

为图割法优化过程中得到的像素p的法向，α_p为像素p处的α扩张系数，

为图割法优化过程中得到的像素q的法向，α_q为像素q处的α扩张系数，t和σ均为0到1之间的小数，对能量函数进行最小化运算，得到所需的表面法线方向n。

Images