2D and 3D Fourier Based Discrete Radon Transform Amir Averbuch With Ronald Coifman – Yale University Dave Donoho – Stanford University Moshe Israeli – Technion, Israel Yoel Shkolnisky – Yale University

Research Activities • Polar processing (Radon, MRI, diffraction tomography, polar processing, image processing) • Dimensionality reduction (hyperspectral processing, segmentation and sub-pixel segmentation, remote sensing, performance monitoring, data mining) • Wavelet and frames (error correction, compression) • Scientific computation (prolate spheroidal wave functions) • XML (fast Xpath, handheld devices, compression) • Nano technology (modeling nano batteries, controlled drug release, material science simulations) interdisciplinary research with material science, medicine, biochemistry, life sciences

Participants: previous and current • Dr Yosi Keller – Gibbs Professor, Yale • Yoel Shkolnisky - Gibbs Professor, Yale • Tamir Cohen – submitted his Ph.D

• Shachar Harrusi – Ph.D student • Ilya Sedelnikov - Ph.D student • Neta Rabin - Ph.D student • Alon Shekler - Ph.D student • Yossi Zlotnick - Ph.D student • Nezer Zaidenberg - Ph.D student • Zur Izhakian – Ph.D student

Computerized tomography

CT-Basics

CT-Basics

Typical CT Images 7

CT Scanner 8

CT-Basics " Projection " : 

along f line

( ρ )

d

ρ

y Ra X ay X -R a y X -R a y X -R ay X-R ay X-R X -R y a X -R a y X -R ay X Ra y X-Ra y

CT-Basics

Introduction – CT Scanning • X-Ray from source to detector • Each ray reflects the total absorption along the ray • Produce projection from many angles • Reconstruct original image

2D Continuous Radon Transform The 2D continuous Radon is defined as

Rf

(  ,

s

)  

L f

(

x

,

y

)

du

        

f

(

x

,

y

)  (

x

cos  

y

sin  

s

)

dxdy

g(s,  s y 0 s  x f(x,y)

2D Continuous Fourier Slice Theorem 2D Fourier slice theorem where

f

ˆ

F

( (  1

Rf

,  )(  2 ) ,    )        

f f

ˆ (  (

x

, cos  ,  sin is the 2D continuous Fourier transform of

f.

 )

y

)

e

 2 

i

(

x

 1 

y

 2 )

dxdy

The1D Fourier transform with respect to

s

 equal to a central slice, at angle

θ

, of the 2D Fourier transform of the function

f(x,y)

.

Discretization Guidelines We will look for both 2D and 3D definitions of the discrete Radon transform with the following properties: • Algebraic exactness • Geometric fidelity • Rapid computation algorithm • Invertibility • Parallels with continuum theory

2D Discrete Radon Transform Definition y  L x • Summation along straight lines with θ<45° • Trigonometric interpolation at non grid points

2D Discrete Radon Transform – Formal Definition For a line

y



sx



t

(

s

 1 ) we define

Radon

({

y



sx



t

},

I

) 

u n

/    2 

n

/ 1 ~ 1

I

2 (

u

,

su



t

) where and

m



D m

2

n

(

t

)  1 .

I

~ 1  (

u m

,

y

) 

n

/ 2   1 sin( sin(  

v t t

 / ) 

n

/

m

) 2

I

(

u

,

v

)

D m

(

y



v

) is the Dirichlet kernel with

2D Definition - Illustration y

Radon

({

y



sx



t

},

I

) 

u n

/    2 

n

/ 1

I

~ 1 2 (

u

,

su



t

)

I

~ 1 (

u

,

y

) 

v n

/   2   1

I n

/ 2 (

u

,

v

)

D m

(

y



v

)

D m

(

t

) 

m

sin( sin(  

t t

/ )

m

)

m

 2

n

 1 L x

2D Discrete Radon Definition – Cont.

For a line

x



sy



t

(

s

 1 ) we define

Radon

({

x



sy



t

},

I

) 

v n

/   2  

n

/ 1 ~

I

2 2 (

sv



t

,

v

) where and

I

~ 2 (

x

,

v

) 

u n

/   2   1

n

/

I

2 (

u

,

v

)

D m

(

x



u

)

D m

(

t

) 

m

sin( sin(  

t t

/ )

m

)

m

 2

n

 1

2D Definition - Illustration y

Radon

({

x



sy



t

},

I

) 

v n

/   2  

n

/ 1 ~

I

2 2 (

sv



t

,

v

)

I

~ 2 (

x

,

v

) 

u n

/   2   1

I n

/ 2 (

u

,

v

)

D m

(

x



u

)

D m

(

t

) 

m

sin( sin(  

t t

/ )

m

)

m

 2

n

 1 L x

Selection of the Parameter

t t=n-1

y x • •

Radon({y=sx+t},I)

• Sum over all lines with non trivial projections.



n



t



n t=-n

• Same arguments for basically vertical lines.

Selection of the Parameter

m

• Periodic interpolation kernel.

• Points out of the grid are interpolated as points inside the grid.

• Summation over broken line.

• Wraparound effect.

Selection of the Parameter

m

– Cont.

• Pad the image prior to using trigonometric interpolation.

• Equivalent to elongating the kernel.

• No wraparound over true samples of

I.

• • Summation over true geometric lines.

Required: m

 2

n

 1

The Translation Operator The translation operator (

T

  )

u

translates the vector  

i n

  

n



i D m

(

u



i

  ) using trigonometric interpolation.

Example: translation of a vector with   1    T 

The Shearing Operator For the slope

θ

of a basically horizontal line: (

S

1 

I

)(

u

,

v

)  (

T



u

tan 

I

(

u

,  ))

v

For the slope

θ

of a basically vertical line: (

S

 2

I

)(

u

,

v

)  (

T



v

cot 

I

(  ,

v

))

u

Motivation: The shearing operator translates the samples along an inclined line into samples along horizontal/vertical line.

I The Shearing Operator Illustration y y x I x

I The Shearing Operator Illustration y y x I x

Alternative Definition of the Discrete Radon Transform

Radon

({

y



sx



t

},

I

) 

u n

/   2   1

n

/ ( 2

S

 1

I

 1 )(

u

,

t

)

Radon

({

x



sy



t

},

I

) 

v n

/   2   1

n

/ ( 2

S

 2

I

 2 )(

t

,

v

)

I

 1

I

 2 the x-axis respectively

I

2D Discrete Fourier Slice Theorem Using the alternative discrete Radon definition we prove:

F

(

R



I

)(

k

)    

I

ˆ

I

ˆ (  (

k

,

s

1

k



s

,

k

2

k

)

s

1 )

s

2   tan   cot    [    [  / / 4 ,  4 , 3  / 4 ] / 4 ] where

I

ˆ (  1 ,  2 ) 

u n

/   2  1

n

/ 2    

n

/ 2

v n

 1 /

I

2 (

u

,

v

)

e

 2 

i

(  1

u

  2

v

) /

m m

 2

n

 1

Discretization of θ • The discrete Radon transform was defined for a continuous set of angles.

• For the discrete set Θ  1  2     / 2   arctan  arctan 2

l

/

n

  2

l

|

l

/ 

n

 |

l

  

n

/   2 

n l

/  2 

n

/

l

2  

n

/ 2     1   2 the discrete Radon transform is discrete in both Θ and

t

.

• For the set Θ, the Radon transform is rapidly computable and invertible.

Θ 2 y Illustration of Θ y x x Θ 1

Fourier Slice Theorem Revisited For  For   

F F

arctan(  ( (

R



I

2

l

)(

k

) / 

n

),

I

ˆ (    2

lk

 2 /

n

,

k

) 

PPI

1 /

R

 2

I

 )( arctan(

k

) 

I

ˆ ( 2

k l

/

n

), ,  2

lk

/ 

n

 )  1 

PPI

2 (

k

, (

k

,

l

)

l

) where

I

ˆ (  1 ,  2 ) 

u n

 /  2

n

 1 / 2

v n

/   2  1  

n

/

I

2 (

u

,

v

)

e

 2 

i

(  1

u

  2

v

) /

m

.

We define the pseudo-polar Fourier transform:

PPI

1 (

k

,

l

)

PPI

2 (

k

,

l

)  

I

ˆ (  2

lk

/

I

ˆ (

k

,  2

lk n

,

k

) /

n

)

The Pseudo-Polar Grid

PPI

1 (

k

,

l

)

PPI

2 (

k

,

l

)  

I

ˆ (  2

lk

/

I

ˆ (

k

,  2

lk n

,

k

) /

n

) special pointset called the pseudo-polar grid.

The pseudo-polar grid is defined by

P

1     2

lk

/

n

,

k

 | 

n

/ 2 

l



n

/ 2 , 

n



k



n



P

2   

k

,  2

lk

/

n

 | 

n

/ 2 

l



n

/ 2 , 

n



k



n



P



P

1 

P

2

The Pseudo-Polar Grid Illustration y y x x

P

2

P

1

The Pseudo-Polar Grid Illustration y x

The Fractional Fourier Transform The fractional Fourier transform is defined as (

F



n

 1

X

)(  ) 

u n

/    2

n

/ 2

X

(

u

)

e

 2 

i



u

/(

n

 1 )    , 

n

/ 2   Can be computed for any  using

O

(

n

log

n

) 

n

/ 2 .

operations.

We can use the fractional Fourier transform to compute samples of the Fourier transform at any spacing.

Resampling in the Frequency Domain resampling operator

G k

,

n



F m

 

F n

 1 1 0.8

0.6

0.4

0.2

0 -0.2

-0.4

-0.6

-0.8

-1 -3 -2 -1 0 1 2 3

G k,n

1 0.8

0.6

0.4

0.2

0 -0.2

-0.4

-0.6

-0.8

-1 -3 -2 -1 0 1 2 3

2D Discrete Radon Algorithm (

i

1 , 2 ) operations by using 1D Fourier transform.

O

(

n

log

n

)  y  y  x G k,n  x

2D Discrete Radon Algorithm – Cont.

Description: (

PP 1 I

) 1. Pad both ends of the y-direction of the image

I

and compute the 2D DFT of the padded image. The results are placed in

I’

.

2. Resample each row

k

with

α = 2k/n

.

in

I’

using the operator

G k,n

3. Flip each row around its center.

Papers http://www.math.tau.ac.il/~amir http://pantheon.yale.edu/~yk253/ • • • • • • • • • • • • • • • • • • Optical Snow Analysis using the 3D-Xray Transform, submitted.

Fast and Accurate Polar Fourier Transform, submitted.

Discrete diffraction tomography, submitted.

2D Fourier Based Discrete Radon Transform, submitted.

Algebraically accurate 3-D rigid registration, IEEE Trans. on Signal Proessing.

Algebraically Accurate Volume Registration using Euler's Theorem and the 3-D Pseudo-Polar FFT, submitted.

Fast Slant Stack: A notion of Radon Transform for Data in a Cartesian Grid which is Rapidly Computible, Algebraically Exact, Geometrically Faithful and Invertible, SIAM Scientific Computing.

Pseudo-polar based estimation of large translations, rotations and scalings in images, IEEE Trans. on Image Processing.

The Angular Difference Function and its application to Image Registration, IEEE PAMI. 3D Discrete X-Ray Transform, Applied and Computational Harmonic Analysis 3D Fourier Based Discrete Radon Transform, Applied and Computational Harmonic Analysis Digital Implementation of Ridgelet Packets, Beyond wavelets – chapter in book.

Multidimensional discrete Radon transform, chapter in book.

The pseudopolar FFT and its Applications, Research Report A signal processing approach to symmetry detection, IEEE Trans. on Image Processing Fast and accurate pseudo-polar protein docking, submitted.