pro aoradnorm, img, outimg, medimg, radimg,  $
               medbox = MEDBOX, drad = DRAD, cbox = CBOX, rmin = RMIN, $
               xcen = xcen, ycen = ycen,  $
               head = head,auto=auto

;+
; PURPOSE:
; Given an image (presumably an AO image), create an unsharp masked
; image. then compute a radial noise profile for it and divide by the
; profile.
;
; NOTES:
; This is a quick & dirty way to "flatten" the image and
; make the speckle noise look constant at all radii.
;
; uses /iterstat for RADPLOTF.PRO to avoid effect of bad pixels
; still need to be careful in the center as things don't work so well
; there (few numbers of pixels)
;
; INPUTS
;  img     original mosaic
;
; OUTPUTS
;  outimg  radial noise-normalized image
;  medimg  unsharp masked mosaic
;  radimg  image with radial profile used
; KEYWORDS
;  auto		don't ask the user any questions
;    
; uses RADPLOTF.PRO
;
; HISTORY: Written by M. Liu (IfA/UH) 05/21/01
; 06/11/01 MCL: now uses ITERSTAT.PRO to compute std dev
; 2002-07-02 MDP: Added /auto
;-


if not(keyword_set(medbox)) then MEDBOX = 11
if not(keyword_set(drad)) then DRAD = 2
if not(keyword_set(cbox)) then CBOX = 11
if not(keyword_set(rmin)) then RMIN = 50
if n_params() lt 2 then begin
    print, 'pro aoradnorm, img, outimg, medimg, radimg,'
    print, '               [medbox=', strc(MEDBOX), '], [drad=', strc(DRAD), +$
      '], [cbox=', strc(CBOX), '], [rmin=', strc(RMIN), ']'
    return
endif    


; initialize
sz = size(img)
if sz(0) ne 2 then  $
  message, 'input image is not a 2-d array!'
imbad = (img eq img)


; create unsharp masked image
medimg = imsub(img, filter_image(img, median = MEDBOX))
;stop
medbad = (medimg eq medimg)


;------------------------------------------------------------
; interactively find center if needed
if not(keyword_set(xcen)) and not(keyword_set(ycen)) then begin
    message, 'no image center passed, running IREG', /info
    ireg, img, xcen, ycen, fwhm = FWHM, cbox = CBOX, badpix = imbad
endif

; print, 'please use mouse button to indicate the center'
; cursor, xcur, ycur, /data, /wait
; FWHM = 4
; iregf, img, 1B, imbad,  $
;   xcur, ycur, xcen, ycen, CBOX, FWHM

; ; overplot centroid location (data coords)
; oplot, [xcen], [ycen], psym = 1, syms = 6, color = 1
; tvbox, cbox, xcen, ycen, 0, /data

; ; show zoom of centroid region
; win, 2
; cut = imcut(medimg, 4*CBOX+1, round(xcen), round(ycen), /fix, /silent)
; bpcut = imcut(medbad, 4*CBOX+1,  $
;               round(xcen), round(ycen), /fix, /silent)
; nbad = total(bpcut eq 0)
; if (nbad gt 400) then  $
;   fixpix, cut, bpcut, cutfix, /silent  $
; else  $
;   cutfix = cut
; display2, cutfix, /tv, xtit = '4*CBOX + 1', ytit = '4*CBOX + 1'
; if not(getyn('satisfied with centering?')) then stop


;------------------------------------------------------------
; extract radial profile
; assumes image has been sky subtracted
radplotf, medimg, xcen, ycen, radout,  $
  inrad = 0.25, outrad = max([xcen, ycen, sz(1)-xcen, sz(2)-ycen]), $
  drad = DRAD, badpix = medbad, sky = 1e-20, /silent, /iter

; fix central pixel of radial profile, which will have sigma=0
radout(0, 7) = radout(1, 7)

; compute noise floor, excluding a region around the star
dist_circle, ddarr, sz(1:2), xcen, ycen
w = where(medimg eq medimg and temporary(ddarr) ge RMIN)
iterstat, medimg(w), istat, /silent, /nobad
sigmin = istat(3)
message, 'noise floor (in DN) = '+strc(sigmin), /info

; plot radial noise profile and the noise floor
win, 0, 600
plot, radout(*, 0), radout(*, 7), /ylog,  $
  xtit = 'radius (pixels)', ytit = 'sigma (DN)',  $
  ps = 1, chars = 1.5, sym = 0.5, /xs
plots, !x.crange, [sigmin, sigmin], line = 1

; join the radial profile with the noise floor
sig = radout(*, 7)
w = where(sig le sigmin)
sig(min(w):*) = sigmin
oplot, radout(*, 0), sig

legend, ['radial floor', 'noise floor', 'final'],  $
  line = [1, 5, 0], box = 0, chars = 1.5, /right

if not(getyn('satisfied with noise profile?',noquery=auto)) then stop

;------------------------------------------------------------
; construct radial profile using RADGEN.PRO 
radgen, radimg, radout(*, 0), sig, xcen, ycen, imsize = sz(1:2)

; and normalize
outimg = medimg / radimg

;stop

if n_elements(head) gt 0 then begin
    sxaddhist, 'AORADNORM.PRO: '+systime()+userid(), head
    sxaddhist, 'flattened image using radial noise profile', head
    sxaddhist, '  median filtering box = '+strc(MEDBOX), head
    sxaddhist, '  centering box = '+strc(CBOX), head
    sxaddhist, '  radial bin size = '+strc(DRAD), head
    sxaddhist, '  object center '+printcoo(xcen, ycen), head
endif

;wdelete, 0

end