Sparse Decompositions in Dictionaries for Interferometric Image - - PowerPoint PPT Presentation

Sparse Decompositions in Dictionaries for Interferometric Image Reconstruction

AIP 2009, Vienna David Mary

1

Overview

1 Radio/Optical interferometry : an undetermined inverse problem 2 Sparsity and compressiblity of real images 3 On the sparse approximation problem in redundant dictionaries 4 Practical recovery algorithms in transformed redundant dictionaries 5 Simulation results

2

• 1. Interferometry in Astronomy

Some examples in radio / optical wavelengths VLA VLTI

3

• 1. Interferometry in Astronomy

figure from [S. Meimon, 2005]

Pupil :{T1 T2 T3} Transfer Function (MTF)

• The pupil is diluted
• Transfer function = autocorrelation of the pupil
• Interferometry provides data related to samples of the Fourier spectrum of

the scene ⇒ uncomplete sampling of the Fourier spectrum

• Assumption in the following : Fourier phases are available (this is the case

in radio interferometry, but currently not in optical interferometry)

4

• 1. Model : an underdetermined problem
• Underdetermined system :

y = Fν x (+ n) , x ∈ R+N : image of interest (unknown) Fν : Fourier transform of x restricted to the set {ν} of probed frequencies y ∈ CM : data points : Fourier spectrum of x sampled at frequencies {ν}

• N > M → Infinity of solutions in general
• Underdetermination can be decreased or even removed if x is compressible
• r sparse
• Sparsity is expressed via representation bases or dictionaries (e.g. union of

bases)

5

• 2. Example : 2D DCT matrix

Best non-linear approximation in DCT : snr = 17.99 dB ; in I : snr = 0.34 dB

6

• 2. Sparsity, compressibility and approximation
• x ∈ RN is S-sparse if supp(x) = ||x||0 = S
• x is compressible if it is well approximated by a few vectors of some repre-

sentation dictionary D = {Di}i=1,...,T with coefficients ui : x ≈

• few i

uiDi

• Best M-term approximation of x in D :

find support Λ (|Λ| = M) minimizing ||x − xΛ||2 = ||x −

• i∈Λ

uiDi||2

• Traditional source coding uses orthonormal bases in which best M-term ap-

proximation of x is obtained by retaining M = |Λ| largest coefficients xΛ =

• i∈Λ⇔ui>T

uiDi

• Highly irregular signals require dictionaries of vectors larger than bases

7

• 2. Redundant dictionaries
• Large variety of available representations, e.g.
• Direct space (W = I)
• Discrete Cosine Transform
• Wavelets
• Curvelets [Starck 2002]
• Bandlets [Peyré & Mallat 2005]
• The choice of a representation is made w.r.t. a class of signals

and for images depends on the existence of fast operators

• Can concatenate representations Wi in a dictionary of T > N vectors, then

x ≈ Du with D = [W1 W2 . . . WK], and u ∈ RT is sparse

• Redundant dictionaries increase the richness of the geometrical a priori

but make the best M-term approximation problem difficult

• Best M-term approximation of x in D requires solving

LO(L, x, Λ) = ||x −

• i∈Γ

uiDi||2 + L2||u||0

8

• 2. Sparsity Promoting Reconstruction
• Model :

y = Fν D

Dν

u (+ n), with Du = x.

• Dν : transformed operator
• The transformed dictionary has vectors {FνDi}i=1,...,T.
• Reconstruction approach :
• Choose D
• Find a few {˜

ui} that best approximate y

• The reconstructed image is ˜

x = D˜ u

9

• 3. Recovery of S-sparse signal
• Assume

y = Dν u0 with ||u0||0 = S.

• To find u0, try to find u∗ = argminu||u||0 s.t. Dνu = y
• The solution of this problem is unique if any set of at least 2S columns of

Dν are independent ⇒ Crucial point : columns dependencies (cf mutual coherence, RIP , ERC,...)

• The direct minimisation of |u0| is generally untractable (e.g. C26

100 ≈ N)

• One alternative : try to solve the relaxed problem :

minu ||u||p s.t. Dνu = y with ||u||p = (

• i

|ui|p)1/p for some 0 < p ≤ 1

10

• 3. Best approximation of compressible signals
• Model : y = Dνu + r
• First approach (here) : Synthesis prior

⇒ Find the best M approximation vectors in Dν by minimizing on u : Ls

O(L, x, Λ) = ||y −

• i∈Γ

uiDνi||2 + L2||u||0

• The coefficients {˜

ui}i∈˜

Λ can be estimated by a sparse denoising estimation

algorithm

• The synthesised image is then ˜

x =

i∈˜ Λ ˜

uiDi

• This corresponds to a sparse approximation of the data : ˜

y =

i∈˜ Λ ˜

uiDνi

11

• 3. Best approximation of compressible signals
• Model : y = Dνu + r
• Second approach : analysis prior

⇒ do not focus on y but on x ≈

• few i

uiDi

• To find the best M approximation vectors in D, minimize on x′ :

La

O = ||y − Fνx′||2 + L2||D∗x′||0

instead of Ls

O = ||y − FνDu||2 + L2||u||0

• The solution minimizing La

O and Ls O are not the same for redundant dictio-

naries [Elad, Milfar & Rubinstein 2007]

• The synthesis coefficients {˜

ui}i∈˜

Λ can also be estimated by sparse denoi-

sing estimation algorithms (e.g. MCA [Bobin et al. 2007])

12

• 4. Practical Recovery Algorithms
• Direct minimization of La

0 and Ls 0 generally untractable

• Can try minimization of lp, 0 < p < 1 (IRLS,..) but risk of local minima
• Many approaches to minimise l1 or other sparsity promoting functionals
• Greedy pursuits
• The following only implements solutions focused on the synthesis problem :

Ls

0 = ||y − FνDu||2 + L2||u||0

13

• 4. Matching pursuit [Mallat & Zhang 1993]
• Algorithm :
• Initialisation : Choose D. m = 0. R0y = y. Compute {< y, Dνi >}i∈Γ
• Best match : Find Dνm = argmax < Rmy, Dνi >
• Update : Rm+1y = Rmy − PDνm(Rmy) = Rmy− <Rmy, Dνm> Dνm
• Criterion stop : Compare ||Rm+1y|| to noise, or analysis coefficients to

threshold

• Back projection : Reestimate ˜

uii∈˜

Λ by minimizing ||y −

• i∈˜

Λ

u′

iDνi||2

• Can be run for several dictionaries

Here : D include combinations of Diracs, wavelets and DCT bases

• The CLEAN algorithm [Hoegbom 1974] is equivalent to MP with D = I

14

• 4. Orthogonal Matching pursuit
• The residual Rm+1y is orthogonal to Dνm but not {Dνi}i=0,...m−1
• The update step Rm+1y = Rmy− <Rmy, Dνm> Dνm reintroduces com-

ponents of the {Dνi}i=0,...m−1

• Orthogonalisation in OMP : the update step is

Rm+1y = Rmy − PV (Rmy) where V is the space generated by {Dνi}i=0,...m−1.

• Same stopping criteria as MP

15

• 4. Iterative Soft Thresholding Algorithm
• Goal : solve the constrainted problem

˜ u = argminu||u||1 s.t. ||y − Dνu||2 ≤ ǫ by minimizing its Langrangian formulation La

1 = ||y − FνDu||2 + L2||u||1

• For each ǫ, there exists L making both problems equivalent
• Algorithm :
• Initialisation : m = 0. Choose u0 = 0.
• Gradient step : ˜

u

′m = ˜

um + 1

τ D∗ ν(y − Dν˜

um)

• Soft thresholding : ˜

uim+1 = ρL/τ( ˜ ui

′m)

• Criterion stop : A fixed tolerance on ||˜

um+1 − ˜ um||1

• Back projection : Reestimate ˜

uii∈˜

Λ by minimizing ||y −

• i∈˜

Λ

˜ uiDνi||2

• To reach ǫ, iteratively update the threshold [Chambolle 2004]
• Many acceleration algorithms exist (e.g. FISTA,GPSR,SPARSA,GPAS)

16

• 5. Simulations

17

• 1. Dictionary = Dirac basis (data snr =29.3 dB)

Original Image MTF Equivalent PSF Peudo Inv Rec snr = 0.90884dB

18

(data snr =29.3 dB)

Original Image Peudo Inv Rec snr = 0.90884dB CLEAN !=0.1 snr = 12.4508dB minl1 I snr = 16.8484dB

19

• 2. Dictionary = [Dirac Wavelet] bases (data snr =39 dB)

Original Image MTF Equivalent PSF Peudo Inv Rec snr = 9.1512dB

20

(data snr =39 dB)

21

The reconstructed image on the left is the sum of two components : one image sparse in the Dirac basis plus one image sparse in the Wavelet (not DCT) basis

• 3. Dictionary = [Dirac DCT] bases (data snr =43.1 dB)

Original Image MTF Equivalent PSF Peudo Inv Rec snr = 15.1911dB

22

(data snr =43.1 dB)

minl1 [DCT I] snr = 43.8411dB I component DCT component

The reconstructed image on the left is the sum of two components : one image sparse in the Dirac basis plus one image sparse in the DCT basis

23

• 3. Real image (Abell cluster). Dictionary = [Dirac Wavelet] bases

(data snr =37.5 dB)

Original Image MTF Equivalent PSF Peudo Inv Rec snr = 11.9241dB

24

(data snr =37.5 dB)

minl1 [W I] snr = 14.9395dB I component W component

25

• 6. Conclusions
• Modeled the astronomical image reconstruction (deconvolution) problem as

a sparse approximation problem in redundant dictionaries

• l1 minimization generally improves the reconstruction w.r.t. standard greedy

approaches (MP (CLEAN), OMP)

• The considered redundant dictionaries generally improve on single repre-

sentation bases

• Model allows to include efficient and flexible prior geometrical information

for the reconstruction

• Improvements : include other representations bases in the dictionary (cur-

velets, bandlets), compare synthesis and analysis approaches, compare to

• ther regularized deconvolution methods (e.g. l2 − l1)

26