CN108399629B - 一种人工复眼相机的图像阵列光流估计方法 - Google Patents

一种人工复眼相机的图像阵列光流估计方法 Download PDF

Info

Publication number
CN108399629B
CN108399629B CN201810108789.0A CN201810108789A CN108399629B CN 108399629 B CN108399629 B CN 108399629B CN 201810108789 A CN201810108789 A CN 201810108789A CN 108399629 B CN108399629 B CN 108399629B
Authority
CN
China
Prior art keywords
sub
eye
optical flow
image
flow estimation
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.)
Expired - Fee Related
Application number
CN201810108789.0A
Other languages
English (en)
Other versions
CN108399629A (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.)
Southwest Jiaotong University
Original Assignee
Southwest Jiaotong University
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 Southwest Jiaotong University filed Critical Southwest Jiaotong University
Priority to CN201810108789.0A priority Critical patent/CN108399629B/zh
Publication of CN108399629A publication Critical patent/CN108399629A/zh
Application granted granted Critical
Publication of CN108399629B publication Critical patent/CN108399629B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/20Analysis of motion
    • G06T7/269Analysis of motion using gradient-based methods

Landscapes

  • Engineering & Computer Science (AREA)
  • Multimedia (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开了一种人工复眼相机的图像阵列光流估计方法,根据所用人工复眼的结构特点,对阵列图像中包含运动目标的图像进行检测,随后构建光流估计能量函数,并利用变分法对能量函数进行优化,完成人工复眼图像的光流估计。本发明与现有技术相比,其显著优点在于:能有效地对人工复眼图像序列进行光流估计,并且由于子眼图像间空间相关的约束项,能够对相邻子眼对应区域进行约束,获得一致的光流估计结果,解决人工复眼相机中子眼视场角小导致的图像边缘光流估计不准确的问题。

Description

一种人工复眼相机的图像阵列光流估计方法
技术领域
本发明涉及计算机视觉技术领域,特别是一种人工复眼相机的图像阵列光流估计方法。
背景技术
视觉是大部分生物认识世界、感知世界的一个重要手段,昆虫复眼作为一种广泛存在的视觉系统,具有体积小、近距离清晰成像、视场大、对运动目标敏感等特点。虽然传统的单目相机已经能够获得高清的图像,但单目相机的体积以及焦距都很难进一步缩小。人工复眼相机作为一种新型的相机类型,通过模拟昆虫复眼成像原理,能在保持与单目相机相比拟的成像效果的同时具有更小的体积。因此人工复眼相机在运动目标估计、三维重建等方面具有重要的研究价值。其中运动估计是计算机视觉领域的一个重要的研究内容,光流估计是运动估计研究中一个重要的研究方向。
虽然目前已有大量的光流估计方法,但在人工复眼相机图像的应用中仍然存在许多问题。人工复眼相机由于其多个小视场角的低分辨率成像通道的特点,虽然相邻子眼图像间存在交叠,但由于缺少有效的支撑信息,采用目前最优的光流估计方法对人工复眼图像进行光流估计仍会存在图像边缘估计结果较差、以及相邻子眼对应区域光流场不一致的问题。因此,发明一种适用于人工复眼相机的光流估计方法具有十分重要的意义,有利于将人工复眼相机应用到实际环境中。
发明内容
本发明目的是针对人工复眼相机提供一种图像阵列光流估计方法,该方法根据所用人工复眼的结构特点,对阵列图像中包含运动目标的图像进行检测,随后构建光流估计能量函数,并利用变分法对能量函数进行优化,完成人工复眼图像的光流估计。
实现本发明目的的技术方案如下:
一种人工复眼相机的图像阵列光流估计方法,包括步骤
S1、获取人工复眼图像:令人工复眼相机具有M×N个子眼,组成一个子眼阵列,i为二维向量表示表示阵列中的一个子眼,则在时间t和t+1分别捕捉获得M×N个图像Ii(x,y,t)和Ii(x,y,t+1);
S2、检测每个子眼图像是否包含运动目标:子眼i在时间t和t+1的图像Ii(x,y,t)和Ii(x,y,t+1)关于时间t的绝对差值dt(x,y)公式为:
dt(x,y)=|Ii(x,y,t)-Ii(x,y,t+1)|
如果dt(x,y)大于阈值τ,则像素(x,y)是运动目标的投影点;如果子眼图像中运动目标的投影点数量大于某一常数Mon,则该子眼捕捉图像中包含运动目标;
S3、令x=(x,y,t)T表示像素(x,y)在时间t的坐标,w=(u,v,1)T表示时间t和t+1间的位移向量;其中u和v分别表示水平和垂直位移量,则光流估计能量函数E(u,v)计算方法为:
S31、对任意子眼i,计算图像数据匹配项
Figure BDA00015685611900000210
Figure BDA0001568561190000021
其中
Figure BDA0001568561190000022
表示空间梯度,Ψ是惩罚函数,
Figure BDA0001568561190000023
ε为一极小正常数以保证惩罚函数为凸函数,γ是权重,Ω表示图像域;
S32、计算平滑项
Figure BDA00015685611900000211
Figure BDA0001568561190000024
其中
Figure BDA0001568561190000025
当只有两张图像,则
Figure BDA0001568561190000026
变为▽;
S33、进一步利用相邻子眼图像的对应信息,计算空间对应项
Figure BDA00015685611900000212
Figure BDA0001568561190000027
其中Ne(i)表示子眼i的4-最近邻子眼,xj表示x在子眼j的对应像素,wj表示其位移向量,δ(xj)是狄拉克函数;如果xj存在则为1,否则为0;g(xj)表示权重函数;
S34、构建光流估计能量函数E(u,v):
Figure BDA0001568561190000028
其中Ei(u,v)表示子眼图像Ii的能量函数,
Figure BDA0001568561190000029
其中α和β分别表示平滑项和空间对应项的权重;
S4、利用变分法对能量函数进行最小化:
S41、首先定义如下的缩写:
Figure BDA0001568561190000031
其中,z替换了t;
S42、根据欧拉-拉格朗日方程和S41的缩写,则Ei(u,v)的欧拉-拉格朗日方程为:
Figure BDA0001568561190000032
Figure BDA0001568561190000033
S43、对光流的估计,采用标准的由粗到精的图像金字塔原理,设每个子眼图像包含r层,其中原图像是第r层,其中每一层降采样因子为ρ,从第1层到第r层,利用固定点迭代方法计算对应的w;设wk=(uk,vk,1)T,k=0,1,...,表示在k次迭代的位移向量,其中在第1层的第0次迭代初始w0=(0,0,1)T,第2-r层的初始位移向量来自前一层的估计结果;因此第k+1次迭代时,wk+1由下式获得,
Figure BDA0001568561190000034
Figure BDA0001568561190000035
另外,设uk+1=uk+Δuk,vk+1=vk+Δvk,以及利用
Figure BDA0001568561190000036
的近似表示,
Figure BDA0001568561190000041
以及将式(1)中数据项、平滑项和空间对应项的导数缩写为以下形式,
Figure BDA0001568561190000042
则将式(2)代入式(3),可得
Figure BDA0001568561190000043
将(2)和(4)代入(1),可得,
Figure BDA0001568561190000044
Figure BDA0001568561190000045
S44、再次采用内部固定点迭代方法,设内部迭代初始值为Δuk,0=0,Δvk,0=0,则在第l+1次迭代,在位移向量增量Δuk,l+1和Δvk,l+1的系统方程为,
Figure BDA0001568561190000046
Figure BDA0001568561190000047
线性系统方程通过逐次超松弛方法求解;
S5、光流估计:对每个子眼的每一层图像,外部固定点迭代k次,每次迭代后使用5×5的加权中值滤波移除噪声和保持运动边界,随后传到下一层,直到第r层,最后获得每个子眼光流估计结果(u,v)。
本发明与现有技术相比,其显著优点在于:能有效地对人工复眼图像序列进行光流估计,并且由于子眼图像间空间相关的约束项,能够对相邻子眼对应区域进行约束,获得一致的光流估计结果,解决人工复眼相机中子眼视场角小导致的图像边缘光流估计不准确的问题。
附图说明
图1是本发明实施例中采用人工复眼相机图像阵列光流估计方法的实验效果图。
其中1是t时刻采集图像,2是t+1时刻图像,3是本发明方法对1和2中方框区域的图像的光流估计结果,4是采用原变分法的光流估计结果。
具体实施方式
下面通过实施例对本发明作进一步详细描述。
一种人工复眼相机的图像阵列光流估计方法,该方法根据所用人工复眼的结构特点,构建光流估计能量函数,有效对相邻子眼图像对应区域进行约束,并利用变分法对光流进行估计,完成人工复眼图像光流估计。本实施例以一种特定的人工复眼相机(电子簇眼,Electronic cluster eye)为实验对象,对目标深度信息进行估计。实例中所用电子簇眼成像通道数为17*13,由于失真影响,只采用其中心13*13个通道,本发明需采取如下步骤:
S1、利用电子簇眼连续采集两帧图像,并提取电子簇眼的13*13个子眼图像Ii(x,y,t)和Ii(x,y,t+1)。提取图像包括图像大小和颜色通道等信息,每个子眼图像修正后分辨率为110*110像素。
S2、检测每个子眼图像是否包含运动目标,假设子眼i在时间t和t+1的图像Ii(x,y,t)和Ii(x,y,t+1)关于时间t的绝对差值dt(x,y)公式为:
dt(x,y)=|Ii(x,y,t)-Ii(x,y,t+1)|,如果dt(x,y)大于阈值τ=0.045,则像素(x,y)被认为是运动目标的投影点。如果子眼图像中运动目标的投影点数量大于某一常数Mon=30,则认为该子眼捕捉图像中包含运动目标;S3、首先假设x=(x,y,t)T表示像素在时间t的坐标,w=(u,v,1)T表示时间t和t+1间的位移向量,其中u和v分别表示水平和垂直位移量。则光流估计能量函数E(u,v)具体计算方法为:
S31、对任意子眼i,按下式计算图像数据匹配项
Figure BDA0001568561190000051
Figure BDA0001568561190000052
其中
Figure BDA0001568561190000053
表示空间梯度,Ψ是惩罚函数,
Figure BDA0001568561190000054
ε=0.001,γ=20;
S32、按以下公式计算平滑项
Figure BDA0001568561190000061
Figure BDA0001568561190000062
其中
Figure BDA0001568561190000063
如果只有两张图像,则
Figure BDA0001568561190000064
就变为
Figure BDA0001568561190000065
S33、进一步利用相邻子眼图像的对应信息,按以下公式计算空间对应项
Figure BDA0001568561190000066
Figure BDA0001568561190000067
其中Ne(i)表示子眼i的4-最近邻子眼,xj表示x在子眼j的对应像素,wj则表示其位移向量,δ(xj)是狄拉克函数,如果xj存在则为1,否则为0。g(xj)表示以子眼图像中心为中心点的高斯分布权重函数;
S34、按下列公式构建光流估计能量函数E(u,v):
Figure BDA0001568561190000068
其中Ei(u,v)表示子眼图像Ii的能量函数,计算方法为:
Figure BDA0001568561190000069
其中α=1,β=2。
S4、利用变分法对能量函数进行最小化:
S41、首先定义如下的缩写:
Figure BDA00015685611900000610
其中,z替换了t。
S42、根据欧拉-拉格朗日方程和S41的缩写,则Ei(u,v)的欧拉-拉格朗日方程为
Figure BDA00015685611900000611
Figure BDA00015685611900000612
S43、对光流的估计,本发明采用标准的由粗到精的图像金字塔原理,假设每个子眼图像包含r=4层,其中原图像是第r层,其中每一层降采样因子为ρ=0.8,从第1层到第r层,利用固定点迭代方法计算对应的w,假设wk=(uk,vk,1)T,k=0,1,...,表示在k次迭代的位移向量,其中在第1层的第0次迭代初始w0=(0,0,1)T,第2-r层的初始位移向量来自前一层的估计结果。因此第k+1次迭代时,wk+1可以由下式获得,
Figure BDA0001568561190000071
Figure BDA0001568561190000072
另外,假设uk+1=uk+Δuk,vk+1=vk+Δvk,以及利用
Figure BDA0001568561190000073
的近似表示,
Figure BDA0001568561190000074
以及将式(1)中数据项、平滑项和空间对应项的导数缩写为以下形式,
Figure BDA0001568561190000075
则将式(2)代入式(3),可得
Figure BDA0001568561190000076
将(2)和(4)代入(1)可得,
Figure BDA0001568561190000081
Figure BDA0001568561190000082
S44、为了估计每一次迭代后的位移向量增量,我们再次采用内部固定点迭代方法,假设内部迭代初始值为Δuk,0=0,Δvk,0=0,则在第l+1次迭代,在位移向量增量Δuk,l+1和Δvk,l+1的系统方程为,
Figure BDA0001568561190000083
Figure BDA0001568561190000084
于是该线性系统方程可以通过逐次超松弛方法求解。
S5、对每个字眼的每一层图像,外部固定点迭代k=3次,每次迭代后使用5×5的加权中值滤波移除噪声和保持运动边界,随后传到下一层,直到第r层,最后获得每个子眼光流估计结果(u,v)。
对实验结果与目前现有变分法光流估计进行比较如下:
利用middlebury数据库中的图像对本发明和原变分法进行比较,表1所示middlebury图像数据集中的光流估计结果的平均角点误差(EPE)和平均角度误差(AAE),其值越小越好。
表1 平均角点误差和平均角度误差结果
Figure BDA0001568561190000085
由表1可以看出,本发明具有更好的效果。图1所示为本发明实施例和现有变分法对光流估计结果的定性比较。由图1可以看出,本发明估计的光流结果具有更好的一致性效果。

Claims (1)

1.一种人工复眼相机的图像阵列光流估计方法,其特征在于,包括步骤
S1、获取人工复眼图像:令人工复眼相机具有M×N个子眼,组成一个子眼阵列,i为二维向量表示表示阵列中的一个子眼,则在时间t和t+1分别捕捉获得M×N个图像Ii(x,y,t)和Ii(x,y,t+1);
S2、检测每个子眼图像是否包含运动目标:子眼i在时间t和t+1的图像Ii(x,y,t)和Ii(x,y,t+1)关于时间t的绝对差值dt(x,y)公式为:
dt(x,y)=|Ii(x,y,t)-Ii(x,y,t+1)|
如果dt(x,y)大于阈值τ,则像素(x,y)是运动目标的投影点;如果子眼图像中运动目标的投影点数量大于某一常数Mon,则该子眼捕捉图像中包含运动目标;
S3、令x=(x,y,t)T表示像素(x,y)在时间t的坐标,w=(u,v,1)T表示时间t和t+1间的位移向量;其中u和v分别表示水平和垂直位移量,则光流估计能量函数E(u,v)计算方法为:
S31、对任意子眼i,计算图像数据匹配项
Figure FDA0002379283940000011
Figure FDA0002379283940000012
其中
Figure FDA0002379283940000013
表示空间梯度,Ψ是惩罚函数,
Figure FDA0002379283940000014
ε为一极小正常数以保证惩罚函数为凸函数,γ是权重,Ω表示图像域;
S32、计算平滑项
Figure FDA0002379283940000015
Figure FDA0002379283940000016
当只有两张图像,则
Figure FDA0002379283940000017
变为
Figure FDA0002379283940000018
S33、进一步利用相邻子眼图像的对应信息,计算空间对应项
Figure FDA0002379283940000019
Figure FDA00023792839400000110
其中Ne(i)表示子眼i的4-最近邻子眼,xj表示x在子眼j的对应像素,wj表示其位移向量,δ(xj)是狄拉克函数;如果xj存在则为1,否则为0;g(xj)表示权重函数;
S34、构建光流估计能量函数E(u,v):
Figure FDA00023792839400000111
其中Ei(u,v)表示子眼图像Ii的能量函数,
Figure FDA0002379283940000021
其中α和β分别表示平滑项和空间对应项的权重;
S4、利用变分法对能量函数进行最小化:
S41、首先定义如下的缩写:
Figure FDA0002379283940000022
其中,z替换了t;
S42、根据欧拉-拉格朗日方程和S41的缩写,则Ei(u,v)的欧拉-拉格朗日方程为:
Figure FDA0002379283940000023
Figure FDA0002379283940000024
S43、对光流的估计,采用标准的由粗到精的图像金字塔原理,设每个子眼图像包含r层,其中原图像是第r层,其中每一层降采样因子为ρ,从第1层到第r层,利用固定点迭代方法计算对应的w;设wk=(uk,vk,1)T,k=0,1,...,表示在k次迭代的位移向量,其中在第1层的第0次迭代初始w0=(0,0,1)T,第2-r层的初始位移向量来自前一层的估计结果;因此第k+1次迭代时,wk+1由下式获得,
Figure FDA0002379283940000031
Figure FDA0002379283940000032
另外,设uk+1=uk+Δuk,vk+1=vk+Δvk,以及利用
Figure FDA0002379283940000033
的近似表示,
Figure FDA0002379283940000034
Figure FDA0002379283940000035
Figure FDA0002379283940000036
以及将式(1)中数据项、平滑项和空间对应项的导数缩写为以下形式,
Figure FDA0002379283940000037
Figure FDA00023792839400000313
Figure FDA0002379283940000038
则将式(2)代入式(3),可得
Figure FDA00023792839400000310
Figure FDA00023792839400000311
Figure FDA00023792839400000312
将(2)和(4)代入(1),可得,
Figure FDA0002379283940000041
Figure FDA0002379283940000042
S44、再次采用内部固定点迭代方法,设内部迭代初始值为Δuk,0=0,Δvk,0=0,则在第l+1次迭代,在位移向量增量Δuk,l+1和Δvk,l+1的系统方程为,
Figure FDA0002379283940000043
Figure FDA0002379283940000044
该系统方程通过逐次超松弛方法求解;
S5、光流估计:对每个子眼的每一层图像,外部固定点迭代k次,每次迭代后使用5×5的加权中值滤波移除噪声和保持运动边界,随后传到下一层,直到第r层,最后获得每个子眼光流估计结果(u,v)。
CN201810108789.0A 2018-02-05 2018-02-05 一种人工复眼相机的图像阵列光流估计方法 Expired - Fee Related CN108399629B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810108789.0A CN108399629B (zh) 2018-02-05 2018-02-05 一种人工复眼相机的图像阵列光流估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810108789.0A CN108399629B (zh) 2018-02-05 2018-02-05 一种人工复眼相机的图像阵列光流估计方法

Publications (2)

Publication Number Publication Date
CN108399629A CN108399629A (zh) 2018-08-14
CN108399629B true CN108399629B (zh) 2020-06-02

Family

ID=63095293

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810108789.0A Expired - Fee Related CN108399629B (zh) 2018-02-05 2018-02-05 一种人工复眼相机的图像阵列光流估计方法

Country Status (1)

Country Link
CN (1) CN108399629B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110738610A (zh) * 2019-09-19 2020-01-31 天津大学 一种复眼图像像质提升方法
CN116342660B (zh) * 2022-09-08 2023-12-08 中北大学 一种多尺度分析融合加权滤波的仿生复眼光流场估计方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1697982A (zh) * 2003-03-14 2005-11-16 C.R.F.阿西安尼顾问公司 用于探测障碍、特别用于自主导航的有源电光装置
CN103247058A (zh) * 2013-05-13 2013-08-14 北京工业大学 一种基于误差分布式多层网格的快速光流场计算方法
WO2015116640A1 (en) * 2014-01-29 2015-08-06 Shazly Tarek A Eye and head tracking device
CN107451952A (zh) * 2017-08-04 2017-12-08 追光人动画设计(北京)有限公司 一种全景视频的拼接融合方法、设备以及系统
CN107507228A (zh) * 2017-06-15 2017-12-22 清华大学 一种基于光流的3d视觉生成方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7929800B2 (en) * 2007-02-06 2011-04-19 Meadow William D Methods and apparatus for generating a continuum of image data
US8233054B2 (en) * 2006-09-25 2012-07-31 Pioneer Corporation Scenery imaging apparatus, scenery imaging method, scenery imaging program, and computer-readable recording medium

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1697982A (zh) * 2003-03-14 2005-11-16 C.R.F.阿西安尼顾问公司 用于探测障碍、特别用于自主导航的有源电光装置
CN103247058A (zh) * 2013-05-13 2013-08-14 北京工业大学 一种基于误差分布式多层网格的快速光流场计算方法
WO2015116640A1 (en) * 2014-01-29 2015-08-06 Shazly Tarek A Eye and head tracking device
CN107507228A (zh) * 2017-06-15 2017-12-22 清华大学 一种基于光流的3d视觉生成方法
CN107451952A (zh) * 2017-08-04 2017-12-08 追光人动画设计(北京)有限公司 一种全景视频的拼接融合方法、设备以及系统

Also Published As

Publication number Publication date
CN108399629A (zh) 2018-08-14

Similar Documents

Publication Publication Date Title
CN106651938B (zh) 一种融合高分辨率彩色图像的深度图增强方法
Yue et al. Image super-resolution: The techniques, applications, and future
CN107025632B (zh) 一种图像超分辨率重建方法及系统
CN108986136B (zh) 一种基于语义分割的双目场景流确定方法及系统
CN110503680B (zh) 一种基于非监督的卷积神经网络单目场景深度估计方法
CN109190581B (zh) 图像序列目标检测识别方法
US9692939B2 (en) Device, system, and method of blind deblurring and blind super-resolution utilizing internal patch recurrence
CN103824273B (zh) 基于复合运动和自适应非局部先验的超分辨率重建方法
KR100985805B1 (ko) 적응적인 칼만필터를 이용한 영상 안정화 장치 및 방법
JP2007257287A (ja) 画像レジストレーション方法
JP5107409B2 (ja) 動き領域の非線形スムージングを用いた動き検出方法およびフィルタリング方法
Su et al. Super-resolution without dense flow
CN109658361B (zh) 一种顾及运动估计误差的运动场景超分辨率重建方法
AU2005280315A1 (en) Method and system for motion correction in a sequence of images
Rochefort et al. An improved observation model for super-resolution under affine motion
CN108399629B (zh) 一种人工复眼相机的图像阵列光流估计方法
CN111145094A (zh) 一种基于表面法向引导与图拉普拉斯先验约束的深度图增强方法
CN112927251A (zh) 基于形态学的场景稠密深度图获取方法、系统及装置
CN111861888A (zh) 图像处理方法、装置、电子设备及存储介质
CN104376544B (zh) 一种基于多区域尺度放缩补偿的非局部超分辨率重建方法
CN110580715A (zh) 一种基于照度约束和格网变形的图像对齐方法
Al Ismaeil et al. Real-time non-rigid multi-frame depth video super-resolution
Amiri et al. A fast video super resolution for facial image
Mohan Adaptive super-resolution image reconstruction with lorentzian error norm
Cho et al. Depth map up-sampling using cost-volume filtering

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
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20200602

Termination date: 20210205