第4章 Python 数字图像处理(DIP) - 频率域滤波11 - 使用高通滤波器锐化图像

使用高通滤波器锐化图像

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

$$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 image
    param: 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 image
    param: radius: input, the radius of the lowest value, greater value, bigger blocker out range, if the radius is 0, then all
                   value is 0
    return 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 = radius
    kernel = D.copy()
    
    kernel[D > D0] = 1
    kernel[D <= D0] = 0
    
    return kernel

def gauss_high_pass_filter(source, center, radius=5):
    """
    create gaussian high pass filter 
    param: source: input, source image
    param: 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 image
    param: radius: input, the radius of the lowest value, greater value, bigger blocker out range, if the radius is 0, then all
                   value is 1
    return 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 = radius
    
    if 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 image
    param: 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 image
    param: radius: input, the radius of the lowest value, greater value, bigger blocker out range, if the radius is 0, then all
                   value is 0
    param: n: input, float, the order of the filter, if n is small, then the BHPF will be close to GHPF, and more smooth from low
              frequency to high freqency.if n is large, will close to IHPF
    return a [0, 1] value filter
    """  
    epsilon = 1e-8
    M, 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 = radius
    kernel = (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 cm

center = img_ori.shape

radius = 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 image
    param: mode:     input, str, is numpy pad parameter, default is 'constant'. for more detail please refere to Numpy pad
    param: **params: can accept multiply parameters, params['filter'] to choose which filter to use, if filter='butterworth', 
                     will need n=1 following, value vary
    param: radius:   the cut-off frequency size, default is 10
    return 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 plt
 
img_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 FFT
    param: fshift: input, FFT shift to center
    param: 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 image
    param: 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 image
    param: radius: input, the radius of the lowest value, greater value, bigger blocker out range, if the radius is 0, then all
                   value is 1
    return 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 = radius
    
    kernel = 1 - np.exp(- c * (D**2)/(D0**2))   
    kernel = (rh - rl) * kernel + rl
    return 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()