网站做缓存吗,江西最新新闻事件今天,王也为什么这么受欢迎,wordpress发表的文章百度抓取失败目录使用高通滤波器锐化图像由低通滤波器得到理想、高斯和巴特沃斯高通滤波器指纹增强频域中的拉普拉斯钝化掩蔽、高提升滤波和高频强调滤波同态滤波使用高通滤波器锐化图像
由低通滤波器得到理想、高斯和巴特沃斯高通滤波器 HHP(u,v)1−HLP(u,v)(4.118)H_{HP}(u, v) 1 - H_{…
目录使用高通滤波器锐化图像由低通滤波器得到理想、高斯和巴特沃斯高通滤波器指纹增强频域中的拉普拉斯钝化掩蔽、高提升滤波和高频强调滤波同态滤波使用高通滤波器锐化图像
由低通滤波器得到理想、高斯和巴特沃斯高通滤波器
HHP(u,v)1−HLP(u,v)(4.118)H_{HP}(u, v) 1 - H_{LP}(u, v) \tag{4.118}HHP(u,v)1−HLP(u,v)(4.118) 理想高通 H(u,v){0,D(u,v)≤D01,D(u,v)D0(4.119)H(u, v) \begin{cases} 0, D(u, v) \leq D_0 \\1, D(u, v) D_0\end{cases} \tag{4.119}H(u,v){0,1,D(u,v)≤D0D(u,v)D0(4.119) 高斯高通 H(u,v)1−e−D2(u,v)/2D02(4.120)H(u,v) 1 - e^{-D^2(u,v) / 2D_0^2} \tag{4.120}H(u,v)1−e−D2(u,v)/2D02(4.120) 巴特沃斯高通 H(u,v)11[D0/D(u,v)]2n(4.121)H(u,v) \frac{1} {1 [D_0 / D(u,v)]^{2n}} \tag{4.121}H(u,v)1[D0/D(u,v)]2n1(4.121)
def idea_high_pass_filter(source, center, radius5):create idea high pass filter param: source: input, source imageparam: center: input, the center of the filter, where is the lowest value, (0, 0) is top left corner, source.shape[:2] is center of the source imageparam: radius: input, the radius of the lowest value, greater value, bigger blocker out range, if the radius is 0, then allvalue is 0return a [0, 1] value filterM, N source.shape[1], source.shape[0]u np.arange(M)v np.arange(N)u, v np.meshgrid(u, v)D np.sqrt((u - center[1]//2)**2 (v - center[0]//2)**2)D0 radiuskernel D.copy()kernel[D D0] 1kernel[D D0] 0return kerneldef gauss_high_pass_filter(source, center, radius5):create gaussian high pass filter param: source: input, source imageparam: center: input, the center of the filter, where is the lowest value, (0, 0) is top left corner, source.shape[:2] is center of the source imageparam: radius: input, the radius of the lowest value, greater value, bigger blocker out range, if the radius is 0, then allvalue is 1return a [0, 1] value filterM, N source.shape[1], source.shape[0]u np.arange(M)v np.arange(N)u, v np.meshgrid(u, v)D np.sqrt((u - center[1]//2)**2 (v - center[0]//2)**2)D0 radiusif D0 0:kernel np.ones(source.shape[:2], dtypenp.float)else:kernel 1 - np.exp(- (D**2)/(D0**2)) return kerneldef butterworth_high_pass_filter(img, center, radius5, n1):create butterworth high pass filter param: source: input, source imageparam: center: input, the center of the filter, where is the lowest value, (0, 0) is top left corner, source.shape[:2] is center of the source imageparam: radius: input, the radius of the lowest value, greater value, bigger blocker out range, if the radius is 0, then allvalue is 0param: n: input, float, the order of the filter, if n is small, then the BHPF will be close to GHPF, and more smooth from lowfrequency to high freqency.if n is large, will close to IHPFreturn a [0, 1] value filter epsilon 1e-8M, N img.shape[1], img.shape[0]u np.arange(M)v np.arange(N)u, v np.meshgrid(u, v)D np.sqrt((u - center[1]//2)**2 (v - center[0]//2)**2)D0 radiuskernel (1 / (1 (D0 / (D epsilon))**(2*n)))return kernel# 理想、高斯、巴特沃斯高通传递函数
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cmcenter img_ori.shaperadius 100
IHPF idea_high_pass_filter(img_ori, img_ori.shape, radiusradius)
GHPF gauss_high_pass_filter(img_ori, img_ori.shape, radiusradius)
BHPF butterworth_high_pass_filter(img_ori, img_ori.shape, radiusradius, n2)filters [IHPF, GHPF, BHPF]
# 用来绘制3D图
M, N img_ori.shape[1], img_ori.shape[0]
u np.arange(M)
v np.arange(N)
u, v np.meshgrid(u, v)fig plt.figure(figsize(15, 15))for i in range(len(filters)):ax_1 fig.add_subplot(3, 3, i*3 1, projection3d)plot_3d(ax_1, u, v, eval(filters[i]))ax_2 fig.add_subplot(3, 3, i*3 2)ax_2.imshow(eval(filters[i]),gray), ax_2.set_title(filters[i]), ax_2.set_xticks([]), ax_2.set_yticks([])h_1 eval(filters[i])[img_ori.shape[0]//2:, img_ori.shape[1]//2]ax_3 fig.add_subplot(3, 3, i*3 3)ax_3.plot(h_1), ax_3.set_xticks([0, radius//2]), ax_3.set_yticks([0, 1]), ax_3.set_xlim([0, 320]), ax_3.set_ylim([0, 1.2])
plt.tight_layout()
plt.show()# 频率域理想、高通、巴特沃斯高通滤波器传递函数对应的空间核函数
# 这里简单用1 - 低通滤波器而生成
img_temp np.zeros([1000, 1000])radius 15
# IHPF idea_high_pass_filter(img_temp, img_temp.shape, radiusradius)
# GHPF gauss_high_pass_filter(img_temp, img_temp.shape, radiusradius)
# BHPF butterworth_high_pass_filter(img_temp, img_temp.shape, radiusradius, n2)
# filters [IHPF, GHPF, BHPF]ILPF idea_low_pass_filter(img_temp, img_temp.shape, radiusradius)
GLPF gauss_low_pass_filter(img_temp, img_temp.shape, radiusradius)
BLPF butterworth_low_pass_filter(img_temp, img_temp.shape, radiusradius, n2)
filters [ILPF, GLPF, BLPF]fig plt.figure(figsize(15, 10))
for i in range(len(filters)):# 这是显示空间域的核ax fig.add_subplot(2, 3, i1)spatial, spatial_s frequen2spatial(eval(filters[i]))ax.imshow(1 - spatial_s, gray), ax.set_xticks([]), ax.set_yticks([])# 这里显示是对应的空间域核水平扫描线的灰度分布ax fig.add_subplot(2, 3, i4)hx spatial[:, 500]hx centralized_2d(hx.reshape(1, -1)).flatten()ax.plot(1 - hx), ax.set_xticks([]), ax.set_yticks([])
plt.tight_layout()
plt.show()def hpf_test(img_ori, modeconstant, radius10, **params):param: img_ori: input imageparam: mode: input, str, is numpy pad parameter, default is constant. for more detail please refere to Numpy padparam: **params: can accept multiply parameters, params[filter] to choose which filter to use, if filterbutterworth, will need n1 following, value varyparam: radius: the cut-off frequency size, default is 10return image with negetive and positive value, you might need to normalize the image or clip the value to certain value to showM, N img_ori.shape[:2]# 填充fp pad_image(img_ori, modemode)# 中心化fp_cen centralized_2d(fp)# 正变换fft np.fft.fft2(fp_cen)# 滤波器if params[filter] gauss:H gauss_high_pass_filter(fp, centerfp.shape, radiusradius)elif params[filter] idea:H idea_high_pass_filter(fp, centerfp.shape, radiusradius)elif params[filter] butterworth:H butterworth_high_pass_filter(fp, centerfp.shape, radiusradius, n params[n])# 滤波HF fft * H# 反变换ifft np.fft.ifft2(HF)# 去中心化gp centralized_2d(ifft.real)# 还回来与原图像的大小g gp[:M, :N]# 把负值截断
# g np.clip(g, 0, g.max())dst g
# dst np.uint8(normalize(g) * 255)return dst# 理想、高斯、巴特沃斯高通滤波器字符测试模式
# 图像具有负值与正值这里把负值都置为0
img_ori cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0448(a)(characters_test_pattern).tif, -1)
radius [60, 160]
filters [idea, gauss, butterworth]
fig plt.figure(figsize(17, 10))
for i in range(len(radius)):for j in range(len(filters)):ax fig.add_subplot(2, 3, 3*ij1)if filters[j] butterworth:img hpf_test(img_ori, modereflect, radiusradius[i], filterfilters[j], n2)else:img hpf_test(img_ori, modereflect, radiusradius[i], filterfilters[j])img np.clip(img, 0, img.max())ax.imshow(img, gray, vmin0, vmax255), ax.set_xticks([]), ax.set_yticks([])ax.set_title(radius str(radius[i]) , filters[j])
plt.tight_layout()
plt.show()# 理想、高斯、巴特沃斯高通滤波器字符测试模式
# 图像具有负值与正值这里把图像重新标定显示
img_ori cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0448(a)(characters_test_pattern).tif, -1)
radius [160]
filters [idea, gauss, butterworth]
fig plt.figure(figsize(17, 10))
for i in range(len(radius)):for j in range(len(filters)):ax fig.add_subplot(2, 3, 3*ij1)img hpf_test(img_ori, filters[j], modereflect, radiusradius[i])img np.uint8(normalize(img) * 255)ax.imshow(img, gray), ax.set_xticks([]), ax.set_yticks([])ax.set_title(radius str(radius[i]) , filters[j])
plt.tight_layout()
plt.show()从上面的图形可以看到高通滤波器在边缘和边界的检测非常重要。
但理想高通滤波器出现了振铃效应。
# 频率域滤波过程
img_ori cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0448(a)(characters_test_pattern).tif, -1)
M, N img_ori.shape[:2]
fig plt.figure(figsize(15, 15))# 这里对原图像进行pad以得到更好的显示
img_ori_show np.ones([img_ori.shape[0]*2, img_ori.shape[1]*2]) * 255
img_ori_show[:img_ori.shape[0], img_ori.shape[1]:] img_ori
ax_1 fig.add_subplot(3, 3, 1)
imshow_img(img_ori_show, ax_1)fp pad_image(img_ori, modeconstant)
ax_2 fig.add_subplot(3, 3, 2)
imshow_img(fp, ax_2)fp_cen centralized_2d(fp)
fp_cen_show np.clip(fp_cen, 0, 255) # 负数的都截断为0
ax_3 fig.add_subplot(3, 3, 3)
ax_3.imshow(fp_cen_show, cmapgray), ax_3.set_xticks([]), ax_3.set_yticks([])fft np.fft.fft2(fp_cen)
spectrum spectrum_fft(fft)
spectrum np.log(1 spectrum)
ax_4 fig.add_subplot(3, 3, 4)
imshow_img(spectrum, ax_4)H gauss_high_pass_filter(fp, centerfp.shape, radius300)
ax_5 fig.add_subplot(3, 3, 5)
imshow_img(H, ax_5)HF fft * H
HF_spectrum np.log(1 spectrum_fft(HF))
ax_6 fig.add_subplot(3, 3, 6)
imshow_img(HF_spectrum, ax_6)ifft np.fft.ifft2(HF)
gp centralized_2d(ifft.real)
ax_7 fig.add_subplot(3, 3, 7)
imshow_img(gp, ax_7)# 最终结果有黑边
g gp[:M, :N]
g np.clip(g, 0, g.max())
ax_8 fig.add_subplot(3, 3, 8)
imshow_img(g, ax_8)plt.tight_layout()
plt.show()指纹增强
# 使用高能滤波和阈值处理增强图像
img_finger cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0457(a)(thumb_print).tif, 0) #直接读为灰度图像
img_back hpf_test(img_finger, modereflect, radius50, filterbutterworth, n4)
img_thres np.clip(img_back, 0, 1)plt.figure(figsize(24, 8))
plt.subplot(1, 3, 1),plt.imshow(img_finger,gray),plt.title(origial), plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 2),plt.imshow(img_back,gray),plt.title(After GHPF), plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 3),plt.imshow(img_thres,gray),plt.title(Threshold), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()# 使用高能滤波和阈值处理增强图像
import cv2
import numpy as np
import matplotlib.pyplot as pltimg_finger cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0457(a)(thumb_print).tif, 0) #直接读为灰度图像plt.figure(figsize(15, 12))
plt.subplot(221),plt.imshow(img_finger,gray),plt.title(origial)#--------------------------------
fft np.fft.fft2(img_finger)
fft_shift np.fft.fftshift(fft)
amp_img np.abs(np.log(1 np.abs(fft_shift)))
plt.subplot(222),plt.imshow(amp_img,gray),plt.title(FFT)#--------------------------------
BHPF butterworth_high_pass_filter(img_finger, img_finger.shape, radius30, n4)
plt.subplot(223),plt.imshow(BHPF,gray),plt.title(BHPF)#--------------------------------
f1shift fft_shift * BHPF
f2shift np.fft.ifftshift(f1shift) #对新的进行平移
img_back np.fft.ifft2(f2shift)#出来的是复数无法显示
img_new np.abs(img_back)#调整大小范围便于显示
img_new (img_new-np.amin(img_new))/(np.amax(img_new)-np.amin(img_new))
plt.subplot(224),plt.imshow(img_new,gray),plt.title(After BHPF)plt.tight_layout()
plt.show()#调整大小范围便于显示
plt.figure(figsize(15, 10))plt.subplot(121),plt.imshow(img_finger,gray),plt.title(Original)img_back_thred np.clip(img_back.real, 0, 1)plt.subplot(122),plt.imshow(np.abs(img_back_thred),gray),plt.title(Thredhold after BHPF)plt.tight_layout()
plt.show()频域中的拉普拉斯
拉普拉斯可用如下滤波器传递函数在频率域中实现
H(u,v)−4π2(u2v2)(4.123)H(u,v) -4\pi^2 (u^2 v^2) \tag{4.123}H(u,v)−4π2(u2v2)(4.123)
或者关于频率矩形的中心
H(u,v)−4π2[(u−M/2)2(v−N/2)2]−4π2D2(u,v)(4.124)H(u,v) -4\pi^2 [(u - M/2)^2 (v-N/2)^2] -4\pi^2 D^2(u,v) \tag{4.124}H(u,v)−4π2[(u−M/2)2(v−N/2)2]−4π2D2(u,v)(4.124)
D(u,v)[(u−M/2)2(v−N/2)2]1/2D(u,v) [(u - M/2)^2 (v-N/2)^2]^{1/2}D(u,v)[(u−M/2)2(v−N/2)2]1/2
∇2f(x,y)J−1[H(u,v)F(u,v)](4.125)\nabla^2 f(x, y) \mathfrak{J}^{-1}[H(u, v)F(u, v)] \tag{4.125}∇2f(x,y)J−1[H(u,v)F(u,v)](4.125) g(x,y)f(xy)c∇2f(x,y)(4.126)g(x, y) f(x y) c\nabla^2 f(x, y) \tag{4.126}g(x,y)f(xy)c∇2f(x,y)(4.126)
其中c−1c-1c−1因为H(u,v)H(u, v)H(u,v)是负的。 用式(4.125)计算∇2f(x,y)\nabla^2 f(x, y)∇2f(x,y)这个因子的幅度要比fff的最大值要大几个数量级所以要将fff与∇2f(x,y)\nabla^2 f(x, y)∇2f(x,y)的差限定在可比的范围内。
方法 计算f(x,y)f(x, y)f(x,y)的DFT之前将f(x,y)f(x, y)f(x,y)的值归一化到[0, 1]并将∇2f(x,y)\nabla^2 f(x, y)∇2f(x,y)除以它的最大值进行限定在近似区间[-1, 1]。
g(x,y)J−1{F(u,v)−H(u,v)F(u,v)}J−1{[1−H(u,v)]F(u,v)}J−1{[14π2D2]F(u,v)}(4.127)\begin{aligned} g(x, y) \mathfrak{J}^{-1}\big\{ F(u, v) - H(u, v)F(u, v) \big\} \\ \mathfrak{J}^{-1}\big\{ [1 - H(u, v)]F(u, v) \big\} \mathfrak{J}^{-1}\big\{ [1 4\pi^2 D^2]F(u, v) \big\} \end{aligned}\tag{4.127}g(x,y)J−1{F(u,v)−H(u,v)F(u,v)}J−1{[1−H(u,v)]F(u,v)}J−1{[14π2D2]F(u,v)}(4.127)
所以一般用式(4.125)公式先计算滤波后的结果再进行处理。
def lapalacian_filter(source):这效果跟原来的不一样有特改进M, N source.shape[1], source.shape[0]u np.arange(M)v np.arange(N)u, v np.meshgrid(u, v)D np.sqrt((u - M//2)**2 (v - N//2)**2)kernel -4 * (np.pi**2 ) * (D**2)kernel_normal kernel # (kernel - kernel.min()) / (kernel.max() - kernel.min())return kernel_normal# 频域中的拉普拉斯滤波过程
img_ori cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0458(a)(blurry_moon).tif, -1)
img_norm img_ori / 255
M, N img_ori.shape[:2]
fp pad_image(img_norm, modeconstant)
fp_cen centralized_2d(fp)
fft np.fft.fft2(fp_cen)
H lapalacian_filter(fp)
HF fft * H
ifft np.fft.ifft2(HF)
gp centralized_2d(ifft.real)
g gp[:M, :N]# 标定因子
g1 g / g.max()
# 归一化
img img_ori / 255.
# 根据公式进行增强
img_res img - g1
img_res np.clip(img_res, 0, 1)plt.figure(figsize(13, 8))
plt.subplot(1,2, 1), plt.imshow(img,gray),plt.title(ORI), plt.xticks([]), plt.yticks([])
plt.subplot(1,2, 2), plt.imshow(img_res,gray),plt.title(Transform), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()钝化掩蔽、高提升滤波和高频强调滤波
gmask(x,y)f(x,y)−fLP(x,y)(4.128)g_{mask}(x, y) f(x, y) - f_{LP}(x,y) \tag{4.128}gmask(x,y)f(x,y)−fLP(x,y)(4.128)
fLP(x,y)J−1[HLP(u,v)F(u,v)](4.129)f_{LP}(x,y) \mathfrak{J}^{-1} \big[ H_{LP}(u, v) F(u, v) \big] \tag{4.129}fLP(x,y)J−1[HLP(u,v)F(u,v)](4.129)
g(x,y)f(x,y)kgmask(x,y)(3.56)g(x,y) f(x,y) k g_{mask}(x, y) \tag{3.56}g(x,y)f(x,y)kgmask(x,y)(3.56)
权值k≥0k \ge 0k≥0k1k 1k1时它是钝化掩蔽k1k 1k1时这个过程称为高提升滤波选择k≤1k \leq 1k≤1可以减少钝化模板的贡献。
g(x,y)J−1{[1k[1−HLP(u,v)]]F(u,v)}(4.131)g(x,y) \mathfrak{J}^{-1} \Big\{\big[1 k[1 - H_{LP}(u, v)]\big] F(u, v) \Big\} \tag{4.131}g(x,y)J−1{[1k[1−HLP(u,v)]]F(u,v)}(4.131) g(x,y)J−1{[1kHHP(u,v)]F(u,v)}(4.132)g(x,y) \mathfrak{J}^{-1} \Big\{\big[1 kH_{HP}(u, v)\big] F(u, v) \Big\} \tag{4.132}g(x,y)J−1{[1kHHP(u,v)]F(u,v)}(4.132)
[1kHHP(u,v)]\big[1 kH_{HP}(u, v)\big][1kHHP(u,v)]称为高频强调滤波器传递函数。
高频强调滤波 g(x,y)J−1{[k1k2HHP(u,v)]F(u,v)}(4.133)g(x,y) \mathfrak{J}^{-1} \Big\{\big[k_1 k_2H_{HP}(u, v)\big] F(u, v) \Big\} \tag{4.133}g(x,y)J−1{[k1k2HHP(u,v)]F(u,v)}(4.133)
k1≥0k_1 \ge 0k1≥0偏移传递函数的值以便使直流项不为零而k2≥0k_2 \ge 0k2≥0控制高频的贡献
def frequency_filter(fshift, h):Frequency domain filtering using FFTparam: fshift: input, FFT shift to centerparam: h: input, filter, value range [0, 1]return an unnormalized reconstruct image, value may have negative and positive value# 滤波res fshift * h# 反变换ifft np.fft.ifft2(res)# 取实数部分gp centralized_2d(ifft.real)return gpdef cdf_interp(img):global histogram equalization used numpyhist, bins np.histogram(img, bins256)cdf hist.cumsum()cdf 255 * cdf / cdf[-1] # 这个方法好像合理点灰度值会比较小img_dst np.interp(img.flatten(), bins[:-1], cdf)img_dst np.uint8(img_dst.reshape(img.shape))return img_dst# 高频强调滤波 直方图均衡化
# 从图像来看确实优于单独使用其中一种方法得到的结果但直接对原图进行直方图均衡化也能得到不错的结果。
img_ori cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0459(a)(orig_chest_xray).tif, -1)
M, N img_ori.shape[:2]# 填充
fp pad_image(img_ori, modeconstant)
fp_cen centralized_2d(fp)
fft np.fft.fft2(fp_cen)# 滤波器
GHPF gauss_high_pass_filter(fp, fp.shape, radius70)
# 滤波后的图像
img_new frequency_filter(fft, GHPF)
img_new img_new[:M, :N]#高频增强滤波
k_1 0.5
k_2 0.75
res fft * (k_1 k_2 * GHPF)# 反变换
ifft np.fft.ifft2(res)# 取实数部分
gp centralized_2d(ifft.real)
g gp[:M, :N]
g np.clip(g, 0, g.max())temp np.uint8(normalize(g) * 255)
img_res cdf_interp(temp)hist_ori cdf_interp(img_ori)plt.figure(figsize(13, 14))
plt.subplot(3, 2, 1), plt.imshow(img_ori,gray),plt.title(ORI), plt.xticks([]), plt.yticks([])
plt.subplot(3, 2, 2), plt.imshow(img_new,gray),plt.title(GHPF), plt.xticks([]), plt.yticks([])
plt.subplot(3, 2, 3), plt.imshow(g,gray),plt.title(High frequency emphasis filtering), plt.xticks([]), plt.yticks([])
plt.subplot(3, 2, 4), plt.imshow(img_res,gray),plt.title(Histogram equalization), plt.xticks([]), plt.yticks([])
plt.subplot(3, 2, 6), plt.imshow(hist_ori,gray),plt.title(Histogram equalization ORI), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()同态滤波
照射-反射模型可以赶通告个频率域程序这个程序通过灰度范围压缩和对比度增强来同时改善图像的外观。图像f(x,y)f(x, y)f(x,y)可以表示为其照射分量i(x,y)i(x, y)i(x,y)和反射分量r(x,y)r(x, y)r(x,y) f(x,y)i(x,y)r(x,y)(4.134)f(x, y) i(x, y) r(x, y) \tag{4.134}f(x,y)i(x,y)r(x,y)(4.134)
因为乘积的傅里叶变换不是变换的乘积 J[f(x,y)]≠J[i(x,y)]J[r(x,y)](4.135)\mathfrak{J}\big[f(x, y)\big] \neq \mathfrak{J}\big[i(x, y)\big] \mathfrak{J}\big[r(x, y)\big] \tag{4.135}J[f(x,y)]J[i(x,y)]J[r(x,y)](4.135)
但我们可以自然对数变换得到 J[z(x,y)]J[lnf(x,y)]J[lni(x,y)]J[lnr(x,y)](4.137)\mathfrak{J}\big[z(x, y)\big] \mathfrak{J}\big[\text{ln}\;f(x, y)\big] \mathfrak{J}\big[\text{ln}\; i(x, y)\big] \mathfrak{J}\big[\text{ln}\; r(x, y)\big] \tag{4.137}J[z(x,y)]J[lnf(x,y)]J[lni(x,y)]J[lnr(x,y)](4.137) 或者是 Z(u,v)Fi(u,v)Fr(u,v)(4.138)Z(u, v) F_i(u, v) F_r(u, v) \tag{4.138}Z(u,v)Fi(u,v)Fr(u,v)(4.138)
在频率域可以用滤波器传递函数H(u,v)H(u, v)H(u,v)对Z(u,v)Z(u, v)Z(u,v)滤波则有 S(u,v)H(u,v)Z(u,v)H(u,v)Fi(u,v)H(u,v)Fr(u,v)(4.139)S(u, v) H(u, v)Z(u, v) H(u, v)F_i(u, v) H(u, v)F_r(u, v) \tag{4.139}S(u,v)H(u,v)Z(u,v)H(u,v)Fi(u,v)H(u,v)Fr(u,v)(4.139)
空间域中滤波后的图像是 s(x,y)J−1[S(u,v)]J−1[H(u,v)Fi(u,v)]J−1[H(u,v)Fr(u,v)](4.140)s(x, y) \mathfrak{J}^{-1}[S(u, v)] \mathfrak{J}^{-1}[H(u, v)F_i(u, v)] \mathfrak{J}^{-1}[H(u, v)F_r(u, v)] \tag{4.140}s(x,y)J−1[S(u,v)]J−1[H(u,v)Fi(u,v)]J−1[H(u,v)Fr(u,v)](4.140)
所以有 i′(x,y)J−1[H(u,v)Fi(u,v)](4.141)i(x, y) \mathfrak{J}^{-1}[H(u, v)F_i(u, v)] \tag{4.141}i′(x,y)J−1[H(u,v)Fi(u,v)](4.141) r′(x,y)J−1[H(u,v)Fr(u,v)](4.142)r(x, y) \mathfrak{J}^{-1}[H(u, v)F_r(u, v)] \tag{4.142}r′(x,y)J−1[H(u,v)Fr(u,v)](4.142) ⇒\Rightarrow⇒ s(x,y)i′(x,y)r′(x,y)(4.143)s(x, y) i(x, y) r(x, y) \tag{4.143}s(x,y)i′(x,y)r′(x,y)(4.143)
因为z(x,y)z(x, y)z(x,y)是通过输入图像的自然对数形成的所以可以通过取滤波后的结果的指数这一反处理来形成输出图像 g(x,y)es(x,y)ei′(x,y)er′(x,y)(4.144)g(x, y) e^{s(x, y)} e^{i(x, y)} e^{r(x, y)} \tag{4.144}g(x,y)es(x,y)ei′(x,y)er′(x,y)(4.144)
这种方法是同态系统的一种特殊情况。可用同态滤波器传递函数H(u,v)H(u, v)H(u,v)分别对照射分量和反射分量进行操作。
图像的照射分量通常由慢空间变化来表征而反射分量往往倾向于突变特别是在不同目标的连接处。根据这此特征我们可以将图像取对数后的傅里叶变换的低频与照射联系起来而将高频与反射联系起来。尽管这些联系只是粗略的挖但它们可以用于图像滤波。
同态滤波的步骤
f(x,y)⇒ln⇒DFT⇒H(u,v)⇒(DFT)−1⇒exp⇒g(x,y)f(x, y) \Rightarrow \text{ln} \Rightarrow \text{DFT} \Rightarrow H(u, v) \Rightarrow (\text{DFT})^{-1} \Rightarrow \text{exp} \Rightarrow g(x, y)f(x,y)⇒ln⇒DFT⇒H(u,v)⇒(DFT)−1⇒exp⇒g(x,y)
使用同态滤波器可更好的控制照射分量和反射分量。这种控制要求规定滤波器传递函数H(u,v)H(u, v)H(u,v)以便以不同的控制方法来影响傅里叶变换的低频和高频分量。
采用形式上稍有变化的GHPF函数可得到同态函数 H(u,v)(γH−γL)[1−e−cD2(u,v)/D02]γL(4.147)H(u, v) (\gamma_H - \gamma_L)\Big[1 - e^{-cD^2(u,v)/D_0^2}\Big] \gamma_L \tag{4.147}H(u,v)(γH−γL)[1−e−cD2(u,v)/D02]γL(4.147) 若选择的γL\gamma_LγL和γH\gamma_HγH满足γL1且γH≥1\gamma_L 1且\gamma_H \ge 1γL1且γH≥1则滤波器函数将衰减低频照射的贡献而放大高频反射的贡献。最终结果是同时进行动态范围的压缩和对比度的增强。
def gauss_homomorphic_filter(source, center, rl0.5, rh1, c1, radius5):create gaussian high pass filter param: source: input, source imageparam: center: input, the center of the filter, where is the lowest value, (0, 0) is top left corner, source.shape[:2] is center of the source imageparam: radius: input, the radius of the lowest value, greater value, bigger blocker out range, if the radius is 0, then allvalue is 1return a [rl, rh] value filterM, N source.shape[1], source.shape[0]u np.arange(M)v np.arange(N)u, v np.meshgrid(u, v)D np.sqrt((u - center[1]//2)**2 (v - center[0]//2)**2)D0 radiuskernel 1 - np.exp(- c * (D**2)/(D0**2)) kernel (rh - rl) * kernel rlreturn kernel# 高斯同态滤波器与传递函数的径向剖面图
img_temp np.zeros([1000, 1000])
GHF gauss_homomorphic_filter(img_temp, img_temp.shape, rl0.6, rh1, c0.4, radius50)hx GHF[500:, 500].flatten()fig plt.figure(figsize(15, 5))
ax_1 fig.add_subplot(1, 2, 1)
ax_1.imshow(GHF, gray), ax_1.set_xticks([]), ax_1.set_yticks([])ax_2 fig.add_subplot(1, 2, 2)
ax_2.plot(hx), #ax_3.set_xticks([]), ax_3.set_yticks([])plt.tight_layout()
plt.show()# 高斯同态滤波器
img_ori cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0462(a)(PET_image).tif, -1)
M, N img_ori.shape[:2]# 填充
fp pad_image(img_ori, modeconstant)
fp_cen centralized_2d(fp)
fft np.fft.fft2(fp_cen)# 滤波器
GHF gauss_homomorphic_filter(fp, fp.shape, rl0.5, rh3, c5, radius20)# 滤波后的图像
img_new frequency_filter(fft, GHF)
img_new img_new[:M, :N]
# print(img_new.min(), img_new.max())
img_res np.uint8(normalize((img_new)) * 255)plt.figure(figsize(13, 14))
plt.subplot(1, 2, 1), plt.imshow(img_ori,gray),plt.title(ORI), plt.xticks([]), plt.yticks([])
plt.subplot(1, 2, 2), plt.imshow(img_res,gray),plt.title(G-Homomorphic), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()# 高斯同态滤波器根据推导过程应用但得到不好的图像然后上面的直接应用的话就算改变不同的参数变化也不是很大
img_ori cv2.imread(DIP_Figures/DIP3E_Original_Images_CH04/Fig0462(a)(PET_image).tif, -1)
M, N img_ori.shape[:2]
img_ln np.log(img_ori1e-8)# 填充
fp pad_image(img_ln, modeconstant)
fp_cen centralized_2d(fp)
fft np.fft.fft2(fp_cen)# 滤波器
GHF gauss_homomorphic_filter(fp, fp.shape, rl0.8, rh1, c3, radius20)# 滤波后的图像
img_new frequency_filter(fft, GHF)
img_new img_new[:M, :N]
img_new np.exp(img_new)img_new np.clip(img_new, 0, img_new.max())
img_res np.uint8(img_new / img_new.max() * 255)plt.figure(figsize(13, 14))
plt.subplot(1, 2, 1), plt.imshow(img_ori,gray),plt.title(ORI), plt.xticks([]), plt.yticks([])
plt.subplot(1, 2, 2), plt.imshow(img_new,gray),plt.title(G-Homomorphic), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()
