Transcript f(x)

Fourier
Transform in
Image
Processing
Review - Image as a Function
• We can think of an image as a function, f,
• f: R2  R
– f (x, y) gives the intensity at position (x, y)
– Realistically, we expect the image only to be defined over a
rectangle, with a finite range:
• f: [a,b]x[c,d]  [0,1]
• A color image is just three functions pasted together.
• We can write this as a “vector-valued” function:
 r ( x, y ) 


f ( x, y )  g ( x, y )


 b( x, y ) 
Image as a Function
• function, f,
• f: R2  R
– f (x, y) gives the
intensity at position
(x, y)
Linear
Transformations
on Images
Transformations on Images
• Define a new image g in terms of an existing
image f
– We can transform:
1. either the domain
2. or the range of f
• Range transformation:
1. What kinds of operations on the image can this
transformation perform?
2. What t transformations are good for what?
Operations that change the domain of image
1. Some operations preserve the range but change
the domain of f
2. Here is an example that changes the domain :
–
What kinds of operations can this perform?
3. Still other operations operate on both the domain
and the range of f .
1.
2.
In general we are interested in all these transformations
Specifically we want them to have some nice (like linear)
properties
Linear Shift Invariant Systems (LSIS)
Linearity:
f1
g1
f1  f 2
g2
f2
g1  g2
Shift invariance:
f x  a 
a
g x  a 
a
Example of LSIS: ideal lens
Linear Shift Invariant Systems (LSIS)
Defocused image ( g ) is a
processed version of the focused
image ( f )
g
f
Ideal lens is a LSIS f x 
LSIS
g x 
Linearity: Brightness variation
Shift invariance: Scene movement
(this is not valid for lenses with non-linear distortions)
A Reminder of
a Convolution
concept
LSIS is Doing a Convolution
Linear Shift Invariant
Systems (LSIS)
convolution is linear and shift invariant
g x  

 f  hx   d
g  f h

f  
1. Continuous convolution
defined by integral
2. Discrete convolution
from previous lectures
defined by finite sums
We negate
(mirror) the kernel
of convolution h
x
h 

h  
kernel h


Convolution – two more Examples
Here we show two more
examples of convolution
of functions f and g
We use two nonidea approximation
of Dirac’s Delta
function as a signal
model
g x  

 f  hx   d
g  f h

f
g
f g
From Eric Weinstein’s Math World
Convolution – Another Example
g x  

 f  hx   d

a x 
b x 
1
-1
1
-1
1
1
c  a b
Here we convolve
signals a and b into
convolved signal c
c x 
1
-2
g  f h
-1
1
2
Convolution Kernel – Impulse Response
f
g
h
g  f h
• We ask question “What h will give us g = f ?”
Dirac Delta Function (Unit Impulse)
Sifting property:
We convolve f with Dirac Delta:




 f x xdx   f 0 xdx
1
 x 
2
 
 0

 f 0  x dx  f 0

And we get this

g x    f   x   d  f x 


    hx   d  hx 

Therefore g=f when we
convolve with delta
Point Spread Function
Now we apply information from last
slide to optical system
scene
Optical
System
image
• Ideally, the optical system should be a Dirac delta function. (no distortion)
• However, optical systems are never ideal.
 x 
point source
Optical
System
• Point spread function of Human Eyes
PSF x 
point spread function
PSF stands for Point
Spread Function
If the PSF of human eye is different
(not a good approximation of delta)
then we have eye troubles, related to
this PSF type
Point Spread Function
PSF stands for
Point Spread
Function
normal vision
myopia
hyperopia
Tell me what your PSF is and I will
tell you what is your eye problem.
astigmatism
Images by Richmond Eye Associates
Properties of Convolution
Convolution
operator
• Commutative
• Associative
a b  b  a
a  b c  a  b  c
This is unlike matrix
multiplication and
Kronecker (tensor)
products
• Cascade system
f
h1
 f
 f
g
h2
h1  h2
h2  h1
We used in cascading filters
in previous lectures. Now
we have a theoretical base
g
g
These are general
properties of convolution
, not related to signal
representation
• So now we know what is convolution and
why it is so important.
• However signals for convolution can be
represented in various ways.
• So the question is, “how to represent
signals?”
– 1D
– 2D
– 3D?
How to Represent Signals?
• Option 1: Taylor series represents any function using
polynomials.
• Polynomials are not the best - unstable and not very
physically meaningful. (they still have some applications)
• Easier to talk about “signals” in terms of its
“frequencies”
(how fast/often signals change, etc).
Mr. Jean Baptiste
Fourier and
Fourier Transform
idea
Jean Baptiste Joseph Fourier (1768-1830)
• Had crazy idea (1807):
•
Any periodic function
can be rewritten as a
weighted 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 true!
– called Fourier Series
– Possibly the greatest tool
used in Engineering
A Sum of Sinusoids
• Our building block:
Asin( x   
• Add enough of them to
get any signal f(x) you
want!
• How many degrees of
freedom?
• What does each control?
• Which one encodes the
coarse vs. fine structure of
the signal?
Math Review - Complex
numbers
• Real numbers:
1
-5.2

• Complex numbers
4.2 + 3.7i
9.4447 – 6.7i
-5.2 (-5.2 + 0i)
i  1
We often denote in EE i by j
Math Review - Complex
numbers
• Complex numbers
4.2 + 3.7i
9.4447 – 6.7i
-5.2 (-5.2 + 0i)
• Amplitude
A = | Z | = √(a2 + b2)
• Phase
 =  Z = tan-1(b/a)
• General Form
Z = a + bi
Re(Z) = a
Im(Z) = b
Real and imaginary parts
Math Review – Complex Numbers
• Polar Coordinate
Z = a + bi
• Amplitude
A = √(a2 + b2)
• Phase
 = tan-1(b/a)

A

a
b

Math Review – Complex
Numbers and Cosine Waves
• Cosine wave has three properties
– Frequency
– Amplitude
– Phase
• Complex number has two properties
– Amplitude
– Wave
• Complex numbers to represent cosine waves at varying frequency
– Frequency 1:
– Frequency 2:
– Frequency 3:
Z1 = 5 +2i
Z2 = -3 + 4i
Z3 = 1.3 – 1.6i
Simple but great idea !!
Fourier Transform Idea
• We want to understand the frequency  of our signal. So, let’s
reparameterize the signal by  instead of x:
Fourier
Transform
f(x)
F()
• For every  from 0 to infinity, F() holds the amplitude A and
phase  of the corresponding sine
Asin( x   
Complex number trick!
F ( )  R( )  iI ( )
– How can F hold both?
A   R( )  I ( )
2
F()
2
Inverse Fourier
Transform
I ( )
  tan
R( )
1
f(x)
All info
preserved
Math Review Periodic
Functions –
more illustrations
to our formalism
Math Review - Periodic Functions
If there is some a, for a function f(x), such that
f(x) = f(x + na)
then function is periodic with the period a
a
0
2a
3a
Math Review - Attributes of
cosine wave
5
3
f(x) = cos (x)
1
-10
-5
-1 0
5
10
5
10
5
10
-3
-5
5
4
Amplitude
3
f(x) = 5 cos (x)
2
1
0
-10
-5
-1 0
-2
-3
-4
-5
5
Phase
f(x) = 5 cos (x + 3.14)
3
1
-10
-5
-1 0
-3
-5
Math Review - Attributes of
cosine wave
5
4
3
Amplitude
f(x) = 5 cos (x)
2
1
0
-10
-5
-1 0
5
10
5
10
5
10
-2
-3
-4
-5
5
Phase
3
f(x) = 5 cos (x + 3.14)
1
-10
-5
-1 0
-3
-5
5
Frequency f(x) = 5 cos (3 x + 3.14)
3
1
-10
-5
-1 0
-3
-5
Math Review - Attributes of
cosine wave
5
3
f(x) = cos (x)
1
-10
-5
-1 0
5
-3
-5
f(x) = A cos (kx + )
Amplitude, Frequency, Phase
This math is base of Fourier Analysis and Fourier Synthesis
10
Time and Frequency
• example : g(t) = sin(2 f t) + (1/3)sin(2(3f) t)
Our periodic signal is a mixture of two
frequencies
Before we go to using
complex numbers and
formulas let us get intuition
about time and frequency
Time and Frequency
• example : g(t) = sin(2 f t) + (1/3)sin(2(3f) t)
=
+
Frequency Spectra
• example : g(t) = sin(2pf t) + (1/3)sin(2p(3f) t)
=
+
In time domain
In frequency domain
Frequency Spectra
• Usually, frequency is more interesting than the phase
1.
2.
We will decompose it to frequencies,
we will add one more coefficient based base
function in each slide
Frequency Spectra
First base function
=
=
second base
function
+
Sum of first and
second base
functions is already
not bad
approximation
Frequency Spectra
=
=
+
Sum of first , second
and third base
functions is even
better approximation
Frequency Spectra
=
=
+
Frequency Spectra
=
=
+
Frequency Spectra
=
=
+
Frequency Spectra

1
= A sin(2 kt )
k 1 k
In time domain
In frequency
domain
Frequency Spectra
Observe the
negative values
of domain that
appear in
spectrum
Fourier
Transform is
Just a change of
basis
Important concept
Fourier Transform is Just a change of basis
M * f(x) = F()
*
.
.
.
M can be a matrix, a
decision diagram, a butterfly,
a system (physical or
simulated)
=
IFT: Just a change of basis
M-1
* F() = f(x)
*
.
.
.
Inverse Fourier
Transform is similarly
only a change of basis
=
Fourier Transform is an invertible
operator
Image
Fourier Transform
FT
v2 will display image or its transform, Matlab examples.
Fourier Transform is an invertible
operator
Image
Fourier Transform
Ny ⁄ 2
Ny
y
F(kx,ky)
f(x,y)
0
Nx
x
{F(kx,ky)} = f(x,y)}
{f(x,y)}
= F(kx,ky)
Similarly as for Hadamard-Walsh or
Reed-Muller Transforms we can talk
about transform pairs.
Nx ⁄ 2
Continuous
Fourier
Transform
Continuous Fourier Transform
F (s)  

f ( x)e i 2xs dx

f(x)
= F(s)

f ( x)   F ( s )e i 2xs ds


F ( s )   f ( x)e ixs dx

1
f ( x) 
2
Euler’s Formula
1
F ( s) 
2
1
f ( x) 
2









F ( s )eixs ds
f ( x)e ixs dx
F ( s)eixs ds
Various normalizations of pairs
Some Conventions for Fourier
Transforms
• Image Domain
• Fourier Domain
–
–
–
–
–
• Forward Transform
• f(x,y,z)
• g(x)
• F
Reciprocal space
Fourier Space
K-space
Frequency Space
Spectral domain
• Reverse Transform,
Inverse Transform
• F(kx,ky,kz)
• G(s)
• F
Fourier
Analysis
Decompose f(x) into a series of cosine
waves that when summed reconstruct f(x)
5
4
Fourier Analysis in 1D. Audio
signals
Amplitude Only
3
2
1
0
-1 0
200
400
600
800
1000
1200
1400
-2
5
-3
-4
10
15
(Hz)
-5
5
4
3
2
1
0
-1 0
-2
200
400
600
800
1000
1200
1400
5
10
-3
-4
-5
(Hz)
15
5
4
Fourier Analysis in 1D. Audio
signals
3
2
1
0
-1 0
200
400
600
800
-2
1000
1200
1400
5
10
-3
-4
-5
Your ear performs fourier analysis.
(Hz)
15
Fourier Analysis in 1D.
Spectrum Analyzer.
iTunes performs fourier analysis.
1.
2.
We know it from our Hi-Fi equipment
The same can be done with standard images or radar images or any other images, in any
number of dimensions
Fourier
Synthesis
Summing cosine waves reconstructs
the original function
Example: Fourier Synthesis of
Boxcar Function
Boxcar function
Periodic Boxcar
Can this function be reproduced with cosine waves?
k=1. One cycle per period
A1·cos(2kx + 1)
k=1
 A ·cos(2kx +  )
1
k
k=1
k
k=2. Two cycles per period
A2·cos(2kx + 2)
k=2
 A ·cos(2kx +  )
2
k
k=1
k
k=3. Three cycles per period
A3·cos(2kx + 3)
k=3
 A ·cos(2kx +  )
3
k
k=1
k
Fourier Synthesis. N Cycles
A3·cos(2kx + 3)
k=3
 A ·cos(2kx +  )
N
k
k=1
k
Example: Fourier Synthesis of a
2D Function
An image is two dimensional
data.
Intensities as a function of x,y
White pixels represent the
highest intensities.
Greyscale image of iris
128x128 pixels
Fourier Synthesis of a 2D
Function
F(2,3)
1. We select coefficients, multiply them by weights and add
2. We can create artificial images
Fourier Filters
• Fourier Filters in 1D data (sound) change
the image by changing which frequencies
of cosine waves go into the image
• Represented by 1D spectral profile
• 2D Specteral Profile is rotationally
symmetrized 1D profile
Low and High Frequency Image
Filters
• Low frequency terms
– Close to origin in Fourier Space
– Changes with great spatial extent (like ice
gradient), or particle size
• High frequency terms
– Closer to edge in Fourier Space
– Necessary to represent edges or highresolution features
Frequency-based Filters
1. Low-pass Filter (blurs)
–
Restricts data to low-frequency components
2. High-pass Filter (sharpens)
–
Restricts data to high-frequency-components
3. Band-pass Filter
–
Restrict data to a band of frequencies
4. Band-stop Filter
–
Suppress a certain band of frequencies
We will show many
practical examples
soon
Cutoff
Low-pass
Filter
Image is blurred
Sharp features are lost
Ringing artifacts
Filter
characteristics
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
0.2
0.4
0.6
0.8
1
Butterworth Low-pass Filter
Flat in the pass-band
Zero in the stop-band
No ringing
Gaussian Low-pass Filter
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
0.2
0.4
0.6
0.8
1
Butterworth High-pass Filter
• Note the loss of solid
densities
How the filter looks in 2D
unprocessed
bandpass
lowpass
highpass
spectra
Convolution
again
Our goal is now to approach the
most important results of
spectral processing
Convolution
Convolution of some function f(x) with
some kernel g(x)
Continuous
Discrete
*
=
Convolution in 2D
image
kernel
result
x
x

=
x
x
x
x
=
x
x
x

x
x
x x
x
x
x
Microscope Point-SpreadFunction is Convolution
Point-Spread-Function
Convolution Theorem
fg=
{FG}
f=
FG
G
Convolution in image domain
Is equivalent to multiplication in
fourier domain
Can be derived from
above because FG/G
= F and f = F-1 (F)
Lowpass Filtering by Convolution
fg=
1
{FG}
0.9
0.8
0.7
•
0.6
Camera shake
0.5
0.4
0.3
We showed examples of this using
standard convolution, but it is done
in modern equipment through
digitally calculated FFT and IFFT
0.2
0.1
0
0
0.2
0.4
0.6
0.8
1
2D gaussian is
example of kernel
2D gaussian as a
matrix of a kernel
Lens Performs Fourier
Transform
Contrast Theory
Power spectrum
PS = F2(s) CTF2(s) Env2(s) + N2(s)
Incoherant average of transform
obs(x) = f(x)  psf(x)  env(x) + n(x)
observed image
noise
f(x) for true particle
envelope function
point-spread function
Gibbs Ringing
• 5 waves
• 25 waves
• 125 waves
Even with 125 waves we still
have some small Gibbs Ringing
Fourier
Transform
Math
Properties
Fourier Transform – more formally
Represent the signal as an infinite weighted sum
of an infinite number of sinusoids

F u    f x e
i 2ux

dx
Note: e  cos k  i sin k
ik
i  1
Arbitrary function
Single Analytic Expression
Spatial Domain (x)
Frequency Domain (u)
(Frequency Spectrum F(u))
Inverse Fourier Transform (IFT)

f x    F u e

i 2ux
dx
Fourier Transform
• Also, defined as:
Fourier
Transform

F u    f x eiuxdx

Note: e  cos k  i sin k
ik
• Inverse Fourier Transform (IFT)
1
f x  
2



F u eiuxdx
i  1
Fourier
Transform
Pairs (I)
Fourier
Transform
Note that these are derived using
angular frequency ( e  iux )
Fourier Transform Pairs (I)
Fourier Transform
Note that these are derived using
angular frequency ( e  iux )
Fourier Transform Pairs (I)
Note that these are
derived using
angular
frequency ( e  iux )
Fourier Transform
and Convolution
Fourier Transform and Convolution
Let
g  f h
Here we derive the most important
formula of image processing and
other applications of Fourier
Transform

Then Gu    g x ei 2uxdx

From definition




 
f  hx   ei 2ux ddx
  f  e
   f  e



i 2u
 


i 2u
d
Add and
subtract tau

d hx   ei 2u  x  dx
 hx'e


i 2ux'

dx'

Separate
domains
 F u H u 
Convolution in spatial domain
 Multiplication in frequency domain
MAIN THEOREM OF
MODERN IMAGE
PROCESSING:
Convolution in spatial domain
 Multiplication in frequency domain
Fourier Transform and Convolution
Spatial Domain (x)
Frequency Domain (u)
g  f h
g  fh
G  FH
G  F H
So, we can find g(x) by Fourier transform
g

IFT
G
f

FT

F
h
FT

H
Properties of
Fourier
Transform
Properties of Fourier Transform 1
Spatial Domain (x)
Frequency Domain (u)
c1 f x  c2 g x
c1F u   c2Gu 
Scaling
f ax
1 u
F 
a a
Shifting
f x  x0 
e i 2ux0 F u 
Symmetry
F x 
f  u 
Conjugation
f  x 
Convolution
f x   g x 
F   u 
Differentiation
d n f x 
dx n
i 2u n F u 
Linearity
F u Gu 
Note that these are derived using
frequency ( e  i 2ux )
Properties of Fourier Transform 2
Star will be used for conjugation in
complex plane
Substitute f instead of g in this formula to
get the upper formula
For instance, when f is R
then F is RE, IO,. etc
1.
2.
3.
4.
5.
6.
RO – real odd
IE – imaginary even
IO – imaginary odd
RE – Real even
CE – Complex even
CO – Complex odd
We explained even-odd
properties for Walsh, we will
explain again for other
transforms
Example use of Fourier: Smoothing/Blurring
1.
We want a smoothed function of f(x)
g x  f x hx
But we do this in
spectral domain
•2. Let us use a Gaussian kernel
h x 
 1 x2 
1
hx  
exp 
2
2

2 


• 3. Then
 1

2
H u   exp  2u   2 
 2

Gu   F u H u 

x
H u 
1
2
u
H(u) attenuates high frequencies in F(u) (Low-pass Filter)!
Our Images are
Discrete
so discrete Fourier applies
Image as a Discrete Function
Digital Images
The scene is
– projected on a 2D plane,
– sampled on a regular grid, and each sample is
– quantized (rounded to the nearest integer)
f i, j   Quantize  f i, j
Image as a matrix
Review
Fourier Transform is invertible operator
Math Review
Periodic functions
Amplitude, Phase and
Frequency
Complex number
Amplitude and Phase
Fourier Filters
Low-pass
High-pass
Band-pass
Band-stop
Convolution Theorem
Fourier Analysis (Forward Transform)
Decomposition of periodic signal
into cosine waves
Fourier Synthesis (Inverse Transform)
Summation of cosine waves into
multi-frequency waveform
Fourier Transforms in 1D, 2D, 3D, ND
Image Analysis
Image (real-valued)
Transform (complex-valued,
amplitude plot)
Deconvolute by Division in
Fourier Space
All Fourier Filters can be
expressed as real-space
Convolution Kernels
Lens does Fourier transforms
Diffraction Microscopy
Where we are in this class?
•
•
•
•
•
We defined continuous Fourier Transform
We defined discrete Fourier Transform DFT
DFT is calculated by matrix multiplication
It is slow
We will need factorizations, Kroneckers,
recursions and butterflies to calculate it
faster.
• Only then FT will become a really useful
tool, FFT.
Further Reading
• Wikipedia
• Mathworld
• The Fourier Transform and its
Applications. Ronald Bracewell
Filtering with EMAN2
LowPass Filters
filtered=image.process(‘filter.lowpass.guass’, {‘sigma’:0.10})
filtered=image.process(‘filter.lowpass.butterworth’, {‘low_cutoff_frequency’:0.10, ‘high_cutoff_frequency’:0.35})
filtered=image.process(‘filter.lowpass.tanh’, {‘cutoff_frequency’:0.10, ‘falloff’:0.2})
HighPass Filters
filtered=image.process(‘filter.highpass.guass’, {‘sigma’:0.10})
filtered=image.process(‘filter.highpass.butterworth’, {‘low_cutoff_frequency’:0.10, ‘high_cutoff_frequency’:0.35})
filtered=image.process(‘filter.highpass.tanh’, {‘cutoff_frequency’:0.10, ‘falloff’:0.2})
BandPass Filters
filtered=image.process(‘filter.bandpass.guass’, {‘center’:0.2,‘sigma’:0.10})
filtered=image.process(‘filter.bandpass.butterworth’, {‘low_cutoff_frequency’:0.10, ‘high_cutoff_frequency’:0.35})
filtered=image.process(‘filter.bandpass.tanh’, {‘cutoff_frequency’:0.10, ‘falloff’:0.2})
Sources of slides
1. S. Narasimhan
2. Mike Marsh
National Center for Macromolecular
Imaging
Baylor College of Medicine