330 likes | 347 Views
Multiscale Analysis of Photon-Limited Astronomical Images. Rebecca Willett. Photon-limited astronomical imaging. NG2997. Saturn. Error performance of R-L algorithm with regularization. Richardson-Lucy performance on Saturn deblurring. Error performance of standard R-L algorithm.
E N D
Multiscale Analysis of Photon-Limited Astronomical Images Rebecca Willett
Photon-limited astronomical imaging NG2997 Saturn
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
Main question: how to best perform Poisson intensity estimation?
Test data Rosetta (Starck) Saturn
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
Wavelet thresholding Saturn image Wavelet coefficients of Saturn image Wavelet coefficient magnitude Sorted wavelet index Approximation using wavelet coeffs. > 0.3
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
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
Haarwavelets Wavelet thresholding results
Variance stabilizing transforms Anscombe 1948
Haarwavelets Anscombe transform results
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
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
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
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
MMI - Complexity Regularization Kolaczyk & Nowak, 2004
MMI - Complexity Regularization pruning = aggregation = data fusion = robustness to noise
penalty(prior) likelihood |P| Partitions selection Complexity penalized estimator: set of all possible partitions
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!
Platelet estimation Donoho, Ann. Stat. ‘99 Willett & Nowak, IEEE-TMI ‘03
Platelet theory No other method can do significantly better asymptotically for this (smoother) class of images! Willett & Nowak, submitted to IEEE-Info.Th. ‘05
á 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
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
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
Wavelet thresholding Observations Wavelets + Anscombe Corrected thresholds MMI - MAP MMI - Complexity Reg. Platelets A trous
Wavelet thresholding Observations Wavelets + Anscombe Corrected thresholds MMI - MAP MMI - Complexity Reg. Platelets A trous