CN103186888B - 一种去除ct图像噪声的方法及装置 - Google Patents
一种去除ct图像噪声的方法及装置 Download PDFInfo
- Publication number
- CN103186888B CN103186888B CN201110462903.8A CN201110462903A CN103186888B CN 103186888 B CN103186888 B CN 103186888B CN 201110462903 A CN201110462903 A CN 201110462903A CN 103186888 B CN103186888 B CN 103186888B
- Authority
- CN
- China
- Prior art keywords
- image
- noise
- calculating
- module
- value
- 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
- 238000000034 method Methods 0.000 title claims abstract description 61
- 238000001914 filtration Methods 0.000 claims abstract description 55
- 230000002708 enhancing effect Effects 0.000 claims abstract description 14
- 238000009792 diffusion process Methods 0.000 claims abstract description 13
- 210000001519 tissue Anatomy 0.000 claims description 43
- 238000004458 analytical method Methods 0.000 claims description 17
- 238000010586 diagram Methods 0.000 claims description 14
- 239000002245 particle Substances 0.000 claims description 10
- 101100273916 Schizosaccharomyces pombe (strain 972 / ATCC 24843) wip1 gene Proteins 0.000 claims description 9
- 238000005070 sampling Methods 0.000 claims description 9
- 238000012935 Averaging Methods 0.000 claims description 6
- 238000005260 corrosion Methods 0.000 claims description 6
- 230000007797 corrosion Effects 0.000 claims description 6
- 238000009499 grossing Methods 0.000 claims description 6
- 210000004872 soft tissue Anatomy 0.000 claims description 6
- 230000009466 transformation Effects 0.000 claims description 6
- 238000000605 extraction Methods 0.000 claims description 4
- 230000008520 organization Effects 0.000 claims description 4
- 238000003706 image smoothing Methods 0.000 claims description 3
- 230000000877 morphologic effect Effects 0.000 claims description 3
- 230000000007 visual effect Effects 0.000 claims description 3
- 230000010339 dilation Effects 0.000 claims 2
- 101100460203 Schizosaccharomyces pombe (strain 972 / ATCC 24843) new2 gene Proteins 0.000 claims 1
- 230000003628 erosive effect Effects 0.000 claims 1
- 238000002591 computed tomography Methods 0.000 description 25
- 239000000284 extract Substances 0.000 description 2
- 230000005855 radiation Effects 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 125000004122 cyclic group Chemical group 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
Landscapes
- Image Processing (AREA)
- Image Analysis (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
本发明名称为“一种去除CT图像噪声的方法及装置”。本发明涉及图像去噪技术领域，公开了一种去除CT图像噪声的方法及装置，该方法包括估计图像的组织权重、估计图像的噪声等级、计算去噪参数、对图像进行各向异性扩散滤波、对滤波输出的图像进行边缘增强、对图像进行细节增强及对比度修正、裁切图像并输出结果的步骤。本装置包括用于估计图像的组织权重的模块，用于估计图像的噪声等级的模块，用于计算去噪参数的模块，用于对图像进行各向异性扩散滤波的模块，用于对滤波输出的图像进行边缘增强的模块，用于对图像进行细节增强及对比度修正的模块，用于裁切图像并输出结果的模块。本方法及装置能够在有效去除CT图像高频噪声的同时，保持图像的边缘及原有的对比度。
Description
技术领域
本发明涉及图像去噪技术领域，特别是一种去除CT图像噪声的方法及装置。
背景技术
计算机断层扫描(CT)作为一种成熟的在临床上普遍认可的检查方法，它的辐射剂量问题已经引起了人们越来越多的关注。在目前的CT中，常通过降低管电流、管电压以减少CT剂量，但这样将会增加图像的块状和粒状噪声，降低图像质量，影响临床医生的诊断。
当前提高低剂量CT图像质量的方法主要分为基于投影空间数据和基于图像空间数据两大类。各向异性滤波器和小波的多尺度分析即为基于图像空间数据的非线性自适应滤波降噪方法，各向异性滤波器能够使CT图像中的低频段小梯度噪声得到有效遏制，然而CT的高频段噪声通常具有较大的梯度，因而通常被各向异性扩散滤波器视为边缘而不能去除。同时，这种降噪方法通常会在滤除噪声的同时，损失图像的原有边缘和低对比度分辨力。针对上述现有技术存在的问题，本发明提出的技术方案，能够提取并滤除CT高频段噪声，在降噪的同时，能保持图像边缘和原有的对比度。
发明内容
针对现有技术不能有效地去除CT图像的高频噪声，并且会在去噪的同时损失图像边缘和对比度的问题，本发明提供的一种去除CT图像噪声的方法，包括如下步骤：
步骤一：估计图像的组织权重；
步骤二：估计图像的噪声等级；
步骤三：计算去噪参数；
步骤四：对图像进行各向异性扩散滤波；
步骤五：对滤波输出的图像进行边缘增强；
步骤六：对图像进行细节增强及对比度修正；
步骤七：裁切图像并输出结果。
所述步骤一包括如下步骤：
a)配置初始参数：组织均值Tissue、软组织第一门限值TissueUp、软组织第二门限值TissueDown、组织降噪权重SRange、平滑处理窗宽度K；
b)在显示视野DFOV范围内，对输入的图像做平滑宽度为K的图像平滑处理，得到每个像素点的平滑值smoothimage(i，j)；
c)算出输入图像的每个像素点的组织权重：
其中，x＝abs(smoothimage(i，j)-Tissue)，abs()表示取绝对值。
所述步骤二包括如下步骤：
a)对输入的图像进行降采样；
b)选取估算区域：在降采样得到的图像上选取一块方形区域作为估算区域EstimationRegion；
c)计算估算区域的组织权重：对估算区域内的每个像素值用所述估计图像的组织权重的方法估算其组织权重；
d)对估算区域加权：用组织权重对估算区域内的每一个像素点进行加权运算，得到加权后的图像像素值EstimationMap；
e)计算图像的内外包络：对加权后的图像像素值分别进行基于形态学的开运算和闭运算，分别得到图像的内包络和外包络；
f)计算噪声分布图：对内、外包络求平均值得到平滑分布图SmoothMap，并计算出噪声分布图NoiseMap＝abs(EstimationMap-SmoothMap)，abs()表示求绝对值；
g)直方图分析：对该噪声分布图做直方图分析，得到噪声等级估算值NoiseLevel。
所述估算区域位于降采样得到的图像的中心位置，为正方形区域。
所述基于形态学的开运算依次包括计算滤波核、对所述加权后的图像像素值EstimationMap进行腐蚀操作和对腐蚀操作得到的结果进行膨胀操作的步骤。
所述的计算滤波核进一步包括如下步骤：
a)求出所述加权后的图像像素值EstimationMap的最大值MaxValue和平均值MeanValue；
b)设定滤波核大小KernelSize、滤波核范围Range以及滤波参数Delta；
c)计算出每个像素点对应的滤波核初始值：
其中，(i，j)为每
个像素点值，FWSize＝(KernelSize+1)/2；
d)计算每个像素点对应的滤波核的准确值：
Kernel(i，j)＝Kernel_init(i，j)*max(Kernel_init(i，j))*TH，其中，TH＝(MaxValue-MeanValue)/Range，max()表示取最大。
所述基于形态学的闭运算依次包括计算滤波核、对所述加权后的图像像素值EstimationMap进行膨胀操作和对膨胀操作的结果进行腐蚀操作的步骤。
所述的计算滤波核进一步包括如下步骤：
a)求出所述加权后的图像像素值EstimationMap的最大值MaxValue和平均值MeanValue；
b)设定滤波核大小KernelSize、滤波核范围Range以及滤波参数Delta；
c)计算出每个像素点对应的滤波核初始值：
其中，(i，j)为每个像素点值，FWSize＝(KernelSize+1)/2；
d)计算每个像素点对应的滤波核的准确值：
Kernel(i，j)＝Kernel_init(i，j)*max(Kernel_init(i，j))*TH，其中，TH＝(MaxValue-MeanValue)/Range，max()表示取最大。
所述的直方图分析包括如下步骤：
a)计算出所述噪声分布图NoiseMap中的像素最大值MaxValue和最小值MinValue，进而得到图像的动态范围Range＝MaxValue-MinValue；
b)将图像的动态范围分成N份小区间，每一份小区间的步长Step＝Range/N；
c)统计每个小区间内的像素点数Hist(i)，对于第i个区间而言，就是对[MinValue+Step*(i-1)，Minvalue+Step*i]区间内的像素点数进行统计；
d)对Hist(i)进行概率分布统计，即计算出公式中i的值，M表示噪声分布图NoiseMap中有效噪声的比例；
e)计算噪声分布图NoiseMap中有效噪声点的最大值：Threshold＝(I-1)*Step+MinValue，其中的I即为4)步算出的i的值；
f)对噪声分布图NoiseMap中所有小于有效噪声点的最大值Threshold的值求平均值，得到图像的噪声等级NoiseLevel＝mean(NoiseMap＜Threshold)。
所述的噪声分布图NoiseMap中有效噪声的比例M＝0.95。
所述步骤三包括如下步骤：
a)设定参数K和fwidth；
b)算出每个像素点位置校正值(x，y)表示把图像的中心点作为坐标原点时，每个像素点的坐标；
c)根据所述步骤一得到的组织权重、步骤二得到的噪声等级和位置校正值，计算出每个像素点的降噪参数：
所述步骤四依次包括如下步骤：
a)计算出应当被去除的噪声值：
其中，
Diff＝ds*InputImage(x+m，y+n)-InputImage(x，y)，m＝[12313123]，n＝[11122333]，ds＝[r，1，r，1，1，r，1，r]，InputImage(x，y)为输入图像的像素值，.*表示点乘运算，j为各向异性扩散滤波系数，NoiseScale为所述步骤三得到的降噪参数；
b)计算并输出去噪以后的图像像素值：
DeNoiseImage＝InputImage-NoiseReduce。
所述b)步完成之后，还包含如下步骤：将b)步得到的去噪以后的图像像素值DeNoiseImage赋值给InputImage，更换不同的各向异性扩散滤波系数值j，然后返回a)步，在至少依此方法循环执行三次后，停止所有步骤并输出去噪后的图像像素值DeNoiseImage。
所述步骤五包括如下步骤：
a)去除大颗粒噪声：用带通滤波器对图像Image进行带通滤波，得到图像中较大颗粒的噪声Noise以及去除该噪声以后的图像Image＝Image-Noise；
b)提取图图像边缘：更改带通滤波器的滤波频段，对去除噪声以后的图像进行边缘提取Edge＝BandPass(Image)；
c)计算有效边缘：对边缘进行直方图分析，得到有效的边缘UsefulEdge；
d)边缘增强：边缘增强以后的图像EnhancedImage＝Image+Scale*UsefulEdge，Scale表示对边缘进行增强的倍数。
所述步骤六包括如下步骤：
a)计算图像的第一层的背景信息和细节信息：对图像做哈尔小波变换，得到：
[背景信息Background1细节信息Foreground1]＝HarrTransform(Image)；
b)去除背景信息中的小噪声：对提取到的背景信息进行各向异性滤波，得到去除小噪声后新的背景信息BackgroundNew1＝AnisotropicDiffusionFilter(Background1)；
c)计算新的背景信息中的第二层细节信息和背景信息：对新的背景信息BackgroundNew1做哈尔小波变换，得到第二层的背景信息和细节信息：[背景信息Background2细节信息Foreground2]＝HarrTransform(BackgroundNew1)；
d)去除第二层背景信息中的小噪声：对第二层的背景信息进行各向异性滤波，得到去除小噪声后新的第二层背景信息：
BackgroundNew2＝AnisotropicDiffusionFilter(Background2)；
e)对第二层的图像细节进行尺度为K2的增强：增强后的第二层图像，Enhance2＝K2*Foreground2+BackgroundNew2；
f)对增强后的第二层图像做哈尔小波反变换：
输出Out2＝HarrInverse(Enhance2)
g)计算第一层图像细节增强的对比度修正因子：
输出
h)对第一层的图像细节进行尺度为K1的增强及对比度修正：增强后的
第一层图像：
Enhance1＝Correctfactor*K1*Foreground1+Background1
i)对增强后的第一层图像做哈尔小波反变换：其输出Out1＝HarrInverse(Enhance1)。
所述步骤七中的裁切图像是指截取图像中的显示视野可见区域，即：DFOV区域。
本发明提供的一种去除CT图像噪声的装置，包括：用于估计图像的组织权重的模块，用于估计图像的噪声等级的模块，用于计算去噪参数的模块，用于对图像进行各向异性扩散滤波的模块，用于对滤波输出的图像进行边缘增强的模块，用于对图像进行细节增强及对比度修正的模块，用于裁切图像并输出结果的模块。
所述的用于估计图像的噪声等级的模块进一步包括：用于对输入的图像进行降采样的模块；用于选取估算区域的模块，其在降采样得到的图像上选取一块方形区域作为估算区域EstimationRegion；用于估算区域的组织权重的模块，其对估算区域内的每个像素值用所述估计图像的组织权重的方法估算其组织权重；用于对估算区域加权的模块，其用组织权重对估算区域内的每一个像素点进行加权运算，得到加权后的图像像素值EstimationMap；用于计算图像的内外包络的模块，其对加权后的图像像素值分别进行基于形态学的开运算和闭运算，分别得到图像的内包络和外包络；用于计算噪声分布图的模块，其对内、外包络求平均值得到平滑分布图SmoothMap，并计算出噪声分布图NoiseMap＝abs(EstimationMap-SmoothMap)，abs()表示求绝对值；用于直方图分析的模块，其对该噪声分布图做直方图分析，得到噪声等级估算值NoiseLevel。
所述的用于对滤波输出的图像进行边缘增强的模块进一步包括：用于去除大颗粒噪声的模块，其用带通滤波器对图像Image进行带通滤波，得到图像中较大颗粒的噪声Noise以及去除该噪声以后的图像Image＝Image-Noise；用于提取图图像边缘的模块，其更改带通滤波器的滤波频段，对去除噪声以后的图像进行边缘提取Edge＝BandPass(Image)；用于计算有效边缘的模块，其对边缘进行直方图分析，得到有效的边缘UsefulEdge；用于边缘增强的模块，其边缘增强以后的图像EnhancedImage＝Image+Scale*UsefulEdge，Scale表示对边缘进行增强的倍数。
所述的用于对图像进行细节增强及对比度修正的模块进一步包括用于计算图像的第一层的背景信息和细节信息的模块，其对图像做哈尔小波变换，得到：[Background1Foreground1]＝HarrTransform(Image)；用于去除背景信息中的小噪声的模块，其对提取到的背景信息进行各向异性滤波，得到去除小噪声后新的背景信息BackgroundNew1＝AnisotropicDiffusionFilter(Background1)；用于计算新的背景信息中的第二层细节信息和背景信息的模块，其对新的背景信息BackgroundNew1做哈尔小波变换，得到第二层的背景信息和细节信息：[Background2 Foreground2]＝HarrTransform(BackgroundNew1)；用于去除第二层背景信息中的小噪声的模块，其对第二层的背景信息进行各向异性滤波，得到去除小噪声后新的第二层背景信息BackgroundNew2＝AnisotropicDiffusionFilter(Background2)；用于对第二层的图像细节进行尺度为K2的增强的模块，其增强后的第二层图像Enhance2＝K2*Foreground2+BackgroundNew2；用于对增强后的第二层图像做哈尔小波反变换的模块，其输出Out2＝HarrInverse(Enhance2)；用于计算第一层图像细节增强的对比度修正因子的模块，其输出 用于对第一层的图像细节进行尺度为K1的增强及对比度修正的模块，其增强后的第一层图像：Enhance1＝Correctfactor*K1*Foreground1+Background1；用于对增强后的第一层图像做哈尔小波反变换的模块，其输出Out1＝HarrInverse(Enhance1)。
与现有技术相比，本发明提供的一种去除CT图像噪声的方法及装置，其有益的技术效果是：
1、在去噪的同时能保持图像的边缘和细节；
2、在去噪的同时不会削弱图像的对比度；
3、不会引入新的假影；
4、在保持图像质量不变的情况下减少了对病员的辐射剂量。
附图说明
图1是本发明提供的一种去除CT图像噪声的方法的总体流程图；
图2是本发明提供的一种去除CT图像噪声的方法中的估计图像的噪声等级的具体流程图；
图3是本发明提供的一种去除CT图像噪声的方法中的对滤波输出的图像进行边缘增强的具体流程图；
图4是本发明提供的一种去除CT图像噪声的方法中的对图像进行细节增强及对比度修正的具体流程图；
图5是本发明提供的一种去除CT图像噪声的装置的示意性框图；
图6是用本发明提供的一种去除CT图像噪声的方法及装置处理一幅CT图像之前和之后的对比图，图6(a)是原始输入的CT图像，图6(b)是经过本方法/装置处理以后的CT图像。
具体实施方式
下面结合附图和具体实施方式对本发明作进一步说明。
本实施例披露了一种去除CT图像噪声的方法，包括如下步骤：
1、估计图像的组织权重
先配置初始参数：组织均值Tissue、软组织第一门限值TissueUp、软组织第二门限值TissueDown、组织降噪权重SRange、均值滤波窗宽度K；
然后在显示视野(DFOV)范围内，对输入的图像做平滑宽度为K的图像平滑处理，得到每个像素点的平滑值：
其中InputImage(i，j)表示原始的输入图像的像素值；
最后算出输入图像的每个像素点的组织权重：
其中，x＝abs(smoothimage(i，j)-Tissue)，abs()表示取绝对值。
2、估计图像的噪声等级，本步骤包括如下步骤：
1)对原始的输入图像进行降采样，比如按照50％的比例进行降采样。
2)在降采样以后的图像上选取一块区域，作为噪声等级估计的估计区域，比如，选取位于正中心位置的正方形区域。
3)用第1步中描述的估计图像的组织权重的方法，对该选定的估计区域计算出每个像素点的组织权重，并用该组织权重对每个像素点的值进行加权，得到加权后的图像像素值EstimationMap。
4)对上一步得到的加权后的图像像素值分别进行形态学滤波运算中的开运算和闭运算，分别得到图像的内包络InProfile和外包络OutProfile。
其中，开运算使用的伪代码如下：
闭运算使用的伪代码如下：
其中，膨胀操作Dalite(f)＝max{f(s+x，t+y)+kernel(x，y)|(s+x，t+y)∈f，(x，y)∈kernel}，腐蚀操作Erode(f)＝min{f(s+x，t+y)-kernel(x，y)|(s+x，t+y)∈f，(x，y)∈kernel}。
上述伪代码中，滤波核大小KernelSize、滤波核范围Range以及滤波参数Delta均根据实际的工程应用情况来设定。
5)对上一步得到的内、外包络求平均值得到平滑分布图SmoothMap＝(InProfile+OutProfile)/2，并计算出噪声分布图NoiseMap＝abs(EstimationMap-SmoothMap)，abs()表示求绝对值。
6)对该噪声分布图做直方图分析，得到噪声等级估算值NoiseLevel。
具体步骤包括：
a)计算出所述噪声分布图NoiseMap中的像素最大值MaxValue和最小值MinValue，进而得到图像的动态范围Range＝MaxValue-MinValue；
b)将图像的动态范围分成N份小区间，每一份小区间的步长Step＝Range/N；
c)统计每个小区间内的像素点数Hist(i)，对于第i个区间而言，就是对[MinValue+Step*(i-1)，Minvalue+Step*i]区间内的像素点数进行统计；
d)对Hist(i)进行概率分布统计，即计算出下面公式中RateHist(i)＝M时i的值：M表示噪声分布图NoiseMap中有效噪声的比例，按照大量实验经验，我们认为在噪声图像中，大约95％的点为有效噪声点，剩下的值为误提取的边界或者结构等信息，因此本实施例取M＝0.95；
e)计算噪声分布图NoiseMap中有效噪声点的最大值：
Threshold＝(I-1)*Step+MinValue，其中的I即为4)步算出的i的值；
f)对噪声分布图NoiseMap中所有小于有效噪声点的最大值Threshold的值求平均值，得到图像的噪声等级NoiseLevel＝mean(NoiseMap＜Threshold)。
3、计算去噪参数
在这里，我们定义的去噪参数：
NoiseScale＝(噪声等级*位置校正值)/组织权重
其中，组织权重来自第1步的输出，噪声等级来自第2步的输出，位置校正值(x，y)表示把图像的中心点作为坐标原点时，每个像素点的坐标，其中的参数K和fwidth是根据CT图像扫描场景预先配置好的。
4、对图像进行各向异性扩散滤波
本实施例对图像做三次各向异性扩散滤波，每次的区别在于滤波器系数K的取值不同，这三次滤波该用伪代码表示如下：
其中，
InputImage表示原始的输入图像，，Diff＝ds*InputImage(x+m，y+n)-InputImage(x，y)，m＝[12313123]，n＝[11122333]，ds＝[r，1，r，1，1，r，1，r]，NoiseScale就是第3步的输出结果。
5、对滤波输出的图像进行边缘增强，具体包括：
a)去除大颗粒噪声：用带通滤波器对图像Image进行带通滤波，得到图像中较大颗粒的噪声Noise以及去除该噪声以后的图像Image＝Image-Noise；
b)提取图图像边缘：更改带通滤波器的滤波频段，对去除噪声以后的图像进行边缘提取Edge＝BandPass(Image)；
c)计算有效边缘：对边缘进行直方图分析，得到有效的边缘UsefulEdge；
d)边缘增强：边缘增强以后的图像EnhancedImage＝Image+Scale*UsefulEdge，Scale表示对边缘进行增强的倍数。
6、对图像进行细节增强及对比度修正，具体包括：
a)计算图像的第一层的背景信息和细节信息：对图像做哈尔小波变换，得到：[背景信息Background1细节信息Foreground1]＝HarrTransform(Image)；
b)去除背景信息中的小噪声：对提取到的背景信息进行各向异性滤波，得到去除小噪声后新的背景信息BackgroundNew1＝AnisotropicDiffusionFilter(Background1)；
c)计算新的背景信息中的第二层细节信息和背景信息：对新的背景信息BackgroundNew1做哈尔小波变换，得到第二层的背景信息和细节信息：[背景信息Background2细节信息Foreground2]＝HarrTransform(BackgroundNew1)；
d)去除第二层背景信息中的小噪声：对第二层的背景信息进行各向异性滤波，得到去除小噪声后新的第二层背景信息：
BackgroundNew2＝AnisotropicDiffusionFilter(Background2)；
e)对第二层的图像细节进行尺度为K2的增强：增强后的第二层图像，
Enhance2＝K2*Foreground2+BackgroundNew2；
f)对增强后的第二层图像做哈尔小波反变换：
输出Out2＝HarrInverse(Enhance2)
g)计算第一层图像细节增强的对比度修正因子：
输出
h)对第一层的图像细节进行尺度为K1的增强及对比度修正：增强后的第一层图像：
Enhance1＝Correctfactor*K1*Foreground1+Background1
i)对增强后的第一层图像做哈尔小波反变换：其输出Out1＝HarrInverse(Enhance1)。
7、裁切图像并输出结果：截取图像中的显示视野可见区域，即：DFOV区域作为最终结果输出。
Claims (18)
1.一种去除CT图像噪声的方法，包括如下步骤：
步骤一：估计图像的组织权重，包括：
a)配置初始参数：组织均值Tissue、软组织第一门限值TissueUp、软组织第二门限值TissueDown、组织降噪权重SRange、平滑处理窗宽度K；
b)在显示视野范围内，对输入的图像做平滑宽度为K的图像平滑处理，得到每个像素点的平滑值smoothimage(i，j)；
c)算出输入图像的每个像素点的组织权重：
<mrow>
<mi>T</mi>
<mi>i</mi>
<mi>s</mi>
<mi>s</mi>
<mi>u</mi>
<mi>e</mi>
<mi>W</mi>
<mi>e</mi>
<mi>i</mi>
<mi>g</mi>
<mi>h</mi>
<mi>t</mi>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mrow>
<mi>x</mi>
<mo>></mo>
<mi>T</mi>
<mi>i</mi>
<mi>s</mi>
<mi>s</mi>
<mi>u</mi>
<mi>e</mi>
<mi>U</mi>
<mi>p</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>-</mo>
<mi>T</mi>
<mi>i</mi>
<mi>s</mi>
<mi>s</mi>
<mi>u</mi>
<mi>e</mi>
<mi>D</mi>
<mi>o</mi>
<mi>w</mi>
<mi>n</mi>
<mo>)</mo>
<mo>*</mo>
<mi>K</mi>
<mi>K</mi>
<mo>+</mo>
<mi>S</mi>
<mi>R</mi>
<mi>a</mi>
<mi>n</mi>
<mi>g</mi>
<mi>e</mi>
</mrow>
</mtd>
<mtd>
<mrow>
<mi>x</mi>
<mo>&Element;</mo>
<mrow>
<mo>(</mo>
<mi>T</mi>
<mi>i</mi>
<mi>s</mi>
<mi>s</mi>
<mi>u</mi>
<mi>e</mi>
<mi>D</mi>
<mi>o</mi>
<mi>w</mi>
<mi>n</mi>
<mo>,</mo>
<mi>T</mi>
<mi>i</mi>
<mi>s</mi>
<mi>s</mi>
<mi>u</mi>
<mi>e</mi>
<mi>U</mi>
<mi>p</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>S</mi>
<mi>R</mi>
<mi>a</mi>
<mi>n</mi>
<mi>g</mi>
<mi>e</mi>
</mrow>
</mtd>
<mtd>
<mrow>
<mi>x</mi>
<mo><</mo>
<mi>T</mi>
<mi>i</mi>
<mi>s</mi>
<mi>s</mi>
<mi>u</mi>
<mi>e</mi>
<mi>D</mi>
<mi>o</mi>
<mi>w</mi>
<mi>n</mi>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
其中，
x＝abs(smoothimage(i，j)-Tissue)，abs()表示取绝对值；
步骤二：估计图像的噪声等级；
步骤三：计算去噪参数，所述去噪参数根据所述组织权重与所述噪声等级计算得到；
步骤四：利用所述去噪参数对图像进行各向异性扩散滤波；
步骤五：对滤波输出的图像进行边缘增强；
步骤六：对图像进行细节增强及对比度修正；
步骤七：裁切图像并输出结果。
2.根据权利要求1所述的一种去除CT图像噪声的方法，其特征是：所述步骤二包括如下步骤：
a)对输入的图像进行降采样；
b)选取估算区域：在降采样得到的图像上选取一块方形区域作为估算区域EstimationRegion；
c)计算估算区域的组织权重：对估算区域内的每个像素值，算出每个像素点的组织权重；
d)对估算区域加权：用组织权重对估算区域内的每一个像素点进行加权运算，得到加权后的图像像素值EstimationMap；
e)计算图像的内外包络：对加权后的图像像素值分别进行基于形态学的开运算和闭运算，分别得到图像的内包络和外包络；
f)计算噪声分布图：对内、外包络求平均值得到平滑分布图SmoothMap，并计算出噪声分布图
NoiseMap＝abs(EstimationMap-SmoothMap)，abs()表示求绝对值；
g)直方图分析：对该噪声分布图做直方图分析，得到噪声等级估算值NoiseLevel。
3.根据权利要求2所述的一种去除CT图像噪声的方法，其特征是：所述估算区域位于降采样得到的图像的中心位置，为正方形区域。
4.根据权利要求2所述的一种去除CT图像噪声的方法，其特征是：所述基于形态学的开运算依次包括计算滤波核、对所述加权后的图像像素值EstimationMap进行腐蚀操作和对腐蚀操作得到的结果进行膨胀操作的步骤。
5.根据权利要求4所述的一种去除CT图像噪声的方法，其特征是：所述的计算滤波核进一步包括如下步骤：
a)求出所述加权后的图像像素值的最大值MaxValue和平均值MeanValue；
b)设定滤波核大小KernelSize、滤波核范围Range以及滤波参数Delta；
c)计算出每个像素点对应的滤波核初始值：
其中，(i，j)为像素点坐标值，FWSize＝(KernelSize+1)/2；
d)计算每个像素点对应的滤波核的准确值：
Kernel(i，j)＝Kernel_init(i，j)*max(Kernel_init(i，j))*TH，其中，
TH＝(MaxValue-MeanValue)/Range，max(Kernel_init(i，j))指的是针对所有像素点遍历像素点坐标值后取的最大值。
6.根据权利要求2所述的一种去除CT图像噪声的方法，其特征是：所述基于形态学的闭运算依次包括计算滤波核、对所述加权后的图像像素值进行膨胀操作和对膨胀操作的结果进行腐蚀操作的步骤。
7.根据权利要求6所述的一种去除CT图像噪声的方法，其特征是：所述的计算滤波核按照如下步骤：
a)求出所述加权后的图像像素值的最大值MaxValue和平均值MeanValue；
b)设定滤波核大小KernelSize、滤波核范围Range以及滤波参数Delta；
c)计算出每个像素点对应的滤波核初始值：
其中，(i，j)为像素点坐标值，FWSize＝(KernelSize+1)/2；
d)计算每个像素点对应的滤波核的准确值：
Kernel(i，j)＝Kernel_init(i，j)*max(Kernel_init(i，j))*TH，其中，
TH＝(MaxValue-MeanValue)/Range，max(Kernel_init(i，j))指的是针对所有像素点遍历像素点坐标值后取的最大值。
8.根据权利要求2所述的一种去除CT图像噪声的方法，其特征是：所述的直方图分析包括如下步骤：
a)计算出所述噪声分布图NoiseMap中的像素最大值MaxValue1和最小值MinValue，进而得到图像的动态范围Range1＝MaxValue1-MinValue；
b)将图像的动态范围分成N份小区间，每一份小区间的步长Step＝Range1/N；
c)统计每个小区间内的像素点数Hist(t)，对于第t个区间而言，就是对[MinValue+Step*(t-1)，Minvalue+Step*t]区间内的像素点数进行统计；
d)对Hist(t)进行概率分布统计，即计算出公式中t的值，M表示噪声分布图NoiseMap中有效噪声的比例，所述M为是经过试验经验得到的估计值；
e)计算噪声分布图NoiseMap中有效噪声点的最大值：
Threshold＝(I-1)*Step+MinValue，其中的I表示上一步算出的t的值；
f)对噪声分布图NoiseMap中所有小于有效噪声点的最大值Threshold的值求平均值，得到图像的噪声等级NoiseLevel＝mean(NoiseMap＜Threshold)。
9.根据权利要求8所述的一种去除CT图像噪声的方法，其特征是：所述的噪声分布图NoiseMap中有效噪声的比例M＝0.95。
10.根据权利要求1所述的一种去除CT图像噪声的方法，其特征是：所述步骤三包括如下步骤：
a)设定参数Q和fwidth，参数Q和fwidth是根据CT图像扫描场景预先配置好的；
b)算出每个像素点位置校正值(x，y)表示把图像的中心点作为坐标原点时，每个像素点的坐标；
c)根据所述步骤一得到的组织权重、步骤二得到的噪声等级、以及所述每个像素点位置校正值，计算出每个像素点的去噪参数：
11.根据权利要求1所述的一种去除CT图像噪声的方法，其特征是：所述步骤四依次包括如下步骤：
a)计算出应当被去除的噪声值：
<mrow>
<mi>N</mi>
<mi>o</mi>
<mi>i</mi>
<mi>s</mi>
<mi>e</mi>
<mi>Re</mi>
<mi>d</mi>
<mi>u</mi>
<mi>c</mi>
<mi>e</mi>
<mo>=</mo>
<msubsup>
<mi>&Sigma;</mi>
<mrow>
<mi>p</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mn>8</mn>
</msubsup>
<mi>D</mi>
<mi>i</mi>
<mi>f</mi>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mi>p</mi>
<mo>)</mo>
</mrow>
<mo>*</mo>
<msup>
<mi>e</mi>
<mrow>
<mo>-</mo>
<mfrac>
<mrow>
<mi>D</mi>
<mi>i</mi>
<mi>f</mi>
<mi>f</mi>
<msup>
<mrow>
<mo>(</mo>
<mi>p</mi>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
<mrow>
<msup>
<mi>NoiseScale</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>.</mo>
<mo>*</mo>
<mi>q</mi>
</mrow>
</msup>
<mo>;</mo>
</mrow>
其中，
Diff(p)＝ds(p)*InputImage(x+m(p)，y+n(p))-InputImage(x，y)
m＝[1，2，3，1，3，1，2，3]，n＝[1，1，1，2，2，3，3，3]，ds＝[r，1，r，1，1，r，1，r]，
InputImage(x，y)为输入图像的像素值，.*表示点乘运算，p对应m、n、ds的第p个序列值，q为各向异性扩散滤波系数，NoiseScale为所述步骤三得到的去噪参数；
b)计算并输出去噪以后的图像像素值：
DeNoiseImage＝InputImage-NoiseReduce。
12.根据权利要求11所述的一种去除CT图像噪声的方法，其特征是：所述b)步完成之后，还包含如下步骤：将b)步得到的去噪以后的图像像素值DeNoiseImage赋值给InputImage，更换不同的各向异性扩散滤波系数值q，然后返回a)步，在至少依此方法循环执行三次后，停止所有步骤并输出去噪后的图像像素值DeNoiseImage。
13.根据权利要求1所述的一种去除CT图像噪声的方法，其特征是：所述步骤五包括如下步骤：
a)去除大颗粒噪声：用带通滤波器对图像Image进行带通滤波，得到图像中较大颗粒的噪声Noise以及去除该噪声以后的图像Filtered_Image＝Image-Noise；
b)提取图像边缘：更改带通滤波器的滤波频段，对去除噪声以后的图像进行边缘提取Edge＝BandPass(Filtered_Image)；
c)计算有效边缘：对边缘进行直方图分析，得到有效的边缘UsefulEdge；
d)边缘增强：边缘增强以后的图像EnhancedImage＝Filtered_Image+Scale*UsefulEdge，Scale表示对边缘进行增强的倍数。
14.根据权利要求1所述的一种去除CT图像噪声的方法，其特征是：所述步骤六包括如下步骤：
a)计算图像的第一层的背景信息和细节信息：对图像做哈尔小波变换，得到：[Background1，Foreground1]＝HarrTransform(Image)，其中Background1表示背景信息，Foreground1表示细节信息；
b)去除背景信息中的小噪声：对提取到的背景信息进行各向异性滤波，得到去除小噪声后新的背景信息BackgroundNew1＝AnisotropicDiffusionFilter(Background1)；
c)计算新的背景信息中的第二层细节信息和背景信息：对新的背景信息BackgroundNew1做哈尔小波变换，得到第二层的背景信息和细节信息：[Background2，Foreground2]＝HarrTransform(BackgroundNew1)，其中Background2表示背景信息，Foreground2表示细节信息；
d)去除第二层背景信息中的小噪声：对第二层的背景信息进行各向异性滤波，得到去除小噪声后新的第二层背景信息：
BackgroundNew2＝AnisotropicDiffusionFilter(Background2)；
e)对第二层的图像细节进行尺度为K2的增强：增强后的第二层图像，
Enhance2＝K2*Foreground2+BackgroundNew2；
f)对增强后的第二层图像做哈尔小波反变换：
输出Out2＝HarrInverse(Enhance2)
g)计算第一层图像细节增强的对比度修正因子：
输出
h)对第一层的图像细节进行尺度为K1的增强及对比度修正：
增强后的第一层图像：
Enhance1＝Correctfactor*K1*Foreground1+Background1
i)对增强后的第一层图像做哈尔小波反变换：其输出Out1＝HarrInverse(Enhance1)。
15.根据权利要求1所述的一种去除CT图像噪声的方法，其特征是：所述步骤七中的裁切图像是指截取图像中的显示视野可见区域，即：DFOV区域。
16.一种去除CT图像噪声的装置，包括：用于估计图像的组织权重的模块，用于估计图像的噪声等级的模块，该模块包括用于对输入的图像进行降采样的模块；用于选取估算区域的模块，其在降采样得到的图像上选取一块方形区域作为估算区域EstimationRegion；用于估算区域的组织权重的模块，其对估算区域内的每个像素值用所述估计图像的组织权重的方法估算其组织权重；用于对估算区域加权的模块，其用组织权重对估算区域内的每一个像素点进行加权运算，得到加权后的图像像素值EstimationMap；用于计算图像的内外包络的模块，其对加权后的图像像素值分别进行基于形态学的开运算和闭运算，分别得到图像的内包络和外包络；用于计算噪声分布图的模块，其对内、外包络求平均值得到平滑分布图SmoothMap，并计算出噪声分布图NoiseMap＝abs(EstimationMap-SmoothMap)，abs()表示求绝对值；用于直方图分析的模块，其对该噪声分布图做直方图分析，得到噪声等级估算值NoiseLevel，根据所述组织权重与所述噪声等级计算去噪参数的模块，利用所述去噪参数对图像进行各向异性扩散滤波的模块，用于对滤波输出的图像进行边缘增强的模块，用于对图像进行细节增强及对比度修正的模块，用于裁切图像并输出结果的模块。
17.根据权利要求16所述的一种去除CT图像噪声的装置，所述的用于对滤波输出的图像进行边缘增强的模块进一步包括：用于去除大颗粒噪声的模块，其用带通滤波器对图像Image进行带通滤波，得到图像中较大颗粒的噪声Noise以及去除该噪声以后的图像Filtered_Image＝Image-Noise；用于提取图像边缘的模块，其更改带通滤波器的滤波频段，对去除噪声以后的图像进行边缘提取Edge＝BandPass(Filtered_Image)；用于计算有效边缘的模块，其对边缘进行直方图分析，得到有效的边缘UsefulEdge；用于边缘增强的模块，其边缘增强以后的图像EnhancedImage＝Filtered_Image+Scale*UsefulEdge，Scale表示对边缘进行增强的倍数。
18.根据权利要求16所述的一种去除CT图像噪声的装置，所述的用于对图像进行细节增强及对比度修正的模块进一步包括用于计算图像的第一层的背景信息和细节信息的模块，其对图像做哈尔小波变换，得到：[Background1，Foreground1]＝HarrTransform(Image)，其中Background1表示背景信息，Foreground1表示细节信息；用于去除背景信息中的小噪声的模块，其对提取到的背景信息进行各向异性滤波，得到去除小噪声后新的背景信息BackgroundNew1＝AnisotropicDiffusionFilter(Background1)；用于计算新的背景信息中的第二层细节信息和背景信息的模块，其对新的背景信息BackgroundNew1做哈尔小波变换，得到第二层的背景信息和细节信息：[Background2，Foreground2]＝HarrTransform(BackgroundNew1)，其中Background2表示背景信息，Foreground2表示细节信息；用于去除第二层背景信息中的小噪声的模块，其对第二层的背景信息进行各向异性滤波，得到去除小噪声后新的第二层背景信息BackgroundNew2＝AnisotropicDiffusionFilter(Background2)；用于对第二层的图像细节进行尺度为K2的增强的模块，其增强后的第二层图像Enhance2＝K2*Foreground2+BackgroundNew2；用于对增强后的第二层图像做哈尔小波反变换的模块，其输出Out2＝HarrInverse(Enhance2)；用于计算第一层图像细节增强的对比度修正因子的模块，其输出用于对第一层的图像细节进行尺度为K1的增强及对比度修正的模块，其增强后的第一层图像：
Enhance1＝Correctfactor*K1*Foreground1+Background1；用于对增强后的第一层图像做哈尔小波反变换的模块，其输出Out1＝HarrInverse(Enhance1)。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201110462903.8A CN103186888B (zh) | 2011-12-30 | 2011-12-30 | 一种去除ct图像噪声的方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201110462903.8A CN103186888B (zh) | 2011-12-30 | 2011-12-30 | 一种去除ct图像噪声的方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN103186888A CN103186888A (zh) | 2013-07-03 |
| CN103186888B true CN103186888B (zh) | 2017-11-21 |
Family
ID=48678046
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201110462903.8A Active CN103186888B (zh) | 2011-12-30 | 2011-12-30 | 一种去除ct图像噪声的方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN103186888B (zh) |
Families Citing this family (18)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104036457B (zh) * | 2013-09-02 | 2015-03-18 | 上海联影医疗科技有限公司 | 一种图像降噪方法 |
| DE102014206720A1 (de) * | 2014-04-08 | 2015-10-08 | Siemens Aktiengesellschaft | Rauschreduktion in Tomogrammen |
| US10169848B2 (en) * | 2014-04-23 | 2019-01-01 | Koninklijke Philips N.V. | Restoration of low contrast structure in de-noise image data |
| CN104331863B (zh) * | 2014-10-15 | 2018-01-09 | 中国航空工业集团公司洛阳电光设备研究所 | 一种图像滤波去噪方法 |
| CN104376564B (zh) * | 2014-11-24 | 2018-04-24 | 西安工程大学 | 基于各向异性高斯方向导数滤波器提取图像粗边缘的方法 |
| CN105812621A (zh) * | 2014-12-31 | 2016-07-27 | 乐视致新电子科技（天津）有限公司 | 自适应调节图像噪声水平方法以及装置 |
| CN106157253B (zh) * | 2015-04-17 | 2019-09-03 | 瑞昱半导体股份有限公司 | 图像处理装置与图像处理方法 |
| CN106157275A (zh) * | 2015-04-21 | 2016-11-23 | 杭州与盟医疗技术有限公司 | 基于热扩散方程的多尺度ct图像去噪算法 |
| CN106296633B (zh) * | 2015-05-22 | 2019-08-02 | 浙江大学 | 一种基于多尺度图像域双能分解算法 |
| CN106994021B (zh) * | 2016-01-22 | 2022-10-14 | 通用电气公司 | 一种计算ct影像上的噪声的方法及装置 |
| EP3503528B1 (en) * | 2017-12-21 | 2020-10-07 | Axis AB | Determination of a contrast value for a digital image |
| CN108810320B (zh) * | 2018-06-01 | 2020-11-24 | 深圳市商汤科技有限公司 | 图像质量提升方法和装置 |
| CN110514672A (zh) * | 2018-09-06 | 2019-11-29 | 永康市缘匠贸易有限公司 | 玻璃噪声自适应去除系统 |
| CN109978782A (zh) * | 2019-03-18 | 2019-07-05 | 新疆大学 | 一种数字全息再现图像中散斑噪声的去除方法及装置 |
| CN110796615B (zh) * | 2019-10-18 | 2023-06-02 | 浙江大华技术股份有限公司 | 一种图像去噪方法、装置以及存储介质 |
| CN111353958B (zh) * | 2020-02-28 | 2023-07-25 | 北京东软医疗设备有限公司 | 图像处理方法、装置及系统 |
| CN112508823B (zh) * | 2021-02-02 | 2021-09-07 | 数坤(北京)网络科技股份有限公司 | 一种影像降噪方法、装置以及计算机可读存储介质 |
| CN118014884B (zh) * | 2024-04-09 | 2024-08-06 | 福州海狸家口腔科技有限公司 | 一种牙冠三维模型的边缘去噪方法、介质和设备 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1614636A (zh) * | 2003-11-06 | 2005-05-11 | Ge医疗系统环球技术有限公司 | 图像重构方法和x射线ct装置 |
| CN101320468A (zh) * | 2007-06-07 | 2008-12-10 | 株式会社东芝 | 数据处理装置、医用诊断装置、数据处理方法以及医用诊断方法 |
-
2011
- 2011-12-30 CN CN201110462903.8A patent/CN103186888B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1614636A (zh) * | 2003-11-06 | 2005-05-11 | Ge医疗系统环球技术有限公司 | 图像重构方法和x射线ct装置 |
| CN101320468A (zh) * | 2007-06-07 | 2008-12-10 | 株式会社东芝 | 数据处理装置、医用诊断装置、数据处理方法以及医用诊断方法 |
Non-Patent Citations (5)
| Title |
|---|
| Contrast Enhancement of Medical Images Using Multiscale Edge Representation;Jian Lu et al.;《Optical engineering》;19940701;第33卷(第7期);第2151-2161页 * |
| Nonlinear Anisotropic Filtering of MRI Data;Guido Gerig et al.;《IEEE TRANSACTIONS ON MEDICAL IMAGING》;19920630;第11卷(第2期);第221-232页 * |
| 医学CT图像的三维可视化技术研究;周宴;《中国优秀硕士学位论文全文数据库》;20051215(第8期);第2.4.2-2.4.3节，第3.1节 * |
| 医学CT图像的三维可视化技术研究;周宴;《中国优秀硕士学位论文全文数据库》;20051215(第8期);第C030-21页 * |
| 基于边缘增强算子的各向异性滤波方法;朱丽琪 等;《电脑知识与技术》;20110930;第7卷(第25期);第6214-6217页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN103186888A (zh) | 2013-07-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN103186888B (zh) | 一种去除ct图像噪声的方法及装置 | |
| US8542944B2 (en) | Method and apparatus for multi-scale based dynamic range compression and noise suppression | |
| JP5851301B2 (ja) | 画像処理装置、画像処理方法、及び磁気共鳴イメージング装置 | |
| EP1909227B1 (en) | Method of and apparatus for minimizing ringing artifacts in an input image | |
| CN102737382B (zh) | 一种前列腺超声图像自动精确分割方法 | |
| CN101499164B (zh) | 一种基于单张低分辨率图像的图像插值重构方法 | |
| EP2588374A2 (en) | Non-linear resolution reduction for medical imagery | |
| EP3072104B1 (en) | Image de-noising method | |
| Jeevakala | Sharpening enhancement technique for MR images to enhance the segmentation | |
| CN104077743A (zh) | 一种x射线图像中滤线栅伪影的抑制方法及装置 | |
| CN107680057A (zh) | 超声图像增强的方法及装置 | |
| CN107437251A (zh) | 头部mri影像颅骨剥离模块 | |
| CN105809633A (zh) | 去除颜色噪声的方法及装置 | |
| CN111429362A (zh) | 一种内窥镜荧光图像的血管增强方法 | |
| CN115375574A (zh) | 基于区域自适应的多尺度非局部低剂量ct图像去噪方法 | |
| CN114708352A (zh) | 一种基于无监督学习的低剂量ct图像恢复方法及系统 | |
| Boudraa et al. | Image contrast enhancement based on 2D Teager-Kaiser operator | |
| CN110443760B (zh) | 图像降噪方法及计算机可读存储介质 | |
| CN104504659A (zh) | 一种基于提升小波变换的快速iso去噪方法及系统 | |
| Nugroho et al. | Performance analysis of filtering techniques for speckle reduction on breast ultrasound images | |
| EP2328124B1 (en) | Method of enhancing the contrast of spatially-localized phenomena in an image | |
| Kazmi et al. | Medical image denoising based on adaptive thresholding in contourlet domain | |
| CN112785520A (zh) | Ct图像伪影去除方法 | |
| Snoeren et al. | Thickness correction of mammographic images by anisotropic filtering and interpolation of dense tissue | |
| Pancholi et al. | A Review of Noise Reduction Filtering Techniques for MRI Images |
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 |