Fourier Reconstruction from Non-Uniform Spectral Data Anne Gelb - - PowerPoint PPT Presentation

Fourier Reconstruction from Non-Uniform Spectral Data

Anne Gelb School of Mathematical and Statistical Sciences Arizona State University anne.gelb@asu.edu

February Fourier Talks February 17 2011 Research supported in part by National Science Foundation grants DMS 0510813 and DMS 0652833 (FRG).

Introduction Non-Uniform Data Current Methods Alternate Approaches

1 Introduction

Magnetic Resonance Imaging Sampling Patterns in MR Imaging Challenges in Cartesian Reconstruction

Spectral Reprojection

2 The Non-Uniform Data Problem

Problem Formulation The Non-harmonic Kernel Reconstruction Results using the Non-harmonic Kernel

3 Current Methods

Reconstruction Methods Error Characteristics

4 Alternate Approaches

Spectral reprojection Incorporating Edge Information Spectral reprojection for Fourier Frames

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

MR Imaging Process

An acquired MR signal can be written as (2D slice) S(kx, ky) = Z Z ρ(x, y) e−i(kx x+ky y) dx dy ρ(x, y) is a measure of the concentration of MR relevant nuclei (spins) and kx = γ 2π Z t Gx(τ)dτ, ky = γ 2π Z t Gy(τ)dτ

Gx and Gy are the applied gradient fields We denote the signal acquisition space k = (kx, ky) as “k-space”

20 40 60 80 100 120 50 100 0.2 0.4 0.6 0.8 1 ωkx Fourier Coefficients ωky |ˆ f (ωkx, ωky)|

(a) Fourier transform

−10 −8 −6 −4 −2 2 4 6 8 10 −10 −8 −6 −4 −2 2 4 6 8 10 kx ky

(b) k-space sampling trajectory

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

Sampling Patterns

−10 −8 −6 −4 −2 2 4 6 8 10 −10 −8 −6 −4 −2 2 4 6 8 10 kx ky

(c) Cartesian Sampling

−0.5 −0.4 −0.3 −0.2 −0.1 0.1 0.2 0.3 0.4 0.5 −0.5 −0.4 −0.3 −0.2 −0.1 0.1 0.2 0.3 0.4 0.5 kx ky Spiral scan trajectory

(d) Non-Cartesian Sampling – Spiral Imaging

Figure: MR Imaging Sampling Patterns1

1Spiral sampling pattern courtesy Dr. Jim Pipe, Barrow Neurological Institute, Phoenix,

Arizona

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

Advantages and Disadvantages of each Pattern

Cartesian Imaging Advantages Well understood and often used in equipment Simple reconstruction procedure Computationally efficient through use

• f the standard FFT

Disadvantages Susceptible to undersampling – causes aliased images Gradient field waveforms can mean sharp transitions in collection Non-Cartesian Imaging Advantages Less susceptible to aliasing artifacts – aliased images are of “diagnostic quality” Easier to generate gradient field waveforms Disadvantages Complex reconstruction procedure (if possible!) Computationally demanding – FFT algorithm no longer applies

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

Advantages and Disadvantages of each Pattern

Cartesian Imaging Advantages Well understood and often used in equipment Simple reconstruction procedure Computationally efficient through use

• f the standard FFT

Disadvantages Susceptible to undersampling – causes aliased images Gradient field waveforms can mean sharp transitions in collection Non-Cartesian Imaging Advantages Less susceptible to aliasing artifacts – aliased images are of “diagnostic quality” Easier to generate gradient field waveforms Disadvantages Complex reconstruction procedure (if possible!) Computationally demanding – FFT algorithm no longer applies

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

Challenges in Cartesian Reconstruction

The Gibbs Phenomenon affects the reconstruction of piecewise-smooth functions.

• ccurs when data is sampled uniformly or non-uniformly.

Two important consequences:

Gibbs ringing artifact – presence of non-physical oscillations in the vicinity

• f discontinuities.

Reduced rate of convergence to first order even in smooth regions of the reconstruction.

Why is this important?

Oscillations cause post-processing problems in tasks like segmentation, edge detection etc. Filtering oscillations may cause loss of important information. The reduced order of convergence means more Fourier coefficients are needed to obtain a good reconstruction.

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

The Gibbs Phenomenon – An Example

SNf(x) =

N

X

k=−N

ˆ f(k)eikx, ˆ f(k) = 1 2π Z π

−π

f(x)e−ikxdx

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −1 −0.5 0.5 1 1.5 2 x SNf(x) Fourier Reconstruction True Fourier

(a) Reconstruction

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −7 −6 −5 −4 −3 −2 −1 x Log10 |e| Log Error in Reconstruction

(b) Reconstruction error

Figure: Gibbs Phenomenon, N = 32

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

The Gibbs Phenomenon – An Example

SNf(x) =

N

X

k=−N

ˆ f(k)eikx, ˆ f(k) = 1 2π Z π

−π

f(x)e−ikxdx

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −1 −0.5 0.5 1 1.5 2 x SNf(x) Fourier Reconstruction True Fourier

(a) Reconstruction

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −7 −6 −5 −4 −3 −2 −1 x Log10 |e| Log Error in Reconstruction N=32 N=64 N=256

(b) Error does not go away by increasing N

Figure: Gibbs Phenomenon

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

Shepp Logan MRI

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

Reconstruction from Fourier Data by Spectral Reprojection (Gottlieb, Shu

• et. al.)

Assume that Fourier coefficients { ˆ fk}N

k=−N are given for a piecewise

analytic function f(x) in [−1, 1] Use edge detection to determine sub-intervals of smoothness [a, b] Compute the Fourier partial sum in [a, b]: SNf(x) = PN

k=−N ˆ

fkeikπx Reproject SNf(x) onto a new basis for x ∈ [a, b]: PM ` SNf(x) ´ → f(x) Many other algorithms available (Banerjee & Geer, Driscoll & Fornberg, Eckhoff, Jung & Shizgal, Solomonoff, Tadmor & Tanner, more ...)

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

Spectral Reprojection (Gottlieb & Shu)

Goal: approximate a function that is smooth but not periodic for ξ ∈ [−1, 1] (after linear transformation from x ∈ [a, b]) We could use an orthogonal polynomial expansion, PMf(x(ξ)) =

M

X

l=0

alψl(ξ), al = 1 γl < f, ψl >ω but we have the Fourier expansion SNf(x(ξ)) =

N

X

k=−N

ˆ fkeikπx(ξ) Can we use SNf to approximate al? What should PM look like?

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

Spectral Reprojection (Gottlieb & Shu)

Determine appropriate orthogonal polynomial basis {ψl(ξ)}M

l=0 on [−1, 1]

Build Fourier approximation inside smooth region [a, b]: SNf(x(ξ)) =

N

X

k=−N

ˆ fkeikπx(ξ), x(ξ) = ǫξ + δ Expansion coefficients for ψl(ξ) are approximated by bl =

1 γl < Snf, ψl >ω≈ al

Reprojection in smooth region [a, b]: Pm(SNf(x(ξ))) =

m

X

l=0

blψl(ξ)

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

Spectral Reprojection (Gottlieb & Shu)

Definition: A Gibbs Complementary Reprojection Basis, {ψl}N

l=0, has the

properties:

1 For an analytic function on [−1, 1], the function’s expansion in the

• rthogonal reprojection basis is exponentially convergent, Pmf(x) → f(x)

exponentially for smooth f(x).

2 The projection of the high modes in the original basis on the low modes in

the new basis is exponentially small. If these conditions are met, then Pm(SNf(x(ξ))) → f(x) exponentially for x ∈ [a, b]

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

Error Analysis For Spectral Reprojection

Err(m, N, f, ω) := ||f − Pm(SNf)|| ≤ ||f − Pmf|| + ||Pm(f − SNf)|| = Trunc(m, f, ω) + Proj(m, N, f, ω) Trunc(m, f, ω) measures the convergence properties of the reprojection basis for m expansion coefficients (reflects first Gibbs complementary basis property)

Converges exponentially for ω(x) ≥ 0 m = βN expansion terms, 0 < β < 1, resolves the function

Proj(m, N, f, ω) measures the near orthogonality of the reprojection space Pm and the space containing the information about the function that is not known, I − SN (reflects second Gibbs complementary basis property)

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

Look at Projection Error

The projection error corresponds to the decay rate of the coefficients bl: Proj(m, N, f, ω) = Pm(f − SNf) =

m

X

l=0

ψl(ξ) 1 γl f − SNf, ψlω =

m

X

l=0

ψl(ξ) 1 γl Z 1

−1

ω(y)ψl(y)(f(x(y)) − SNf(x(y)))dy =

m

X

l=0

1 γl X

|k|>N

ˆ fkψl(ξ) Z 1

−1

eiπkx(y)ω(y)ψl(y)dy

• Error will be small if the weight function ω is chosen so that the

corresponding integral is very small.

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

Reprojection using the Gegenbauer polynomial basis (Gottlieb, Shu, et. al)

Define the weight function ω as: ωλ(ξ) = (1 − ξ2)λ− 1

2 . Large λ ensures

that R 1

−1 eiπkx(y)ω(y)ψl(y)dy is small. Subsequently, the errors from the

boundaries do not enter the approximation Corresponding reprojection basis, ψl(ξ) for ωλ(ξ), are the Gegenbauer polynomials, Cλ

l (ξ) with < Cλ l , Cλ m >ωλ= 1 γl δlm

Note that Chebyshev (λ = 0) and Legendre (λ = 1

2) do not make good

reprojection bases.

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

Gegenbauer Reconstruction

If λ = λ(N) then the reprojection coefficients bl := 1 γl Z 1

−1

SNf(x(ξ))Cλ

l (ξ)(1 − ξ2)λ(N)− 1

2 dξ,

decay exponentially. Can rewrite implementation to avoid quadrature (use FFT) For many imaging applications, small m and λ work well Robust with respect to noise in imaging data (Archibald & Gelb)

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

Shepp Logan Phantom (Archibald & Gelb)

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Magnetic Resonance Imaging Sampling Patterns Challenges in Cartesian Reconstruction

Improvement of image quality for segmentation (Archibald, Chen, Gelb & Renaut)

Gray matter segmented prob- ability maps 256 × 256 randomly generated MNI digital brain 9% noise level non-uniform tissue intensity

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Problem Formulation The Non-harmonic Kernel Examples

Outline

1 Introduction

Magnetic Resonance Imaging Sampling Patterns in MR Imaging Challenges in Cartesian Reconstruction

Spectral Reprojection

2 The Non-Uniform Data Problem

Problem Formulation The Non-harmonic Kernel Reconstruction Results using the Non-harmonic Kernel

3 Current Methods

Reconstruction Methods Error Characteristics

4 Alternate Approaches

Spectral reprojection Incorporating Edge Information Spectral reprojection for Fourier Frames

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Problem Formulation The Non-harmonic Kernel Examples

Problem Formulation

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Template function

(a) Template Function

−30 −20 −10 10 20 30 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 ω | ˆ f(ω)| Fourier Transform

(b) Fourier Transform

Given these coefficients, can we/how do we reconstruct the function? What accuracy can we achieve given a finite (usually small) number of coefficients? Computational issue – no FFT available

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Problem Formulation The Non-harmonic Kernel Examples

Problem Formulation

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Template function

(c) Template Function

−30 −20 −10 10 20 30 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 ωk | ˆ f(ωk)| Fourier Transform Fourier transform Acquired samples

(d) Fourier Coefficients, N = 32

Given these coefficients, can we/how do we reconstruct the function? What accuracy can we achieve given a finite (usually small) number of coefficients? Computational issue – no FFT available

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Problem Formulation The Non-harmonic Kernel Examples

Problem Formulation

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Template function

(e) Template Function

−30 −20 −10 10 20 30 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 ωk | ˆ f(ωk)| Fourier Transform Fourier transform Acquired samples

(f) Fourier Coefficients, N = 32

Given these coefficients, can we/how do we reconstruct the function? What accuracy can we achieve given a finite (usually small) number of coefficients? Computational issue – no FFT available

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Problem Formulation The Non-harmonic Kernel Examples

Problem Formulation

Let f be defined on R and supported in (−π, π) with Fourier transform ˆ f(ω) = 1 2π Z π

−π

f(x)e−iωxdx, ω ∈ R Objective Recover f given a finite number of its non-harmonic Fourier coefficients, ˆ f(ωk), k = −N, ..., N ωk not necessarily ∈ Z We are particularly interested in sampling patterns with variable sampling density The underlying function f is piecewise-smooth

−40 −30 −20 −10 10 20 30 40 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 ω Sampling scheme −0.5 −0.4 −0.3 −0.2 −0.1 0.1 0.2 0.3 0.4 0.5 −0.5 −0.4 −0.3 −0.2 −0.1 0.1 0.2 0.3 0.4 0.5 kx ky Spiral scan trajectory

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Problem Formulation The Non-harmonic Kernel Examples

Problem Formulation

Let f be defined on R and supported in (−π, π) with Fourier transform ˆ f(ω) = 1 2π Z π

−π

f(x)e−iωxdx, ω ∈ R Objective Recover f given a finite number of its non-harmonic Fourier coefficients, ˆ f(ωk), k = −N, ..., N ωk not necessarily ∈ Z We are particularly interested in sampling patterns with variable sampling density The underlying function f is piecewise-smooth

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) −3 −2 −1 1 2 3 −1.5 −1 −0.5 0.5 1 1.5 x f(x)

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Problem Formulation The Non-harmonic Kernel Examples

Sampling Patterns

Jittered Sampling: ωk = k ± τk, τk ∼ U[0, θ], k = −N, −(N − 1), ..., N Log Sampling: |ωk| is logarithmically distributed between 10−v and N, with v > 0 and 2N + 1 being the total number of samples.

2 4 6 8 10 12 14 16 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 ω

(k) Jittered Sampling

2 4 6 8 10 12 14 16 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 ω

(l) Log Sampling

Figure: Non-uniform Sampling Schemes (right half plane), N = 16

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Problem Formulation The Non-harmonic Kernel Examples

The Non-harmonic Reconstruction Kernel

Standard (harmonic) Fourier reconstruction: SNf(x) = X

|k|≤N

ˆ f(k)eikx = (f ∗ DN)(x) where DN(x) = X

|k|≤N

eikx is the Dirichlet kernel. The non-harmonic Fourier reconstruction: SN ˜ f(x) = X

|k|≤N

ˆ f(ωk)eiωkx = (f ∗ AN)(x) where AN(x) = X

|k|≤N

eiωkx is the non-harmonic kernel. The non-harmonic kernels do not constitute an

• rthogonal basis for span {eikx, |k| ≤ N}

0.5 1 1.5 2 2.5 3 −40 −20 20 40 60 80 100 120 140 x DN(x) Dirichlet

Figure: The Dirichlet Kernel plotted on the right half plane, N = 64

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Problem Formulation The Non-harmonic Kernel Examples

The Non-harmonic Reconstruction Kernel

Standard (harmonic) Fourier reconstruction: SNf(x) = X

|k|≤N

ˆ f(k)eikx = (f ∗ DN)(x) where DN(x) = X

|k|≤N

eikx is the Dirichlet kernel. The non-harmonic Fourier reconstruction: SN ˜ f(x) = X

|k|≤N

ˆ f(ωk)eiωkx = (f ∗ AN)(x) where AN(x) = X

|k|≤N

eiωkx is the non-harmonic kernel. The non-harmonic kernels do not constitute an

• rthogonal basis for span {eikx, |k| ≤ N}

0.5 1 1.5 2 2.5 3 −40 −20 20 40 60 80 100 120 140 x AN(x) Non−harmonic Dirichlet

(a) Jittered Sampling

0.5 1 1.5 2 2.5 3 −40 −20 20 40 60 80 100 120 140 x AN(x) Non−harmonic Dirichlet

(b) Log Sampling

Figure: Non-harmonic Kernel, N = 64

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Problem Formulation The Non-harmonic Kernel Examples

The Non-harmonic Reconstruction Kernel

Standard (harmonic) Fourier reconstruction: SNf(x) = X

|k|≤N

ˆ f(k)eikx = (f ∗ DN)(x) where DN(x) = X

|k|≤N

eikx is the Dirichlet kernel. The non-harmonic Fourier reconstruction: SN ˜ f(x) = X

|k|≤N

ˆ f(ωk)eiωkx = (f ∗ AN)(x) where AN(x) = X

|k|≤N

eiωkx is the non-harmonic kernel. The non-harmonic kernels do not constitute an

• rthogonal basis for span {eikx, |k| ≤ N}

5 10 15 20 25 30 35 −0.4 −0.2 0.2 0.4 0.6 0.8 1 x

(a) Jittered Sampling

5 10 15 20 25 30 35 0.75 0.8 0.85 0.9 0.95 1 x

(b) Log Sampling

Figure: Autocorrelation plot of the kernels

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Problem Formulation The Non-harmonic Kernel Examples

Non-harmonic Kernels

0.5 1 1.5 2 2.5 3 −20 −10 10 20 30 40 50 60 70 x AN(x) Dirichlet

(a) Dirichlet, N = 32

0.5 1 1.5 2 2.5 3 −20 −10 10 20 30 40 50 60 70 x AN(x) Non−harmonic

(b) Jittered, N = 32

0.5 1 1.5 2 2.5 3 10 20 30 40 50 60 70 x AN(x) Non−harmonic

(c) Log, N = 32

0.5 1 1.5 2 2.5 3 −100 −50 50 100 150 200 250 300 x AN(x) Dirichlet

(d) Dirichlet, N = 128

0.5 1 1.5 2 2.5 3 −100 −50 50 100 150 200 250 300 x AN(x) Non−harmonic

(e) Jittered, N = 128

0.5 1 1.5 2 2.5 3 50 100 150 200 250 300 x AN(x) Non−harmonic

(f) Log, N = 128

0.5 1 1.5 2 2.5 3 −200 −100 100 200 300 400 500 600 x AN(x) Dirichlet

(g) Dirichlet, N = 256

0.5 1 1.5 2 2.5 3 −200 −100 100 200 300 400 500 600 x AN(x) Non−harmonic

(h) Jittered, N = 256

0.5 1 1.5 2 2.5 3 50 100 150 200 250 300 350 400 450 500 550 x AN(x) Non−harmonic

(i) Log, N = 256

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Problem Formulation The Non-harmonic Kernel Examples

Reconstruction Examples

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) True Gridding

(a) “Jittered” Sampling

−4 −3 −2 −1 1 2 3 4 −10 10 20 30 40 50 x f(x) True Non−equispaced sum

(b) “Log” Sampling (c) “Spiral” Sampling

Figure: Non-harmonic Fourier sum Reconstruction,N = 128

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Outline

1 Introduction

Magnetic Resonance Imaging Sampling Patterns in MR Imaging Challenges in Cartesian Reconstruction

Spectral Reprojection

2 The Non-Uniform Data Problem

Problem Formulation The Non-harmonic Kernel Reconstruction Results using the Non-harmonic Kernel

3 Current Methods

Reconstruction Methods Error Characteristics

4 Alternate Approaches

Spectral reprojection Incorporating Edge Information Spectral reprojection for Fourier Frames

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Conventional Reconstruction Methods

Several approaches available to perform reconstruction Convolutional Gridding – most popular Uniform Resampling Iterative Methods – “Fix” the quadrature rule while evaluating the non-harmonic sum SN ˜ f(x) =

N

X

k=−N

αk ˆ f(ωk)eiωkx – αk are density compensation factors e.g., αk = ωk+1 − ωk – Evaluated using a “non-uniform” FFT

−10 −8 −6 −4 −2 2 4 6 8 10 −1 1 2 3 4 5 6 f^ Trapezoidal Rule − Unequal subintervals ω ← → ← → ∆p ∆q ωq ωq+1 ωp ωp+1 f^(ωp+1) f^(ωp) f^(ωq) f^(ωq+1)

Figure: Evaluating the non-uniform Fourier sum

Although there are distinct difference in methodology and computational cost, reconstruction accuracy is similar in most schemes. We will look at convolutional gridding and uniform resampling to obtain an intuitive understanding of the problems in reconstruction.

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Conventional Reconstruction Methods

Several approaches available to perform reconstruction Convolutional Gridding – most popular Uniform Resampling Iterative Methods

−30 −20 −10 10 20 30 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 ωk | ˆ f(ωk)| Fourier Transform Fourier transform Acquired samples

(a) Non-harmonic modes

−150 −100 −50 50 100 150 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 k F(k) Fourier Coefficients True Interpolated

(b) Obtain uniform modes

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Reconstruction True URS

(c) (Filtered) Fourier sum

Figure: Uniform Resampling

Although there are distinct difference in methodology and computational cost, reconstruction accuracy is similar in most schemes. We will look at convolutional gridding and uniform resampling to obtain an intuitive understanding of the problems in reconstruction.

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Conventional Reconstruction Methods

Several approaches available to perform reconstruction Convolutional Gridding – most popular Uniform Resampling Iterative Methods

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Reconstruction (2 iterations) True CG

(a) after iteration 2

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Reconstruction (4 iterations) True CG

(b) after iteration 4

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Reconstruction (15 iterations) True CG

(c) after iteration 15

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Reconstruction (45 iterations) True CG

(d) after iteration 45

Figure: Iterative Reconstruction

Although there are distinct difference in methodology and computational cost, reconstruction accuracy is similar in most schemes. We will look at convolutional gridding and uniform resampling to obtain an intuitive understanding of the problems in reconstruction.

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Conventional Reconstruction Methods

Several approaches available to perform reconstruction Convolutional Gridding – most popular Uniform Resampling Iterative Methods

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Reconstruction (2 iterations) True CG

(a) after iteration 2

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Reconstruction (4 iterations) True CG

(b) after iteration 4

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Reconstruction (15 iterations) True CG

(c) after iteration 15

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Reconstruction (45 iterations) True CG

(d) after iteration 45

Figure: Iterative Reconstruction

Although there are distinct difference in methodology and computational cost, reconstruction accuracy is similar in most schemes. We will look at convolutional gridding and uniform resampling to obtain an intuitive understanding of the problems in reconstruction.

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Convolutional Gridding

Instead of approximating f by

N

X

k=−N

αk ˆ f(ωk)eiωkx, we compute SN ˜ g(x) =

N

X

k=−N

ˆ ˜ g(k)eikx

1 Map the non-uniform modes to a uniform grid via convolution. The new

coefficients on the uniform grid are therefore given by ˆ ˜ g(k) = ˆ f ∗ ˆ φ ˛ ˛ ˛

ω=k ≈

X

m st. |k−ωm|≤q

αm ˆ f(ωm)ˆ φ(k − ωm)

2 Compute a (filtered) Fourier partial sum. 3 “Compensate” for the mapping operation (divide by φ(x)). 4 Choose the interpolating function φ to be essentially bandlimited, i.e.,

ˆ φ(ω) ≈ 0 |ω| > q, q ∈ R, small φ(x) ≈ 0 |x| > π φ(x) = 0 x ∈ [−π, π]

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Convolutional Gridding

Instead of approximating f by

N

X

k=−N

αk ˆ f(ωk)eiωkx, we compute SN ˜ g(x) =

N

X

k=−N

ˆ ˜ g(k)eikx

1 Map the non-uniform modes to a uniform grid via convolution. The new

coefficients on the uniform grid are therefore given by ˆ ˜ g(k) = ˆ f ∗ ˆ φ ˛ ˛ ˛

ω=k ≈

X

m st. |k−ωm|≤q

αm ˆ f(ωm)ˆ φ(k − ωm)

2 Compute a (filtered) Fourier partial sum. 3 “Compensate” for the mapping operation (divide by φ(x)). 4 Choose the interpolating function φ to be essentially bandlimited, i.e.,

ˆ φ(ω) ≈ 0 |ω| > q, q ∈ R, small φ(x) ≈ 0 |x| > π φ(x) = 0 x ∈ [−π, π]

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Convolutional Gridding

Instead of approximating f by

N

X

k=−N

αk ˆ f(ωk)eiωkx, we compute SN ˜ g(x) =

N

X

k=−N

ˆ ˜ g(k)eikx

1 Map the non-uniform modes to a uniform grid via convolution. The new

coefficients on the uniform grid are therefore given by ˆ ˜ g(k) = ˆ f ∗ ˆ φ ˛ ˛ ˛

ω=k ≈

X

m st. |k−ωm|≤q

αm ˆ f(ωm)ˆ φ(k − ωm)

2 Compute a (filtered) Fourier partial sum. 3 “Compensate” for the mapping operation (divide by φ(x)). 4 Choose the interpolating function φ to be essentially bandlimited, i.e.,

ˆ φ(ω) ≈ 0 |ω| > q, q ∈ R, small φ(x) ≈ 0 |x| > π φ(x) = 0 x ∈ [−π, π]

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Convolutional Gridding

Instead of approximating f by

N

X

k=−N

αk ˆ f(ωk)eiωkx, we compute SN ˜ g(x) =

N

X

k=−N

ˆ ˜ g(k)eikx

1 Map the non-uniform modes to a uniform grid via convolution. The new

coefficients on the uniform grid are therefore given by ˆ ˜ g(k) = ˆ f ∗ ˆ φ ˛ ˛ ˛

ω=k ≈

X

m st. |k−ωm|≤q

αm ˆ f(ωm)ˆ φ(k − ωm)

2 Compute a (filtered) Fourier partial sum. 3 “Compensate” for the mapping operation (divide by φ(x)). 4 Choose the interpolating function φ to be essentially bandlimited, i.e.,

ˆ φ(ω) ≈ 0 |ω| > q, q ∈ R, small φ(x) ≈ 0 |x| > π φ(x) = 0 x ∈ [−π, π]

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Convolutional Gridding

−30 −25 −20 −15 −10 −5 5 10 0.05 0.1 ω (f ∗φ)(ω)|ω=−

10

Gridding f^(ω) φ^(k−ω)

Figure: Gridding: ˆ ˜ g = ˆ f ∗ ˆ φ

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Convolutional Gridding

−3 −2 −1 1 2 3 −4 −3 −2 −1 1 2 3 4 x g(x) FFT Reconstruction

Figure: Fourier Reconstruction of g(x) = f(x)φ(x)

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Convolutional Gridding

−3 −2 −1 1 2 3 −4 −3 −2 −1 1 2 3 4 x f(x) Compensation g(x) φ(x) f(x)

Figure: Compensation f(x) = g(x)/φ(x)

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Reconstruction Examples – Convolutional Gridding

−3 −2 −1 1 2 3 −1 −0.5 0.5 1 1.5 x f(x) True Gridding reconstruction

(a) Test Function reconstruction

−3 −2 −1 1 2 3 −0.2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x f(x) True Gridding reconstruction

(b) Cross-section of a brain scan

Figure: Gridding reconstruction, N = 128 (processed by a fourth order exponential filter)

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Uniform Resampling

Reconstruction is accomplished in two steps:

1 recover equispaced coefficients ˆ

f(k)

2 partial Fourier reconstruction using the FFT algorithm

Since f is compactly supported, we use the sampling theorem to relate ˆ f(ωk) and ˆ f(k). ˆ f(ω) =

∞

X

p=−∞

ˆ f(p) sinc(ω − p), ω ∈ R, p ∈ N To recover ˆ f(k), we have to invert the above system, i.e., solve Ax = b, Aij = sinc(ωi − j), b = n ˆ f(ωk)

• N

k=−N ,

x = n ˆ f(p)

• M

p=−M

Any number of methods to do so - iterative methods, pseudoinverse-based methods with regularization ... The condition number of the sinc matrix depends on the sampling pattern’s deviation from equispaced nodes.

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Representative Results – Uniform Resampling

−150 −100 −50 50 100 150 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 k F(k) Fourier Coefficients True Interpolated

(a) Recovered Fourier coefficients

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Reconstruction True URS

(b) Filtered reconstruction

Figure: URS solution, N = 128

Solved a square 128 × 128 system Inverted the system by computing the pseudoinverse Pseudoinverse was computed using TSVD, with a SVD threshold of 10−5

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Convolutional Gridding Error

Theorem (Convolutional Gridding Error (Viswanathan, 2009)) Let ˆ g = ˆ f ∗ ˆ φ denote the true gridding coefficients and ˆ ˜ g denote the approximate gridding coefficients. Let ∆k be the maximum distance between sampling points and dk :=

1 ∆k be the minimum sample density in the q-vicinity

• f k. Then, the gridding error at mode k is bounded by

e(k) ≤ C ·

1 d2

k , k = −N, ..., N for some positive constant C. Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Convolutional Gridding Error

Physical space reconstruction error e(x) ≈ g(x) − SN ˜ g(x) = g(x) − SNg(x) + SNg(x) − SN ˜ g(x) = X

|k|>N

ˆ g(k)eikx | {z } standard Fourier truncation + X

|k|≤N

“ ˆ g(k) − ˆ ˜ g(k) ” eikx | {z } gridding SNg suffers from Gibbs; the maximum error occurs in the vicinity of a jump (≈ 1.09 of the jump value). There is also a reduced rate of convergence with g − SNg = O (N). Gridding (sampling) error |SNg(x) − SN ˜ g(x)| = ˛ ˛ ˛ ˛ ˛ ˛ X

|k|≤N

“ ˆ g(k) − ˆ ˜ g(k) ” eikx ˛ ˛ ˛ ˛ ˛ ˛ ≤ X

|k|≤N

˛ ˛ ˛ˆ g(k) − ˆ ˜ g(k) ˛ ˛ ˛ ≤ C X

|k|≤N

1 d2

k

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Error Plots

−150 −100 −50 50 100 150 −4 −3.5 −3 −2.5 −2 −1.5 −1 k log10 |e|

(a) Error bound, log sampling

20 40 60 80 100 120 140 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 N H(ωk,N)

(b) |SNg(x) − SN ˜ g(x)| vs N

Figure: Error Plots

. Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Error Plots

30 40 50 60 70 80 90 100 110 120 −0.03 −0.02 −0.01 0.01 0.02 0.03 0.04 0.05 0.06 k F(k) Fourier Coefficients True Interpolated

(a) Fourier coefficients – High modes

20 40 60 80 100 120 140 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 N H(ωk,N)

(b) |SNg(x) − SN ˜ g(x)| vs N

Figure: Error Plots

. Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Error vs Sampling Density

The reconstruction error is e(x) ≈ X

|k|>N

ˆ g(k)eikx + X

|k|≤N

“ ˆ g(k) − ˆ ˜ g(k) ” eikx 1st term decreases as N increases 2nd term increases as N increases

−150 −100 −50 50 100 150 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 k | ˆ f(k)|

(a) Fourier coefficients

−150 −100 −50 50 100 150 −15 −10 −5 k Log | E | Error in Fourier Coefficients

(b) Coefficient error

Figure: Error in gridding coefficients

For a given sampling trajectory and function, there is a critical value Ncrit beyond which adding coefficients does not improve the accuracy. While filtering decreases the error, the underlying problem is not solved.

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Error vs Sampling Density

The reconstruction error is e(x) ≈ X

|k|>N

ˆ g(k)eikx + X

|k|≤N

“ ˆ g(k) − ˆ ˜ g(k) ” eikx 1st term decreases as N increases 2nd term increases as N increases

−150 −100 −50 50 100 150 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 k | ˆ f(k)|

(a) Fourier coefficients

−150 −100 −50 50 100 150 −15 −10 −5 k Log | E | Error in Fourier Coefficients

(b) Coefficient error

Figure: Error in uniform re-sampling

For a given sampling trajectory and function, there is a critical value Ncrit beyond which adding coefficients does not improve the accuracy. While filtering decreases the error, the underlying problem is not solved.

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Error vs Sampling Density

The reconstruction error is e(x) ≈ X

|k|>N

ˆ g(k)eikx + X

|k|≤N

“ ˆ g(k) − ˆ ˜ g(k) ” eikx 1st term decreases as N increases 2nd term increases as N increases

−150 −100 −50 50 100 150 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 k | ˆ f(k)|

(a) Fourier coefficients

−150 −100 −50 50 100 150 −12 −10 −8 −6 −4 −2 2 k Log | Enorm | Normalized Error in Fourier Coefficients

(b) Coefficient error

Figure: Error in uniform re-sampling

For a given sampling trajectory and function, there is a critical value Ncrit beyond which adding coefficients does not improve the accuracy. While filtering decreases the error, the underlying problem is not solved.

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Reconstruction Methods Error Characteristics

Error in Recovered Coefficients – Uniform Resampling

30 40 50 60 70 80 90 100 110 120 −0.03 −0.02 −0.01 0.01 0.02 0.03 0.04 0.05 0.06 k F(k) Fourier Coefficients True Interpolated

(a) Recovered Fourier coefficients

−150 −100 −50 50 100 150 −15 −10 −5 k Log | E | Error in Fourier Coefficients

(b) Error in recovered coefficients

Figure: URS solution, N = 128

The condition number of the resampling matrix is directly related to the error in the coefficients.

. Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Spectral reprojection Incorporating Edge Information Spectral reprojection fo

Outline

1 Introduction

Magnetic Resonance Imaging Sampling Patterns in MR Imaging Challenges in Cartesian Reconstruction

Spectral Reprojection

2 The Non-Uniform Data Problem

Problem Formulation The Non-harmonic Kernel Reconstruction Results using the Non-harmonic Kernel

3 Current Methods

Reconstruction Methods Error Characteristics

4 Alternate Approaches

Spectral reprojection Incorporating Edge Information Spectral reprojection for Fourier Frames

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Spectral reprojection Incorporating Edge Information Spectral reprojection fo

Gibbs Phenomenon, non-uniform case

We still have the “Gibbs phenomenon” – non-physical oscillations at discontinuities, and a reduced rate of convergence (first order). Hence, we require a large number of coefficients to get acceptable reconstructions. However, by formulation of the sampling scheme and recovery procedure, the coefficients recovered at large ω are inaccurate. = ⇒ we need more coefficients, but the coefficients we get are inaccurate! Spectral Reprojection – High frequency modes of f in the original basis have exponentially small contributions on the low modes in the new basis

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Spectral reprojection Incorporating Edge Information Spectral reprojection fo

Reducing the Impact of the High Mode Coefficients using Spectral Reprojection (Viswanathan, Cochran, Gelb, & Renaut, 2010)

Spectral reprojection expansion coefficients: 1 γλ

l

< SN ˜ g, Cλ

l >ω(λ)= 1

γλ

l

Z 1

−1

(1 − η2)λ−1/2Cλ

l (η)

X

|k|≤N

ˆ ˜ g(k)eiπkηdη

1 2 3 4 5 6 7 8 9 10 11 −5 −4.5 −4 −3.5 −3 −2.5 −2 −1.5 −1 −0.5 k log 10 |ˆ g(k)|

Figure: Decay of Gegenbauer coefficients

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Spectral reprojection Incorporating Edge Information Spectral reprojection fo

Spectral Reprojection Error of Gridded Coefficients

Err(m, N, g, ω) := ||g − Pm(SN ˜ g)|| ≤ ||g − Pmg|| + ||Pm(g − SNg)|| + ||PmSNg − PmSN ˜ g|| ||PmSNg − PmSN ˜ g||∞ ≤ Const Pm

l=0

P

|k|≤N 1 d2

k

˛ ˛ ˛

Cλ

l (1)

γλ

l

R 1

−1(1 − η2)λ−1/2Cλ l (η)eiπkηdη

˛ ˛ ˛ ˛ ˛ ˛ ˛ Cλ

l (1)

γλ

l

Z 1

−1

(1 − η2)λ−1/2Cλ

l (η)eiπkηdη

˛ ˛ ˛ ˛ ≤ Γ(λ)(l + λ)Γ(l + 2λ) l!Γ(2λ) „ 2 π|k| «λ The corresponding gridding error does not increase rapidly with N

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Spectral reprojection Incorporating Edge Information Spectral reprojection fo

Spectral Reprojection – Contribution from Gridding Error

20 40 60 80 100 120 140 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 N log10 H(dk,N)

(a) Gridding Error λ = 2

20 40 60 80 100 120 140 0.12 0.14 0.16 0.18 0.2 0.22 0.24 N log10 H(dk,N)

(b) Gridding Error λ = 3

20 40 60 80 100 120 140 0.13 0.135 0.14 0.145 0.15 0.155 0.16 0.165 0.17 0.175 N log10 H(dk,N)

(c) Gridding Error λ = 4

20 40 60 80 100 120 140 0.13 0.1305 0.131 0.1315 0.132 0.1325 0.133 0.1335 0.134 N 10 k

(d) Gridding Error λ = 8

Figure: Contribution from Gridding Error in Spectral Reprojection

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Spectral reprojection Incorporating Edge Information Spectral reprojection fo

Gegenbauer Reconstruction - Results

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) True Fourier (256 coeffs.) Gegenbauer(121 coeffs.)

(a) Reconstruction

−3 −2 −1 1 2 3 −4.5 −4 −3.5 −3 −2.5 −2 −1.5 −1 −0.5 x log10 |e| Filtered Fourier, 256 coeffs. Gegenbauer, 121 coeffs.

(b) Reconstruction error

Figure: Gegenbauer reconstruction from log spectral samples

Filtered Fourier (second-order exponential) reconstruction uses 256 coefficients Gegenbauer reconstruction uses 2N + 1 = 121 coefficients Edges detected using the concentration method (Gelb & Tadmor) with thresholding and non-linear post-processing. Parameters (m, λ) chosen proportional to size of reconstruction interval.

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Spectral reprojection Incorporating Edge Information Spectral reprojection fo

Error Plots

20 30 40 50 60 70 80 90 100 110 120 1 2 3 4 5 6 N f − PmSNf 2 2nd order exp. 4th order exp. 8th order exp. Gegenbauer

(a) 2-norm error

20 30 40 50 60 70 80 90 100 110 120 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 N f − PmSNf ∞ 2nd order exp. 4th order exp. 8th order exp. Gegenbauer

(b) Maximum-norm error

Figure: Error Plots – Filtered and Gegenbauer Reconstruction

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Spectral reprojection Incorporating Edge Information Spectral reprojection fo

Incorporating Edge Information into Reconstruction (Viswanathan, 2010)

Obtain edge information by Using the concentration method on the recovered coefficients Sσ

N[ ˜

f](x) = i

N

X

k=−N

ˆ ˜ f(k) sgn(k) σ „|k| N « eikx Solving for the jump function directly from the non-harmonic Fourier data Let ˘ f(ωk) = i sgn(ωk) ˆ f(ωk) min

p

F{Wp}|ωk− ˘ f(ωk) 2

2+λ p 1

−3 −2 −1 1 2 3 −1 −0.5 0.5 1 1.5 x f(x) Function Jump function −3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Jump function

Figure: Edge Detection using the iterative formulation

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Spectral reprojection Incorporating Edge Information Spectral reprojection fo

Methods Incorporating Edge Information

Compute the high frequency modes using the relation

Compare

ˆ f(k) = X

p∈P

[f](ζp) e−ikζp 2πik

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Reconstruction True With edge information

(a) Reconstruction - Using edge informa- tion

30 40 50 60 70 80 90 100 110 120 130 −0.02 −0.01 0.01 0.02 0.03 k F(k) Fourier Coefficients True With edge information

(b) The high modes - Using edge infor- mation

Figure: Reconstruction of a test function using edge information

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Spectral reprojection Incorporating Edge Information Spectral reprojection fo

Spectral reprojection for Fourier Frames (Gelb & Hines, 2010)

Finite frame approximation: Assume {eiωnx}n∈Z is a frame for L2(I). Given frame coefficients { ˆ f(ωn)}N

n=−N = {< f(x), eiωnx >}N n=−N

Compute TNf = PN

n=−N ˆ

f(ωn)S−1eiωnx Filtering does not improve convergence: T σ

Nf = PN n=−N σ(ωn) ˆ

f(ωn)S−1eiωnx f Spectral reprojection yields exponential convergence. Theorem (Gelb and Hines (2010)): Let {eiωnx}n∈Z be a frame generated by a balanced sampling sequence. Suppose we are given the first 2N + 1 frame coefficients of f ∈ L2[−1, 1]. If λ = αN and M = βN for 0 < α < 1 and 0 < β < 1, then there exists a constant 0 < q < 1 such that ||P λ

M(f − TNf)||∞ ≤ cN 2qN

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Spectral reprojection Incorporating Edge Information Spectral reprojection fo

Spectral reprojection for Fourier Frames (Gelb & Hines, 2010)

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0.5 1 1.5 x f(x) TN f(x) Pµ

M T N f(x)

(a) Reconstruction with Gegenbauer frame algorithm

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −8 −7 −6 −5 −4 −3 −2 −1 x Error in Gegenbauer reconstruction N = 64 N = 128 N = 256

(b) pointwise error with Gegenbauer frame algorithm

Figure: Gegenbauer frame reconstruction

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Spectral reprojection Incorporating Edge Information Spectral reprojection fo

Spectral reprojection for Fourier Frames (Gelb & Hines, 2010)

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0.5 1 1.5 x f(x) TN f(x) Pµ

M T N f(x)

(a) Reconstruction with Gegenbauer frame algorithm

5 10 15 20 25 30 35 40 −20 −18 −16 −14 −12 −10 −8 −6 −4 −2 m Gegenbauer coeff, N = 256 µ = 8 µ = 16 µ = 24 µ = 32

(b) coefficient decay

Figure: Gegenbauer frame reconstruction

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Spectral reprojection Incorporating Edge Information Spectral reprojection fo

Spectral reprojection for Fourier Frames (Gelb & Hines, 2010)

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −6 −5 −4 −3 −2 −1 1 x Error in frame reconstruction N = 64 N = 128 N = 256

(a) pointwise error with frame recon- struction

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −8 −7 −6 −5 −4 −3 −2 −1 x Error in Gegenbauer reconstruction N = 64 N = 128 N = 256

(b) pointwise error with Gegenbauer frame algorithm

Figure: Gegenbauer frame reconstruction

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Spectral reprojection Incorporating Edge Information Spectral reprojection fo

Summary

Problem: Fourier reconstruction of piecewise smooth functions from non-uniform coefficients – MRI, SAR Conventional reconstruction methods from MRI community (density compensation, uniform resampling, iterative methods) Error analysis and characteristics related to sampling density Spectral reprojection mitigates the impact of the error from the high frequency coefficients Fourier based edge information can help to obtain better approximation to high frequency coefficients Fourier frames may provide a better alternative in reconstruction since no interpolation is needed Special thanks to Adityavikram Viswanathan, ASU EE PhD 2010, currently a postdoctoral fellow at ACM. California Institute of Technology.

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data

Introduction Non-Uniform Data Current Methods Alternate Approaches Spectral reprojection Incorporating Edge Information Spectral reprojection fo

Special thanks

Special thanks to the organizers and to my sons (for letting me come)

Figure: Sam and Josh Bagatell

Anne Gelb Fourier Reconstruction from Non-Uniform Spectral Data