NAME: BOXAVE PURPOSE: Box-average a 1 or 2 dimensional array. This procedure differs from the intrinsic REBIN function in the follow 2 ways: (1) the box size parameter is specified rather that the output array size (2) for INTEGER arrays, BOXAVE computes intermediate steps using REAL*4 arithmetic. This is considerably slower than REBIN but avoids integer truncation CALLING SEQUENCE: result = BOXAVE( Array, Xsize,[ Ysize ] ) INPUTS: ARRAY - Two dimensional input Array to be box-averaged. Array may be one or 2 dimensions and of any type except character. OPTIONAL INPUTS: XSIZE - Size of box in the X direction, over which the array is to be averaged. If omitted, program will prompt for this parameter. YSIZE - For 2 dimensional arrays, the box size in the Y direction. If omitted, then the box size in the X and Y directions are assumed to be equal OUTPUT: RESULT - Output array after box averaging. If the input array has dimensions XDIM by YDIM, then RESULT has dimensions XDIM/NBOX by YDIM/NBOX. The type of RESULT is the same as the input array. However, the averaging is always computed using REAL arithmetic, so that the calculation should be exact. If the box size did not exactly divide the input array, then then not all of the input array will be boxaveraged. PROCEDURE: BOXAVE boxaverages all points simultaneously using vector subscripting REVISION HISTORY: Written, W. Landsman, October 1986 Call REBIN for REAL*4 and REAL*8 input arrays, W. Landsman Jan, 1992
(See /usr/local/idl/lib/zastron/image/boxave.pro)
NAME: CONVOLVE PURPOSE: Convolution of image with Point Spread Function (PSF), default is to compute using product of Fourier transforms. CALLING SEQUENCE: imconv = convolve( image1, psf, FT_PSF = psf_FT ) or correl = convolve( image1, image2, /CORREL ) or correl = convolve( image, /AUTO ) INPUTS: image = the image to be convolved with psf psf = the Point Spread Function, (size < or = to size of image). OPTIONAL INPUT KEYWORDS: FT_PSF = passes the Fourier transform of the PSF, so that it can be re-used for the next time function is called. /ROTATE : convolution with rotated PSF, same as /CORREL if using FFT. /CORRELATE uses the conjugate of the Fourier transform of PSF, to compute the cross-correlation of image and PSF. /AUTO_CORR computes the auto-correlation function of image. /NO_FT overrides the use of FFT, using the IDL function convol() instead. (then PSF is rotated by 180 degrees to give equivalent result) METHOD: When using FFT, PSF is centered & expanded to size of image. HISTORY: written 1992 Frank Varosi STX @ NASA/GSFC
(See /usr/local/idl/lib/zastron/image/convolve.pro)
NAME: CORREL_IMAGES PURPOSE: Computes the 2-D cross-correlation function of two images for a range of (x,y) shifting by pixels of one image relative to the other. CALLING SEQUENCE: Result = correl_images( image_A, image_B ) INPUTS: image_A, image_B = the two images of interest. KEYWORDS: XSHIFT = the + & - shift to be applied in X direction, default=7. YSHIFT = the Y direction + & - shifting, default=7. XOFFSET_B = initial X pixel offset of image_B relative to image_A. YOFFSET_B = Y pixel offset, defaults are (0,0). REDUCTION = optional reduction factor causes computation of Low resolution correlation of bin averaged images, thus faster. Can be used to get approximate optimal (x,y) offset of images, and then called for successive lower reductions in conjunction with CorrMat_Analyze until REDUCTION=1, getting offset up to single pixel. MAGNIFICATION = option causes computation of high resolution correlation of magnified images, thus much slower. Shifting distance is automatically = 2 + Magnification, and optimal pixel offset should be known and specified. Optimal offset can then be found to fractional pixels using CorrMat_Analyze( correl_images( ) ). /NUMPIX causes the number of pixels for each correlation to be saved in a second image, concatenated to the correlation image, so Result is fltarr( Nx, Ny, 2 ). /MONITOR causes the progress of computation to be briefly printed. OUTPUTS: Result is the cross-correlation function, given as a matrix. PROCEDURES USED: nint - Nearest Integer function PROCEDURE: Loop over all possible (x,y) shifts, compute overlap and correlation for each shift. Correlation set to zero when there is no overlap. MODIFICATION HISTORY: Written, July,1991, Frank Varosi, STX @ NASA/GSFC
(See /usr/local/idl/lib/zastron/image/correl_images.pro)
NAME: CORREL_OPTIMIZE PURPOSE: Find the optimal (x,y) pixel offset of image_B relative to image_A, by means of maximizing the correlation function. CALLING SEQUENCE: correl_optimize, image_A, image_B, xoffset_optimum, yoffset_optimum INPUTS: image_A, image_B = the two images of interest. KEYWORDS: XOFF_INIT = initial X pixel offset of image_B relative to image_A, YOFF_INIT = Y pixel offset, (default offsets are 0 and 0). MAGNIFICATION = option to determine offsets up to fractional pixels, (example: MAG=2 means 1/2 pixel accuracy, default=1). /NUMPIX: sqrt( sqrt( # pixels )) used as correlation weighting factor. /MONITOR causes the progress of computation to be briefly printed. /PRINT causes the results of analysis to be printed. OUTPUTS: xoffset_optimum = optimal X pixel offset of image_B relative to image_A. yoffset_optimum = optimal Y pixel offset. CALLS: function correl_images( image_A, image_B ) pro corrmat_analyze PROCEDURE: The combination of function correl_images( image_A, image_B ) and corrmat_analyze of the result is used to obtain the (x,y) offset yielding maximal correlation. The combination is first executed at large REDUCTION factors to speed up computation, then zooming in recursively on the optimal (x,y) offset by factors of 2. Finally, the MAGNIFICATION option (if specified) is executed to determine the (x,y) offset up to fractional pixels. MODIFICATION HISTORY: Written, July,1991, Frank Varosi, STX @ NASA/GSFC
(See /usr/local/idl/lib/zastron/image/correl_optimize.pro)
NAME: CORRMAT_ANALYZE PURPOSE: Analyze the 2-D cross-correlation function of two images and find the optimal(x,y) pixel offsets. Intended for use with function CORREL_IMAGES. CALLING SEQUENCE: corrmat_analyze, correl_mat, xoffset_optimum, yoffset_optimum, [...] INPUTS: correl_mat = the cross-correlation matrix of 2 images. (as computed by function CORREL_IMAGES( imA, imB ) ). NOTE: If correl_mat(*,*,1) is the number of pixels for each correlation, (the case when /NUMPIX was specified in call to CORREL_IMAGES) then sqrt( sqrt( # pixels )) is used as correlation weighting factor. KEYWORDS: XOFF_INIT = initial X pixel offset of image_B relative to image_A. YOFF_INIT = Y pixel offset, (both as specified to correl_images). REDUCTION = reduction factor used in call to CORREL_IMAGES. MAGNIFICATION = magnification factor used in call to CORREL_IMAGES, this allows determination of offsets up to fractions of a pixel. PLATEAU_THRESH = threshold used for detecting plateaus in the cross-correlation matrix near maximum, (default=0.01), used only if MAGNIFICATION > 1 /PRINT causes the result of analysis to be printed. OUTPUTS: xoffset_optimum = optimal X pixel offset of image_B relative to image_A. yoffset_optimum = optimal Y pixel offset. max_corr = the maximal correlation corresponding to optimal offset. edge = 1 if maximum is at edge of correlation domain, otherwise=0. plateau = 1 if maximum is in a plateua of correlation function, else=0. FUNCTION CALLS: function nint( number ) nearest integer function PROCEDURE: Find point of maximum cross-correlation and calc. corresponding offsets. If MAGNIFICATION > 1: the correl_mat is checked for plateau near maximum, and if found, the center of plateau is taken as point of maximum cross-correlation. MODIFICATION HISTORY: Written, July-1991, Frank Varosi, STX @ NASA/GSFC
(See /usr/local/idl/lib/zastron/image/corrmat_analyze.pro)
NAME: DIST_CIRCLE PURPOSE: Form a square array in which the value of each element is its distance to a specified center. Useful for circular aperture photometry. CALLING SEQUENCE: DIST_CIRCLE, IM, N, XCEN, YCEN INPUTS: N = either a scalar specifying the size of the N x N square output array, or a 2 element vector specifying the size of the N x M rectangular output array. XCEN,YCEN = Scalars designating the X,Y pixel center. These need not be integers, and need not be located within the output image OUTPUTS: IM - N by N (or M x N) floating array in which the value of each pixel is equal to its distance to XCEN,YCEN EXAMPLE: Total the flux in a circular aperture within 3' of a specified RA and DEC on an 512 x 512 image IM, with a header H. IDL> adxy, H, RA, DEC, x, y ;Convert RA and DEC to X,Y IDL> getrot, H, rot, cdelt ;CDELT gives plate scale deg/pixel IDL> cdelt = cdelt*3600. ;Convert to arc sec/pixel IDL> dist_circle, circle, 512, x, y ;Create a distance circle image IDL> circle = circle*abs(cdelt(0)) ;Distances now given in arcseconds IDL> good = where(circle LT 180) ;Within 3 arc minutes IDL> print,total( IM(good) ) ;Total pixel values within 3' RESTRICTIONS: The speed of DIST_CIRCLE decreases and the the demands on virtual increase as the square of the output dimensions. Users should dimension the output array as small as possible, and re-use the array rather than re-calling DIST_CIRCLE MODIFICATION HISTORY: Adapted from DIST W. Landsman March 1991 Allow a rectangular output array W. Landsman June 1994
(See /usr/local/idl/lib/zastron/image/dist_circle.pro)
NAME:
DIST_ELLIPSE
PURPOSE:
Form an array in which the value of each element is equal to the
semi-major axis of the ellipse of specified center, axial ratio, and
position angle, which passes through that element. Useful for
elliptical aperture photometry.
CALLING SEQUENCE:
DIST_ELLIPSE, IM, N, XC, YC, RATIO, POS_ANG
INPUTS:
N = either a scalar specifying the size of the N x N square output
array, or a 2 element vector specifying the size of the
M x N rectangular output array.
XC,YC - Scalars giving the position of the ellipse center. This does
not necessarily have to be within the image
RATIO - Scalar giving the ratio of the major to minor axis. This
should be greater than 1 for postion angle to have its
standard meaning.
OPTIONAL INPUTS:
POS_ANG - Position angle of the major axis, measured counter-clockwise
from the Y axis. For an image in standard orientation
(North up, East left) this is the astronomical position angle.
OUTPUT:
IM - REAL*4 elliptical mask array, of size M x N. THe value of each
pixel is equal to the semi-major axis of the ellipse of center
XC,YC, axial ratio RATIO, and position angle POS_ANG, which
passes through the pixel.
EXAMPLE:
Total the flux in a elliptical aperture with a major axis of 3', an
axial ratio of 2.3, and a position angle of 25 degrees centered on
a specified RA and DEC. The image array, IM is 200 x 200, and has
an associated FITS header H.
ADXY, H, ra, dec, x, y ;Get X and Y corresponding to RA and Dec
GETROT, H, rot, cdelt ;CDELT gives plate scale degrees/pixel
cdelt = abs( cdelt)*3600. ;CDELT now in arc seconds/pixel
DIST_ELLIPSE, ell, 200, x, y, 2.3, 25 ;Create a elliptical image mask
ell = ell*cdelt(0) ;Distances now given in arcseconds
good = where( ell lt 180 ) ;Within 3 arc minutes
print,total( im(good) ) ;Total pixel values within 3'
RESTRICTIONS:
The speed of DIST_ELLIPSE decreases and the the demands on virtual
increase as the square of the output dimensions. Users should
dimension the output array as small as possible, and re-use the
array rather than re-calling DIST_ELLIPSE
REVISION HISTORY:
Written W. Landsman April, 1991
Somewhat faster algorithm August, 1992
Allow rectangular output array June, 1994
(See /usr/local/idl/lib/zastron/image/dist_ellipse.pro)
NAME:
FILTER_IMAGE
PURPOSE:
Computes the average and/or median of pixels in moving box,
replacing center pixel with the computed average and/or median,
(using the IDL smooth or median functions).
The main reason for using this function is the option to
also process the pixels at edges and corners of image, and,
to apply iterative smoothing simulating convolution with gaussian.
CALLING EXAMPLE:
Result = filter_image( image, SMOOTH=box_width, /MEDIAN, /ALL )
INPUTS:
image = 2-D image (matrix)
KEYWORDS:
MEDIAN = width of square moving box for median, in # pixels.
/MEDIAN means use box width= 3 pixels for median filter.
SMOOTH = width of square moving box for averaging, in # pixels.
/SMOOTH means use box width= 3 pixels for smoothing.
/ALL_PIXELS causes the edges of image to be filtered as well.
/ITERATE means apply smooth(*,3) iteratively for (box_width-1)/2 times,
approaching convolution with a Gaussian FWHM=sqrt( box_width ).
Note that /ALL_PIXELS is automatically applied,
giving better results in the iteration limit.
(also, MEDIAN is ignored when /ITER is specified).
RESULT:
Smoothed or median filtered image is result of function.
If both SMOOTH and MEDIAN are specified, median filter is applied first.
PROCEDURE:
Create a larger image by reflecting the edges outward,
then call the IDL median and/or smooth function on the larger image,
and just return the central part (the orginal size image).
MODIFICATION HISTORY:
Written, Sep.1991, Frank Varosi, STX @ NASA/GSFC
Mod, Nov.1991, Frank Varosi, added /ITERATE option.
(See /usr/local/idl/lib/zastron/image/filter_image.pro)
NAME IMLIST PURPOSE: Display pixel values on an image surrounding a specified X,Y center. IMLIST is similar to TVLIST but the center pixel is supplied directly by the user, rather than being read off of the image display CALLING SEQUENCE: IMLIST, Image, Xc, Yc, [ TEXTOUT = , DX = , DY = ,WIDTH = ,DESCRIP = ] INPUTS: Image - Two-dimensional array containing the image Xc - X pixel value at which to center the display, integer scalar Yc - Y pixel value at which to center the display, integer scalar OPTIONAL INPUTS KEYWORDS: TEXTOUT - Scalar number (1-5) or string which determines output device. (see TEXTOPEN) The following dev/file is opened for output. textout=1 TERMINAL using /more option textout=2 TERMINAL without /more option textout=3.prt textout=4 laser.tmp textout=5 user must open file textout = filename (default extension of .prt) DX -Integer scalar giving the number of pixels inthe X direction to be displayed. If omitted then DX = 18 for byte images, and DX = 14 for integer images. IMLIST will display REAL data with more significant figures if more room is availble to print. DY - Same as DX, but in Y direction. If omitted, then DY = DX WIDTH - Integer scalar giving the character width of the output device. Default is 80 characters. DESCRIP = Scalar string which will be written as a description over the output pixel values. If DESCRIP is not supplied, and the output device specified by TEXTOUT is not a terminal, then the user will be prompted for a description. OUTPUTS: None. PROCEDURE: Corresponding region of image is then displayed at the terminal. If necessary, IMLIST will divide all pixel values in a REAL*4 image by a (displayed) factor of 10 to make a pretty format. SYSTEM VARIABLES: If the keyword TEXTOUT is not supplied, then the non-standard system variable !TEXTOUT will be read. RESTRICTIONS: IMLIST may not be able to correctly format all pixel values if the dynamic range of the values near the center pixel is very large EXAMPLE: Display the pixel values of an image array IM in the vicinity of 254,111 IDL> imlist, IM, 254, 111 PROCEDURES USED TEXTOPEN, F_FORMAT, TEXTCLOSE, DATATYPE REVISION HISTORY: Written, W. Landsman June, 1991 Added DESCRIP keyword W. Landsman December, 1991
(See /usr/local/idl/lib/zastron/image/imlist.pro)
NAME:
MAX_ENTROPY
PURPOSE:
Deconvolution of data by Maximum Entropy analysis, given the
instrument point spread response function (spatially invariant psf).
Data can be an observed image or spectrum, result is always positive.
Default is convolutions using FFT (faster when image size = power of 2).
CALLING SEQUENCE:
for i = 1,Niter do begin
Max_Entropy, data, psf, deconv, Lmult, FT_PSF = psf_ft
INPUTS:
data = observed image or spectrum, should be mostly positive,
with mean sky (background) near zero.
psf = Point Spread Function of instrument (response to point source,
must sum to unity).
deconv = result of previous call to Max_Entropy,
multipliers = the Lagrange multipliers of max.entropy theory
(on first call, set = 0, giving flat first result).
OUTPUTS:
deconv = deconvolution result of one more iteration by Max_Entropy.
multipliers = the Lagrange multipliers saved for next iteration.
KEYWORDS:
FT_PSF = passes (out/in) the Fourier transform of the PSF,
so that it can be reused for the next time procedure is called,
/NO_FT overrides the use of FFT, using the IDL function convol() instead.
/LINEAR switches to Linear convergence mode, much slower than the
default Logarithmic convergence mode.
LOGMIN = minimum value constraint for taking Logarithms (default=1.e-9).
EXTERNAL CALLS:
function convolve( image, psf ) for convolutions using FFT or otherwise.
METHOD:
Iteration with PSF to maximize entropy of solution image with
constraint that solution convolved with PSF gives data image.
Based on paper by Hollis, Dorband, Yusef-Zadeh, Ap.J. Feb.1992,
which refers to Agmon, Alhassid, Levine, J.Comp.Phys. 1979.
HISTORY:
written Feb.1992 Frank Varosi STX @ NASA/GSFC
(See /usr/local/idl/lib/zastron/image/max_entropy.pro)
NAME:
MAX_LIKELIHOOD
PURPOSE:
Deconvolution of an observed image (or spectrum), given the
instrument point spread response function (spatially invariant psf).
Performs iteration based on the Maximum Likelihood solution for
the restoration of a blurred image (or spectrum) with additive noise.
Maximum Likelihood formulation can assume Poisson noise statistics
or Gaussian additive noise, yielding two types of iteration.
CALLING SEQUENCE:
for i=1,Niter do begin
Max_Likelihood, data, psf, deconv, FT_PSF=psf_ft
INPUTS:
data = observed image or spectrum, should be mostly positive,
with mean sky (background) near zero.
psf = Point Spread Function of the instrument,
(response to a point source, must sum to unity).
deconv = result of previous call to Max_Likelihood,
(on first call, initial guess default = average of data).
OUTPUTS:
deconv = deconvolution result of one more iteration by Max_Likelihood.
KEYWORDS:
/GAUSSIAN causes max-likelihood iteration for Gaussian additive noise
to be used, otherwise the default is Poisson statistics.
CONSTRAINT = used only when /GAUSS specified, and then the image is
constrained to be > CONSTRAINT, otherwise image > 0.
FT_PSF = passes (out/in) the Fourier transform of the PSF,
so that it can be reused for the next time procedure is called,
/NO_FT overrides the use of FFT, using the IDL function convol() instead.
EXTERNAL CALLS:
function convolve( image, psf ) for convolutions using FFT or otherwise.
METHOD:
Maximum Likelihood solution is a fixed point of an iterative eq.
(derived by setting partial derivatives of Likelihood to zero).
Poisson noise case was derived by Richardson(1972) & Lucy(1974).
Gaussian noise case is similar with subtraction instead of division.
HISTORY:
written 1992 Frank Varosi STX @ NASA/GSFC
(See /usr/local/idl/lib/zastron/image/max_likelihood.pro)
NAME: psf_Gaussian PURPOSE: Return a point spread function having Gaussian profiles, as either a 1D vector, a 2D image, or 3D volumetric-data. CALLING: psf = psf_Gaussian( NPIXEL=, FWHM= , [/NORMALIZE, /ST_DEV, ) or: psf = psf_Gaussian( parameters, NPIXEL = ) REQUIRED KEYWORDS: NPIXEL = number pixels for each dimension, specify as an array, or just one number to make all sizes equal. OPTIONAL KEYWORDS: NDIMEN = dimension of result: 1 (vector), 2 (image), or 3 (volume), default = 2 (an image result). FWHM = the desired Full-Width Half-Max (pixels) in each dimension, specify as an array, or single number to make all the same. CENTROID = pixels numbers of PSF maximum ( 0.5 is center of a pixel ), default is exact center of requested vector/image/volume. STDEV = optional way to specify width by standard deviation param. XY_CORREL = scalar between 0 and 1 specifying correlation coefficient Use this keyword, for example, to specify an elliptical gaussian oriented at an angle to the X,Y axis /NORMALIZE causes resulting PSF to be normalized so Total( psf ) = 1. INPUTS (optional): parameters = an NDIMEN by 3 array giving for each dimension: [ maxval, center, stdev ], overrides other keywords. EXAMPLE: Create a 31 x 31 array containing a normalized centered gaussian with an X FWHM = 4.3 and a Y FWHM = 3.6 IDL> array = PSF_GAUSSIAN( Npixel=31, FWHM=[4.3,3.6], /NORMAL EXTERNAL CALLS: function Gaussian HISTORY: Written, Frank Varosi NASA/GSFC 1991.
(See /usr/local/idl/lib/zastron/image/psf_gaussian.pro)
NAME: RINTER PURPOSE: Cubic interpolation of an image at a set of reference points points. Optionally obtain the X and Y derivatives at these points. Used by the DAOPHOT PSF sequence but can also be for general cubic interpolation. CALLING SEQUENCE: Z = RINTER( P, X, Y, [ DFDX, DFDY ] ) INPUTS: P - Two dimensional data array, X - Either an N element vector or an N x M element array, containing X subscripts where cubic interpolation is desired. Y - Either an N element vector or an N x M element arra, containing Y subscripts where cubic interpolation is desired. OUTPUT: Z - Result = interpolated vector or array. If X and Y are vectors, then so is Z, but if X and Y are arrays then Z will be also. If P is DOUBLE precision, then so is Z, otherwise Z is REAL. OPTIONAL OUTPUT: DFDX - Vector or Array, (same size and type as Z), containing the derivatives with respect to X DFDY - Array containing derivatives with respect to Y EXAMPLE: suppose P is a 256 x 256 element array and X = FINDGEN(50)/2. + 100. and Y = X. Then Z will be a 50 element array, containing the cubic interpolated points. SIDE EFFECTS: can be time consuming. RESTRICTION: Interpolation is not possible at positions outside the range of the array (including all negative subscripts), or within 2 pixel units of the edge. No error message is given but values of the output array are meaningless at these positions. PROCEDURE: invokes CUBIC interpolation algorithm to evaluate each element in Z at virtual coordinates contained in X and Y with the data in P. COMMON BLOCKS: If repeated interpolation of the same array is to occur, then one can save time by initializing the common block RINTER REVISION HISTORY: March 1988 written W. Landsman STX Co. Checked for IDL Version 2, J. Isensee, September, 1990 Corrected call to HISTOGRAM, W. Landsman November 1990
(See /usr/local/idl/lib/zastron/image/rinter.pro)
NAME:
sigma_filter
PURPOSE:
Computes the mean and standard deviation of pixels in a box centered at
each pixel of the image, but excluding the center pixel. If the center
pixel value exceeds some # of standard deviations from the mean, it is
replaced by the mean in box. Note option to process pixels on the edges.
CALLING SEQUENCE:
Result = sigma_filter( image, box_width, N_sigma=(#), /ALL,/MON )
INPUTS:
image = 2-D image (matrix)
box_width = width of square filter box, in # pixels (default = 3)
KEYWORDS:
N_sigma = # standard deviations to define outliers, floating point,
recommend > 2, default = 3. For gaussian statistics:
N_sigma = 1 smooths 35% of pixels, 2 = 5%, 3 = 1%.
RADIUS = alternative to specify box radius, so box_width = 2*radius+1.
/ALL_PIXELS causes computation to include edges of image,
/KEEP causes opposite effect: pixels with values outside of specified
deviation are not changed, pixels within deviation are smoothed.
/ITERATE causes sigma_filter to be applied recursively (max = 20 times)
until no more pixels change (only allowed when N_sigma >= 2).
/MONITOR prints information about % pixels replaced.
Optional Outputs:
N_CHANGE = # of pixels changed (replaced with neighborhood mean).
VARIANCE = image of pixel neighborhood variances * (N_sigma)^2,
DEVIATION = image of pixel deviations from neighborhood means, squared.
CALLS:
function filter_image( )
PROCEDURE:
Compute mean over moving box-cars using smooth, subtract center values,
compute variance using smooth on deviations from mean,
check where pixel deviation from mean is within variance of box,
replace those pixels in smoothed image (mean) with orignal values,
return the resulting partial mean image.
MODIFICATION HISTORY:
Written, 1991, Frank Varosi and Dan Gezari NASA/GSFC
F.V.1992, added optional keywords /ITER,/MON,VAR=,DEV=,N_CHANGE=.
(See /usr/local/idl/lib/zastron/image/sigma_filter.pro)