CN111340902A - 光相位调制方法及任意位置和形状照射的空间光调制方法 - Google Patents

光相位调制方法及任意位置和形状照射的空间光调制方法 Download PDF

Info

Publication number
CN111340902A
CN111340902A CN201911277585.0A CN201911277585A CN111340902A CN 111340902 A CN111340902 A CN 111340902A CN 201911277585 A CN201911277585 A CN 201911277585A CN 111340902 A CN111340902 A CN 111340902A
Authority
CN
China
Prior art keywords
amplitude
distribution
complex amplitude
algorithm
modulation method
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
Application number
CN201911277585.0A
Other languages
English (en)
Other versions
CN111340902B (zh
Inventor
李鹏程
闻栋
陆锦玲
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hust-Suzhou Institute For Brainsmatics
Original Assignee
Hust-Suzhou Institute For Brainsmatics
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Hust-Suzhou Institute For Brainsmatics filed Critical Hust-Suzhou Institute For Brainsmatics
Priority to CN201911277585.0A priority Critical patent/CN111340902B/zh
Publication of CN111340902A publication Critical patent/CN111340902A/zh
Application granted granted Critical
Publication of CN111340902B publication Critical patent/CN111340902B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • GPHYSICS
    • G03PHOTOGRAPHY; CINEMATOGRAPHY; ANALOGOUS TECHNIQUES USING WAVES OTHER THAN OPTICAL WAVES; ELECTROGRAPHY; HOLOGRAPHY
    • G03HHOLOGRAPHIC PROCESSES OR APPARATUS
    • G03H1/00Holographic processes or apparatus using light, infrared or ultraviolet waves for obtaining holograms or for obtaining an image from them; Details peculiar thereto
    • G03H1/04Processes or apparatus for producing holograms
    • G03H1/08Synthesising holograms, i.e. holograms synthesized from objects or objects from holograms
    • G03H1/0866Digital holographic imaging, i.e. synthesizing holobjects from holograms
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02DCLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
    • Y02D30/00Reducing energy consumption in communication networks
    • Y02D30/70Reducing energy consumption in communication networks in wireless communication networks

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Optical Modulation, Optical Deflection, Nonlinear Optics, Optical Demodulation, Optical Logic Elements (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开了一种光相位调制方法及任意位置和形状照射的空间光调制方法,其中,所述光相位调制方法包括,傅里叶变换步骤;第一幅值替换步骤;傅里叶逆变换步骤;第二幅值替换步骤和将上述步骤进行迭代计算。根据本发明的光相位调制方法可以通过计算全息的方法获得所需的光场强度分布,从而可以对光进行任意位置及形状的多点照射,根据算法中误差衰减的过程,增加背景值来提高迭代算法的性能和图像均一性,通过掩模的优化迭代算法,抑制支持域中背景部分的幅值分布,从而消除其对光场分布的影响,通过空间滤波的方法滤除支持域外的背景,有着更高的光能利用率。

Description

光相位调制方法及任意位置和形状照射的空间光调制方法
技术领域
本发明涉及成像技术领域,尤其是涉及一种光相位调制方法及任意位置和形状照射 的空间光调制方法。
背景技术
对于光学调控系统而言,选取的往往是视场内(比如大脑皮层)的某个点、或某些点、或者某些区域进行光照射。相对于整个视场范围而言,这些区域的面积相对较小。 目前一种是通过DMD(Digital Micromirror Device)一类的光强调制器件,调控输入 光束使得其光照范围限制在刺激区域处,这种方法光利用率较低。另一种是通过计算全 息方法提高光能利用率,通过计算全息的方法,直接利用计算机模拟光波衍射和干涉的 过程,得到数字化的相息图,主要包括振幅型计算全息图和相位型计算全息图,其中, 相位型计算全息图的方法,可以通过优化迭代算法逆向获得输入平面(调制平面)的相 位分布,常见的迭代算法有傅里叶迭代算法、菲涅尔迭代算法、模拟退火算法等。其中, 傅里叶迭代算法利用了快速傅里叶变化将光强分布转换到频域中的频谱分布,并在频域 面上与像面上进行幅值约束,进而逐步降低误差。
但是,目前的相位傅里叶迭代算法,随着迭代的增加,误差衰减逐渐减少,甚至迭代前后的误差几乎不变,算法停滞,进而造成光调制面的光利用率不理想,有待进一步 优化。
发明内容
本发明旨在至少解决现有技术中存在的技术问题之一。
本发明第一方面实施例提供了一种光相位调制方法,该光相位调制方法可以进一步 增大误差衰减,提高空间光调制光利用率。
本发明第二方面实施例提供了一种任意位置和形状照射的空间光调制方法。
为解决上述问题,本发明第一方面实施例的光相位调制方法包括:赋予输入光束一 初始相位分布φm,n,定义输入光束复振幅分布的幅值|Um,n|;傅里叶变换步骤,对输 入光束复振幅分布um,n进行傅里叶变换以获得输入光束在频域上的复振幅分布
Figure RE-GDA0002498205000000021
第一幅值替换步骤,将所述输入光束在频域上的复振幅分布的幅值
Figure RE-GDA0002498205000000022
替换为目标光场复振幅分布的幅值|Au,v|,获得幅值替换后的复振幅分布f′u,v=|Au,v|exp(iφu,v),其中,所述目标光场复振幅分布的幅值包括目标图像及背 景值,且所述背景值不为零;傅里叶逆变换步骤,对幅值替换后的复振幅分布f′u,v进行 傅里叶逆变换,获得所述输入光束在空域上的复振幅分布
Figure RE-GDA0002498205000000023
第二幅值替换步骤,将所述输入光束在空域上的复振幅分布的幅值
Figure RE-GDA0002498205000000024
替换为所述输 入光束的复振幅分布的幅值|Um,n|,获得空域上的新的输入光束复振幅分布 g′m,n=|Um,n|exp(iφ′m,n);以所述傅里叶变换步骤、所述第一幅值替换步骤、所述 傅里叶逆变换步骤、所述第二幅值替换步骤进行迭代计算,获得对输入光束进行相位调 制所需的相位分布函数φ′m,n
根据本发明实施例的光相位调制方法,通过在第一幅值替换步骤中,增加目标图像 及背景值,其中,背景值不为零,使得背景区域可以参与误差的衰减过程,从而可以增大迭代算法的误差衰减,进而可以使得迭代算法获得的相息图有着更高的均一性和更小的误差分布。
在一些实施例中,当所述输入光束在频域上的复振幅分布与目标光场的复振幅分布 小于规定的误差值,或者,迭代计算的次数达到预设次数时,确定所述傅里叶逆变换步骤所获得的输入光束复振幅分布中的相位分布为对输入光束进行相位调制所需的相位 分布函数。在一些实施例中,将所述目标图像添加在新的背景板上成为新的目标图像, 初始的所述目标图像的尺寸为nx×ny,所述新的背景板的尺寸为Nx×Ny,其中,Nx>nx, Ny>ny,在所述新的目标图像上制作图像掩模,所述图像掩模的尺寸与初始的所述目标 图像的尺寸相当,其中,在所述第一幅值替换步骤中,对应于所述图像掩模支持域内和 所述图像掩模支持域外采用不同的背景值。
在一些实施例中,所述第一幅值替换步骤包括:将所述输入光束在频域上的复振幅 分布的幅值替换为目标光场复振幅分布的幅值,获得幅值替换后的复振幅分布,所述目标光场复振幅分布的幅值为所述包含背景值的新的目标图像,其中,对应于所述图像掩 模支持域内,所述背景值为零,对应于所述图像掩模支持域外,所述背景值不为零。
在一些实施例中,在所述图像掩模支持域内,所述目标光场复振幅分布的幅值|Vu,v| 满足以下约束条件:
Figure RE-GDA0002498205000000031
其中,E为单位矩阵,
Figure RE-GDA0002498205000000032
为目标光场复振幅分布的幅值,
Figure RE-GDA0002498205000000033
为迭代计算中归一 化的输入光束频域复振幅分布的的幅值,Wk和Wk-1为迭代结果,k为迭代次数。
在一些实施例中,所述光相位调制方法还包括:在进行迭代计算时,保留所述图像掩模支持域外幅值的迭代计算结果,仅对所述图像掩模支持域内的目标图像部分进行幅值约束,其中,在迭代过程中,所述第一幅值替换步骤在频域内对所述图像掩模支持域 内的幅值A′n采用以下约束条件:A′n=A0·mask+An·(1-mask),其中,A0为目标图像 幅值,mask为图像掩模幅值,An为第n次迭代的幅值结果。
在一些实施例中,所述方法还包括:函数初始化步骤,赋予初始输入的输入光束初始相 位分布,并生成二维矩阵,以及根据所述输入光束确定所述输入光束的复振幅函数分布或者 输入光束采用平面波函数分布或高斯光束波面函数分布。
在一些实施例中,所述光相位调制方法根据噪声幅值尺度设置所述目标图像的背景 值。
进一步地,所述目标图像的背景值的取值范围为0.1倍-0.2倍。
为解决上述问题,本发明第二方面实施例的任意位置和形状照射的空间光调制方法, 包括:根据所述的光相位调制方法获取对输入光束进行光相位调制的调制相位分布函数; 根据所述调制相位分布函数,利用相位型空间光调制器对输入光束进行调制,在衍射平 面获得任意位置及形状的目标光强分布。
根据本发明实施例的任意位置和形状照射的空间光调制方法,通过采用所述的光相 位调制方法获得对输入光束进行光相位调制的调制相位分布函数,可以获得更加精准的 血流成像,利于疾病诊断。
本发明的附加方面和优点将在下面的描述中部分给出,部分将从下面的描述中变得 明显,或通过本发明的实践了解到。
附图说明
本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明 显和容易理解,一个或多个实施例通过与之对应的附图进行示例性说明,这些示例性说明和附图并不构成对实施例的限定,附图中具有相同参考数字标号的元件示为类似的元件,附图不构成比例限制,并且其中:
图1是傅里叶迭代算法的流程图;
图2中(a)和(b)为针对图1的傅里叶迭代算法的误差分析原理图;
图3中(a)和(b)为针对图1的傅里叶迭代算法的不同迭代次数的成像效果和误 差变化的示意图;
图4中(a)和(b)为相关技术中G-S算法和加权G-S算法的不同迭代次数与衍射 效率变化的曲线图;
图5为相关技术中G-S算法与加权G-S算法不同迭代次数下的均一性变化的曲线图;
图6中(a)、(b)和(c)为针对图1的傅里叶迭代算法的误差分析原理图;
图7为针对图1的傅里叶迭代算法的误差分析原理图;
图8为针对图1的傅里叶迭代算法不同迭代次数的误差衰减曲线图;
图9中(a)、(b)、(c)和(d)为G-S算法迭代过程中幅角差变化的示意图;
图10是根据本发明一个实施例的光相位调制方法的流程图;
图11中(a)和(b)为不同背景空间占比的目标图像的示意图,图11中(c)为针 对图11中(a)和(b)的相息图的误差的曲线图;
图12中(a)为增加背景值G-S算法和原始G-S算法的均一性的曲线图;图12中(b)为增加背景值G-S算法和原始G-S算法的误差变化的曲线图;
图13中(a)和(b)为添加不同背景值的傅里叶算法的均一性和误差变化的曲线图;
图14中(a)和(b)为根据本发明的一个实施例的添加背景和制作掩模的示意图;
图15中(a)、(b)和(c)为根据本发明的一个实施例的增加掩模后的迭代算法 迭代结果的示意图;
图16为添加掩模G-S算法和原始G-S算法的幅值的分布图;
图17中(a)、(b)和(c)为添加掩模G-S算法和原始G-S算法的衍射效率、误 差和均一性变化的曲线图;
图18中(a)、(b)和(c)为添加不同背景值并添加掩模的G-S算法的衍射效率、 误差和均一性变化的曲线图;
图19中(a)、(b)和(c)为添加掩模的加权G-S算法和未添加掩模的加权G-S 算法的衍射效率、误差和均一性变化的曲线图;
图20为不同傅里叶迭代算法下的图像对比度的比较图;
图21为增加掩模G-S算法和背景值保留G-S算法的幅值的分布图;
图22是根据本发明的一个实施例的血流成像空间光调制方法的流程图。
具体实施方式
为了能够更加详尽地了解本发明实施例的特点与技术内容,下面结合附图对本发明 实施例的实现进行详细阐述,所附附图仅供参考说明之用,并非用来限定本发明实施例。 在以下的技术描述中,为方便解释起见,通过多个细节以提供对所披露实施例的充分理 解。然而,在没有这些细节的情况下,一个或多个实施例仍然可以实施。在其它情况下,为简化附图,熟知的结构和装置可以简化展示。
在血流成像系统中,对于通过计算全息生成选择性空间光照射图案的方法,迭代算 法的目标是使经过调制的光束在衍射面上生成的衍射图案与目标图案一致,迭代算法能 够实现的条件是算法的进行应该沿着误差减小的方向进行。本发明实施例的方法对光相 位调整迭代算法进行改进。
下面先对目前采用的傅里叶变换算法(简称G-S算法)进行介绍,进而对目前的迭代算法进行评估。
图1是傅里叶变换算法的流程图,如图1所示具体步骤包括:
步骤①-初始化光束复振幅分布um,n,即赋予输入光束一初始相位分布φm,n,由计算机随机生成二维矩阵,m、n为目标图像的像素数目,定义输入光束复振幅分布的幅值 |Um,n|:输入光束复振幅分布um,n的幅值|Um,n|可由实际的光束分布决定,也可定义 为平面波函数或者高斯光束波面函数。
步骤②-对输入光束复振幅分布um,n进行离散傅里叶变换(Discrete FourierTransform,简称DFT),得到输入光束频域上的复振幅分布
Figure RE-GDA0002498205000000051
步骤③-在频域上对所得的复振幅分布fu,v进行约束,使用目标图案约束计算的复振 幅分布fu,v,即用目标光场复振幅分布的幅值(目标光强分布)|Au,v|替换计算所得复 振幅分布fu,v的幅值
Figure RE-GDA0002498205000000052
进而得到f′u,v=|Au,v|exp(iφu,v),使得得到的复振幅分 布f′u,v中的强度分布|Au,v|为目标衍射光强分布。
步骤④-对f′u,v进行离散傅里叶逆变换(DFT),得到输入光束在空域上的复振幅分布
Figure RE-GDA0002498205000000053
步骤⑤-误差判定或次数判定,当计算所得的输入光束在空域上的复振幅分布与目标 图案光场的复振幅分布函数误差小于设定值ε时,或者迭代次数大于最大迭代次数n时, 即跳出迭代并输出④中的相位分布φ′m,n,否则继续进入循环迭代运算。
步骤⑥-在空域上对复振幅分布进行约束得到新的输入光束复振幅分布,使用输入光 束的空间幅度进行约束,即使用输入光束的复振幅分布的幅值|Um,n|替换
Figure RE-GDA0002498205000000061
得到 空域上的新的输入光束复振幅分布g′m,n=|Um,n|exp(iφ′m,n)。
步骤⑦输出相位分布:进行多次迭代后,当频域上的复振幅分布的幅值
Figure RE-GDA0002498205000000062
与期 望所得的目标衍射光强分布|Au,v|达到相似的预期值,
Figure RE-GDA0002498205000000063
则为目标振幅分布,这时,空域上的复振幅分布的相位值φ′m,n即为对输入光束进行相位调制所需的相位分布函数,即目标相位分布函数。
定义误差函数
Figure RE-GDA0002498205000000064
为第i次迭代的误差,m、n为空间采样 数,显然,当εi=0时,
Figure RE-GDA0002498205000000065
即为目标振幅分布,对应的φi′就是相位分布。为了能获得所需的相位分布,需要迭代算法是收敛的,即:εi+1≤εi
如图2所示,(a)和(b)分别为G-S算法空域和频域的误差分析效果图,由于傅 里叶变换对的幺正性,根据帕赛瓦尔定理,即∫∫|f(x,y)|2dx dy=∫∫|F(u,v)|2du dv, 对于空域和频域对应的两参考点
Figure RE-GDA0002498205000000066
Figure RE-GDA0002498205000000067
Figure RE-GDA0002498205000000068
|t1|为输入光 束的复振幅分布,|t2|为目标衍射光强分布。
根据G-S算法对|t1|赋予了初始相位
Figure RE-GDA0002498205000000069
经傅里叶变换后在频域内得到
Figure RE-GDA00024982050000000610
其 中|h2|的分布与|t2|是不一致的,但有
Figure RE-GDA00024982050000000611
通过保留h2的幅 角、限制幅度来约束迭代算法,因此使得|h2|+|d2|=|h′2|=|t2|,其中
Figure RE-GDA00024982050000000612
将h′2进行逆傅里叶变换为空间域,得到
Figure RE-GDA00024982050000000613
对于
Figure RE-GDA00024982050000000614
有g1=F(h′2)=F(h2+d2)=F(h2)+F(d2)=t1+c1,其中,c1=F(d2),且
Figure RE-GDA00024982050000000615
Figure RE-GDA00024982050000000616
同样的,通过保留g1的幅角、限制幅度来约束迭代算法,使得|g1|+|e1|= |g′1|=|t1|,对于空域,有
Figure RE-GDA00024982050000000617
根据矢量的三角运算可知,|e1|≤|c1|,所以 有
Figure RE-GDA00024982050000000618
即ε′≤ε。再次进行傅里叶运算和幅度限制后,有ε″≤ε′≤ε,即随着迭代的进行,误差函数始终是减小的,满足εi+1≤εi,即迭代算法 是收敛的。
如图3所示,为G-S算法的收敛性分析,图3中(a)为随着迭代次数n(n=1,2,3,4)时所得的衍射图样,即目标图案,经过两次迭代后,所得相位已基本符合需求,因此G-S 算法本质是一种误差快速收敛算法;图3中(b)为随着迭代次数n(n=20)的算法误差 函数
Figure RE-GDA0002498205000000071
从结果中可已看出,误差函数ε随着迭代次数的增加迅速 下降,当迭代次数达到5次后,随着迭代次数的增加,误差函数ε基本保持不变。通过 在迭代的过程中加入上一次迭代的结果作为加权,改善算法的搜索过程,可以提高了迭 代算法的精确度。
改进的加权G-S算法主要是在原有的迭代算法中的步骤③中,使用一个新的加权幅 值|Vuv|代替|Auv|对复振幅分布进行约束,其中:
Figure RE-GDA0002498205000000072
其中,E为单位矩阵,
Figure RE-GDA0002498205000000073
Figure RE-GDA0002498205000000074
为归一化的幅值分布。改进的加权算法通过Wk这一迭代结果与目标的比值控制迭代的进程,并优化算法的结果。
迭代算法可以通过光束利用率来评估,即光束的衍射效率η,定义为:
Figure RE-GDA0002498205000000075
其中,P衍射图案中的光强分布目标区域,即所需的光强分布区域。
如图4中(b)所示,为标准G-S算法与改进的加权G-S算法在图4中(a)进行相 位迭代运算的衍射效率,在迭代次数达到10次后,算法的衍射效率趋于平缓,标准G-S 算法的衍射效率稍高于加权G-S算法(70%~67%),可见,改进的加权G-S算法能够优 化迭代效果。
另一方面可以通过调制光场的均一性进行评估,在光调控研究中,对于光斑的均匀 性有较高的要求,比如激光技工需要对光能量进行均匀控制,生物医学研究中刺激光斑的均匀能保证刺激强度的一致。均一性uniformity定义为:
Figure RE-GDA0002498205000000076
如图5所示,对于标准G-S算法与改进的加权G-S算法,对比其结果衍射图案的均一性,图4中(a)为目标图案,图6为结果。从图中可见,标准的G-S算法在迭代次 数达到20次左右时,其调制光场均一性已基本达到极限,约为0.65,而改进的加权G-S 算法随着迭代次数的增加,调制光场均一性迅速增加,迭代次数10次左右均一性即超 过了标准G-S算法,并能随着迭代次数的增加而继续增加,最终稳定在0.95以上。结 果表明,加权G-S算法相比标准的G-S相位迭代算法对于光场的调制有着更高的均一性。
由于误差的衰减依赖于误差矢量与原始矢量的空间关系,对于Nx×Ny的二维图像,误差的衰减是在全局范围内的统计,对算法过程中相幅矢量的空间分布进行统计,可以 更好的评估迭代算法的效果。
在迭代过程中,使用离散型傅里叶变换(DFT)进行空域和频域的相互转化:
Figure RE-GDA0002498205000000081
为了便于分析,使用一维的离散傅里叶变换(DFT)得到:
Figure RE-GDA0002498205000000082
根据欧拉公式,
Figure RE-GDA0002498205000000083
对于输出面X(k),其相幅矢量 的实部和虚部可以分别表示为:
Figure RE-RE-GDA0002498205000000084
其中,A(n)为输入面上第n个像素上相幅矢量的幅值,A(n)为常数,对于G-S算法,输入面上会施加约束。φn为第n个像素上相幅矢量的相位角,在算法中,初始输入相位 采用随机相位,φn满足区间[0,2π]上的均匀分布。
假设输入面的相幅矢量满足:
(1)所有的振幅和相位相互统计独立;
(2)当m≠n时,第n个单元的振幅和相位与第m个单元的振幅相位相互独立。
根据中心极限定理,N·R(n)~N(nμ,nσ2)、N·I(n)~N(nμ,nσ2)。其中,
Figure RE-RE-GDA0002498205000000085
其中,取A(n)=1,因此,DFT变换后的频域内的相幅矢量的实部和虚部均满足正态分布N(0,σ2/N),频域内实部和虚部的联合概率密度可以表示为:
Figure RE-GDA0002498205000000086
根据R2+I2=A2,对上式作积分变换,得到频域上相幅矢量幅值和相位的联合概率密度分布:
Figure RE-GDA0002498205000000087
因此,频域上相幅矢量幅值和相位各自的概率密度函数为:
Figure RE-GDA0002498205000000091
由上式可知,频域内的相幅矢量的幅值满足瑞利分布,相位满足均匀分布。
如图6所示,输入面大小为500×600,选用幅值为1的约束条件,初始相位为程序生成的随机数,满足[0,2π]内的均匀分布,进行DFT变换到频域内后,其所得的相幅矢 量的幅值和相位的分布分别如图6中(a)(b)所示,(c)为迭代200次后衍射面的 幅值分布,幅值分布分为两部分,分别对应着背景区域与目标图像区域。
根据G-S算法收敛性的分析可知,迭代算法是通过三角形两边之和大于第三边的来 实现误差的减小,这使得算法的误差衰减依赖于初始相幅与误差相幅的夹角,如图7中角α,当夹角α均匀分布时,且ε较小时,有ε′=ε·cosα,对整个像面进行统计,有误差 的衰减系数k:
Figure RE-GDA0002498205000000092
可见,误差ε随着迭代次数的增加以
Figure RE-GDA0002498205000000093
衰减,假设α在迭代过程中始终均匀分布,经过一定次数的迭代后,则有ε→0。
如图8所示,为500次相位迭代算法过程中误差的变化,从图中可见,在20次迭代运算后,误差衰减逐渐变小,趋于不变。即在迭代运算过程中,夹角α并不是均匀分布 的,而是逐渐减小,从而使得误差的衰减速度下降。
如图9所示,对夹角α进行测量,随着迭代次数的增加,α主要集中在0π、π与2π处,此时,迭代前后的误差几乎没有发生变化,算法停滞。
进一步地,对迭代算法中误差的衰减进行分析可知,在目标区域,误差通过三角形原理逐步减小,而在背景区域,由于幅值约束为0,迭代前后误差仍保持不变,该背景 区域未参与误差衰减过程,在一定程度上限制了算法误差的衰减,当背景在目标图像中 的占比越高时,这种误差衰减的限制就越大,因此需要对上面的迭代算法进行改进,以 提高成像的光利用率和减小迭代算法的误差。
下面参考附图描述本发明第一方面实施例的光相位调制方法,该算法通过增加目标 图像的背景值,可以增大迭代算法的误差衰减,进而可以使得迭代算法获得的相息图有着更高的均一性和更小的误差分布。
图10是根据本发明一个实施例的一种光相位调制方法的流程图。
如图10所示,光相位调制方法包括步骤S1、步骤S2、步骤S3、步骤S4和步骤S5。
步骤S1-傅里叶变换,对输入光束复振幅分布进行傅里叶变换以获得输入光束在频 域上的复振幅分布。
具体地,在迭代算法的初始输入,还包括函数初始化步骤,赋予初始输入的输入光束初始相位分布,并由计算机设备随机生成二维矩阵,以及根据实际输入光束的分布确 定输入光束复振幅分布的幅值,也可以将输入光束的复振幅分布幅值定义为平面波函数 分布或高斯光束波面函数分布。以及,对初始化的输入光束复振幅分布进行离散傅里叶 变换,得到频域上的复振幅分布。
步骤S2-第一幅值替换,将输入光束在频域上的复振幅分布的幅值替换为目标光场 复振幅分布的幅值,获得幅值替换后的复振幅分布,其中,目标光场复振幅分布的幅值包括目标图像及背景值,且背景值不为零,即加入背景值。
步骤S3-傅里叶逆变换,对幅值替换后的复振幅分布进行傅里叶逆变换,获得输入光束在空域上的复振幅分布。
步骤S4-第二幅值替换,将输入光束在空域上的复振幅分布的幅值替换为输入光束 的复振幅分布的幅值,获得空域上的新的输入光束复振幅分布函数。
步骤S5-迭代计算,以傅里叶变换步骤、第一幅值替换步骤、傅里叶逆变换步骤、第二幅值替换步骤进行迭代计算。
由于目前采用的迭代算法,在背景区域幅值约束为0,迭代前后误差保持不变,即背景区域未参与误差衰减过程,在一定程度上限制了算法误差的衰减,当背景在目标图 像中的占比越高时,这种误差衰减的限制就越大。如图11所示,其中(a)和(b)图 像尺寸均为500×500,(a)的背景空间占比为77.7%,(b)的背景空间占比为35.5%。 (c)为目标图像通过500次相位迭代算法后重建图像的误差变化。如图11中(c)所 示,背景空间占比较小的(b)经过迭代运算后获得的相息图有着更小的误差。
由于不同的目标图像有着不同的背景空间占比,而背景空间占比越小的图像经过迭 代后有着更小的误差,因此,背景空间占比的改变对迭代算法的误差有着很重要的影响。 但是对于固定的目标图像往往无法改变其背景空间占比,因此,在本发明实施例中,通过在迭代计算中目标光场复振幅分布的幅值中添加背景值,使背景区域的目标光场复振幅分布的幅值不为0,在迭代的幅值约束过程中,背景区域就能参与误差的衰减过程, 从而使图像的整体误差减小。
如图12所示,通过增加背景的G-S算法与原始G-S算法的比较结果,其中(a)为 算法均一性的比较,(b)为算法的误差衰减的比较,如图所示,通过添加背景后,迭 代算法获得的相息图有着更高的均一性和更小的误差分布。
根据本发明实施例的光相位调制方法,在迭代算法中的第一幅值替换步骤中,通过 目标光场复振幅分布的幅值中添加背景值,使得背景值参与迭代,可以增大迭代算法的误差衰减,进而可以使得迭代算法获得的相息图有着更高的均一性和更小的误差分布。 经过多次迭代,使输入光束在频域上的复振幅分布更加接近于目标复振幅分布。
在一些实施例中,可以根据噪声幅值尺度设置目标图像的背景值。
如图13所示,在不同背景值下的算法效果,其中,背景值分别取目标图像幅值的0.01倍、0.1倍、0.2倍。如图13所示,随着背景值的增加,迭代算法的均一性随之增 加。但是过低的背景值并不会明显增加迭代算法的误差衰减和重建图像的均一性,这是 由于误差的幅值范围相比背景值较大,迭代前后误差变化并不显著。另外,目标光场复 振幅分布的幅值中背景值的增加,也会使衍射图案背景区域有着相应的光强分布,这在 一定程度上会降低成像的对比度。
选取过低的背景值会使背景区域在误差衰减过程中贡献下降,而过高的背景值虽然 会提高算法的误差衰减,但同时也会增加生成图像中的背景值,对光调控的实验结果产生不利的影响。因此,为了获得较好的算法结果,应将目标图像的背景值设置与噪声幅 值尺度相当,合理折中选择背景的取值,在本发明实施例中,背景值选取目标图像的幅 值的0.1倍-0.2倍为最优范围,例如背景值选取目标图像幅值的0.1倍、0.11倍、0.13 倍、0.14倍、0.16倍、0.18倍、0.2倍等。
在一些实施例中,为了减少因引入图像背景值而造成的给光场分布应用带来不利的 影响,可以通过图像掩模的方式,在掩模支持域内和支持域外进行不同的幅值约束,从而在不增加背景值的条件下,降低迭代算法的误差。
在实施例中,为了在目标图像中不引入新的背景值,首先将目标图像添加在新的背 景板上成为新的目标图像,初始的目标图像的尺寸为nx×ny,新的背景板的尺寸为 Nx×Ny,其中,Nx>nx,Ny>ny,在新的目标图像上制作图像掩模,如图14所示, 在迭代运算过程中,掩模支持域内的幅值约束以目标图像为基准,而掩模支持域外的背 景区域,则使用非零的幅值作为约束进行迭代,及在第一幅值替换步骤中,对应于图像 掩模支持域内和图像掩模支持域外采用不同的背景值。通过这种方法,可以不增加目标 图像的背景,并获得更小的迭代误差。
在一些实施例中,在第一幅值替换步骤中,将输入光束在频域上的复振幅分布的幅 值替换为目标光场复振幅分布的幅值,获得幅值替换后的复振幅分布,目标光场复振幅分布的幅值为包含背景值的新的目标图像,其中,对应于图像掩模支持域内,背景值为 零,对应于图像掩模支持域外,背景值不为零。
例如,如图15所示,取
Figure RE-GDA0002498205000000121
掩模支持域外的背景幅值取0.1,通过500次 迭代运算后得到的重建图像。其中(a)为全幅图像,其中,中间区域为支持域区域, (b)展示了支持域区内的迭代结果,(c)为加权G-S算法在500次迭代运算下的恢复 结果。如图所示,掩模支持下的G-S算法获得了背景区域相较加权G-S算法要更加干净。
图16是掩模支持域内的幅值分布情况,如图16所示,其中背景区域的幅值取值要小于G-S迭代算法的背景幅值范围,且目标区域幅值分布更为集中。
图17是将改进的相位迭代算法的衍射效率、误差函数和重建图像均一性与G-S算法 进行比较后的效果图,如图17所示,改进的相位迭代算法相比G-S算法衍射效率从70%上升至90%,平均误差从0.012降低至0.002,重建的图像均一性也从0.7升高至0.9, 根据本发明的实施例,通过该方法能整体上提高重建图像的质量。
图18是对掩模支持域外背景幅值的不同取值进行的比较效果图,如图18所示,背景值分别取值目标图像幅值的0.01倍、0.1倍、0.2倍、1倍。在500次迭代运算后, 重建图像的衍射效率、均一性及误差函数均随着背景值的增加而改善,因此,可以通过 改变不同的掩模背景值来改善迭代算法的收敛情况。但在实际应用中,掩模支持域外的 部分由于背景值得添加,在调制后会将一部分光能量分散到该区域,较大的掩模背景值 虽可以提高目标图像的质量,但会在一定程度上降低系统整体的光利用率,因此掩模背 景值不宜选取太高,经过试验发现想要效果达到最优效果,掩模背景值最好大于0.1。
在一些实施例中,利用加权G-S算法中对幅值约束的思想,在掩模迭代算法中对幅值约束进行改进。在图像掩模支持域内,目标光场复振幅分布的幅值|Vu,v|满足以下约 束条件:
Figure RE-GDA0002498205000000122
其中,E为单位矩阵,
Figure RE-GDA0002498205000000123
为目标光场复振幅分布的幅值,和
Figure RE-GDA0002498205000000124
为迭代计算中 归一化的输入光束频域复振幅分布的幅值,WkWk和Wk-1为迭代结果,k为迭代次数。 在支持域外,幅值取设定的背景值为目标图像幅值的0.1倍。
如图19所示,为图像掩模算法对加权G-S算法的优化,通过图像掩模的引入,新的算法不仅保持了加权G-S算法能获得高均一性图像的同时,也提高了算法的衍射效率, 降低了误差。另一方面,在掩模支持下的加权G-S算法在图像的衍射效率、误差衰减以 及均一性上收敛速度均更快,且在图像的均一性上避免了加权G-S算法在初始阶段的振 荡过程,算法更加稳定。
进一步地,在迭代计算过程中,假设图像掩模支持域外的区域没有误差,即在迭代过程中保留支持域外幅值的迭代运算结果,仅对图像掩模支持域内的目标图像部分进行幅值约束,这样可以进一步的减小迭代过程中的误差。
在迭代过程中,所述第一幅值替换步骤在频域内对所述图像掩模支持域内的幅值A′n采用以下约束条件:A′n=A0·mask+An·(1-mask),其中,A0为目标图像,mask为 图像掩模,An为第n次迭代的结果。为了更好地比较不同迭代算法的恢复效果,定义图 像的对比度β:
Figure RE-GDA0002498205000000131
其中,对比度β为迭代结果中目标区域与背景区域幅值均值之比。
在迭代算法中,越高的对比度意味着更好地重建图像质量,有着更低的背景噪声。如图20所示,为不同算法在迭代过程中对比度的变化,随着迭代的增加,加权G-S算 法和标准G-S算法仅能达到10左右的对比度,而掩模支持下的G-S算法和加权G-S算 法均能提高一定的图像对比度。而保留迭代过程中掩模外背景的G-S算法,则能获得60 以上的图像对比度,大大提高了重建图像的质量。
图21为500次迭代后目标图像区域的幅值分布的示意图,如图21所示,相对于一般的图像掩模G-S算法,背景幅值继续被压缩,集中在一个很小的范围内,目标区域幅 值相比掩模G-S算法有些许下降,但幅值分布更为集中,说明图像均一性更高。同样的, 掩模支持域内的幅值同样可以使用加权约束的方法以获得更高的图像均一性。
根据本发明的实施例,通过图像掩模的方法,可以有效的抑制结果图像中支持域内 背景区域的幅值分布,降低了迭代算法的噪声,而支持域外的区域则可使用光阑,使中间区域的光束通过空间滤波的方法滤除,不会对成像带来影响。
下面参考附图描述本发明第二方面实施例的一种任意位置和形状照射的空间光调制 方法,该方法根据目标光场相位分布函数对输入光束进行调制,可以减小光束成像过程 中的误差,提升光能利用率和图像衍射效率。
图22是根据本发明的一个实施例的一种任意位置和形状照射的空间光调制方法的 流程图。
如图22所示,光相位调制方法包括步骤S1和步骤S2。
步骤S10-根据上面实施例的光相位调制方法获取对输入光束进行光相位调制的调 制相位分布函数。
其中,光相位调制方法的实现过程可以参照上面实施例的说明,其中,当输入光束在频域上的复振幅分布与目标光场的复振幅分布小于规定的误差值,或者,迭代计算的 次数达到预设次数时,确定的傅里叶逆变换步骤所获得的输入光束复振幅分布中的相位 分布即为对输入光束进行相位调制所需的相位分布函数。
步骤S20-根据调制相位分布函数,利用相位型空间光调制器对输入光束进行调制, 在衍射平面获得任意位置及形状的目标光强分布。
具体地,在实施例中,可以采用激光光源、准直扩束镜、空间光调制器、偏振片、 傅里叶透镜及数据处理装置组成任意位置和形状照射的空间光调制系统,其中,数据处 理装置执行上面实施例的光相位调制方法,获得所需的调制相位分布函数,进而通过相 位型空间光调制器对输入光束进行调制,即可在衍射平面上获得目标光强分布,其中, 可以参照需求根据调制相位分布函数在衍射平面的任意位置获得任意形状的目标光强 分布。
根据本发明实施例的血流成像空间光调制方法,通过采用上面实施例的光相位调制 方法计算全息光场分布,获得所需的光场强度分布,可以增加光场分布的均匀性,提高光能利用率。
在本说明书的描述中,参考术语“一个实施例”、“一些实施例”、“示意性实施 例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描 述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说 明书中,对上述术语的示意性表述不一定指的是相同的实施例或示例。
尽管已经示出和描述了本发明的实施例,本领域的普通技术人员可以理解:在不脱 离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型,本发明的范围由权利要求及其等同物限定。

Claims (10)

1.一种光相位调制方法,其特征在于,包括:赋予输入光束的复振幅分布um,n,包括其幅值分布|Um,n|,及其初始相位分布φm,n
傅里叶变换步骤,对输入光束复振幅分布um,n进行傅里叶变换以获得输入光束在频域上的复振幅分布
Figure RE-FDA0002395722060000011
第一幅值替换步骤,将所述输入光束在频域上的复振幅分布的幅值
Figure RE-FDA0002395722060000012
替换为目标光场复振幅分布的幅值|Au,v|,获得幅值替换后的复振幅分布f′u,v=|Au,v|exp(iφu,v),其中,所述目标光场复振幅分布的幅值包括目标图像及背景值,且所述背景值不为零;
傅里叶逆变换步骤,对幅值替换后的复振幅分布f′u,v进行傅里叶逆变换,获得所述输入光束在空域上的复振幅分布
Figure RE-FDA0002395722060000013
第二幅值替换步骤,将所述输入光束在空域上的复振幅分布的幅值
Figure RE-FDA0002395722060000014
替换为所述输入光束的复振幅分布的幅值|Um,n|,获得空域上的新的输入光束复振幅分布g′m,n=|Um,n|exp(iφ′m,n);
以所述傅里叶变换步骤、所述第一幅值替换步骤、所述傅里叶逆变换步骤、所述第二幅值替换步骤进行迭代计算,获得对输入光束进行相位调制所需的相位分布函数φ′m,n
2.根据权利要求1所述的光相位调制方法,其特征在于,为获得对输入光束进行相位调制所需的相位分布函数:
在迭代过程中,当所述输入光束在频域上的复振幅分布与目标光场的复振幅分布小于规定的误差值,或者,迭代计算的次数达到预设次数时,确定所述傅里叶逆变换步骤所获得的输入光束在空域上的复振幅分布中的相位分布φ′m,n为对输入光束进行相位调制所需的相位分布函数。
3.根据权利要求1所述的光相位调制方法,其特征在于,将所述目标图像添加在新的背景板上成为新的目标图像,初始的目标图像的尺寸为nx×ny,新的背景板的尺寸为Nx×Ny,其中,Nx>nx,Ny>ny,在新的目标图像上制作图像掩模,所述图像掩模的尺寸与初始的所述目标图像的尺寸相当,其中,在所述第一幅值替换步骤中,对应于所述图像掩模支持域内和所述图像掩模支持域外采用不同的背景值。
4.根据权利要求3所述的光相位调制方法,其特征在于,所述第一幅值替换步骤包括:
将所述输入光束在频域上的复振幅分布的幅值替换为目标光场复振幅分布的幅值,获得幅值替换后的复振幅分布,所述目标光场复振幅分布的幅值为包含背景值的新的目标图像,其中,对应于所述图像掩模支持域内,所述背景值为零,对应于所述图像掩模支持域外,所述背景值不为零。
5.根据权利要求4所述的光相位调制方法,其特征在于,在所述图像掩模支持域内,所述目标光场复振幅分布的幅值|Vu,v|满足以下约束条件:
Figure RE-FDA0002395722060000021
其中,
E为单位矩阵,
Figure RE-FDA0002395722060000022
为目标光场复振幅分布的幅值,
Figure RE-FDA0002395722060000023
为迭代计算中归一化的输入光束频域复振幅分布的幅值,Wk和wk-1为迭代结果,k为迭代次数。
6.根据权利要求4所述的光相位调制方法,其特征在于,所述光相位调制方法还包括:
在进行迭代计算时,保留所述图像掩模支持域外幅值的迭代计算结果,仅对所述图像掩模支持域内的目标图像部分进行幅值约束,其中,在迭代过程中,所述第一幅值替换步骤在频域内对所述图像掩模支持域内的幅值A'n采用以下约束条件:
A′n=A0·mask+An·(1-mask),其中,
A0为目标图像幅值,mask为图像掩模幅值,An为第n次迭代的幅值结果。
7.根据权利要求1所述的光相位调制方法,其特征在于,还包括:
函数初始化步骤,赋予初始输入光束初始相位分布,并生成二维矩阵,以及根据所述输入光束确定所述输入光束的复振幅函数分布或者输入光束采用平面波函数分布或高斯光束波面函数分布。
8.根据权利要求1所述的光相位调制方法,其特征在于,其中,根据噪声幅值设置所述目标图像的背景值。
9.根据权利要求7所述的光相位调制方法,其特征在于,所述目标图像的背景值的取值范围为目标图像幅值的0.1倍-0.2倍。
10.一种任意位置和形状照射的空间光调制方法,其特征在于,包括:
根据权利要求1-9任一项所述的光相位调制方法获取对输入光束进行光相位调制的相位分布函数φ′m,n
根据所述相位分布函数,利用相位型空间光调制器对输入光束进行调制,在衍射平面获得任意位置及形状的目标光强分布。
CN201911277585.0A 2019-12-11 2019-12-11 光相位调制方法及任意位置和形状照射的空间光调制方法 Active CN111340902B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911277585.0A CN111340902B (zh) 2019-12-11 2019-12-11 光相位调制方法及任意位置和形状照射的空间光调制方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911277585.0A CN111340902B (zh) 2019-12-11 2019-12-11 光相位调制方法及任意位置和形状照射的空间光调制方法

Publications (2)

Publication Number Publication Date
CN111340902A true CN111340902A (zh) 2020-06-26
CN111340902B CN111340902B (zh) 2023-08-18

Family

ID=71183268

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911277585.0A Active CN111340902B (zh) 2019-12-11 2019-12-11 光相位调制方法及任意位置和形状照射的空间光调制方法

Country Status (1)

Country Link
CN (1) CN111340902B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113485011A (zh) * 2021-07-27 2021-10-08 南通智能感知研究院 一种基于模拟退火法的软边光阑整形相位优化方法
WO2023109428A1 (zh) * 2021-12-16 2023-06-22 中国科学院上海高等研究院 一种全光场成像相机及其成像方法及全光场成像装置

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120116703A1 (en) * 2009-04-24 2012-05-10 Nicolas Pavillon Method and apparatus for enhanced spatial bandwidth wavefronts reconstructed from digital interferograms or holograms
CN104035313A (zh) * 2014-05-26 2014-09-10 北京理工大学 一种精确强度调制全息显示的全息图编码方法
CN105651399A (zh) * 2016-01-15 2016-06-08 中国科学院上海光学精密机械研究所 时域相位恢复全光纤激光脉冲弱相位测量装置和测量方法
US20180217629A1 (en) * 2015-09-28 2018-08-02 Cambridge Enterprise Limited Optical system for performing complex fourier transforms
CN109565334A (zh) * 2016-08-22 2019-04-02 三菱电机株式会社 用于检测光纤传输网络中的波长分离器处的失谐的处理和光纤传输网络

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120116703A1 (en) * 2009-04-24 2012-05-10 Nicolas Pavillon Method and apparatus for enhanced spatial bandwidth wavefronts reconstructed from digital interferograms or holograms
CN104035313A (zh) * 2014-05-26 2014-09-10 北京理工大学 一种精确强度调制全息显示的全息图编码方法
US20180217629A1 (en) * 2015-09-28 2018-08-02 Cambridge Enterprise Limited Optical system for performing complex fourier transforms
CN105651399A (zh) * 2016-01-15 2016-06-08 中国科学院上海光学精密机械研究所 时域相位恢复全光纤激光脉冲弱相位测量装置和测量方法
CN109565334A (zh) * 2016-08-22 2019-04-02 三菱电机株式会社 用于检测光纤传输网络中的波长分离器处的失谐的处理和光纤传输网络

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
李兴龙: "无阈值窗口傅里叶变换滤波方法及涡旋光在变形测量中的应用" *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113485011A (zh) * 2021-07-27 2021-10-08 南通智能感知研究院 一种基于模拟退火法的软边光阑整形相位优化方法
WO2023109428A1 (zh) * 2021-12-16 2023-06-22 中国科学院上海高等研究院 一种全光场成像相机及其成像方法及全光场成像装置

Also Published As

Publication number Publication date
CN111340902B (zh) 2023-08-18

Similar Documents

Publication Publication Date Title
US10949954B2 (en) Imaging system and method for imaging objects with reduced image blur
US11137719B2 (en) Methods, systems, and computer readable media for improved digital holography and display incorporating same
US8254714B2 (en) Methods and systems for designing electromagnetic wave filters and electromagnetic wave filters designed using same
CN111340902A (zh) 光相位调制方法及任意位置和形状照射的空间光调制方法
CN114690404B (zh) 一种基于超表面的多衍射级独立光波场调控方法
CN108885425A (zh) 用于编码重建三维物体的复值信号的方法和装置
Alsaka et al. A comparison of iterative Fourier transform algorithms for image quality estimation
Ye et al. Robust phase retrieval with green noise binary masks
CN114002931A (zh) 基于深度学习加速计算的大视场全息投影方法及系统
EP3535554A1 (en) Wavefront sensor and method of reconstructing distorted wavefronts
CN112697751B (zh) 多角度照明无透镜成像方法、系统及装置
Whitehead et al. Fresnel diffractive imaging: Experimental study of coherence and curvature
CN114355743A (zh) 一种多衍射级独立光波场调控的全息编码方法
CN109814258B (zh) 提高光束光强均匀度的复振幅整形方法
Konyaev Computer simulation of adaptive optics for laser systems in atmospheric applications
Hu et al. Mixed multi-region design of diffractive optical element for projection exposure system
Guo et al. Speckle-reduced diffractive optical elements beam shaping with regional padding algorithm
Younis et al. Evolving Deformable Mirror Control to Generate Partially Coherent Light Fields
CN114494596B (zh) 一种计算全息图的生成方法及电子设备
Paterson Diffractive optical elements with spiral phase dislocations
Gladrow Digital phase-only holography using deep conditional generative models
CN116309073B (zh) 一种基于深度学习的低对比度条纹sim重建方法和系统
Zhang et al. Diffraction characteristics of a pixelated spatial light modulator for wave field reconstruction
Zavarzin et al. THE INFLUENCE OF THE APERTURE STOP FRACTAL SHAPE OF AN OPTICAL SYSTEM ON THE ILLUMINANCE DISTRIBUTION
Tervonen et al. Pulse-frequency-modulated high-frequency-carrier diffractive elements for pattern projection

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