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,$ subpixel=subpixel, INTERP_TYPE=INTERP_TYPE,$ pieces=pieces,xsize=xsize,ysize=ysize,extremes=extremes ;+ ; PURPOSE: ; 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 ; subpixel positions will be handled via ; interpolated shifts. (MDP Aug 2001) ; ; 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 (like badpix, but a 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] ; /subpixel Perform subpixel shifts on the images before ; registering. Requires Nyquist sampled images or better! ; ; CAVEATS: ; - 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 ; 08/15/2001 (MDP): Added sub-pixel interpolation with keyword SUBPIXEL ; 06/25/2002 (MDP): Re-added the pieces variable, for measurement of ; noise at each pixel for polarimetry. ; 08/06/2002 (MDP): Made median keyword work again. ; 2004-3-16 (MDP): Fixed bug with exptimes vs exptimes0 ; 2004-4-12 (MDP): Fixed bug introduced in previous bugfix; made /setsky work ; with exposure times. Fixed bug in subpixel code when ; some images are marked not good. ; ; -> modify /setsky so will take information in the weight map (wtmap) ;- ;########################################################################### ; ; LICENSE ; ; This software is OSI Certified Open Source Software. ; OSI Certified is a certification mark of the Open Source Initiative. ; ; Copyright © 1995-2003 by Michael Liu & Marshall Perrin ; ; This software is provided "as-is", without any express or ; implied warranty. In no event will the authors be held liable ; for any damages arising from the use of this software. ; ; Permission is granted to anyone to use this software for any ; purpose, including commercial applications, and to alter it and ; redistribute it freely, subject to the following restrictions: ; ; 1. The origin of this software must not be misrepresented; you must ; not claim you wrote the original software. If you use this software ; in a product, an acknowledgment in the product documentation ; would be appreciated, but is not required. ; ; 2. Altered source versions must be plainly marked as such, and must ; not be misrepresented as being the original software. ; ; 3. This notice may not be removed or altered from any source distribution. ; ; For more information on Open Source Software, visit the Open Source ; web site: http://www.opensource.org. ; ;########################################################################### ; empty spaces of the mosiac and mosaic pieces ; get filled with BADVAL BADVAL = -1e6 ;if not(keyword_set(INTERP_TYPE)) then INTERP_TYPE="I" ; default is regular IDL interpolation ; Default interp type changed. MDP 2005-6-30 if not(keyword_set(INTERP_TYPE)) then INTERP_TYPE="F" ; default is now FOURIER interpolation if ~(keyword_set(subpixel)) then begin message,/info," Using PIXEL level registration and shifting" endif else begin message,/info,"Using SUBPIXEL level registration and shifting" message,/info," with interpolation method: "+INTERP_TYPE endelse ; 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 = 4000 if arg_present(pieces) or keyword_set(median) then use_pieces=1 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) if sz[0] eq 2 then sz[3]=1 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 ; 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 ; exposure map (or array) if keyword_set(exptimes0) then begin message, '* exposure times passed as 1-d array *', /info if n_elements(exptimes0) 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!' endif exptimes = exptimes0 endif else $ exptimes = fltarr(sz(3))+1.0 ; find extreme shifts if ~(keyword_set(extremes)) then begin xmax = round(max(xobj)) & xmin = round(min(xobj)) ymax = round(max(yobj)) & ymin = round(min(yobj)) endif else begin xmax = extremes[2] & xmin = extremes[0] ymax = extremes[3] & ymin = extremes[1] endelse ;-------------------------------------------------- ; (1) make mosaic (with trimming) if keyword_set(trim) eq 1 then message,"Trimming version no longer supported. Too bad!" ; -------------------------------------------------- ; (2) make mosaic (no trimming) ; (assumes all frames are equal exposure time) ; ; OK, let's change that assumption. The code will ; now assume that all frames have been *scaled* to the ; same intensity scale. The 'exptime' array should ; contain the original exposure times. The 'weights' ; array should contain the weights for combining the ; images. If 'weights' is not set, it will be set ; equal to 'exptime'. If both are not set, both are ; set equal to 1 and none of this matters. ; ; Output will be as follows. ; out mosaic constructed as a weighted average ; of the good pixels in the input images. ; outexp output exposure map, consisting of the ; total integration time at each pixel ; sumweight the total of the weights at each ; pixel, used to normalize the weighted ; average when computing out ; ; NaN pixels are fixed before subpixel shifting. ; -------------------------------------------------- ; ; (2a) determine final mosaic size ; if not(keyword_set(xsize)) then xs = round(sz(1) + xmax - xmin) else xs=xsize if not(keyword_set(ysize)) then ys = round(sz(2) + ymax - ymin) else ys=ysize ; sanity check if (xs gt MAXSIZE) or (ys gt MAXSIZE) then begin message,"Final array bigger than MAXSIZE!" endif ; initialize out = fltarr(xs,ys) outexp = fltarr(xs,ys) if keyword_set(weights) then sumweights = fltarr(xs,ys) ; pieces should have ntot elements, not ngood. ; Pieces which are marked not good will be all NANs ; Changed 2006-02-22 MDP if keyword_set(use_pieces) then pieces = fltarr(xs,ys,ntot)+!values.f_nan 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) x0 = xmax-round(xobj) x1 = x0 + sz(1) - 1 ;y0 = round(ymax-yobj) y0 = ymax-round(yobj) ; Subtle change in definition of y0 and x0 made on 2004-08-28 ; by Marshall Perrin. This is to prevent a bug when yobj is exactly ; something.5000000 which always rounds up, so y0 was getting set ; one integer too high for consistency with ysize. y1 = y0 + sz(2) - 1 if keyword_set(subpixel) then begin fdx=(xmax-xobj) - round(xmax-xobj) ; sub-pixel portions of shifts fdy=(ymax-yobj) - round(ymax-yobj) endif xf[*, wg] = transpose([[x0], [x1]]) yf[*, wg] = 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) tmp = imgs0[*, *, ii] ; 0 is bad, so test for what's good to create the badmask: imbad = (tmp ne BADVAL and finite(tmp)) if keyword_set(masks0) then begin bp = badpix * (masks0[*, *, ii] ne 0) * imbad endif else $ bp = badpix * imbad ; record median of the image mm[i] = median(tmp(where(bp ne 0))) ; clean up any NaNs present in the image - needed to ; prevent NaNs from spreading throughout the totals... if keyword_set(subpixel) then begin wnf = where(finite(tmp,/NAN),wnfcount) if (wnfcount gt 0) then begin message,/info,"Cleaning up NaNs before shifting/mosaicing." fixpix,tmp,0,tmp0,/NaN,/silent,/quick tmp = temporary(tmp0) endif endif ;tmp = tmp*bp ; do this *before* subpixel shifting -MDP ; NO! The above comment is completely wrong I think now (2004-10-15) ; because we want all bad pixels to have been cleaned (e.g. by ; ircal_fixpix) before we do the subpixel shifting. That prevents ; the "badness" from corrupting adjacent pixels. ; perform sub-pixel shifting on the temporary buffer ; there better not be any NANs in this or we're screwed. ; ; At this point, we only do the *fractional part* of the subpixel ; shifts. This will always be < 0.5 pixel either way. ; By doing it this way, we can use the original bad pixel mask ; for the nearest pixel after we do the shifting, and it's ; still kinda sorta correct. ; ; EXCEPT what about the edges of the array? Those may get messed ; up... ; tmp0=tmp if keyword_set(subpixel) then begin if not(keyword_set(INTERP_TYPE)) then begin tmp=fftshift(tmp,fdx[i],fdy[i]) ; my fft shift endif else begin tmp=image_shift(tmp,fdx[i],fdy[i],INTERP_TYPE=INTERP_TYPE) ; starfinder shift. ; modify the bad pixel mask to block out bad edges image_shift_maskedges,bp,fdx[i],fdy[i] ;atv,[tmp,tmp0],/bl ;stop endelse ;bp=image_shift(float(bp),fdx[i],fdy[i],INTERP_TYPE=INTERP_TYPE) ; starfinder shift. ;pixmask_fillinholes,bp ;bp = fix(bp) ;print,fdx[i],fdy[i] endif tmp=tmp*bp ; moved outside above if block so this always applies. bpmask = fltarr(sz[1],sz[2]) wbad = where(bp eq 0,badct) if badct gt 0 then bpmask[wbad]=!values.f_nan ;atv,[[[tmp0]],[[tmp]]],/bl ;stop ; 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 ; use individual weights at each pixel out[x0[i]:x1[i], y0[i]:y1[i]] = $ out[x0[i]:x1[i], y0[i]:y1[i]] + (tmp)*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] if keyword_set(use_pieces) then pieces[x0(i):x1(i),y0(i):y1(i),ii] = tmp/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] if keyword_set(use_pieces) then pieces[x0(i):x1(i),y0(i):y1(i),ii] = tmp/exptimes[ii] endelse endif else begin if (keyword_set(weights)) then begin ; use weights for each exposure message,"Weights doesn't necessarily work yet." out[x0[i]:x1[i], y0[i]:y1[i]] = $ out[x0[i]:x1[i], y0[i]:y1[i]] + temporary(tmp)*weights[ii] b = fix(bp ne 0) ; in this case we normalize by the weights, not by exposure time sumweights[x0[i]:x1[i],y0[i]:y1[i]] = $ sumweights[x0[i]:x1[i],y0[i]:y1[i]] + b * weights[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]] + b * expmap[*, *, ii] endif else begin outexp[x0[i]:x1[i],y0[i]:y1[i]] = $ outexp[x0[i]:x1[i],y0[i]:y1[i]] + b * exptimes[ii] endelse endif else begin ; uniform weighting, just using exposure time or expmap and bad pixel ; mask out[x0[i]:x1[i], y0[i]:y1[i]] +=tmp if keyword_set(expmap) then begin outexp[x0[i]:x1[i],y0[i]:y1[i]] += fix(bp ne 0) * expmap[*, *, ii] if keyword_set(use_pieces) then pieces[x0(i):x1(i),y0(i):y1(i),ii] = tmp/expmap[*,*,ii]+bpmask endif else begin outexp(x0(i):x1(i),y0(i):y1(i)) += fix(bp ne 0) * exptimes[ii] if keyword_set(use_pieces) then pieces[x0(i):x1(i),y0(i):y1(i),ii] = tmp/exptimes[ii]+bpmask endelse 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 if keyword_set(weights) or keyword_set(wtmap) then $ message,"Sky level overlap code doesn't work with weights! Better fix this..." overavg(i, j) = $ median((imgs0(x0i:x1i, y0i:y1i, ii))(ww))/exptimes[ii] overavg(j, i) = $ median((imgs0(x0j:x1j, y0j:y1j, jj))(ww))/exptimes[jj] ; 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) ; check the case where the sky values are bogus. ; This may happen if you try the degenerate case of ; mosaicing only one image. wnan = where(~finite(sky),nanct) if nanct gt 0 then sky[wnan]=0 ; add sky offsets to each image and print it if not(keyword_set(silent)) then begin print,format='($," adding sky offsets: ")' endif ;pieces0 = pieces ;out0 = out 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)] += sky[i]*bp if keyword_set(use_pieces) then begin pieces[xf(0, ii):xf(1, ii), yf(0, ii):yf(1, ii),i] += sky[i]*bp/exptimes[i] endif if not(keyword_set(silent)) then print,format='($,A,"(",A6,") ")',strc(ii),strc(sky(i)) endfor if not(keyword_set(silent)) then begin print 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 and finite(imgs0(*, *, ii))) 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 if keyword_set(weights) then begin ; normalize by the weights message,"Normalizing weighted average.",/info w = where(sumweights ne 0.0, nexp) if (nexp gt 0) then out(w) = out(w) / (sumweights(w)) wbad = where(sumweights eq 0.0,nw) endif else begin message,"Normalizing average by exposure time.",/info ; divide by exposure map w = where(outexp ne 0.0, nexp) if (nexp gt 0) then out(w) = out(w) / (outexp(w)) wbad = where(outexp eq 0.0,nw) endelse ; set pixels w/o any good images to Not-a-Number if nw gt 0 then out(wbad) = !values.F_NAN endif else begin ;---------------------------------------------------------------------- ;---------------------------------------------------------------------- message, 'NREJ and MEDIAN probably dont work; this code not updated yet' 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) = !values.F_NAN outexp(i,j) = 0 endelse endfor endfor endelse end