pro find, image, x, y, flux, sharp, roundness, hmin, fwhm, roundlim, sharplim,$ PRINT = print, SILENT=silent, MONITOR= monitor ;+ ; NAME: ; FIND ; PURPOSE: ; Find positive brightness perturbations (i.e stars) in an image ; EXPLANATION: ; Also returns centroids and shape parameters (roundness & sharpness). ; Adapted from 1991 version of DAOPHOT, but does not allow for bad pixels ; and uses a slightly different centroid algorithm. ; ; Modified in March 2008 to use marginal Gaussian fits to find centroids ; CALLING SEQUENCE: ; FIND, image, [ x, y, flux, sharp, round, hmin, fwhm, roundlim, sharplim ; PRINT= , /SILENT ] ; ; INPUTS: ; image - 2 dimensional image array (integer or real) for which one ; wishes to identify the stars present ; ; OPTIONAL INPUTS: ; FIND will prompt for these parameters if not supplied ; ; hmin - Threshold intensity for a point source - should generally ; be 3 or 4 sigma above background RMS ; fwhm - FWHM to be used in the convolve filter ; sharplim - 2 element vector giving low and high cutoff for the ; sharpness statistic (Default: [0.2,1.0] ). Change this ; default only if the stars have significantly larger or ; or smaller concentration than a Gaussian ; roundlim - 2 element vector giving low and high cutoff for the ; roundness statistic (Default: [-1.0,1.0] ). Change this ; default only if the stars are significantly elongated. ; ; OPTIONAL INPUT KEYWORDS: ; /MONITOR - Normally, FIND will display the results for each star ; only if no output variables are supplied. Set /MONITOR ; to always see the result of each individual star. ; /SILENT - set /SILENT keyword to suppress all output display ; PRINT - if set and non-zero then FIND will also write its results to ; a file find.prt. Also one can specify a different output file ; name by setting PRINT = 'filename'. ; ; OPTIONAL OUTPUTS: ; x - vector containing x position of all stars identified by FIND ; y- vector containing y position of all stars identified by FIND ; flux - vector containing flux of identified stars as determined ; by a Gaussian fit. Fluxes are NOT converted to magnitudes. ; sharp - vector containing sharpness statistic for identified stars ; round - vector containing roundness statistic for identified stars ; ; NOTES: ; (1) The sharpness statistic compares the central pixel to the mean of ; the surrounding pixels. If this difference is greater than the ; originally estimated height of the Gaussian or less than 0.2 the height of the ; Gaussian (for the default values of SHARPLIM) then the star will be ; rejected. ; ; (2) More recent versions of FIND in DAOPHOT allow the possibility of ; ignoring bad pixels. Unfortunately, to implement this in IDL ; would preclude the vectorization made possible with the CONVOL function ; and would run extremely slowly. ; ; (3) Modified in March 2008 to use marginal Gaussian distributions to ; compute centroid. (Formerly, find.pro determined centroids by locating ; where derivatives went to zero -- see cntrd.pro for this algorithm. ; This was the method used in very old (~1984) versions of DAOPHOT. ) ; As discussed in more detail in the comments to the code, the centroid ; computation here is the same as in IRAF DAOFIND but differs slightly ; from the current DAOPHOT. ; PROCEDURE CALLS: ; GETOPT() ; REVISION HISTORY: ; Written W. Landsman, STX February, 1987 ; ROUND now an internal function in V3.1 W. Landsman July 1993 ; Change variable name DERIV to DERIVAT W. Landsman Feb. 1996 ; Use /PRINT keyword instead of TEXTOUT W. Landsman May 1996 ; Changed loop indices to type LONG W. Landsman Aug. 1997 ; Replace DATATYPE() with size(/TNAME) W. Landsman Nov. 2001 ; Fix problem when PRINT= filename W. Landsman October 2002 ; Fix problems with >32767 stars D. Schlegel/W. Landsman Sep. 2004 ; Fix error message when no stars found S. Carey/W. Landsman Sep 2007 ; Rewrite centroid computation to use marginal Gaussians W. Landsman ; Mar 2008 ; Added Monitor keyword, /SILENT now suppresses all output ; W. Landsman Nov 2008 ;- ; On_error,2 ;Return to caller compile_opt idl2 npar = N_params() if npar EQ 0 then begin print,'Syntax - FIND, image,' + $ '[ x, y, flux, sharp, round, hmin, fwhm, roundlim, sharplim' print,' PRINT= , /SILENT ]' return endif ;Determine if hardcopy output is desired doprint = keyword_set( PRINT) silent = keyword_set( SILENT ) if N_elements(monitor) EQ 0 then $ monitor = (not silent) and (not arg_present(flux) ) maxbox = 13 ;Maximum size of convolution box in pixels ; Get information about the input image type = size(image) if ( type[0] NE 2 ) then message, $ 'ERROR - Image array (first parameter) must be 2 dimensional' n_x = type[1] & n_y = type[2] message, NoPrint=Silent, $ 'Input Image Size is '+strtrim(n_x,2) + ' by '+ strtrim(n_y,2),/INF if ( N_elements(fwhm) NE 1 ) then $ read, 'Enter approximate FWHM: ', fwhm radius = 0.637*FWHM > 2.001 ;Radius is 1.5 sigma radsq = radius^2 nhalf = fix(radius) < (maxbox-1)/2 ; nbox = 2*nhalf + 1 ;# of pixels in side of convolution box middle = nhalf ;Index of central pixel lastro = n_x - nhalf lastcl = n_y - nhalf sigsq = ( fwhm/2.35482 )^2 mask = bytarr( nbox, nbox ) ;Mask identifies valid pixels in convolution box g = fltarr( nbox, nbox ) ;g will contain Gaussian convolution kernel dd = indgen(nbox-1) + 0.5 - middle ;Constants need to compute ROUND dd2 = dd^2 row2 = (findgen(Nbox)-nhalf)^2 for i = 0, nhalf do begin temp = row2 + i^2 g[0,nhalf-i] = temp g[0,nhalf+i] = temp endfor mask = fix(g LE radsq) ;MASK is complementary to SKIP in Stetson's Fortran good = where( mask, pixels) ;Value of c are now equal to distance to center ; Compute quantities for centroid computations that can be used for all stars g = exp(-0.5*g/sigsq) ; In fitting Gaussians to the marginal sums, pixels will arbitrarily be ; assigned weights ranging from unity at the corners of the box to ; NHALF^2 at the center (e.g. if NBOX = 5 or 7, the weights will be ; ; 1 2 3 4 3 2 1 ; 1 2 3 2 1 2 4 6 8 6 4 2 ; 2 4 6 4 2 3 6 9 12 9 6 3 ; 3 6 9 6 3 4 8 12 16 12 8 4 ; 2 4 6 4 2 3 6 9 12 9 6 3 ; 1 2 3 2 1 2 4 6 8 6 4 2 ; 1 2 3 4 3 2 1 ; ; respectively). This is done to desensitize the derived parameters to ; possible neighboring, brighter stars. xwt = fltarr(nbox,nbox) wt = nhalf - abs(findgen(nbox)-nhalf ) + 1 for i=0,nbox-1 do xwt[0,i] = wt ywt = transpose(xwt) sgx = total(g*xwt,1) p = total(wt) sgy = total(g*ywt,2) sumgx = total(wt*sgy) sumgy = total(wt*sgx) sumgsqy = total(wt*sgy*sgy) sumgsqx = total(wt*sgx*sgx) vec = nhalf - findgen(nbox) dgdx = sgy*vec dgdy = sgx*vec sdgdxs = total(wt*dgdx^2) sdgdx = total(wt*dgdx) sdgdys = total(wt*dgdy^2) sdgdy = total(wt*dgdy) sgdgdx = total(wt*sgy*dgdx) sgdgdy = total(wt*sgx*dgdy) c = g*mask ;Convolution kernel now in c sumc = total(c) sumcsq = total(c^2) - sumc^2/pixels sumc = sumc/pixels c[good] = (c[good] - sumc)/sumcsq c1 = exp(-.5*row2/sigsq) sumc1 = total(c1)/nbox sumc1sq = total(c1^2) - sumc1 c1 = (c1-sumc1)/sumc1sq message,/INF,Noprint=Silent, $ 'RELATIVE ERROR computed from FWHM ' + strtrim(sqrt(total(c[good]^2)),2) if N_elements(hmin) NE 1 then read, $ 'Enter minimum value above background for threshold detection: ',hmin if N_elements(sharplim) NE 2 then begin print,'Enter low and high cutoffs, press [RETURN] for defaults:' GETSHARP: ans = '' read, 'Image Sharpness Statistic (DEFAULT = 0.2,1.0): ', ans if ans EQ '' then sharplim = [0.2,1.0] else begin sharplim = getopt(ans,'F') if N_elements(sharplim) NE 2 then begin message, 'ERROR - Expecting 2 scalar values',/CON goto, GETSHARP endif endelse GETROUND: ans = '' read, 'Image Roundness Statistic [DEFAULT = -1.0,1.0]: ',ans if ans EQ '' then roundlim = [-1.,1.] else begin roundlim = getopt( ans, 'F' ) if N_elements( roundlim ) NE 2 then begin message,'ERROR - Expecting 2 scalar values',/CON goto, GETROUND endif endelse endif message,'Beginning convolution of image', /INF, NoPrint=Silent h = convol(float(image),c) ;Convolve image with kernel "c" h[0:nhalf-1,*] = 0 & h[n_x-nhalf:n_x-1,*] = 0 h[*,0:nhalf-1] = 0 & h[*,n_y-nhalf:n_y-1] = 0 message,'Finished convolution of image', /INF, NoPrint=Silent mask[middle,middle] = 0 ;From now on we exclude the central pixel pixels = pixels -1 ;so the number of valid pixels is reduced by 1 good = where(mask) ;"good" identifies position of valid pixels xx= (good mod nbox) - middle ;x and y coordinate of valid pixels yy = fix(good/nbox) - middle ;relative to the center offset = yy*n_x + xx SEARCH: ;Threshold dependent search begins here index = where( h GE hmin, nfound) ;Valid image pixels are greater than hmin if nfound EQ 0 then begin ;Any maxima found? message,'ERROR - No maxima exceed input threshold of ' + $ string(hmin,'(F9.1)'),/CON goto,FINISH endif for i= 0L, pixels-1 do begin stars = where (h[index] GE h[index+offset[i]], nfound) if nfound EQ 0 then begin ;Do valid local maxima exist? message,'ERROR - No maxima exceed input threshold of ' + $ string(hmin,'(F9.1)'),/CON goto,FINISH endif index = index[stars] endfor ix = index mod n_x ;X index of local maxima iy = index/n_x ;Y index of local maxima ngood = N_elements(index) message,/INF,Noprint=Silent, $ strtrim(ngood,2)+' local maxima located above threshold' nstar = 0L ;NSTAR counts all stars meeting selection criteria badround = 0L & badsharp=0L & badcntrd=0L if (npar GE 2) or (doprint) then begin ;Create output X and Y arrays? x = fltarr(ngood) & y = x endif if (npar GE 4) or (doprint) then begin ;Create output flux,sharpness arrays? flux = x & sharp = x & roundness = x endif if doprint then begin ;Create output file? if ( size(print,/TNAME) NE 'STRING' ) then file = 'find.prt' $ else file = print message,'Results will be written to a file ' + file,/INF,Noprint=Silent openw,lun,file,/GET_LUN printf,lun,' Program: FIND '+ systime() printf,lun,format='(/A,F7.1)',' Threshold above background:',hmin printf,lun,' Approximate FWHM:',fwhm printf,lun,format='(2(A,F6.2))',' Sharpness Limits: Low', $ sharplim[0], ' High',sharplim[1] printf,lun,format='(2(A,F6.2))',' Roundness Limits: Low', $ roundlim[0],' High',roundlim[1] printf,lun,format='(/A,i6)',' No of sources above threshold',ngood endif if (not SILENT) and MONITOR then $ print,format='(/8x,a)',' STAR X Y FLUX SHARP ROUND' ; Loop over star positions; compute statistics for i = 0L,ngood-1 do begin temp = float(image[ix[i]-nhalf:ix[i]+nhalf,iy[i]-nhalf:iy[i]+nhalf]) d = h[ix[i],iy[i]] ;"d" is actual pixel intensity ; Compute Sharpness statistic sharp1 = (temp[middle,middle] - (total(mask*temp))/pixels)/d if ( sharp1 LT sharplim[0] ) or ( sharp1 GT sharplim[1] ) then begin badsharp = badsharp + 1 goto, REJECT ;Does not meet sharpness criteria endif ; Compute Roundness statistic dx = total( total(temp,2)*c1) dy = total( total(temp,1)*c1) if (dx LE 0) or (dy LE 0) then begin badround = badround + 1 goto, REJECT ;Cannot compute roundness endif around = 2*(dx-dy) / ( dx + dy ) ;Roundness statistic if ( around LT roundlim[0] ) or ( around GT roundlim[1] ) then begin badround = badround + 1 goto,REJECT ;Does not meet roundness criteria endif ; ; Centroid computation: The centroid computation was modified in Mar 2008 and ; now differs from DAOPHOT which multiplies the correction dx by 1/(1+abs(dx)). ; The DAOPHOT method is more robust (e.g. two different sources will not merge) ; especially in a package where the centroid will be subsequently be ; redetermined using PSF fitting. However, it is less accurate, and introduces ; biases in the centroid histogram. The change here is the same made in the ; IRAF DAOFIND routine (see ; http://iraf.net/article.php?story=7211&query=daofind ) ; sd = total(temp*ywt,2) sumgd = total(wt*sgy*sd) sumd = total(wt*sd) sddgdx = total(wt*sd*dgdx) hx = (sumgd - sumgx*sumd/p) / (sumgsqy - sumgx^2/p) ; HX is the height of the best-fitting marginal Gaussian. If this is not ; positive then the centroid does not make sense if (hx LE 0) then begin badcntrd = badcntrd + 1 goto, REJECT endif skylvl = (sumd - hx*sumgx)/p dx = (sgdgdx - (sddgdx-sdgdx*(hx*sumgx + skylvl*p)))/(hx*sdgdxs/sigsq) if abs(dx) GE nhalf then begin badcntrd = badcntrd + 1 goto, REJECT endif xcen = ix[i] + dx ;X centroid in original array ; Find Y centroid sd = total(temp*xwt,1) sumgd = total(wt*sgx*sd) sumd = total(wt*sd) sddgdy = total(wt*sd*dgdy) hy = (sumgd - sumgy*sumd/p) / (sumgsqx - sumgy^2/p) if (hy LE 0) then begin badcntrd = badcntrd + 1 goto, REJECT endif skylvl = (sumd - hy*sumgy)/p dy = (sgdgdy - (sddgdy-sdgdy*(hy*sumgy + skylvl*p)))/(hy*sdgdys/sigsq) if abs(dy) GE nhalf then begin badcntrd = badcntrd + 1 goto, REJECT endif ycen = iy[i] +dy ;Y centroid in original array ; This star has met all selection criteria. Print out and save results if monitor then $ print,FORM = '(12x,i5,2f7.1,f9.1,2f9.2)', $ nstar, xcen, ycen, d, sharp1, around if (npar GE 2) or (doprint) then begin x[nstar] = xcen & y[nstar] = ycen endif if ( npar GE 4 ) or (doprint) then begin flux[nstar] = d & sharp[nstar] = sharp1 & roundness[nstar] = around endif nstar = nstar+1 REJECT: endfor nstar = nstar-1 ;NSTAR is now the index of last star found if doprint then begin printf,lun,' No. of sources rejected by SHARPNESS criteria',badsharp printf,lun,' No. of sources rejected by ROUNDNESS criteria',badround printf,lun,' No. of sources rejected by CENTROID criteria',badcntrd endif if (not SILENT) and (MONITOR) then begin print,' No. of sources rejected by SHARPNESS criteria',badsharp print,' No. of sources rejected by ROUNDNESS criteria',badround print,' No. of sources rejected by CENTROID criteria',badcntrd endif if nstar LT 0 then return ;Any stars found? if (npar GE 2) or (doprint) then begin x=x[0:nstar] & y = y[0:nstar] endif if (npar GE 4) or (doprint) then begin flux= flux[0:nstar] & sharp=sharp[0:nstar] roundness = roundness[0:nstar] endif if doprint then begin printf,lun, $ format = '(/8x,a)',' STAR X Y FLUX SHARP ROUND' for i = 0L, nstar do $ printf,lun,format='(12x,i5,2f8.2,f9.1,2f9.2)', $ i+1, x[i], y[i], flux[i], sharp[i], roundness[i] free_lun, lun endif FINISH: if SILENT or (not MONITOR) then return print,form='(A,F8.1)',' Threshold above background for this pass was',hmin ans = '' read,'Enter new threshold or [RETURN] to exit: ',ans ans = getopt(ans,'F') if ans GT 0. then begin hmin = ans goto, SEARCH endif return end