Blind galaxy survey images Deconvolution with Shape Constraint - - PowerPoint PPT Presentation

Jean-Luc Starck

http://jstarck.cosmostat.org Collaborators:

Blind galaxy survey images Deconvolution with Shape Constraint

Morgan Schmitz Fred Nole

Fadi Nammour

• Part I: Introduction to the Blind Deconvolution Problem
• Part II: Point Spread Fonction (PSF) Field Recovery
• Part III: Shape Measurements & Deconvolution

Weak Gravitational Lensing & Blind Deconvolution

Weak Lensing

Observer Gravitational lens Background galaxies

Euclid Mission

• Medium-class mission in ESA’s Cosmic Vision program
• 6 year survey, launch in Q4 2022

Optimized for weak lensing (and galaxy clustering) — 15,000 deg2 survey area — over 1.5 billion galaxies — redshifts out to z = 2

Euclid ESA Space Mission

2022

probe the properties and nature of dark energy,

dark matter, gravity and distinguish their effects decisively Gains in space:

Stable data: homogeneous data set over the whole sky Systematics are small, understood and controlled Homogeneity : Selection function perfectly controlled

Galaxies

Detection + Classification stars/galaxies Galaxies Stars

Motivation for spatial observations

Convolution Operator + Sampling

Space Variant PSF

Euclid PSF Wavelength Dependency

• PSF Modeling

➡Undersampling ➡Space dependency ➡Wavelength dependency ➡Time dependency

Shape Measurements

Galaxies are convolved by an asymetric PSF + Images are undersampled

Shape measurements must be deconvolved

Methods: Moments (KSB), Shapelets, Forward-Fitting, Bayesian estimation, etc

State of the art: PSFextractor [E. Bertin, 2011]

where h is Lanczos interpolation kernel

Regularization

• 1. Forward Modelling approach (FM) leaded by Lance Miller
• Model the exit pupil using pre-defined Zemax modes.
• Fit to all stars.

➡ PRO: No other existing method can achieve the very strong requirement on the PSF. ➡ CON: But hard to validate since i) the simulations are done with the same model, and ii) the model may change (vibrations at launch).

• 2. Need a data driven approach
• Validate the FM solution on real data.
• All surveys ((KIDS, CFHTLenS, DES) have used the data driven approach.
• HST: TinyTim HST modelling software (Krist 1995) for the Hubble Space Telescope not

as good as a data driven approach (Jee et al, PASP, 2007; Hoffmann and Anderson, Instrument Science Report ACS 2017-8, 2017).

• What does not exist can be invented.
• Combination of both approaches could lead to the optimal solution.

State of the art: PSF Modeling

–Introduction to the Blind Deconvolution Problem –Euclid PSF Field Measurement

Monochromatic PSF: Matrix Factorization + Laplacian Graph + Sparsity Polychromatic PSF and Optimal Transport

–Shape Measurements & Deconvolution

Blind Deconvolution

• F. M. Mboula, J.-L. Starck, S. Ronayette, K. Okumura, and J. Amiaux, Super-resolution method using sparse regularization for

point-spread function recovery. A&A, 575, id.A86, 2015.

• F. Ngole, J.-L Starck, et al, “Constraint matrix factorization for space variant PSFs field restoration”, Inverse Problems, 2016.
• M.A. Schmitz, J.-L. Starck and F.M. Ngolè, "Euclid Point Spread Function field recovery through interpolation on a Graph

Laplacian", submitted, 2018.

• F.M. Ngolè Mboula, and J.-L. Starck, "PSFs field learning based on Optimal Transport distances", SIAM

Journal on Imaging Sciences, 10, 3, pp. 979-1004, 2017. DOI: 10.1137/16M1093677.

• M. A. Schmitz et a; "Wasserstein Dictionary Learning:Optimal Transport-based unsupervised

representation learning", SIAM Journal on Imaging Sciences, 11, 2018.

Decimation Random shifting Noise

Datapoints Pixels K

. . .

N

matrix

Degradation

Pixels K

. . .

N/4

matrix

Flat to 1-D vector Observed PSF Datapoints

. . . . . . . . . . . . . . . . . .

“True” PSF (punctual star convolved with true PSF) Datapoint Datapoint

The PSF Field Recovery Problem

: true PSFs at different positions of the FOV. : decimation and shifting. : observed PSFs. X ∈ RN×K

+

min

X

Y MX 2 X

matrix M

The PSF Super-resolution Challenge

Y = HX + N

Inverse Problems & Regularization

Physical Knowledge on X (ex: Gaussian Random Field, etc). ==> Gaussian smoothing, Wiener reconstruction, etc ==> Log normal distribution prior Knowledge on the histogram of X in pixel space or in another one. ==> Positivity constraint, sparsity constraint, etc. X properties are understood through a representative data set. ==> Machine Learning

min

X ||Y − HX||2

C(X) +

ai,k = coefficient corresponding to contribution

• f the i-th vector to the k-th PSF.

PSF(k) = xk =

r

X

i=1

ai,ksi

si = ith vector (2D image)

Eigen PSF

Joint estimations of super-resolved PSFs at stars positions

➡

Low rank constraint: Constraint the PSFs to be a linear combination of the eigenvectors PSFs:

Constraints

a1,i

a2,i

aK,i

ak,i

a.,i a.,i ui ∈ Ustars

f( uk1 uk2 2) ➡

Proximity constraint on the ai,k coefficients: the closer are the stars, the more the coefficients of the linear combination are similar.

X

r

k Φsi k1

➡

Smoothness constraint on each si

➡

Positivity constraint (xk > 0)

PSF Laplacian Graph

Laplacian Graph P = D - J, where D is the degree matrix and J is the adjacent matrix of the graph

P tP = V ΛV t

Degree matrix= diagonal matrix with information about the degree of each vertex—that is, the number of edges attached to each vertex. Adjacent matrix A: element Aij is one when there is an edge from vertex i to vertex j, and zero when there is no edge.

si are ”eigen PSF”

Matrix Factorization

PSF(k) = xk =

r

X

i=1

ai,ksi

S = [s1, ..., sp] http://www.cosmostat.org/software/rca/

Each eigen PSF is sparse in a wavelet dictionary

X

r

k Φsi k1

MODEL:

Matrix Factorization

http://www.cosmostat.org/software/rca/

Sparsity on the eigen PSF: the PSF should have a sparse representation in an appropriate basis Positivity Constraint Smoothness of the PSF field constraint : the

smaller the difference between two PSFs’ positions ui , uj , the smaller the difference between their estimated representations

Data fidelity term

min

S,α

1 2Y M(SαV )2

2 + r

• i=1

wi Φsi1 + ι+(SαV ) + ιΩ(α)

The PSF Interpolation Challenge

ui ∈ Ustars 1 uk1 uk2 ek

2

Numerical Experiments

With undersampling (upsampling factor of 2)

Center Local

Corner

Obs Ref RCA PSFEX

Data: 500 Euclid-like PSFs (Zemax), field observed with different SNRs

Theese PSFs account for mirrors polishing imperfections, manufacturing and alignments errors and thermal stability of the telescope.

Numerical Experiments

Stars SNR

Multiplicative Bias

• M.A. Schmitz, J.-L. Starck and F.M. Ngolè, "Euclid Point Spread Function field recovery through

interpolation on a Graph Laplacian", submitted, 2018.

20% Improvement

PSF Wavelength Dependency

• G. Monge, “Mémoire sur la théorie des déblais et des remblais”, 1781

Optimal Transport

The Monge-Kantorovich problem

Barycenter OT Linear interpolation Plausible average distribution in presence of shifts and changing of shape

Optimal Transport & PSF

Euclidean barycenter

,

Optimal Transport

If instead we replace the Euclidean distance by the Wasserstein distance..

Wasserstein barycenter

,

Optimal Transport

1

Wasserstein barycenter

,

ATTRACTIVE TOOL, BUT COSTLY TO COMPUTE: ==> fixed using a regularisation (Sinkhorn iterations)

• M. A. Schmitz et a; "Wasserstein Dictionary Learning:Optimal Transport-based unsupervised representation learning", SIAM Journal on Imaging

Sciences, 11, 2018.

Optimal Transport

Back to the PSF estimation problem: … … … … … … …

How to break the degeneracy?

1.0 …

Hypothesis: monochromatic PSFs of a intermediary wavelength are Wasserstein barycenters of the extremes.

Optimal Transport

Hypothesis: monochromatic PSFs of a intermediary wavelength are Wasserstein barycenters of the extremes. Experiment results:

Simulated Euclid- like PSFs (ground truth) OT Modeling

linear projection of wavelengths into [0,1].

• RCA

λ

Simple constraint on dictionary: probability distribution. Eigen-PSF extreme right Eigen-PSF extreme left Reconstructed PSFs Ground-truth

Very Preliminary Results

• M. A. Schmitz et al, in preparation, 2018

Numerical Experiments

• Part I: Weak Lensing and Euclid
• Part II: Point Spread Fonction (PSF) Field Recovery
• Part III: Shape Measurements & Deconvolution

Weak Gravitational Lensing: The prospects of Euclid

Astronomical Image Deconvolution

y = Hx + n

Observed Image PSF Convolution True Image Noise

Sandard deconvolution framework: H is huge !!!

argmin

X

1 2kY HXk2

2 + kΦtXkp

s.t. X 0

Sandard deconvolution framework:

argmin

X

1 2kY H(X)k2

2 + λkΦtXkp

s.t. X 0

Big Astronomical Image Deconvolution

Y = H(X) + N

[y0, y1, …, yn] [H0x0, H1x1, …, Hnxn] [n0, n1, …, nn]

Object Oriented Deconvolution For each galaxy, we use the PSF related to its center pixel:

min

α kY HΦαk2 + λkαkp p

➡ Wavelet Denoising introduces a bias. ➡ Fitting a galaxy model directly on significant wavelet coefficients introduces a bias. ➡ Need to find a way to preserve the ellipticity during the restoration process

Inverse Problems: Sparse Recovery

Measuring Galaxies shapes

a b θ

: Galaxy

X

e(X) = 1 − b

a

2 1 + b

a

2 exp(2iθ) = e1(X) + ie2(X)

The complex ellipticity of a galaxy image uses quadrupole moments

e = γ + eg < e > γ

and

X ∈ Mn(R)

µs,t(X) = X

i,j

X[i, j](i − ic)s(j − jc)t

a b θ a b θ

From Ellipticity to Moments

is an unbiased estimator of “ellipticity” (Schneidez & Seitz 1994)

γ1(X) = µ2,0(X) − µ0,2(X) µ2,0(X) + µ0,2(X)

γ2(X) = 2µ1,1(X) µ2,0(X) + µ0,2(X)

γ = γ1 + iγ2 = µ2,0 − µ0,2 + 2µ1,1 µ2,0 + µ0,2

• <latexit sha1_base64="(nul)">(nul)</latexit><latexit sha1_base64="(nul)">(nul)</latexit><latexit sha1_base64="(nul)">(nul)</latexit><latexit sha1_base64="(nul)">(nul)</latexit>

=

Observed galaxy Matched Gaussian window Windowed galaxy

Y

W

W Y

Moments Windowing

Does not work in presense of noise! Need to apply a window function.

• “HSM” (or adaptive moments, Hirata & Seljak 2003; Mandelbaum et al 2004): match the gaussian window to the object
• Handling the PSF: Kaiser, Squires & Broadhurst 1995 (KSB).

High Accuracy Requirements

˜ γ = (1 + m)γ + c γ = γ1 + iγ2 = |γ|ei2Φ m = m1 + im2 = |m|ei2Φm c = c1 + ic2 = |c|ei2Φc σm < 2 × 10−3 σc < 5 × 10−4

Requirements for the ESA Euclid Space Telescope Additive bias Multiplicative bias

Criteria for a good method

➡ Fast Calculation. ➡ Level of bias (multiplicative and additive bias). ➡ Capacity to measure the bias (number of simulations, depth, resolution, etc). ➡ Robustness to calibration errors (PSF measurement errors, detectors effects, background subtraction, etc).

X = image properties

e1(X) = µ2,0(X) − µ0,2(X) µ2,0(X) + µ0,2(X) , e2(X) = 2µ1,1 µ2,0(X) + µ0,2(X)

X ∈ Mn(R)

µs,t(X) = X

i,j

X[i, j](i − ic)s(j − jc)t

Moments using Inner Products

e1(X) = hX, U5ihX, U3i hX, U1i2 + hX, U2i2 hX, U5ihX, U3i hX, U1i2 hX, U2i2

, e2(X) = 2(hX, U6ihX, U3i hX, U1ihX, U2i) hX, U5ihX, U3i hX, U1i2 hX, U2i2

U1 =      1 1 . . . 1 2 2 . . . 2 . . . . . . ... . . . n n . . . n     

U2 =      1 2 . . . n 1 2 . . . n . . . . . . ... . . . 1 2 . . . n     

U3 =      1 1 . . . 1 1 1 . . . 1 . . . . . . ... . . . 1 1 . . . 1     

U4 =      12 + 12 12 + 22 . . . 12 + n2 22 + 12 22 + 22 . . . 22 + n2 . . . . . . ... . . . n2 + 12 n2 + 22 . . . n2 + n2     

U5 =      12 − 12 12 − 22 . . . 12 − n2 22 − 12 22 − 22 . . . 22 − n2 . . . . . . ... . . . n2 − 12 n2 − 22 . . . n2 − n2     

Matrices of indices combination

U6 =      1 2 . . . n 2 4 . . . 2n . . . . . . ... . . . n 2n . . . n2     

e1(X) = hX, U5ihX, U3i hX, U1i2 + hX, U2i2 hX, U5ihX, U3i hX, U1i2 hX, U2i2 , e2(X) = 2(hX, U6ihX, U3i hX, U1ihX, U2i) hX, U5ihX, U3i hX, U1i2 hX, U2i2

Idea : Conserving the inner products is equivalent to conserving the ellipticity parameters. Advantage : These inner products are linear functions of the image we want to restore so it is mathematically easier to work with. Question : How to formulate a constraint using these inner products? What would be the additional contribution of this constraint regarding the data fidelity term?

Ellipticities as inner products

The Shape Contraint : The constraint looks just like a data fidelity term expressed in the moments space. Idea :

• Remove noise as much as possible inside the constraint in a inverse pb.

M(X) =

6

X

i=1

µi hX ⇤ H Y, Uii2

The Moment Shape Constraint

M(X) =

6

X

i=1

µi hW (X ⇤ H Y ), Uii2

• <latexit sha1_base64="(nul)">(nul)</latexit><latexit sha1_base64="(nul)">(nul)</latexit><latexit sha1_base64="(nul)">(nul)</latexit><latexit sha1_base64="(nul)">(nul)</latexit>

=

Observed galaxy Matched Gaussian window Windowed galaxy

Y

W

W Y

Formulation with a Gaussian window :

Constraint with a Gaussian Window

Multi-Window Constraints

Rather than trying to fit a Gaussian window to the data, why not trying all windows in all directions and all scales ? This can be done easily using a Curvelet like Decomposition.

Contrast Enhancement

Formulation with shearlets : Aleternative form :

• : adjoint of the shearlet transform operator

Ψ∗

j

• : Curvelet/Shearlet transform (Kutyniok & Labate, 2012)

Ψj

M(X) = 1 6

6

X

i=1

1 N

N

X

j=1

µij ⌦ X ∗ H − Y, Ψ∗

j(Ui)

↵2

M(X) = 1 6

6

X

i=1

1 N

N

X

j=1

µij hΨj(X ⇤ H Y ), Uii2

Advantage : is a constant.

Ψ∗

j(Ui)

Measuring the Multi-Scale Moments

• Number of bands

• : galaxy observation

Y ∈ Mn(R)

• : noise standard deviation in

σ

Y

• : trade-off parameter between the data fidelity term and the moments constraint

γ

• : weighting tensor

λκ−σMAD

• : starlet transform operator

Φ

• : moments constraint

M(.)

L(X) = 1 2σ2 kX ⇤ H Y k2

F +

γ 2σ2 M(X) + kλκ−σMAD ΦXk0

• : Point Spread Function

H ∈ Mn(R)

Functional to Minimise

Denoising Experiment: Gaussian window v/s shearlets

300 GALSIM galaxies

Gaussian window v/s shearlets

Experiment: Deconvolution

100 GALSIM galaxies + Optical GALSIM PSF

Experiment: Deconvolution

Deconvolution : GALSIM galaxies + MeerKat PSF

1024 GALSIM galaxies + Radio PSF (MeerKat, CASA)

Weak Lensing and Shape Measurements ➡ Serious mathematical challenges to well control systematics. ➡ Need the best methods.

Euclid PSF: RCA is new method which uses all available PSFs to derive the PSF field.

➡ Use Optimal Transport, Graph Laplacian and Sparsity. Moments Constraint: ➡ A new approach to preserve the ellipticity information during the restoration process. ➡ Shape constraints using shearlets outperform the standard Gaussian windowing. Perspectives: ➡ Radio: Is this shape constraint reconstruction method competitive with galaxy model fitting in Fourier space ? ➡ Optical: Compare the constraint denoising method to the Gaussian windowing.

Conclusions