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 : { T 1 T 2 T 3 } 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 ∈ C M : 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 or 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 ∈ R N 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 = { D i } i =1 ,...,T with coefficients u i : � x ≈ u i D i few i • Best M -term approximation of x in D : find support Λ ( | Λ | = M ) minimizing || x − x Λ || 2 = || x − � u i D i || 2 i ∈ Λ • Traditional source coding uses orthonormal bases in which best M -term ap- proximation of x is obtained by retaining M = | Λ | largest coefficients � x Λ = u i D i i ∈ Λ ⇔ u i >T • 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 W i in a dictionary of T > N vectors, then x ≈ Du with D = [ W 1 W 2 . . . W K ] , and u ∈ R T 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 � u i D i || 2 + L 2 || u || 0 L O ( L, x, Λ) = || x − i ∈ Γ 8

2. Sparsity Promoting Reconstruction • Model : y = F ν D u (+ n ) , Du = x . with � �� � D ν � �� � D ν : transformed operator • The transformed dictionary has vectors { F ν D i } i =1 ,...,T . • Reconstruction approach : • Choose D • Find a few { ˜ u i } that best approximate y • The reconstructed image is ˜ x = D ˜ u 9

3. Recovery of S -sparse signal • Assume with || u 0 || 0 = S . y = D ν u 0 • To find u 0 , try to find u ∗ = argmin u || u || 0 s.t. D ν u = y • The solution of this problem is unique if any set of at least 2 S columns of D ν are independent ⇒ Crucial point : columns dependencies (cf mutual coherence, RIP , ERC,...) • The direct minimisation of | u 0 | is generally untractable (e.g. C 26 100 ≈ N ) • One alternative : try to solve the relaxed problem : | u i | p ) 1 /p for some 0 < p ≤ 1 � min u || u || p s . t . D ν u = y with || u || p = ( i 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 : u i Dν i || 2 + L 2 || u || 0 � L s O ( L, x, Λ) = || y − i ∈ Γ • The coefficients { ˜ u i } i ∈ ˜ Λ can be estimated by a sparse denoising estimation algorithm x = � • The synthesised image is then ˜ Λ ˜ u i D i i ∈ ˜ y = � • This corresponds to a sparse approximation of the data : ˜ Λ ˜ u i D νi 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 ≈ u i D i few i • To find the best M approximation vectors in D , minimize on x ′ : O = || y − F ν x ′ || 2 + L 2 || D ∗ x ′ || 0 L a instead of O = || y − F ν Du || 2 + L 2 || u || 0 L s • The solution minimizing L a O and L s O are not the same for redundant dictio- naries [Elad, Milfar & Rubinstein 2007] • The synthesis coefficients { ˜ u i } 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 L a 0 and L s 0 generally untractable • Can try minimization of l p , 0 < p < 1 (IRLS,..) but risk of local minima • Many approaches to minimise l 1 or other sparsity promoting functionals • Greedy pursuits • The following only implements solutions focused on the synthesis problem : 0 = || y − F ν Du || 2 + L 2 || u || 0 L s 13

4. Matching pursuit [Mallat & Zhang 1993] • Algorithm : • Initialisation : Choose D . m = 0 . R 0 y = y . Compute { < y, Dν i > } i ∈ Γ • Best match : Find Dν m = argmax < R m y, Dν i > • Update : R m +1 y = R m y − P Dν m ( R m y ) = R m y − <R m y, Dν m > Dν m • Criterion stop : Compare || R m +1 y || to noise, or analysis coefficients to threshold � u ′ i Dν i || 2 • Back projection : Reestimate ˜ Λ by minimizing || y − u i i ∈ ˜ i ∈ ˜ Λ • 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 R m +1 y is orthogonal to Dν m but not { Dν i } i =0 ,...m − 1 • The update step R m +1 y = R m y − <R m y, Dν m > Dν m reintroduces com- ponents of the { Dν i } i =0 ,...m − 1 • Orthogonalisation in OMP : the update step is R m +1 y = R m y − P V ( R m y ) 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 = argmin u || u || 1 s . t . || y − D ν u || 2 ≤ ǫ ˜ by minimizing its Langrangian formulation 1 = || y − F ν Du || 2 + L 2 || u || 1 L a • For each ǫ , there exists L making both problems equivalent • Algorithm : • Initialisation : m = 0 . Choose u 0 = 0 . ′ m = ˜ u m + 1 τ D ∗ u m ) • Gradient step : ˜ u ν ( y − D ν ˜ ′ m ) u im +1 = ρ L/τ ( ˜ • Soft thresholding : ˜ u i u m +1 − ˜ u m || 1 • Criterion stop : A fixed tolerance on || ˜ � u i D νi || 2 • Back projection : Reestimate ˜ Λ by minimizing || y − u i i ∈ ˜ ˜ i ∈ ˜ Λ • 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 • l 1 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 other regularized deconvolution methods (e.g. l 2 − l 1 ) 26