pyimcom.imdestripe

Program to remove correlated noise stripes from RST images.

Classes

Sca_img: Class defining an SCA image object
Parameters: Class holding the destriping parameters for a given mosaic
Cost_models: Class holding the cost function models. Currently only quadratic is supported

Functions

write_to_file: Function to write some text to an output file
save_fits: Save a 2D image to a FITS file with locking, retries, and atomic rename.
apply_object_mask: Apply a bright object mask to an image.
quadratic: Quadratic cost function f(x) = x^2
absolute: Absolute cost function f(x) = |x|
huber_loss: Huber loss cost function
quad_prime: Derivative of quadratic cost function f’(x) = 2x
abs_prime: Derivative of absolute cost function f’(x) = sign(x)
huber_prime: Derivative of Huber loss cost function
get_scas: Function to get a list of all SCA images and their WCSs for this mosaic
interpolate_image_bilinear: Interpolate values from a “reference” SCA image onto a “target” SCA coordinate grid
transpose_interpolate: Interpolate backwards from image_A to image_B space.
transpose_par: Sum up the values of an image across rows
get_effective_gain: retrieve the effective gain and n_eff of the image. valid only for already-interpolated images
get_ids: Take an SCA label and parse it out to get the Obsid and SCA id strings.
save_snapshot: Save restart state to pickle file.
get_neighbors: Get the neighboring SCAs for each SCA in the mosaic
residual_function: Calculate the residual image, = grad(epsilon)
residual_function_single: Helper function to calculate residuals for a single SCA
cost_function_single: Helper function to calculate cost for a single SCA
linear_search_general: Perform a linear search to find the optimal step size alpha along a given direction
linear_search_quadratic: Calculate optimal step size alpha along a given direction, for quadratic cost function
conjugate_gradient: Perform the conjugate gradient algorithm to minimize the cost function

Attributes

`testoutputs`
`filters`
`t0_global`
`use_output_float`
`tempdir`
`img_full_output`
`img_full_output`
`profiler`

Classes

`Cost_models`	Class holding the cost function models. This is a dictionary of functions
`Sca_img`	Class defining an SCA image object.
`Parameters`	Class holding the parameters for a given mosaic. This can be the destriping parameters, or additional

Functions

`write_to_file`(text[, filename])	Function to write some text to an output file
`save_fits`(image, filename[, dir, overwrite, s, ...])	Save a 2D image to a FITS file with locking, retries, and atomic rename.
`apply_object_mask`(image[, mask, threshold_m, ...])	Apply a bright object mask to an image.
`quadratic`(x)	Quadratic cost function f(x) = x^2
`absolute`(x)	Absolute cost function f(x) = \|x\|
`huber_loss`(x, d)	Huber loss cost function
`quad_prime`(x)	Derivative of quadratic cost function f'(x) = 2x
`abs_prime`(x)	Derivative of absolute cost function f'(x) = sign(x)
`huber_prime`(x, d)	Derivative of Huber loss cost function
`get_scas`(filter_, obsfile, cfg[, indata_type])	Function to get a list of all SCA images and their WCSs for this mosaic
`interpolate_image_bilinear`(image_B, image_A, ...[, mask])	Interpolate values from a "reference" SCA image onto a "target" SCA coordinate grid
`transpose_interpolate`(image_A, wcs_A, image_B, ...)	Interpolate backwards from image_A to image_B space.
`transpose_par`(img)	Sum up the values of an image across rows
`get_effective_gain`(sca[, tempdir])	Retrieve the effective gain and n_eff of the image. valid only for already-interpolated images
`get_ids`(sca)	Take an SCA label and parse it out to get the Obsid and SCA id strings.
`save_snapshot`(p, grad, epsilon, psi, direction, ...)	Save restart state to pickle file.
`get_neighbors`(scalist, ov_mat[, overlap_thresh])	Get a dictionary of overlapping SCAs for each SCA in the mosaic
`residual_function`(psi, f_prime, scalist, wcslist, ...)	Calculate the residual image, = grad(epsilon)
`residual_function_single`(k, sca_a, wcs_a, psi_a, ...)	Calculate the residual for a single SCA image
`cost_function_single`(j, sca_a, p, f, scalist, ...)	Calculate the cost function for a single SCA image
`cost_function`(p, f, thresh, workers, scalist, ...[, ...])	Calculate the cost function with the current de-striping parameters.
`linear_search_general`(p, direction, f, f_prime, ...[, ...])	Linear search via combination bisection and secant methods for parameters that minimize the function
`linear_search_quadratic`(p, direction, f, f_prime, ...)	For the quadratic cost function, direct calculation of alpha that minimizes the function
`conjugate_gradient`(p, f, f_prime, thresh, workers, ...)	Algorithm to use conjugate gradient descent to optimize the parameters for destriping.
`main`([cfg_file, overlaponly])	Main function to run destriping via conjugate gradient descent.

Module Contents

testoutputs[source]

filters = ['W146', 'F184', 'H158', 'J129', 'Y106', 'Z087', 'R062', 'PRSM', 'DARK', 'GRSM', 'K213'][source]

t0_global[source]

use_output_float[source]

tempdir = ''[source]

img_full_output[source]

img_full_output[source]

class Cost_models(cfg)[source]

Class holding the cost function models. This is a dictionary of functions

model[source]

class Sca_img(obsid, scaid, cfg, tempdir=tempdir, interpolated=False, add_objmask=True, indata_type='fits')[source]

Class defining an SCA image object.

Parameters:

scaid (Str) – the SCA id
obsid (Str) – the observation id
cfg (Config object) – built from the configuration file
interpolated (Bool) – True if you want the interpolated version of this SCA and not the original. Default False
add_objmask (Bool) – True if you want to apply the permanent pixel mask and a bright object mask. Default True

image

the SCA image (4088x4088)

Type:: 2D np array

shape

the shape of the image

Type:: Tuple

w

the astropy.wcs object associated with this SCA

Type:: WCS object

obsid[source]

observation ID of this SCA image

Type:: Str

scaid[source]

SCA ID (position on focal plane) of this SCA image

Type:: Str

mask[source]

The full pixel mask that is used on this image. Is correct only after applying masks to image

Type:: 2D np array

g_eff

Effective gain in each pixel of the image

Type:: 2D np array

params_subtracted[source]

True if parameters have been subtracted from this image.

Type:: Bool

cfg[source]

the configuration object passed in at initialization

Type:: Config object

apply_noise()[source]: Apply the appropriate lab noise frame to the SCA image

apply_permanent_mask()[source]: Apply the SCA permanent pixel mask to the image

apply_asdf_mask()[source]: Apply the SCA ASDF file mask to the image

apply_all_mask()[source]: Apply the full SCA mask to the image

subtract_parameters()[source]: Subtract a given set of parameters from self.image; updates self.image, self.params_subtracted

get_coordinates()[source]: Create an array of ra, dec coords for the image

make_interpolated()[source]: Construct a version of this SCA interpolated from other, overlapping ones. Writes the interpolated image out to the disk, to be read/used later

obsid[source]

scaid[source]

mask[source]

params_subtracted = False[source]

cfg[source]

apply_noise()[source]: Add detector noise to self.image

apply_permanent_mask()[source]: Apply permanent pixel mask. Updates self.image and self.mask

apply_asdf_mask()[source]: Apply ASDF file mask. Updates self.image and self.mask

get_permanent_mask()[source]: Apply permanent pixel mask. Updates self.image and self.mask

apply_all_mask()[source]: Apply permanent pixel mask. Updates self.image in-place

subtract_parameters(p, j)[source]

Subtract a set of parameters from the SCA image. Updates self.image and self.params_subtracted

Parameters:

p (Parameters object) – containing params of current iteration
j (int) – the index of the SCA image into all_scas list

get_coordinates(pad=0.0)[source]

Create an array of ra, dec coords for the image

Parameters:: pad (Float64) – N pixels of padding to add to the array. Default 0.0
Returns:: 1D arrays of ra, dec coordinates for each pixel in the image
Return type:: ra, dec; 1D np.arrays of length (height*width)

make_interpolated(ind, scalist, neighbors, tempdir=tempdir, params=None, N_eff_min=0.5)[source]

Construct a version of this SCA interpolated from other, overlapping ones. Writes the interpolated image out to the disk, to be read/used later The N_eff_min parameter requires some minimum effective coverage, otherwise masks that pixel.

Parameters:

ind (int) – index of this SCA in all_scas list
scalist (List of Str) – the list of all SCAs in this mosaic
neighbors (Dict) – dictionary where keys are SCA indices and values are lists of indices of overlapping SCAs
tempdir (Str) – directory to write temporary files to
params (Parameters object) – parameters to be subtracted from contributing SCAs; default None
N_eff_min (float) – Effective coverage needed for a pixel to contribute to the interpolation

class Parameters(cfg, scalist=[])[source]

Class holding the parameters for a given mosaic. This can be the destriping parameters, or additional parameters that need to be the same shape and have the same methods

Parameters:

cfg (Config object) – built from the configuration file
scalist (list of Strings) – the list of SCAs in this mosaic

model[source]

Which destriping model to use, which then specifies the number of parameters per row. Must be a key of the model_params dict

Type:: Str

n_rows[source]

Number of rows in the image, or number of rows to fit ds model over

Type:: Int

params_per_row[source]

Number of parameters per row, set by model_params[model]

Type:: Int

params[source]

The actual array of parameters.

Type:: 2D np array

current_shape[source]

The current shape (1D or 2D) of SCA params

Type:: Str

scalist[source]

the list of SCAs in this mosaic, format: filter_obsid_scaid

Type:: list of Strings

params_2_images()[source]: Reshape params into a 2D array, with one row per SCA

forward_par()[source]: Reshape one row of params array (one SCA) into a 2D array by projection along rows

model[source]

n_rows[source]

params_per_row[source]

params[source]

current_shape = '2D'[source]

scalist = [][source]

params_2_images()[source]: Reshape flattened parameters into 2D array with 1 row per sca and n_rows (in image) * params_per_row entries

forward_par(sca_i)[source]

Takes one SCA row (n_rows) from the params and casts it into 2D (n_rows x n_rows)

Parameters:: sca_i (Int) – Index of which SCA to recast into 2D
Return type:: 2D np.array, the image of SCA_i’s parameters

write_to_file(text, filename='destripe_out.txt')[source]

Function to write some text to an output file

Parameters:

text (Str) – The text to print
filename (Str) – Filename to write out to. Default ‘destripe_out.txt’

save_fits(image, filename, dir=None, overwrite=True, s=False, header=None, retries=3)[source]

Save a 2D image to a FITS file with locking, retries, and atomic rename.

Parameters:

image (np.ndarray) – 2D array to write.
filename (str) – Output filename without extension.
dir (str) – Directory to save into.
overwrite (bool) – Whether to overwrite the final target file.
s (bool) – Whether to print status messages.
header (fits.Header or None) – Optional FITS header.
retries (int) – Number of write retry attempts if write fails.

apply_object_mask(image, mask=None, threshold_m=0, threshold_c=0.3, inplace=False)[source]

Apply a bright object mask to an image.

Parameters:

image (2D numpy array) – the image to be masked.
mask (2D boolean array, optional) – the pre-existing object mask. Default: None
threshold_m (float) – factor to multiply with the median for thresholding.
threshold_c (float) – constant to add to the threshold.
inplace (Bool) – Whether to modify the input image directly.

Returns:

image_out (2D np.array) – the masked image.
neighbor_mask (2D np.array) – the mask applied

quadratic(x)[source]: Quadratic cost function f(x) = x^2

absolute(x)[source]: Absolute cost function f(x) = |x|

huber_loss(x, d)[source]: Huber loss cost function

quad_prime(x)[source]: Derivative of quadratic cost function f’(x) = 2x

abs_prime(x)[source]: Derivative of absolute cost function f’(x) = sign(x)

huber_prime(x, d)[source]: Derivative of Huber loss cost function

get_scas(filter_, obsfile, cfg, indata_type='fits')[source]

Function to get a list of all SCA images and their WCSs for this mosaic

Parameters:

filter (Str) – which filter to use for this run. Options: Y106, J129, H158, F184, K213
obsfile (Str) – prefix / path to the SCA images
cfg (Config object) – built from the configuration file
indata_type (Str) – input data type: ‘fits’ or ‘asdf’. Default ‘fits’

Returns:

all_scas (list of strings) – list of all the SCAs in this mosiac
all_wcs (list of WCS objects) – the WCS object for each SCA in all_scas (same order)

interpolate_image_bilinear(image_B, image_A, interpolated_image, mask=None)[source]

Interpolate values from a “reference” SCA image onto a “target” SCA coordinate grid Uses pyimcom_croutines.bilinear_interpolation(

float* image, float* g_eff, int rows, int cols, float* coords, int num_coords, float* interpolated_image)

Parameters:

image_B (SCA object) – the image to be interpolated
image_A (SCA object) – the image whose grid you are interpolating B onto
interpolated_image (2D np array) – all zeros with shape of Image A. Updated in place to be the interpolation of img. B onto A’s grid
mask (2D np array, optional) – if provided, this mask is interpolated instead of image_B.image

transpose_interpolate(image_A, wcs_A, image_B, original_image)[source]

Interpolate backwards from image_A to image_B space. Uses bilinear_transpose(

float* image, int rows, int cols, float* coords, int num_coords, float* original_image)

image_A2D np array
the already-interpolated gradient image

wcs_Awcs.WCS object
image A’s WCS object

image_BSCA object
the image we’re interpolating the gradient back onto

original_image2D np array
the gradient image re-interpolated into image B space Updated in place

transpose_par(img)[source]

Sum up the values of an image across rows

Parameters:

img (2D np.array) –
Input array

Returns
-------- – 1D np.array, the sum across each row of I

get_effective_gain(sca, tempdir=tempdir)[source]

Retrieve the effective gain and n_eff of the image. valid only for already-interpolated images

Parameters:

sca (Str) – format like “<prefix>_<obsid>_<scaid>” describing which SCA to get the effective gain for
tempdir (Str) – directory where the effective gain files are stored

Returns:

g_eff (memmap 2D np.array) – the effective gain in each pixel
N_eff (memmap 2D np.array) – how many image “B”s contributed to that interpolated image

get_ids(sca)[source]

Take an SCA label and parse it out to get the Obsid and SCA id strings.

Parameters:

sca (Str) – The sca name from all_scas list, formatted like <obsid>_<scaid>

Returns:

obsid (Str) – the observation ID
scaid (Str) – the SCA ID (position in focal plane)

save_snapshot(p, grad, epsilon, psi, direction, grad_prev, direction_prev, cg_model, tol, thresh, norm_0, cost_model, i, restart_file)[source]

Save restart state to pickle file.

Parameters:

p (Parameters object) – current parameters
grad (2D np array) – current gradient
epsilon (2D np array) – current cost
psi (3D np array) – current residuals
direction (2D np array) – current CG direction
grad_prev (2D np array) – previous gradient
direction_prev (2D np array) – previous CG direction
cg_model (Str) – which CG model is being used
tol (Float) – tolerance for convergence
thresh (Float) – threshold for Huber loss cost function
norm_0 (Float) – initial norm of the gradient
cost_model (Str) – which cost function is being used
i (Int) – current iteration number
restart_file (Str) – path to the restart pickle file

get_neighbors(scalist, ov_mat, overlap_thresh=0.1)[source]

Get a dictionary of overlapping SCAs for each SCA in the mosaic

Parameters:

scalist (List of Str) – the list of all SCAs in this mosaic
ov_mat (2D np array) – the overlap matrix for all SCAs in this mosaic
overlap_thresh (Float) – minimum overlap fraction to consider two SCAs as neighbors; default 0.1

Returns:

neighbors – dictionary where keys are SCA indices and values are lists of indices of overlapping SCAs

Return type:

Dict

residual_function(psi, f_prime, scalist, wcslist, neighbors, thresh, workers, cfg, extrareturn=False)[source]

Calculate the residual image, = grad(epsilon)

Parameters:

psi (3D np array) – the image difference array (I_A - J_A) (N_SCA, 4088, 4088)
f_prime (Function) – the derivative of the cost function f in the future this should be set by default based on what you pass for f
scalist (List of Str) – the list of all SCAs in this mosaic
wcslist (List of WCS objects) – the WCS object for each SCA in scalist (same order)
neighbors (Dict) – dictionary where keys are SCA indices and values are lists of indices of overlapping SCAs
thresh (Float) – threshold for Huber loss cost function; default None
workers (Int) – number of parallel workers to use
cfg (Config object) – the configuration for this run
extrareturn (Bool) –

if True, return residual terms 1 and 2 separately; Default False
in addition to full residuals. returns resids, resids1, resids2

Returns:

resids – with one row per SCA and one col per parameter

Return type:

2D np array

residual_function_single(k, sca_a, wcs_a, psi_a, f_prime, scalist, neighbors, thresh, cfg)[source]

Calculate the residual for a single SCA image

Parameters:

k (Int) – index of this SCA in scalist
sca_a (Str) – the SCA label, formatted like <obsid>_<scaid>
wcs_a (wcs.WCS object) – the WCS object for this SCA
psi_a (2D np array) – the difference image I_A - J_A for this SCA
f_prime (Function) – the derivative of the cost function f
scalist (List of Str) – the list of all SCAs in this mosaic
neighbors (Dict) – dictionary where keys are SCA indices and values are lists of indices of overlapping SCAs
thresh (Float) – threshold for Huber loss cost function; default None
cfg (Config object) – the configuration for this run

Returns:

k (Int) – index of this SCA in scalist
term_1 (1D np array) – the first residual term for this SCA
term_2_list (List of tuples) – list of (j, term_2) tuples containing value for term 2 for each SCA j that overlaps with this one

cost_function_single(j, sca_a, p, f, scalist, neighbors, thresh, cfg)[source]

Calculate the cost function for a single SCA image

Parameters:

j (Int) – index of this SCA in scalist
sca_a (Str) – the SCA label, formatted like <obsid>_<scaid>
p (Parameters object) – the current parameters for de-striping
f (Function) – the cost function form
scalist (List of Str) – the list of all SCAs in this mosaic
neighbors (Dict) – dictionary where keys are SCA indices and values are lists of indices of overlapping SCAs
thresh (Float) – threshold for Huber loss cost function; default None
cfg (Config object) – the configuration for this run

Returns:

j (Int) – index of this SCA in scalist
psi (2D np array) – the difference image I_A - J_A for this SCA
local_epsilon (Float) – the cost function value for this SCA

cost_function(p, f, thresh, workers, scalist, neighbors, cfg, tempdir=tempdir)[source]

Calculate the cost function with the current de-striping parameters.

Parameters:

p (parameters object) – the current parameters for de-striping
f (st) – keyword for function dictionary options; should also set an f_prime
thresh (Float) – threshold for Huber loss cost function; default None
workers (Int) – number of parallel workers to use
scalist (List of Str) – the list of all SCAs in this mosaic
neighbors (Dict) – dictionary where keys are SCA indices and values are lists of indices of threshold-overlapping SCAs
cfg (Config object) – the configuration for this run
tempdir (Str) – directory to store temporary files

Returns:

epsilon (int) – the total cost function summed over all images
psi (D np array) – the difference images I_A-J_A

linear_search_general(p, direction, f, f_prime, cost_model, epsilon_current, psi_current, grad_current, thresh, workers, scalist, wcslist, neighbors, cfg, n_iter=100, tol=10**-4)[source]

Linear search via combination bisection and secant methods for parameters that minimize the function: d_epsilon/d_alpha in the given direction . Note alpha = depth of step in direction

Parameters:

p (params object) – the current de-striping parameters
direction (2D np array) – direction of conjugate gradient search
f (function) – cost function form
f_prime (function) – derivative of cost function form
cost_model (Str) – which cost function is being used; options: ‘quadratic’, ‘huber_loss’
epsilon_current (float) – current cost function value
psi_current (3D np array) – current difference images (I_A - J_A)
grad_current (2D np array) – current gradient AKA current residuals
thresh (float) – threshold for Huber loss cost function; default None
workers (Int) – number of parallel workers to use
scalist (List of Str) – the list of all SCAs in this mosaic
wcslist (List of WCS objects) – the WCS object for each SCA in scalist (same order)
neighbors (Dict) – dictionary where keys are SCA indices and values are lists of indices of overlapping SCAs
cfg (Config object) – the configuration for this run
n_iter (int) – number of iterations at which to stop searching
tol (float) – absolute value of d_cost at which to converge

Returns:

best_p (parameters object) – containing the best parameters found via search
best_psi (3D numpy array) – the difference images made from images with the best_p params subtracted off

linear_search_quadratic(p, direction, f, f_prime, grad_current, thresh, workers, scalist, wcslist, neighbors, cfg)[source]

For the quadratic cost function, direct calculation of alpha that minimizes the function: d_epsilon/d_alpha in the given direction . Note alpha = depth of step in direction

Finds the best alpha, computes the new parameters and diff image, and prints the new cost and convergence criteria

Parameters:

p (params object) – the current de-striping parameters
direction (2D np array) – direction of conjugate gradient search
f (function) – cost function form
f_prime (function) – derivative of cost function form
grad_current (2D np array) – current gradient AKA current residuals
thresh (float or None) – threshold for Huber loss cost function
workers (Int) – number of parallel workers to use
scalist (List of Str) – the list of all SCAs in this mosaic
wcslist (List of WCS objects) – the WCS object for each SCA in scalist (same order)
neighbors (Dict) – dictionary where keys are SCA indices and values are lists of indices of overlapping SCAs
cfg (Config object) – the configuration for this run

Returns:

new_p (parameters object) – containing the new parameters found via direct calculation
new_psi (3D numpy array) – the difference images made from images with the new_p params subtracted off
new_resids (2D np array) – the new residuals calculated with new_p
new_epsilon (float) – the new cost function value calculated with new_p

conjugate_gradient(p, f, f_prime, thresh, workers, scalist, wcslist, neighbors, restart_file=None, time_limit=None, cfg=None, of='destripe_out.txt')[source]

Algorithm to use conjugate gradient descent to optimize the parameters for destriping. Direction is updated using Fletcher-Reeves method

Parameters:

p (parameters object) – containing initial parameters guess
f (function) – functional form to use for cost function
f_prime (function) – the derivative of f. KL: eventually f should dictate f prime
thresh (float or None) – threshold for Huber loss cost function; default None
workers (Int) – number of parallel workers to use
scalist (List of Str) – the list of all SCAs in this mosaic
wcslist (List of WCS objects) – the WCS object for each SCA in scalist (same order)
neighbors (Dict) – dictionary where keys are SCA indices and values are lists of indices of overlapping SCAs
restart_file (Str or None) – if not None, path to pickle file containing restart state
time_limit (int or None) – if not None, how much time to elapse before stopping (minutes)
cfg (config object) – containing all config parameters
of (Str) – output file to write log messages to

Returns:

p – the best fit parameters for destriping the SCA images

Return type:

params object

main(cfg_file=None, overlaponly=False)[source]

Main function to run destriping via conjugate gradient descent.

Parameters:

cfg_file (str, optional) – Configuration file (if not provided, reads from command line arguments).
overlaponly (bool, optional) – Only compute the overlap matrix, then stop.

Returns:

Prefix for destriped images. Add f”_{obsid}_{sca}.fits” to get the full file name.

Return type:

str

profiler[source]