پردازش تصویر در OpenCV

تبدیل فوریه (Fourier Transform)

۱۰ دقیقه مطالعهOpenCV 5.xبه‌روزرسانی: ۹ آبان ۱۴۰۵

هدف#

در این بخش، یاد خواهیم گرفت

چگونه تبدیل فوریه تصاویر را با استفاده از OpenCV بیابیم
چگونه از توابع FFT موجود در Numpy استفاده کنیم
برخی از کاربردهای تبدیل فوریه
توابع زیر را خواهیم دید: cv.dft(), cv.idft() و غیره

تئوری#

تبدیل فوریه (Fourier Transform) برای تحلیل ویژگی‌های فرکانسی (Frequency Characteristics) فیلترهای مختلف استفاده می‌شود. برای تصاویر، از تبدیل فوریه گسسته دوبعدی (2D Discrete Fourier Transform - DFT) برای یافتن حوزه فرکانس (Frequency Domain) استفاده می‌شود. یک الگوریتم سریع به نام تبدیل فوریه سریع (Fast Fourier Transform - FFT) برای محاسبه DFT استفاده می‌شود. جزئیات بیشتر درباره این موارد را می‌توانید در هر کتاب پردازش تصویر یا پردازش سیگنال بیابید. لطفاً به بخش منابع اضافی_ مراجعه کنید.

برای یک سیگنال سینوسی (Sinusoidal Signal)، $x (t) = A sin (2 π f t)$ ، می‌توانیم بگوییم $f$ فرکانس سیگنال است و اگر حوزه فرکانسی آن گرفته شود، یک پیک در $f$ مشاهده می‌کنیم. اگر سیگنال نمونه‌برداری شود تا یک سیگنال گسسته تشکیل دهد، همان حوزه فرکانسی را به دست می‌آوریم، اما در بازه $[- π, π]$ یا $[0, 2 π]$ (یا $[0, N]$ برای DFT با N نقطه) تناوبی است. می‌توانید یک تصویر را به عنوان یک سیگنال که در دو جهت نمونه‌برداری شده در نظر بگیرید. بنابراین گرفتن تبدیل فوریه در هر دو جهت X و Y نمایش فرکانسی تصویر را به شما می‌دهد.

به طور شهودی‌تر، برای سیگنال سینوسی، اگر دامنه (Amplitude) خیلی سریع در زمان کوتاه تغییر کند، می‌توان گفت این یک سیگنال با فرکانس بالا (High Frequency) است. اگر به آرامی تغییر کند، یک سیگنال با فرکانس پایین (Low Frequency) است. می‌توانید همین ایده را به تصاویر تعمیم دهید. دامنه در کجای تصاویر به شدت تغییر می‌کند؟ در نقاط لبه (Edge) یا نویزها (Noises). بنابراین می‌توانیم بگوییم، لبه‌ها و نویزها محتوای فرکانس بالا در یک تصویر هستند. اگر تغییر زیادی در دامنه وجود نداشته باشد، مؤلفه فرکانس پایین (Low Frequency Component) است. (برخی پیوندها به منابع اضافی_ اضافه شده که تبدیل فرکانس را به صورت شهودی با مثال توضیح می‌دهد.)

اکنون خواهیم دید که چگونه تبدیل فوریه را بیابیم.

تبدیل فوریه در Numpy#

ابتدا خواهیم دید که چگونه تبدیل فوریه را با استفاده از Numpy بیابیم. Numpy یک بسته FFT برای انجام این کار دارد. np.fft.fft2() تبدیل فرکانسی را به ما می‌دهد که یک آرایه مختلط (Complex Array) خواهد بود. اولین آرگومان آن تصویر ورودی است که باید در مقیاس خاکستری باشد. آرگومان دوم اختیاری است که اندازه آرایه خروجی را تعیین می‌کند. اگر از اندازه تصویر ورودی بزرگتر باشد، تصویر ورودی قبل از محاسبه FFT با صفر پُر (Padded) می‌شود. اگر کوچکتر از تصویر ورودی باشد، تصویر ورودی برش می‌خورد. اگر هیچ آرگومانی ارسال نشود، اندازه آرایه خروجی با ورودی یکسان خواهد بود.

حال وقتی نتیجه را گرفتید، مؤلفه فرکانس صفر (DC Component) در گوشه بالا چپ قرار خواهد داشت. اگر می‌خواهید آن را به مرکز منتقل کنید، باید نتیجه را به اندازه $\frac{N}{2}$ در هر دو جهت جابجا کنید. این کار به سادگی با تابع np.fft.fftshift() انجام می‌شود. (تحلیل آسان‌تر است.) وقتی تبدیل فرکانسی را یافتید، می‌توانید طیف قدر (Magnitude Spectrum) را بیابید.

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt

img = cv.imread('messi5.jpg', cv.IMREAD_GRAYSCALE)
assert img is not None, "file could not be read, check with os.path.exists()"
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.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

نتیجه به صورت زیر است:

تبدیل فوریه (Fourier Transform) — تصویر ۱

می‌بینید، می‌توانید ناحیه سفیدتری را در مرکز مشاهده کنید که نشان می‌دهد محتوای فرکانس پایین بیشتر است.

بنابراین تبدیل فرکانسی را یافتید. اکنون می‌توانید برخی عملیات را در حوزه فرکانس انجام دهید، مانند فیلتر کردن فرکانس بالا (High Pass Filtering) و بازسازی تصویر، یعنی یافتن IDFT معکوس. برای این کار کافی است فرکانس‌های پایین را با ماسک کردن با یک پنجره مستطیلی به اندازه 60x60 حذف کنید. سپس جابجایی معکوس را با np.fft.ifftshift() اعمال کنید تا مؤلفه DC دوباره در گوشه بالا چپ قرار گیرد. سپس FFT معکوس را با استفاده از تابع np.ifft2() بیابید. نتیجه دوباره یک عدد مختلط خواهد بود. می‌توانید قدر مطلق آن را بگیرید.

rows, cols = img.shape
crow, ccol = rows//2, cols//2
fshift[crow-30:crow+31, ccol-30:ccol+31] = 0
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.real(img_back)

plt.subplot(131),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(img_back, cmap = 'gray')
plt.title('Image after HPF'), plt.xticks([]), plt.yticks([])
plt.subplot(133),plt.imshow(img_back)
plt.title('Result in JET'), plt.xticks([]), plt.yticks([])

plt.show()

نتیجه به صورت زیر است:

تبدیل فوریه (Fourier Transform) — تصویر ۲

نتیجه نشان می‌دهد که فیلترگذاری فرکانس بالا (High Pass Filtering) یک عملیات تشخیص لبه است. این همان چیزی است که در فصل گرادیان‌های تصویری دیدیم. همچنین نشان می‌دهد که بیشتر داده‌های تصویر در ناحیه فرکانس پایین طیف وجود دارند. به هر حال دیدیم که چگونه DFT، IDFT و غیره را در Numpy بیابیم. حالا بیایید ببینیم چگونه این کار را در OpenCV انجام دهیم.

اگر به دقت به نتیجه نگاه کنید، به خصوص آخرین تصویر در رنگ JET، می‌توانید برخی مصنوعات (Artifacts) را مشاهده کنید (یک نمونه با پیکان قرمز مشخص کرده‌ام). این نشان‌دهنده ساختارهایی شبیه امواج در آنجا است، که به آن اثر زنگوله‌ای (Ringing Effects) می‌گویند. این مشکل ناشی از پنجره مستطیلی است که برای ماسک کردن استفاده کردیم. این ماسک به شکل sinc تبدیل می‌شود که این مشکل را ایجاد می‌کند. بنابراین از پنجره مستطیلی برای فیلترگذاری استفاده نمی‌شود. گزینه بهتر، پنجره‌های گاوسی (Gaussian Windows) هستند.

تبدیل فوریه در OpenCV#

OpenCV توابع cv.dft() و cv.idft() را برای این منظور فراهم می‌کند. نتیجه مشابه قبلی را برمی‌گرداند، اما با دو کانال. کانال اول بخش حقیقی نتیجه و کانال دوم بخش موهومی نتیجه را خواهد داشت. تصویر ورودی باید ابتدا به np.float32 تبدیل شود. خواهیم دید که چگونه این کار را انجام دهیم.

import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt

img = cv.imread('messi5.jpg', cv.IMREAD_GRAYSCALE)
assert img is not None, "file could not be read, check with os.path.exists()"

dft = cv.dft(np.float32(img),flags = cv.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)

magnitude_spectrum = 20*np.log(cv.magnitude(dft_shift[:,:,0],dft_shift[:,:,1]))

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

بنابراین، اکنون باید DFT معکوس انجام دهیم. در جلسه قبلی یک HPF ایجاد کردیم، این بار خواهیم دید که چگونه محتوای فرکانس بالا را از تصویر حذف کنیم، یعنی یک LPF (فیلتر فرکانس پایین) به تصویر اعمال کنیم. این عملاً تصویر را تار می‌کند. برای این کار، ابتدا یک ماسک با مقدار بالا (۱) در فرکانس‌های پایین ایجاد می‌کنیم، یعنی محتوای LF را عبور می‌دهیم و ۰ در ناحیه HF.

rows, cols = img.shape
crow, ccol = rows//2, cols//2

# create a mask first, center square is 1, remaining all zeros
mask = np.zeros((rows,cols,2),np.uint8)
mask[crow-30:crow+30, ccol-30:ccol+30] = 1

# apply mask and inverse DFT
fshift = dft_shift*mask
f_ishift = np.fft.ifftshift(fshift)
img_back = cv.idft(f_ishift)
img_back = cv.magnitude(img_back[:,:,0],img_back[:,:,1])

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(img_back, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

نتیجه را ببینید:

تبدیل فوریه (Fourier Transform) — تصویر ۳

بهینه‌سازی عملکرد DFT#

عملکرد محاسبه DFT برای برخی اندازه‌های آرایه بهتر است. سریع‌ترین حالت زمانی است که اندازه آرایه توان عدد ۲ باشد. آرایه‌هایی که اندازه آن‌ها حاصل ضرب ۲ها، ۳ها و ۵ها است نیز به طور کارآمدی پردازش می‌شوند. بنابراین اگر نگران عملکرد کد خود هستید، می‌توانید اندازه آرایه را به هر اندازه بهینه‌ای (با پُر کردن با صفر) قبل از یافتن DFT تغییر دهید. برای OpenCV باید به صورت دستی صفر اضافه کنید. اما برای Numpy، اندازه جدید محاسبه FFT را مشخص می‌کنید و به طور خودکار صفر اضافه می‌کند.

پس چگونه این اندازه بهینه را بیابیم؟ OpenCV یک تابع به نام cv.getOptimalDFTSize() برای این منظور فراهم می‌کند. این تابع برای هر دو cv.dft() و np.fft.fft2() قابل استفاده است. بیایید عملکرد آن‌ها را با استفاده از دستور جادویی IPython به نام %timeit بررسی کنیم.

In [15]: img = cv.imread('messi5.jpg', cv.IMREAD_GRAYSCALE)
In [16]: assert img is not None, "file could not be read, check with os.path.exists()"
In [17]: rows,cols = img.shape
In [18]: print("{} {}".format(rows,cols))
342 548

In [19]: nrows = cv.getOptimalDFTSize(rows)
In [20]: ncols = cv.getOptimalDFTSize(cols)
In [21]: print("{} {}".format(nrows,ncols))
360 576

می‌بینید، اندازه (342،548) به (360، 576) تغییر یافته است. حالا بیایید آن را با صفر پُر کنیم (برای OpenCV) و عملکرد محاسبه DFT آن‌ها را بیابیم. می‌توانید این کار را با ایجاد یک آرایه صفر بزرگ جدید و کپی کردن داده‌ها به آن انجام دهید، یا از cv.copyMakeBorder() استفاده کنید.

nimg = np.zeros((nrows,ncols))
nimg[:rows,:cols] = img

یا:

right = ncols - cols
bottom = nrows - rows
bordertype = cv.BORDER_CONSTANT #just to avoid line breakup in PDF file
nimg = cv.copyMakeBorder(img,0,bottom,0,right,bordertype, value = 0)

اکنون مقایسه عملکرد DFT تابع Numpy را محاسبه می‌کنیم:

In [22]: %timeit fft1 = np.fft.fft2(img)
10 loops, best of 3: 40.9 ms per loop
In [23]: %timeit fft2 = np.fft.fft2(img,[nrows,ncols])
100 loops, best of 3: 10.4 ms per loop

این نشان‌دهنده سرعت ۴ برابری است. اکنون همین کار را با توابع OpenCV امتحان می‌کنیم.

In [24]: %timeit dft1= cv.dft(np.float32(img),flags=cv.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 13.5 ms per loop
In [27]: %timeit dft2= cv.dft(np.float32(nimg),flags=cv.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 3.11 ms per loop

همچنین سرعت ۴ برابری را نشان می‌دهد. می‌توانید ببینید که توابع OpenCV تقریباً ۳ برابر سریع‌تر از توابع Numpy هستند. این را می‌توان برای FFT معکوس نیز آزمایش کرد، که به عنوان تمرین برای شما باقی می‌ماند.

چرا لاپلاسی یک فیلتر فرکانس بالا است؟#

یک سؤال مشابه در یک انجمن پرسیده شد. سؤال این بود که چرا لاپلاسی یک فیلتر فرکانس بالا (High Pass Filter - HPF) است؟ چرا Sobel یک HPF است؟ و غیره. اولین پاسخ داده شده به آن بر اساس تبدیل فوریه بود. کافی است تبدیل فوریه لاپلاسی را برای اندازه بزرگتری از FFT بگیرید. آن را تحلیل کنید:

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt

# simple averaging filter without scaling parameter
mean_filter = np.ones((3,3))

# creating a gaussian filter
x = cv.getGaussianKernel(5,10)
gaussian = x*x.T

# different edge detecting filters
# scharr in x-direction
scharr = np.array([[-3, 0, 3],
                   [-10,0,10],
                   [-3, 0, 3]])
# sobel in x direction
sobel_x= np.array([[-1, 0, 1],
                   [-2, 0, 2],
                   [-1, 0, 1]])
# sobel in y direction
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.xticks([]), plt.yticks([])

plt.show()

نتیجه را ببینید:

تبدیل فوریه (Fourier Transform) — تصویر ۴

از تصویر می‌توانید ببینید که هر کرنل (Kernel) کدام ناحیه فرکانسی را مسدود می‌کند و کدام ناحیه را عبور می‌دهد. از این اطلاعات می‌توانیم بگوییم چرا هر کرنل یک HPF یا یک LPF است.

منابع اضافی#

An Intuitive Explanation of Fourier Theory نوشته Steven Lehar
Fourier Transform در HIPR
What does frequency domain denote in case of images?