Viewing contents of file '../idllib/contrib/buie/findsrc.pro'
;+
; NAME:
; findsrc
; PURPOSE:
; Automatic source detection and photometry from a digital image.
; DESCRIPTION:
;
; CATEGORY:
; CCD data processing
; CALLING SEQUENCE:
; findsrc,file
; INPUTS:
; file - Name of image file to search for sources.
; OPTIONAL INPUT PARAMETERS:
;
; KEYWORD INPUT PARAMETERS:
;
; EXTLIST - If image is a multi-extension FITS image, this list will
; force the reduction of only the extension numbers listed.
; The default is to do all the extensions, one at a time.
;
; GAIN - Gain of image, in photons/DN, default=1.0
;
; GAP - This is a number used to avoid looking at pixels near the
; object. It should be set to a value that is roughly equal
; to the FWHM of a typical stellar image. Default=2 pixels.
;
; KEYLIST - Name of a file containing a correspondence list. This list
; associates a set of standard names with the actual keyword
; names found in a FITS file header. If this keyword is
; omitted, a default list is used, as if a file with the
; following contents had been supplied:
; AIRMASS K AIRMASS
; DATE K DATE-OBS
; DATETMPL T DD-MM-YYYY
; EXPDELTA V 0.0
; EXPTIME K EXPTIME
; FILTER K FILTERS
; FILENAME K CCDFNAME
; OBJECT K OBJECT
; UT K UT
; The middle column is a flag. It may be K, for Keyword,
; T, for Template, or V, for Value. If it is V, the contents
; of the third field on that line should make sense for the
; name in the first field.
;
; MAXPHOTSIG- Maximum DN value for a useful signal. Any source with a peak
; above this level is passed over. Default=60000.0 DN
;
; NODISPLAY - Flag, when set will suppress all image display allowing program
; to be run in background or batch mode. This will be somewhat
; faster as well. The display steps take a small but non-trivial
; amount of time.
;
; OBJRAD - Radius of object aperture, in pixels, for photometry extraction.
; Default=GAP
;
; PATH - Optional path for original image directory.
; If not specified, the current directory is used.
;
; SIGTHRESH - Sigma threshold for source detection. Anything brighter
; than this many sigma above sky will be considered a source.
; Default = 2.5
;
; WINDOW - Size of region to average over in each direction, default=6
;
; OUTPUTS:
;
; KEYWORD OUTPUT PARAMETERS:
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; MODIFICATION HISTORY:
; 98/03/11, Written by Marc W. Buie, Lowell Observatory
; 98/03/22, MWB, added OBJRAD keyword
; 98/03/23, MWB, added EXTLIST keyword
;
;-
PRO findsrc_detect,a,sz,sigthresh,x0,dx,y0,dy,detect
fmt='($,a)'
for i=0,sz-1 do begin
if i eq 0 then begin
avg=shift(a,x0,y0)
endif else begin
avg=temporary(avg)+shift(a,x0+dx*i,y0+dy*i)
endelse
endfor
avg=temporary(avg)/sz
print,'.',format=fmt
for i=0,sz-1 do begin
if i eq 0 then begin
std=(shift(a,x0,y0)-avg)^2
endif else begin
std=temporary(std)+(shift(a,x0+dx*i,y0+dy*i)-avg)^2
endelse
endfor
std=sqrt(temporary(std)/(sz-1))
std = (a-avg)/temporary(std)
z=where(std gt sigthresh,count)
if count ne 0 then detect[z]=ishft(detect[z],1)
END
PRO findsrc_phot,fn,exttag,disp,detect,nx,ny,gain,a,binfac,exptime, $
sz,gap,objrad,maxphotsig,sigthresh,nobj
if disp then setwin,1
z=where(detect,count)
; z=where(detect eq 1,count)
if count ne 0 then begin
x = z mod nx
y = z / nx
zg=where(a[x,y] lt maxphotsig,countg)
if countg gt 0 then begin
if countg gt 0 and countg ne n_elements(xc) then begin
x = x[zg]
y = y[zg]
endif
basphote,gain,a,exptime,x,y,objrad,objrad+1.0,objrad+5.0,/nolog,/silent, $
xcen=xc,ycen=yc,mag=mag,err=err,fwhm=fwhm,max=max, $
skymean=sky,skyerr=skyerr,boxmrad=-1*objrad
zg = where(fwhm gt 1.0 and mag lt 30.0 and $
max-sky gt sigthresh*skyerr, countg)
if countg gt 0 and countg ne n_elements(xc) then begin
xc = xc[zg]
yc = yc[zg]
fwhm = fwhm[zg]
mag = mag[zg]
err = err[zg]
endif
endif
nobj=countg
if countg eq 0 then begin
print,'FINDSRC: Warning, no valid objects found!'
xc = 0.
yc = 0.
fwhm = 0.
mag = 99.999
err = 99.999
endif
if disp and nobj ne 0 then $
plots,xc,yc,psym=8,/device,symsize=2.0,color=120
;
; idx=0
; xc=fltarr(count)
; yc=fltarr(count)
; fwhm=fltarr(count)
; mag=fltarr(count)
; err=fltarr(count)
; for i=0,count-1 do begin
; x = z[i] mod nx
; y = z[i] / nx
; basphote,gain,a,exptime,x,y,objrad,objrad+1.0,objrad+5.0,/nolog,/silent, $
; xcen=xc0,ycen=yc0,mag=mag0,err=err0,fwhm=fwhm0,max=max, $
; skymean=sky0,skyerr=skyerr0
; if fwhm0 gt 1.0 and mag0 lt 30.0 and max lt maxphotsig and $
; max-sky0 gt sigthresh*skyerr0 then begin
; xc[idx]=xc0
; yc[idx]=yc0
; fwhm[idx]=fwhm0
; mag[idx]=mag0
; err[idx]=err0
; if disp then plots,xc0,yc0,psym=8,/device,symsize=2.0,color=120
; idx=idx+1
; endif
; endfor
;
; data=[[xc[0:idx-1]],[yc[0:idx-1]],[fwhm[0:idx-1]],[mag[0:idx-1]],[err[0:idx-1]]]
; nobj=idx
data=[[xc*binfac],[yc*binfac],[fwhm*binfac],[mag],[err]]
mkhdr,hdr,data
sxaddpar,hdr,'SIGTHRSH',sigthresh,' Sigma threshold for source detection'
sxaddpar,hdr,'GAP',gap*binfac,' Object gap size, is approximately the FWHM'
sxaddpar,hdr,'OBJRAD',objrad*binfac,' Object aperture radius for photometry'
sxaddpar,hdr,'SIGWSIZE',sz,' Sigma window size for source detection'
sxaddpar,hdr,'GAIN',gain/binfac,' Gain of CCD in e-/DN'
sxaddpar,hdr,'MAXSIG',maxphotsig,' Saturated above this DN level'
writefits,fn+'.src'+exttag,data,hdr
robomean,fwhm*binfac,3.0,0.5,avgfwhm
print,'('+strcompress(string(avgfwhm,format='(f10.1)')+')',/remove_all), $
format='($,a)'
endif
END
pro findsrc,fn,PATH=path,KEYLIST=keylist,GAP=gap,WINDOW=sz,OBJRAD=objrad, $
SIGTHRESH=sigthresh,GAIN=in_gain,MAXPHOTSIG=maxphotsig,NODISPLAY=nodisplay, $
EXTLIST=extlist,BINFAC=binfac
IF badpar(fn,7,0,CALLER='FINDSRC: (file) ') THEN return
IF badpar(keylist,[7,0],0,CALLER='FINDSRC: (KEYLIST) ',DEFAULT='[[DEFAULT]]') THEN return
IF badpar(path,[0,7],0, $
CALLER='FINDSRC: (PATH) ',DEFAULT='') THEN RETURN
if path ne '' then path=addslash(path)
IF badpar(gap,[0,2,3],0,caller='FINDSRC: (GAP) ',default=2) THEN return
IF badpar(sz,[0,2,3],0,caller='FINDSRC: (WINDOW) ',default=6) THEN return
IF badpar(sigthresh,[0,2,3,4,5],0,caller='FINDSRC: (SIGTHRESH) ',default=2.5) THEN return
IF badpar(in_gain,[0,2,3,4,5],0,caller='FINDSRC: (GAIN) ',default=1.0) THEN return
IF badpar(maxphotsig,[0,2,3,4,5],0,caller='FINDSRC: (MAXPHOTSIG) ',default=60000.0) THEN return
IF badpar(objrad,[0,2,3,4,5],0,caller='FINDSRC: (OBJRAD) ',default=gap) THEN return
if badpar(extlist,[0,1,2,3],[0,1],CALLER='FINDSRC: (EXTLIST) ', $
default=-1) then return
IF badpar(binfac,[0,2,3],0,caller='FINDSRC: (BINFAC) ',default=1) THEN return
gain=in_gain*(binfac^2)
loadkeys,keylist,hdrlist
disp = (!d.name eq 'X' or !d.name eq 'PS') and not keyword_set(nodisplay)
fmt='($,a)'
; Check header of image to see if it is a multi-extension image.
hdr=headfits(path+fn)
numext=sxpar(hdr,'NEXTEND')
IF numext eq 0 THEN BEGIN
extlist=0
ENDIF ELSE BEGIN
IF extlist[0] eq -1 THEN BEGIN
extlist=indgen(numext)+1
ENDIF ELSE BEGIN
IF max(extlist) gt numext THEN BEGIN
print,'FINDSRC: Input extension list is incompatible with the number of extensions'
print,'in the file. This file had ',numext,' extensions and the largest item in'
print,'your list is ',max(extlist)
print,'Aborting.'
return
ENDIF ELSE IF min(extlist) le 0 THEN BEGIN
print,'FINDSRC: Input extension list is invalid. You have one or more values less'
print,'than or equal to zero.'
print,'Aborting.'
return
ENDIF
ENDELSE
ENDELSE
numext=n_elements(extlist)
FOR ext=0,numext-1 DO BEGIN
IF extlist[ext] eq 0 THEN BEGIN
extstr = ''
exttag = ''
ENDIF ELSE BEGIN
extstr = strb36(extlist[ext])
exttag = 'x'+extstr
ENDELSE
print,fn+exttag+' '+systime(),format=fmt
a=0.
a=readfits(path+fn,hdr,exten_no=extlist[ext],/silent)
if binfac ne 1 then begin
arrsz=size(a)
nx=(arrsz[1]/binfac)*binfac
ny=(arrsz[2]/binfac)*binfac
print,' RB',format=fmt
a=rebin(a[0:nx-1,0:ny-1],nx/binfac,ny/binfac)
endif
parsekey,hdr,hdrlist,info
exptime=info.exptime
arrsz=size(a)
nx=arrsz[1]
ny=arrsz[2]
IF arrsz[3] ne 4 THEN a = float(a)
IF disp THEN BEGIN
setwin,0,xsize=nx,ysize=ny
tvscl,a
ENDIF
if ext eq 0 then detect=replicate(1B,nx,ny) else detect[*] = 1B
print,' UP',format=fmt
findsrc_detect,a,sz,sigthresh,0,0,-gap,-1,detect
print,'DOWN',format=fmt
findsrc_detect,a,sz,sigthresh,0,0,gap,1,detect
print,'RIGHT',format=fmt
findsrc_detect,a,sz,sigthresh,-gap,-1,0,0,detect
print,'LEFT',format=fmt
findsrc_detect,a,sz,sigthresh,gap,1,0,0,detect
print,' THRSH',format=fmt
; requires detection in 3 or 4 directions.
detect=(temporary(detect) and 24B) < 1B
; zz=where(detect le 2,count)
; if count ne 0 then detect[zz]=0
IF disp THEN BEGIN
print,' DISP',format=fmt
setwin,1,xsize=nx,ysize=ny
skysclim,a,lowval,hival,meanval,sigma
lowval=meanval-8*sigma
hival=meanval+16*sigma
tv,bytscl(a,min=lowval,max=hival,top=!d.n_colors-1)
setwin,2,xsize=nx,ysize=ny
tvscl,detect
ENDIF
detect=median(detect,2)
detect[0:4,*]=0
detect[*,0:4]=0
detect[nx-5:nx-1,*]=0
detect[*,ny-5:ny-1]=0
IF disp THEN BEGIN
setwin,3,xsize=nx,ysize=ny
tvscl,detect
ENDIF
print,' LOC',format=fmt
if disp then setusym,-1
for j=0,ny-1 do begin
z=where(detect[*,j],count)
; z=where(detect[*,j] gt 1,count)
if count ne 0 then begin
for i=0,count-1 do begin
; if detect[z[i],j] gt 1 then begin
if detect[z[i],j] then begin
xmax0 = z[i]
ymax0 = j
boxm,a,xmax0,ymax0,gap,gap,xmax,ymax,/nocheck
repeat begin
i0=max([0,xmax0-gap])
i1=min([nx-1,xmax0+gap])
j0=max([0,ymax0-gap])
j1=min([ny-1,ymax0+gap])
detect[i0:i1,j0:j1]=0B
xmax0=xmax
ymax0=ymax
boxm,a,xmax0,ymax0,gap,gap,xmax,ymax,/nocheck
endrep until xmax0 eq xmax and ymax0 eq ymax
detect[xmax,ymax]=1B
IF disp THEN plots,xmax,ymax,psym=8,/device,symsize=2.0,color=64
endif
endfor
endif
endfor
print,' PHOT',format=fmt
findsrc_phot,fn,exttag,disp,detect,nx,ny,gain,a,binfac,exptime, $
sz,gap,objrad,maxphotsig,sigthresh,nobjs
if disp then setusym,1
print,' '+systime(),nobjs
ENDFOR ; extension loop
end