;+
; NAME:
; imgclean
;
; PURPOSE:
; This procedure is designed to remove cosmic ray hits (CRs) from an image,
; especially WFPC images. The procedure's effectiveness varies from image
; to image. Often, it does an excellent job, but sometimes
; the operational parameters must be modified (via keywords) to tweak
; the performance of the software.
; It should not be assumed that this procedure will automatically do a good
; job! Keep a close eye on before and after images. Occasionally, a faint
; star gets obliterated, especially if there is a cosmic ray around 3 pixels
; from it. Naturally, this procedure is only recommended when a CR split is
; not possible, as that will often yield better results.
; Run IMGclean twice on the same image to achieve best cleaning results
; especially on images with very dense CR's.
; Still in a somewhat primitive state, this procedure still need some
; improvement to the algorithm. Suggestions welcome to deutsch@stsci.edu.
; Use the CLEANEXAMINE main program to examine the results of the IMGclean
; run.
; Currently, IMGclean is tweaked for PC images. I haven't worked with
; WF images in a while, so you may to diddle with the parameters for best
; effect.
;
; CALLING SEQUENCE:
; imgclean,IMAGE,HEADER,[CR_ARRAY],[optional keywords]'
; imgclean,/help
;
; INPUT/OUTPUT PARAMETERS:
; IMAGE This parameter supplies the image in which the CRs are to be
; removed. This image is MODIFIED! The output cleaned image
; is returned via this variable.
; HEADER This parameter is the FITS header. The header will be modified
; with added HISTORY notes and keywords detailing the settings
; used on the last IMGclean run.
;
; OPTIONAL OUTPUT:
; CR_ARRAY This array is a BYTE array of the same size as the input image
; containing some information of which pixels the software
; considered, and what decisions it made as far as CR or STAR.
; DN Meaning
; 70 Low scan pixel not removed
; 120 Hi scan pixel not removed (i.e. suspected star)
; 255 Hi scan pixel removed (i.e. cosmic ray)
; 200 Neighboring Low scan pixel removed
;
; OPTIONAL INPUT:
; HELP Displays quick keyword defaults values
; TESTRUN If set, dumps a file with stats on each object considered
; STAR_PSF_SENS Sensitivity of the STAR checking.
; The higher the value, the more CRs will be left in the
; image because they are suspected as being stars. The
; lower the value, the more faint stars will be removed.
; Default=.35 Typical=.35, .25, .30, .40.
; -Value: .35 is probably the optimum value for doing a
; good job with CR's but leaving virtually faint
; scrungy stars for PC frames.
; -Value: .25 is probably the optimum value for doing
; the best possible job of cleaning CR's at the
; expense of removing a few more faint scrungy stars
; than is proper. If faint, scrungy stars are not
; important, this may be the best setting (again, PC).
; SKY_VAL_SAMP_SZ Size of the boxes used when calculating the local
; background. The X,Y image size must be a multiple of the
; box size. Use smaller boxes for steeper images.
; Default values are 10 or 16.
; HI_SCAN_HEIGHT Pixel value above the sky value in sky RMS units that is
; flagged as a suspected CR. Default=6. Typical=5,6,7
; LO_SCAN_HEIGHT Pixel value above the sky value in sky RMS units that is
; removed if it is adjacent to a positively identified CR,
; provided FINE_CLEAN=1. Default=1.5. Typical=2,2.5,3
; FINE_CLEAN Set (=1) this keyword (Flag) if pixels above LO_SCAN_HEIGHT
; neighboring pixels identified as CRs are to be removed.
; PHASE1_ITER Number of iterations for Phase 1. Should always be 2
; because 1 iteration tends to leave more CRs behind.
; QPHOT If set, quick photometry is run on stars using qphot.pro
; NOT SUPPORTED!!!!
; NOCLEAN If set, all objects are assumed to be stars and nothing
; is checked or removed. Used in conjunction with QPHOT
; on an already cleaned image.
;
; HISTORY:
; 02-SEP-90 Version 1.0 written by Eric W. Deutsch (EWD)
; 27-JUN-91 Version 1.1 released with several minor modifications including
; changing the operational parameters to keywords. (EWD)
; 27-JAN-92 Version 1.2 released with the additions of changing the
; operational parameters to keywords and adding better on-line
; help. (EWD)
; 07-MAY-92 Version 1.21 released: proper header added and /HELP keyword
; added. EWD.
; 10-MAY-92 Version 1.22: Fixed "Donutting" problem and changed default
; STAR_PSF_SENS value to .50 instead of .66 since it was removing
; too many actual stars. EWD.
; 24-JUL-92 Version 1.3: Better Star/CR determination algorithm added.
; This algorithm is less Kludgy, but sometimes gives poorer
; results than the v1.22 release. Needs Work! EWD.
; 08-DEC-92 Version 1.31: Tidied things up a bit, but nothing major. EWD.
; 09-DEC-92 Version 1.32: Fixed a rather major bug in the focus pixel
; identifier and fixed the logic for negative ratios. This
; now shows a major improvement! The major task now is to add
; a special identifier for streak CR hits. EWD
; 04-JAN-93 Version 1.40: Changed the default sensitivity to 0.35 which
; make IMGclean more lenient as far as stars go, but added
; another simple-minded neighboring pixels check which catches
; A LOT of CR's. Basically, I say that if more than 2 of the
; 8 neighboring pixels aren't above the low scan hight, then
; it MUST be a cosmic ray. This now shows another big im-
; provement in number of CR's removed as well as star removed.
; This also does a better job at removing long streaks, but
; an algorithm specifically for longs streaks is needful and
; will be implemented soon. Header updated. Eric W. Deutsch
; 05-JAN-93 Version 1.41: Little Faster. Fixed FITS header updating. EWD
; 06-JAN-93 Added a few more comments and text to program header. EWD
; 28-FEB-93 Version 1.42: Added /QPHOT and /NOCLEAN options and code so that
; IMGclean can do some rudimentary quick aperture photometry. EWD
;-
pro skyline,origline,skyv,rms,interactive
;+
; NAME:
; skyline
;
; PURPOSE:
; Determine sky value and RMS for an image.
; Called by IMGCLEAN.
;
; HISTORY:
; Written by Eric Deutsch as part of the IMGCLEAN suite.
; Header added by PP/ACC April 16, 1997
;
;-
;
;;;if (n_params(0) lt 4) then interactive=0
tmp=where(origline ne 0)
if (tmp(0) eq -1) then begin
; print,'[SKYLINE]: Warning: piece of sky is all 0.00's!'
skyv=0.0 & rms=1.0
return
endif
line=origline(where(origline ne 0))
els=n_elements(line)
av1=avg(line)
rms1=sigma(line)
;;;if (interactive eq 1) then begin
if keyword_set(interactive) then begin
print,'1st Iteration: AVG=',strn(av1),',RMS=',strn(rms1)
plot,indgen(els),line
oplot,[0,els],[av1,av1]
oplot,[0,els],[av1+rms1,av1+rms1]
ch=get_kbrd(1)
endif
atmp=where(line lt av1+3*rms1)
chk=size(atmp)
if (chk(0) eq 0) then begin
a=line
goto,SKIP1
endif
a=line(atmp)
av1=avg(a)
rms1=sigma(a)
;;;if (interactive eq 1) then begin
if keyword_set(interactive) then begin
print,'2nd Iteration: AVG=',strn(av1),',RMS=',strn(rms1)
plot,indgen(els),a
oplot,[0,els],[av1,av1]
oplot,[0,els],[av1+rms1,av1+rms1]
ch=get_kbrd(1)
endif
SKIP1:
atmp=where((a lt av1+2*rms1) and (a gt av1-2*rms1))
chk=size(atmp)
if (chk(0) eq 0) then goto,DONE
a=a(atmp)
av1=avg(a)
rms1=sigma(a)
;;;if (interactive eq 1) then begin
if keyword_set(interactive) then begin
print,'3nd Iteration: AVG=',strn(av1),',RMS=',strn(rms1)
plot,indgen(els),a
oplot,[0,els],[av1,av1]
oplot,[0,els],[av1+rms1,av1+rms1]
ch=get_kbrd(1)
endif
DONE:
skyv=1.0
; print,'Skyvalue=',strn(av1),', RMS=',strn(rms1)
skyv=av1 & rms=rms1
return
end
FUNCTION SIGMA,ARRAY,N_PAR,DIMENSION
;+
; NAME:
; SIGMA
; PURPOSE:
; Calculate the standard deviation value of an array, or calculate the
; standard deviation over one dimension of an array as a function of all
; the other dimensions.
; CALLING SEQUENCE:
; RESULT = SIGMA(ARRAY)
; RESULT = SIGMA(ARRAY,N_PAR)
; RESULT = SIGMA(ARRAY,N_PAR,DIMENSION)
; INPUTS:
; ARRAY = Input array. May be any type except string.
; OPTIONAL INPUT PARAMETERS:
; N_PAR = Number of parameters. Default value is zero. The number of
; degrees of freedom is N_ELEMENTS(ARRAY) - N_PAR. The value
; of sigma varies as one over the square root of the number of
; degrees of freedom.
; DIMENSION = Optional dimension to do standard deviation over.
; OUTPUTS:
; The standard deviation value of the array when called with one
; parameter.
;
; If DIMENSION is passed, then the result is an array with all the
; dimensions of the input array except for the dimension specified,
; each element of which is the standard deviation of the corresponding
; vector in the input array.
;
; For example, if A is an array with dimensions of (3,4,5), then the
; command B = SIGMA(A,N,1) is equivalent to
;
; B = FLTARR(3,5)
; FOR J = 0,4 DO BEGIN
; FOR I = 0,2 DO BEGIN
; B(I,J) = SIGMA( A(I,*,J) , N )
; ENDFOR
; ENDFOR
; RESTRICTIONS:
; Dimension specified must be valid for the array passed; otherwise the
; input array is returned as the output array.
; PROCEDURE:
; Uses the function AVG.
; When DIMENSION is passed, then the function SUM is used.
; MODIFICATION HISTORY:
; William Thompson Applied Research Corporation
; July, 1986 8201 Corporate Drive
; Landover, MD 20785
; Converted to Version 2 July, 1990
;-
ON_ERROR,2
S = SIZE(ARRAY)
IF S(0) EQ 0 THEN MESSAGE,'*** Variable must be an array, name= ARRAY'
;
IF N_PARAMS(0) EQ 3 THEN BEGIN
IF ((DIMENSION GE 0) AND (DIMENSION LT S(0))) THEN BEGIN
SIG = SQRT( AVG(ARRAY^2,DIMENSION) - AVG(ARRAY,DIMENSION)^2 )
N = S(DIMENSION+1)
END ELSE $
MESSAGE,'*** Dimension out of range, name= ARRAY'
END ELSE BEGIN
DIFF = ARRAY - AVG(ARRAY)
SIG = SQRT( AVG( DIFF*DIFF ))
N = N_ELEMENTS(ARRAY)
ENDELSE
;
IF N_PARAMS(0) GE 2 THEN SIG = SIG * SQRT( N / ((N - N_PAR) > 1.) )
;
RETURN,SIG
END
pro imgclean,img,h,cr,$
SKY_VAL_SAMP_SZ = SKY_VAL_SAMP_SZ, FINE_CLEAN = FINE_CLEAN, $
HI_SCAN_HEIGHT = HI_SCAN_HEIGHT, LO_SCAN_HEIGHT = LO_SCAN_HEIGHT, $
STAR_PSF_SENS = STAR_PSF_SENS, PHASE1_ITER = PHASE1_ITER, $
HELP = HELP, TESTRUN = TESTRUN, NOCLEAN = NOCLEAN $
, interactive=interactive
;QPHOT = QPHOT qphot not supported
if (n_elements(help) eq 1) then goto,PRNT_HELP
if (n_params(0) lt 1) then begin
print,'Call> IMGCLEAN,image_array,header,[returned_CR_array],[several keywords]'
print,'e.g.> IMGCLEAN,img1,h1,cr1'
print,"Type> IMGCLEAN,/HELP for parameter listing/default values"
return
endif
s=size(img)
if (s(0) ne 2) then goto,PRNT_HELP
if (n_elements(SKY_VAL_SAMP_SZ) eq 0) then SKY_VAL_SAMP_SZ=10.
if (n_elements(HI_SCAN_HEIGHT) eq 0) then HI_SCAN_HEIGHT=6.
if (n_elements(LO_SCAN_HEIGHT) eq 0) then LO_SCAN_HEIGHT=1.5
if (n_elements(STAR_PSF_SENS) eq 0) then STAR_PSF_SENS=.35
if (n_elements(FINE_CLEAN) eq 0) then FINE_CLEAN=1
if (n_elements(PHASE1_ITER) eq 0) then PHASE1_ITER=2
if (n_elements(testrun) eq 0) then testrun=0
if (n_elements(qphot) eq 0) then qphot=0
if (n_elements(noclean) eq 0) then noclean=0
if (testrun eq 1) then openw,5,'imgclean.dmp'
if (qphot ge 1) then begin
openw,8,'qphot.dmp'
qphot,img,-999,-999
endif
if keyword_set(interactive) then begin
ICversion='1.41'
print,'[IMGclean '+ICversion+'] MESSAGE: Released 05-JAN-1993. '
print,' Please report problems or bugs to deutsch@stsci.edu.'
print,' Suggestions for improvement are welcome.'
print,' Include example images (or portions thereof) if possible.'
endif
s=size(img)
NAXIS1=s(1) & NAXIS2=s(2)
cr=bytarr(NAXIS1,NAXIS2)
flag=2
; yboxes=100 & xboxes=100 & flag=0 ; manual override
while (flag ge 2) do begin
xboxes=NAXIS1/SKY_VAL_SAMP_SZ & yboxes=NAXIS1/SKY_VAL_SAMP_SZ
if (xboxes eq fix(xboxes)) and (yboxes eq fix(yboxes)) then flag=0
if (flag ge 14) then flag=1
if (flag ge 2) then begin
if (SKY_VAL_SAMP_SZ eq 16) and ((flag and 4) eq 0) then begin
SKY_VAL_SAMP_SZ=10. & flag=flag+4 & endif
if (SKY_VAL_SAMP_SZ eq 10) and ((flag and 8) eq 0) then begin $
SKY_VAL_SAMP_SZ=16. & flag=flag+8 & endif
if (SKY_VAL_SAMP_SZ ne 10) and (SKY_VAL_SAMP_SZ ne 16) then $
SKY_VAL_SAMP_SZ=10.
endif
endwhile
if (flag eq 1) then begin
print,'[IMGClean] WARNING: Axes are not multiples of either 10 or 16.'
print,' Using boxsize of ',strn(fix(SKY_VAL_SAMP_SZ)), $
'. Edges of image may not be cleaned.'
xboxes=NAXIS1/SKY_VAL_SAMP_SZ & yboxes=NAXIS1/SKY_VAL_SAMP_SZ
endif
xboxes=fix(xboxes) & yboxes=fix(yboxes)
xboxsz=NAXIS1/xboxes & yboxsz=NAXIS2/yboxes
hipix=0L & lopix=0L & phase=0
PH1ITER:
if keyword_set(interactive) then begin
print,'PHASE 1(Iteration ',strn(phase+1),'): Scanning for High Pixels (Stars and CRs)....'
print,' High Scan height: ',strn(HI_SCAN_HEIGHT),'*RMS'
print,' Low Scan height: ',strn(LO_SCAN_HEIGHT),'*RMS'
endif
for y=0,yboxes-1-phase do begin
for x=0,xboxes-1-phase do begin
tmp=1 & crtmp=1
tmp=extrac(img,(x+phase/2.)*xboxsz,(y+phase/2.)*yboxsz,xboxsz,yboxsz)
skyline,tmp,skyv,rms,interactive
lo=where(tmp gt skyv+rms*LO_SCAN_HEIGHT)
chk=size(lo)
if (chk(0) ne 0) then begin
crtmp=cr((x+phase/2.)*xboxsz:(x+phase/2.+1)*xboxsz-1, $
(y+phase/2.)*yboxsz:(y+phase/2.+1)*yboxsz-1)
crtmp(lo)=70
cr((x+phase/2.)*xboxsz:(x+phase/2.+1)*xboxsz-1, $
(y+phase/2.)*yboxsz:(y+phase/2.+1)*yboxsz-1)=crtmp
lopix=lopix+n_elements(lo)
endif
high=where(tmp gt skyv+rms*HI_SCAN_HEIGHT)
chk=size(high)
if (chk(0) ne 0) then begin
crtmp=cr((x+phase/2.)*xboxsz:(x+phase/2.+1)*xboxsz-1, $
(y+phase/2.)*yboxsz:(y+phase/2.+1)*yboxsz-1)
crtmp(high)=120
cr((x+phase/2.)*xboxsz:(x+phase/2.+1)*xboxsz-1, $
(y+phase/2.)*yboxsz:(y+phase/2.+1)*yboxsz-1)=crtmp
hipix=hipix+n_elements(high)
endif
endfor
if keyword_set(interactive) then begin
if (y/5 eq y/5.) then $
print,strn(fix((y+1)*yboxsz*100./NAXIS2),length=3),'% complete:', $
strn(hipix,length=5),' high pixels found',lopix,' low-level pixels found'
endif
endfor
lopix=0L & hipix=0L
if (PHASE1_ITER eq 2) and (phase eq 0) then begin
phase=1 & goto,PH1ITER & endif
no_crpix=0 & no_cr=0 & no_starpix=0 & no_star=0 & t2=0
no_adj_CR=0 & no_adj_CR_pix=0
if keyword_set(interactive) then begin
print,'PHASE 2: Identifying High Pixels (Stars or CRs)....'
endif
for y=0,NAXIS2-1 do begin
band=where(cr(*,y) eq 120)
chk=size(band)
if (chk(0) ne 0) then begin
for i=0,n_elements(band)-1 do begin
current=lonarr(20) & t1=0 & flag=0
current(t1)=band(i) & t1=t1+1 & starti=i
while (flag eq 0) do begin
if (i eq n_elements(band)-1) or (t1 eq 19) then flag=2
if (flag eq 0) then begin
if (band(i) eq band(i+1)-1) then begin
i=i+1 & current(t1)=band(i) & t1=t1+1
endif else begin
flag=1
endelse
endif
endwhile
if (flag ne 0) then begin
cent=starti+t1/2
starchck,img,band(cent),y,ratio,radius,skyv,rms,test=testrun, $
crim=cr,sens=STAR_PSF_SENS,qphot=qphot,noclean=noclean $
, interactive=interactive
if (ratio lt STAR_PSF_SENS) then begin
is_a_CR=0
no_starpix=no_starpix+t1 & no_star=no_star+1
endif else begin
is_a_CR=1
no_crpix=no_crpix+t1 & no_cr=no_cr+1
tmp=extrac(img,band(cent)-7,y-7,15,15)
skyline,tmp,skyv,rms,interactive
endelse
seed=skyv
for t2=starti,i do begin
if (is_a_CR eq 1) then begin
cr(band(t2),y)=255
n_rnd=n_elements(band(t2)) & rnd=fltarr(n_rnd)
for rndi=0,n_rnd-1 do rnd(rndi)=randomu(seed)*rms*1.5-rms*.5
img(band(t2),y)=skyv+rnd
endif
endfor
if (is_a_CR eq 1) and (band(starti) gt 0) and (band(i) lt NAXIS1-1) $
and (y gt 0) and (y lt NAXIS2-1) then begin
subcr=extrac(cr,band(starti)-1,y-1,band(i)-band(starti)+3,3)
lowcr=where(subcr eq 70)
chk=size(lowcr)
if (chk(0) ne 0) then begin
subimg=extrac(img,band(starti)-1,y-1,band(i)-band(starti)+3,3)
n_rnd=n_elements(lowcr) & rnd=fltarr(n_rnd)
for rndi=0,n_rnd-1 do rnd(rndi)=randomu(seed)*rms*1.5-rms*.5
subimg(lowcr)=skyv+rnd
subcr(lowcr)=200
img(band(starti)-1:band(i)+1,y-1:y+1)=subimg
cr(band(starti)-1:band(i)+1,y-1:y+1)=subcr
no_adj_CR=no_adj_CR+1
no_adj_CR_pix=no_adj_CR_pix+n_elements(lowcr)
endif
endif
endif
;;;if (skyv gt 200) then print,' Large Sky!: ',skyv,rms
if (skyv gt 200 and keyword_set(interactive)) then begin
print,' Large Sky!: ',skyv,rms
endif
endfor
endif
if keyword_set(interactive) then begin
if (y/25. eq y/25) then $
print,strn(fix((y+1)*100./NAXIS2),length=2),'% complete', $
': (Ln-Obj,Pix): CRs=',vect([no_cr,no_crpix]), $
' AdjCRs=',vect([no_adj_CR,no_adj_CR_pix]), $
' Stars=',vect([no_star,no_starpix])
endif
endfor
if keyword_set(interactive) then begin
print,'*** Final: ', $
': (Ln-Obj,Pix): CRs=',vect([no_cr,no_crpix]), $
' AdjCRs=',vect([no_adj_CR,no_adj_CR_pix]), $
' Stars=',vect([no_star,no_starpix])
endif
if (testrun eq 1) then close,5
if (qphot eq 1) then close,8
s=size(h)
if (n_elements(s) lt 3) then return
if (s(2) ne 7) then return
on_error,2 ; return if no SX routines available
tmp=sxpar(h,'IC_CRSRM') & iter='_'
if (!ERR ge 0) then begin
i=2 & next=9
while (i lt 8) do begin
tmp=sxpar(h,'IC'+strn(i)+'CRSRM')
if (!ERR lt 0) then begin & next=i & i=9 & endif
i=i+1
endwhile
iter=strn(next)
endif
sxaddhist,'[IMGclean '+ICversion+'] '+!stime+': IC'+iter+' keywords added',h
sxaddpar,h,'IC'+iter+'CRSRM',no_crpix,' Number of CR pixels removed'
sxaddpar,h,'IC'+iter+'ADJRM',no_adj_CR_pix,' Number of pixels adjacent to CRs removed'
sxaddpar,h,'IC'+iter+'STRPX',no_starpix,' Pixels part of a star (not removed)'
sxaddpar,h,'IC'+iter+'SMPSZ',SKY_VAL_SAMP_SZ,' Size of Box for local Sky Value'
sxaddpar,h,'IC'+iter+'HISCN',HI_SCAN_HEIGHT,' RMSs above sky for HI classif.'
sxaddpar,h,'IC'+iter+'LOSCN',LO_SCAN_HEIGHT,' RMSs above sky for LO classif.'
sxaddpar,h,'IC'+iter+'PSFSN',STAR_PSF_SENS,' Sensitivity of Star/CR Discrimination'
sxaddpar,h,'IC'+iter+'FINCL',FINE_CLEAN,' Adjacent Pixels removed? (1=Yes)'
sxaddpar,h,'IC'+iter+'PH1IT',PHASE1_ITER,' Iterations of Phase 1 (Find)'
return
PRNT_HELP:
print,'Call> IMGCLEAN,image_array,header,[returned_CR_array],[optional keywords]
print,'Keywords:
print,' SKY_VAL_SAMP_SZ Default: 10.0 pixels squared'
print,' HI_SCAN_HEIGHT Default: 6.0 RMS units above local sky'
print,' LO_SCAN_HEIGHT Default: 1.5 RMS units above local sky'
print,' STAR_PSF_SENS Default: 0.35 height to volume ratio'
print,' FINE_CLEAN Default: 1 (TRUE) (Remove highish adjacent pixels'
print,' PHASE1_ITER Default: 2 iterations (each iteration slightly different).'
return
end