pro mosf, imgs0, xobj0, yobj0, out, outexp, sky, xf, yf, $ badpix=badpix, masks=masks0, $ wtmap = wtmap, expmap = expmap, exptimes = exptimes0, $ trim=trim, $ nrej=nrej, median=median, $ setsky=setsky, minsky=minsky, submed=submed, $ goodlist = goodlist0, $ maxsize = MAXSIZE silent = silent ;+ ; routine to shift and stack a set of images to create a final mosaic, ; using masks to exclude bad pixels or cosmic rays ; ; INPUTS: ; imgs 3d array of images (all the same size) ; xobj, yobj position of registration object in each frame ; (need not be integer, but non-integer shifts ; will be rounded) ; ; OUTPUTS: ; out resulting mosaic ; outexp exposure map for resulting mosaic ; sky sky values determined for each frame by /setsky ; xf lower + upper x coords for each image's position ; in the mosiac (2 by N array where N is number ; of images) -> for debugging ; yf y coords for each image's position in the mosaic ; pieces indiv shifted (and maybe trimmed) images ; which go into the mosiac - buffered by BADVAL, ; bad pixels are turned either to BADVAL, ; and with sky offset from /setsky added ; ; KEYWORD INPUTS: ; badpix badpix 2d image (0=bad, not 0=good) (optional) ; masks masks for each individual image (3d array) ; wtmap map of weights for combining images (3-d array) ; final mosaic will be scaled to wtmap=1.0 ; (for images of diff int. times, optimal ; weights are the int. times) ; expmap map of exposure times for combining images (3-d array) ; final mosaic will be scaled to expmap=1.0 ; exptimes 1-d array of exp times for combining images ; final mosaic will be scaled to expmap=1.0 ; /trim only mosaic the overlapping regions ; /setsky iteratively adds a constant to frames to ; minimize the 'seams' of the mosaic ; (adds considerable time but often necessary, ; will crash unless ALL the images overlap) ; nrej number of images to reject at high/low ends ; at a given pixel location before combining, ; i.e., a quick cosmic ray deterrent ; (if not enough images at a gievn pixel, ; then just averages) ; /median take the median value of the good images ; at each pixel ; minval min values to stop iterations of /setsky ; (0.001 default) ; /submed set median of each image to zero before ; combining ; goodlist logical list of which images to use (0=bad,!0=good) ; very useful if passing a large number of ; images so don't have to pass, eg, imgs(*,*,10:150) ; ; CAVEATS: ; - does only integer shifts ; - due to crappy programming, need to be VERY CAREFUL with ; variables. most of the working variables used to make the ; mosaic (e.g. xobj, yobj, sky) can have fewer entries than the ; input variables (e.g. imgs, xobj0, yobj0, exptimes0, etc.). ; this occurs when some images are flagged as bad, via ; 'goodlist' or shifts=-999. ; the mapping btwn the two is kept in the 'wg' index ; variable. when doing FOR loops, the 'i' index loops over only the ; good subset while the 'ii' index loops over all images ; ; HISTORY: ; Written by M. C. Liu (UCB): 12/14/94 (spun off of 'ireg.pro') ; 06/05/95 (MCL): added high/low pixel rejection ability ; 07/27/95 (MCL): improved bad pixel treatment, ; corrected bug in pixel rejection scheme ; 08/20/95 (MCL): adjusted algorithm for greater speed, ; made sure rounding of offsets is done properly ; 04/16/96 (MCL): able to pass masks for the individual frames ; badpix and masks set to BYTE type ; exposure map (outexp) set to INT type ; changed the assignment of BADVAL to blank areas and ; also 'eq BADVAL' to 'le/gt BADVAL' (round-off worries) ; made sky offsets total up to 0 ; added /submed ; changed /setsky so it uses the median of the ; overlapping regions instead of the average ; to reduce sensitivity to unflagged bad pixels ; 11/30/96 (MCL): corrected subtle error in determining final image ; size using rounding when the max image offset ends ; in 0.5 (only when using /trim) ; 12/04/96 (MCL): big changes from MOSF.PRO - less memory intensive! ; developed under the name MOSFQUICK.PRO ; eliminated the 'pieces' variable ; disabled (for now) all the filtering schemes (nrej, etc) ; for combining the mosaic ; added 'goodlist' feature ; changed 'le/gt BADVAL' back to 'ne BADVAL' ; increased MINSKY from 1e-5 (excessive) to 1e-3 DN ; checked against MOSF.PRO, ok to w/in roundoff errors ; 02/10/97 (MCL): added /silent ; 02/28/97 (MCL): replaced existing version of MOSF.PRO with MOSFQUICK.PRO ; previous version renamed MOSF_SLOW.PRO (still useful) ; 04/24/98 (MCL): added 'expmap', a 3d-array of pixel exposure time ; for int time weighting when combining diff mosaics ; final mosaic will be normalized to expmap=1.0 ; NOTE: /setsky doesn't handle the expmap (yet) ; 10/25/98 (MCL): renamed 'expmap' parm into 'wtmap' (map of weights) ; b/c this is treated like weights, not diff int. times ; *NOTE: therefore, past use of 'expmap' gave *WRONG* ; phot, but did do exposure time (variance) wts! ; added new 'expmap' variable which does behave right ; small bug fix: /submed wouldn't work if no /setsky ; 01/29/99 (MCL): added MAXSIZE limit on output array size (sanity check) ; 03/13/99 (MCL): added 'exptimes' - untested!!! ; 09/29/00 (MCL): added kludge fix for /setsky if some images don't overlap ; now identifies images to be ignored (shift=-999) ; very convenient, should have done this earlier ; *** placed under RCS, version 1.1 *** ; 10/02/00 (MCL): bug fix - wasn't using masks when some shifts=-999 ; *** RCS, version 1.2 (05/10/01) *** ; 05/10/01 (MCL): better memory management for 'wtmap' and 'expmap' ; (code is slightly uglier, but takes less memory) ; *** RCS, version 1.3 *** ; 05/12/01 (MCL): MAJOR CHANGES ; 1. much better memory usage, no longer makes extra ; copies of the input images & masks! however, the coding ; is slightly uglier (use of both 'i' and 'ii' indicies!) ; 2. recognizes BADVAL pixels in input images, no need to pass masks! ; fairly extensive testing against previous version 1.3 (at least for ; most commonly used options: didn't test /trim, /submed) ; minor changes: increased use of MESSSAGE, instead of PRINT ; 05/13/01 (MCL): found bug - I think "shifts=-999" doesn't work when user ; doesn't pass 'goodlist'. never encountered this before! ; tried to fix bug, but got bogged down, then tried to ; return to previous tested version. think I did this ; successfully (re-ran tests from yesterday) ; for now, temporarily disable use of 'shifts=-999' ; -> still good to be cautious with this ... <- ; *** RCS version 1.4 *** ; 05/13/01 (MCL): above bug is more extensive than I thought - ; previous version can't use 'goodlist' at all! ; fixed implementation of 'goodlist' (should *always* ; use this instead of passing arrays like im[*,*,10:20]) ; disabled shifts=-999 for now ; 05/24/01 (MCL): adjusted disabling of shifts=-999, still turned off ; but now routine will run if 'goodlist' is passed ; 07/03/01 (MCL): big fix - /submed was adding median value, instead ; of subtracting it from the indiv images! ; ** version 1.6 ** ; ; -> modify /setsky so will take information in the weight map (wtmap) ;- ; empty spaces of the mosiac and mosaic pieces ; get filled with BADVAL BADVAL = -1e6 ; for /setsky, the largest sky offset must be larger ; than MINSKY to continue iterating on the sky value if not(keyword_set(minsky)) then MINSKY = 0.001 ; max dimension length of output array if not(keyword_set(MAXSIZE)) then MAXSIZE = 2000 if not(keyword_set(nrej)) then nrej = 0 if n_params() lt 4 then begin print, 'mosf, imgs, xobj, yobj, out, outexp, sky, xf, yf,' print, ' [badpix=], [masks=], [wtmap=], [expmap=], [exptimes=]' print, ' [trim],' print, ' [nrej=0], [median],' print, ' [setsky], [minsky=', strc(minsky), '], [submed],' print, ' [goodlist=],' print, ' [maxsize=],'+strc(MAXSIZE) print, ' [silent]' retall endif ; sanity checks sz = size(imgs0) ntot = sz(3) ;n = sz(3) if not(keyword_set(badpix)) then begin if not(keyword_set(silent)) then $ message,'no bad pixel mask passed', /info badpix = bytarr(sz(1),sz(2)) + 1B endif if total(abs((size(badpix))[1:2] - sz[1:2])) ne 0 then $ message, 'badpix mask and images sizes disagree!' ; innitialize goodlist if keyword_set(goodlist0) then begin if n_elements(goodlist0) ne sz(3) then $ message, 'goodlist does not have right number of entries!' goodlist = goodlist0 endif else $ goodlist = bytarr(sz(3))+1B ; check for images to be ignored (shift = -999 and/or goodlist = 0) wg = where(xobj0 ne -999 and goodlist ne 0, ng) wb = where(xobj0 eq -999 or goodlist eq 0, nb) w = where(xobj0 eq -999, nw) if (nw ge 1) and n_elements(goodlist) eq 0 then $ message, 'shifts=-999 is disabled b/c of possible bug' ; 05/13/01 ;message, strc(nb)+' out of the '+strc(ng)+ $ ;' input images will be ignored (shift = -999)', /info if not(keyword_set(silent)) then $ message, 'mosaicking '+ strc(ng)+ ' frames', /info ; initialize: need to be CAREFUL with variables ; most of the working variables will have 'ng' items (i.e. # of good images) ; i.e. potentially fewer entries than the input variables ; (e.g. xobj0, yobj0, etc.) ; the mapping btwn these two is kept in the 'wg' index variable xobj = xobj0(wg) yobj = yobj0(wg) mm = fltarr(ng) ; image medians ;imgs = imgs0(*, *, wg) ;if keyword_set(masks0) then masks = masks0(*, *, wg) ; however, (xf,yf) will have the same # of items as original variables ; so we can more easily use with the input shifts: change made 05/13/01 xf = intarr(2,ntot) + BADVAL yf = intarr(2,ntot) + BADVAL ;xf = intarr(2,ng) ;yf = intarr(2,ng) ; exposure map (or array) if keyword_set(exptimes0) then begin message, '* exposure times passed as 1-d array *', /info if n_elements(exptimes) ne sz(3) then begin message, '# of exposure times not the same as # of images!' stop endif if n_elements(expmap) ne 0 then begin message, 'cannot pass both expmap and exptimes!', /info stop endif exptimes = exptimes0(wg) ;expmap = fltarr(sz(1), sz(2), sz(3)) ;for i = 0, sz(3)-1 do $ ; expmap(*, *, i) = float(exptimes(i)) ;stop endif else $ exptimes = fltarr(sz(3))+1.0 ;endif else if not(keyword_set(expmap)) then begin ; if not(keyword_set(silent)) then $ ; print,'* exposure map passed *' ; expmap = fltarr(sz(1),sz(2), sz(3)) + 1.0 ;endif ;; weight map ;if keyword_set(wtmap) eq 0 then begin ; if not(keyword_set(silent)) then $ ; print,'* no weight map passed *' ; wtmap = fltarr(sz(1),sz(2), sz(3)) + 1 ;endif ; find extreme shifts xmax = max(xobj) & xmin = min(xobj) ymax = max(yobj) & ymin = min(yobj) ; ;-------------------------------------------------- ; (1) make mosaic (with trimming) ;-------------------------------------------------- ; if keyword_set(trim) eq 1 then begin xs = sz(1) - round(xmax - xmin) ys = sz(2) - round(ymax - ymin) ;xs = round(sz(1) - (xmax - xmin)) ;ys = round(sz(2) - (ymax - ymin)) ; sanity check if (xs gt MAXSIZE) or (ys gt MAXSIZE) then begin message, 'output array gtrater than MAXSIZE='+strc(MAXSIZE), /info print, 'output array dimensions = ', printcoo(xs, ys) stop endif ; initialize out = fltarr(xs,ys) outexp = fltarr(xs,ys) pieces = fltarr(xs,ys,ng) for i=0,ng-1 do begin ii = wg(i) ;if (goodlist(i) ne 0) then begin if not(keyword_set(silent)) then $ print,format='($,A," ")',strc(i) xf(*,ii) = [0,xs-1] yf(*,ii) = [0,ys-1] x0 = round(xobj(i) - xmin) x1 = xs + x0 - 1 y0 = round(yobj(i) - ymin) y1 = ys + y0 - 1 tmp = imgs0(x0:x1, y0:y1, ii) imbad = tmp ne BADVAL if keyword_set(masks0) then begin bp = badpix(x0:x1, y0:y1) * (masks0(x0:x1, y0:y1, ii) ne 0) * imbad ; bp = badpix(x0:x1, y0:y1) * (masks0(x0:x1, y0:y1, i) ne 0) $ ; * (imgs(x0:x1, y0:y1, i) ne BADVAL) endif else begin bp = badpix(x0:x1, y0:y1) * imbad ; bp = badpix(x0:x1, y0:y1) * (imgs(x0:x1, y0:y1, i) ne BADVAL) endelse wbad = where(bp eq 0, nbad) if (nbad gt 0) then tmp(wbad) = 0.0 out = out + temporary(tmp)*wtmap(x0:x1, y0:y1, ii) outexp = temporary(outexp) + $ fix(bp ne 0) * wtmap(x0:x1, y0:y1, ii) * expmap(x0:x1, y0:y1, ii) ;endif endfor if not(keyword_set(silent)) then $ print ; -------------------------------------------------- ; (2) make mosaic (no trimming) ; (assumes all frames are equal exposure time) ; -------------------------------------------------- endif else begin ; ; (2a) determine final mosaic size ; xs = round(sz(1) + xmax - xmin) ys = round(sz(2) + ymax - ymin) ; sanity check if (xs gt MAXSIZE) or (ys gt MAXSIZE) then begin message, 'output array gtrater than MAXSIZE='+strc(MAXSIZE), /info print, 'output array dimensions = ', printcoo(xs, ys) stop endif ; initialize out = fltarr(xs,ys) outexp = fltarr(xs,ys) if not(keyword_set(silent)) then $ message,' output image will be '+strc(xs)+' x '+strc(ys), /info ; ; (2b) stack images into the output image in the proper places ; if not(keyword_set(silent)) then $ print, form = '($, " stacking images into mosaic: ")' x0 = round(xmax-xobj) x1 = x0 + sz(1) - 1 y0 = round(ymax-yobj) y1 = y0 + sz(2) - 1 xf(*, wg) = transpose([[x0], [x1]]) yf(*, wg) = transpose([[y0], [y1]]) ;xf = transpose([[x0], [x1]]) ;yf = transpose([[y0], [y1]]) for i=0,ng-1 do begin ii = wg(i) ;if (goodlist(i) ne 0) then begin if not(keyword_set(silent)) then $ print,format='($,A," ")',strc(i) ;x0 = round(xmax - xobj(i)) ;y0 = round(ymax - yobj(i)) ;x1 = x0 + sz(1) - 1 ;y1 = y0 + sz(2) - 1 ;xf(*,i) = [x0,x1] ;yf(*,i) = [y0,y1] tmp = imgs0(*, *, ii) imbad = tmp ne BADVAL if keyword_set(masks0) then begin bp = badpix * (masks0(*, *, ii) ne 0) * imbad ;bp = badpix * (masks0(*, *, i) ne 0) * (imgs(*, *, i) ne BADVAL) $ endif else $ bp = badpix * imbad ;bp = badpix * (imgs(*, *, i) ne BADVAL) ; record median of the image mm(i) = median(tmp(where(bp ne 0))) ; add to mosaic ; this is pretty kludgy, but it does save memory! ; tested against last version when no wtmap nor expmap -> ok if keyword_set(wtmap) then begin out(x0(i):x1(i), y0(i):y1(i)) = $ out(x0(i):x1(i), y0(i):y1(i)) + temporary(tmp)*bp*wtmap(*, *, ii) if keyword_set(expmap) then begin outexp(x0(i):x1(i),y0(i):y1(i)) = $ outexp(x0(i):x1(i),y0(i):y1(i)) $ + fix(bp ne 0) * wtmap(*, *, ii) * expmap(*, *, ii) endif else begin outexp(x0(i):x1(i),y0(i):y1(i)) = $ outexp(x0(i):x1(i),y0(i):y1(i)) $ + fix(bp ne 0) * wtmap(*, *, ii) * exptimes(ii) endelse endif else begin out(x0(i):x1(i), y0(i):y1(i)) = $ out(x0(i):x1(i), y0(i):y1(i)) + temporary(tmp)*bp if keyword_set(expmap) then begin outexp(x0(i):x1(i),y0(i):y1(i)) = $ outexp(x0(i):x1(i),y0(i):y1(i)) + fix(bp ne 0) * expmap(*, *, ii) endif else begin outexp(x0(i):x1(i),y0(i):y1(i)) = $ outexp(x0(i):x1(i),y0(i):y1(i)) + fix(bp ne 0) * exptimes(ii) endelse endelse ;out(x0:x1, y0:y1) = out(x0:x1, y0:y1) + tmp*bp ;outexp(x0:x1,y0:y1) = outexp(x0:x1,y0:y1) $ ; + fix(bp ne 0) ;endif endfor if not(keyword_set(silent)) then $ print ; ; (2c) if desired, set the relative sky level to avoid edges ; be sure not to include bad pixels when calculating offsets ; *** this does not take 'wtmap' and 'expmap' into account so far *** ; sky = fltarr(ng) if keyword_set(setsky) then begin if not(keyword_set(silent)) then begin message,'setting optimum sky levels', /info print, ' iterations stop once d(sky) <= '+strc(MINSKY) endif ; ; for each image, see which frames it overlaps with ; and calculate the sky offsets ; over = fltarr(ng,ng) ; which pairs of images overlap ; if do overlap, then has # of overlapping pixels overavg = fltarr(ng,ng) ; overavg(x,y) is the average pixel value ; of image x in the overlapping region ; between x and y if not(keyword_set(silent)) then $ print,format='($," determining overlap regions: ")' for i=0,ng-1 do begin ii = wg(i) if not(keyword_set(silent)) then $ print,format='($,A," ")',strc(wg(i)) bpi = 0 imbad = imgs0(*, *, ii) if keyword_set(masks0) then begin bpi = badpix * (masks0(*, *, ii) ne 0) * imbad ; bpi = badpix * (masks0(*, *, ii) ne 0)$ ; * (imgs(*, *, ii) ne BADVAL) endif else begin bpi = badpix * imbad ; bpi = badpix * (imgs(*, *, ii) ne BADVAL) endelse for j=i,ng-1 do begin ; for j=0,ng-1 do begin jj = wg(j) if (ii ne jj) then begin ; find number of overlapping pixels for each ; pair of images, avoiding bad pixels; ; 'overlap' is a logical map of the good (1) ; and bad (0) pixels for a given overlap region xmx = max([xobj(i), xobj(j)]) xmn = min([xobj(i), xobj(j)]) ymx = max([yobj(i), yobj(j)]) ymn = min([yobj(i), yobj(j)]) ; find size of the overlap region ; crash if the images don't overlap oxs = sz(1) - round(xmx - xmn) oys = sz(2) - round(ymx - ymn) x0i = round(xobj(i) - xmn) x1i = oxs + x0i - 1 y0i = round(yobj(i) - ymn) y1i = oys + y0i - 1 x0j = round(xobj(j) - xmn) x1j = oxs + x0j - 1 y0j = round(yobj(j) - ymn) y1j = oys + y0j - 1 bpj = 0 overlap = 0 imbad = imgs0(*, *, jj) ne BADVAL if keyword_set(masks0) then begin bpj = badpix * (masks0(*, *, jj) ne 0) * imbad ;bpj = badpix * (masks0(*, *, jj) ne 0) $ ; * (imgs(*, *, jj) ne BADVAL) endif else begin bpj = badpix * imbad ;bpj = badpix * (imgs(*, *, jj) ne BADVAL) endelse ; determine # of overlapping pixels ; be sure to check if there is no overlap ; (a bit of a kludge) if (x0i lt 0) or (x0i ge sz(1)) or $ (x1i lt 0) or (x1i ge sz(1)) or $ (y0i lt 0) or (y0i ge sz(2)) or $ (y1i lt 0) or (y1i ge sz(2)) or $ (y0i lt 0) or (y1i ge sz(2)) then begin message, 'image pair = ('+strc(i)+','+strc(j)+ $ ') do not overlap', /info overlap = 0 endif else $ overlap = bpi(x0i:x1i, y0i:y1i) * bpj(x0j:x1j, y0j:y1j) ww = where(overlap ne 0, nover) over(i,j) = nover over(j,i) = nover ; and if they do overlap, find overavg(i,j), ; multiplication by 'overlap' insures we use ; only good pixels if (nover gt 0) then begin overavg(i, j) = $ median((imgs0(x0i:x1i, y0i:y1i, ii))(ww)) overavg(j, i) = $ median((imgs0(x0j:x1j, y0j:y1j, jj))(ww)) ; tmp = pieces(*, *, i) * overlap ; w = where(overlap ne 0, nover) ; if (nover gt 0) then overavg(i, j) = median(tmp(w)) ; overavg(i,j) = total(pieces(*,*,i)*overlap) / over(i,j) ; print,format='($,A," ")',strc(round(over(i,j))) endif endif endfor endfor if not(keyword_set(silent)) then $ print ; ; now calculate iteratively the sky offsets for each image, weighting ; by # of pixels in the overlapping region for each pair of images ; if not(keyword_set(silent)) then begin print,' calculating sky offsets:' print,format='($," (skymax,niter)= ")' ; print,format='($," calculating sky offsets, skymax = ")' endif dsky = fltarr(ng) skymax = 1.0 niter = 0 ; number of iterations skym0 = 1e4 while (skymax gt MINSKY) and (niter lt 2*ng*ng) do begin niter = niter + 1 for i=0,ng-1 do begin skysum = 0.0 nsum = 0.0 for j=0,ng-1 do begin if over(i,j) ne 0 then begin skysum = skysum + $ 0.5*(overavg(j,i)-overavg(i,j)) * over(i,j) nsum = nsum + over(i,j) endif endfor dsky(i) = skysum / nsum sky(i) = sky(i) + dsky(i) endfor for i=0,ng-1 do begin for j=0,ng-1 do begin overavg(i,j) = overavg(i,j) + dsky(i) endfor endfor skymax = max(abs(dsky)) if (skymax lt skym0/10.) then begin skym0 = skymax if not(keyword_set(silent)) then $ print, format = '($,"(",A,",",A,") ")', $ strc(skym0), strc(niter) endif endwhile if not(keyword_set(silent)) then $ print ; normalize sky offsets so that total shift = 0 sky = sky - total(sky)/float(ng) ; add sky offsets to each image and print it if not(keyword_set(silent)) then $ print,format='($," adding sky offsets: ")' for i=0,ng-1 do begin ii = wg(i) imbad = imgs0(*, *, ii) ne BADVAL if keyword_set(masks0) then begin bp = badpix * (masks0(*, *, ii) ne 0) * imbad endif else begin bp = badpix * imbad endelse out(xf(0, ii):xf(1, ii), yf(0, ii):yf(1, ii)) = $ out(xf(0, ii):xf(1, ii), yf(0, ii):yf(1, ii)) + $ sky(i)*bp if not(keyword_set(silent)) then $ print,format='($,A,"(",A6,") ")',strc(ii),strc(sky(i)) endfor if not(keyword_set(silent)) then begin print if not(keyword_set(silent)) then $ print,' required ',strc(niter),' iterations' endif endif ; end of /setsky endelse ; ; (2d) or if desired, set the median of each image to zero ; if keyword_set(submed) then begin if not(keyword_set(silent)) then begin print, '-- setting medians of images to zero --' print, format = '($," medians: ")' endif sky = mm for i = 0, ng-1 do begin ii = wg(i) bpi = 0 imbad = imgs0(*, *, ii) ne BADVAL if keyword_set(masks0) then $ bpi = badpix * (masks0(*, *, ii) ne 0) * imbad $ else $ bpi = badpix * imbad ii = wg(i) out(xf(0, ii):xf(1, ii), yf(0, ii):yf(1, ii)) = $ out(xf(0, ii):xf(1, ii), yf(0, ii):yf(1, ii)) - $ sky(i)*bpi if not(keyword_set(silent)) then $ print, format = '($,A,"(",A6,") ")', strc(ii), strc(sky(i)) endfor endif ; ;---------------------------------------------------------------------- ; (3) now add the images with proper spatial & sky offsets, and ; applying a filter (tossing highest & lowest pixels) if desired ;---------------------------------------------------------------------- ; if not(keyword_set(silent)) then $ message, 'combining images', /info if not(keyword_set(nrej) or keyword_set(median)) then begin ; for i=0,n-1 do begin ; out = out + pieces(*,*,i) * float(pieces(*,*,i) gt BADVAL) ; endfor ; divide by exposure map w = where(outexp ne 0.0, nexp) if (nexp gt 0) then out(w) = out(w) / (outexp(w)) ; set pixels w/o any good images to BADVAL w = where(outexp eq 0.0,nw) if nw gt 0 then out(w) = BADVAL endif else begin ;---------------------------------------------------------------------- ;---------------------------------------------------------------------- message, 'this portion not modified yet', /info ; stop ;---------------------------------------------------------------------- ;---------------------------------------------------------------------- if keyword_set(median) then $ if not(keyword_set(silent)) then $ print,' applying median filter' $ else $ if not(keyword_set(silent)) then $ print,' applying rejection filter, nrej=',strc(nrej) for j=0,ys-1 do begin for i=0,xs-1 do begin cut = pieces(i,j,*) ww = where(cut gt BADVAL,npix) ; apply rejection or median filtering if (npix gt 2*nrej) then begin pix = cut(ww) if keyword_set(median) then begin out(i,j) = median(pix) endif else begin ss = sort(pix) pix = pix(ss) out(i,j) = total(pix(nrej:npix-nrej-1)) / (npix-2*nrej) endelse outexp(i,j) = (npix-2*nrej) ; if not enough images for filtering, just average endif else if (npix gt 0) then begin out(i,j) = total(cut(ww)) / npix outexp(i,j) = npix ; fill empty pixel with BADVAL endif else begin out(i,j) = BADVAL outexp(i,j) = 0 endelse endfor endfor endelse end