function rdif_hi, A, B, rawdif = rawdif, hdr = hdr, noshift = noshift, med = med, min = min, $
	side=side, tel=tel, nreg=nreg, local=local, scale=scale, scl=scl, pos=pos, split=split, $
	xout = xout, yout = yout, model=model, mfilt=mfilt,box=box, debug=debug, DEPROJECT=deproject
;
; $Id: rdif_hi.pro,v 1.11 2010/05/11 16:08:03 nathan Exp $
;
;   Subtracts A from B after appropriate shift
;
; Example Algorithm to weaken star field in running differences
; A --> Filename (with Full Path) Raw image at t = 0
; B --> Filename (with Full Path) Raw image at t = t + Delta t
;
; CALLING SEQUENCE:
; files = findfile('/net/earth/earth2/secchi/lz/hi_1/L0/a/img/20070204/????????_??????_s4*.fts')
; img = rdif_hi(files(0), files(1), xout=xout, yout=yout, rawdif=rawdif)
; tvscl, img>(-100)<100  (for images taken 2 hours apart)
; tvscl, img>(-200)<200  (for images taken 4 hours apart)
;
; INPUT KEYWORDS:
;   BOX = [x1,x2,y1,y2]     Box to use for cross_correlation (defaults to area
;   	    near center of fov offset toward sun) based on 512x512 FFV
;   HDR = Header to use to getbkgimg, if input MODEL= is undefined
;   /NOSHIFT	Returns zero image
;   MED=    	Value to which to normalize result instead of median of model, only used
;   	    	if ratio of medians of 1st and 2nd image is LT 0.97.
;   MIN=    	Not used
;   SIDE=   	'a' or 'b'
;   TEL=    	'hi_1' or 'hi_2'
;   NREG=   	purpose unclear
;   /LOCAL  	Use sccreadfits instead of secchi_prep if input are filenames
;   SCALE=  	Factor by which to scale input before test_crosscorrelation
;   SCL=    	Factor by which to scale output of test_crosscorr before shift and diff
;   POS=    	Arraw containing x,y of start corner of x1,y1,x2,y2 of subfield to use
;   /SPLIT  	Not used
;   MODEL=  	Background model image to use (optional with HDR defined or filename input)
;   MFILT=  	Apply median filter to result
;   /DEBUG  	Print info messages and display images in process
;
; OUTPUT (Optional)
;   HDR =   header from input filename1 (if any)
;   xout=, yout=    Amount of shift 
;   rawdif= Standard running difference image
;
; $Log: rdif_hi.pro,v $
; Revision 1.11  2010/05/11 16:08:03  nathan
; remove non-finite values from test region before crosscorr
;
; Revision 1.10  2009/07/31 13:37:42  cooper
; changed format of info printed
;
; Revision 1.9  2009/04/23 19:41:21  nathan
; renamed fshift -> fshift2 to not conflict with yohkhoh/ucon/idl/sakao/fshift.pro
;
; Revision 1.8  2009/04/23 19:38:27  nathan
; add /DEPROJECT keyword
;
; Revision 1.7  2009/01/15 22:27:04  nathan
; added keyword explanations
;
; Revision 1.6  2008/12/10 16:16:21  nathan
; added /debug; use secchi_prep IF input are filenames
;
; Revision 1.5  2007/12/04 16:25:08  nathan
; added BOX=[x1,x2,y1,y2] functionality for rdif_hi.pro
;
; Revision 1.4  2007/11/30 16:08:01  sheeley
; added /match to getbkgimg call
;
; Revision 1.3  2007/11/29 18:33:39  nathan
; version from adam/star_rem/ on Nov 7, 2007
;

; Extract the day of the image to look for the daily model

	COMMON DIR_COMMON, mvi_dir, jmap_dir, median_dir, diff_dir, image_dir

	IF(NOT(keyword_set(SIDE))) THEN side = 'a'
	side = STRLOWCASE(side)
	IF(NOT(keyword_set(TEL))) THEN tel = 'hi_1'
	tel = STRLOWCASE(tel)
	IF(tel EQ 'hi1' OR tel EQ 'hi2') THEN tel = 'hi_' + strmid(tel,2)
	IF(~keyword_set(NREG)) THEN nreg = 1
	local = keyword_set(LOCAL)
	min = keyword_set(MIN)
	shift = ~keyword_set(NOSHIFT)

	filename = DATATYPE(A) EQ 'STR'
if keyword_set(DEBUG) THEN help,local

	IF(DATATYPE(A) EQ 'STR') THEN BEGIN
		IF(~local) THEN secchi_prep, A, hdr, img1, /NOCALFAC, /CALIMG_OFF $
		ELSE img1 = sccreadfits(A, hdr)
;		img1 = sccreadfits(A, hdr)
	ENDIF ELSE img1 = A

	IF(DATATYPE(B) EQ 'STR') THEN BEGIN
		IF(~local) THEN secchi_prep, B, hdr, img2, /NOCALFAC, /CALIMG_OFF $
		ELSE img2 = sccreadfits(B, hdr)
;		img2 = sccreadfits(B, hdr)
	ENDIF ELSE img2 = B

	IF(datatype(model) EQ 'UND') THEN model = scc_getbkgimg(hdr,outhdr=medhdr,/match)

    	s = size(img1)
    	nx1 = s[1]
    	nx2 = s[2]
	IF(~keyword_set(POS)) THEN pos = [0,0]
	IF(n_elements(POS) EQ 2) THEN BEGIN
		pos = [pos, nx1, nx2]
		model2 = model[pos[0]:pos[0]+nx1-1, pos[1]:pos[1]+nx2-1]
	ENDIF ELSE BEGIN
		img1 = img1[pos[0]:pos[0]+pos[2]-1, pos[1]:pos[1]+pos[3]-1]
		img2 = img2[pos[0]:pos[0]+pos[2]-1, pos[1]:pos[1]+pos[3]-1]
		model2 = model[pos[0]:pos[0]+pos[2]-1, pos[1]:pos[1]+pos[3]-1]
	ENDELSE
	IF(~keyword_set(MED)) THEN med = median(model2)

	med1 = median(img1)
	med2 = median(img2)
	med12 = med1 / med2
	IF(med12 GE 1) THEN med12 = 1./med12
	IF(med12 GE 0.97) THEN BEGIN
		med1 = 1.
		med2 = 1.
	ENDIF ELSE BEGIN
		print, 'median normalization'
		med1 = med / med1
		med2 = med / med2
	ENDELSE
if keyword_set(DEBUG) THEN help,med12,med

	IF(tel EQ 'hi_2') THEN BEGIN
		IF keyword_set(DEPROJECT) and keyword_set(HDR) THEN BEGIN
		    wcs = FITSHEAD2WCS(hdr)
		    dsun=hdr.dsun_obs
	    	    inp=fltarr(2,nx1,nx2)
		    inp[0,*,*]=float(img1*med1) - model2
		    inp[1,*,*]=float(img2*med2) - model2
		    coord = WCS_GET_COORD(wcs) ; array of coordinates for all pixels in original image
		    PLANE_PROJECT, inp, coord, outp, LIMIT=0, SPHER=1, SCEN=scen, RANGE=rad, $
			PIXEL_MAP=proj_map, /INTERP, NEW_COORD=proj_coord, DSUN_OBS=dsun, debug=debug
		    ;Make sure outp has even columns
		    sz2=size(outp)
		    IF (sz2[3] mod 2) THEN coln=sz2[3]-2 ELSE coln=sz2[3]-1
		    frame1=reform(outp[0,*,0:coln])
		    frame2=reform(outp[1,*,0:coln])
		ENDIF ELSE BEGIN
		    frame1 = float(img1*med1) - model2
		    frame2 = float(img2*med2) - model2
		ENDELSE
		szf=size(frame1)
		nx1=szf[1] & nx2=szf[2]
		dum1 = frame1 ;alog(frame1)		; to make background gradient smoother when taking 'derivative'
		dum2 = frame2 ;alog(frame2)		; to make background gradient smoother when taking 'derivative'
		dum_ant = frame1
		dum_act = frame2
	ENDIF ELSE BEGIN
		frame1 = float(img1*med1)
		szf=size(frame1)
		nx1=szf[1] & nx2=szf[2]
		dum1 = alog(frame1)
		frame2 = float(img2*med2)
		dum2 = alog(frame2)
		dum_ant = frame1 - model2
		dum_act = frame2 - model2
	ENDELSE
    	
	IF shift THEN BEGIN
		ant = dum1 - shift(dum1,0,3)	; test for best shifting
		act = dum2 -shift(dum2,0,3)	; test for best shifting
	ENDIF

        rawdif = (dum_act) - (dum_ant)                    ; to compare 
	IF keyword_set(DEBUG) THEN BEGIN
	    tvscl,rebin(rawdif,nx1/2,nx2/2)>(-1000)<1000
    	    print,'ant:'
	    maxmin,ant
	    print,'Act:'
	    maxmin,act
	    wait,3
	ENDIF
	IF(~keyword_set(SCALE)) THEN scale = 1
	IF(~keyword_set(SCL)) THEN scl = 1
	IF(~keyword_set(SPLIT)) THEN split = 1
	x = 0.0
	y = 0.0
	regx = 0
	regy = 0
	x1 = 0
	y1 = 0
	;s = size(img1)
	;nx1 = s[1]
	;nx2 = s[2]
	a2 = fltarr(nx1, nx2)
	shifts = fltarr(nreg, nreg, 2)
	IF shift THEN BEGIN
		FOR i=0, nreg-1 DO BEGIN
		FOR j=0, nreg-1 DO BEGIN
			x2 = floor((i+1) * float(nx1) / nreg - 1)
			y2 = floor((j+1) * float(nx2) / nreg - 1)
			
			reg1 = ant(x1:x2, y1:y2)
			reg2 = act(x1:x2, y1:y2)
			regx = x2 - x1 + 1
			regy = y2 - y1 + 1

			IF(scale NE 1) THEN BEGIN
				reg1 = rebin(temporary(reg1), scale*regx, scale*regy, /sample)
				reg2 = rebin(temporary(reg2), scale*regx, scale*regy, /sample)
			ENDIF
;print,'cross correllate'
; old #s a = (700:1000,300:700), b = (0:300,300:700) nrs 10/2/07
                        IF keyword_set(BOX) THEN BEGIN
			    box1=box*2
			    r1=box1[0]
			    r2=box1[1]
			    s1=box1[2]
			    s2=box1[3]
			ENDIF ELSE IF side EQ 'a' THEN BEGIN
			    IF tel EQ 'hi_2' THEN BEGIN
    	    	    	    ; h2 box for side A
				r1=230*2 ; 512->1024 means x2
                        	r2=450*2
                        	s1=130*2
                        	s2=400*2
			    ENDIF ELSE BEGIN
			    ; tel eq hi_1
			    	r1=0
				r2=500
				s1=700
				s2=1000
			    ENDELSE
			ENDIF ELSE BEGIN
			    IF tel EQ 'hi_2' THEN BEGIN
    	    	    	    ; h2 box for side B
                        	r1=139*2
                        	r2=311*2
                        	s1=044*2
                        	s2=473*2
			    ENDIF ELSE BEGIN
			    ; tel eq hi_1
			    	r1=500
				r2=1000
				s1=700
				s2=1000
			    ENDELSE
			ENDELSE
			   
			;h=test_crosscorr(reg1,reg2,x,y, /gauss)
			IF(nreg GT 1 OR nx1 LT 1024 OR nx2 LT 1024) THEN BEGIN
			    ccreg1=reg1
			    ccreg2=reg2
			ENDIF ELSE BEGIN
	    			print,'Using x=',STRTRIM(STRING(r1/2),2),':',STRTRIM(STRING(r2/2),2),$
				' y=',STRTRIM(STRING(s1/2),2),':',STRTRIM(STRING(s2/2),2),' for cross correlation.'
				ccreg1=reg1[r1:r2,s1:s2]
				ccreg2=reg2[r1:r2,s1:s2]
    			ENDELSE
			
			wf=where(finite(ccreg1),complement=wi1)
			wf=where(finite(ccreg2),complement=wi2)
			IF keyword_set(DEBUG) THEN help,wi1,wi2
			IF wi1[0] GT 0 THEN ccreg1[wi1]=0
			IF wi2[0] GT 0 THEN ccreg2[wi2]=0
			h=test_crosscorr(ccreg1,ccreg2,x,y, /gauss, verbose=debug)

;print, 'scaling by ', scl, ' for shift'
			x = x * float(scl)/scale
			y = y * float(scl)/scale

			reg1 = (dum_ant)(x1:x2, y1:y2)
			reg2 = (dum_act)(x1:x2, y1:y2)
			IF(scl NE 1) THEN BEGIN
				reg1 = rebin(temporary(reg1), scl*regx, scl*regy, /sample)
				reg2 = rebin(temporary(reg2), scl*regx, scl*regy, /sample)
			ENDIF
			d = fshift2(reg2,x,y) - reg1  ;fdiff(reg1, reg2, x,y) ;fshift better?
			IF(keyword_set(MFILT)) THEN BEGIN
				IF(mfilt NE scl) THEN d = congrid(temporary(d), mfilt*regx, mfilt*regy, /interp)
				d = median(temporary(d), 7*mfilt)
			ENDIF
			diff = rebin(d, regx, regy)

			a2(x1:x2, y1:y2) = diff

;			print, x,y
        		xout = x
		        yout = y
			shifts[i,j,0] = x
			shifts[i,j,1] = y
			y1 = y2+1
			IF(y1 GE nx2) THEN y1 = 0
;stop
			undefine, d  
			undefine, diff 
			undefine, reg1 
			undefine, reg2 
		ENDFOR
			x1 = x2+1
			IF(x1 GE nx1) THEN x1 = 0
		ENDFOR
	ENDIF

     ;; ant = act
     ;; dum_ant = dum_act
;stop
undefine, model2 
    
    return, a2
    
    ; a2 = d = fshift(frame2,x,y) - frame1
    
    ; x,y from  frame1, frame2 (shift frame2 to match frame1)
    
 end