;+
; NAME:
;	correl_optimize
; PURPOSE:
;	Assuming the images are two detections of the same scene,
;	find the optimal (x,y) pixel offset of image_B relative to image_A,
;	by means of maximizing the cross-correlation function.
;	Note that bad pixels or extremely noisy data will bias the
;	cross-correlation function and give possibly erroneous results.
;	Small amount of noise is no problem, pre-filtering is recommended,
;	returned variable max_correl can be used to decide validity.
;	Cross-correlation is computed directly by default, but if
;	keyword /FFT is set it will be computed using FFTs of images.
; CALLING:
;	correl_optimize, image_A, image_B, xoffset_optimum, yoffset_optimum 
; INPUTS:
;	image_A, image_B = the two images to be analyzed, the offset of
;			image_B relative to image_A is determined.
; 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).
;	MAX_REDUCTION = maximum reduction factor allowed, otherwise
;		the first reduction rebins image into about 8 x 8 pixels.
;	MAGNIFICATION = option to determine offsets up to fractional pixels,
;			(example: MAG=2 means 1/2 pixel accuracy, default=1).
;	/NUMPIX : (# overlap pixels)^(1/4) used as correlation weighting factor.
;	/MONITOR causes the progress of computation to be briefly printed.
;	/PRINT causes the results of analysis to be printed to standard output.
;	/FFT : use FFT to compute correlation of images, usually faster,
;		but may not be always accurate because periodic nature of FT
;		may cause confusion (in this case /NUMPIX has no effect).
;		But if data is well inside image boundaries it is accurate.
; OUTPUTS:
;	xoffset_optimum = optimal X pixel offset of image_B relative to image_A.
;	yoffset_optimum = optimal Y pixel offset.
;	max_correl = 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 plateau of correlation function, else=0.
; EXTERNAL CALLS:
;	function  correl_images( image_A, image_B )
;	function  convolve( image_A, image_B,/CORREL )	if /FFT is set.
;	pro  corrmat_analyze
;	IDL system variable !DEBUG should be defined.
; PROCEDURE:
;	The combination of function correl_images( image_A, image_B ) and
;	pro corrmat_analyze is used to obtain the (x,y) offset
;	yielding maximal correlation. The combination is first executed at
;	large REDUCTION factor (to speed up computation) and with
;	all possible shifts to find the approximate optimal (x,y) offsets.
;	The reduction factor is then decreased by factors of 2 and
;	the x & y shifts are kept small, thus zooming into the optimal offsets.
;	Finally, the MAGNIFICATION option (if specified)
;	is executed to determine the (x,y) offset up to fractional pixels.
;	If /FFT is set then just use function convolve( image_A,image_B,/CORR )
;	and pro corrmat_analyze to obtain the optimal (x,y) offset,
;	rebinning the images first if magnification is greater than one.
; NOTES:
;	Recommend that the images be normalized so that maximim=1,
;	and thresholded at some uniform percent of max=1.
;	See for example function modify_image.
; HISTORY:
;	Written: Frank Varosi, NASA/GSFC, 1991.
;-

pro correl_optimize, image_A, image_B, xoffset_optimum, yoffset_optimum, $
					max_correl, edge, plateau,       $
					XOFF_INIT = xoff_init, $
					YOFF_INIT = yoff_init, $
					PRINT=print, MONITOR=monitor, $
					NUMPIX=numpix, MAGNIFICATION=Magf, $
					MAX_REDUCTION = maxreducf, FFT=use_FFT
	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
		return
	   endif

	if N_elements( xoff_init ) NE 1 then xoff_init=0
	if N_elements( yoff_init ) NE 1 then yoff_init=0
	xoff = xoff_init
	yoff = yoff_init
	if N_elements( Magf ) NE 1 then Magf=1

;This is option to use FFT and then return:

	if keyword_set( use_FFT ) then begin

	   xoff = ( simA(1) - simB(1) )/2.
	   yoff = ( simA(2) - simB(2) )/2.
	   simA = Magf * simA
	   simB = Magf * simB
	   norm = 1/( Magf^2 * sqrt( total( image_A^2 ) * total( image_B^2 ) ) )

	   if N_elements( image_B ) GT N_elements( image_A ) then begin

		corrmat_analyze, convolve( rebin( image_B, simB(1), simB(2) ), $
					   rebin( image_A, simA(1), simA(2) ), $
					/CORREL ) * norm, MAG=Magf, $
					xoffset_optimum, yoffset_optimum, $
					max_correl, edge, plateau, PR=print, $
					XOFF=-xoff, YOFF=-yoff

		xoffset_optimum = -xoffset_optimum
		yoffset_optimum = -yoffset_optimum

	    endif else begin

		corrmat_analyze, convolve( rebin( image_A, simA(1), simA(2) ), $
					   rebin( image_B, simB(1), simB(2) ), $
					/CORREL ) * norm, MAG=Magf, $
					xoffset_optimum, yoffset_optimum, $
					max_correl, edge, plateau, PR=print, $
					XOFF=xoff, YOFF=yoff
	     endelse

	   return
	 endif

;This is normal method of direct correlations:

	reducf = min( [simA(1:2),simB(1:2)] ) / 8	;Bin average to about
							; 8 by 8 pixel image.
	reducf = 2^round_off( aLog( reducf )/aLog(2) )	; make it power of 2.
	if N_elements( maxreducf ) EQ 1 then reducf = reducf < maxreducf

	xsiz = simA(1) > simB(1)
	ysiz = simA(2) > simB(2)
	xshift = xsiz
	yshift = ysiz		;shift over the whole images first correlation.

	while (reducf GT 1) do begin

		corrmat_analyze, correl_images( image_A, image_B, NUM=numpix, $
						XOFF=xoff, YOFF=yoff, $
						XS=xshift, YS=yshift, $
					       REDUC=reducf, MONIT=monitor ), $
			xoff, yoff, XOFF=xoff, YOFF=yoff, PR=print, REDUC=reducf

		xshift = 2*reducf
		yshift = 2*reducf	;shift over coarse pixel grid to refine
		reducf = reducf/2	; correlations at increasing resolution.
	  endwhile

	if keyword_set( monitor ) then print,"Magnification = 1"

	corrmat_analyze, correl_images( image_A, image_B, XSHIFT=3, YSHIFT=3, $
				XOFF=xoff,YOFF=yoff, MON=monitor,NUM=numpix ),$
				xoffset_optimum, yoffset_optimum, max_correl, $
				edge, plateau, XOFF=xoff, YOFF=yoff, PR=print

	if (Magf GE 2) then begin

		xoff = xoffset_optimum		;refine offsets to
		yoff = yoffset_optimum		; fractional pixels.

		corrmat_analyze, correl_images( image_A, image_B, XOFF=xoff, $
					YOFF=yoff, MAG=Magf, MON=monitor ),  $
					xoffset_optimum, yoffset_optimum, $
					max_correl, edge, plateau,  $
					XOFF=xoff,YOFF=yoff,PR=print, MAG=Magf
	   endif
end