微网站和普通网站区别,网站建设服务后所有权归谁,找网络公司建网站每年收维护费,广州汽车网站建设专栏地址#xff1a;『youcans 的图像处理学习课』 文章目录#xff1a;『youcans 的图像处理学习课 - 总目录』 【youcans 的 OpenCV 学习课】10. 图像复原与重建
图像复原是对图像退化过程建模#xff0c;并以图像退化的先验知识来恢复退化的图像。
图像增强是一种主观处… 专栏地址『youcans 的图像处理学习课』 文章目录『youcans 的图像处理学习课 - 总目录』 【youcans 的 OpenCV 学习课】10. 图像复原与重建
图像复原是对图像退化过程建模并以图像退化的先验知识来恢复退化的图像。
图像增强是一种主观处理而图像复原是一种客观处理。
本文提供上述各种算法的完整例程和运行结果。 文章目录【youcans 的 OpenCV 学习课】10. 图像复原与重建1. 图像退化/复原处理模型1.1 图像退化模型2. 噪声模型2.1 高斯噪声Gauss Noise例程 9.1高斯噪声Gauss Noise2.2 瑞利噪声 Rayleigh Noise例程 9.2瑞利噪声 Rayleigh Noise2.3 爱尔兰噪声 Ireland Noise例程 9.3伽马噪声 Gamma Noise2.4 指数噪声 Exponential Noise例程 9.4指数噪声 Exponent Noise2.5 均匀噪声 Uniform Noise例程 9.5均匀噪声 Uniform Noise2.6 椒盐噪声 Salt-pepper Noise例程 9.6椒盐噪声 Salt-pepper Noise2.7 图像噪声小结例程 9.7各种噪声模型的直方图分布在这里插入图片描述3. 仅噪声存在的空间滤波图像复原3.1 算术平均滤波器Arithmentic mean filter例程 9.8算术平均滤波器3.2 几何均值滤波器Geometric mean filter例程 9.9几何均值滤波器3.3 谐波平均滤波器Harmonic mean filter例程 9.10谐波平均滤波器3.4 反谐波平均滤波器Invert-harmonic mean filter例程 9.11反谐波平均滤波器3.5 统计排序滤波器例程 9.12统计排序滤波器3.6 修正阿尔法均值滤波器Modified alpha-mean filter例程 9.13修正阿尔法均值滤波器3.7 自适应局部降噪滤波器例程 9.14自适应局部降噪滤波器3.8 自适应中值滤波器Adaptive median filter例程 9.15自适应中值滤波器4. 频率域滤波图像复原4.1 陷波滤波器**Notch Filter**例程 9.16陷波带阻滤波器的传递函数例程 9.17陷波带阻滤波器消除周期噪声干扰4.2 最优陷波滤波器5. 估计退化函数5.1 观察法估计退化函数5.2 试验法估计退化函数5.3 模型法估计退化函数例程 9.18运动模糊退化模型例程 9.19湍流模糊退化模型6. 退化图像复原6.1 退化图像的逆滤波Inverse filter例程 9.20湍流模糊退化图像的逆滤波6.2 退化图像的维纳滤波Wiener filter例程 9.21: 维纳滤波 Wiener filter6.3 约束最小二乘方滤波Constrained Least Squares Filtering例程 9.21: 约束最小二乘方滤波6.4 几何均值滤波例程 9.22: 几何均值滤波7. 投影重建图像7.1 计算机断层成像Computed tomography7.2 投影和雷登变换Radon transform例程 9.23雷登变换正弦图7.2 雷登变换反投影重建图像Radon transform back projection例程 9.24雷登变换反投影重建图像7.3 滤波反投影重建图像 Filtered back projection例程 9.25滤波反投影重建图像1. 图像退化/复原处理模型
1.1 图像退化模型
图像退化
图像退化是指图像在形成、记录、处理和传输过程中由于 成像系统、记录设备、传输介质和处理方法的不完善导致图像质量下降。
把图像退化表示为退化算子 H\mathcal{H}H退化算子与一个加性噪声项 η(x,y)\eta (x,y)η(x,y) 共同对输入图像 f(x,y)f(x,y)f(x,y) 进行运算生成退化图像 g(x,y)g(x,y)g(x,y)。
图像复原
图像复原是指利用退化过程的先验知识消除或降低在图像获取、传输过程中造成的图像品质下降恢复图像的本来面目。具体地要根据退化的原因分析引起退化的环境因素建立相应的数学模型并沿着使图像降质的逆过程恢复图像。因此图象复原是一个求逆问题就是把退化模型化采用相反的过程处理以复原出原始图像。
已知退化图像 g(x,y)g(x,y)g(x,y)关于退化算子 H\mathcal{H}H 与加性噪声项 η(x,y)\eta (x,y)η(x,y) 的概括信息通过图像复原处理得到输入图像 f(x,y)f(x,y)f(x,y) 的估计 f^(x,y)\hat{f}(x,y)f^(x,y) 尽可能接近原图像。
一般来说关于退化算子 H\mathcal{H}H 与加性噪声项 η(x,y)\eta (x,y)η(x,y) 的信息越多得到的估计 f^(x,y)\hat{f}(x,y)f^(x,y) 就越接近原图像 f(x,y)f(x,y)f(x,y) 。
由于引起退化的因素众多、性质不同为了描述图像退化过程所建立的数学模型多种多样而恢复的质量标准也存在差异性因此图像复原是一个复杂的数学过程相关的方法和技术也各不相同。
退化模型
若退化算子 H\mathcal{H}H 是一个线性位置不变算子则空间域的退化图像为 g(x,y)(h⋆f)(x,y)η(x,y)g(x,y) (h \star f)(x,y) \eta(x,y) g(x,y)(h⋆f)(x,y)η(x,y) 频率域的等效公式为 G(u,v)H(u,v)F(u,v)N(u,v)G(u,v) H(u,v) F(u,v) N(u,v) G(u,v)H(u,v)F(u,v)N(u,v) 2. 噪声模型
数字图像中的噪声源主要来自图像获取和传输过程。在获取图像时光照水平和传感器温度影响图像中的噪声。在传输图像时传输信道中的干扰对图像产生污染。
白噪声
白光等比例地包含可见光谱中的所有频率。当噪声的傅里叶谱是常量时称为白噪声。
2.1 高斯噪声Gauss Noise
高斯噪声中空间域和频率域中都很方便进行数学处理因而得到了广泛的应用。
高斯噪声的概率密度函数为 p(z)12πσe−(z−zˉ)2/2σ2p(z) \frac{1}{\sqrt{2\pi} \sigma} e^{-(z-\bar{z})^2/2\sigma ^2} p(z)2πσ1e−(z−zˉ)2/2σ2
高斯噪声的均值为 zˉ\bar{z}zˉ标准差为σ2\sigma ^2σ2。 例程 9.1高斯噪声Gauss Noise # 9.1高斯噪声 (GaussNoise)img cv2.imread(../images/Fig0503.tif, 0) # flags0 读取为灰度图像# img np.ones([256, 256]) * 128mu, sigma 0.0, 20.0noiseGause np.random.normal(mu, sigma, img.shape)imgGaussNoise img noiseGauseimgGaussNoise np.uint8(cv2.normalize(imgGaussNoise, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]plt.figure(figsize(9, 3))plt.subplot(131), plt.title(Origin), plt.axis(off)plt.imshow(img, gray, vmin0, vmax255)plt.subplot(132), plt.title(GaussNoise), plt.axis(off)plt.imshow(imgGaussNoise, gray)plt.subplot(133), plt.title(Gray Hist)histNP, bins np.histogram(imgGaussNoise.flatten(), bins255, range[0, 255], densityTrue)plt.bar(bins[:-1], histNP[:])plt.tight_layout()plt.show()2.2 瑞利噪声 Rayleigh Noise
瑞利噪声的概率密度函数为 p(z){2/b∗(z−a)e−(z−a)2/b,z≥a0,zap(z) \begin{cases} 2/b * (z-a) e^{-(z-a)^2 /b} , z \ge a\\ 0, z a \end{cases} p(z){2/b∗(z−a)e−(z−a)2/b0,z≥a,za
瑞利噪声的均值和标准差为 zˉaπb/4σ2b(4−π)/4\bar{z} a \sqrt{\pi b/4} \\ \sigma ^2 b(4-\pi)/4 zˉaπb/4σ2b(4−π)/4
瑞利噪声概率密度分布到原点的距离及密度的基本形状右偏常用于倾斜形状直方图的建模。 例程 9.2瑞利噪声 Rayleigh Noise # # 9.2瑞利噪声 (RayleighNoise)img cv2.imread(../images/Fig0503.tif, 0) # flags0 读取为灰度图像# img np.ones([256, 256]) * 128a 30.0noiseRayleigh np.random.rayleigh(a, sizeimg.shape)imgRayleighNoise img noiseRayleighimgRayleighNoise np.uint8(cv2.normalize(imgRayleighNoise, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]plt.figure(figsize(9, 3))plt.subplot(131), plt.title(Origin), plt.axis(off)plt.imshow(img, gray, vmin0, vmax255)plt.subplot(132), plt.title(RayleighNoise), plt.axis(off)plt.imshow(imgRayleighNoise, gray)plt.subplot(133), plt.title(Gray Hist)histNP, bins np.histogram(imgRayleighNoise.flatten(), bins255, range[0, 255], densityTrue)plt.bar(bins[:-1], histNP[:])plt.tight_layout()plt.show()2.3 爱尔兰噪声 Ireland Noise
爱尔兰噪声的概率密度函数为 p(z){abzb−1(b−1)!e−az,z≥00,z0p(z) \begin{cases} \frac{a^b z^{b-1}}{(b-1)!} e^{-az} , z \ge 0\\ 0, z 0 \end{cases} p(z){(b−1)!abzb−1e−az0,z≥0,z0
爱尔兰噪声的均值和标准差为 zˉb/aσ2b/a2\bar{z} b/a \\ \sigma ^2 b/a^2 zˉb/aσ2b/a2
当标准差的分母 a2a^2a2为伽马函数时称为伽马噪声。 例程 9.3伽马噪声 Gamma Noise # # 9.3伽马噪声 (Gamma Noise)img cv2.imread(../images/Fig0503.tif, 0) # flags0 读取为灰度图像# img np.ones([256, 256]) * 128a, b 10.0, 2.5noiseGamma np.random.gamma(shapeb, scalea, sizeimg.shape)imgGammaNoise img noiseGammaimgGammaNoise np.uint8(cv2.normalize(imgGammaNoise, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]plt.figure(figsize(9, 3))plt.subplot(131), plt.title(Origin), plt.axis(off)plt.imshow(img, gray, vmin0, vmax255)plt.subplot(132), plt.title(Gamma noise), plt.axis(off)plt.imshow(imgGammaNoise, gray)plt.subplot(133), plt.title(Gray hist)histNP, bins np.histogram(imgGammaNoise.flatten(), bins255, range[0, 255], densityTrue)plt.bar(bins[:-1], histNP[:])plt.tight_layout()plt.show()2.4 指数噪声 Exponential Noise
指数噪声的概率密度函数为 p(z){ae−az,z≥00,z0p(z) \begin{cases} ae^{-az} , z \ge 0\\ 0, z 0 \end{cases} p(z){ae−az0,z≥0,z0
指数噪声的均值和标准差为 zˉ1/aσ21/a2\bar{z} 1/a \\ \sigma ^2 1/a^2 zˉ1/aσ21/a2
显然指数噪声是爱尔兰噪声在 b1 时的特殊情况。 例程 9.4指数噪声 Exponent Noise # # 9.4指数噪声 (Exponential noise)img cv2.imread(../images/Fig0503.tif, 0) # flags0 读取为灰度图像# img np.ones([256, 256]) * 128a 10.0noiseExponent np.random.exponential(scalea, sizeimg.shape)imgExponentNoise img noiseExponentimgExponentNoise np.uint8(cv2.normalize(imgExponentNoise, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]plt.figure(figsize(9, 3))plt.subplot(131), plt.title(Origin), plt.axis(off)plt.imshow(img, gray, vmin0, vmax255)plt.subplot(132), plt.title(Exponential noise), plt.axis(off)plt.imshow(imgExponentNoise, gray)plt.subplot(133), plt.title(Gray hist)histNP, bins np.histogram(imgExponentNoise.flatten(), bins255, range[0, 255], densityTrue)plt.bar(bins[:-1], histNP[:])plt.tight_layout()plt.show()2.5 均匀噪声 Uniform Noise
均匀噪声的概率密度函数为 p(z){1/(b−a),a≤z≤b0,otherp(z) \begin{cases} 1/(b-a) , a \le z \le b\\ 0, other \end{cases} p(z){1/(b−a)0,a≤z≤b,other
均匀噪声的均值和标准差为 zˉ(ab)/2σ2(b−a)2/12\bar{z} (ab)/2 \\ \sigma ^2 (b-a)^2/12 zˉ(ab)/2σ2(b−a)2/12 例程 9.5均匀噪声 Uniform Noise # # 9.5均匀噪声 (Uniform noise)img cv2.imread(../images/Fig0503.tif, 0) # flags0 读取为灰度图像# img np.ones([256, 256]) * 128mean, sigma 10, 100a 2 * mean - np.sqrt(12 * sigma) # a -14.64b 2 * mean np.sqrt(12 * sigma) # b 54.64noiseUniform np.random.uniform(a, b, img.shape)imgUniformNoise img noiseUniformimgUniformNoise np.uint8(cv2.normalize(imgUniformNoise, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]plt.figure(figsize(9, 3))plt.subplot(131), plt.title(Origin), plt.axis(off)plt.imshow(img, gray, vmin0, vmax255)plt.subplot(132), plt.title(Uniform noise), plt.axis(off)plt.imshow(imgUniformNoise, gray)plt.subplot(133), plt.title(Gray hist)histNP, bins np.histogram(imgUniformNoise.flatten(), bins255, range[0, 255], densityTrue)plt.bar(bins[:-1], histNP[:])plt.tight_layout()plt.show()2.6 椒盐噪声 Salt-pepper Noise
椒盐噪声的概率密度函数为 p(z){Ps,z2k−1Pp,z01−(PsPp),zVp(z) \begin{cases} Ps , z 2^k -1 \\ Pp , z 0 \\ 1-(PsPp), z V \end{cases} p(z)⎩⎪⎨⎪⎧PsPp1−(PsPp),z2k−1,z0,zV
当 Ps、Pp 都不为 0 时噪声值是白色的 (2k−1)(2^k-1)(2k−1) 或黑色的 (0)(0)(0)就像盐粒或胡椒粒那样随机地分布在整个图像上因此称为椒盐噪声也称为双极冲击噪声。当 Ps 或 Pp 为 0 时称为单极冲击噪声。
椒盐噪声的均值和标准差为 zˉ(0)PpK(1−Ps−Pp)(2k−1)Psσ2(0−zˉ)2Pp(K−zˉ)2(1−Ps−Pp)(2k−1)2Ps\bar{z} (0)PpK(1-Ps-Pp)(2^k-1)Ps \\ \sigma ^2 (0-\bar{z})^2 Pp(K-\bar{z})^2(1-Ps-Pp)(2^k-1)^2Ps zˉ(0)PpK(1−Ps−Pp)(2k−1)Psσ2(0−zˉ)2Pp(K−zˉ)2(1−Ps−Pp)(2k−1)2Ps
像素被白色盐粒、黑色胡椒粒污染的概率 P 称为噪声密度 PPsPpP Ps Pp PPsPp 例如Ps0.05Pp0.02则噪声密度 P0.07表示图像中约 5% 的像素被盐粒噪声污染约 2% 的像素被胡椒粒噪声污染噪声密度为 7%即图像中 7% 的像素被椒盐噪声污染。 例程 9.6椒盐噪声 Salt-pepper Noise # # 9.6椒盐噪声 (Salt-pepper)img cv2.imread(../images/Fig0503.tif, 0) # flags0 读取为灰度图像ps, pp 0.05, 0.02mask np.random.choice((0, 0.5, 1), sizeimg.shape[:2], p[pp, (1-ps-pp), ps])imgChoiceNoise img.copy()imgChoiceNoise[mask1] 255imgChoiceNoise[mask0] 0plt.figure(figsize(9, 3))plt.subplot(131), plt.title(Origin), plt.axis(off)plt.imshow(img, gray, vmin0, vmax255)plt.subplot(132), plt.title(Choice noise), plt.axis(off)plt.imshow(imgChoiceNoise, gray)plt.subplot(133), plt.title(Gray hist)histNP, bins np.histogram(imgChoiceNoise.flatten(), bins255, range[0, 255], densityTrue)plt.bar(bins[:-1], histNP[:])plt.tight_layout()plt.show()2.7 图像噪声小结
上述噪声模型为各种图像噪声的建模提供了有效的工具。例如
高斯噪声常用于电子电路及传感器噪声光照不足和/或高温引起等因素所导致噪声的建模
瑞利噪声常用于距离成像中的噪声建模
指数噪声和伽马噪声常用于激光成像中的噪声建模
冲激噪声常用于在成像期间的快速瞬变如开关故障的噪声建模
均匀噪声是对实际情况最基本描述也是数字图像处理中各种随机数发生器的基础。 例程 9.7各种噪声模型的直方图分布 # # 9.7:各种噪声模型的直方图分布img cv2.imread(../images/Fig0503.tif, 0) # flags0 读取为灰度图像mu, sigma 0.0, 20.0noiseGause np.random.normal(mu, sigma, img.shape)imgGaussNoise img noiseGauseimgGaussNoise np.uint8(cv2.normalize(imgGaussNoise, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]a 30.0noiseRayleigh np.random.rayleigh(a, sizeimg.shape)imgRayleighNoise img noiseRayleighimgRayleighNoise np.uint8(cv2.normalize(imgRayleighNoise, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]a, b 10.0, 2.5noiseGamma np.random.gamma(shapeb, scalea, sizeimg.shape)imgGammaNoise img noiseGammaimgGammaNoise np.uint8(cv2.normalize(imgGammaNoise, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]a 10.0noiseExponent np.random.exponential(scalea, sizeimg.shape)imgExponentNoise img noiseExponentimgExponentNoise np.uint8(cv2.normalize(imgExponentNoise, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]mean, sigma 10, 100a 2 * mean - np.sqrt(12 * sigma) # a -14.64b 2 * mean np.sqrt(12 * sigma) # b 54.64noiseUniform np.random.uniform(a, b, img.shape)imgUniformNoise img noiseUniformimgUniformNoise np.uint8(cv2.normalize(imgUniformNoise, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]ps, pp 0.05, 0.02mask np.random.choice((0, 0.5, 1), sizeimg.shape[:2], p[pp, (1-ps-pp), ps])imgChoiceNoise img.copy()imgChoiceNoise[mask1] 255imgChoiceNoise[mask0] 0plt.figure(figsize(12, 8))plt.subplot(231), plt.title(Gauss noise)histNP, bins np.histogram(imgGaussNoise.flatten(), bins255, range[0, 255], densityTrue)plt.bar(bins[:-1], histNP[:])plt.subplot(232), plt.title(Rayleigh noise)histNP, bins np.histogram(imgRayleighNoise.flatten(), bins255, range[0, 255], densityTrue)plt.bar(bins[:-1], histNP[:])plt.subplot(233), plt.title(Gamma noise)histNP, bins np.histogram(imgGammaNoise.flatten(), bins255, range[0, 255], densityTrue)plt.bar(bins[:-1], histNP[:])plt.subplot(234), plt.title(Exponential noise)histNP, bins np.histogram(imgExponentNoise.flatten(), bins255, range[0, 255], densityTrue)plt.bar(bins[:-1], histNP[:])plt.subplot(235), plt.title(Uniform noise)histNP, bins np.histogram(imgUniformNoise.flatten(), bins255, range[0, 255], densityTrue)plt.bar(bins[:-1], histNP[:])plt.subplot(236), plt.title(Salt-pepper noise)histNP, bins np.histogram(imgChoiceNoise.flatten(), bins255, range[0, 255], densityTrue)plt.bar(bins[:-1], histNP[:])plt.tight_layout()plt.show()3. 仅噪声存在的空间滤波图像复原
当一幅图像中唯一存在的退化是噪声时退化模型简化为 g(x,y)f(x,y)η(x,y)G(u,v)F(u,v)N(u,v)g(x,y) f(x,y) \eta(x,y) \\ G(u,v) F(u,v) N(u,v) g(x,y)f(x,y)η(x,y)G(u,v)F(u,v)N(u,v) 当仅存在加性随机噪声时可以采用空间滤波方法来估计原图像 f(x,y)f(x,y)f(x,y)即对退化图像 g(x,y)g(x,y)g(x,y) 去除噪声。
空间滤波方法在《7. 空间域图像滤波》中进行了详细介绍本章简要讨论空间滤波的降噪性能。
3.1 算术平均滤波器Arithmentic mean filter
算术平均滤波器是最简单的均值滤波器与空间域滤波中的盒式滤波器相同。
令 SxySxySxy 表示中心在点 (x,y)(x,y)(x,y) 、大小为 m∗nm*nm∗n 的矩形子窗口邻域的一组坐标算术平均滤波器在由 SxySxySxy 定义的区域中计算被污染图像 g(x,y)g(x,y)g(x,y) 的平均值。
复原的图像 f^\hat{f}f^ 在点 (x,y)(x,y)(x,y) 的值是使用 SxySxySxy 定义的邻域中的像素算出的算术平均值即 f^(x,y)1mn∑(r,c)∈Sxyg(r,c)\hat{f}(x,y) \frac{1}{mn} \sum _{(r,c) \in Sxy} g(r,c) f^(x,y)mn1(r,c)∈Sxy∑g(r,c) 这一操作可以通过一个大小为 m∗nm*nm∗n 的空间滤波器核来实现核的所有系数都是 1/mn1/mn1/mn。
均值滤波平滑图像中的局部变化可以降低图像中的噪声但会模糊图像。 例程 9.8算术平均滤波器 # 9.8: 算术平均滤波器 (Arithmentic mean filter)# 参见例程 1.70图像的低通滤波 (盒式滤波器核)img cv2.imread(../images/Fig0507b.tif, 0) # flags0 读取为灰度图像kSize (3,3)kernalMean np.ones(kSize, np.float32) / (kSize[0]*kSize[1]) # 生成归一化盒式核imgConv1 cv2.filter2D(img, -1, kernalMean) # cv2.filter2D 方法kSize (5,5)imgConv3 cv2.boxFilter(img, -1, kSize) # cv2.boxFilter 方法 (默认normalizeTrue)plt.figure(figsize(9, 6))plt.subplot(131), plt.axis(off), plt.title(Original)plt.imshow(img, cmapgray, vmin0, vmax255)plt.subplot(132), plt.axis(off), plt.title(filter2D(kSize[3,3]))plt.imshow(imgConv1, cmapgray, vmin0, vmax255)plt.subplot(133), plt.axis(off), plt.title(boxFilter(kSize[5,5]))plt.imshow(imgConv3, cmapgray, vmin0, vmax255)plt.tight_layout()plt.show()3.2 几何均值滤波器Geometric mean filter
使用几何均值滤波器复原图像复原图像 f^\hat{f}f^ 在点 (x,y)(x,y)(x,y) 的值是邻域中的像素的几何平均值 f^(x,y)[∏(r,c)∈Sxyg(r,c)]1/mn\hat{f}(x,y) \begin{bmatrix} \prod _{(r,c) \in Sxy} g(r,c) \end{bmatrix} ^{1/mn} f^(x,y)[∏(r,c)∈Sxyg(r,c)]1/mn 几何均值滤波器实现的平滑与算术平均滤波器相当但损失的图像细节更少。 例程 9.9几何均值滤波器 # 9.9: 几何均值滤波器 (Geometric mean filter)img cv2.imread(../images/Fig0507b.tif, 0) # flags0 读取为灰度图像img_h img.shape[0]img_w img.shape[1]# 算术平均滤波 (Arithmentic mean filter)kSize (3,3)kernalMean np.ones(kSize, np.float32) / (kSize[0]*kSize[1]) # 生成归一化盒式核imgAriMean cv2.filter2D(img, -1, kernalMean)# 几何均值滤波器 (Geometric mean filter)m, n 3, 3order 1/(m*n)kernalMean np.ones((m,n), np.float32) # 生成盒式核hPad int((m-1) / 2)wPad int((n-1) / 2)imgPad np.pad(img.copy(), ((hPad, m-hPad-1), (wPad, n-wPad-1)), modeedge)imgGeoMean img.copy()for i in range(hPad, img_h hPad):for j in range(wPad, img_w wPad):prod np.prod(imgPad[i-hPad:ihPad1, j-wPad:jwPad1]*1.0)imgGeoMean[i-hPad][j-wPad] np.power(prod, order)plt.figure(figsize(9, 6))plt.subplot(131), plt.axis(off), plt.title(Original)plt.imshow(img, cmapgray, vmin0, vmax255)plt.subplot(132), plt.axis(off), plt.title(Arithmentic mean filter)plt.imshow(imgAriMean, cmapgray, vmin0, vmax255)plt.subplot(133), plt.axis(off), plt.title(Geometric mean filter)plt.imshow(imgGeoMean, cmapgray, vmin0, vmax255)plt.tight_layout()plt.show() 3.3 谐波平均滤波器Harmonic mean filter
谐波平均滤波器的复原图像 f^\hat{f}f^ 由下式给出 f^(x,y)mn∑(r,c)∈Sxy1/g(r,c)\hat{f}(x,y) \frac{mn}{\sum _{(r,c) \in Sxy} 1/{g(r,c)}} f^(x,y)∑(r,c)∈Sxy1/g(r,c)mn 谐波平均滤波器既能处理盐粒噪声白色噪点又能处理类似于高斯噪声的其他噪声但不能处理胡椒噪声黑色噪点。
例程 9.10谐波平均滤波器 # 9.10: 谐波平均滤波器 (Harmonic mean filter)img cv2.imread(../images/Fig0507b.tif, 0) # flags0 读取为灰度图像img_h img.shape[0]img_w img.shape[1]# 算术平均滤波 (Arithmentic mean filter)kSize (3, 3)kernalMean np.ones(kSize, np.float32) / (kSize[0]*kSize[1]) # 生成归一化盒式核imgAriMean cv2.filter2D(img, -1, kernalMean)# 谐波平均滤波器 (Harmonic mean filter)m, n 3, 3order m * nkernalMean np.ones((m,n), np.float32) # 生成盒式核hPad int((m-1) / 2)wPad int((n-1) / 2)imgPad np.pad(img.copy(), ((hPad, m-hPad-1), (wPad, n-wPad-1)), modeedge)epsilon 1e-8imgHarMean img.copy()for i in range(hPad, img_h hPad):for j in range(wPad, img_w wPad):sumTemp np.sum(1.0 / (imgPad[i-hPad:ihPad1, j-wPad:jwPad1] epsilon))imgHarMean[i-hPad][j-wPad] order / sumTempplt.figure(figsize(9, 6))plt.subplot(131), plt.axis(off), plt.title(Original)plt.imshow(img, cmapgray, vmin0, vmax255)plt.subplot(132), plt.axis(off), plt.title(Arithmentic mean filter)plt.imshow(imgAriMean, cmapgray, vmin0, vmax255)plt.subplot(133), plt.axis(off), plt.title(Harmonic mean filter)plt.imshow(imgHarMean, cmapgray, vmin0, vmax255)plt.tight_layout()plt.show() 3.4 反谐波平均滤波器Invert-harmonic mean filter
反谐波平均滤波器适用于降低或消除椒盐噪声复原图像 f^\hat{f}f^ 由下式给出 f^(x,y)∑(r,c)∈Sxyg(r,c)Q1∑(r,c)∈Sxyg(r,c)Q\hat{f}(x,y) \frac {\sum _{(r,c) \in Sxy} {g(r,c)}^{Q1}} {\sum _{(r,c) \in Sxy} {g(r,c)}^Q} f^(x,y)∑(r,c)∈Sxyg(r,c)Q∑(r,c)∈Sxyg(r,c)Q1 Q 称为滤波器的阶数Q 取正整数时可以消除胡椒噪声Q 取负整数时可以消除盐粒噪声但不能同时消除这两种噪声。使用反谐波平均滤波器必须知道噪声是亮噪声还是暗噪声据此选择 Q 的正负号。
反谐波平均滤波器当 Q0Q0Q0 时简化为算术平均滤波器当 Q−1Q-1Q−1 时简化为谐波平均滤波器。 例程 9.11反谐波平均滤波器 # 9.11: 反谐波平均滤波器 (Inv-harmonic mean filter)img cv2.imread(../images/Fig0508a.tif, 0) # flags0 读取为灰度图像img_h img.shape[0]img_w img.shape[1]m, n 3, 3order m * nkernalMean np.ones((m,n), np.float32) # 生成盒式核hPad int((m-1) / 2)wPad int((n-1) / 2)imgPad np.pad(img.copy(), ((hPad, m-hPad-1), (wPad, n-wPad-1)), modeedge)Q 1.5 # 反谐波平均滤波器 阶数epsilon 1e-8imgHarMean img.copy()imgInvHarMean img.copy()for i in range(hPad, img_h hPad):for j in range(wPad, img_w wPad):# 谐波平均滤波器 (Harmonic mean filter)sumTemp np.sum(1.0 / (imgPad[i-hPad:ihPad1, j-wPad:jwPad1] epsilon))imgHarMean[i-hPad][j-wPad] order / sumTemp# 反谐波平均滤波器 (Inv-harmonic mean filter)temp imgPad[i-hPad:ihPad1, j-wPad:jwPad1] epsilonimgInvHarMean[i-hPad][j-wPad] np.sum(np.power(temp, (Q1))) / np.sum(np.power(temp, Q) epsilon)plt.figure(figsize(9, 6))plt.subplot(131), plt.axis(off), plt.title(Original)plt.imshow(img, cmapgray, vmin0, vmax255)plt.subplot(132), plt.axis(off), plt.title(Harmonic mean filter)plt.imshow(imgHarMean, cmapgray, vmin0, vmax255)plt.subplot(133), plt.axis(off), plt.title(Invert harmonic mean)plt.imshow(imgInvHarMean, cmapgray, vmin0, vmax255)plt.tight_layout()plt.show() 3.5 统计排序滤波器
统计排序滤波器是空间滤波器其响应是基于滤波器邻域中的像素值的顺序排序结果决定了滤波器的输出。
统计排序包括中值滤波器、最大值滤波器、最小值滤波器、中点滤波器和修正阿尔法均值滤波器。
中值滤波器用预定义的像素邻域中的灰度中值来代替像素的值与线性平滑滤波器相比能有效地降低某些随机噪声且模糊度要小得多。
f^(x,y)median(r,c)∈Sxy{g(r,c)}\hat{f}(x,y) {median} _{(r,c) \in Sxy} \{g(r,c)\} f^(x,y)median(r,c)∈Sxy{g(r,c)}
由于需要排序操作中值滤波消耗的运算时间很长。 OpenCV 提供了 cv.medianBlur 函数实现中值滤波算法详见《例程 1.73图像的非线性滤波—中值滤波器》。
最大值滤波器用预定义的像素邻域中的灰度最大值来代替像素的值可用于找到图像中的最亮点或用于消弱与明亮区域相邻的暗色区域也可以用来降低胡椒噪声。
f^(x,y)max(r,c)∈Sxy{g(r,c)}\hat{f}(x,y) \max _{(r,c) \in Sxy} \{g(r,c)\} f^(x,y)(r,c)∈Sxymax{g(r,c)}
最小值滤波器用预定义的像素邻域中的灰度最小值来代替像素的值可用于找到图像中的最暗点或用于削弱与暗色区域相邻的明亮区域也可以用来降低盐粒噪声。
f^(x,y)min(r,c)∈Sxy{g(r,c)}\hat{f}(x,y) \min _{(r,c) \in Sxy} \{g(r,c)\} f^(x,y)(r,c)∈Sxymin{g(r,c)}
中点滤波器用预定义的像素邻域中的灰度的最大值与最小值的均值来代替像素的值注意中点的取值与中值常常是不同的。中点滤波器是统计排序滤波器与平均滤波器的结合适合处理随机分布的噪声例如高斯噪声、均匀噪声。
f^(x,y)[max(r,c)∈Sxy{g(r,c)}min(r,c)∈Sxy{g(r,c)}]/2\hat{f}(x,y) [\max _{(r,c) \in Sxy} \{g(r,c)\} \min _{(r,c) \in Sxy} \{g(r,c)\}]/2 f^(x,y)[(r,c)∈Sxymax{g(r,c)}(r,c)∈Sxymin{g(r,c)}]/2 例程 9.12统计排序滤波器 # 9.12: 统计排序滤波器 (Statistical sorting filter)img cv2.imread(../images/Fig0508a.tif, 0) # flags0 读取为灰度图像img_h img.shape[0]img_w img.shape[1]m, n 3, 3kernalMean np.ones((m, n), np.float32) # 生成盒式核# 边缘填充hPad int((m-1) / 2)wPad int((n-1) / 2)imgPad np.pad(img.copy(), ((hPad, m-hPad-1), (wPad, n-wPad-1)), modeedge)imgMedianFilter np.zeros(img.shape) # 中值滤波器imgMaxFilter np.zeros(img.shape) # 最大值滤波器imgMinFilter np.zeros(img.shape) # 最小值滤波器imgMiddleFilter np.zeros(img.shape) # 中点滤波器for i in range(img_h):for j in range(img_w):# # 1. 中值滤波器 (median filter)pad imgPad[i:im, j:jn]imgMedianFilter[i, j] np.median(pad)# # 2. 最大值滤波器 (maximum filter)pad imgPad[i:im, j:jn]imgMaxFilter[i, j] np.max(pad)# # 3. 最小值滤波器 (minimum filter)pad imgPad[i:im, j:jn]imgMinFilter[i, j] np.min(pad)# # 4. 中点滤波器 (middle filter)pad imgPad[i:im, j:jn]imgMiddleFilter[i, j] int(pad.max()/2 pad.min()/2)plt.figure(figsize(9, 7))plt.subplot(221), plt.axis(off), plt.title(median filter)plt.imshow(imgMedianFilter, cmapgray, vmin0, vmax255)plt.subplot(222), plt.axis(off), plt.title(maximum filter)plt.imshow(imgMaxFilter, cmapgray, vmin0, vmax255)plt.subplot(223), plt.axis(off), plt.title(minimum filter)plt.imshow(imgMinFilter, cmapgray, vmin0, vmax255)plt.subplot(224), plt.axis(off), plt.title(middle filter)plt.imshow(imgMiddleFilter, cmapgray, vmin0, vmax255)plt.tight_layout()plt.show() 程序说明
需要说明的是本例程和各种滤波器图像处理结果是为了说明滤波器实现的编程方法和程序运行结果。图中一些滤波器图像处理的效果较差并不能全面反映该滤波器的性能只能说明该滤波器不适合处理某些类型的噪声。关于统计滤波器在处理不同噪声的选择和比较可以参考冈萨雷斯《数字图像处理第四版》第五章的相关内容。 3.6 修正阿尔法均值滤波器Modified alpha-mean filter
修正阿尔法均值滤波器也属于统计排序滤波器其思想类似于比赛中去掉最高分和最低分后计算平均分的方法。
令 SxySxySxy 表示中心在点 (x,y)(x,y)(x,y) 、大小为 m∗nm*nm∗n 的矩形子窗口邻域的一组坐标修正阿尔法均值滤波器在由 SxySxySxy 定义的邻域中删除 d 个最低灰度值和 d 个最高灰度值计算剩余像素 gR(r,c)g_R(r,c)gR(r,c) 的算术平均值作为输出结果即
f^(x,y)1mn−2d∑(r,c)∈SRgR(r,c)\hat{f}(x,y) \frac{1}{mn-2d} \sum _{(r,c) \in S_R} g_R(r,c) f^(x,y)mn−2d1(r,c)∈SR∑gR(r,c) d 的取值范围是 [0,mn/2−1][0, mn/2-1][0,mn/2−1]。选择 d 的大小对图像处理的效果影响很大当 d0d0d0 时简化为算术平均滤波器当 dmn/2−1dmn/2-1dmn/2−1 简化为中值滤波器。d 取其它值时适合于处理多种混合噪声如高斯噪声和椒盐噪声。 例程 9.13修正阿尔法均值滤波器 # 9.13: 修正阿尔法均值滤波器 (Modified alpha-mean filter)img cv2.imread(../images/Fig0507b.tif, 0) # flags0 读取为灰度图像img_h img.shape[0]img_w img.shape[1]m, n 5, 5kernalMean np.ones((m, n), np.float32) # 生成盒式核# 边缘填充hPad int((m-1) / 2)wPad int((n-1) / 2)imgPad np.pad(img.copy(), ((hPad, m-hPad-1), (wPad, n-wPad-1)), modeedge)imgAlphaFilter0 np.zeros(img.shape)imgAlphaFilter1 np.zeros(img.shape)imgAlphaFilter2 np.zeros(img.shape)for i in range(img_h):for j in range(img_w):# 邻域 m * npad imgPad[i:im, j:jn]padSort np.sort(pad.flatten()) # 对邻域像素按灰度值排序d 1sumAlpha np.sum(padSort[d:m*n-d-1]) # 删除 d 个最大灰度值, d 个最小灰度值imgAlphaFilter0[i, j] sumAlpha / (m*n-2*d) # 对剩余像素进行算术平均d 2sumAlpha np.sum(padSort[d:m*n-d-1])imgAlphaFilter1[i, j] sumAlpha / (m*n-2*d)d 4sumAlpha np.sum(padSort[d:m*n-d-1])imgAlphaFilter2[i, j] sumAlpha / (m*n-2*d)plt.figure(figsize(9, 7))plt.subplot(221), plt.axis(off), plt.title(Original)plt.imshow(img, cmapgray, vmin0, vmax255)plt.subplot(222), plt.axis(off), plt.title(Modified alpha-mean(d1))plt.imshow(imgAlphaFilter0, cmapgray, vmin0, vmax255)plt.subplot(223), plt.axis(off), plt.title(Modified alpha-mean(d2))plt.imshow(imgAlphaFilter1, cmapgray, vmin0, vmax255)plt.subplot(224), plt.axis(off), plt.title(Modified alpha-mean(d4))plt.imshow(imgAlphaFilter2, cmapgray, vmin0, vmax255)plt.tight_layout()plt.show() 3.7 自适应局部降噪滤波器
前述滤波器直接应用到图像处理并未考虑图像本身的特征。自适应滤波器的特性根据 m∗nm*nm∗n 矩形邻域 SxySxySxy 定义的滤波区域内的图像的统计特性变化。通常自适应滤波器的性能优于前述的滤波器但滤波器的复杂度和计算量也更大。
均值和方差是随机变量最基础的统计量。在图像处理中均值是像素邻域的平均灰度方差是像素邻域的图像对比度。
令 SxySxySxy 表示中心在点 (x,y)(x,y)(x,y) 、大小为 m∗nm*nm∗n 的矩形子窗口邻域滤波器在由 SxySxySxy 定义的邻域操作。
噪声图像在点 (x,y)(x,y)(x,y) 的值 g(x,y)g(x,y)g(x,y)噪声的方差 ση2\sigma^2_{\eta}ση2 由噪声图像估计SxySxySxy 中像素的局部平均灰度为 zˉSxy\bar{z}_{Sxy}zˉSxy灰度的局部方差为 σSxy2\sigma^2_{Sxy}σSxy2。
使用自适应局部降噪滤波器复原的图像 f^\hat{f}f^ 在点 (x,y)(x,y)(x,y) 的值由如下自适应表达式描述 f^(x,y)g(x,y)−ση2σSxy2[g(x,y)−zSxyˉ]\hat{f}(x,y) g(x,y) - \frac{\sigma^2_{\eta}}{\sigma^2_{Sxy}} [g(x,y)-\bar{z_{Sxy}}] f^(x,y)g(x,y)−σSxy2ση2[g(x,y)−zSxyˉ]
例程 9.14自适应局部降噪滤波器 # 9.14: 自适应局部降噪滤波器 (Adaptive local noise reduction filter)img cv2.imread(../images/Fig0507b.tif, 0) # flags0 读取为灰度图像hImg img.shape[0]wImg img.shape[1]m, n 5, 5imgAriMean cv2.boxFilter(img, -1, (m, n)) # 算术平均滤波# 边缘填充hPad int((m-1) / 2)wPad int((n-1) / 2)imgPad np.pad(img.copy(), ((hPad, m-hPad-1), (wPad, n-wPad-1)), modeedge)# 估计原始图像的噪声方差 sigmaEtamean, stddev cv2.meanStdDev(img)sigmaEta stddev ** 2print(sigmaEta)# 自适应局部降噪epsilon 1e-8imgAdaLocal np.zeros(img.shape)for i in range(hImg):for j in range(wImg):pad imgPad[i:im, j:jn] # 邻域 Sxy, m*ngxy img[i,j] # 含噪声图像的像素点zSxy np.mean(pad) # 局部平均灰度sigmaSxy np.var(pad) # 灰度的局部方差rateSigma min(sigmaEta / (sigmaSxy epsilon), 1.0) # 加性噪声假设sigmaEta/sigmaSxy 1imgAdaLocal[i, j] gxy - rateSigma * (gxy - zSxy)plt.figure(figsize(9, 6))plt.subplot(131), plt.axis(off), plt.title(Original)plt.imshow(img, cmapgray, vmin0, vmax255)plt.subplot(132), plt.axis(off), plt.title(Arithmentic mean filter)plt.imshow(imgAriMean, cmapgray, vmin0, vmax255)plt.subplot(133), plt.axis(off), plt.title(Adaptive local filter)plt.imshow(imgAdaLocal, cmapgray, vmin0, vmax255)plt.tight_layout()plt.show() 3.8 自适应中值滤波器Adaptive median filter
中值滤波器的窗口尺寸是固定大小不变的不能同时兼顾去噪和保护图像的细节在噪声的密度较小时的性能较好当噪声概率较高时的性能就会劣化。
自适应中值滤波器根据预先设定的条件在滤波的过程中动态改变滤波器的窗口尺寸大小进一步地根据条件判断当前像素是否噪声由此决定是否用邻域中值替换当前像素。
自适应中值滤波器可以处理较大概率的脉冲噪声平滑非脉冲噪声尽可能保护图像细节信息避免图像边缘的细化或者粗化。
令 SxySxySxy 表示中心在点 (x,y)(x,y)(x,y) 、大小为 m∗nm*nm∗n 的矩形子窗口邻域滤波器在由 SxySxySxy 定义的邻域操作。
Zmin、Zmax、ZmedZmin、Zmax、ZmedZmin、Zmax、Zmed 表示 SxySxySxy 中的最小灰度值、最大灰度值和灰度值的中值ZxyZxyZxy 是点 (x,y)(x,y)(x,y) 的灰度值SmaxSmaxSmax 是允许的最大窗口尺寸。
自适应中值滤波器分为两个过程 Step A A1 Zmed - ZminA2 Zmed - Zmax如果 A10 且 A20则跳转到 Step B否则增大窗口尺寸如果增大后的尺寸 ≤ Smax则重复 A否则输出 Zmed Step B B1 Zxy - ZminB2 Zxy - Zmax如果 B10 且 B20则输出 Zxy否则输出 Zmed
StepA 的目的是确定当前窗口内得到的中值 Zmed 是否是噪声。如果 ZminZmedZmax则中值 Zmed 不是噪声转到StepB。如果ZminZxyZmax则 Zxy 不是噪声滤波器输出 Zxy。如果不满足上述条件则判定 Zxy 是噪声输出中值 Zmed。
如果在StepA中Zmed 不满足条件 ZminZmedZmax则可判断得到的中值 Zmed 是噪声。这时要增大滤波器的窗口尺寸在更大的范围内寻找一个非噪声点的中值。
因此如果图像中噪声的概率较低自适应中值滤波器可以使用较小的窗口尺寸以提高计算速度反之如果噪声的概率较高则需要增大滤波器的窗口尺寸以改善滤波效果。 例程 9.15自适应中值滤波器 # # 9.15: 自适应中值滤波器 (Adaptive median filter)img cv2.imread(../images/Fig0514a.tif, 0) # flags0 读取为灰度图像hImg img.shape[0]wImg img.shape[1]smax 7 # 允许最大窗口尺寸m, n smax, smaximgAriMean cv2.boxFilter(img, -1, (m, n)) # 算术平均滤波# 边缘填充hPad int((m-1) / 2)wPad int((n-1) / 2)imgPad np.pad(img.copy(), ((hPad, m-hPad-1), (wPad, n-wPad-1)), modeedge)imgMedianFilter np.zeros(img.shape) # 中值滤波器imgAdaMedFilter np.zeros(img.shape) # 自适应中值滤波器for i in range(hPad, hPadhImg):for j in range(wPad, wPadwImg):# 1. 中值滤波器 (Median filter)ksize 3k int(ksize/2)pad imgPad[i-k:ik1, j-k:jk1] # 邻域 Sxy, m*nimgMedianFilter[i-hPad, j-wPad] np.median(pad)# 2. 自适应中值滤波器 (Adaptive median filter)ksize 3k int(ksize/2)pad imgPad[i-k:ik1, j-k:jk1]zxy img[i-hPad][j-wPad]zmin np.min(pad)zmed np.median(pad)zmax np.max(pad)if zmin zmed zmax:if zmin zxy zmax:imgAdaMedFilter[i-hPad, j-wPad] zxyelse:imgAdaMedFilter[i-hPad, j-wPad] zmedelse:while True:ksize ksize 2if zmin zmed zmax or ksize smax:breakk int(ksize / 2)pad imgPad[i-k:ik1, j-k:jk1]zmed np.median(pad)zmin np.min(pad)zmax np.max(pad)if zmin zmed zmax or ksize smax:if zmin zxy zmax:imgAdaMedFilter[i-hPad, j-wPad] zxyelse:imgAdaMedFilter[i-hPad, j-wPad] zmedplt.figure(figsize(9, 6))plt.subplot(131), plt.axis(off), plt.title(Original)plt.imshow(img, cmapgray, vmin0, vmax255)plt.subplot(132), plt.axis(off), plt.title(Median filter)plt.imshow(imgMedianFilter, cmapgray, vmin0, vmax255)plt.subplot(133), plt.axis(off), plt.title(Adaptive median filter)plt.imshow(imgAdaMedFilter, cmapgray, vmin0, vmax255)plt.tight_layout()plt.show()4. 频率域滤波图像复原
通过频率域滤波可以有效分析并滤除周期噪声其理论基础是傅里叶变换后周期噪声在对应周期干扰的频率显示为集中突发的能量因此可以使用选择性滤波器来分离滤除噪声。
4.1 陷波滤波器**Notch Filter**
陷波滤波器阻止或通过预定的频率矩形邻域中的频率可以很好地复原被周期性噪声干扰的图像。在《5.2 陷波滤波器Notch Filter》已进行介绍并给出了例程。
陷波滤波器可以在某一个频率点迅速衰减输入信号以达到阻碍此频率信号通过的滤波效果的滤波器。
陷波带阻滤波器的传递函数是中心平移到陷波中心的各个高通滤波器的乘积 HNR(u,v)∏k1QHk(u,v)H−k(u,v)H_{NR}(u,v) \prod_{k1}^Q H_k(u,v) H_{-k}(u,v) HNR(u,v)k1∏QHk(u,v)H−k(u,v)
其中滤波器的距离计算公式为 Dk(u,v)(u−M/2−uk)2(v−N/2−vk)2D−k(u,v)(u−M/2uk)2(v−N/2vk)2D_k(u,v) \sqrt{(u-M/2-u_k)^2 (v-N/2-v_k)^2} \\ D_{-k}(u,v) \sqrt{(u-M/2u_k)^2 (v-N/2v_k)^2} Dk(u,v)(u−M/2−uk)2(v−N/2−vk)2D−k(u,v)(u−M/2uk)2(v−N/2vk)2 例如具有 3个陷波对的 n 阶巴特沃斯陷波带阻滤波器为 HNR(u,v)∏k13[11[D0k/Dk(u,v)]n][11[D−k/Dk(u,v)]n]H_{NR}(u,v) \prod_{k1}^3 [\frac{1}{1[D_{0k}/D_k(u,v)]^n}] [\frac{1}{1[D_{-k}/D_k(u,v)]^n}] HNR(u,v)k1∏3[1[D0k/Dk(u,v)]n1][1[D−k/Dk(u,v)]n1]
例程 9.16陷波带阻滤波器的传递函数 # 9.16: 陷波带阻滤波器的传递函数def ideaBondResistFilter(shape, radius10, w5): # 理想带阻滤波器u, v np.meshgrid(np.arange(shape[1]), np.arange(shape[0]))D np.sqrt((u - shape[1]//2)**2 (v - shape[0]//2)**2)D0 radiushalfW w/2kernel np.piecewise(D, [DD0halfW, DD0-halfW], [1, 0])kernel 1 - kernel # 带阻return kerneldef butterworthNRFilter(img, radius10, uk60, vk80, n1): # 巴特沃斯陷波带阻滤波器M, N img.shape[1], img.shape[0]u, v np.meshgrid(np.arange(M), np.arange(N))Dkm np.sqrt((u - M//2 - uk)**2 (v - N//2 - vk)**2) # D_kDkp np.sqrt((u - M//2 uk)**2 (v - N//2 vk)**2) # D_-kD0 radiusn2 n * 2epsilon 1e-6kernel (1 / (1 (D0 / (Dkmepsilon))**n2)) * (1 / (1 (D0 / (Dkpepsilon))**n2))return kerneldef gaussNRFilter(img, radius10, uk60, vk80): # 高斯陷波带阻滤波器M, N img.shape[1], img.shape[0]u, v np.meshgrid(np.arange(M), np.arange(N))Dkm np.sqrt((u - M//2 - uk)**2 (v - N//2 - vk)**2) # D_kDkp np.sqrt((u - M//2 uk)**2 (v - N//2 vk)**2) # D_-kD0 radiuskernel (1 - np.exp(-(Dkm**2)/(D0**2))) * (1 - np.exp(-(Dkp**2)/(D0**2)))return kerneldef ideaNRFilter(img, radius10, uk60, vk80): # 高斯陷波带阻滤波器M, N img.shape[1], img.shape[0]u, v np.meshgrid(np.arange(M), np.arange(N))Dkm np.sqrt((u - M//2 - uk)**2 (v - N//2 - vk)**2) # D_kDkp np.sqrt((u - M//2 uk)**2 (v - N//2 vk)**2) # D_-kD0 radiusk1 Dkm.copy()k1[DkmD0] 1k1[DkmD0] 0k2 Dkp.copy()k2[DkpD0] 1k2[DkpD0] 0kernel k1 * k2return kernel# 理想、高斯、巴特沃斯陷波带阻滤波器传递函数img np.zeros([128, 128])shape img.shaperadius 15INRF ideaNRFilter(img, radiusradius, uk20, vk30) # (uk,vk) 陷波中心GNRF gaussNRFilter(img, radiusradius, uk20, vk30)BNRF butterworthNRFilter(img, radiusradius, uk20, vk30, n2)filters [INRF, GNRF, BNRF]u, v np.mgrid[-1:1:2.0/shape[0], -1:1:2.0/shape[1]]fig plt.figure(figsize(10, 8))for i in range(3):nrFilter eval(filters[i]).copy()ax1 fig.add_subplot(3, 3, 3*i1)ax1.imshow(nrFilter, gray)ax1.set_title(filters[i]), ax1.set_xticks([]), ax1.set_yticks([])ax2 plt.subplot(3,3,3*i2, projection3d)ax2.set_title(transfer function)ax2.plot_wireframe(u, v, nrFilter, rstride2, linewidth0.5, colorc)ax2.set_xticks([]), ax2.set_yticks([]), ax2.set_zticks([])ax3 plt.subplot(3,3,3*i3)profile nrFilter[:, 30]ax3.plot(profile), ax3.set_title(profile), ax3.set_xticks([]), ax3.set_yticks([])plt.show()例程 9.17陷波带阻滤波器消除周期噪声干扰 # 9.17: 陷波带阻滤波器消除周期噪声干扰def butterworthNRFilter(img, radius10, uk10, vk10, n2): # 巴特沃斯陷波带阻滤波器M, N img.shape[1], img.shape[0]u, v np.meshgrid(np.arange(M), np.arange(N))Dm np.sqrt((u - M//2 - uk)**2 (v - N//2 - vk)**2)Dp np.sqrt((u - M//2 uk)**2 (v - N//2 vk)**2)D0 radiusn2 2 * nkernel (1 / (1 (D0 / (Dm 1e-6))**n2)) * (1 / (1 (D0 / (Dp 1e-6))**n2))return kernel# (1) 读取原始图像img cv2.imread(../images/Fig0505a.tif, flags0) # flags0 读取为灰度图像imgFloat32 np.float32(img) # 将图像转换成 float32rows, cols img.shape[:2] # 图片的高度和宽度fig plt.figure(figsize(9, 6))plt.subplot(231), plt.title(Original image), plt.axis(off), plt.imshow(img, cmapgray)# (2) 中心化, centralized 2d array f(x,y) * (-1)^(xy)mask np.ones(img.shape)mask[1::2, ::2] -1mask[::2, 1::2] -1fImage imgFloat32 * mask # f(x,y) * (-1)^(xy)# (3) 快速傅里叶变换rPadded cv2.getOptimalDFTSize(rows) # 最优 DFT 扩充尺寸cPadded cv2.getOptimalDFTSize(cols) # 用于快速傅里叶变换dftImage np.zeros((rPadded, cPadded, 2), np.float32) # 对原始图像进行边缘扩充dftImage[:rows, :cols, 0] fImage # 边缘扩充下侧和右侧补0cv2.dft(dftImage, dftImage, cv2.DFT_COMPLEX_OUTPUT) # 快速傅里叶变换dftAmp cv2.magnitude(dftImage[:,:,0], dftImage[:,:,1]) # 傅里叶变换的幅度谱 (rPad, cPad)dftAmpLog np.log(1.0 dftAmp) # 幅度谱对数变换以便于显示dftAmpNorm np.uint8(cv2.normalize(dftAmpLog, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]plt.subplot(232), plt.axis(off), plt.title(DFT spectrum)plt.imshow(dftAmpNorm, cmapgray)plt.arrow(445, 370, 25, 30, width5, length_includes_headTrue, shapefull) # 在图像上加上箭头plt.arrow(550, 490, -25, -30, width5, length_includes_headTrue, shapefull) # 在图像上加上箭头# (4) 构建陷波带阻滤波器 传递函数BRFilter butterworthNRFilter(dftImage, radius15, uk25, vk16, n3) # 巴特沃斯陷波带阻滤波器, 处理周期噪声plt.subplot(233), plt.axis(off), plt.title(Butterworth notch resist filter)plt.imshow(BRFilter, cmapgray)# (5) 在频率域修改傅里叶变换: 傅里叶变换 点乘 陷波带阻滤波器dftFilter np.zeros(dftImage.shape, dftImage.dtype) # 快速傅里叶变换的尺寸(优化尺寸)for i in range(2):dftFilter[:rPadded, :cPadded, i] dftImage[:rPadded, :cPadded, i] * BRFilter# 频域滤波傅里叶变换的傅里叶谱nrfDftAmp cv2.magnitude(dftFilter[:, :, 0], dftFilter[:, :, 1]) # 傅里叶变换的幅度谱nrfDftAmpLog np.log(1.0 nrfDftAmp) # 幅度谱对数变换以便于显示nrfDftAmpNorm np.uint8(cv2.normalize(nrfDftAmpLog, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]plt.subplot(234), plt.axis(off), plt.title(BNRF DFT Spectrum)plt.imshow(nrfDftAmpNorm, cmapgray)# (6) 对频域滤波傅里叶变换 执行傅里叶逆变换并只取实部idft np.zeros(dftAmp.shape, np.float32) # 快速傅里叶变换的尺寸(优化尺寸)cv2.dft(dftFilter, idft, cv2.DFT_REAL_OUTPUT cv2.DFT_INVERSE cv2.DFT_SCALE)# (7) 中心化, centralized 2d array g(x,y) * (-1)^(xy)mask2 np.ones(dftAmp.shape)mask2[1::2, ::2] -1mask2[::2, 1::2] -1idftCen idft * mask2 # g(x,y) * (-1)^(xy)plt.subplot(235), plt.axis(off), plt.title(g(x,y)*(-1)^(xy))plt.imshow(idftCen, cmapgray)# (8) 截取左上角大小和输入图像相等idftCenClip np.clip(idftCen, 0, 255) # 截断函数将数值限制在 [0,255]imgFiltered idftCenClip.astype(np.uint8)imgFiltered imgFiltered[:rows, :cols]plt.subplot(236), plt.axis(off), plt.title(BNRF filtered image)plt.imshow(imgFiltered, cmapgray)plt.tight_layout()plt.show()print(image.shape:{}.format(img.shape))print(imgFloat32.shape:{}.format(imgFloat32.shape))print(dftImage.shape:{}.format(dftImage.shape))print(dftAmp.shape:{}.format(dftAmp.shape))print(idft.shape:{}.format(idft.shape))print(dftFilter.shape:{}.format(dftFilter.shape))print(imgFiltered.shape:{}.format(imgFiltered.shape))程序说明
傅立叶变换的频谱反映能量分布。通过中心化移频到原点以后傅里叶变换的频谱图是以原点为中心对称分布的。中心化不仅可以清晰地看出图像频率分布还可以分离出有周期性规律的干扰信号。
本例程中的图像受到正弦噪声的干扰从中心化的频谱图可以看出除原点以外还存在一对对称分布的亮点箭头指示处这就干扰噪声产生的。因此在该亮点位置设计陷波带阻滤波器可以消除干扰正弦噪声。 4.2 最优陷波滤波器
当存在多个干扰分量时陷波滤波器的设计比较困难而且容易过多地滤除图像信息。特别是当干扰分量不是单频率突发能量而是携带干扰模式信息的宽裙摆在正常的变换背景中就难以检测和分离。
最优陷波滤波方法在最小化复原估计的局部方差方面是最优的该方法首先分离干扰模式的各个主要贡献然后从被污染图像中减去该模式的一个可变加权部分。 5. 估计退化函数
估计图像复原中所用的退化函数主要有三种方法观察法、试验法和数学建模方法。
5.1 观察法估计退化函数
为了降低噪声的影响选择一个信号内容很强的区域如一个高对比度区域。
令gs(x,y)g_s(x,y)gs(x,y) 表示观察子图像f^s(x,y)\hat{f}_s(x,y)f^s(x,y) 表示处理后的子图像则 Hs(u,v)Gs(u,v)f^s(u,v)H_s(u,v) \frac{G_s(u,v)}{\hat{f}_s(u,v)} Hs(u,v)f^s(u,v)Gs(u,v) 由该函数的特性可以根据位置不变假设来推断完整的退化函数 H(u,v)H(u,v)H(u,v)。 5.2 试验法估计退化函数
如果设备与获取退化图像的设备相似原则上就可能获得退化模型的精确估计——如果能使用原设备当然更好。在各种系统设置下获取类似于退化图像的图像然后再该系统设置条件下对一个冲激小光点成像就得到退化的冲激响应。 H(u,v)G(u,v)AH(u,v) \frac{G(u,v)}{A} H(u,v)AG(u,v) G(u,v)G(u,v)G(u,v) 是观察图像的傅里叶变换A 是描述冲激强度的常量。 5.3 模型法估计退化函数
分析导致退化的原因根据基本原理提出退化模型如湍流导致的模糊、匀速运动导致的模糊可以基于模型更加准确地估计退化函数。
下面以运动模糊和大气湍流模型为例采用退化模型对图像的退化建模。 例程 9.18运动模糊退化模型
运动模糊是相机物体背景间相对运动造成的效果通常由于长时间曝光或场景内的物体快速移动导致在摄影中可以借助移动镜头追踪移动的物体来避免。
对匀速线性运动模糊建模假设图像 f(x,y)f(x,y)f(x,y) 做平面运动运动在 x、y 方向的时变分量分别为 x0(t)at/Tx_0(t)at/Tx0(t)at/T、y0(t)bt/Ty_0(t)bt/Ty0(t)bt/T。记录介质上任意点的总曝光量是瞬时曝光量的积分可以建立运动模糊退化函数模型 H(u,v)Tπ(uavb)sin[π(uavb)]e−jπ(uavb)H(u,v) \frac{T}{\pi (uavb)} sin[\pi(uavb)]e^{-j \pi (uavb)} H(u,v)π(uavb)Tsin[π(uavb)]e−jπ(uavb) # 9.18: 运动模糊退化图像 (Motion blur degradation)def motionBlur(image, degree10, angle45):image np.array(image)center (degree/2, degree/2) # 旋转中心M cv2.getRotationMatrix2D(center, angle, 1) # 无损旋转kernel np.diag(np.ones(degree) / degree) # 运动模糊内核kernel cv2.warpAffine(kernel, M, (degree, degree))blurred cv2.filter2D(image, -1, kernel) # 图像卷积blurredNorm np.uint8(cv2.normalize(blurred, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]return blurredNorm# 运动模糊图像img cv2.imread(../images/Fig0526a.tif, 0) # flags0 读取为灰度图像imgBlur1 motionBlur(img, degree30, angle45)imgBlur2 motionBlur(img, degree40, angle45)imgBlur3 motionBlur(img, degree60, angle45)plt.figure(figsize(9, 6))plt.subplot(131), plt.title(degree20), plt.axis(off), plt.imshow(imgBlur1, gray)plt.subplot(132), plt.title(degree40), plt.axis(off), plt.imshow(imgBlur2, gray)plt.subplot(133), plt.title(degree60), plt.axis(off), plt.imshow(imgBlur3, gray)plt.tight_layout()plt.show()例程 9.19湍流模糊退化模型
湍流是自然界中普遍存在的一种复杂的流动现象。物体通过湍流大气成像时受到湍流效应的影响出现光强闪烁、光束方向漂移、光束宽度扩展及接收面上相位的起伏造成图像模糊和抖动甚至扭曲变形。
Hufnagel and Stanley 根据大气湍流的物理特性提出一种退化模型
H(u,v)e−k(u2v2)5/6H(u,v) e^{-k (u^2v^2)^{5/6}} H(u,v)e−k(u2v2)5/6
k 是湍流常数反映湍流的强烈程度。 # 9.19: 湍流模糊退化模型 (turbulence blur degradation model)def getDegradedImg(image, Huv): # 根据退化模型生成退化图像rows, cols image.shape[:2] # 图片的高度和宽度# (1) 中心化, centralized 2d array f(x,y) * (-1)^(xy)mask np.ones((rows, cols))mask[1::2, ::2] -1mask[::2, 1::2] -1imageCen image * mask# (2) 快速傅里叶变换dftImage np.zeros((rows, cols, 2), np.float32)dftImage[:, :, 0] imageCencv2.dft(dftImage, dftImage, cv2.DFT_COMPLEX_OUTPUT) # 快速傅里叶变换 (rows, cols, 2)# (4) 构建 频域滤波器传递函数:Filter np.zeros((rows, cols, 2), np.float32) # (rows, cols, 2)Filter[:, :, 0], Filter[:, :, 1] Huv, Huv# (5) 在频率域修改傅里叶变换: 傅里叶变换 点乘 滤波器传递函数dftFilter dftImage * Filter# (6) 对修正傅里叶变换 进行傅里叶逆变换并只取实部idft np.ones((rows, cols), np.float32) # 快速傅里叶变换的尺寸cv2.dft(dftFilter, idft, cv2.DFT_REAL_OUTPUT cv2.DFT_INVERSE cv2.DFT_SCALE) # 只取实部# (7) 中心化, centralized 2d array g(x,y) * (-1)^(xy)mask2 np.ones(dftImage.shape[:2])mask2[1::2, ::2] -1mask2[::2, 1::2] -1idftCen idft * mask2 # g(x,y) * (-1)^(xy)# (8) 截取左上角大小和输入图像相等imgDegraded np.uint8(cv2.normalize(idftCen, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]# print(image.shape, dftFilter.shape, imgDegraded.shape)return imgDegradeddef turbulenceBlur(img, k0.001): # 湍流模糊传递函数# H(u,v) exp(-k(u^2v^2)^5/6)M, N img.shape[1], img.shape[0]u, v np.meshgrid(np.arange(M), np.arange(N))radius (u - M//2)**2 (v - N//2)**2kernel np.exp(-k * np.power(radius, 5/6))return kernel# 读取原始图像img cv2.imread(../images/Fig0525a.tif, 0) # flags0 读取为灰度图像# 生成湍流模糊图像HBlur1 turbulenceBlur(img, k0.001) # 湍流模糊传递函数imgBlur1 getDegradedImg(img, HBlur1) # 生成湍流模糊图像HBlur2 turbulenceBlur(img, k0.0025)imgBlur2 getDegradedImg(img, HBlur2)plt.figure(figsize(9, 6))plt.subplot(131), plt.title(origin), plt.axis(off), plt.imshow(img, gray)plt.subplot(132), plt.title(turbulence blur(k0.001)), plt.axis(off), plt.imshow(imgBlur1, gray)plt.subplot(133), plt.title(turbulence blur(k0.0025)), plt.axis(off), plt.imshow(imgBlur2, gray)plt.tight_layout()plt.show() 6. 退化图像复原
图像复原是对图像退化的过程进行估计并补偿退化过程造成的失真以便获得未经退化的原始图像或原始图像的最优估值从而改善图像质量的一种方法。
典型的图像复原方法是根据图像退化的先验知识建立退化模型以退化模型为基础采用滤波等手段进行处理使复原后的图像符合一定的准则达到改善图像质量的目的。
因此图像复原是沿着质量降低的逆过程来重现真实的原始图像通过去模糊函数而去除图像模糊。
6.1 退化图像的逆滤波Inverse filter
图像退化表示为退化算子 H\mathcal{H}H 退化函数可以用观察法、试验法或建模法估计则通过逆滤波就可以直接实现图像复原。用退化图像的傅里叶变换除以退化函数的傅里叶变换得到原始图像的傅里叶变换估计 F^(u,v)G(u,v)H(u,v)\hat{F}(u,v) \frac{G(u,v)}{H(u,v)} F^(u,v)H(u,v)G(u,v)
但是由于实际上退化图像是退化算子与加性噪声项共同作用的结果由此得到 F^(u,v)F(u,v)N(u,v)H(u,v)\hat{F}(u,v) F(u,v) \frac{N(u,v)}{H(u,v)} F^(u,v)F(u,v)H(u,v)N(u,v)
这表明即使获得退化函数 H(u,v)H(u,v)H(u,v) 的估计由于噪声项是未知的因此也不能准确地复原原始图像。
进一步地如果退化函数为 0 或很小则噪声项的影响将非常严重信噪比低。这时需要将频率限制到原点附近进行分析可以减少遇到零值的可能性。 例程 9.20湍流模糊退化图像的逆滤波
如前所述通过湍流退化模型可以得到退化图像。使用该退化模型进行逆滤波退化函数与生成退化图像所用的退化函数相反 H(u,v)e−k[(u−M/2)2(v−N/2)2]5/6H(u,v) e^{-k [(u-M/2)^2(v-N/2)^2]^{5/6}} H(u,v)e−k[(u−M/2)2(v−N/2)2]5/6
但是直接使用退化模型 H(u,v) 逆滤波的结果D0full很差用理想低通滤波器对退化模型 H(u,v) 在半径 D0 之外截止后则视觉效果较好。 # 9.20: 湍流模糊退化图像的逆滤波def turbulenceBlur(img, k0.001): # 湍流模糊传递函数: H(u,v) exp(-k(u^2v^2)^5/6)M, N img.shape[1], img.shape[0]u, v np.meshgrid(np.arange(M), np.arange(N))radius (u - M//2)**2 (v - N//2)**2kernel np.exp(-k * np.power(radius, 5/6))return kerneldef getDegradedImg(image, Huv, eps): # 根据退化模型生成退化图像# (1) 傅里叶变换, 中心化fft np.fft.fft2(image.astype(np.float32)) # 傅里叶变换fftShift np.fft.fftshift(fft) # 将低频分量移动到频域图像中心# (2) 在频率域修改傅里叶变换: 傅里叶变换 点乘 滤波器传递函数fftShiftFilter fftShift * Huv # Guv Fuv * Huv# (3) 对修正傅里叶变换 进行傅里叶逆变换逆中心化invShift np.fft.ifftshift(fftShiftFilter) # 将低频分量逆转换回图像四角imgIfft np.fft.ifft2(invShift) # 逆傅里叶变换返回值是复数数组imgDegraded np.uint8(cv2.normalize(np.abs(imgIfft), None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]return imgDegradeddef ideaLPFilter(img, radius10): # 理想低通滤波器M, N img.shape[1], img.shape[0]u, v np.meshgrid(np.arange(M), np.arange(N))D np.sqrt((u - M//2)**2 (v - N//2)**2)kernel np.zeros(img.shape[:2], np.float32)kernel[D radius] 1return kerneldef inverseFilter(image, Huv, D0): # 根据退化模型逆滤波# (1) 傅里叶变换, 中心化fft np.fft.fft2(image.astype(np.float32)) # 傅里叶变换fftShift np.fft.fftshift(fft) # 将低频分量移动到频域图像中心# (2) 在频率域修改傅里叶变换: 傅里叶变换 点乘 滤波器传递函数if D00:fftShiftFilter fftShift / Huv # Guv Fuv / Huvelse:lpFilter ideaLPFilter(image, radiusD0)fftShiftFilter fftShift / Huv * lpFilter # Guv Fuv / Huv# (3) 对修正傅里叶变换 进行傅里叶逆变换逆中心化invShift np.fft.ifftshift(fftShiftFilter) # 将低频分量逆转换回图像四角imgIfft np.fft.ifft2(invShift) # 逆傅里叶变换返回值是复数数组imgRebuild np.uint8(cv2.normalize(np.abs(imgIfft), None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]return imgRebuild# 读取原始图像img cv2.imread(../images/Fig0525a.tif, 0) # flags0 读取为灰度图像# 生成湍流模糊图像HTurb turbulenceBlur(img, k0.0025)imgBlur np.abs(getDegradedImg(img, HTurb, 0.0))print(imgBlur.max(), imgBlur.min())# # 逆滤波imgRebuild inverseFilter(imgBlur, HTurb, 480) # Huv 全滤波器imgRebuild1 inverseFilter(imgBlur, HTurb, D040) # 在半径 D0 之外 Huv 截止imgRebuild2 inverseFilter(imgBlur, HTurb, D070)imgRebuild3 inverseFilter(imgBlur, HTurb, D0100)plt.figure(figsize(9, 7))plt.subplot(231), plt.title(origin), plt.axis(off), plt.imshow(img, gray)plt.subplot(232), plt.title(turbulence blur), plt.axis(off), plt.imshow(imgBlur, gray)plt.subplot(233), plt.title(inverse filter(D0full)), plt.axis(off), plt.imshow(imgRebuild, gray)plt.subplot(234), plt.title(inverse filter(D040)), plt.axis(off), plt.imshow(imgRebuild1, gray)plt.subplot(235), plt.title(inverse filter(D070)), plt.axis(off), plt.imshow(imgRebuild2, gray)plt.subplot(236), plt.title(inverse filter(D0100)), plt.axis(off), plt.imshow(imgRebuild3, gray)plt.tight_layout()plt.show() 6.2 退化图像的维纳滤波Wiener filter
逆滤波对加性噪声特别敏感使得恢复的图像几乎不可用例程9.20中的退化图像未加入噪声项 。
最小均方误差滤波用来去除含有噪声的模糊图像其目标是寻找未污染图像的一个估计使均方误差最小 e2E{(f−f^)2}e^2 E\{ (f - \hat{f})^2 \} e2E{(f−f^)2} 最小均方差滤波由 Wiener [1942] 首先提出是最早提出的线性图像复原方法因此称为维纳滤波。
信噪比 SNR(f)S(f)/N(f)SNR(f) S(f)/N(f)SNR(f)S(f)/N(f)是信息承载信号功率未退化的原图像水平与噪声功率水平的测度。低噪声图像的信噪比较高高噪声图像的信噪比较低。
将复原后的图像视为“信号”而复原图像与原图像的差视为“噪声”则可以将空间域的信噪比定义为 SNR∑x0M−1∑y0N−1f^(x,y)2/∑x0M−1∑y0N−1[f(x,y)−f^(x,y)]2SNR \sum_{x0}^{M-1} \sum_{y0}^{N-1} \hat{f}(x,y)^2 / \sum_{x0}^{M-1} \sum_{y0}^{N-1} [f(x,y) -\hat{f}(x,y)]^2 SNRx0∑M−1y0∑N−1f^(x,y)2/x0∑M−1y0∑N−1[f(x,y)−f^(x,y)]2 最小均方误差滤波的传递函数描述为 G(f)H∗(f)S(f)∣H(f)∣2S(f)N(f)H∗(f)∣H(f)∣2S(f)N(f)/S(f)G(f) \frac{H^* (f) S(f)}{|H(f)|^2 S(f) N(f)} \frac{H^* (f)}{|H(f)|^2 S(f) N(f)/S(f)} G(f)∣H(f)∣2S(f)N(f)H∗(f)S(f)∣H(f)∣2S(f)N(f)/S(f)H∗(f) 式中 G(f)、H(f)G(f)、H(f)G(f)、H(f) 是 g 和 h 在频率域的傅里叶变换S(f)、N(f)S(f)、N(f)S(f)、N(f) 是信号 x(t)、噪声 n(t) 的功率谱。
用信噪比将上式进一步表示为 G(f)1H(f)∣H(f)∣2∣H(f)∣21/SNR(f)G(f) \frac{1}{H(f)} \frac{|H(f)|^2}{|H(f)|^2 1/SNR(f)} G(f)H(f)1∣H(f)∣21/SNR(f)∣H(f)∣2 当处理白噪声时∣N(u,v)∣2|N(u,v)|^2∣N(u,v)∣2 是一个常数 但未退化图像和噪声的功率谱通常未知或不能估计则可用下式近似 F^(u,v)1H(u,v)∣H(u,v)∣2∣H(u,v)∣2KG(u,v)\hat{F}(u,v) \frac{1}{H(u,v)} \frac{|H(u,v)|^2}{|H(u,v)|^2 K} G(u,v) F^(u,v)H(u,v)1∣H(u,v)∣2K∣H(u,v)∣2G(u,v) K 是加到 ∣H(u,v)∣2|H(u,v)|^2∣H(u,v)∣2的所有项上的一个规定常数。 例程 9.21: 维纳滤波 Wiener filter # 9.21: 退化图像的维纳滤波 (Wiener filter)def getMotionDsf(shape, angle, dist):xCenter (shape[0] - 1) / 2yCenter (shape[1] - 1) / 2sinVal np.sin(angle * np.pi / 180)cosVal np.cos(angle * np.pi / 180)PSF np.zeros(shape) # 点扩散函数for i in range(dist): # 将对应角度上motion_dis个点置成1xOffset round(sinVal * i)yOffset round(cosVal * i)PSF[int(xCenter - xOffset), int(yCenter yOffset)] 1return PSF / PSF.sum() # 归一化def makeBlurred(image, PSF, eps): # 对图片进行运动模糊fftImg np.fft.fft2(image) # 进行二维数组的傅里叶变换fftPSF np.fft.fft2(PSF) epsfftBlur np.fft.ifft2(fftImg * fftPSF)fftBlur np.abs(np.fft.fftshift(fftBlur))return fftBlurdef inverseFilter(image, PSF, eps): # 逆滤波fftImg np.fft.fft2(image)fftPSF np.fft.fft2(PSF) eps # 噪声功率这是已知的考虑epsilonimgInvFilter np.fft.ifft2(fftImg / fftPSF) # 计算F(u,v)的傅里叶反变换imgInvFilter np.abs(np.fft.fftshift(imgInvFilter))return imgInvFilterdef wienerFilter(input, PSF, eps, K0.01): # 维纳滤波K0.01fftImg np.fft.fft2(input)fftPSF np.fft.fft2(PSF) epsfftWiener np.conj(fftPSF) / (np.abs(fftPSF)**2 K)imgWienerFilter np.fft.ifft2(fftImg * fftWiener)imgWienerFilter np.abs(np.fft.fftshift(imgWienerFilter))return imgWienerFilter# 读取原始图像img cv2.imread(../images/Fig0526a.tif, 0) # flags0 读取为灰度图像hImg, wImg img.shape[:2]# 不含噪声的运动模糊PSF getMotionDsf((hImg, wImg), 45, 100) # 运动模糊函数imgBlurred np.abs(makeBlurred(img, PSF, 1e-6)) # 生成不含噪声的运动模糊图像imgInvFilter inverseFilter(imgBlurred, PSF, 1e-6) # 逆滤波imgWienerFilter wienerFilter(imgBlurred, PSF, 1e-6) # 维纳滤波# 带有噪声的运动模糊scale 0.05 # 噪声方差noisy imgBlurred.std() * np.random.normal(loc0.0, scalescale, sizeimgBlurred.shape) # 添加高斯噪声imgBlurNoisy imgBlurred noisy # 带有噪声的运动模糊imgNoisyInv inverseFilter(imgBlurNoisy, PSF, scale) # 对添加噪声的模糊图像进行逆滤波imgNoisyWiener wienerFilter(imgBlurNoisy, PSF, scale) # 对添加噪声的模糊图像进行维纳滤波plt.figure(figsize(9, 7))plt.subplot(231), plt.title(blurred image), plt.axis(off), plt.imshow(imgBlurred, gray)plt.subplot(232), plt.title(inverse filter), plt.axis(off), plt.imshow(imgInvFilter, gray)plt.subplot(233), plt.title(Wiener filter), plt.axis(off), plt.imshow(imgWienerFilter, gray)plt.subplot(234), plt.title(blurred image with noisy), plt.axis(off), plt.imshow(imgBlurNoisy, gray)plt.subplot(235), plt.title(inverse filter), plt.axis(off), plt.imshow(imgNoisyInv, gray)plt.subplot(236), plt.title(Wiener filter), plt.axis(off), plt.imshow(imgNoisyWiener, gray)plt.tight_layout()plt.show()程序说明
对于不含噪声的运动模糊图像在已知运动模糊退化模型和参数的前提下使用逆滤波可以很好地复原退化图像逆滤波的性能优于维纳滤波。但是考虑实际退化图像往往含有一定水平的加性噪声此时即使已知退化模型逆滤波的后的噪声几乎掩盖了图像内容而维纳滤波的结果则较好。 6.3 约束最小二乘方滤波Constrained Least Squares Filtering
维纳滤波建立在退化函数和信噪比已知的前提下这在实践中并不容易满足。
约束最小二乘方滤波仅要求噪声方差和均值的知识或估计这些参数通常可以由一幅给定的退化图像算出因而具有更为广泛的应用。而且维纳滤波是以最小化一个统计准则为基础因此是平均意义上的最优而约束最小二乘方滤波则对每幅图像都会产生最优估计。
与维纳滤波相比最小二乘方滤波对于高噪声和中等噪声处理效果更好对于低噪声二者处理结果基本相同。当手工选择参数以取得更好的视觉效果时约束最小二乘方滤波的效果有可能比维纳滤波的效果更好。
约束最小二乘方滤波的核心是解决退化函数对噪声的敏感性问题。降低噪声敏感的一种方法是以平滑度量的最佳复原为基础如图像的二阶导数拉普拉斯算子其数学描述是求一个带约束条件的准则函数的最小值 minC∑x0M−1∑y0N−1[▽2f(x,y)]2s.t.:∣∣g−Hf^∣∣2∣∣η∣∣2min \ C \sum_{x0}^{M-1} \sum_{y0}^{N-1} [\triangledown ^2 f(x,y)]^2 \\ s.t.: ||g - H \hat{f}||^2 ||\eta||^2 min Cx0∑M−1y0∑N−1[▽2f(x,y)]2s.t.:∣∣g−Hf^∣∣2∣∣η∣∣2 其中∣∣a∣∣2||a||^2∣∣a∣∣2是欧几里德范数f^\hat{f}f^ 是未退化图像的估计$ \triangledown ^2 $ 是拉普拉斯算子。求准则函数 C 的最小值便得到最好的平滑效果即退化图像的最佳复原。
该约束最小二乘方问题的解是 F^(u,v)[H∗(u,v)∣H(u,v)∣2γ∣P(u,v)∣2]G(u,v)\hat{F}(u,v) [\frac{H^*(u,v)}{|H(u,v)|^2 \gamma |P(u,v)|^2}] \ G(u,v) F^(u,v)[∣H(u,v)∣2γ∣P(u,v)∣2H∗(u,v)] G(u,v) 式中γ\gammaγ 是一个必须满足约束条件的参数P(u,v)P(u,v)P(u,v) 是函数 p(x,y)p(x,y)p(x,y) 的傅里叶变换 p(x,y)[0−10−14−10−10]p(x,y) \begin{bmatrix} 0 -1 0 \\ -1 4 -1\\ 0 -1 0 \end{bmatrix} p(x,y)⎣⎡0−10−14−10−10⎦⎤ 注意函数 P(u,v)P(u,v)P(u,v) 和 H(u,v)H(u,v)H(u,v) 的大小必须相等。如果 H(u,v)H(u,v)H(u,v) 的大小为 M∗NM*NM∗N则 p(x,y) 必须嵌入 M∗NM*NM∗N 零阵列的中心并保持偶对称。 例程 9.21: 约束最小二乘方滤波 # 9.21: 约束最小二乘方滤波def getMotionDsf(shape, angle, dist):xCenter (shape[0] - 1) / 2yCenter (shape[1] - 1) / 2sinVal np.sin(angle * np.pi / 180)cosVal np.cos(angle * np.pi / 180)PSF np.zeros(shape) # 点扩散函数for i in range(dist): # 将对应角度上motion_dis个点置成1xOffset round(sinVal * i)yOffset round(cosVal * i)PSF[int(xCenter - xOffset), int(yCenter yOffset)] 1return PSF / PSF.sum() # 归一化def makeBlurred(image, PSF, eps): # 对图片进行运动模糊fftImg np.fft.fft2(image) # 进行二维数组的傅里叶变换fftPSF np.fft.fft2(PSF) epsfftBlur np.fft.ifft2(fftImg * fftPSF)fftBlur np.abs(np.fft.fftshift(fftBlur))return fftBlurdef wienerFilter(input, PSF, eps, K0.01): # 维纳滤波K0.01fftImg np.fft.fft2(input)fftPSF np.fft.fft2(PSF) epsfftWiener np.conj(fftPSF) / (np.abs(fftPSF)**2 K)imgWienerFilter np.fft.ifft2(fftImg * fftWiener)imgWienerFilter np.abs(np.fft.fftshift(imgWienerFilter))return imgWienerFilterdef getPuv(image):h, w image.shape[:2]hPad, wPad h - 3, w - 3pxy np.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]])pxyPad np.pad(pxy, ((hPad//2, hPad - hPad//2), (wPad//2, wPad - wPad//2)), modeconstant)fftPuv np.fft.fft2(pxyPad)return fftPuvdef leastSquareFilter(image, PSF, eps, gamma0.01): # 约束最小二乘方滤波fftImg np.fft.fft2(image)fftPSF np.fft.fft2(PSF)conj fftPSF.conj()fftPuv getPuv(image)# absConj np.abs(fftPSF) ** 2Huv conj / (np.abs(fftPSF)**2 gamma * (np.abs(fftPuv)**2))ifftImg np.fft.ifft2(fftImg * Huv)ifftShift np.abs(np.fft.fftshift(ifftImg))imgLSFilter np.uint8(cv2.normalize(np.abs(ifftShift), None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]return imgLSFilter# # 读取原始图像img cv2.imread(../images/Fig0526a.tif, 0) # flags0 读取为灰度图像hImg, wImg img.shape[:2]# 带有噪声的运动模糊PSF getMotionDsf((hImg, wImg), 45, 100) # 运动模糊函数imgBlurred np.abs(makeBlurred(img, PSF, 1e-6)) # 生成不含噪声的运动模糊图像scale 0.01 # 噪声方差noisy imgBlurred.std() * np.random.normal(loc0.0, scalescale, sizeimgBlurred.shape) # 添加高斯噪声imgBlurNoisy imgBlurred noisy # 带有噪声的运动模糊imgWienerFilter wienerFilter(imgBlurNoisy, PSF, scale, K0.01) # 对含有噪声的模糊图像进行维纳滤波imgLSFilter leastSquareFilter(imgBlurNoisy, PSF, scale, gamma0.01) # 约束最小二乘方滤波plt.figure(figsize(9, 7))plt.subplot(231), plt.title(blurred image (dev0.01)), plt.axis(off), plt.imshow(imgBlurNoisy, gray)plt.subplot(232), plt.title(Wiener filter), plt.axis(off), plt.imshow(imgWienerFilter, gray)plt.subplot(233), plt.title(least square filter), plt.axis(off), plt.imshow(imgLSFilter, gray)scale 0.1 # 噪声方差noisy imgBlurred.std() * np.random.normal(loc0.0, scalescale, sizeimgBlurred.shape) # 添加高斯噪声imgBlurNoisy imgBlurred noisy # 带有噪声的运动模糊imgWienerFilter wienerFilter(imgBlurNoisy, PSF, scale, K0.01) # 维纳滤波imgLSFilter leastSquareFilter(imgBlurNoisy, PSF, scale, gamma0.1) # 约束最小二乘方滤波plt.subplot(234), plt.title(blurred image (dev0.1)), plt.axis(off), plt.imshow(imgBlurNoisy, gray)plt.subplot(235), plt.title(Wiener filter), plt.axis(off), plt.imshow(imgWienerFilter, gray)plt.subplot(236), plt.title(least square filter), plt.axis(off), plt.imshow(imgLSFilter, gray)plt.tight_layout()plt.show()程序说明
对于含有不同水平噪声的运动模糊图像最小二乘方滤波的图像复原效果都较好特别是对于高噪声和中等噪声的处理效果更好。
对于最小二乘方滤波可以交互地或循环地调整参数 γ\gammaγ 以取得更好的复原效果例如参数 γ\gammaγ 可以从较小的取值递增以获得最优估计。 6.4 几何均值滤波
几何均值滤波器是维纳滤波的推广其传递函数由括号内幂次分别为 α\alphaα 和 1−α1-\alpha1−α 的两个表达式组成 F^(u,v)[H∗(u,v)∣H(u,v)∣2]α[H∗(u,v)∣H(u,v)∣2β[Sη(u,v)/Sf(u,v)]]1−αG(u,v)\hat{F}(u,v) \Bigg[ \frac{H^*(u,v)}{|H(u,v)|^2} \Bigg] ^{\alpha} \ \Bigg[\frac{H^*(u,v)}{|H(u,v)|^2 \beta [S_{\eta} (u,v) /S_f(u,v)]}\Bigg]^{1-\alpha}\ G(u,v) F^(u,v)[∣H(u,v)∣2H∗(u,v)]α [∣H(u,v)∣2β[Sη(u,v)/Sf(u,v)]H∗(u,v)]1−α G(u,v) 式中$ \alpha$ 和 β\betaβ 是非负的实常数。
当 α1/2\alpha1/2α1/2 时几何均值滤波器是幂次相同的两个量的乘积这就是几何均值的含义。
当 α1/2,β1\alpha1 / 2, \beta1α1/2,β1 时称为频谱均衡滤波器。当 α1\alpha1α1 时简化为逆滤波器当 α0\alpha0α0 时简化为带参数的维纳滤波器并在 β1\beta1β1 时成为标准维纳滤波器。 例程 9.22: 几何均值滤波 # 9.22: 几何均值滤波器def getMotionDsf(shape, angle, dist):xCenter (shape[0] - 1) / 2yCenter (shape[1] - 1) / 2sinVal np.sin(angle * np.pi / 180)cosVal np.cos(angle * np.pi / 180)PSF np.zeros(shape) # 点扩散函数for i in range(dist): # 将对应角度上motion_dis个点置成1xOffset round(sinVal * i)yOffset round(cosVal * i)PSF[int(xCenter - xOffset), int(yCenter yOffset)] 1return PSF / PSF.sum() # 归一化def makeBlurred(image, PSF, eps): # 对图片进行运动模糊fftImg np.fft.fft2(image) # 进行二维数组的傅里叶变换fftPSF np.fft.fft2(PSF) epsfftBlur np.fft.ifft2(fftImg * fftPSF)fftBlur np.abs(np.fft.fftshift(fftBlur))return fftBlurdef wienerFilter(input, PSF, eps, K0.01): # 维纳滤波K0.01fftImg np.fft.fft2(input)fftPSF np.fft.fft2(PSF) epsfftWiener np.conj(fftPSF) / (np.abs(fftPSF)**2 K)imgWienerFilter np.fft.ifft2(fftImg * fftWiener)imgWienerFilter np.abs(np.fft.fftshift(imgWienerFilter))return imgWienerFilterdef geometricMeanFilter(image, PSF, eps, K1, alpha1, beta1): # 几何均值滤波器fftImg np.fft.fft2(image)fftPSF np.fft.fft2(PSF)conj fftPSF.conj()squarePSF (fftPSF * conj).realHuv np.power(conj / (squarePSF), alpha) * np.power(conj / (squarePSF beta * K), 1-alpha)ifftImg np.fft.ifft2(fftImg * Huv)ifftShift np.abs(np.fft.fftshift(ifftImg))imgGMFilter np.uint8(cv2.normalize(np.abs(ifftShift), None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]return imgGMFilter# # 读取原始图像img cv2.imread(../images/Fig0526a.tif, 0) # flags0 读取为灰度图像hImg, wImg img.shape[:2]# 带有噪声的运动模糊PSF getMotionDsf((hImg, wImg), 45, 100) # 运动模糊函数imgBlurred np.abs(makeBlurred(img, PSF, 1e-6)) # 生成不含噪声的运动模糊图像scale 0.01 # 噪声方差noisy imgBlurred.std() * np.random.normal(loc0.0, scalescale, sizeimgBlurred.shape) # 添加高斯噪声imgBlurNoisy imgBlurred noisy # 带有噪声的运动模糊imgWienerFilter wienerFilter(imgBlurNoisy, PSF, scale, K0.01) # 对含有噪声的模糊图像进行维纳滤波imgGMFilter geometricMeanFilter(imgBlurNoisy, PSF, scale, K0.01, alpha0.5, beta1) # 约束最小二乘方滤波plt.figure(figsize(9, 7))plt.subplot(231), plt.title(blurred image (dev0.01)), plt.axis(off), plt.imshow(imgBlurNoisy, gray)plt.subplot(232), plt.title(Wiener filter), plt.axis(off), plt.imshow(imgWienerFilter, gray)plt.subplot(233), plt.title(geometric mean filter), plt.axis(off), plt.imshow(imgGMFilter, gray)scale 0.1 # 噪声方差noisy imgBlurred.std() * np.random.normal(loc0.0, scalescale, sizeimgBlurred.shape) # 添加高斯噪声imgBlurNoisy imgBlurred noisy # 带有噪声的运动模糊imgWienerFilter wienerFilter(imgBlurNoisy, PSF, scale, K0.01) # 维纳滤波imgGMFilter geometricMeanFilter(imgBlurNoisy, PSF, scale, K0.01, alpha0, beta1) # 约束最小二乘方滤波plt.subplot(234), plt.title(blurred image (dev0.1)), plt.axis(off), plt.imshow(imgBlurNoisy, gray)plt.subplot(235), plt.title(Wiener filter), plt.axis(off), plt.imshow(imgWienerFilter, gray)plt.subplot(236), plt.title(geometric mean filter), plt.axis(off), plt.imshow(imgGMFilter, gray)plt.tight_layout()plt.show()程序说明 7. 投影重建图像
7.1 计算机断层成像Computed tomography
图像重建Image Reconstruction的基本思想就是通过探测物体的投影数据重建物体的实际内部构造。
断层成像就是要获得物体内部的截面图像。利用 X射线、超声波等射线穿透被遮挡物体的透视投影图可以计算恢复物体的断层图进而可以利用断层图或直接从二维透视投影图重建物体的形状和内部结构。其原理是射线在穿过不同组织时的吸收率不同在成像面上得到不同的投射强度由此反演求得内部分布的图像。
X 射线、CT 技术就是应用断层重建的医学诊断方法。计算机断层扫描CT已经发展成为一种非常成功和必不可少的医学诊断工具被认为是 X 射线发现以来医学影像领域最伟大的发明。
X 射线计算机断层成像的基本原理是使用 X 射线从多个不同方向和角度对物体进行扫描通过反投影算法获取物体内部结构的切片堆叠这些切片就可以得到人体的三维表示。
投影重建还应用于地矿探测接收不同地层和矿体反射的超声波 按照超声波在媒质的透射率和反射规律对透射投影图进行分析计算就可以恢复重建地下的矿体形状。 7.2 投影和雷登变换Radon transform
雷登变换是三维重建的数学基础。
一条直线 yaxbyaxbyaxb 的法线束为 xcosθysinθρx cos \theta y sin \theta \rho xcosθysinθρ 沿该直线对函数 f(x,y)f(x,y)f(x,y) 进行积分这个积分结果就是雷登变换后的函数 Rf\mathcal{R}fRf 在这条直线上的值 g(ρ,θ)∫−∞∞∫−∞∞f(x,y)δ(xcosθysinθ−ρ)dxdyg(\rho, \theta) \int ^{\infty}_{-\infty} \int ^{\infty}_{-\infty} f(x,y) \ \delta(x cos \theta y sin \theta - \rho) dxdy g(ρ,θ)∫−∞∞∫−∞∞f(x,y) δ(xcosθysinθ−ρ)dxdy 雷登变换的离散形式 g(ρ,θ)∑x0M−1∑y0N−1f(x,y)δ(xcosθysinθ−ρ)g(\rho, \theta) \sum ^{M-1}_{x0} \sum^{N-1}_{y0} f(x,y) \ \delta(x cos \theta y sin \theta - \rho) g(ρ,θ)x0∑M−1y0∑N−1f(x,y) δ(xcosθysinθ−ρ)
雷登变换将二维空间 xyxyxy 映射到了另一个二维空间 αs\alpha sαs称为直线空间。
雷登变换常被用于医学影像处理利用反雷登变换可以进行三维图像重建。 例程 9.23雷登变换正弦图 # # 9.23: 雷登变换正弦图from scipy import ndimagedef discreteRadonTransform(image, steps):channels image.shape[0]res np.zeros((channels, channels), dtypenp.float32)for s in range(steps):rotation ndimage.rotate(image, -s * 180/steps, reshapeFalse).astype(np.float32)res[:, s] sum(rotation)return res# # 读取原始图像plt.figure(figsize(9, 7))fileImg [Fig0534a.tif, Fig0534b.tif, Fig0534c.tif] # 原始图像 文件名for i in range(len(fileImg)):img cv2.imread(../images/fileImg[i], flags0) # flags0 读取为灰度图像imgRadon discreteRadonTransform(img, img.shape[0]) # Radon 变换print(img.shape, imgRadon.shape)plt.subplot(2, len(fileImg), i 1), plt.axis(off) # 绘制原始图像plt.title(origin image), plt.imshow(img, gray)plt.subplot(2, len(fileImg), i len(fileImg) 1), plt.axis(off) # 绘制 sinogram 图plt.title(Radon transform), plt.imshow(imgRadon, gray)plt.tight_layout()plt.show() 7.2 雷登变换反投影重建图像Radon transform back projection
空间点 X 通过摄像机 P 被作用到图像平面的图像点 m PX 获得采集到的图像数据这种投影关系称为摄像机的正向投影 forward projection 简称投影。
反向投影是针对图像平面的基本几何元素而言的图像平面点 m 的反投影是指在摄像机 P 的作用下具有像点 m 的所有空间点的集合。
反投影一个点形成部分图像的过程是将直线 L(ρj,θk)L(\rho_j,\theta_k)L(ρj,θk) 复制到图像上直线上每点的灰度值是 g(ρj,θk)g(\rho_j,\theta_k)g(ρj,θk)。遍历投影信号中的每个点得到 fθ(x,y)g(ρ,θ)g(xcosθysinθ)f_{\theta}(x,y) g(\rho,\theta) g(x cos \theta ysin \theta) fθ(x,y)g(ρ,θ)g(xcosθysinθ) 对所有反投影图像积分得到最终的图像 f(x,y)∫0πfθ(x,y)dθf(x,y) \int_0^{\pi}f_{\theta}(x,y)d\theta f(x,y)∫0πfθ(x,y)dθ 其离散形式为
f(x,y)∑θ0πfθ(x,y)f(x,y) \sum_{\theta0}^{\pi}f_{\theta}(x,y) f(x,y)θ0∑πfθ(x,y)
反投影的图像有时称为层图可以理解为投影图像的一个近似。 例程 9.24雷登变换反投影重建图像 # 9.24: 雷登变换反投影重建图像from scipy import ndimagedef discreteRadonTransform(image, steps): # 离散雷登变换channels image.shape[0]resRadon np.zeros((channels, channels), dtypenp.float32)for s in range(steps):rotation ndimage.rotate(image, -s * 180/steps, reshapeFalse).astype(np.float32)resRadon[:, s] sum(rotation)return resRadondef inverseRadonTransform(image, steps): # 雷登变换反投影channels image.shape[0]res np.zeros((steps, channels, channels))for s in range(steps):expandDims np.expand_dims(image[:, s], axis0)repeat expandDims.repeat(channels, axis0)res[s] ndimage.rotate(repeat, s * 180/steps, reshapeFalse).astype(np.float32)invRadon np.sum(res, axis0)return invRadon# 读取原始图像img1 cv2.imread(../images/Fig0534a.tif, 0) # flags0 读取为灰度图像img2 cv2.imread(../images/Fig0534c.tif, 0)# 雷登变换imgRadon1 discreteRadonTransform(img1, img1.shape[0]) # Radon 变换imgRadon2 discreteRadonTransform(img2, img2.shape[0])# 雷登变换反投影imgInvRadon1 inverseRadonTransform(imgRadon1, imgRadon1.shape[0])imgInvRadon2 inverseRadonTransform(imgRadon2, imgRadon2.shape[0])plt.figure(figsize(9, 7))plt.subplot(231), plt.axis(off), plt.title(origin image), plt.imshow(img1, gray) # 绘制原始图像plt.subplot(232), plt.axis(off), plt.title(Radon transform), plt.imshow(imgRadon1, gray) # 绘制 sinogram 图plt.subplot(233), plt.axis(off), plt.title(Inv-Radon transform), plt.imshow(imgInvRadon1, gray)plt.subplot(234), plt.axis(off), plt.title(origin image), plt.imshow(img2, gray)plt.subplot(235), plt.axis(off), plt.title(Radon transform), plt.imshow(imgRadon2, gray)plt.subplot(236), plt.axis(off), plt.title(Inv-Radon transform), plt.imshow(imgInvRadon2, gray)plt.tight_layout()plt.show()7.3 滤波反投影重建图像 Filtered back projection
直接由正弦图得到反投影图像会存在严重的模糊这是早期 CT 系统所存在的问题。
傅立叶中心切片定理表明投影的一维傅立叶变换是得到投影区域的二维傅立叶变换的切片。滤波反投影重建算法在反投影前将每一个采集投影角度下的投影进行卷积处理从而改善点扩散函数引起的形状伪影有效地改善了重建的图像质量。
以 f(x,y)f(x,y)f(x,y) 表示需要重建的图像F(u,v)F(u,v)F(u,v) 的二维傅里叶反变换是 f(x,y)∫−∞∞∫−∞∞F(u,v)ej2π(uxvy)dudvf(x,y) \int^{\infty}_{-\infty} \int^{\infty}_{-\infty} F(u,v) e^{j 2 \pi (uxvy)} dudv f(x,y)∫−∞∞∫−∞∞F(u,v)ej2π(uxvy)dudv 利用傅里叶切片定理可以得到 f(x,y)∫0π[∫−∞∞∣ω∣G(ω,θ)ej2πωρdω]dθf(x,y) \int^{\pi}_{0} \big[ \int^{\infty}_{-\infty} |\omega|G(\omega,\theta) e^{j 2 \pi \omega \rho} d\omega \big] d\theta f(x,y)∫0π[∫−∞∞∣ω∣G(ω,θ)ej2πωρdω]dθ 括号 [] 内部是一个一维傅里叶反变换可以认为这是一个一维滤波器的传递函数。由于 ∣ω∣|\omega|∣ω∣ 是一个不可积的斜坡函数Slope function可以通过对斜坡加窗进行限制典型地如汉明窗Hamming window、韩窗Hann window。
该式也可以使用空间卷积来实现 f(x,y)∫0π[∫−∞∞∣ω∣G(ω,θ)ej2πωρdω]dθ∫0π[∫−∞∞g(ρ,θ)s(xcosθysinθ−ρ)dρ]dθf(x,y) \int^{\pi}_{0} \big[ \int^{\infty}_{-\infty} |\omega|G(\omega,\theta) e^{j 2 \pi \omega \rho} d\omega \big] d\theta \\ \int^{\pi}_{0} \big[ \int^{\infty}_{-\infty} g(\rho,\theta) s(xcos \theta ysin \theta - \rho) d\rho \big] d\theta \\ f(x,y)∫0π[∫−∞∞∣ω∣G(ω,θ)ej2πωρdω]dθ∫0π[∫−∞∞g(ρ,θ)s(xcosθysinθ−ρ)dρ]dθ 这表明将对应的投影 g(ρ,θ)g(\rho, \theta)g(ρ,θ) 与斜坡滤波器传递函数 s(ρ)s(\rho)s(ρ) 的傅里叶反变换进行卷积可以得到角度 θ\thetaθ 的各个反投影整个反投影图像可以通过对所有反投影图像积分得到。 例程 9.25滤波反投影重建图像
SL 滤波器是使用 Sinc 函数对斜坡滤波器进行截断产生其离散形式为 hSL(nδ)−2π2δ2(4∗n2−1),n0,±1,±2,...h_{SL}(n \delta) - \frac{2}{\pi ^2 \delta ^2 (4*n^2 -1)}, \ n0,\pm 1,\pm 2,... hSL(nδ)−π2δ2(4∗n2−1)2, n0,±1,±2,... 利用投影值和 SL 滤波器进行卷积然后再进行反投影就可以实现图像重建。 # 9.25: 滤波反投影重建图像from scipy import ndimagedef discreteRadonTransform(image, steps): # 离散雷登变换channels image.shape[0]resRadon np.zeros((channels, channels), dtypenp.float32)for s in range(steps):rotation ndimage.rotate(image, -s * 180/steps, reshapeFalse).astype(np.float32)resRadon[:, s] sum(rotation)return resRadondef inverseRadonTransform(image, steps): # 雷登变换反投影重建图像channels image.shape[0]backproject np.zeros((steps, channels, channels))for s in range(steps):expandDims np.expand_dims(image[:, s], axis0)repeat expandDims.repeat(channels, axis0)backproject[s] ndimage.rotate(repeat, s * 180/steps, reshapeFalse).astype(np.float32)invRadon np.sum(backproject, axis0)return invRadondef SLFilter(N, d): # SL 滤波器, Sinc 函数对斜坡滤波器进行截断filterSL np.zeros((N,))for i in range(N):filterSL[i] - 2 / (np.pi**2 * d**2 * (4 * (i-N/2)**2 - 1))return filterSLdef filterInvRadonTransform(image, steps): # 滤波反投影重建图像channels len(image[0])backproject np.zeros((steps, channels, channels)) # 反投影filter SLFilter(channels, 1) # SL 滤波器for s in range(steps):value image[:, s] # 投影值valueFiltered np.convolve(filter, value, same) # 投影值和 SL 滤波器进行卷积filterExpandDims np.expand_dims(valueFiltered, axis0)filterRepeat filterExpandDims.repeat(channels, axis0)backproject[s] ndimage.rotate(filterRepeat, s * 180/steps, reshapeFalse).astype(np.float32)filterInvRadon np.sum(backproject, axis0)return filterInvRadon# 读取原始图像img1 cv2.imread(../images/Fig0534a.tif, 0) # flags0 读取为灰度图像img2 cv2.imread(../images/Fig0534c.tif, 0)# 雷登变换imgRadon1 discreteRadonTransform(img1, img1.shape[0]) # Radon 变换imgRadon2 discreteRadonTransform(img2, img2.shape[0])# 雷登变换反投影重建图像imgInvRadon1 inverseRadonTransform(imgRadon1, imgRadon1.shape[0])imgInvRadon2 inverseRadonTransform(imgRadon2, imgRadon2.shape[0])# 滤波反投影重建图像imgFilterInvRadon1 filterInvRadonTransform(imgRadon1, imgRadon1.shape[0])imgFilterInvRadon2 filterInvRadonTransform(imgRadon2, imgRadon2.shape[0])plt.figure(figsize(9, 7))plt.subplot(231), plt.axis(off), plt.title(origin image), plt.imshow(img1, gray) # 绘制原始图像plt.subplot(232), plt.axis(off), plt.title(inverse Radon trans), plt.imshow(imgInvRadon1, gray)plt.subplot(233), plt.axis(off), plt.title(filter inv-Radon trans), plt.imshow(imgFilterInvRadon1, gray)plt.subplot(234), plt.axis(off), plt.title(origin image), plt.imshow(img2, gray)plt.subplot(235), plt.axis(off), plt.title(inverse Radon trans), plt.imshow(imgInvRadon2, gray)plt.subplot(236), plt.axis(off), plt.title(filter inv-Radon trans), plt.imshow(imgFilterInvRadon2, gray)plt.tight_layout()plt.show()版权声明
注本文中部分原始图片来自 Rafael C. Gonzalez “Digital Image Processing, 4th.Ed.”特此致谢。
版权声明 youcansxupt 原创作品转载必须标注原文链接(https://blog.csdn.net/youcans/article/details/123192478) Copyright 2022 youcans, XUPT 欢迎关注 『youcans 的 OpenCV 学习课』 系列持续更新 【youcans 的 OpenCV 学习课】1. 安装与环境配置 【youcans 的 OpenCV 学习课】2. 图像读取与显示 【youcans 的 OpenCV 学习课】3. 图像的创建与修改 【youcans 的 OpenCV 学习课】4. 图像的叠加与混合 【youcans 的 OpenCV 学习课】5. 图像的几何变换 【youcans 的 OpenCV 学习课】6. 灰度变换与直方图处理 【youcans 的 OpenCV 学习课】7. 空间域图像滤波 【youcans 的 OpenCV 学习课】8. 频率域图像滤波上 【youcans 的 OpenCV 学习课】9. 频率域图像滤波下 【youcans 的 OpenCV 学习课】10. 图像复原与重建
