function remove_cr,im1,h1,crnt_im,h2, again, shrink, ALL = all, INIT=init, N_CR=n_cr, PF=pf
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