;+ ; 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,osz,binfac,info, $ 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,info.exptime,x,y,objrad,objrad+1.0,objrad+5.0,/nolog,/silent, $ xcen=xc,ycen=yc,mag=mag,err=err,fwhm=fwhm,max=max,flux=flux, $ skymean=sky,skyerr=skyerr,skysig=skysig,boxmrad=-1*objrad crsdet = (max-sky)/(flux/gain*info.exptime) snr = (max-sky)/skysig zg = where(fwhm gt 1.0 and mag lt 30.0 and $ (crsdet ge 0.0 and crsdet le 0.25 and finite(crsdet)) and $ snr ge sigthresh, 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] snr = snr[zg] sky = sky[zg] skysig= skysig[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 snr = 0.000 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 robomean,fwhm*binfac,3.0,0.5,avgfwhm robomean,sky,3.0,0.5,avgsky robomean,skysig,3.0,0.5,skysg zo=where(a gt avgsky+skysg*50.0,countz) obscura1= float(countz)/float(n_elements(a)) zo=where(a gt avgsky+skysg*5.0,countz) obscura2= float(countz)/float(n_elements(a)) data=[[xc*binfac],[yc*binfac],[fwhm*binfac],[mag],[err],[snr]] mkhdr,hdr,data sxaddpar,hdr,'OBJECT',info.object,' Object name' sxaddpar,hdr,'AIRMASS',info.airmass,' Airmass of observation' sxaddpar,hdr,'XSIZE',osz[0],' Original x-size of image' sxaddpar,hdr,'YSIZE',osz[1],' Original y-size of image' 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^2),' Gain of CCD in e-/DN' sxaddpar,hdr,'BINFAC',binfac,' Software binning factor used' sxaddpar,hdr,'EXPTIME',info.exptime,' Exposure time in seconds' sxaddpar,hdr,'MAXSIG',maxphotsig,' Saturated above this DN level' sxaddpar,hdr,'MEANFWHM',avgfwhm,' Mean FWHM in pixels' sxaddpar,hdr,'SKYLEVEL',avgsky,' Average sky signal level counts/pixel' sxaddpar,hdr,'SKYSIGMA',skysg,' Standard deviation of the sky signal' sxaddpar,hdr,'OBSCURA1',obscura1,' Fraction of image obscured by 50*skysig bright pixels' sxaddpar,hdr,'OBSCURA2',obscura2,' Fraction of image obscured by 5*skysig bright pixels' writefits,fn+'.src'+exttag,data,hdr print,'('+strn(avgfwhm,format='(f10.1)')+') s'+ $ strn(avgsky,format='(i5)',length=5)+' +- '+ $ strn(skysg,format='(i4)',length=4)+' '+ $ strn(obscura1,format='(f5.3)')+' '+ $ strn(obscura2,format='(f5.3)'), $ format='($,a)' ; send information to log file fnlog = 'info'+exttag+'.log' tag = fn+exttag info = strn(info.object,padtype=1,length=10) + ' ' + $ strn(info.exptime,length=6,format='(f6.1)') + ' ' + $ strn(info.airmass,length=4,format='(f4.2)') + ' ' + $ strn(avgfwhm,length=5,format='(f5.2)') + ' ' + $ strn(avgsky,length=5,format='(i5)') + ' '+ $ strn(skysg,length=4,format='(i4)') + ' '+ $ strn(nobj,length=5,format='(i5)') + ' ' + $ strn(obscura1,format='(f5.3)') + ' ' + $ strn(obscura2,format='(f5.3)') repwrite,fnlog,tag,tag+info 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. if exists(path+fn+'.fits') then ft='.fits' else ft='' hdr=headfits(path+fn+ft) 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) fits_read,path+fn+ft,a,hdr,exten_no=extlist[ext] 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 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,[arrsz[1],arrsz[2]], $ binfac,info,sz,gap,objrad,maxphotsig,sigthresh,nobjs if disp then setusym,1 print,' '+systime(),nobjs ENDFOR ; extension loop end