空間周波数フィルタリング

基本関数

min max scaler

import cv2
import numpy as np
import math

def _min_max(x):
    return (x - x.min()) / (x.max() - x.min())

FFT

def make_FFT(img):
    # FFT
    dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
    # ゼロ周波数の成分を中心に移動
    dft_shift = np.fft.fftshift(dft)
    # パワースペクトル
    magnitude_spectrum = 20 * np.log((cv2.magnitude(dft_shift[:, :, 0], dft_shift[:, :, 1])) + 1)

    return dft_shift, magnitude_spectrum

Innverse FFT

def make_IFFT(fshift):
    #fshift is not magnitude
    # Inverse FFT
    f_ishift = np.fft.ifftshift(fshift)
    img_back = cv2.idft(f_ishift)
    img_back = cv2.magnitude(img_back[:, :, 0], img_back[:, :, 1])

    return img_back

Masking parts

def make_mask(img, dft_shift, freq, LowPass=False):

    #円形のマスク画像を作ります
    rows, cols = img.shape
    crow, ccol = int(rows / 2), int(cols / 2)

    # mask radius r
    if rows >= cols:
        r = int(0.5 * freq * cols)
    else:
        r = int(0.5 * freq * rows)

    center = [crow, ccol]
    x, y = np.ogrid[:rows, :cols]

    #direct
    mask_area = np.sqrt((x - center[0]) ** 2 + (y - center[1]) ** 2) <= r

    if LowPass == False:
        # 内側 0, 外側 1
        mask = np.ones((rows, cols, 2), np.uint8)
        mask[mask_area] = 0
    elif LowPass == True:
        # 内側 1, 外側 0
        mask = np.zeros((rows, cols, 2), np.uint8)
        mask[mask_area] = 1

    #マスク処理

    fshift = dft_shift * mask
    mag = 20 * np.log((cv2.magnitude(fshift[:, :, 0], fshift[:, :, 1])) + 1)

    return fshift, mag

テストコード

img = cv2.imread('./image/messi.jpg', 0)
print(img.max())
print(img.min())

#input image -> fft
dft_shift, img_mag = make_FFT(img=img)


#low pass filter
f_shift1, fshift_mask_mag1 = make_mask(img = img,dft_shift=dft_shift, freq=0.3, LowPass=True)
f_shift2, fshift_mask_mag2 = make_mask(img = img,dft_shift=dft_shift, freq=0.5, LowPass=True)

#band pass filter
f_shift3 = f_shift2 - f_shift1
fshift_mask_mag3 = fshift_mask_mag2 - fshift_mask_mag1

#IFFT
img_back1 = make_IFFT(fshift=f_shift1)
img_back2 = make_IFFT(fshift=f_shift2)
img_back3 = make_IFFT(fshift=f_shift3)


#Plot
cv2.imshow('original image', img)
cv2.imshow('fft image', (255*_min_max(img_mag)).astype(np.uint8))

cv2.imshow('masked image1', (255*_min_max(fshift_mask_mag1)).astype(np.uint8))
cv2.imshow('masked image2', (255*_min_max(fshift_mask_mag2)).astype(np.uint8))
cv2.imshow('masked image3', (255*_min_max(fshift_mask_mag3)).astype(np.uint8))

cv2.imshow('ifft masked_image1', (255*_min_max(img_back1)).astype(np.uint8))
cv2.imshow('ifft masked_image2', (255*_min_max(img_back2)).astype(np.uint8))
cv2.imshow('ifft masked_image3', (255*_min_max(img_back3)).astype(np.uint8))

cv2.waitKey(0)

結果

Original Image

Frequency Domain

BandPassFiltering

BandPassFiltered Image