The HSM Module

Routines for adaptive moment estimation and PSF correction.

This module contains code for estimation of second moments of images, and for carrying out PSF correction using a variety of algorithms. The algorithms are described in Hirata & Seljak (2003), and were tested/characterized using real data in Mandelbaum et al. (2005). Note that these routines for moment measurement and shear estimation are not accessible via config, only via python. There are a number of default settings for the code (often governing the tradeoff between accuracy and speed) that can be adjusting using an optional hsmparams argument as described below.

The moments that are estimated are “adaptive moments” (see the first paper cited above for details); that is, they use an elliptical Gaussian weight that is matched to the image of the object being measured. The observed moments can be represented as a Gaussian sigma and a Shear object representing the shape.

The PSF correction includes several algorithms, three that are re-implementations of methods originated by others and one that was originated by Hirata & Seljak:

One from Kaiser, Squires, & Broadhurst (1995), “KSB”
One from Bernstein & Jarvis (2002), “BJ”
One that represents a modification by Hirata & Seljak (2003) of methods in Bernstein & Jarvis (2002), “LINEAR”
One method from Hirata & Seljak (2003), “REGAUSS” (re-Gaussianization)

These methods return shear (or shape) estimators, which may not in fact satisfy conditions like \(|e|<=1\), and so they are represented simply as e1/e2 or g1/g2 (depending on the method) rather than using a Shear object, which IS required to satisfy \(|e|<=1\).

These methods are all based on correction of moments, but with different sets of assumptions. For more detailed discussion on all of these algorithms, see the relevant papers above.

Users can find a listing of the parameters that can be adjusted using the hsmparams keyword, along with default values, under galsim.hsm.HSMParams below.

Shape Measurement Functions

galsim.hsm.FindAdaptiveMom(object_image, weight=None, badpix=None, guess_sig=5.0, precision=1e-06, guess_centroid=None, strict=True, check=True, round_moments=False, hsmparams=None, use_sky_coords=False)[source]

Measure adaptive moments of an object.

This method estimates the best-fit elliptical Gaussian to the object (see Hirata & Seljak 2003 for more discussion of adaptive moments). This elliptical Gaussian is computed iteratively by initially guessing a circular Gaussian that is used as a weight function, computing the weighted moments, recomputing the moments using the result of the previous step as the weight function, and so on until the moments that are measured are the same as those used for the weight function. FindAdaptiveMom can be used either as a free function, or as a method of the Image class.

By default, this routine computes moments in pixel coordinates, which generally use (x,y) for the coordinate variables, so the underlying second moments are Ixx, Iyy, and Ixy. If the WCS is (at least approximately) just a PixelScale, then this scale can be applied to convert the moments’ units from pixels to arcsec. The derived shapes are unaffected by the pixel scale.

However, there is also an option to apply a non-trivial WCS, which may potentially rotate and/or shear the (x,y) moments to the local sky coordinates, which generally use (u,v) for the coordinate variables. These coordinates are measured in arcsec and are oriented such that +v is towards North and +u is towards West. In this case, the returned values are all in arcsec, and are based instead on Iuu, Ivv, and Iuv. To enable this feature, use use_sky_coords=True. See also the method ShapeData.applyWCS for more details.

Note

The application of the WCS implicitly assumes that the WCS is locally uniform across the size of the object being measured. This is normally a very good approximation for most applications of interest.

Like EstimateShear, FindAdaptiveMom works on Image inputs, and fails if the object is small compared to the pixel scale. For more details, see EstimateShear.

Example:

>>> my_gaussian = galsim.Gaussian(flux=1.0, sigma=1.0)
>>> my_gaussian_image = my_gaussian.drawImage(scale=0.2, method='no_pixel')
>>> my_moments = galsim.hsm.FindAdaptiveMom(my_gaussian_image)

or:

>>> my_moments = my_gaussian_image.FindAdaptiveMom()

Assuming a successful measurement, the most relevant pieces of information are my_moments.moments_sigma, which is |det(M)|^(1/4) (= sigma for a circular Gaussian) and my_moments.observed_shape, which is a Shear. In this case, my_moments.moments_sigma is precisely 5.0 (in units of pixels), and my_moments.observed_shape is consistent with zero.

Methods of the Shear class can be used to get the distortion e, the shear g, the conformal shear eta, and so on.

As an example of how to use the optional hsmparams argument, consider cases where the input images have unusual properties, such as being very large. This could occur when measuring the properties of a very over-sampled image such as that generated using:

>>> my_gaussian = galsim.Gaussian(sigma=5.0)
>>> my_gaussian_image = my_gaussian.drawImage(scale=0.01, method='no_pixel')

If the user attempts to measure the moments of this very large image using the standard syntax,

>>> my_moments = my_gaussian_image.FindAdaptiveMom()

then the result will be a GalSimHSMError due to moment measurement failing because the object is so large. While the list of all possible settings that can be changed is accessible in the docstring of the HSMParams class, in this case we need to modify max_amoment which is the maximum value of the moments in units of pixel^2. The following measurement, using the default values for every parameter except for max_amoment, will be successful:

>>> new_params = galsim.hsm.HSMParams(max_amoment=5.0e5)
>>> my_moments = my_gaussian_image.FindAdaptiveMom(hsmparams=new_params)

Parameters:

object_image – The Image for the object being measured.
weight – The optional weight image for the object being measured. Can be an int or a float array. Currently, GalSim does not account for the variation in non-zero weights, i.e., a weight map is converted to an image with 0 and 1 for pixels that are not and are used. Full use of spatial variation in non-zero weights will be included in a future version of the code. [default: None]
badpix – The optional bad pixel mask for the image being used. Zero should be used for pixels that are good, and any nonzero value indicates a bad pixel. [default: None]
guess_sig – Optional argument with an initial guess for the Gaussian sigma of the object (in pixels). [default: 5.0]
precision – The convergence criterion for the moments. [default: 1e-6]
guess_centroid – An initial guess for the object centroid (useful in case it is not located at the center, which is used if this keyword is not set). The convention for centroids is such that the center of the lower-left pixel is (image.xmin, image.ymin). [default: object_image.true_center]
strict – Whether to require success. If strict=True, then there will be a GalSimHSMError exception if shear estimation fails. If set to False, then information about failures will be silently stored in the output ShapeData object. [default: True]
check – Check if the object_image, weight and badpix are in the correct format and valid. [default: True]
round_moments – Use a circular weight function instead of elliptical. [default: False]
hsmparams – The hsmparams keyword can be used to change the settings used by FindAdaptiveMom when estimating moments; see HSMParams documentation for more information. [default: None]
use_sky_coords – Whether to convert the measured moments to sky_coordinates. Setting this to true is equivalent to running applyWCS(object_image.wcs, image_pos=object_image.true_center) on the result. [default: False]

Returns:

a ShapeData object containing the results of moment measurement.

galsim.hsm.EstimateShear(gal_image, PSF_image, weight=None, badpix=None, sky_var=0.0, shear_est='REGAUSS', recompute_flux='FIT', guess_sig_gal=5.0, guess_sig_PSF=3.0, precision=1e-06, guess_centroid=None, strict=True, check=True, hsmparams=None)[source]

Carry out moments-based PSF correction routines.

Carry out PSF correction using one of the methods of the HSM package (see references in docstring for file hsm.py) to estimate galaxy shears, correcting for the convolution by the PSF.

This method works from Image inputs rather than GSObject inputs, which provides additional flexibility (e.g., it is possible to work from an Image that was read from file and corresponds to no particular GSObject), but this also means that users who wish to apply it to compound GSObject classes (e.g., Convolve) must take the additional step of drawing their GSObject into Image instances.

This routine assumes that (at least locally) the WCS can be approximated as a PixelScale, with no distortion or non-trivial remapping. Any non-trivial WCS gets completely ignored.

Note

There is not currently an option to use_sky_coords as we have for FindAdaptiveMom. We would welcome a PR adding this feature if someone wants it for their science case.

Note that the method will fail if the PSF or galaxy are too point-like to easily fit an elliptical Gaussian; when running on batches of many galaxies, it may be preferable to set strict=False and catch errors explicitly, as in the second example below.

This function has a number of keyword parameters, many of which a typical user will not need to change from the default.

Example:

Typical application to a single object:

>>> galaxy = galsim.Gaussian(flux=1.0, sigma=1.0)
>>> galaxy = galaxy.shear(g1=0.05, g2=0.0)  # shears the Gaussian by (0.05, 0) using the
>>>                                         # |g| = (a-b)/(a+b) definition
>>> psf = galsim.Kolmogorov(flux=1.0, fwhm=0.7)
>>> final = galsim.Convolve(galaxy, psf)
>>> final_image = final.drawImage(scale=0.2)
>>> final_epsf_image = psf.drawImage(scale=0.2)
>>> result = galsim.hsm.EstimateShear(final_image, final_epsf_image)

After running the above code, result.observed_shape is a Shear object with a value of (0.0438925349133, -2.85747392701e-18) and result.corrected_e1, result_corrected_e2 are (0.09934103488922119, -3.746108423463568e-10), compared with the expected (0.09975, 0) for a perfect PSF correction method.

The code below gives an example of how one could run this routine on a large batch of galaxies, explicitly catching and tracking any failures:

>>> n_image = 100
>>> n_fail = 0
>>> for i=0, range(n_image):
>>>     #...some code defining this_image, this_final_epsf_image...
>>>     result = galsim.hsm.EstimateShear(this_image, this_final_epsf_image, strict=False)
>>>     if result.error_message != "":
>>>         n_fail += 1
>>> print "Number of failures: ", n_fail

Parameters:

gal_image – The Image of the galaxy being measured.
PSF_image – The Image for the PSF.
weight – The optional weight image for the galaxy being measured. Can be an int or a float array. Currently, GalSim does not account for the variation in non-zero weights, i.e., a weight map is converted to an image with 0 and 1 for pixels that are not and are used. Full use of spatial variation in non-zero weights will be included in a future version of the code.
badpix – The optional bad pixel mask for the image being used. Zero should be used for pixels that are good, and any nonzero value indicates a bad pixel.
sky_var – The variance of the sky level, used for estimating uncertainty on the measured shape. [default: 0.]
shear_est – A string indicating the desired method of PSF correction: ‘REGAUSS’, ‘LINEAR’, ‘BJ’, or ‘KSB’. The first three options return an e-type distortion, whereas the last option returns a g-type shear. [default: ‘REGAUSS’]
recompute_flux – A string indicating whether to recompute the object flux, which should be ‘NONE’ (for no recomputation), ‘SUM’ (for recomputation via an unweighted sum over unmasked pixels), or ‘FIT’ (for recomputation using the Gaussian + quartic fit). [default: ‘FIT’]
guess_sig_gal – Optional argument with an initial guess for the Gaussian sigma of the galaxy (in pixels). [default: 5.]
guess_sig_PSF – Optional argument with an initial guess for the Gaussian sigma of the PSF (in pixels). [default: 3.]
precision – The convergence criterion for the moments. [default: 1e-6]
guess_centroid – An initial guess for the object centroid (useful in case it is not located at the center, which is used if this keyword is not set). The convention for centroids is such that the center of the lower-left pixel is (image.xmin, image.ymin). [default: gal_image.true_center]
strict – Whether to require success. If strict=True, then there will be a GalSimHSMError exception if shear estimation fails. If set to False, then information about failures will be silently stored in the output ShapeData object. [default: True]
check – Check if the object_image, weight and badpix are in the correct format and valid. [default: True]
hsmparams – The hsmparams keyword can be used to change the settings used by EstimateShear when estimating shears; see HSMParams documentation for more information. [default: None]

Returns:

a ShapeData object containing the results of shape measurement.

HSM output

class galsim.hsm.ShapeData(image_bounds=galsim.BoundsI(), moments_status=-1, observed_shape=galsim.Shear(0j), moments_sigma=-1.0, moments_amp=-1.0, moments_centroid=galsim.PositionD(x=0.0, y=0.0), moments_rho4=-1.0, moments_n_iter=0, correction_status=-1, corrected_e1=-10.0, corrected_e2=-10.0, corrected_g1=-10.0, corrected_g2=-10.0, meas_type='None', corrected_shape_err=-1.0, correction_method='None', resolution_factor=-1.0, psf_sigma=-1.0, psf_shape=galsim.Shear(0j), error_message='')[source]

A class to contain the outputs of the HSM shape and moments measurement routines.

The ShapeData class contains the following information about moment measurement (from either EstimateShear or FindAdaptiveMom:

image_bounds: a BoundsI object describing the image.
moments_status: the status flag resulting from moments measurement; -1 indicates no attempt to measure, 0 indicates success.
observed_shape: a Shear object representing the observed shape based on adaptive moments.
observed_e1, observed_e2: floats representing the e1 and e2 components respectively of the observed_shape.
moments_sigma: size sigma=(det M)^(1/4) from the adaptive moments, in units of pixels; -1 if not measured. (If FindAdaptiveMom is called with use_sky_coords=True, then the units will be arcsec.)
moments_amp: total image intensity for best-fit elliptical Gaussian from adaptive moments. Normally, this field is simply equal to the image flux (for objects that follow a Gaussian light distribution, otherwise it is something approximating the flux). However, if the image was drawn using drawImage(method='sb') then moments_amp relates to the flux via flux=(moments_amp)*(pixel scale)^2.
moments_centroid: a PositionD object representing the centroid based on adaptive moments, in units of pixels. The indexing convention is defined with respect to the BoundsI object defining the bounds of the input Image, i.e., the center of the lower left pixel is (image.xmin, image.ymin). An object drawn at the center of the image should generally have moments_centroid equal to image.true_center. (If FindAdaptiveMom is called with use_sky_coords=True, then the units will be arcsec, measured in sky coordinates with respect to the image center.)
moments_rho4: the weighted radial fourth moment of the image (dimensionless).
moments_n_iter: number of iterations needed to get adaptive moments, or 0 if not measured.

If EstimateShear was used, then the following fields related to PSF-corrected shape will also be populated:

correction_status: the status flag resulting from PSF correction; -1 indicates no attempt to measure, 0 indicates success.
corrected_e1, corrected_e2, corrected_g1, corrected_g2: floats representing the estimated shear after removing the effects of the PSF. Either e1/e2 or g1/g2 will differ from the default values of -10, with the choice of shape to use determined by the correction method (since the correction method determines what quantity is estimated, either the shear or the distortion). After a measurement is made, the type of shape measurement is stored in the ShapeData structure as ‘meas_type’ (a string that equals either ‘e’ or ‘g’).
corrected_shape_err: shape measurement uncertainty sigma_gamma per component. The estimate of the uncertainty will only be non-zero if an estimate of the sky variance was passed to EstimateShear.
correction_method: a string indicating the method of PSF correction (will be “None” if PSF-correction was not carried out).
resolution_factor: Resolution factor R_2; 0 indicates object is consistent with a PSF, 1 indicates perfect resolution.
psf_sigma: size sigma=(det M)^(1/4) of the PSF from the adaptive moments, in units of pixels; -1 if not measured.
psf_shape: a Shear object representing the observed PSF shape based on adaptive moments.
error_message: a string containing any error messages from the attempt to carry out PSF-correction.

The ShapeData object can be initialized completely empty, or can be returned from the routines that measure object moments (FindAdaptiveMom) and carry out PSF correction (EstimateShear).

applyWCS(wcs, image_pos)[source]

Convert moments in pixel coordinates to moments in sky coordinates.

Natively, the HSM algorithm computes second moments in (x,y) coordinates in the sensor coordinate system. However, many applications of second moments for shape measurements need the measurements in sky coordinates. This method converts the moments to (u,v) coordinates, where +v is towards North and +u is towards West.

The values that are different from the original are:

moments_sigma is changed to be in units of arcsec.
observed_shape is changed to be in (u,v) coordinates.
moments_centroid is changed to be in (u,v) coordinates relative to image_pos.

Note

This currently only works for the measurements from FindAdaptiveMom. If the input ShapeData instance has any values set from EstimateShear, they will not be present in the return value. We would welcome a PR adding the ability to work on corrected shapes if someone wants it for their science case.

Parameters:

wcs – The WCS to apply.
image_pos – The position in image coordinates (x,y) of the object whose moments have been measured.

HSM parameters

class galsim.hsm.HSMParams(nsig_rg=3.0, nsig_rg2=3.6, max_moment_nsig2=0, regauss_too_small=1, adapt_order=2, convergence_threshold=1e-06, max_mom2_iter=400, num_iter_default=-1, bound_correct_wt=0.25, max_amoment=8000.0, max_ashift=15.0, ksb_moments_max=4, ksb_sig_weight=0.0, ksb_sig_factor=1.0, failed_moments=-1000.0)[source]

A collection of parameters that govern how the HSM functions operate.

HSMParams stores a set of numbers that determine how the moments/shape estimation routines make speed/accuracy tradeoff decisions and/or store their results.

Parameters:

nsig_rg – A parameter used to optimize convolutions by cutting off the galaxy profile. In the first step of the re-Gaussianization method of PSF correction, a Gaussian approximation to the pre-seeing galaxy is calculated. If nsig_rg > 0, then this approximation is cut off at nsig_rg sigma to save computation time in convolutions. [default: 3.0]
nsig_rg2 – A parameter used to optimize convolutions by cutting off the PSF residual profile. In the re-Gaussianization method of PSF correction, a ‘PSF residual’ (the difference between the true PSF and its best-fit Gaussian approximation) is constructed. If nsig_rg2 > 0, then this PSF residual is cut off at nsig_rg2 sigma to save computation time in convolutions. [default: 3.6]
max_moment_nsig2 – No longer used. (Now calculated based on convergence_threshold.)
regauss_too_small – A parameter for how strictly the re-Gaussianization code treats small galaxies. If this parameter is 1, then the re-Gaussianization code does not impose a cut on the apparent resolution before trying to measure the PSF-corrected shape of the galaxy; if 0, then it is stricter. Using the default value of 1 prevents the re-Gaussianization PSF correction from completely failing at the beginning, before trying to do PSF correction, due to the crudest possible PSF correction (Gaussian approximation) suggesting that the galaxy is very small. This could happen for some usable galaxies particularly when they have very non-Gaussian surface brightness profiles – for example, if there’s a prominent bulge that the adaptive moments attempt to fit, ignoring a more extended disk. Setting a value of 1 is useful for keeping galaxies that would have failed for that reason. If they later turn out to be too small to really use, this will be reflected in the final estimate of the resolution factor, and they can be rejected after the fact. [default: 1]
adapt_order – The order to which circular adaptive moments should be calculated for KSB method. This parameter only affects calculations using the KSB method of PSF correction. Warning: deviating from default value of 2 results in code running more slowly, and results have not been significantly tested. [default: 2]
convergence_threshold – Accuracy (in x0, y0, and sigma, each as a fraction of sigma) when calculating adaptive moments. [default: 1.e-6]
max_mom2_iter – Maximum number of iterations to use when calculating adaptive moments. This should be sufficient in nearly all situations, with the possible exception being very flattened profiles. [default: 400]
num_iter_default – Number of iterations to report in the output ShapeData structure when code fails to converge within max_mom2_iter iterations. [default: -1]
bound_correct_wt – Maximum shift in centroids and sigma between iterations for adaptive moments. [default: 0.25]
max_amoment – Maximum value for adaptive second moments before throwing exception. Very large objects might require this value to be increased. [default: 8000]
max_ashift – Maximum allowed x / y centroid shift (units: pixels) between successive iterations for adaptive moments before throwing exception. [default: 15]
ksb_moments_max – Use moments up to ksb_moments_max order for KSB method of PSF correction. [default: 4]
ksb_sig_weight – The width of the weight function (in pixels) to use for the KSB method. Normally, this is derived from the measured moments of the galaxy image; this keyword overrides this calculation. Can be combined with ksb_sig_factor. [default: 0.0]
ksb_sig_factor – Factor by which to multiply the weight function width for the KSB method (default: 1.0). Can be combined with ksb_sig_weight. [default: 1.0]
failed_moments – Value to report for ellipticities and resolution factor if shape measurement fails. [default: -1000.]

Note

Parameters that are given in units of pixels should still be in pixels, even if one is calling FindAdaptiveMom with the option use_sky_coords=True.

After construction, all of the above are available as read-only attributes.

static check(hsmparams, default=None)[source]

Checks that hsmparams is either a valid HSMParams instance or None.

In the former case, it returns hsmparams, in the latter it returns default (HSMParams.default if no other default specified).