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
- 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.
- 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).
- 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.
- Calculate the DFT of the image obtained in step 3, i.e F ( u , v ) F(u, v) F(u,v).
- 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).
- 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).
- 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.
- 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
- They are a Fourier transform pair, both of which are Gaussian functions and real functions
- 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π σ1Ae−2π2σ12x2−2π σ2Be−2π2σ22x2(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()