# 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.
#
__all__ = [ 'GSObject' ]
import numpy as np
import math
import copy
from . import _galsim
from .gsparams import GSParams
from .position import _PositionD, _PositionI, Position, parse_pos_args
from ._utilities import lazy_property
from .errors import GalSimError, GalSimRangeError, GalSimValueError, GalSimIncompatibleValuesError
from .errors import GalSimFFTSizeError, GalSimNotImplementedError, convert_cpp_errors, galsim_warn
from .image import Image, ImageD, ImageF, ImageCD, ImageCF
from .shear import Shear, _Shear
from .angle import Angle
from .bounds import BoundsI, _BoundsI
from .random import BaseDeviate
from .sensor import Sensor
from .random import PoissonDeviate
from .wcs import BaseWCS, PixelScale
[docs]class GSObject:
"""Base class for all GalSim classes that represent some kind of surface brightness profile.
A GSObject is not intended to be constructed directly. Normally, you would use whatever
derived class is appropriate for the surface brightness profile you want::
>>> gal = galsim.Sersic(n=4, half_light_radius=4.3)
>>> psf = galsim.Moffat(beta=3, fwhm=2.85)
>>> conv = galsim.Convolve([gal,psf])
All of these classes are subclasses of GSObject, so you should see those docstrings for
more details about how to construct the various profiles. Here we discuss attributes and
methods that are common to all GSObjects.
GSObjects are always defined in sky coordinates. So all sizes and other linear dimensions
should be in terms of some kind of units on the sky, arcsec for instance. Only later (when
they are drawn) is the connection to pixel coordinates established via a pixel scale or WCS.
(See the documentation for galsim.BaseWCS for more details about how to specify various kinds
of world coordinate systems more complicated than a simple pixel scale.)
For instance, if you eventually draw onto an image that has a pixel scale of 0.2 arcsec/pixel,
then the normal thing to do would be to define your surface brightness profiles in terms of
arcsec and then draw with ``pixel_scale=0.2``. However, while arcsec are the usual choice of
units for the sky coordinates, if you wanted, you could instead define the sizes of all your
galaxies and PSFs in terms of radians and then use ``pixel_scale=0.2/206265`` when you draw
them.
**Transforming methods**:
The GSObject class uses an "immutable" design[1], so all methods that would potentially modify
the object actually return a new object instead. This uses pointers and such behind the
scenes, so it all happens efficiently, but it makes using the objects a bit simpler, since
you don't need to worry about some function changing your object behind your back.
In all cases below, we just give an example usage. See the docstrings for the methods for
more details about how to use them.::
>>> obj = obj.shear(shear) # Apply a shear to the object.
>>> obj = obj.dilate(scale) # Apply a flux-preserving dilation.
>>> obj = obj.magnify(mu) # Apply a surface-brightness-preserving magnification.
>>> obj = obj.rotate(theta) # Apply a rotation.
>>> obj = obj.shift(dx,dy) # Shft the object in real space.
>>> obj = obj.transform(dudx,dudy,dvdx,dvdy) # Apply a general jacobian transformation.
>>> obj = obj.lens(g1,g2,mu) # Apply both a lensing shear and magnification.
>>> obj = obj.withFlux(flux) # Set a new flux value.
>>> obj = obj * ratio # Scale the surface brightness profile by some factor.
**Access Methods**:
There are some access methods and properties that are available for all GSObjects.
Again, see the docstrings for each method for more details.::
>>> obj.flux
>>> obj.centroid
>>> obj.nyquist_scale
>>> obj.stepk
>>> obj.maxk
>>> obj.has_hard_edges
>>> obj.is_axisymmetric
>>> obj.is_analytic_x
>>> obj.is_analytic_k
>>> obj.xValue(x,y) or obj.xValue(pos)
>>> obj.kValue(kx,ky) os obj.kValue(kpos)
Most subclasses have additional methods that are available for values that are particular to
that specific surface brightness profile. e.g. ``sigma = gauss.sigma``. However, note
that class-specific methods are not available after performing one of the above transforming
operations.::
>>> gal = galsim.Gaussian(sigma=5)
>>> gal = gal.shear(g1=0.2, g2=0.05)
>>> sigma = gal.sigma # This will raise an exception.
It is however possible to access the original object that was transformed via the
``original`` attribute.::
>>> sigma = gal.original.sigma # This works.
No matter how many transformations are performed, the ``original`` attribute will contain the
_original_ object (not necessarily the most recent ancestor).
**Drawing Methods**:
The main thing to do with a GSObject once you have built it is to draw it onto an image.
There are two methods that do this. In both cases, there are lots of optional parameters.
See the docstrings for these methods for more details.::
>>> image = obj.drawImage(...)
>>> kimage = obj.drawKImage(...)
There two attributes that may be available for a GSObject.
Attributes:
original: This was mentioned above as a way to access the original object that has
been transformed by one of the transforming methods.
noise: Some types, like `RealGalaxy`, set this attribute to be the intrinsic noise that
is already inherent in the profile and will thus be present when you draw the
object. The noise is propagated correctly through the various transforming
methods, as well as convolutions and flux rescalings. Note that the ``noise``
attribute can be set directly by users even for GSObjects that do not naturally
have one. The typical use for this attribute is to use it to whiten the noise in
the image after drawing. See `BaseCorrelatedNoise` for more details.
**GSParams**:
All GSObject classes take an optional ``gsparams`` argument, so we document that feature here.
For all documentation about the specific derived classes, please see the docstring for each
one individually.
The ``gsparams`` argument can be used to specify various numbers that govern the tradeoff
between accuracy and speed for the calculations made in drawing a GSObject. The numbers are
encapsulated in a class called `GSParams`, and the user should make careful choices whenever
they opt to deviate from the defaults. For more details about the parameters and their default
values, please see the docstring of the `GSParams` class.
For example, let's say you want to do something that requires an FFT larger than 4096 x 4096
(and you have enough memory to handle it!). Then you can create a new `GSParams` object with a
larger ``maximum_fft_size`` and pass that to your GSObject on construction::
>>> gal = galsim.Sersic(n=4, half_light_radius=4.3)
>>> psf = galsim.Moffat(beta=3, fwhm=2.85)
>>> conv = galsim.Convolve([gal,psf])
>>> im = galsim.Image(1000,1000, scale=0.02) # Note the very small pixel scale!
>>> im = conv.drawImage(image=im) # This uses the default GSParams.
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "galsim/gsobject.py", line 1666, in drawImage
added_photons = prof.drawFFT(draw_image, add)
File "galsim/gsobject.py", line 1877, in drawFFT
kimage, wrap_size = self.drawFFT_makeKImage(image)
File "galsim/gsobject.py", line 1802, in drawFFT_makeKImage
raise GalSimFFTSizeError("drawFFT requires an FFT that is too large.", Nk)
galsim.errors.GalSimFFTSizeError: drawFFT requires an FFT that is too large.
The required FFT size would be 12288 x 12288, which requires 3.38 GB of memory.
If you can handle the large FFT, you may update gsparams.maximum_fft_size.
>>> big_fft_params = galsim.GSParams(maximum_fft_size=12300)
>>> conv = galsim.Convolve([gal,psf],gsparams=big_fft_params)
>>> im = conv.drawImage(image=im) # Now it works (but is slow!)
>>> im.write('high_res_sersic.fits')
Note that for compound objects such as `Convolution` or `Sum`, not all `GSParams` can be
changed when the compound object is created. In the example given here, it is possible to
change parameters related to the drawing, but not the Fourier space parameters for the
components that go into the `Convolution`. To get better sampling in Fourier space,
for example, the ``gal`` and/or ``psf`` should be created with ``gsparams`` that have a
non-default value of ``folding_threshold``. This statement applies to the threshold and
accuracy parameters.
"""
_gsparams_opt = { 'minimum_fft_size' : int,
'maximum_fft_size' : int,
'folding_threshold' : float,
'stepk_minimum_hlr' : float,
'maxk_threshold' : float,
'kvalue_accuracy' : float,
'xvalue_accuracy' : float,
'realspace_relerr' : float,
'realspace_abserr' : float,
'integration_relerr' : float,
'integration_abserr' : float,
'shoot_accuracy' : float,
'allowed_flux_variation' : float,
'range_division_for_extrema' : int,
'small_fraction_of_flux' : float
}
def __init__(self):
raise NotImplementedError("The GSObject base class should not be instantiated directly.")
# Note: subclasses are expected to define the following attributes or properties:
#
# Required for all profiles:
#
# _flux (the object's flux, natch)
# _gsparams (use GSParams.check(None) if you just want the default)
# _stepk (the sampling in k space necessary to avoid folding of image in x space)
# _maxk (the value of k beyond which aliasing can be neglected)
# _has_hard_edges (true if should use real-space convolution with another hard edge profile)
# _is_axisymmetric (true if f(x,y) = f(r))
# _is_analytic_x (true if _xValue and _drawReal are implemented)
# _is_analytic_k (true if _kValue and _drawKImage are implemented)
#
# Required and typically class attributes, but there are defaults in each case:
#
# _req_params (dict of required config parameters: name : type, default: {})
# _opt_params (dict of optional config parameters: name : type, default: {})
# _single_params (list of dicts for parameters where exactly one of several is required,
# default: [])
# _takes_rng (bool specifying whether rng is an input parameter, default: False)
#
# Optional
#
# _centroid (default = PositionD(0,0), which is often the right value)
# _positive_flux (default = _flux + _negative_flux)
# _negative_flux (default = 0; note: this should be absolute value of the negative flux)
# _max_sb (default 1.e500, which in this context is equivalent to "unknown")
# _noise (default None)
#
# In addition, subclasses should typically define most of the following methods.
# The default in each case is to raise a NotImplementedError, so if you cannot implement one,
# you may simply not define it.
#
# _xValue(self, pos)
# _kValue(self, kpos)
# _drawReal(self, image)
# _shoot(self, photons, rng):
# _drawKImage(self, image)
#
# Required for real-space convolution
#
# _sbp which must be an attribute or property providing a C++-layer SBProfile instance.
#
# Note that most objects don't need to implement real-space convolution, so use of a C++-layer
# SBProfile sub-class is usually only an implementation detail to improve efficiency.
#
# TODO: For now, _sbp is also required for transformations, but this is expected to be
# addressed in a future PR.
@property
def flux(self):
"""The flux of the profile.
"""
return self._flux
@property
def gsparams(self):
"""A `GSParams` object that sets various parameters relevant for speed/accuracy trade-offs.
"""
return self._gsparams
@property
def maxk(self):
"""The value of k beyond which aliasing can be neglected.
"""
return self._maxk
@property
def stepk(self):
"""The sampling in k space necessary to avoid folding of image in x space.
"""
return self._stepk
@property
def nyquist_scale(self):
"""The pixel spacing that does not alias maxk.
"""
return math.pi / self.maxk
@property
def has_hard_edges(self):
"""Whether there are any hard edges in the profile, which would require very small k
spacing when working in the Fourier domain.
"""
return self._has_hard_edges
@property
def is_axisymmetric(self):
"""Whether the profile is axially symmetric; affects efficiency of evaluation.
"""
return self._is_axisymmetric
@property
def is_analytic_x(self):
"""Whether the real-space values can be determined immediately at any position without
requiring a Discrete Fourier Transform.
"""
return self._is_analytic_x
@property
def is_analytic_k(self):
"""Whether the k-space values can be determined immediately at any position without
requiring a Discrete Fourier Transform.
"""
return self._is_analytic_k
@property
def centroid(self):
"""The (x, y) centroid of an object as a `PositionD`.
"""
return self._centroid
@lazy_property
def _centroid(self):
# Most profiles are centered at 0,0, so make this the default.
return _PositionD(0,0)
@property
def positive_flux(self):
"""The expectation value of flux in positive photons.
Some profiles, when rendered with photon shooting, need to shoot both positive- and
negative-flux photons. For such profiles, this method returns the total flux
of the positive-valued photons.
For profiles that don't have this complication, this is equivalent to getFlux().
It should be generally true that ``obj.positive_flux - obj.negative_flux`` returns the same
thing as ``obj.flux``. Small difference may accrue from finite numerical accuracy in
cases involving lookup tables, etc.
"""
return self._positive_flux
@property
def negative_flux(self):
"""Returns the expectation value of flux in negative photons.
Some profiles, when rendered with photon shooting, need to shoot both positive- and
negative-flux photons. For such profiles, this method returns the total absolute flux
of the negative-valued photons (i.e. as a positive value).
For profiles that don't have this complication, this returns 0.
It should be generally true that ``obj.positive_flux - obj.negative_flux`` returns the same
thing as ``obj.flux``. Small difference may accrue from finite numerical accuracy in
cases involving lookup tables, etc.
"""
return self._negative_flux
@lazy_property
def _positive_flux(self):
# The usual case.
return self.flux + self._negative_flux
@lazy_property
def _negative_flux(self):
# The usual case.
return 0.
@lazy_property
def _flux_per_photon(self):
# The usual case.
return 1.
def _calculate_flux_per_photon(self):
# If negative_flux is overriden, then _flux_per_photon should be overridden as well
# to return this calculation.
posflux = self.positive_flux
negflux = self.negative_flux
eta = negflux / (posflux + negflux)
return 1.-2.*eta
@property
def max_sb(self):
"""An estimate of the maximum surface brightness of the object.
Some profiles will return the exact peak SB, typically equal to the value of
obj.xValue(obj.centroid). However, not all profiles (e.g. Convolution) know how to
calculate this value without just drawing the image and checking what the maximum value is.
Clearly, this would be inefficient, so in these cases, some kind of estimate is returned,
which will generally be conservative on the high side.
This routine is mainly used by the photon shooting process, where an overestimate of
the maximum surface brightness is acceptable.
Note, for negative-flux profiles, this will return the absolute value of the most negative
surface brightness. Technically, it is an estimate of the maximum deviation from zero,
rather than the maximum value. For most profiles, these are the same thing.
"""
return self._max_sb
@lazy_property
def _max_sb(self):
# The way this is used, overestimates are conservative.
# So the default value of 1.e500 will skip the optimization involving the maximum sb.
return 1.e500
@property
def noise(self):
"""An estimate of the noise already in the profile.
Some profiles have some noise already in their definition. E.g. those that come from
observations of galaxies in real data. In GalSim, `RealGalaxy` objects are an example of
this. In these cases, the noise attribute gives an estimate of the Noise object that
would generate noise consistent with that already in the profile.
It is permissible to attach a noise estimate to an existing object with::
>>> obj.noise = noise # Some BaseNoise instance
"""
return self._noise
@noise.setter
def noise(self, n):
# We allow the user to set the noise with obj.noise = n
self._noise = n
@lazy_property
def _noise(self):
# Most profiles don't have any noise.
return None
# a few definitions for using GSObjects as duck-typed ChromaticObjects
@property
def separable(self): return True
@property
def interpolated(self): return False
@property
def deinterpolated(self): return self
@property
def sed(self):
return sed.SED(self.flux, 'nm', '1')
@property
def SED(self):
from .deprecated import depr
depr('obj.SED', 2.5, 'obj.sed')
return self.sed
@property
def spectral(self): return False
@property
def dimensionless(self): return True
@property
def wave_list(self): return np.array([], dtype=float)
@property
def redshift(self): return getattr(self, '_redshift', 0.)
# Also need these methods to duck-type as a ChromaticObject
def evaluateAtWavelength(self, wave):
return self
def _approxWavelength(self, wave):
return wave, self
# Make op+ of two GSObjects work to return an Add object
# Note: we don't define __iadd__ and similar. Let python handle this automatically
# to make obj += obj2 be equivalent to obj = obj + obj2.
[docs] def __add__(self, other):
"""Add two GSObjects.
Equivalent to Add(self, other)
"""
return _sum.Add([self, other])
# op- is unusual, but allowed. It subtracts off one profile from another.
[docs] def __sub__(self, other):
"""Subtract two GSObjects.
Equivalent to Add(self, -1 * other)
"""
return _sum.Add([self, (-1. * other)])
# Make op* work to adjust the flux of an object
[docs] def __mul__(self, other):
"""Scale the flux of the object by the given factor.
obj * flux_ratio is equivalent to obj.withScaledFlux(flux_ratio)
It creates a new object that has the same profile as the original, but with the
surface brightness at every location scaled by the given amount.
You can also multiply by an `SED`, which will create a `ChromaticObject` where the `SED`
acts like a wavelength-dependent ``flux_ratio``.
"""
return self.withScaledFlux(other)
[docs] def __rmul__(self, other):
"""Equivalent to obj * other. See `__mul__` for details."""
return self.__mul__(other)
# Likewise for op/
[docs] def __div__(self, other):
"""Equivalent to obj * (1/other). See `__mul__` for details."""
return self * (1. / other)
__truediv__ = __div__
def __neg__(self):
return -1. * self
# Some calculations that can be done for all GSObjects.
[docs] def calculateHLR(self, size=None, scale=None, centroid=None, flux_frac=0.5):
"""Returns the half-light radius of the object.
If the profile has a half_light_radius attribute, it will just return that, but in the
general case, we draw the profile and estimate the half-light radius directly.
This function (by default at least) is only accurate to a few percent, typically.
Possibly worse depending on the profile being measured. If you care about a high
precision estimate of the half-light radius, the accuracy can be improved using the
optional parameter scale to change the pixel scale used to draw the profile.
The default scale is half the Nyquist scale, which were found to produce results accurate
to a few percent on our internal tests. Using a smaller scale will be more accurate at
the expense of speed.
In addition, you can optionally specify the size of the image to draw. The default size is
None, which means `drawImage` will choose a size designed to contain around 99.5% of the
flux. This is overkill for this calculation, so choosing a smaller size than this may
speed up this calculation somewhat.
Also, while the name of this function refers to the half-light radius, in fact it can also
calculate radii that enclose other fractions of the light, according to the parameter
``flux_frac``. E.g. for r90, you would set flux_frac=0.90.
The default scale should usually be acceptable for things like testing that a galaxy
has a reasonable resolution, but they should not be trusted for very fine grain
discriminations.
Parameters:
size: If given, the stamp size to use for the drawn image. [default: None,
which will let `drawImage` choose the size automatically]
scale: If given, the pixel scale to use for the drawn image. [default:
0.5 * self.nyquist_scale]
centroid: The position to use for the centroid. [default: self.centroid]
flux_frac: The fraction of light to be enclosed by the returned radius.
[default: 0.5]
Returns:
an estimate of the half-light radius in physical units
"""
try:
# It there is a half_light_radius attribute, use that.
return self.half_light_radius
except (AttributeError, GalSimError):
# Otherwise, or (e.g. with Airy where it is only implemented for obscuration=0)
# if there is an error trying to use it, then keep going with this calculation.
pass
if scale is None:
scale = self.nyquist_scale * 0.5
if centroid is None:
centroid = self.centroid
# Draw the image. Note: need a method that integrates over pixels to get flux right.
# The offset is to make all the rsq values different to help the precision a bit.
offset = _PositionD(0.2, 0.33)
im = self.drawImage(nx=size, ny=size, scale=scale, offset=offset, dtype=float)
center = im.true_center + offset + centroid/scale
return im.calculateHLR(center=center, flux=self.flux, flux_frac=flux_frac)
[docs] def calculateMomentRadius(self, size=None, scale=None, centroid=None, rtype='det'):
"""Returns an estimate of the radius based on unweighted second moments.
The second moments are defined as:
Q_ij = int( I(x,y) i j dx dy ) / int( I(x,y) dx dy )
where i,j may be either x or y.
If I(x,y) is a Gaussian, then T = Tr(Q) = Qxx + Qyy = 2 sigma^2. Thus, one reasonable
choice for a "radius" for an arbitrary profile is sqrt(T/2).
In addition, det(Q) = sigma^4. So another choice for an arbitrary profile is det(Q)^1/4.
This routine can return either of these measures according to the value of the ``rtype``
parameter. ``rtype='trace'`` will cause it to return sqrt(T/2). ``rtype='det'`` will cause
it to return det(Q)^1/4. And ``rtype='both'`` will return a tuple with both values.
Note that for the special case of a Gaussian profile, no calculation is necessary, and
the ``sigma`` attribute will be used in both cases. In the limit as scale->0, this
function will return the same value, but because finite pixels are drawn, the results
will not be precisely equal for real use cases. The approximation being made is that
the integral of I(x,y) i j dx dy over each pixel can be approximated as
int(I(x,y) dx dy) * i_center * j_center.
This function (by default at least) is only accurate to a few percent, typically.
Possibly worse depending on the profile being measured. If you care about a high
precision estimate of the radius, the accuracy can be improved using the optional
parameters size and scale to change the size and pixel scale used to draw the profile.
The default is to use the the Nyquist scale for the pixel scale and let `drawImage`
choose a size for the stamp that will enclose at least 99.5% of the flux. These
were found to produce results accurate to a few percent on our internal tests.
Using a smaller scale and larger size will be more accurate at the expense of speed.
The default parameters should usually be acceptable for things like testing that a galaxy
has a reasonable resolution, but they should not be trusted for very fine grain
discriminations. For a more accurate estimate, see galsim.hsm.FindAdaptiveMom.
Parameters:
size: If given, the stamp size to use for the drawn image. [default: None,
which will let `drawImage` choose the size automatically]
scale: If given, the pixel scale to use for the drawn image. [default:
self.nyquist_scale]
centroid: The position to use for the centroid. [default: self.centroid]
rtype: There are three options for this parameter:
- 'trace' means return sqrt(T/2)
- 'det' means return det(Q)^1/4
- 'both' means return both: (sqrt(T/2), det(Q)^1/4)
[default: 'det']
Returns:
an estimate of the radius in physical units (or both estimates if rtype == 'both')
"""
if rtype not in ('trace', 'det', 'both'):
raise GalSimValueError("Invalid rtype.", rtype, ('trace', 'det', 'both'))
if hasattr(self, 'sigma'):
if rtype == 'both':
return self.sigma, self.sigma
else:
return self.sigma
if scale is None:
scale = self.nyquist_scale
if centroid is None:
centroid = self.centroid
# Draw the image. Note: need a method that integrates over pixels to get flux right.
im = self.drawImage(nx=size, ny=size, scale=scale, dtype=float)
center = im.true_center + centroid/scale
return im.calculateMomentRadius(center=center, flux=self.flux, rtype=rtype)
[docs] def calculateFWHM(self, size=None, scale=None, centroid=None):
"""Returns the full-width half-maximum (FWHM) of the object.
If the profile has a fwhm attribute, it will just return that, but in the general case,
we draw the profile and estimate the FWHM directly.
As with `calculateHLR` and `calculateMomentRadius`, this function optionally takes size and
scale values to use for the image drawing. The default is to use the the Nyquist scale
for the pixel scale and let `drawImage` choose a size for the stamp that will enclose at
least 99.5% of the flux. These were found to produce results accurate to well below
one percent on our internal tests, so it is unlikely that you will want to adjust
them for accuracy. However, using a smaller size than default could help speed up
the calculation, since the default is usually much larger than is needed.
Parameters:
size: If given, the stamp size to use for the drawn image. [default: None,
which will let `drawImage` choose the size automatically]
scale: If given, the pixel scale to use for the drawn image. [default:
self.nyquist_scale]
centroid: The position to use for the centroid. [default: self.centroid]
Returns:
an estimate of the full-width half-maximum in physical units
"""
if hasattr(self, 'fwhm'):
return self.fwhm
if scale is None:
scale = self.nyquist_scale
if centroid is None:
centroid = self.centroid
# Draw the image. Note: draw with method='sb' here, since the fwhm is a property of the
# raw surface brightness profile, not integrated over pixels.
# The offset is to make all the rsq values different to help the precision a bit.
offset = _PositionD(0.2, 0.33)
im = self.drawImage(nx=size, ny=size, scale=scale, offset=offset, method='sb', dtype=float)
# Get the maximum value, assuming the maximum is at the centroid.
if self.is_analytic_x:
Imax = self.xValue(centroid)
else:
im1 = self.drawImage(nx=1, ny=1, scale=scale, method='sb', offset=-centroid/scale)
Imax = im1(1,1)
center = im.true_center + offset + centroid/scale
return im.calculateFWHM(center=center, Imax=Imax)
[docs] def xValue(self, *args, **kwargs):
"""Returns the value of the object at a chosen 2D position in real space.
This function returns the surface brightness of the object at a particular position
in real space. The position argument may be provided as a `PositionD` or `PositionI`
argument, or it may be given as x,y (either as a tuple or as two arguments).
The object surface brightness profiles are typically defined in world coordinates, so
the position here should be in world coordinates as well.
Not all `GSObject` classes can use this method. Classes like Convolution that require a
Discrete Fourier Transform to determine the real space values will not do so for a single
position. Instead a GalSimError will be raised. The xValue() method is available if and
only if ``obj.is_analytic_x == True``.
Users who wish to use the xValue() method for an object that is the convolution of other
profiles can do so by drawing the convolved profile into an image, using the image to
initialize a new `InterpolatedImage`, and then using the xValue() method for that new
object.
Parameters:
position: The position at which you want the surface brightness of the object.
Returns:
the surface brightness at that position.
"""
pos = parse_pos_args(args,kwargs,'x','y')
return self._xValue(pos)
[docs] def _xValue(self, pos):
"""Equivalent to `xValue`, but ``pos`` must be a `galsim.PositionD` instance
Parameters:
pos: The position at which you want the surface brightness of the object.
Returns:
the surface brightness at that position.
"""
raise NotImplementedError("%s does not implement xValue"%self.__class__.__name__)
[docs] def kValue(self, *args, **kwargs):
"""Returns the value of the object at a chosen 2D position in k space.
This function returns the amplitude of the fourier transform of the surface brightness
profile at a given position in k space. The position argument may be provided as a
`Position` argument, or it may be given as kx,ky (either as a tuple or as two arguments).
Technically, kValue() is available if and only if the given obj has ``obj.is_analytic_k
== True``, but this is the case for all `GSObject` classes currently, so that should never
be an issue (unlike for `xValue`).
Parameters:
position: The position in k space at which you want the fourier amplitude.
Returns:
the amplitude of the fourier transform at that position.
"""
kpos = parse_pos_args(args,kwargs,'kx','ky')
return self._kValue(kpos)
[docs] def _kValue(self, kpos): # pragma: no cover (all our classes override this)
"""Equivalent to `kValue`, but ``kpos`` must be a `galsim.PositionD` instance.
"""
raise NotImplementedError("%s does not implement kValue"%self.__class__.__name__)
[docs] def withGSParams(self, gsparams=None, **kwargs):
"""Create a version of the current object with the given `GSParams`.
You may either provide a `GSParams` instance::
>>> gsparams = galsim.GSParams(folding_threshold=1.e-4, maxk_threshold=1.e-4)
>>> obj = obj.withGSParams(gsparams)
or one or more named parameters as keyword arguments::
>>> obj = obj.withGSParams(folding_threshold=1.e-4, maxk_threshold=1.e-4)
.. note::
The latter style will leave all non-named parameters at their current
values. It only updates the named parameters to the given values.
"""
# Note to developers: objects that wrap other objects should override this in order
# to apply the new gsparams to the components.
# This implementation relies on getstate/setstate clearing out any _sbp or similar
# attribute that depends on the details of gsparams. If there are stored calculations
# aside from these, you should also clear them as well, or update them.
if gsparams == self.gsparams: return self
ret = copy.copy(self)
ret._gsparams = GSParams.check(gsparams, self.gsparams, **kwargs)
return ret
[docs] def withFlux(self, flux):
"""Create a version of the current object with a different flux.
This function is equivalent to ``obj.withScaledFlux(flux / obj.flux)``.
It creates a new object that has the same profile as the original, but with the
surface brightness at every location rescaled such that the total flux will be
the given value. Note that if ``flux`` is an `SED`, the return value will be a
`ChromaticObject` with specified `SED`.
Parameters:
flux: The new flux for the object.
Returns:
the object with the new flux
"""
return self.withScaledFlux(flux / self.flux)
[docs] def withScaledFlux(self, flux_ratio):
"""Create a version of the current object with the flux scaled by the given ``flux_ratio``.
This function is equivalent to ``obj.withFlux(flux_ratio * obj.flux)``. Indeed, withFlux()
is implemented in terms of this one.
It creates a new object that has the same profile as the original, but with the
surface brightness at every location scaled by the given amount. If ``flux_ratio`` is an
`SED`, then the returned object is a `ChromaticObject` with the `SED` multiplied by
its current ``flux``.
Note that in this case the ``flux`` attribute of the `GSObject` being scaled gets
interpreted as being dimensionless, instead of having its normal units of [photons/s/cm^2].
The photons/s/cm^2 units are (optionally) carried by the `SED` instead, or even left out
entirely if the `SED` is dimensionless itself (see discussion in the `ChromaticObject`
docstring). The `GSObject` ``flux`` attribute *does* still contribute to the
`ChromaticObject` normalization, though. For example, the following are equivalent::
>>> chrom_obj = gsobj.withScaledFlux(sed * 3.0)
>>> chrom_obj2 = (gsobj * 3.0).withScaledFlux(sed)
An equivalent, and usually simpler, way to effect this scaling is::
>>> obj = obj * flux_ratio
Parameters:
flux_ratio: The ratio by which to rescale the flux of the object when creating a new
one.
Returns:
the object with the new flux.
"""
# Prohibit non-SED callable flux_ratio here as most likely an error.
if hasattr(flux_ratio, '__call__') and not isinstance(flux_ratio, sed.SED):
raise TypeError('callable flux_ratio must be an SED.')
if flux_ratio == 1:
return self
else:
return transform.Transform(self, flux_ratio=flux_ratio)
[docs] def expand(self, scale):
"""Expand the linear size of the profile by the given ``scale`` factor, while preserving
surface brightness.
e.g. ``half_light_radius`` <-- ``half_light_radius * scale``
This doesn't correspond to either of the normal operations one would typically want to do to
a galaxy. The functions dilate() and magnify() are the more typical usage. But this
function is conceptually simple. It rescales the linear dimension of the profile, while
preserving surface brightness. As a result, the flux will necessarily change as well.
See dilate() for a version that applies a linear scale factor while preserving flux.
See magnify() for a version that applies a scale factor to the area while preserving surface
brightness.
Parameters:
scale: The factor by which to scale the linear dimension of the object.
Returns:
the expanded object.
"""
return transform.Transform(self, jac=[scale, 0., 0., scale])
[docs] def dilate(self, scale):
"""Dilate the linear size of the profile by the given ``scale`` factor, while preserving
flux.
e.g. ``half_light_radius`` <-- ``half_light_radius * scale``
See expand() and magnify() for versions that preserve surface brightness, and thus
changes the flux.
Parameters:
scale: The linear rescaling factor to apply.
Returns:
the dilated object.
"""
# equivalent to self.expand(scale) * (1./scale**2)
return transform.Transform(self, jac=[scale, 0., 0., scale], flux_ratio=scale**-2)
[docs] def magnify(self, mu):
"""Create a version of the current object with a lensing magnification applied to it,
scaling the area and flux by ``mu`` at fixed surface brightness.
This process applies a lensing magnification mu, which scales the linear dimensions of the
image by the factor sqrt(mu), i.e., ``half_light_radius`` <--
``half_light_radius * sqrt(mu)`` while increasing the flux by a factor of mu. Thus,
magnify() preserves surface brightness.
See dilate() for a version that applies a linear scale factor while preserving flux.
See expand() for a version that applies a linear scale factor while preserving surface
brightness.
Parameters:
mu: The lensing magnification to apply.
Returns:
the magnified object.
"""
return self.expand(math.sqrt(mu))
[docs] def shear(self, *args, **kwargs):
"""Create a version of the current object with an area-preserving shear applied to it.
The arguments may be either a `Shear` instance or arguments to be used to initialize one.
For more details about the allowed keyword arguments, see the `Shear` docstring.
The shear() method precisely preserves the area. To include a lensing distortion with
the appropriate change in area, either use shear() with magnify(), or use lens(), which
combines both operations.
Parameters:
shear: The `Shear` to be applied. Or, as described above, you may instead supply
parameters do construct a shear directly. eg. ``obj.shear(g1=g1,g2=g2)``.
Returns:
the sheared object.
"""
if len(args) == 1:
if kwargs:
raise TypeError("Error, gave both unnamed and named arguments to GSObject.shear!")
if not isinstance(args[0], Shear):
raise TypeError("Error, unnamed argument to GSObject.shear is not a Shear!")
shear = args[0]
elif len(args) > 1:
raise TypeError("Error, too many unnamed arguments to GSObject.shear!")
elif len(kwargs) == 0:
raise TypeError("Error, shear argument is required")
else:
shear = Shear(**kwargs)
return transform.Transform(self, shear.getMatrix())
[docs] def _shear(self, shear):
"""Equivalent to `GSObject.shear`, but without the overhead of sanity checks or other
ways to input the ``shear`` value.
Also, it won't propagate any noise attribute.
Parameters:
shear: The `Shear` to be applied.
Returns:
the sheared object.
"""
return transform._Transform(self, shear.getMatrix())
[docs] def lens(self, g1, g2, mu):
"""Create a version of the current object with both a lensing shear and magnification
applied to it.
This `GSObject` method applies a lensing (reduced) shear and magnification. The shear must
be specified using the g1, g2 definition of shear (see `Shear` for more details).
This is the same definition as the outputs of the PowerSpectrum and NFWHalo classes, which
compute shears according to some lensing power spectrum or lensing by an NFW dark matter
halo. The magnification determines the rescaling factor for the object area and flux,
preserving surface brightness.
Parameters:
g1: First component of lensing (reduced) shear to apply to the object.
g2: Second component of lensing (reduced) shear to apply to the object.
mu: Lensing magnification to apply to the object. This is the factor by which
the solid angle subtended by the object is magnified, preserving surface
brightness.
Returns:
the lensed object.
"""
shear = Shear(g1=g1, g2=g2)
return transform.Transform(self, shear.getMatrix() * math.sqrt(mu))
def _lens(self, g1, g2, mu):
"""Equivalent to `GSObject.lens`, but without the overhead of some of the sanity checks.
Also, it won't propagate any noise attribute.
Parameters:
g1: First component of lensing (reduced) shear to apply to the object.
g2: Second component of lensing (reduced) shear to apply to the object.
mu: Lensing magnification to apply to the object. This is the factor by which
the solid angle subtended by the object is magnified, preserving surface
brightness.
Returns:
the lensed object.
"""
shear = _Shear(g1 + 1j*g2)
return transform._Transform(self, shear.getMatrix() * math.sqrt(mu))
[docs] def rotate(self, theta):
"""Rotate this object by an `Angle` ``theta``.
Parameters:
theta: Rotation angle (`Angle` object, positive means anticlockwise).
Returns:
the rotated object.
"""
if not isinstance(theta, Angle):
raise TypeError("Input theta should be an Angle")
s, c = theta.sincos()
return transform.Transform(self, jac=[c, -s, s, c])
[docs] def shift(self, *args, **kwargs):
"""Create a version of the current object shifted by some amount in real space.
After this call, the caller's type will be a `GSObject`.
This means that if the caller was a derived type that had extra methods or properties
beyond those defined in `GSObject` (e.g. `Gaussian.sigma`), then these methods are no
longer available.
Note: in addition to the dx,dy parameter names, you may also supply dx,dy as a tuple,
or as a `Position` object.
The shift coordinates here are sky coordinates. `GSObject` profiles are always defined in
sky coordinates and only later (when they are drawn) is the connection to pixel coordinates
established (via a pixel_scale or WCS). So a shift of dx moves the object horizontally
in the sky (e.g. west in the local tangent plane of the observation), and dy moves the
object vertically (north in the local tangent plane).
The units are typically arcsec, but we don't enforce that anywhere. The units here just
need to be consistent with the units used for any size values used by the `GSObject`.
The connection of these units to the eventual image pixels is defined by either the
``pixel_scale`` or the ``wcs`` parameter of `GSObject.drawImage`.
Note: if you want to shift the object by a set number (or fraction) of pixels in the
drawn image, you probably want to use the ``offset`` parameter of `GSObject.drawImage`
rather than this method.
Parameters:
dx: Horizontal shift to apply.
dy: Vertical shift to apply.
Alternatively, you may supply a single parameter as a `Position` instance, rather than
the two components separately if that is more convenient.
Parameter:
offset: The shift to apply, given as PositionD(dx,dy) or PositionI(dx,dy)
Returns:
the shifted object.
"""
offset = parse_pos_args(args, kwargs, 'dx', 'dy')
return transform.Transform(self, offset=offset)
[docs] def _shift(self, dx, dy):
"""Equivalent to `shift`, but without the overhead of sanity checks or option
to give the shift as a PositionD.
Also, it won't propagate any noise attribute.
Parameters:
dx: The x-component of the shift to apply
dy: The y-component of the shift to apply
Returns:
the shifted object.
"""
new_obj = transform._Transform(self, offset=(dx,dy))
return new_obj
[docs] def atRedshift(self, redshift):
"""Create a version of the current object with a different redshift.
For regular GSObjects, this method doesn't do anything aside from setting a ``redshift``
attribute with the given value. But this allows duck typing with ChromaticObjects
where this function will adjust the SED appropriately.
Returns:
the object with the new redshift
"""
ret = copy.copy(self)
ret._redshift = redshift
return ret
# Make sure the image is defined with the right size and wcs for drawImage()
def _setup_image(self, image, nx, ny, bounds, add_to_image, dtype, center, odd=False):
# If image is given, check validity of nx,ny,bounds:
if image is not None:
if bounds is not None:
raise GalSimIncompatibleValuesError(
"Cannot provide bounds if image is provided", bounds=bounds, image=image)
if nx is not None or ny is not None:
raise GalSimIncompatibleValuesError(
"Cannot provide nx,ny if image is provided", nx=nx, ny=ny, image=image)
if dtype is not None and image.array.dtype != dtype:
raise GalSimIncompatibleValuesError(
"Cannot specify dtype != image.array.dtype if image is provided",
dtype=dtype, image=image)
# Resize the given image if necessary
if not image.bounds.isDefined():
# Can't add to image if need to resize
if add_to_image:
raise GalSimIncompatibleValuesError(
"Cannot add_to_image if image bounds are not defined",
add_to_image=add_to_image, image=image)
N = self.getGoodImageSize(1.0)
if odd: N += 1
bounds = _BoundsI(1,N,1,N)
image.resize(bounds)
# Else use the given image as is
# Otherwise, make a new image
else:
# Can't add to image if none is provided.
if add_to_image:
raise GalSimIncompatibleValuesError(
"Cannot add_to_image if image is None", add_to_image=add_to_image, image=image)
# Use bounds or nx,ny if provided
if bounds is not None:
if nx is not None or ny is not None:
raise GalSimIncompatibleValuesError(
"Cannot set both bounds and (nx, ny)", nx=nx, ny=ny, bounds=bounds)
if not bounds.isDefined():
raise GalSimValueError("Cannot use undefined bounds", bounds)
image = Image(bounds=bounds, dtype=dtype)
elif nx is not None or ny is not None:
if nx is None or ny is None:
raise GalSimIncompatibleValuesError(
"Must set either both or neither of nx, ny", nx=nx, ny=ny)
image = Image(nx, ny, dtype=dtype)
if center is not None:
image.shift(_PositionI(np.floor(center.x+0.5-image.true_center.x),
np.floor(center.y+0.5-image.true_center.y)))
else:
N = self.getGoodImageSize(1.0)
if odd: N += 1
image = Image(N, N, dtype=dtype)
if center is not None:
image.setCenter(_PositionI(np.ceil(center.x), np.ceil(center.y)))
return image
def _local_wcs(self, wcs, image, offset, center, use_true_center, new_bounds):
# Get the local WCS at the location of the object.
if wcs._isUniform:
return wcs.local()
elif image is None:
bounds = new_bounds
else:
bounds = image.bounds
if not bounds.isDefined():
raise GalSimIncompatibleValuesError(
"Cannot provide non-local wcs with automatically sized image",
wcs=wcs, image=image, bounds=new_bounds)
elif center is not None:
obj_cen = center
elif use_true_center:
obj_cen = bounds.true_center
else:
obj_cen = bounds.center
# Convert from PositionI to PositionD
obj_cen = _PositionD(obj_cen.x, obj_cen.y)
# _parse_offset has already turned offset=None into PositionD(0,0), so it is safe to add.
obj_cen += offset
return wcs.local(image_pos=obj_cen)
def _parse_offset(self, offset):
if offset is None:
return _PositionD(0,0)
elif isinstance(offset, Position):
return _PositionD(offset.x, offset.y)
else:
# Let python raise the appropriate exception if this isn't valid.
return _PositionD(offset[0], offset[1])
def _parse_center(self, center):
# Almost the same as _parse_offset, except we leave it as None in that case.
if center is None:
return None
elif isinstance(center, Position):
return _PositionD(center.x, center.y)
else:
# Let python raise the appropriate exception if this isn't valid.
return _PositionD(center[0], center[1])
def _get_new_bounds(self, image, nx, ny, bounds, center):
if image is not None and image.bounds.isDefined():
return image.bounds
elif nx is not None and ny is not None:
b = BoundsI(1,nx,1,ny)
if center is not None:
b = b.shift(_PositionI(np.floor(center.x+0.5)-b.center.x,
np.floor(center.y+0.5)-b.center.y))
return b
elif bounds is not None and bounds.isDefined():
return bounds
else:
return BoundsI()
def _adjust_offset(self, new_bounds, offset, center, use_true_center):
# Note: this assumes self is in terms of image coordinates.
if center is not None:
if new_bounds.isDefined():
offset += center - new_bounds.center
else:
# Then will be created as even sized image.
offset += _PositionD(center.x-np.ceil(center.x), center.y-np.ceil(center.y))
elif use_true_center:
# For even-sized images, the SBProfile draw function centers the result in the
# pixel just up and right of the real center. So shift it back to make sure it really
# draws in the center.
# Also, remember that numpy's shape is ordered as [y,x]
dx = offset.x
dy = offset.y
shape = new_bounds.numpyShape()
if shape[1] % 2 == 0: dx -= 0.5
if shape[0] % 2 == 0: dy -= 0.5
offset = _PositionD(dx,dy)
return offset
def _determine_wcs(self, scale, wcs, image, default_wcs=None):
# Determine the correct wcs given the input scale, wcs and image.
if wcs is not None:
if scale is not None:
raise GalSimIncompatibleValuesError(
"Cannot provide both wcs and scale", wcs=wcs, scale=scale)
if not isinstance(wcs, BaseWCS):
raise TypeError("wcs must be a BaseWCS instance")
if image is not None: image.wcs = None
elif scale is not None:
wcs = PixelScale(scale)
if image is not None: image.wcs = None
elif image is not None and image.wcs is not None:
wcs = image.wcs
# If the input scale <= 0, or wcs is still None at this point, then use the Nyquist scale:
if wcs is None or (wcs._isPixelScale and not wcs.scale > 0):
if default_wcs is None:
wcs = PixelScale(self.nyquist_scale)
else:
wcs = default_wcs
return wcs
def _prepareDraw(self):
# Do any work that was postponed until drawImage.
pass
[docs] def drawImage(self, image=None, nx=None, ny=None, bounds=None, scale=None, wcs=None, dtype=None,
method='auto', area=1., exptime=1., gain=1., add_to_image=False,
center=None, use_true_center=True, offset=None,
n_photons=0., rng=None, max_extra_noise=0.,
poisson_flux=None, sensor=None, photon_ops=None, n_subsample=3, maxN=None,
save_photons=False, bandpass=None, setup_only=False,
surface_ops=None):
"""Draws an `Image` of the object.
The drawImage() method is used to draw an `Image` of the current object using one of several
possible rendering methods (see below). It can create a new `Image` or can draw onto an
existing one if provided by the ``image`` parameter. If the ``image`` is given, you can
also optionally add to the given `Image` if ``add_to_image = True``, but the default is to
replace the current contents with new values.
**Providing an input image**:
Note that if you provide an ``image`` parameter, it is the image onto which the profile
will be drawn. The provided image *will be modified*. A reference to the same image
is also returned to provide a parallel return behavior to when ``image`` is ``None``
(described above).
This option is useful in practice because you may want to construct the image first and
then draw onto it, perhaps multiple times. For example, you might be drawing onto a
subimage of a larger image. Or you may want to draw different components of a complex
profile separately. In this case, the returned value is typically ignored. For example::
>>> im1 = bulge.drawImage()
>>> im2 = disk.drawImage(image=im1, add_to_image=True)
>>> assert im1 is im2
>>> full_image = galsim.Image(2048, 2048, scale=pixel_scale)
>>> b = galsim.BoundsI(x-32, x+32, y-32, y+32)
>>> stamp = obj.drawImage(image = full_image[b])
>>> assert (stamp.array == full_image[b].array).all()
**Letting drawImage create the image for you**:
If drawImage() will be creating the image from scratch for you, then there are several ways
to control the size of the new image. If the ``nx`` and ``ny`` keywords are present, then
an image with these numbers of pixels on a side will be created. Similarly, if the ``bounds``
keyword is present, then an image with the specified bounds will be created. Note that it
is an error to provide an existing `Image` when also specifying ``nx``, ``ny``, or
``bounds``. In the absence of ``nx``, ``ny``, and ``bounds``, drawImage will decide a good
size to use based on the size of the object being drawn. Basically, it will try to use an
area large enough to include at least 99.5% of the flux.
.. note::
This value 0.995 is really ``1 - folding_threshold``. You can change the value of
``folding_threshold`` for any object via `GSParams`.
You can set the pixel scale of the constructed image with the ``scale`` parameter, or set
a WCS function with ``wcs``. If you do not provide either ``scale`` or ``wcs``, then
drawImage() will default to using the Nyquist scale for the current object.
You can also set the data type used in the new `Image` with the ``dtype`` parameter that has
the same options as for the `Image` constructor.
**The drawing "method"**:
There are several different possible methods drawImage() can use for rendering the image.
This is set by the ``method`` parameter. The options are:
auto
This is the default, which will normally be equivalent to 'fft'. However,
if the object being rendered is simple (no convolution) and has hard edges
(e.g. a Box or a truncated Moffat or Sersic), then it will switch to
'real_space', since that is often both faster and more accurate in these
cases (due to ringing in Fourier space).
fft
The integration of the light within each pixel is mathematically equivalent
to convolving by the pixel profile (a `Pixel` object) and sampling the result
at the centers of the pixels. This method will do that convolution using
a discrete Fourier transform. Furthermore, if the object (or any component
of it) has been transformed via shear(), dilate(), etc., then these
transformations are done in Fourier space as well.
real_space
This uses direct integrals (using the Gauss-Kronrod-Patterson algorithm)
in real space for the integration over the pixel response. It is usually
slower than the 'fft' method, but if the profile has hard edges that cause
ringing in Fourier space, it can be faster and/or more accurate. If you
use 'real_space' with something that is already a Convolution, then this
will revert to 'fft', since the double convolution that is required to also
handle the pixel response is far too slow to be practical using real-space
integrals.
phot
This uses a technique called photon shooting to render the image.
Essentially, the object profile is taken as a probability distribution
from which a finite number of photons are "shot" onto the image. Each
photon's flux gets added to whichever pixel the photon hits. This process
automatically accounts for the integration of the light over the pixel
area, since all photons that hit any part of the pixel are counted.
Convolutions and transformations are simple geometric processes in this
framework. However, there are two caveats with this method: (1) the
resulting image will have Poisson noise from the finite number of photons,
and (2) it is not available for all object types (notably anything that
includes a Deconvolution).
no_pixel
Instead of integrating over the pixels, this method will sample the profile
at the centers of the pixels and multiply by the pixel area. If there is
a convolution involved, the choice of whether this will use an FFT or
real-space calculation is governed by the ``real_space`` parameter of the
Convolution class. This method is the appropriate choice if you are using
a PSF that already includes a convolution by the pixel response. For
example, if you are using a PSF from an observed image of a star, then it
has already been convolved by the pixel, so you would not want to do so
again. Note: The multiplication by the pixel area gets the flux
normalization right for the above use case. cf. ``method = 'sb'``.
sb
This is a lot like 'no_pixel', except that the image values will simply be
the sampled object profile's surface brightness, not multiplied by the
pixel area. This does not correspond to any real observing scenario, but
it could be useful if you want to view the surface brightness profile of an
object directly, without including the pixel integration.
The 'phot' method has a few extra parameters that adjust how it functions. The total
number of photons to shoot is normally calculated from the object's flux. This flux is
taken to be given in photons/cm^2/s, so for most simple profiles, this times area * exptime
will equal the number of photons shot. (See the discussion in Rowe et al, 2015, for why
this might be modified for `InterpolatedImage` and related profiles.) However, you can
manually set a different number of photons with ``n_photons``. You can also set
``max_extra_noise`` to tell drawImage() to use fewer photons than normal (and so is faster)
such that no more than that much extra noise is added to any pixel. This is particularly
useful if you will be subsequently adding sky noise, and you can thus tolerate more noise
than the normal number of photons would give you, since using fewer photons is of course
faster. Finally, the default behavior is to have the total flux vary as a Poisson random
variate, which is normally appropriate with photon shooting. But you can turn this off with
``poisson_flux=False``. It also defaults to False if you set an explicit value for
``n_photons``.
Given the periodicity implicit in the use of FFTs, there can occasionally be artifacts due
to wrapping at the edges, particularly for objects that are quite extended (e.g., due to
the nature of the radial profile). See `GSParams` for parameters that you can use to reduce
the level of these artifacts, in particular ``folding_threshold`` may be helpful if you see
such artifacts in your images.
Setting the offset:
The object will by default be drawn with its nominal center at the center location of the
image. There is thus a qualitative difference in the appearance of the rendered profile
when drawn on even- and odd-sized images. For a profile with a maximum at (0,0), this
maximum will fall in the central pixel of an odd-sized image, but in the corner of the four
central pixels of an even-sized image. There are three parameters that can affect this
behavior. First, you can specify any arbitrary pixel position to center the object using
the ``center`` parameter. If this is None, then it will pick one of the two potential
"centers" of the image, either ``image.true_center`` or ``image.center``. The latter is
an integer position, which always corresponds to the center of some pixel, which for even
sized images won't (cannot) be the actual "true" center of the image. You can choose which
of these two centers you want to use with the ``use_true_center`` parameters, which
defaults to False. You can also arbitrarily offset the profile from the image center with
the ``offset`` parameter to handle any aribtrary offset you want from the chosen center.
(Typically, one would use only one of ``center`` or ``offset`` but it is permissible to use
both.)
Setting the overall normalization:
Normally, the flux of the object should be equal to the sum of all the pixel values in the
image, less some small amount of flux that may fall off the edge of the image (assuming you
don't use ``method='sb'``). However, you may optionally set a ``gain`` value, which
converts between photons and ADU (so-called analog-to-digital units), the units of the
pixel values in real images. Normally, the gain of a CCD is in electrons/ADU, but in
GalSim, we fold the quantum efficiency into the gain as well, so the units are photons/ADU.
Another caveat is that, technically, flux is really in units of photons/cm^2/s, not photons.
So if you want, you can keep track of this properly and provide an ``area`` and ``exptime``
here. This detail is more important with chromatic objects where the `SED` is typically
given in erg/cm^2/s/nm, so the exposure time and area are important details. With achromatic
objects however, it is often more convenient to ignore these details and just consider the
flux to be the total number of photons for this exposure, in which case, you would leave the
area and exptime parameters at their default value of 1.
On return, the image will have an attribute ``added_flux``, which will be set to the total
flux added to the image. This may be useful as a sanity check that you have provided a
large enough image to catch most of the flux. For example::
>>> obj.drawImage(image)
>>> assert image.added_flux > 0.99 * obj.flux
The appropriate threshold will depend on your particular application, including what kind
of profile the object has, how big your image is relative to the size of your object,
whether you are keeping ``poisson_flux=True``, etc.
The following code snippet illustrates how ``gain``, ``exptime``, ``area``, and ``method``
can all influence the relationship between the ``flux`` attribute of a `GSObject` and
both the pixel values and ``.added_flux`` attribute of an `Image` drawn with
``drawImage()``::
>>> obj = galsim.Gaussian(fwhm=1)
>>> obj.flux
1.0
>>> im = obj.drawImage()
>>> im.added_flux
0.9999630988657515
>>> im.array.sum()
0.99996305
>>> im = obj.drawImage(exptime=10, area=10)
>>> im.added_flux
0.9999630988657525
>>> im.array.sum()
99.996315
>>> im = obj.drawImage(exptime=10, area=10, method='sb', scale=0.5, nx=10, ny=10)
>>> im.added_flux
0.9999973790505298
>>> im.array.sum()
399.9989
>>> im = obj.drawImage(exptime=10, area=10, gain=2)
>>> im.added_flux
0.9999630988657525
>>> im.array.sum()
49.998158
Using a non-trivial sensor:
Normally the sensor is modeled as an array of pixels where any photon that hits a given
pixel is accumulated into that pixel. The final pixel value then just reflects the total
number of pixels that hit each sensor. However, real sensors do not (quite) work this way.
In real CCDs, the photons travel some distance into the silicon before converting to
electrons. Then the electrons diffuse laterally some amount as they are pulled by the
electric field toward the substrate. Finally, previous electrons that have already been
deposited will repel subsequent electrons, both slowing down their descent, leading to
more diffusion, and pushing them laterally toward neighboring pixels, which is called
the brighter-fatter effect.
Users interested in modeling this kind of effect can supply a ``sensor`` object to use
for the accumulation step. See `SiliconSensor` for a class that models silicon-based CCD
sensors.
Some related effects may need to be done to the photons at the surface layer before being
passed into the sensor object. For instance, the photons may need to be given appropriate
incidence angles according to the optics of the telescope (since this matters for where the
photons are converted to electrons). You may also need to give the photons wavelengths
according to the `SED` of the object. Such steps are specified in a ``photon_ops``
parameter, which should be a list of any such operations you wish to perform on the photon
array before passing them to the sensor. See `FRatioAngles` and `WavelengthSampler` for
two examples of such photon operators.
Since the sensor deals with photons, it is most natural to use this feature in conjunction
with photon shooting (``method='phot'``). However, it is allowed with FFT methods too.
But there is a caveat one should be aware of in this case. The FFT drawing is used to
produce an intermediate image, which is then converted to a `PhotonArray` using the
factory function `PhotonArray.makeFromImage`. This assigns photon positions randomly
within each pixel where they were drawn, which isn't always a particularly good
approximation.
To improve this behavior, the intermediate image is drawn with smaller pixels than the
target image, so the photons are given positions closer to their true locations. The
amount of subsampling is controlled by the ``n_subsample`` parameter, which defaults to 3.
Larger values will be more accurate at the expense of larger FFTs (i.e. slower and using
more memory).
Parameters:
image: If provided, this will be the image on which to draw the profile.
If ``image`` is None, then an automatically-sized `Image` will be
created. If ``image`` is given, but its bounds are undefined (e.g. if
it was constructed with ``image = galsim.Image()``), then it will be
resized appropriately based on the profile's size [default: None].
nx: If provided and ``image`` is None, use to set the x-direction size of
the image. Must be accompanied by ``ny``.
ny: If provided and ``image`` is None, use to set the y-direction size of
the image. Must be accompanied by ``nx``.
bounds: If provided and ``image`` is None, use to set the bounds of the image.
scale: If provided, use this as the pixel scale for the image.
If ``scale`` is None and ``image`` is given, then take the provided
image's pixel scale.
If ``scale`` is None and ``image`` is None, then use the Nyquist scale.
If ``scale <= 0`` (regardless of ``image``), then use the Nyquist scale.
If ``scale > 0`` and ``image`` is given, then override ``image.scale``
with the value given as a keyword. [default: None]
wcs: If provided, use this as the wcs for the image (possibly overriding any
existing ``image.wcs``). At most one of ``scale`` or ``wcs`` may be
provided. [default: None]
dtype: The data type to use for an automatically constructed image. Only
valid if ``image`` is None. [default: None, which means to use
numpy.float32]
method: Which method to use for rendering the image. See discussion above
for the various options and what they do. [default: 'auto']
area: Collecting area of telescope in cm^2. [default: 1.]
exptime: Exposure time in s. [default: 1.]
gain: The number of photons per ADU ("analog to digital units", the units of
the numbers output from a CCD). [default: 1]
add_to_image: Whether to add flux to the existing image rather than clear out
anything in the image before drawing.
Note: This requires that ``image`` be provided and that it have defined
bounds. [default: False]
center: The position on the image at which to place the nominal center of the
object (usually the centroid, but not necessarily). [default: None]
use_true_center: If ``center`` is None, then the object is normally centered at the
true center of the image (using the property image.true_center).
If you would rather use the integer center (given by image.center),
set this to ``False``. [default: True]
offset: The offset in pixel coordinates at which to center the profile being
drawn relative to either ``center`` (if given) or the center of the
image (either the true center or integer center according to the
``use_true_center`` parameter). [default: None]
n_photons: If provided, the number of photons to use for photon shooting.
If not provided (i.e. ``n_photons = 0``), use as many photons as
necessary to result in an image with the correct Poisson shot
noise for the object's flux. For positive definite profiles, this
is equivalent to ``n_photons = flux``. However, some profiles need
more than this because some of the shot photons are negative
(usually due to interpolants).
[default: 0]
rng: If provided, a random number generator to use for photon shooting,
which may be any kind of `BaseDeviate` object. If ``rng`` is None, one
will be automatically created. [default: None]
max_extra_noise: If provided, the allowed extra noise in each pixel when photon
shooting. This is only relevant if ``n_photons=0``, so the number of
photons is being automatically calculated. In that case, if the image
noise is dominated by the sky background, then you can get away with
using fewer shot photons than the full ``n_photons = flux``.
Essentially each shot photon can have a ``flux > 1``, which increases
the noise in each pixel. The ``max_extra_noise`` parameter specifies
how much extra noise per pixel is allowed because of this approximation.
A typical value for this might be ``max_extra_noise = sky_level / 100``
where ``sky_level`` is the flux per pixel due to the sky. Note that
this uses a "variance" definition of noise, not a "sigma" definition.
[default: 0.]
poisson_flux: Whether to allow total object flux scaling to vary according to
Poisson statistics for ``n_photons`` samples when photon shooting.
[default: True, unless ``n_photons`` is given, in which case the default
is False]
sensor: An optional `Sensor` instance, which will be used to accumulate the
photons onto the image. [default: None]
photon_ops: A list of operators that can modify the photon array that will be
applied in order before accumulating the photons on the sensor.
[default: None]
n_subsample: The number of sub-pixels per final pixel to use for fft drawing when
using a sensor. The sensor step needs to know the sub-pixel positions
of the photons, which is lost in the fft method. So using smaller
pixels for the fft step keeps some of that information, making the
assumption of uniform flux per pixel less bad of an approximation.
[default: 3]
maxN: Sets the maximum number of photons that will be added to the image
at a time. (Memory requirements are proportional to this number.)
.. note::
Using this parameter will not necessarily produce identical
results as when not using it due to potentially different order
of various random number generations in either the photon_ops,
or the sensor, or (for method='fft') the conversion of the FFT
image to photons.
[default: None, which means no limit]
save_photons: If True, save the `PhotonArray` as ``image.photons``. Only valid if
method is 'phot' or sensor is not None. [default: False]
bandpass: This parameter is ignored, but is allowed to enable duck typing
eqivalence between this method and the ChromaticObject.drawImage
method. [default: None]
setup_only: Don't actually draw anything on the image. Just make sure the image
is set up correctly. This is used internally by GalSim, but there
may be cases where the user will want the same functionality.
[default: False]
Returns:
the drawn `Image`.
"""
if surface_ops is not None:
from .deprecated import depr
depr('surface_ops', 2.3, 'photon_ops')
photon_ops = surface_ops
# Check that image is sane
if image is not None and not isinstance(image, Image):
raise TypeError("image is not an Image instance", image)
# Make sure (gain, area, exptime) have valid values:
if not gain > 0.:
raise GalSimRangeError("Invalid gain <= 0.", gain, 0., None)
if not area > 0.:
raise GalSimRangeError("Invalid area <= 0.", area, 0., None)
if not exptime > 0.:
raise GalSimRangeError("Invalid exptime <= 0.", exptime, 0., None)
if method not in ('auto', 'fft', 'real_space', 'phot', 'no_pixel', 'sb'):
raise GalSimValueError("Invalid method name", method,
('auto', 'fft', 'real_space', 'phot', 'no_pixel', 'sb'))
# Check that the user isn't convolving by a Pixel already. This is almost always an error.
if method == 'auto' and isinstance(self, convolve.Convolution):
if any([ isinstance(obj, box.Pixel) for obj in self.obj_list ]):
galsim_warn(
"You called drawImage with ``method='auto'`` "
"for an object that includes convolution by a Pixel. "
"This is probably an error. Normally, you should let GalSim "
"handle the Pixel convolution for you. If you want to handle the Pixel "
"convolution yourself, you can use method=no_pixel. Or if you really meant "
"for your profile to include the Pixel and also have GalSim convolve by "
"an _additional_ Pixel, you can suppress this warning by using method=fft.")
# Some parameters are only relevant for method == 'phot'
if method != 'phot':
if n_photons != 0.:
raise GalSimIncompatibleValuesError(
"n_photons is only relevant for method='phot'",
method=method, sensor=sensor, n_photons=n_photons)
if max_extra_noise != 0.:
raise GalSimIncompatibleValuesError(
"max_extra_noise is only relevant for method='phot'",
method=method, sensor=sensor, max_extra_noise=max_extra_noise)
if poisson_flux is not None:
raise GalSimIncompatibleValuesError(
"poisson_flux is only relevant for method='phot'",
method=method, sensor=sensor, poisson_flux=poisson_flux)
# If using photon ops with fft, then need a sensor.
# Note: "if photon_ops" so photon_ops being [] or () is the same as None.
if photon_ops and method != 'phot' and sensor is None:
sensor = Sensor()
if photon_ops is None:
# Easier to just make it an empty tuple, rather than deal with None below.
photon_ops = ()
# Some parameters are only relevant for either phot or when using a sensor.
if method != 'phot' and sensor is None:
if rng is not None:
raise GalSimIncompatibleValuesError(
"rng is only relevant for method='phot' or when using a sensor",
method=method, sensor=sensor, rng=rng)
if maxN != None:
raise GalSimIncompatibleValuesError(
"maxN is only relevant for method='phot' or when using a sensor",
method=method, sensor=sensor, maxN=maxN)
if save_photons:
raise GalSimIncompatibleValuesError(
"save_photons is only valid for method='phot' or when using a sensor",
method=method, sensor=sensor, save_photons=save_photons)
else:
# If we want to save photons, it doesn't make sense to limit the number per shoot call.
if save_photons and maxN is not None:
raise GalSimIncompatibleValuesError(
"Setting maxN is incompatible with save_photons=True")
# Do any delayed computation needed by fft or real_space drawing.
if method != 'phot':
self._prepareDraw()
# Figure out what wcs we are going to use.
wcs = self._determine_wcs(scale, wcs, image)
# Make sure offset and center are PositionD, converting from other formats (tuple, array,..)
# Note: If None, offset is converted to PositionD(0,0), but center will remain None.
offset = self._parse_offset(offset)
center = self._parse_center(center)
# Determine the bounds of the new image for use below (if it can be known yet)
new_bounds = self._get_new_bounds(image, nx, ny, bounds, center)
# Get the local WCS, accounting for the offset correctly.
local_wcs = self._local_wcs(wcs, image, offset, center, use_true_center, new_bounds)
# Account for area and exptime.
flux_scale = area * exptime
# For surface brightness normalization, also scale by the pixel area.
if method == 'sb':
flux_scale /= local_wcs.pixelArea()
# Only do the gain here if not photon shooting, since need the number of photons to
# reflect that actual photons, not ADU.
if gain != 1 and method != 'phot' and sensor is None:
flux_scale /= gain
# Determine the offset, and possibly fix the centering for even-sized images
offset = self._adjust_offset(new_bounds, offset, center, use_true_center)
# Convert the profile in world coordinates to the profile in image coordinates:
prof = local_wcs.profileToImage(self, flux_ratio=flux_scale, offset=offset)
if offset != _PositionD(0,0):
local_wcs = local_wcs.shiftOrigin(offset)
# If necessary, convolve by the pixel
if method in ('auto', 'fft', 'real_space'):
if method == 'auto':
real_space = None
elif method == 'fft':
real_space = False
else:
real_space = True
prof_no_pixel = prof
prof = convolve.Convolve(prof, box.Pixel(scale=1.0, gsparams=self.gsparams),
real_space=real_space, gsparams=self.gsparams)
# Make sure image is setup correctly
image = prof._setup_image(image, nx, ny, bounds, add_to_image, dtype, center)
image.wcs = wcs
if setup_only:
image.added_flux = 0.
return image
# Making a view of the image lets us change the center without messing up the original.
imview = image._view()
imview._shift(-image.center) # equiv. to setCenter(0,0), but faster
imview.wcs = PixelScale(1.0)
orig_center = image.center # Save the original center to pass to sensor.accumulate
if method == 'phot':
added_photons, photons = prof.drawPhot(imview, gain, add_to_image,
n_photons, rng, max_extra_noise, poisson_flux,
sensor, photon_ops, maxN,
orig_center, local_wcs)
else:
# If not using phot, but doing sensor, then make a copy.
if sensor is not None:
if imview.dtype in (np.float32, np.float64):
dtype = None
else:
dtype = np.float64
draw_image = imview.real.subsample(n_subsample, n_subsample, dtype=dtype)
draw_image._shift(-draw_image.center) # eqiv. to setCenter(0,0)
if method in ('auto', 'fft', 'real_space'):
# Need to reconvolve by the new smaller pixel instead
pix = box.Pixel(scale=1.0/n_subsample, gsparams=self.gsparams)
prof = convolve.Convolve(prof_no_pixel, pix, real_space=real_space,
gsparams=self.gsparams)
elif n_subsample != 1:
# We can't just pull off the pixel-free version, so we need to deconvolve
# by the original pixel and reconvolve by the smaller one.
dec = convolve.Deconvolve(box.Pixel(scale=1.0, gsparams=self.gsparams))
pix = box.Pixel(scale=1.0/n_subsample, gsparams=self.gsparams)
prof = convolve.Convolve(prof, dec, pix, gsparams=self.gsparams)
add = False
if not add_to_image: imview.setZero()
else:
draw_image = imview
add = add_to_image
if prof.is_analytic_x:
added_photons = prof.drawReal(draw_image, add)
else:
added_photons = prof.drawFFT(draw_image, add)
if sensor is not None:
if maxN is not None and maxN < added_photons:
niter = int(np.ceil(added_photons / maxN))
draw_image /= niter
else:
niter = 1
added_photons = 0
resume = False
for it in range(niter):
photons = pa.PhotonArray.makeFromImage(draw_image, rng=rng)
for op in photon_ops:
op.applyTo(photons, local_wcs, rng)
if imview.dtype in (np.float32, np.float64):
added_photons += sensor.accumulate(photons, imview, orig_center,
resume=resume)
resume = True # Resume from this point if there are any further iterations.
else:
# Need a temporary
im1 = ImageD(bounds=imview.bounds)
added_photons += sensor.accumulate(photons, im1, orig_center)
imview += im1
image.added_flux = added_photons / flux_scale
if save_photons:
image.photons = photons
return image
[docs] def drawReal(self, image, add_to_image=False):
"""
Draw this profile into an `Image` by direct evaluation at the location of each pixel.
This is usually called from the `drawImage` function, rather than called directly by the
user. In particular, the input image must be already set up with defined bounds. The
profile will be drawn centered on whatever pixel corresponds to (0,0) with the given
bounds, not the image center (unlike `drawImage`). The image also must have a `PixelScale`
wcs. The profile being drawn should have already been converted to image coordinates via::
>>> image_profile = original_wcs.toImage(original_profile)
Note that the image produced by ``drawReal`` represents the profile sampled at the center
of each pixel and then multiplied by the pixel area. That is, the profile is NOT
integrated over the area of the pixel. This is equivalent to method='no_pixel' in
`drawImage`. If you want to render a profile integrated over the pixel, you can convolve
with a `Pixel` first and draw that.
Parameters:
image: The `Image` onto which to place the flux. [required]
add_to_image: Whether to add flux to the existing image rather than clear out
anything in the image before drawing. [default: False]
Returns:
The total flux drawn inside the image bounds.
"""
if image.wcs is None or not image.wcs._isPixelScale:
raise GalSimValueError("drawReal requires an image with a PixelScale wcs", image)
if image.dtype in (np.float64, np.float32) and not add_to_image and image.iscontiguous:
self._drawReal(image)
return image.array.sum(dtype=float)
else:
# Need a temporary
if image.dtype in (np.complex128, np.int32, np.uint32):
im1 = ImageD(bounds=image.bounds, scale=image.scale)
else:
im1 = ImageF(bounds=image.bounds, scale=image.scale)
self._drawReal(im1)
if add_to_image:
image += im1
else:
image._copyFrom(im1)
return im1.array.sum(dtype=float)
[docs] def _drawReal(self, image, jac=None, offset=(0.,0.), flux_scaling=1.):
"""A version of `drawReal` without the sanity checks or some options.
This is nearly equivalent to the regular ``drawReal(image, add_to_image=False)``, but
the image's dtype must be either float32 or float64, and it must have a c_contiguous array
(``image.iscontiguous`` must be True).
"""
raise NotImplementedError("%s does not implement drawReal"%self.__class__.__name__)
[docs] def getGoodImageSize(self, pixel_scale):
"""Return a good size to use for drawing this profile.
The size will be large enough to cover most of the flux of the object. Specifically,
at least (1-gsparams.folding_threshold) (i.e. 99.5% by default) of the flux should fall
in the image.
Also, the returned size is always an even number, which is usually desired in practice.
Of course, if you prefer an odd-sized image, you can add 1 to the result.
Parameters:
pixel_scale: The desired pixel scale of the image to be built.
Returns:
N, a good (linear) size of an image on which to draw this object.
"""
# Start with a good size from stepk and the pixel scale
Nd = 2. * math.pi / (pixel_scale * self.stepk)
# Make it an integer
# (Some slop to keep from getting extra pixels due to roundoff errors in calculations.)
N = int(math.ceil(Nd*(1.-1.e-12)))
# Round up to an even value
N = 2 * ((N+1) // 2)
return N
[docs] def drawFFT_makeKImage(self, image):
"""
This is a helper routine for drawFFT that just makes the (blank) k-space image
onto which the profile will be drawn. This can be useful if you want to break
up the calculation into parts for extra efficiency. E.g. save the k-space image of
the PSF so drawing many models of the galaxy with the given PSF profile can avoid
drawing the PSF each time.
Parameters:
image: The `Image` onto which to place the flux.
Returns:
(kimage, wrap_size), where wrap_size is either the size of kimage or smaller if
the result should be wrapped before doing the inverse fft.
"""
# Start with what this profile thinks a good size would be given the image's pixel scale.
N = self.getGoodImageSize(image.scale)
# We must make something big enough to cover the target image size:
image_N = max(np.max(np.abs((image.bounds._getinitargs()))) * 2,
np.max(image.bounds.numpyShape()))
N = max(N, image_N)
# Round up to a good size for making FFTs:
N = image.good_fft_size(N)
# Make sure we hit the minimum size specified in the gsparams.
N = max(N, self.gsparams.minimum_fft_size)
dk = 2.*np.pi / (N * image.scale)
maxk = self.maxk
if N*dk/2 > maxk:
Nk = N
else:
# There will be aliasing. Make a larger image and then wrap it.
Nk = int(np.ceil(maxk/dk)) * 2
if Nk > self.gsparams.maximum_fft_size:
raise GalSimFFTSizeError("drawFFT requires an FFT that is too large.", Nk)
bounds = _BoundsI(0,Nk//2,-Nk//2,Nk//2)
if image.dtype in (np.complex128, np.float64, np.int32, np.uint32):
kimage = ImageCD(bounds=bounds, scale=dk)
else:
kimage = ImageCF(bounds=bounds, scale=dk)
return kimage, N
[docs] def drawFFT_finish(self, image, kimage, wrap_size, add_to_image):
"""
This is a helper routine for drawFFT that finishes the calculation, based on the
drawn k-space image.
It applies the Fourier transform to ``kimage`` and adds the result to ``image``.
Parameters:
image: The `Image` onto which to place the flux.
kimage: The k-space `Image` where the object was drawn.
wrap_size: The size of the region to wrap kimage, which must be either the same
size as kimage or smaller.
add_to_image: Whether to add flux to the existing image rather than clear out
anything in the image before drawing.
Returns:
The total flux drawn inside the image bounds.
"""
# Wrap the full image to the size we want for the FT.
# Even if N == Nk, this is useful to make this portion properly Hermitian in the
# N/2 column and N/2 row.
bwrap = _BoundsI(0, wrap_size//2, -wrap_size//2, wrap_size//2-1)
kimage_wrap = kimage._wrap(bwrap, True, False)
# Perform the fourier transform.
breal = _BoundsI(-wrap_size//2, wrap_size//2+1, -wrap_size//2, wrap_size//2-1)
real_image = Image(breal, dtype=float)
with convert_cpp_errors():
_galsim.irfft(kimage_wrap._image, real_image._image, True, True)
# Add (a portion of) this to the original image.
temp = real_image.subImage(image.bounds)
if add_to_image:
image += temp
else:
image._copyFrom(temp)
added_photons = temp.array.sum(dtype=float)
return added_photons
[docs] def drawFFT(self, image, add_to_image=False):
"""
Draw this profile into an `Image` by computing the k-space image and performing an FFT.
This is usually called from the `drawImage` function, rather than called directly by the
user. In particular, the input image must be already set up with defined bounds. The
profile will be drawn centered on whatever pixel corresponds to (0,0) with the given
bounds, not the image center (unlike `drawImage`). The image also must have a `PixelScale`
wcs. The profile being drawn should have already been converted to image coordinates via::
>>> image_profile = original_wcs.toImage(original_profile)
Note that the `Image` produced by drawFFT represents the profile sampled at the center
of each pixel and then multiplied by the pixel area. That is, the profile is NOT
integrated over the area of the pixel. This is equivalent to method='no_pixel' in
`drawImage`. If you want to render a profile integrated over the pixel, you can convolve
with a `Pixel` first and draw that.
Parameters:
image: The `Image` onto which to place the flux. [required]
add_to_image: Whether to add flux to the existing image rather than clear out
anything in the image before drawing. [default: False]
Returns:
The total flux drawn inside the image bounds.
"""
if image.wcs is None or not image.wcs._isPixelScale:
raise GalSimValueError("drawPhot requires an image with a PixelScale wcs", image)
kimage, wrap_size = self.drawFFT_makeKImage(image)
self._drawKImage(kimage)
return self.drawFFT_finish(image, kimage, wrap_size, add_to_image)
[docs] def _calculate_nphotons(self, n_photons, poisson_flux, max_extra_noise, rng):
"""Calculate how many photons to shoot and what flux_ratio (called g) each one should
have in order to produce an image with the right S/N and total flux.
This routine is normally called by `drawPhot`.
Returns:
n_photons, g
"""
# For profiles that are positive definite, then N = flux. Easy.
#
# However, some profiles shoot some of their photons with negative flux. This means that
# we need a few more photons to get the right S/N = sqrt(flux). Take eta to be the
# fraction of shot photons that have negative flux.
#
# S^2 = (N+ - N-)^2 = (N+ + N- - 2N-)^2 = (Ntot - 2N-)^2 = Ntot^2(1 - 2 eta)^2
# N^2 = Var(S) = (N+ + N-) = Ntot
#
# So flux = (S/N)^2 = Ntot (1-2eta)^2
# Ntot = flux / (1-2eta)^2
#
# However, if each photon has a flux of 1, then S = (1-2eta) Ntot = flux / (1-2eta).
# So in fact, each photon needs to carry a flux of g = 1-2eta to get the right
# total flux.
#
# That's all the easy case. The trickier case is when we are sky-background dominated.
# Then we can usually get away with fewer shot photons than the above. In particular,
# if the noise from the photon shooting is much less than the sky noise, then we can
# use fewer shot photons and essentially have each photon have a flux > 1. This is ok
# as long as the additional noise due to this approximation is "much less than" the
# noise we'll be adding to the image for the sky noise.
#
# Let's still have Ntot photons, but now each with a flux of g. And let's look at the
# noise we get in the brightest pixel that has a nominal total flux of Imax.
#
# The number of photons hitting this pixel will be Imax/flux * Ntot.
# The variance of this number is the same thing (Poisson counting).
# So the noise in that pixel is:
#
# N^2 = Imax/flux * Ntot * g^2
#
# And the signal in that pixel will be:
#
# S = Imax/flux * (N+ - N-) * g which has to equal Imax, so
# g = flux / Ntot(1-2eta)
# N^2 = Imax/Ntot * flux / (1-2eta)^2
#
# As expected, we see that lowering Ntot will increase the noise in that (and every
# other) pixel.
# The input max_extra_noise parameter is the maximum value of spurious noise we want
# to allow.
#
# So setting N^2 = Imax + nu, we get
#
# Ntot = flux / (1-2eta)^2 / (1 + nu/Imax)
# g = (1 - 2eta) * (1 + nu/Imax)
#
# Returns the total flux placed inside the image bounds by photon shooting.
#
flux = self.flux
if flux == 0.0:
return 0, 1.0
# The _flux_per_photon property is (1-2eta)
# This factor will already be accounted for by the shoot function, so don't include
# that as part of our scaling here. There may be other adjustments though, so g=1 here.
eta_factor = self._flux_per_photon
mod_flux = flux / (eta_factor * eta_factor)
g = 1.
# If requested, let the target flux value vary as a Poisson deviate
if poisson_flux:
# If we have both positive and negative photons, then the mix of these
# already gives us some variation in the flux value from the variance
# of how many are positive and how many are negative.
# The number of negative photons varies as a binomial distribution.
# <F-> = eta * Ntot * g
# <F+> = (1-eta) * Ntot * g
# <F+ - F-> = (1-2eta) * Ntot * g = flux
# Var(F-) = eta * (1-eta) * Ntot * g^2
# F+ = Ntot * g - F- is not an independent variable, so
# Var(F+ - F-) = Var(Ntot*g - 2*F-)
# = 4 * Var(F-)
# = 4 * eta * (1-eta) * Ntot * g^2
# = 4 * eta * (1-eta) * flux
# We want the variance to be equal to flux, so we need an extra:
# delta Var = (1 - 4*eta + 4*eta^2) * flux
# = (1-2eta)^2 * flux
absflux = abs(flux)
mean = eta_factor*eta_factor * absflux
pd = PoissonDeviate(rng, mean)
pd_val = pd() - mean + absflux
ratio = pd_val / absflux
g *= ratio
mod_flux *= ratio
if n_photons == 0.:
n_photons = abs(mod_flux)
if max_extra_noise > 0.:
gfactor = 1. + max_extra_noise / abs(self.max_sb)
n_photons /= gfactor
g *= gfactor
# Make n_photons an integer.
iN = int(n_photons + 0.5)
return iN, g
[docs] def makePhot(self, n_photons=0, rng=None, max_extra_noise=0., poisson_flux=None,
photon_ops=(), local_wcs=None, surface_ops=None):
"""
Make photons for a profile.
This is equivalent to drawPhot, except that the photons are not placed onto
an image. Instead, it just returns the PhotonArray.
.. note::
The (x,y) positions returned are in the same units as the distance units
of the GSObject being rendered. If you want (x,y) in pixel coordinates, you
should call this function for the profile in image coordinates::
>>> photons = image.wcs.toImage(obj).makePhot()
Or if you just want a simple pixel scale conversion from sky coordinates to image
coordinates, you can instead do
>>> photons = obj.makePhot()
>>> photons.scaleXY(1./pixel_scale)
Parameters:
n_photons: If provided, the number of photons to use for photon shooting.
If not provided (i.e. ``n_photons = 0``), use as many photons as
necessary to result in an image with the correct Poisson shot
noise for the object's flux. For positive definite profiles, this
is equivalent to ``n_photons = flux``. However, some profiles need
more than this because some of the shot photons are negative
(usually due to interpolants). [default: 0]
rng: If provided, a random number generator to use for photon shooting,
which may be any kind of `BaseDeviate` object. If ``rng`` is None, one
will be automatically created, using the time as a seed.
[default: None]
max_extra_noise: If provided, the allowed extra noise in each pixel when photon
shooting. This is only relevant if ``n_photons=0``, so the number of
photons is being automatically calculated. In that case, if the image
noise is dominated by the sky background, then you can get away with
using fewer shot photons than the full ``n_photons = flux``.
Essentially each shot photon can have a ``flux > 1``, which increases
the noise in each pixel. The ``max_extra_noise`` parameter specifies
how much extra noise per pixel is allowed because of this approximation.
A typical value for this might be ``max_extra_noise = sky_level / 100``
where ``sky_level`` is the flux per pixel due to the sky. Note that
this uses a "variance" definition of noise, not a "sigma" definition.
[default: 0.]
poisson_flux: Whether to allow total object flux scaling to vary according to
Poisson statistics for ``n_photons`` samples when photon shooting.
[default: True, unless ``n_photons`` is given, in which case the default
is False]
photon_ops: A list of operators that can modify the photon array that will be
applied in order before accumulating the photons on the sensor.
[default: ()]
local_wcs: The local wcs in the original image. [default: None]
Returns:
- a `PhotonArray` with the data about the photons.
"""
if surface_ops is not None:
from .deprecated import depr
depr('surface_ops', 2.3, 'photon_ops')
photon_ops = surface_ops
# Make sure the type of n_photons is correct and has a valid value:
if not n_photons >= 0.:
raise GalSimRangeError("Invalid n_photons < 0.", n_photons, 0., None)
if poisson_flux is None:
# If n_photons is given, poisson_flux = False
poisson_flux = (n_photons == 0.)
# Check that either n_photons is set to something or flux is set to something
if n_photons == 0. and self.flux == 1.:
galsim_warn(
"Warning: drawImage for object with flux == 1, area == 1, and "
"exptime == 1, but n_photons == 0. This will only shoot a single photon.")
Ntot, g = self._calculate_nphotons(n_photons, poisson_flux, max_extra_noise, rng)
try:
photons = self.shoot(Ntot, rng)
except (GalSimError, NotImplementedError) as e:
raise GalSimNotImplementedError(
"Unable to draw this GSObject with photon shooting. Perhaps it "
"is a Deconvolve or is a compound including one or more "
"Deconvolve objects.\nOriginal error: %r"%(e))
if g != 1.:
photons.scaleFlux(g)
for op in photon_ops:
op.applyTo(photons, local_wcs, rng)
return photons
[docs] def drawPhot(self, image, gain=1., add_to_image=False,
n_photons=0, rng=None, max_extra_noise=0., poisson_flux=None,
sensor=None, photon_ops=(), maxN=None, orig_center=_PositionI(0,0),
local_wcs=None, surface_ops=None):
"""
Draw this profile into an `Image` by shooting photons.
This is usually called from the `drawImage` function, rather than called directly by the
user. In particular, the input image must be already set up with defined bounds. The
profile will be drawn centered on whatever pixel corresponds to (0,0) with the given
bounds, not the image center (unlike `drawImage`). The image also must have a `PixelScale`
wcs. The profile being drawn should have already been converted to image coordinates via::
>>> image_profile = original_wcs.toImage(original_profile)
Note that the image produced by `drawPhot` represents the profile integrated over the
area of each pixel. This is equivalent to convolving the profile by a square `Pixel`
profile and sampling the value at the center of each pixel, although this happens
automatically by the shooting algorithm, so you do not need to manually convolve by
a `Pixel` as you would for `drawReal` or `drawFFT`.
Parameters:
image: The `Image` onto which to place the flux. [required]
gain: The number of photons per ADU ("analog to digital units", the units of
the numbers output from a CCD). [default: 1.]
add_to_image: Whether to add to the existing images rather than clear out
anything in the image before drawing. [default: False]
n_photons: If provided, the number of photons to use for photon shooting.
If not provided (i.e. ``n_photons = 0``), use as many photons as
necessary to result in an image with the correct Poisson shot
noise for the object's flux. For positive definite profiles, this
is equivalent to ``n_photons = flux``. However, some profiles need
more than this because some of the shot photons are negative
(usually due to interpolants). [default: 0]
rng: If provided, a random number generator to use for photon shooting,
which may be any kind of `BaseDeviate` object. If ``rng`` is None, one
will be automatically created, using the time as a seed.
[default: None]
max_extra_noise: If provided, the allowed extra noise in each pixel when photon
shooting. This is only relevant if ``n_photons=0``, so the number of
photons is being automatically calculated. In that case, if the image
noise is dominated by the sky background, then you can get away with
using fewer shot photons than the full ``n_photons = flux``.
Essentially each shot photon can have a ``flux > 1``, which increases
the noise in each pixel. The ``max_extra_noise`` parameter specifies
how much extra noise per pixel is allowed because of this approximation.
A typical value for this might be ``max_extra_noise = sky_level / 100``
where ``sky_level`` is the flux per pixel due to the sky. Note that
this uses a "variance" definition of noise, not a "sigma" definition.
[default: 0.]
poisson_flux: Whether to allow total object flux scaling to vary according to
Poisson statistics for ``n_photons`` samples when photon shooting.
[default: True, unless ``n_photons`` is given, in which case the default
is False]
sensor: An optional `Sensor` instance, which will be used to accumulate the
photons onto the image. [default: None]
photon_ops: A list of operators that can modify the photon array that will be
applied in order before accumulating the photons on the sensor.
[default: ()]
maxN: Sets the maximum number of photons that will be added to the image
at a time. (Memory requirements are proportional to this number.)
[default: None, which means no limit]
orig_center: The position of the image center in the original image coordinates.
[default: (0,0)]
local_wcs: The local wcs in the original image. [default: None]
Returns:
(added_flux, photons) where:
- added_flux is the total flux of photons that landed inside the image bounds, and
- photons is the `PhotonArray` that was applied to the image.
"""
if surface_ops is not None:
from .deprecated import depr
depr('surface_ops', 2.3, 'photon_ops')
photon_ops = surface_ops
# Make sure the type of n_photons is correct and has a valid value:
if not n_photons >= 0.:
raise GalSimRangeError("Invalid n_photons < 0.", n_photons, 0., None)
if poisson_flux is None:
# If n_photons is given, poisson_flux = False
poisson_flux = (n_photons == 0.)
# Check that either n_photons is set to something or flux is set to something
if n_photons == 0. and self.flux == 1.:
galsim_warn(
"Warning: drawImage for object with flux == 1, area == 1, and "
"exptime == 1, but n_photons == 0. This will only shoot a single photon.")
# Make sure the image is set up to have unit pixel scale and centered at 0,0.
if image.wcs is None or not image.wcs._isPixelScale:
raise GalSimValueError("drawPhot requires an image with a PixelScale wcs", image)
if sensor is None:
sensor = Sensor()
elif not isinstance(sensor, Sensor):
raise TypeError("The sensor provided is not a Sensor instance")
Ntot, g = self._calculate_nphotons(n_photons, poisson_flux, max_extra_noise, rng)
if gain != 1.:
g /= gain
# total flux falling inside image bounds, this will be returned on exit.
added_flux = 0.
if maxN is None:
maxN = Ntot
if not add_to_image: image.setZero()
# Nleft is the number of photons remaining to shoot.
Nleft = Ntot
photons = None # Just in case Nleft is already 0.
resume = False
while Nleft > 0:
# Shoot at most maxN at a time
thisN = min(maxN, Nleft)
try:
photons = self.shoot(thisN, rng)
except (GalSimError, NotImplementedError) as e:
raise GalSimNotImplementedError(
"Unable to draw this GSObject with photon shooting. Perhaps it "
"is a Deconvolve or is a compound including one or more "
"Deconvolve objects.\nOriginal error: %r"%(e))
if g != 1. or thisN != Ntot:
photons.scaleFlux(g * thisN / Ntot)
if image.scale != 1.:
photons.scaleXY(1./image.scale) # Convert x,y to image coords if necessary
for op in photon_ops:
op.applyTo(photons, local_wcs, rng)
if image.dtype in (np.float32, np.float64):
added_flux += sensor.accumulate(photons, image, orig_center, resume=resume)
resume = True # Resume from this point if there are any further iterations.
else:
# Need a temporary
im1 = ImageD(bounds=image.bounds)
added_flux += sensor.accumulate(photons, im1, orig_center)
image += im1
Nleft -= thisN
return added_flux, photons
[docs] def shoot(self, n_photons, rng=None):
"""Shoot photons into a `PhotonArray`.
Parameters:
n_photons: The number of photons to use for photon shooting.
rng: If provided, a random number generator to use for photon shooting,
which may be any kind of `BaseDeviate` object. If ``rng`` is None, one
will be automatically created, using the time as a seed.
[default: None]
Returns:
A `PhotonArray`.
"""
photons = pa.PhotonArray(n_photons)
if n_photons == 0:
# It's ok to shoot 0, but downstream can have problems with it, so just stop now.
return photons
if rng is None:
rng = BaseDeviate()
self._shoot(photons, rng)
return photons
[docs] def _shoot(self, photons, rng):
"""Shoot photons into the given `PhotonArray`.
This is the backend implementation of `shoot` once the `PhotonArray` has been constructed.
Parameters:
photons: A `PhotonArray` instance into which the photons should be placed.
rng: A `BaseDeviate` instance to use for the photon shooting,
"""
raise NotImplementedError("%s does not implement shoot"%self.__class__.__name__)
[docs] def applyTo(self, photon_array, local_wcs=None, rng=None):
"""Apply this surface brightness profile as a convolution to an existing photon array.
This method allows a GSObject to duck type as a PhotonOp, so one can apply a PSF
in a photon_ops list.
Parameters:
photon_array: A `PhotonArray` to apply the operator to.
local_wcs: A `LocalWCS` instance defining the local WCS for the current photon
bundle in case the operator needs this information. [default: None]
rng: A random number generator to use to effect the convolution.
[default: None]
"""
p1 = pa.PhotonArray(len(photon_array))
if photon_array.hasAllocatedWavelengths():
p1._wave = photon_array._wave
if photon_array.hasAllocatedPupil():
p1._pupil_u = photon_array._pupil_u
p1._pupil_v = photon_array._pupil_v
if photon_array.hasAllocatedTimes():
p1._time = photon_array._time
obj = local_wcs.toImage(self) if local_wcs is not None else self
obj._shoot(p1, rng)
photon_array.convolve(p1, rng)
[docs] def drawKImage(self, image=None, nx=None, ny=None, bounds=None, scale=None,
add_to_image=False, recenter=True, bandpass=None, setup_only=False):
"""Draws the k-space (complex) `Image` of the object, with bounds optionally set by input
`Image` instance.
Normalization is always such that image(0,0) = flux. Unlike the real-space `drawImage`
function, the (0,0) point will always be one of the actual pixel values. For even-sized
images, it will be 1/2 pixel above and to the right of the true center of the image.
Another difference from `drawImage` is that a wcs other than a simple pixel scale is not
allowed. There is no ``wcs`` parameter here, and if the images have a non-trivial wcs (and
you don't override it with the ``scale`` parameter), a TypeError will be raised.
Also, there is no convolution by a pixel. This is just a direct image of the Fourier
transform of the surface brightness profile.
Parameters:
image: If provided, this will be the `Image` onto which to draw the k-space
image. If ``image`` is None, then an automatically-sized image will be
created. If ``image`` is given, but its bounds are undefined, then it
will be resized appropriately based on the profile's size.
[default: None]
nx: If provided and ``image`` is None, use to set the x-direction size of
the image. Must be accompanied by ``ny``.
ny: If provided and ``image`` is None, use to set the y-direction size of
the image. Must be accompanied by ``nx``.
bounds: If provided and ``image`` is None, use to set the bounds of the image.
scale: If provided, use this as the pixel scale, dk, for the images.
If ``scale`` is None and ``image`` is given, then take the provided
images' pixel scale (which must be equal).
If ``scale`` is None and ``image`` is None, then use the Nyquist scale.
If ``scale <= 0`` (regardless of ``image``), then use the Nyquist scale.
[default: None]
add_to_image: Whether to add to the existing images rather than clear out
anything in the image before drawing.
Note: This requires that ``image`` be provided and that it has defined
bounds. [default: False]
recenter: Whether to recenter the image to put k = 0 at the center (True) or to
trust the provided bounds (False). [default: True]
bandpass: This parameter is ignored, but is allowed to enable duck typing
eqivalence between this method and the ChromaticObject.drawImage
method. [default: None]
setup_only: Don't actually draw anything on the image. Just make sure the image
is set up correctly. This is used internally by GalSim, but there
may be cases where the user will want the same functionality.
[default: False]
Returns:
an `Image` instance (created if necessary)
"""
# Make sure provided image is complex
if image is not None:
if not isinstance(image, Image):
raise TypeError("Provided image must be galsim.Image", image)
if not image.iscomplex:
raise GalSimValueError("Provided image must be complex", image)
# Possibly get the scale from image.
if image is not None and scale is None:
# Grab the scale to use from the image.
# This will raise a TypeError if image.wcs is not a PixelScale
scale = image.scale
# The input scale (via scale or image.scale) is really a dk value, so call it that for
# clarity here, since we also need the real-space pixel scale, which we will call dx.
if scale is None or not scale > 0:
dk = self.stepk
else:
dk = float(scale)
if image is not None and image.bounds.isDefined():
dx = np.pi/( max(image.array.shape) // 2 * dk )
elif scale is None or not scale > 0:
dx = self.nyquist_scale
else:
# Then dk = scale, which implies that we need to have dx smaller than nyquist_scale
# by a factor of (dk/stepk)
dx = self.nyquist_scale * dk / self.stepk
# If the profile needs to be constructed from scratch, the _setup_image function will
# do that, but only if the profile is in image coordinates for the real space image.
# So make that profile.
if image is None or not image.bounds.isDefined():
real_prof = PixelScale(dx).profileToImage(self)
dtype = np.complex128 if image is None else image.dtype
image = real_prof._setup_image(image, nx, ny, bounds, add_to_image, dtype,
center=None, odd=True)
else:
# Do some checks that setup_image would have done for us
if bounds is not None:
raise GalSimIncompatibleValuesError(
"Cannot provide bounds if image is provided", bounds=bounds, image=image)
if nx is not None or ny is not None:
raise GalSimIncompatibleValuesError(
"Cannot provide nx,ny if image is provided", nx=nx, ny=ny, image=image)
# Can't both recenter a provided image and add to it.
if recenter and image.center != _PositionI(0,0) and add_to_image:
raise GalSimIncompatibleValuesError(
"Cannot use add_to_image=True unless image is centered at (0,0) or recenter=False",
recenter=recenter, image=image, add_to_image=add_to_image)
# Set the center to 0,0 if appropriate
if recenter and image.center != _PositionI(0,0):
image._shift(-image.center)
# Set the wcs of the images to use the dk scale size
image.scale = dk
if setup_only:
return image
if not add_to_image and image.iscontiguous:
self._drawKImage(image)
else:
im2 = Image(bounds=image.bounds, dtype=image.dtype, scale=image.scale)
self._drawKImage(im2)
image += im2
return image
[docs] def _drawKImage(self, image, jac=None): # pragma: no cover (all our classes override this)
"""A version of `drawKImage` without the sanity checks or some options.
Equivalent to ``drawKImage(image, add_to_image=False, recenter=False, add_to_image=False)``,
but without the option to create the image automatically.
The input image must be provided as a complex `Image` instance (dtype=complex64 or
complex128), and the bounds should be set up appropriately (e.g. with 0,0 in the center if
so desired). This corresponds to recenter=False for the normal `drawKImage`. And, it must
have a c_contiguous array (image.iscontiguous must be True).
Parameters:
image: The `Image` onto which to draw the k-space image. [required]
"""
raise NotImplementedError("%s does not implement drawKImage"%self.__class__.__name__)
# Derived classes should define the __eq__ function
def __ne__(self, other): return not self.__eq__(other)
# Put these at the bottom to avoid circular import errors.
from . import convolve
from . import sum as _sum
from . import transform
from . import box
from . import photon_array as pa
from . import sed