第4章 Python 数字图像处理(DIP) - 频率域滤波9 - 频率域滤波基础、频率域的滤波过程、低通、高通

频率域滤波基础

频率域的其他特性

• 频率域中的滤波过程如下：
  • 首先修改傅里叶变换以在到特定目的
  • 然后计算IDFT，返回到空间域

# 频率域中的其他特性
img = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0429(a)(blown_ic).tif', -1)

# FFT
img_fft = np.fft.fft2(img.astype(np.float32))

# 非中心化的频谱
spectrum = spectrum_fft(img_fft)
spectrum = np.uint8(normalize(spectrum) * 255)

# 中心化
fshift = np.fft.fftshift(img_fft)            # 将变换的频率图像四角移动到中心

# 中心化后的频谱
spectrum_fshift = spectrum_fft(fshift)
spectrum_fshift_n = np.uint8(normalize(spectrum_fshift) * 255)

# 对频谱做对数变换
spectrum_log = np.log(1 + spectrum_fshift)

plt.figure(figsize=(15, 8))
plt.subplot(1, 2, 1), plt.imshow(img, cmap='gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(1, 2, 2), plt.imshow(spectrum_log, cmap='gray'), plt.title('FFT Shift log'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

频率域滤波基础知识

已知大小为$P\times Q$像素的一幅数字图像$f(x,y)$，则我们感兴趣取的基本滤波公式是：
g ( x , y ) = R e a l { J − 1 [ H ( u , v ) F ( u , v ) ] } (4.104) g(x, y) = Real\{ \mathfrak{J}^{-1}[H(u, v) F(u, v)] \} \tag{4.104} g(x,y)=Real{J−1[H(u,v)F(u,v)]}(4.104)

${\mathfrak{J}}^{-1}$是IDFT，$H(u,v)$是滤波器传递函数（经常称为滤波器或滤波器函数），$F(u,v)$是输入图像$f(x,y)$的DFT，$g(x,y)$是滤波后的图像。

中以化：用$(-1{)}^{x+y}$乘以输入图像来完成的

• 低通滤波器

  • 会衰减高频而通过低频的函数，会模糊图像

• 高通滤波器

  • 将使图像的细节更清晰，但会降低图像的对比度

def centralized_2d_loop(arr):
    """
    centralized 2d array $f(x, y) (-1)^{x + y}$
    """    
    #==================中心化=============================
    height, width = arr.shape[:2]
    dst = np.zeros([height, width])
    
    for h in range(height):
        for w in range(width):
            dst[h, w] = arr[h, w] * ((-1) ** (h + w))
    return dst

def centralized_2d_index(arr):
    """
    centralized 2d array $f(x, y) (-1)^{x + y}$, about 35 times faster than loop
    """
    
    def get_1_minus1(img):
        """
        get 1 when image index is even, -1 while index is odd, same shape as input image, need this array to multiply with input image
        to get centralized image for DFT
        Parameter: img: input, here we only need img shape to create the array
        return such a [[1, -1, 1], [-1, 1, -1]] array, example is 3x3
        """
        M, N = img.shape
        m, n = img.shape
        # 这里判断一个m, n是否是偶数，如果是偶数的话，结果不正确，需要变换为奇数
        if m % 2 == 0 :
            m += 1
        if n % 2 ==0 :
            n += 1
        temp = np.arange(m * n).reshape(m, n)
        dst = np.piecewise(temp, [temp % 2 == 0, temp % 2 == 1], [1, -1])
        dst = dst[:M, :N]
        return dst

    #==================中心化=============================
    mask = get_1_minus1(arr)
    dst = arr * mask

    return dst

def centralized_2d(arr):
    """
    centralized 2d array $f(x, y) (-1)^{x + y}$, about 4.5 times faster than index, 160 times faster than loop,
    """
    
    def get_1_minus1(img):
        """
        get 1 when image index is even, -1 while index is odd, same shape as input image, need this array to multiply with input image
        to get centralized image for DFT
        Parameter: img: input, here we only need img shape to create the array
        return such a [[1, -1, 1], [-1, 1, -1]] array, example is 3x3
        """
        dst = np.ones(img.shape)
        dst[1::2, ::2] = -1
        dst[::2, 1::2] = -1
        return dst

    #==================中心化=============================
    mask = get_1_minus1(arr)
    dst = arr * mask

    return dst

a = np.arange(25).reshape(5, 5)
a

array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19],
       [20, 21, 22, 23, 24]])

centralized_2d(a)

array([[  0.,  -1.,   2.,  -3.,   4.],
       [ -5.,   6.,  -7.,   8.,  -9.],
       [ 10., -11.,  12., -13.,  14.],
       [-15.,  16., -17.,  18., -19.],
       [ 20., -21.,  22., -23.,  24.]])

def idea_high_pass_filter(source, radius=5):
    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)
    D0 = radius
    
    kernel = D.copy()
    
    kernel[D > D0] = 1
    kernel[D <= D0] = 0
    
    kernel_normal = (kernel - kernel.min()) / (kernel.max() - kernel.min())
    return kernel_normal

# 频率域滤波基础知识
img = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0429(a)(blown_ic).tif', -1)

# FFT
img_fft = np.fft.fft2(img.astype(np.float32))

# 中心化
fshift = np.fft.fftshift(img_fft)            # 将变换的频率图像四角移动到中心

# 滤波器
h = idea_high_pass_filter(img, radius=1)

# 滤波
res = fshift * h

# 先反中以化，再反变换
f1shift = np.fft.ifftshift(res)
ifft = np.fft.ifft2(f1shift)

# 取实数部分
recon = np.real(ifft)
# 小于0都被置为0
# recon = np.where(recon > 0, recon, 0)
recon = np.clip(recon, 0, recon.max())
recon = np.uint8(normalize(recon) * 255)

plt.figure(figsize=(15, 8))
plt.subplot(1, 2, 1), plt.imshow(img, cmap='gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(1, 2, 2), plt.imshow(recon, cmap='gray'), plt.title('FFT Shift log'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

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))   
        kernel = (kernel - kernel.min()) / (kernel.max() - kernel.min())
    return kernel

def plot_3d(ax, x, y, z):
    ax.plot_surface(x, y, z, antialiased=True, shade=True)
    ax.view_init(30, 60), ax.grid(b=False), ax.set_xticks([]), ax.set_yticks([]), ax.set_zticks([])
    
def imshow_img(img, ax, cmap='gray'):
    ax.imshow(img, cmap=cmap), ax.set_xticks([]), ax.set_yticks([])

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 a uint8 [0, 255] reconstruct image
    """
    # 滤波
    res = fshift * h

    # 先反中以化，再反变换
    f1shift = np.fft.ifftshift(res)
    ifft = np.fft.ifft2(f1shift)

    # 取实数部分
    recon = np.real(ifft)
    # 小于0都被置为0
    # recon = np.where(recon > 0, recon, 0)
    recon = np.clip(recon, 0, recon.max())
    recon = np.uint8(normalize(recon) * 255)
    
    return recon

# 高斯核滤波器与对应的结果
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm

# 显示多个3D高斯核函数
k = 3
sigma = 0.8
X = np.linspace(-k, k, 200)
Y = np.linspace(-k, k, 200)
x, y = np.meshgrid(X, Y)
z = np.exp(- ((x)**2 + (y)**2)/ (2 * sigma**2))

fig = plt.figure(figsize=(15, 10))
ax_1 = fig.add_subplot(2, 3, 1, projection='3d')
plot_3d(ax_1, x, y, z)

ax_2 = fig.add_subplot(2, 3, 2, projection='3d')
plot_3d(ax_2, x, y, 1 - z)

# 这里只显示不明显，
a = 1
x = a + x
y = a + y
ax_3 = fig.add_subplot(2, 3, 3, projection='3d')
plot_3d(ax_3, x, y, 1 - z)

img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0431(d)(blown_ic_crop).tif', -1)

# FFT
img_fft = np.fft.fft2(img_ori.astype(np.float32))

# 中心化
fshift = np.fft.fftshift(img_fft)            # 将变换的频率图像四角移动到中心

# 滤波器
radius = 20
center = img_ori.shape[:2]
h_high_1 = gauss_high_pass_filter(img_ori, center, radius=radius)
h_low = 1 - h_high_1
center = (1026, 872)
h_high_2 = gauss_high_pass_filter(img_ori, center, radius=radius)

img_low = frequency_filter(fshift, h_low)
img_high_1 = frequency_filter(fshift, h_high_1)
img_high_2 = frequency_filter(fshift, h_high_2)

ax_4 = fig.add_subplot(2, 3, 4)
imshow_img(img_low, ax_4)

ax_5 = fig.add_subplot(2, 3, 5)
imshow_img(img_high_1, ax_5)

ax_6 = fig.add_subplot(2, 3, 6)
imshow_img(img_high_2, ax_6)

plt.tight_layout()
plt.show()

下面展示的是函数未填充时，会产生交叠错误；而函数填充0后，产生我们期望的结果，但是此时的滤波器应该是原来的2倍大小，如果是还是原来的大小的话，那么图像会比之前更模糊的。

这里采用的方法是对图像进行零填充，再直接在频率域中创建期望的滤波器传递函数，该函数的到底会给你个填充零后的图像的大小相同。这个明显降低交叠错误，只是会出现振铃效应。

# 未填充与填充0的结果
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0432(a)(square_original).tif', -1)

#================未填充===================
# FFT
img_fft = np.fft.fft2(img_ori.astype(np.float32))

# 中心化
fshift = np.fft.fftshift(img_fft)            # 将变换的频率图像四角移动到中心

# 滤波器
radius = 50
center = img_ori.shape[:2]
h_high_1 = gauss_high_pass_filter(img_ori, center, radius=radius)
h_low = 1 - h_high_1

img_low = frequency_filter(fshift, h_low)

fig = plt.figure(figsize=(15, 5))
ax_1 = fig.add_subplot(1, 3, 1)
imshow_img(img_ori, ax_1)

ax_2 = fig.add_subplot(1, 3, 2)
imshow_img(img_low, ax_2)

#================填充0====================
fp = np.zeros([img_ori.shape[0]*2, img_ori.shape[1]*2])
fp[:img_ori.shape[0], :img_ori.shape[1]] = img_ori

img_fft = np.fft.fft2(fp.astype(np.float32))
fshift = np.fft.fftshift(img_fft)   

hp = 1 - gauss_high_pass_filter(fp, fp.shape, radius=100)

img_low = frequency_filter(fshift, hp)
img_dst = img_low[:img_ori.shape[0], :img_ori.shape[1]]

ax_3 = fig.add_subplot(1, 3, 3)
imshow_img(img_dst, ax_3)

plt.tight_layout()
plt.show()

另一个合理的做法是：构建一个与未填充图像大小相同的函数，计算该函数的IDFT得到相应的空间表示，在空间域中对这一表示填充，然后计算其DFT返回到频率域。这也会出现振铃效应

def centralized(arr):
    """
    centralized 1d array $f(x) (-1)^{n}$
    """
    u = arr.shape[0]
    dst = np.zeros(arr.shape)
    
    # 中心化
    for n in range(u):
        dst[n] = arr[n] * ((-1) ** n)
    return dst

# 频率域反变换到空间域，再填充，再DFT，会出现振铃效应
h = np.zeros([256])
h[:9] = 1
h_c = np.fft.fftshift(h)

fig = plt.figure(figsize=(15, 10))
grid = plt.GridSpec(2, 3, wspace=0.2, hspace=0.1)
ax_1 = fig.add_subplot(grid[0, 0])
ax_1.plot(h_c), ax_1.set_xticks([0, 128, 255]), ax_1.set_yticks(np.arange(-0.2, 1.4, 0.2)), ax_1.set_xlim([0, 255])

# 用公式法的中以化才能得到理想的结果
h_pad = np.pad(h, (0, 256), mode='constant', constant_values=0)
ifshift = centralized_2d(h_pad.reshape(1, -1)).flatten() # 也可以用centralized_2d这个函数
ifft = np.fft.ifft(ifshift)
ifft = ifft.real * 2
ifft[:128] = 0
ifft[384:] = 0
ax_2 = fig.add_subplot(grid[0, 1:3])
ax_2.plot(ifft), ax_2.set_xticks([0, 128, 256, 384, 511]), ax_2.set_yticks(np.arange(-0.01, 0.05, 0.01)), ax_2.set_xlim([0, 511])

f = ifft[128:384]
ax_3 = fig.add_subplot(grid[1, 0])
ax_3.plot(f), ax_3.set_xticks([0, 128, 255]), ax_3.set_yticks(np.arange(-0.01, 0.05, 0.01)), ax_3.set_xlim([0, 255])

# 已经中心化，截取部分构成原函数，但出现的效果不是太好，但可以看出振铃效应
d = 120
f_ = np.zeros([h.shape[0]])
f_[:256-d] = f[d:]
f_pad = np.pad(f_, (0, 256), mode='constant', constant_values=0)
f_centralized = centralized(f_pad)
fft = np.fft.fft(f_centralized)
ax_4 = fig.add_subplot(grid[1, 1:3])
ax_4.plot(fft.real), ax_4.set_xticks([0, 128, 256, 384, 511]), ax_4.set_yticks(np.arange(-0.2, 1.4, 0.2)), ax_4.set_xlim([0, 511])

plt.show()

零相移滤波器

• 相角计算为复数的实部与虚部之比的反正切。因为$H(u,v)$同是乘以$R,I$，所以当这个比率形成后，它会被抵消。对实部和虚部的影响相同且对相角无影响的滤波器。

# 相角对反变换的影响
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0431(d)(blown_ic_crop).tif', -1)

#================未填充===================
# FFT
img_fft = np.fft.fft2(img_ori.astype(np.float32))

# 中心化
fshift = np.fft.fftshift(img_fft)            

# 频谱与相角
spectrum = spectrum_fft(fshift)
phase = phase_fft(fshift)

# 相角乘以-1
phase_1 = phase * -1
new_fshift = spectrum * np.exp(phase_1 * -1j)
new_fshift = np.fft.ifftshift(new_fshift)
recon_1 = np.fft.fft2(new_fshift)

# 相角乘以0.5，需要乘以一个倍数，以便可视化
phase_2 = phase * 0.25
new_fshift = spectrum * np.exp(phase_2 * -1j)
new_fshift = np.fft.ifftshift(new_fshift)
recon_2 = np.fft.fft2(new_fshift)
recon_2 = np.uint8(normalize(recon_2.real) * 4000)

fig = plt.figure(figsize=(15, 5))
ax_1 = fig.add_subplot(1, 3, 1)
imshow_img(img_ori, ax_1)

ax_2 = fig.add_subplot(1, 3, 2)
imshow_img(recon_1.real, ax_2)

ax_3 = fig.add_subplot(1, 3, 3)
ax_3.imshow(recon_2, cmap='gray')

plt.tight_layout()
plt.show()

频率域滤波步骤小结

1. 已知一幅大小为$M\times N$的输入图像$f(x,y)$，分别得到填充零后的尺寸$P$和$Q$，即$P=2M$和$Q=2N$。
2. 使用零填充、镜像填充或复制填充，形成大小为$P\times Q$的填充后的图像$f_P(x,y)$。
3. 将$f_P(x,y)$乘以$(-1{)}^{x+y}$，使用傅里叶变换位于$P\times Q$大小的频率矩形的中心。
4. 计算步骤3得到的图像的DFT，即$F(u,v)$。
5. 构建一个实对称滤波器传递函数$H(u,v)$，其大小为$P\times Q$，中心中$(P/2,Q/2)$处。
6. 采用对应像素相乘得到$G(u,v)=H(u,v)F(u,v)$。
7. 计算$G(u,v)$的IDFT得到滤波后的图像（大小为$P\times Q$），$g_P(x,y)=\left\{\text{Real}\right[{\mathfrak{J}}^{-1}[G(u,v)]\left]\right\}(-1{)}^{x+y}$。
8. 最后，从$g_P(x,y)$的左上象限提取一个大小为$M\times N$的区域，得到与输入图像大小相同的滤波后的结果$g(x,y)$。

def pad_image(img, mode='constant'):
    """
    pad image into PxQ shape, orginal is in the top left corner.
    param: img: input image
    param: mode: input, str, is numpy pad parameter, default is 'constant'. for more detail please refere to Numpy pad
    return PxQ shape image padded with zeros or other values
    """
    dst = np.pad(img, ((0, img.shape[0]), (0, img.shape[1])), mode=mode)
    return dst    

# 频率域滤波过程
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0431(d)(blown_ic_crop).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 = 1 - gauss_high_pass_filter(fp, center=fp.shape, radius=70)
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]
ax_8 = fig.add_subplot(3, 3, 8)
imshow_img(g, ax_8)

plt.tight_layout()
plt.show()

def frequency_filter(fft, h):
    """
    Frequency domain filtering using FFT
    param: fshift: input, FFT shift to center
    param: h: input, filter, value range [0, 1]
    return a uint8 [0, 255] reconstruct image
    """
    # 滤波
    res = fft * h

    # 先反中以化，再反变换
    ifft = np.fft.ifft2(res)

    # 取实数部分
    recon = centralized_2d(ifft.real)
    # 小于0都被置为0
    recon = np.clip(recon, 0, recon.max())
    recon = np.uint8(normalize(recon) * 255)
    
    return recon

# 频率域滤波
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0431(d)(blown_ic_crop).tif', -1)
M, N = img_ori.shape[:2]

# 填充
fp = pad_image(img_ori, mode='reflect')
# 中心化
fp_cen = centralized_2d(fp)
# 正变换
fft = np.fft.fft2(fp_cen)
# 滤波器
H = 1 - gauss_high_pass_filter(fp, center=fp.shape, radius=100)
# 滤波
HF = fft * H
# 反变换
ifft = np.fft.ifft2(HF)
# 去中心化
gp = centralized_2d(ifft.real)
# 还回来与原图像的大小
g = gp[:M, :N]
fig = plt.figure(figsize=(15, 15))
ax_8 = fig.add_subplot(3, 3, 8)
imshow_img(g, ax_8)
plt.tight_layout()
plt.show()

空间域和频率域滤波之间的对应关系

空间域中的滤波器核与频率域傅里叶变换对：
$$h(x,y)\iff H(u,v)\text{\hspace{1em}(4.106)}$$

$h(x,y)$是空间核。这个核 可以由一个频率域滤波器对一个冲激的响应得到，因此有时我们称$h(x,y)$是$H(u,v)$的冲激响应。

一维频率域高斯传递函数
$$H(u)=A{e}^{-u^2/2{\sigma }^2}\text{\hspace{1em}(4.107)}$$
空间域的核 由$H(u)$的傅里叶反变换得到
$$h(x)=\sqrt{2\pi }\sigma A{e}^{-2{\pi }^2{\sigma }^2{x}^2}\text{\hspace{1em}(4.108)}$$

这是两个很重要的公式

1. 它们是一个傅里叶变换对，两者都是高斯函数和实函数
2. 两个函数的表现相反

• 高斯差

  • 用一个常量减去低通函数，可以构建一个高通滤波器。
  • 因为使用所谓的高斯差（涉及两个低通函数）就可更好地控制滤波函数的形状

• 频率域 $A\ge B,{\sigma }_1>{\sigma }_2$
  $$H(u)=A{e}^{-u^2/2{\sigma }_1^2}-B{e}^{-u^2/2{\sigma }_2^2}\text{\hspace{1em}(4.109)}$$

• 空间域
  $$h(x)=\sqrt{2\pi }{\sigma }_{1}A{e}^{-2{\pi }^2{\sigma }_1^2{x}^2}-\sqrt{2\pi }{\sigma }_{2}B{e}^{-2{\pi }^2{\sigma }_2^2{x}^2}\text{\hspace{1em}(4.110)}$$

# 空间域与频率域的滤波器对应关系
u = np.linspace(-3, 3, 200)

A = 1
sigma = 0.2
Hu = A * np.exp(-np.power(u, 2) / (2 * np.power(sigma, 2)))

x = u
hx = np.sqrt(2 * np.pi) * sigma * A * np.exp(-2 * np.pi**2 * sigma**2 * x**2)

fig = plt.figure(figsize=(10, 10))
ax_1 = setup_axes(fig, 221)
ax_1.plot(u, Hu), ax_1.set_title('H(u)', loc='center', y=1.05), ax_1.set_ylim([0, 1]), ax_1.set_yticks([])

B = 1
sigma_2 = 4
Hu_2 = B * np.exp(-np.power(u, 2) / (2 * np.power(sigma_2, 2)))
Hu_diff = Hu_2 - Hu

ax_2 = setup_axes(fig, 222)
ax_2.plot(x, Hu_diff), ax_2.set_title('H(u)', loc='center', y=1.05), ax_2.set_ylim([0, 1.1]), ax_2.set_yticks([])

ax_3 = setup_axes(fig, 223)
ax_3.plot(x, hx), ax_3.set_title('h(x)', loc='center', y=1.05), ax_3.set_ylim([0, 0.8]), ax_3.set_yticks([])

hx_2 = np.sqrt(2 * np.pi) * sigma_2 * B * np.exp(-2 * np.pi**2 * sigma_2**2 * x**2)
hx_diff = hx_2 - hx

ax_4 = setup_axes(fig, 224)
ax_4.plot(x, hx_diff), ax_4.set_title('h(x)', loc='center', y=1.05), ax_4.set_ylim([-1, 10]), ax_4.set_yticks([])

plt.tight_layout()
plt.show()

# 空间核得到频率域
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0438(a)(bld_600by600).tif', -1)

fp = pad_zero(img_ori)
fp_cen = centralized_2d(fp)
fft = np.fft.fft2(fp_cen)
spectrum = spectrum_fft(fft)
spectrum = np.log(1 + spectrum)

plt.figure(figsize=(18, 15))
plt.subplot(1, 2, 1), plt.imshow(img, cmap='gray'), plt.title('No Log')
plt.subplot(1, 2, 2), plt.imshow(spectrum, cmap='gray'), plt.title('FFT')
plt.tight_layout()
plt.show()

注：频率域的Sobel滤波器还没有解决

# 频域
center = (fp.shape[0], fp.shape[1] + 300)
GHPF = gauss_high_pass_filter(fp, center, radius=200)
GHPF = (1 - np.where(GHPF > 1e-6, GHPF, 1)) * -2

center = (fp.shape[0], fp.shape[1] - 300)
GHPF_1 = gauss_high_pass_filter(fp, center, radius=200)
GHPF_1 = (1 - np.where(GHPF_1 > 1e-6, GHPF_1, 1)) * 2

sobel_f = GHPF + GHPF_1

plt.figure(figsize=(15, 5))

plt.subplot(1, 3, 1), plt.imshow(GHPF, 'gray')
plt.subplot(1, 3, 2), plt.imshow(GHPF_1, 'gray')
plt.subplot(1, 3, 3), plt.imshow(sobel_f, 'gray')

plt.tight_layout()
plt.show()

# 频域
fig = plt.figure(figsize=(15, 5))
HF = fft * sobel_f
HF_spectrum = np.log(1 + spectrum_fft(HF))
ax_6 = fig.add_subplot(1, 3, 1)
imshow_img(HF_spectrum, ax_6)

ifft = np.fft.ifft2(HF)
gp = centralized_2d(ifft.real)
ax_7 = fig.add_subplot(1, 3, 2)
imshow_img(gp, ax_7)

# 最终结果有黑边
g = gp[:M, :N]
ax_8 = fig.add_subplot(1, 3, 3)
imshow_img(g, ax_8)

# plt.figure(figsize=(15, 5))

# plt.subplot(1, 3, 1), plt.imshow(GHPF, 'gray')
# plt.subplot(1, 3, 2), plt.imshow(GHPF_1, 'gray')
# plt.subplot(1, 3, 3), plt.imshow(sobel_f, 'gray')

plt.tight_layout()
plt.show()

空间域的Sobel滤波器

def kernel_seperate(kernel):
    """
    这里的分离核跟之前不太一样
    get kernel seperate vector if separable
    param: kernel: input kerel of the filter
    return two separable vector c, r
    """
    rank = np.linalg.matrix_rank(kernel)

    if rank == 1:

    #-----------------Numpy------------------
    # here for sobel kernel seperate
        c = kernel[:, -1]
        r = kernel[-1, :]

    else:
        print('Not separable')

#     c = np.array(c).reshape(-1, 1)
#     r = np.array(r).reshape(-1, 1)
    
    return c, r

def separate_kernel_conv2D(img, kernel):
    """
    separable kernel convolution 2D, if 2D kernel not separable, then canot use this fuction. in the fuction. first need to
    seperate the 2D kernel into two 1D vectors.
    param: img: input image you want to implement 2d convolution with 2d kernel
    param: kernel: 2D separable kernel, such as Gaussian kernel, Box Kernel
    return image after 2D convolution
    """
    m, n = kernel.shape[:2]
    pad_h = int((m - 1) / 2)
    pad_w = int((n - 1) / 2)
    w_1, w_2 = kernel_seperate(kernel)
#     w_1 = np.array([-1, 0, 1])
#     w_2 = np.array([1, 2, 1])
    convovle = np.vectorize(np.convolve, signature='(x),(y)->(k)')
    r_2 = convovle(img, w_2)
    r_r = np.rot90(r_2)
    r_1 = convovle(r_r, w_1)
    r_1 = r_1.T
    r_1 = r_1[:, ::-1]
    img_dst = img.copy().astype(np.float)
    img_dst = r_1[pad_h:-pad_h, pad_w:-pad_w]
    return img_dst

# Sobel梯度增强边缘
img_ori = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH04/Fig0438(a)(bld_600by600).tif", 0)

sobel_x = np.zeros([3, 3], np.int)
sobel_x[0, :] = np.array([-1, -2, -1])
sobel_x[2, :] = np.array([1, 2, 1])

sobel_y = np.zeros([3, 3], np.int)
sobel_y[:, 0] = np.array([-1, -2, -1])
sobel_y[:, 2] = np.array([1, 2, 1])

gx = separate_kernel_conv2D(img_ori, kernel=sobel_x)
gy = separate_kernel_conv2D(img_ori, kernel=sobel_y)

# gx = cv2.filter2D(img_ori, ddepth=-1, kernel=sobel_x)
# gy = cv2.filter2D(img_ori, ddepth=-1, kernel=sobel_y)

# 先对gx gy做二值化处理再应用下面的公式
img_sobel = np.sqrt(gx**2 + gy**2)   

# 二值化后，平方根的效果与绝对值很接近
img_sobel = abs(gx) + abs(gy)
img_sobel = np.uint8(normalize(img_sobel) * 255)

plt.figure(figsize=(15, 12))
plt.subplot(1, 2, 1), plt.imshow(img_ori,   'gray', vmax=255), plt.title("Original"), plt.xticks([]), plt.yticks([])
plt.subplot(1, 2, 2), plt.imshow(img_sobel, 'gray', vmax=255), plt.title("Sobel"), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()