# 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 NIRSPEC Data Analysis (if verbose options are used in dia() or dialin()) numpy, matplotlib, sys 2008-11-19 08:54 IJC: Created @ UCLA 2012-04-02 14:26 IJMC: Updated documentation; code unchanged.
dia.centroid(a)[source]

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)[source]

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 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 (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. ```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') ``` 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)[source]
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 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 (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. 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)[source]

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 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) (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. I – inverse matrix 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 dsamulti()
dia.dsamulti(mod, obs, N, **kw)[source]

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 w=None – L x M array containing the statistical weights of ‘obs’ meansub=False – subtract a weighted mean from obs (along axis=1) 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 dsa()
dia.rconvolve1d(a, b, mode='valid', extend='nearest')[source]

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 This is “reversed” from the canonical definition of the convolution. dsa()
dia.rconvolve2d(a, b, mode='valid', extend=0, verbose=False)[source]

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 This is “reversed” from the canonical definition of the convolution. Note that I could also implement an FFT convolution algorithm, for increased speed.

#### Previous topic

“Kernel density estimate” statistics

#### Next topic

Fourier imaging routines