## Abstract

The shapelets method for image analysis is based upon the decomposition of localized objects into a series of orthogonal components with convenient mathematical properties. We extend the ‘Cartesian shapelet’ formalism from earlier work, and construct ‘polar shapelet’ basis functions that separate an image into components with explicit rotational symmetries. These frequently provide a more compact parametrization, and can be interpreted in an intuitive way. Image manipulation in shapelet space is simplified by the concise expressions for linear coordinate transformations, and shape measures (including object photometry, astrometry and galaxy morphology estimators) take a naturally elegant form. Particular attention is paid to the analysis of astronomical survey images, and we test shapelet techniques upon real data from the *Hubble Space Telescope*. We present a practical method to automatically optimize the quality of an arbitrary shapelet decomposition in the presence of observational noise, pixelization and a point spread function. A central component of this procedure is the adaptive choice of the scale size and the truncation order of the shapelet expansion. A complete software package to perform shapelet image analysis is made available on the World Wide Web.

## Introduction

In the shapelets formalism (Refregier 2003; hereafter Shapelets I), individual objects in an image are decomposed into weighted sums of orthogonal basis functions. The particular set of basis functions has been chosen to be mathematically convenient for image manipulation and analysis. In astronomical images, it also provides a compact representation for the shapes of galaxies of all morphological types. Refregier & Bacon (2003; hereafter Shapelets II) showed how these properties could be used to measure the slight distortions in galaxy shapes due to weak gravitational lensing. The elegant behaviour of shapelets under a Fourier transform also enabled Chang & Refregier (2002) to reconstruct images from interferometric observations. Massey et al. (2004) used shapelets to simulate realistic astronomical images containing galaxies with complex morphologies. A classification scheme for galaxy morphologies using principal-component analysis of the shapelet basis set was discussed in that paper and applied to the Sloan Digital Sky Survey by Kelly & McKay (2004). A method similar to shapelets has also been independently suggested by Bernstein & Jarvis (2002; hereafter BJ02).

In this paper, we expand upon the earlier work of Shapelets I, Shapelets II and BJ02, developing the formalism of ‘polar shapelets’, in which an image is decomposed into components with explicit rotational symmetries. Whilst the original Cartesian shapelets remain useful in certain situations, the polar shapelets, which are separable in *r* and θ, frequently provide a more elegant and intuitive form. We find estimators of the flux, position and size of an object, that form naturally from groups of its polar shapelet coefficients. We calculate the behaviour of a polar shapelet model during linear coordinate transformations. We also improve the basic shapelet decomposition by incorporating treatments of pixelization, observational noise and point-spread functions, and optimizing the overall quality of image reconstruction while maximizing data compression. To test our method upon real data, we extract isolated galaxies from the *Hubble Deep Fields* (hereafter HDFs; Williams et al. 1996, 1998). These deep, high-resolution images from the *Hubble Space Telescope* (*HST*) provide typical examples of the irregular shapes of distant galaxies.

A complete idl software package to perform the image decomposition and shape analyses presented in this paper can be downloaded from .

This paper is laid out as follows. In Section 2, we introduce the Cartesian and polar shapelet basis functions, and their relation to each other. In Section 3, we investigate the qualitative effects of varying the shapelet scale size β and set quantitative goals for the optimization of this choice. In Section 4, we develop practical techniques to cope with the effects of pixelization, seeing and noise in real data. We then demonstrate various applications of shapelets: in Section 5, we illustrate the manipulation of images in terms of their changing polar shapelet coefficients under coordinate transforms. In Section 6, we construct basic shape estimators for a shapelet model, including flux, centroid and size measures. In Section 7, we develop more advanced shape measures that can be used to quantitatively distinguish galaxies of various morphological types. We conclude in Section 8.

## Shapelets Formalism

### Cartesian basis functions

The shapelet image decomposition method was introduced in Shapelets I, and a related method has been independently suggested by BJ02. The idea is to express the surface brightness of an object *f*(*x*, *y*) as a linear sum of orthogonal two-dimensional (2D) functions,

*H*(

_{n}*x*) is a Hermite polynomial of the order of

*n*, and β is the shapelet scale size. They are shown in Fig. 1. These Gauss-Hermite polynomials form a complete and orthonormal basis set; this ensures that the shapelet coefficients for any image can be simply and uniquely determined by evaluating the ‘overlap integral’

In practice, a shapelet expansion (equation 1) must be truncated at a finite order *n*_{1}+*n*_{2}≤*n*_{max}. The array of shapelet coefficients is sparse for typical galaxy morphologies, which therefore can be accurately modelled using only a few coefficients. As shown in Shapelets I, data compression ratios as high as 60:1 can be achieved for well resolved *HST* images. Note, however, that our choice of Gauss-Hermite basis functions was not governed by the physics of galaxy morphology and evolution but by the mathematics of image manipulation. As we shall see throughout this paper, a shapelet parametrization is mathematically convenient for many tasks common in astronomy and other sciences.

### Polar shapelet basis functions

Polar shapelets were introduced in Shapelets I as an orthogonal transformation of the Cartesian basis states, and were independently proposed by BJ02. They have all the useful properties of Cartesian shapelets, and a similar Gaussian weighting function with a given scale size β. However, the polar shapelet basis functions are instead separable in *r* and θ. This renders many operations more intuitive, and makes polar shapelet coefficients easy to interpret in terms of their explicit rotational symmetries.

The polar shapelet basis functions χ_{n,m}(*r*, θ; β) are also parametrized by two integers, *n* and *m*, and a smooth function *f*(*r*, θ) in polar coordinates may be decomposed into

*f*

_{n,m}are again given by the ‘overlap integral’

BJ02 showed that the ‘polar Hermite polynomials’ , which were described in Shapelets I, are related to associated Laguerre polynomials

for*n*

_{r}>

*n*

_{l}by where

*n*

_{l}and

*n*

_{r}are any non-negative integers. In this paper we shall instead prefer the simpler

*n*,

*m*notation, where

*n*=

*n*

_{r}+

*n*

_{l}and

*m*=

*n*

_{r}−

*n*

_{l}. In this scheme,

*n*can be any non-negative integer, and

*m*can be any integer between −

*n*and

*n*in steps of 2. We truncate the series at

*n*≤

*n*

_{max}. Although the only allowed states are those with

*n*and

*m*both even or both odd, this condition will not be written explicitly alongside every summation for the sake of brevity.

As plotted in Fig. 2, the dimensionful polar shapelet basis functions are therefore

These are different from the Laguerre expansion used by BJ02 in two ways. Those are the complex conjugate of our basis functions: i.e. their*m*is equivalent to our −

*m*. The Laguerre expansion in BJ02 is also normalized by one less factor of β. This dimensionality ensures that, as in the case of Cartesian shapelets, the polar shapelet basis functions are orthonormal and complete (see, e.g., Wünsche 1998) where δ is the Kronecker delta and the asterisk denotes complex conjugation. Only those basis functions with

*m*= 0 contain net flux.

Fig. 3 demonstrates the polar shapelet decomposition of a galaxy found in the HDF. The original image (middle left-hand panel) agrees well with the reconstruction using *n*_{max}= 20. The top panel shows the modulus of the polar shapelet coefficients as a function of the *n* and *m* indices. The dominant coefficients have small values of both indices, demonstrating the compactness of a polar shapelet representation, and further improved prospects for data compression. The bottom panel shows the reconstruction of the galaxy using only coefficients with given values of *n* or ∣*m*∣, thus highlighting the contributions of terms with specific rotational symmetries. The off-central bulge is captured in the ∣*m*∣= 1 coefficients and the main spiral arms in the ∣*m*∣= 2 coefficients. The spiral arms can also be seen as the rotation of the *n*-only reconstructions with increasing radius. The fainter spiral arms appear as an interplay of the ∣*m*∣= 4 and 5 coefficients.

### Conversion between Cartesian and polar shapelets

Cartesian shapelets are real functions, but polar shapelet basis functions χ_{n,m} and coefficients *f*_{n,m} are both complex. However, their symmetries

*f*(

**), such as the surface brightness of an image. Equations (9) and (12) imply that**

*x**f*(

**) is real if and only if Coefficients with**

*x**m*= 0 are thus wholly real. All polar shapelet coefficients are paired with their complex conjugate on the other side of the line

*m*= 0. Therefore, even though the polar shapelet coefficients

*f*

_{n,m}are generally complex, the number of independent parameters in the shapelet decomposition of a real function is conserved from the Cartesian case.

A set of Cartesian shapelet coefficients with *n*_{1}+*n*_{2}≤*n*_{max} can be transformed, into polar shapelet representation with *n*≤*n*_{max}, using

## Choice of Shapelet Scale Size

A shapelet decomposition requires values for the scale size β and for the centre of the basis functions *x*_{c} to be specified in advance. Choosing the centre is relatively easy: there are many methods well known in the astronomical literature to accurately determine astrometry from the flux-weighted moments of objects. However, the selection of β is a less well-posed problem. In this section, we shall first use some properties of polar shapelets to describe the effect that the choice of the scale size has upon a shapelet decomposition. We shall then set quantitative criteria for the selection of β in arbitrary galaxy images that we can implement in a practical algorithm.

Note that the selection of β will be linked to the selection of *n*_{max}. As shown in Shapelets I Section 2.4, these two parameters determine the maximum extent θ_{max} and finest resolution θ_{min} that can be successfully captured by a shapelet model. If *n*_{max}→∞, any object can be represented using any scale size β. However, if the shapelet expansion is truncated at finite *n*_{max}, the shape information needs to be more efficiently contained within fewer coefficients. It is clearly desirable in this situation to select a scale size β that compresses information, and lets us store the smallest possible number of coefficients.

### Radial profiles

Our discussion can be simplified by initially considering only the radial profile of an object, thus reducing the task to a one-dimensional problem. Let us consider an object with surface brightness *f*(** x**). The radial profile of the object is its brightness averaged in concentric rings about its centre, i.e.

*m*= 0 basis functions are invariant under rotations. These are given by The first few rotationally invariant basis functions are written explicitly in Table 1.

As a concrete example, we consider the decomposition of galaxy images from the Hubble Deep Fields (Williams et al. 1996, 1998). The mean radial profile of spiral galaxies is typically exponential, , with some characteristic radius scale *r*_{0}. Fig. 4 shows the shapelet reconstruction of an exponential radial profile using various values of β, with *n*_{max}= 20 and the integral in equation (5) calculated numerically.

As can be seen in the top panel of Fig. 4, the quality of the reconstruction depends on the choice of β. For small values (β ≲ 0.4*r*_{0}) the reconstruction is oscillatory and cuts off the profile at large radii (*r*≳ 1.5*r*_{0}). For large values (β≳ 0.8*r*_{0}), the reconstruction fails to reproduce the cusp at small radii (*r* ≲ 0.4*r*_{0}) and exceeds the true profile at *r*≃ 0.6*r*_{0}. However, for intermediate values (0.5*r*_{0} < β < 1.1*r*_{0}), the reconstruction is good throughout the range 0.1*r*_{0} ≲ *r* ≲ 2.8*r*_{0}. This range can of course be expanded by including more shapelet coefficients of higher order. As *n*_{max}→∞, the input model can be recovered with arbitrary precision using *any* value of β.

The corresponding behaviour in shapelet space is apparent in the bottom panel of Fig. 4. The *f*_{n,0} coefficients can be thought as the profile of the galaxy in shapelet space or the ‘shapelet profile’. For low values of β the shapelet profile is very flat, showing that the power is distributed almost evenly throughout all orders. For β= 0.5*r*_{0}, the coefficients *a*_{n,m} are seen numerically to be ∝ (*n* + 1)^{−2}. This will be an important result for the convergence of shape estimators formed from series of shapelet coefficients in Section 6. Convergence is fastest at β≈ 0.8*r*_{0}, with *a*_{n,m}∝ (*n* + 1)^{−5/2}. For higher values of β, the signs of *a*_{n,m} begin to alternate and the convergence once again falls below ∝ (*n* + 1)^{−2} at β≈ 1.1*r*_{0}.

Fig. 5 demonstrates the importance of a proper choice of the parameters β and *n*_{max} for the practical decomposition of a spiral galaxy in the HDF. Its spiral arms complicate measurement, but its radial profile is approximately an exponential with a scalelength of *r*_{0}≈ 3 arcsec (12 pixel). The left-hand column shows the growth in complexity of a shapelet model using increasing *n*_{max}. Note, in particular, the rotation of the core ellipticity as *n*_{max} is increased from 2 to 8 and higher-order moments are included to resolve the spiral arms. In this column, β is allowed to vary in order to minimize the least-squares difference between the model and the HDF image, shown at the bottom. The middle column shows shapelet decompositions at fixed *n*_{max}= 20, with varying β. The residuals are plotted in the right-hand column. As in Fig. 4, we find that the best overall image reconstruction uses 0.5*r*_{0} ≲β ≲ 0.7*r*_{0}. This is perhaps at the low end of the range suggested by Fig. 4 because of the extra high-frequency detail contained in the spiral arms.

By experimentation we have found a fairly wide range of β values that produce a faithful shapelet reconstruction. The information is then concentrated into the few lowest shapelet states, with fast convergence to the final model, and truncation is possible at a computationally acceptable value of *n*_{max}. We shall now consider ways to formalize this process, and hone our choice of *x*_{c}, β and *n*_{max} using quantitative criteria.

### Existing optimization methods

Methods in the literature that face similar choices suggest several distinct philosophies for the quantitative selection of parameters equivalent to *x*_{c}, β and *n*_{max}. The suggestions, outlined below, differ both in the goals set for an ideal decomposition and the method used to achieve it.

- (i)
Shapelets I suggested a geometrical argument using θ

_{min}, θ_{max}: the minimum [point spread function (PSF) or pixel] and maximum (entire image) sizes on which information exists. This could be iterated using functional rules on*x*_{c}and*r*^{2}_{f}as defined by shapelet coefficients. However, the coefficients change as a function of*n*_{max}, and it is not clear what the rules should be. - (ii)
Van der Marel & Franx (1993) fit one-dimensional (1D) Gauss-Hermite polynomials to spectral lines. They arbitrarily fix

*n*_{max}= 6, probably finding this sufficient because their spectra have relatively high signal-to-noise (S/N) ratios and their lines have a nearly Gaussian profile. Parameters equivalent to*x*_{c}and β are obtained from the best-fitting Gaussian. This also determines*f*_{0}and in 1D is equivalent to constraining*f*_{1}=*f*_{2}= 0, i.e. the derivatives of the Gaussian with respect to*x*_{c}and β. The number of variables is reduced and the problem rendered tractable. Unfortunately, this does not help us in 2D because while both*a*_{1±1}can be forced to zero by varying*x*_{c}, no unique recipe can be found for setting the three*n*= 2 states using only one value β. - (iii)
Van der Marel et al. (1994) relax the constraint on

*f*_{1}. This is an improvement as*f*_{1}is only the first term of an expression for the centroid, expanded using*all*odd*f*in equation (51). Without the higher-order corrections, the true object centroid is moved slightly from the origin: amongst other things rendering rotations and shear operations more complicated. Instead, they set_{n}*x*_{c}from the theoretical rest wavelength of a line. Unfortunately, astrometric calibration clearly cannot be performed with such accuracy. Nor has the*n*= 2 problem been solved. - (iv)
Kaiser, Squires & Broadhurst (1995) combine fitting with a stand-alone object detection algorithm, hfindpeaks. Translated into shapelet language, their approach is roughly equivalent to placing

*x*_{c}at data peaks and then finding a width β such that the signal-to-noise ratio ν in*f*_{00}is maximized. - (v)
Bernstein & Jarvis (2002) propose a similar approach. They prescribe β by requiring

*f*_{20}= 0, while moving*x*_{c}to ensure*f*_{1±1}= 0. Higher coefficients are then determined afterwards by linear decomposition. To first order, this β constraint is equivalent to that for hfindpeaks. This β is generally larger than values chosen by our χ^{2}method below, and it can be several times larger for a high signal-to-noise ratio object containing lots of substructure such as the galaxy in Fig. 3. This method may indeed prescribe the optimal decomposition for weak lensing as the shear signal in the quadrupole moments becomes concentrated in one number; however, a predisposition towards particular states often leads to poor overall image reconstruction and PSF deconvolution, so it is not necessarily ideal for all applications. - (vi)
Kelly & McKay (2004) were able to set a fixed

*physical*scale of β= 2 kpc for galaxies in the Sloan Digital Sky Survey, where photometric redshifts were available. However, galaxies have a broad distribution of physical sizes, and it may in fact become more difficult to interpret a shapelet model derived using this method. - (vii)
Marshall (in preparation) describes a fully Bayesian approach to applying the shapelet transform in the context of image reconstruction. Here,

*x*_{c}, β and*n*_{max}are varied in order to maximize the evidence (the probability of observing the data, marginalized over all shapelet coefficients). At high S/N ratio, this method gives a value of β which approaches the same as that from our χ^{2}method below, but otherwise tends to prefer a fractionally larger β, conservatively eliminating some ‘noise’ in favour of a smoother image reconstruction. However, this is computationally slow, a serious issue when analysing large images.

### Optimization of image reconstruction

We shall adopt a choice of β and *n*_{max} that is suitable for many applications, including overall image reconstruction and PSF deconvolution. Different models will be quantitatively compared via the overall reconstruction residual

*f*

_{obs}(

**) is the observed image and**

*x**f*

_{rec}(

**; β) is the reconstructed image from the shapelet model.**

*x**V*is the covariance matrix between pixel values, i.e. its diagonal elements are the noise variance in each pixel. We will need to know this a priori, or estimate it from blank areas of the image.

*n*

_{pixels}is the number of pixels in the observed image and

*n*

_{coeffs}is the number of shapelet coefficients used in the model. The residual itself has variance (Lupton 1993)

An example of typical χ^{2}_{r} isocontours on an *n*_{max} versus β plane is shown in Fig. 6 for an elliptical galaxy from the HDF. The horizontal trough is present for all galaxies (and many other isolated objects). This demonstrates that there is indeed an optimum β for the reconstruction of this image. As one might expect, it is roughly independent of *n*_{max}, but decreases very slightly as more coefficients are added. By increasing *n*_{max}→∞, the reconstruction can be improved to arbitrary precision. However, stopping at the χ^{2}_{r}= 1 contour produces a model where the residual is consistent with the noise. Additional coefficients would just model the background noise and should be excluded.

The form of these typical contours thus suggests a unique location in parameter space. We will choose β and *n*_{max} so that the model lies at the intersection of the trough and the χ^{2}_{r}= 1 contour, i.e. at the leftmost point on the contour. To achieve this, we set quantitative goals of

^{2}

_{r}= 1 or flattens out The first constraint ensures that the scale size is well suited to efficiently model the image. The second ensures that the shapelet centre matches the object centroid. The third guarantees that sufficient coefficients (

*n*

_{max}) are included to model an object, but with truncation that ‘smooths over’ observational noise. A flatness constraint is also included (in the right-hand side of equation 22). This is particularly important for galaxies with a near neighbour or for very faint objects that have noisy and fragmented χ

^{2}

_{r}contours. In these cases, including additional shapelet coefficients may not significantly improve a fit, so the series is truncated early.

We apply extra geometrical constraints to the minimum θ_{min} and maximum θ_{max} scales of the decomposition, to prevent the model from containing features smaller than the pixel scale or extending off the edge of an image, where it would be unconstrained.

### Automatic optimization algorithm

Satisfying the three conditions (20)–(22) would ensure that a shapelet decomposition uses the optimum values of *n*_{max}, β and *x*_{c}. It is easy to determine the values of these parameters once the entire *n*_{max} versus β plane has been examined, as in Fig. 6. However, this is a slow process, so we need a practical algorithm to more efficiently explore this parameter space, and to iterate rapidly towards our targets. The numerical implementation of this iteration will inevitably be non-trivial, because it combines both minimization and root finding, in a space with one axis discrete. Here we describe a code that we have developed to repeatedly decompose an object into shapelets, test the residual, and improve the decomposition parameters. Its stepwise approach is shown in Fig. 6, and the full code can be downloaded from the World Wide Web.

Objects are first detected in an image using sextractor (Bertin & Arnouts 1996), a friends-of-friends peak-finding algorithm. After experimenting on various data sets, we have found the results of sextractor highly sensitive to input settings. To avoid reliance upon these settings, we use sextractor as sparingly as possible. We set low detection thresholds in order to obtain a complete catalogue, and filter out false detections later. We use the measurement of the FWHM of each object to make an initial guess at β, and also to set the size of the fixed, circular ‘postage stamp’ region that is extracted around each object. We aim for a postage stamp large enough to contain the entire object, but small enough to isolate it from its neighbours and to make the routine computationally efficient. We then use the sextractor segmentation map to identify pixels in the postage stamp but well away from any object. These are used to estimate the background noise level, or to locally renormalize the pixel weight map. Within reasonable limits, the process is stable with respect to such parameters and we shall not be too concerned as to the exact sextractor settings.

Using constant *n*_{max}= 2 for speed, β is varied in order to minimize χ^{2}_{r} and satisfys the criterion in equation (20), via a 1D version of the *Numerical Recipes*amoeba routine: crawling vertically in Fig. 6. During each step of this iteration, the centroid is simultaneously shifted to re-zero the series in equation (51) in the shapelet coefficients and thus satisfy the criterion in equation (21). As the calculation of the centroid is independent of β for isolated objects (see Section 7), this part of the iteration is both stable and fast. Fig. 6 also shows the additional geometrical constraints of θ_{min} > 0.2 pixel and θ_{max} not falling off the edge of the postage stamp. These act as hard boundaries to the region of parameter space that the amoeba is allowed to explore.

Once the optimum β has been found, *n*_{max} is increased until the criterion in equation (22) is satisfied: crawling horizontally in Fig. 6. The increases are performed in steps of two, because even *n* states frequently improve the fit more than odd *n* states (primarily due to the addition of a new *f*_{n,0} circular state). The value of *n*_{max} is fine-tuned to the exact best value at the end. If two values of *n*_{max} both allow a decomposition with χ^{2}_{r}= 1 ± 1σ, the lower value is taken.

If the object warranted more coefficients than the initial guess of *n*_{max}= 2, β and *x*_{c} are again readjusted at the new *n*_{max}, using our 1D amoeba routine. Another *n*_{max} search then starts back at *n*_{max}= 2 and increases again in steps of two. The algorithm terminates when either the horizontal or vertical search returns to the same value as it started. All three conditions in equations (20)–(22) have then been met. Computation time for each object increases ∝*n*^{4}_{max}. On a single 2-GHz processor, our algorithm takes ∼ 45 min to process all of the 3596 objects detected in the HDF North.

A selection of reasonably bright HDF galaxies is shown with their shapelet models in Fig. 7. The right-hand column shows the reconstruction residuals, which are consistent with noise even for irregular galaxy morphologies. A comparison of their shapelet-based shape estimators to traditional sextractor measurements is shown in Fig. 8.

## Decomposition of Real Data

### Least-squares fitting

Unlike the continuous, analytic formalism presented in Section 2, real images are complicated by pixelization, PSF convolution and noise. In order to incorporate these effects, we shall first adopt a somewhat different approach to shapelet decomposition than the overlap integrals (3) and (5). We shall instead fit shapelet coefficients to the data using a least-squares method. As the model *f*_{rec}(** x**) in equation (18) is linear in the shapelet coefficients, we can solve for the minimum χ

^{2}

_{r}solution (18) exactly. We obtain (see Lupton 1993; Chang & Refregier 2002)

*f*_{n,m}is a vector of the derived shapelet coefficients,

*f*_{x,y}is the surface brightness in each pixel arranged as a data vector,

*V*is the covariance matrix between pixel values and

*M*is a matrix of each shapelet basis function evaluated in each pixel. A fit achieving χ

^{2}

_{r}= 1 has successfully modelled all significant spatial variation in the image and removed observational noise.

If the noise per pixel is known, 1σ confidence limits can be derived on all of the assigned coefficients using this fitting method (Lupton 1993). If a complete pixel noise map is available (e.g. from multiple exposures stacked using drizzle software —Fruchter & Hook 2002), it can be used to down-weight noisy pixels where cosmic rays or hot/cold pixels were present in some of the exposures. Although the code available on the World Wide Web simply uses a diagonal matrix for *V* that contains only the noise level in each pixel, the method is, in general, able to use the full covariance matrix that contains the amount of covariance between different pixels. In real data, the flux in adjacent pixels is indeed slightly correlated because of convolution with the PSF and also because of additional aliasing effects introduced by drizzle. If this effect is important, the pixel-to-pixel covariances could be estimated from empty regions of an image and included in the calculation. In particular, this may have a small improvement on statistics measured from very small objects (cf. Massey et al. 2004).

A constant background level can also be removed using this method, by adding an undetermined constant to the set of basis functions. Poor flat-fielding or local background gradients near a bright object can also be fitted and removed by adding a plane with a variable slope. Although these functions are not strictly orthogonal, the procedure works well in practice as long as there are sufficient pixels around the fitted object that contain only background noise.

### PSF deconvolution

All real images are inevitably seen after convolution with a point spread function. In astronomy, this is typically caused by atmospheric turbulence or ‘seeing’ (for ground-based observatories), aperture diffraction at the primary mirror and imperfect telescope tracking or optics. The combination of such effects can be measured from the size and shape of stars observed in an image (because these distant objects would be point-like in the absence of a PSF), and can be fitted with a shapelet model in the same way as the galaxies. Shapelets I presented the matrix operation for convolving an image with a Gaussian PSF in shapelet space. Shapelets II extended this derivation to a general PSF and demonstrated PSF *de*convolution via matrix inversion. However, the inversion of the PSF matrix is potentially slow and may be numerically unstable. Our least-squares fitting method will allow us to elegantly sidestep this process by convolving the basis functions with the PSF model in advance, then fitting this new basis set to the data. The returned shapelet model, reconstructed using the unconvolved basis functions, will be automatically deconvolved from the PSF.

The formalism for convolution in shapelet space is presented in Shapelets I Section 4 and involves three separate scale sizes for three separate objects: α for the unconvolved model, β for the PSF and γ for the convolved model (there are also corresponding values of *n*^{α}_{max}, *n*^{β}_{max} and *n*^{γ}_{max}). We assume that β is known. We can optimize α as in Section 3.3. However, the choice of γ is a matter entirely internal to the fitting procedure. Just as before, if *n*^{γ}_{max}→∞, any γ will work (but this time without increasing the number of external free parameters in the model). In practice, however, it is still necessary to truncate this series somewhere. Note that γ^{2}=α^{2}+β^{2} was incorrectly suggested as a ‘natural choice’ for this parameter in Shapelets I. Another choice would be γ=α, which, with *n*^{γ}_{max}=*n*^{α}_{max}, makes the convolution matrix *P*_{n,m} symmetric and thus simplifies its calculation.

The optimum values for γ and *n*^{γ}_{max} are, in fact, obtained from an argument concerning the information present in shapelet coefficients. A shapelet model contains information only between minimum and maximum scales

During convolution, θ^{α}_{min} and θ^{β}_{min} add in quadrature to produce θ^{γ}_{min}; similarly for θ^{γ}_{max}. γ and *n*^{γ}_{max} should therefore be chosen to precisely capture the information contained on all of these new scales. Writing (*n*^{α}_{max}+ 1) as *N*_{α}, etc. for brevity, we find

The PSFs of cameras on board the *Hubble Space Telescope* are well known and stable. Fig. 9 shows an oversampled tinytim (Krist & Hook 1997) model of the *Wide-Field Planetary Camera 2* (WFPC2) PSF, raytraced through an engineering model, plus charge diffusion to simulate photon capture within the CCD cameras. This is easy to model with shapelets, except for the fact that its steep cusp and extended wings are intrinsically ill-matched to the Gaussian around which shapelets are constructed. The PSF is shown in the figure beside a shapelet decomposition up to *n*^{PSF}_{max}= 15. This is sufficient to accurately capture the core and the first two diffraction rings, which are already more than two orders of magnitude below the maximum, but does not extend to the four faint diffraction spikes or far into the low-level wings (note that the colour scales are logarithmic). In principle, this could be further extended at a cost to processor time by using more shapelet coefficients.

Fig. 10 demonstrates successful PSF deconvolution. A galaxy from the HDF is convolved with the WFPC2 PSF (in real space). This is treated as the observed image, and deconvolved from the PSF using a shapelet fit. The resulting reconstruction is in good agreement with the original galaxy image, as can be seen from its small residual. Note that the optimum scale size β for the model is slightly lower when PSF deconvolution is performed. This reflects the need to capture finer details.

### Pixelization

Real image data are typically stored in discrete pixels. To link this to the analytic shapelet formalism, one must either smooth the data or pixellate the shapelet basis functions. Smoothing the data requires an arbitrary interpolation scheme to be defined, and resampling the data on to smaller pixels can be very slow. A better approach is to leave the data alone and to discretize the smooth shapelet basis functions. This reduces the integrals in equations (3) and (5) to sums over pixel values, which are fast to compute. However, they are no longer analytically exact. We therefore need to define a discretization scheme that keeps the basis functions as orthogonal as possible, and the integrals as accurate as possible.

As pointed out by Berry, Hobson & Withington (2004), one cannot simply adopt the value of basis functions at the centre of each pixel. Basis functions that contain oscillations on scales smaller than the pixel size are sampled in an essentially random manner. Their discrete versions are then neither representative of the analytic function nor orthogonal. Degeneracies are introduced between shapelet coefficients during the decomposition that destabilize the inversion of coefficient matrices in the reconstructed model, and bias quantities such as the flux of an object. Fortunately, this is rarely a problem in practical cases, because we can choose *n*_{max} and β in advance to isolate only those basis functions that contain oscillations on scales larger than the pixel (or seeing) size. Under these conditions, Berry et al. (2004) show that the shapelet basis functions are indeed orthogonal.

We suggest an even safer alternative here. The Cartesian basis functions are separable in *x* and *y*, and may be analytically integrated within rectangular pixels. This is exactly the same process undergone by photons arriving at a CCD, where the smooth function of a real scene becomes binned into digital squares. Once we have convolved the basis functions with the PSF, and integrated them within pixels, they can be suitably matched to the data.

To integrate the 2D Cartesian basis functions, first consider the 1D basis functions from Shapelets I,

Integrating by parts and using two well-known identities (see, e.g., Boas, chapter 12) and one can obtain the recurrence relation Finally, note that This supplies all the necessary integrals. As the 2D Cartesian basis functions are separable in*x*and

*y*, it is easy to extend this derivation to integrate within square CCD pixels: where, if there is no ‘dead zone’ around the edge of a pixel, (

*b*

_{1}−

*a*

_{1}) × (

*b*

_{2}−

*a*

_{2}) is the angular size of a pixel. A missing pixel border, due for instance to electronics which is unresponsive to light, can be included by altering the limits on the integral.

We can either use this result to obtain a model in Cartesian shapelet space, which can later be converted to a polar shapelet representation using equation (14), or we can integrate the polar shapelet basis functions within pixels using the same equations. This integration is a particularly important advance for small galaxies or for shapelet basis functions at high *n*, that can contain oscillations smaller than a single pixel.

The symmetries of polar shapelets can also be used to integrate models within circular apertures using equations (46)–(49).

## Coordinate Transformation

Image manipulation via linear transformations is simple in shapelet space. As in Shapelets I, let us consider an infinitesimal coordinate transformation ** x**→ (1 +

**Ψ**)

**+**

*x***ε**, where

**ε**={ε

_{1}, ε

_{2}} is a displacement and

**Ψ**is a 2 × 2 matrix parametrized as

_{i}correspond to infinitesimal rotations, dilations, translations and shears.

An image transforms as *f*(** x**) →

*f*′(

**) ≃**

*x**f*(

**−**

*x***Ψ**−

*x***ε**), which can be written as

^{−1}. These transformations can be viewed as a mapping of

*f*

_{n,m}coefficients in shapelet space. For example, a finite rotation is so a rotation through 180° can be written as

An (infinitesimal) dilation can be performed in polar shapelet space by mapping the shapelet coefficients as

The shapelet model may require more coefficients after this transformation. Note that this dilation operation increases both the flux and the image area by a factor of 1 + 2κ, thus conserving surface brightness. To instead perform a dilation that conserves the total*flux*, divide the right-hand side of equation (39) by this factor. To first order, this is In Section 6, we shall ensure that shape estimators for a shapelet model are independent of the scalefactor chosen for the decomposition by ensuring that the estimators are unchanged under this mapping.

Rather than these first-order approximations, finite dilations can be performed to all orders using the rescaling matrix in the appendix of Shapelets I. This is identical to the convolution matrix, but the image is convolved with a δ-function.

Shears and translations can be performed using

and with the translation specified in units of β.Other image manipulations can also be represented as mappings of shapelet coefficients. Changes of flux by a factor *B* are trivially implemented by the mapping

*x*-axis using Combining this with the rotation operator allows reflections to be performed in any axis.

The actions of these operators are demonstrated upon a real galaxy image in Fig. 11.

## Object Shape Measurement

The above symmetries of the polar shapelet basis functions can be used to identify combinations of shapelet coefficients that measure the flux (photometry), centroid position (astrometry) and the size of an object. Similar weighted combinations of Cartesian shapelet coefficients were found in Shapelets I, but we find the interpretation of polar shapelets more intuitive, and the expressions below are usually simpler than their Cartesian equivalents. For example, the rotationally invariant part of an object is isolated into its *m*= 0 coefficients. The linear offset of an object from the origin is described by its *m*=±1 coefficients and the ellipticity of an object by its *m*=±2 coefficients. In the latter cases, the magnitude of the coefficients indicate an amplitude and the phases a direction.

### Photometry

Practical measurements of the flux of an object usually introduce a Gaussian or top-hat weight function in order to limit contamination from surrounding noise and nearby objects. The flux of a shapelet model inside a circular aperture can be calculated using only the coefficients with *m*= 0. All other coefficients correspond to basis functions with positive and negative regions that cancel out under integration around θ. From equations (16) and (17) for the radial profile of an object, we find that

However, the imposition of an integration boundary is unnecessary with shapelets because the model is analytic and noise-free. In the limit of *R*→∞, we obtain a simple expression for the total flux in a shapelet model

This extrapolation to large radii does rely upon the faithful representation of an object by a shapelet expansion, and the removal of its noise via series truncation. Such truncation restricts the completeness of the basis functions, and a weight function (constructed from a combination of the allowed basis functions) akin to a ‘prior probability’ is subtly implicit inside our fitting procedure. However, a fitting method such as ours can beat a direct, pixel-by-pixel measurement. Our fit is able to include flux from the extended wings of an object, by integrating it over a large area, even when the signal lies beneath the noise level in any individual pixel. The wings of galaxies in Fig. 7 are indeed well captured by the shapelet models.

### Astrometry

It can similarly be shown that the unweighted centroid (*x*_{c}, *y*_{c}) is

*n*, because only these have the

*m*=±1 coefficients that possess the desired rotational symmetries.

### Size

Measures for the size and ellipticity of an object can be derived from the unweighted quadrupole moments of the object,

The rms radius*R*of an object is given by Integrals (46)–(49) can also be used to calculate Petrosian radii that enclose a specified fraction of the total flux within a circular aperture.

### Ellipticity

The unweighted ellipticity of an object can also be calculated from its quadrupole moments.

where the complex ellipticity notation of Blandford et al. (1991), with ɛ= ∣*e*∣ cos 2θ+ i ∣

*e*∣ sin 2θ, arises here naturally.

## Galaxy Morphology Classification

The shapelet decomposition of an object captures its entire structure, and useful information is frequently found in coefficients of higher order than those considered above. In particular, galaxy morphologies are well known to provide an indication of their physical properties, local environment and formation history. The classical ‘Hubble sequence’ of morphological types has been recently improved by several shape estimators that attempt to classify galaxies in a more quantitative manner, which correlates directly with the physical properties of interest (e.g. Simard 1998; Bershady, Jangren & Conselice 2000; van den Bergh 2002).

It is possible to manufacture such morphology diagnostics from weighted combinations of shapelet coefficients. Introducing shapelets to this field allows a measurement to take advantage of our robust treatment of noise, pixelization and PSF deconvolution. The shapelet expressions for existing shape measures are frequently elegant; and the natural symmetries in shapelets also suggest new diagnostics.

### General scale-invariant quantities

One approach is to consider general shape estimators *Q*, formed from a linear combination of shapelet coefficients, as an extension of the previous section

*w*

_{n,m}are arbitrary weights and the exponent

*s*sets the dimension of the estimator. These are also linear in the surface brightness of a galaxy. We initially restrict ourselves to using those combinations which are independent of β to at least first order. This ensures that the choice of the scalefactor does not affect the final result, and is also equivalent to invariance under object dilations (40). We can then impose further constraints that the estimator must be independent of or linearly dependent upon the various other operations described in Section 5. Setting ∂

*Q*/∂β= 0 and using the result that it is easy to show that we require Note that all quantities so formed mix coefficients with only one value of ∣

*m*∣. This can be chosen to give

*Q*the desired properties under rotation. Any term on the right-hand side should be ignored if it refers to non-existent states with negative

*n*. The normalization of the first term in each series,

*w*

_{m,m}, is arbitrary: this can be set to ensure independence to changes of object flux.

Setting (*s*, *m*) = (1, 0), (2, 1), (3, 0) and (3, 2) recovers the flux *F*, centroid *x*_{c}, rms square radius *R*^{2} and ellipticity ɛ, up to the normalization factor of *F*^{−1} for the latter three quantities. This proves that these are indeed the *only*β-invariant linear quantities with such dimensionality and rotational symmetries. Furthermore, as equations (50), (51), (54) and (55) describe for *unweighted* moments, they must in fact be independent of β to all orders.

All of these basic shape estimators converge for any galaxy with a shapelet spectrum steeper than *n*^{−2}. This includes both spiral galaxies with an exponential profile, and elliptical galaxies with a ‘de Vaucouleurs’ profile, as long as *n*_{max} is kept sufficiently low to prevent the high-*n* coefficients from modelling background noise at large radii. The flux and centroid estimators converge most rapidly, so are least sensitive to the choice of *n*_{max}. The error on these series due to truncation can be calculated using any of a range of methods for generic Taylor series in, for example, Boas (1983).

### Concentration

We can extend this sequence by raising *s* further. For example, setting *s*= 5 and *m*= 0 gives the 2D unweighted kurtosis of the image, producing an estimate of the concentration of the object. Unfortunately, such a high value of *s* yields a series of shapelet coefficients that does not converge for galaxies with a de Vaucouleurs or exponential radial profile.

We have also noted that a ratio of the two existing shapelet scale sizes, *R* and β, also works rather well as a concentration index (although this estimator is not independent of β). Further work will need to be performed to calibrate this estimator to the physical properties of galaxies.

An alternative approach is to mimic pre-existing, and pre-calibrated, morphology diagnostics. Bershady et al. (2000) define a concentration index

where*r*

_{80}and

*r*

_{20}are the radii of circular apertures containing 80 and 20 per cent of the total flux of the object. This correlates well with a Hubble-type galaxy (Bershady et al. 2000) and also its mass (Conselice, Gallagher & Wyse 2002). Integrals (46)–(49) can be evaluated for various values of

*R*, to find

*r*

_{80}and

*r*

_{20}, and thus calculate this quantity for a shapelet model.

### Asymmetry

Conselice, Bershady & Jangren (2000) define an asymmetry index

where the superscript denotes an image rotated through 180°. A term dealing with the background noise and sky level has been omitted here, as these are automatically dealt with during the shapelet decomposition process in Section 4. The asymmetry correlates with the star formation rate of a galaxy (Conselice et al. 2000), and high asymmetry values often indicate recent galaxy interactions or mergers.In a shapelet expansion, all of the information concerning galaxy asymmetry is contained in coefficients with odd *m* (and *n*). Using the orthonormality condition (9) and rotating via equation (38), we find the simple form

Estimators of asymmetry under rotations of 120° or 90° can also be formed from sums of shapelet coefficients with *m* not divisible by 3 or 4, respectively.

### Chirality

A quadratic combination of shapelet coefficients can be used to describe the ‘chirality’ or ‘handedness’ of an object. One dimensionless estimator χ_{∣m∣} can be formed for each value of ∣*m*∣, to trace the relative rotation of those coefficients, with increasing *n*. This is roughly equivalent to tracing the rotation of the isophotes of a galaxy with increasing radius. For example, the galaxy shown in Fig. 3 has two prominent spiral arms that unwind in a clockwise sense, so it has a large, positive value of χ_{2}.

We require that the chirality estimators should be invariant under global rotation of the object, invariant under changes of flux, invariant to first order under changes of β and to flip sign when the object is mirror-imaged. These conditions uniquely specify

where*w*

_{m,m +2}= 1 and thus mixing all coefficients with the same value of

*m*.

This estimator has yet to be calibrated against physical quantities. However, this approach ought to be able to automatically distinguish between elliptical galaxies and spiral galaxies in way that mimics a visual classification, and could also be adapted as a function of *n*_{max} to find bars in the cores of spiral galaxies.

## Conclusions

We have extended the formalism of shapelets for image analysis from basis functions separable in Cartesian to polar coordinates. Cartesian shapelets are convenient for the initial object decomposition. In particular, we have shown that they can analytically be integrated inside a square boundary, thus facilitating the pixelization of the smooth basis functions. On the other hand, polar shapelets decompose an object into components with explicit rotational symmetries, and often have a more direct physical interpretation. In addition, they yield more compact representations of typical galaxy images, as terms with low orders of rotational symmetry tend to dominate.

We have quantitatively investigated the effects of the choice of the shapelet scale size parameter, β. For most objects in astronomical images, one scale size is clearly optimal for high-quality image reconstruction, data compression and the fast convergence of shape estimators. We have developed a practical algorithm to find this value of β for arbitrary objects in real images, plus optimum values for the shapelet centre *x*_{c} and truncation order *n*_{max}. This algorithm can also take into account observational effects including noise, pixelization and PSF deconvolution.

We have then described a number of applications of polar shapelets. Shapelet models can be rotated, enlarged and sheared by simple analytic operations. As the shapelet basis functions are invariant under Fourier transform, analytic convolutions and deconvolutions (e.g. from a PSF) are also easy to perform. Linear combinations of the polar shapelet coefficients of an object produce elegant expressions for its flux (photometry), position (astrometry) and size. We also showed how other combinations of shapelet coefficients can be used to distinguish between morphological types of galaxies. A complete idl software package to perform the image decomposition and shapelet image analysis is publicly available on the World Wide Web at .

### Acknowledgments

The authors thank David Bacon, Gary Bernstein, Sarah Bridle, Tzu-Ching Chang, Mark Coffey, Chris Conselice, Phil Marshall and Jean-Luc Starck for invaluable insights and constructive conversations. Astute suggestions from the referee helped shape this paper into a more logical progression.

## References

### Appendix

#### Appendix: Laguerre polynomials

Different conventions have been used to define the Laguerre polynomials, especially before the 1960s. The *p*! term is omitted from equation (6) in many older books, and caution must be observed with the resulting relations. Several useful recursion relations can be derived to simplify their calculation (e.g. Boas 1983, chapter 12), which we gather here, using our convention, for future reference: