Edge Detection from Spectral Data With Applications to PSF - - PowerPoint PPT Presentation

Edge Detection from Spectral Data

With Applications to PSF Estimation and Fourier Reconstruction

Aditya Viswanathan

School of Electrical, Computer and Energy Engineering, Arizona State University aditya.v@asu.edu

Oct 22 2009

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 1 / 41

• Prof. Anne Gelb is with the School of

Mathematical and Statistical Sciences, Arizona State University

• Prof. Douglas Cochran is with the School of

Electrical, Computer and Energy Engineering, Arizona State University

• Prof. Rosemary Renaut is with the School of

Mathematical and Statistical Sciences, Arizona State University

• Dr. Wolfgang Stefan is with the Department
• f Computational and Applied Mathematics,

Rice University Research supported in part by National Science Foundation grants DMS 0510813 and DMS 0652833 (FRG).

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 2 / 41

Introduction Motivation

Why do we need edge data?

Solve PDE’s with shocks more accurately. Identify tissue boundaries in medical images and segment them. Reconstruct piecewise-analytic functions from Fourier and other spectral expansion coefficients with uniform and exponential convergence properties.

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 −4 −2 2 4 x (x pi) Un(x,t) PDE Solution 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 −10 −5 5 Jump Function Approximation [f](x) x True Solution tf = 0 tf = 50, Exp. Filtering, p=2 tf = 50, Exp. Filtering, p=8 tf = 0 tf = 50, Exp. Filtering, p=2 tf = 50, Exp. Filtering, p=8

(a) PDE solution with shock discontinuity (b) MRI tissue bound- ary detection

−3 −2 −1 1 2 3 −2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 2.5 Function Reconstruction x g(x) True Fourier reconstruction

(c) Fourier reconstruc- tion showing Gibbs

−3 −2 −1 1 2 3 −2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 2.5 Function Reconstruction x g(x) True

• Alt. reconstruction

(d) Alt. reconstruction

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 3 / 41

Introduction Motivation

Application – Magnetic Resonance Imaging

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

(e) Acquired Fourier Samples

−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

(f) Sampling Trajectory

Reconstructed phantom

(g) Reconstructed Image

Figure: MR Imaginga

aSampling pattern courtesy Dr. Jim Pipe, Barrow Neurological Institute,

Phoenix, Arizona

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 4 / 41

Introduction Problem Statement

Problem Statement

Objective: To recover location, magnitude and sign of jump discontinuities from a finite number of spectral coefficients Assumptions: f is 2π-periodic and piecewise-smooth function in [−π, π). Its jump function is defined as [f](x) := f(x+) − f(x−) It has Fourier series coefficients ˆ f(k) = 1 2π Z π

−π

f(x)e−ikxdx , k ∈ [−N, N] A jump discontinuity is a local feature; i.e., the jump function at any point x only depends on the values of f at x+ and x−. However, ˆ f is a global representation; i.e., ˆ f(k) are obtained using values of f over the entire domain [−π, π).

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 5 / 41

Concentration Method

Outline

1 Introduction

Motivation Problem Statement

2 Jump detection using the Concentration method

The Concentration Method Concentration Factors Statistical Analysis of the Concentration Method Iterative Formulations

3 Applications

PSF Estimation in Blurring Problems Applications to Fourier Reconstruction

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 6 / 41

Concentration Method The Concentration Method

Getting Jump Data from Fourier Coefficients

Let f contain a single jump at x = ζ. ˆ f(k) = 1 2π Z π

−π

f(x)e−ikxdx = 1 2π Z ζ−

−π

f(x)e−ikxdx + Z π

ζ+ f(x)e−ikxdx

! = 1 2π f(x) e−ikx −ik ˛ ˛ ˛ ˛

ζ− −π

− Z ζ−

−π

f ′(x)e−ikx −ik dx + f(x) e−ikx −ik ˛ ˛ ˛ ˛

π ζ+

− Z π

ζ+ f ′(x) e−ikx

−ik dx « = 1 2π f(ζ−)e−ikζ − f(−π)eikπ −ik − Z ζ−

−π

f ′(x)e−ikx −ik dx + f(π)e−ikπ − f(ζ+)e−ikζ −ik − Z π

ζ+ f ′(x)e−ikx

−ik dx «

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 7 / 41

Concentration Method The Concentration Method

Getting Jump Data from Fourier Coefficients

Let f contain a single jump at x = ζ. ˆ f(k) = 1 2π Z π

−π

f(x)e−ikxdx = 1 2π Z ζ−

−π

f(x)e−ikxdx + Z π

ζ+ f(x)e−ikxdx

! = 1 2π f(x) e−ikx −ik ˛ ˛ ˛ ˛

ζ− −π

− Z ζ−

−π

f ′(x)e−ikx −ik dx + f(x) e−ikx −ik ˛ ˛ ˛ ˛

π ζ+

− Z π

ζ+ f ′(x) e−ikx

−ik dx « = 1 2π f(ζ−)e−ikζ − f(−π)eikπ −ik − Z ζ−

−π

f ′(x)e−ikx −ik dx + f(π)e−ikπ − f(ζ+)e−ikζ −ik − Z π

ζ+ f ′(x)e−ikx

−ik dx «

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 7 / 41

Concentration Method The Concentration Method

Getting Jump Data from Fourier Coefficients

Let f contain a single jump at x = ζ. ˆ f(k) = 1 2π Z π

−π

f(x)e−ikxdx = 1 2π Z ζ−

−π

f(x)e−ikxdx + Z π

ζ+ f(x)e−ikxdx

! = 1 2π f(x) e−ikx −ik ˛ ˛ ˛ ˛

ζ− −π

− Z ζ−

−π

f ′(x)e−ikx −ik dx + f(x) e−ikx −ik ˛ ˛ ˛ ˛

π ζ+

− Z π

ζ+ f ′(x) e−ikx

−ik dx « = 1 2π f(ζ−)e−ikζ − f(−π)eikπ −ik − Z ζ−

−π

f ′(x)e−ikx −ik dx + f(π)e−ikπ − f(ζ+)e−ikζ −ik − Z π

ζ+ f ′(x)e−ikx

−ik dx «

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 7 / 41

Concentration Method The Concentration Method

Getting Jump Data from Fourier Coefficients

ˆ f(k) = 1 2π f(ζ−)e−ikζ − f(−π)eikπ −ik − Z ζ−

−π

f ′(x)e−ikx −ik dx + f(π)e−ikπ − f(ζ+)e−ikζ −ik − Z π

ζ+ f ′(x)e−ikx

−ik dx « = ` f(ζ+) − f(ζ−) ´ e−ikζ 2πik + f(−π)eikπ − f(π)e−ikπ 2πik + O „ 1 k2 « Since f is periodic, f(−π) = f(π) and the second term vanishes. ˆ f(k) = [f](ζ) e−ikζ 2πik + O „ 1 k2 «

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 8 / 41

Concentration Method The Concentration Method

Getting Jump Data from Fourier Coefficients

ˆ f(k) = 1 2π f(ζ−)e−ikζ − f(−π)eikπ −ik − Z ζ−

−π

f ′(x)e−ikx −ik dx + f(π)e−ikπ − f(ζ+)e−ikζ −ik − Z π

ζ+ f ′(x)e−ikx

−ik dx « = ` f(ζ+) − f(ζ−) ´ e−ikζ 2πik + f(−π)eikπ − f(π)e−ikπ 2πik + O „ 1 k2 « Since f is periodic, f(−π) = f(π) and the second term vanishes. ˆ f(k) = [f](ζ) e−ikζ 2πik + O „ 1 k2 «

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 8 / 41

Concentration Method The Concentration Method

Getting Jump Data from Fourier Coefficients

ˆ f(k) = 1 2π f(ζ−)e−ikζ − f(−π)eikπ −ik − Z ζ−

−π

f ′(x)e−ikx −ik dx + f(π)e−ikπ − f(ζ+)e−ikζ −ik − Z π

ζ+ f ′(x)e−ikx

−ik dx « = ` f(ζ+) − f(ζ−) ´ e−ikζ 2πik + f(−π)eikπ − f(π)e−ikπ 2πik + O „ 1 k2 « Since f is periodic, f(−π) = f(π) and the second term vanishes. ˆ f(k) = [f](ζ) e−ikζ 2πik + O „ 1 k2 «

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 8 / 41

Concentration Method The Concentration Method

Extracting Jump Information

Let us compute a ‘filtered’ partial Fourier sum of the form SN[f](x) =

N

X

k=−N

„iπk N « ˆ f(k)eikx SN[f](x) =

N

X

k=−N

„iπk N « ˆ f(k)eikx =

N

X

k=−N

„iπk N « » [f](ζ) e−ikζ 2πik + O „ 1 k2 «– eikx = [f](ζ) 1 2N

N

X

k=−N

eik(x−ζ) +

N

X

k=−N

O „ 1 k « eik(x−ζ) First term is a scaled (by the jump value) Dirac delta localized at x = ζ (the jump location) The second term is a manifestation of the global nature of Fourier data

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 9 / 41

Concentration Method The Concentration Method

Concentration Method (Gelb, Tadmor)

Generalized conjugate partial Fourier sum: Sσ

N[f](x) = i N

X

k=−N

ˆ f(k) sgn(k) σ „|k| N « eikx where, σk,N(η) = σ( |k|

N ) are known as concentration factors. Concentration factors have

to satisfy certain properties in order to be admissible. These include,

1 N

X

k=1

σ „ k N « sin(kx) be odd

2

σ(η) η ∈ C2(0, 1)

3

Z 1

ǫ

σ(η) η → −π, ǫ = ǫ(N) > 0 being small Under these conditions, we have the following relation (concentration property) Sσ

N[f](x) = [f](x) + O(ǫ),

ǫ = ǫ(N) > 0 being small

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 10 / 41

Concentration Method The Concentration Method

Illustration of the Concentration Method

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 11 / 41

Concentration Method Concentration Factors

Classical Concentration Factors

Factor Expression Trigonometric σT (η) = π sin(α η) Si(α) Si(α) = Z α sin(x) x dx Polynomial σP (η) = −p π ηp p is the order of the factor Exponential σE(η) = C η exp „ 1 α η (η − 1) « C - normalizing constant; α - order C = π R 1− 1

N 1 N

exp “

1 α τ (τ−1)

” dτ

Table: Examples of concentration factors

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 12 / 41

Concentration Method Concentration Factors

Concentration Factors - Visualization in Fourier Space

−80 −60 −40 −20 20 40 60 80 0.5 1 1.5 2 2.5 3 3.5 4 k sig(eta) Concentration Factors Trigonometric Polynomial Exponential

Figure: Envelopes of the Concentration Factors in Fourier Space

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 13 / 41

Concentration Method Concentration Factors

Concentration Factors - Typical Responses

−3 −2 −1 1 2 3 −3.5 −3 −2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 SN[f](x) x f SN[f]

(a) Trigonometric Factor

−3 −2 −1 1 2 3 −3.5 −3 −2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 SN[f](x) x f SN[f]

(b) Polynomial Factor

−3 −2 −1 1 2 3 −3.5 −3 −2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 SN[f](x) x f SN[f]

(c) Exponential Factor

Figure: Characteristic Responses of the Concentration Factors to a Sawtooth Function

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 14 / 41

Concentration Method Concentration Factors

Examples

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) f SN[f]

(a) Trigonometric Factor

−3 −2 −1 1 2 3 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 x f(x) f SN[f]

(b) Exponential Factor

Figure: Jump Response, N = 128

For images, we apply the concentration method to each dimension separately Sσ

N[f](x(¯

y)) = i

N

X

l=−N

sgn(l) σ „ |l| N « ·

N

X

k=−N

ˆ fk,l ei(kx+l¯

y)

where the overbar represents the dimension held constant.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 15 / 41

Concentration Method Concentration Factors

Examples

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) f SN[f]

(a) Trigonometric Factor

−3 −2 −1 1 2 3 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 x f(x) f SN[f]

(b) Trigonometric Factor

Figure: Jump Response, N = 128

For images, we apply the concentration method to each dimension separately Sσ

N[f](x(¯

y)) = i

N

X

l=−N

sgn(l) σ „ |l| N « ·

N

X

k=−N

ˆ fk,l ei(kx+l¯

y)

where the overbar represents the dimension held constant.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 15 / 41

Concentration Method Concentration Factors

Examples

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) f SN[f]

(a) Trigonometric Factor

−3 −2 −1 1 2 3 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 x f(x) f SN[f]

(b) Trigonometric Factor

Figure: Jump Response, N = 128

For images, we apply the concentration method to each dimension separately Sσ

N[f](x(¯

y)) = i

N

X

l=−N

sgn(l) σ „ |l| N « ·

N

X

k=−N

ˆ fk,l ei(kx+l¯

y)

where the overbar represents the dimension held constant.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 15 / 41

Concentration Method Concentration Factors

Concentration Factor Design (Viswanathan, Gelb)

Choose a template function f, eg. the sawtooth function f(x) = 

x−π 2

x < 0

x+π 2

x > 0 [f](x) =  −π x = 0 else Design σ such that Sσ

N[f](x) = i N

X

k=−N

ˆ f(k) sgn(k) σ „|k| N « eikx closely matches [f]. Typical problem formulation min

σ

φ0(σ) subject to φm(σ) = cm, ψn(σ) ≤ cn where cm, cn are constants and the objective and constraints are typically norm measures or the conjugate partial sum evaluated at certain points/intervals in the domain.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 16 / 41

Concentration Method Concentration Factors

Equivalence of Constraints and Admissibility Conditions

1 N

X

k=1

σ „ k N « sin(kx) be odd ⇐ ⇒ conjugate kernel real, odd; hence σ(−η) = σ(η) with σ(0) = 0.

2

σ(η) η ∈ C2(0, 1) ⇐ ⇒ include σ2 or similar smoothness metric in objective.

3

Z 1

ǫ

σ(η) η → −π, ǫ = ǫ(N) > 0 ⇐ ⇒ Sσ

N[f](ζ) = [f](ζ)

ζ being a point of discontinuity. (2) is not really an equivalence, but this turns out be practically inconsequential.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 17 / 41

Concentration Method Concentration Factors

Examples

Problem Formulation 1 min

σ

[f] − Sσ

N[f] 2

subject to Sσ

N[f](0) = −π

Problem Formulation 2 min

σ

[f] − Sσ

N[f] 1

subject to Sσ

N[f](0) = −π

Problem Formulation 3 min

σ

[f] − Sσ

N[f] 1

subject to Sσ

N[f](0) = −π

|Sσ

N[f](x)| ≤ 10−3, |x| ∈ (0.25, 3)

−4 −3 −2 −1 1 2 3 4 −4 −3 −2 −1 1 2 x [f](x) f SN[f]

(a) Response

−80 −60 −40 −20 20 40 60 80 0.5 1 1.5 2 2.5 3 3.5 k σ(k)

(b) Conc. Factor

Figure: Problem Formulation 1, N = 64

−4 −3 −2 −1 1 2 3 4 −4 −3 −2 −1 1 2 x [f](x) f SN[f]

(a) Response

−80 −60 −40 −20 20 40 60 80 0.5 1 1.5 2 2.5 3 k σ(k)

(b) Conc. Factor

Figure: Problem Formulation 2, N = 64

−4 −3 −2 −1 1 2 3 4 −4 −3 −2 −1 1 2 x [f](x) f SN[f]

(a) Response

−80 −60 −40 −20 20 40 60 80 0.5 1 1.5 2 2.5 3 3.5 k σ(k)

(b) Conc. Factor

Figure: Problem Formulation 3, N = 64

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 18 / 41

Concentration Method Concentration Factors

The Missing Coefficient Problem

Let us assume that mid-range Fourier coefficients are missing (eg., modes |30 − 40| for N = 64). Since the response and kernel are real, we assume that both ˆ f(±p) are missing. Use of the standard concentration factors results in spurious oscillations in smooth regions.

−4 −3 −2 −1 1 2 3 4 −3 −2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x [f](x) f SN[f]

(a) Response – Trigonometric factor

−4 −3 −2 −1 1 2 3 4 −2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x [f](x) f SN[f]

(b) Response – Exponential factor

Figure: Jump approximation with Fourier modes |30 − 40| missing, N = 64

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 19 / 41

Concentration Method Concentration Factors

Design for Missing Spectral Coefficients

Choose a template function f, eg. the sawtooth function f(x) = 

x−π 2

x < 0

x+π 2

x > 0 [f](x) =  −π x = 0 else Explicitly set ˆ f(k) ˛ ˛ ˛

30≤|k|≤40 = 0.

Design σ such that Sσ

N[f](x) = i N

X

k=−N

ˆ f(k) sgn(k) σ „|k| N « eikx closely matches [f]. Use the standard problem formulation min

σ

φ0(σ) subject to φm(σ) = cm, ψn(σ) ≤ cn where cm, cn are constants and the objective and constraints are typically norm measures or the conjugate partial sum evaluated at certain points/intervals in the domain.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 20 / 41

Concentration Method Concentration Factors

An Example

min

σ

[f] − Sσ

N[f] 1

subject to Sσ

N[f](0) = −π

σ(η) ≥ 0 |Sσ

N[f](x)| ≤ 10−3, |x| ∈ (0.25, 3)

−4 −3 −2 −1 1 2 3 4 −4 −3 −2 −1 1 2 x [f](x) f SN[f]

(a) Response

−80 −60 −40 −20 20 40 60 80 0.5 1 1.5 2 2.5 3 k σ(k)

(b) Conc. Factor

Figure: Jump detection with missing spectral data, N = 64

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 21 / 41

Concentration Method Statistics

Edge detection in noisy data (Viswanathan, Cochran, Gelb, Cates)

Typically, building an edge detector out of the jump approximation requires comparing the jump function SN[f] against a threshold However, the jump approximations have unique characteristics - mainlobe width, sidelobe oscillations etc. Under these circumstances, when deciding whether a point is an edge, it is advantageous to take into consideration measurements in the vicinity of the point In the presence of noise, we have to additionally weight the measurements by a covariance matrix The final form of the detector (assuming additive white Gaussian noise in Fourier space) is a weighted inner product of the form M T C−1

V Y > γ

where Y is the vector of noisy jump function measurements, M is a template or jump response (noiseless) to a single step edge, CV is the covariance matrix and γ is a threshold.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 22 / 41

Concentration Method Statistics

The Details

Assume zero-mean, additive complex white Gaussian noise ˆ g(k) = ˆ f(k) + ˆ v(k) ˆ v(k) ∼ N[0, ρ2] Concentration method is linear, i.e., Sσ

N[g](x) = Sσ N[f](x) + Sσ N[v](x)

Mean: E [ Sσ

N[g](x) ] = Sσ N[f](x)

Covariance: (Cv)xa,xb

p,q

= ρ2 X

l

σp( |l| N )σq( |l| N )eil(xa−xb) The detection problem is H0 : Y = V ∼ N[0, CV] H1 : Y = M + V ∼ N[M, CV] Solve using Neyman-Pearson Lemma → H1 : Pr(Y|H1) Pr(Y|H0) > γ Detector is a generalized matched filter → H1 : M T C−1

V Y > γ

M T C−1

V M is the “SNR” metric and governs detector performance

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 23 / 41

Concentration Method Statistics

A Representative Result

−4 −3 −2 −1 1 2 3 4 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Detected jumps in the test function Noisy Function True Jumps Detected Jumps

Figure: Results - Edge Detection with Noisy Data, N = 128, ρ2 = 7.5, 3−point Trigonometric detector

performs well with the exception of false alarms in the vicinity of an edge false alarms can be addressed using statistical sidelobe mitigation methods Use of multiple concentration factors is possible and indeed encouraged; we use those combinations of concentration factors which yield high SNR

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 24 / 41

Concentration Method Iterative Formulations

In Search of Precise Jump Locations (Stefan, Viswanathan, Gelb, Renaut)

Jump approximations using the concentration method are computed using a Fourier partial sum. Since the jump function is piecewise-analytic, this approximation suffers from convergence issues. In particular, the jump has a minimum, non-zero resolution and there are spurious responses in smooth regions. Non Linear Post-processing

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

(a) The minmod response

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

(b) Enhancement of scales

Figure: Non-linear post-processing of the jump approximation

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 25 / 41

Concentration Method Iterative Formulations

Problem Formulation

In setting up our iterative formulations, We will take inspiration from sparsity enforcing regularization routines and their iterative solutions (Tadmor and Zou). We will exploit the characteristic responses of the concentration factors in our problem formulation. We may write Sσ

N[f](x) = (f ∗ Kσ N)(x) ≈ ([f] ∗ W σ N)(x)

Sσ

N[f] is the jump approximation computed using the concentration method and

concentration factor σ(η). W σ

N is the characteristic response to a unit jump using the concentration factor σ(η).

Iterative Problem Formulation min

p

Wp − SN[f] 2

2 + λ p 1

where W is a Toeplitz matrix containing shifted replicates of the characteristic response.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 26 / 41

Concentration Method Iterative Formulations

Representative Examples

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

Figure: Jump detection – Iterative Formulation (N = 40, Exponential factor)

Works well with even a small number of measurements Works with missing coefficients and non-harmonic Fourier coefficients, i.e., ˆ f(ωk) := 1 2π Z π

−π

f(x)e−iωkxdx Can be extended to blurred data (g = f ∗ h + n) by modifying the problem formulation min

p

H · W · p − SN[g] 2

2 + λ p 1

where H is a Toeplitz matrix containing shifted replicates of the blur or psf.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 27 / 41

Concentration Method Iterative Formulations

Representative Examples

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Blurred function Edges Blur −3 −2 −1 1 2 3 −1 −0.5 0.5 1 1.5 x f(x) Function Blurred function Edges Blur

Figure: Jump detection – Iterative Formulation (N = 128, Exponential factor, Gaussian Blur)

Works well with even a small number of measurements Works with missing coefficients and non-harmonic Fourier coefficients, i.e., ˆ f(ωk) := 1 2π Z π

−π

f(x)e−iωkxdx Can be extended to blurred data (g = f ∗ h + n) by modifying the problem formulation min

p

H · W · p − SN[g] 2

2 + λ p 1

where H is a Toeplitz matrix containing shifted replicates of the blur or psf.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 27 / 41

Applications

Outline

1 Introduction

Motivation Problem Statement

2 Jump detection using the Concentration method

The Concentration Method Concentration Factors Statistical Analysis of the Concentration Method Iterative Formulations

3 Applications

PSF Estimation in Blurring Problems Applications to Fourier Reconstruction

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 28 / 41

Applications PSF Estimation

Convolutional Blurring (Viswanathan, Stefan, Cochran, Gelb, Renaut)

Often, the input data we observe has been subjected to a distorting physical process during measurement, transmission or instrumentation. A large class of these distortions can be explained using a convolutional blurring model. Corrective actions often require an accurate estimate of the distortion or blur. The convolutional blurring model can be written as g = f ∗ h + n f is the true function h is the blur or the point-spread function (psf) n is noise g is the observed function Let f ∈ L2(−π, π) be piecewise-smooth. We estimate the psf starting with 2N + 1 blurred Fourier coefficients ˆ g(k), k = −N, ..., N.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 29 / 41

Applications PSF Estimation

PSF Estimation using Edge Detection

Given the blurring model g = f ∗ h + n Principle: Apply a linear edge detector. We shall assume that the edge detector can be written as a convolution with an appropriate kernel Sσ

N[g] = T (f ∗ h + n)

= (f ∗ h + n) ∗ K = f ∗ h ∗ K + n ∗ K = (f ∗ K) ∗ h + n ∗ K ≈ [f] ∗ h + ˜ n Hence, we observe shifted and scaled replicates of the psf.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 30 / 41

Applications PSF Estimation

PSF Estimation using Edge Detection

Given the blurring model g = f ∗ h + n Principle: Apply a linear edge detector. We shall assume that the edge detector can be written as a convolution with an appropriate kernel Sσ

N[g] = T (f ∗ h + n)

= (f ∗ h + n) ∗ K = f ∗ h ∗ K + n ∗ K = (f ∗ K) ∗ h + n ∗ K ≈ [f] ∗ h + ˜ n Hence, we observe shifted and scaled replicates of the psf.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 30 / 41

Applications PSF Estimation

PSF Estimation using Edge Detection

Given the blurring model g = f ∗ h + n Principle: Apply a linear edge detector. We shall assume that the edge detector can be written as a convolution with an appropriate kernel Sσ

N[g] = T (f ∗ h + n)

= (f ∗ h + n) ∗ K = f ∗ h ∗ K + n ∗ K = (f ∗ K) ∗ h + n ∗ K ≈ [f] ∗ h + ˜ n Hence, we observe shifted and scaled replicates of the psf.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 30 / 41

Applications PSF Estimation

PSF Estimation using Edge Detection

Given the blurring model g = f ∗ h + n Principle: Apply a linear edge detector. We shall assume that the edge detector can be written as a convolution with an appropriate kernel Sσ

N[g] = T (f ∗ h + n)

= (f ∗ h + n) ∗ K = f ∗ h ∗ K + n ∗ K = (f ∗ K) ∗ h + n ∗ K ≈ [f] ∗ h + ˜ n Hence, we observe shifted and scaled replicates of the psf.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 30 / 41

Applications PSF Estimation

Example (No Noise)

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

(a) True function

−3 −2 −1 1 2 3 −0.2 0.2 0.4 0.6 0.8 1 1.2 x h(x)

(b) Gaussian Blur PSF

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

(c) Blurred Function

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Blur after edge detection True blur (normalized)

(d) After applying edge detection

Figure: Function subjected to Gaussian blur, N = 128

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 31 / 41

Applications PSF Estimation

Example (No Noise)

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

(a) True function

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

(b) Motion Blur PSF

−3 −2 −1 1 2 3 −1.5 −1 −0.5 0.5 1 1.5 x g(x) Blurred Function

(c) Blurred Function

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Blur after edge detection True blur (normalized)

(d) After applying edge detection

Figure: Function subjected to motion blur, N = 128

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 31 / 41

Applications PSF Estimation

Representative Examples

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Blurred, noisy function Blur after edge detection True blur (normalized)

(a) Noisy blur estimation

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Blurred, noisy function Blur after edge detection True blur (normalized)

(b) After low-pass filtering

Figure: Function subjected to Gaussian blur, N = 128

Noise distribution – ˆ n ∼ CN(0,

1.5 (2N+1)2 )

Second picture subjected to low-pass (Gaussian) filtering It is conceivable that the parameter estimation can take into account the effect of Gaussian filtering

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 32 / 41

Applications PSF Estimation

Representative Examples

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Blurred, noisy function Blur after edge detection True blur (normalized)

(a) Noisy blur estimation

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) Function Blurred, noisy function Blur after edge detection True blur (normalized)

(b) After TV denoising

Figure: Function subjected to Motion blur, N = 128

Cannot perform conventional low-pass filtering since blur is piecewise-smooth We compute the noisy blur estimate SN[g] ≈ [f] ∗ h + n ∗ KN Denoising problem formulation min

p

p − SN[g] 2

2 + λ Lp 1

where L is a derivative matrix.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 32 / 41

Applications Fourier Reconstruction

Fourier Reconstruction of Piecewise-smooth functions

SNf(x) =

N

X

k=−N

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

−π

f(x)e−ikxdx Consequences of Gibbs Non-uniform convergence – presence of non-physical oscillations in the vicinity of discontinuities Reduced order of convergence – first order convergence even in smooth regions of the reconstruction

−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

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 33 / 41

Applications Fourier Reconstruction

Filtered Fourier Reconstructions

Filtering helps to ameliorate the effects of Gibbs, but does not eliminate it. In fact, it introduces a smearing artifact in the vicinity of a discontinuity.

−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) Filtered Reconstruction

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −16 −14 −12 −10 −8 −6 −4 −2 x Log10 |e| Log Error in Reconstruction

(b) Reconstruction error

Figure: Exponentially Filtered Reconstruction, p = 2, N = 64

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 34 / 41

Applications Fourier Reconstruction

Spectral Re-projection (Gottlieb and Shu)

Spectral re-projection schemes were formulated to resolve the Gibbs phenomenon. They involve reconstructing the function using an alternate basis, Ψ (known as a Gibbs complementary basis). Reconstruction is performed using the rapidly converging series f(x) ≈

m

X

l=0

clψl(x), where cl = fN, ψlw ψl2

w

, fN is the Fourier expansion of f Reconstruction is performed in each smooth interval. Hence, we require jump discontinuity locations High frequency modes of f have exponentially small contributions on the low modes in the new basis

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 35 / 41

Applications Fourier Reconstruction

Spectral Re-projection – An Example

−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 f(x) Gegenbauer reconstruction

(a) Reconstruction

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −16 −14 −12 −10 −8 −6 −4 −2 x Log |e| Log error N=32 N=64 N=256

(b) Reconstruction error

Figure: Spectral Reprojection Reconstructions

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 36 / 41

Applications Fourier Reconstruction

Non-harmonic Fourier Reconstruction (Viswanathan, Gelb, Cochran, Renaut)

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) 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

Figure: A motivating example

Fourier samples violate the quadrature rule for discrete Fourier expansion Computational issue – no FFT available Mathematical issue – given these coefficients, can we/how do we reconstruct the function? Resolution – what accuracy can we achieve given a finite (usually small) number of coefficients?

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 37 / 41

Applications Fourier Reconstruction

Non-harmonic Fourier Reconstruction (Viswanathan, Gelb, Cochran, Renaut)

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) 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 ωk | ˆ f(ωk)| Fourier Transform Fourier transform Acquired samples

(b) Fourier Coefficients, N = 32

Figure: A template example

Fourier samples violate the quadrature rule for discrete Fourier expansion Computational issue – no FFT available Mathematical issue – given these coefficients, can we/how do we reconstruct the function? Resolution – what accuracy can we achieve given a finite (usually small) number of coefficients?

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 37 / 41

Applications Fourier Reconstruction

Non-harmonic Fourier Reconstruction (Viswanathan, Gelb, Cochran, Renaut)

−3 −2 −1 1 2 3 −2 −1.5 −1 −0.5 0.5 1 1.5 2 x f(x) 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 ωk | ˆ f(ωk)| Fourier Transform Fourier transform Acquired samples

(b) Fourier Coefficients, N = 32

Figure: A template example

Fourier samples violate the quadrature rule for discrete Fourier expansion Computational issue – no FFT available Mathematical issue – given these coefficients, can we/how do we reconstruct the function? Resolution – what accuracy can we achieve given a finite (usually small) number of coefficients?

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 37 / 41

Applications Fourier Reconstruction

Reconstruction Procedure

1 recover equispaced coefficients ˆ

f(k)

2 partial Fourier reconstruction using the FFT algorithm

Equispaced coefficients are obtained by inverting a model derived from application of the sampling theorem.

−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 coeffi- cients

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

(b) Error in recovered coeffi- cients

−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 reconstruction

Figure: Reconstruction result, N = 128

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 38 / 41

Applications Fourier Reconstruction

Reconstruction Procedure

1 recover equispaced coefficients ˆ

f(k)

2 partial Fourier reconstruction using the FFT algorithm

Equispaced coefficients are obtained by inverting a model derived from application of the sampling theorem.

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 coeffi- cients

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

(b) Error in recovered coeffi- cients

−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 reconstruction

Figure: Reconstruction result, N = 128

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 38 / 41

Applications Fourier Reconstruction

Incorporating Edge Information

Compute the high frequency modes using the relation ˆ 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

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 39 / 41

Summary

Summary

We summarized the concentration method of edge detection. We looked at routines to design concentration factors. We showed iterative routines for accurate detection of edges. We briefly surveyed the statistical properties of concentration edge detection. We looked at applications of edge detection to

Point spread function estimation in blurring problems. Fourier reconstruction of piecewise-smooth functions – spectral re-projection. Non-harmonic Fourier reconstruction.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 40 / 41

Summary

References

Concentration Method

1 A. Gelb and E. Tadmor, Detection of Edges in Spectral Data, in Appl. Comp.

Harmonic Anal., 7 (1999), pp. 101–135.

2 A. Gelb and E. Tadmor, Detection of Edges in Spectral Data II. Nonlinear

Enhancement, in SIAM J. Numer. Anal., Vol. 38, 4 (2000), pp. 1389–1408. Statistical Formulation

1 A. Viswanathan, D. Cochran, A. Gelb and D. Cates, Detection of Signal

Discontinuities from Noisy Fourier Data, in Conf. Record of the 42nd Asilomar Conf.

• n Signals, Systems and Comp., Oct. 2008

Iterative Formulations

1 E. Tadmor and J. Zou, Novel edge detection methods for incomplete and noisy

spectral data, in J. Fourier Anal. Appl., Vol. 14, 5 (2008), pp. 744–763. Spectral Re-projection

1 D. Gottlieb and C.W. Shu, On the Gibbs phenomenon and its resolution, in

SIAM Review (1997), pp. 644–668.

Aditya Viswanathan (Arizona State University) Edge Detection from Spectral Data Oct 22 2009 41 / 41