MEM for HETGS Analysis

Notes / Publications:

• MSSL Workshop, Oct. 2002
  -- see "Extended Source Analysis for Grating Spectrometers" by D. Dewey.
• note_010319.txt and fig_010319.gif ( and .ps )

References:

• "Inverse Problems and the Use of A Priori Information",
  "Maximum Entropy Image Restoration",
  in Numerical Recipes in Fortran, 2nd Edition, as sections 18.4 and 18.7,
  W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, Cambridge, 1992.
• "Deconvolution for Real and Synthetic Apertures" by T.J. Cornwell,
  in Astronomical Data Analysis Software and Systems I,
  Ed.s D.M. Worrall, C. Biemesderfer, and J. Barnes, ASP Conf. Series Vol.25, p.163, 1992.
• "Maximum Entropy Image Restoration in Astronomy",
  R. Narayan and R. Nityananda, ARA&A Vol.24, p.127-170, 1986.
  PDF version
• "Image Restoration by a Powerful Maximum Entropy Method",
  S.F. Burch, S.F. Gull, and J. Skilling, Comp. Vision Graphics and Image Proc., Vol.23, pp.113-128, 1983.
• "A Simple Maximum Entropy Deconvolution Algorithm",
  T.J. Cornwell and K.F. Evans, Astron. Astrophys. Vol.143, pp.77-83, 1985.
• "Doppler imaging from artificial data",
  J.B. Rice and K.G. Strassmeier, A&A Suppl. Ser. Vol.147, pp.151-168, 2000.
  (Briefly mentions MEM and Tikhonov (2nd derriv. smoothness) criteria)

Software:

• IDL astro-lib's max_entropy.pro
  • The procedure me_test.pro calls the IDL max_entropy.pro routine to deconvolve some simple test data... The routine is for a spatially invariant PSF and doesn't operate through explicit minimization... It could be a nice way to sharpen up a zeroth-order image to be used as a default map...
  • Screen output and ouptut plots for me_test are shown here for the cases of noise-free data and (one realization of) data with Poisson noise added:

    Case	Screen	Plots
    Noise-free	Screen	Plots
    Poisson noise	Screen	Plots

  • The outputs show that with each iteration the Chi-squared agreement between data and convolved model decreases, while the three measures of entropy ( "- f ln(f)", "ln(I)", and "sqrt(I)" ) monotonically decrease from the initial maximum entropy solution (flat!). The simulated data shows i) an isolated narrow peak, ii) a plateu feature, and iii) a He-like "R-I-F" feature.
• PIXON Original IDL Code:
  • The original IDL pixon code is available from the pixon link below. Very briefly, the pixon method can be explained with the following diagram:

                             pixon map
                                 |
                                 V
      pseudo_image   -->/ pixon map smoothing /-->   deconvolved image  --->
                     ( pixon smoothing function )
    
          --->   deconvolved image   -->/ PSF application /-->   modelled data
                                           ( PSF array )
    

    The software carries out forward folding and adjusts the pseudo-image pixels to get the modelled data to have the lowest Chi-squared fit to the actual data.

    The "pixon map smoothing" provides a constraint on the set of images which can be produced from the pseudo-image pixels. The pixon map assigns to each pixel in the deconvolved image a smoothing size from a discrete set, going from 1 to some maximum size. Then a set of smoothed pseudo images is created - each smoothed by a simple function (e.g., an upside down parabola) scaled by the smoothing size. Then each pixel in the deconvolved image is assigned its value from the smoothed pseudo-image corresponding to its smoothing scale as given in the pixon map.

    After each Chi-square minimization, the pixon map is re-evaluated and another minimization carried out until the pixon map and resulting images converge.

    There are two main pieces to the algorithm:
    i) calculate the pixon map from the current data (and its noise level) and deconvolved image, and
    ii) for a given pixon map minimize the Chi-squared by adjusting the pseudo-image pixel values.

    The minimization is performed with a conjugate-gradient algorithm written in IDL:

    • minf_conj_grad.pro
    • which uses: minf_parabol_d.pro and minf_bracket.pro.

  • For convenience in running the code some local changes were made
    • local notes
  • The pixon code was tested on a simple 1D example. using the procedure pxn_test.pro.
  • It appears that the PSF array has to be the same size as the data(image) array with the non-zero PSF values centered in the array.
  • Because pixon saves arrays in common it is important to reset various arrays before starting the "real" iterations after doing the pxnmap=1.0 test solution.
  • Ouput plots at some iterations are given here:
    • Iteration "0" - "maximum likelyhood solution", no pixon constraints, i.e., the pixon map is identically 1.0 so that the mapping from pseudo-image to deconvolved image is the identity.
    • Iteration 1 - the first iteration starts with mostly large values of the pixon map (scale sizes) so the deconvolved image has to be mostly smooth...
    • Iteration 9 - that looks better:

    PIXON output (w/strange one pixel offset)

  • There appears to be a one pixel offset between the deconvolved image and the input data.
  • For comparison, the 1D simulated data (with a different noise realization) was processed with max_entropy.pro giving these iteration plots. The pitfalls of MEM show up here in the poor plateau reconstruction in the 16 to 23 range and the overly jagged reconstruction of the tall plateau at 25-30 (even more so at iteration 50):
    

    max_entropy output

Related and Random:

• Pixons
• Computer Vision Bibliography