![[OpenCV] 04-11. Image Transforms in OpenCV](https://img1.daumcdn.net/thumb/R750x0/?scode=mtistory2&fname=https%3A%2F%2Ft1.daumcdn.net%2Fcfile%2Ftistory%2F99CB754A5D78FF3A27)
< Fourier Transform >
이번 장에서는
- OpenCV를 사용하여 이미지의 Fourier Transform을 찾을 것이다.
- Numpy에서 FFT함수를 이용하기 위해
- Fourier Transform의 기능들
- 다음의 함수를 볼 것 이다 : cv2.dft(), cv2.idft() 등등
Theory
Fourier Transform은 다양한 필터의 주파수(frequency) 특성을 분석하는 데 사용된다. 이미지의 경우, 주파수 영역을 찾기 위해 2D Discrete Fourier Transfrom(DFT)를 사용한다. Fast Fourier Transform(FFT)이라는 고속 알고리즘이 이미지 처리 또는 신호 처리에서 확인해 볼 수 있다.
사인곡선의 신호의 경우,
더 직관적으로, 사인곡석의 신호의 경우, 짧은 시간에 진폭이 그렇게 빠르게 변화하면 고주파 신호라고 할 수 있다. 천천히 변화하면 저주파 신호이다. 같은 아이디어를 이미지로 확장할 수 있다. 이미지에서 진폭이 크게 변하는 곳은 어디일까? 가장자리 부분이나 노이즈가 있는 부분일 것이다. 그래서 다음과 같이, 가장자리와 노이즈는 이미지에서 고주파 부분이라고 할 수 있다. 진폭에 큰 변화가 없으면 저주파 성분인 것이다.
이제 Fourier Transform을 어떻게 찾을지 볼 것이다.
Fourier Transfrom in Numpy
먼저 Numpy를 사용하여 Fourier Transform을 찾는 법을 볼 것이다. Numpy는 이를 하기 위한 FFT 패키지가 있다. np.fft.fft2()는 복잡한 배열로 주파수 변환을 제공한다. 첫 번째 인자는 입력 이미지이고, 흑백스케일이다. 두 번째 인자는 임의로써 결과 배열의 크기를 정한다. 만약에 입력 이미지의 크기보다 크다면 입력 이미지는 FFT 연산이전에 제로-패딩을 거친다. 만약에 입력 이미지보다 작을 경우, 입력 이미지는 잘려진다(cropped). 아무 인자도 안넘겨진다면, 결과 배열의 크기는 입력과 동일할 것이다.
이제 결과를 받았으면, 0 주파수 성분은 왼쪽 위에 위치할 것이다. 만약 이를 가운데로 가지고오고 싶으면, 두 방향으로
import numpy as np
import cv2
import matplotlib.pyplot as plt
img = cv2.imread('./images/face2.jpg',0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
magnitude_spectrum = 20*np.log(np.abs(fshift))
plt.subplot(121), plt.imshow(img,cmap='gray')
plt.title("Input Image"), plt.axis('off')
plt.subplot(122), plt.imshow(magnitude_spectrum,cmap='gray')
plt.title("Magnitude Spectrum"), plt.axis('off')
plt.show()
결과는 아래와 같다.
가운데의 흰색 영역을 볼 수 있는데, 이는 저주파수 함량이 더 많다는 것을 말한다.
주파수 변환을 찾았으니 이제 고역 통과 필터링, 이미지 재구성, 즉 역(inverse) DFT와 같은 주파수 영역의 기능을 할 수 있다. 이를 위해 60x60 크기의 직사각형 창으로 마스킹해서 낮은 주파수를 제거하기만 하면 된다. 그런 다음 np.fft.ifftshift() 함수를 사용하여 역방향 이동을 적용하여 DC 성분이 다시 왼쪽 상단에 오도록 한다. 그런 다음 np.ifft2() 함수를 사용하여 역 FFT를 찾는다. 결과는 다시 복잡한 숫자가 될 것이다. 이제 절대값을 취하면 된다.
rows,cols = img.shape
# 소수점으로 떨어지는 것을 방지하기 위함 : round
crow,ccol = round(rows/2), round(cols/2)
# 60x60 크기로 창 만들기
fshift[crow-30:crow+60, ccol-30:ccol+30] = 0
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)
plt.figure(figsize=(12,8))
plt.subplot(131),plt.imshow(img,cmap='gray')
plt.title("Input Image"), plt.axis('off')
plt.subplot(132),plt.imshow(img_back,cmap='gray')
plt.title("Image after HPF"), plt.axis('off')
plt.subplot(133),plt.imshow(img_back)
plt.title("Result in JET"), plt.axis('off')
plt.show()
결과는 아래와 같다.
결과를 보면 HPF(High Pass Filtering)가 가장자리 검출 역할을 한다는 것을 알 수 있다. 이는 Image Gradient 장에서 봤었던 것이다. 이는 또한 대부분의 이미지 데이터가 스펙트럼의 저주파 영역에 존재한다는 것을 보여준다. 어쨌든 Numpy에서 DFT, IDFT 등을 찾는 법을 보았다. 이제 OpenCV에서 어떻게 하는지 보자!
결과를 자세히 봤다면, 마지막 결과인 JET 색버전에서 몇 가지 인공물(artifacts)를 볼 수 있다. 이는 구조가 거기 있듯이 잔물결을 보여준다. 그리고 이를 Ringing Effects라고 부른다. 이는 우리가 마스킹하는데 사용한 직사각형 창에 의해 발생한다. 이 마스크는 이 문제를 일으키는 sinc 모양으로 변환된다. 그래서 직사각형 창은 필터링에 적용되지 않는다. 더 나은 방법은 가우시안 창을 이용하는 것이다.
Fourier Transform in OpenCV
OpenCV는 이를 위한 cv2.dft()와 cv2.idft() 함수를 제공한다. 이는 이전과 같은 결과를 돌려주지만, 두 개의 채널을 가지고 있다. 처음 채널은 결과의 실제 부분을 담당하고, 두 번째 채널은 결과의 상상의 부분을 맡는다. 입력 이미지는 np.float32로 변환되어져야 한다. 어떻게 하는지 봐보자!
img = cv2.imread('./images/face2.jpg',0)
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]))
plt.subplot(121), plt.imshow(img,cmap='gray')
plt.title("Input Image"), plt.axis('off')
plt.subplot(122), plt.imshow(magnitude_spectrum,cmap='gray')
plt.title("Magnitude Spectrum"), plt.axis('off')
plt.show()
결과는 당연히 같다!
** 한 번에 크기와 위상을 반환하는 cv2.cartToPolar()를 사용할 수 있다.
이제 역 DFT를 구해야한다, 이전 파트에서는 HPF를 생성했고, 이번에는 이미지에서 고주파 부분을 제거하는 방법, 즉 영상에 LPF를 적용하는 방법을 살펴보자. 이는 사실 이미지를 흐리게 한다. 이를 위해 먼저 저주파에서는 높은 값(1)으로 마스크를 생성한다. 즉 LF를 통과하고 HF 지역에서는 0을 통과한다.
rows, cols = img.shape
crow,ccol = round(rows/2), round(cols/2)
# 마스크를 먼저 생성하고, 가운데 네모를 1로 나머지를 0으로 한다.
mask = np.zeros((rows,cols,2),np.uint8)
mask[crow-30:crow+30,ccol-30:ccol+30] = 1
# 마스크와 역 DFT를 적용한다.
fshift = dft_shift*mask
f_ishift = np.fft.ifftshift(fshift)
img_back = cv2.idft(f_ishift)
img_back = cv2.magnitude(img_back[:,:,0],img_back[:,:,1])
plt.figure(figsize=(12,8))
plt.subplot(121), plt.imshow(img,cmap='gray')
plt.title("Input Image"), plt.axis('off')
plt.subplot(122), plt.imshow(img_back,cmap='gray')
plt.title("Magnitude Spectrum"), plt.axis('off')
plt.show()
결과를 보자.
** 보통, OpenCV의 cv2.dft()와 cv2.idft() 함수는 Numpy 보다 빠르다. 하지만 Numpy 함수는 좀 더 유저친화적이다. 성능에 관한 세부적인 내용은 아래 섹션을 보자.
Performance Optimization of DFT
DFT 연산의 성능은 일부 배열 크기에서 더 좋다. 이는 배열의 크기가 2배 일 때 더 빠르다. 배열의 크기가 2배, 3배, 5배인 경우도 효율적이다. 코드의 성능에 대해 걱정이라면, 배열의 크기를 DFT를 찾기 이전에 어떠한 최적의 크기(제로-패딩으로)로 수정하면 된다. OpenCV에서는, 직접 제로패딩을 해야한다. 하지만 Numpy에서는, FFT 계산의 새 크기를 지정하면 자동으로 0이 패딩된다.
그러면 최적의 크기를 어떻게 찾을 수 있을까? OpenCV는 이를 위한 함수인, cv2.getOptimalDFTSize()를 제공한다. 이는 cv2.dft()와 np.fft.fft2() 모두에 적용된다. IPython의 마법같은 명령어인 %timeit을 사용해서 성능을 확인해보자.
img = cv2.imread('./images/face2.jpg',0)
rows,cols = img.shape
print(rows,cols)
# 501 398
nrows = cv2.getOptimalDFTSize(rows)
ncols = cv2.getOptimalDFTSize(cols)
print(nrows,ncols)
# 512 400
보면, 크기 (501, 398)이 (512, 400)으로 수정된다. 이제 제로 패딩을 해보고 DFT 연산 성능을 찾아보자. 커다란 0 배열을 만들어도 되고 데이터를 복사하고 cv2.copyMakeBorder()를 사용해도 된다.
nimg = np.zeros((nrows,ncols))
nimg[:rows,:cols] = img
또는
right = ncols - cols
bottom = nrows - rows
bordertype = cv2.BORDER_CONSTANT
nimg = cv2.copyMakeBorder(img,0,bottom,0,right,bordertype, value=0)
이제 Numpy 함수의 DFT 성능을 비교해보자.
>>>%timeit fft1 = np.fft.fft2(img)
71.4 ms ± 1.62 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
>>>%timeit fft2 = np.fft.fft2(img,[nrows,ncols])
8.1 ms ± 330 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
이는 4배 빠른 것을 보여준다. 이를 OpenCV 함수로도 해볼 것이다.
>>>%timeit dft1= cv2.dft(np.float32(img),flags=cv2.DFT_COMPLEX_OUTPUT)
2.13 ms ± 101 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
>>>%timeit dft2= cv2.dft(np.float32(nimg),flags=cv2.DFT_COMPLEX_OUTPUT)
1.16 ms ± 69.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
이 또한 4배 빠른 것을 보여준다. OpenCV 함수가 Numpy 함수보다 3배정도 빠른 것을 볼 수 있다. 이는 역 FFT로 해도 비슷하다.
Why Laplacian is a High Pass Filter?
포럼에서도 비슷한 질문이 있었다. 그 질문은, 왜 Laplacian이 HPF인가? 왜 Sobel이 HPF인가? 이었다. 그리고 그것에 대한 첫 번째 대답은 Fourier Transform에 관한 것이다. 좀 더 큰 크기의 FFT를 위해 Laplacian의 Fourier Transform을 해봐라!
# 조정 파라미터 없는 간단한 평균 필터
mean_filter = np.ones((3,3))
# 가우시안 필터 만들기
x = cv2.getGaussianKernel(3,3)
gaussian = x*x.T
## 다른 가장자리 검출 필터들
# x 방향으로의 Scharr
scharr = np.array([[-3,0,3],
[-10,0,10],
[-3,0,3]])
# x 방향으로의 Sobel
sobel_x = np.array([[-1,0,1],
[-2,0,2],
[-1,0,1]])
# y 방향으로의 Sobel
sobel_y = np.array([[-1,-2,-1],
[0,0,0],
[1,2,1]])
# laplacian
laplacian = np.array([[0,1,0],
[1,-4,1],
[0,1,0]])
filters = [mean_filter, gaussian, laplacian, sobel_x, sobel_y, scharr]
filter_name = ['mean_filter', 'gaussian','laplacian', 'sobel_x', \
'sobel_y', 'scharr_x']
fft_filters = [np.fft.fft2(x) for x in filters]
fft_shift = [np.fft.fftshift(y) for y in fft_filters]
mag_spectrum = [np.log(np.abs(z)+1) for z in fft_shift]
for i in range(6):
plt.subplot(2,3,i+1), plt.imshow(mag_spectrum[i],cmap='gray')
plt.title(filter_name[i]), plt.axis('off')
plt.show()
아래 결과를 보자.
이미지에서, 각 커널이 차단하는 주파수 영역과 통과되는 영역을 볼 수 있다. 이 정보로부터, 우리는 왜 각 커널이 HPF인지 LPF인지 말할 수 있다.