Arbitrary Profiles

If none of the above classes seem appropriate, it is possible to define any arbitrary profile by an image and then interpolating between the pixel positions to define the surface brightness at any arbitrary location. Similarly, one can define the image in Fourier space, or use shapelets, which are a different complete basis set, instead of using pixels.

Interpolated Images

class galsim.InterpolatedImage(image, x_interpolant=None, k_interpolant=None, normalization='flux', scale=None, wcs=None, flux=None, pad_factor=4.0, noise_pad_size=0, noise_pad=0.0, rng=None, pad_image=None, calculate_stepk=True, calculate_maxk=True, use_cache=True, use_true_center=True, depixelize=False, offset=None, gsparams=None, _force_stepk=0.0, _force_maxk=0.0, hdu=None)[source]

Bases: GSObject

A class describing non-parametric profiles specified using an Image, which can be interpolated for the purpose of carrying out transformations.

The InterpolatedImage class is useful if you have a non-parametric description of an object as an Image, that you wish to manipulate / transform using GSObject methods such as GSObject.shear, GSObject.magnify, GSObject.shift, etc. Note that when convolving an InterpolatedImage, the use of real-space convolution is not recommended, since it is typically a great deal slower than Fourier-space convolution for this kind of object.

There are three options for determining the flux of the profile.

You can simply specify a flux value explicitly.
If you set normalization = 'flux', the flux will be taken as the sum of the pixel values in the input image. This corresponds to an image that was drawn with drawImage(method='no_pixel'). This is the default if flux is not given.
If you set normalization = 'sb', the pixel values are treated as samples of the surface brightness profile at each location. This corresponds to an image drawn with drawImage(method='sb'). The flux is then the sum of the pixels in the input image multiplied by the pixel area.

You can also use images that were drawn with one of the pixel-integrating methods (‘auto’, ‘phot’, ‘fft’, or ‘real_space’), or real images where nature has integrated over the pixel for you. However, the resulting profile will by default include a convolution by a Pixel (this is equivalent to integration over the pixel). This is often acceptable, and the resulting InterpolatedImage object can be convolved by other profiles as usual. One just needs to remember to draw the final convolved profile using method='no_pixel' to avoid including the pixel convolution a second time. In particular, one should not use it in conjunction with photon shooting, for the same reason.

However, if you want to try to remove the effect of the pixel and have the InterpolatedImage model the pre-pixelized profile, then you can set depixelize=True. This will call Image.depixelize on the image automatically to try to remove the effect of the pixelization. We recommend using a Lanczos interpolant with this option for best results. (Higher order tends to work better here.) This step can be rather slow and memory-demanding, so use this with caution. But if used, the resulting profile represents the true underlying profile, without the pixel convolution. It can therefore be rotated, sheared, etc. And when rendering, one should use the methods that do involve integration over the pixel: auto, phot, etc.

Warning

Input images that are undersampled and/or noisy may not necessarily work well with the depixelize=True option. Users should treat this option with some care and validate that the results are sufficiently accurate for your particular use case.

If the input Image has a scale or wcs associated with it, then there is no need to specify one as a parameter here. But if one is provided, that will override any scale or wcs that is native to the Image.

The user may optionally specify an interpolant, x_interpolant, for real-space manipulations (e.g., shearing, resampling). If none is specified, then by default, a Quintic interpolant is used. The user may also choose to specify two quantities that can affect the Fourier space convolution: the k-space interpolant (k_interpolant) and the amount of padding to include around the original images (pad_factor). The default values for x_interpolant, k_interpolant, and pad_factor were chosen based on the tests of branch #389 to reach good accuracy without being excessively slow. Users should be particularly wary about changing k_interpolant and pad_factor from the defaults without careful testing. The user is given complete freedom to choose interpolants and pad factors, and no warnings are raised when the code is modified to choose some combination that is known to give significant error. More details can be found in http://arxiv.org/abs/1401.2636, especially table 1, and in comment https://github.com/GalSim-developers/GalSim/issues/389#issuecomment-26166621 and the following comments.

The user can choose to pad the image with a noise profile if desired. To do so, specify the target size for the noise padding in noise_pad_size, and specify the kind of noise to use in noise_pad. The noise_pad option may be a Gaussian random noise of some variance, or a Gaussian but correlated noise field that is specified either as a BaseCorrelatedNoise instance, an Image (from which a correlated noise model is derived), or a string (interpreted as a filename of an image to use for deriving a CorrelatedNoise). The user can also pass in a random number generator to be used for noise generation. Finally, the user can pass in a pad_image for deterministic image padding.

By default, the InterpolatedImage recalculates the Fourier-space step and number of points to use for further manipulations, rather than using the most conservative possibility. For typical objects representing galaxies and PSFs this can easily make the difference between several seconds (conservative) and 0.04s (recalculated). However, the user can turn off this option, and may especially wish to do so when using images that do not contain a high S/N object - e.g., images of noise fields.

Example:

>>> interpolated_image = galsim.InterpolatedImage(
        image, x_interpolant=None, k_interpolant=None, normalization='flux', scale=None,
        wcs=None, flux=None, pad_factor=4., noise_pad_size=0, noise_pad=0., use_cache=True,
        pad_image=None, rng=None, calculate_stepk=True, calculate_maxk=True,
        use_true_center=True, offset=None, hdu=None)

Initializes interpolated_image as an InterpolatedImage instance.

For comparison of the case of padding with noise or zero when the image itself includes noise, compare im1 and im2 from the following code snippet (which can be executed from the examples/ directory):

>>> image = galsim.fits.read('data/147246.0_150.416558_1.998697_masknoise.fits')
>>> int_im1 = galsim.InterpolatedImage(image)
>>> int_im2 = galsim.InterpolatedImage(image, noise_pad='data/blankimg.fits')
>>> im1 = galsim.ImageF(1000,1000)
>>> im2 = galsim.ImageF(1000,1000)
>>> int_im1.drawImage(im1, method='no_pixel')
>>> int_im2.drawImage(im2, method='no_pixel')

Examination of these two images clearly shows how padding with a correlated noise field that is similar to the one in the real data leads to a more reasonable appearance for the result when re-drawn at a different size.

Parameters:

image – The Image from which to construct the object. This may be either an Image instance or a string indicating a fits file from which to read the image. In the latter case, the hdu kwarg can be used to specify a particular HDU in that file.
x_interpolant – Either an Interpolant instance or a string indicating which real-space interpolant should be used. Options are ‘nearest’, ‘sinc’, ‘linear’, ‘cubic’, ‘quintic’, or ‘lanczosN’ where N should be the integer order to use. [default: galsim.Quintic()]
k_interpolant – Either an Interpolant instance or a string indicating which k-space interpolant should be used. Options are ‘nearest’, ‘sinc’, ‘linear’, ‘cubic’, ‘quintic’, or ‘lanczosN’ where N should be the integer order to use. We strongly recommend leaving this parameter at its default value; see text above for details. [default: galsim.Quintic()]
normalization –
Two options for specifying the normalization of the input Image:
- ”flux” or “f” means that the sum of the pixels is normalized to be equal to the total flux.
- ”surface brightness” or “sb” means that the pixels sample the surface brightness distribution at each location.
This is overridden if you specify an explicit flux value. [default: “flux”]
scale – If provided, use this as the pixel scale for the Image; this will override the pixel scale stored by the provided Image, in any. If scale is None, then take the provided image’s pixel scale. [default: None]
wcs – If provided, use this as the wcs for the image. At most one of scale or wcs may be provided. [default: None]
flux – Optionally specify a total flux for the object, which overrides the implied flux normalization from the Image itself. [default: None]
pad_factor – Factor by which to pad the Image with zeros. We strongly recommend leaving this parameter at its default value; see text above for details. [default: 4]
noise_pad_size – If provided, the image will be padded out to this size (in arcsec) with the noise specified by noise_pad. This is important if you are planning to whiten the resulting image. You want to make sure that the noise-padded image is larger than the postage stamp onto which you are drawing this object. [default: None]
noise_pad –
Noise properties to use when padding the original image with noise. This can be specified in several ways:
1. as a float, which is interpreted as being a variance to use when padding with uncorrelated Gaussian noise;
2. as a galsim.BaseCorrelatedNoise, which contains information about the desired noise power spectrum - any random number generator passed to the rng keyword will take precedence over that carried in an input BaseCorrelatedNoise instance;
3. as an Image of a noise field, which is used to calculate the desired noise power spectrum; or
4. as a string which is interpreted as a filename containing an example noise field with the proper noise power spectrum (as an Image in the first HDU).
It is important to keep in mind that the calculation of the correlation function that is internally stored within a BaseCorrelatedNoise object is a non-negligible amount of overhead, so the recommended means of specifying a correlated noise field for padding are (b) or (d). In the case of (d), if the same file is used repeatedly, then the use_cache keyword (see below) can be used to prevent the need for repeated CorrelatedNoise initializations. [default: 0, i.e., pad with zeros]
use_cache – Specify whether to cache noise_pad read in from a file to save having to build a CorrelatedNoise object repeatedly from the same image. [default: True]
rng – If padding by noise, the user can optionally supply the random noise generator to use for drawing random numbers as rng (may be any kind of BaseDeviate object). Such a user-input random number generator takes precedence over any stored within a user-input BaseCorrelatedNoise instance (see the noise_pad parameter). If rng=None, one will be automatically created, using the time as a seed. [default: None]
pad_image –
Image to be used for deterministically padding the original image. This can be specified in two ways:
- as an Image; or
- as a string which is interpreted as a filename containing an image to use (in the first HDU).
The pad_image scale or wcs is ignored. It uses the same scale or wcs for both the image and the pad_image. The user should be careful to ensure that the image used for padding has roughly zero mean. The purpose of this keyword is to allow for a more flexible representation of some noise field around an object; if the user wishes to represent the sky level around an object, they should do that after they have drawn the final image instead. [default: None]
calculate_stepk – Specify whether to perform an internal determination of the extent of the object being represented by the InterpolatedImage; often this is useful in choosing an optimal value for the stepsize in the Fourier space lookup table. If you know a priori an appropriate maximum value for stepk, then you may also supply that here instead of a bool value, in which case the stepk value is still calculated, but will not go above the provided value. [default: True]
calculate_maxk – Specify whether to perform an internal determination of the highest spatial frequency needed to accurately render the object being represented by the InterpolatedImage; often this is useful in choosing an optimal value for the extent of the Fourier space lookup table. If you know a priori an appropriate maximum value for maxk, then you may also supply that here instead of a bool value, in which case the maxk value is still calculated, but will not go above the provided value. [default: True]
use_true_center – Similar to the same parameter in the GSObject.drawImage function, this sets whether to use the true center of the provided image as the center of the profile (if use_true_center=True) or the nominal center given by image.center (if use_true_center=False) [default: True]
depixelize – Whether to try to remove the effect of the pixelization. If this is True, then drawing this profile with method=’auto’ should render an image equivalent to the input image. If this is False (the default), then you would need to draw with method=’no_pixel’ to get an equivalent image. See discussion above. [default: False]
offset – The location in the input image to use as the center of the profile. This should be specified relative to the center of the input image (either the true center if use_true_center=True, or the nominal center if use_true_center=False). [default: None]
gsparams – An optional GSParams argument. [default: None]
hdu – When reading in an Image from a file, this parameter can be used to select a particular HDU in the file. [default: None]
_force_stepk – Override the normal stepk calculation (using gsparams.folding_threshold) and force stepk to the given value. [default: 0]
_force_maxk – Override the normal maxk calculation (using gsparams.maxk_threshold) and force maxk to the given value. This option in particular can help reduce FFT artifacts in a manner that is currently unobtainable by lowering maxk_threshold. [default: 0]

property image: The underlying Image being interpolated.

property k_interpolant: The Fourier-space Interpolant for this profile.

withGSParams(gsparams=None, **kwargs)[source]

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.

property x_interpolant: The real-space Interpolant for this profile.

galsim._InterpolatedImage(image, x_interpolant=galsim.Quintic(gsparams=galsim.GSParams(128, 8192, 0.005, 5.0, 0.001, 1e-05, 1e-05, 1, 0.0001, 1e-06, 1e-06, 1e-08, 1e-05)), k_interpolant=galsim.Quintic(gsparams=galsim.GSParams(128, 8192, 0.005, 5.0, 0.001, 1e-05, 1e-05, 1, 0.0001, 1e-06, 1e-06, 1e-08, 1e-05)), use_true_center=True, offset=None, gsparams=None, force_stepk=0.0, force_maxk=0.0)[source]

Approximately equivalent to InterpolatedImage, but with fewer options and no sanity checks.

Some notable reductions in functionality relative to InterpolatedImage:

There are no padding options. The image must be provided with all padding already applied.
The stepk and maxk values will not be calculated. If you want to use values for these other than the default, you may provide them as force_stepk and force_maxk. Otherwise stepk ~= 2pi / image_size and maxk ~= x_interpolant.krange() / pixel_scale.
The flux is just the flux of the image. It cannot be rescaled to a different flux value.
The input image must have a defined wcs.
The image is not recentered to have its center at (0,0). The returned profile will be centered wherever the (0,0) location is in the image, possibly with an offset governed by offset and use_true_center. If you want to mimic the behavior of the regular InterpolatedImage initializer, you can call image.setCenter(0,0) before calling this function.

Parameters:

image – The Image from which to construct the object.
x_interpolant – An Interpolant instance for real-space interpolation [default: Quintic]
k_interpolant – An Interpolant instance for k-space interpolation [default: Quintic]
use_true_center – Whether to adjust the offset by the difference between the integer center and the true center. For odd-sized images, this does nothing, but for even-sized dimensions, it adjusts the offset by -0.5. [default: True]
offset – The location in the input image to use as the center of the profile relative to position (0,0). [default: None]
gsparams – An optional GSParams argument. [default: None]
force_stepk – A stepk value to use rather than the default value. [default: 0.]
force_maxk – A maxk value to use rather than the default value. [default: 0.]

Returns:

an InterpolatedImage instance

Interpolated Fourier-space Images

class galsim.InterpolatedKImage(kimage=None, k_interpolant=None, stepk=None, gsparams=None, real_kimage=None, imag_kimage=None, real_hdu=None, imag_hdu=None)[source]

Bases: GSObject

A class describing non-parametric profiles specified by samples of their complex Fourier transform.

The InterpolatedKImage class is useful if you have a non-parametric description of the Fourier transform of the profile (provided as either a complex Image or two images giving the real and imaginary parts) that you wish to manipulate / transform using GSObject methods such as GSObject.shear, GSObject.magnify, GSObject.shift, etc. Note that neither real-space convolution nor photon-shooting of InterpolatedKImages is currently implemented. Please submit an issue at http://github.com/GalSim-developers/GalSim/issues if you require either of these use cases.

The images required for creating an InterpolatedKImage are precisely those returned by the GSObject.drawKImage method. The a and b objects in the following command will produce essentially equivalent images when drawn with the GSObject.drawImage method:

>>> a = returns_a_GSObject()
>>> b = galsim.InterpolatedKImage(a.drawKImage())

The input kimage must have dtype=numpy.complex64 or dtype=numpy.complex128, which are also known as ImageCF and ImageCD objects respectively. The only wcs permitted is a simple PixelScale (or OffsetWCS), in which case kimage.scale is used for the stepk value unless overridden by the stepk initialization argument.

Furthermore, the complex-valued Fourier profile given by kimage must be Hermitian, since it represents a real-valued real-space profile. (To see an example of valid input to InterpolatedKImage, you can look at the output of GSObject.drawKImage).

The user may optionally specify an interpolant, k_interpolant, for Fourier-space manipulations (e.g., shearing, resampling). If none is specified, then by default, a Quintic interpolant is used. The Quintic interpolant has been found to be a good compromise between speed and accuracy for real-and Fourier-space interpolation of objects specified by samples of their real-space profiles (e.g., in InterpolatedImage), though no extensive testing has been performed for objects specified by samples of their Fourier-space profiles (e.g., this class).

Example:

>>> interpolated_kimage = galsim.InterpolatedKImage(kimage, k_interpolant=None, stepk=0.,
                                                    gsparams=None)

Initializes interpolated_kimage as an InterpolatedKImage instance.

Parameters:

kimage – The complex Image corresponding to the Fourier-space samples.
k_interpolant – Either an Interpolant instance or a string indicating which k-space interpolant should be used. Options are ‘nearest’, ‘sinc’, ‘linear’, ‘cubic’, ‘quintic’, or ‘lanczosN’ where N should be the integer order to use. [default: galsim.Quintic()]
stepk – By default, the stepk value (the sampling frequency in Fourier-space) of the profile is set by the scale attribute of the supplied images. This keyword allows the user to specify a coarser sampling in Fourier- space, which may increase efficiency at the expense of decreasing the separation between neighboring copies of the DFT-rendered real-space profile. (See the GSParams docstring for the parameter folding_threshold for more information). [default: kimage.scale]
gsparams – An optional GSParams argument. [default: None]
real_kimage – Optionally, rather than provide kimage, you may provide the real and imaginary parts separately. These separate real-valued images may be strings, in which case they refer to FITS files from which to read the images. [default: None]
imag_kimage – The imaginary image corresponding to real_kimage. [default: None]
real_hdu – When reading in real_kimage from a file, this parameter can be used to select a particular HDU in the file. [default: None]
imag_hdu – When reading in imag_kimage from a file, this parameter can be used to select a particular HDU in the file. [default: None]

property k_interpolant: The Fourier-space Interpolant for this profile.

property kimage: The underlying Image being interpolated.

withGSParams(gsparams=None, **kwargs)[source]

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.

galsim._InterpolatedKImage(kimage, k_interpolant, gsparams)[source]

Approximately equivalent to InterpolatedKImage, but with fewer options and no sanity checks.

Parameters:

kimage – The complex Image corresponding to the Fourier-space samples.
k_interpolant – An Interpolant instance indicating which k-space interpolant should be used.
gsparams – An optional GSParams argument. [default: None]

Shapelet Decomposition

class galsim.Shapelet(sigma, order, bvec=None, gsparams=None)[source]

Bases: GSObject

A class describing polar shapelet surface brightness profiles.

This class describes an arbitrary profile in terms of a shapelet decomposition. A shapelet decomposition is an eigenfunction decomposition of a 2-d function using the eigenfunctions of the 2-d quantum harmonic oscillator. The functions are Laguerre polynomials multiplied by a Gaussian. See Bernstein & Jarvis, 2002 or Massey & Refregier, 2005 for more detailed information about this kind of decomposition. For this class, we follow the notation of Bernstein & Jarvis.

The decomposition is described by an overall scale length, sigma, and a vector of coefficients, b. The b vector is indexed by two values, which can be either (p,q) or (N,m). In terms of the quantum solution of the 2-d harmonic oscillator, p and q are the number of quanta with positive and negative angular momentum (respectively). Then, N=p+q, m=p-q.

The 2D image is given by (in polar coordinates):

\[I(r,\theta) = \frac{1}{\sigma^2} \sum_{pq} b_{pq} \psi_{pq}(r/\sigma, \theta)\]

where \(\psi_{pq}\) are the shapelet eigenfunctions, given by:

\[\psi_pq(r,\theta) = \frac{(-)^q}{\sqrt{\pi}} \sqrt{\frac{q!}{p!}} r^m \exp(i m \theta) \exp(-r^2/2) L_q^{(m)}(r^2)\]

and \(L_q^{(m)}(x)\) are generalized Laguerre polynomials.

The coeffients \(b_{pq}\) are in general complex. However, we require that the resulting \(I(r,\theta)\) be purely real, which implies that \(b_{pq} = b_{qp}^*\) (where \({}^*\) means complex conjugate). This further implies that \(b_{pp}\) (i.e. \(b_{pq}\) with \(p==q\)) is real.

Make a blank Shapelet instance with all \(b_{pq} = 0.\):
```
>>> shapelet = galsim.Shapelet(sigma, order)
```

Make a Shapelet instance using a given vector for the \(b_{pq}\) values.:

>>> order = 2
>>> bvec = [ 1, 0, 0, 0.2, 0.3, -0.1 ]
>>> shapelet = galsim.Shapelet(sigma, order, bvec)

We use the following order for the coefficients, where the subscripts are in terms of p,q.

[ b00 Re(b10) Im(b10) Re(b20) Im(b20) b11 Re(b30) Im(b30) Re(b21) Im(b21) … ]

i.e. we progressively increase N, and for each value of N, we start with m=N and go down to m=0 or 1 as appropriate. And since m=0 is intrinsically real, it only requires one spot in the list.

Parameters:

sigma – The scale size in the standard units (usually arcsec).
order – The order of the shapelet decomposition. This is the maximum N=p+q included in the decomposition.
bvec – The initial vector of coefficients. [default: None, which means to use all zeros]
gsparams – An optional GSParams argument. [default: None]

property bvec: The vector of shapelet coefficients

classmethod fit(sigma, order, image, center=None, normalization='flux', gsparams=None)[source]

Fit for a shapelet decomposition of a given image.

The optional normalization parameter mirrors the parameter of the InterpolatedImage class. The following sequence should produce drawn images that are approximate matches to the original image:

>>> image = [...]
>>> shapelet = galsim.FitShapelet(sigma, order, image, normalization='sb')
>>> im2 = shapelet.drawImage(image=im2, scale=image.scale, method='sb')
>>> shapelet = galsim.FitShapelet(sigma, order, image, normalization='flux')
>>> im3 = shapelet.drawImage(image=im3, scale=image.scale, method='no_pixel')

Then im2 and im3 should be as close as possible to image for the given sigma and order. Increasing the order can improve the fit, as can having sigma match the natural scale size of the image. However, it should be noted that some images are not well fit by a shapelet for any (reasonable) order.

Parameters:

sigma – The scale size in the standard units (usually arcsec).
order – The order of the shapelet decomposition. This is the maximum N=p+q included in the decomposition.
image – The Image for which to fit the shapelet decomposition
center – The position in pixels to use for the center of the decomposition. [default: image.true_center]
normalization – The normalization to assume for the image. [default: “flux”]
gsparams – An optional GSParams argument. [default: None]

Returns:

the fitted Shapelet profile

getNM(N, m)[source]

Return the coefficient according to N,m rather than p,q where N=p+q and m=p-q.

Parameters:

N – The value of N=p+q to get.
m – The value of m=p-q to get.

Returns:

a tuple (Re(b_pq), Im(b_pq))

getPQ(p, q)[source]

Return the (p,q) coefficient.

Parameters:

p – The p index to get.
q – The q index to get.

Returns:

a tuple (Re(b_pq), Im(b_pq))

property order: The shapelet order.

property sigma: The scale size, sigma.

classmethod size(order)[source]: The size of the shapelet vector.