CN113077380B - 一种栅格数据快速插值方法 - Google Patents

一种栅格数据快速插值方法 Download PDF

Info

Publication number
CN113077380B
CN113077380B CN202110383204.8A CN202110383204A CN113077380B CN 113077380 B CN113077380 B CN 113077380B CN 202110383204 A CN202110383204 A CN 202110383204A CN 113077380 B CN113077380 B CN 113077380B
Authority
CN
China
Prior art keywords
interpolation
point
single cross
interpolated
points
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
CN202110383204.8A
Other languages
English (en)
Other versions
CN113077380A (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.)
Individual
Original Assignee
Individual
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 Individual filed Critical Individual
Priority to CN202110383204.8A priority Critical patent/CN113077380B/zh
Publication of CN113077380A publication Critical patent/CN113077380A/zh
Application granted granted Critical
Publication of CN113077380B publication Critical patent/CN113077380B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00Geometric image transformations in the plane of the image
    • G06T3/04Context-preserving transformations, e.g. by using an importance map

Landscapes

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

Abstract

本发明提供一种栅格数据快速插值方法,包括以下步骤:S1:把待插值图像与插值区域加载到内存;S2:在插值区域内,进行沿着水平垂直方向的单十字内插或加速单十字内插;S3:在插值区域内,进行沿着对角线方向的单十字内插或加速单十字内插;S4:把步骤S2‑S3的预测值逐点求平均,进行双十字内插,得到待插点预测值。本发明在可接受的插值误差下极大提高了计算机处理的实时性,在低性能终端上仅用单核CPU实现了高速内插,还可以在脚本语言上实现实时插值。

Description

一种栅格数据快速插值方法
技术领域
本发明涉及一种遥感影像的自图像插值方法,尤其是涉及一种栅格数据快速插值方法。
背景技术
图像插值算法被广泛应用在航空拍摄、影像预览、数字产品制作等场景中,解决图像拼接、离线插值预览及自动化产品等任务。在现实情况下,受到传输速率、硬件条件以及不同平台的影响,往往需要在受限的算力之下运行插值算法。例如在航空拍摄中,为了能实时获取多个摄像头的拼接影像,需要在性能受限设备上运行插值算法。又或是在网页端浏览影像时,会在图像传输完成之前先用上一级的影像插值出预览图,待传输结束后替换。传统插值算法依赖于GPU进行并行加速,也不适合运行在性能较低的脚本语言上(如网页端上的JavaScript),存在一定的局限性。
基于栅格的数据具有很强的空间关联结构,每一个数据周围的都存在八个相邻的其他数据,而这些相邻数据在内存中的位置可以直接计算获得,若能充分利用这种栅格数据间的空间关联性能够简化索引与遍历数据点的过程。
根据地理学第一定律,任何事物都是与其他事物相关的,而相近的事物关联更紧密。在栅格数据上的未知量可以通过邻近的已知量来表示,这是加权平均插值法的核心思想。通用的搜索邻近点的方法有逐项遍历、构建KD-Tree等,通常都有极大的计算量。栅格数据的空间关联结构极大地简化了寻找邻近已知量的搜索过程,传统方法是借助由小到大递增的偏移量算出邻近点的内存地址,直到获得足够多的邻近点即可结束搜索过程。
IDW(反距离权重插值算法)因其直观、高速、可并行等特点,被广泛应用于图像插值中,但过于依赖GPU加速降低了算法的兼容性与跨平台性。为了在无法调用GPU的终端上使用IDW算法,提出一种IDW简化算法,仅利用栅格数据上待插点的十字方向近邻来插值。为进一步提升十字内插算法速度,提出一种仅遍历两次像元的优化算法。实验结果显示,该算法有效降低计算时间,在低性能终端上可实时运算,适合运用在对实时性有要求的场合。
发明内容
本发明提供了一种栅格数据快速插值方法,解决了插值运算的问题,其技术方案如下所述:
一种栅格数据快速插值方法,包括以下步骤:
S1:把待插值图像与插值区域加载到内存;
S2:在插值区域内,进行沿着水平垂直方向的单十字内插或加速单十字内插;
S3:在插值区域内,进行沿着对角线方向的单十字内插或加速单十字内插;
S4:把步骤S2-S3的预测值逐点求平均,进行双十字内插,得到待插点预测值。
进一步的,步骤S1中,插值区域不能与图像边界直接接触。
进一步的,步骤S2和S3中,针对密集成块的待插值区域,所述加速单十字内插是通过当前待插点的上、下、左、右四个测量值直接计算出预测值,所有待插点中与当前待插点在同一待插连通区域,并处于当前待插点的十字方向的一系列待插点都共享一半的相同测量值。
密集成块的待插值区域进行单十字内插或加速单十字内插时,其过程如下所述:
1)若栅格数据的数据量大小为N*M,且N≤M,设M为图像宽度,加速单十字内插需要开辟长度为2M的整型缓存数组来存放垂直方向上的两个测量值,前M个整型称为垂向起始点,用来缓存对应栅格横坐标的上方向的测量值纵坐标,后M个整型称为垂直终止点,用来缓存对应栅格横坐标的下方向的测量值纵坐标;
2)先初始化缓存数组并置为零,按行依次从左到右进行遍历:当从第一个测量值遍历到待插点时,记录第一个测量值的坐标为区间头,继续向后遍历直到指针从待插点移动到第二个测量值为止,记录第二个测量值的坐标为区间尾,对这个区间内每个点依次进行单十字内插;
3)对于任一待插点,寻找对应横坐标的垂直终止点,若垂直终止点的纵坐标小于待插点的纵坐标,则直到满足垂直终止点的纵坐标大于待插点为止,更新垂直起始点与垂直终止点,至此得到了该待插点的上、下、左、右四个点的坐标,分别对应起始点、终止点、区间头、区间尾,之后就是按照反距离权重插值的算法计算出预测值。
进一步的,步骤S3中,通过改变指针步进方向,把图片视作旋转45°区域,以水平垂直十字内插方法为基础,进行对角线方向的单十字内插。
进一步的,步骤S4中,使用上两步单十字内插或加速单十字内插的预测值作平均,即得到待插点的双十字内插预测值。
步骤S2和S3中,根据待插点密度选择单十字内插或加速单十字内插。
本发明在可接受的插值误差下极大提高了计算机处理的实时性,在低性能终端上仅用单核CPU实现了高速内插,还可以在脚本语言上实现实时插值。
传统的插值算法过于依赖GPU等设备进行并行运算,移植到低性能终端上会导致性能无法接受。而本发明仅使用单核CPU即可进行插值运算。对于密集成块待插点区域,加速单十字内插算法在搜索近邻点的时间复杂度上比扩大包围盒搜索算法快两个量级。并且,十字内插算法可以稳定计算不同大小的栅格数据插值,双十字内插法则在精度评估上略优于4近邻IDW方法。
附图说明
图1是所述单十字内插算法的示意图;
图2是所述加速单十字内插的示意图;
图3是更新垂直起始点与终止点的示意图;
图4是栅格图按照顺时针旋转的示意图;
图5是采用扩大包围盒搜索时的复杂度示意图;
图6是不同栅格大小与待插点密度下的性能对比示意图;
图7是本发明的流程示意图。
具体实施方式
本发明在IDW方法基础上进行改进,利用栅格数据上待插点的十字方向近邻来插值,简称为单十字内插算法。根据栅格数据遍历的特点,仅利用栅格数据上待插点的水平垂直十字方向与对角线十字方向近邻点共计八个近邻点来快速完成两趟十字内插算法,简称为双十字内插算法。并提出一种加速单十字内插算法,仅遍历两次栅格数据即可完成一趟单十字内插算法。
如图7所示,本发明的主要实施方案可以概括为四个步骤:
S1:加载待插值图像与插值区域(掩膜)到内存:把待插值图像与插值区域(掩膜)加载到内存,插值区域不能与图像边界直接接触。
十字内插法需要读取内存中的图像,因此需要预先把图像从其他平台、设备或者驱动器加载到内存。若图像规定了缺失像素点的缺省值,则可以通过判断像素点是否为缺省值,从而在内存中直接生成掩膜。由于内插的特性,不能给边界像素进行插值,需要注意插值区域不能触碰图像边界。
S2:进行沿着水平垂直方向的单十字内插或加速单十字内插:在寻找栅格数据邻近点时,一种有代表性的IDW方法,是扩大包围盒搜索算法。设置以待插点为中心的包围盒半径为r,沿着盒边缘遍历数据点,直到找到n个测量值为止,逐渐扩大半径。该算法对每一个测量值计算欧式距离,需要进行开平方操作。
考虑到传统算法的性能瓶颈在于寻找邻近点时对已经遍历过的点多次遍历的冗余过程,以及求欧式距离时的指数运算,简化版的IDW算法仅考虑十字方向的测量值,加速了遍历过程,而且待插点十字方向的点的欧式距离天然地不需要开平方操作,可以简单地通过偏移量做乘法运算获得距离值,从而减少算法在求解权重时运算的耗时。
如图1所示的单十字内插算法,白底格点是测量值;灰底格点代表待插点;以黑色实心圆形符号标记的格点代表当前待插点;待插点由用户指定,以掩膜的方式作为参数传入。按像素遍历掩膜,考虑遍历指针所处位置是待插点的情况:在待插图像中,当前待插点指向的四个以黑色实心三角形符号标记的格点代表着插值所参考的上、下、左、右(或斜上下)四个测量值。十字内插算法在搜索遍历点时,偏移量仅单独作用在某一个坐标值上,相比于扩大包围盒搜索算法,简化了程序逻辑的复杂程度。十字方向的点自然地平分四个象限,采样相对均匀,而反距离权重本身解决了在均值内插时十字长短差距大时导致的插值误差增大的问题。
针对密集成块的待插值区域,现在介绍一种在密集区域加速单十字内插算法的思路。考虑一般十字内插算法,通过待插点的上、下、左、右(或斜上下)四个测量值可以直接计算出预测值。而且与待插点在同一待插连通区域,并处于待插点的十字方向的一系列待插点都共享一半的相同测量值。可以得到加速优化思路,即在遍历时候缓存在十字方向的那些待插点与测量值,充分利用遍历得到的先验知识,减少遍历过程。
步骤S2和S3中,根据待插点密度选择单十字内插或加速单十字内插,当密度小于50%时,选择单十字内插,当密度大于50%时,选择加速单十字内插。
若栅格数据的数据量大小为N*M(N≤M),设M为图像宽度,加速单十字内插算法需要开辟长度为2M的整型缓存数组来存放垂直方向上的两个测量值,前M个整型称为垂向起始点,用来缓存对应栅格横坐标的上方向的测量值纵坐标,相似的,后M个整型称为垂直终止点,用来缓存下方向的测量值纵坐标。
算法开始时,先初始化缓存数组并置为零,按行依次从左到右进行遍历。当从第一个测量值遍历到待插点时,记录第一个测量值的坐标为区间头,继续向后遍历直到指针从待插点移动到第二个测量值为止,第二个测量值的坐标标记为区间尾。区间头与区间尾内仅包含待插点,对这个待插区间内每个待插点依次进行十字内插。对于待插区间内任一待插点,寻找对应横坐标的垂直终止点,若终止点的纵坐标小于待插点的纵坐标,则直到满足终止点的纵坐标大于待插点为止,更新垂直起始点与终止点,至此得到了待插点的上(起始点)、下(终止点)、左(区间头)、右(区间尾)四个点的坐标,之后就是按照反距离权重插值的算法计算出预测值。
内插点的八邻域可以都是内插点,与本方法不冲突。对于一个给定的待插点,本方法取的是该点的水平方向与垂直方向中非待插的最近测量点作为插值依据,与该点八邻域点无关。在本段所述待插区间中,区间的待插点实际上共享水平方向上的区间头和尾这两个测量值,也是以此作为加速原理。
加速单十字内插算法如下图2所示,以黑框空心三角形符号标记的格点代表着缓存数组,每一列都会缓存两个坐标,给其跨越的待插点提供上、下测量值参考。
更新垂直起始点与终止点的过程即单列向下遍历,以终止点的坐标为遍历的开始位置,向纵向寻找直到找到从测量值到待插点的跳变过程时,更新起始点,相似的寻找到从待插点到测量值的跳变时更新终止点。更新过程如下图3所示,图中虚线代表着一次更新后垂直起始点与终止点的位置变化状态。
S3:沿着对角线方向的单十字内插或加速单十字内插:若把栅格图按照顺时针或逆时针旋转45°,栅格数据依然满足进行十字内插算法的条件。在实现时只需要调整指针步进的偏移,改变行列迭代的方式为对角方向,而不需要在内存里旋转栅格图像。图像顺时针旋转45°后的遍历过程,相当于图4所示的对角遍历过程。
通过改变指针步进方向,把图片视作旋转45°区域,以水平垂直十字内插方法为基础,进行对角线方向的单十字内插。
S4:双十字内插:使用上两步单十字内插的预测值作平均,即得到待插点的双十字内插预测值。在同一图像的不同处理中,待插区域是相同的,不同处理也仅是对每个待插点计算一个预测值。作平均就是把同一像素点的每个预测值取平均。
具体来说,在上两步中已经把未旋转的栅格图像与旋转的栅格图像分别进行十字内插算法,将得到的预测值逐点对应算平均,这样就考虑到了空间位置8个近邻的点。双十字内插算法则总会选取上、下、左、右、斜上、斜下方共8个点来计算插值,具有空间适应性。
S5:算法时间复杂度分析:
时间复杂度能够衡量一个算法的性能,实现算法的实时性则要求具备较低的时间复杂度。在对实时性要求比较高的场合,需要分析最坏情况下的时间复杂度。在待插点数量一定的情况下,单十字内插算法与使用4近邻的IDW的耗时区别仅在于搜索近邻点上,因此最坏的情况即是除了栅格数据最外层一圈格点之外,其余区域都是待插点。
若栅格数据的数据量大小为N*M(N≤M),设M为图像宽度。对于使用4近邻的IDW算法,扩大包围盒搜索会依照半径遍历待插点周围一圈的邻近点,直到搜索出4个测量值为止,半径以步长1递增。扩大包围盒搜索对于每个插值点的遍历近邻数如图5所示,可以由外向内逐层划分为c个圈,第i圈所包含的待插点个数为ci,具体计算如式1和式2所示。
ci=2*(N-2i+M-2i)-4 (2)
对于第i圈的每个待插点,搜索近邻直到触及栅格最外层的格点为止,停止搜索时的半径即为当前圈数,则可以得出搜索格点数量Ti与当前圈数i的关系式为:
Ti=(2i+1)2-1 (3)
则扩大包围盒搜索一幅N*M(N≤M)大小的栅格数据的时间复杂度为:
O(SIDW)=O(MN3) (4)
下方公式是来自上方公式的时间复杂度表示。大写O通常表示渐进上界,即算法时间复杂度。时间复杂度反应算法最坏情况所用时间。下面公式都是同样表示。
当使用十字内插算法时,每一个待插点搜索的格点数量都为N+M,可得十字内插算法的时间复杂度为:
Scross=(M-2)(N-2)*(N+M)
O(Scross)=O(M2N) (5)
而使用加速单十字内插算法时,由于缓存了测量点,对任意的栅格数据仅需要遍历两次即可完成搜索:
Scross_fast=2MN
O(Scross_fast)=O(MN) (6)
在一般插值任务中,通常栅格数据的长宽N与M的大小不会相差很大,可以把N近似看成M,比较三个时间复杂度可得:
O(MN)≈O(M2)<O(M2N)≈O(M3)<O(MN3)≈O(M4)此时加速单十字内插算法比未加速的版本快一个量级,比传统IDW算法快了两个量级,加速效果非常显著。
S6:性能实验:为对比三种算法的性能情况,缩放试验区域至长宽均为104、5×103、103、5×102,随机标记其内部100%、90%、50%、10%、1%的点为待插点,运行三种插值算法。
图6展示了不同栅格大小与待插点密度下的性能对比。图中纵坐标为对数坐标轴,利于在量级差距较大的时候进行对比分析。由图可见两种十字内插算法在相同场景下效率都比扩大包围盒搜索算法高。加速单十字内插算法在栅格数据大小固定时效率也相差不大,在密集区域的加速效率远超过其他两种算法,在填补大面积空洞区域时其余两种算法性能近乎于不可接受。在待插点密度小于50%的稀疏情况下,未加速的十字内插算法的性能相对突出。

Claims (5)

1.一种栅格数据快速插值方法,包括以下步骤:
S1:把待插值图像与插值区域加载到内存;
S2:在插值区域内,进行沿着水平垂直方向的单十字内插或加速单十字内插;
S3:在插值区域内,进行沿着对角线方向的单十字内插或加速单十字内插;
S4:把步骤S2-S3的预测值逐点求平均,进行双十字内插,得到待插点预测值;
其中,步骤S2和S3中,根据待插点密度选择单十字内插或加速单十字内插;当插值密度小于50%时,选择单十字内插,在插值区域内,进行沿着水平垂直方向的单十字内插;然后进行沿着对角线方向的单十字内插;当插值密度大于50%时,选择加速单十字内插,在插值区域内,进行沿着水平垂直方向的加速单十字内插;然后进行沿着对角线方向的加速单十字内插;
所述单十字内插,是指在待插图像中,当前待插点指向的四个格点代表着插值所参考的上、下、左、右四个测量值或斜上下四个测量值,四个测量值能够直接计算出预测值;
针对密集成块的待插值区域,所述加速单十字内插是通过当前待插点的上、下、左、右四个测量值直接计算出预测值,所有待插点中与当前待插点在同一待插连通区域,并处于当前待插点的十字方向的一系列待插点都共享一半的相同测量值。
2.根据权利要求1所述的栅格数据快速插值方法,其特征在于:步骤S1中,插值区域不能与图像边界直接接触。
3.根据权利要求1所述的栅格数据快速插值方法,其特征在于:针对密集成块的待插值区域进行单十字内插或加速单十字内插的过程如下所述:
1)若栅格数据的数据量大小为N*M,且N≤M,设M为图像宽度,加速单十字内插需要开辟长度为2M的整型缓存数组来存放垂直方向上的两个测量值,前M个整型称为垂向起始点,用来缓存对应栅格横坐标的上方向的测量值纵坐标,后M个整型称为垂直终止点,用来缓存对应栅格横坐标的下方向的测量值纵坐标;
2)先初始化缓存数组并置为零,按行依次从左到右进行遍历:当从第一个测量值遍历到待插点时,记录第一个测量值的坐标为区间头,继续向后遍历直到指针从待插点移动到第二个测量值为止,记录第二个测量值的坐标为区间尾,对这个区间内每个点依次进行单十字内插;
3)对于任一待插点,寻找对应横坐标的垂直终止点,若垂直终止点的纵坐标小于待插点的纵坐标,则直到满足垂直终止点的纵坐标大于待插点为止,更新垂直起始点与垂直终止点,至此得到了该待插点的上、下、左、右四个点的坐标,分别对应起始点、终止点、区间头、区间尾,之后就是按照反距离权重插值的算法计算出预测值。
4.根据权利要求1所述的栅格数据快速插值方法,其特征在于:步骤S3中,通过改变指针步进方向,把图片视作旋转45°区域,以水平垂直十字内插方法为基础,进行对角线方向的单十字内插。
5.根据权利要求1所述的栅格数据快速插值方法,其特征在于:步骤S4中,使用上两步单十字内插或加速单十字内插的预测值作平均,即得到待插点的双十字内插预测值。
CN202110383204.8A 2021-04-09 2021-04-09 一种栅格数据快速插值方法 Active CN113077380B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110383204.8A CN113077380B (zh) 2021-04-09 2021-04-09 一种栅格数据快速插值方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110383204.8A CN113077380B (zh) 2021-04-09 2021-04-09 一种栅格数据快速插值方法

Publications (2)

Publication Number Publication Date
CN113077380A CN113077380A (zh) 2021-07-06
CN113077380B true CN113077380B (zh) 2023-12-08

Family

ID=76615826

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110383204.8A Active CN113077380B (zh) 2021-04-09 2021-04-09 一种栅格数据快速插值方法

Country Status (1)

Country Link
CN (1) CN113077380B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102096902A (zh) * 2010-11-30 2011-06-15 上海高清数字科技产业有限公司 考虑十字处理的图像处理方法、数字图像缩放方法及装置
CN102917198A (zh) * 2012-09-04 2013-02-06 上海高清数字科技产业有限公司 一种去隔行技术中的内插方法及系统
CN106570936A (zh) * 2016-11-14 2017-04-19 河海大学 一种基于栅格dem数据的等距离权重内插加密方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109993693B (zh) * 2017-12-29 2023-04-25 澜至电子科技(成都)有限公司 用于对图像进行插值的方法和装置

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102096902A (zh) * 2010-11-30 2011-06-15 上海高清数字科技产业有限公司 考虑十字处理的图像处理方法、数字图像缩放方法及装置
CN102917198A (zh) * 2012-09-04 2013-02-06 上海高清数字科技产业有限公司 一种去隔行技术中的内插方法及系统
CN106570936A (zh) * 2016-11-14 2017-04-19 河海大学 一种基于栅格dem数据的等距离权重内插加密方法

Also Published As

Publication number Publication date
CN113077380A (zh) 2021-07-06

Similar Documents

Publication Publication Date Title
US10587806B2 (en) Method and apparatus for motion estimation
US8811749B1 (en) Determining correspondence between image regions
US8861869B2 (en) Determining correspondence between image regions
CN103167218B (zh) 一种基于非局部性的超分辨率重建方法和设备
WO2022141178A1 (zh) 图像处理方法及装置
CN110992246B (zh) 影像的金字塔分层切片方法
JP4982471B2 (ja) 信号処理方法、装置およびプログラム
CN111091567B (zh) 医学图像配准方法、医疗设备及存储介质
JP5166156B2 (ja) 解像度変換装置、方法およびプログラム
Zhang et al. A novel tracking method based on improved FAST corner detection and pyramid LK optical flow
US8300015B2 (en) Method of detecting the movement of an entity equipped with an image sensor and device for implementing same
CN115578255A (zh) 一种基于帧间亚像素块匹配的超分辨率重建方法
CN113077380B (zh) 一种栅格数据快速插值方法
CN115546027A (zh) 图像缝合线确定方法、装置以及存储介质
RU2747965C1 (ru) Обработка окклюзий для frc c помощью глубокого обучения
JP2020098455A (ja) 物体識別システム、物体識別方法、並びに、画像識別プログラム
CN107995488A (zh) 一种基于直线扩展收缩模式的图像缩放方法
CN115841496A (zh) Sem图像的轮廓提取方法、装置、计算机设备及存储介质
CN113012165A (zh) 一种基于gpu的各向异性视频过分割方法
CN115115520A (zh) 一种具有提高图像清晰度和分辨率的图像成像方法
CN113497886B (zh) 视频处理方法、终端设备及计算机可读存储介质
US8483505B2 (en) Rendering piece-wise smooth image from colour values along paths
JP2009064193A (ja) 解像度変換装置及び方法
CN107274430A (zh) 一种对象运动轨迹预测方法和装置
CN108876889B (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