CN103679652A - 一种大幅提高成像质量的图像复原系统 - Google Patents

一种大幅提高成像质量的图像复原系统 Download PDF

Info

Publication number
CN103679652A
CN103679652A CN201310629499.8A CN201310629499A CN103679652A CN 103679652 A CN103679652 A CN 103679652A CN 201310629499 A CN201310629499 A CN 201310629499A CN 103679652 A CN103679652 A CN 103679652A
Authority
CN
China
Prior art keywords
image
module
data
row
mtf
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.)
Granted
Application number
CN201310629499.8A
Other languages
English (en)
Other versions
CN103679652B (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.)
Beijing Institute of Space Research Mechanical and Electricity
Original Assignee
Beijing Institute of Space Research Mechanical and Electricity
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 Beijing Institute of Space Research Mechanical and Electricity filed Critical Beijing Institute of Space Research Mechanical and Electricity
Priority to CN201310629499.8A priority Critical patent/CN103679652B/zh
Publication of CN103679652A publication Critical patent/CN103679652A/zh
Application granted granted Critical
Publication of CN103679652B publication Critical patent/CN103679652B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Image Processing (AREA)
  • Facsimile Image Signal Circuits (AREA)

Abstract

本发明一种大幅提高成像质量的图像复原系统,包括:图像数据存储模块、刃边法计算边缘扩展函数模块、线扩展函数计算模块、归一化MTF函数计算模块、MTFC数字滤波器系数计算存储模块、图像数据边界处理模块、行向图像数字滤波模块、行向滤波图像抑噪模块、行向滤波图像存储模块、列向图像数字滤波模块、列向滤波图像抑噪模块、复原后图像存储模块。本发明达到了在轨实时处理遥感图像,增强了图像清晰度并有效抑制噪声。

Description

一种大幅提高成像质量的图像复原系统
技术领域
本发明涉及一种大幅提高成像质量的图像复原系统,属于图像处理技术领域。
背景技术
目前遥感图像处理大部分是在地面完成,但是由于受传输带宽的限制,遥感图像从卫星传送回地面前一般要进行3:1--8:1的有损压缩,地面接收后再解压缩得到的图像存在一定程度的失真,这在一定程度上影响后续的图像增强处理。所以,需在图像压缩前针对原始的14bit全色图像数据进行数字滤波等一系列图像处理,这样相对于图像传回地面再处理减少了有损压缩带来的像质损失,而且实时性高,图像在传回地面前便已处理完成,有助于提升全色图像零级图像像质。同时,在地面设备对图像的处理中,进行调制传递函数补偿(MTFC),可使成像系统MTF得到提升,具有良好的成像性能。美国研制的1m以下的分辨率的商业卫星Ikonos-2、QuickBird-2、OrbView-3等都成功地应用了地面MTFC技术,使像质得到显著的提高。其中Ikonos-2在轨测量的MTF为0.03~0.05,经过MTFC之后能够提高到0.17;OrbView-3在轨测量的MTF为0.1±0.01,经过MTFC之后能提高到0.15±0.01。
对退化图像进行复原最常用的方法是有约束复原方法,包括约束最小二乘方滤波、维纳滤波等。该方法通常已知各退化过程的精确MTF曲线进行拟合并构造MTF矩阵。
综上所述,上述方法在实际应用中存在不足之处,以下重要因素制约了采用以上方法的有效性:
(1)图像在星上做有损压缩处理会部分损失图像信息,图像在地面需针对此问题做图像增强处理。
(2)地面MTFC技术可以在图像传回地面后做处理,再给客户使用,使用步骤繁琐。
(3)星上MTFC补偿在频率域进行,占用资源以及处理速度。对于超过1Gbits的图像进行傅立叶变换、反变换和变换域乘积处理,其运算量非常大,会严重制约系统的处理速度,实时性较差。
(4)MTF参数无法完整、精确地获得。且大气条件因时、因地而改变导致其MTF也不确定。各个退化过程中的误差传递造成理论求解与实际情况差异较大。
(5)常规图像复原的方法,对图像增强有益,但增强图像细节同时放大噪声。
发明内容
本发明技术解决问题:克服现有技术实时性差,系统运算量大,复原精度有限的问题,提供一种大幅提高成像质量的图像复原系统,达到在轨实时处理遥感图像,增强了图像清晰度并有效抑制噪声。
本发明的技术解决方案:一种大幅提高成像质量的图像复原系统,包括:图像数据存储模块、刃边法计算边缘扩展函数模块、线扩展函数计算模块、归一化MTF函数计算模块、MTFC数字滤波器系数计算存储模块、图像数据边界处理模块、行向图像数字滤波模块、行向滤波图像抑噪模块、行向滤波图像存储模块、列向图像数字滤波模块、列向滤波图像抑噪模块、复原后图像存储模块:首先,图像数据存储模块接收、存储遥感相机CCD焦面电路输出的原始图像数据;再将原始图像数据输出至刃边法计算边缘扩展函数模块和图像数据边界处理模块;刃边法边缘扩展函数模块使用原始图像数据计算得到相机边缘扩展函数(LSF);线扩展函数计算模块使用LSF计算得到单行图像线扩展函数(ESF);归一化MTF函数计算模块使用ESF计算得到相机系统MTF序列曲线函数。MTFC数字滤波器系数计算存储模块使用MTF计算、存储并输出行向、列向数字滤波器系数。图像数据边界处理模块对原始图像数据做边界处理;行向图像数字滤波模块使用行向数字滤波器系数对边界处理后数据滤波;行向滤波图像抑噪模块对行向滤波后数据做抑噪处理;行向滤波图像存储模块存储并输出行向滤波抑噪后图像数据;列向图像数字滤波模块对行向滤波抑噪后图像数据做列向滤波处理;列向滤波图像抑噪模块对列向滤波后数据做抑噪处理;列向滤波抑噪处理后数据存入复原后图像存储模块。
所述刃边法计算边缘扩展函数模块实现如下:
(1)选取含有地面边缘纹理(冰水分隔线、建筑物边缘等)的某遥感图像做边缘检测。
(2)抠取具有刀刃形状的目标子图,保证刀刃边缘两侧的灰度分布较均匀无其他因素干扰。
(3)根据边缘成像的灰度分布拟合出边缘扩展函数曲线。
:所述线扩展函数计算模块实现如下:
(1)对边缘扩展函数曲线序列数据做简单差分,即计算x(n)-x(n-1),并取max(x(n)-x(n-1)),即可检测得到最大斜率位置。
(2)对图像的每一行数据点进行三次样条插值,通常选取插值分辨率为0.05,即在每两个数据点之间插入20个数据点,插值结果为一条灰度曲线,即为该行的边缘扩展函数ESF。
(3)对各行重排数据得到的ESF求平均,以最小化误差及降低噪声。ESF曲线与原始灰度分布对比如图2所示。
(4)对ESF做差分运算LSF(n)=ESF(n)-ESF(n-1),得到线扩展函数LSF曲线。
所述归一化MTF函数计算模块实现如下:
(1)对线扩展函数曲线LSF做离散傅里叶变换后取各分量的模MTF(n)=|DFT(LSF(n))|。
(2)以变换后的直流分量即第一个MTF值为标准做归一化处理norm_MTF(n)=MTF(n)/MTF(1),即得到MTF曲线序列如图4所示。
(3)使用归一化频率轴,Nyquist频率为截止频率的一半
Nyquist_frequency=(whole_data_size×resolution)/2+1=(Number_of_trimmed_pixel)/2+1
其中:Nyquist_frequency为Nyquist频率点,whole_data_size为包含所有插值点的截取后LSF亚像素级宽度,resolution为插值分辨率,Number_of_trimmed_pixel为截取后的LSF像素级宽度。
所述MTFC数字滤波器系数计算存储模块实现如下:
(1)假定成像系统为线性移不变系统,通用图像退化模型来表示:f(x,y)*h(x,y)+n(x,y)=g(x,y)。在这个模型中,一个作用在输入图像f(x,y)上的系统h(x,y)与加性噪声n(x,y)的联合作用导致产生退化图像g(x,y)。
(2)根据最小二乘原理,采用被系统h(x,y)模糊的图像与观测到的图像g(x,y)之间的差在均方意义最小的准则下来估计结果;
f ^ ( x , y ) = F - 1 [ F ^ ( u , v ) ] = F - 1 [ G ( u , v ) H ( u , v ) ]
在频域下,G(u,v)、F(u,v)、H(u,v)分别是g(x,y)、f(x,y)、h(x,y)的傅立叶变换,系统H(u,v)作为滤波函数与F(u,v)的乘积是退化图像g(x,y)的傅里叶变换;
(3)把系统函数H(u,v)看作复数,对幅值作归一化。将H(u,v)替换为k·MTF(u,v),k为系统函数归一化系数。重写上式为:
f ^ ( x , y ) = F - 1 [ 1 k · MTF ( u , v ) G ( u , v ) ]
(4)为了达到提升MTF的效果,对上式进行修改为
F ^ ( u , v ) = G ( u , v ) H ( u , v ) = [ 1 r ′ ( u , v ) k · MTF ( u , v ) ] G ( u , v )
其中,r′(u,v)是实现在频域下不同频率段的补偿,从而有针对性提升MTF。逆滤波复原局限于中低频区域进行,且适用于信噪比较高的情况;
(5)再令 r ( u , v ) = 1 r ′ ( u , v ) k · MTF ( u , v ) , F ^ ( u , v ) = G ( u , v ) [ r ( u , v ) ] α , 其中α为补偿因子,恒大于等于0,此时MTFC补偿函数为: MTFC ( u , v ) = [ r ( u , v ) ] α = [ 1 r ′ ( u , v ) k · MTF ( u , v ) ] α . 通过调整r′及α可以实现对MTF的补偿;
(6)选择重点频段对r′进行补偿和运用不同的补偿因子进行优化。r′选用凸曲线比凹曲线或直线对图像的补偿效果更好。α进一步优化补偿效果,在图像清晰度和噪声之间寻求平衡,信噪比、熵的增加和均方梯度相对最大能达到补偿效果和控制噪声的目的;
(7)计算输出MTFC数字滤波器系数后存入PROM。
所述图像数据边界处理模块实现如下:
(1)图像数据按次序输入,假设一行CCD像元为8190个,图像共有1024行数据,使用D1,D2…,D8190分别代表第1、2至第8190个像元的输出值。每个像元输入时,带有图像数据有效使能;
(2)使用图像数据效使能激发计数器1,输出结果n1
(3)当1≤n1≤5或8186≤n1≤8190时,使用此判断有效使能激发计数器2,输出结果n2
(4)将n2作为选择器输入,datain表示为选择器输入数据,dataout表示选择器输出数据,条件对应关系如下:
if n2=0dataout=D5
if n2=1dataout=D4
if n2=2dataout=D3
if n2=3dataout=D2
if n2=4dataout=D1
if n2=5dataout=datain
if n2=6dataout=D8190
if n2=7dataout=D8189
if n2=8dataout=D8188
if n2=9dataout=D8187
if n2=10dataout=D8186
最后得到加边界的图像数据序列为:
D5,D4,D3,D2,D1,D1,D2,D3,…,D8190,D8190,D8189,D8180,D8187,D8186。在原始图像线阵数据一行的开始和结束分别填补5个像元数据,为这行数据开始5个和结束5个像元数据作对称处理。
所述行向图像数字滤波模块实现如下:
MTFC数字滤波算法图像水平(线阵)方向和垂直(锥扫)方向的复原分开进行,先水平(线阵)方向上的卷积、抑噪,再垂直方向上的卷积、抑噪;
水平(线阵)方向卷积系数为11个。每个卷积核系数的范围为[-8,8],16bit存储在Prom中,第1位为符号位,第2-4位为整数位,第5-16位为小数位;
算法为:
dataout=k1·Dn-5+k2·Dn-4+k3·Dn-3+k4·Dn-2+k5·Dn-1+k6·Dn+k7·Dn+1+k8·Dn+2+k9·Dn+3+k10·Dn+4+k11·Dn+5
其中,n=1,2,…8190。k1,k2,…k11为行向滤波器系数。由原始图像数据边界处理模块计算,D-5,D-4,D-3,D-2,D-1分别对应图像数据D5,D4,D3,D2,D1。D8191,D8192,D8193,D8194,D8195分别对应图像数据D8190,D8189,D8188,D8187,D8186
所述行向滤波图像抑噪模块实现如下:
(1)假设某一像元原始灰度值为i0,卷积运算后结果为i1,取i2=i0+k(i1-i0)为最终复原值。当k=1时,则复原后MTF为全额补偿,若k=0.5,则MTF为半额补偿,若k=0,则为零补偿;
(2)取Δ=K||i1-i0||,其中K为抑噪阈值参数。
(3)当Δ≥i0时,k=1,即i2=i1
(4)当Δ<i0时,k=Δ/2n,即i2=i0+Δ(i1-i0)/2n
n为像元量化位数。其中此处除法可用移位代替
最终值复原值为i2
k为浮点数,具体实现时可以不计算k,而直接计算i2,从而避免浮点数的运算。
所述行向滤波图像存储模块实现如下:行向滤波处理和抑噪处理后同时十列数据图像数据缓存到RAM中,以进行列方向数据处理。
本发明与现有技术相比,有如下优点:
(1)本发明无有损压缩处理,减少损失图像信息。
(2)本发明使用MTFC图像滤波抑噪算法模块,在轨实时对原始图像数据做复原处理,对地面顾客以及星上其他数据系统提供直观更清晰,细节更丰富图像。
(3)本发明在空间域做MTFC滤波补偿处理克服了频域处理速度低的问题,可实现MTFC的在轨计算机的实时快速运算。
(4)本发明使用MTF在轨测量计算模块在轨测量相机系统退化过程,可针对MTF退化曲线特征,在保证原始图像信噪比未受较大影响的条件下,增强了图像灰度梯度、点锐度、边缘对比度,对于MTF的中高频域增强。
(5)本发明添加抑噪算法的MTFC图像滤波抑噪算法模块,使复原后图像信噪比与原始图像信噪比相比降低很小,硬件实现结果与理论计算结果相差<5%,开启后信噪比下降<1dB。
附图说明
图1为本发明大幅提高成像质量的图像复原系统组成模块结构图;
图2为ESF曲线与原始灰度分布对比;
图3为LSF曲线;
图4为MTF曲线(频率标准化后);
图5为图像退化模型;
图6为PROM参数存储格式;
图7为图像数据边界处理模块流程图;
图8为行(列)向滤波图像抑噪模块SG模块设计;
图9为滤波器幅频响应曲线;
图10为相机分系统全链路的MTF曲线;
图11为原始图像及使用抑噪方案复原后图像,其中a为仿真模糊图像(SNR=39dB,b为复原后图像。
具体实施方式
如图1所示,本发明系统组成为:图像数据存储模块1、刃边法计算边缘扩展函数模块2、线扩展函数计算模块3、归一化MTF函数计算模块4、MTFC数字滤波器系数计算存储模块5、原始图像数据边界处理模块6、行向图像数字滤波模块7、行向滤波图像抑噪模块8、行向滤波图像存储模块9、列向图像数字滤波模块10、列向滤波图像抑噪模块11、复原图像存储模块12。
1、图像数据存储模块。
CCD输出焦面模拟信号后,经AD转换数字量化,缓存入FPGA中BlockRam并输出图像数字信号。
2、刃边法计算边缘扩展函数模块。
(1)选取含有地面边缘纹理(冰水分隔线、建筑物边缘等)的某遥感图像做边缘检测。
(2)抠取具有刀刃形状的目标子图,保证刀刃边缘两侧的灰度分布较均匀无其他因素干扰。
(3)根据边缘成像的灰度分布拟合出边缘扩展函数曲线。图像中某行灰度值(DN)分布如下:
表1灰度分布序列
Figure BDA0000426795670000101
3、线扩展函数计算模块。
对边缘扩展函数曲线一次求导,得出线扩展函数曲线。
(1)对边缘扩展函数曲线序列数据做简单差分,即计算x(n)-x(n-1),并取max(x(n)-x(n-1)),即可检测得到最大斜率位置。理论上该最大斜率点(拐点)即为像素级边缘点。
(2)刃边法假设图像边缘位于同一直线,对图像逐行确定的亚像元位置进行最小二乘拟合,以获得亚像素级边缘中心线位置。具体做法为,对图像的每一行数据点进行三次样条插值,通常选取插值分辨率为0.05,即在每两个数据点之间插入20个数据点,插值结果为一条灰度曲线,即为该行的边缘扩展函数ESF。
(3)对各行重排数据得到的ESF求平均,以最小化误差及降低噪声。ESF曲线与原始灰度分布对比如图2所示。
(4)对ESF做差分运算LSF(n)=ESF(n)-ESF(n-1),得到线扩展函数LSF曲线如图3所示。
4、归一化MTF函数计算模块。
(1)对线扩展函数曲线LSF做离散傅里叶变换后取各分量的模MTF(n)=|DFT(LSF(n))|。
(2)以变换后的直流分量即第一个MTF值为标准做归一化处理norm_MTF(n)=MTF(n)MTF(1),即得到MTF曲线序列如图4所示。
(3)使用归一化频率轴,Nyquist频率为截止频率的一半
Nyquist_frequency=(whole_data_size×resolution)/2+1=(Number_of_trimmed_pixel)/2+1
其中:Nyquist_frequency为Nyquist频率点,whole_data_size为包含所有插值点的截取后LSF亚像素级宽度,resolution为插值分辨率,Number_of_trimmed_pixel为截取后的LSF像素级宽度。
5、MTFC数字滤波器系数计算存储模块。
从成像过程来看,图像复原是针对引起图像退化的原因,建立退化模型,并根据模型得到退化前的图像。采用MTFC复原遥感图像,是在频域下对退化模型进行补偿,得到近似的退化前的图像。
(1)假定成像系统为线性移不变系统,且噪声是加性的,在线性代数范畴内,成像系统可以用一个通用图像退化模型来表示,如图5所示:f(x,y)*h(x,y)+n(x,y)=g(x,y)。
在这个模型中,一个作用在输入图像f(x,y)上的系统h(x,y)与加性噪声n(x,y)的联合作用导致产生退化图像g(x,y)。根据这个模型恢复图像,即在给定g(x,y)和h(x,y)的基础上得到对f(x,y)的某个近似的过程(这里假设已知n(x,y)的统计特性)。
(2)根据最小二乘原理,采用被系统h(x,y)模糊的图像与观测到的图像g(x,y)之间的差在均方意义最小的准则下来估计结果。
f ^ ( x , y ) = F - 1 [ F ^ ( u , v ) ] = F - 1 [ G ( u , v ) H ( u , v ) ]
在频域下,G(u,v)、F(u,v)、H(u,v)分别是g(x,y)、f(x,y)、h(x,y)的傅立叶变换,系统H(u,v)作为滤波函数与F(u,v)的乘积是退化图像g(x,y)的傅里叶变换。将
Figure BDA0000426795670000121
反傅里叶变换是一个逆滤波过程,可得到复原图像。
(3)把系统函数H(u,v)看作复数,表示为
Figure BDA0000426795670000122
中,|H|是幅值,
Figure BDA0000426795670000123
为相位。对幅值作归一化,使其零频率的幅值为1,称此归一化的幅值为调制传递函数(Modulation Transfer Function,MTF)。将H(u,v)替换为k·MTF(u,v),k为系统函数归一化系数。重写上式为:
f ^ ( x , y ) = F - 1 [ 1 k &CenterDot; MTF ( u , v ) G ( u , v ) ]
(4)为了达到提升MTF的效果,对上式进行修改为
F ^ ( u , v ) = G ( u , v ) H ( u , v ) = [ 1 r &prime; ( u , v ) k &CenterDot; MTF ( u , v ) ] G ( u , v )
其中,r′(u,v)是实现在频域下不同频率段的补偿,从而有针对性提升MTF。实际H(u,v)随频率增大快速衰减,在低频段就衰减很多,而噪声则大部分在高频,衰减速度较慢。为避免H(u,v)值太小,导致结果出现不定值,因此逆滤波复原常局限于中低频区域进行,且适用于信噪比较高的情况。
(5)再令 r ( u , v ) = 1 r &prime; ( u , v ) k &CenterDot; MTF ( u , v ) . 为进一步提升MTF,理论上r(u,v)应该恒大于1,但为了适应不同遥感图像MTFC复原,需对r(u,v)进行不同程度的衰减或放大。可使用指数方式,即:
Figure BDA0000426795670000127
其中α为补偿因子,恒大于等于0,此时MTFC补偿函数为: MTFC ( u , v ) = [ r ( u , v ) ] &alpha; = [ 1 r &prime; ( u , v ) k &CenterDot; MTF ( u , v ) ] &alpha; . 通过调整r′及α可以实现对MTF的补偿。
(6)选择重点频段对r′进行补偿和运用不同的补偿因子进行优化。r′选用凸曲线比凹曲线或直线对图像的补偿效果更好。α可以进一步优化补偿效果,在图像清晰度和噪声之间寻求平衡。信噪比、熵的增加和均方梯度相对最大能达到补偿效果和控制噪声的目的。
(7)计算输出MTFC数字滤波器系数后存入PROM。PROM中参数存储格式如图6所示。
MTFC数字滤波算法图像水平(线阵)方向和垂直(锥扫)方向的复原分开进行,先水平(线阵)方向上的卷积、抑噪,再垂直方向上的卷积、抑噪。
经实验分析,选用卷积核大小为N=11时,滤波器频率形状丰富,且在轨计算机的运算量在硬件限制内。水平(线阵)方向和垂直(锥扫)方向上卷积系数各11个。每个卷积核系数的范围为[-8,8],16bit存储,第1位为符号位,第2-4位为整数位,第5-16位为小数位。
6、原始图像数据边界处理模块
原始图像数据边界处理模块首先对图像数据进行边界处理。
原始图像数据边界处理模块SG模块设计见图7,步骤如下:
(1)图像数据按次序输入,假设一行CCD像元为8190个,图像共有1024行数据。使用D1,D2…,D8190分别代表第1、2至第8190个像元的输出值。每个像元输入时,带有图像数据有效使能。
(2)使用图像数据有效使能激发计数器1,输出结果n1
(3)当1≤n1≤5或8186≤n1≤8190时,使用此判断有效使能激发计数器2,输出结果n2
(4)将n2作为选择器输入,datain表示为选择器输入数据,dataout表示选择器输出数据,条件对应关系如下:
if n2=0dataout=D5
if n2=1dataout=D4
if n2=2dataout=D3
if n2=3dataout=D2
if n2=4dataout=D1
if n2=5dataout=datain
if n2=6dataout=D8190
if n2=7dataout=D8189
if n2=8dataout=D8188
if n2=9dataout=D8187
if n2=10dataout=D8186
最后得到加边界的图像数据序列为:
D5,D4,D3,D2,D1,D1,D2,D3,…,D8190,D8190,D8189,D8180,D8187,D8186。在原始图像线阵数据一行的开始和结束分别填补5个像元数据,为这行数据开始5个和结束5个像元数据作对称处理。
7、行向图像数字滤波模块。
行向图像数字滤波模块接收来自MTFC数字滤波器系数计算存储模块获得的行卷积参数对图像数据进行行方向卷积运算并输出给行向滤波图像抑噪模块。
算法为:
dataout=k1·Dn-5+k2·Dn-4+k3·Dn-3+k4·Dn-2+k5·Dn-1+k6·Dn+k7·Dn+1+k8·Dn+2+k9·Dn+3+k10·Dn+4+k11·Dn+5
其中,n=1,2,…8190。k1,k2,…k11为行向滤波器系数。由原始图像数据边界处理模块计算,D-5,D-4,D-3,D-2,D-1分别对应图像数据D5,D4,D3,D2,D1。D8191,D8192,D8193,D8194,D8195分别对应图像数据D8190,D8189,D8188,D8187,D8186。
8、行(列)向滤波图像抑噪模块。
图像经过单纯高通滤波器,卷积后图像信噪比严重下降,需要进行抑噪处理。抑噪思路是通过对补偿像元邻近区域图像的高频信息丰富程度进行判断结果来选择MTF是否全额补偿。对于判断结果高于抑噪阈值的像元(即该像元附近高频信息较丰富),MTF全额补偿,对判断结果低于抑噪阈值的像元(即该像元附近高频信息较少),MTF不全额补偿,只进行适当补偿。
另外,当不进行抑噪时,卷积后图像信噪比下降严重,需要进行抑噪处理并需确定一种计算量小而抑噪效果好并对图像锐化程度影响小的算法。一种抑噪思路是通过对补偿像元邻近区域图像的高频信息丰富程度进行判断结果来选择MTF是否全额补偿。对于判断结果高于抑噪阈值的像元(即该像元附近高频信息较丰富),MTF全额补偿,对判断结果低于抑噪阈值的像元(即该像元附近高频信息较少),MTF不全额补偿,只进行适当补偿。
根据上述思路给出以下抑噪方案:
(1)假设某一像元原始灰度值为i0,卷积运算后结果为i1,取i2=i0+k(i1-i0)为最终复原值。k为补偿系数,当k=1时复原后MTF为全额补偿;若k=0.5,MTF为半额补偿;若k=0,则为零补偿。
(2)取Δ=K||i1-i0||,其中K为抑噪阈值参数。
(3)当Δ≥i0时,k=1,即i2=i1
(4)当Δ<i0时,k=Δ/2n,即i2=i0+Δ(i1-i0)/2n
n为像元量化位数。其中此处除法可用移位代替。最终值复原值为i2。k为浮点数,具体实现时不计算k,而直接计算i2,从而避免浮点数的运算。行(列)向滤波图像抑噪模块SG模块设计见图8。完成抑噪方案(1)、(2)、(3)、(4)步骤。
行(列)向滤波图像抑噪模块接收来自行(列)向图像数字滤波模块计算后的图像数据,进行抑噪处理,输出处理后图像数据。抑噪阈值参数K存储在Prom中为正整数,范围为[1,1024],用16bit存储。
9、行向滤波图像存储模块。
原始图像数据为串行数据,传完线阵CCD一行数据后传下一行。列向图像数字滤波操作需要同时十列数据做运算。行向滤波图像存储模块主要将行向滤波处理和抑噪处理后图像数据存储到RAM中,以进行列方向数据处理。本模块使用10个blockram缓存十行图像数据。
10、列向图像数字滤波模块。
列向图像数字滤波模块接收来自MTFC数字滤波器系数计算存储模块获得的列向滤波参数对行向滤波图像存储模块输出图像数据进行列方向滤波运算并输出至列向滤波图像抑噪模块。
列向图像数字滤波模块实现的算法为:
dataout=k12·Dl+k13·Dl+1+k14·Dl+2+k15·Dl+3+k16·Dl+4+k17·Dl+5+k18·Dl+6+k19·Dl+7+k20·Dl+8+k21·Dl+9+k22·Dl+10
其中,l=1,2,…1024,代表图像第l行。k12,k13,…k21为列向滤波器系数。
SG模型功能测试使用matlab的M文件实现,Corr_MTFC.m程序可对MTFC算法关状态时存储的原始图像在matlab环境下做MTFC运算结果,与硬件MTFC算法打开后复原图像作比较。可确定硬件计算的精度。
2、图像复原结果分析
使用刃边法计算型号1遥感器图像MTF结果如图4所示。在Nyquist频率处的MTF值计算结果为0.098。
选取11阶数字滤波器系数为:
[0.00342  -0.00024  -0.03101  -0.09399  -0.160161.5625  -0.16016  -0.09399  -0.031006  -0.00024  0.00342]
MTFC补偿滤波器幅频响应曲线如图9所示。
经补偿后,相机分系统全链路的MTF曲线形状提升如图10所示。
外景成像某全色图像中某块较均匀区域图像进行仿真,不抑噪时复原前后图像信噪比由43.976dB下降为36.976dB,下降了7.000dB。抑噪后,复原后信噪比为43.974dB,信噪比下降0.002dB。
图像复原结果如图11所示。原始图像及其复原后图像客观评价如表2所示。
表2原始图像及其复原后图像客观评价
灰度梯度 点锐度 边缘能量 方差
原始图像 162.02 68.01 8.48 43.03
复原后图像 374.66 100.47 29.23 47.83
本发明未详细阐述部分属于本领域公知技术。

Claims (9)

1.一种大幅提高成像质量的图像复原系统,其特征在于包括:图像数据存储模块、刃边法计算边缘扩展函数模块、线扩展函数计算模块、归一化MTF函数计算模块、MTFC数字滤波器系数计算存储模块、图像数据边界处理模块、行向图像数字滤波模块、行向滤波图像抑噪模块、行向滤波图像存储模块、列向图像数字滤波模块、列向滤波图像抑噪模块、复原后图像存储模块:首先,图像数据存储模块接收、存储遥感相机CCD焦面电路输出的原始图像数据;再将原始图像数据输出至刃边法计算边缘扩展函数模块和图像数据边界处理模块;刃边法边缘扩展函数模块使用原始图像数据计算得到相机边缘扩展函数(LSF);线扩展函数计算模块使用LSF计算得到单行图像线扩展函数(ESF);归一化MTF函数计算模块使用ESF计算得到相机系统MTF序列曲线函数。MTFC数字滤波器系数计算存储模块使用MTF计算、存储并输出行向、列向数字滤波器系数。图像数据边界处理模块对原始图像数据做边界处理;行向图像数字滤波模块使用行向数字滤波器系数对边界处理后数据滤波;行向滤波图像抑噪模块对行向滤波后数据做抑噪处理;行向滤波图像存储模块存储并输出行向滤波抑噪后图像数据;列向图像数字滤波模块对行向滤波抑噪后图像数据做列向滤波处理;列向滤波图像抑噪模块对列向滤波后数据做抑噪处理;列向滤波抑噪处理后数据存入复原后图像存储模块。
2.根据权利要求1所述的一种大幅提高成像质量的图像复原系统,其特征在于:所述刃边法计算边缘扩展函数模块实现如下:
(1)选取含有地面边缘纹理的某遥感图像做边缘检测;
(2)抠取具有刀刃形状的目标子图,保证刀刃边缘两侧的灰度分布较均匀无其他因素干扰;
(3)根据边缘成像的灰度分布拟合出边缘扩展函数曲线。
3.根据权利要求1所述的一种大幅提高成像质量的图像复原系统,其特征在于:所述线扩展函数计算模块实现如下:
(1)对边缘扩展函数曲线序列数据做简单差分,即计算x(n)-x(n-1),并取max(x(n)-x(n-1)),即可检测得到最大斜率位置;
(2)对图像的每一行数据点进行三次样条插值,通常选取插值分辨率为0.05,即在每两个数据点之间插入20个数据点,插值结果为一条灰度曲线,即为该行的边缘扩展函数ESF;
(3)对各行重排数据得到的ESF求平均,以最小化误差及降低噪声;
(4)对ESF做差分运算LSF(n)=ESF(n)-ESF(n-1),得到线扩展函数LSF曲线。
4.根据权利要求1所述的一种大幅提高成像质量的图像复原系统,其特征在于:所述归一化MTF函数计算模块实现如下:
(1)对线扩展函数曲线LSF做离散傅里叶变换后取各分量的模MTF(n)=|DFT(LSF(n))|;
(2)以变换后的直流分量即第一个MTF值为标准做归一化处理norm_MTF(n)=MTF(n)/MTF(1),即得到MTF曲线序列;
(3)使用归一化频率轴,Nyquist频率为截止频率的一半
Nyquist_frequency=(whole_data_size×resolution)/2+1=(Number_of_trimmed_pixel)/2+1
其中:Nyquist_frequency为Nyquist频率点,whole_data_size为包含所有插值点的截取后LSF亚像素级宽度,resolution为插值分辨率,Number_of_trimmed_pixel为截取后的LSF像素级宽度。
5.根据权利要求1所述的一种大幅提高成像质量的图像复原系统,其特征在于:所述MTFC数字滤波器系数计算存储模块实现如下:
(1)假定成像系统为线性移不变系统,通用图像退化模型来表示:f(x,y)*h(x,y)+n(x,y)=g(x,y)。在这个模型中,一个作用在输入图像f(x,y)上的系统h(x,y)与加性噪声n(x,y)的联合作用导致产生退化图像g(x,y);
(2)根据最小二乘原理,采用被系统h(x,y)模糊的图像与观测到的图像g(x,y)之间的差在均方意义最小的准则下来估计结果;
f ^ ( x , y ) = F - 1 [ F ^ ( u , v ) ] = F - 1 [ G ( u , v ) H ( u , v ) ]
在频域下,G(u,v)、F(u,v)、H(u,v)分别是g(x,y)、f(x,y)、h(x,y)的傅立叶变换,系统H(u,v)作为滤波函数与F(u,v)的乘积是退化图像g(x,y)的傅里叶变换;
(3)把系统函数H(u,v)看作复数,对幅值作归一化,将H(u,v)替换为k·MTF(u,v),k为系统函数归一化系数,重写上式为:
f ^ ( x , y ) = F - 1 [ 1 k &CenterDot; MTF ( u , v ) G ( u , v ) ]
(4)为了达到提升MTF的效果,对上式进行修改为
F ^ ( u , v ) = G ( u , v ) H ( u , v ) = [ 1 r &prime; ( u , v ) k &CenterDot; MTF ( u , v ) ] G ( u , v )
其中,r′(u,v)是实现在频域下不同频率段的补偿,从而有针对性提升MTF。逆滤波复原局限于中低频区域进行,且适用于信噪比较高的情况;
(5)再令 r ( u , v ) = 1 r &prime; ( u , v ) k &CenterDot; MTF ( u , v ) , F ^ ( u , v ) = G ( u , v ) [ r ( u , v ) ] &alpha; , 其中α为补偿因子,恒大于等于0,此时MTFC补偿函数为: MTFC ( u , v ) = [ r ( u , v ) ] &alpha; = [ 1 r &prime; ( u , v ) k &CenterDot; MTF ( u , v ) ] &alpha; . 通过调整r′及α实现对MTF的补偿;
(6)选择重点频段对r′进行补偿和运用不同的补偿因子进行优化,r′选用凸曲线比凹曲线或直线对图像的补偿效果更好,α进一步优化补偿效果,在图像清晰度和噪声之间寻求平衡,信噪比、熵的增加和均方梯度相对最大能达到补偿效果和控制噪声的目的;
(7)计算输出MTFC数字滤波器系数后存入PROM。
6.根据权利要求1所述的一种大幅提高成像质量的图像复原系统,其特征在于:所述图像数据边界处理模块实现如下:
(1)图像数据按次序输入,假设一行CCD像元为8190个,图像共有1024行数据,使用D1,D2…,D8190分别代表第1、2至第8190个像元的输出值。每个像元输入时,带有图像数据有效使能;
(2)使用图像数据效使能激发计数器1,输出结果n1
(3)当1≤n1≤5或8186≤n1≤8190时,使用此判断有效使能激发计数器2,输出结果n2
(4)将n2作为选择器输入,datain表示为选择器输入数据,dataout表示选择器输出数据,条件对应关系如下:
if n2=0dataout=D5
if n2=1dataout=D4
if n2=2dataout=D3
if n2=3dataout=D2
if n2=4dataout=D1
if n2=5dataout=datain
if n2=6dataout=D8190
if n2=7dataout=D8189
if n2=8dataout=D8188
if n2=9dataout=D8187
if n2=10dataout=D8186
最后得到加边界的图像数据序列为:
D5,D4,D3,D2,D1,D1,D2,D3,…,D8190,D8190,D8189,D8180,D8187,D8186。在原始图像线阵数据一行的开始和结束分别填补5个像元数据,为这行数据开始5个和结束5个像元数据作对称处理。
7.根据权利要求1所述的一种大幅提高成像质量的图像复原系统,其特征在于:所述行向图像数字滤波模块实现如下:
MTFC数字滤波算法图像水平即方向和垂直即方向的复原分开进行,先水平方向上的卷积、抑噪,再垂直方向上的卷积、抑噪;
水平方向卷积系数为11个。每个卷积核系数的范围为[-8,8],16bit存储在Prom中,第1位为符号位,第2-4位为整数位,第5-16位为小数位;
算法为:
dataout=k1·Dn-5+k2·Dn-4+k3·Dn-3+k4·Dn-2+k5·Dn-1+k6·Dn+k7·Dn+1+k8·Dn+2+k9·Dn+3+k10·Dn+4+k11·Dn+5
其中,n=1,2,…8190。k1,k2,…k11为行向滤波器系数。由原始图像数据边界处理模块计算,D-5,D-4,D-3,D-2,D-1分别对应图像数据D5,D4,D3,D2,D1。D8191,D8192,D8193,D8194,D8195分别对应图像数据D8190,D8189,D8188,D8187,D8186
8.根据权利要求1所述的一种大幅提高成像质量的图像复原系统,其特征在于:所述行向滤波图像抑噪模块实现如下:
(1)假设某一像元原始灰度值为i0,卷积运算后结果为i1,取i2=i0+k(i1-i0)为最终复原值。当k=1时,则复原后MTF为全额补偿,若k=0.5,则MTF为半额补偿,若k=0,则为零补偿;
(2)取Δ=K||i1-i0||,其中K为抑噪阈值参数。
(3)当Δ≥i0时,k=1,即i2=i1
(4)当Δ<i0时,k=Δ/2n,即i2=i0+Δ(i1-i0)/2n
n为像元量化位数。其中此处除法可用移位代替
最终值复原值为i2
k为浮点数,具体实现时可以不计算k,而直接计算i2,从而避免浮点数的运算。
9.根据权利要求1所述的一种大幅提高成像质量的图像复原系统,其特征在于:所述行向滤波图像存储模块实现如下:行向滤波处理和抑噪处理后同时十列数据图像数据缓存到RAM中,以进行列方向数据处理。
CN201310629499.8A 2013-11-29 2013-11-29 一种大幅提高成像质量的图像复原系统 Active CN103679652B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201310629499.8A CN103679652B (zh) 2013-11-29 2013-11-29 一种大幅提高成像质量的图像复原系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201310629499.8A CN103679652B (zh) 2013-11-29 2013-11-29 一种大幅提高成像质量的图像复原系统

Publications (2)

Publication Number Publication Date
CN103679652A true CN103679652A (zh) 2014-03-26
CN103679652B CN103679652B (zh) 2017-04-19

Family

ID=50317108

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201310629499.8A Active CN103679652B (zh) 2013-11-29 2013-11-29 一种大幅提高成像质量的图像复原系统

Country Status (1)

Country Link
CN (1) CN103679652B (zh)

Cited By (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104537646A (zh) * 2014-12-12 2015-04-22 南京理工大学 遥感图像的多角度自动mtf估计方法
CN104732531A (zh) * 2015-03-11 2015-06-24 中国空间技术研究院 一种高分辨率遥感图像信噪比曲线自适应获取方法
CN104905784A (zh) * 2015-06-25 2015-09-16 南京大学 一种减小生物电采集波形失真的方法
CN105069313A (zh) * 2015-08-24 2015-11-18 北京理工大学 基于相位非线性重采样拟合刃边的在轨mtf估计方法
CN105354837A (zh) * 2015-10-13 2016-02-24 哈尔滨工程大学 一种基于人体解剖结构的x光消化道图像清晰度评价方法
CN106127699A (zh) * 2016-06-16 2016-11-16 哈尔滨理工大学 一种道路监控随机运动模糊图像快速复原仿真系统
CN106157253A (zh) * 2015-04-17 2016-11-23 瑞昱半导体股份有限公司 图像处理装置与图像处理方法
CN106373097A (zh) * 2016-08-29 2017-02-01 合肥康胜达智能科技有限公司 一种图像处理方法
CN107424583A (zh) * 2017-07-03 2017-12-01 威创集团股份有限公司 异形图像的显示数据处理方法和系统
CN107802269A (zh) * 2017-10-30 2018-03-16 温州医科大学 实时单次快照多频解调空间频域成像方法
CN107854872A (zh) * 2017-12-11 2018-03-30 北京北排装备产业有限公司 一种非金属链条式刮泥机水下监测装置及其使用方法
CN108652630A (zh) * 2018-04-13 2018-10-16 廉江市人民医院 一种神经外科脑部手术辅助仪
CN108765286A (zh) * 2018-05-18 2018-11-06 浙江大学 一种辐射调制函数保真的数字图像插值算法
CN108961163A (zh) * 2018-06-28 2018-12-07 长光卫星技术有限公司 一种高分辨率卫星影像超分辨重建方法
CN109035139A (zh) * 2018-06-28 2018-12-18 长光卫星技术有限公司 一种高分辨率卫星影像调制传递函数补偿方法
CN110460765A (zh) * 2018-05-02 2019-11-15 展讯通信(上海)有限公司 空间频率响应曲线的获取方法、装置及电子设备
CN111610001A (zh) * 2020-05-25 2020-09-01 中国科学院长春光学精密机械与物理研究所 一种宽幅遥感图像mtf地面模拟测试装置

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6154574A (en) * 1997-11-19 2000-11-28 Samsung Electronics Co., Ltd. Digital focusing method and apparatus in image processing system
US20060013479A1 (en) * 2004-07-09 2006-01-19 Nokia Corporation Restoration of color components in an image model
CN101635050A (zh) * 2009-06-26 2010-01-27 武汉大学 一种图像复原方法
CN102692273A (zh) * 2012-05-31 2012-09-26 中国资源卫星应用中心 一种干涉型高光谱成像仪的mtf在轨检测方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6154574A (en) * 1997-11-19 2000-11-28 Samsung Electronics Co., Ltd. Digital focusing method and apparatus in image processing system
US20060013479A1 (en) * 2004-07-09 2006-01-19 Nokia Corporation Restoration of color components in an image model
CN101635050A (zh) * 2009-06-26 2010-01-27 武汉大学 一种图像复原方法
CN102692273A (zh) * 2012-05-31 2012-09-26 中国资源卫星应用中心 一种干涉型高光谱成像仪的mtf在轨检测方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
ZHAI GUOFANG ET AL.: "Design of high-speed multi-channel TDICCD image acquisition and processing system based on FPGA", 《INTERNATIONAL SYMPOSIUM ON PHOTOELECTRONIC DETECTION AND IMAGING 2013: IMAGING SPECTROMETER TECHNOLOGIES AND APPLICATIONS》 *
吴昀昭 等: "北京一号小卫星MTF在轨测量与图像复原", 《遥感应用》 *
柴雅琼 等: "改进的MTF遥感影像复原算法研究", 《遥感应用》 *

Cited By (27)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104537646A (zh) * 2014-12-12 2015-04-22 南京理工大学 遥感图像的多角度自动mtf估计方法
CN104537646B (zh) * 2014-12-12 2017-06-27 南京理工大学 遥感图像的多角度自动mtf估计方法
CN104732531B (zh) * 2015-03-11 2017-07-07 中国空间技术研究院 一种高分辨率遥感图像信噪比曲线自适应获取方法
CN104732531A (zh) * 2015-03-11 2015-06-24 中国空间技术研究院 一种高分辨率遥感图像信噪比曲线自适应获取方法
CN106157253B (zh) * 2015-04-17 2019-09-03 瑞昱半导体股份有限公司 图像处理装置与图像处理方法
CN106157253A (zh) * 2015-04-17 2016-11-23 瑞昱半导体股份有限公司 图像处理装置与图像处理方法
CN104905784A (zh) * 2015-06-25 2015-09-16 南京大学 一种减小生物电采集波形失真的方法
CN105069313A (zh) * 2015-08-24 2015-11-18 北京理工大学 基于相位非线性重采样拟合刃边的在轨mtf估计方法
CN105069313B (zh) * 2015-08-24 2017-11-10 北京理工大学 基于相位非线性重采样拟合刃边的在轨mtf估计方法
CN105354837B (zh) * 2015-10-13 2017-12-19 哈尔滨工程大学 一种基于人体解剖结构的x光消化道图像清晰度评价方法
CN105354837A (zh) * 2015-10-13 2016-02-24 哈尔滨工程大学 一种基于人体解剖结构的x光消化道图像清晰度评价方法
CN106127699A (zh) * 2016-06-16 2016-11-16 哈尔滨理工大学 一种道路监控随机运动模糊图像快速复原仿真系统
CN106127699B (zh) * 2016-06-16 2018-10-19 哈尔滨理工大学 一种道路监控随机运动模糊图像快速复原仿真系统
CN106373097A (zh) * 2016-08-29 2017-02-01 合肥康胜达智能科技有限公司 一种图像处理方法
CN107424583A (zh) * 2017-07-03 2017-12-01 威创集团股份有限公司 异形图像的显示数据处理方法和系统
WO2019085114A1 (zh) * 2017-10-30 2019-05-09 温州医科大学 实时单次快照多频解调空间频域成像方法
CN107802269A (zh) * 2017-10-30 2018-03-16 温州医科大学 实时单次快照多频解调空间频域成像方法
CN107854872A (zh) * 2017-12-11 2018-03-30 北京北排装备产业有限公司 一种非金属链条式刮泥机水下监测装置及其使用方法
CN108652630A (zh) * 2018-04-13 2018-10-16 廉江市人民医院 一种神经外科脑部手术辅助仪
CN110460765A (zh) * 2018-05-02 2019-11-15 展讯通信(上海)有限公司 空间频率响应曲线的获取方法、装置及电子设备
CN110460765B (zh) * 2018-05-02 2020-12-18 展讯通信(上海)有限公司 空间频率响应曲线的获取方法、装置及电子设备
CN108765286A (zh) * 2018-05-18 2018-11-06 浙江大学 一种辐射调制函数保真的数字图像插值算法
CN108765286B (zh) * 2018-05-18 2020-06-30 浙江大学 一种辐射调制函数保真的数字图像插值算法
CN109035139A (zh) * 2018-06-28 2018-12-18 长光卫星技术有限公司 一种高分辨率卫星影像调制传递函数补偿方法
CN108961163A (zh) * 2018-06-28 2018-12-07 长光卫星技术有限公司 一种高分辨率卫星影像超分辨重建方法
CN111610001A (zh) * 2020-05-25 2020-09-01 中国科学院长春光学精密机械与物理研究所 一种宽幅遥感图像mtf地面模拟测试装置
CN111610001B (zh) * 2020-05-25 2021-09-21 中国科学院长春光学精密机械与物理研究所 一种宽幅遥感图像mtf地面模拟测试装置

Also Published As

Publication number Publication date
CN103679652B (zh) 2017-04-19

Similar Documents

Publication Publication Date Title
CN103679652A (zh) 一种大幅提高成像质量的图像复原系统
CN102609904B (zh) 双变量非局部平均滤波x射线图像消噪方法
CN104657945B (zh) 复杂背景下多尺度时空联合滤波的红外小目标检测方法
CN101727658B (zh) 图像处理方法及装置
CN101882304B (zh) 一种sar图像自适应去噪和特征增强方法
CN106952228A (zh) 基于图像非局部自相似性的单幅图像的超分辨率重建方法
CN101847257B (zh) 基于非局部均值与多级定向图像的图像降噪方法
CN101699509B (zh) 一种利用气象数据进行大气模糊遥感影像恢复的方法
CN103136734B (zh) 一种凸集投影超分辨率图像重建时边缘晕轮效应的抑制方法
CN103530896A (zh) 一种红外图像的图像压缩和细节增强方法
CN102572499B (zh) 基于小波变换多分辨率预测的无参考图像质量评价方法
CN109146797A (zh) 一种基于Lp伪范数与交叠组稀疏的脉冲噪声古籍图像修复方法
CN109377464B (zh) 一种红外图像的双平台直方图均衡化方法及其应用系统
CN103279935A (zh) 基于map算法的热红外遥感图像超分辨率重建方法及系统
CN103473755B (zh) 基于变化检测的sar图像稀疏去噪方法
CN106169181A (zh) 一种图像处理方法及系统
CN102332153A (zh) 基于核回归的图像压缩感知重构方法
CN104867122A (zh) 一种红外自适应非均匀性校正及细节增强级联处理方法
CN109949256B (zh) 一种基于傅里叶变换的天文图像融合方法
CN104680485A (zh) 一种基于多分辨率的图像去噪方法及装置
CN103971345A (zh) 一种基于改进双边滤波的图像去噪方法
CN102789634A (zh) 一种获取光照均一化图像的方法
CN102339460B (zh) 卫星图像自适应复原方法
CN103606130A (zh) 红外退化图像自适应复原方法
CN115409872B (zh) 一种水下摄像机图像优化方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant