Transcript The Pseudo-Polar FFT (PPFFT)

Summer School at Inzell
Direct Exact Inverse Pseudo-Polar
FFT and Radon Transform Using
Orthogonalizing Weights
Ofer Levi
Department of Industrial Engineering and Management
Ben-Gurion University of the Negev
Beer Sheva, Israel
On sabbatical leave at the Institute for Computational and
Mathematical Engineering, Stanford University
April 8, 2015
My research interests
Mathematical and statistical modeling of
physical and biological processes and systems
 Inverse problems solving using Sparse
Representation and Compressed Sensing
 High dimensional image analysis and
processing
 Large scale Optimization, Parallel and
Distributed Computing
 Numerical Analysis and Matrix Computation

Lecture Structure

Background
• Rectilinear DFT: General Definition and Properties
• Polar DFT: Importance and difficulties
• The Pseudo-Polar FFT (PPFFT): Definition and
Properties
• The Fast Slant-Stack Algorithm

Direct and Exact Inverse PPFFT (IPPFFT)

Conclusions
3
Background – 1D DFT
1D DFT: General Definition and Properties
F : C n  C n , F ( x)  xˆ
xˆ k  1n
e ix  cos( x )  i sin( x )
n 1
 i 2 jk / n
x
e
, k  [0, n  1]
 j
j 0
Matrix-vector notation
xˆ  Ax, Ak j 
1
e i 2 jk / n
n
A* A  I , Ak j  A jk
Reconstruction (IDFT)
 xk 
1
n
n 1
i 2 jk / n
ˆ
x
e
,
 j
j  [ 0, n  1]
k 0
4
Background – 1D DFT
1D DFT: Computability
Direct evaluation of the 1D DFT costs o(n2)
n 1
xˆ k  1n  x j e i 2 jk / n , k  [0, n  1]
j 0
FFT – an n·log(n) DFT Algorithm
w j  x2 j
z j  x2 j 1
 n 
j  0,  1
 2 
 i 2 k / n
xˆ k  wˆ k  e
zˆk
 i 2 k / n
 ˆ
x n  wˆ k  e
zˆk
k
 n 
 1
 2 
k   0,
2
5
Example – Spectral Decomposition
+
+
=
6
Example – Spectral Decomposition
Time Domain
Frequency Domain
7
Example - Denoising
Time Domain
Frequency Domain
8
Background – 2D DFT
2D DFT – Cartesian Grid
F : C mn  C mn , F ( x)  xˆˆ
xˆˆ k1k2
1

mn
m 1 n 1
i 2 ( j1k1 / m  j2 k 2 / n )
x
e
,
  j1 j2
j2  0 j1  0
k1  [ 0 , m  1] , k 2  [ 0 , n  1]
Direct evaluation → o(m2n2)
xˆˆ k1k 2
1

m
 1



j2  0  n
m 1
n 1
x
j1  0
j1 j 2
e
 i 2 ( j1k1 / m )
  i 2 ( j 2 k 2 / m )
e


9
Background – 2D DFT
2D FFT
1. Apply 1D FFT for each column
n times mlog(m)
2. Apply 1D FFT for each row
m times nlog(n)
A total of o(Nlog(N)), N=mn
Same for 2D IFFT and for higher dimensions
10
Background – 2D DFT
Applications of rectilinear 2D FFT

Spectral Analysis

Compression, Denoising

Trigonometric Interpolation – Shift Property

Fast Convolution/Correlation - n·log(n) instead of n2
Convolution Theorem -
Z  X * Y  Zˆ  Xˆ  Yˆ
11
Background – Polar DFT
Polar DFT
xˆˆ k1k 2
1

n
n 1 n 1
x
j2  0 j1  0
j1 j2
e
 i 2 ( j1k1  j2 k 2 ) / n
,
k 1  r cos( ) , k 2  r sin( ) , r  [  n / 2 , n / 2  1],   [ 0 ,  ]
ˆxˆ  1
r
n
n 1 n 1
x
j2  0 j1  0
j1 j 2
e
 i 2 r( j1cos( )  j 2 sin( )) / n
Difficulties:
1 – Impossible to separate to series of 1D FFTs
2 – Non Orthogonal (Ill Conditioned)
12
Background – Polar DFT
Polar DFT
Direct Polar DFT is impractical –
o(N2) and no direct inverse
 Common Solution – Interpolation to and from
Cartesian Grid with Oversampling


Trade off between time and accuracy
13
Background – Polar DFT
Importance of Polar DFT
Accurate Rotation (Shift in Polar Coordinates)
 Rigid body Registration
 MRI – Reconstruction from Polar Grid
 CT – Reconstruction from Projections

Projection-Slice Theorem
R ( , r )    f ( x, y ) r  x cos( )  y sin( ) dxdy
y x
ˆˆ
R (r )  F f ( r cos( ), r sin( ))
1
1
14
Pseudo-Polar FFT
The Pseudo-Polar FFT (PPFFT)
(Averbuch et. al.)
Basically Vertical
- Concentric squares
Basically
Horizontal
- Equally sloped lines
A total of 4n2 grid points

ml
 nn nn   
ll
ml

BV  k12  ;;ll[[nn, ,nn11]], ,kk21 ; ;mm , , 11 
BH
22
nn
 22 22   

15
Pseudo-Polar FFT
The Pseudo Polar FFT
ˆxˆ BH
lm
BH
xˆˆlm
1

n
1

n
n 1 n 1
x
j 2  0 j1  0
 1



j2  0  n
n 1
j1 j 2
e
n 1
x
j1  0
j1 j 2
e
 i l ( j1n  2 j 2 m ) / n 2
 i 2 j1l /( 2 n )
 i 2 j2 m / n
l
e
; 

n

Fractional FFT – o(nlog(n))
xˆ k 
1
n
n 1
 n , n  1
 i 2 jk / n
x
e
,
k


j
 2 2 
j 0
16
Fractional FFT Algorithm
(D. Bailey and P. Swarztrauber 1990)
n 1
xˆ k   x j e i 2 jk / n ; k  [0, n  1]
j 0
 x j ei j  / n ; j  [0, n  1]
1. y j  

0
; j  [ n,2n  1]
2
 e i j  / n ; j  [0, n  1]
2. z j   i ( j  2 n ) 2  / n
e
; j  [ n,2n  1]
2
3. wk 

2 n 1
y z
j 0
4. xˆk  e
j k j
ik 2 / n
; k  [0, n  1]
wk ; k  [0, n  1]
17
Some basic facts about Toeplitz Matrices
T has a Toeplitz structure if Tjk=f(j-k), i.e. T has constant diag’s
Any n by n Toeplitz matrix T can be expanded into a 2n by 2n
circulant Matrix C as follows:
T S 
C

S
T


Where S is also Toeplitz
A circulant Matrix C can be diagonalyzed using the DFT
Matrix F as follows: C=FDF-1, D=Diag(v) where v is the
Fourier transform of the first column of C
Tx can be computed in nlog(n)
flops by doing the following
T S   x  Tx 
 
    

 S T  0   Sx 
This procedure is very similar to the FFFT Algorithm!
18
FFFT and Structured Matrices
FFFT in Matrix-vector notation xˆ

 Vx, Vk j 
1
e i 2 jk / n
n
V is symmetric Vandermonde
What is the structure of V ?
Reminder: if V=V(a) then Vjk=ajk
V=Vt => Vjk=ajk
Theorem : A symmetric Vandermonde Matrix V can be
decomposed as V=DTD where T is Toeplitz
Proof: If V is symmetric Vandermonde then there
exist a unique scalar β such that Vjk= β-2jk
Define Dj= βjk and T=DVD
 T jk  
( j k ) 2
19
Pseudo-Polar FFT
The PPFFT Algorithm
1. 1D FFT for each 0-padded column
n times 2nlog(2n)
2. Apply Fractional FFT for each row
With α=l/n 2n times nlog(n)
A total of o(Nlog(N)), N=n2
Repeat the same procedure for the transposed image
matrix to compute the BH coefficients
20
Pseudo-Polar FFT
The PPFFT – Matrix notation
A  C 4n
Alm, j1 j2 
BH
1
n
e
2
n
2
,
 il ( j1n  2 j2 m ) / n 2
 A BH 
A   BV 
A 
, Alm, j1 j2 
BV
1
n
e
 il ( 2 j1m  j2 n ) / n 2

A can be implicitly applied in O(Nlog(N)) operations

Denote the Adjoint PPFFT by A*
A* can be also implicitly applied in o(Nlog(N))
21
Pseudo-Polar FFT
Inverse PPFFT
x  arg min ( Au  y 2 )
u
Solve A* Ax  A* y
Use CGLS or LSQR for the Normal Equations
A problem: A is ill conditioned, k(A) is proportional to n
Solution: Solve
A*W 2 A~x  A * W 2 y
W is diagonal when each diagonal element is the grid
point radius of the corresponding PPFFT coefficient
If a zero residual solution exists then
~x  x
22
Pseudo-Polar FFT
Weighted PPFFT
- Each coefficient is multiplied
by its grid point radius
- The weights compensates for
the non-uniform grid sampling
- Experimental result:
K ( A*W 2 A)  1.2
The weighted IPPFFT converges within 4-5 iterations
23
Pseudo-Polar FFT
Applications of the PPFFT





Fast and Accurate interpolation to PFT
Fast and Accurate interpolation to Spiral FT
Fast Slant-stack – An nlog(n) Radon Transform
Fast and Accurate Rotations
3D Pseudo Spherical FFT and Radon
24
Direct IPPFFT
Direct Inverse PPFFT
Theorem:
There exists W, a Real Positive Diagonal Matrix
such that:
A*W 2 A  I
If y belongs to the image of A then:
x AW y
*
2
The elements of W can be rapidly pre-calculated for any given n
25
Direct IPPFFT
Problem formulation
A W A  I, W  R
*
2
4 n 2 4 n 2
, A R
4 n 2 n 2
- 4n4 constraints
- 4n2 variables (W is diagonal)
The problem is over-determined
There is no solution for an arbitrary A
and a diagonal W
26
Direct IPPFFT
Finding the Ideal weights
A*W 2 A  I  Mx  b
M C
4 n 4 4 n 2
, x  diag (W 2 ), b  vec( I )
1. Experiments showed that the over-determined system is
solvable for n≤8
~ ~
Mx  b  Mx  b
System Reduction
~
2 n 2 4 n 2 ~
4n2
M C
,b R
2. The under-determined reduced system is solvable.
The weights could be computed numerically for n≤32
3. The reduced system could be solved by a fast
iterative FFT based solver within o(n2log(n)) operations
27
Direct IPPFFT
The ideal Weights
28
Direct IPPFFT
Fast solver for the ideal weights
Define: PPFFTos1,os2 = PPFFT with over-sampling ration = os1·os2
os1n slopes and os2n radiuses
Aos ,os : C nn  C os1 nos2 n
1
2
The Conventional PPFFT can be denoted as PPFFT2,2
xos1 ,os2  diag (Wos21 ,os2 )
Define:
If x exists it satisfies -

n2  n
2
~
~ 1 ; j 
b  Rn , bj  
2
0 ; otherwise

*
os1 / 2 , os2 / 2 os1 , os2
A
x
~
b
29
Direct IPPFFT
Existence of ideal weights
The system
*
os1 / 2 ,os2 / 2
A
*
os1 / 2 , os2 / 2 os1 , os2
A
x
~
b
has a solution if
is invertible
There is no solution for PPFFT2,2 in its standard form
*
since A1,1 is singular
For a slightly modify PP grid solution exists. It can be
verified using the Vandermonde structure of A and the
fact that the modified grid has distict set of points
30
Direct IPPFFT
Modified PP Grid
31
Error analysis
x1  arg min ( Au  y 2 )
u
x2  arg min ( W ( Au  y ) 2 )
u
y  ~y  r , Ax  ~y
x1  ( A* A) 1 A* y  ( A* A) 1 A* ~y  ( A* A) 1 A*r 
( A* A) 1 A* Ax  ( A* A) 1 A*r  x  ( A* A) 1 A*r
x2  A*W 2 y  A*W 2 ~y  A*W 2 r  A*W 2 Ax  A*W 2 r  x  A*W 2 r
A*W 2  ( A* A) 1 A*

x  x2  x  x1
32
Conclusions

Matrix approach can be valuable for better
understanding and analysis of discrete signal and image
transformations

The Pseudo-Polar Fourier Transform (PPFFT) posses
several attractive computational properties

The PPFFT combines the best properties of both the
Cartesian FT and the Polar FT

Can be generalized to higher dimensions

Provides fast and accurate discrete Radon Transform

Direct inverse
33
Thank You!
34
The Slow Slant-Stack Transform
1. Zero-pad the image from the left and the right (BV lines)
2. Shear the padded image in a θ angle using trigonometric
interpolation via FFT
3. Sum the sheared image array column-wise to get the θ-projection
35
The Fast Slant-Stack Algorithm
n / 2 1
~
SS ( y  sx  t , I )   I (u , su  t )
u  n / 2
~
I (u , y ) 
n / 2 1
 I (u, v) D
2n
u  n / 2
Dl (t )  e
i t / l

( y  v)
Backprojections of
SS coefficients
sin( t )
l  sin( t / l )
Dl(t) is an interpolating kernel
~
 n n 
I (u , v )  I (u , v ) , u , v    ,  1
 2 2 
36
The Fast Slant-Stack transform
Projection-Slice Theorem
For a given n by n discrete image
1. Compute the PPFFT coefficients
2. Apply 1D IFFT to each vector of
“same slope” coefficients in the
PPFFT coefficients array
Slant-Stack Lines
2
2


BH   y  sx  t ; s  1,1  ,  ,1  , t   n, n  1,  , n  1
n
n


2
2


BV   x  sy  t ; s  1,1  ,  ,1  , t   n, n  1,  , n  1
n
n


Uniform horizontal
/vertical spacing in
each projection !
37
Conclusions
- A direct solution for an important inverse problem
- A direct result – Direct inverse for the Fast Slant-Stack
- Can be generalized to higher dimensions
- A new methodology for LS problem – pre-computation
of weights
- Various Applications
38