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 linearlyvarying kernel 

REQUIREMENTS:  Numerical analysis routines (for dsamulti()) NIRSPEC Data Analysis (if verbose options are used in dia() or dialin()) numpy, matplotlib, sys 
HISTORY:  20081119 08:54 IJC: Created @ UCLA 20120402 14:26 IJMC: Updated documentation; code unchanged. 
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
Computational tool for Difference Image Analysis (DIA)
INPUTS: 
I – Current image to be analysed.


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

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 20081118 11:12 IJC: Extrapolating from my 1D algorithm 20120403 08:19 IJMC: Added example to documentation. 20130425 15:28 IJMC: Optimized a bit for speed; 2x boost, but still slow. 
INPUTS: 
I – Current image to be analysed.


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

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 20081119 21:50 IJC: Trying it out 
Computational tool for Difference Spectral Analysis (DSA)
INPUTS: 
I – Current spectrum to be analysed. Nk – number of pixels in the desired convolution kernel 

OPTIONS: 
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

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: 
Run LSD using a given set of spectra and a model deltafunction spectrum.
INPUTS: 
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 varianceuncertainties on ‘lsdbox’ chisq – length M array with the chisquared values from the LSD 
SEE_ALSO: 
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: 
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: 