Transcript Fast Polar Fourier Transform Michael Elad* Scientific Computing and Computational Mathematics Stanford University FoCM Conference, August 2002 Image and Signal Processing Workshop IMA - Minneapolis * Joint.

Fast Polar Fourier Transform
Michael Elad*
Scientific Computing and Computational Mathematics
Stanford University
FoCM Conference, August 2002
Image and Signal Processing Workshop
IMA - Minneapolis
* Joint work with Dave Donoho (Stanford-Statistics),
Amir Averbuch (TAU-CS-Israel), and Ronald Coifman (Yale-Math)
-1-
Collaborators
Dave Donoho
Amir Averbuch
Ronald Coifman
Statistics Department
Stanford
CS Department
Tel-Aviv University
Math. Department
Yale
-2-
Fast Polar Fourier Transform
 FFT is one of top 10 algorithms of 20th century.
 We'll extend utility of FFT algorithms to new
class of settings in image processing.
 Create a tool which is:

Free of emotional involvement, &

Freely available over the internet.
 Current Stage:

We have the tool, and its analysis,

Have not demonstrated applications yet.
-3-
Agenda
1. ThinkingThinking
Polar – Continuum
Polar – Continuum
2. Thinking Polar – Discrete
3. Current State-Of-The-Art
Background
&
Motivation
4. Our Approach - General
5. The Pseudo-Polar Fast Transform
6. From Pseudo-Polar to Polar
7. Algorithm Analysis
New
Approach
and its
Results
8. Conclusions
-4-
1. Thinking Polar - Continuum
 For today f(x,y) function of (x,y)2
 Continuous Fourier Transform
ˆf u, v   f x , y   f x , y  exp ixu  iyvdxdy

 Polar coordinates: u=r·cos() , v=r·sin()
~
f r ,   ˆf r  cos(), r  sin() 

  f x, y  exp ixr  cos()  iy  sin()dxdy
 Important Processes easy to continuum polar domain.
-5-
1. Thinking Polar - Continuum
v
r
-5
0
-4
-3
1
r
-2

-1
u
2
3
0
1
4
2
5
3
4
5
6
-5
-4
-3
-2
-1
0
1
2
3
4
5
0
1
2
3
4
5
6
7

-6-
1. Thinking Polar - Continuum
Natural Operations: 1. Rotation
Using the polar coordinates, rotation is simply a shift in
the angular variable.
 Q0 – planar rotation by 0 degrees


 Rotation f x, y   f Q x, y
0
0
 In polar coordinates – shift in angular variable
~
~
f0 r ,   f r ,   0 
-7-
1. Thinking Polar - Continuum
Natural Operations: 2. Scaling
Using the polar coordinates, 2D scaling is simply a 1D
scaling in the radial variable.
 S – planar scaling by a factor 
 Scaling f x, y   f Sx, y
 In polar coordinates – 1D scale in radial variable
~
~
f r ,   Const  f r , 
 Log-Polar – shift in the radial variable.
-8-
1. Thinking Polar - Continuum
Natural Operations: 3. Registration
Using the polar coordinates, rotation+shift registration
simply amounts to correlations.


 f(x,y) and g(x,y): f x, y   g Q x, y  x0 , y0
0
 Goal: recover x 0 , y 0 , 0 .
~
 Angular cross-correlation between f r ,  and ~
g r , 
– Estimate 0.
 Rotation & cross-correlation on regular Fourier transform
gives the shift.
-9-
1. Thinking Polar - Continuum
Natural Operations: 4. Tomography
Using the polar coordinates, we obtain a method to
obtain the Inverse Radon Transform.
 Radon Transform:
Rf t ,    f ( x , y)x cos()  y sin()  t dxdy
 Goal: Given Rf(t,), recover f.
~
 Projection-Slice-Theorem: 1Rf t ,   f r , .
~
 Reconstruction: Rf  f  ˆf  f .
- 10 -
1. Thinking Polar - Continuum
More Natural Operations
New orthonormal bases:
 Ridgelets,

 Curvelets,
 Fourier Integral operations,
 Ridgelet packets.
Analysis of textures.
r
Analysis of singularities.
More …
- 11 -
Agenda
1. Thinking Polar – Continuum
2.
Thinking Polar – Discrete
3. Current State-Of-The-Art
4. Our Approach - General
5. The Pseudo-Polar Fast Transform
6. From Pseudo-Polar to Polar
7. Algorithm Analysis
8. Conclusions
- 12 -
2. Thinking Polar - Discrete
 Certain procedures very important to digitize
 Rotation,
 Scaling,
 Registration,
 Tomography, and
 …
 These look so easy in continuous theory – Can’t we
use it in the discrete domain?
 Why not Polar-FFT?
- 13 -
2. Thinking Polar - Discrete
The FFT Miracles
 1D Discrete Fourier Transform
 Uniformly sampled in time and frequency – FFT.
 Complexity – O(5Nlog2N) instead of O(N2).
 2D Discrete Fourier Transform
 Cartesian grid in space and frequency – Separability
 Only 1D-FFT operations.
 Smart memory management.
- 14 -
2. Thinking Polar - Discrete
2D DFT – Cartesian Grid

N
1
2n1 N  2
y
2
N
x 




2

n
N
2  n ,n   N
 y
1 2
2

-
2
N
x
-
- 15 -
2D FFT Complexity
 Complexity: O(10N2log2N)
instead of O(N4).
 Important Feature: All
operations are 1D
– leading to
efficient cache usage
Cartesian N-by-N
Data
5N2logN
1D FFT to columns
5N2logN
1D FFT to rows
Cartesian
2D-FFT 10N2logN
- 16 -
2. Thinking Polar - Discrete
Discrete Polar Coordinates?
y
Choice of grid?
NS 1

n1  r
r 

NS

r  n1  0

NSr

NS 1

2n2  
,  

NS
  n2  0

Resulting with NS
rays with NSr elements

-
x
on each:
For S=Sr=1, we
have N2 grid points.
2
NS 
-
- 17 -
2. Thinking Polar - Discrete
Grid Problematics
y

 Grid spacing?
 Fate of corners?

-
x
 No X-Y separability !!
-
- 18 -
2. Thinking Polar - Discrete
Polar FFT - Current Status
 Current widespread belief - There cannot be a
fast method for DFT on the polar grid. See e.g.
The DFT: an owner’s manual, Briggs and
Henson, SIAM, 1995, p. 284.
 Consequence of Non-existence:
 Continuous Fourier – vague inspiration only.
 Fourier in implementations widely deprecated.
 Fragmentation: each field special algorithm.
- 19 -
Agenda
1. Thinking Polar – Continuum
2. Thinking Polar – Discrete
3.
Current State-Of-The-Art
4. Our Approach - General
5. The Pseudo-Polar Fast Transform
6. From Pseudo-Polar to Polar
7. Algorithm Analysis
8. Conclusions
- 20 -
3. Current State-Of-The-Art
 Assessing T: Unequally-spaced FFT (USFFT)
 Data in Cartesian set.
 Approximate transform in non-Cartesian set.
 Oriented to 1D – not 2D and definitely not Polar.
 Assessing T †: For Tomography
 Data in Polar coordinates in frequency.
 Approximate inverse transform to Cartesian grid.
 Inspired by the projection-slice-theorem.
- 21 -
3. Current State-of-the-Art
USFFT - Rational
y
3
+ Destination Polar
grid
O Critically sampled
Cartesian grid
2
1
0
o Over-sampled
Cartesian grid
-1
-2
-3
-3
-2
-1
0
1
2
3
x
- 22 -
3. Current State-of-the-Art
USFFT - Detailed
 Over-sample Cartesian grid.
 Rapidly evaluate FT:
 Values F.
 Possibly - partial derivatives.
 Associate Cartesian neighbors to each
polar grid point.
 Approximate interpolation.
- 23 -
3. Current State-of-the-Art
Our Reading of Literature
 Boyd (1992) – Over-sampling and interpolation
by Euler sum or Langrangian interpolation.
 Dutt-Rokhlin (1993,1995) - Over-sampling and
interpolation by the Fast-Multipole method.
 Anderson-Dahleh (1996) – Over-sampling and
obtaining the partial derivatives, and then
interpolating by Taylor series.
 Ware (1998) – Survey on USFFT methods.
- 24 -
3. Current State-of-the-Art
USFFT Problematics
 Several involved parameters:
 Over-sampling factor,
 Method of interpolation, and
 Order of interpolation.
 Good accuracy calls for extensive over-sampling.
 Correspondence overhead: spoils vectorizability of
algorithm - causes high cache misses.
 Emotionally involved.
- 25 -
Agenda
1. Thinking Polar – Continuum
2. Thinking Polar – Discrete
3. Current State-Of-The-Art
4.
Our Approach - General
5. The Pseudo-Polar Fast Transform
6. From Pseudo-Polar to Polar
7. Algorithm Analysis
8. Conclusions
- 26 -
4. Our Approach - General
We propose a
Fast Polar Fourier Transform
with the following features:
 Low complexity – O(Const·N2log2N)
 Vectorizability – 1D operations only
 Non-Expansiveness – Factor 2 (or 4) only
 Stability – via Regularization
 Accuracy – 2 orders of magnitude over USFFT methods
- 27 -
Our Strategy
Fast and Exact
Fourier Trans.
on a polar-like
grid
1D
interpolations
to the polar
grid
Pseudo
Polar
Grid
- 28 -
Agenda
1. Thinking Polar – Continuum
2. Thinking Polar – Discrete
3. Current State-Of-The-Art
4. Our Approach - General
5.
The Pseudo-Polar Fast Transform
6. From Pseudo-Polar to Polar
7. Algorithm Analysis
8. Conclusions
- 29 -
5. The Pseudo-Polar FFT
 Developed by Averbuch, Coifman, Donoho,
Israeli, and Waldén (1998).
 Basic idea: A “Polar-Like” grid that enables
 EXACT Fourier transform,
 FAST computation,
 1D operations only.
 Applications: Tomography, image processing,
Ridgelets, and more.
- 30 -
5. The Pseudo-Polar FFT
The Pseudo-Polar Skeleton
NSr equi-spaced
concentric squares,
NSt ‘equi-spaced’
(not in angle)
We separate our
treatment to
basically vertical
and basically
horizontal lines.
y
3
2
1
0
-1
-2
-3
-3
-2
-1
0
1
2
3
x
- 31 -
5. The Pseudo-Polar FFT
Fast Fourier Transform
 The destination samples
are uniformly sampled
vertically,
 Per each row, destination
samples are uniformly
sampled horizontally,
 Fractional Fourier
Transform is the answer
(Chirp-Z), with complexity:
O(20Nlog2N).
[Why?]
y
3
2
1
0
-1
-2
-3
-3
-2
-1
0
1
2
3
x
- 32 -
5. The Pseudo-Polar FFT
PP-FFT versus 2D-FFT
Cartesian N-by-N
Data
PP-FFT
5N2logN
1D FFT to columns
20N2logN
1D FRFFT to rows
PP-FFT
25N2logN vertical
2D-FFT
5N2logN
1D FFT to columns
5N2logN
1D FFT to rows
Cartesian
2D-FFT 10N2logN
- 33 -
5. The Pseudo-Polar FFT
The PP-FFT - Properties
 Exact in exact arithmetic.
 No parameters involved !!
 Complexity - O(50·N2log2N) versus O(N4).
 1D operations only.
 For the chosen grid (Sr=St=2) - 5.
- 34 -
Agenda
1. Thinking Polar – Continuum
2. Thinking Polar – Discrete
3. Current State-Of-The-Art
4. Our Approach - General
5. The Pseudo-Polar Fast Transform
6.
From Pseudo-Polar to Polar
7. Algorithm Analysis
8. Conclusions
- 35 -
6. From Pseudo-Polar to Polar
Fast and Exact
Fourier Trans.
on a polar-like
grid
2 stages of 1D
interpolations
to get to the
polar grid
- 36 -
6. From Pseudo-Polar to Polar
The Interpolation Stages
y
3
The original PseudoPolar Grid
Warping to equi-spaced
angles
Warping each ray to
have the same step
2
1
0
-1
-2
-3
-3
-2
-1
0
1
2
3
x
- 37 -
6. From Pseudo-Polar to Polar
First Interpolation Stage
y
Rotation of the rays to
have equi-spaced angles
(S-Pseudo-Polar grid):
 Every row is a
trigonometric polynomial of
order N,
 FRFT on over-sampled
array and 1D interpolation,
 Very accurate.
3
2
1
0
-1
-2
-3
-3
-2
-1
0
1
2
3
x
- 38 -
6. From Pseudo-Polar to Polar
The Required Warping
New x
Basically vertical lines:
3
NS t
2
1


2
2m


,



 y
x
y
NSr
NSt

,m  NS t
2
y  
2
1
0
NS t
-1
1

 m   2
 
x  y  tan
 2NSt  m NS t

-2
2
[Why?]
-3
-3
-2
-1
0
1
2
3
Original x
- 39 -
6. From Pseudo-Polar to Polar
The Actual Interpolation
Cartesian
Data
N-by-N
PP-FFT
5N2logN
S-PP-FFT
5N2logN
1D FFT to columns
20N2logN
1D FRFFT to rows
1D FFT to columns
20N2S·log(NS)
1D Over-sampled (S) FRFFT
to rows
O{N2}
25N2logN
1D Interpolation
PP-FFT Vertical
(20S+5)N2logN
S-PP-FFT Vertical
- 40 -
6. From Pseudo-Polar to Polar
Second Interpolation Stage
y
 As opposed to the previous
step, the rays are not
trigonometric polynomials
of order N,
 We proved that the rays
are band-limited (smooth)
functions,
 Over-sampling and
interpolation is expected to
perform well.
3
2
1
0
-1
-2
-3
-3
-2
-1
0
1
2
3
x
- 41 -
6. From Pseudo-Polar to Polar
Over-Sampling Along Rays
 Over-sampling along
rays cannot be done
by taking the 1D ray
and over-sampling it.
 Sr>1:
 More concentric squares.
 Sr longer 1D-FFT’s at the
beginning of the algorithm.
 Sr times FRFFT operations.
y
3
2
1
0
-1
-2
-3
-3
-2
-1
0
1
2
3
x
- 42 -
6. From Pseudo-Polar to Polar
The Actual Interpolation
Cartesian
Data
N-by-N
N·5(NSr)·log(NSr)
Full Polar FFT
O{40SrStN2logN}
1DFFT to over-sampled columns
NSr·20(NSt)·log(NSt)
1D Over-sampled (S) FRFFT
to rows
O{(NSr)·N}
Polar-FFT Vertical
O{N·N}
1D Interpolation
1D Interpolation
S-PP-FFT Vertical
- 43 -
6. From Pseudo-Polar to Polar
To Summarize
We propose a
Fast Polar Fourier Transform
with the following features:
 Low complexity – O(Const·N2log2N)
 Vectorizability – 1D operations only
 Non-Expansiveness – Factor 2 (or 4) only
 Stability – via Regularization
 Accuracy – 2 orders of magnitude over USFFT methods
- 44 -
Agenda
1. Thinking Polar – Continuum
2. Thinking Polar – Discrete
3. Current State-Of-The-Art
4. Our Approach - General
5. The Pseudo-Polar Fast Transform
6. From Pseudo-Polar to Polar
7.
Algorithm Analysis
8. Conclusions
- 45 -
7. Algorithm Analysis
We have a code performing the Polar-FFT in Matlab:
Y=Polar_FFT(X);
OR
Y=Polar_FFT(X,St,Sr);
Where: X
– Input array of N-by-N samples
St,Sr – Over-sampling factors in the approximations
Y
– Polar-FFT result as an 2N-by-2N array with rows being
the rays and columns being the concentric circles.
- 46 -
7. Algorithm Analysis
The Implementation
 The current Polar-FFT code implements Taylor
method for the first interpolation stage and spline
ONLY (no-derivatives) for the second stage.
 For comparison, we demonstrate the performance
of the USFFT method with over-sampling S and
interpolation based on the Taylor interpolation
(found to be the best).
- 47 -
7. Algorithm Analysis
Error for Specific Signal
||Approximation error||1
2
10
• Input is random
32-by-32 array,
• USFFT method
uses one
parameter
whereas there
are two for the
up-sampling in
the Polar-FFT.
• Thumb rule:
Sr·St= S2.
Thumb rule: Sr=4St
1
10
0
10
St=1
-1
Taylor USFFT
10
-2
St=2
10
-3
10
St=4
St=3
10
20
30
40
St·Sr
or S2
50
60
70
80
90
100
- 48 -
7. Algorithm Analysis
Error For Specific Signals
 Similar curves obtained of ||*|| and ||*||2 norms.
 Similar behavior is found for other signals.
 Conclusion: For the proper choice of St and Sr, we
get 2-orders-of-magnitude improvement in the
accuracy comparing to the best USFFT method.
 Further improvement should be achieved for Taylor
interpolation in the second stage.
 US-FFT takes longer due the 2D operations.
- 49 -
7. Algorithm Analysis
The Transform as a Matrix
Approximate
All the involved
transformations (accurate
and approximate) are
linear - they can be
represented as a matrix of
size 4N2-by-N2.
Ya=Ax
Or
Yt=Tx
True
- 50 -
7. Algorithm Analysis
Regularization of T/A
 An input signal of N-by-N is transformed to an
array or 2N-by-2N.
 For N=16, T size is 1024-by-256, and 60,000
(bad for inversion).
 Adding the assumption that the Frequency
corners should be zeroed, we obtain
y  TPolar x
 TPolar 
y
T
 x   0
 
 Corner 
and the condition number becomes 5 !!!
- 51 -
7. Algorithm Analysis
Discarding the Corners?
 If the given signal was
sampled at 1.4142 the
Nyquist Rate, the
corners can be
removed.
 If it is not, oversampling it can be done
by FFT.




2
2

-


2

2
-
- 52 -
7. Algorithm Analysis
Error Analysis – Worst Signal
Approximation error is :
Worst error :x worst , e
2
worst
A
Polar FFT
  Arg / Max
A
 T x  ePolar FFT
Polar FFT
 TPolar x
2
2
2
x2
x
2
Worst relative error : x rworst , erworst
  Arg / Max
x
A
2
Polar FFT
 TPolar x 2
2
TPolar x 2
- 53 -
7. Algorithm Analysis
Worst Signal - Results
N=16  TC 1024256, S=Sr=St=4
USFFT
worst signal (abs.
Value) =3.469
Polar-FFT
worst signal (abs.
Value) =0.0319
The worst case
signal in the freq.
Domain (abs. and
shifted)
The worst case
signal in the freq.
Domain (abs. and
shifted)
- 54 -
7. Algorithm Analysis
Relative Worst Signal - Results
Same parameters: N=16  TC 1024256, S=Sr=St=4
USFFT
worst signal (abs.
Value) =0.0613
Polar-FFT
worst signal (abs.
Value) =0.0023
The worst case
signal in the freq.
Domain (abs. and
shifted)
The worst case
signal in the freq.
Domain (abs. and
shifted)
- 55 -
7. Algorithm Analysis
Comparing Approximations
 Solve for the eigenvalue/vectors of the matrix
A
Polar FFT
 TPolar H APolar FFT  TPolar 
and obtained k , xk k 1 ( k in ascending order).
N2
 Compare to AUSFFT by computing
2
k  A USFFT  TPolar xk 2
using the above eigenvectors and compare to k.
- 56 -
7. Algorithm Analysis
Comparing Approximations - Results
2
Mean Squared Error
10
[N=32]
0
10
USFFT
-2
10
-4
10
-6
10
Polar-FFT
-8
Eigen-space
of the PolarFFT
10
-10
10
0
200
400
600
800
1000
1200
- 57 -
Agenda
1. Thinking Polar – Continuum
2. Thinking Polar – Discrete
3. Current State-Of-The-Art
4. Our Approach - General
5. The Pseudo-Polar Fast Transform
6. From Pseudo-Polar to Polar
7. Algorithm Analysis
8.
Conclusions
- 58 -
8. Conclusions
 We have proposed a fast, accurate, stable, and
reliable Polar Discrete-Fourier-Transform.
 By this we extend utility of FFT algorithms to
new class of settings in image processing.
 Future plans:

Extend the analysis and improve further,

Demonstrate applications,

Publish the code for the procedure and some
applications over the internet.
- 59 -
Beyond this slides are
the appendix or
redundant slides
- 60 -
3. Current State-of-the-Art
USFFT for T†
 Over-sample Polar grid (and possibly
partial derivatives).
 Associate polar neighbors to each
Cartesian grid point.
 Approximate interpolation to get the
Cartesian grid values.
 Perform the Cartesian 2D Inverse-FFT.
- 61 -
3. Current State-of-the-Art
Our Reading of Literature
Direct Fourier method with over-sampling and
interpolation (also called gridding) – see
 Natterer (1985).
 Jackson, Meyer, Nishimura and Macovski (1991).
 Schomberg and Timmer (1995).
 Choi and Munson (1998).
- 62 -
A. The Fractional Fourier Transform
The Pseudo-Polar Sampling
Basically vertical lines:
NS r
1

2  2
y 

NS

r   NS r
2
NS t
2
1
y
3
2
1
0


2m
y 
x 
NSt

m  NS t
2
-1
For St=Sr=1, we have
N2 grid points
-3
-2
-3
-2
-1
0
1
2
3
x
- 63 -
A. The Fractional Fourier Transform
The Pseudo-Polar FT – Stage 1
F  x , y  


N1
N1
 
k1 0
k 2 0
f k1 ,k 2  exp ik1x  ik 2 y  
N1
N1
k1 0
k 2 0


2m
f k1 ,k 2  exp ik1
y  ik 2y  
NS t


 
N1

k1 0


2m
exp ik1
y 
NS t


N1

k 2 0
f k1 ,k 2  exp ik 2 y 
ˆf k1 , 
This part is obtained by 1D-FFT along the rows !!
- 64 -
A. The Fractional Fourier Transform
The Pseudo-Polar FT – Stage 2
N1
2y 

ˆ
F x , y  F[m, ]   f k1 ,  exp  ik1m

NSt 

k1  0


 This summation takes columns of ˆ
f k1 ,  (being equispaced 1D signals) and computes Fourier transform of it.
 The destination grid points are also 1D equi-spaced in
the frequency domain, BUT THEY ARE NOT IN THE
RANGE [-,], but rather [-y,y].
 This task is called Fractional Fourier/Chirp-Z Transform.
- 65 -
A. The Fractional Fourier Transform
Fractional Fourier Transform
N1
 2km 
F[m]   f k  exp  i
 
N


k 0
 For =1 we get the ordinary 1D-FFT,
 For =-1 we get the ordinary 1D-IFFT,
 There exists a Fast Fractional Fourier Transform with the
complexity of O(20·Nlog2N), based on 1D-FFT operations.
See: Fast fractional Fourier transforms and applications, by D. H. Bailey and P. N.
Swarztrauber, SIAM Review, 1991, and also Bluestein, Rabiner, and Rader (1960’s).
- 66 -
A. The Fractional Fourier Transform
FR-FFT Detailed
N 1
 2km 
F[m]   f k  exp  i
  
N


k 0
N1

m2
i
 N 1
e N 
f
k 2
i

e N

  k  m2  k 2  m2

  f k  exp  i
  
N


k 0

 k 
k 0
Post
Multiplication
 k  m2 
 exp  i

N


Pre-Multiplication
Convolution
[Back]
- 67 -
B. From Pseudo-Polar to Polar
Interpolation As 1D Operation
F  x , y  

N1
N1
k1 0
k 2 0
 
N1

k1 0

N1

k1 0
f k1 ,k 2  exp ik1x  ik 2y  

 m  
exp ik1 tan 
 y 
 2NS t  

N1

k 2 0
f k1 ,k 2  exp ik 2y  

 m   ˆ
exp ik1 tan 
 y  f k1 , 

 2NS t  
It is a 1D operation – But it is not the Fractional-FFT.
Can be computed by over-sampled FRFFT and interpolation.
- 68 -