Image Deconvolution by S03OBzbO


									Image Deconvolution

Applied to CCD astronomical images
     with astrometric purposes

                      Maria Teresa Merino Espasa
                      Universidad de Barcelona (Spain)
   Introduction to Image Restoration
   Adaptative deconvolution with wavelets
   Applications of our deconvolution algorithm with
    astrometric purposes:
          Increase of faint detections in Baker-Nunn images for faint
           asteroid discovering
          Quest aplication to increase astrometrical resolution to detect
           gravitational lenses candidates
          Work in progress to increse faint detection + increase
           astrometrical resolution

   Conclusions
   1. Introduction to image

Presentation of our algorithm for MLE deconvolution in
linear space-invariant PSF systems with complete basic
CCD images statistics (Poisson + Gaussian)
                     Image restoration
Our detected data images are not exactly the radiation from sources of
interest. It’s important to consider that may be affected by:
             atmosphere effects
             optical distorsions and aberrations
             different pixel sensitivity
             additional noise (in emision, travelling and/or detection)
             …

Image Restoration is the removal or reduction of those degradations that were
incurred while the digital image was being obtained to bring the image toward
what it would have been if it had been recorded without degradation.

In astronomy, the most usual objectives are:
                  deblurring
                  removal atmospheric seeing radiation
                  image sharpening
                  fusion of data from diferent instruments
                  increase S/N of sources
                  astrometric resolution improvement
                  …
Linear image formation system with additive noise
                                                               p(x,y): projection (measured) data
 p( x, y)  g ( x, y)  n( x, y)                               g(x,y): blurred image
                                                               n(x,y): additive noise
   g ( x, y )    h( x, y; x1 , y1 )·d ( x1 , y1 )dx1dy1     h(x,y): blurring effects
                                                               d(x,y): source+background radiation

  When space-invariant PSF:
      g ( x, y )    h( x  x1 , y  y1 )·d ( x1 , y1 )dx1dy1  h( x, y )  d ( x, y )

                 p( x, y)  h( x, y)  d ( x, y)  n( x, y)

  In Fourier
  Space:          P(u, v)  H (u, v)·D(u, v)  N (u, v)
          Deconvolution ill-posed problem
Deconvolution: Image restoration based on some knowledge of the
degradation processes effects and/or the expected source distribution using
models of the statistical and blurring effects of the image formation system.

                              P(u, v) N (u, v)
Why not….. D(u, v)                   
                              H (u, v) H (u, v)
                                           1. at those points the division
Because:                                      operation is undefined or results in
   H(u,v) has zero o near zero                meaningless values !!!!
   values in most of (u,v) range           2. for points having very small
                                              |H(u,v)|, although the division can
                                              be done, the noise will be amplified
Conclusion:                                   to an intorelable extent !!!!

  Image deconvolution is an inverse solution ill-posed problem !!!
    All the algorithms which try to solve it should be concerned about obtaining a
    solution with the three following properties:     • existence
                                                      • uniqueness
                                                      • stability
            Deconvolution algorithms
1. Image formation, degradation and recovering model basic
               Classical approach:
                          •   Direct inverse filter in Fourier space
                          •   Linear Regularized Filters (Wiener)
                          •   Minimum root mean squares
                          •   CLEAN
               Bayessian approach:
                          •   Fixed Model: MLE, MAP
                          •   Variable Model: Pixons

2. Regularization hypotesis (positiveness, energy conservation,
   frequency bounders, …)

3. Convergence iterative algorithm (gradient methods, EM,
   successive substitutions,…)

4. Iteration Stopping protocol (a posteriori, feasibility tests,
   crossvalidation maximum,…)
The Bayesian Approach                                       (fixed models)
                                       P(a|p): “a posterior probability”, radiation a
                                       given the recorded data p
             ( p | a )  ( a )       P(a): “a prior probability” of the source of
( a | p )                            radiation
                  ( p )               P(p|a): likelihood
                                       P(p): Normalizing constant

Bayes Hypotesis: The most probable image a(x,y), given measured
data p(x,y) is obtained by maximizing P(a|p) with a suposed P(p|a)
image formation model (and maybe also a P(a) “a priory” knowledge).

   The most usual deconvolution algorithms division regarding
   how they maximize P(a|p) are:
      • MAP (Maximum A Posteriori), maximizes P(p|a)·P(a)
      • MLE (Maximum Likelihood Estimator), maximizes P(p|a) with P(a)=cte
      • MEM (Maximum Entropy Method) is MAP with p(a) given by maximum
        entropy hypotesis (=minimum data uncertainty)
                                       G+P_MLE deconvolution
                                     • sucessive sustitutions convergence method
    Algorithm                                                     B             B D f a          D            D
                                     • energy conservation:  qi ai                         p j   b j
                                                                                         ji i
    basics:                                                     i 1          i 1 j 1 Cj      j 1         j 1
                                     • Poisson photon gathering+Gauss readout
                                                   hj  hj 
                                                              k                                 ( p j k )2

                                                       
                                                                                                                         B f
                                                                                         1                       h j   ij ·ai  b j
                                        P k hj  e                   P( p j | k )            e     2 2

                                                         k!                             2                            i 1 C j

                                                                                                                        
                                                                                        D                  f ji p ' j   
                                                                               (s) 1
                                                             ( s 1)
                                                        ai                ai                                         
                                                                                   qi                 f jl al  C j b j 
                                                                                         j 1
                                                                                               l 1
                                                                                                                                    D

                ( k  p j )    2
                                      (h j ) k        Iterative and nonlinear algorithm                                           ( p        bj )
          ke
                     2    2
                                                                                                                                    j 1
                                                                                                             (s) 
                                        k!                                                                                                                    m
         k 0                                                                                                                                            
p' j           ( k  p j ) 2
                                                                 NOTE: This resulting algorithm,                          B      D            f ji p' j    
                                     (h j )   k
                                                                without flat and background, is
                                                                                                                         a i  B                          

         e         2 2
                                      k!
                                                                 similar to MLE Poisson
                                                                 (Richardson Lucy) computation
                                                                                                                        i 1      j 1
                                                                                                                                  l 1  f jl al  C j b j 
         k 0                                                  structure if p’=p and K=1                                                                 
2. Adaptative Deconvolution
       with Wavelets

Presentation of our multiresolution adaptative
deconvolution algorithm based on wavelets
            Wavelets descomposition
The wavelet transform replaces the Fourier transform’s sinusoidal waves
by a family generated by translations and dilations of a function called
Mother Wavelet.
                                                                  - b space and    a   scale
                                                                  -  a ,b ( x )
       WT ( f (a, b))  a       2
                                       f ( x)  ( x b ) dx

     Discrete wavelet transform (DWT) algorithm specific properties:
      DWT is computed with à trous algorithm.
      Mother wavelet function  is derived from a B3 cubic spline scaling function.
      This is a dyadic decomposition (scale  2i , i=1,...,n).
      No subsampling is applied  wavelet planes have same size as original image.

       General properties of wavelets descomposition:
        Spatial and frequential content is almost decoupled.
        Better noise vs.signal discrimination than Fourier transform.
        Good processing flexibility.
     Example of Wavelets descomposition
c0                          w1                                     w2

w3                          w4                                     c4

        Original                                                Residual plane
        image                       Wavelet planes                at scale 4
          co = w1            +      w2    +   w3     +   w4    +    c4

                   High frequency                             Low frequency
 Adaptative deconvolution with wavelets
To deconvolve the image, this is decomposed in wavelet planes at every
iteration of the method. The signal detection masks only applies
deconvolution to those pixels with significant signal.
                  Objectives of using wavelets: Multiescale deconvolution
                  Objectives of adaptative deconvolution:
                            - Do not amplify noise.
                            - To detect and deconvolute only those features with signal.
 Signal detection mask:
                      3                   

                                                                                                               w                             
                         v, j   v, j                                                                                      w vh,(js )
                                          r                                                                                 p'                   2

              1  exp  2                  
                                                                                
                                                                                                                           v, j
                                             if                    vr, j  0                   v, j     t
    mv , j
                           2( vr, j ) 2      
                                                              v, j
                                                                                      with
                                                                                                            noise std. desv. at
                           0                     if        v, j       r
                                                                           v, j    0               v, j 
                                                                                                               wavelet plane n
                                                                                                               sorrounding pixel j

                                                                                                                    
                                                                     Nw
                                                               f ji  w vh,(js )  mv , j w v , jj  w vh,(js ) 
                                                      D
                                            (s) 1
G+P_AWMLE : ai ( s 1)                 ai 
                                                qi
                                                                   v
                                                       j 1
                                                                             f jl al  C j b j
                                                                                                               
                                                                      l 1                                     
    3. Applications of our
 deconvolution algorithm with
    astrometric purposes:

Show astrometric objectives, aplications and
deconvolution artifacts and effects
  Deconvolution of Baker-Nunn (RAO) images
 Data overview
   -50cm - f/1 Baker-Nunn instrument!!!
   -very large FOV 4.2ºx4.2º,
   -moderately undersampled data:
   -significant object blending,
   -bright stars cause overblooming,
   -PSF is low correlated with seeing
   and dominated by systematic
   distortions across FOV, which can be
   -due to large pixel size, background
   level is dominated by Poisson noise.

                                            0.8degx0.8deg subframe: note object blending and
Deconvolution objective:                              blooming around bright stars

 Increase limit magnitude to detect faint
                                                     RAO: Rothney Astrophysical Observatory
 NEOs with minimum false detections                  University of Calgary (Canada)
    Protocol for Baker-Nunn images
       1) CCD image calibration: bias, dark and flatfield correction


  2) Compute astrometric plate transformation with reference catalogue
                                         wcstools.sua2, sextractor,iraf.ccmap

3) Perform aperture photometry of selected stars (bright but not saturated)
                                        iraf.phot, iraf.pstselect

                    4) PSF fitting with list created in 3)
                                        iraf.psf, iraf.seepsf

                        5) Image deconvolution

  6) sextractor object detection and matching with USNOA2.0 catalogue

             7) Compute limiting magnitude gain and rest of
        Deconvolution of Baker-Nunn images

  Benefits of image deconvolution

  1) Increase of SNR
Limiting magnitude gain

2) Increase of resolution
   Object deblending

                            Original image: in blue matched   Deconvoluted image 60 iterations:
                                       detections             yellow and red are new matched
        Deconvolution of Baker-Nunn images

Matched vs. non-matched detections

 Evolution of raw and matched                Adaptative        wavelet    deconvolution
 detections with the number of iterations.   introduces small number of non-matched
 Asymptotic converge translated into         detections.
 stable number of detections for large       In contrast to classical Richardson-Lucy
 number of iterations.                       algorithms, false detections increase
                                             linearly, not exponentially.
    Deconvolution of Baker-Nunn images

Where do the non-matched detections come from?

   Left: original image with explanation of the origin of non-matched detections displayed in red in the right
                        frame, which corresponds to a 60 iteration deconvoluted image.
           Deconvolution of Baker-Nunn images
Limiting magnitude gain
Image            # matched detections
Original                   1559
180 iterations

But deriving magnitude gain directly from
ratio of the number of detections is not
allways reliable, due to intrinsic magnitude
distribution of studied stellar field.

From estimating the completeness of
magnitude histograms with respect to
      USNOA2.0 distribution:

                                                 USNOA2.0 R magnitude distribution for the catalogue itself
                                               (top), deconvoluted image 180 iterations (middle) and original

       mR  0.6
                                                                      image (bottom)
             Deconvolution of QUEST images
                   QUEST= QUasar Ecuatorial Survey Team

Data overview
 QUEST images to deconvolve                 WIYN images to evaluate results
      V filter, frames from 2 CCDs             2 MiniMosaic WIYN
      1”/pixel,                                0.141”/pixel
      2 diferent fields in 2 single CCDs       Covers aprox. the same field
      Limiting magnitude (S/N>10)              6+38 QSO candidates detected by
       aprox. 19.2                               varability criteria in those fields

Deconvolution objective:
  Detect a subsample of QSO candidates likely to be lensed by:
     1. improving the resolution (deblend close companions)
     2. improving the limiting magnitude
                Deconvolution of QUEST images

Aplication of image deconvolution to 6 QSO candidates in the field. Each panel includes original QUEST
image, deconvolved QUEST (400 iterations) and high resolution WIYN image, repectively. Those in green
correspond to objecte present in all three images, in black those only resolved by WIYN and in red those
resolved in both deconvolved QUEST and WIYN.
              Algorithm properties
   Powerful technique for increasing the number of useful science
    objects from the faint part of magnitude distribution.

   Enhances spatial resolution which translates into better object
    deblending discrimination.

   Generalizable for whatever observing facility employing a CCD.

   Cheap: only requires fast computers. The algorithm is highly
    parallelizable to be run in distributed computers.

   Low speed convergence and heavy CPU and RAM usage.

   Systematic usage best suited for astrometric surveys with
    moderate data throughput.
   Selective usage is suited for all kind of images, either high or
    low resolution.
                  Work in progress
   Use multi frame reduction for removing cosmics and hot pixels.

   Improve PSF fitting to decrease false detections around bright stars.

   Compare internal astrometric accuracy of original and deconvoluted

   Accurate determination of CCD flats, gain and readout noise for better
    noise distribution modelling when deconvoluting.

   Define strategy to assess space variant PSF across the FOV

   Further evaluation of astrometric improvements and
    photometric effects.
                         To sum up:
 In this presentation we have presented our selected algorithm:
Gausian and Poison Adaptative Wavelets Maximum Likehood Estimator

 This algorithm has been proved with two diferent astrometrical type
of data and diferent scientific goals
 Those examples have shown how our algorithm can produce in the
particular case of astrometry of point sources:
                 • recover faint objects
                 • improve resolution

                      Final Conclusions
    Deconvolution can be a powerful technique in astrometry projects.
    But it’s crucial to select the proper algorithm carefully considering:
                       1. your data
                       2. your scientifical objectives

To top