Summer School at Inzell Direct Exact Inverse Pseudo-Polar

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

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 • The Fast Slant-Stack Algorithm Polar DFT: Importance and difficulties The Pseudo-Polar FFT (PPFFT): Definition and Properties Direct and Exact Inverse PPFFT (IPPFFT) Conclusions 3

Background – 1 D DFT: General Definition and Properties Matrix-vector notation Reconstruction (IDFT) 4

Background – 1 D DFT: Computability Direct evaluation of the 1 D DFT costs o(n 2) FFT – an n·log(n) DFT Algorithm 5

Example – Spectral Decomposition + + = 6

Example – Spectral Decomposition Time Domain Frequency Domain 7

Example - Denoising Time Domain Frequency Domain 8

Background – 2 D DFT – Cartesian Grid Direct evaluation → o(m 2 n 2) 9

Background – 2 D DFT 2 D FFT 1. Apply 1 D FFT for each column n times mlog(m) 2. Apply 1 D FFT for each row m times nlog(n) A total of o(Nlog(N)), N=mn Same for 2 D IFFT and for higher dimensions 10

Background – 2 D DFT Applications of rectilinear 2 D FFT Spectral Analysis Compression, Denoising Trigonometric Interpolation – Shift Property Fast Convolution/Correlation - n·log(n) instead of n 2 Convolution Theorem - 11

Background – Polar DFT Difficulties: 1 – Impossible to separate to series of 1 D FFTs 2 – Non Orthogonal (Ill Conditioned) 12

Background – Polar DFT Direct Polar DFT is impractical – o(N 2) 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 14

Pseudo-Polar FFT The Pseudo-Polar FFT (PPFFT) (Averbuch et. al. ) Basically Vertical - Concentric squares Basically Horizontal - Equally sloped lines A total of 4 n 2 grid points 15

Pseudo-Polar FFT The Pseudo Polar FFT Fractional FFT – o(nlog(n)) 16

Fractional FFT Algorithm (D. Bailey and P. Swarztrauber 1990) 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 2 n by 2 n circulant Matrix C as follows: 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 This procedure is very similar to the FFFT Algorithm! 18

FFFT and Structured Matrices FFFT in Matrix-vector notation What is the structure of V ? V is symmetric Vandermonde 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= β-2 jk Define Dj= βjk and T=DVD 19

Pseudo-Polar FFT The PPFFT Algorithm 1. 1 D FFT for each 0 -padded column n times 2 nlog(2 n) 2. Apply Fractional FFT for each row With α=l/n 2 n times nlog(n) A total of o(Nlog(N)), N=n 2 Repeat the same procedure for the transposed image matrix to compute the BH coefficients 20

Pseudo-Polar FFT The PPFFT – Matrix notation 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 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 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: 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 3 D Pseudo Spherical FFT and Radon 24

Direct IPPFFT Direct Inverse PPFFT Theorem: There exists W, a Real Positive Diagonal Matrix such that: If y belongs to the image of A then: The elements of W can be rapidly pre-calculated for any given n 25

Direct IPPFFT Problem formulation - 4 n 4 constraints - 4 n 2 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 1. Experiments showed that the over-determined system is solvable for n≤ 8 System Reduction 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(n 2 log(n)) operations 27

Direct IPPFFT The ideal Weights 28

Direct IPPFFT Fast solver for the ideal weights Define: PPFFTos 1, os 2 = PPFFT with over-sampling ration = os 1·os 2 os 1 n slopes and os 2 n radiuses The Conventional PPFFT can be denoted as PPFFT 2, 2 Define: If x exists it satisfies 29

Direct IPPFFT Existence of ideal weights The system has a solution if is invertible There is no solution for PPFFT 2, 2 in its standard form since 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 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 Backprojections of SS coefficients Dl(t) is an interpolating kernel 36

The Fast Slant-Stack transform Projection-Slice Theorem For a given n by n discrete image 1. Compute the PPFFT coefficients 2. Apply 1 D IFFT to each vector of “same slope” coefficients in the PPFFT coefficients array Slant-Stack Lines 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