pro aoradnorm, img, outimg, radprof, $ medimg, radimg, $ medbox = MEDBOX, drad = DRAD, cbox = CBOX, rmin = RMIN, $ xcen = xcen, ycen = ycen, $ head = head, $ outfile = outfile, $ nofloor = nofloor, nodisp = nodisp ;+ ; 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. 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 ; radprof radial profile with avg flux & sigma in annulus ([3,n] array) ; ; 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 ; 08/07/01 MCL: writes measurements to outfile ; 09/28/01 MCL: now returns 'radprof' data, changed order of output variables ; computes 1-sigma flux ratio vs radius, saved to output ; (uses max pixel in input image, so use with caution) ; added /nodisp and /nofloor ;- BADVAL = -1e6 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, radprof,' print, ' [medbox=', strc(MEDBOX), '], [drad=', strc(DRAD), +$ '], [cbox=', strc(CBOX), '], [rmin=', strc(RMIN), ']' print, ' [xcen=], [ycen=],' print, ' [head=],' print, ' [outfile=]' print, ' [nofloor], [nodisp]' return endif ; initialize sz = size(img) if sz(0) ne 2 then $ message, 'input image is not a 2-d array!' imbad = (img ne BADVAL) ; create unsharp masked image medimg = imsub(img, filter_image(img, median = MEDBOX)) medbad = (medimg ne BADVAL) ;------------------------------------------------------------ ; 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 ne BADVAL and temporary(ddarr) ge RMIN) iterstat, medimg(w), istat, /silent, /nobad sigmin = istat(3) message, 'noise floor (in DN) = '+strc(sigmin), /info ; join the radial profile with the noise floor sig = radout(*, 7) if not(keyword_set(nofloor)) then begin w = where(sig le sigmin) sig(min(w):*) = sigmin endif ; plot radial noise profile and the noise floor if not(keyword_set(nodisp)) then begin win, 0, 450 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 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?')) then stop wdelete, 0 endif ;------------------------------------------------------------ ; construct radial profile using RADGEN.PRO radgen, radimg, radout(*, 0), sig, xcen, ycen, imsize = sz(1:2) ; and normalize outimg = imdiv(medimg, radimg) ; compute 1-sigma flux ratio relative to max pixel ; use with caution! (will be bogus for saturated images, of course) dflux = max(img)/sig ; store output array radprof = [[radout(*, 0)], [radout(*, 6)], [sig], [dflux]] if keyword_set(outfile) then begin if filecheck(outfile, /over) then begin openw, unit0, outfile, /get_lun printf, unit0, '# AORADNORM.PRO: ', systime(), userid() printf, unit0, '# flattened image using radial noise profile' printf, unit0, '# median filtering box = '+strc(MEDBOX) printf, unit0, '# centering box = '+strc(CBOX) printf, unit0, '# radial bin size = '+strc(DRAD) printf, unit0, '# object center '+printcoo(xcen, ycen) printf, unit0, '# max pixel in input image = ', strc(max(img)), ' DN' printf, unit0, '#' printf, unit0, '# r(pix) avg(flux) sig(flux) d(flux),1-sigma' printf, unit0 free_lun, unit0 printarray, radprof, out = outfile, /silent endif endif ; update history, if passed 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 end