The Image class

class galsim.Image(*args, **kwargs)

A class for storing image data along with the pixel scale or WCS information

The Image class encapsulates all the relevant information about an image including a NumPy array for the pixel values, a bounding box, and some kind of WCS that converts between pixel coordinates and world coordinates. The NumPy array may be constructed by the Image class itself, or an existing array can be provided by the user.

This class creates shallow copies unless a deep copy is explicitly requested using the copy method. The main reason for this is that it allows users to work directly with and modify subimages of larger images (for example, to successively draw many galaxies into one large image). For other implications of this convention, see the description of initialization instructions below.

In most applications with images, we will use (x,y) to refer to the coordinates. We adopt the same meaning for these coordinates as most astronomy applications do: ds9, SAOImage, SExtractor, etc. all treat x as the column number and y as the row number. However, this is different from the default convention used by numpy. In numpy, the access is by [row_num,col_num], which means this is really [y,x] in terms of the normal x,y values. Users are typically insulated from this concern by the Image API, but if you access the numpy array directly via the array attribute, you will need to be careful about this difference.

There are 6 data types that the Image can use for the data values. These are numpy.uint16, numpy.uint32, numpy.int16, numpy.int32, numpy.float32, and numpy.float64. If you are constructing a new Image from scratch, the default is numpy.float32, but you can specify one of the other data types.

There are several ways to construct an Image: (Optional arguments are shown with their default values after the = sign.)

Image(ncol, nrow, dtype=numpy.float32, init_value=0, xmin=1, ymin=1, ...)

This constructs a new image, allocating memory for the pixel values according to the number of columns and rows. You can specify the data type as dtype if you want. The default is numpy.float32 if you don’t specify it. You can also optionally provide an initial value for the pixels, which defaults to 0. The optional xmin,ymin allow you to specify the location of the lower-left pixel, which defaults to (1,1). Reminder, with our convention for x,y coordinates described above, ncol is the number of pixels in the x direction, and nrow is the number of pixels in the y direction.

Image(bounds, dtype=numpy.float32, init_value=0, ...)

This constructs a new image, allocating memory for the pixel values according to a given Bounds object. Particularly, the bounds should be a BoundsI instance. You can specify the data type as dtype if you want. The default is numpy.float32 if you don’t specify it. You can also optionally provide an initial value for the pixels, which defaults to 0.

Image(array, xmin=1, ymin=1, make_const=False, copy=False ...)

This views an existing NumPy array as an Image, where updates to either the image or the original array will affect the other one. The data type is taken from array.dtype, which must be one of the allowed types listed above. You can also optionally set the origin xmin, ymin if you want it to be something other than (1,1).

You can also optionally force the Image to be read-only with make_const=True, though if the original NumPy array is modified then the contents of Image.array will change.

If you want to make a copy of the input array, rather than just view the existing array, you can force a copy with:

>>> image = galsim.Image(array, copy=True)

Image(image, dtype=image.dtype, copy=True)

This creates a copy of an Image, possibly changing the type. e.g.:

>>> image_float = galsim.Image(64, 64) # default dtype=numpy.float32
>>> image_double = galsim.Image(image_float, dtype=numpy.float64)

You can see a list of valid values for dtype in galsim.Image.valid_dtypes. Without the dtype argument, this is equivalent to image.copy(), which makes a deep copy. If you want a copy that shares data with the original, see the view method.

If you only want to enforce the image to have a given type and not make a copy if the array is already the correct type, you can use, e.g.:

>>> image_double = galsim.Image(image, dtype=numpy.float64, copy=False)

You can specify the ncol, nrow, bounds, array, or image parameters by keyword argument if you want, or you can pass them as simple arg as shown aboves, and the constructor will figure out what they are.

The other keyword arguments (shown as … above) relate to the conversion between sky coordinates, which is how all the GalSim objects are defined, and the pixel coordinates. There are three options for this:

scale

You can optionally specify a pixel scale to use. This would normally have units arcsec/pixel, but it doesn’t have to be arcsec. If you want to use different units for the physical scale of your galsim objects, then the same unit would be used here.

wcs

A WCS object that provides a non-trivial mapping between sky units and pixel units. The scale parameter is equivalent to wcs=PixelScale(scale). But there are a number of more complicated options. See the WCS class for more details.

None

If you do not provide either of the above, then the conversion is undefined. When drawing onto such an image, a suitable pixel scale will be automatically set according to the Nyquist scale of the object being drawn.

After construction, you can set or change the scale or wcs with:

>>> image.scale = new_scale
>>> image.wcs = new_wcs

Note that image.scale will only work if the WCS is a PixelScale. Once you set the wcs to be something non-trivial, then you must interact with it via the wcs attribute. The image.scale syntax will raise an exception.

There are also two read-only attributes:

>>> image.bounds
>>> image.array

The array attribute is a NumPy array of the Image’s pixels. The individual elements in the array attribute are accessed as image.array[y,x], matching the standard NumPy convention, while the Image class’s own accessor uses either (x,y) or [x,y].

That is, the following are equivalent:

>>> ixy = image(x,y)
>>> ixy = image[x,y]
>>> ixy = image.array[y,x]
>>> ixy = image.getValue(x,y)

Similarly, for setting individual pixel values, the following are equivalent:

>>> image[x,y] = new_ixy
>>> image.array[y,x] = new_ixy
>>> image.setValue(x,y,new_ixy)

__getitem__(*args)

Return either a subimage or a single pixel value.

For example,:

>>> subimage = im[galsim.BoundsI(3,7,3,7)]
>>> value = im[galsim.PositionI(5,5)]
>>> value = im[5,5]

__setitem__(*args)

Set either a subimage or a single pixel to new values.

For example,:

>>> im[galsim.BoundsI(3,7,3,7)] = im2
>>> im[galsim.PositionI(5,5)] = 17.
>>> im[5,5] = 17.

_wrap(bounds, hermx, hermy)

A version of wrap without the sanity checks.

Equivalent to image.wrap(bounds, hermitian='x' if hermx else 'y' if hermy else False)

_view()

Equivalent to view, but without some of the sanity checks and extra options.

_shift(delta)

Equivalent to shift, but without some of the sanity checks and delta must be a PositionI instance.

Parameters:

delta – The amount to shift as a PositionI.

__call__(*args, **kwargs)

Get the pixel value at given position

The arguments here may be either (x, y) or a PositionI instance. Or you can provide x, y as named kwargs.

_getValue(x, y)

Equivalent to getValue, except there are no checks that the values fall within the bounds of the image.

_setValue(x, y, value)

Equivalent to setValue except that there are no checks that the values fall within the bounds of the image, and the coordinates must be given as x, y.

Parameters:

• x – The x coordinate of the pixel to set.

• y – The y coordinate of the pixel to set.

• value – The value to set the pixel to.

_addValue(x, y, value)

Equivalent to addValue except that there are no checks that the values fall within the bounds of the image, and the coordinates must be given as x, y.

Parameters:

• x – The x coordinate of the pixel to add to.

• y – The y coordinate of the pixel to add to.

• value – The value to add to this pixel.

_fill(value)

Equivalent to fill, except that there are no checks that the bounds are defined.

_invertSelf()

Equivalent to invertSelf, except that there are no checks that the bounds are defined.

FindAdaptiveMom(weight=None, badpix=None, guess_sig=5.0, precision=1e-06, guess_centroid=None, strict=True, check=True, round_moments=False, hsmparams=None, use_sky_coords=False)

Measure adaptive moments of an object.

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

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

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

Note

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

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

Example:

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

or:

>>> my_moments = my_gaussian_image.FindAdaptiveMom()

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

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

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

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

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

>>> my_moments = my_gaussian_image.FindAdaptiveMom()

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

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

Parameters:

• object_image – The Image for the object being measured.

• weight – The optional weight image for the object being measured. Can be an int or a float array. Currently, GalSim does not account for the variation in non-zero weights, i.e., a weight map is converted to an image with 0 and 1 for pixels that are not and are used. Full use of spatial variation in non-zero weights will be included in a future version of the code. [default: None]

• badpix – The optional bad pixel mask for the image being used. Zero should be used for pixels that are good, and any nonzero value indicates a bad pixel. [default: None]

• guess_sig – Optional argument with an initial guess for the Gaussian sigma of the object (in pixels). [default: 5.0]

• precision – The convergence criterion for the moments. [default: 1e-6]

• guess_centroid – An initial guess for the object centroid (useful in case it is not located at the center, which is used if this keyword is not set). The convention for centroids is such that the center of the lower-left pixel is (image.xmin, image.ymin). [default: object_image.true_center]

• strict – Whether to require success. If strict=True, then there will be a GalSimHSMError exception if shear estimation fails. If set to False, then information about failures will be silently stored in the output ShapeData object. [default: True]

• check – Check if the object_image, weight and badpix are in the correct format and valid. [default: True]

• round_moments – Use a circular weight function instead of elliptical. [default: False]

• hsmparams – The hsmparams keyword can be used to change the settings used by FindAdaptiveMom when estimating moments; see HSMParams documentation for more information. [default: None]

• use_sky_coords – Whether to convert the measured moments to sky_coordinates. Setting this to true is equivalent to running applyWCS(object_image.wcs, image_pos=object_image.true_center) on the result. [default: False]

Returns:

a ShapeData object containing the results of moment measurement.

addNoise(noise)

Add noise to the image according to a supplied noise model.

Parameters:

noise – The noise (BaseNoise) model to use.

addNoiseSNR(noise, snr, preserve_flux=False)

Adds noise to the image in a way that achieves the specified signal-to-noise ratio.

The specified snr (signal-to-noise ratio, or \(S/N\)) can be achieved either by scaling the flux of the object while keeping the noise level fixed, or the flux can be preserved and the noise variance changed. This choice is specified using the parameter preserve_flux, where False means the former and True means the latter.

The definition of \(S/N\) is equivalent to the one used by Great08. We take the signal to be a weighted sum of the pixel values:

\[\begin{split}S &= \frac{\sum W(x,y) I(x,y)}{\sum W(x,y)} \\ N^2 = Var(S) &= \frac{\sum W(x,y)^2 Var(I(x,y))}{(\sum W(x,y))^2}\end{split}\]

and assume that \(Var(I(x,y))\) is constant, \(V \equiv Var(I(x,y))\). We then assume that we are using a matched filter for \(W\), so \(W(x,y) = I(x,y)\). Then a few things cancel and we find that

\[(S/N)^2 = \frac{\sum I(x,y)^2}{V}\]

and therefore, for a given \(I(x,y)\) and \(S/N\) (aka snr)

\[V = \frac{\sum I(x,y)^2}{(S/N)^2}\]

Note

For noise models such as PoissonNoise and CCDNoise, the assumption of constant \(Var(I(x,y))\) is only approximate, since the flux of the object adds to the Poisson noise in those pixels. Thus, the real \(S/N\) on the final image will be slightly lower than the target snr value, and this effect will be larger for brighter objects.

This function relies on BaseNoise.getVariance to determine how much variance the noise model will add. Thus, it will not work for noise models that do not have a well-defined variance, such as VariableGaussianNoise.

Parameters:

• noise – The noise (BaseNoise) model to use.

• snr – The desired signal-to-noise after the noise is applied.

• preserve_flux – Whether to preserve the flux of the object (True) or the variance of the noise model (False) to achieve the desired snr. [default: False]

Returns:

the variance of the noise that was applied to the image.

addReciprocityFailure(exp_time, alpha, base_flux)

Accounts for the reciprocity failure and includes it in the original Image directly.

Reciprocity, in the context of photography, is the inverse relationship between the incident flux (I) of a source object and the exposure time (t) required to produce a given response (p) in the detector, i.e., p = I*t. At very low (also at high) levels of incident flux, deviation from this relation is observed, leading to reduced sensitivity at low flux levels. The pixel response to a high flux is larger than its response to a low flux. This flux-dependent non- linearity is known as ‘Reciprocity Failure’ and is known to happen in photographic films since 1893. Interested users can refer to http://en.wikipedia.org/wiki/Reciprocity_(photography)

CCDs are not known to suffer from this effect. HgCdTe detectors that are used for near infrared astrometry, although to an extent much lesser than the photographic films, are found to exhibit reciprocity failure at low flux levels. The exact mechanism of this effect is unknown and hence we lack a good theoretical model. Many models that fit the empirical data exist and a common relation is

\[\frac{p_R}{p} = \left(1 + \alpha \log_{10}\left(\frac{p}{t}\right) - \alpha \log_{10}\left(\frac{p^\prime}{t^\prime}\right)\right)\]

where \(t\) is the exposure time (in units of seconds), \(p\) is the pixel response (in units of electrons) and \(p_R\) is the response if the reciprocity relation fails to hold. \(p^\prime/t^\prime\) is the count rate (in electrons/second) corresponding to the photon flux (base flux) at which the detector is calibrated to have its nominal gain. alpha is the parameter in the model, measured in units of per decade and varies with detectors and the operating temperature. The functional form for the reciprocity failure is motivated empirically from the tests carried out on H2RG detectors.

See for reference Fig. 1 and Fig. 2 of http://arxiv.org/abs/1106.1090. Since \(p_R/p\) remains close to unity over a wide range of flux, we convert this relation to a power law by approximating \(p_R/p \approx 1 + \log(p_R/p)\). This gives a relation that is better behaved than the logarithmic relation at low flux levels.

\[\frac{p_R}{p} = \left(\frac{p/t}{p^\prime/t^\prime}\right)^\frac{\alpha}{\log(10)}.\]

Because of how this function is defined, the input image must have non-negative pixel values for the resulting image to be well-defined. Negative pixel values result in ‘nan’s. The image should be in units of electrons, or if it is in ADU, then the value passed to exp_time should be the exposure time divided by the nominal gain. The image should include both the signal from the astronomical objects as well as the background level. The addition of nonlinearity should occur after including the effect of reciprocity failure.

Parameters:

• exp_time – The exposure time (t) in seconds, which goes into the expression for reciprocity failure given in the docstring.

• alpha – The alpha parameter in the expression for reciprocity failure, in units of ‘per decade’.

• base_flux – The flux (\(p^\prime/t^\prime\)) at which the gain is calibrated to have its nominal value.

addValue(*args, **kwargs)

Add some amount to the pixel value at given (x,y) position

The arguments here may be either (x, y, value) or (pos, value) where pos is a PositionI. Or you can provide x, y, value as named kwargs.

This is equivalent to self[x,y] += rhs

applyIPC(IPC_kernel, edge_treatment='extend', fill_value=None, kernel_nonnegativity=True, kernel_normalization=True)

Applies the effect of interpixel capacitance to the Image instance.

In NIR detectors, the quantity that is sensed is not the charge as in CCDs, but a voltage that relates to the charge present within each pixel. The voltage read at a given pixel location is influenced by the charges present in the neighboring pixel locations due to capacitive coupling of sense nodes.

This interpixel capacitance is approximated as a linear effect that can be described by a 3x3 kernel that is convolved with the image. The kernel must be an Image instance and could be intrinsically anisotropic. A sensible kernel must have non-negative entries and must be normalized such that the sum of the elements is 1, in order to conserve the total charge. The (1,1) element of the kernel is the contribution to the voltage read at a pixel from the electrons in the pixel to its bottom-left, the (1,2) element of the kernel is the contribution from the charges to its left and so on.

The argument ‘edge_treatment’ specifies how the edges of the image should be treated, which could be in one of the three ways:

1. ‘extend’: The kernel is convolved with the zero-padded image, leading to a larger intermediate image. The central portion of this image is returned. [default]

2. ‘crop’: The kernel is convolved with the image, with the kernel inside the image completely. Pixels at the edges, where the center of the kernel could not be placed, are set to the value specified by ‘fill_value’. If ‘fill_value’ is not specified or set to ‘None’, then the pixel values in the original image are retained. The user can make the edges invalid by setting fill_value to numpy.nan.

3. ‘wrap’: The kernel is convolved with the image, assuming periodic boundary conditions.

The size of the image array remains unchanged in all three cases.

Parameters:

• IPC_kernel – A 3x3 Image instance that is convolved with the Image instance

• edge_treatment – Specifies the method of handling edges and should be one of ‘crop’, ‘extend’ or ‘wrap’. See above for details. [default: ‘extend’]

• fill_value – Specifies the value (including nan) to fill the edges with when edge_treatment is ‘crop’. If unspecified or set to ‘None’, the original pixel values are retained at the edges. If edge_treatment is not ‘crop’, then this is ignored.

• kernel_nonnegativity – Specify whether the kernel should have only non-negative entries. [default: True]

• kernel_normalization – Specify whether to check and enforce correct normalization for the kernel. [default: True]

applyNonlinearity(NLfunc, *args)

Applies the given non-linearity function (NLfunc) on the Image instance directly.

This routine can transform the image in a non-linear manner specified by the user. However, the typical kind of non-linearity one sees in astronomical images is voltage non-linearity, also sometimes known as ‘classical non-linearity’, refers to the non-linearity in charge-to-voltage conversion process. This arises as charge gets integrated at the junction capacitance of the pixel node. Voltage non-linearity decreases signals at higher signal levels, causing the attenuation of brighter pixels. The image should include both the signal from the astronomical objects as well as the background level. Other detectors effects such as dark current and persistence (not currently included in GalSim) would also occur before the inclusion of nonlinearity.

The argument NLfunc is a callable function (for example a lambda function, a galsim.LookupTable, or a user-defined function), possibly with arguments that need to be given as subsequent arguments to the applyNonlinearity function (after the NLfunc argument). NLfunc should be able to take a 2d NumPy array as input, and return a NumPy array of the same shape. It should be defined such that it outputs the final image with nonlinearity included (i.e., in the limit that there is no nonlinearity, the function should return the original image, NOT zero). The image should be in units of electrons when this routine is being used to generate classical non-linearity. When used for other purposes, the units can be in electrons or in ADU, as found appropriate by the user.

Examples:

>>> f = lambda x: x + (1.e-7)*(x**2)
>>> img.applyNonlinearity(f)

>>> f = lambda x, beta1, beta2: x - beta1*x*x + beta2*x*x*x
>>> img.applyNonlinearity(f, 1.e-7, 1.e-10)

On calling the method, the Image instance img is transformed by the user-defined function f with beta1 = 1.e-7 and beta2 = 1.e-10.

Parameters:

• NLfunc – The function that maps the input image pixel values to the output image pixel values.

• *args – Any subsequent arguments are passed along to the NLfunc function.

applyPersistence(imgs, coeffs)

Applies the effects of persistence to the Image instance.

Persistence refers to the retention of a small fraction of the signal after resetting the imager pixel elements. The persistence signal of a previous exposure is left in the pixel even after several detector resets. This effect is most likely due to charge traps in the material. Laboratory tests on the Roman Space Telescope CMOS detectors show that if exposures and readouts are taken in a fixed cadence, the persistence signal can be given as a linear combination of prior pixel values that can be added to the current image.

This routine takes in a list of Image instances and adds them to Image weighted by the values passed on to ‘coeffs’. The pixel values of the Image instances in the list must correspond to the electron counts before the readout. This routine does NOT keep track of realistic dither patterns. During the image simulation process, the user has to queue a list of previous Image instances (imgs) outside the routine by inserting the latest image in the beginning of the list and deleting the oldest image. The values in ‘coeffs’ tell how much of each Image is to be added. This usually remains constant in the image generation process.

Parameters:

• imgs – A list of previous Image instances that still persist.

• coeffs – A list of floats that specifies the retention factors for the corresponding Image instances listed in ‘imgs’.

property array

The underlying numpy array.

bin(nx, ny)

Bin the image pixels in blocks of nx x ny pixels.

This returns a new image that is a binned version of the current image. Adjacent pixel values in nx x ny blocks are added together to produce the flux in each output pixel.

If the current number of pixels in each direction is not a multiple of nx, ny, then the last pixel in each direction will be the sum of fewer than nx or ny pixels as needed.

See also subsample, which is the opposite of this.

If the wcs is a Jacobian (or simpler), the output image will have its wcs set properly. But if the wcs is more complicated, the output wcs would be fairly complicated to figure out properly, so we leave it as None. The user should set it themselves if required.

Parameters:

• nx – The number of adjacent pixels in the x direction to add together into each output pixel.

• ny – The number of adjacent pixels in the y direction to add together into each output pixel.

Returns:

a new Image

property bounds

The bounds of the Image.

calculateFWHM(center=None, Imax=0.0)

Returns the full-width half-maximum (FWHM) of a drawn object.

This method is equivalent to GSObject.calculateFWHM when the object has already been drawn onto an image. Note that the profile should be drawn using a method that does not integrate over pixels, so either ‘sb’ or ‘no_pixel’. Also, if there is a significant amount of noise in the image, this method may not work well.

If the image has a wcs other than a PixelScale, an AttributeError will be raised.

Parameters:

• center – The position in pixels to use for the center, r=0. [default: self.true_center]

• Imax – The maximum surface brightness. [default: max(self.array)] Note: If Imax is provided, and the maximum pixel value is larger than this value, Imax will be updated to use the larger value.

Returns:

an estimate of the full-width half-maximum in physical units defined by the pixel scale.

calculateHLR(center=None, flux=None, flux_frac=0.5)

Returns the half-light radius of a drawn object.

This method is equivalent to GSObject.calculateHLR when the object has already been been drawn onto an image. Note that the profile should be drawn using a method that integrates over pixels and does not add noise. (The default method=’auto’ is acceptable.)

If the image has a wcs other than a PixelScale, an AttributeError will be raised.

Parameters:

• center – The position in pixels to use for the center, r=0. [default: self.true_center]

• flux – The total flux. [default: sum(self.array)]

• 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 defined by the pixel scale.

calculateMomentRadius(center=None, flux=None, rtype='det')

Returns an estimate of the radius based on unweighted second moments of a drawn object.

This method is equivalent to GSObject.calculateMomentRadius when the object has already been drawn onto an image. Note that the profile should be drawn using a method that integrates over pixels and does not add noise. (The default method=’auto’ is acceptable.)

If the image has a wcs other than a PixelScale, an AttributeError will be raised.

Parameters:

• center – The position in pixels to use for the center, r=0. [default: self.true_center]

• flux – The total flux. [default: sum(self.array)]

• 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 defined by the pixel scale (or both estimates if rtype == ‘both’).

calculate_fft()

Performs an FFT of an Image in real space to produce a k-space Image.

Note: the image will be padded with zeros as needed to make an image with bounds that look like BoundsI(-N/2, N/2-1, -N/2, N/2-1).

The input image must have a PixelScale wcs. The output image will be complex (an ImageCF or ImageCD instance) and its scale will be 2pi / (N dx), where dx is the scale of the input image.

Returns:

an Image instance with the k-space image.

calculate_inverse_fft()

Performs an inverse FFT of an Image in k-space to produce a real-space Image.

The starting image is typically an ImageCD, although if the Fourier function is real valued, then you could get away with using an ImageD or ImageF.

The image is assumed to be Hermitian. In fact, only the portion with x >= 0 needs to be defined, with f(-x,-y) taken to be conj(f(x,y)).

Note: the k-space image will be padded with zeros and/or wrapped as needed to make an image with bounds that look like BoundsI(0, N/2, -N/2, N/2-1). If you are building a larger k-space image and then wrapping, you should wrap directly into an image of this shape.

The input image must have a PixelScale wcs. The output image will be real (an ImageD instance) and its scale will be 2pi / (N dk), where dk is the scale of the input image.

Returns:

an Image instance with the real-space image.

property center

The current nominal center (xcen,ycen) of the image as a PositionI instance.

In terms of the rows and columns, xcen is the x value for the central column, and ycen is the y value of the central row. For even-sized arrays, there is no central column or row, so the convention we adopt in this case is to round up. For example:

>>> im = galsim.Image(numpy.array(range(16),dtype=float).reshape((4,4)))
>>> im.center
galsim.PositionI(x=3, y=3)
>>> im(im.center)
10.0
>>> im.setCenter(56,72)
>>> im.center
galsim.PositionI(x=56, y=72)
>>> im(im.center)
10.0

static clear_depixelize_cache()

Release the cached solver used by depixelize to make repeated calls more efficient.

property conjugate

Return the complex conjugate of an image.

This works for real or complex. For real images, it acts the same as view.

Note that for complex images, this is not a conjugate view into the original image. So changing the original image does not change the conjugate (or vice versa).

copy()

Make a copy of the Image

copyFrom(rhs)

Copy the contents of another image

depixelize(x_interpolant)

Return a depixelized version of the image.

Specifically, this function creates an image that could be used with InterpolatedImage with the given x_interpolant, which when drawn with method=auto would produce the current image.

>>> alt_image = image.depixelize(x_interpolant)
>>> ii = galsim.InterpolatedImage(alt_image, x_interpolant=x_interpolant)
>>> image2 = ii.drawImage(image.copy(), method='auto')

image2 will end up approximately equal to the original image.

Warning

This function is fairly expensive, both in memory and CPU time, so it should only be called on fairly small images (~100x100 or smaller typically). The memory requirement scales as Npix^2, and the execution time scales as Npix^3.

However, the expensive part of the calculation is independent of the image values. It only depends on the size of the image and interpolant being used. So this part of the calculation is cached and reused if possible. If you make repeated calls to depixelize using the same image size and interpolant, it will be much faster after the first call.

If you need to release the cache (since it can be a non-trivial amount of memory), you may do so using Image.clear_depixelize_cache.

Parameters:

x_interpolant – The Interpolant to use in the InterpolatedImage to describe how the profile should be interpolated between the pixel centers.

Returns:

an Image representing the underlying profile without the pixel convolution.

property dtype

The dtype of the underlying numpy array.

fill(value)

Set all pixel values to the given value

Parameter:

value: The value to set all the pixels to.

flip_lr()

Return a version of the image flipped left to right.

Note: The returned image will have an undefined wcs. If you care about the wcs, you will need to set it yourself.

flip_ud()

Return a version of the image flipped top to bottom.

Note: The returned image will have an undefined wcs. If you care about the wcs, you will need to set it yourself.

getValue(x, y)

This method is a synonym for im(x,y). It is a bit faster than im(x,y), since GalSim does not have to parse the different options available for __call__. (i.e. im(x,y) or im(pos) or im(x=x,y=y))

Parameters:

• x – The x coordinate of the pixel to get.

• y – The y coordinate of the pixel to get.

get_pixel_centers()

A convenience function to get the x and y values at the centers of the image pixels.

Returns:

(x, y), each of which is a numpy array the same shape as self.array

classmethod good_fft_size(input_size)

Round the given input size up to the next higher power of 2 or 3 times a power of 2.

This rounds up to the next higher value that is either 2^k or 3*2^k. If you are going to be performing FFTs on an image, these will tend to be faster at performing the FFT.

property imag

Return the imaginary part of an image.

This is a property, not a function. So write im.imag, not im.imag().

This works for real or complex. For real images, the returned array is read-only and all elements are 0.

invertSelf()

Set all pixel values to their inverse: x -> 1/x.

Note: any pixels whose value is 0 originally are ignored. They remain equal to 0 on the output, rather than turning into inf.

property iscomplex

Whether the Image values are complex.

property isconst

Whether the Image is constant. I.e. modifying its values is an error.

property iscontiguous

Indicates whether each row of the image is contiguous in memory.

Note: it is ok for the end of one row to not be contiguous with the start of the next row. This just checks that each individual row has a stride of 1.

property isinteger

Whether the Image values are integral.

property ncol

The number of columns in the image

property nrow

The number of rows in the image

property origin

Return the origin of the image. i.e. the (x,y) position of the lower-left pixel.

In terms of the rows and columns, this is the (x,y) coordinate of the first column, and first row of the array. For example:

>>> im = galsim.Image(numpy.array(range(16),dtype=float).reshape((4,4)))
>>> im.origin
galsim.PositionI(x=1, y=1)
>>> im(im.origin)
0.0
>>> im.setOrigin(23,45)
>>> im.origin
galsim.PositionI(x=23, y=45)
>>> im(im.origin)
0.0
>>> im(23,45)
0.0
>>> im.bounds
galsim.BoundsI(xmin=23, xmax=26, ymin=45, ymax=48)

property outer_bounds

The bounds of the outer edge of the pixels.

Equivalent to galsim.BoundsD(im.xmin-0.5, im.xmax+0.5, im.ymin-0.5, im.ymax+0.5)

quantize()

Rounds the pixel values in an image to integer values, while preserving the type of the data.

At certain stages in the astronomical image generation process, detectors effectively round to the nearest integer. The exact stage at which this happens depends on the type of device (CCD vs. NIR detector). For example, for H2RG detectors, quantization happens in two stages: first, when detecting a certain number of photons, corresponding to the sum of background and signal multiplied by the QE and including reciprocity failure. After this, a number of other processes occur (e.g., nonlinearity, IPC, read noise) that could result in non-integer pixel values, only rounding to an integer at the stage of analog-to-digital conversion.

Because we cannot guarantee that quantization will always be the last step in the process, the quantize() routine does not actually modify the type of the image to ‘int’. However, users can easily do so by doing:

>>> image.quantize()
>>> int_image = galsim.Image(image, dtype=int)

property real

Return the real part of an image.

This is a property, not a function. So write im.real, not im.real().

This works for real or complex. For real images, it acts the same as view.

replaceNegative(replace_value=0)

Replace any negative values currently in the image with 0 (or some other value).

Sometimes FFT drawing can result in tiny negative values, which may be undesirable for some purposes. This method replaces those values with 0 or some other value if desired.

Parameters:

replace_value – The value with which to replace any negative pixels. [default: 0]

resize(bounds, wcs=None)

Resize the image to have a new bounds (must be a BoundsI instance)

Note that the resized image will have uninitialized data. If you want to preserve the existing data values, you should either use subImage (if you want a smaller portion of the current Image) or make a new Image and copy over the current values into a portion of the new image (if you are resizing to a larger Image).

Parameters:

• bounds – The new bounds to resize to.

• wcs – If provided, also update the wcs to the given value. [default: None, which means keep the existing wcs]

rot_180()

Return a version of the image rotated 180 degrees.

Note: The returned image will have an undefined wcs. If you care about the wcs, you will need to set it yourself.

rot_ccw()

Return a version of the image rotated 90 degrees counter-clockwise.

Note: The returned image will have an undefined wcs. If you care about the wcs, you will need to set it yourself.

rot_cw()

Return a version of the image rotated 90 degrees clockwise.

Note: The returned image will have an undefined wcs. If you care about the wcs, you will need to set it yourself.

property scale

The pixel scale of the Image. Only valid if the wcs is a PixelScale.

If the WCS is either not set (i.e. it is None) or it is a PixelScale, then it is permissible to change the scale with:

>>> image.scale = new_pixel_scale

setCenter(*args, **kwargs)

Set the center of the image to the given (integral) (xcen, ycen)

The arguments here may be either (xcen, ycen) or a PositionI instance. Or you can provide xcen, ycen as named kwargs.

In terms of the rows and columns, xcen is the new x value for the central column, and ycen is the new y value of the central row. For even-sized arrays, there is no central column or row, so the convention we adopt in this case is to round up. For example:

>>> im = galsim.Image(numpy.array(range(16),dtype=float).reshape((4,4)))
>>> im(1,1)
0.0
>>> im(4,1)
3.0
>>> im(4,4)
15.0
>>> im(3,3)
10.0
>>> im.setCenter(0,0)
>>> im(0,0)
10.0
>>> im(-2,-2)
0.0
>>> im(1,-2)
3.0
>>> im(1,1)
15.0
>>> im.setCenter(234,456)
>>> im(234,456)
10.0
>>> im.bounds
galsim.BoundsI(xmin=232, xmax=235, ymin=454, ymax=457)

setOrigin(*args, **kwargs)

Set the origin of the image to the given (integral) (x0, y0)

The arguments here may be either (x0, y0) or a PositionI instance. Or you can provide x0, y0 as named kwargs.

In terms of the rows and columns, x0 is the new x value for the first column, and y0 is the new y value of the first row. For example:

>>> im = galsim.Image(numpy.array(range(16),dtype=float).reshape((4,4)))
>>> im(1,1)
0.0
>>> im(4,1)
3.0
>>> im(1,4)
12.0
>>> im(4,4)
15.0
>>> im.setOrigin(0,0)
>>> im(0,0)
0.0
>>> im(3,0)
3.0
>>> im(0,3)
12.0
>>> im(3,3)
15.0
>>> im.setOrigin(234,456)
>>> im(234,456)
0.0
>>> im.bounds
galsim.BoundsI(xmin=234, xmax=237, ymin=456, ymax=459)

setSubImage(bounds, rhs)

Set a portion of the full image to the values in another image

This is equivalent to self[bounds] = rhs

setValue(*args, **kwargs)

Set the pixel value at given (x,y) position

The arguments here may be either (x, y, value) or (pos, value) where pos is a PositionI. Or you can provide x, y, value as named kwargs.

This is equivalent to self[x,y] = rhs

setZero()

Set all pixel values to zero.

shift(*args, **kwargs)

Shift the pixel coordinates by some (integral) dx,dy.

The arguments here may be either (dx, dy) or a PositionI instance. Or you can provide dx, dy as named kwargs.

In terms of columns and rows, dx means a shift in the x value of each column in the array, and dy means a shift in the y value of each row. In other words, the following will return the same value for ixy. The shift function just changes the coordinates (x,y) used for that pixel:

>>> ixy = im(x,y)
>>> im.shift(3,9)
>>> ixy = im(x+3, y+9)

subImage(bounds)

Return a view of a portion of the full image

This is equivalent to self[bounds]

subsample(nx, ny, dtype=None)

Subdivide the image pixels into nx x ny sub-pixels.

This returns a new image that is a subsampled version of the current image. Each pixel’s flux is split (uniformly) into nx x ny smaller pixels.

See also bin, which is the opposite of this. Note that subsample(nx,ny) followed by bin(nx,ny) is essentially a no op.

If the wcs is a Jacobian (or simpler), the output image will have its wcs set properly. But if the wcs is more complicated, the output wcs would be fairly complicated to figure out properly, so we leave it as None. The user should set it themselves if required.

Parameters:

• nx – The number of sub-pixels in the x direction for each original pixel.

• ny – The number of sub-pixels in the y direction for each original pixel.

• dtype – Optionally provide a dtype for the return image. [default: None, which means to use the same dtype as the original image]

Returns:

a new Image

symmetrizeNoise(noise, order=4)

Impose N-fold symmetry (where N=``order`` is an even integer >=4) on the noise in a square image assuming that the noise currently in the image can be described by the BaseCorrelatedNoise object noise. See BaseCorrelatedNoise.symmetrizeImage for more details of how this method works.

Parameters:

• noise – The BaseCorrelatedNoise model to use when figuring out how much noise to add to make the final noise have symmetry at the desired order.

• order – Desired symmetry order. Must be an even integer larger than 2. [default: 4]

Returns:

the theoretically calculated variance of the combined noise fields in the updated image.

transpose()

Return the tranpose of the image.

Note: The returned image will have an undefined wcs. If you care about the wcs, you will need to set it yourself.

property true_center

The current true center of the image as a PositionD instance.

Unline the nominal center returned by im.center, this value may be half-way between two pixels if the image has an even number of rows or columns. It gives the position (x,y) at the exact center of the image, regardless of whether this is at the center of a pixel (integer value) or halfway between two (half-integer). For example:

>>> im = galsim.Image(numpy.array(range(16),dtype=float).reshape((4,4)))
>>> im.center
galsim.PositionI(x=3, y=3)
>>> im.true_center
galsim.PositionI(x=2.5, y=2.5)
>>> im.setCenter(56,72)
>>> im.center
galsim.PositionI(x=56, y=72)
>>> im.true_center
galsim.PositionD(x=55.5, y=71.5)
>>> im.setOrigin(0,0)
>>> im.true_center
galsim.PositionD(x=1.5, y=1.5)

view(scale=None, wcs=None, origin=None, center=None, make_const=False, dtype=None, contiguous=False)

Make a view of this image, which lets you change the scale, wcs, origin, etc. but view the same underlying data as the original image.

If you do not provide either scale or wcs, the view will keep the same wcs as the current Image object.

Parameters:

• scale – If provided, use this as the pixel scale for the image. [default: None]

• wcs – If provided, use this as the wcs for the image. [default: None]

• origin – If provided, use this as the origin position of the view. [default: None]

• center – If provided, use this as the center position of the view. [default: None]

• make_const – Make the view’s data array immutable. [default: False]

• dtype – If provided, ensure that the output has this dtype. If the original Image is a different dtype, then a copy will be made. [default: None]

• contiguous – If provided, ensure that the output array is contiguous. [default: False]

whitenNoise(noise)

Whiten the noise in the image assuming that the noise currently in the image can be described by the BaseCorrelatedNoise object noise. See BaseCorrelatedNoise.whitenImage for more details of how this method works.

Parameters:

noise – The BaseCorrelatedNoise model to use when figuring out how much noise to add to make the final noise white.

Returns:

the theoretically calculated variance of the combined noise fields in the updated image.

wrap(bounds, hermitian=False)

Wrap the values in a image onto a given subimage and return the subimage.

This would typically be used on a k-space image where you initially draw a larger image than you want for the FFT and then wrap it onto a smaller subset. This will cause aliasing of course, but this is often preferable to just using the smaller image without wrapping.

For complex images of FFTs, one often only stores half the image plane with the implicit understanding that the function is Hermitian, so im(-x,-y) == im(x,y).conjugate(). In this case, the wrapping needs to work slightly differently, so you can specify that your image is implicitly Hermitian with the hermitian argument. Options are:

hermitian=False

(default) Normal non-Hermitian image.

hermitian=’x’

Only x>=0 values are stored with x<0 values being implicitly Hermitian. In this case im.bounds.xmin and bounds.xmin must be 0.

hermitian=’y’

Only y>=0 values are stored with y<0 values being implicitly Hermitian. In this case im.bounds.ymin and bounds.ymin must be 0.

Also, in the two Hermitian cases, the direction that is not implicitly Hermitian must be symmetric in the image’s bounds. The wrap bounds must be almost symmetric, but missing the most negative value. For example,:

>>> N = 100
>>> im_full = galsim.ImageCD(bounds=galsim.BoundsI(0,N/2,-N/2,N/2), scale=dk)
>>> # ... fill with im[i,j] = FT(kx=i*dk, ky=j*dk)
>>> N2 = 64
>>> im_wrap = im_full.wrap(galsim.BoundsI(0,N/2,-N2/2,N2/2-1, hermitian='x')

This sets up im_wrap to be the properly Hermitian version of the data appropriate for passing to an FFT.

Note that this routine modifies the original image (and not just the subimage onto which it is wrapped), so if you want to keep the original pristine, you should call wrapped_image = image.copy().wrap(bounds).

Parameters:

• bounds – The bounds of the subimage onto which to wrap the full image.

• hermitian – Whether the image is implicitly Hermitian and if so, whether it is the x or y values that are not stored. [default: False]

Returns:

the subimage, image[bounds], after doing the wrapping.

write(file_name=None, dir=None, hdu_list=None, clobber=True, compression='auto')

Write a single image to a FITS file.

Write the Image instance image to a FITS file, with details depending on the arguments. This function can be called directly as galsim.fits.write(image, ...), with the image as the first argument, or as an Image method: image.write(...).

Parameters:

• image – The Image to write to file. Per the description of this method, it may be given explicitly via galsim.fits.write(image, ...) or the method may be called directly as an image method, image.write(...). Note that if the image has a ‘header’ attribute containing a FitsHeader, then the FitsHeader is written to the header in the PrimaryHDU, followed by the WCS as usual.

• file_name – The name of the file to write to. [Either file_name or hdu_list is required.]

• dir – Optionally a directory name can be provided if file_name does not already include it. [default: None]

• hdu_list – An astropy.io.fits.HDUList. If this is provided instead of file_name, then the Image is appended to the end of the HDUList as a new HDU. In that case, the user is responsible for calling either hdu_list.writeto(...) or galsim.fits.writeFile(...) afterwards. [Either file_name or hdu_list is required.]

• clobber – Setting clobber=True will silently overwrite existing files. [default: True]

• compression –

  Which compression scheme to use (if any). Options are:

  • None or ‘none’ = no compression

  • ’rice’ = use rice compression in tiles (preserves header readability)

  • ’gzip’ = use gzip to compress the full file

  • ’bzip2’ = use bzip2 to compress the full file

  • ’gzip_tile’ = use gzip in tiles (preserves header readability)

  • ’hcompress’ = use hcompress in tiles (only valid for 2-d images)

  • ’plio’ = use plio compression in tiles (only valid for pos integer data)

  • ’auto’ = determine the compression from the extension of the file name (requires file_name to be given):

    • ’.fz’ => ‘rice’

    • ’.gz’ => ‘gzip’

    • ’.bz2’ => ‘bzip2’

    • otherwise None

  [default: ‘auto’]

property xmax

Alias for self.bounds.xmax.

property xmin

Alias for self.bounds.xmin.

property ymax

Alias for self.bounds.ymax.

property ymin

Alias for self.bounds.ymin.

galsim._Image(array, bounds, wcs)

Equivalent to Image(array, bounds, wcs), but without the overhead of sanity checks, and the other options for how to provide the arguments.

galsim.ImageF(*args, **kwargs)

Alias for galsim.Image(…, dtype=numpy.float32)

galsim.ImageD(*args, **kwargs)

Alias for galsim.Image(…, dtype=numpy.float64)

galsim.ImageI(*args, **kwargs)

Alias for galsim.Image(…, dtype=numpy.int32)

galsim.ImageS(*args, **kwargs)

Alias for galsim.Image(…, dtype=numpy.int16)

galsim.ImageUI(*args, **kwargs)

Alias for galsim.Image(…, dtype=numpy.uint32)

galsim.ImageUS(*args, **kwargs)

Alias for galsim.Image(…, dtype=numpy.uint16)

galsim.ImageCF(*args, **kwargs)

Alias for galsim.Image(…, dtype=numpy.complex64)

galsim.ImageCD(*args, **kwargs)

Alias for galsim.Image(…, dtype=numpy.complex128)