ISOPHOT Interactive Analysis (PIA): Astrolib

IMAGE


List of Routines


Routine Descriptions

BOXAVE

[Next Routine] [List of Routines]
 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)


CONVOLVE

[Previous Routine] [Next Routine] [List of Routines]
 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)


CORREL_IMAGES

[Previous Routine] [Next Routine] [List of Routines]
 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)


CORREL_OPTIMIZE

[Previous Routine] [Next Routine] [List of Routines]
 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)


CORRMAT_ANALYZE

[Previous Routine] [Next Routine] [List of Routines]
 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)


DIST_CIRCLE

[Previous Routine] [Next Routine] [List of Routines]
 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)


DIST_ELLIPSE

[Previous Routine] [Next Routine] [List of Routines]
 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)


FILTER_IMAGE

[Previous Routine] [Next Routine] [List of Routines]
 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)


IMLIST

[Previous Routine] [Next Routine] [List of Routines]
 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)


MAX_ENTROPY

[Previous Routine] [Next Routine] [List of Routines]
 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)


MAX_LIKELIHOOD

[Previous Routine] [Next Routine] [List of Routines]
 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)


PSF_GAUSSIAN

[Previous Routine] [Next Routine] [List of Routines]
 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)


RINTER

[Previous Routine] [Next Routine] [List of Routines]
 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)


SIGMA_FILTER

[Previous Routine] [List of Routines]
 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)