CN114399448B - 一种基于非下采样剪切波变换的多偏振信息选通融合方法 - Google Patents

一种基于非下采样剪切波变换的多偏振信息选通融合方法 Download PDF

Info

Publication number
CN114399448B
CN114399448B CN202111388499.4A CN202111388499A CN114399448B CN 114399448 B CN114399448 B CN 114399448B CN 202111388499 A CN202111388499 A CN 202111388499A CN 114399448 B CN114399448 B CN 114399448B
Authority
CN
China
Prior art keywords
image
source image
frequency
fusion
neighborhood
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
Application number
CN202111388499.4A
Other languages
English (en)
Other versions
CN114399448A (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.)
XiAn Institute of Optics and Precision Mechanics of CAS
Original Assignee
XiAn Institute of Optics and Precision Mechanics of CAS
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 XiAn Institute of Optics and Precision Mechanics of CAS filed Critical XiAn Institute of Optics and Precision Mechanics of CAS
Priority to CN202111388499.4A priority Critical patent/CN114399448B/zh
Publication of CN114399448A publication Critical patent/CN114399448A/zh
Application granted granted Critical
Publication of CN114399448B publication Critical patent/CN114399448B/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
    • G06T5/00Image enhancement or restoration
    • G06T5/50Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/20Image enhancement or restoration using local operators
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20212Image combination
    • G06T2207/20221Image fusion; Image merging

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Image Processing (AREA)

Abstract

本发明涉及偏振图像融合技术,具体涉及一种基于非下采样剪切波变换的多偏振信息选通融合方法。解决了现有图像融合时偏振角图像参与度不高无法保留有效信息以及偏振角图像噪声滤除较差,造成融合图像信噪比较低的技术问题。本发明方法包括以下步骤:S1)利用非下采样剪切波变换对可见光强度源图像、偏振度源图像、偏振角源图像进行分解得到低频分量A0,B0,C0及高频方向分量
Figure DDA0003367953570000011
S2)利用加权平均法对源图像的低频分量进行融合得到低频融合分量;S3)对源图像的高频方向分量进行选通,取邻域能量最大的像素点灰度值得到融合图像的高频方向分量;S4)对低频融合分量和融合图像高频方向分量进行非下采样剪切波变换的逆变换,得到融合图像。

Description

一种基于非下采样剪切波变换的多偏振信息选通融合方法
技术领域
本发明涉及偏振图像融合技术,具体涉及一种基于非下采样剪切波变换的多偏振信息选通融合方法。
背景技术
物体表面的反射光和散射光会包含其自身特性的偏振信息,通过这些偏振信息,可以反应被测物体的轮廓、纹理、粗糙度等物理信息。偏振度图像及偏振角图像分别根据目标表面偏振度及偏振角进行成像,可实现对目标边缘特征的增强,并对耀光、过曝等影响普通可见光强度图像的因素有较强的抑制作用,有着独特的目标探测能力。
为实现偏振角图像、偏振度图像、可见光强度图像信息的互补以及后续数据利用的便利,需要将上述三种图像进行融合为一。但偏振角图像、偏振度图像中噪声影响较大,尤其偏振角图像,噪声影响较为显著,使得融合过程中想要同时实现滤除噪声和保留源图像信息比较困难。
现有偏振图像融合算法可分为两类:第一类是只利用可见光强度图像和偏振度图像进行融合,噪声影响较大的偏振角图像不参与融合。这类融合方法得到的融合图像信噪比较高,但也同样丢失了偏振角图像中包含的有效信息。第二类是将偏振角图像、偏振度图像、可见光强度图像三者进行融合。这类方法可实现三种源图像信息的融合,但在噪声抑制及信息保留两种目标之间未实现有效的均衡。要么未对偏振角图像中较强的噪声进行有效处理,造成融合图像信噪比较低,要么抑制了噪声,同时保留的有效信息也较少。
发明内容
本发明目的是解决现有偏振角图像、偏振度图像、可见光强度图像三者进行融合时,在噪声抑制及信息保留两种目标之间未实现有效的均衡,造成融合图像质量较差的技术问题,而提供一种基于非下采样剪切波变换的多偏振信息选通融合方法。该方法利用非下采样剪切波变换进行滤波分解,通过邻域方差大小对源图像窗口邻域进行选通,判定是否噪声较大,是否参与融合,可有效抑制源图像中的噪声并保留其有效信息。该方法可对图像各分量进行精细融合调控,实现了多偏振信息的高质量融合。
本发明技术方案为:
一种基于非下采样剪切波变换的多偏振信息选通融合方法,其特殊之处在于,包括以下步骤:
S1)利用非下采样剪切波变换对可见光强度源图像、偏振度源图像、偏振角源图像进行多级多方向分解,得到可见光强度源图像、偏振度源图像、偏振角源图像的低频分量A0,B0,C0及高频方向分量
Figure BDA0003367953550000021
k为非下采样剪切波变换的分解级数,l为各分解级下的方向级数;
所述k为3~4;
所述l为4的倍数;
S2)利用加权求和方法对步骤S1)得到的可见光强度源图像、偏振度源图像、偏振角源图像的低频分量灰度值进行融合得到低频融合分量R0
R0=λ1A02B03C0
λ1,λ2,λ3分别为可见光强度源图像、偏振度源图像、偏振角源图像的低频分量在融合图像中的权重;
所述λ1,λ2,λ3均不为0,且λ123=1,λ2=λ3,λ1为0.5~0.6;
S3)对步骤S1)得到的可见光强度源图像、偏振度源图像、偏振角源图像的高频方向分量中各像素点计算邻域窗口ω1的邻域方差大小,并根据邻域方差大小进行选通;若邻域方差大于阈值,则对应像素点不参与后续融合;若邻域方差小于阈值,则对相应的像素点再取邻域窗口ω2,并计算该邻域窗口ω2的能量,选取可见光强度源图像、偏振度源图像、偏振角源图像中邻域窗口ω2能量最大的像素点灰度值作为融合图像的相应高频方向分量像素点灰度值;各像素点灰度值确定后,便得到融合图像高频方向分量
Figure BDA0003367953550000031
Figure BDA0003367953550000032
Figure BDA0003367953550000033
分别为可见光强度源图像高频方向分量、偏振度源图像高频方向分量、偏振角源图像高频方向分量;
F{}为高频融合算法;
Figure BDA0003367953550000034
为经过高频算法选通后得到的融合图像的高频方向分量;
所述邻域窗口ω1为21~31个像素;
所述邻域窗口ω2为3~5个像素;
所述阈值为图像整体邻域方差平均值的1.2~1.5倍;
S4)对步骤S2)得到的低频融合分量R0及步骤S3)得到的融合图像的高频方向分量
Figure BDA0003367953550000035
进行非下采样剪切波变换的逆变换,得到最终的融合图像R。
进一步地,步骤S3)中,方差Var计算公式为:
Figure BDA0003367953550000036
MV为邻域灰度均值;
(i,j)为邻域坐标;
g(i,j)为图像(i,j)处的灰度值;
M×N为窗口大小;
邻域能量Ar计算公式为:
Figure BDA0003367953550000041
ω(p,q)为邻域窗口各像素权重;
g(i+p,j+q)为图像(i+p,j+q)处的灰度值;
窗口横向宽度为2m+1,纵向宽度为2n+1。
进一步地,步骤S1)中,k=4;l在各分解级上分别为8,8,16,16;
步骤S2)中,λ1,λ2,λ3分别取0.5,0.25,0.25。
进一步地,步骤S3)中,所述邻域窗口ω1为21×21;
所述阈值为图像整体邻域方差平均值的1.2倍;
所述邻域窗口ω2为3×3,邻域窗口各像素权重(ω(p,q))为1/16*[1 2 1;2 4 2;1 2 1]。
本发明的有益效果:
1、该方法利用非下采样剪切波变换进行滤波分解,通过邻域方差大小对源图像高频方向分量窗口邻域进行选通,判定是否噪声较大,是否参与融合,可有效抑制源图像中的噪声并保留其有效信息。
2、该方法可对各源图像高频方向分量进行精细融合调控,实现了多偏振信息的高质量融合。
3、该方法能够较好的保留三种源图像的重要特征信息实现了融合图像的高质量表达。
附图说明
图1为本发明实施例的偏振度源图像;
图2为本发明实施例的偏振角源图像;
图3为本发明实施例的可见光强度源图像;
图4为本发明实施例的融合图像。
具体实施方式
一种基于非下采样剪切波变换的多偏振信息选通融合方法,包括以下步骤:
S1)利用非下采样剪切波变换对可见光强度源图像、偏振度源图像、偏振角源图像进行多级多方向分解,得到可见光强度源图像、偏振度源图像、偏振角源图像的低频分量A0,B0,C0及高频方向分量
Figure BDA0003367953550000051
k为非下采样剪切波变换的分解级数,l为各分解级下的方向级数。分解级数太少,融合不够精细;分解级数太多,计算量太大,因此k一般为3~4。同样分解方向太少,不利用融合分析;分解方向太多,计算量太大,l为4的倍数,一般取8或16。本实施例中分解级数k=4,在各分解级上的方向级数l分别为8,8,16,16。可见光强度源图像、偏振度源图像、偏振角源图像均被分解为1幅低频分量图像及48幅高频方向分量图像。
S2)利用加权求和方法对步骤S1)得到的可见光强度源图像、偏振度源图像、偏振角源图像的低频分量灰度值进行融合得到低频融合分量R0
R0=λ1A02B03C0
λ1,λ2,λ3分别为可见光强度源图像、偏振度源图像、偏振角源图像的低频分量在融合图像中的权重;λ1,λ2,λ3均不为0,且λ123=1,λ2=λ3,λ1为0.5~0.6。本实施例λ1,λ2,λ3分别取0.5,0.25,0.25。
S3)先对步骤S1)得到的可见光强度源图像、偏振度源图像、偏振角源图像的高频方向分量中各像素点计算邻域窗口ω1的邻域方差大小,并根据邻域方差大小进行选通;若邻域方差大于阈值,则对应像素点不参与后续融合;若邻域方差小于阈值,则对相应的像素点再取邻域窗口ω2,并计算该邻域窗口ω2的能量,选取可见光强度源图像、偏振度源图像、偏振角源图像中邻域窗口ω2能量最大的像素点灰度值作为融合图像的相应高频方向分量像素点灰度值;各像素点灰度值确定后,便得到融合图像高频方向分量
Figure BDA0003367953550000061
Figure BDA0003367953550000062
Figure BDA0003367953550000063
为各源图像高频方向分量;
F{}为高频融合算法;
Figure BDA0003367953550000064
为经过高频算法选通后得到的融合图像高频方向分量。
方差Var计算公式为:
Figure BDA0003367953550000065
MV为邻域灰度均值;
(i,j)为邻域坐标;
g(i,j)为图像(i,j)处的灰度值;
M×N为窗口大小。
邻域能量Ar计算公式为:
Figure BDA0003367953550000066
ω(p,q)为邻域窗口各像素权重,权重大小一般以窗口中间像素为主,兼顾窗口边缘;
g(i+p,j+q)为图像(i+p,j+q)处的灰度值;
窗口横向宽度为2m+1,纵向宽度为2n+1。
阈值为图像整体邻域方差平均值的1.2~1.5倍,阈值较大,各源图像噪声部分易引入到融合图像中,阈值较小,各源图像的有效信息可能被滤除,本实施例的选通阈值设为图像整体邻域方差平均值的1.2倍。为区分像素点区域是否噪声较大,窗口ω1的取值一般在21~31个像素之间,而窗口ω2一般取3~5个像素,其是在窗口ω1内且噪声较小的像素点区域部分,选取较小领域内质量最好的像素值进行融合。本实施例中邻域窗口ω1为21×21,计算能量的邻域窗口ω2为3×3,m为1,n为1。邻域窗口ω2的各像素权重为1/16*[1 2 1;2 4 2;1 2 1]。
S4)对步骤S2)得到的低频融合分量R0及步骤S3)得到的融合图像的高频方向分量
Figure BDA0003367953550000071
进行非下采样剪切波变换的逆变换,得到最终的融合图像R。
图1、图2和图3分别为室外同一远景的偏振度源图像、偏振角源图像和可见光强度源图像。如图1所示,偏振度源图像由于玻璃镜面的强偏振反射,偏振程度显著,造成图像局部亮度高,整体亮度较低,近处树木细节已无法分辨,但对楼房、窗口等目标的边缘反应较好。如图2所示,偏振角源图像噪声较大,树木处噪声尤其显著,但对远处电线及楼宇窗口可清晰感知。如图3所示,可见光强度源图像,远处天空出现过曝情况,电线处已经无法分辨,但近处树木边缘较清晰。
如图4所示,经过本发明融合算法处理后得到融合图像,融合图像对偏振角源图像中的噪声进行了有效选通抑制,而且树木、楼房、窗口、远处电线等目标的有效信息经过融合均得到有效保留,证明本方法可有效实现对源图像的噪声滤除及信息保留。
此外,本发明方法将可见光强度源图像、偏振度源图像、偏振角源图像进行了多层多方向分解,最终各源图像均被分解为1幅低频分量图像及48幅高频方向分量图像,实现了精细融合操作。

Claims (4)

1.一种基于非下采样剪切波变换的多偏振信息选通融合方法,其特征在于,包括以下步骤:
S1)利用非下采样剪切波变换对可见光强度源图像、偏振度源图像、偏振角源图像进行多级多方向分解,得到可见光强度源图像、偏振度源图像、偏振角源图像的低频分量A0,B0,C0及高频方向分量
Figure FDA0003367953540000011
k为非下采样剪切波变换的分解级数,l为各分解级下的方向级数;
所述k为3~4;
所述l为4的倍数;
S2)利用加权求和方法对步骤S1)得到的可见光强度源图像、偏振度源图像、偏振角源图像的低频分量灰度值进行融合得到低频融合分量R0
R0=λ1A02B03C0
λ1,λ2,λ3分别为可见光强度源图像、偏振度源图像、偏振角源图像的低频分量在融合图像中的权重;
所述λ1,λ2,λ3均不为0,且λ123=1,λ2=λ3,λ1为0.5~0.6;
S3)对步骤S1)得到的可见光强度源图像、偏振度源图像、偏振角源图像的高频方向分量中各像素点计算邻域窗口ω1的邻域方差大小,并根据邻域方差大小进行选通;若邻域方差大于阈值,则对应像素点不参与后续融合;若邻域方差小于阈值,则对相应的像素点再取邻域窗口ω2,并计算该邻域窗口ω2的能量,选取可见光强度源图像、偏振度源图像、偏振角源图像中邻域窗口ω2能量最大的像素点灰度值作为融合图像的相应高频方向分量像素点灰度值;各像素点灰度值确定后,便得到融合图像高频方向分量
Figure FDA0003367953540000012
Figure FDA0003367953540000021
Figure FDA0003367953540000022
分别为可见光强度源图像高频方向分量、偏振度源图像高频方向分量、偏振角源图像高频方向分量;
F{}为高频融合算法;
Figure FDA0003367953540000023
为经过高频融合 算法选通后得到的融合图像的高频方向分量;
所述邻域窗口ω1为21~31个像素;
所述邻域窗口ω2为3~5个像素;
所述阈值为图像整体邻域方差平均值的1.2~1.5倍;
S4)对步骤S2)得到的低频融合分量R0及步骤S3)得到的融合图像的高频方向分量
Figure FDA0003367953540000024
进行非下采样剪切波变换的逆变换,得到最终的融合图像R。
2.根据权利要求1所述的一种基于非下采样剪切波变换的多偏振信息选通融合方法,其特征在于:
步骤S3)中,方差Var计算公式为:
Figure FDA0003367953540000025
MV为邻域灰度均值;
(i,j)为邻域坐标;
g(i,j)为图像(i,j)处的灰度值;
M×N为窗口大小;
邻域能量Ar计算公式为:
Figure FDA0003367953540000026
ω(p,q)为邻域窗口各像素权重;
g(i+p,j+q)为图像(i+p,j+q)处的灰度值;
窗口横向宽度为2m+1,纵向宽度为2n+1。
3.根据权利要求2所述的一种基于非下采样剪切波变换的多偏振信息选通融合方法,其特征在于:步骤S1)中,k=4;l在各分解级上分别为8,8,16,16;
步骤S2)中,λ1,λ2,λ3分别取0.5,0.25,0.25。
4.根据权利要求3所述的一种基于非下采样剪切波变换的多偏振信息选通融合方法,其特征在于:步骤S3)中,所述邻域窗口ω1为21×21;
所述阈值为图像整体邻域方差平均值的1.2倍;
所述邻域窗口ω2为3×3,邻域窗口各像素权重(ω(p,q))为1/16*[1 2 1;2 4 2;1 21]。
CN202111388499.4A 2021-11-22 2021-11-22 一种基于非下采样剪切波变换的多偏振信息选通融合方法 Active CN114399448B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111388499.4A CN114399448B (zh) 2021-11-22 2021-11-22 一种基于非下采样剪切波变换的多偏振信息选通融合方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111388499.4A CN114399448B (zh) 2021-11-22 2021-11-22 一种基于非下采样剪切波变换的多偏振信息选通融合方法

Publications (2)

Publication Number Publication Date
CN114399448A CN114399448A (zh) 2022-04-26
CN114399448B true CN114399448B (zh) 2023-04-11

Family

ID=81225849

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111388499.4A Active CN114399448B (zh) 2021-11-22 2021-11-22 一种基于非下采样剪切波变换的多偏振信息选通融合方法

Country Status (1)

Country Link
CN (1) CN114399448B (zh)

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103530862A (zh) * 2013-10-30 2014-01-22 重庆邮电大学 基于nsct的邻域特性区域化的红外与微光图像融合方法
CN104732504A (zh) * 2015-01-23 2015-06-24 天津大学 基于压缩感知和wbct变换的图像融合方法
CN105139367A (zh) * 2015-07-27 2015-12-09 中国科学院光电技术研究所 一种基于非下采样剪切波的可见光偏振图像融合方法
WO2018120936A1 (en) * 2016-12-27 2018-07-05 Zhejiang Dahua Technology Co., Ltd. Systems and methods for fusing infrared image and visible light image
CN108492274A (zh) * 2018-04-03 2018-09-04 中国人民解放军国防科技大学 一种长波红外偏振特征提取与融合的图像增强方法
CN109658371A (zh) * 2018-12-05 2019-04-19 北京林业大学 红外图像与可见光图像的融合方法、系统及相关设备
CN109816618A (zh) * 2019-01-25 2019-05-28 山东理工大学 一种基于自适应阈值的区域能量光子计数图像融合算法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103530862A (zh) * 2013-10-30 2014-01-22 重庆邮电大学 基于nsct的邻域特性区域化的红外与微光图像融合方法
CN104732504A (zh) * 2015-01-23 2015-06-24 天津大学 基于压缩感知和wbct变换的图像融合方法
CN105139367A (zh) * 2015-07-27 2015-12-09 中国科学院光电技术研究所 一种基于非下采样剪切波的可见光偏振图像融合方法
WO2018120936A1 (en) * 2016-12-27 2018-07-05 Zhejiang Dahua Technology Co., Ltd. Systems and methods for fusing infrared image and visible light image
CN108492274A (zh) * 2018-04-03 2018-09-04 中国人民解放军国防科技大学 一种长波红外偏振特征提取与融合的图像增强方法
CN109658371A (zh) * 2018-12-05 2019-04-19 北京林业大学 红外图像与可见光图像的融合方法、系统及相关设备
CN109816618A (zh) * 2019-01-25 2019-05-28 山东理工大学 一种基于自适应阈值的区域能量光子计数图像融合算法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
Infrared polarization image fusion via multi-scale sparse representation and pulse coupled neural network;Jiajia Zhang 等;《Optical Sensing and Imaging Technology》;20191218;全文 *
基于NSST变换的微光偏振图像融合算法;沈薛晨等;《激光杂志》;20191010(第02期);全文 *
基于区域能量的NSST域偏振图像融合算法;朱达荣等;《安徽建筑大学学报》;20170815(第04期);全文 *
基于可见光偏振成像的目标探测技术研究;刘征;《知网硕士电子期刊》;20160815;全文 *
基于局部能量匹配的红外偏振图像融合;杨风暴等;《红外技术》;20160504(第04期);全文 *
基于边缘信息的偏振图像融合算法及评价;张晶晶等;《光电工程》;20071115(第11期);全文 *

Also Published As

Publication number Publication date
CN114399448A (zh) 2022-04-26

Similar Documents

Publication Publication Date Title
CN109767439B (zh) 一种自适应窗口的多尺度差异与双边滤波的目标检测方法
CN105976330B (zh) 一种嵌入式雾天实时视频稳像方法
Khan et al. Localization of radiance transformation for image dehazing in wavelet domain
CN102800074B (zh) 基于轮廓波变换的sar图像变化检测差异图生成方法
CN111986120A (zh) 一种基于帧累加和多尺度Retinex的低光照图像增强优化方法
CN103578084A (zh) 基于亮通道滤波的彩色图像增强方法
Riaz et al. Single image dehazing via reliability guided fusion
CN107403134B (zh) 基于局部梯度三边的图域多尺度红外弱小目标检测方法
CN103679173A (zh) 图像显著区域检测方法
CN110675340A (zh) 一种基于改进的非局部先验的单幅图像去雾方法及介质
CN109961408B (zh) 基于nsct和块匹配滤波的光子计数图像去噪方法
CN113421205B (zh) 一种结合红外偏振成像的小目标检测方法
Bansal et al. A review of image restoration based image defogging algorithms
CN114332079A (zh) 基于图像处理的塑料餐盒裂缝检测方法、装置及介质
CN114677336A (zh) 一种基于红外图像的幕墙面板损伤识别方法
CN114399449B (zh) 一种基于均值滤波分解的形态选通偏振图像融合方法
CN110599422A (zh) 一种基于边缘保护的加权平均椒盐噪声降噪算法
CN110111280A (zh) 一种多尺度梯度域引导滤波的低照度图像增强算法
CN114399448B (zh) 一种基于非下采样剪切波变换的多偏振信息选通融合方法
KR102618580B1 (ko) 대기광 추정 및 레티넥스 기반의 야간 저조도 영상 복원 방법, 이를 수행하는 장치 및 컴퓨터 프로그램
Zhang et al. A new image filtering method: Nonlocal image guided averaging
CN107203976B (zh) 一种基于噪点检测的自适应非局部均值去噪方法及系统
CN112070765B (zh) 一种基于双边滤波结合改进的otsu的布匹检测方法
CN109961413B (zh) 大气光方向优化估计的图像去雾迭代算法
CN113052833A (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