CN105574893B - 一种基于投影梯度的稳像方法 - Google Patents

一种基于投影梯度的稳像方法 Download PDF

Info

Publication number
CN105574893B
CN105574893B CN201510958272.7A CN201510958272A CN105574893B CN 105574893 B CN105574893 B CN 105574893B CN 201510958272 A CN201510958272 A CN 201510958272A CN 105574893 B CN105574893 B CN 105574893B
Authority
CN
China
Prior art keywords
mrow
msub
curve
image block
projection
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
CN201510958272.7A
Other languages
English (en)
Other versions
CN105574893A (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.)
Zhejiang University ZJU
Original Assignee
Zhejiang University ZJU
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 Zhejiang University ZJU filed Critical Zhejiang University ZJU
Priority to CN201510958272.7A priority Critical patent/CN105574893B/zh
Publication of CN105574893A publication Critical patent/CN105574893A/zh
Application granted granted Critical
Publication of CN105574893B publication Critical patent/CN105574893B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10016Video; Image sequence

Landscapes

  • Control Of Indicators Other Than Cathode Ray Tubes (AREA)
  • Transforming Electric Information Into Light Information (AREA)

Abstract

本发明公开了一种基于投影梯度的稳像方法。针对室内拍摄的高帧频视频中帧间亮度不均以及高噪声的情况,引入图像投影梯度,并利用阈值判别函数作为灰度位移函数和梯度位移函数的权重函数,其中判别函数的输入为参考帧与当前帧总灰度值的相对灰度差,最后获得加权图像位移矢量。通过使用投影梯度以及阈值判断,本发明提出的方法能很好地解决灰度投影在帧间亮度不均以及图像高噪声时误检问题,从而得到更好的稳像效果。

Description

一种基于投影梯度的稳像方法
技术领域
本发明属于图像处理领域,尤其涉及一种基于投影梯度的稳像方法。
背景技术
灰度投影算法常用于视频稳像以及图像匹配,相较于其它算法,灰度投影的优势在于算法简单、时间复杂度低,具有很好的实时性。以灰度投影为基础,之后出现了分区灰度投影、多方向投影等一系列改进型灰度投影算法,进一步增强了算法的实时性和鲁棒性。
通常情况下,用于图像匹配和稳像的视频序列为自然光下拍摄的视频或低帧频(30Hz)的视频。由于自然光无闪烁,室内照明灯闪烁频率一般为50-70Hz,因此此时拍摄的视频中各帧图像亮度恒定,或亮度变化很小。然而现实情况中,有时会遇到需要拍摄室内照明或夜间灯光照明下的高速视频,这种情况随着摄影设备的轻型化和高速化越发普遍,比如如今市场上的主流手机都具备300Hz以上的高速摄影能力。对于闪烁灯光下的高速摄影,视频图像各帧之间会产生亮度不均的问题,如附图1所示,并且该问题会随着拍摄帧频与照明闪烁频率差异变大而更加明显。此外,视频帧频的增加也会导致图像噪声增大。
传统的稳像与匹配算法均无法解决帧间亮度不均以及噪声大带来的配准误差:基于灰度投影的算法依赖于图像灰度,亮度不均会导致参考帧与当前帧的灰度投影存在y方向(灰度方向)上的位移,如附图2(a)所示,而亮度归一化处理也会因为大噪声的存在效果受到影响,通常位移矢量误差在0.5个像素左右,如附图2(b)所示;基于光流的匹配算法由于光流方程建立在亮度恒定的基础上,所以无法处理亮度非恒定的视频图像;而基于特征的匹配算法,尽管对亮度和噪声不敏感,但相对高昂的时间和计算复杂度限制了算法的实时性。
梯度作为信号的基本特征之一,反应了信号的变化率,具有与信号本身大小无关的特性,可解决信号幅值平移的匹配问题。对于图像来说,图像梯度反映了图像的边缘信息,即高频信息。但图像梯度在突出边缘信息的同时,也丢失了低频信息,这使得梯度投影曲线相较于灰度投影曲线波动剧烈,如附图3所示,因而更易产生匹配误差。此外,计算图像梯度也增加了算法的时间,降低了算法的实时性。
发明内容
本发明的目的是在视频图像的匹配和稳像中,针对室内拍摄的高帧频视频中帧间亮度不均以及高噪声图像,提出一种基于投影梯度的稳像方法,利用图像投影梯度,通过相对灰度差的判别函数对灰度投影位移函数和梯度位移函数进行加权,获得加权位移量。
本发明的目的是通过以下技术方案来实现的:一种基于投影梯度的稳像方法,该方法包括以下步骤:
(1)图像分块,具体是:
(1.1)将参考图像与当前图像分为对应的K个大小为N×N的图像块。
(2)计算图像块的灰度投影曲线,具体是:
(2.1)水平方向上的灰度投影为所有行的像素值在列方向上的累积;垂直方向上的灰度投影为所有列像素值在行方向上的累积;
(2.2)去除方差小于方差阈值的图像块的灰度投影曲线。
(3)将灰度投影曲线中的第n+1项减第n-1项,得到投影梯度曲线;
(4)互相关运算,具体是:
(4.1)利用互相关公式对参考帧与当前帧对应图像块的水平灰度投影曲线和垂直投影曲线分别做互相关计算,得到互相关曲线,根据互相关曲线的谷值确定当前帧相对于参考帧的灰度位移矢量Mkp
(4.2)利用互相关公式对参考帧与当前帧对应图像块的水平投影梯度曲线和垂直投影梯度曲线分别做互相关计算,得到互相关曲线,根据互相关曲线的谷值确定当前帧相对于参考帧的梯度位移矢量Mkg
(5)计算判别函数,具体是:
(5.1)计算参考帧与当前帧对应图像块总灰度的相对灰度差rk
其中,Ikc(i,j)为当前帧的第k个图像块第i行第j列的灰度值,Ikr(i,j)为参考帧的第k个图像块第i行第j列的的灰度值;
(5.2)将步骤(5.1)中计算所得rk作为判别函数λ(r)的输入:
其中,α为坡度参数,β为阈值参数。
(6)计算位移矢量Mk,公式如下:
Mk=(1-λ(r))Mkp+λ(r)Mkg (3)
去除图像块位移矢量Mk中由于步骤2.2中的误差和局部运动造成的偏差位移矢量,将剩余位移矢量求平均得到当前帧与参考帧的位移矢量M。
其中,l1为偏差位移矢量的个数,l2为步骤2.2中方差小于方差阈值的图像块个数。
(7)将当前帧向参考帧移动-M,实现当前帧的稳像。
进一步地,步骤1.1中,所述图像块大小N为128。
进一步地,步骤5.2中,所述α=250,β=10。
本发明的有益效果:本发明针对室内拍摄的高帧频视频中帧间亮度不均以及高噪声的情况,提出了一种基于投影梯度的稳像方法。首先,本发明引入图像投影梯度用于处理亮度不等的问题,使用投影梯度而非梯度投影的优势在于投影梯度不需要对二维图像求梯度,而只需对一维的灰度投影曲线求梯度,因此不会增加额外的计算量,保证了算法的实时性。其次,针对梯度曲线波动大、对噪声敏感的问题,采用分块处理的方式,从而减少了曲线波动带来的相关运算误差。最后,利用相对灰度差的判别函数来判断当前帧和参考帧图像是否存在亮度不等的问题,判别函数为一平滑函数,并以权重的形式对灰度位移矢量和梯度位移矢量加权求和,从而降低了位移矢量的绝对误差。本发明提出的方法在不增加算法计算时间的前提下能很好地解决灰度投影在帧间亮度不均以及图像高噪声时误检问题,从而得到更好的稳像效果。
附图说明
图1(a)为室内照明条件下的高速视频参考帧。
图1(b)为室内照明条件下的高速视频当前帧。
图2(a)为图1(a)所示图像垂直方向灰度投影曲线。
图2(b)为图1(a)所示图像亮度归一化后垂直方向灰度投影曲线。
图3为图1(a)所示图像的垂直梯度投影曲线。
图4为本发明方法流程图。
图5(a)为参考图像块。
图5(b)为当前图像块。
图6(a)为图5所示图像块垂直方向灰度投影曲线。
图6(b)为图5所示图像块垂直方向投影梯度曲线。
图7(a)为图5所示图像块垂直方向灰度投影曲线的SAD曲线。
图7(b)为图5所示图像块垂直方向投影梯度曲线的SAD曲线。
图8为判别函数。
具体实施方式
以下结合附图对本发明作进一步说明。
如图4所示,本发明提供的一种基于投影梯度的稳像方法,主要包括图像分块、计算投影曲线、互相关运算、阈值计算和位移矢量计算等几个步骤。
步骤1.图像分块
1-1为了减小互相关运算时投影梯度曲线波动带来的估计误差,将参考图像与当前图像分为对应的K个不重叠图像块,图像块大小为N×N,如附图5所示。N取值太大将无法降低曲线波动造成的位移矢量估计误差,取值太小则会使步骤4中的相对灰度差估计不准确,因此这里N=128。
步骤2.计算图像块的灰度投影曲线
2-1:对当前帧图像和参考帧图像进行预处理,将图像灰度值投影成为2个独立的一维曲线,如附图6(a)所示,其投影公式为:
Gkcurv(j),Gkcurh(i)分别为当前帧第k个图像块的第j列和第i行的投影值,xkcur(i,j)为当前帧第k个图像块的第i行、第j列的像素值。参考帧图像的投影曲线与之类似。
为防止均匀图像块造成的误差,去除方差小于方差阈值的投影曲线,具体是:对投影曲线计算中心化后的方差,若方差小于设定的方差阈值,则认为图像块为均匀图像块,其投影曲线可信度较低,因此去除。
步骤3.计算灰度投影的梯度曲线
3-1利用梯度曲线解决帧间图像亮度不等的问题。为避免计算图像梯度带来的额外计算、保证算法实时性,这里使用投影梯度而非梯度投影,从而仅需对已求得的一维灰度投影曲线G求梯度,而不需要先计算二维图像梯度再进行梯度投影。
定义投影梯度曲线为灰度投影曲线G中第n+1项减第n-1项,并计算得到投影梯度曲线,如图6(b)所示,计算公式为:
Dkcurv(j)=Gkcurv(j+1)-Gkcurv(j-1) (3)
Dkcurh(i)=Gkcurh(i+1)-Gkcurh(i-1) (4)
Dkcurv,Dkcurh分别为当前帧图像块垂直方向和水平方向灰度投影的梯度曲线。参考帧图像的投影梯度曲线计算过程与之类似。容易证明投影梯度与梯度投影成正比,(3)式可进一步写为:
令GDkcurv为对应图像块的垂直梯度投影,则有:
其中Ikcur为当前帧第k图像块,为卷积符号,Pv为垂直Prewitt算子:
由式(5)和式(6)可知,投影梯度与梯度投影间仅存在一个倍数。此外,容易证明对其它边缘算子该关系仍成立,如Sobel算子、Laplace算子等。
步骤4.互相关运算
分别对当前帧与参考帧的行、列灰度投影曲线做互相关计算,即计算投影曲线中对应元素差的绝对值SAD(Sum of absolute differences)。相关曲线计算公式为:
式中,-m≤ω≤m;Ck(ω)为当前帧图像块k相对参考帧图像块k偏移为ω时的相关量;Gkr为参考帧水平方向或者垂直方向第j个灰度投影;Gkc为当前帧水平方向或者垂直方向第j个灰度投影;m为搜索范围最大值;N为图像块大小。
根据两条相关曲线的谷值确定当前帧相对于参考帧的灰度位移矢量Mkp,灰度位移矢量可通过其它改进型投影方法获取,以提高准确率。投影梯度曲线相关运算过程与之类似,如附图7所示。
步骤5.计算判别函数
首先计算参考帧与当前帧对应图像块总灰度的相对灰度差rk
定义判别函数λ为:
其中,α为坡度参数,用于调节函数曲线过度区域的范围,β为阈值参数,用于控制阈值大小,这里α=250,β=10,如附图8所示。计算所得rk作为判别函数λ(r)的输入。
步骤6.计算位移矢量
分别将灰度位移矢量Mkp和梯度位移矢量Mkg乘以大小为(1-λ)和λ的权重然后相加,其和为当前帧与参考帧对应图像块的位移矢量Mk,即当两帧图像总像素和差别大于阈值时,结果更偏向于梯度位移矢量,反之偏向于灰度位移矢量。
Mk=(1-λ(r))Mkp+λ(r)Mkg (11)
在视频稳像中,需要去除计算得到的图像块位移矢量Mk中由于误差和局部运动造成的偏差位移矢量,将剩余位移矢量求平均得到当前帧与参考帧的位移矢量M。
其中,l1为偏差位移矢量的个数,l2为步骤2中方差小于方差阈值的图像块个数。
步骤7.将当前帧向参考帧移动-M,实现当前帧的稳像。
表1为直接使用灰度投影、亮度归一化、以及本发明提出的方法在亮度恒定情况和高噪声亮度非恒定情况下的位移矢量平均准确率。
表1位移矢量准确率比较

Claims (3)

1.一种基于投影梯度的稳像方法,其特征在于,该方法包括以下步骤:
(1)图像分块,具体是:
(1.1)将参考图像与当前图像分为对应的K个大小为N×N的图像块;
(2)计算图像块的灰度投影曲线,具体是:
(2.1)水平方向上的灰度投影为所有行的像素值在列方向上的累积;垂直方向上的灰度投影为所有列像素值在行方向上的累积;
(2.2)去除方差小于方差阈值的图像块的灰度投影曲线;
(3)将灰度投影曲线中的第n+1项减第n-1项,得到投影梯度曲线;
(4)互相关运算,具体是:
(4.1)利用互相关公式对参考帧与当前帧对应图像块的水平灰度投影曲线和垂直灰度投影曲线分别做互相关计算,得到互相关曲线,根据互相关曲线的谷值确定当前帧相对于参考帧的灰度位移矢量Mkp
(4.2)利用互相关公式对参考帧与当前帧对应图像块的水平投影梯度曲线和垂直投影梯度曲线分别做互相关计算,得到互相关曲线,根据互相关曲线的谷值确定当前帧相对于参考帧的梯度位移矢量Mkg
(5)计算判别函数,具体是:
(5.1)计算参考帧与当前帧对应图像块总灰度的相对灰度差rk
<mrow> <msub> <mi>r</mi> <mi>k</mi> </msub> <mo>=</mo> <mo>|</mo> <mfrac> <mrow> <munder> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </munder> <msub> <mi>I</mi> <mrow> <mi>k</mi> <mi>c</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>-</mo> <munder> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </munder> <msub> <mi>I</mi> <mrow> <mi>k</mi> <mi>r</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> </mrow> <mrow> <munder> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </munder> <msub> <mi>I</mi> <mrow> <mi>k</mi> <mi>r</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>|</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
其中,Ikc(i,j)为当前帧的第k个图像块第i行第j列的灰度值,Ikr(i,j)为参考帧的第k个图像块第i行第j列的的灰度值;
(5.2)将步骤(5.1)中计算所得rk作为判别函数λ(r)的输入:
<mrow> <mi>&amp;lambda;</mi> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>1</mn> <mo>+</mo> <mi>exp</mi> <mrow> <mo>(</mo> <mi>&amp;beta;</mi> <mo>-</mo> <mi>&amp;alpha;</mi> <mi>r</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
其中,α为坡度参数,β为阈值参数;
(6)计算图像块位移矢量Mk,公式如下:
Mk=(1-λ(r))Mkp+λ(r)Mkg (3)
去除图像块位移矢量Mk中由于步骤(2.2)中的误差和局部运动造成的偏差位移矢量,将剩余位移矢量求平均得到当前帧与参考帧的位移矢量M;
<mrow> <mi>M</mi> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mi>K</mi> <mo>-</mo> <msub> <mi>l</mi> <mn>1</mn> </msub> <mo>-</mo> <msub> <mi>l</mi> <mn>2</mn> </msub> </mrow> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mi>k</mi> <mrow> <mi>K</mi> <mo>-</mo> <msub> <mi>l</mi> <mn>1</mn> </msub> <mo>-</mo> <msub> <mi>l</mi> <mn>2</mn> </msub> </mrow> </munderover> <msub> <mi>M</mi> <mi>k</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>
其中,l1为偏差位移矢量的个数,l2为步骤(2.2)中方差小于方差阈值的图像块个数;
(7)将当前帧向参考帧移动-M,实现当前帧的稳像。
2.根据权利要求1所述一种基于投影梯度的稳像方法,其特征在于,步骤(1.1)中,所述图像块大小N为128。
3.根据权利要求1所述一种基于投影梯度的稳像方法,其特征在于,步骤(5.2)中,所述α=250,β=10。
CN201510958272.7A 2015-12-18 2015-12-18 一种基于投影梯度的稳像方法 Expired - Fee Related CN105574893B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510958272.7A CN105574893B (zh) 2015-12-18 2015-12-18 一种基于投影梯度的稳像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510958272.7A CN105574893B (zh) 2015-12-18 2015-12-18 一种基于投影梯度的稳像方法

Publications (2)

Publication Number Publication Date
CN105574893A CN105574893A (zh) 2016-05-11
CN105574893B true CN105574893B (zh) 2018-02-09

Family

ID=55884980

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510958272.7A Expired - Fee Related CN105574893B (zh) 2015-12-18 2015-12-18 一种基于投影梯度的稳像方法

Country Status (1)

Country Link
CN (1) CN105574893B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106357956A (zh) * 2016-08-25 2017-01-25 哈尔滨工程大学 一种uuv横垂荡时的光视觉图像稳像方法
CN108805832B (zh) * 2018-05-29 2022-02-08 重庆大学 适于隧道环境特性的改进灰度投影稳像方法
CN110581933A (zh) * 2019-07-23 2019-12-17 国网山东省电力公司电力科学研究院 一种视频优化处理系统、处理方法、终端设备及存储介质
CN114298984B (zh) * 2021-12-16 2023-12-08 苏州凌云光工业智能技术有限公司 一种屏幕贯穿线的检测方法、装置、电子设备及存储介质

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103475802A (zh) * 2013-09-26 2013-12-25 中国矿业大学 一种电子稳像方法
CN103955930A (zh) * 2014-04-28 2014-07-30 中国人民解放军理工大学 基于灰度积分投影互相关函数特征的运动参数估计方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103620621B (zh) * 2011-06-30 2017-10-24 诺基亚技术有限公司 用于利用积分梯度投影的面部跟踪的方法和装置

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103475802A (zh) * 2013-09-26 2013-12-25 中国矿业大学 一种电子稳像方法
CN103955930A (zh) * 2014-04-28 2014-07-30 中国人民解放军理工大学 基于灰度积分投影互相关函数特征的运动参数估计方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Projection gradient method for energy functional minimization with a constraint and its application to computing the ground state of spin-orbit-coupled Bose-Einstein condensates;Hanquan Wang et al;《Computer Physics Communication》;20140528;第2803-2808页 *
一种基于灰度投影的电子稳像改进算法;张永祥 等;《微电子学与计算机》;20081130;第25卷(第11期);第212-215页 *
基于梯度投影法的电子稳像算法;齐蕴光 等;《计算机工程》;20120331;第38卷(第6期);第230-232、235页 *

Also Published As

Publication number Publication date
CN105574893A (zh) 2016-05-11

Similar Documents

Publication Publication Date Title
CN105574893B (zh) 一种基于投影梯度的稳像方法
US8280165B2 (en) System and method for segmenting foreground and background in a video
CN103455984B (zh) 一种Kinect深度图像获取方法与装置
CN105100807B (zh) 一种基于运动矢量后处理的帧率提升方法
CN104639920B (zh) 基于单帧双次曝光模式的宽动态融合方法
CN105611116B (zh) 一种全局运动矢量估计方法及监控视频稳像方法及装置
CN105931213B (zh) 基于边缘检测和帧差法的高动态范围视频去鬼影的方法
JP4940164B2 (ja) 撮像装置及び撮像方法
JP2011199716A (ja) 画像処理装置、および画像処理方法、並びにプログラム
CN104809716B (zh) 一种确定图像偏移量的方法及设备
CN106210712A (zh) 一种图像坏点检测及处理方法
TWI425445B (zh) 用於判斷一移動平台之自我運動量的方法以及偵測系統
CN108009997B (zh) 一种调整图像对比度的方法及装置
CN108805832B (zh) 适于隧道环境特性的改进灰度投影稳像方法
CN106210448B (zh) 一种视频图像抖动消除处理方法
WO2012027904A1 (zh) 从照片合成肖像素描的系统和方法
CN112866581A (zh) 摄像头自动曝光补偿方法、装置及电子设备
CN106296732B (zh) 一种复杂背景下的运动目标追踪方法
CN104079800A (zh) 一种视频监控中视频图像的抗抖动方法
CN104036471B (zh) 一种图像噪声估值方法及图像噪声估值装置
CN109711267A (zh) 一种行人重识别、行人运动轨迹生成方法及装置
JP2011258179A (ja) 画像処理装置及びその制御方法
TWI394097B (zh) 移動物體的偵測方法以及偵測系統
CN106454061A (zh) 电子装置及影像处理方法
TWI608447B (zh) 立體影像深度圖產生裝置及方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into 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: 20180209

Termination date: 20181218