CN102607455B - 基于光学显微镜和变化光照的微观形貌三维测量方法 - Google Patents
基于光学显微镜和变化光照的微观形貌三维测量方法 Download PDFInfo
- Publication number
- CN102607455B CN102607455B CN 201210045676 CN201210045676A CN102607455B CN 102607455 B CN102607455 B CN 102607455B CN 201210045676 CN201210045676 CN 201210045676 CN 201210045676 A CN201210045676 A CN 201210045676A CN 102607455 B CN102607455 B CN 102607455B
- Authority
- CN
- China
- Prior art keywords
- normal direction
- pixel
- lambda
- sigma
- gbr
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Landscapes
- Length Measuring Devices By Optical Means (AREA)
Abstract
本发明公开了一种基于光学显微镜和变换光照的微观形貌三维测量方法。该方法根据光学显微镜在固定视点不同光照方向下拍摄得到的多张微观图像重建微观结构表面的三维形貌。其过程如下:首先采用UPS方法和表面可积性约束获得具有GBR歧义的表面反射率和法线方向;然后用基于最小熵的GBR消歧法得到无歧义的表面反射率和法线方向;为了减少噪声的影响,在马尔可夫随机场模型下,利用图割法进一步优化表面法线方向;最后,根据表面法线方向,采用积分的方法重建微观结构表面的三维形貌。
Description
技术领域
本发明属于精密测量领域,具体涉及一种基于光学显微镜和变化光照的微观形貌三维测量方法。本发明将所提方法称为MPS方法。
背景技术
微观形貌三维测量是一个非常重要的研究领域,目前国内外的研究机构及企业已开发出多种原理的形貌形貌三维测量设备,使用比较广泛的有聚焦显微镜(Confocal microscopy),白光干涉仪(white-light interferometry),显微光栅投影(microscopic fringe projection)及显微立体视觉(lightmicroscope stereo vision)等。其中,聚焦显微镜和白光干涉仪可以测量纳米尺度的三维数据;而显微光栅投影和光学显微镜能够快速、稳定地测量范围在1cm2-几个cm2内的微观物体。上述设备大多结构复杂且价格昂贵,不适合用于很多对精度要求不高、仅需要微观三维形貌的领域(如:缺陷检测、皮肤医学、金属腐蚀学等)。因此,开发一种易于实施且成本低廉的微观形貌测量方法具有非常明确的现实应用需求。
光度立体视觉技术(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直方图{ai},i=1,2,...,h,h表示直方图中的采样数目,h为小于等于m的正整数;ai表示数值分布在第i个区间内的采样点数目,ai≤m;采用最大释然估计算子,得到最小化目标函数
采用模拟退火法进行最小化优化,确定参数u,v,λ;
作为上述技术方案的进一步改进,第4步优选的实现过程为:
设(p,q)为一对相邻像素,X=N,其中N=[n1,n2,...,nm]为被测对象各像素的法线方向;Y=I,其中I=[I1,I2,...,If]为拍摄的一组图像,则光度立体视觉技术的马尔可夫随机场模型表示为:
其中,
其中, 为图割法优化过程中得到的像素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个像素处的图像灰度为:
Iij=ρini Tsj (1)
其中ρi和ni为第i个像素处的表面反射率和法线方向,sj包含了拍摄第j幅图像时光源的强度和方向。为了简化分析,将公式(1)表达为矩阵形式,设I为一个m×f的矩阵,其中m表示每幅图像的像素数,f表示图像的数目(要求f>3)。因此,I中每一列代表一幅图像。
I=BS (2)
其中I(m×f)存储所有图像像素的灰度值;B(m×3)包含了每个像素处物体表面的反射率ρi和法线方向ni,第i行为第i个像素处物体表面反射率与表面法线方向的乘积,即:bi=ρini;S(3×f)包含了所有图像拍摄时的光源强度和方向,第j列为拍摄第j幅图像时的光源sj。
如果光源强度和方向已知,可以利用已标定的光度立体视觉方法,直接使用最小二乘法由公式(2)计算物体表面的反射率和法线方向B。然而,当光源强度和方向未知时,通过SVD分解(即奇异值分解)方法,只能获取具有3×3线性变换的B和S:
其中,U和V代表二个相互正交矩阵,D代表一对角矩阵,为伪法线方向(pseudo-normals),为伪光源强度与方向(pseudo-lights);P3和Q3均为3×3矩阵,他们满足P3Q3=D3,D3为3×3的对角矩阵,其对角元素为IT的最大三个特征值的平方根;E是V的前三列,F是U的前三列。P3和Q3存在歧义: 其中A为任意未知的可逆矩阵。这种歧义是公式(2)所固有的,必须进行消除。现有的研究已证明,借助表面可积性约束,可以将这种九个参数的线性模糊(或线性歧义)A转变为三个参数的广义浅浮雕(generalized bas relief,GBR)变换。
GBR变换可表示为如下形式:
u,v,λ为GBR变换的三个参数(这三个参数没有特定含义,用于表明G是一个下三角矩阵)。假设一个朗伯表面的高度为z=f(x,y)(其中x,y为某像素的二维坐标),表面反射率和法线方向分别为ρ(x,y)和n(x,y),GBR变换对他们具有如下影响:
(3)GBR消歧
本发明所提的GBR消歧法基于如下观察:绝大部分人造或自然物体,其表面反射率是由少量的占用主导地位的反射率组成的,GBR变换将使得这类物体的表面反射率的概率密度函数(Probability Density Function,PDF)变得更加平滑,从而增加表面反射率的熵。因此,本发明将GBR消歧问题转化为表面反射率的熵的最小化问题,通过获取能够具有最小熵的u,v,λ,来实现GBR消歧。
设为使用UPS方法得到的GBR变换后的表面反射率与表面法线方向的乘积,为待求的真实的表面反射率。很明显,ρi是GBR参数u,v,λ的函数。本发明采用直方图近似计算待测物体表面反射率的概率密度函数,假设从m个采样点(每个像素即为一个采样点)中生成一个h-bin直方图{ai},i=1,2,...,h,h表示直方图中的采样数目,h为小于等于m的正整数;ai表示数值分布在第i个区间内的采用点数目,ai≤m。采用最大释然估计算子,能够得到如下最小化目标函数:
但是在实际测量过程中,由于拍摄得到的微观形貌的图像中存在严重的图像噪声,采用上述方法很难得到非常准确的GBR参数,计算的表面法线方向存在误差,会影响后续三维重建的精度。因此,需要进一步优化表面法线方向。
(4)表面法线方向优化
本发明在推导光度立体视觉技术的马尔可夫随机场模型(MarkovRandom Field,MRF)的基础上,采用图割法对表面法线方向进行优化。
马尔可夫网络(Markov network)包括两种节点:X和Y,X为隐藏变量集,Y为观测变量集。其后验概率分布P(X|Y)可表示为:
其中代表局部信息,φ(xp,xq)表示兼容性函数,(p,q)为一对相邻像素。设:X=N,其中N=[n1,n2,...,nm]为被测对象各像素的法线方向;Y=I,其中I=[I1,I2,...,If]为拍摄的一组图像。则光度立体视觉技术的马尔可夫随机场模型可表示为:
其中,为图像中像素p的匹配代价函数(matching cost function),为邻域内法线方向的兼容性函数(compatibility function)。最佳概率根据上述公式的最大后验概率进行计算:
其中,函数D和V为计划使用图割法进行最小化的能量函数,分别代表图割法中的数据单元和平滑单元,其中数据单元D对应公式(8)中的匹配代价函数平滑单元V对应公式(8)中的兼容性函数。E(N)是表示根据最大后验概率计算得到的最佳概率。本发明中,log均是指以e为底的对数。
本发明构建的数据单元为:
如果保持光源强度和方向不变,则图像中某像素p的灰度主要取决于该像素的表面发射率和法线方向的乘积(即非单位化的法向 因此,能够使用非单位化的法向与图割法优化过程中得到的法向之间的差值来替代公式(10),即:数据单元可表示为:
本发明构建的平滑单元为:
其中为图割法优化过程中得到的像素q的法向,αq为像素q处的α扩张系数,t是正则化系数,σ是一个用来控制洛仑兹函数形状的参数。t和σ的数值根据被测物体的不连续性进行设置,t和σ均为0到1之间的小数,被测物体越平滑t和σ的取值越小。
构建出上述能量函数后,即可使用Kolmogorov开发的程序进行最小化计算,从而得到更优的表面法线方向n。
(5)表面三维重建
使用上述过程计算得到物体表面法线方向n后,即可使用现有的多种方法重建出微观形貌表面的三维数据。
设某点的高度值为z(x,y),其法线方向为n(x,y)=(r,t,-1)T,其中r=zx,t=zy,zx,zy为z在x方向和y方向的偏微分,可表示为:
zx=z(x+1,y)-z(x,y)=-nx/nz,
(13)
zy=z(x,y+1)-z(x,y)=-ny/nz.
其中(nx,ny,nz)为计算得到的法线方向。上述方程为过约束的线性方程,可以使用最小二乘法进行求解,得到被测物体的高度值z(x,y)。
上述具体实施方式的内容只是本发明的一个具体实例,本发明的内容并不局限于上述实例的内容。本领域一般技术人员根据本发明公开的内容,可以采用其它多种具体实施方式实施本发明,因此,凡是采用本发明的设计结构和思路,做一些简单的变化或更改的设计,都落入本发明保护的范围。
Claims (3)
1.一种基于光学显微镜和变换光照的微观形貌三维测量方法,其特征在于,该方法包括下述步骤:
第1步在固定视点不同光照方向下拍摄多幅被测微观物体的微观图像;
第2步采用UPS方法和表面可积性约束获得具有GBR歧义的表面反射率和法线方向;
第3步用基于最小熵的GBR消歧法得到无歧义的表面反射率和法线方向;
第4步在马尔可夫随机场模型下,利用图割法优化表面法线方向;
第5步根据表面法线方向,采用积分的方法重建微观结构表面的三维形貌。
3.根据权利要求1所述的微观形貌三维测量方法,其特征在于,第4步具体包括下述步骤:
设(p,q)为一对相邻像素,X=N,其中N=[n1,n2,...,nm]为被测对象各像素的法线方向;Y=I,其中I=[I1,I2,...,If]为拍摄的一组图像,则光度立体视觉技术的马尔可夫随机场模型表示为:
其中,
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN 201210045676 CN102607455B (zh) | 2012-02-27 | 2012-02-27 | 基于光学显微镜和变化光照的微观形貌三维测量方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN 201210045676 CN102607455B (zh) | 2012-02-27 | 2012-02-27 | 基于光学显微镜和变化光照的微观形貌三维测量方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN102607455A CN102607455A (zh) | 2012-07-25 |
CN102607455B true CN102607455B (zh) | 2013-12-25 |
Family
ID=46525062
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN 201210045676 Active CN102607455B (zh) | 2012-02-27 | 2012-02-27 | 基于光学显微镜和变化光照的微观形貌三维测量方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN102607455B (zh) |
Families Citing this family (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107677216B (zh) * | 2017-09-06 | 2019-10-29 | 西安交通大学 | 一种基于光度立体视觉的多个磨粒三维形貌同步获取方法 |
CN107767418A (zh) * | 2017-10-25 | 2018-03-06 | 梧州学院 | 一种低倍显微物镜识别与定位方法 |
CN108303045A (zh) * | 2018-02-01 | 2018-07-20 | 北京科技大学 | 一种表面粗糙度测量方法与装置 |
CN108460773B (zh) * | 2018-02-28 | 2022-03-18 | 哈尔滨工程大学 | 一种基于偏移场水平集的声纳图像分割方法 |
CN108458658A (zh) * | 2018-04-20 | 2018-08-28 | 南京航空航天大学 | 一种基于光照反射模型的微细孔三维形貌测量装置及方法 |
CN113318913A (zh) * | 2021-02-01 | 2021-08-31 | 北京理工大学 | 一种基于未标定光度立体视觉的胶点三维重建方法 |
CN114061488B (zh) * | 2021-11-15 | 2024-05-14 | 华中科技大学鄂州工业技术研究院 | 一种物体测量方法、***以及计算机可读存储介质 |
CN115185028A (zh) * | 2022-07-25 | 2022-10-14 | 深圳博升光电科技有限公司 | 零级光强可调的二维分束衍射光栅制备方法 |
CN117059512B (zh) * | 2023-10-13 | 2024-01-26 | 苏州瑞霏光电科技有限公司 | 一种晶圆的三维表面微观测量方法及*** |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JPS60200111A (ja) * | 1984-03-26 | 1985-10-09 | Hitachi Ltd | 3次元物体認識装置 |
CN1963390A (zh) * | 2005-11-08 | 2007-05-16 | 李啸炜 | 一种便捷高效的三维测量方法 |
CN100495441C (zh) * | 2007-07-25 | 2009-06-03 | 西北师范大学 | 基于数字图像技术构建金属表面层三维宏观形貌的方法 |
-
2012
- 2012-02-27 CN CN 201210045676 patent/CN102607455B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN102607455A (zh) | 2012-07-25 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN102607455B (zh) | 基于光学显微镜和变化光照的微观形貌三维测量方法 | |
CN102800096B (zh) | 一种摄像机参数的鲁棒性估计算法 | |
CN108541322A (zh) | 三维视觉测量数据的处理方法和装置 | |
CN106504276A (zh) | 非局部立体匹配算法的组合匹配代价算法和视差联合填充算法 | |
Zhou et al. | A novel way of understanding for calibrating stereo vision sensor constructed by a single camera and mirrors | |
CN103325120A (zh) | 一种快速自适应支持权值双目视觉立体匹配方法 | |
CN103454276B (zh) | 一种基于动态序列图像的织物形态风格评价方法 | |
CN103353388A (zh) | 一种具摄像功能的双目体式显微成像***标定方法及装置 | |
CN103530880A (zh) | 基于投影高斯网格图案的摄像机标定方法 | |
CN101833666A (zh) | 一种对离散点云数据微分几何量的估计方法 | |
Nayak et al. | A modified triangle box-counting with precision in error fit | |
CN103308000B (zh) | 基于双目视觉的曲线物体测量方法 | |
CN107240132B (zh) | 一种利用灰度共生矩阵计算分形维数的方法 | |
CN102346857B (zh) | 人脸图像光照参数和去光照图的高精度同时估计方法 | |
CN116310219A (zh) | 一种基于条件扩散模型的三维脚型生成方法 | |
Willenborg et al. | Integration of semantic 3D city models and 3D mesh models for accuracy improvements of solar potential analyses | |
Yang et al. | A dual-platform laser scanner for 3D reconstruction of dental pieces | |
CN114723884A (zh) | 三维人脸重建方法、装置、计算机设备及存储介质 | |
CN104813217A (zh) | 用于设计能够估计景深的无源单通道成像器的方法 | |
CN106815871A (zh) | 一种扫描电子显微镜成像***的建模方法 | |
CN103954238B (zh) | 一种基于高斯函数的光纤干涉条纹图像背景光补偿方法 | |
CN115601423A (zh) | 一种双目视觉场景下基于边缘增强的圆孔位姿测量方法 | |
CN108364274A (zh) | 微纳米尺度下光学图像的无损清晰重构方法 | |
Wang et al. | Shape-from-focus reconstruction using block processing followed by local heat-diffusion-based refinement | |
Takimoto et al. | Shape reconstruction from multiple RGB-D point cloud registration |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant |