pro Imgcorf, im1, im2, xshift, yshift, cr, $
             nofit=nofit, verbose=verbose, display=display, $
             badpix=badpix 

; use  cross correlation to find the offset between two images
;
; INPUTS:
;   im1, im2	images (must be same size)
;
; OUTPUTS:
;   xshift	x shift between the two images (see below)
;   yshift	y shift between the two images
;   crosscr	cross correlation image between the two frames
;
; KEYWORD INPUTS:
;   /nofit	don't use Goddard CNTRD to find shift, just
;		use max pixel
;   /verbose	print results
;   /display	use windows 0-3 to display results
;   badpix      bad pixel map
;
; NOTES:
; - the shift convention is defined so that for an object
;   with positions x1 and x2 in the images, the result is
;   xshift = x2 - x1
; - unlike James' 'offcor' routine, there is no accomodation
;   for the DC offset result though an error flag is printed
;   if the resulting shift is 0.  
;
; USES:
;   crossc, wrap, ndshift, cntrd
;
; HISTORY:
;   Written by MCL: 7/30/95 (based on 'offcor' routine by JRG)
;


; fwhm parameter used in finding centroid
FITFWHM = 5.0

if n_params() lt 4 then begin
	print,'pro imgcorf,im1,im2,xshift,yshift,(crosscor),'
	print,'    [nofit],[display],[verbose],[badpix=]'
	retall
endif


sz = size(im1)
if total(size(im1)-size(im2)) ne 0.0 then $
	message,'images are not same size!'

if keyword_set(badpix) then begin
    fixpix, im1, badpix, im1fix
    fixpix, im2, badpix, im2fix
endif else begin
    im1fix =  im1
    im2fix =  im2
endelse

; display if desired
if keyword_set(display) then begin
	win,0,title='IDL (0): image 1'
	display2,im1fix,/asp
	win,1,title='IDL (1): image 2'
	display2,im2fix,/asp
endif


; find cross correlation using Fourier product
cr = crossc(im1fix,im2fix)


; find peak in cross correlation, being wary of multiple peaks
found =  0
repeat begin
    w = where(cr eq max(cr), nmax)
    whereis, cr, w, xmax, ymax
    if nmax gt 1 then begin
        message, 'multiple maxima in cross correlation', /info
        xmax = round(total(xmax)/n_elements(xmax)) 
        ymax = round(total(ymax)/n_elements(ymax)) 
    endif 
    
;   if found DC component, set it to zero and continue
    xshift = sz(1)/2 - xmax
    yshift = sz(2)/2 - ymax
    if (abs(xshift) le 0.5) and (abs(yshift) le 0.5) then begin
        message, 'found DC component - masking and trying again', $
          /info
        map = bytarr(sz(1), sz(2)) + 1B
        map(xmax, ymax) =  0
        fixpix,  cr, map, cr, /silent

    endif else $
      found =  1

endrep until (found eq 1)


; try to find shift to fractional pixel if desired
if keyword_set(nofit) then begin
	xcen=xmax
	ycen=ymax
endif else begin
	cntrd,cr,xmax,ymax,xcen,ycen,FITFWHM
	if xcen eq -1 then begin
		xcen = xmax
		ycen = ymax
	endif
endelse


; display cross correlation image and region close to peak
if keyword_set(display) then begin
	win,2,title='IDL (2): cross corr'
;	tvh,cr,/silent
	display2,cr,/asp
	xc = round(xcen)
	yc = round(ycen)
;	tvcircle,10,xc,yc,/data
	win,3,title='IDL (3): cross corr peak'
	x0 = (xc-15) > 0
	x1 = (xc+15) < (sz(1)-1)
	y0 = (yc-15) > 0
	y1 = (yc+15) < (sz(2)-1)
	display2,cr(x0:x1,y0:y1), $
		findgen(x1-x0+1)+x0,findgen(y1-y0+1)+y0, $
		/asp
	oplot,[xcen],[ycen],psym=1,color=0,thick=1
endif


; I think this is ok if image has an odd number of pixels on a side
xshift = sz(1)/2 - xcen
yshift = sz(2)/2 - ycen


; output
if keyword_set(verbose) then begin
	print,xmax,ymax,xcen,ycen,xshift,yshift
endif

if (abs(xshift) lt 1) and (abs(yshift) lt 1) then $
  message,'No shift found - recovered DC component only',/info




end