find.pro

basic median filter simulation

    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