function cor_imag, image_A, image_B, XSHIFT = x_shift,	$
					  YSHIFT = y_shift, 	$
					  XOFFSET_B = x_offset, $
					  YOFFSET_B = y_offset, $
					  REDUCTION = reducf,	$
					  MAGNIFICATION = Magf, $
					  NUMPIX=numpix, MONITOR=monitor
;+
; NAME:
;	cor_imag
; 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 = cor_imag( 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 cor_anal
;			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 cor_anal ( cor_imag( ) ).
;
;	/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
;-
	simA = size( image_A )
	simB = size( image_B )

	if (simA(0) LT 2) OR (simB(0) LT 2) then begin
		message,"first two arguments must be images",/INFO,/CONTIN
		return,[-1]
	   endif

	if N_elements( x_offset ) NE 1 then x_offset=0
	if N_elements( y_offset ) NE 1 then y_offset=0

	if N_elements( x_shift ) NE 1 then x_shift = 7
	if N_elements( y_shift ) NE 1 then y_shift = 7
	x_shift = abs( x_shift )
	y_shift = abs( y_shift )

	if keyword_set( reducf ) then begin

		reducf = fix( reducf ) > 1
		if keyword_set( monitor ) then $
				print,"Reduction = ",strtrim( reducf, 2 )
		simA = simA/reducf
		LA = simA * reducf -1	;may have to drop edges of images.
		simB = simB/reducf
		LB = simB * reducf -1

		imtmp_A = Rebin( image_A( 0:LA(1), 0:LA(2) ), simA(1), simA(2) )
		imtmp_B = Rebin( image_B( 0:LB(1), 0:LB(2) ), simB(1), simB(2) )

		xoff = nint ( x_offset/reducf )
		yoff = nint ( y_offset/reducf )
		xs = x_shift/reducf
		ys = y_shift/reducf

		return, cor_imag( imtmp_A, imtmp_B, XS=xs,YS=ys,$
							XOFF=xoff, YOFF=yoff, $
						MONITOR=monitor, NUMPIX=numpix )

	  endif else if keyword_set( Magf ) then begin

		Magf = fix( Magf ) > 1
		if keyword_set( monitor ) then $
				print,"Magnification = ",strtrim( Magf, 2 )
		simA = simA*Magf
		simB = simB*Magf

		imtmp_A = rebin( image_A, simA(1), simA(2), /SAMPLE )
		imtmp_B = rebin( image_B, simB(1), simB(2), /SAMPLE )

		xoff = nint( x_offset*Magf )
		yoff = nint( y_offset*Magf )

		return, cor_imag( imtmp_A, imtmp_B, XS=Magf+2, YS=Magf+2,$
							XOFF=xoff, YOFF=yoff, $
						MONITOR=monitor, NUMPIX=numpix)
	   endif

	Nx = 2 * x_shift + 1
	Ny = 2 * y_shift + 1
	if keyword_set( numpix ) then Nim=2 else Nim=1

	correl_mat = fltarr( Nx, Ny, Nim )

	xs = nint( x_offset ) - x_shift
	ys = nint( y_offset ) - y_shift

	sAx = simA(1)-1
	sAy = simA(2)-1
	sBx = simB(1)-1
	sBy = simB(2)-1

	for y = 0, Ny-1 do begin	;compute correlation for each y,x shift.

	    yoff = ys + y
	    yAmin = yoff > 0
	    yAmax = sAy < (sBy + yoff)
	    yBmin = (-yoff) > 0
	    yBmax = sBy < (sAy - yoff)		;Y overlap

	    if (yAmax GT yAmin) then begin

	       for x = 0, Nx-1 do begin

		   xoff = xs + x
		   xAmin = xoff > 0
		   xAmax = sAx < (sBx + xoff)
		   xBmin = (-xoff) > 0
		   xBmax = sBx < (sAx - xoff)		;X overlap

		   if (xAmax GT xAmin) then begin

			im_ov_A = image_A( xAmin:xAmax, yAmin:yAmax )
			im_ov_B = image_B( xBmin:xBmax, yBmin:yBmax )
			Npix = N_elements( im_ov_A )

			if N_elements( im_ov_B ) NE Npix then begin
			       message,"overlap error: # pixels NE",/INFO,/CONT
			       print, Npix, N_elements( im_ov_B )
			       stop
			   endif
			
			im_ov_A = im_ov_A - total( im_ov_A )/Npix
			im_ov_B = im_ov_B - total( im_ov_B )/Npix
			totAA = total( im_ov_A * im_ov_A )
			totBB = total( im_ov_B * im_ov_B )

			correl_mat(x,y) = total( im_ov_A * im_ov_B ) / $
							sqrt( totAA * totBB )

			if keyword_set( numpix ) then correl_mat(x,y,1) = Npix
		     endif

	          endfor
		endif

		if keyword_set( monitor ) then print, Ny-y, FORM="($,i3)"
	  endfor

	if keyword_set( monitor ) then print," "

return, correl_mat
end