PaperlayerView exact revision on arXiv

Source-linked AI summary

Imfit: A Fast, Flexible New Program for Astronomical Image Fitting

Peter Erwin

arXiv:1408.1097v2astro-ph.IMastro-ph.GA

TL;DR

IMFIT addresses the need for accurate, flexible galaxy image fitting when existing programs are difficult to extend and low-count Gaussian approximations can bias parameters. It combines an open object-oriented architecture with multiple components, minimizers, and statistical models. The paper shows that Poisson MLE statistics produce essentially unbiased fitted parameters where Gaussian χ2 fits can remain biased, even with Gaussian read noise.

  • Problem

    Existing image-fitting programs are often difficult to extend, while Gaussian approximations to Poisson statistics can bias fitted parameters in low-count images.

  • Method

    IMFIT uses an open object-oriented design with user-extensible image components, multiple minimization algorithms, and Gaussian or Poisson-based fitting statistics.

  • Results

    Poisson MLE statistics yield essentially unbiased parameter values, including when Gaussian read noise is present, whereas χ2 fits show directional parameter biases.

  • Takeaways & Limitations

    Poisson-based fitting is appropriate for low-count galaxy images where Gaussian χ2 approximations can bias derived parameters.

  • Takeaways & Limitations

    Levenberg-Marquardt can become trapped in local minima and requires nonnegative least-squares objectives, while bootstrap confidence intervals are time-consuming.

Abstract

from arXiv · show

I describe a new, open-source astronomical image-fitting program called Imfit, specialized for galaxies but potentially useful for other sources, which is fast, flexible, and highly extensible. A key characteristic of the program is an object-oriented design which allows new types of image components (2D surface-brightness functions) to be easily written and added to the program. Image functions provided with Imfit include the usual suspects for galaxy decompositions (Sersic, exponential, Gaussian), along with Core-Sersic and broken-exponential profiles, elliptical rings, and three components which perform line-of-sight integration through 3D luminosity-density models of disks and rings seen at arbitrary inclinations. Available minimization algorithms include Levenberg-Marquardt, Nelder-Mead simplex, and Differential Evolution, allowing trade-offs between speed and decreased sensitivity to local minima in the fit landscape. Minimization can be done using the standard chi^2 statistic (using either data or model values to estimate per-pixel Gaussian errors, or else user-supplied error images) or Poisson-based maximum-likelihood statistics; the latter approach is particularly appropriate for cases of Poisson data in the low-count regime. I show that fitting low-S/N galaxy images using chi^2 minimization and individual-pixel Gaussian uncertainties can lead to significant biases in fitted parameter values, which are avoided if a Poisson-based statistic is used; this is true even when Gaussian read noise is present.

1. INTRODUCTION

Galaxy structure is complex enough that direct 2D image fitting is often preferable to 1D profiles, while existing tools remain difficult to extend. IMFIT addresses these limitations with an open, extensible design, multiple fitting methods, and parallel computation.

  • Galaxy decompositions are needed to measure structures such as envelopes, nuclei, disks, scale lengths, half-light radii, and Sérsic indices.
  • 1D profiles discard most image information and can confuse non-axisymmetric structures with axisymmetric components.
  • Direct 2D fitting constructs full model images and compares them pixel-by-pixel with observations, while multi-Gaussian expansions complicate component interpretation.
  • Existing programs differ in supported components and minimizers, with GALFIT offering broad component flexibility and fast Levenberg-Marquardt optimization.
  • Most image-fitting programs are closed or difficult to extend with entirely new image functions.
  • Low-count images can produce biased fitted parameters under Gaussian approximations to Poisson statistics.
  • IMFIT provides an open-source, object-oriented framework for user-designed components, multiple fitting algorithms and statistics, and parallel execution across CPU cores.

2. GENERAL OUTLINE OF THE PROGRAM

IMFIT reads data, optional PSF, noise, and mask images, then builds a model object from a configuration file. A user-selected nonlinear minimizer repeatedly generates, convolves, and evaluates model images against the data.

  • IMFIT reads the data image, optional PSF, noise, and mask images, plus a configuration file specifying models, initial values, and parameter limits.
  • The ModelObject class stores data structures and image-function instances, while χ2 fitting constructs a noise image and combines weights with any mask.
  • Each nonlinear minimizer updates parameters, passes them to the model object, optionally applies PSF convolution, and receives a fit statistic from pixel-wise comparison.

3. CONSTRUCTING THE MODEL IMAGE

IMFIT builds model images from readable configuration files containing function blocks and parameterized image functions. Its extensible architecture supports component summation, PSF convolution, model-image generation, and component-flux analysis.

  • Configuration File: Configuration files define function blocks containing shared centers and one or more components such as bulges, disks, bars, rings, or nuclei.
  • Configuration File: Each image function specifies initial parameter values, optional bounds, or fixed parameters, and all parameters form the minimization vector.
  • Image Functions: An image function maps setup parameters and pixel coordinates to a computed intensity, while component intensities are summed into the model image.
  • Image Functions: IMFIT’s main program is independent of individual function implementations, allowing components ranging from constant backgrounds to line-of-sight integrations through 3D luminosity models.
  • PSF Convolution: Model images can be PSF-convolved using normalized FITS PSFs and FFT-based convolution, with image expansion used to reduce edge effects.
  • Makeimage: MAKEIMAGE reuses IMFIT’s model-construction code and configuration format to save complete or component-separated model images without fitting.
  • Makeimage: Complex simulations can be time-consuming because even simple fits typically require dozens to hundreds of function evaluations.

4. THE FITTING PROCESS

IMFIT fits model parameters by minimizing likelihood-based statistics that compare predicted and observed pixel values. It supports Gaussian and Poisson formulations, with Poisson maximum-likelihood statistics avoiding biases that can affect low-count Gaussian χ2 fits.

  • Likelihood Framework: The fitting process searches for parameters θ that make predicted pixel values m_i match observed values d_i under an assumed statistical model.
  • Likelihood Framework: Maximum-likelihood fitting converts the product of per-pixel probabilities into a negative log-likelihood sum suitable for minimization.
  • Statistical Models: The exact combination of Poisson photon statistics and Gaussian read noise produces an impractical infinite-series likelihood for rapid fitting.
  • Statistical Models: The Gaussian χ2 approximation replaces Poisson variance with σ = √m and is commonly considered appropriate at roughly 20 or more counts per pixel.
  • Statistical Models: IMFIT supports data-based and model-based Gaussian error estimates, user-supplied error images, and assumptions that final per-pixel errors remain Gaussian.
  • Poisson Statistics: For pure Poisson data, the Cash statistic and PMLR provide Poisson maximum-likelihood alternatives based on model and observed counts.
  • Poisson Statistics: PMLR has the same best-fitting parameters as C, but its nonnegative values enable fast least-squares minimization and support χ2-like goodness-of-fit interpretation at large N.
  • Results: Gaussian χ2 fitting can bias parameters even above the nominal count threshold, whereas Poisson MLE statistics avoid these biases in artificial and real galaxy images, including with Gaussian read noise.

4.2. Implementation: Specifying Per-Pixel Errors and Masking

IMFIT constructs per-pixel uncertainties for χ2 fitting from data or model intensities, user-supplied error maps, and masking information. Poisson-based statistics instead use mask-derived binary weights.

  • IMFIT estimates per-pixel Gaussian errors from data or model values when no separate error image is supplied.The estimates use the Gaussian approximation to Poisson statistics.
  • Data-based error estimates incorporate A/D gain, pre-subtracted sky background, and read noise in the input image units.
  • Mask values are converted to zi = 1 for valid pixels and zi = 0 for bad pixels, then combined with variances as wi = zi/σ2_i.
  • Users can provide FITS error maps containing Gaussian errors, variances, or precomputed pixel weights.
  • For Cash or PMLR minimization, the weight map is based directly on the binary mask, while model and data counts are represented by mi and di.

4.3. Minimization Algorithms

IMFIT offers three minimizers that trade speed against robustness to local minima and objective-function constraints. Levenberg–Marquardt is fastest, while Nelder–Mead and Differential Evolution provide broader applicability or greater resistance to local minima at higher computational cost.

  • Levenberg-Marquardt: Levenberg–Marquardt is the default fast minimizer and supports parameter bounds, with numerical differentiation simplifying the addition of new image functions.
  • Nelder-Mead simplex: Nelder–Mead requires an initial guess, handles objective functions that may be negative, is less prone to local minima than L-M, and is roughly an order of magnitude slower.
  • Differential Evolution: Differential Evolution begins with randomly generated parameter vectors, does not require initial guesses, and is least likely among IMFIT’s algorithms to become trapped in local minima.
  • Recommendations: For exploratory or large-scale fitting, the paper recommends L-M; N-M or DE are alternatives when local minima are a concern.
  • Comparison: 37: In a 256×256-pixel IC 3478 fit, N-M took ∼4 times and DE ∼60 times as long as L-M, while all three reached the same solution.
  • Outputs and statistics: IMFIT reports best-fit parameters, fit statistics, and optionally model and residual FITS images; reduced χ2 is available for χ2 and PMLR, but not directly for Cash.
  • Model comparison: AIC and BIC compare models fit to the same data, with a difference of at least ∼6 usually required to regard one model as clearly superior or inferior.

5. CONFIDENCE INTERVALS FOR FITTED

IMFIT estimates fitted-parameter uncertainty through L-M Hessian intervals or bootstrap resampling. Bootstrap distributions can provide asymmetric confidence intervals and agree with L-M estimates in the illustrated simple case, but require many repeated fits.

  • L-M confidence intervals: L-M can produce approximate 1-σ confidence intervals by inverting the Hessian matrix computed during minimization.
  • Limitations: The non-L-M minimizers do not compute confidence intervals, and rerunning L-M on their solutions does not work when minimizing Cash.
  • Bootstrap resampling: Bootstrap resampling repeatedly refits data images generated by sampling unmasked pixels with replacement, using the original best-fit parameters as starting values.
  • Bootstrap resampling: After n bootstrap iterations, the parameter distribution yields asymmetric 68% confidence intervals and a standard deviation; 68% corresponds to ±1-σ for near-Gaussian distributions.
  • Illustrated result: In the illustrated 500-round Sérsic bootstrap, parameter distributions were approximately Gaussian and had dispersions similar to the original L-M estimates.
  • Limitation: Bootstrap resampling is time-consuming because it effectively reruns the fit hundreds to thousands of times, although original best-fit starting values save some time.

6. IMAGE FUNCTIONS

IMFIT implements image functions as subclasses of a common abstract FunctionObject interface. This object-oriented structure makes adding new image functions relatively simple, while the documentation and command-line tools expose available functions and parameters.

  • Each image function subclasses the abstract FunctionObject base class and conforms to its interface.
  • Adding a new image function requires writing header and implementation files, registering the function, and recompiling IMFIT.
  • IMFIT documents function parameters and can print available function names and parameter lists with the --list-parameters flag.

6.1. 2D Components

IMFIT’s 2D components combine common analytic galaxy profiles with flexible geometry, specialized core and ring models, and broken-exponential structures. Their parameters control shape, scale lengths, transitions, and isophotal forms.

  • General geometry: IMFIT converts ellipticity to axis ratio and position angle to a rotated coordinate system before evaluating each component’s intensity profile.Pixel or subpixel coordinates are centered, rotated, and converted to a scaled radius used by the relevant one-dimensional profile.
  • Basic profiles: The program includes basic 2D background, Gaussian, Moffat, exponential, and generalized-ellipse components for modeling image structure and telescope PSFs.The Moffat shape parameter β controls wing strength, while β →∞ approaches a Gaussian.
  • Sérsic profiles: Sérsic components use effective radius re and index n to represent profiles ranging from Gaussian-like to exponential and de Vaucouleurs-like forms.The profile is equivalent to a Gaussian at n = 0.5, an exponential at n = 1, and a de Vaucouleurs profile at n = 4.
  • Specialized profiles: Core-Sérsic components join an inner power law to an outer Sérsic profile, with α controlling the sharpness of the transition at break radius rb.The model is designed for core galaxies and uses Ib to set the intensity at the break radius.
  • Specialized profiles: Broken-exponential components join inner and outer exponential zones at Rb, allowing different scale lengths and adjustable transition sharpness.Low α produces a gradual break, whereas high α produces an abrupt transition; the profile is intended for truncated or antitruncated disks.

6.2. 3D Components

IMFIT’s 3D components project disk and ring luminosity-density models into 2D images through numerical line-of-sight integration. They support arbitrary inclinations and flexible radial and vertical structures.

  • 3D projection framework: Three IMFIT components perform line-of-sight integration through 3D luminosity-density models to produce 2D surface-brightness projections.The models assume a symmetry plane inclined relative to the line of sight, with inclination ranging from face-on at i = 0° to edge-on at i = 90°.
  • 3D projection framework: For each image pixel, IMFIT transforms coordinates into a frame aligned with the line of nodes and integrates luminosity density along the line-of-sight coordinate s.The sky plane is defined at s = 0, and each integration point is mapped into the component’s native coordinates.
  • Disk models: ExponentialDisk3D uses numerical integration through an axisymmetric disk with an exponential radial profile and a sech2/n vertical profile.Its parameters include radial scale length h, vertical scale height z0, vertical-shape index n, and central luminosity density J0.
  • Disk models: ExponentialDisk3D produces the same model at inclinations of 75°, 85°, and 89°, with the 89° projection nearly matching the analytic EdgeOnDisk model.This comparison connects the numerical 3D projection to the perfectly edge-on analytic solution.
  • Disk models: BrokenExponential3D extends the disk model by replacing its radial density profile with the broken-exponential form while retaining the vertical structure.Its radial parameters are inherited from the two-dimensional BrokenExponential function.
  • Ring models: GaussianRing3D projects an elliptical ring with Gaussian radial density and exponential vertical structure at arbitrary inclination.The ring’s in-plane major axis is handled by rotating the native disk coordinates before line-of-sight integration.

7. PROGRAMMING NOTES

IMFIT is implemented in C++ with open-source numerical libraries and parallelized model-image computation. Its multithreading substantially reduces fitting time, including when PSF convolution is required.

  • Implementation: IMFIT is written in standard C++ and uses CFITSIO and FFTW, with NLopt and the GNU Scientific Library available as optional recommended libraries.The project also uses SCons for builds and CxxTest for unit testing.
  • Parallel computation: OpenMP distributes model-image computation across processor cores because image generation is usually the slowest part of fitting.This supports shared-memory systems, including multicore CPUs.
  • Parallel computation: FFTW accelerates Fourier transforms used for PSF convolution, while OpenMP parallelizes the broader model-image calculation.Both forms of parallelism can be combined during a fit.
  • Performance: Including a 35 × 35-pixel PSF, fitting a 500 × 500-pixel image drops from ∼280s without multithreading to ∼50s with FFT and OpenMP threading.Using multithreading only for the FFT reduces the time to ∼120s.
  • Performance: Users can reduce or disable multithreading through a command-line option when they do not want a fit to use all available CPU cores.

8. SAMPLE APPLICATIONS

The sample applications show Imfit modeling increasingly complex galaxy structures and demonstrate systematic low-count biases in Gaussian-χ2 fits that Poisson maximum-likelihood fitting avoids.

  • PGC 35772: Disk, Bar, and Ring: The two-component fit reduced residuals, with an elongated, misaligned Sérsic component modeling the galaxy’s bar.
  • PGC 35772: Disk, Bar, and Ring: Adding a GaussianRing to the exponential-plus-Sérsic model improved residuals inside and around the ring and produced the smallest AIC.
  • IC 5176: Thin and Thick Disks: The general 3D disk model took approximately 29 minutes versus 3m20s for the analytic 2D approach, while yielding ΔAIC ≈2305.
  • Simple Model Galaxy Images: The χ2 biases decreased with increasing counts but remained present, while Poisson MLE fits stayed unbiased even with moderate Gaussian read noise.
  • Simple Model Galaxy Images: In low-S/N images, data-based χ2 underestimated Sérsic n and re by 12.3% and 15.4%, respectively, while overestimating Ie by 34%.
  • Simple Model Galaxy Images: Model-based χ2 produced smaller biases in the opposite direction, whereas Poisson C fits differed from the input parameters by less than 0.1%.
  • Scope of Biases: For typical optical and near-infrared images, χ2 and Poisson MLE differences can often be ignored, but low-background u-band images can make χ2 bias significant.

10. SUMMARY

IMFIT is an open-source program for flexible 2D astronomical image fitting, with extensible image components and multiple fitting statistics. Its comparisons show that Poisson maximum-likelihood fitting avoids parameter biases found with χ2 minimization in low-count images.

  • IMFIT provides a fast, flexible, open-source framework for fitting galaxy or other astronomical images.Its design emphasizes extensibility through easily added 2D image functions.
  • Models combine multiple 2D image functions, including Gaussian, Moffat, exponential, Sérsic, broken-exponential, edge-on-disk, Core-Sérsic, and ring profiles.
  • Fits can minimize standard χ2 using estimated or user-supplied per-pixel variances, or use Poisson-based maximum-likelihood estimators for low-count images.The Poisson approach includes Cash and likelihood-ratio statistics.
  • Poisson MLE produces essentially unbiased Sérsic parameter values, even when Gaussian read noise is present.By contrast, χ2-based fits show directional biases in n, re, total luminosity, and Ie.
  • IMFIT recommends Poisson MLE instead of χ2 minimization when possible, especially below approximately 100 photoelectrons per pixel.

APPENDIX

The appendix compares Levenberg–Marquardt and bootstrap estimates of parameter uncertainties for image fits.

  • The appendix compares Levenberg–Marquardt estimates with bootstrap estimates.
  • The comparison concerns parameter uncertainties.
  • The setting is image fitting.

Parameter Estimates for Fits to PGC 35772

For PGC 35772, bootstrap and Levenberg–Marquardt uncertainty estimates are similar for a simple Sérsic model but diverge as model complexity increases.

  • For the Sérsic model, bootstrap uncertainties had mean and median ratios to L-M uncertainties of 1.38 and 0.81.
  • For the Sérsic + Exponential model, bootstrap-to-L-M uncertainty ratios had mean and median values of 0.48 and 0.51.
  • For the Sérsic + GaussianRing + Exponential model, bootstrap-to-L-M uncertainty ratios had mean and median values of 0.61 and 0.41.
  • Bootstrap uncertainties were typically about half the L-M uncertainties for the more complex models.

Parameter Estimates for Multiple Exposures of Elliptical Galaxies

Comparisons across repeated and differently exposed elliptical-galaxy images indicate that both L-M and bootstrap uncertainty estimates can underestimate parameter variability, with bootstrap estimates performing inconsistently across galaxies.

  • The exposure-comparison test expected approximately 68% of points below σ15s/|∆x| = 1 and 32% above it when uncertainties were accurate.
  • Bootstrap estimates exceeded unity for 50% of deviations, compared with 1/8 for L-M estimates.
  • Repeated-image comparisons found essentially all bootstrap and L-M estimates within ±1, indicating that both estimates tended to be too small, particularly for re.
  • Bootstrap estimates performed better for NGC 3379 but worse for NGC 3377.
  • Overall, both L-M and bootstrap estimates should probably be treated as underestimates of the true parameter uncertainties.
Loading 1408.1097v2…