# Copyright (c) 2012-2023 by the GalSim developers team on GitHub
# https://github.com/GalSim-developers
#
# This file is part of GalSim: The modular galaxy image simulation toolkit.
# https://github.com/GalSim-developers/GalSim
#
# GalSim is free software: redistribution and use in source and binary forms,
# with or without modification, are permitted provided that the following
# conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice, this
# list of conditions, and the disclaimer given in the accompanying LICENSE
# file.
# 2. Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions, and the disclaimer given in the documentation
# and/or other materials provided with the distribution.
#
import numpy as np
from . import _galsim
from .position import PositionD
from .bounds import BoundsI
from .shear import Shear
from .image import Image, ImageI, ImageF, ImageD
from .errors import GalSimValueError, GalSimHSMError, GalSimIncompatibleValuesError
[docs]class ShapeData:
"""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`).
"""
def __init__(self, image_bounds=BoundsI(), moments_status=-1,
observed_shape=Shear(), moments_sigma=-1.0, moments_amp=-1.0,
moments_centroid=PositionD(), moments_rho4=-1.0, moments_n_iter=0,
correction_status=-1, corrected_e1=-10., corrected_e2=-10.,
corrected_g1=-10., corrected_g2=-10., meas_type="None",
corrected_shape_err=-1.0, correction_method="None",
resolution_factor=-1.0, psf_sigma=-1.0,
psf_shape=Shear(), error_message=""):
# Avoid empty string, which can caus problems in C++ layer.
if error_message == "": error_message = "None"
if not isinstance(image_bounds, BoundsI):
raise TypeError("image_bounds must be a BoundsI instance")
# The others will raise an appropriate TypeError from the call to _galsim.ShapeData
# when converting to int, float, etc.
self._data = _galsim.ShapeData(
image_bounds._b, int(moments_status), observed_shape.e1, observed_shape.e2,
float(moments_sigma), float(moments_amp), moments_centroid._p,
float(moments_rho4), int(moments_n_iter), int(correction_status),
float(corrected_e1), float(corrected_e2), float(corrected_g1), float(corrected_g2),
str(meas_type), float(corrected_shape_err), str(correction_method),
float(resolution_factor), float(psf_sigma), psf_shape.e1, psf_shape.e2,
str(error_message))
@property
def image_bounds(self): return BoundsI(self._data.image_bounds)
@property
def moments_status(self): return self._data.moments_status
@property
def observed_e1(self):
return self._data.observed_e1
@property
def observed_e2(self):
return self._data.observed_e2
@property
def observed_shape(self):
return Shear(e1=self.observed_e1, e2=self.observed_e2)
@property
def moments_sigma(self): return self._data.moments_sigma
@property
def moments_amp(self): return self._data.moments_amp
@property
def moments_centroid(self): return PositionD(self._data.moments_centroid)
@property
def moments_rho4(self): return self._data.moments_rho4
@property
def moments_n_iter(self): return self._data.moments_n_iter
@property
def correction_status(self): return self._data.correction_status
@property
def corrected_e1(self): return self._data.corrected_e1
@property
def corrected_e2(self): return self._data.corrected_e2
@property
def corrected_g1(self): return self._data.corrected_g1
@property
def corrected_g2(self): return self._data.corrected_g2
@property
def meas_type(self): return self._data.meas_type
@property
def corrected_shape_err(self): return self._data.corrected_shape_err
@property
def correction_method(self): return self._data.correction_method
@property
def resolution_factor(self): return self._data.resolution_factor
@property
def psf_sigma(self): return self._data.psf_sigma
@property
def psf_shape(self):
return Shear(e1=self._data.psf_e1, e2=self._data.psf_e2)
@property
def error_message(self):
# We use "None" in C++ ShapeData to indicate no error messages to avoid problems on
# (some) Macs using zero-length strings. Here, we revert that back to "".
if self._data.error_message == "None":
return ""
else:
return self._data.error_message
def __repr__(self):
s = 'galsim.hsm.ShapeData('
if self.image_bounds.isDefined(): s += 'image_bounds=%r, '%self.image_bounds
if self.moments_status != -1: s += 'moments_status=%r, '%self.moments_status
# Always include this one:
s += 'observed_shape=%r'%self.observed_shape
if self.moments_sigma != -1: s += ', moments_sigma=%r'%self.moments_sigma
if self.moments_amp != -1: s += ', moments_amp=%r'%self.moments_amp
if self.moments_centroid != PositionD():
s += ', moments_centroid=%r'%self.moments_centroid
if self.moments_rho4 != -1: s += ', moments_rho4=%r'%self.moments_rho4
if self.moments_n_iter != 0: s += ', moments_n_iter=%r'%self.moments_n_iter
if self.correction_status != -1: s += ', correction_status=%r'%self.correction_status
if self.corrected_e1 != -10.: s += ', corrected_e1=%r'%self.corrected_e1
if self.corrected_e2 != -10.: s += ', corrected_e2=%r'%self.corrected_e2
if self.corrected_g1 != -10.: s += ', corrected_g1=%r'%self.corrected_g1
if self.corrected_g2 != -10.: s += ', corrected_g2=%r'%self.corrected_g2
if self.meas_type != 'None': s += ', meas_type=%r'%self.meas_type
if self.corrected_shape_err != -1.:
s += ', corrected_shape_err=%r'%self.corrected_shape_err
if self.correction_method != 'None': s += ', correction_method=%r'%self.correction_method
if self.resolution_factor != -1.: s += ', resolution_factor=%r'%self.resolution_factor
if self.psf_sigma != -1.: s += ', psf_sigma=%r'%self.psf_sigma
if self.psf_shape != Shear(): s += ', psf_shape=%r'%self.psf_shape
if self.error_message != "": s += ', error_message=%r'%self.error_message
s += ')'
return s
def __eq__(self, other):
return (self is other or
(isinstance(other,ShapeData) and self._getinitargs() == other._getinitargs()))
def __ne__(self, other): return not self.__eq__(other)
def __hash__(self): return hash(("galsim.hsm.ShapeData", self._getinitargs()))
def _getinitargs(self):
return (self.image_bounds, self.moments_status, self.observed_shape,
self.moments_sigma, self.moments_amp, self.moments_centroid, self.moments_rho4,
self.moments_n_iter, self.correction_status, self.corrected_e1, self.corrected_e2,
self.corrected_g1, self.corrected_g2, self.meas_type, self.corrected_shape_err,
self.correction_method, self.resolution_factor, self.psf_sigma,
self.psf_shape, self.error_message)
def __getstate__(self):
return self._getinitargs()
def __setstate__(self, state):
self.__init__(*state)
[docs] def applyWCS(self, wcs, image_pos):
"""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.
"""
jac = wcs.jacobian(image_pos=image_pos)
scale, shear, theta, flip = jac.getDecomposition()
# Fix moments_sigma
moments_sigma = self.moments_sigma * scale
# Fix observed_shape
shape = self.observed_shape
# First the flip, if any.
if flip:
shape = Shear(g1 = -shape.g1, g2 = shape.g2)
# Next the rotation
shape = Shear(g = shape.g, beta = shape.beta + theta)
# Finally the shear
observed_shape = shear + shape
# Fix moments_centroid
moments_centroid = jac.toWorld(self.moments_centroid) - jac.toWorld(image_pos)
return ShapeData(image_bounds=self.image_bounds,
moments_status=self.moments_status,
observed_shape=observed_shape,
moments_sigma=moments_sigma,
moments_amp=self.moments_amp,
moments_centroid=moments_centroid,
moments_rho4=self.moments_rho4,
moments_n_iter=self.moments_n_iter,
error_message=self.error_message)
# The other values are reset to the defaults, since they are
# results from EstimateShear.
[docs]class HSMParams:
"""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.
"""
def __init__(self, nsig_rg=3.0, nsig_rg2=3.6, max_moment_nsig2=0, regauss_too_small=1,
adapt_order=2, convergence_threshold=1.e-6, max_mom2_iter=400,
num_iter_default=-1, bound_correct_wt=0.25, max_amoment=8000., max_ashift=15.,
ksb_moments_max=4, ksb_sig_weight=0.0, ksb_sig_factor=1.0, failed_moments=-1000.):
if max_moment_nsig2 != 0:
from .deprecated import depr
depr('max_moment_nsig2', 2.4, '', 'This parameter is no longer used.')
self._nsig_rg = float(nsig_rg)
self._nsig_rg2 = float(nsig_rg2)
self._regauss_too_small = int(regauss_too_small)
self._adapt_order = int(adapt_order)
self._convergence_threshold = float(convergence_threshold)
self._max_mom2_iter = int(max_mom2_iter)
self._num_iter_default = int(num_iter_default)
self._bound_correct_wt = float(bound_correct_wt)
self._max_amoment = float(max_amoment)
self._max_ashift = float(max_ashift)
self._ksb_moments_max = int(ksb_moments_max)
self._ksb_sig_weight = float(ksb_sig_weight)
self._ksb_sig_factor = float(ksb_sig_factor)
self._failed_moments = float(failed_moments)
self._make_hsmp()
def _make_hsmp(self):
self._hsmp = _galsim.HSMParams(*self._getinitargs())
def _getinitargs(self):
# TODO: For now, leave 3rd param as unused max_moment_nsig2.
# Remove it at version 3.0 to avoid changing C++ API yet.
return (self.nsig_rg, self.nsig_rg2, 0., self.regauss_too_small,
self.adapt_order, self.convergence_threshold, self.max_mom2_iter,
self.num_iter_default, self.bound_correct_wt, self.max_amoment, self.max_ashift,
self.ksb_moments_max, self.ksb_sig_weight, self.ksb_sig_factor,
self.failed_moments)
@property
def nsig_rg(self): return self._nsig_rg
@property
def nsig_rg2(self): return self._nsig_rg2
@property
def max_moment_nsig2(self): return 0.
@property
def regauss_too_small(self): return self._regauss_too_small
@property
def adapt_order(self): return self._adapt_order
@property
def convergence_threshold(self): return self._convergence_threshold
@property
def max_mom2_iter(self): return self._max_mom2_iter
@property
def num_iter_default(self): return self._num_iter_default
@property
def bound_correct_wt(self): return self._bound_correct_wt
@property
def max_amoment(self): return self._max_amoment
@property
def max_ashift(self): return self._max_ashift
@property
def ksb_moments_max(self): return self._ksb_moments_max
@property
def ksb_sig_weight(self): return self._ksb_sig_weight
@property
def ksb_sig_factor(self): return self._ksb_sig_factor
@property
def failed_moments(self): return self._failed_moments
[docs] @staticmethod
def check(hsmparams, default=None):
"""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).
"""
if hsmparams is None:
return default if default is not None else HSMParams.default
elif not isinstance(hsmparams, HSMParams):
raise TypeError("Invalid HSMParams: %s"%hsmparams)
else:
return hsmparams
def __repr__(self):
return ('galsim.hsm.HSMParams(' + 14*'%r,' + '%r)')%self._getinitargs()
def __eq__(self, other):
return (self is other or
(isinstance(other, HSMParams) and self._getinitargs() == other._getinitargs()))
def __ne__(self, other):
return not self.__eq__(other)
def __hash__(self):
return hash(('galsim.hsm.HSMParams', self._getinitargs()))
def __getstate__(self):
d = self.__dict__.copy()
del d['_hsmp']
return d
def __setstate__(self, d):
self.__dict__ = d
self._make_hsmp()
# We use the default a lot, so make it a class attribute.
HSMParams.default = HSMParams()
# A helper function that checks if the weight and the badpix bounds are
# consistent with that of the image, and that the weight is non-negative.
def _checkWeightAndBadpix(image, weight=None, badpix=None):
# Check that the weight and badpix, if given, are sensible and compatible
# with the image.
if weight is not None:
if weight.bounds != image.bounds:
raise GalSimIncompatibleValuesError(
"Weight image does not have same bounds as the input Image.",
weight=weight, image=image)
# also make sure there are no negative values
if np.any(weight.array < 0):
raise GalSimValueError("Weight image cannot contain negative values.", weight)
if badpix is not None and badpix.bounds != image.bounds:
raise GalSimIncompatibleValuesError(
"Badpix image does not have the same bounds as the input Image.",
badpix=badpix, image=image)
# A helper function for taking input weight and badpix Images, and returning a weight Image in the
# format that the C++ functions want
def _convertMask(image, weight=None, badpix=None):
# Convert from input weight and badpix images to a single mask image needed by C++ functions.
# This is used by EstimateShear() and FindAdaptiveMom().
# if no weight image was supplied, make an int array (same size as gal image) filled with 1's
if weight is None:
mask = ImageI(bounds=image.bounds, init_value=1)
else:
# if weight is an ImageI, then we can use it as the mask image:
if weight.dtype == np.int32:
if not badpix:
mask = weight
else:
# If we need to mask bad pixels, we'll need a copy anyway.
mask = ImageI(weight)
# otherwise, we need to convert it to the right type
else:
mask = ImageI(bounds=image.bounds, init_value=0)
mask.array[weight.array > 0.] = 1
# if badpix image was supplied, identify the nonzero (bad) pixels and set them to zero in weight
# image; also check bounds
if badpix is not None:
mask.array[badpix.array != 0] = 0
# if no pixels are used, raise an exception
if not np.any(mask.array):
raise GalSimHSMError("No pixels are being used!")
# finally, return the Image for the weight map
return mask
# A simpler helper function to force images to be of type ImageF or ImageD
def _convertImage(image):
# Convert the given image to the correct format needed to pass to the C++ layer.
# This is used by EstimateShear() and FindAdaptiveMom().
# if weight is not of type float/double, convert to float/double
if (image.dtype == np.int16 or image.dtype == np.uint16):
image = ImageF(image)
elif (image.dtype == np.int32 or image.dtype == np.uint32):
image = ImageD(image)
return image
[docs]def 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=1.0e-6, guess_centroid=None,
strict=True, check=True, hsmparams=None):
"""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.
"""
# prepare inputs to C++ routines: ImageF or ImageD for galaxy, PSF, and ImageI for weight map
gal_image = _convertImage(gal_image)
PSF_image = _convertImage(PSF_image)
hsmparams = HSMParams.check(hsmparams)
if check:
_checkWeightAndBadpix(gal_image, weight=weight, badpix=badpix)
weight = _convertMask(gal_image, weight=weight, badpix=badpix)
if guess_centroid is None:
guess_centroid = gal_image.true_center
try:
result = ShapeData()
_galsim.EstimateShearView(result._data,
gal_image._image, PSF_image._image, weight._image,
float(sky_var), shear_est.upper(), recompute_flux.upper(),
float(guess_sig_gal), float(guess_sig_PSF), float(precision),
guess_centroid._p, hsmparams._hsmp)
return result
except RuntimeError as err:
if (strict == True):
raise GalSimHSMError(str(err)) from None
else:
return ShapeData(error_message = str(err))
[docs]def FindAdaptiveMom(object_image, weight=None, badpix=None, guess_sig=5.0, precision=1.0e-6,
guess_centroid=None, strict=True, check=True, round_moments=False, hsmparams=None,
use_sky_coords=False):
"""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.
"""
# prepare inputs to C++ routines: ImageF or ImageD for galaxy, PSF, and ImageI for weight map
object_image = _convertImage(object_image)
hsmparams = HSMParams.check(hsmparams)
if check:
_checkWeightAndBadpix(object_image, weight=weight, badpix=badpix)
weight = _convertMask(object_image, weight=weight, badpix=badpix)
if guess_centroid is None:
guess_centroid = object_image.true_center
try:
result = ShapeData()
_galsim.FindAdaptiveMomView(result._data,
object_image._image, weight._image,
float(guess_sig), float(precision), guess_centroid._p,
bool(round_moments), hsmparams._hsmp)
if use_sky_coords:
result = result.applyWCS(object_image.wcs, image_pos=object_image.true_center)
return result
except RuntimeError as err:
if (strict == True):
raise GalSimHSMError(str(err)) from None
else:
return ShapeData(error_message = str(err))
# make FindAdaptiveMom a method of Image class
Image.FindAdaptiveMom = FindAdaptiveMom