Fundamentals of frequency domain filtering

Other characteristics in frequency domain

• The filtering process in the frequency domain is as follows:
  • First, the Fourier transform is modified to achieve a specific purpose
  • Then calculate the IDFT and return to the spatial domain

# Other characteristics in frequency domain
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))

# Decentralized spectrum
spectrum = spectrum_fft(img_fft)
spectrum = np.uint8(normalize(spectrum) * 255)

# Centralization
fshift = np.fft.fftshift(img_fft)            # Move the four corners of the transformed frequency image to the center

# Centralized spectrum
spectrum_fshift = spectrum_fft(fshift)
spectrum_fshift_n = np.uint8(normalize(spectrum_fshift) * 255)

# Logarithmic transformation of spectrum
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()

Basic knowledge of frequency domain filtering

Known size is P × Q P \times Q P × A digital image of Q pixels f ( x , y ) f(x, y) f(x,y), then the basic filtering formula we are interested in is:
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)

J − 1 \mathfrak{J}^{-1} J − 1 is IDFT, H ( u , v ) H(u, v) H(u,v) is the filter transfer function (often referred to as filter or filter function), F ( u , v ) F(u, v) F(u,v) is the input image f ( x , y ) f(x, y) DFT of f(x,y), g ( x , y ) g(x, y) g(x,y) is the filtered image.

Zhongyihua: use ( − 1 ) x + y (-1)^{x + y} (− 1)x+y multiplied by the input image

• low pass filter

  • It will attenuate the high frequency and blur the image through the function of low frequency

• High pass filter

  • It will make the details of the image clearer, but it will reduce the contrast of the image

def centralized_2d_loop(arr):
    """
    centralized 2d array $f(x, y) (-1)^{x + y}$
    """    
    #==================Centralization=============================
    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
        # Here, judge whether a m, n is even. If it is even, the result is incorrect and needs to be transformed into an odd number
        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

    #==================Centralization=============================
    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

    #==================Centralization=============================
    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

# Basic knowledge of frequency domain filtering
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))

# Centralization
fshift = np.fft.fftshift(img_fft)            # Move the four corners of the transformed frequency image to the center

# wave filter
h = idea_high_pass_filter(img, radius=1)

# wave filtering
res = fshift * h

# First reverse to transform, and then reverse transformation
f1shift = np.fft.ifftshift(res)
ifft = np.fft.ifft2(f1shift)

# Real part
recon = np.real(ifft)
# Less than 0 is set to 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
    """
    # wave filtering
    res = fshift * h

    # First reverse to transform, and then reverse transformation
    f1shift = np.fft.ifftshift(res)
    ifft = np.fft.ifft2(f1shift)

    # Real part
    recon = np.real(ifft)
    # Less than 0 is set to 0
    # recon = np.where(recon > 0, recon, 0)
    recon = np.clip(recon, 0, recon.max())
    recon = np.uint8(normalize(recon) * 255)
    
    return recon

# Gaussian kernel filter and corresponding results
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm

# Display multiple 3D Gaussian kernels
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)

# It's just not obvious here,
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))

# Centralization
fshift = np.fft.fftshift(img_fft)            # Move the four corners of the transformed frequency image to the center

# wave filter
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()

The following shows that overlapping errors will occur when the function is not filled; After the function is filled with 0, it will produce the desired result, but the filter should be twice the original size. If it is still the original size, the image will be more blurred than before.

The method used here is to fill the image with zero, and then directly create the desired filter transfer function in the frequency domain. The function will give you an image with the same size after filling zero. This significantly reduces overlapping errors, but there is a ringing effect.

# Unfilled and 0 filled results
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0432(a)(square_original).tif', -1)

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

# Centralization
fshift = np.fft.fftshift(img_fft)            # Move the four corners of the transformed frequency image to the center

# wave filter
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)

#================Fill 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()

Another reasonable approach is to construct a function with the same size as the unfilled image, calculate the IDFT of the function to obtain the corresponding spatial representation, fill this representation in the spatial domain, and then calculate its DFT to return to the frequency domain. This will also have a ringing effect

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

# If the frequency domain is inversely transformed to the spatial domain, then filled, and then DFT, the ringing effect will appear
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])

# The ideal result can be obtained by the neutralization of formula method
h_pad = np.pad(h, (0, 256), mode='constant', constant_values=0)
ifshift = centralized_2d(h_pad.reshape(1, -1)).flatten() # You can also use centralized_2d this function
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])

# It has been centralized, and the intercepted part constitutes the original function, but the effect is not very good, but the ringing effect can be seen
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()

Zero phase shift filter

• The phase angle is calculated as the arctangent of the ratio of the real part to the imaginary part of the complex number. because H ( u , v ) H(u, v) H(u,v) is also multiplied by R , I R,I R. I, so when this ratio is formed, it will be offset. A filter that has the same effect on the real part and imaginary part and has no effect on the phase angle.

# Influence of phase angle on inverse transformation
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0431(d)(blown_ic_crop).tif', -1)

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

# Centralization
fshift = np.fft.fftshift(img_fft)            

# Spectrum and phase angle
spectrum = spectrum_fft(fshift)
phase = phase_fft(fshift)

# Phase angle multiplied by - 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)

# If the phase angle is multiplied by 0.5, it needs to be multiplied by a multiple for visualization
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()

Summary of filtering steps in frequency domain

1. A known size is M × N M\times N M × Input image of N f ( x , y ) f(x, y) f(x,y), respectively obtain the size after filling zero P P P and Q Q Q. Namely P = 2 M P=2M P=2M and Q = 2 N Q=2N Q=2N.
2. Use zero fill, mirror fill, or duplicate fill to form a size of P × Q P\times Q P × Q filled image f P ( x , y ) f_P(x, y) fP​(x,y).
3. take f P ( x , y ) f_P(x, y) fP (x,y) times ( − 1 ) x + y (-1)^{x+y} (− 1)x+y, using Fourier transform at P × Q P\times Q P × Q is the size of the center of the frequency rectangle.
4. Calculate the DFT of the image obtained in step 3, i.e F ( u , v ) F(u, v) F(u,v).
5. A transfer function of real symmetric filter is constructed H ( u , v ) H(u, v) H(u,v), whose size is P × Q P\times Q P × Q. Center ( P / 2 , Q / 2 ) (P/2, Q/2) (P/2,Q/2).
6. Obtained by multiplying the corresponding pixels G ( u , v ) = H ( u , v ) F ( u , v ) G(u, v) = H(u, v)F(u, v) G(u,v)=H(u,v)F(u,v).
7. calculation G ( u , v ) G(u, v) The IDFT of G(u,v) obtains the filtered image (size: P × Q P\times Q P×Q)， g P ( x , y ) = { Real [ J − 1 [ G ( u , v ) ] ] } ( − 1 ) x + y g_P(x, y) = \Big\{\text{Real} \big[\mathfrak{J}^{-1}[G(u, v)] \big]\Big\}(-1)^{x+y} gP​(x,y)={Real[J−1[G(u,v)]]}(−1)x+y.
8. Finally, from g P ( x , y ) g_P(x, y) The upper left quadrant of gP (x,y) extracts a size of M × N M\times N M × N, the filtered result with the same size as the input image is obtained g ( x , y ) g(x, y) 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    

# Frequency domain filtering process
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))

# Here, the original image is pad to get better display
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) # Negative numbers are truncated to 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)

# The end result has a black edge
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
    """
    # wave filtering
    res = fft * h

    # First reverse to transform, and then reverse transformation
    ifft = np.fft.ifft2(res)

    # Real part
    recon = centralized_2d(ifft.real)
    # Less than 0 is set to 0
    recon = np.clip(recon, 0, recon.max())
    recon = np.uint8(normalize(recon) * 255)
    
    return recon

# Frequency domain filtering
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0431(d)(blown_ic_crop).tif', -1)
M, N = img_ori.shape[:2]

# fill
fp = pad_image(img_ori, mode='reflect')
# Centralization
fp_cen = centralized_2d(fp)
# Positive transformation
fft = np.fft.fft2(fp_cen)
# wave filter
H = 1 - gauss_high_pass_filter(fp, center=fp.shape, radius=100)
# wave filtering
HF = fft * H
# Inverse transformation
ifft = np.fft.ifft2(HF)
# Decentralization
gp = centralized_2d(ifft.real)
# Back with the size of the original image
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()

The corresponding relationship between spatial domain and frequency domain filtering

Filter core in spatial domain and Fourier transform pair in frequency domain:
h ( x , y ) ⇔ H ( u , v ) (4.106) h(x, y) \Leftrightarrow H(u, v) \tag{4.106} h(x,y)⇔H(u,v)(4.106)

h ( x , y ) h(x, y) h(x,y) is the space kernel. This kernel can be obtained from the response of a frequency domain filter to an impulse, so it is sometimes called h ( x , y ) h(x, y) h(x,y) is H ( u , v ) H(u, v) Impulse response of H(u,v).

Gaussian transfer function in one-dimensional frequency domain
H ( u ) = A e − u 2 / 2 σ 2 (4.107) H(u) = A e^{-u^2/2\sigma^2} \tag{4.107} H(u)=Ae−u2/2σ2(4.107)
Kernel of spatial domain H ( u ) H(u) The inverse Fourier transform of H(u) is obtained
h ( x ) = 2 π σ A e − 2 π 2 σ 2 x 2 (4.108) h(x) = \sqrt{2\pi}\sigma A e^{-2\pi^2\sigma^2 x^2} \tag{4.108} h(x)=2π ​σAe−2π2σ2x2(4.108)

These are two very important formulas

1. They are a Fourier transform pair, both of which are Gaussian functions and real functions
2. The performance of the two functions is opposite

• Gaussian difference

  • A high pass filter can be constructed by subtracting the low-pass function from a constant.
  • Because the shape of the filter function can be better controlled by using the so-called Gaussian difference (involving two low-pass functions)

• Frequency domain A ≥ B , σ 1 > σ 2 A \ge B, \sigma_1 > \sigma_2 A≥B,σ1​>σ2​
  H ( u ) = A e − u 2 / 2 σ 1 2 − B e − u 2 / 2 σ 2 2 (4.109) H(u) = A e^{-u^2/2\sigma_1^2} - B e^{-u^2/2\sigma_2^2}\tag{4.109} H(u)=Ae−u2/2σ12​−Be−u2/2σ22​(4.109)

• Spatial domain
  h ( x ) = 2 π σ 1 A e − 2 π 2 σ 1 2 x 2 − 2 π σ 2 B e − 2 π 2 σ 2 2 x 2 (4.110) 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} \tag{4.110} h(x)=2π ​σ1​Ae−2π2σ12​x2−2π ​σ2​Be−2π2σ22​x2(4.110)

# Corresponding relationship between filters in spatial domain and frequency domain
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()

# Spatial kernel frequency domain
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()

Note: the Sobel filter in frequency domain has not been solved

# frequency domain
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()

# frequency domain
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)

# The end result has a black edge
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 filter in spatial domain

def kernel_seperate(kernel):
    """
    The separated nuclei here are not the same as before
    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 gradient enhanced edge
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)

# First binarize GX and Gy, and then apply the following formula
img_sobel = np.sqrt(gx**2 + gy**2)   

# After binarization, the effect of square root is very close to the absolute value
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()