Update image-processing lab: implementation of fft

Author: heqin-zhu

计算机图像学/labs/dft.py

@@ -0,0 +1,119 @@
+''' mbinary
+#########################################################################
+# File : fft.py
+# Author: mbinary
+# Mail: [email protected]
+# Blog: https://mbinary.xyz
+# Github: https://github.com/mbinary
+# Created Time: 2019-06-11  12:48
+# Description:
+#########################################################################
+'''
+import numpy as np
+
+
+def _fft_n2(a, invert):
+    '''O(n^2)'''
+    N = len(a)
+    w = np.arange(N)
+    i = 2j if invert else -2j
+    m = w.reshape((N, 1)) * w
+    W = np.exp(m * i * np.pi / N)
+    return np.concatenate(np.dot(W, a.reshape((N,
+                                               1))))  # important, cannot use *
+
+
+def _fft(a, invert=False):
+    '''recursion version'''
+    N = len(a)
+    if N == 1:
+        return [a[0]]
+    elif N & (N - 1) == 0:  # O(nlogn),  2^k
+        even = _fft(a[::2], invert)
+        odd = _fft(a[1::2], invert)
+        i = 2j if invert else -2j
+        factor = np.exp(i * np.pi * np.arange(N // 2) / N)
+        prod = factor * odd
+        return np.concatenate([even + prod, even - prod])
+    else:
+        return _fft_n2(a, invert)
+
+
+def _fft2(a, invert=False):
+    ''' iteration version'''
+
+    def rev(x):
+        ret = 0
+        for i in range(r):
+            ret <<= 1
+            if x & 1:
+                ret += 1
+            x >>= 1
+        return ret
+
+    N = len(a)
+    if N & (N - 1) == 0:  # O(nlogn),  2^k
+        r = int(np.log(N))
+        c = np.array(a,dtype='complex')
+        i = 2j if invert else -2j
+        w = np.exp(i * np.pi / N)
+        for h in range(r - 1, -1, -1):
+            p = 2**h
+            z = w**(N / p / 2)
+            for k in range(N):
+                if k % p == k % (2 * p):
+                    c[k], c[k + p] = c[k] + c[k + p], c[k] * z**(k % p)
+
+        return np.asarray([c[rev(i)] for i in range(N)])
+    else:  # O(n^2)
+        return _fft_n2(a, invert)
+
+
+def fft(a):
+    '''fourier[a]'''
+    n = len(a)
+    if n == 0:
+        raise Exception("[Error]: Invalid length: 0")
+    return _fft(a)
+
+
+def ifft(a):
+    '''invert fourier[a]'''
+    n = len(a)
+    if n == 0:
+        raise Exception("[Error]: Invalid length: 0")
+    return _fft(a, True) / n
+
+
+def fft2(arr):
+    return np.apply_along_axis(fft, 0,
+                               np.apply_along_axis(fft, 1, np.asarray(arr)))
+
+
+def ifft2(arr):
+    return np.apply_along_axis(ifft, 0,
+                               np.apply_along_axis(ifft, 1, np.asarray(arr)))
+
+
+def test(n=128):
+    print('\nsequence length:', n)
+    print('fft')
+    li = np.random.random(n)
+    print(np.allclose(fft(li), np.fft.fft(li)))
+
+    print('ifft')
+    li = np.random.random(n)
+    print(np.allclose(ifft(li), np.fft.ifft(li)))
+
+    print('fft2')
+    li = np.random.random(n * n).reshape((n, n))
+    print(np.allclose(fft2(li), np.fft.fft2(li)))
+
+    print('ifft2')
+    li = np.random.random(n * n).reshape((n, n))
+    print(np.allclose(ifft2(li), np.fft.ifft2(li)))
+
+
+if __name__ == '__main__':
+    for i in range(1, 4):
+        test(i * 16)

计算机图像学/labs/lab1.py

@@ -103,7 +103,7 @@ def show(img, s='opencv'):
 
     cmap = mpl.cm.gray # mpl.cm.gray_r   'gray'
 
-    plt.figure(figsize=(10, 10))
+    plt.figure(figsize=(8, 8))
 
     plt.subplot(321), plt.imshow(img,cmap=cmap), plt.title('origin'),plt.xticks([]), plt.yticks([])
     plt.subplot(322), plt.imshow(img2,cmap=cmap), plt.title(f'tran k={k},b={b}'),plt.xticks([]), plt.yticks([])

计算机图像学/labs/lab2.py

@@ -142,7 +142,7 @@ def find_median(arr):
     img3 = median_filter(noised_img)
 
     cmap = mpl.cm.gray
-    plt.figure(figsize=(10, 10))
+    plt.figure(figsize=(8, 8))
 
     plt.subplot(221), plt.xticks([]), plt.yticks([])
     plt.imshow(img, cmap=cmap)

计算机图像学/labs/lab3.py

@@ -46,7 +46,7 @@ def prewitt(img):
     img2 = roberts(img)
     img3 = prewitt(img)
     cmap = mpl.cm.gray
-    plt.figure(figsize=(10, 10))
+    plt.figure(figsize=(8, 8))
 
     plt.subplot(221),plt.xticks([]), plt.yticks([])
     plt.imshow(img, cmap=cmap)

计算机图像学/labs/lab4.py

@@ -4,36 +4,33 @@
 import numpy as np
 import matplotlib as mpl
 from matplotlib import pyplot as plt
+import dft
 
 
-def fft(img):
-    f = np.fft.fft2(img)
-    # fshift = np.fft.fftshift(f)
-    # magnitude_spectrum
-    mag = np.abs(f)
-    imag = np.fft.ifft2(mag)
-    phase = np.angle(f)
-    iphase = np.fft.ifft2(phase)
-    return mag, phase, imag, iphase
 
 
 if __name__ == '__main__':
     path = sys.argv[1]
     img = cv2.imread(path, cv2.IMREAD_GRAYSCALE)
 
-    mag, phase, imag, iphase = fft(img)
+    f = dft.fft2(img)
+    # magnitude_spectrum
+    mag = np.abs(f)
+    invert_magnitude = dft.ifft2(mag)
     mag = np.log(mag + 1)
-    imag = np.real(imag)
+    phase = np.angle(f)
+    iphase = dft.ifft2(phase)
+    invert_magnitude = np.real(invert_magnitude)
     iphase = np.real(iphase)
 
     cmap = mpl.cm.gray
-    plt.figure(figsize=(10, 10))
+    plt.figure(figsize=(6,6))
 
     plt.subplot(221), plt.imshow(img, cmap=cmap)
     plt.title('Input Image'), plt.xticks([]), plt.yticks([])
     plt.subplot(222), plt.imshow(mag, cmap=cmap)
     plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
-    plt.subplot(223), plt.imshow(imag, cmap=cmap)
+    plt.subplot(223), plt.imshow(invert_magnitude, cmap=cmap)
     plt.title('idft [ magnitude] '), plt.xticks([]), plt.yticks([])
     plt.subplot(224), plt.imshow(iphase, cmap=cmap)
     plt.title('idft [ phase ]'), plt.xticks([]), plt.yticks([])