#### Transcript 5_frequency_domain

```CS 691 Computational
Photography
Instructor: Gianfranco Doretto
Frequency Domain
Overview
• Frequency domain analysis
• Sampling and reconstruction
Linear image transformations
• In analyzing images, it’s often useful to
make a change of basis.


F  Uf
transformed image
Vectorized image
Fourier transform, or
Wavelet transform, or
Steerable pyramid transform
Self-inverting transforms
Same basis functions are used for the inverse transform


1
f U F


U F
U transpose and complex conjugate
“Gala Contemplating the Mediterranean Sea,
which at 30 meters becomes the portrait
of Abraham Lincoln”, 1976
A nice set of basis
Teases away fast vs. slow changes in the image.
This change of basis has a special name…
Jean Baptiste Joseph Fourier (1768-1830)
manner in which the author arrives at these
equations is not exempt of difficulties and...his
Any univariate function can analysis to integrate them still leaves something to be
be rewritten as a weighted desired on the score of generality and even rigour.
sum of sines and cosines of
different frequencies.
• Don’t believe it?
– Neither did
Lagrange, Laplace,
Poisson and other
big wigs
– Not translated into
English until 1878!
• But it’s (mostly)
true!
– called Fourier Series
– there are some
Laplace
Lagrange
Legendre
A sum of sines
Our building block:
Asin(x   
Add enough of them to get
any signal f(x) you want!
Frequency Spectra
• example : g(t) = sin(2πf t) + (1/3)sin(2π(3f) t)
=
+
Frequency Spectra
Frequency Spectra
=
=
+
Frequency Spectra
=
=
+
Frequency Spectra
=
=
+
Frequency Spectra
=
=
+
Frequency Spectra
=
=
+
Frequency Spectra

1
= A sin(2 kt )
k 1 k
Fourier basis for image analysis
Intensity Image
Frequency Spectra
http://sharp.bu.edu/~slehar/fourier/fourier.html#filtering
Signals can be composed
+
=
http://sharp.bu.edu/~slehar/fourier/fourier.html#filtering
More: http://www.cs.unm.edu/~brayer/vision/fourier.html
Fourier Transform
• Fourier transform stores the magnitude and phase at
each frequency
– Magnitude encodes how much signal there is at a particular
frequency
– Phase encodes spatial information (indirectly)
– For mathematical convenience, this is often notated in terms of
real and complex numbers
Amplitude:
A   R( )  I ( )
Euler’s formula:
2
2
I ( )
Phase:   tan
R( )
1
Computing the Fourier Transform
Continuous
Discrete
k = -N/2..N/2
Fast Fourier Transform (FFT): NlogN
Fourier Transform pairs
Phase and Magnitude
• Fourier transform of a
real function is complex
– difficult to plot, visualize
– instead, we can think of the
phase and magnitude of
the transform
• Phase is the phase of the
complex transform
• Magnitude is the
magnitude of the complex
transform
• Curious fact
– all natural images have
transform
– hence, phase seems to
matter, but magnitude
largely doesn’t
• Demonstration
– Take two pictures, swap
the phase transforms,
compute the inverse - what
does the result look like?
This is the
magnitude
transform
of the
cheetah
pic
This is the
phase
transform
of the
cheetah
pic
This is the
magnitude
transform
of the
zebra pic
This is the
phase
transform
of the
zebra pic
Reconstructio
n with zebra
phase,
cheetah
magnitude
Reconstruction
with cheetah
phase, zebra
magnitude
The Convolution Theorem
• The Fourier transform of the convolution of two
functions is the product of their Fourier
transforms
• The inverse Fourier transform of the product of
two Fourier transforms is the convolution of
the two inverse Fourier transforms
• Convolution in spatial domain is equivalent to
multiplication in frequency domain!
Properties of Fourier Transforms
• Linearity
• Fourier transform of a real signal is
• The energy of the signal is the same as
the energy of its Fourier transform
See Szeliski Book (3.4)
Filtering in spatial domain
1
0
-1
2
0
-2
1
0
-1
*
=
Filtering in frequency domain
FFT
FFT
=
Inverse FFT
Play with FFT in Matlab
• Filtering with fft
im = ... % “im” should be a gray-scale floating point image
[imh, imw] = size(im);
fftsize = 1024; % should be order of 2 (for speed) and include padding
im_fft = fft2(im, fftsize, fftsize); % 1) fft im with padding
hs = 50; % filter half-size
fil = fspecial('gaussian', hs*2+1, 10);
fil_fft = fft2(fil, fftsize, fftsize); % 2) fft fil, pad to same size as image
im_fil_fft = im_fft .* fil_fft; % 3) multiply fft images
im_fil = ifft2(im_fil_fft); % 4) inverse fft2
im_fil = im_fil(1+hs:size(im,1)+hs, 1+hs:size(im, 2)+hs); % 5) remove padding
• Displaying with fft
figure(1), imagesc(log(abs(fftshift(im_fft)))), axis image, colormap jet
Convolution versus FFT
• 1-d FFT: O(NlogN) computation time,
where N is number of samples.
• 2-d FFT: 2N(NlogN), where N is number of
pixels on a side
• Convolution: K N2, where K is number of
samples in kernel
• Say N=210, K=100. 2-d FFT: 20 220, while
convolution gives 100 220
Why is the Fourier domain
particularly useful?
• It tells us the effect of linear convolutions.
• There is a fast algorithm for performing the
DFT, allowing for efficient signal filtering.
• The Fourier domain offers an alternative
domain for understanding and
manipulating the image.
coefficient
Analysis of our simple filters
1.0
0
Pixel offset
original
M 1
F [m]   f [k ]e
Filtered
(no change)
 km 
i 

M


k 0
1
1.0
0
constant
coefficient
Analysis of our simple filters

1.0
0
Pixel offset
original
shifted
M 1
F [m]   f [k ]e
k 0
e
i
m
M
 km 
i 

M 
Constant
magnitude,
1.0 linearly shifted
phase
0
coefficient
Analysis of our simple filters
0.3
0
Pixel offset
original
blurred
M 1
F [m]   f [k ]e
 km 
i 

M


k 0
1
 m  
 1  2 cos
 
3
 M 
1.0
0
Low-pass
filter
Analysis of our simple filters
2.0
0.33
0
0
sharpened
original
M 1
F [m]   f [k ]e
 km 
i 

M


high-pass filter
2.3
k 0
1
 m  
 2  1  2 cos
 
3
 M 
1.0
0
Playing with the DFT of an image
Can change spectrum, then reconstruct
Low and High Pass filtering
Why does the Gaussian give a nice smooth image, but
the square filter give edgy artifacts?
Gaussian
Box filter
Gaussian Filter
Box Filter
Overview
• Frequency domain analysis
• Sampling and reconstruction
Sampling and Reconstruction
Sampled representations
• How to store and compute with continuous
functions?
• Common scheme for representation: samples
[FvDFH fig.14.14b / Wolberg]
– write down the function’s values at many points
Reconstruction
• Making samples back into a continuous function
[FvDFH fig.14.14b / Wolberg]
– for output (need realizable method)
– for analysis or processing (need mathematical method)
– amounts to “guessing” what the function did in between
1D Example: Audio
low
high
frequencies
Sampling in digital audio
• Recording: sound to analog to samples to disc
• Playback: disc to samples to analog to sound
again
– how can we be sure we are filling in the gaps
correctly?
Sampling and Reconstruction
• Simple example: a sign wave
Undersampling
• What if we “missed” things between the
samples?
• Simple example: undersampling a sine wave
– unsurprising result: information is lost
Undersampling
• What if we “missed” things between the
samples?
• Simple example: undersampling a sine wave
– unsurprising result: information is lost
– surprising result: indistinguishable from lower
frequency
Undersampling
• What if we “missed” things between the
samples?
• Simple example: undersampling a sine wave
– unsurprising result: information is lost
– surprising result: indistinguishable from lower
frequency
– also was always indistinguishable from higher
frequencies
– aliasing: signals “traveling in disguise” as other
frequencies
Aliasing in video
Aliasing in images
What’s happening?
Input signal:
Plot as image:
x = 0:.05:5; imagesc(sin((2.^x).*x))
Alias!
Not enough samples
Fourier Transform of a Sampled Signal
Sampling period T
Sampling frequency F = 1/T
T 2T 3T …
-F
0
F
…
Fourier Transform of a Sampled Signal
Sampling period T
Sampling frequency F = 1/T
T
2T
3T …
-2F
-F
0
F
2F
…
Nyquist-Shannon Sampling Theorem
• When sampling a signal at discrete intervals, the sampling
frequency must be F > 2 x fmax
• fmax = max frequency of the input signal
• This will allow to reconstruct the original perfectly from the
sampled version
• Do not know fmax or cannot sample at that rate? Prefiltering! We
lose something but still better than aliasing!!!!!
v
v
v
good
Image halfsizing
This image is too big to
fit on the screen. How
can we reduce it?
How to generate a halfsized version?
Image sub-sampling
1/8
1/4
Throw away every other row and
column to create a 1/2 size image
- called image sub-sampling
Image sub-sampling
1/2
1/4
Aliasing! What do we do?
(2x zoom)
1/8
(4x zoom)
Gaussian (lowpass) pre-filtering
G 1/8
G 1/4
Gaussian 1/2
Solution: filter the image, then subsample
• Filter size should double for each ½ size reduction. Why?
Subsampling with Gaussian pre-filtering
Gaussian 1/2
G 1/4
G 1/8
Compare with...
1/2
1/4
(2x zoom)
1/8
(4x zoom)
Algorithm for downsampling by factor of 2
2. Apply low-pass filter
im_blur = imfilter(image, fspecial(‘gaussian’, 7, 1))
3. Sample every other pixel
im_small = im_blur(1:2:end, 1:2:end);
Slide Credits
• This set of sides also contains
contributions
kindly made available by the following
authors
– Alexei Efros
– Frédo Durand
– Bill Freeman
– Steve Seitz
– Derek Hoiem
– Steve Marschner
```