Transcript Document
Multiscale Analysis of Photon-Limited
Astronomical Images
Rebecca Willett
Photon-limited astronomical imaging
NG2997
Saturn
Richardson-Lucy performance on
Saturn deblurring
Error performance of
standard R-L algorithm
Error performance of R-L
algorithm with regularization
Iteration Number
QuickTime™ and a
decompressor
are needed to see this picture.
Main question: how to best perform
Poisson intensity estimation?
Test data
Saturn
Rosetta (Starck)
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
Saturn image
Wavelet
coefficients of
Saturn image
Approximation using
wavelet coeffs. > 0.3
Wavelet coefficient magnitude
Wavelet thresholding
QuickTime™ and a
TIFF (LZW) decompressor
are needed to see this picture.
Sorted wavelet index
Noisy Saturn
image
Wavelet
coefficients of
Noisy Saturn
image
Estimate using wavelet
coeffs. > 0.3
Noise wavelet coefficient magnitude
Wavelet thresholding for denoising
QuickTime™ and a
TIFF (LZW) decompressor
are needed to see this picture.
Sorted wavelet index
Translation invariance
1. Approximate with
Haar wavelets as
on previous slide
1. Shift image by 1/3 in
each direction
2. approximate as before
3. shift back by 1/3
Avoid this difficulty by averaging over all different possible shifts;
this can be done quickly with undecimated (redundant) wavelets
Wavelet thresholding results
Haar
wavelets
Variance stabilizing transforms
QuickTime™ and a
TIFF (LZW) decompressor
are needed to see this picture.
Anscombe 1948
QuickTime™ and a
TIFF (LZW) decompressor
are needed to see this picture.
QuickTime™ and a
TIFF (LZW) decompressor
are needed to see this picture.
QuickTime™ and a
TIFF (LZW) decompressor
are needed to see this picture.
Anscombe transform results
Haar
wavelets
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
Corrected Haar threshold results
Multiplicative Multiscale Innovation Models
(aka Bayesian Multiscale Models)
0,0,0 X0,0,0
1,0,0 X1,0,0
1,1,0 X1,1,0
1,0,1 X1,0,1
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
Timmermann & Nowak, 1999
Kolaczyk, 1999
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,0 X1,0,0 1,1,0 X1,1,0 1,0,1 X1,0,1 1,1,1 X1,1,1
MMI-MAP estimation results
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
0,0,0 X0,0,0
1,0,0 X1,0,0 1,1,0 X1,1,0 1,0,1 X1,0,1 1,1,1 X1,1,1
MMI - Complexity Regularization
Kolaczyk & Nowak, 2004
MMI - Complexity Regularization
pruning = aggregation = data fusion = robustness to noise
Partitions selection
penalty
(prior)
likelihood
|P|
Complexity penalized estimator:
set of all possible partitions
MMI-Complexity regularization results
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
Platelet results
á 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 1989
Starck 2002
Intensity estimation with á trous wavelets
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 á trous coefficient
Use level-dependent,
wavelet-dependent,
location-dependent
thresholds
(result on next slide)
á trous results
Truth
Observations; 1.74
Wavelet thresholding; 0.325
Wavelets + Anscombe; 0.465
Corrected thresholds; 0.198
MMI - MAP; 0.245
MMI - Complexity Reg.; 0.173
Platelets; 0.163
Observations
Wavelet thresholding
Wavelets + Anscombe
Corrected thresholds
MMI - MAP
MMI - Complexity Reg.
Platelets
A trous
Observations
Wavelet thresholding
Wavelets + Anscombe
Corrected thresholds
MMI - MAP
MMI - Complexity Reg.
Platelets
A trous
Method
Speed
Effectiveness
Wavelet
thresholding
Wavelets +
Anscombe
Fast
Poor
Fast
Poor
Corrected
thresholds
MMI-MAP
Fast
Medium
Fast
Medium
MMI-EMC2
Medium
High; significance maps!
MMI-Complexity
regularization
Platelets
Fast
High
A trous
Medium
Medium-slow High
High