CN114612513B - 一种基于fpga的图像金字塔光流值计算方法及系统 - Google Patents
一种基于fpga的图像金字塔光流值计算方法及系统 Download PDFInfo
- Publication number
- CN114612513B CN114612513B CN202210240902.7A CN202210240902A CN114612513B CN 114612513 B CN114612513 B CN 114612513B CN 202210240902 A CN202210240902 A CN 202210240902A CN 114612513 B CN114612513 B CN 114612513B
- Authority
- CN
- China
- Prior art keywords
- optical flow
- layer
- frame
- image pyramid
- feature 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
Links
- 238000004364 calculation method Methods 0.000 title claims abstract description 49
- 230000003287 optical effect Effects 0.000 claims abstract description 147
- 238000012795 verification Methods 0.000 claims abstract description 25
- 230000002457 bidirectional effect Effects 0.000 claims abstract description 24
- 238000000034 method Methods 0.000 claims description 66
- 230000033001 locomotion Effects 0.000 claims description 31
- 230000008569 process Effects 0.000 claims description 19
- 238000004422 calculation algorithm Methods 0.000 claims description 18
- 238000013461 design Methods 0.000 claims description 14
- 238000012216 screening Methods 0.000 claims description 13
- 238000000605 extraction Methods 0.000 claims description 10
- 239000011159 matrix material Substances 0.000 claims description 10
- 238000010586 diagram Methods 0.000 description 13
- 238000003860 storage Methods 0.000 description 12
- 238000012545 processing Methods 0.000 description 10
- 230000005540 biological transmission Effects 0.000 description 8
- 238000004590 computer program Methods 0.000 description 7
- 230000001133 acceleration Effects 0.000 description 6
- 230000000694 effects Effects 0.000 description 6
- 230000008901 benefit Effects 0.000 description 5
- 230000006870 function Effects 0.000 description 5
- 230000004044 response Effects 0.000 description 4
- 230000008859 change Effects 0.000 description 3
- 238000001514 detection method Methods 0.000 description 3
- 238000006073 displacement reaction Methods 0.000 description 3
- 238000005457 optimization Methods 0.000 description 3
- 238000004519 manufacturing process Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000001934 delay Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/20—Analysis of motion
- G06T7/246—Analysis of motion using feature-based methods, e.g. the tracking of corners or segments
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/04—Architecture, e.g. interconnection topology
- G06N3/045—Combinations of networks
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/20—Analysis of motion
- G06T7/269—Analysis of motion using gradient-based methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20016—Hierarchical, coarse-to-fine, multiscale or multiresolution image processing; Pyramid transform
-
- Y—GENERAL 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
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
- Y02D10/00—Energy efficient computing, e.g. low power processors, power management or thermal management
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Multimedia (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Biophysics (AREA)
- Evolutionary Computation (AREA)
- Artificial Intelligence (AREA)
- Biomedical Technology (AREA)
- Health & Medical Sciences (AREA)
- Computational Linguistics (AREA)
- Data Mining & Analysis (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Molecular Biology (AREA)
- Computing Systems (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Physics (AREA)
- Software Systems (AREA)
- Image Analysis (AREA)
Abstract
本发明公开了一种基于FPGA的图像金字塔光流值计算方法及系统,将第K帧和第K+1帧图像作为输入,经异步FIFO实现位宽和时钟的切换后写入DDR中;以第K帧特征点的坐标作为输入,从DDR中读取以第K帧特征点为中心的连续两帧图像所需最大像素窗口,隐式的创建图像金字塔;对读取的最大像素窗口进行双线性插值,筛选满足特征验证条件的特征点;对筛选后的特征点进行光流的双向估计与迭代求解,经特征点光流的双向估计与迭代求解,得到图像金字塔第1层中第K帧特征点在第K+1帧中的坐标,将图像金字塔第1层中第K+1帧特征点的坐标作为第K帧特征点经光流计算最终得到的结果。本发明使用金字塔的逐层估计,并在每层金字塔中使用了双向估计与迭代求解,大大提升了光流估计的精度。
Description
技术领域
本发明属于FPGA平台算法硬件加速技术领域,具体涉及一种基于FPGA的图像金字塔光流值计算方法及系统。
背景技术
近年来,随着各种硬件平台的迭代更新,而原有算法不能满足实时性要求的前提下,算法硬件化加速获得了人们的研究兴趣。光流法是计算机视觉中一个基础算法,可广泛用于运动检测、运动估计及视频分析等领域,由于匹配计算的高度复杂性,对于实际应用环境来说通常是一项非常具有挑战性的任务,难以实现处理的高速度和高精度。同时,为了避免两帧物体运动位移较大情况下算法的误差影响,建立图像金字塔对光流来说不可或缺,这需要大量的计算能力和存储区域。因此,为了实现高分辨率的实时处理,设计特定应用的光流系统变得至关重要。
光流法最大的问题是计算复杂、速度慢,限制了它在实际系统尤其是嵌入式系统中的应用。基于FPGA计算光流可以大大缩短开发周期,在并行数据处理和流水线体系结构的设计上具有较大优势,能够使光流定位达到又快又准的效果。现有的基于FPGA加速的光流法硬件设计主要分为两类,第一类是在原有光流法基础上省去了图像金字塔的建立,另一类是在原有光流法基础上省去了每层光流计算中的双向估计与迭代过程。
现有一种划分子图像的方法用来减小所需的电路尺寸,但忽略了实现图像金字塔的重要性。光流法的假设中包括小运动这一前提,在两帧间物体运动快速时,省去金字塔容易导致大量特征点跟踪失败,对物体旋转及较大几何变换敏感,误匹配率高。
现有一种改进的基于KLT特征跟踪器的无人机返回框架,并集成了多种优化方法设计了硬件加速架构。虽然实现了图像金字塔的创建,但仅单纯依靠金字塔结构将图像的大运动分解成符合光流计算假设条件的小运动,忽略光流计算的双向估计和迭代过程,精度难免造成损失。
显然,这两类设计虽然在时间运行上得到一定程度的改善,但在精度上会逊色于OpenCV的光流法。
综上所述,尽管基于FPGA加速的光流法硬件设计已经存在多种解决方案,但是仍缺少一种兼顾光流金字塔实时性及保证特征点跟踪的精度的方法。
发明内容
本发明所要解决的技术问题在于针对上述现有技术中的不足,提供一种基于FPGA的图像金字塔光流值计算方法及系统,基于FPGA平台,使用预读取特征点像素窗口的方法减少图像金字塔建立的时间和存储空间,改进光流模型,保留每层计算中的双向估计与迭代过程,优化方法的固有并行性来加速计算,兼顾了运行效率与精度。
本发明采用以下技术方案:
一种基于FPGA的图像金字塔光流值计算方法,包括以下步骤:
S1、将第K帧和第K+1帧图像作为输入,经异步FIFO实现位宽和时钟的切换后写入DDR中;
S2、以第K帧特征点的坐标作为输入,从步骤S1写入存储第K帧和第K+1帧图像的DDR中读取以第K帧特征点为中心的连续两帧图像所需最大像素窗口,隐式的创建图像金字塔;
S3、对步骤S2读取的最大像素窗口进行双线性插值,将双线性插值过后的窗口进行梯度求解,并对相应的特征点进行KLT特征验证,筛选满足特征验证条件的特征点;
S4、对步骤S3筛选后的特征点进行光流的双向估计与迭代求解,得到图像金字塔第3层中第K帧特征点在第K+1帧中的坐标,以坐标初始化图像金字塔第2层中第K+1帧特征点的坐标,继续经步骤S3对特征点光流的双向估计与迭代求解,得到图像金字塔第2层中第K帧特征点在第K+1帧中的坐标,以坐标初始化图像金字塔第1层中第K+1帧特征点的坐标,经步骤S3和步骤S4中的特征点光流的双向估计与迭代求解,得到图像金字塔第1层中第K帧特征点在第K+1帧中的坐标,将图像金字塔第1层中第K+1帧特征点的坐标作为第K帧特征点经光流计算最终得到的结果。
具体的,步骤S1中,以32bit为一个单元存储,将4个8bit的像素数据写入DDR中。
具体的,步骤S2中,以特征点Pk、Pk+1为中心、像素窗口左上角坐标为基地址,分别读取第K帧、第K+1帧图像在各层中所需的像素窗口,隐式的创建图像金字塔,通过采用预读取像素窗口的方式实现图像金字塔的多层读取。
具体的,步骤S3中,最大像素窗口的尺寸Lk为:
Lk=max(Li,Ld)
其中,Li,Ld分别为8×8和10×10的尺寸。
具体的,步骤S3中,对相应的特征点进行KLT特征验证具体为:
对梯度窗口加和得到像素点的特征提取矩阵G,根据特征提取矩阵G的特征值λ与给定阈值λt进行比较,如果λ>λt,特征点满足要求,设置输出状态属性status为1,否则设置为0。
进一步的,特征提取矩阵G及特征值判定方程f(λ)具体为:
其中,a,b,c分别为图像在以特征点为中心的7×7领域内的像素窗口中水平方向梯度的平方和,图像在以特征点为中心的7×7领域内的像素窗口中水平方向梯度与垂直方向梯度乘积的加和,图像在以特征点为中心的7×7领域内的像素窗口中垂直方向梯度的平方和,Ix为图像在特征点水平方向的梯度,Iy为图像在特征点垂直方向的梯度。
具体的,步骤S4中,基于梯度信息构造运动模型确定光流估计的目标,利用光流法估计图像运动中的光流值,借助图像金字塔结构将大运动分解为若干连续小运动的同时,使用了双向估计与迭代过程,每层图像金字塔只读取一次在偏移量范围内的最大像素窗口,当层迭代求解完成后,将此时的特征点Pk+1更新值作为对应层光流估计所得的最终结果,当上一层金字塔的光流估计完成后,下层光流估计时将上层估计出的光流值作为初始光流值初始化本层图像,再在本层图像中根据光流算法进行光流估计,当图像金字塔自顶向下的光流全部计算完成后,最底层的特征点Pk+1更新值即为光流估计所得的最终结果,并重新计算此时的时间梯度窗口Wt表示为输出误差属性error。
进一步的,迭代过程中,当某一层图像金字塔迭代次数达到预先界定最大值;或者偏移量超过约束范围;或者计算的光流值变化量小于设定阈值时,跳出迭代,若对应层为最底层图像金字塔,直接输出光流估计结果,否则进入下一层图像金字塔继续进行计算。
更进一步的,第L层的初始光流值PL和图像总光流d表示为:
PL-1=2(PL+dL)
d=P1+d1
跟踪成功的特征点满足条件如下:
其中,PL-1为第L-1层的初始光流值,dL为第L层的光流计算结果,P1为第1层的初始光流值,d1为第1层的光流计算结果,status为输出状态属性。
第二方面,本发明实施例提供了一种基于FPGA的图像金字塔光流值计算系统,包括:
输入模块,将第K帧和第K+1帧图像作为输入,经异步FIFO实现位宽和时钟的切换后写入DDR中;
读取模块,以第K帧特征点的坐标作为输入,从输入模块写入存储第K帧和第K+1帧图像的DDR中读取以第K帧特征点为中心的连续两帧图像所需最大像素窗口,隐式的创建图像金字塔;
筛选模块,对读取模块读取的最大像素窗口进行双线性插值,将双线性插值过后的窗口进行梯度求解,并对相应的特征点进行KLT特征验证,筛选满足特征验证条件的特征点;
设计模块,对筛选模块筛选后的特征点进行光流的双向估计与迭代求解,得到图像金字塔第3层中第K帧特征点在第K+1帧中的坐标,以坐标初始化图像金字塔第2层中第K+1帧特征点的坐标,经筛选模块特征点光流的双向估计与迭代求解,得到图像金字塔第2层中第K帧特征点在第K+1帧中的坐标,以坐标初始化图像金字塔第1层中第K+1帧特征点的坐标,再经设计模块的特征点光流的双向估计与迭代求解,得到图像金字塔第1层中第K帧特征点在第K+1帧中的坐标,将图像金字塔第1层中第K+1帧特征点的坐标作为第K帧特征点经光流计算最终得到的结果。
与现有技术相比,本发明至少具有以下有益效果:
本发明一种基于FPGA的图像金字塔光流值计算方法,使用DDR进行突发传输,做到连续的突发传输,随机读取图像中的像素窗口而无需进行整张图像的降采样来创建光流金字塔,减少读写访问次数和存储空间,且对于每一个特征点的光流计算,只需进行两次预读取(第K、K+1帧)即可完成;使用预读取特征点像素窗口的方法减少图像金字塔建立的时间和存储空间,改进光流模型,保留每层计算中的双向估计与迭代过程,利用优化方法的固有并行性加速计算,兼顾运行效率与精度,利用图像序列中像素在时间域上的变化以及相邻帧之间的相关性找到相邻帧之间存在的对应关系,从而计算出相邻帧之间物体的运动信息,表征世界中物体的真实运动。
进一步的,因硬件平台是基于Xilinx Zynq 7z100,其ARM架构的CPU cortex-A9是32bit的,所以在一个周期内,以32bit为一个单元存储4个8bit的像素数据比存储1个8bit的像素数据传输速率要快,提高了一定的效率。
进一步的,采用预读取像素窗口的方式节省了创建图像金字塔时降采样的时间开销,并且节省了对于图像金字塔中不同尺度图像存储的空间开销,大大提升了效率。
进一步的,在每层图像金字塔的光流计算中,因双线性插值、梯度解算和保证在偏移量范围内的移动都需要不同尺寸的像素窗口,为了避免重复读取,时间资源的消耗,提前计算最大的像素窗口尺寸,只需读取一次最大的像素窗口即可完成光流的计算。
进一步的,对特征点进行KLT特征验证一方面选择两个特征值不太小的点,可以排除噪声的部分影响,另一方面选择特征值差别不太大的点,认为这是角点,能够更好的表征图像的性质。
进一步的,以往KLT特征点的验证方式是通过开根号进行运算,由于在FPGA中,开根号的逻辑资源开销大,且计算缓慢,所以通过特征提取矩阵G及特征值判定方程f(λ)更改验证条件将原本复杂的开根号运算变为简单的加减法运算及乘法运算。
进一步的,光流法的假设条件针对的是小运动,如果运动速度较快时,算法的误差会很大,依靠图像金字塔的结构将大运动分解为若干小运动,从而在不同尺度的各层图像金字塔中求解光流进行叠加来减小误差。
进一步的,将光流法构建成一个非线性优化问题,采用逐次搜索的迭代算法,找到新的近似点,如此迭代反复,直到满足预定的精度要求或达到迭代次数,从而提升光流值的精度。
进一步的,通过分层计算光流值进行叠加的方式,从第L层(最顶层)的初始光流值PL逐步计算得到图像总光流d,可以将大运动分解为小运动进行光流计算,减小误差。设置跟踪成功的特征点必须满足KLT的验证条件且跟踪后的点与原始点的像素差小于一定阈值,可以提升光流运算结果的精度。
综上所述,本发明方法平均处理一帧的时间为10.82ms,即平均每秒可处理至少90帧图像,完全满足实时性的要求。
下面通过附图和实施例,对本发明的技术方案做进一步的详细描述。
附图说明
图1为本发明流程图;
图2为第K帧8×8像素窗口在各层降采样过程图;
图3为方向梯度求解过程图;
图4为第K+1帧各层预读取像素窗口过程图;
图5为本发明光流求解硬件实现流程图;
图6为本发明与对比算法的光流效果图;
图7为本发明的仿真时序图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
在本发明的描述中,需要理解的是,术语“包括”和“包含”指示所描述特征、整体、步骤、操作、元素和/或组件的存在,但并不排除一个或多个其它特征、整体、步骤、操作、元素、组件和/或其集合的存在或添加。
还应当理解,在本发明说明书中所使用的术语仅仅是出于描述特定实施例的目的而并不意在限制本发明。如在本发明说明书和所附权利要求书中所使用的那样,除非上下文清楚地指明其它情况,否则单数形式的“一”、“一个”及“该”意在包括复数形式。
还应当进一步理解,在本发明说明书和所附权利要求书中使用的术语“和/或”是指相关联列出的项中的一个或多个的任何组合以及所有可能组合,并且包括这些组合,例如,A和/或B,可以表示:单独存在A,同时存在A和B,单独存在B这三种情况。另外,本文中字符“/”,一般表示前后关联对象是一种“或”的关系。
应当理解,尽管在本发明实施例中可能采用术语第一、第二、第三等来描述预设范围等,但这些预设范围不应限于这些术语。这些术语仅用来将预设范围彼此区分开。例如,在不脱离本发明实施例范围的情况下,第一预设范围也可以被称为第二预设范围,类似地,第二预设范围也可以被称为第一预设范围。
取决于语境,如在此所使用的词语“如果”可以被解释成为“在……时”或“当……时”或“响应于确定”或“响应于检测”。类似地,取决于语境,短语“如果确定”或“如果检测(陈述的条件或事件)”可以被解释成为“当确定时”或“响应于确定”或“当检测(陈述的条件或事件)时”或“响应于检测(陈述的条件或事件)”。
在附图中示出了根据本发明公开实施例的各种结构示意图。这些图并非是按比例绘制的,其中为了清楚表达的目的,放大了某些细节,并且可能省略了某些细节。图中所示出的各种区域、层的形状及它们之间的相对大小、位置关系仅是示例性的,实际中可能由于制造公差或技术限制而有所偏差,并且本领域技术人员根据实际所需可以另外设计具有不同形状、大小、相对位置的区域/层。
本发明提供了一种基于FPGA的图像金字塔光流值计算方法,基于FPGA平台,使用预读取特征点像素窗口的方法减少了图像金字塔建立的时间和存储空间,改进光流模型,保留了每层计算中的双向估计与迭代过程,优化算法的固有并行性来加速计算,兼顾了运行效率与精度,本发明方法平均处理一帧的时间为10.82ms,平均每秒可处理至少90帧图像,完全满足实时性的要求。
请参阅图1,本发明一种基于FPGA的图像金字塔光流值计算方法,包括以下步骤:
S1、连续两帧图像(第K、K+1帧)作为输入经过异步FIFO实现位宽和时钟的切换,通过控制电路写入DDR(外部存储器)中;
将连续两帧图像(第K、K+1帧)数据流经过异步FIFO实现位宽和时钟的切换,并起到数据缓冲的作用写入DDR中,对图像进行整合,采用乒乓操作,工作频率为120MHz。因硬件平台是基于Xilinx Zynq 7z100,其ARM架构的CPU cortex-A9是32bit的,所以以32bit为一个单元存储4个8bit的像素数据,方便步骤S2进行像素窗口的读取,提高了一定的效率。
S2、以第K帧特征点坐标作为输入,从步骤S1的DDR中读取以特征点为中心的连续两帧所需最大像素窗口;
创建图像金字塔需要对原图进行降采样操作,直接对原图进行处理得到各层图像进行存储不仅效率低,还需要大量的逻辑资源,得益于步骤S1中DDR突发传输的优势,采用预读取像素窗口的方式实现图像金字塔的多层读取。
请参阅图3,对于第K帧图像,每层像素窗口需经过步骤S3的双线性插值和梯度求解得到方向梯度窗口Wx、Wy,因计算要求Wx、Wy尺寸大小为7×7,故每层读取的原像素窗口需以输入的特征点Pk为中心,尺寸大小为10×10。
请参阅图2,针对不同层的10×10像素窗口,使用固定位置取值进行降采样,图2以8×8的像素窗口为例,经过1次降采样可得到4×4的像素窗口,经过2次降采样可得到2×2的像素窗口。
故为了得到上述在3层中的10×10像素窗口,需要遍历以点Pk为中心,尺寸分别为40×40、20×20、10×10的像素窗口,即第K帧图像需要遍历的最大像素窗口是以点Pk为中心,尺寸为40×40的。
请参阅图4,对于第K+1帧图像,每层像素窗口需经过步骤S3的双线性插值,步骤S4的迭代求解,因计算要求其中的时间梯度窗口Wt尺寸大小为7×7,故每层读取的原像素窗口需以特征点Pk+1初始坐标(Pk+1初始坐标与Pk一致)为中心,尺寸大小分别为24×24、16×16、12×12,
参见上述第K帧降采样方式,故需要遍历以点Pk+1为中心,尺寸分别为24×24、32×32、48×48的像素窗口,即第K+1帧图像需要遍历的最大像素窗口是以点Pk+1为中心,尺寸为48×48的。
综上,以特征点Pk、Pk+1为中心、像素窗口左上角坐标为基地址,分别读取第K帧、第K+1帧图像在各层中所需的像素窗口,隐式的创建图像金字塔。
S3、对步骤S2中的像素窗口进行双线性插值,保证光流法的精度,将插值过后的窗口进行梯度求解,并对相应的特征点进行KLT特征验证,判断当前特征点是否为一个适合跟踪的特征点;
在得到步骤S2每层的特征点领域像素窗口后,需对第K、K+1帧每层的领域像素窗口进行双线性插值,对插值后的第K帧窗口进行梯度求解。因两帧间进行计算的方向梯度窗口Wx、Wy,时间梯度窗口Wt尺寸为7×7,以某一层为例说明,经过以下处理获得经过双线性插值、梯度求解所需要的7×7窗口。
双线性插值根据特征点的浮点数值,计算四个坐标二次线性插值的系数组成一个2×2卷积核IW对一个尺寸为Li=8×8的像素窗口Wi做卷积,得到更加精确的7×7领域像素窗口。在做卷积操作之前,由于在FPGA中,浮点数运算的逻辑资源开销大,且计算较慢,为了高效的计算效率,对所有浮点数使用定点数的形式进行计算,下式IW皆为定点数,IW表示为:
a=float(x)-int(x)
b=float(y)-int(y)
iw00=(1-a)*(1-b)
iw01=a*(1-b)
iw10=(1-a)*b
iw11=1-iw00-iw01-iw10
其中,Scharr算子具有比Sobel算子更加良好的效果,且抗噪能力很好。
请参阅图3,求解水平和垂直梯度时,先用卷积核IW对一个尺寸Ld=10×10的像素窗口Wd做卷积得到一个精确的9×9领域像素窗口,再使用3×3的水平和垂直Sobel算子对9×9领域像素窗口做卷积分别得到7×7的水平梯度窗口Wx和垂直梯度窗口Wy。
综上,第K帧所需要读取的最大像素窗口为Wk,Wk的尺寸为Lk,对于第K帧图像,所需要预读取的是以特征点Pk为中心、尺寸为Lk的像素窗口。
Lk表示为:
Lk=max(Li,Ld)=10×10
经过双线性插值和梯度解算过后,需要对输入的特征点Pk进行KLT验证,通过对梯度窗口加和得到像素点的特征提取矩阵G,根据矩阵G的特征值λ与给定阈值λt进行比较,判断该特征点是否满足要求,排除噪声的影响。特征提取矩阵G及特征值判定方程f(λ)表示为:
若KLT验证成功,则设置输出状态属性status为1,否则设置为0。以往KLT特征点的验证方式是通过开根号进行运算,由于在FPGA中,开根号的逻辑资源开销大,且计算缓慢,故更改验证条件表示为:
其中,a,b,c分别为图像在以特征点为中心的7×7领域内的像素窗口中水平方向梯度的平方和,图像在以特征点为中心的7×7领域内的像素窗口中水平方向梯度与垂直方向梯度乘积的加和,图像在以特征点为中心的7×7领域内的像素窗口中垂直方向梯度的平方和,Ix为图像在特征点水平方向的梯度,Iy为图像在特征点垂直方向的梯度。
S4、对步骤S3中经过筛选过后的特征点实现光流的双向估计与迭代求解,完成特征点的光流值和相关状态属性的计算。
对S3中满足KTL验证要求的特征点进行光流估计,光流估计需要提前满足三个假设条件:
亮度恒定;小运动;空间一致。
基于梯度信息构造运动模型,则光流估计的目标表示为:
其中,I、J分别为第K、K+1帧图像,d=[dxdy]T为光流的估计值,整理上式,可得:
∫∫W[J(x)-I(x)+gTd]g(x)w(x)dx=0
其中,
-[∫∫Wg(x)g(x)Tw(x)dx]d=∫∫W[J(x)-I(x)]g(x)w(x)dx
Gd=e
为了利用光流法能够精确的估计出图像运动中的光流值,在借助图像金字塔结构将大运动分解为若干连续小运动的同时,使用了双向估计与迭代过程,这得益于步骤S2预读取特征点像素窗口模块的方式使得FPGA资源的使用在可控范围内。在迭代过程中,特征点Pk+1是不断更新的,导致以特征点Pk+1为中心的时间梯度窗口Wt也需要不断变化,但Wt的尺寸Lt=7×7是一直保持不变的。由于第K+1帧预读取领域像素窗口是以Pk+1初始坐标(Pk+1的初始坐标与Pk一致)为中心的,为了避免频繁的DDR读取降低算法运行效率,约束对于每一个特征点Pk+1,只读取一次在偏移量范围内的最大像素窗口,即预读取的像素窗口一经读取就不再改变,为了满足时间梯度窗口Wt在预读取的像素窗口不会产生越界,第K+1帧各层预读取像素窗口过程如图4所示。
从图4可知,约束在第3层特征点Pk+1的最大偏移量是2pixel,在第2层特征点Pk+1的最大偏移量是4pixel,在第1层特征点Pk+1的最大偏移量是8pixel,这里的偏移量是指特征点Pk+1的最大变化量,在图像金字塔各层光流估计的迭代过程中,跳出循环有3种情况,一是迭代次数达到预先界定最大值;二是一旦偏移量超过约束范围;三是计算的光流值变化量小于某一阈值,就会跳出迭代,认为此时的Pk+1更新值为对应层光流估计所得的最终结果。
表1为各层偏移量及两帧所预读取的像素窗口大小汇总。
当上一层金字塔的光流估计完成后,下层光流估计时将上层估计出的光流值作为初始光流值(假设最顶层图像初始光流值为0)初始化本层图像,再在本层图像中根据光流算法进行光流估计。
第L层的初始光流值PL和图像总光流d表示为:
PL-1=2(PL+dL)
d=P1+d1
当3层金字塔光流计算完成后,重新计算此时的时间梯度窗口Wt表示为输出误差属性error,它是衡量此时Pk和Pk+1两个特征点的误差指标。跟踪成功的特征点需满足下式:
。
本发明再一个实施例中,提供一种基于FPGA的图像金字塔光流值计算系统,该系统能够用于实现上述基于FPGA的图像金字塔光流值计算方法,具体的,该基于FPGA的图像金字塔光流值计算系统包括输入模块、读取模块、筛选模块以及设计模块。
其中,输入模块,将第K帧和第K+1帧图像作为输入,经异步FIFO实现位宽和时钟的切换后写入DDR中;
读取模块,以第K帧特征点的坐标作为输入,从输入模块写入存储第K帧和第K+1帧图像的DDR中读取以第K帧特征点为中心的连续两帧图像所需最大像素窗口,隐式的创建图像金字塔;
筛选模块,对读取模块读取的最大像素窗口进行双线性插值,将双线性插值过后的窗口进行梯度求解,并对相应的特征点进行KLT特征验证,筛选满足特征验证条件的特征点;
设计模块,对筛选模块筛选后的特征点进行光流的双向估计与迭代求解,得到图像金字塔第3层中第K帧特征点在第K+1帧中的坐标,以坐标初始化图像金字塔第2层中第K+1帧特征点的坐标,经筛选模块特征点光流的双向估计与迭代求解,得到图像金字塔第2层中第K帧特征点在第K+1帧中的坐标,以坐标初始化图像金字塔第1层中第K+1帧特征点的坐标,再经设计模块的特征点光流的双向估计与迭代求解,得到图像金字塔第1层中第K帧特征点在第K+1帧中的坐标,将图像金字塔第1层中第K+1帧特征点的坐标作为第K帧特征点经光流计算最终得到的结果。
光流就是利用图像序列中像素在时间域上的变化以及相邻帧之间的相关性来找到相邻帧之间存在的对应关系,从而计算出相邻帧之间物体的运动信息,表征世界中物体的真实运动。
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明一部分实施例,而不是全部的实施例。通常在此处附图中的描述和所示的本发明实施例的组件可以通过各种不同的配置来布置和设计。因此,以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围,而是仅仅表示本发明的选定实施例。基于本发明中的实施例,本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
实施例
请参阅图5,基于Xilinx Zynq 7z100平台,工作频率为120MHz,将分辨率为376×240的K、K+1两帧连续图像缓存入DDR中,并将预先在K帧中通过Fast角点提取的特征点作为输入进行光流和坐标初始化,通过AXI总线预读取两帧所需像素窗口,经由双线性插值、方向和时间梯度求解、KLT特征验证、光流估计等模块,进行光流值及相关状态属性的求解,最终将光流值及相关状态属性通过block RAM输出至计算机端进行数据通信,以便对数据进行处理与分析。
请参阅图6,其中,图6(a)与图6(b)分别为FPGA端和OpenCV端图像光流匹配效果图。图6(c)与图6(d)分别为FPGA端与OpenCV端的x、y方向上的特征点坐标对比曲线,可知本发明的精度与OpenCV端图像光流匹配效果基本一致。图7为完成一次完整的全局迭代光流金字塔计算的仿真时序图,可知本发明处理一帧的时间为10.82ms,且在仿真波形图的输出结果与OpenCV端基本一致。
为了进一步评估精度,选取200组分辨率为376×240的连续两帧K、K+1图像,并预先对所有第K帧图像进行Fast角点提取170左右个特征点,将这200组特征点序列和图像对数据作为输入计算本发明的光流结果,同时将相同数据作为OpenCV函数calcOpticalFlowPyrLK()的输入计算光流结果,精度评估如表2所示。
表2本发明与对比算法的精度对比
从实验结果可知,本发明和OpenCV的光流算法所能跟踪到的特征点数量基本一致,对于跟踪成功的特征点,在200组数据的验证下,与OpenCV光流结果对比,本发明的光流结果在x方向特征点坐标的平均位移是0.15个pixel,y方向特征点坐标的平均位移是0.06个pixel,由此可知本发明的精度与OpenCV的基本一致。
表3本发明的硬件资源占比
表4本发明与对比算法的性能对比
光流算法硬件系统的FPGA资源消耗如表3所示,可知本发明算法的总体硬件资源占比开销在可接受范围之内。对比本发明与基于FPGA的先进光流设计性能,运行速率如表4所示。尽管FAOFE在速率上实现了更快的加速,但本发明的实现达到了90fps以上,这明显高于普遍接受的实时性能(25fps),且本发明实现了图像金字塔的创建以及在光流估计中的迭代求解,在保证精度的前提下速率也大大提升。
综上所述,本发明一种基于FPGA的图像金字塔光流值计算方法及系统,具有以下优势:
1、使用DDR能够进行突发传输,对同一行中相邻的存储单元连续进行数据传输,只要指定起始列地址与突发长度,就会依次的自动对后面相应数量的存储单元进行读/写操作而不再需要控制器连续的提供列地址,除了第一笔数据的传输需要若干个周期延迟外,其后每个数据只需一个周期即可获得,而只要控制好两段突发读取命令的间隔周期,即可做到连续的突发传输,随机读取图像中的像素窗口而无需进行整张图像的降采样来创建光流金字塔,减少了读写访问次数和存储空间,且对于每一个特征点的光流计算,只需进行两次预读取(第K、K+1帧)即可完成。
2、对像素窗口做了预处理,使用双线性插值的方法提高了像素值的精度,进而提高了方向梯度窗口Wx、Wy、时间梯度窗口Wt和最终光流估计的精度。
3、光流求解过程使用了金字塔的逐层估计,并在每层金字塔中使用了双向估计与迭代求解,大大提升了光流估计的精度。
4、Pipeline并行性,尽管在光流求解过程中数据依赖性存在于不同的阶段,不同阶段中的操作必须一个接一个的进行计算,但使用Pipeline无需等待前一阶段的所有结果可用后再开始下一阶段。可以使用流水线的方式进行同时工作。
本领域内的技术人员应明白,本申请的实施例可提供为方法、系统、或计算机程序产品。因此,本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且,本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器,使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中,使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品,该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上,使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理,从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上内容仅为说明本发明的技术思想,不能以此限定本发明的保护范围,凡是按照本发明提出的技术思想,在技术方案基础上所做的任何改动,均落入本发明权利要求书的保护范围之内。
Claims (10)
1.一种基于FPGA的图像金字塔光流值计算方法,其特征在于,包括以下步骤:
S1、将第K帧和第K+1帧图像作为输入,经异步FIFO实现位宽和时钟的切换后写入DDR中;
S2、以第K帧特征点的坐标作为输入,从步骤S1写入存储第K帧和第K+1帧图像的DDR中读取以第K帧特征点为中心的连续两帧图像所需最大像素窗口,隐式的创建图像金字塔;
S3、对步骤S2读取的最大像素窗口进行双线性插值,将双线性插值过后的窗口进行梯度求解,并对相应的特征点进行KLT特征验证,筛选满足特征验证条件的特征点;
S4、对步骤S3筛选后的特征点进行光流的双向估计与迭代求解,得到图像金字塔第3层中第K帧特征点在第K+1帧中的坐标,以坐标初始化图像金字塔第2层中第K+1帧特征点的坐标,继续经步骤S3对特征点光流的双向估计与迭代求解,得到图像金字塔第2层中第K帧特征点在第K+1帧中的坐标,以坐标初始化图像金字塔第1层中第K+1帧特征点的坐标,经步骤S3的特征点光流的双向估计与迭代求解,得到图像金字塔第1层中第K帧特征点在第K+1帧中的坐标,将图像金字塔第1层中第K+1帧特征点的坐标作为第K帧特征点经光流计算最终得到的结果。
2.根据权利要求1所述的基于FPGA的图像金字塔光流值计算方法,其特征在于,步骤S1中,以32bit为一个单元存储,将4个8bit的像素数据写入DDR中。
3.根据权利要求1所述的基于FPGA的图像金字塔光流值计算方法,其特征在于,步骤S2中,以特征点Pk、Pk+1为中心、像素窗口左上角坐标为基地址,分别读取第K帧、第K+1帧图像在各层中所需的像素窗口,隐式的创建图像金字塔,通过采用预读取像素窗口的方式实现图像金字塔的多层读取。
4.根据权利要求1所述的基于FPGA的图像金字塔光流值计算方法,其特征在于,步骤S3中,最大像素窗口的尺寸Lk为:
Lk=max(Li,Ld)
其中,Li,Ld分别为8×8和10×10的尺寸。
5.根据权利要求1所述的基于FPGA的图像金字塔光流值计算方法,其特征在于,步骤S3中,对相应的特征点进行KLT特征验证具体为:
对梯度窗口加和得到像素点的特征提取矩阵G,根据特征提取矩阵G的特征值λ与给定阈值λt进行比较,如果λ>λt,特征点满足要求,设置输出状态属性status为1,否则设置为0。
6.根据权利要求5所述的基于FPGA的图像金字塔光流值计算方法,其特征在于,特征提取矩阵G及特征值判定方程f(λ)具体为:
其中,a,b,c分别为图像在以特征点为中心的7×7领域内的像素窗口中水平方向梯度的平方和,图像在以特征点为中心的7×7领域内的像素窗口中水平方向梯度与垂直方向梯度乘积的加和,图像在以特征点为中心的7×7领域内的像素窗口中垂直方向梯度的平方和,Ix为图像在特征点水平方向的梯度,Iy为图像在特征点垂直方向的梯度。
7.根据权利要求1所述的基于FPGA的图像金字塔光流值计算方法,其特征在于,步骤S4中,基于梯度信息构造运动模型确定光流估计的目标,利用光流法估计图像运动中的光流值,借助图像金字塔结构将大运动分解为若干连续小运动的同时,使用了双向估计与迭代过程,每层图像金字塔只读取一次在偏移量范围内的最大像素窗口,当层迭代求解完成后,将此时的特征点Pk+1更新值作为对应层光流估计所得的最终结果,当上一层金字塔的光流估计完成后,下层光流估计时将上层估计出的光流值作为初始光流值初始化本层图像,再在本层图像中根据光流算法进行光流估计,当图像金字塔自顶向下的光流全部计算完成后,最底层的特征点Pk+1更新值即为光流估计所得的最终结果,并重新计算此时的时间梯度窗口Wt表示为输出误差属性error。
8.根据权利要求7所述的基于FPGA的图像金字塔光流值计算方法,其特征在于,迭代过程中,当某一层图像金字塔迭代次数达到预先界定最大值;或者偏移量超过约束范围;或者计算的光流值变化量小于设定阈值时,跳出迭代,若对应层为最底层图像金字塔,直接输出光流估计结果,否则进入下一层图像金字塔继续进行计算。
9.根据权利要求8所述的基于FPGA的图像金字塔光流值计算方法,其特征在于,第L层的初始光流值PL和图像总光流d表示为:
PL-1=2(PL+dL)
d=P1+d1
跟踪成功的特征点满足条件如下:
其中,PL-1为第L-1层的初始光流值,dL为第L层的光流计算结果,P1为第1层的初始光流值,d1为第1层的光流计算结果,status为输出状态属性。
10.一种基于FPGA的图像金字塔光流值计算系统,其特征在于,包括:
输入模块,将第K帧和第K+1帧图像作为输入,经异步FIFO实现位宽和时钟的切换后写入DDR中;
读取模块,以第K帧特征点的坐标作为输入,从输入模块写入存储第K帧和第K+1帧图像的DDR中读取以第K帧特征点为中心的连续两帧图像所需最大像素窗口,隐式的创建图像金字塔;
筛选模块,对读取模块读取的最大像素窗口进行双线性插值,将双线性插值过后的窗口进行梯度求解,并对相应的特征点进行KLT特征验证,筛选满足特征验证条件的特征点;
设计模块,对筛选模块筛选后的特征点进行光流的双向估计与迭代求解,得到图像金字塔第3层中第K帧特征点在第K+1帧中的坐标,以坐标初始化图像金字塔第2层中第K+1帧特征点的坐标,经筛选模块特征点光流的双向估计与迭代求解,得到图像金字塔第2层中第K帧特征点在第K+1帧中的坐标,以坐标初始化图像金字塔第1层中第K+1帧特征点的坐标,再经设计模块的特征点光流的双向估计与迭代求解,得到图像金字塔第1层中第K帧特征点在第K+1帧中的坐标,将图像金字塔第1层中第K+1帧特征点的坐标作为第K帧特征点经光流计算最终得到的结果。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210240902.7A CN114612513B (zh) | 2022-03-10 | 2022-03-10 | 一种基于fpga的图像金字塔光流值计算方法及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210240902.7A CN114612513B (zh) | 2022-03-10 | 2022-03-10 | 一种基于fpga的图像金字塔光流值计算方法及系统 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114612513A CN114612513A (zh) | 2022-06-10 |
CN114612513B true CN114612513B (zh) | 2024-03-01 |
Family
ID=81862232
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202210240902.7A Active CN114612513B (zh) | 2022-03-10 | 2022-03-10 | 一种基于fpga的图像金字塔光流值计算方法及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114612513B (zh) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117876206B (zh) * | 2024-01-17 | 2024-07-23 | 电子科技大学 | 基于可重构计算阵列的金字塔光流加速电路 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104978728A (zh) * | 2014-04-08 | 2015-10-14 | 南京理工大学 | 一种光流法的图像匹配系统 |
CN108776971A (zh) * | 2018-06-04 | 2018-11-09 | 南昌航空大学 | 一种基于分层最近邻域的变分光流确定方法及系统 |
CN114140502A (zh) * | 2021-12-10 | 2022-03-04 | 江南大学 | 一种基于嵌入式gpu的图像光流计算方法、装置以及设备 |
Family Cites Families (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP3785225B1 (en) * | 2018-04-24 | 2023-09-13 | Snap Inc. | Efficient parallel optical flow algorithm and gpu implementation |
CN109753940B (zh) * | 2019-01-11 | 2022-02-22 | 京东方科技集团股份有限公司 | 图像处理方法及装置 |
-
2022
- 2022-03-10 CN CN202210240902.7A patent/CN114612513B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104978728A (zh) * | 2014-04-08 | 2015-10-14 | 南京理工大学 | 一种光流法的图像匹配系统 |
CN108776971A (zh) * | 2018-06-04 | 2018-11-09 | 南昌航空大学 | 一种基于分层最近邻域的变分光流确定方法及系统 |
CN114140502A (zh) * | 2021-12-10 | 2022-03-04 | 江南大学 | 一种基于嵌入式gpu的图像光流计算方法、装置以及设备 |
Non-Patent Citations (2)
Title |
---|
基于局部光流约束的角点匹配算法;蒋晓瑜;姚军;宋小杉;汪熙;;光学技术;20100315(02);全文 * |
蒋晓瑜 ; 姚军 ; 宋小杉 ; 汪熙 ; .基于局部光流约束的角点匹配算法.光学技术.2010,(02),全文. * |
Also Published As
Publication number | Publication date |
---|---|
CN114612513A (zh) | 2022-06-10 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US10776688B2 (en) | Multi-frame video interpolation using optical flow | |
CN113286194B (zh) | 视频处理方法、装置、电子设备及可读存储介质 | |
US11475542B2 (en) | Neural network system with temporal feedback for adaptive sampling and denoising of rendered sequences | |
US11557022B2 (en) | Neural network system with temporal feedback for denoising of rendered sequences | |
JP4964159B2 (ja) | ビデオのフレームのシーケンスにおいてオブジェクトを追跡するコンピュータに実装される方法 | |
US20190180469A1 (en) | Systems and methods for dynamic facial analysis using a recurrent neural network | |
CN111914698B (zh) | 图像中人体的分割方法、分割系统、电子设备及存储介质 | |
CN110120065B (zh) | 一种基于分层卷积特征和尺度自适应核相关滤波的目标跟踪方法及系统 | |
US11494879B2 (en) | Convolutional blind-spot architectures and bayesian image restoration | |
CN108717571B (zh) | 一种用于人工智能的加速方法和装置 | |
CN102509071B (zh) | 光流计算系统和方法 | |
CN111008040A (zh) | 缓存装置及缓存方法、计算装置及计算方法 | |
CN114612513B (zh) | 一种基于fpga的图像金字塔光流值计算方法及系统 | |
CN111161306A (zh) | 一种基于运动注意力的视频目标分割方法 | |
Li et al. | High throughput hardware architecture for accurate semi-global matching | |
CN112419372B (zh) | 图像处理方法、装置、电子设备及存储介质 | |
CN111382759A (zh) | 一种像素级分类方法、装置、设备及存储介质 | |
Holzer et al. | Multilayer adaptive linear predictors for real-time tracking | |
CN115019148A (zh) | 一种目标检测方法 | |
CN108764182B (zh) | 一种优化的用于人工智能的加速方法和装置 | |
KR20210051707A (ko) | 드론 정지 비행을 위한 영상 기반 특징점 추적 장치 및 그 방법 | |
CN103413326A (zh) | Fast approximated SIFT算法中特征点检测方法及装置 | |
CN112907750A (zh) | 一种基于卷积神经网络的室内场景布局估计方法及系统 | |
US11861811B2 (en) | Neural network system with temporal feedback for denoising of rendered sequences | |
CN114359222B (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 |