The Fast Octa-Polar Fourier Transform and its expansion to an accurate discrete radon transform

1
The Fast Octa-Polar Fourier Transform and its
expansion to an accurate discrete radon transform

• Ofer Levi

Department of Industrial Engineering and
Management
Ofir Harari and Ronen Peretz
Department of Mathematics Ben-Gurion University
of the Negev Beer Sheva, Israel
January 12, 2014
2
Outline

• Main result and the relation to Radon Transform
• Rectilinear DFT General Definition and
  Properties
• The Pseudo-Polar FFT (PPFFT)
• The FFFT and its relation to structured matrices
• The Octa-Polar FFT
• Summary and on-going work

3
2D Fast Fourier Transforms
The computational cost for a single 2D Fourier
coefficient is O(n2)
Are there any other 2D meshes that allow a fast
simultaneous computation ?
4
Main Result a new, almost-polar FFTThe
Octa-Polar FFT
5
Importance of Polar DFT
Background Polar DFT

• Accurate Rotation (Shift in Polar Coordinates)

• MRI Reconstruction from Polar Grid

• CT Reconstruction from Projections

6
Relation between Radon and Fourier Transforms
An easy exercise !
7
Relation between Radon and Fourier Transforms
Serious problems - Exact Polar DFT takes O(n4)
flops! The
Transform is highly ill-conditioned
8
Approximated Polar DFTs
Fast but inaccurate and require non sequential
memory access
9
Approximated Polar DFTs
Improved PP-FFT based interpolation (Elad et. Al.)
10
1D DFT General Definition and Properties
Background 1D DFT
Matrix-vector notation
Reconstruction (IDFT)
11
1D DFT Computability
Background 1D DFT

• Direct evaluation of the 1D DFT costs o(n2)

FFT an nlog(n) DFT Algorithm
12
Example Spectral Decomposition
13
Example Spectral Decomposition
14
Example - Denoising
15
2D DFT Cartesian Grid
Background 2D DFT
Direct evaluation ? o(m2n2)
16
2D complex exponents
17
2D FFT
Background 2D DFT

• 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)), Nmn
Same for 2D IFFT and for higher dimensions
18
Applications of rectilinear 2D FFT
Background 2D DFT

• Spectral Analysis

• Compression, Denoising

• Trigonometric Interpolation Shift Property

• Fast Convolution/Correlation - nlog(n) instead
  of n2

19
Polar DFT
Background Polar DFT
Difficulties 1 Impossible to separate to
series of 1D FFTs 2 Non Orthogonal (Ill
Conditioned)
20
Polar DFT
Background Polar DFT

• Direct Polar DFT is impractical
• o(n4) and no direct inverse

• Common Solution Interpolation to and from
  Cartesian Grid with Oversampling

• Trade off between time and accuracy

21
The Pseudo-Polar FFT (PPFFT)(Donoho et. al.)
Pseudo-Polar FFT
Basically Vertical

• Concentric squares
• Equally sloped lines

A total of 4n2 grid points
22
The Pseudo Polar FFT
Pseudo-Polar FFT
23
Fractional FFT Algorithm (D. Bailey and P.
Swarztrauber 1990)
24
Some basic facts about Toeplitz Matrices
T has a Toeplitz structure if Tjkf(j-k), i.e. T
has constant diags
A circulant Matrix C can be diagonalyzed using
the DFT Matrix F as follows CFDF-1, DDiag(v)
where v is the Fourier transform of the first
column of C
This procedure is very similar to the FFFT
Algorithm!
25
FFFT and Structured Matrices
V is symmetric Vandermonde
What is the structure of V ?
Reminder if VV(a) then Vjkajk VVt gt Vjkajk
Theorem A symmetric Vandermonde Matrix V can be
decomposed as VDTD where T is Toeplitz
Proof If V is symmetric Vandermonde then there
exist a unique scalar ß such that Vjk ß-2jk
Define Dj ßjk and TDVD
26
The PPFFT Algorithm
Pseudo-Polar FFT

• 1D FFT for each 0-padded column
• n times 2nlog(2n)

2. Apply Fractional FFT for each row With al/n
2n times nlog(n)
A total of o(Nlog(N)), Nn2
Repeat the same procedure for the transposed
image matrix to compute the BH coefficients
27
The PPFFT Matrix notation
Pseudo-Polar FFT

• 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))

28
Inverse PPFFT
Pseudo-Polar FFT
Use CGLS or LSQR for the Normal Equations
A problem A is ill conditioned, k(A) is
proportional to n
Solution Solve
W is diagonal when each diagonal element is the
grid point radius of the corresponding PPFFT
coefficient
If a zero residual solution exists then
29
Weighted PPFFT
Pseudo-Polar FFT

• Each coefficient is multiplied by its grid
  point radius

• The weights compensates for the non-uniform
  grid sampling

The weighted IPPFFT converges within 4-5
iterations
30
The Slow Slant-Stack Transform
31
The Fast Slant-Stack Algorithm
32
The Fast Slant-Stack transform

• For a given n by n discrete image
• Compute the PPFFT coefficients
• Apply 1D IFFT to each vector of same slope
  coefficients in the PPFFT coefficients array

Uniform horizontal /vertical spacing in each
projection !
33

• The Octa-Polar FFT

Rays slopes are sampled uniformly, and points
are equi-spaced along each ray.
N/S and E/W are treated similarly to BH BV in the
PPFFT
34
Treating The NW/SW and NE/SE grid points sets
1. 45 degrees rotation by reordering and embading
in a big zeros matrix
2. Applying rectangular FFFT both vertically and
horizontally
After an appropriate change of variables to the
indices
35
Summary and Future research

• The Octa-Polar grid a new almost-polar exact FFT
  
• Can be expanded to a new efficient DRT
• Provides much better approximation to PFT by
  using Octagons instead of Squres
• - Expansion to 3D
• Finding other non-uniform meshes for which fast
  algorithm exist
• Preconditioning the inverse solver
• Error Analysis
• Testing on real raw data