← Back to branch summary

~paparazzi-uav/paparazzi/v5.0-manual

« back to all changes in this revision

Viewing changes to sw/ext/opencv_bebop/opencv/doc/py_tutorials/py_imgproc/py_transforms/py_fourier_transform/py_fourier_transform.markdown

  • Committer: Paparazzi buildbot
  • Date: 2016-05-18 15:00:29 UTC
  • Revision ID: felix.ruess+docbot@gmail.com-20160518150029-e8lgzi5kvb4p7un9
Manual import commit 4b8bbb730080dac23cf816b98908dacfabe2a8ec from v5.0 branch.

Show diffs side-by-side

added added

removed removed

Lines of Context:
sw/ext/opencv_bebop/opencv/doc/py_tutorials/py_imgproc/py_transforms/py_fourier_transform/py_fourier_transform.markdown
 
1
Fourier Transform {#tutorial_py_fourier_transform}
 
2
=================
 
3
 
 
4
Goal
 
5
----
 
6
 
 
7
In this section, we will learn
 
8
    -   To find the Fourier Transform of images using OpenCV
 
9
    -   To utilize the FFT functions available in Numpy
 
10
    -   Some applications of Fourier Transform
 
11
    -   We will see following functions : **cv2.dft()**, **cv2.idft()** etc
 
12
 
 
13
Theory
 
14
------
 
15
 
 
16
Fourier Transform is used to analyze the frequency characteristics of various filters. For images,
 
17
**2D Discrete Fourier Transform (DFT)** is used to find the frequency domain. A fast algorithm
 
18
called **Fast Fourier Transform (FFT)** is used for calculation of DFT. Details about these can be
 
19
found in any image processing or signal processing textbooks. Please see Additional Resources_
 
20
section.
 
21
 
 
22
For a sinusoidal signal, \f$x(t) = A \sin(2 \pi ft)\f$, we can say \f$f\f$ is the frequency of signal, and
 
23
if its frequency domain is taken, we can see a spike at \f$f\f$. If signal is sampled to form a discrete
 
24
signal, we get the same frequency domain, but is periodic in the range \f$[- \pi, \pi]\f$ or \f$[0,2\pi]\f$
 
25
(or \f$[0,N]\f$ for N-point DFT). You can consider an image as a signal which is sampled in two
 
26
directions. So taking fourier transform in both X and Y directions gives you the frequency
 
27
representation of image.
 
28
 
 
29
More intuitively, for the sinusoidal signal, if the amplitude varies so fast in short time, you can
 
30
say it is a high frequency signal. If it varies slowly, it is a low frequency signal. You can extend
 
31
the same idea to images. Where does the amplitude varies drastically in images ? At the edge points,
 
32
or noises. So we can say, edges and noises are high frequency contents in an image. If there is no
 
33
much changes in amplitude, it is a low frequency component. ( Some links are added to
 
34
Additional Resources_ which explains frequency transform intuitively with examples).
 
35
 
 
36
Now we will see how to find the Fourier Transform.
 
37
 
 
38
Fourier Transform in Numpy
 
39
--------------------------
 
40
 
 
41
First we will see how to find Fourier Transform using Numpy. Numpy has an FFT package to do this.
 
42
**np.fft.fft2()** provides us the frequency transform which will be a complex array. Its first
 
43
argument is the input image, which is grayscale. Second argument is optional which decides the size
 
44
of output array. If it is greater than size of input image, input image is padded with zeros before
 
45
calculation of FFT. If it is less than input image, input image will be cropped. If no arguments
 
46
passed, Output array size will be same as input.
 
47
 
 
48
Now once you got the result, zero frequency component (DC component) will be at top left corner. If
 
49
you want to bring it to center, you need to shift the result by \f$\frac{N}{2}\f$ in both the
 
50
directions. This is simply done by the function, **np.fft.fftshift()**. (It is more easier to
 
51
analyze). Once you found the frequency transform, you can find the magnitude spectrum.
 
52
@code{.py}
 
53
import cv2
 
54
import numpy as np
 
55
from matplotlib import pyplot as plt
 
56
 
 
57
img = cv2.imread('messi5.jpg',0)
 
58
f = np.fft.fft2(img)
 
59
fshift = np.fft.fftshift(f)
 
60
magnitude_spectrum = 20*np.log(np.abs(fshift))
 
61
 
 
62
plt.subplot(121),plt.imshow(img, cmap = 'gray')
 
63
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
 
64
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
 
65
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
 
66
plt.show()
 
67
@endcode
 
68
Result look like below:
 
69
 
 
70
![image](images/fft1.jpg)
 
71
 
 
72
See, You can see more whiter region at the center showing low frequency content is more.
 
73
 
 
74
So you found the frequency transform Now you can do some operations in frequency domain, like high
 
75
pass filtering and reconstruct the image, ie find inverse DFT. For that you simply remove the low
 
76
frequencies by masking with a rectangular window of size 60x60. Then apply the inverse shift using
 
77
**np.fft.ifftshift()** so that DC component again come at the top-left corner. Then find inverse FFT
 
78
using **np.ifft2()** function. The result, again, will be a complex number. You can take its
 
79
absolute value.
 
80
@code{.py}
 
81
rows, cols = img.shape
 
82
crow,ccol = rows/2 , cols/2
 
83
fshift[crow-30:crow+30, ccol-30:ccol+30] = 0
 
84
f_ishift = np.fft.ifftshift(fshift)
 
85
img_back = np.fft.ifft2(f_ishift)
 
86
img_back = np.abs(img_back)
 
87
 
 
88
plt.subplot(131),plt.imshow(img, cmap = 'gray')
 
89
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
 
90
plt.subplot(132),plt.imshow(img_back, cmap = 'gray')
 
91
plt.title('Image after HPF'), plt.xticks([]), plt.yticks([])
 
92
plt.subplot(133),plt.imshow(img_back)
 
93
plt.title('Result in JET'), plt.xticks([]), plt.yticks([])
 
94
 
 
95
plt.show()
 
96
@endcode
 
97
Result look like below:
 
98
 
 
99
![image](images/fft2.jpg)
 
100
 
 
101
The result shows High Pass Filtering is an edge detection operation. This is what we have seen in
 
102
Image Gradients chapter. This also shows that most of the image data is present in the Low frequency
 
103
region of the spectrum. Anyway we have seen how to find DFT, IDFT etc in Numpy. Now let's see how to
 
104
do it in OpenCV.
 
105
 
 
106
If you closely watch the result, especially the last image in JET color, you can see some artifacts
 
107
(One instance I have marked in red arrow). It shows some ripple like structures there, and it is
 
108
called **ringing effects**. It is caused by the rectangular window we used for masking. This mask is
 
109
converted to sinc shape which causes this problem. So rectangular windows is not used for filtering.
 
110
Better option is Gaussian Windows.
 
111
 
 
112
Fourier Transform in OpenCV
 
113
---------------------------
 
114
 
 
115
OpenCV provides the functions **cv2.dft()** and **cv2.idft()** for this. It returns the same result
 
116
as previous, but with two channels. First channel will have the real part of the result and second
 
117
channel will have the imaginary part of the result. The input image should be converted to
 
118
np.float32 first. We will see how to do it.
 
119
@code{.py}
 
120
import numpy as np
 
121
import cv2
 
122
from matplotlib import pyplot as plt
 
123
 
 
124
img = cv2.imread('messi5.jpg',0)
 
125
 
 
126
dft = cv2.dft(np.float32(img),flags = cv2.DFT_COMPLEX_OUTPUT)
 
127
dft_shift = np.fft.fftshift(dft)
 
128
 
 
129
magnitude_spectrum = 20*np.log(cv2.magnitude(dft_shift[:,:,0],dft_shift[:,:,1]))
 
130
 
 
131
plt.subplot(121),plt.imshow(img, cmap = 'gray')
 
132
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
 
133
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
 
134
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
 
135
plt.show()
 
136
@endcode
 
137
 
 
138
@note You can also use **cv2.cartToPolar()** which returns both magnitude and phase in a single shot
 
139
 
 
140
So, now we have to do inverse DFT. In previous session, we created a HPF, this time we will see how
 
141
to remove high frequency contents in the image, ie we apply LPF to image. It actually blurs the
 
142
image. For this, we create a mask first with high value (1) at low frequencies, ie we pass the LF
 
143
content, and 0 at HF region.
 
144
 
 
145
@code{.py}
 
146
rows, cols = img.shape
 
147
crow,ccol = rows/2 , cols/2
 
148
 
 
149
# create a mask first, center square is 1, remaining all zeros
 
150
mask = np.zeros((rows,cols,2),np.uint8)
 
151
mask[crow-30:crow+30, ccol-30:ccol+30] = 1
 
152
 
 
153
# apply mask and inverse DFT
 
154
fshift = dft_shift*mask
 
155
f_ishift = np.fft.ifftshift(fshift)
 
156
img_back = cv2.idft(f_ishift)
 
157
img_back = cv2.magnitude(img_back[:,:,0],img_back[:,:,1])
 
158
 
 
159
plt.subplot(121),plt.imshow(img, cmap = 'gray')
 
160
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
 
161
plt.subplot(122),plt.imshow(img_back, cmap = 'gray')
 
162
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
 
163
plt.show()
 
164
@endcode
 
165
See the result:
 
166
 
 
167
![image](images/fft4.jpg)
 
168
 
 
169
@note As usual, OpenCV functions **cv2.dft()** and **cv2.idft()** are faster than Numpy
 
170
counterparts. But Numpy functions are more user-friendly. For more details about performance issues,
 
171
see below section.
 
172
 
 
173
Performance Optimization of DFT
 
174
===============================
 
175
 
 
176
Performance of DFT calculation is better for some array size. It is fastest when array size is power
 
177
of two. The arrays whose size is a product of 2’s, 3’s, and 5’s are also processed quite
 
178
efficiently. So if you are worried about the performance of your code, you can modify the size of
 
179
the array to any optimal size (by padding zeros) before finding DFT. For OpenCV, you have to
 
180
manually pad zeros. But for Numpy, you specify the new size of FFT calculation, and it will
 
181
automatically pad zeros for you.
 
182
 
 
183
So how do we find this optimal size ? OpenCV provides a function, **cv2.getOptimalDFTSize()** for
 
184
this. It is applicable to both **cv2.dft()** and **np.fft.fft2()**. Let's check their performance
 
185
using IPython magic command %timeit.
 
186
@code{.py}
 
187
In [16]: img = cv2.imread('messi5.jpg',0)
 
188
In [17]: rows,cols = img.shape
 
189
In [18]: print rows,cols
 
190
342 548
 
191
 
 
192
In [19]: nrows = cv2.getOptimalDFTSize(rows)
 
193
In [20]: ncols = cv2.getOptimalDFTSize(cols)
 
194
In [21]: print nrows, ncols
 
195
360 576
 
196
@endcode
 
197
See, the size (342,548) is modified to (360, 576). Now let's pad it with zeros (for OpenCV) and find
 
198
their DFT calculation performance. You can do it by creating a new big zero array and copy the data
 
199
to it, or use **cv2.copyMakeBorder()**.
 
200
@code{.py}
 
201
nimg = np.zeros((nrows,ncols))
 
202
nimg[:rows,:cols] = img
 
203
@endcode
 
204
OR:
 
205
@code{.py}
 
206
right = ncols - cols
 
207
bottom = nrows - rows
 
208
bordertype = cv2.BORDER_CONSTANT #just to avoid line breakup in PDF file
 
209
nimg = cv2.copyMakeBorder(img,0,bottom,0,right,bordertype, value = 0)
 
210
@endcode
 
211
Now we calculate the DFT performance comparison of Numpy function:
 
212
@code{.py}
 
213
In [22]: %timeit fft1 = np.fft.fft2(img)
 
214
10 loops, best of 3: 40.9 ms per loop
 
215
In [23]: %timeit fft2 = np.fft.fft2(img,[nrows,ncols])
 
216
100 loops, best of 3: 10.4 ms per loop
 
217
@endcode
 
218
It shows a 4x speedup. Now we will try the same with OpenCV functions.
 
219
@code{.py}
 
220
In [24]: %timeit dft1= cv2.dft(np.float32(img),flags=cv2.DFT_COMPLEX_OUTPUT)
 
221
100 loops, best of 3: 13.5 ms per loop
 
222
In [27]: %timeit dft2= cv2.dft(np.float32(nimg),flags=cv2.DFT_COMPLEX_OUTPUT)
 
223
100 loops, best of 3: 3.11 ms per loop
 
224
@endcode
 
225
It also shows a 4x speed-up. You can also see that OpenCV functions are around 3x faster than Numpy
 
226
functions. This can be tested for inverse FFT also, and that is left as an exercise for you.
 
227
 
 
228
Why Laplacian is a High Pass Filter?
 
229
------------------------------------
 
230
 
 
231
A similar question was asked in a forum. The question is, why Laplacian is a high pass filter? Why
 
232
Sobel is a HPF? etc. And the first answer given to it was in terms of Fourier Transform. Just take
 
233
the fourier transform of Laplacian for some higher size of FFT. Analyze it:
 
234
@code{.py}
 
235
import cv2
 
236
import numpy as np
 
237
from matplotlib import pyplot as plt
 
238
 
 
239
# simple averaging filter without scaling parameter
 
240
mean_filter = np.ones((3,3))
 
241
 
 
242
# creating a guassian filter
 
243
x = cv2.getGaussianKernel(5,10)
 
244
gaussian = x*x.T
 
245
 
 
246
# different edge detecting filters
 
247
# scharr in x-direction
 
248
scharr = np.array([[-3, 0, 3],
 
249
                   [-10,0,10],
 
250
                   [-3, 0, 3]])
 
251
# sobel in x direction
 
252
sobel_x= np.array([[-1, 0, 1],
 
253
                   [-2, 0, 2],
 
254
                   [-1, 0, 1]])
 
255
# sobel in y direction
 
256
sobel_y= np.array([[-1,-2,-1],
 
257
                   [0, 0, 0],
 
258
                   [1, 2, 1]])
 
259
# laplacian
 
260
laplacian=np.array([[0, 1, 0],
 
261
                    [1,-4, 1],
 
262
                    [0, 1, 0]])
 
263
 
 
264
filters = [mean_filter, gaussian, laplacian, sobel_x, sobel_y, scharr]
 
265
filter_name = ['mean_filter', 'gaussian','laplacian', 'sobel_x', \
 
266
                'sobel_y', 'scharr_x']
 
267
fft_filters = [np.fft.fft2(x) for x in filters]
 
268
fft_shift = [np.fft.fftshift(y) for y in fft_filters]
 
269
mag_spectrum = [np.log(np.abs(z)+1) for z in fft_shift]
 
270
 
 
271
for i in xrange(6):
 
272
    plt.subplot(2,3,i+1),plt.imshow(mag_spectrum[i],cmap = 'gray')
 
273
    plt.title(filter_name[i]), plt.xticks([]), plt.yticks([])
 
274
 
 
275
plt.show()
 
276
@endcode
 
277
See the result:
 
278
 
 
279
![image](images/fft5.jpg)
 
280
 
 
281
From image, you can see what frequency region each kernel blocks, and what region it passes. From
 
282
that information, we can say why each kernel is a HPF or a LPF
 
283
 
 
284
Additional Resources
 
285
--------------------
 
286
 
 
287
-#  [An Intuitive Explanation of Fourier
 
288
    Theory](http://cns-alumni.bu.edu/~slehar/fourier/fourier.html) by Steven Lehar
 
289
2.  [Fourier Transform](http://homepages.inf.ed.ac.uk/rbf/HIPR2/fourier.htm) at HIPR
 
290
3.  [What does frequency domain denote in case of images?](http://dsp.stackexchange.com/q/1637/818)
 
291
 
 
292
Exercises
 
293
---------