To perform image restoration (and many other useful image
processing algorithms) in a computer, we need a Fourier
Transform (FT) that is discrete and two-dimensional.
Fkl=Fuv|
u=2πkN
,
v=2πlN
F
k
l
u
2
k
N
v
2
l
N
F
u
v
(1)
for
k=0…N−1
k
0
…
N
1
and
l=0…N−1
l
0
…
N
1
.
Fuv=∑m∑nfmne−(ium)e−(ivm)
F
u
v
m
n
f
m
n
u
m
v
m
(2)
Fkl=∑m=0N−1∑n=0N−1fmne(−i)2πkmNe(−i)2πlnN
F
k
l
m
N
1
0
n
N
1
0
f
m
n
2
k
m
N
2
l
n
N
(3)
where the above equation (
Equation 3)
has finite support for an
NNx
NN
image.
As with our regular fourier transforms, the 2D DFT also has
an inverse transform that allows us to reconstruct an image
as a weighted combination of complex sinusoidal basis
functions.
fmn=1N2∑k=0N−1∑l=0N−1Fklei2πkmNei2πlnN
f
m
n
1
N
2
k
N
1
0
l
N
1
0
F
k
l
2
k
m
N
2
l
n
N
(4)
The regular 2D convolution equation is
gmn=∑k=0N−1∑l=0N−1fklhm−kn−l
g
m
n
k
N
1
0
l
N
1
0
f
k
l
h
m
k
n
l
(5)
Below we go through the steps of convolving two
two-dimensional arrays. You can think of
ff as representing an image and
hh represents a PSF, where
hmn=0
h
m
n
0
for
m∧n>1
m
n
1
and
m∧n<0
m
n
0
.
h=(
h00h01
h10h11
)
h
h
0
0
h
0
1
h
1
0
h
1
1
f=(
f00…f0N−1
⋮⋱⋮
fN−10…fN−1N−1
)
f
f
0
0
…
f
0
N
1
⋮
⋱
⋮
f
N
1
0
…
f
N
1
N
1
Step 1 (Flip hh):
h−m−n=(
h11h100
h01h000
000
)
h
m
n
h
1
1
h
1
0
0
h
0
1
h
0
0
0
0
0
0
(6)
Step 2 (Convolve):
g00=h00f00
g
0
0
h
0
0
f
0
0
(7)
We use the standard 2D convolution equation (
Equation 5) to find the first element of
our convolved image. In order to better understand what is
happening, we can think of this visually. The basic idea is
to take
h−m−n
h
m
n
and place it "on top" of
fkl
f
k
l
, so that just the bottom-right element,
h00
h
0
0
of
h−m−n
h
m
n
overlaps with the top-left element,
f00
f
0
0
, of
fkl
f
k
l
. Then, to get the next element of our convolved
image, we slide the flipped matrix,
h−m−n
h
m
n
, over one element to the right and get the
following result:
g01=h00f01+h01f00
g
0
1
h
0
0
f
0
1
h
0
1
f
0
0
We continue in this fashion to find all of the elements of
our convolved image,
gmn
g
m
n
. Using the above method we define the general
formula to find a particular element of
gmn
g
m
n
as:
gmn=h00fmn+h01fmn−1+h10fm−1n+h11fm−1n−1
g
m
n
h
0
0
f
m
n
h
0
1
f
m
n
1
h
1
0
f
m
1
n
h
1
1
f
m
1
n
1
(8)
What does
HklFkl
H
k
l
F
k
l
produce?
2D Circular Convolution
g
~
mn=IDFTHklFkl=circular convolution in 2D
g
~
m
n
IDFT
H
k
l
F
k
l
circular convolution in 2D
(9)
Due to periodic extension by DFT (Figure 2):
Based on the above solution, we will let
g
~
mn=IDFTHklFkl
g
~
m
n
IDFT
H
k
l
F
k
l
(10)
Using this equation, we can calculate the value for each
position on our final image,
g
~
mn
g
~
m
n
. For example, due to the periodic extension of
the images, when circular convolution is applied we will
observe a
wrap-around effect.
g
~
00=h00f00+h10fN−10+h01f0N−1+h11fN−1N−1
g
~
0
0
h
0
0
f
0
0
h
1
0
f
N
1
0
h
0
1
f
0
N
1
h
1
1
f
N
1
N
1
(11)
Where the last three terms in
Equation 11 are a result of the wrap-around effect caused
by the presence of the images copies located all around it.
If the support of hh is
MMxMM
and ff is
NNxNN,
then we zero pad ff and
hh to
M+N−1
M
N
1
x
M+N−1
M
N
1
(see Figure 3).
Circular Convolution = Regular Convolution
Fkl=∑m=0N−1∑n=0N−1fmne(−i)2πkmNe(−i)2πlnN
F
k
l
m
N
1
0
n
N
1
0
f
m
n
2
k
m
N
2
l
n
N
(12)
where in the above equation,
∑n=0N−1fmne(−i)2πlnN
n
N
1
0
f
m
n
2
l
n
N
is simply a 1D DFT over
nn.
This means that we will take the 1D FFT of each row; if we
have
NN rows, then it will
require
NlogN
N
N
operations per row. We can rewrite this as
Fkl=∑m=0N−1
f
′
mle(−i)2πkmN
F
k
l
m
N
1
0
f
′
m
l
2
k
m
N
(13)
where now we take the 1D FFT of each column, which means
that if we have
NN columns,
then it requires
NlogN
N
N
operations per column.
Therefore the overall complexity of a 2D FFT is
ON2logN
O
N
2
N
where
N2
N
2
equals the number of pixels in the image.
