;+
; NAME:
;	Shifter_2d
; PURPOSE:
; 	Compute X and Y shifts required to coalign
;	a set of IVM raw images
; CATEGORY:
; 	Data analysis
; CALLING SEQUENCE:
;	shifter_2d,file,nn,nimages,shifts, vb
; INPUTS:
; 	NN is the size (=nx=ny) of the images
;	NIMAGES is the number of images in the file
;	vb is verbosity (0=none, 1=start/end, 2=shifts)
; OPTIONAL INPUT PARAMETERS:
; KEYWORD PARAMETERS:
; 	geod - geometry data array - fltarr(nn,nn,nimages)
;	geof - if set, geometry file from which to read 'geod'
; 	BX, BY - size of subimage to use for correlation
;       HX, HY - center of subimage
;       EDGE - If set, use sobel filter before fft
;	NO_LIM -- do not set upper limit of shifts applied.
;       SHREFIDX -- Reference image index into input geom array.
;	SHREFIMG -- Reference image.  Takes precedence over SHREFIDX.
; OUTPUTS:
; 	SHIFTS is the 2 x nimages array of shifts calculated.
; OPTIONAL OUTPUT PARAMETERS:
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
;	Calls xyoff2 to find the displacement of each image in the
;	input file from the reference image, taken to be the central
;	image in the file.
; MODIFICATION HISTORY:
;	GC   Jan 94
;	TRM 1994-10-05 Modified to use xyoff2 rather than xyoff and to use
;               a Hanning window.  Also use only central 128x128.
;	DLM 1995-04-03 Print shifts from here, not in xyoff2. (use ivm_xyoff2)
;       JPW 1996-05-14 Added keywords bx,by,hx,hy, and streamlined code
;       ELW 2009-11-20 Added 'geod' input array;
;			read from 'geof' if set AND filein==1
;       ELW 2010-02-17 Added 'vb' to choose verbosity
;       ELW 2010-05-26 filein parameter, better input logic
;       ELW 2011-04-05 Added/shuffled keyword to allow for an external image
;			reference, rather than just an index.  Useful primarily
;			for timeseries data. SHREFIDX is now what SHREFIMG was.
;-

PRO SHIFTER_2D, NN,nimages,SHIFTS, vb, filein, geod=geod,geof=geof, $
		bx=bx,by=by,hx=hx,hy=hy,edge=edge,detrend=detrend,  $
		no_lim=no_lim, shrefidx=shrefidx, shrefimg=shrefimg

;
; Declarations
;
shifts=fltarr(2,nimages)	;Shift values for all frames

;
; Open file; read reference image (the central one).
;
if (vb ge 1) then print,'Shifter_2d: ',systime()

if filein && keyword_set(geof) then begin
    geod=fltarr(nn,nn,nimages)
    openr,lun, geof, /get_lun
    readu,lun, geod
    close,lun & free_lun,lun
endif else if ~keyword_set(geod) then begin
    print, 'shifter_2d: Input array (geod) not set.  Stopping...'
    stop
endif

; Get reference image from 'shrefimg' or from geometry input
if (n_elements(shrefimg) eq long(nn)*long(nn)) then begin
    ref = shrefimg 		; External ref image passed in via 'shrefimg'
endif else begin
    ; Determine index into geom array, according to whether 'shrefidx' is set
    if (keyword_set(shrefidx)) then begin
	if n_elements(shrefidx) NE 1 then shrefidx = nimages/2L
    endif else begin
	shrefidx=nimages/2L
    endelse
    ref=geod[*,*,shrefidx]  	; Better than choosing the first image 
endelse

if n_elements(bx) ne 1 then bx = 128
if n_elements(by) ne 1 then by = 128
bx = bx < nn
by = by < nn
; hx,hy passed to ivm_xyoff2.  Default is nx/2, ny/2 unless keyword-set here 

for i=0,nimages-1 do begin
   shifts[*,i]=ivm_xyoff2(ref,geod[*,*,i],bx,by,hanning(bx,by),1.,hx,hy,fna,$
                         sobel = edge,detrend=detrend)
   if (vb ge 2) then print,i,':    dx = ',shifts[0,i],'    dy = ',shifts[1,i]
endfor

upshift = 5.  		              ;Maximum shift univerally accepted
if (keyword_set(no_lim)) then upshift = nn/2.
				      ; Start trapping for bad images:
wh_all=where(abs(shifts) ge upshift,count_all)
if (count_all ne 0) then begin
   print,'WARNING!  images ',wh_all / 2,' have large shifts!'
   wh=where(abs(shifts) ge 3*upshift,count)
   if (count ne 0) then begin
      print,'limiting images ',wh / 2,'to zero shift'
      shifts[wh]=0. 
   endif
endif				      
end