Difference Image Analysis

Contents:

Python scripts to perform 1D and 2D deconvolution.

CONTENTS:

dsa() – 1D: Difference Spectral Analysis

dsamulti() – 1D: DSA for sets of multiple spectra

dia() – 2D: Difference Image Analysis

dialin() – 2D DIA, with a linearly-varying kernel

REQUIREMENTS:

Numerical analysis routines (for dsamulti())

NIRSPEC Data Analysis (if verbose options are used in dia() or dialin())

numpy, matplotlib, sys

HISTORY:

2008-11-19 08:54 IJC: Created @ UCLA

2012-04-02 14:26 IJMC: Updated documentation; code unchanged.

dia.centroid(a)

Return the pixel coordinates of the centroid of an array.

(xc, yc) = centroid(a)

Uses the algorithm from: http://www.dfanning.com/tips/centroid.html

dia.dia(r, i, k, w=None, verbose=False, noback=False)

Computational tool for Difference Image Analysis (DIA)

INPUTS:

R – reference image. This should have the highest possible

signal-to-noise and the sharpest PSF.

I – Current image to be analysed.

k – 2D kernel basis mask: 1 for pixels to be used, 0 for

pixels to be ignored

OPTIONS:

w – weights of the pixel values in I; typically (sigma)^-2

noback – do not fit for a variable background; assume constant.

verbose – Print output statements and make a plot or two

OUTPUTS:

(M, K, B, C): M – R, convolved to match I

K – kernel used in convolution

B – background offset

C – chisquared of fit. If no weights were specified, weights

are set to unity for this calculation.

EXAMPLE:

import pylab as py
import dia

# Define reference and blurred images:
ref = py.zeros((10,10))
img = py.zeros((10,10))
ref[3,3] = 1; ref[3,6] = 1;  ref[6,3] = 1
img[2,3] = 1; img[3,2:5] = 1;  img[4,3] = 1
img[2,6] = 1; img[3,5:8] = 1;  img[4,6] = 1
img[5,3] = 1; img[6,2:5] = 1;  img[7,3] = 1

# Add some noise:
img += py.randn(10, 10)/10.

# Define kernel basis
kb = py.ones((3,3))    

# Compute Difference Image Analysis:
m, kern, bkg, chisq = dia.dia(ref, img, kb)

# Display results:
py.figure()
py.subplot(231)
py.imshow(ref)
py.title('Reference image')
py.subplot(232)
py.imshow(img)
py.title('Observed image')
py.subplot(234)
py.imshow(kern)
py.title('Optimal kernel')
py.subplot(235)
py.imshow(m)
py.title('Convolved Reference')
py.subplot(236)
py.imshow(img - m)
py.title('Residuals')

NOTES:

Best results are obtained with proper registration of the images.

Also, beware of edge effects. As a general rule, anything within a kernel width of the edges is suspect.

Based on the 2D Bramich (2008) DIA algorithm

2008-11-18 11:12 IJC: Extrapolating from my 1D algorithm

2012-04-03 08:19 IJMC: Added example to documentation.

2013-04-25 15:28 IJMC: Optimized a bit for speed; 2x boost, but still slow.

dia.dialin(r, i, k, w=None, verbose=True, noback=False)

Computational tool for Difference Image Analysis (DIA) with a

linearly-varying kernel

INPUTS:

R – reference image. This should have the highest possible

signal-to-noise and the sharpest PSF.

I – Current image to be analysed.

k – pixel mask for kernel; 1 for pixels to be used, 0 for

pixels to be ignored

OPTIONS:

w – weights of the pixel values in I; typically (sigma)^-2

noback – do not fit for a variable background; assume constant.

verbose – Print output statements and make a plot or two

OUTPUTS:

(M, K, X, Y, B, C): M – R, convolved to match I

K – kernel used in convolution

B – background offset

C – chisquared of fit. If no weights were specified, weights

are set to unity for this calculation.

NOTES:

Best results are obtained with proper registration of the images.

Also, beware of edge effects. As a general rule, anything within a kernel width of the edges is suspect.

Originally based on the 2D Bramich (2008) DIA algorithm

2008-11-19 21:50 IJC: Trying it out

dia.dsa(r, i, Nk, **kw)

Computational tool for Difference Spectral Analysis (DSA)

INPUTS:

R – reference spectrum. This should have the highest possible

signal-to-noise and the highest spectral resolution.

I – Current spectrum to be analysed.

Nk – number of pixels in the desired convolution kernel

OPTIONS:

w – weights of the pixel values in I; typically (sigma)^-2

(Not HANDELED CORRECTLY?!?!?)

noback – do not fit for a variable background; assume constant.

tol=1e-10 – if matrix determinant is less than tol, use

pseudoinverse rather than straight matrix inversion

verbose – Print output statements and make a plot or two

retinv – return a fourth output, the Least Squares inverse

matrix (False by default)

OUTPUTS:

(M, K, B, C): M – R, convolved to match I

K – kernel used in convolution

B – background offset

C – chisquared of fit. If no weights were specified, weights

are set to unity for this calculation.

OPTIONS:

I – inverse matrix

NOTES:

Best results are obtained with proper registration of the spectra.

Also, beware of edge effects. As a general rule, anything within a kernel width of the edges is suspect.

Also

SEE_ALSO:

dsamulti()

dia.dsamulti(mod, obs, N, **kw)

Run LSD using a given set of spectra and a model delta-function spectrum.

INPUTS:

mod – length L array containing a model delta-functon spectrum,

made with a model linelist and nsdata.linespec

obs – L X M array containing M observations of spectra of length L

N – size of deconvolution kernel to solve for

OPTIONS:

w=None – L x M array containing the statistical weights of ‘obs’

meansub=False – subtract a weighted mean from obs (along axis=1)

OUTPUT:

lsdbox – M x N array of computer deconvolution kernels

varbox – M x N array of statistical variance-uncertainties on ‘lsdbox’

chisq – length M array with the chi-squared values from the LSD

SEE_ALSO:

dsa()

dia.rconvolve1d(a, b, mode='valid', extend='nearest')

Compute a 1D reverse convolution in the style of Bramich 2008.

INPUTS:

‘a’ should be longer than ‘b’ – i.e., ‘b’ is the kernel.

‘extend’ tells how to extend the boundaries – either ‘nearest’-neighbor or a number

NOTES:

This is “reversed” from the canonical definition of the convolution.

SEE_ALSO:

dsa()

dia.rconvolve2d(a, b, mode='valid', extend=0, verbose=False)

Compute a 2D reverse convolution in the style of Bramich 2008.

INPUTS:

‘a’ should be larger than ‘b’ – i.e., ‘b’ is the kernel – and both should be square.

‘extend’ tells how to extend the boundaries – either ‘nearest’-neighbor or a number

NOTES:

This is “reversed” from the canonical definition of the convolution.

Note that I could also implement an FFT convolution algorithm, for increased speed.

SEE_ALSO:

dsa(), rconvolve1d()