pro imaskf, image, mask, new = new, badpix = badpix, $
            cbox = cbox, fwhm = fwhm, maxshift = maxshift, $
            sigscl = sigscl, minval = minval, $
            radius = radius, minrad = minrad, maxrad = maxrad, $
            nocenter = nocenter, maxpix = maxpix, search = search, $
            silent = silent

;+
; semi-interactively make mask frame for *one* given image
;   0 = to be masked
;   1 = not to be masked
; use IMASK.PRO (a wrapper for this routine) to mask a stack of images.
;
; algorithm: use cursor to pick out objects in the frame, then routine
;       uses median and sigma from 'iterstat' along with IDL's SEARCH2D
;       routine to find all connected pixels above a minimum value.
;       It then finds the distance to the farthest pixel from
;       center point (chosen as the max pixel in a box around the
;       cursor) and constructs a circular mask with that distance
;       and center.
;
; to use: display image in window with desired stretch using tvh'
;       (do not REBIN the image! if use DISPLAY to show the image,
;       need to edit the 'cursor' command so it reads in /data coords)
;
; warning: 'mask' must be a TWO-DIMENSIONAL array 
;
; minval: user-specified minval (cancels use of sigscl)
; /new: initialize masks (default is to add to existing masks)
; /search: non-circular mask, just use SEARCH2D results (jagged mask)
;
; can switch between 3 types of masks by choosing 
; the relevant keyword parm to be non-zero, with others set to zero
;   1. fixed radius circular masks (RADIUS parm) - default
;   2. isophotal circular masks (SIGSCL > 0.0)
;   3. purely isophotal masks (SIGSCL < 0.0, or /search)
;
; HISTORY: 
; Written by M. Liu (UCB): 06/03/95 
; 08/23/95(MCL): changed mask array to integer
; 09/24/95(MCL): added 'radius' parameter - use mask of fixed radius
; 09/29/95(MCL): added 'maxrad' - maximum size for a mask
; 04/12/96(MCL): changed mask array to byte
; 05/02/96(MCL): added /nocenter option
; 08/19/96(MCL): user allowed to enter RADIUS using middle button
; 10/31/96(MCL): iterstat measurements done only on the good pixels
; 03/01/97(MCL): flip between SIGSCL & RADIUS by setting them to <=0
;                added 'minrad' (min radius for mask)
; 05/13/97(MCL): now able to handle bad pixel mask
; 10/03/97(MCL): added check to see if IMCENTERF blowing up
; 05/21/98(MCL): fixed bug when passing sigscl=0, and then trying to change
; 05/26/98(MCL): added 'fwhm' parm (previously, FWHM=4 always)
;                centroiding errors now signaled with 'CFLAG' variable
;                mask radius now sets minimum size for 'cbox' parm
;                added 'maxshift' so won't centroid too far from cursor
;                (should reduce sensitivity to bad pixels when centroiding)
; 02/17/99(MCL): added simple option to set /search while using IMASKF!
;                user sets SIGSCL<0.0 to indicate /search
;                (side effect - not possible to set SIGSCL to truly <0)
; 07/20/01(MCL): finally converted from data coords to device coords
;                instead of having a separate window showing the mask,
;                it now draws the mask right on top of the image
;-


if not(keyword_set(cbox)) then cbox = 5
if not(keyword_set(FWHM)) then FWHM = 3
if not(keyword_set(maxshift)) then maxshift = 5
if not(keyword_set(sigscl) or keyword_set(radius)) then sigscl = 0.7
if keyword_set(minrad) eq 0 then minrad = 2.
if keyword_set(maxrad) eq 0 then maxrad = 30.
if keyword_set(radius) then $
    if (radius gt maxrad) then maxrad =  radius
if (keyword_set(image) eq 0) then begin
	print, 'pro imaskf,image,mask,[new],[badpix=],'
        print, '           [cbox='+strc(cbox)+'],[fwhm='+strc(fwhm)+'],[maxshift='+strc(maxshift)+']'
        print, '           [sigscl='+strc(sigscl)+'],[minval=],'
        print, '           [radius=],[minrad='+strc(minrad)+'],[maxrad='+strc(maxrad)+'],'
        print, '           [nocenter],[maxpix],[search],'
        print, '           [silent]'
	retall
endif


; initialize mask options (SIGSCL is default)
if not(keyword_set(sigscl)) then sigscl = 0 $
else radius = 0


; initialize masks if desired
sz = size(image)
if keyword_set(new) or (n_elements(mask) eq 0) then begin
	print,'initializing masks'
	mask = bytarr(sz(1),sz(2)) + 1B
endif


;; set up display
;maskwin = 12
;imwin = !D.WINDOW
;if (maskwin eq imwin) then maskwin=13
;window,maskwin,xs=sz(1)>256,ys=sz(2)>256,title='IMASK: current mask frame'
;display2,mask,/notv
;wset,imwin


;set up for device or data coords
;if keyword_set(dev) then begin
;    dat = 0
;    dev = 1
;endif else begin
;    dat = 1
;    dev = 0
;endelse


; get image stats from evenly distributed subset of pixels
NPIX = 8192
ii = indgen(NPIX) * round((sz(1)*sz(2))/NPIX)
pix = image(ii)
pix = pix(where(mask(ii) ne 0))
iterstat, pix, istat, maxit = 2, /silent
;iterstat, image(where(mask ne 0)), istat, maxit = 2, /silent
imed = istat(2)
isig = istat(3)
minval = imed + abs(sigscl)* isig


; print masking info
if not(keyword_set(radius) or keyword_set(umin)) then begin
    message, 'current SIGSCL = '+strc(sigscl), /info
    message, 'current MINVAL = '+strc(minval), /info
endif else if keyword_set(radius) then begin
    message, 'using user specified RADIUS = '+strc(radius), /info
    minval =  min(image)
endif else begin
    minval = umin
    message, 'using user specified MINVAL = '+strc(minval), /info
endelse
maxval =  max(image)


; print info
print,'* index button to click on location of objects *'
if keyword_set(radius) then $
  print, '* middle button to change RADIUS (set <0 to use SIGSCL) *' $
else $
  print, '* middle button to change SIGSCL (set =0 to use RADIUS) *' 
print,'* ring button to quit '


repeat begin
    
  CFLAG = 1     ; error flag for centroiding

  ; if using mask of fixed radius, 
  ; centering box cannot be larger than the mask size
  if keyword_set(radius) then minbox = 2*RADIUS $
    else minbox = cbox

  cursor,x,y,/data,/wait
  ;cursor,x,y,/dev,/wait
  ;cursor,x,y,/wait, dev = dev, dat = dat

; mask object
  if (!err eq 1) then begin
      if (x ge 0) and (x lt sz(1)) and (y ge 0) and (y lt sz(2)) then begin
          
;         find max pixel location near cursor
          if keyword_set(nocenter) then begin
              xcen = x
              ycen = y
          endif else begin
;              imcenterf, image, x, y, xcen, ycen, /maxpix, /silent
              imcenterf, image, x, y, xcen, ycen, badpix = badpix, $
                cbox = cbox < minbox, bigbox = cbox < minbox, $
                maxpix = maxpix, /silent
              print, '   cursor at ', printcoo(x, y), $
                ', center at ', printcoo(round(xcen), round(ycen))

;	      check centering using Goddard 'cntrd.pro' function 
              if (xcen(0)-FWHM) ge 0 and (xcen(0)+FWHM) lt (sz(1)-1) and $
                (ycen(0)-FWHM) ge 0 and (ycen(0)+FWHM) lt (sz(2)-1) then begin
                  cntrd, image, x, y, xcen2, ycen2, FWHM
                  if (xcen(0)-(xcen2)+ycen(0)-ycen2 gt 1.0) then $
                    print, '  diff w/CNTRD > 1 pix = ', $
                    printcoo(xcen(0)-xcen2, ycen(0)-ycen2)
              endif

              if (n_elements(xcen) gt 1) then xcen = round(avg(xcen))
              if (n_elements(ycen) gt 1) then ycen = round(avg(ycen))
              wait, 0.2

              ; sanity check - sometimes IMCENTERF blows up
              if (xcen gt sz(1)-1) or (xcen lt 0) or $
                (ycen gt sz(2)-1) or (ycen lt 0) then begin
                  message, 'IMCENTERF gives bad centroid of '$
                    +printcoo(xcen, ycen)+', using CNTRD centroid of '$
                    +printcoo(xcen2, ycen2), /info
                  xcen = xcen2
                  ycen = ycen2
              endif
              if (xcen lt 0) or (ycen lt 0) then begin
                  message, 'bad centroiding, ignoring object', /info
                  CFLAG = 0
              endif

              ; check for excessively large shift from cursor location
              dx = xcen-x
              dy = ycen-y
              if (abs(dx) gt maxshift) or (abs(dy) gt maxshift) then $
                CFLAG = -1

          endelse 
          
; 	  create mask if max pixel is above minval,
;	  be careful to select only one pixel!
          if ((image(xcen,ycen) ge minval) or keyword_set(radius)) $
            and (CFLAG eq 1) then begin
              
	      xyouts,xcen,ycen,'M',/dev
;              xyouts,xcen,ycen,'M',dev = dev, dat = dat
              
;	      find all pixels above minval connected to max pixel
              if not(keyword_set(radius)) then begin
                  region = search2d(image,xcen,ycen,minval,maxval)
                  whereis,image,region,rx,ry
                  
;		  find distance from center pixel to farthest masked pixel,
;		  be careful to select only one pixel!
                  dd2 = (float(xcen-rx))^2. + (float(ycen-ry))^2.
                  mm = where(dd2 eq max(dd2))
                  if n_elements(mm) gt 1 then begin
                      mx = rx(round(avg(mm)))
                      my = ry(round(avg(mm)))
                  endif else begin
                      mx = rx(mm(0))
                      my = ry(mm(0))
                  endelse
                  dd = sqrt((xcen-mx)^2. + (ycen-my)^2.)
                  if (dd gt maxrad) then dd = maxrad
                  if (dd lt minrad) then dd = minrad
                  print,'   mask radius = ',strc(ceil(dd))
              endif else $
                dd = radius
              
;	      get a square region of size ceil(dd) whose values
;	      are the distance to the center of the region
              cdd = ceil(dd)
              dist_circle,circ,fix(2*cdd+1),fix(cdd),fix(cdd)
              
; 	      find the pixels in the region within (dd) of the center
              rr = where(circ le total(dd))
              whereis,circ,rr,cx,cy
              
;	      convert these pixel locations to locations in the image
              cx = cx + xcen - total(cdd)
              cy = cy + ycen - total(cdd)
              
;	      mask the region
              if keyword_set(search) then mask(region) = 0 $
              else mask(cx,cy) = 0
              
;	      update display of the mask: removed 07/20/01
	      ;wset,maskwin
	      ;display2,mask,/notv,/noerase
	      ;wset,imwin
              display2, image*mask*badpix
              
;              plots, ps = 3, cx, cy, color = 0, dev = dev, dat = dat
              
          endif else if (image(xcen,ycen) lt minval) then $
            print,'   object below sigma cutoff -> not masked' $
          else if (CFLAG eq 0) then $
            print,'   multiple max pixels near cursor -> no masking' $
          else if (CFLAG eq -1) then $
            print, '  centroid is too far from cursor location, delta = ', $
            printcoo(dx, dy), ' -> no masking'
          
      endif


  ; middle button: adjust mask parameters
  endif else if (!err eq 2) then $
    repeat begin
      if (radius gt 0) or (sigscl eq 0) then begin
          message, 'current RADIUS='+strc(radius), /info
          read, 'enter new RADIUS (<=0 to switch to SIGSCL) = ', radius
          if (radius gt maxrad) then begin
              message, 'warning! exceeded previous MAXRAD', /info
              maxrad = radius
          endif
      endif 
      if (radius le 0) then begin
          sigscl = float(sigscl)
          message, 'current SIGSCL='+strc(abs(sigscl))+$
            ', MINVAL='+strc(minval), /info
          if abs(sigscl gt 0) then  $
            message, 'using circular masks', /info  $
          else  $
            message, 'using isophotal masks (/search)', /info

          read, 'enter new SIGSCL (0 to switch to RADIUS, <0 to toggle isophotal masks) = ', sigscl
          minval = imed + abs(sigscl)*isig
          search = (sigscl lt 0)

          message, 'new SIGSCL='+strc(sigscl)+', MINVAL='+strc(minval), /info
          if not(keyword_set(search)) then  $
            message, 'using circular masks', /info  $
          else  $
            message, 'using isophotal masks (/search)', /info

      endif
  endrep until (radius gt 0) or (abs(sigscl) gt 0)


endrep until (!err eq 4)

;wdelete,maskwin

end