使用高通滤波器锐化图像

由低通滤波器得到理想、高斯和巴特沃斯高通滤波器

$$H_{HP}(u,v)=1-H_{LP}(u,v)\text{\hspace{1em}(4.118)}$$

• 理想高通
  $$H(u,v)=\{\begin{array}{ll}0,& D(u,v)\le D_0\\ 1,& D(u,v)>D_0\end{array}\text{\hspace{1em}(4.119)}$$

• 高斯高通
  $$H(u,v)=1-e^{-D^2(u,v)/2D_0^2}\text{\hspace{1em}(4.120)}$$

• 巴特沃斯高通
  $$H(u,v)=\frac{1}{1+[D_0/D(u,v){]}^{2n}}\text{\hspace{1em}(4.121)}$$

def idea_high_pass_filter(source, center, radius=5):"""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 filter"""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 - center[1]//2)**2 + (v - center[0]//2)**2)D0 = radiuskernel = D.copy()kernel[D > D0] = 1kernel[D <= D0] = 0return kernel

def gauss_high_pass_filter(source, center, radius=5):"""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 filter"""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 - center[1]//2)**2 + (v - center[0]//2)**2)D0 = radiusif D0 == 0:kernel = np.ones(source.shape[:2], dtype=np.float)else:kernel = 1 - np.exp(- (D**2)/(D0**2))   return kernel

def butterworth_high_pass_filter(img, center, radius=5, n=1):"""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, radius=radius)
GHPF = gauss_high_pass_filter(img_ori, img_ori.shape, radius=radius)
BHPF = butterworth_high_pass_filter(img_ori, img_ori.shape, radius=radius, n=2)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, projection='3d')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, radius=radius)
# GHPF = gauss_high_pass_filter(img_temp, img_temp.shape, radius=radius)
# BHPF = butterworth_high_pass_filter(img_temp, img_temp.shape, radius=radius, n=2)
# filters = ['IHPF', 'GHPF', 'BHPF']ILPF = idea_low_pass_filter(img_temp, img_temp.shape, radius=radius)
GLPF = gauss_low_pass_filter(img_temp, img_temp.shape, radius=radius)
BLPF = butterworth_low_pass_filter(img_temp, img_temp.shape, radius=radius, n=2)
filters = ['ILPF', 'GLPF', 'BLPF']fig = plt.figure(figsize=(15, 10))
for i in range(len(filters)):# 这是显示空间域的核ax = fig.add_subplot(2, 3, i+1)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, i+4)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, mode='constant', radius=10, **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 filter='butterworth', will need n=1 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 show"""M, N = img_ori.shape[:2]# 填充fp = pad_image(img_ori, mode=mode)# 中心化fp_cen = centralized_2d(fp)# 正变换fft = np.fft.fft2(fp_cen)# 滤波器if params['filter'] == "gauss":H = gauss_high_pass_filter(fp, center=fp.shape, radius=radius)elif params['filter'] == "idea":H = idea_high_pass_filter(fp, center=fp.shape, radius=radius)elif params['filter'] == "butterworth":H = butterworth_high_pass_filter(fp, center=fp.shape, radius=radius, 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*i+j+1)if filters[j] == 'butterworth':img = hpf_test(img_ori, mode='reflect', radius=radius[i], filter=filters[j], n=2)else:img = hpf_test(img_ori, mode='reflect', radius=radius[i], filter=filters[j])img = np.clip(img, 0, img.max())ax.imshow(img, 'gray', vmin=0, vmax=255), 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*i+j+1)img = hpf_test(img_ori, filters[j], mode='reflect', radius=radius[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, mode='constant')
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, cmap='gray'), 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, center=fp.shape, radius=300)
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, mode='reflect', radius=50, filter='butterworth', n=4)
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, radius=30, n=4)
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{\pi }^2(u^2+v^2)\text{\hspace{1em}(4.123)}$$

或者关于频率矩形的中心

$$H(u,v)=-4{\pi }^2[(u-M/2{)}^2+(v-N/2{)}^2]=-4{\pi }^2{D}^2(u,v)\text{\hspace{1em}(4.124)}$$

$$D(u,v)=[(u-M/2{)}^2+(v-N/2{)}^2{]}^{1/2}$$

$${\nabla }^{2}f(x,y)={\mathfrak{J}}^{-1}[H(u,v)F(u,v)]\text{\hspace{1em}(4.125)}$$
$$g(x,y)=f(xy)+c{\nabla }^{2}f(x,y)\text{\hspace{1em}(4.126)}$$

其中，$c=-1$，因为$H(u,v)$是负的。
用式(4.125)计算${\nabla }^{2}f(x,y)$，这个因子的幅度要比$f$的最大值要大几个数量级，所以要将$f$与${\nabla }^{2}f(x,y)$的差限定在可比的范围内。

• 方法：
  • 计算$f(x,y)$的DFT之前将$f(x,y)$的值归一化到[0, 1]
  • 并将${\nabla }^{2}f(x,y)$除以它的最大值，进行限定在近似区间[-1, 1]。

$$\begin{array}{rl}g(x,y)& ={\mathfrak{J}}^{-1}\{F(u,v)-H(u,v)F(u,v)\}\\ & ={\mathfrak{J}}^{-1}\{[1-H(u,v)]F(u,v)\}={\mathfrak{J}}^{-1}\{[1+4{\pi }^2{D}^2]F(u,v)\}\end{array}\text{\hspace{1em}(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, mode='constant')
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()

钝化掩蔽、高提升滤波和高频强调滤波

$$g_{mask}(x,y)=f(x,y)-f_{LP}(x,y)\text{\hspace{1em}(4.128)}$$

$$f_{LP}(x,y)={\mathfrak{J}}^{-1}[H_{LP}(u,v)F(u,v)]\text{\hspace{1em}(4.129)}$$

$$g(x,y)=f(x,y)+k{g}_{mask}(x,y)\text{\hspace{1em}(3.56)}$$

权值$k\ge 0$，$k=1$时，它是钝化掩蔽，$k>1$时，这个过程称为高提升滤波，选择$k\le 1$可以减少钝化模板的贡献。

$$g(x,y)={\mathfrak{J}}^{-1}\left\{\right[1+k[1-H_{LP}(u,v)]]F(u,v)\}\text{\hspace{1em}(4.131)}$$
$$g(x,y)={\mathfrak{J}}^{-1}\left\{\right[1+k{H}_{HP}(u,v)]F(u,v)\}\text{\hspace{1em}(4.132)}$$

$[1+k{H}_{HP}(u,v)]$称为高频强调滤波器传递函数。

• 高频强调滤波
  $$g(x,y)={\mathfrak{J}}^{-1}\left\{\right[k_1+k_2{H}_{HP}(u,v)]F(u,v)\}\text{\hspace{1em}(4.133)}$$

$k_1\ge 0$偏移传递函数的值，以便使直流项不为零，而$k_2\ge 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 gp

def cdf_interp(img):"""global histogram equalization used numpy"""hist, bins = np.histogram(img, bins=256)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, mode='constant')
fp_cen = centralized_2d(fp)
fft = np.fft.fft2(fp_cen)# 滤波器
GHPF = gauss_high_pass_filter(fp, fp.shape, radius=70)
# 滤波后的图像
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)$可以表示为其照射分量$i(x,y)$和反射分量$r(x,y)$：
$$f(x,y)=i(x,y)r(x,y)\text{\hspace{1em}(4.134)}$$

因为乘积的傅里叶变换不是变换的乘积：
$$\mathfrak{J}[f(x,y)]\ne \mathfrak{J}[i(x,y)]\mathfrak{J}[r(x,y)]\text{\hspace{1em}(4.135)}$$

但我们可以自然对数变换得到：
$$\mathfrak{J}[z(x,y)]=\mathfrak{J}[\text{ln}f(x,y)]=\mathfrak{J}[\text{ln}i(x,y)]\mathfrak{J}[\text{ln}r(x,y)]\text{\hspace{1em}(4.137)}$$
或者是
$$Z(u,v)=F_i(u,v)+F_r(u,v)\text{\hspace{1em}(4.138)}$$

在频率域，可以用滤波器传递函数$H(u,v)$对$Z(u,v)$滤波，则有
$$S(u,v)=H(u,v)Z(u,v)=H(u,v)F_i(u,v)+H(u,v)F_r(u,v)\text{\hspace{1em}(4.139)}$$

空间域中滤波后的图像是：
$$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)]\text{\hspace{1em}(4.140)}$$

所以有
$$i^{\prime }(x,y)={\mathfrak{J}}^{-1}[H(u,v)F_i(u,v)]\text{\hspace{1em}(4.141)}$$
$$r^{\prime }(x,y)={\mathfrak{J}}^{-1}[H(u,v)F_r(u,v)]\text{\hspace{1em}(4.142)}$$
$\Rightarrow$
$$s(x,y)=i^{\prime }(x,y)+r^{\prime }(x,y)\text{\hspace{1em}(4.143)}$$

因为$z(x,y)$是通过输入图像的自然对数形成的，所以可以通过取滤波后的结果的指数这一反处理来形成输出图像：
$$g(x,y)=e^{s(x,y)}=e^{i^{\prime }(x,y)}{e}^{r^{\prime }(x,y)}\text{\hspace{1em}(4.144)}$$

这种方法是同态系统的一种特殊情况。可用同态滤波器传递函数$H(u,v)$分别对照射分量和反射分量进行操作。

图像的照射分量通常由慢空间变化来表征，而反射分量往往倾向于突变，特别是在不同目标的连接处。根据这此特征，我们可以将图像取对数后的傅里叶变换的低频与照射联系起来，而将高频与反射联系起来。尽管这些联系只是粗略的挖，但它们可以用于图像滤波。

• 同态滤波的步骤

$$f(x,y)\Rightarrow \text{ln}\Rightarrow \text{DFT}\Rightarrow H(u,v)\Rightarrow (\text{DFT}{)}^{-1}\Rightarrow \text{exp}\Rightarrow g(x,y)$$

使用同态滤波器可更好的控制照射分量和反射分量。这种控制要求规定滤波器传递函数$H(u,v)$，以便以不同的控制方法来影响傅里叶变换的低频和高频分量。

• 采用形式上稍有变化的GHPF函数可得到同态函数：
  $$H(u,v)=({\gamma }_H-{\gamma }_L)[1-e^{-c{D}^2(u,v)/D_0^2}]+{\gamma }_L\text{\hspace{1em}(4.147)}$$
  若选择的${\gamma }_L$和${\gamma }_H$满足${\gamma }_L<1且{\gamma }_H\ge 1$，则滤波器函数将衰减低频（照射）的贡献，而放大高频（反射）的贡献。最终结果是同时进行动态范围的压缩和对比度的增强。

def gauss_homomorphic_filter(source, center, rl=0.5, rh=1, c=1, radius=5):"""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 filter"""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 - 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, rl=0.6, rh=1, c=0.4, radius=50)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, mode='constant')
fp_cen = centralized_2d(fp)
fft = np.fft.fft2(fp_cen)# 滤波器
GHF = gauss_homomorphic_filter(fp, fp.shape, rl=0.5, rh=3, c=5, radius=20)# 滤波后的图像
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_ori+1e-8)# 填充
fp = pad_image(img_ln, mode='constant')
fp_cen = centralized_2d(fp)
fft = np.fft.fft2(fp_cen)# 滤波器
GHF = gauss_homomorphic_filter(fp, fp.shape, rl=0.8, rh=1, c=3, radius=20)# 滤波后的图像
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()