mccd.utils

UTILS.

These functions include several functions needed by PSF modelling algorithms.

Authors

Morgan Schmitz <https://github.com/MorganSchmitz>

Tobias Liaudat <tobias.liaudat@cea.fr>

apply_transform(data, filters)

Transform data through application of a set of filters.

Parameters

• data (numpy.ndarray) – Data to be transformed. Should be in rca_format, where the image index is contained on last/2nd axis, ie (n_pix,n_pix,n_images).

• filters (numpy.ndarray) – Set of filters. Usually the wavelet transform filters.

acc_sig_maps(shap_im, ker_stack, sig_est, flux_est, flux_ref, upfact, w, sig_data=None)

Apply acc_sig_map() several times.

Calls:

• utils.acc_sig_map()

acc_sig_map(shap_im, ker_stack, sig_est, flux_est, flux_ref, upfact, w, sig_data=None)

Estimate the simga noise maps from the observed data.

Computes the square root of \(\mathcal{F}^{2*}(\hat\sigma^2)(A^\top\odot A^\top)\). See equation (27) in RCA paper (Ngole et al.).

Notes

\(\mathrm{Var}(B)\) has been replaced by the noise level as estimated from the data, and here we do not have the term \(\mu\) (gradient step size in the paper).

return_neighbors(new_pos, obs_pos, vals, n_neighbors)

Find the nearest neighbors locally in one ccd.

rca_format(cube)

Switch from regular format to RCA format.

RCA format: image index is contained on last axis [:,:,it] Regular format: image index is contained on first axis [it,:,:]

reg_format(rca_cube)

Switch from RCA format to regular format.

RCA format: image index is contained on last axis [:,:,it] Regular format: image index is contained on first axis [it,:,:]

decim(im, d, av_en=1, fft=1)

Decimate image to lower resolution.

pairwise_distances(obs_pos)

Compute pairwise distances.

transpose_decim(im, decim_fact, av_en=0)

Apply the transpose of the decimation matrix.

SoftThresholding(data, thresh)

Perform element-wise soft thresholding.

HardThresholding(data, thresh)

Perform element-wise hard thresholding.

kthresholding(x, k)

Apply k-thresholding.

Keep only k highest values and set the rest to 0.

lineskthresholding(mat, k)

Apply k-thresholding to each line of input matrix.

Calls:

• utils.kthresholding()

mad(x, weight=None)

Compute MAD (Median Absolute Deviation).

transform_mask(weights, filt)

Propagate bad pixels to 1st wavelet scale and mask all pixels affected.

Bad pixels are the ones with weight 0.

lanczos(U, n=10, n2=None)

Generate Lanczos kernel for a given shift.

flux_estimate(im, cent=None, sigma=4)

Estimate flux for one image.

Parameters

• im (numpy.ndarray) – Image stamp containing the star.

• cent (numpy.ndarray) – Centroid of the star. If not provided, the centroid is calculated. Default is None.

• sigma (float) – Size of the star in sigma that will be used to calculate the flux and possibly the centroid too. Default is 4.

Returns

flux – Photometric flux value of the star.

Return type

float

Notes

See SPRITE paper (Ngole et al.), section 3.4.1., subsection ‘Photometric flux’.

flux_estimate_stack(stack, cent=None, sigmas=2.0)

Estimate flux for a bunch of images.

Calls:

• utils.flux_estimate()

shift_ker_stack(shifts, upfact, lanc_rad=8)

Generate shifting kernels and rotated shifting kernels.

Calls:

• utils.lanczos()

gen_Pea(distances, e, a)

Compute the graph Laplacian for a given set of parameters.

Parameters

• distances (numpy.ndarray) – Array of pairwise distances

• e (float) – Exponent to which the pairwise distances should be raised.

• a (float) – Constant multiplier along Laplacian’s diagonal.

Returns

Pea – Graph laplacian.

Return type

numpy.ndarray

Notes

Computes \(P_{e,a}\) matrix for given e, a couple. See Equations (16-17) in RCA paper (Ngole et al.). Watch out with the e parameter as it plays a vital role in the graph definition as it is a parameter of the distance that defines the graph’s weights.

select_vstar(eigenvects, R, weights)

Pick best eigenvector from a set of \((e,a)\).

i.e., solve (35) from RCA paper (Ngole et al.).

Parameters

• eigenvects (numpy.ndarray) – Array of eigenvects to be tested over.

• R (numpy.ndarray) – \(R_i\) matrix.

• weights (numpy.ndarray) – Entry-wise weights for \(R_i\).

class GraphBuilder(obs_data, obs_pos, obs_weights, n_comp, n_eigenvects=None, n_iter=3, ea_gridsize=10, distances=None, auto_run=True, verbose=2)

Bases: object

GraphBuilder class.

This class computes the necessary quantities for RCA’s graph constraint.

Parameters

• obs_data (numpy.ndarray) – Observed data.

• obs_pos (numpy.ndarray) – Corresponding positions.

• obs_weights (numpy.ndarray) – Corresponding per-pixel weights.

• n_comp (int) – Number of RCA components.

• n_eigenvects (int) – Maximum number of eigenvectors to consider per \((e,a)\) couple. Default is None; if not provided, all eigenvectors will be considered, which can lead to a poor selection of graphs, especially when data is undersampled. Ignored if VT and alpha are provided.

• n_iter (int) – How many alternations should there be when optimizing over \(e\) and \(a\). Default is 3.

• ea_gridsize (int) – How fine should the logscale grid of \((e,a)\) values be. Default is 10.

• distances (numpy.ndarray) – Pairwise distances for all positions. Default is None; if not provided, will be computed from given positions.

• auto_run (bool) – Whether to immediately build the graph quantities. Default is True.

pick_emax(epsilon=1e-15)

Pick maximum value for e parameter.

From now, we fix the maximum \(e\) to 1 and ignore the old procedure that was giving values that were too big.

Old procedure: Select maximum value of \(e\) for the greedy search over set of \((e,a)\) couples, so that the graph is still fully connected.

select_params(R, e_range, a_range)

Select best graph parameters.

Select \((e,a)\) parameters and best eigenvector for current \(R_i\) matrix.

Parameters

• R (numpy.ndarray) – Current \(R_i\) matrix (as defined in RCA paper (Ngole et al.), sect. 5.5.3.)

• e_range (numpy.ndarray) – List of \(e\) values to be tested.

• a_range (numpy.ndarray) – List of \(a\) values to be tested.

gen_eigenvects(mat)

Compute input matrix’s eigenvectors.

Keep only the n_eigenvects associated with the smallest eigenvalues.

poly_pos(pos, max_degree, center_normalice=True, x_lims=None, y_lims=None, normalice_Pi=True, min_degree=None)

Construct polynomial matrix.

Return a matrix Pi containing polynomials of stars positions up to max_degree.

Defaulting to CFIS CCD limits.

New method: The positions are scaled to the [-0.5, 0.5]x[-0.5, 0.5]. Then the polynomials are constructed with the normalized positions.

Old method: Positions are centred, the polynomials are constructed. Then the polynomials are normalized.

class CentroidEstimator(im, sig=7.5, n_iter=5, auto_run=True, xc=None, yc=None)

Bases: object

Estimate intra-pixel shifts.

It calculates the centroid of the image and compare it with the stamp centroid and returns the proper shift. The star centroid is calculated following an iterative procedure where a matched elliptical gaussian is used to calculate the moments.

Parameters

• im (numpy.ndarray) – Star image stamp.

• sig (float) – Estimated shape of the star in sigma. Default is 7.5.

• n_iter (int) – Max iteration number for the iterative estimation procedure. Default is 5.

• auto_run (bool) – Auto run the intra-pixel shif calculation in the initialization of the class. Default is True.

• xc (float) – First guess of the x component of the star centroid. (optional) Default is None.

• yc (float) – First guess of the y component of the star centroid. (optional) Default is None.

UpdateGrid()

Update the grid where the star stamp is defined.

EllipticalGaussian(e1=0, e2=0)

Compute an elliptical 2D gaussian with arbitrary centroid.

ComputeMoments()

Compute the star moments.

Compute the star image normalized first order moments with the current window function.

estimate()

Estimate the star image centroid iteratively.

return_shifts()

Return intra-pixel shifts.

Intra-pixel shifts are the difference between the estimated centroid and the center of the stamp (or pixel grid).

adjoint_degradation_op(x_i, shift_ker, D)

Apply adjoint of the degradation operator degradation_op.

degradation_op(X, shift_ker, D)

Shift and decimate fine-grid image.

handle_SExtractor_mask(stars, thresh)

Handle Sextractor masks.

Reads SExtracted star stamps, generates MCCD-compatible masks (that is, binary weights), and replaces bad pixels with 0s - they will not be used by MCCD, but the ridiculous numerical values can otherwise still lead to problems because of convolutions.

match_psfs(test_stars, PSFs)

Match psfs.DEPRECATED.

See MCCD.validation_stars instead. Takes as input the test_stars vignets and the PSFs vignets that were outputs from the psf modelling method. The function outputs the PSFs matching the corresponding test stars. This allows to compute the pixel RMSE. Intended to be used with PSFEx validation functions.

Parameters

• test_stars (numpy.ndarray) – reg format (n_stars,n_pix,n_pix)

• PSFs (numpy.ndarray) – reg format (n_stars,n_pix,n_pix)

Returns

deg_PSFs – reg format (n_stars,n_pix,n_pix)

Return type

numpy.ndarray

class NoiseEstimator(img_dim, win_rad)

Bases: object

Noise estimator.

Parameters

• img_dim (tuple of int) – Image size

• win_rad (int) – window radius in pixels

static sigma_mad(x)

Compute an estimation of the standard deviation of a Gaussian distribution using the robust MAD (Median Absolute Deviation) estimator.

estimate_noise(image)

Estimate the noise level of the image.