function remove_cr,im1,h1,crnt_im,h2, again, shrink, ALL = all, INIT=init, N_CR=n_cr, PF=pf
;+
; NAME:
;       REMOVE_CR 
;
; PURPOSE:
;       Removes cosmic rays (and other bright pixels which are not stars). 
;
; PROCEDURE:
;	Compares each pixel of crnt_im with the pixel in im1 where 
;	a star would be if there was a star in that pixel. If the ratio
;	exceeds 'threshold', then that pixel = median value of superpixel
;	in the previous image.
;
; CALLING SEQUENCE:
;       Result = REMOVE_CR(im1,h1,crnt_im,h2,again) 
;
; INPUTS:
;       im1	(num)ARR	previous image
;       h1	STRUCT		header from previous image
;	crnt_im (num)ARR	current image
;	h2	STRUCT		header from current image
;	AGAIN	STRING 		'y' to prompt for ROI
;	shrink	FLOAT		Factor by which the image is shrunk from 1024x1024
;
; KEYWORDS:
;	ALL	remove cosmic rays AND stars
;	INIT	define ROI on first time through
;	SUMSTARS	Sum star images
;	PF	Return point-filtered image
;
; Returned in Common block:
;	N_CR	Returns number of pixels removed
; 	COORDS	Returns coordinates of removed pixels
;
;
; OUTPUTS:
;       crnt_im with cosmic rays removed is returned
;
; PROCEDURE:
; Compute point filtered image. Subtract point filtered image from current image (imdiff). 
; Subtract background image from point filtered image (bkgdiff).
; Compute std dev of each superpixel of bkgdiff. Candidate crays 
; and stars are where imdiff is greater than TWICE the interpolated standard deviation for each 
; superpixel of bkgdiff. For each candidate pixel, check in the previous image in a 5x5 box 
; where that pixel should have been if it was a star, based on a given star motion across 
; the FOV. If the maximum value from this box is less than the value of that pixel in the 
; current image minus the standard deviation of the superpixel from diff, than replace with 
; the median value of the superpixel of the previous image.

;

; MODIFICATION HISTORY:
;               Written 12/23/98, N Rich - NRL/Interferometrics
;	4/12/99 NBR	Use median-filtered image to check for CRs
;	7/14/99 NBR	Fix problem when gaps in im1
;	11/14/00 NBR	Add N_CR keyword, reduce_history common block; use point_filter
;	11/21/00 NBR - Optimized for raw images
;	12/28/00 NBR - Use C3clearmask2.fts for defining no-check area
;	 7/19/01 NBR - Tinkering; using dn/sec images
;
; LIMITATIONS:
;	- may remove comets, planets, etc. if not exempted by setting 
;	  again to 'y' which calls the DEFROI routine
;
version = '@(#)remove_cr.pro	1.3 09/13/10'  ; NRL LASCO IDL LIBRARY
;
;-

COMMON STAR_SUMS, allstarims,numims, bkg, count, coords
COMMON reduce_history, cmnver, prev_a, prev_hdr, zblocks0

cmnver = STRMID(version,4,strlen(version))

mark = [1025,1025]
coords=0	
;x3_pixpersec = 0.00071753D	; C3 pix/sec, 9/10/99
;y3_pixpersec = 0.000009587D	; C3 pix/sec, 9/10/99
;x3_pixpersec = 0.00073841084D	; value for 1999/03/24
;y3_pixpersec = 0.000033187545D	; value for 1999/03/24
;

IF datatype(h1) NE 'STC' THEN h1 = lasco_fitshdr2struct(h1)
IF datatype(h2) NE 'STC' THEN h2 = lasco_fitshdr2struct(h2)
IF datatype(again) EQ 'UND' THEN again = 'n'
IF datatype(shrink) EQ 'UND' THEN shrink = 1.
maxmin,im1

IF keyword_set(ALL) THEN print,'Removing cosmic rays AND stars...' $
	ELSE print,'Removing cosmic rays...'

nc1 = !d.table_size-1	;# of colors available
sz = size(im1)
nxs = sz(1)
nys = sz(2)
sumfac = (h1.lebxsum > h1.sumcol)
mask=replicate(1,nxs,nys)
zeros=where(im1 EQ 0)
IF zeros(0) NE -1 THEN mask(zeros)=1
dmax = 5000
dmin= 0

IF datatype(im2) NE 'BYT' THEN BEGIN
	if h2.detector EQ 'C3' THEN threshold = 0.02 ELSE threshold = 0.03
ENDIF ELSE threshold = 0.3

count=0L

IF nys EQ 1024 and NOT(keyword_set(INIT)) THEN BEGIN
   IF h2.detector EQ 'C3' THEN BEGIN
	;restore,'$imax/c3ctrarea.sav'
	mask = READFITS('$ANCIL_DATA/calib/c3clearmask2.fts')
   ENDIF ELSE BEGIN
	cx1=382
	cx2=641
	cy1=365
	cy2=645
	xvals=[cx1,cx1,cx2,cx2]
	yvals=[cy1,cy2,cy2,cy1]
	mask(cx1:cx2,cy1:cy2)=2		; do not filter image center
   ENDELSE
   ctr_area = where(mask EQ 2 or mask EQ 0)   ; do not filter occulter mask or outer mask area
   boxsize = 15
ENDIF ELSE BEGIN
   IF init GT 1 THEN restore , '/tmp/roi.sav'
   boxsize = 7
ENDELSE

IF keyword_set(INIT) THEN BEGIN
	IF init GT 1 THEN goto, done
	window,0,xsize=nxs,ysize=nys
	tv,bytscl(crnt_im,max=dmax,min=dmin)
	print,string(7B)	; Bell
	print,'Define a region where not to do cosmic ray filtering:'
	print
	ctr_area=defroi(nxs,nys,xvals,yvals)
	init=2		; after first time, notify that INIT was activated
	save,ctr_area,xvals,yvals,filename='/tmp/roi.sav'
	done:
ENDIF
im2=crnt_im>0
im1=im1>0
mask=intarr(nxs,nys)
mask(ctr_area)=1		; do not filter image center

IF again EQ 'n' THEN BEGIN
	wset,0
	tvscl,crnt_im<dmax>dmin
	device, set_graphics=6             ;Set XOR mode
;	polyfill,xvals,yvals ,$
;		/dev,color = nc1,/noclip 	;Complete polygon
	tvscl,mask
	device, set_graphics=3             ;
	print,'Waiting before processing...'
	wait,2
ENDIF ELSE BEGIN

WHILE again NE '' DO BEGIN
	window,0,xsize=nxs,ysize=nys
	tvscl,crnt_im<dmax>dmin
	print,string(7B)	; Bell
	device, set_graphics=6             ;Set XOR mode
	polyfill,xvals,yvals,$
		/dev,color = nc1,/noclip 	;Complete polygon
	device, set_graphics=3             ;
	print,'Define an additional region where not to do cosmic ray filtering:'
	print
	roi=defroi(nxs,nys)
	IF roi(0) GE 0 THEN BEGIN
		mask(roi)=1
		read,'Enter any char and <RETURN> to define another ROI, else just <RETURN>: ',again
	ENDIF ELSE again='n'
ENDWHILE
ENDELSE
IF again EQ '' THEN again='y'

time2 = h2.date_obs+' '+h2.time_obs
t1=utc2tai(str2utc(h1.date_obs+' '+h1.time_obs))
t2=utc2tai(str2utc(time2))
tdiff=t2-t1
;
; ** Begin compute star motion

get_star_motion,h2,x_pixpersec,y_pixpersec
print,x_pixpersec,y_pixpersec
;
; ** End compute star motion
;
xdiff = x_pixpersec*tdiff*shrink
ydiff = y_pixpersec*tdiff*shrink

IF h2.detector EQ 'C3' THEN BEGIN
;	xdiff=x3_pixpersec*tdiff*shrink
;	ydiff=y3_pixpersec*tdiff*shrink
	bacfac = 1.06
	boxsize=15
ENDIF ELSE BEGIN
;   IF nxs EQ 1024 THEN BEGIN		; this is for if distortion correction is applied
;	xdiff=0.00347*tdiff*shrink	; nbr, 4/15/99
;	ydiff=-0.000441*tdiff*shrink
;   ENDIF ELSE BEGIN
;	xdiff=0.00341*tdiff*shrink	; nbr, 4/12/99, for video
;	ydiff=0.000344*tdiff*shrink
;   ENDELSE
   bacfac = 1.12
ENDELSE

IF keyword_set(ALL) or keyword_set(PF) THEN BEGIN
	xdiff=0.001*tdiff
	ydiff=-0.0005*tdiff
	n_iter=2
ENDIF ELSE n_iter=1

xdiff = xdiff/sumfac
ydiff = ydiff/(h1.lebysum > h1.sumrow)
xsp = nxs/16
ysp = nys/16
m1=fltarr(xsp,ysp)
m2=fltarr(xsp,ysp)
stdv = m1
stdv_d=m1
stdv_med2=m1
FOR i=0,xsp-1 DO BEGIN		; ** Compute medians of quarter superpixels, std dev of im1
   FOR j=0,ysp-1 DO BEGIN
	spp = im2(i*16:(i+1)*16-1,j*16:(j+1)*16-1)
	m2(i,j) = median(spp)
	m1(i,j) = median(im1(i*16:(i+1)*16-1,j*16:(j+1)*16-1))
   ENDFOR
ENDFOR
print,'OK, beginning...'
n1=systime(1)
median2t = median(im2)
t0=systime(1)
help,bkg
IF datatype(bkg) EQ 'UND' THEN bkg=getbkgimg(h2,bh)/bh.exptime
tvscl,(im2-bkg)<10>0
;stop

; Compute point filtered image. 
print,'Doing point_filter.....'
med2=filter_image(im2,MEDIAN=5,/ALL)
;point_filter,im2,boxsize,3,n_iter,med2,outpts

help,boxsize
cmnver = cmnver+' boxsize='+TRIM(STRING(boxsize))

print,'med2:',fix(systime(1)-t0),' seconds'

IF keyword_set(PF) THEN BEGIN
   sz=size(outpts)
   n_cr = sz[1]
   help,n_cr
   return,med2
ENDIF ELSE BEGIN
; Subtract point filtered image from current image.

diff2 = im2 - med2

t0=systime(1)
;fdiff2 = filter_image(diff2,median=5,/all)

; Subtract background image from median filtered image
diffbkg=med2-bkg
mediffbkg=median(diffbkg)
diffbkg(ctr_area)=mediffbkg
diffbkg=diffbkg>0
print,'diffbkg:',fix(systime(1)-t0),' seconds'
; Compute std dev of each superpixel of bkg difference image and of point filtered image

FOR i=0,xsp-1 DO BEGIN		
   FOR j=0,ysp-1 DO BEGIN
	spp = diffbkg(i*16:(i+1)*16-1,j*16:(j+1)*16-1)
	result = moment(spp,SDEV=stdv0)
	stdv_d(i,j)=stdv0
	spp = med2(i*16:(i+1)*16-1,j*16:(j+1)*16-1)
	result = moment(spp,sdev=stdv0)
	stdv_med2(i,j)=stdv0
   ENDFOR
ENDFOR
;stdv_d = stdv_d<100

t0=systime(1)
;; PREVIOUSLY: Candidate crays and stars are where (current image minus pointfiltered image) is greater than the standard deviation for each superpixel of bkg difference image
;
;nog = where(diff2 GT rebin(stdv_med2,nxs,nys),num)

zz = where(med2 LE 0)
IF zz[0] NE -1 THEN med2[zz]=1.
nog = where((diff2/med2) GT 0.05,num)

print,'nog:',fix(systime(1)-t0),' seconds'

; For each candidate pixel, check in the previous image in a 5x5 box where that pixel should have been if it was a star, based on a given star motion across the FOV. If the maximum value from this box is less than the value of that pixel in the current image minus the standard deviation of the superpixel from diff, than replace with the value from the median image.

;stop
ind = 0L
t0=systime(1)
starimg = fltarr(nxs,nys)
;tvscl,(im1-bkg)<1000>0

WHILE ind LE num-1 DO BEGIN
   IF mask(nog(ind)) EQ 1 THEN goto, skip
   one2two,nog(ind),[nxs,nys],x,y
   x1 = x-xdiff
   y1 = y-ydiff
   IF x1 GE 2 and x1 LE nxs-3 and y1 GE 2 and y1 LE nys-3 THEN BEGIN 
	sub1 = im1(x1-2.0:x1+2.0,y1-2.0:y1+2.0)
   	maxsub1=max(sub1)
	;IF diff2(nog(ind)) LE 0.04 THEN fac = 1.+diff2(nog(ind)) ELSE fac = bacfac
	boxcol=80
	flag='sta'
	;variat=INTERPOLATE(stdv_d,x/16.,y/16.)
	variat=3*sqrt(med2[x,y]*13.)/13.
	;newval = INTERPOLATE(m2,x/16.,y/16)
	newval= med2[x,y]
   	IF  maxsub1 LT (im2(x,y)-variat) THEN BEGIN	; OK, not a star
	   im2(x,y) =  newval
	   coords=[coords,nog[ind]]
	   boxcol=160
	   flag= 'nos'
   	   count=count+1
	ENDIF ELSE BEGIN
	   IF maxsub1 LT (im2[x,y]-variat/2.) THEN BEGIN
		im2[x,y] = maxsub1
		coords=[coords,nog[ind]]
		boxcol = 240
	 	flag = 'ns1'
   	   	count=count+1
	   ENDIF ELSE starimg(x,y)=crnt_im(x,y)
	ENDELSE
	print,flag,x,y,maxsub1,crnt_im[x,y],variat,newval
	wset,0
	tvbox,5,x,y,color=boxcol
	;wset,2
	;IF flag EQ 'nos' THEN surface,sub1,/lego
	;wait,1
   ENDIF
   skip:
   ind = ind+1
ENDWHILE
wset,0
print,'star check:',fix(systime(1)-t0),' seconds'
count = n_elements(coords)-1
;IF keyword_set(SUMSTARS) THEN BEGIN
   IF datatype(allstarims) EQ 'UND' THEN allstarims=fltarr(1024,1024)
   allstarims = allstarims+starimg
   IF datatype(numims) EQ 'UND' THEN numims=0
   numims = numims + 1
;ENDIF
n2=systime(1)
print,' Remove_cr took',n2-n1,' seconds and removed',count,' out of',num,' pixels tested.'
n_cr=count

array = bytarr(nxs,nys)
array[*]=1
array[coords]=200
;tvscl,array
;tvscl,im2<dmax

wait,4
;stop
return,im2

ENDELSE

end