; ; defringeflat.pro - Removes fringes from an image ; ; Copyright (C) 2005,2006 Patricio Rojo ; ; This program is free software; you can redistribute it and/or ; modify it under the terms of the GNU General Public License ; as published by the Free Software Foundation; either version 2 ; of the License, or (at your option) any later version. ; ; This program is distributed in the hope that it will be useful, ; but WITHOUT ANY WARRANTY; without even the implied warranty of ; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ; GNU General Public License for more details. ; ; You should have received a copy of the GNU General Public License ; along with this program; if not, write to the Free Software ; Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, ; MA 02110-1301, USA. ; ;+ ; NAME: ; DEFRINGEFLAT ; ; PURPOSE: ; This function removes the fringes from a flat field. ; ; CATEGORY: ; Pato's Array Data Reduction Potpourri. ; ; CALLING SEQUENCE: ; ; defringedflat = DEFRINGEFLAT(Flat, Refper) ; or ; defringeflat_version = DEFRINGEFLAT(/QUERY) ; ; ; INPUTS: ; Flat: Two dimensional array containing the flat field to ; defringe. ; Refper: Approximate fringe period (in pixels) at position ; COLREF, ROWREF. ; ; KEYWORD PARAMETERS: ; GENERAL KEYWORDS ; QUERY: If specified it ignores any other parameters or ; arguments and just return a string with the current ; release number of the defringeflat software. ; ; ARRAY KEYWORDS ; BINW: Number of rows included in median average to improve ; signal-to-noise ratio. Default is 10. ; COLREF: Reference image column for which REFPER value is ; given. Because the default is the middle column, this ; parameter should always be set to something different ; in the common case of a discontinuity in this ; column. ; ROWREF: Reference row for which REFPER value is ; given. Default is middle row. ; TRIM: [Deprecated]. Use ITRIM instead. Keyword is going to ; be disabled in a future releas. ; ITRIM: Ignore datapoints outside these boundaries. Format is ; [x0, y0, x1, y1]. ; GOODPIX: Mask of good pixels. Must have same shape as Flat ; array. 1 is good pixel, 0 is bad. ; GPIXRAD: Radius of median filter to apply whenever a good ; pixel mask is supplied. The median value will replace ; the bad value. Default is 5. ; ; TRACE FITTING KEYWORDS ; NOCOI: Ignore cone-of-influence (COI) boundary when fitting ; the trace. If this is not set, then trace fitting ; will stop at the column where the trace passes ; through the COI boundary. ; SILENT: Set this for a silent run without any kind of ; progress output. Error messages will show, though. ; FITNAME: What kind of fit to do. Currently 'GFC', ; 'GFCest', 'GFCnb', 'GVC', and 'true' are ; supported. Default is 'GFC'. See below for an ; explanation of each method under TRACE FITTING METHOD. ; OSAMP: Trace-fit oversampling factor. Default is 10. ; BLINDCONTINUE: Continue the specified number of column blindly ; if a peak was not successfully found. It uses the ; last known trace center to find the next one. If ; none of this extra columns gives a successfull fit ; then it stop looking in that direction. Default is 25 ; pixels. ; KEEPBACK: If set, do not subtract the background from the trace, ; i.e. do not consider the wavelet's background as part ; of the flat but as part of the fringe pattern ; instead. This option also put zeroes outside the ; local minima to avoid contamination from wavelets ; from different amplitudes. This switch has no ; effect on the methods without background. ; TRACELIM: Specify this to have a limit on the trace ; value. It should be a two-element array. The center ; of the trace won't go below TRACELIM[0] or above ; TRACELIM[1]. If the trace tries going beyond these ; values, then a backsearch is started using the ; reference period as a reference. Note that this ; keyword has not been tested thoroughly, please report ; any bug that you might find. ; PERUPDATE: Set this to update the value of REFPER with the ; fitted period in the initial position. Useful when ; processing several flats in a loop. ; SHOWTRACE: Show the fitted trace after having calculated it for ; each row. ; ; WAVELET PARAMETERS KEYWORDS ; WAVEMOTHER: Name of wavelet to use. Currently only Morlet, Paul, ; or DOG, with WAVEPARAM 6, 4, or 2, respectively, are ; supported. ; WAVEPARAM: Parameter of the wavelet, currently only some values ; are supported, see WAVEMOTHER. ; NOPAD: Set this flag not to zero-pad the data before doing ; the wavelet transform. ; ; FRINGE PARAMETERS KEYWORDS ; SMOOTHFRINGE: Set this to patch and smooth (see below) the ; parameters of the trace in the wavelet space before ; reconstructing. ; NSIGMA: Rejection threshold for fitted fringe parameters. If ; any parameter is more than NSIGMA standard deviations ; different from its median value then patch it. It has ; an effect only when SMOOTHFRINGE is ; specified. Default is 1.5. ; RADMEDFILTER: Radius of the median filter to apply when ; discriminating bad parameter values. It has an effect ; only when SMOOTHFRINGE is specified. Default is 4. ; RADSMOOTH: Radius of mean filter to apply when smoothing in the ; parameter space. It has an effect only when ; SMOOTHFRINGE is specified. Default is 9. ; ; OPTIONAL OUTPUT KEYWORD ; OTRIM: Output trim. Same format as ITRIM. This area is given ; as a maximum possible area, it does not include the ; limits set by Cone of Influence (COI). ; FRINGE: Returns the two dimensional fringe that was ; constructed from the wavelet fit. ; PRMFIT: Returns the parametric fit used internally to ; reconstruct the fit. The order of each parameter is ; dependent on the value of FITNAME keyword, it has no ; effect when 'true' is selected. This keyword is ; mostly intended for debugging. ; ; OUTPUTS: ; This function returns the corrected flat (input Flat minus the ; fitted fringe). ; ; SIDE EFFECTS: ; Value of input parameter Refper is overwritten when using ; PERUPDATE option. ; ; RESTRICTIONS: ; Do no operate heavy machinery while using this:D. ; ; PROCEDURE: ; Area outside TRIM is cropped. ; A 1D median average of BINW rows is computed for each of the ; rows. Let's call each of the rows obtained in this manner ; an enhanced row. ; For each of those meanrows, the slope is subtracted ; through a linear fit. Then, a wavelet transform is performed ; according to the WAVEMOTHER and the WAVEPARAM keywords. ; Starting from Refper at position COLREF and ROWREF, the trace ; of prominent periods in the absolute value of the wavelet space ; is followed along the columns. The trace profile is fit by a ; FITNAME function. Repeating this procedure for each row, we ; obtain a fitted set of parameters (e.g., center, width, ; amplitude of a Gaussian) for each pixel in the array. ; If SMOOTHFRINGE is set, then each of the parameters is ; patched (compared against the median value using NSIGMA and ; RADMEDFILTER) and painted (smoothed according to RADSMOOTH). ; Then, the wavelet transform is reconstructed from the ; (smoothed) parameters. Finally, an inverse wavelet transform is ; performed to obtain the fringe pattern, which is then removed ; from the flat. ; ; TRACE FITTING METHODS: ; All of the following procedures are applied to the data points ; within the local minima. Those minimas are found at both sides ; of the profile which contains the extrapolated reference period. ; The different alternatives for the FITNAME keyword are: ; * "GFC": (default). After oversampling, a center is found by ; the highest point in the local profile. Then a ; Gaussian function is fit to the trace's cross ; section but with the center previously fixed. Both ; centers estimates will be exactly the same only ; when the profile is symmetric. Four parameters ; (height, center, width, background) are obtained ; and can be smoothed. A negative background makes no ; sense. Hence, if such is the fit at a particular ; trace position, then the fit is retried on that ; pixel using 'GFCnb'. Using this method will not ; fit every wavelet point perfectly, but is more ; robust against removing non-fringe patterns from ; the image. ; * "GFCnb": Another modification to 'GFC', where there is ; no background parameter (height, center, and width ; only). For real fringes, this and 'GFC' option ; give similar results but for sine-wave generated ; fringes, this option seems to work better. ; * "GVC": Modification of 'GFC'. The center parameter is ; now fitted through a Gaussian and not fixed by the ; highest value. It fallbacks to 'GVCnb' when ; the background is negative. ; * "GVCnb": Modification of 'GVC', where there is no ; background parameter. ; * "GFCest": Yet another modification of 'GFC' fit, but the ; parameters are only estimated. 'Center' is the ; coordinate where the trace peaks, 'height' is the ; trace's value at center, 'width' is one-fifth ; (chosen arbitrarily as a good approximation) of the ; separation between the local lower minima, and ; 'background' is zero. There is not really any ; advantage in choosing this alternative instead of ; the above other than for speed or as a failproof ; alternative. ; * "true": Once the peak of trace's cross section has been ; found, all the points between the local lower ; minima at the sides of the trace are used to ; reconstruct the fringe's wavelet. The main ; disadvantage of this approach is that every single ; feature of the wavelet around the peak period is ; going to be removed from the original image. This ; might include real features from the image that ; were not intended to be cleaned. Smoothing of ; parameters is not available. ; ; DEPENDS ON: ; maxtrace.pro Finds a trace in a 2 dimensional array. ; waverecon.pro Produces an inverse wavelet transform. ; wavelet.pro Do a forward or inverse wavelet ; transform. ; Routine by Torrence and Compo. It is ; available from the DEFRINGEFLAT homepage ; http://www.das.uchile.cl/~pato/sw/ ; The latest version should be available directly ; from the authors at ; http://paos.colorado.edu/research/wavelets/ ; ; ; EXAMPLE: ; Read example flat and defringe it. Start with a reference period ; of 40 pixels at the center of the array, some trimming. Use a ; Gaussian fit, interpolate period in the wavelet transform to 10 ; times more points for improved trace fitting. Show the fitted ; trace after finding it in every row. Patch and smooth each of ; the parameters of the Gaussian fit (center, width, height, and ; background). ; ; IDL> flat = readfits('example/flat.fits', h) ; IDL> func = 'GFC' ; IDL> clean = defringeflat(flat, 40, binw=40, /show, $ ; trim=[10,145,950,950], osamp=10, $ ; /smooth, fitname=func, fringe=fringe) ; IDL> tvscl, flat ; IDL> tvscl, fringe ; IDL> tvscl, clean ; ; MODIFICATION HISTORY: ; Written by: Pato Rojo, Cornell. Jul 2005 - Mar 2006 ; pato@astro.cornell.edu ;- FUNCTION findandpatch_defringe, data, radius, nsigma, noval COMPILE_OPT idl2, hidden ;look for outliers and patch med = MEDIAN(data, radius*2, /DOUBLE) n = SIZE(data, /DIM) - radius - 1 difprm = data - med outliers = WHERE(ABS(difprm) GT $ STDDEV(MEDIAN(difprm[radius:n[0],radius:n[1]], 2))*nsigma $ OR data EQ noval) IF outliers[0] EQ -1 THEN RETURN, data ret = data ret[outliers] = med[outliers] RETURN, ret END ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO pnp_defringe, prmfit, radsmooth, radmedfilter, $ nsigma, silent, noval COMPILE_OPT idl2, hidden ;Patch and paint(smooth) each of the fitted parameters IF ~ silent THEN print, ' Correcting parameters', FORMAT='(a,$)' FOR e = 0, (size(prmfit, /dim))[2]-1 DO BEGIN IF ~ silent THEN print, string(e, FORMAT='(i2)') + ".", FORMAT='(a,$)' ;Patch. Find and fix outliers smoothed = FINDANDPATCH_DEFRINGE(REFORM(prmfit[*, *, e]), $ radmedfilter, nsigma, noval) ;Paint. Smooth parameters IF ~ silent THEN print,".", FORMAT='(a,$)' ind = WHERE(smoothed EQ noval) IF (ind[0] NE -1) THEN smoothed[ind] = !VALUES.F_NAN smoothed = SMOOTH(smoothed, 2*radsmooth+1, /NAN, /EDGE) x = WHERE(~ FINITE(smoothed)) IF x[0] NE -1 THEN smoothed[x] = noval IF ~ silent THEN print, ".", FORMAT='(a,$)' prmfit[*, *, e] = smoothed ENDFOR IF ~ silent THEN print, " " END ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FUNCTION maketrace_defringe, prm, lm, idx, period, fitname, noval COMPILE_OPT idl2, hidden ;Rebuild trace nwa = DBLARR(N_ELEMENTS(idx), N_ELEMENTS(period)) cnt = prm[*,1] nx = (SIZE(prm, /DIM))[0] nprm = (SIZE(prm, /DIM))[1] FOR xi = 0, nx-1 DO IF (cnt[xi] NE noval) THEN BREAK SWITCH fitname OF 'GFCnb': 'GVCnb': BEGIN FOR i = xi, nx-1 DO BEGIN IF (cnt[i] EQ noval) THEN CONTINUE l1 = lm[0,i] l2 = lm[1,i] gauss = prm[i,0] * EXP(-0.5*((period-prm[i,1])/prm[i,2])^2) nwa[i,l1:l2] = gauss[l1:l2] ENDFOR BREAK END 'GFC': 'GVC': 'GFCest': BEGIN FOR i = xi, nx-1 DO BEGIN IF (cnt[i] EQ noval) THEN CONTINUE l1 = lm[0,i] l2 = lm[1,i] gauss = prm[i,0] * EXP(-0.5*((period-prm[i,1])/prm[i,2])^2) + $ prm[i,3] nwa[i,l1:l2] = gauss[l1:l2] ENDFOR BREAK END ELSE: message, "Method '" + fitname + $ "' not enabled in maketrace_defringe" ENDSWITCH RETURN, nwa END ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO buildfringe_defringe, silent, prmfit, locmin, nwp, noval, $ ;border cdelta, psi0, period, scale, fitname, $ ;parameters fringe ;returns COMPILE_OPT idl2, hidden ;rebuild fringe from fitted wavelet nx = (SIZE(prmfit, /DIM))[0] ny = (SIZE(prmfit, /DIM))[1] idx = INDGEN(nx) IF size(flat, /type) EQ 5 THEN fringe = DBLARR(nx, ny) $ ELSE fringe = fltarr(nx, ny) IF keyword_set(locmin) THEN dozero = 1 ELSE dozero = 0 ;reconstruct IF ~ silent THEN print," Reconstructing..." FOR j = 0, ny-1 DO BEGIN ;rebuild trace IF dozero THEN lm = REFORM(locmin[*,*,j]) $ ELSE lm = [0, N_ELEMENTS(period)-1]#REPLICATE(1.0,nx) nwa = MAKETRACE_DEFRINGE(REFORM(prmfit[*,j,*]), lm, $ idx, period, fitname, noval) cnw = COMPLEX(nwa*COS(nwp[*,*,j]), nwa*SIN(nwp[*,*,j])) ;rebuild row fringe[*,j] = WAVERECON(cnw, 1, scale, cdelta, psi0) ENDFOR END ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO tracefit_defringe, flat, binw, osamp, $ ;basics rrefper, colref, rowref, $ ;reference blindcontinue, doshowtrace, keepback, fitname, $ tracelim, silent, nocoi, $ ;flags waveparam, wavemother, pad, noval, $ ;wavelet parameters prmfit, nwp, period, scale, locmin ;returns COMPILE_OPT idl2, hidden ;find fringe from fitted wavelet of data IF ~ silent THEN print,' Fitting trace for each row...' nx = (SIZE(flat, /DIM))[0] ny = (SIZE(flat, /DIM))[1] cbin = binw/2 idx = indgen(nx) w = WAVELET(idx, 1, COI=coi, SCALE=scale, PERIOD=period, $ MOTHER=wavemother, PARAM=waveparam, PAD=pad) nper = N_ELEMENTS(scale) nwp = DBLARR(nx, nper, ny) IF (STRPOS(fitname,'GFC') NE -1 OR STRPOS(fitname,'GVC') NE -1) THEN BEGIN IF ~ silent THEN print, fitname, FORMAT='(a,$)' prmfit = DBLARR(nx, ny, 4) ENDIF ELSE $ MESSAGE, "Fitting method '" + fitname + "' not enabled in tracefit" IF keepback THEN BEGIN locmin = DBLARR(2, nx, ny) PRINT, 'kb', FORMAT='(a,$)' ENDIF PRINT, ':', FORMAT='(a,$)' degree = 2 finrow = ny - binw - 1 steprow = 1 FOR order = 0, 1 DO BEGIN refper = rrefper FOR j = rowref - cbin, finrow, steprow DO BEGIN ;median average, subtract slope, and nullify borders p = MEDIAN(REFORM(flat[*, j:j+binw]), $ DIMENSION=2, /DOUBLE) a = POLY_FIT(idx, p, degree, YFIT=yfit) p = p - yfit ;produce wavelet w = WAVELET(p, 1, COI=coi, SCALE=scale, PERIOD=period, $ MOTHER=wavemother, PARAM=waveparam, PAD=pad) IF ~ silent THEN BEGIN IF steprow LT 0 THEN print, '<', FORMAT='(a,$)' $ ELSE print,'>', FORMAT='(a,$)' ENDIF IF nocoi THEN coi2max = 0 ELSE coi2max = coi ;find trace prm = MAXTRACE(abs(w), period, colref, refper, $ OSAMP=osamp, FIT=fitname, /EXPY, TRACELIM=tracelim, $ NOVAL=noval, COIBORDERS=coi2max, PEAKS=peaks, $ BLINDCONTINUE=blindcontinue, LOCALMIN=localmin) prmfit[*, j+cbin, *] = prm IF keepback THEN locmin[*, *, j+cbin] = localmin ;show if requested IF doshowtrace AND refper NE 0 THEN BEGIN v = ABS(ALOG(refper)/ALOG(10)-FIX(ALOG(refper)/ALOG(10))) ;Use different yrange if it is very close to 10^i if (v GT 0.04 AND v LT 0.95) then $ yrng = 10^(fix(alog(refper)/alog(10)+[1,0])) $ else $ yrng = 10^(fix(alog(refper)/alog(10)+[1,-1]+0.5)) !P.POSITION = [.1,.83,.95,.95] blank = REPLICATE(" ",30) PLOT, p, /YSTYLE, /XSTYLE, XTICKNAME=blank, $ TITLE=STRING('Row: trimx + bin/2 + ', j, FORMAT='(a,i0)'), $ YTICKS=3 !P.POSITION = [.1,.18,.95,.83] xr = !X.CRANGE ttrace = WHERE(prmfit[*,j+cbin,1] NE 0) rpidx = (WHERE(period GT refper))[0] mval = MAX(ABS(w[colref-5:colref+5, (rpidx-5)>0:(rpidx+5)