Multiscale Analysis of Photon-Limited Astronomical Images

1. Multiscale Analysis of Photon-Limited Astronomical Images Rebecca Willett

2. Photon-limited astronomical imaging NG2997 Saturn

3. Error performance of R-L algorithm with regularization Richardson-Lucy performance on Saturn deblurring Error performance of standard R-L algorithm MSE of deconvolvedestimate Iteration Number

4. Main question: how to best perform Poisson intensity estimation?

5. Test data Rosetta (Starck) Saturn

6. Methods reviewed in this talk • Wavelet thresholding • Variance stabilizing transforms • Corrected Haar wavelet thresholds • Multiplicative Multiscale Innovation models • MAP estimation • EMC2 estimation • Complexity Regularization • Platelets • á trous wavelet thresholding

7. Wavelet thresholding Saturn image Wavelet coefficients of Saturn image Wavelet coefficient magnitude Sorted wavelet index Approximation using wavelet coeffs. > 0.3

8. Wavelet thresholding for denoising Noisy Saturn image Wavelet coefficients of Noisy Saturn image Noise wavelet coefficient magnitude Sorted wavelet index Estimate using wavelet coeffs. > 0.3

9. Avoid this difficulty by averaging over all different possible shifts;this can be done quickly with undecimated (redundant) wavelets Translation invariance Approximate with Haar wavelets as on previous slide Shift image by 1/3 in each direction approximate as before shift back by 1/3

10. Haarwavelets Wavelet thresholding results

11. Variance stabilizing transforms Anscombe 1948

12. Haarwavelets Anscombe transform results

13. Kolaczyk’s corrected Haar thresholds Basic idea: Keep wavelet coeffs which correspond to signal;Threshold wavelet coeffs which correspond to noise (or background) If we had Gaussian noise (variance 2) and no signal: (j,k)th Gaussian wavelet coeff. For Poisson noise, design similar bound for background 0 (noise): (j,k)th Poisson wavelet coeff. Threshold becomes: Background intensity level Kolaczyk 1999

14. Corrected Haar threshold results

15. 0,0,0 X0,0,0     1,0,1 X1,0,1 1,0,0 X1,0,0 1,1,0 X1,1,0 1,1,1 X1,1,1 • Recursively subdivide image into squares • Let { denote the ratio between child and parent intensities • Knowing { Knowing { • Estimate {} from empirical estimates based on counts in each partition square Multiplicative Multiscale Innovation Models (aka Bayesian Multiscale Models) Timmermann & Nowak, 1999Kolaczyk, 1999

16. 0,0,0 X0,0,0     1,0,1 X1,0,1 1,0,0 X1,0,0 1,1,0 X1,1,0 1,1,1 X1,1,1 MMI-MAP estimation Basic idea: place Dirichlet prior distribution with parameters {} on {estimate {by maximizing posterior distribution

17. MMI-MAP estimation results

18. 0,0,0 X0,0,0     1,0,1 X1,0,1 1,0,0 X1,0,0 1,1,0 X1,1,0 1,1,1 X1,1,1 MMI-EMC2 • Before (with MMI-MAP): • place Dirichlet prior distribution with parameters {} on { • user sets parameters {} • estimate {by maximizing posterior distribution • Now (with MMI-EMC2): • place hyperprior distribution on parameters {} • user only controls few hyperparameters • prior information about intensity built into hyperprior • use MCMC to draw samples from posterior • Estimate posterior mean • Estimate posterior variance Esch, Connors, Karovska, van Dyk 2004

19. MMI - Complexity Regularization Kolaczyk & Nowak, 2004

20. MMI - Complexity Regularization pruning = aggregation = data fusion = robustness to noise

21. penalty(prior) likelihood |P| Partitions selection Complexity penalized estimator: set of all possible partitions

22. MMI-Complexity regularization results

23. MMI-Complexity regularization theory No other method can do significantly better asymptotically for this class of images! This theory also supports other Haar-wavelet based methods!

24. Platelet estimation Donoho, Ann. Stat. ‘99 Willett & Nowak, IEEE-TMI ‘03

25. Platelet theory No other method can do significantly better asymptotically for this (smoother) class of images! Willett & Nowak, submitted to IEEE-Info.Th. ‘05

26. Platelet results

27. á trous wavelet transform 1. Redefine wavelet as difference between scaling functions at successive levels 2. Compute coeffs. at one level by filtering coeffs at next finer scale 3. This means synthesis (getting image back from wavelet coeffs.) is simple addition Holschneider 1989Starck 2002

28. Method 1(Classical) Compute Anscombe transform of data Perform á trous wavelet thresholding as if iid Gaussian noise (same problems as other Anscombe-based approaches for very few photon counts) Method 2(Starck + Murtagh, 2nd ed., unpublished) Compute variance stabilizing transform of each á trouscoefficient Use level-dependent, wavelet-dependent, location-dependent thresholds (result on next slide) Intensity estimation with á trous wavelets

29. á trous results

30. Observations; 1.74 Truth Wavelet thresholding; 0.325 Wavelets + Anscombe; 0.465 Corrected thresholds; 0.198 MMI - MAP; 0.245 MMI - Complexity Reg.; 0.173 Platelets; 0.163

31. Wavelet thresholding Observations Wavelets + Anscombe Corrected thresholds MMI - MAP MMI - Complexity Reg. Platelets A trous

32. Wavelet thresholding Observations Wavelets + Anscombe Corrected thresholds MMI - MAP MMI - Complexity Reg. Platelets A trous