CN102607460B - 一种应用于三维测量的全局相位滤波方法 - Google Patents

一种应用于三维测量的全局相位滤波方法 Download PDF

Info

Publication number
CN102607460B
CN102607460B CN201210063823.XA CN201210063823A CN102607460B CN 102607460 B CN102607460 B CN 102607460B CN 201210063823 A CN201210063823 A CN 201210063823A CN 102607460 B CN102607460 B CN 102607460B
Authority
CN
China
Prior art keywords
slope
value
phase
point
variable
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
CN201210063823.XA
Other languages
English (en)
Other versions
CN102607460A (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.)
Tianjin Polytechnic University
Original Assignee
Tianjin Polytechnic 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 Tianjin Polytechnic University filed Critical Tianjin Polytechnic University
Priority to CN201210063823.XA priority Critical patent/CN102607460B/zh
Publication of CN102607460A publication Critical patent/CN102607460A/zh
Application granted granted Critical
Publication of CN102607460B publication Critical patent/CN102607460B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Image Processing (AREA)
  • Image Analysis (AREA)

Abstract

本发明属于三维机器视觉领域,涉及一种应用于三维测量的全局相位滤波方法。该方法可以去除全局相位中的噪声干扰,在去除了所有噪声信号的同时,很好的保留了有效的相位信息,使得三维重建的效果更好的包含了边缘等细节信息。滤波后的全局相位用于高精度三维测量中,可以获得更好的三维重建效果。本发明所设计的全局相位滤波方法,除了具有很好的滤波效果之外,也具有更广泛的适用性,不但可以适用于格雷码的编码方式,也可以适用于外差多频以及其他的多频率的编码方式,对于其他类似信号的噪声干扰难题,也有很好的适应性。

Description

一种应用于三维测量的全局相位滤波方法
技术领域
本发明涉及一种应用于三维测量的全局相位滤波方法,更具体的说,本发明涉及一种能够用于高精度三维测量的有效滤除全局相位噪声的新的滤波方法。
背景技术
光学三维测量方法已广泛应用于工业检测、逆向工程、人体扫描、文物保护、服装鞋帽等多个领域,对自由曲面的检测具有速度快、精度高的优势。按照成像照明方式的不同,光学三维测量技术可分为被动三维测量和主动三维测量两大类。在主动三维测量技术中,结构光三维测量技术发展最为迅速,尤其是相位测量轮廓术(Phase Measuring Profilometry,PMP),也被称为相移测量轮廓术(Phase Shifting Profilometry,PSP),是目前三维测量产品中常用的测量方法。相位测量方法是向被测物体上投射固定周期的按照三角函数(正弦或者余弦)规律变化的光亮度图像,此光亮度图像经过大于3步的均匀相移,最好为4-6步均匀相移,向物体投射4-6次光亮度图像,最终完成一个周期的相位移动。物体上面的每个点,经过相移图像的投射后,在图像中会分别获得几个不同的亮度值。此亮度值经过解相运算,会获得唯一的相位值。由于目前采集到的图像的幅面较大,为了提高相位精度,需要向被测物体投射多个周期的相位图,因此,在一副图像中,相同相位值会出现多次。为了在图像中获得唯一的相位值,格雷码方法是常用的辅助解相方法。目前出现的三维测量产品,普遍采用格雷码加相移的光学投射方法,如德国GOM公司的Atos-I型结构光三维测量系统、德国Steinbichler公司的COMMET系列结构光三维测量系统、德国Breuckmann公司的optoTOP系列结构光三维测量系统、北京天远三维科技有限公司的OKIO-II型三维扫描仪、上海数造科技有限公司的3DSS综合型三维扫描仪、天津世纪动力光电科学仪器有限公司的CPOS三维扫描仪等。由于格雷码的编码方法主要靠图像的二值化来进行编码,因此对于物体表面颜色变化较多的情况,一般需要喷涂显影剂才能实现较好的测量效果。多频率的光学投射方法可以解决无法喷涂显影剂的三维测量难题,因为多频率的方法是通过不同频率的波形叠加进行相位解相,因此不是简单的通过图像二值化的操作,对于物体表面颜色不一致的情况,无需喷涂显影剂,就可以实现较好的测量。
但是无论是基于格雷码还是多频率的三维测量方式,最终的目的都是获得在被测空间全场范围内的唯一的全局相位。由于被测场景光学干扰、物体遮挡等多方面的原因,最终获得的全局相位并不是理想的递增线性函数,而是杂乱的波形。如果不对此相位进行滤波处理,在三维重构过程中,会生成大量的噪声点云,并且影响运算速度,也给用户带来很多的麻烦。本发明设计了一种全局相位噪声滤波方法,可以有效滤除全局相位中得噪声干扰,实现更好的三维测量效果。
发明内容
本发明提供一种全局相位滤波方法,该方法能够应用于高精度三维测量中,该方法可以去除全局相位中的噪声干扰,滤波后的全局相位用于高精度三维测量中,可以获得更好的三维重建效果。
所述的全局相位滤波方法是用于高精度三维测量装置中,其硬件系统包括:
用于投射格雷码或者多种频率光信号的光源投射装置,光源投射装置的分辨率为LR×LC
用于精度控制、图像采集和数据处理的计算机;
用于采集图像的彩色摄像机,图像分辨率为CR×CC,相机个数为1-2个;
用于放置所述的光源和所述的摄像机的扫描平台;
对于分辨率为N×M的相位图像,存在N行相位数据,每行M个像素,每一行相位数据去噪的详细操作步骤如下:
步骤1:对M个数据相邻像素进行斜率运算,可以计算出M-1个斜率,即:K1、K2、.....KM-1,开辟t-1个浮点类型的内存空间,用于存储K1~Kt-1,定义计数变量:n、z,编号变量:i、q,将两个点(k点和j点)形成的斜率编号为q(q的初始值为1)并存储于第i(i的初始值为0)个存储空间中,k点和j点的斜率应按照如下公式计算:
K j = y k - y j x k - x j
步骤2:判定q是否大于等于M-1如果为真,执行步骤5,如果为假则将编号i,q分别加1并将斜率Kq存入第i个存储空间,按照如下判定条件并计数到变量n中:
如果|Ki-Kj|>ε(j=i-1,ε为设定的阈值),则重新将Ki代入第0个存储空间,并且将原先计数值n清0,编号q应执行加i+1的操作,将i的值归零,并重复执行步骤2;
如果|Ki-Kj|≤ε(j=i-1,ε为设定的阈值),变量n计数值应加1,则重复执行步骤2直到n与t-2相等为止,然后执行步骤3;
步骤3:根据已经找出的t个(t-1个斜率)连续光滑点作基准向后按照公式|Ki-Kj|≤ε进行延伸判定,编号q应执行加1操作,而此时Kj指的是刚算出t-1个斜率中的最后一个斜率,Ki在M-1个斜率中应指的是编号为q的斜率;
当判定出|Ki-Kj|>ε时将第q个点的幅值置0并开始向变量z中计数;
如果满足|Ki-Kj|≤ε条件时,计数变量z清0;
当计数值达到人为给定a(a是出现连续噪声的个数,并且a的值不能大于等于t,经测试当a大于等于t时,会出现相当多的噪点)值时,执行步骤4,否则执行步骤3,直到q大于等于M-1为止;
步骤4:将编号q执行减去a-1的操作,重新执行步骤2;
步骤5:删除已开辟的内存空间,去噪结束。
附图说明
图1:被测物体相位图
图2:512行相位图
(a)512行局部相位图
(b)去噪之前512行全局相位图
(c)去噪之后512行全局相位图
图3:相位滤波算法流程图
图4:属于编码M的像素示意图
(a)去除噪声干扰前,属于编码M的像素示意图
(b)去除噪声干扰后,属于编码M的像素示意图
图5:三维重建效果
(a)去除噪声干扰前的三维重建效果
(b)去除噪声干扰后的三维重建效果
具体实施方式
相位测量的原理是是向被测物体投射周期性的光波(正弦或者余弦函数波),此周期函数经过N步相移,完成一个周期的移动,N的值应大于3,最好是4-8步。
以一种频率的光信号为例,某一个像素点(x,y),实际采集到的光亮度Ir的公式如公式(1)所示:
Ir(x,y)=a(x,y)+bsin(θ(x,y)+δ(x,y))                   公式(1)
其中:
Ir(x,y):为(x,y)点采集到光亮度;
a(x,y):为(x,y)点的背景光亮度;
θ(x,y):为(x,y)点的相位值;
b:为(x,y)点投射光的最大亮度值;
δ(x,y):为(x,y)点的正弦相移值;
以四步相移为例,对于某一个像素点(x,y),通过四次光投射,采集到的图像灰度值分别为Irp1(x,y)、Irp2(x,y)、Irp3(x,y)、Irp4(x,y)。
因此相位值θ(x,y)可通过公式(2)获得:
θ ( x , y ) = arctan ( I rp 1 ( x , y ) - I rp 3 ( x , y ) I rp 2 ( x , y ) - I rp 4 ( x , y ) ) 公式(2)
由公式(2)所得到的相位值为(x,y)点在一个正弦周期内的局部相位值,要想得到该点在整个相位图中得全局相位值,需要知道该点所在的相位周期在整个相位图中的编码。假设该点所在周期的编码值为M(x,y),则无论采用格雷码方式,还是多频率的测量方法,最终的全局相位的表达式都可表示为公式(3)所示。
θG(x,y)=θi(x,y)+2π×M(x,y)                     公式(3)
其中:θG(x,y)为(x,y)点的全局相位值;
θi(x,y)为(x,y)点根据相移公式计算得到的局部相位值;
M(x,y)为(x,y)点所在的相位周期在全局相位中的编码值。
以一个模型手的测量为例,选取一副相移图像如图1所示。
由于物体表面颜色、物体表面形状以及背景光亮度的干扰,使得最后的全局相位并不是理想的均匀递增的线性函数变化规律,而是会存在很多的噪声干扰。取图1第512行的局部相位和全局相位,分别如图2(a)和2(b)所示。从图2(b)可以看出,最终的全局相位存在很多偏离直线的噪声干扰点。这些点必须去除才能保证三维测量的准确性和可靠性。
由于理想的合成相位应该是线性递增的函数,但是如果使用简单的斜率滤波方法的话,在不同测量之间,斜率的值可能会有所变化,故而采用相邻斜率差的方法,程序先计算相邻的两两点间直线斜率,然后将这两个斜率相减取绝对值。如果得出的结果在阈值范围内则认为这三个点连成的线是光滑的,反之,认为这三个点中至少有一个是噪点。程序会事先查找一定数量的点,这些点符合上述规则,并且以这些点作为基准进行前后查找,如果查找出前面存在断点,从断点开始的前面所有点将会被置为噪点。而后面如果出现了断点,则从断点的前一点处向后重新查找一定数量符合规则的新点作为基准,向后查找符合规则的点,直到所有点全部搜索完毕为止。至于判断断点的准则经过调试,认为判断噪声点的个数应小于等于作为基准点的数量时,为曲线出现断点,需要重新找新的基准。
去噪的原理是建立在连续条件:和光滑条件:的基础之上。具体执行过程如下:
假设图像的分辨率为N×M(即N行M列,则每行图像含有M个像素),设y1、y2、......、yM为x1、x2、......、xM一一对应的映像,x1、x2、......、xM是连续的自然数,它们都是离散的点。要想按照连续条件和光滑条件两个方法去除噪声,就必须求出两个点(k点和j点)的斜率即:由于待处理的全部是离散点,因此连续条件与光滑条件需要根据相邻两点的斜率值进行更改。首先是光滑条件,假设3个相邻的点k、j和i点,计算相继三个点的两个斜率: 根据光滑条件的规则,对于离散点,光滑条件可定义为|K1-K2|≤ε(ε为设定的光滑阈值)。定义好光滑的条件后,我们将根据光滑条件进行去噪。随着噪点的消去,定义连续就非常重要了。如果不定义好连续,光滑也就无从谈起。由于待处理的数据是图像中每行的数据,由于前面假设图像的分辨率为M×N,则每行存在M个需要去除噪声的相位值,这些相位数据的横坐标X为图像中像素点的横向坐标值,x值的变化是从1-M;Y坐标值为相位值,该值为浮点类型数据,理论上Y值应该是个递增的线性函数,但是因为背景干扰等原因,会有很多杂散噪声附着在Y值上,使得Y值被噪声干扰的部分为散乱的数据点,而正常数据部分仍然为局部的线性递增函数。非噪声点所对应的x值应当是递增(递减),并且增(减)量全部相同即:
xj-xk≡a  (a是常数,j-k≡1)
为了方便起见本说明全令a≡1。
对于分辨率为N×M的相位图像,存在N行相位数据,每行M个像素,每一行相位数据去噪的详细操作步骤如下:
步骤1:对M个数据相邻像素进行斜率运算,可以计算出M-1个斜率,即:K1、K2、.....KM-1,开辟t-1个浮点类型的内存空间,用于存储K1~Kt-1,定义计数变量:n、z,编号变量:i、q,将两个点(k点和j点)形成的斜率编号为q(q的初始值为1)并存储于第i(i的初始值为0)个存储空间中,k点和j点的斜率应按照如下公式计算:
K j = y k - y j x k - x j
步骤2:判定q是否大于等于M-1如果为真,执行步骤5,如果为假则将编号i、q分别加1并将斜率Kq存入第i个存储空间,按照如下判定条件并计数到变量n中:
如果|Ki-Kj|>ε(j=i-1,ε为设定的阈值),则重新将Ki代入第0个存储空间,并且将原先计数值n清0,编号q应执行加i+1的操作,将i的值归零,并重复执行步骤2;
如果|Ki-Kj|≤ε(j=i-1,ε为设定的阈值),变量n计数值应加1,则重复执行步骤2直到n与t-2相等为止,然后执行步骤3;
步骤3:根据已经找出的t个(t-1个斜率)连续光滑点作基准向后按照公式|Ki-Kj|≤ε进行延伸判定,编号q应执行加1操作,而此时Kj指的是刚算出t-1个斜率中的最后一个斜率,Ki在M-1个斜率中应指的是编号为q的斜率;
当判定出|Ki-Kj|>ε时将第q个点的幅值置0并开始向变量z中计数;
如果满足|Ki-Kj|≤ε条件时,计数变量z清0;
当计数值达到人为给定a(a是出现连续噪声的个数,并且a的值不能大于等于t,经测试当a大于等于t时,会出现相当多的噪点)值时,执行步骤4,否则执行步骤3,直到q大于等于M-1为止;
步骤4:将编号q执行减去a-1的操作,重新执行步骤2;
步骤5:删除已开辟的内存空间,去噪结束。
经过噪声滤除后的全局相位图如图2(c)所示。
噪声滤除算法示意图如图3所示。
取编码M=26的一个周期的相位在灰度图中的对应关系,在去除噪声之前,如图4(a)所示,从图中可以看到除了本应该属于编码M的像素之外,很多毛糙的噪声也含在里面。经过相位滤波后,编码M的噪声全部去除了,如图4(b)所示。
使用本发明所设计的滤波方法之前,三维重建的效果很差,如图5(a)所示。经过本发明的滤波方法之后,三维重建效果很光滑,如图5(b)所示。
本发明与已有的噪声滤波方法的最大区别是:
(1)本发明是针对最终的全局相位进行噪声滤波,而以往的滤波方法,都是注重前端的噪声去除方法,使得最终的相位精度无法得到有效控制。
(2)本发明所设计的全局相位滤波方法,具有更广泛的适用性,不但可以适用于格雷码的编码方式,也可以适用于多频率的编码方式。对于其他类似信号的噪声干扰难题,也有适应性。
综上所述,本发明所述的相位滤波方法的优点是:
(1)由于本发明所述的相位滤波方法,是建立在全局相位的基础之上,因此比已有的局部相位滤波方式具有更好的处理效果,有效消除了最后环节出现的噪声干扰。
(2)已有的相位滤波方法,在去除噪声干扰的同时,很多有效的相位信息也被去除,而本发明所介绍的相位滤波方法,在去除了所有噪声信号的同时,很好的保留了有效的相位信息,使得三维重建的效果更好的包含了边缘等细节信息。
(3)使用本发明的相位滤波方法,可以避免杂散噪声在三维重建时的数据干扰和误判断,得到平滑的三维重建数据,减少用户的人工干预工作量,提高工作效率。
(4)由于本发明所介绍的相位滤波方法,是一种适应性广泛的去除噪声干扰的方法,对于目前三维重建领域中所使用的格雷码方法、多频率方法等具有较好的适用性。除此之外,在其余信号的噪声干扰问题上,本发明所使用的方法也有较好的适用性。
以上示意性的对本发明及其实施方式进行了描述,该描述没有局限性,附图中所示的也只是本发明的实施方式之一。所以,如果本领域的普通技术人员受其启示,在不脱离本发明创造宗旨的情况下,采用其它形式的同类部件或其它形式的各部件布局方式,不经创造性的设计出与该技术方案相似的技术方案与实施例,均应属于本发明的保护范围。

Claims (1)

1.一种应用于三维测量的全局相位滤波方法,其特征是,包括下列步骤:
步骤1:对于分辨率为N×M的相位图像,存在N行相位数据,每行M个像素,对M个数据相邻像素进行斜率运算,可以计算出M-1个斜率,即:K1、K2、......KM-1,开辟t-1个浮点类型的内存空间,用于存储K1~Kt-1,定义计数变量:n、z,编号变量:i、q,将k点和j点两个点形成的斜率编号为q,q的初始值为1,并存储于第i个存储空间,i的初始值为0,k点和j点的斜率应按照如下公式计算:
K j = y k - y j x k - x j
步骤2:判定q是否大于等于M-1如果为真,执行步骤5,如果为假则将编号i、q分别加1并将斜率Kq存入第i个存储空间,按照如下判定条件并计数到变量n中:
如果|Ki-Kj|>ε,j=i-1,ε为设定的阈值,则重新将Ki代入第0个存储空间,并且将原先计数值n清0,编号q应执行加i+1的操作,将i的值归零,并重复执行步骤2;
如果|Ki-Kj|≤ε,j=i-1,ε为设定的阈值,变量n计数值应加1,则重复执行步骤2直到n与t-2相等为止,然后执行步骤3;
步骤3:根据已经找出的t个(t-1个斜率)连续光滑点作基准向后按照公式|Ki-Kj|≤ε进行延伸判定,编号q应执行加1操作,而此时Kj指的是刚算出t-1个斜率中的最后一个斜率,Ki在M-1个斜率中应指的是编号为q的斜率;
当判定出|Ki-Kj|>ε时将第q个点的幅值置0并开始向变量z中计数;
如果满足|Ki-Kj|≤ε条件时,计数变量z清0;
当计数值z达到人为给定a值时,a是出现连续噪声的个数,并且a的值不能大于等于t,执行步骤4,否则执行步骤3,直到q大于等于M-1为止;
步骤4:将编号q执行减去a-1的操作,重新执行步骤2;
步骤5:删除已开辟的内存空间,去噪结束。
CN201210063823.XA 2012-03-13 2012-03-13 一种应用于三维测量的全局相位滤波方法 Expired - Fee Related CN102607460B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201210063823.XA CN102607460B (zh) 2012-03-13 2012-03-13 一种应用于三维测量的全局相位滤波方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201210063823.XA CN102607460B (zh) 2012-03-13 2012-03-13 一种应用于三维测量的全局相位滤波方法

Publications (2)

Publication Number Publication Date
CN102607460A CN102607460A (zh) 2012-07-25
CN102607460B true CN102607460B (zh) 2014-07-30

Family

ID=46525067

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201210063823.XA Expired - Fee Related CN102607460B (zh) 2012-03-13 2012-03-13 一种应用于三维测量的全局相位滤波方法

Country Status (1)

Country Link
CN (1) CN102607460B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1628410A (zh) * 2002-05-28 2005-06-15 Rf信息公司 二次奈奎斯特斜率滤波器
CN1921557A (zh) * 2005-08-26 2007-02-28 株式会社日立制作所 图像信号处理装置及图像信号处理方法
GB2457215A (en) * 2007-03-07 2009-08-12 Nikolaos Kokkas Automatic 3D Modelling
CN201361622Y (zh) * 2009-01-19 2009-12-16 中国计量学院 基于斜率滤波算法的快速重量分选系统
CN102222325A (zh) * 2011-06-17 2011-10-19 电子科技大学 一种基于相位特性的物体轮廓检测方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1628410A (zh) * 2002-05-28 2005-06-15 Rf信息公司 二次奈奎斯特斜率滤波器
CN1921557A (zh) * 2005-08-26 2007-02-28 株式会社日立制作所 图像信号处理装置及图像信号处理方法
GB2457215A (en) * 2007-03-07 2009-08-12 Nikolaos Kokkas Automatic 3D Modelling
CN201361622Y (zh) * 2009-01-19 2009-12-16 中国计量学院 基于斜率滤波算法的快速重量分选系统
CN102222325A (zh) * 2011-06-17 2011-10-19 电子科技大学 一种基于相位特性的物体轮廓检测方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
刘淑娟等.基于斜率的自适应中值滤波算法.《计算机应用》.2012,第32卷(第3期),全文.
基于斜率的自适应中值滤波算法;刘淑娟等;《计算机应用》;20120301;第32卷(第3期);全文 *
洪微.限斜率滤波.《微型计算机与应用系统设计》.1986,正文第20页第8.4.2节. *

Also Published As

Publication number Publication date
CN102607460A (zh) 2012-07-25

Similar Documents

Publication Publication Date Title
Feng et al. Robust dynamic 3-D measurements with motion-compensated phase-shifting profilometry
CN101881605B (zh) 基于相位编码技术的光学三维测量方法
CN102589478B (zh) 一种应用于多频率三维测量的全局相位解相方法
CN109341589B (zh) 一种光栅图像投影方法、三维重建方法及三维重建系统
CN110174079B (zh) 一种基于四步相移编码型面结构光的三维重建方法
CN108519064B (zh) 一种应用于多频率三维测量的反光抑制方法
CN111207692B (zh) 一种改进的分段阶梯相位编码三维测量方法
CN101975558A (zh) 基于彩色光栅投影的快速三维测量方法
CN102519390A (zh) 三编码周期灰度梯形相移结构光三维信息获取方法
CN102853783A (zh) 一种高精度多波长三维测量方法
CN111207694B (zh) 一种双步相移法结合相位编码的三维测量方法
CN104482877A (zh) 动态物体三维成像中的运动补偿方法与系统
CN103267496A (zh) 一种基于小波变换的改进窗口傅里叶三维测量法
Liu et al. Background and amplitude encoded fringe patterns for 3D surface-shape measurement
CN109631798A (zh) 一种基于π相移方法的三维面形垂直测量方法
CN112504165A (zh) 一种基于双边滤波优化的复合立体相位展开方法
CN102881043A (zh) 一种应用于明暗反差大场景的高精度三维重建方法
Rodriguez-Vera et al. Phase locked loop profilometry
CN106840039B (zh) 一种基于复合光栅投影的快速三维测量方法
CN108534715A (zh) 一种应用于不同波长三维测量的展开相位解相方法
CN102607460B (zh) 一种应用于三维测量的全局相位滤波方法
CN102519396B (zh) 三个灰度对称线性编码周期的采样点三维信息获取方法
CN107421467A (zh) 结构光三维测量方法及系统
CN112212806A (zh) 一种基于相位信息导向的立体相位展开方法
CN114252026B (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
C14 Grant of patent or utility model
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20140730

Termination date: 20150313

EXPY Termination of patent right or utility model