# 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__ = [ 'BaseCorrelatedNoise', 'CorrelatedNoise', 'UncorrelatedNoise',
'getCOSMOSNoise', 'CovarianceSpectrum', ]
import numpy as np
import os
import math
from .image import Image
from .random import BaseDeviate
from .gsparams import GSParams
from .wcs import PixelScale, BaseWCS, compatible
from . import utilities
from . import fits
from . import meta_data
from .errors import GalSimError, GalSimValueError, GalSimRangeError, GalSimUndefinedBoundsError
from .errors import GalSimIncompatibleValuesError, galsim_warn
from .interpolant import Linear
from .position import _PositionD
from .angle import radians
from .random import GaussianDeviate
from .box import Pixel
from .convolve import Convolve, AutoCorrelate, AutoConvolve
from .interpolatedimage import InterpolatedImage, InterpolatedKImage
def whitenNoise(self, noise):
# This will be inserted into the Image class as a method. So self = image.
"""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.
"""
return noise.whitenImage(self)
def symmetrizeNoise(self, noise, order=4):
# This will be inserted into the Image class as a method. So self = image.
"""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.
"""
return noise.symmetrizeImage(self, order=order)
# Now inject whitenNoise and symmetrizeNoise as methods of the Image class.
Image.whitenNoise = whitenNoise
Image.symmetrizeNoise = symmetrizeNoise
###
# Now a standalone utility function for generating noise according to an input (square rooted)
# Power Spectrum
#
def _generate_noise_from_rootps(rng, shape, rootps):
# Utility function for generating a NumPy array containing a Gaussian random noise field with
# a user-specified power spectrum also supplied as a NumPy array.
# shape is the shape of the output array, needed because of the use of Hermitian symmetry to
# increase inverse FFT efficiency using the np.fft.irfft2 function (gets sent
# to the kwarg s of np.fft.irfft2)
# rootps is a NumPy array containing the square root of the discrete power spectrum ordered
# in two dimensions according to the usual DFT pattern for np.fft.rfft2 output
# Returns a NumPy array (contiguous) of the requested shape, filled with the noise field.
# Quickest to create Gaussian rng each time needed, so do that here...
# Note sigma scaling: 1/sqrt(2) needed so <|gaussvec|**2> = product(shape)
# shape needed because of the asymmetry in the 1/N^2 division in the NumPy FFT/iFFT
gd = GaussianDeviate(rng, sigma=np.sqrt(.5 * shape[0] * shape[1]))
# Fill a couple of arrays with this noise
gvec_real = utilities.rand_arr((shape[0], shape[1]//2+1), gd)
gvec_imag = utilities.rand_arr((shape[0], shape[1]//2+1), gd)
# Prepare a complex vector upon which to impose Hermitian symmetry
gvec = gvec_real + 1J * gvec_imag
# Now impose requirements of Hermitian symmetry on random Gaussian halfcomplex array, and ensure
# self-conjugate elements (e.g. [0, 0]) are purely real and multiplied by sqrt(2) to compensate
# for lost variance, see https://github.com/GalSim-developers/GalSim/issues/563
# First do the bits necessary for both odd and even shapes:
gvec[-1:shape[0]//2:-1, 0] = np.conj(gvec[1:(shape[0]+1)//2, 0])
rt2 = np.sqrt(2.)
gvec[0, 0] = rt2 * gvec[0, 0].real
# Then make the changes necessary for even sized arrays
if shape[1] % 2 == 0: # x dimension even
gvec[-1:shape[0]//2:-1, shape[1]//2] = np.conj(gvec[1:(shape[0]+1)//2, shape[1]//2])
gvec[0, shape[1]//2] = rt2 * gvec[0, shape[1]//2].real
if shape[0] % 2 == 0: # y dimension even
gvec[shape[0]//2, 0] = rt2 * gvec[shape[0]//2, 0].real
# Both dimensions even
if shape[1] % 2 == 0:
gvec[shape[0]//2, shape[1]//2] = rt2 * gvec[shape[0]//2, shape[1]//2].real
# Finally generate and return noise using the irfft
return np.fft.irfft2(gvec * rootps, s=shape)
###
# Then we define the CorrelatedNoise, which generates a correlation function by estimating it
# directly from images:
#
def _cf_periodicity_dilution_correction(cf_shape):
"""Return an array containing the correction factor required for wrongly assuming periodicity
around noise field edges in an DFT estimate of the discrete correlation function.
Uses the result calculated by MJ on GalSim Pull Request #366.
See https://github.com/GalSim-developers/GalSim/pull/366.
Returns a 2D NumPy array with the same shape as the input parameter tuple ``cf_shape``. This
array contains the correction factor by which elements in the naive CorrelatedNoise estimate of
the discrete correlation function should be multiplied to correct for the erroneous assumption
of periodic boundaries in an input noise field.
Note this should be applied only to correlation functions that have *not* been rolled to place
the origin at the array centre. The convention used here is that the lower left corner is the
[0, 0] origin, following standard FFT conventions (see e.g numpy.fft.fftfreq). You should
therefore only apply this correction before using galsim.utilities.roll2d() to recentre the
image of the correlation function.
"""
# First calculate the Delta_x, Delta_y
deltax, deltay = np.meshgrid( # Remember NumPy array shapes are [y, x]
np.fft.fftfreq(cf_shape[1]) * float(cf_shape[1]),
np.fft.fftfreq(cf_shape[0]) * float(cf_shape[0]))
# Then get the dilution correction
correction = (
cf_shape[1] * cf_shape[0] / (cf_shape[1] - np.abs(deltax)) / (cf_shape[0] - np.abs(deltay)))
return correction
# Free function for returning a COSMOS noise field correlation function
[docs]def getCOSMOSNoise(file_name=None, rng=None, cosmos_scale=0.03, variance=0., x_interpolant=None,
gsparams=None):
"""Returns a representation of correlated noise in the HST COSMOS F814W unrotated science coadd
images.
See http://cosmos.astro.caltech.edu/astronomer/hst.html for information about the COSMOS survey,
and Leauthaud et al (2007) for detailed information about the unrotated F814W coadds used for
weak lensing science.
This function uses a stacked estimate of the correlation function in COSMOS noise fields.
The correlation function was computed by the GalSim team as described in::
GalSim/devel/external/hst/make_cosmos_cfimage.py
The resulting file is distributed with GalSim as::
os.path.join('galsim.meta_data.share_dir', 'acs_I_unrot_sci_20_cf.fits')
Parameters:
file_name: If provided, override the usual location of the file with the given
file name. [default: None]
rng: If provided, a random number generator to use as the random number
generator of the resulting noise object. (may be any kind of
`BaseDeviate` object) [default: None, in which case, one will be
automatically created, using the time as a seed.]
cosmos_scale: COSMOS ACS F814W coadd image pixel scale in the units you are using to
describe `GSObject` instances and image scales in GalSim. [default: 0.03
(arcsec), see below for more information.]
variance: Scales the correlation function so that its point variance, equivalent
to its value at zero separation distance, matches this value.
[default: 0., which means to use the variance in the original COSMOS
noise fields.]
x_interpolant: Forces use of a non-default interpolant for interpolation of the
internal lookup table in real space. See below for more details.
[default: galsim.Linear()]
gsparams: An optional `GSParams` argument. [default: None]
Returns:
a `BaseCorrelatedNoise` instance representing correlated noise in F814W COSMOS images.
The default ``x_interpolant`` is a ``galsim.Linear()``, which uses bilinear interpolation.
The use of this interpolant is an approximation that gives good empirical results without
requiring internal convolution of the correlation function profile by a `Pixel` object when
applying correlated noise to images: such an internal convolution has been found to be
computationally costly in practice, requiring the Fourier transform of very large arrays.
The use of the bilinear interpolants means that the representation of correlated noise will be
noticeably inaccurate in at least the following two regimes:
1. If the pixel scale of the desired final output (e.g. the target image of
`BaseCorrelatedNoise.drawImage`, `BaseCorrelatedNoise.applyTo` or
`BaseCorrelatedNoise.whitenImage`) is small relative to the separation between pixels in
the ``image`` used to instantiate ``cn`` as shown above.
2. If the `BaseCorrelatedNoise` instance ``cn`` was instantiated with an image of scale
comparable to that in the final output, and ``cn`` has been rotated or otherwise transformed
(e.g. via the `BaseCorrelatedNoise.rotate`, `BaseCorrelatedNoise.shear` methods; see below).
Conversely, the approximation will work best in the case where the correlated noise used to
instantiate the ``cn`` is taken from an input image for which ``image.scale`` is smaller than
that in the desired output. This is the most common use case in the practical treatment of
correlated noise when simulating galaxies from COSMOS, for which this function is expressly
designed.
Changing from the default bilinear interpolant is made possible, but not recommended.
In our validation tests, we found that the Linear interpolant usually gave the most
accurate results.
You may also specify a gsparams argument. See the docstring for `GSParams` for more
information about this option.
.. note::
The ACS coadd images in COSMOS have a pixel scale of 0.03 arcsec, and so the pixel scale
``cosmos_scale`` adopted in the representation of of the correlation function takes a
default value ``cosmos_scale = 0.03``
If you wish to use other units, ensure that the input keyword ``cosmos_scale`` takes the
value corresponding to 0.03 arcsec in your chosen system.
**Example**:
The following commands use this function to generate a 300 pixel x 300 pixel image of noise with
HST COSMOS correlation properties (substitute in your own file and path for the ``filestring``)::
>>> rng = galsim.BaseDeviate(123456)
>>> noise = galsim.getCOSMOSNoise(rng=rng)
>>> image = galsim.ImageD(nx, ny, scale=0.03)
>>> image.addNoise(cf)
If your image has some other pixel scale or a complicated WCS, then the applied noise will
have the correct correlations in world coordinates, which may not be what you wanted if you
expected the pixel-to-pixel correlations to match the COSMOS noise profile. However, in
this case, you would want to view your image with the COSMOS pixel scale when you apply
the noise::
>>> image = galsim.Image(nx, ny, wcs=complicated_wcs)
>>> noise = galsim.getCOSMOSNoise(rng=rng)
>>> image.view(wcs=noise.wcs).addNoise(noise)
The FITS file ``out.fits`` should then contain an image of randomly-generated, COSMOS-like noise.
"""
if file_name is None:
file_name = os.path.join(meta_data.share_dir,'acs_I_unrot_sci_20_cf.fits')
return BaseCorrelatedNoise.from_file(file_name, cosmos_scale,
rng=rng, variance=variance,
x_interpolant=x_interpolant,
gsparams=gsparams)
[docs]class CovarianceSpectrum:
"""Class to hold a `ChromaticSum` noise covariance spectrum (which is a generalization of a
power spectrum or equivalently a correlation function).
Analogous to how a `galsim.CorrelatedNoise` object stores the variance and covariance of a
`galsim.Image` object, a `galsim.CovarianceSpectrum` stores the variance and covariance of the
Fourier mode amplitudes in different components of a `ChromaticSum`.
Note that the covariance in question exists between different `SED` components of the
`ChromaticSum`, and not between different Fourier modes, which are assumed to be uncorrelated.
This structure arises naturally for a `ChromaticRealGalaxy` (see devel/modules/CGNotes.pdf for
more details).
Parameters:
Sigma: A dictionary whose keys are tuples numerically indicating a pair of
`ChromaticSum` components whose Fourier mode amplitude covariances are
described by the corresponding `GSObject` values.
SEDs: `SED` instances of associated `ChromaticSum` components.
"""
def __init__(self, Sigma, SEDs):
self.Sigma = Sigma
self.SEDs = SEDs
def __mul__(self, variance_ratio):
return self.withScaledVariance(variance_ratio)
def transform(self, dudx, dudy, dvdx, dvdy):
Sigma = {}
for k, v in self.Sigma.items():
Sigma[k] = v.transform(dudx, dudy, dvdx, dvdy)
return CovarianceSpectrum(Sigma, self.SEDs)
def withScaledVariance(self, variance_ratio):
Sigma = {}
for k, v in self.Sigma.items():
Sigma[k] = v * variance_ratio
return CovarianceSpectrum(Sigma, self.SEDs)
[docs] def toNoise(self, bandpass, PSF, wcs, rng=None):
"""Derive the `CorrelatedNoise` object for the associated `ChromaticSum` when convolved
with ``PSF`` and drawn through ``bandpass`` onto pixels with specified ``wcs``.
Parameters:
bandpass: `Bandpass` object representing filter image is drawn through.
PSF: output chromatic PSF to convolve by.
wcs: WCS of output pixel scale.
rng: Random number generator to forward to resulting CorrelatedNoise object.
Returns:
CorrelatedNoise object.
"""
NSED = len(self.SEDs)
maxk = np.min([PSF.evaluateAtWavelength(bandpass.blue_limit).maxk,
PSF.evaluateAtWavelength(bandpass.red_limit).maxk])
stepk = np.max([PSF.evaluateAtWavelength(bandpass.blue_limit).stepk,
PSF.evaluateAtWavelength(bandpass.red_limit).stepk])
nk = 2*int(np.ceil(maxk/stepk))
PSF_eff_kimgs = np.empty((NSED, nk, nk), dtype=np.complex128)
for i, sed in enumerate(self.SEDs):
# Assume that PSF does not yet include pixel contribution, so add it in.
conv = Convolve(PSF, Pixel(wcs.scale, gsparams=PSF.gsparams)) * sed
PSF_eff_kimgs[i] = conv.drawKImage(bandpass, nx=nk, ny=nk, scale=stepk).array
pkout = np.zeros((nk, nk), dtype=np.float64)
for i in range(NSED):
for j in range(i, NSED):
s = self.Sigma[(i, j)].drawKImage(nx=nk, ny=nk, scale=stepk).array
pkout += (np.conj(PSF_eff_kimgs[i]) * s * PSF_eff_kimgs[j] *
(2 if i != j else 1)).real
pk = Image(pkout + 0j, scale=stepk, dtype=complex) # imag part should be zero
iki = InterpolatedKImage(pk)
iki *= wcs.pixelArea()**2 # determined this empirically
return BaseCorrelatedNoise(rng, iki, wcs)
def __eq__(self, other):
return (self is other or
(isinstance(other, CovarianceSpectrum) and
self.SEDs == other.SEDs and
self.Sigma == other.Sigma))
def __ne__(self, other): return not self.__eq__(other)
def __hash__(self):
return hash(("galsim.CovarianceSpectrum", tuple(self.SEDs),
frozenset(self.Sigma.items())))
def __repr__(self):
return "galsim.CovarianceSpectrum(%r, %r)" % (self.Sigma, self.SEDs)
def __str__(self):
sigma_str = '{' + ', '.join([':'.join((str(k),str(v))) for k,v in self.Sigma.items()]) + '}'
seds_str = '[' + ', '.join([str(s) for s in self.SEDs]) + ']'
return "galsim.CovarianceSpectrum(%s, %s)" % (sigma_str, seds_str)