CN108007385A - 一种大视场弹坑表面形貌成像系统及方法 - Google Patents
一种大视场弹坑表面形貌成像系统及方法 Download PDFInfo
- Publication number
- CN108007385A CN108007385A CN201711496086.1A CN201711496086A CN108007385A CN 108007385 A CN108007385 A CN 108007385A CN 201711496086 A CN201711496086 A CN 201711496086A CN 108007385 A CN108007385 A CN 108007385A
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- mtd
- image
- matrix
- 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.)
- Granted
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 46
- 238000000034 method Methods 0.000 title claims abstract description 34
- 238000012876 topography Methods 0.000 title claims abstract description 33
- 230000000007 visual effect Effects 0.000 title claims abstract description 26
- 239000011159 matrix material Substances 0.000 claims abstract description 54
- 230000005622 photoelectricity Effects 0.000 claims abstract description 33
- 238000005070 sampling Methods 0.000 claims abstract description 21
- 230000001360 synchronised effect Effects 0.000 claims abstract description 11
- 238000012545 processing Methods 0.000 claims abstract description 9
- 238000005259 measurement Methods 0.000 claims description 27
- 230000008569 process Effects 0.000 claims description 12
- 239000000835 fiber Substances 0.000 claims description 11
- 238000006073 displacement reaction Methods 0.000 claims description 9
- 239000013307 optical fiber Substances 0.000 claims description 8
- 230000000694 effects Effects 0.000 claims description 7
- 238000001514 detection method Methods 0.000 claims description 6
- 238000000354 decomposition reaction Methods 0.000 claims description 5
- 238000004364 calculation method Methods 0.000 claims description 4
- 230000008859 change Effects 0.000 claims description 4
- 230000015572 biosynthetic process Effects 0.000 claims description 3
- 230000000638 stimulation Effects 0.000 claims description 3
- 210000004556 brain Anatomy 0.000 claims description 2
- 230000010354 integration Effects 0.000 claims description 2
- 239000000463 material Substances 0.000 claims description 2
- 230000008447 perception Effects 0.000 claims description 2
- 238000002310 reflectometry Methods 0.000 claims description 2
- 230000007704 transition Effects 0.000 claims description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 230000006870 function Effects 0.000 description 2
- 239000000203 mixture Substances 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 241000630665 Hada Species 0.000 description 1
- VYPSYNLAJGMNEJ-UHFFFAOYSA-N Silicium dioxide Chemical compound O=[Si]=O VYPSYNLAJGMNEJ-UHFFFAOYSA-N 0.000 description 1
- XUIMIQQOPSSXEZ-UHFFFAOYSA-N Silicon Chemical compound [Si] XUIMIQQOPSSXEZ-UHFFFAOYSA-N 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 210000001367 artery Anatomy 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 238000000205 computational method Methods 0.000 description 1
- 235000009508 confectionery Nutrition 0.000 description 1
- 230000005611 electricity Effects 0.000 description 1
- 230000004313 glare Effects 0.000 description 1
- 238000002372 labelling Methods 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 230000035800 maturation Effects 0.000 description 1
- 230000015654 memory Effects 0.000 description 1
- 230000000149 penetrating effect Effects 0.000 description 1
- 230000010363 phase shift Effects 0.000 description 1
- 230000002035 prolonged effect Effects 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000009738 saturating Methods 0.000 description 1
- 229910052710 silicon Inorganic materials 0.000 description 1
- 239000010703 silicon Substances 0.000 description 1
- 238000003786 synthesis reaction Methods 0.000 description 1
- 210000003462 vein Anatomy 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B11/00—Measuring arrangements characterised by the use of optical techniques
- G01B11/24—Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures
- G01B11/2433—Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures for measuring outlines by shadow casting
Abstract
本发明公开了一种大视场弹坑表面形貌成像系统及方法,系统包括同步控制系统、结构光投影系统、光电探测系统、光电门控系统以及信号处理系统,利用哈达玛矩阵调制的数字微镜阵列投影模板采样,从初始分辨率图像开始,逐渐获取目标场景由低分辨率到高分辨率的细节信息,用于重构最终分辨率三维图像。本发明对目标场景有自适应性,有效减少了采样次数,缩短了成像时间,适用于快速高分辨率三维成像应用;降低了计算开销,减少了重构所需的时间。
Description
技术领域
本发明属于激光三维成像技术领域,特别是一种大视场弹坑表面形貌成像系统及成像方法。
背景技术
目前,传统的三维成像系统主要采用条纹投影法和点扫描法([1]基于空间正交条纹投影的多分辨三维数字成像方法,CN200410019972.1.[2]一种基于TOF相机的三维扫描系统及扫描方法,CN201510868410.2.)。其中,条纹投影法三维测量过程可分为三个步骤:首先是条纹投影,计算机生成标准正弦条纹图案并由投影仪投射到物体表面,条纹结构光经物体表面调制后由相机采集得到变形的条纹图像;然后是相位计算,从这些变形条纹图像中计算出调制相位;最后是三维重建,通过系统标定从调制相位中重建出物体表面的三维形貌。例如:杨福俊等人([3]3D shape measurement of discontinuous objects basedon single frequency fringe projection in 4-step phase shifting profilemetry)提出了基于单频四步相移条纹投影的不连续物体三维形貌测量方法。首先计算机产生相移步长为π/2的4幅正弦条纹图像,利用液晶投影仪投影到被测物体表面,由CCD相机同步采集获得含有形貌信息的畸变条纹图并保存在计算机中,通过四步相移算法获得条纹相位信息。对于表面形貌不连续的物体,通过4幅畸变条纹图灰度算术和运算并进行二值化处理定位阴影或暗背景,利用二值图像修正解调的相位从而得到正确相位。
点扫描法测量过程是通过激光脉冲逐像素扫描,获得每个像素的TOF从而计算出像素点距离,得到单像素的距离信息的同时,通过扫描装置实现光束方向的控制,确定探测点的水平偏转角度和俯仰偏转角度,从而实现单点的空间定位,对扫描区域重复这个过程,最终获得整个采样区域的三维形貌。
上述条纹投影法由于要求条纹精密,计算方法复杂导致重建效率低下,不适合需要快速成像的场合,并且易受干扰的缺点使其在实际应用的效果受限。点扫描法尽管技术成熟,计算简单,但由于采样时间长,成像速度缓慢,在采集高分辨率图像时成像缓慢的缺点尤为明显,因此不适合用于快速成像。
目前的表面形貌快速成像方法有:Dai Huidong等人提出的AC3DI([4]Adaptivecompressed 3D imaging based on wavelet trees and Hadamard multiplexing with asingle photon counting detector,arXiv:1709.05961v1[cs.CV]15Sep 2017),SunMingjie等人提出的单像素三维成像方法(5]Single-pixel three-dimensional imagingwith time-based depth resolution,Nature Communications,2016,7:12010)等。由于AC3DI使用光子计数的方法成像,为了保证测量时能接收到足够光子,单次测量时间需要足够长,因而成像速度受限。另外,文献[5]提出的方法由于采样次数较多,导致采样总时间长,也难以满足快速成像的要求。
发明内容
本发明的目的在于提供一种简易、可靠、快速、准确、分辨率高、可自适应的大视场弹坑表面形貌成像系统及方法,解决传统表面形貌成像方法受限于计算复杂,成像时间长的问题。
实现本发明目的的技术解决方案为:一种大视场弹坑表面形貌成像系统,包括同步控制系统、结构光投影系统、光电探测系统、光电门控系统以及信号处理系统,所述同步控制系统由同步信号发生器组成,输出同步信号给长脉冲激光光源和短脉冲激光光源;
结构光投影系统包括长脉冲激光光源、分束镜、扩束镜、数字微镜器件DMD和投影透镜,长脉冲激光光源收到同步信号后发射脉冲激光,脉冲激光经过扩束镜后照亮DMD,经过DMD调制的脉冲激光通过投影透镜照射到目标场景区域内;
光电探测系统包括收束镜、窄带滤光片以及由光电雪崩二极管APD和放大器级联构成的单像素探测模块,激光回波在收束镜作用下收束,经过窄带滤光片后被单像素探测模块中的APD探测,经过放大器后将信号传输给光电门控系统;
光电门控系统包括短脉冲激光器、时间延迟模块和光电门控模块,其中时间延迟模块和光电门控模块之间通过多模光纤连接;
信号处理系统由采集卡和计算机PC组成,采集卡对光电门控模块的输出波形进行采样,PC对这些采样值进行处理,生成DMD图案和完成表面形貌成像。
一种大视场弹坑表面形貌成像方法,利用哈达玛矩阵对投影图案进行调制,从初始分辨率深度图像开始,利用多分辨率小波逼近逐渐获取目标场景中的深度细节信息,逐渐提高图像分辨率,最终通过细节信息和低分辨率深度图像得到高分辨率图像。
本发明与现有技术相比,其显著优点:(1)避免了条纹投影法成像所需的大量计算开销,同时避免了点扫描式表面形貌扫描系统中需要精确控制光束方向的要求,使设备简易可靠,稳定性好的特点。(2)采用的多分辨率小波逼近方法自适应采集目标场景重要深度信息,有效减少了测量次数,缩短了成像时间,实现了快速测量。(3)根据组合测量原理,用哈达玛矩阵对DMD投影图案进行调制,增加了图像的信噪比,提高了测量结果的精确度。(4)由于单次采样时间短,且无需调整阈值,减少了测量时间。(5)总采样次数少,计算量小,总耗时短,因此更适合快速成像。(6)采样时间可根据目标区域深度进行调整,对目标场景有自适应性,提高了采集效率。
下面结合附图对本发明作进一步详细描述。
附图说明
图1是本发明大视场弹坑表面形貌成像系统结构总图。
图2是本发明大视场弹坑表面形貌成像系统成像时间延迟模块结构图。
图3是本发明大视场弹坑表面形貌成像系统成像光电门控模块结构图。
图4是本发明大视场弹坑表面形貌成像系统成像方法计算调制图像/强度图像流程图。
图5是本发明大视场弹坑表面形貌成像方法流程图。
具体实施方式
如图1所示,本发明大视场弹坑表面形貌成像系统采用光电雪崩二极管和数字微镜器件组合的单像素相机结构对弹坑成像,获得弹坑的表面形貌,包括同步控制系统、结构光投影系统、光电探测系统、光电门控系统以及信号处理系统。同步控制系统由同步信号发生器h组成,输出同步信号给长脉冲激光光源a和短脉冲激光光源i。结构光投影系统包括长脉冲激光光源a、扩束镜b、数字微镜器件c(digital micromirror device,DMD)以及投影透镜d。长脉冲激光光源收到同步信号后发射脉冲激光,脉冲激光经过扩束镜后照亮DMD。经过DMD调制的脉冲激光通过投影透镜照射到目标场景区域内。光电探测系统包括收束镜e、窄带滤光片f以及由光电雪崩二极管g(avalanche photodiode,APD)和放大器级联构成的单像素探测模块。激光回波在收束镜作用下收束,经过窄带滤光片后被单像素探测模块中的APD探测,经过放大器后将信号传输给光电门控系统。光电门控系统包括短脉冲激光器i、时间延迟模块j和光电门控模块k。其中时间延迟模块和光电门控模块之间通过多模光纤连接,内部结构分别如图2和图3所示。信号处理系统由采集卡l和计算机PCm组成,采集卡l对光电门控模块的输出波形进行采样,PCm对这些采样值进行处理,生成DMD图案和完成表面形貌成像。
如图2所示,所述时间延迟模块包括窄带滤光片a2、位移反射镜组b2、光纤准直镜c2、激光-光纤耦合器d2,短脉冲激光器i接收到同步信号后发射出脉宽为数百ps的短激光脉冲,短激光脉冲进入时间延迟模块的窄带滤光片a2后,被位移反射镜组b2反射,经过光纤准直镜c2,最后从激光-光纤耦合器d2进入光纤。其中位移反射镜组b2由两面位移反射镜组成,可在电机控制下沿图2中的黑色箭头方向垂直位移。
如图3所示,所述光电门控模块包括fiberport耦合器a3、扩束镜b3、光电二极管c3、电阻d3,时间延迟模块的输出光通过光纤进入光电门控模块后,通过fiberport耦合器a3将光纤中的激光耦合成激光束,激光束被扩束镜b3扩大后照在光电二极管c3上,光电二极管c3在偏置电压和激光脉冲作用下形成光电门,探测器信号经过光电门后进入电阻d3,最后将电阻d3两端电压峰值作为输出值传给信号处理系统。
本发明所采用的大视场弹坑表面形貌成像系统工作过程如下:
对于一次成像过程,首先需要利用哈达玛矩阵在低分辨率下采集场景。然后根据小波系数关系预测并标记场景中的边缘区域。再根据标记在更高的分辨率下对边缘区域上采样,并通过哈达玛反变换计算出各像素的深度值。对于平滑区域,无需再次采集。同时,根据深度值中的最值重新分配位移反射镜组的位移区间。重复预测再采集这个过程直到达到目标分辨率。最后根据采集的数据恢复出场景的深度图像和三维形貌。
对于一次采样过程,先由同步信号发生器产生同步信号,长脉冲激光光源和短脉冲激光光源同时发射脉冲激光。来自长脉冲激光光源a的光进入扩束镜b扩束,照亮数字微镜器件d平面。根据PCm生成的调制图案,数字微镜器件d对来自长脉冲激光光源a的脉冲激光进行空间调制,经投影透镜d照射目标场景。目标场景反射调制脉冲激光经过收束镜e和窄带滤光片f滤除杂散光后,通过光电雪崩二极管g测量得到连续的回波信号,并由放大器对其进行放大得到探测器信号。
来自短脉冲激光光源a的光进入时间延迟模块j后,经过窄带滤光片a2后照射到位移反射镜组b2上。位移反射镜组可以沿图2中黑色尖头方向垂直微弱位移从而改变时间延迟模块j中激光的光程。激光离开位移反射镜组b2后经过光纤准直镜c2收束,通过激光-光纤耦合器d2进入光纤并传输到光电门控模块k。从时间延迟模块j传来的激光进入光电门控模块k后,在fiberport耦合器a3作用下耦合为激光束,并被扩束镜b3扩束,最终照射在光电二极管c3上。
当短脉冲激光光源a发射的激光脉冲光程与目标场景反射的调制脉冲激光光程相等时,光电二极管c3导通。由于短脉冲激光光源a发射的激光脉冲脉宽远小于目标场景反射的调制脉冲,并且其光程可根据位移反射镜组b2的位置确定,因此可以得到光电二极管c3导通时对应目标场景反射的调制脉冲激光的光程。采集卡l对光电二极管c3导通时加在电阻d3两端的电压值进行测量得到该光程下组合测量的结果。在PCm控制下多次改变位移反射镜组的位置进行采集即可得到不同光程下的组合测量结果。
PCm通过位移反射镜组的位置得到各采样值对应的飞行时间。将不同TOF的采样值综合得到TOF直方图。直方图中所有采样值纵坐标相加得到强度图像的组合采样值。各采样值纵坐标与对应横坐标相乘再求和得出本次投影的调制图像的组合采样值。本发明所采用的大视场弹坑表面形貌成像方法根据这些测量结果,恢复出目标场景三维信息以及生成下一阶段的DMD投影图案。PCm使用软件Labview控制系统,生成DMD投影图案,最终结果由Matlab生成。
本发明所采用大视场弹坑表面形貌成像方法,利用哈达玛矩阵对投影图案进行调制,从初始分辨率深度图像开始,利用多分辨率小波逼近逐渐获取目标场景中的深度细节信息,逐渐提高图像分辨率,最终通过细节信息和低分辨率深度图像得到高分辨率图像,结合图5,具体步骤如下:
第一步,数字微镜器件d使用哈达玛矩阵调制图案,对目标场景进行组合测量,得到初始分辨率深度图像。为了使系统简易,本发明使用S型哈达玛矩阵代替常规的H型哈达玛矩阵对图案进行调制。H型哈达玛矩阵生成方式为:
其中:
S型哈达玛矩阵由0和1元素组成,可由H型哈达玛矩阵变换构造。移除H型哈达玛矩阵第一行和第一列,将新方阵中的-1元素用0元素代替,即可得到S型哈达玛矩阵。
设初始分辨率为根据S型哈达玛矩阵特性,需要在矩阵最后一列之后再补一列0向量使矩阵的可以被DMD投影,每次测量的像素数为(m/2)-1,对于每一个图案,第m个像素处的元素恒定为0。根据组合测量原理和S型哈达玛矩阵特性,测量所需DMD调制图案个数为m-1。并且,每个DMD调制图案积分时间内获得的飞行时间离散信号记录了图案中各像素到达场景并反射回探测器的飞行时间以及对应的强度。测量完毕后,可由哈达玛反变换求出调制图像和强度图像,通过调制图像和强度图像对应像素相除可获得深度图像,图像的第m个像素的值由邻近像素的深度插值近似,最终得到图2中第一阶段的初始分辨率深度图像。
第二步,由于深度图像在小波域中更具有稀疏性,因此选用上一分辨率的深度图像作为当前分辨率细节像素位置估计的依据,可以减少所需采样次数。具体地,首先对上一分辨率采集的深度图像或边缘区域进行单层小波分解,得到其小波系数矩阵(生成表示当前分辨率采样位置的标记矩阵)。其中对应图像边缘、纹理等细节信息的小波系数称为重要小波系数,这些系数的绝对值相比平滑区域较大,只占自然图像小波系数矩阵中很少一部分。使用少量重要小波系数可以在不改变视觉效果的前提下重建图像。根据上一分辨率小波系数与阈值的关系,可估计出当前分辨率下重要系数的像素位置,即图2中的边缘标记区域。
其中估计过程和阈值的选取如下:对上一分辨率深度图像单层小波分解后,所得的小波系数矩阵分为低频子带(LL)、垂直高频子带(LH)、水平高频子带(HL)、对角高频子带(HH)。根据韦伯定理,若
则大脑可感受到该刺激的变化。式中k为韦伯分数,I为原始刺激强度,ΔI为刺激变化值。因此,将韦伯常数keye作为阈值,当的绝对值大于韦伯常数keye时,则该系数对应的位置在人眼看来是边缘位置。当小于韦伯常数keye时,若满足其绝对值大于keye,则对应位置的像素是边缘位置。
第三步,PC根据第二步估计结果,结合S型哈达玛矩阵生成DMD调制图案,其具体方法如下:
设由上一分辨率的深度图D1小波分解结果所得到的重要细节像素位置矩阵为M(n)。矩阵中值为1的位置表示该位置为图像细节,需要被再次采集;值为0的位置表示该位置处于平滑区域,无需再次采集。设当前分辨率m×m下需要采集的像素所在位置矩阵为M’(n),则M(n)与M’(n)存在如下关系:
M'=U(M)
其中,是克罗内克积运算。
通过上述步骤得到位置矩阵M’(n)后,即可在当前分辨率对k个含有细节的像素位置采样,k应当满足首先生成L阶S型哈达玛矩阵HL,并把矩阵HL中的-1转换为0,其中然后根据M’(n)对矩阵HL逐行按顺序重新排列成一个大小为k×m的用于当前阶段测量的感知矩阵Hj,则H可表示为:
其中,x、y与n存在如下关系:
x=[n/m]+1
y=n mod m
矩阵Hj的每一行取出并重排为一个m×m大小的矩阵后,即为DMD下一阶段采样的投影图案。
第四步,单像素探测模块对被目标反射的DMD投影图案测量,并将测量结果输出到光电探测模块,采集卡对光电探测模块输出采样。一次采样过程中,位移反射镜组多次移动,最终可得到一系列采样值。PC根据采样值可计算出TOF,从而由TOF重建深度图像,其过程如图4所示,具体方法如下:
一次采样过程中,同步控制系统h会根据位移反射镜组b2的总移动次数产生相同数量的同步脉冲。每产生一个同步脉冲,采集卡l就会得到一个采样值,同时位移反射镜组b2产生一个步长的位移。重复这个过程直到产生与位移反射镜组b2的总移动次数相同数量的同步脉冲,最终根据各采样值和对应时间得到TOF直方图。将TOF直方图中各点纵坐标相加,最终得到一个M×1的表示当前阶段测量结果的强度图像向量YI:
式中XI表示大小为N×1的向量化的强度图像。hmn是测量矩阵Hj第m行第n列的元素,表示DMD对应模板对应位置的状态(‘0’或‘1’)。ηmn是使用第m个模板测量时,到达DMD上第n个像素入射场景,经场景反射后的到达探测器的激光强度。
因为深度图像和测量结果之间是非线性关系,所以需要引入调制图像的概念,将深度图像和强度图像关联。调制图像XQ由强度图像XI和深度图像XD对应元素相乘得到,表示边缘的调制图像向量YQ则可以通过测量采集时的飞行时间来对深度进行线性采集,其表达式为:
XQ是一个N×1的调制图像向量,Tmn是飞行时间,C是一个常数,用于表征将飞行时间转换到深度的比率。则边缘信息的强度图像可以通过对YI做哈达玛反变换得以还原出来。随后,通过结合在当前阶段采样并恢复出的边缘强度图像和在上一阶段得到的强度图像X′I即可获得第j阶段采样结束后,对应分辨率下的强度图像XI,即:
同理,恢复深度图像的方法为:
XD=Nz(a)·iH(YQ)·/iH(YI)+U(X′D)·(1-M)
上式Nz(a)为一个零一矩阵,对于iH(YI)中的值为0的位置,Nz(a)对应位置为0,否则为1(为了避免分母为0的情况出现)。X′D为上一阶段采集的的深度图像,XD即为恢复出的当前阶段深度图像,如图2中第二、三阶段的深度图像。同时,根据当前阶段深度图像中的最大值和最小值,计算出下一阶段位移反射镜组所需的位移次数i:
式中XDmax为深度图像最大值,XDmin为深度图像最小值,step为位移反射镜组步长。
第五步,重复第二步到第四步,获取更高分辨率深度图像,直到获得最终分辨率深度图像。最后,PC使用软件matlab根据深度图像生成弹坑表面形貌三维模型,即图5中的三维形貌。
本发明所述成像系统的长脉冲激光光源和短脉冲激光光源频率均为40MHz,长脉冲激光光源脉宽为10ns,短脉冲激光光源脉宽为300ps。所用DMD是Vialux ALP 4.2系列,具有1024×768分辨率和4G片上内存,微镜最高翻转频率约为22.7kHz。投影透镜和收束镜分别选用50mm和35mm尼康标准镜头。单像素探测器选用索雷博的DET-100A硅偏置探测器。窄带滤光片透过波长宽度为10nm。
Claims (9)
1.一种大视场弹坑表面形貌成像系统,其特征在于包括同步控制系统、结构光投影系统、光电探测系统、光电门控系统以及信号处理系统,所述同步控制系统由同步信号发生器(h)组成,输出同步信号给长脉冲激光光源(a)和短脉冲激光光源(i);
结构光投影系统包括长脉冲激光光源(a)、分束镜(b)、扩束镜(c)、数字微镜器件DMD(d)和投影透镜(e),长脉冲激光光源(a)收到同步信号后发射脉冲激光,脉冲激光经过扩束镜(c)后照亮DMD(d),经过DMD(d)调制的脉冲激光通过投影透镜(e)照射到目标场景区域内;
光电探测系统包括收束镜(f)、窄带滤光片(g)以及由光电雪崩二极管APD(h)和放大器级联构成的单像素探测模块,激光回波在收束镜(f)作用下收束,经过窄带滤光片(g)后被单像素探测模块中的APD(h)探测,经过放大器后将信号传输给光电门控系统;
光电门控系统包括短脉冲激光器(i)、时间延迟模块(j)和光电门控模块(k),其中时间延迟模块(j)和光电门控模块(k)之间通过多模光纤连接;
信号处理系统由采集卡(l)和计算机PC(m)组成,采集卡(l)对光电门控模块的输出波形进行采样,PC(m)对这些采样值进行处理,生成DMD图案和完成表面形貌成像。
2.根据权利要求1所述的大视场弹坑表面形貌成像系统,其特征在于所述的时间延迟模块包括窄带滤光片(a2)、位移反射镜组(b2)、光纤准直镜(c2)、激光-光纤耦合器(d2),短脉冲激光器(i)接收到同步信号后发射出脉宽为数百ps的短激光脉冲,短激光脉冲进入时间延迟模块的窄带滤光片(a2)后,被位移反射镜组(b2)反射,经过光纤准直镜(c2),最后从激光-光纤耦合器(d2)进入光纤。
3.根据权利要求1或2所述的大视场弹坑表面形貌成像系统,其特征在于所述的光电门控模块包括fiberport耦合器(a3)、扩束镜(b3)、光电二极管(c3)、电阻(d3),时间延迟模块的输出光通过光纤进入光电门控模块后,通过fiberport耦合器(a3)将光纤中的激光耦合成激光束,激光束被扩束镜(b3)扩大后照在光电二极管(c3)上,光电二极管(c3)在偏置电压和激光脉冲作用下形成光电门,探测器信号经过光电门后进入电阻(d3),最后将电阻(d3)两端电压峰值作为输出值传给信号处理系统。
4.一种大视场弹坑表面形貌成像方法,其特征在于:利用哈达玛矩阵对投影图案进行调制,从初始分辨率深度图像开始,利用多分辨率小波逼近逐渐获取目标场景中的深度细节信息,逐渐提高图像分辨率,最终通过细节信息和低分辨率深度图像得到高分辨率图像。
5.根据权利要求4所述的大视场弹坑表面形貌成像方法,其特征在于具体步骤如下:
第一步,DMD使用哈达玛矩阵调制图案,对目标场景进行组合测量,得到初始分辨率深度图像;
第二步,根据上一分辨率深度图像,估计当前分辨率需要测量的细节像素所在区域,生成表示当前分辨率采样位置的标记矩阵;
第三步,PC根据第二步估计结果,结合S型哈达玛矩阵生成DMD调制图案;
第四步,单像素探测模块对被目标反射的DMD投影图案测量,并将测量结果输出到光电探测模块,采集卡对光电探测模块输出采样。一次采样过程中,位移反射镜组多次移动,最终可得到一系列采样值。PC根据采样值可计算出TOF,从而由TOF重建深度图像;
第五步,重复第二步到第四步,获取更高分辨率深度图像,直到获得最终分辨率深度图像,最后由PC使用软件matlab根据深度图像生成弹坑表面形貌三维模型。
6.根据权利要求5所述的大视场弹坑表面形貌成像方法,其特征在于第一步中,使用S型哈达玛矩阵代替常规的H型哈达玛矩阵对图案进行调制,H型哈达玛矩阵生成方式为:
<mfenced open = "" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>H</mi>
<mrow>
<mn>2</mn>
<mi>n</mi>
</mrow>
</msub>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<msqrt>
<mn>2</mn>
</msqrt>
</mfrac>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>H</mi>
<mi>n</mi>
</msub>
</mtd>
<mtd>
<msub>
<mi>H</mi>
<mi>n</mi>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>H</mi>
<mi>n</mi>
</msub>
</mtd>
<mtd>
<mrow>
<mo>-</mo>
<msub>
<mi>H</mi>
<mi>n</mi>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
</mtd>
<mtd>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>=</mo>
<msup>
<mn>2</mn>
<mi>m</mi>
</msup>
<mo>,</mo>
<mi>m</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mn>3</mn>
<mo>,</mo>
<mn>...</mn>
<mo>)</mo>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
其中:
<mrow>
<msub>
<mi>H</mi>
<mn>2</mn>
</msub>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
S型哈达玛矩阵由0和1元素组成,由H型哈达玛矩阵变换构造,移除H型哈达玛矩阵第一行和第一列,将新方阵中的-1元素用0元素代替,即可得到S型哈达玛矩阵;
设初始分辨率为根据S型哈达玛矩阵特性,在矩阵最后一列之后再补一列0向量使矩阵被DMD投影,每次测量的像素数为(m/2)-1,对于每一个图案,第m个像素处的元素恒定为0,根据组合测量原理和S型哈达玛矩阵特性,测量所需DMD调制图案个数为m-1;每个DMD调制图案积分时间内获得的飞行时间离散信号记录了图案中各像素到达场景并反射回探测器的飞行时间以及对应的强度;测量完毕后,由哈达玛反变换求出调制图像和强度图像,通过调制图像和强度图像对应像素相除可获得深度图像,图像的第m个像素的值由邻近像素的深度插值近似,最终得到第一阶段的初始分辨率深度图像。
7.根据权利要求5所述的大视场弹坑表面形貌成像方法,其特征在于第二步的估计过程和阈值选取为:对上一分辨率深度图像单层小波分解后,所得的小波系数矩阵分为低频子带LL、垂直高频子带LH、水平高频子带HL、对角高频子带HH,根据韦伯定理,若
<mrow>
<mi>k</mi>
<mo>=</mo>
<mfrac>
<mrow>
<mi>&Delta;</mi>
<mi>I</mi>
</mrow>
<mi>I</mi>
</mfrac>
</mrow>
则大脑可感受到该刺激的变化,式中k为韦伯分数,I为原始刺激强度,ΔI为刺激变化值;将韦伯常数keye作为阈值,当的绝对值大于韦伯常数keye时,则该系数对应的位置在人眼看来是边缘位置;
当小于韦伯常数keye时,若 满足其绝对值大于keye,则对应位置的像素是边缘位置。
8.根据权利要求5所述的大视场弹坑表面形貌成像方法,其特征在于第三步的具体过程如下:
设由上一分辨率的深度图D1小波分解结果所得到的重要细节像素位置矩阵为M(n),矩阵中值为1的位置表示该位置为图像细节,需要被再次采集;值为0的位置表示该位置处于平滑区域,无需再次采集;设当前分辨率m×m下需要采集的像素所在位置矩阵为M’(n),则M(n)与M’(n)存在如下关系:
M'=U(M)
<mrow>
<mi>U</mi>
<mrow>
<mo>(</mo>
<mi>X</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>X</mi>
<mo>&CircleTimes;</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
其中,是克罗内克积运算;
通过上述步骤得到位置矩阵M’(n)后,即可在当前分辨率对k个含有细节的像素位置采样,k应当满足首先生成L阶S型哈达玛矩阵HL,并把矩阵HL中的-1转换为0,其中然后根据M’(n)对矩阵HL逐行按顺序重新排列成一个大小为k×m的用于当前阶段测量的感知矩阵Hj,则H可表示为:
<mrow>
<msub>
<mi>H</mi>
<mi>j</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>H</mi>
<mi>L</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
<mtd>
<mrow>
<mi>i</mi>
<mi>f</mi>
<mi> </mi>
<msup>
<mi>M</mi>
<mo>&prime;</mo>
</msup>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mn>1</mn>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mrow>
<mi>i</mi>
<mi>f</mi>
<mi> </mi>
<msup>
<mi>M</mi>
<mo>&prime;</mo>
</msup>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mn>0</mn>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
其中,x、y与n存在如下关系:
x=[n/m]+1
y=n mod m
矩阵Hj的每一行取出并重排为一个m×m大小的矩阵后,即为DMD下一阶段采样的投影图案。
9.根据权利要求5所述的大视场弹坑表面形貌成像方法,其特征在于第三步的具体过程如下:
一次采样过程中,同步控制系统会根据位移反射镜组的总移动次数产生相同数量的同步脉冲,每产生一个同步脉冲,采集卡就会得到一个采样值,同时位移反射镜组产生一个步长的位移,重复这个过程直到产生与位移反射镜组的总移动次数相同数量的同步脉冲,最终根据各采样值和对应时间得到TOF直方图;
将TOF直方图中各点纵坐标相加,最终得到一个M×1的表示当前阶段测量结果的强度图像向量YI:
<mrow>
<msub>
<mi>Y</mi>
<mi>I</mi>
</msub>
<mo>=</mo>
<msub>
<mi>H</mi>
<mi>j</mi>
</msub>
<msub>
<mi>X</mi>
<mi>I</mi>
</msub>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>m</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>M</mi>
</munderover>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>n</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>N</mi>
</munderover>
<msub>
<mi>h</mi>
<mrow>
<mi>m</mi>
<mi>n</mi>
</mrow>
</msub>
<msub>
<mi>&eta;</mi>
<mrow>
<mi>m</mi>
<mi>n</mi>
</mrow>
</msub>
</mrow>
式中XI表示大小为N×1的向量化的强度图像,hmn是测量矩阵Hj第m行第n列的元素,表示DMD对应模板对应位置的状态(‘0’或‘1’),ηmn是使用第m个模板测量时,到达DMD上第n个像素入射场景,经场景反射后的到达探测器的激光强度;
将深度图像和强度图像关联,调制图像XQ由强度图像XI和深度图像XD对应元素相乘得到,表示边缘的调制图像向量YQ则可以通过测量采集时的飞行时间来对深度进行线性采集,其表达式为:
<mrow>
<msub>
<mi>Y</mi>
<mi>Q</mi>
</msub>
<mo>=</mo>
<msub>
<mi>H</mi>
<mi>j</mi>
</msub>
<msub>
<mi>X</mi>
<mi>Q</mi>
</msub>
<mo>=</mo>
<mi>C</mi>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>m</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>M</mi>
</munderover>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>n</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>N</mi>
</munderover>
<msub>
<mi>h</mi>
<mrow>
<mi>m</mi>
<mi>n</mi>
</mrow>
</msub>
<msub>
<mi>&eta;</mi>
<mrow>
<mi>m</mi>
<mi>n</mi>
</mrow>
</msub>
<msub>
<mi>T</mi>
<mrow>
<mi>m</mi>
<mi>n</mi>
</mrow>
</msub>
</mrow>
XQ是一个N×1的调制图像向量,Tmn是飞行时间,C是一个常数,用于表征将飞行时间转换到深度的比率,则边缘信息的强度图像可以通过对YI做哈达玛反变换得以还原出来;随后,通过结合在当前阶段采样并恢复出的边缘强度图像和在上一阶段得到的强度图像X′I即可获得第j阶段采样结束后,对应分辨率下的强度图像XI,即:
<mrow>
<msub>
<mi>X</mi>
<mi>I</mi>
</msub>
<mo>=</mo>
<msubsup>
<mi>X</mi>
<mi>I</mi>
<mi>e</mi>
</msubsup>
<mo>+</mo>
<mi>U</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>X</mi>
<mi>I</mi>
<mo>&prime;</mo>
</msubsup>
<mo>)</mo>
</mrow>
<mo>&CenterDot;</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<msup>
<mi>M</mi>
<mo>&prime;</mo>
</msup>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<msubsup>
<mi>X</mi>
<mi>I</mi>
<mi>e</mi>
</msubsup>
<mo>=</mo>
<mi>i</mi>
<mi>H</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>Y</mi>
<mi>I</mi>
</msub>
<mo>)</mo>
</mrow>
</mrow>
同理,恢复深度图像的方法为:
XD=Nz(a)·iH(YQ)·/iH(YI)+U(X′D)·(1-M)
上式Nz(a)为一个零一矩阵,对于iH(YI)中的值为0的位置,Nz(a)对应位置为0,否则为1,X′D为上一阶段采集的的深度图像,XD即为恢复出的当前阶段深度图像;同时,根据当前阶段深度图像中的最大值和最小值,计算出下一阶段位移反射镜组所需的位移次数i:
<mrow>
<mi>i</mi>
<mo>=</mo>
<mfrac>
<mrow>
<msub>
<mi>X</mi>
<mrow>
<mi>D</mi>
<mi>m</mi>
<mi>a</mi>
<mi>x</mi>
</mrow>
</msub>
<mo>-</mo>
<msub>
<mi>X</mi>
<mrow>
<mi>D</mi>
<mi>min</mi>
</mrow>
</msub>
</mrow>
<mrow>
<mi>s</mi>
<mi>t</mi>
<mi>e</mi>
<mi>p</mi>
</mrow>
</mfrac>
</mrow>
式中XDmax为深度图像最大值,XDmin为深度图像最小值,step为位移反射镜组步长。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201711496086.1A CN108007385B (zh) | 2017-12-31 | 2017-12-31 | 一种大视场弹坑表面形貌成像系统及方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201711496086.1A CN108007385B (zh) | 2017-12-31 | 2017-12-31 | 一种大视场弹坑表面形貌成像系统及方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN108007385A true CN108007385A (zh) | 2018-05-08 |
CN108007385B CN108007385B (zh) | 2023-11-14 |
Family
ID=62049841
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201711496086.1A Active CN108007385B (zh) | 2017-12-31 | 2017-12-31 | 一种大视场弹坑表面形貌成像系统及方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108007385B (zh) |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108895985A (zh) * | 2018-06-19 | 2018-11-27 | 中国科学院合肥物质科学研究院 | 一种基于单像素探测器的物体定位方法 |
WO2022110947A1 (zh) * | 2020-11-25 | 2022-06-02 | Oppo广东移动通信有限公司 | 电子装置的控制方法、电子装置及计算机可读存储介质 |
CN115290004A (zh) * | 2022-10-08 | 2022-11-04 | 南京理工大学 | 一种基于压缩感知与hsi的水下并行单像素成像方法 |
Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102788770A (zh) * | 2012-08-16 | 2012-11-21 | 中国科学技术大学 | 一种基于双波长区分三种气溶胶颗粒的方法 |
CN102865833A (zh) * | 2012-10-17 | 2013-01-09 | 南京理工大学 | 基于等高信息稀疏测量的三维成像装置及方法 |
CN103234479A (zh) * | 2013-04-10 | 2013-08-07 | 浙江大学 | 基于压缩感知脉冲编码调制的三维信息获取装置和方法 |
CN103913228A (zh) * | 2014-04-09 | 2014-07-09 | 辽宁大学 | 编码模板多目标超分辨率时间飞行成像系统及方法 |
EP2755327A1 (en) * | 2011-10-25 | 2014-07-16 | Center For Space Science And Applied Research, Chinese Academy of Sciences | Time-resolved single-photon or ultra-weak light multi-dimensional imaging spectrum system and method |
US9163936B1 (en) * | 2012-05-07 | 2015-10-20 | Physical Optics Corporation | Three-dimensional profilometer |
CN106772430A (zh) * | 2016-12-30 | 2017-05-31 | 南京理工大学 | 基于多分辨率小波逼近的单像素光子计数三维成像系统及方法 |
JP2018516377A (ja) * | 2015-09-10 | 2018-06-21 | テンヨウン サンウェイ(天津)科技有限公司 | マルチラインアレイレーザ光3次元走査システム、及びマルチラインアレイレーザ光3次元走査方法 |
CN207779348U (zh) * | 2017-12-31 | 2018-08-28 | 中国人民解放军陆军工程大学 | 一种大视场弹坑表面形貌成像系统 |
-
2017
- 2017-12-31 CN CN201711496086.1A patent/CN108007385B/zh active Active
Patent Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP2755327A1 (en) * | 2011-10-25 | 2014-07-16 | Center For Space Science And Applied Research, Chinese Academy of Sciences | Time-resolved single-photon or ultra-weak light multi-dimensional imaging spectrum system and method |
US9163936B1 (en) * | 2012-05-07 | 2015-10-20 | Physical Optics Corporation | Three-dimensional profilometer |
CN102788770A (zh) * | 2012-08-16 | 2012-11-21 | 中国科学技术大学 | 一种基于双波长区分三种气溶胶颗粒的方法 |
CN102865833A (zh) * | 2012-10-17 | 2013-01-09 | 南京理工大学 | 基于等高信息稀疏测量的三维成像装置及方法 |
CN103234479A (zh) * | 2013-04-10 | 2013-08-07 | 浙江大学 | 基于压缩感知脉冲编码调制的三维信息获取装置和方法 |
CN103913228A (zh) * | 2014-04-09 | 2014-07-09 | 辽宁大学 | 编码模板多目标超分辨率时间飞行成像系统及方法 |
JP2018516377A (ja) * | 2015-09-10 | 2018-06-21 | テンヨウン サンウェイ(天津)科技有限公司 | マルチラインアレイレーザ光3次元走査システム、及びマルチラインアレイレーザ光3次元走査方法 |
CN106772430A (zh) * | 2016-12-30 | 2017-05-31 | 南京理工大学 | 基于多分辨率小波逼近的单像素光子计数三维成像系统及方法 |
CN207779348U (zh) * | 2017-12-31 | 2018-08-28 | 中国人民解放军陆军工程大学 | 一种大视场弹坑表面形貌成像系统 |
Non-Patent Citations (1)
Title |
---|
马彦鹏;王亚南;王义坤;葛明锋;王雨曦;亓洪兴;舒嵘;: "基于压缩感知的单点探测计算成像技术研究", 光学学报, no. 12 * |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108895985A (zh) * | 2018-06-19 | 2018-11-27 | 中国科学院合肥物质科学研究院 | 一种基于单像素探测器的物体定位方法 |
WO2022110947A1 (zh) * | 2020-11-25 | 2022-06-02 | Oppo广东移动通信有限公司 | 电子装置的控制方法、电子装置及计算机可读存储介质 |
CN115290004A (zh) * | 2022-10-08 | 2022-11-04 | 南京理工大学 | 一种基于压缩感知与hsi的水下并行单像素成像方法 |
CN115290004B (zh) * | 2022-10-08 | 2023-02-28 | 南京理工大学 | 一种基于压缩感知与hsi的水下并行单像素成像方法 |
Also Published As
Publication number | Publication date |
---|---|
CN108007385B (zh) | 2023-11-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN106772430B (zh) | 基于多分辨率小波逼近的单像素光子计数三维成像系统及方法 | |
CN104267407B (zh) | 基于压缩采样的主动成像方法和系统 | |
KR102456875B1 (ko) | 심도 촬상 장치, 방법 및 응용 | |
KR101214702B1 (ko) | 웨이블릿 변환에 기반한 3차원 측정 방법 | |
KR102624984B1 (ko) | ToF 센서와 이를 이용한 3차원 영상 장치 및 3차원 영상 장치의 구동 방법 | |
CN207779348U (zh) | 一种大视场弹坑表面形貌成像系统 | |
CN111351450B (zh) | 基于深度学习的单帧条纹图像三维测量方法 | |
CN1330928C (zh) | 一种采用双波长结构光测量物体轮廓的方法 | |
CN108007385A (zh) | 一种大视场弹坑表面形貌成像系统及方法 | |
CN111047681B (zh) | 基于深度学习的单像素三维端到端重建方法及装置 | |
CN110264540B (zh) | 一种并行单像素成像方法 | |
JP2009523040A (ja) | 1以上のスペクトルを符号化する内視鏡技術によるデータ生成システム及び方法 | |
CN106546513A (zh) | 一种基于正交双视场的三维降水粒子测量与重构装置及方法 | |
CN103234479B (zh) | 基于压缩感知脉冲编码调制的三维信息获取装置和方法 | |
CN104266587A (zh) | 一种三维测量系统及获得真实3d纹理点云数据方法 | |
CN107121682B (zh) | 一种基于相位式激光测距的三维关联成像方法 | |
CN113379818B (zh) | 一种基于多尺度注意力机制网络的相位解析方法 | |
CN113114882B (zh) | 一种高采样效率的傅里叶单像素成像方法 | |
CN104680534B (zh) | 基于单帧复合模板的物体深度信息获取方法 | |
KR101011925B1 (ko) | 삼차원 형상 측정장치 | |
CN106067162A (zh) | 集成成像超分辨率微单元图像阵列采集与重构方法 | |
CN102865833B (zh) | 基于等高信息稀疏测量的三维成像装置及方法 | |
Gutierrez-Barragan et al. | itof2dtof: A robust and flexible representation for data-driven time-of-flight imaging | |
Je et al. | Color-phase analysis for sinusoidal structured light in rapid range imaging | |
CN108447113A (zh) | 一种基于深度学习模型的脉冲式强度关联三维成像方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |