pro DIcal, inFn, DNFN=dnfn, RADFN=radfn, REVRADFN=revradfn, IFFN=iffn, REVIFFN=reviffn,$ VERBOSE=verbose, SILENT=silent, FORCE=force, MODULES=modules, $ CHOSEN_ONLY=chosen_only, BADPIXS=badpixs, DESPIKE=despike, $ FILLGAPS=fillgaps, DENOISE=denoise, GEOM=geom, MTF=mtf, $ DESMEARALG=desmearalg, DARKFN=darkfn, FLATFN=flatfn, $ MAXGAPSIZE=maxgapsize, USEDARKMODEL=usedarkmodel, $ COMPRESSFN=compressfn, ADCFN=adcfn, BADPIXSFN=badpixsfn, $ GEOMFN=geomfn, MTFFN=mtffn, VISCONSTFN=visconstfn, LINDNFN=lindnfn, $ SPECFN=specfn, CONSTMAPFN=constmapfn, GAINFN=gainfn, SPIKESIG=spikesig,$ SPIKEITER=spikeIter, SPIKEBOX=spikebox, SPIKEMED=spikemed, SPIKEALG=spikealg,$ BADPIXSINTERP=badpixsinterp, MISSINGINTERP=missinginterp, $ HLUTFN=hlutfn, COMPRESSMETH=compressmeth, MTFALG=mtfalg, MTFPARAM=mtfparam $ , debug=debugArg dodebug = keyword_set(debugArg) ? 1b : 0b nodebug = 1b - dodebug ;;----------------------------------------------------------------------------- ;; PURPOSE: ;; Performs the calibration on a Deep Impact FITS image. ;; 1. Temperatures and Voltages calibrated ;; 2. Decompression ;; 3. Flag saturated pixels ;; 4. Correct for uneven ADC bit weighting ;; 5. Linearize DN values (IR only) ;; 6. Remove electrical crosstalk ;; 7. Subtract the dark frame (dark found using interpolation) ;; 8. Normalize quadrant gains ;; 9. Divide by the normalized flat field ;; 10. Flag bad pixels ;; 11. Transfer smear removed (VIS only) ;; 12. Image converted from DN units to radiance and/or I/F units ;; 13. Bad Pixels and Gaps are Reclaimed ;; 14. Cosmic Ray Removal ;; 15. Denoising ;; 16. Deconvolution ;; 17. Geometric calibration ;; ;; By default, modules 2,3,5,6,7,9,10,11,12 and 13 are applied ;; ;; For those modules that are applied by default and have associated ;; keywords, if the keyword is explicitly set to 0, the module is not ;; applied. ;; ;; The output FITS files consists of: ;; Primary Data Unit - Calibrated Image ;; 1st Extension - compressed version of the flags image. To extract, ;; use the di_readfits_flagmap routine. ;; 2nd Extension - Wavelength map (IR only) ;; 3rd Extension - Delta Wavelength map (IR only) ;; 4th Extension - A SNR map giving an estimate of the SNR for that pixel ;; ;; CALLING SEQUENCE: ;; DIcal, inFn, RADFN=radfn, IFFN=iffn, /VERBOSE, /SILENT, /FORCE, $ ;; MODULES=modules, /CHOSEN_ONLY, /BADPIXS, /DESPIKE, /FILLGAPS,$ ;; /DENOISE, /GEOM, /MTF, DESMEARALG=desmearalg, /RAD, /IOVERF, $ ;; DARKFN=darkfn, FLATFN=flatfn, MAXGAPSIZE=maxgapsize, $ ;; /USEDARKMODEL ;; ;; REQUIRED INPUTS: ;; inFn - The filename of the FITS image to calibrate ;; ;; OUTPUTS: ;; ;; OPTIONAL INPUT KEYWORDS: ;; DNFN - Filename of the calibrate file in DN ;; RADFN - Filename of the calibrated file in radience units to output ;; REVRADFN - Filename of the reversibly calibrated file in radiance units ;; to output ;; IFFN - Filename of the calibrated file in I/F units to output. RADFN ;; and/or IFFN must be set ;; REVIFFN - Filename of the reversibly calibrated file in I/F units to ;; output ;; VERBOSE - If applied then status messages will be printed ;; SILENT - If flagged, then error message won't be displayed ;; FORCE - If flagged then an error in a calibration step won't abort ;; the calibration ;; MODULES - A activeModules structure specifying which modules to use. ;; This allows complete control of what modules are applied ;; (Any module whose structure element is flagged is applied) ;; CHOSEN_ONLY - Only apply those modules that are explicity chosen are ;; applied. This can be used to apply badPixs, despike, fillGaps, ;; denoise, geom, mtf, desmear, dark frame and flat field ;; operations only. ;; BADPIXS - Applies the module to interpolate over bad pixels. ;; DESPIKE - Applies the module to despike the image. ;; FILLGAPS - Causes gaps in data whose size is less than MAXGAPSIZE ;; (default=9999) to be interpolated over. ;; DENOISE - Applies the module that tries to remove random gaussian noise ;; GEOM - Applies a geometric rubber sheet correction ;; MTF - Applies a convolution filter to compensate for MTF and scattered ;; light ;; DESMEARALG - Determines which desmearing algorithm to use (VIS only) ;; 0 - No desmearing occurs ;; 1 - Remove smear using the POC values (default for modes 1-6,9) ;; 2 - Remove smear using column averaging (default for modes 7&8) ;; DARKFN - Filename of the dark FITS file to sutract from the image. ;; If this is not set, then a dark frame is found based on the ;; USEDARKMODEL keyword. Setting this keyword is equivilant to ;; explicity choosing the the dark subtraction operation. ;; FLATFN - Filename of the flat FITS file to divide the image by. If ;; this is not set, the default flat field is used. Setting this ;; keyword is equivilant to explicity choosing the flat field ;; compensation module ;; MAXGAPSIZE - The maximum size of a data gap to fill if the Fill gap ;; module is run. ;; USEDARKMODEL - If flagged then the dark frame to subtract is found ;; by using a pre-calculated dark model at every pixel. ;; Explicity setting this keyword to any value is equivilant of ;; choosing the dark subtraction module. ;; COMPRESSFN - filename of the decompression LUT to use ;; ADCFN - filename of the bit weighting correction LUT to use ;; BADPIXSFN - filename of the bad pixel map to use ;; GEOMFN - filename of the geometric tie points to use ;; LINDNFN - filename of the linearization polynomial to use ;; MTFFN - filename of the convolution filter to use ;; SPECFN - filename of the spectral map to use ;; CONSTMAPFN - filename of the ST calibration constant map to use ;; VISCONSTFN - filename of the calibration constant to use ;; GAINFN - filename of the gain maps to use ;; SPIKESIG - threshold for the despiking algorithm in distance from mean ;; in standard deviations ;; SPIKEITER - maximum number of iterations for the despking routine ;; SPIKEBOX - How big a box to use to calculate statistics for despiking ;; SPIKEMED - Should despiking use Median? Otherwise it will use mean ;; SPIKEALG - Despiking routine algorithm ;; 1 = Sigma filtering ;; 2 = imgclean (VIS only) ;; BADPIXSINTERP - Should the bad pixels be interpolated over. ;; MISSINGINTERP - Should missing data be interpolated over. ;; HLUTFN - The filename of the H lambda lookup table for IR images ;; COMPRESSMETH - The method used to distribute DN values in decompression ;; 1 - average value ;; 2 - gaussian distribution ;; 3 - random uniform distribution ;; MTFALG - Algorithm used for deconvolution ;; 1 - constrained least squares ;; 2 - Richardson-Lucy with wavelet noise suppression ;; MTFPARAM - Algorithm specific parameter for the deconvolution algorithm ;; ;; EXAMPLE: ;; IDL> DIcal,'in.fit',radFn='rad.fit',/VERBOSE ;; ;; PROCEDURES USED (i.e. called directly!): ;; READFITS - Reads a FITS file ;; WRITEFITS - Writes a FITS file ;; SXPAR - Reads a FITS parameter ;; DIcal_VIS - Calibrates a VIS image ;; DIcal_IR - Calibrates an IR image ;; ;; MODIFICATION HISTORY: ;; 2004-05-25 M. Desnoyer Created ;; 2005-04-05 B. Carcich Modified to do check_FITS,...,/SILENT ;; 2005-04-12 M. Desnoyer Removed call to compile routines ;; 2005-05-26 M. Desnoyer Added new options to despiking on deconvolution ;; ;;----------------------------------------------------------------------------- @dical_flags doverbose = keyword_set(VERBOSE) IF not doverbose then begin ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;BTCarcich2005-04-14 ;;; Use saved quiet state savequiet=!QUIET !QUIET = 1 ENDIF ; Do error handling ;CATCH, error error=0 curStep = 0 IF error NE 0 THEN BEGIN ;; Throw an error if an error causes things to stop CATCH, /CANCEL IF keyword_set(SILENT) THEN $ message, /NOPRINT message, inFn+': '+!ERROR_STATE.MSG ENDIF ;; Make sure an output file is specified dn = keyword_set(dnfn) ioverf = keyword_set(iffn) revIoverf = keyword_set(revIffn) rad = keyword_set(radfn) revRad = keyword_set(revRadfn) IF (rad EQ 0) AND (ioverf EQ 0) AND (dn EQ 0) AND (revIoverf EQ 0) AND (revRad EQ 0) THEN $ message, 'Output filename must be specified', /NONAME ;; Open the FITS file imgIn = readFITS(inFn,hdr,/silent) ;; Make sure we got a valid file IF size(imgIn,/n_dimensions) NE 2 THEN $ message, 'Invalid FITS file', /NONAME $ ELSE $ IF doverbose THEN print, 'File opened' ;; Get the flag map catch, valid IF valid EQ 0 THEN BEGIN flags = di_readfits_flagmap(inFn) ENDIF ELSE BEGIN sz = size(imgIn,/dimensions) flags = intarr(sz[0],sz[1]) message, /reset ENDELSE catch, /cancel ;; Get the detector dtctr = sxpar(hdr,'IMGH021', COUNT=c1) ;; Make sure the items needed were in the FITS header IF (c1 EQ 0) THEN message, 'Invalid FITS header', /NONAME ;; IF activeModules structure was not specified in modules argument, ;; THEN this IF-ENDIF clause sets the defaults by creating one and ;; filling it in to determine which modules are to be applied IF keyword_set(modules) EQ 0 THEN BEGIN ;;; Create a new activeModules structure with all modules turned off; ;;; the code which follows then turns on a default set of modules modules = {activeModules} ;;; Temperature and Voltage re-calibrations are turned off by default modules.tempvolt = 0 ;; Some modules default to be on unless the ,/CHOSEN_ONLY keyword ;; argument was specified, so set NotChosen_only to be its inverse NotChosen_only = keyword_set(CHOSEN_ONLY) EQ 0 ;;; Perform flagging of saturated pixels and crosstalk correction ;;; by default if ,/CHOSEN_ONLY was *NOT* specified modules.flagSat = NotChosen_only modules.xTalk = NotChosen_only ;;; Perform normalization of quadrant gains and/or uneven bit weighting ;;; by default only if their respective filenames were specified modules.gain = keyword_set(GAINFN) modules.bitWeight = keyword_set(adcfn) ;; Modules that require more complex logic ;;; Perform dark subtraction by default if any ;;; of the following are true: modules.dark = NotChosen_only $ ;;; No ,/CHOSEN_ONLY OR keyword_set(DARKFN) $ ;;; Dark file specified OR (n_elements(USEDARKMODEL) GT 0) ;;; Dark Model specified ;;; Perform DN linearization by default if either of the first two ;;; following items are true, ***FOR IR DATA ONLY*** modules.linDN = (NotChosen_only $ ;;; No ,/CHOSEN_ONLY OR keyword_set(lindnfn) $ ;;; Linearization file specfified ) AND (dtctr EQ 1) ;;; IR data ONLY ;;; Perform flat fielding by default ;;; if either of the following are true: modules.flat = NotChosen_only $ ;;; No ,/CHOSEN_ONLY OR keyword_set(FLATFN) ;;; Flat file specified ;;; Perform Desmear algorithm by default if either ;;; - a desmear algorithm number was specified, ;;; - OR this is VIS and ,/CHOSEN_ONLY was *NOT* specified modules.desmear = n_elements(DESMEARALG) eq 1L $ ? (DESMEARALG NE 0) $ : (NotChosen_only AND (dtctr EQ 0)) ;;; ,/BADPIXS ,/FILLGAPS ,/DESPIKE ,/DENOISE ,/GEOM ,/MTF ;;; Perform flagging of bad pixels and/or filling of gaps by default ;;; if either their respective keywords were specified, OR if ;;; ,/CHOSEN_ONLY was *NOT* specified; ;;; also, specify the algorithm number from the keyword value modules.badPixs = n_elements(BADPIXS) eq 1L ? BADPIXS : NotChosen_Only modules.fillGaps = n_elements(FILLGAPS) eq 1L ? FILLGAPS : NotChosen_Only ;;; Perform the rest only if their respective keywords were set; ;;; also, specify the algorithm number from the keyword value modules.despike = n_elements(DESPIKE) eq 1L ? DESPIKE : 0 modules.denoise = n_elements(DENOISE) eq 1L ? DENOISE : 0 modules.geom = n_elements(GEOM) eq 1L ? GEOM : 0 modules.mtf = n_elements(MTF) eq 1L ? MTF : 0 ENDIF ;; Insert all of the calibration FITS keywords with default values dical_insertCalKeywords, hdr sxaddpar, hdr, 'DATE', dical_systime2utc(), ' UTC creation time for this FITS file' ;; For adding to headers later on stars='***********************************************************************' ptr_free,ptr_new(rtnPY,/no_copy) ;; Determine which calibration pipeline to enter based on dectector IF dtctr EQ 0 THEN BEGIN IF doverbose THEN print, 'Calibrating VIS image' imgOut = DIcal_VIS(imgIn,hdr,modules, verbose=verbose, silent=silent,$ force=force, maxgapsize=maxgapsize, desmearalg=desmearalg, $ dn=dn, rad=rad, revrad=revrad, ioverf=ioverf, revioverf=revioverf,$ flatfn=flatfn, darkfn=darkfn, usedarkmodel=usedarkmodel, $ compressfn=compressfn, adcfn=adcfn, badpixsfn=badpixsfn, $ geomfn=geomfn, mtffn=mtffn, visconstfn=visconstfn, flags=flags,$ gainfn=gainfn, spikesig=spikesig, spikeiter=spikeiter, spikealg=spikealg, $ spikebox=spikebox, spikemed=spikemed, badpixsinterp=badpixsinterp, $ compressmeth=compressmeth, missinginterp=missinginterp, mtfalg=mtfalg,$ mtfparam=mtfparam, snr=snr, rtnPY=rtnPY) ENDIF ELSE IF dtctr EQ 1 THEN BEGIN IF doverbose THEN print, 'Calibrating IR image' imgOut = DIcal_IR(imgIn,hdr,modules, verbose=verbose, silent=silent,$ force=force, maxgapsize=maxgapsize, rad=rad, revrad=revrad,$ ioverf=ioverf, revioverf=revioverf, dn=dn, darkfn=darkfn, $ usedarkmodel=usedarkmodel, compressfn=compressfn, adcfn=adcfn,$ badpixsfn=badpixsfn, mtffn=mtffn, gainfn=gainfn,$ lindnfn=lindnfn, specfn=specfn, constmapfn=constmapfn, $ spikesig=spikesig, spikeiter=spikeiter, spikebox=spikebox, $ spikemed=spikemed, spikealg=spikealg, flags=flags, specMap=specMap, $ badpixsinterp=badpixsinterp, hlutfn=hlutfn, flatfn=flatfn, $ compressmeth=compressmeth, missinginterp=missinginterp, mtfalg=mtfalg,$ mtfparam=mtfparam, snr=snr, debug=debugArg) ;; Make the headers for the extensions mkhdr, specHdr0, float(specMap[*,*,0]), /image sxaddpar, specHdr0, 'BSCALE', 1.0, ' Data scaling factor' sxaddpar, specHdr0, 'BZERO', 0.0, ' Data offset value' sxaddpar, specHdr0, 'BUNIT', 'MICROMETER', ' Physical Units of Data' sxaddpar, specHdr0, 'COMMENT', stars specHdr1 = specHdr0 sxaddpar, specHdr0, 'COMMENT', '** This extension is a wavelength map for the' sxaddpar, specHdr0, 'COMMENT', '** primary image array' sxaddpar, specHdr0, 'COMMENT', stars sxaddpar, specHdr1, 'COMMENT', '** This extension is a spectral bandwidth map for the' sxaddpar, specHdr1, 'COMMENT', '** primary image array' sxaddpar, specHdr1, 'COMMENT', stars sxaddpar, specHdr0, 'EXTNAME', 'WAVELENGTH', $ ' Unique name for this image extension' sxaddpar, specHdr1, 'EXTNAME', 'DELTA_WAVELENGTH', $ ' Unique name for this image extension' ENDIF ELSE message, 'Invalid detector ID', /NONAME ;; Create the header for the SNR map mkhdr, snrHdr, float(snr), /image sxaddpar, snrHdr, 'COMMENT', stars sxaddpar, snrHdr, 'COMMENT', '** This extension is a signal-to-noise ratio map for the' sxaddpar, snrHdr, 'COMMENT', '** primary image array' sxaddpar, snrHdr, 'COMMENT', stars sxaddpar, snrHdr, 'BSCALE', 1.0, ' Data scaling factor' sxaddpar, snrHdr, 'BZERO', 0.0, ' Data offset value' sxaddpar, snrHdr, 'BUNIT', 'N/A', ' Data are unitless' sxaddpar, snrHdr, 'EXTNAME', 'SNR', ' Unique name for this image extension' ;;;BTCarcich_201105 Added ProfileY extension if n_elements(rtnPY) gt 2L then begin mkhdr, pyHdr, float(rtnPY), /image sxaddpar, pyHdr, 'COMMENT', stars sxaddpar, pyHdr, 'COMMENT', '** This extension records the de-stripe correction as a 2xN pixel' sxaddpar, pyHdr, 'COMMENT', '** array, where N is the number of rows in the primary image (e.g.,' sxaddpar, pyHdr, 'COMMENT', '** 1024 or 512). The extension has units of DN. The first column' sxaddpar, pyHdr, 'COMMENT', '** of this de-stripe extension is the one-dimensional array of' sxaddpar, pyHdr, 'COMMENT', '** corrections subtracted from image columns <=M, where M is the' sxaddpar, pyHdr, 'COMMENT', '** middle column of the primary image (i.e., 512 for a 1024 column' sxaddpar, pyHdr, 'COMMENT', '** image, 1-based indexing). The second column of this de-stripe' sxaddpar, pyHdr, 'COMMENT', '** extension was subtracted from image columns numbered >M.' sxaddpar, pyHdr, 'COMMENT', '** Note that the serial over-clock columns that border the primary' sxaddpar, pyHdr, 'COMMENT', '** image array are not modified by the de-stripe procedure.' sxaddpar, pyHdr, 'COMMENT', '**' sxaddpar, pyHdr, 'COMMENT', '** If the de-stripe correction was not applied then all values in this' sxaddpar, pyHdr, 'COMMENT', '** extension are set to 0.' sxaddpar, pyHdr, 'COMMENT', stars sxaddpar, pyHdr, 'BUNIT', 'DN', 'Data unit' sxaddpar, pyHdr, 'EXTNAME', 'DESTRIPE', 'Unique name for this image extension' endif ;;BTCarcich_20081126 Ensure SOC & POC pixels do not change rtn = getPOCrows( hdr[*,0], COMPLEMENT=imgRows) rtn = getBIAScols( hdr[*,0], COMPLEMENT=imgCols) ;; Create a list of the files to create fns=['','', ''] IF dn THEN fns = [[fns], [dnfn, 'DN', 'dn']] IF rad THEN fns = [[fns], [radfn, 'RAD', 'r']] IF revrad THEN fns = [[fns], [revradfn, 'RADREV', 'rr']] IF ioverf THEN fns = [[fns], [iffn, 'IF', 'i']] IF revioverf THEN fns = [[fns], [reviffn, 'IFREV', 'ir']] n = n_elements(fns[0,*]) - 1 ;; Write the output FITS file(s) FOR i=1L, n DO BEGIN idata = i-1L fn = fns[0,i] caltype = fns[1,i] caltypeShort = fns[2,i] lclhdr = hdr[*,idata] ;; Make modifications to the header fxaddpar, lclhdr, 'CALTYPE', caltype $ , ' REVersible, DN=data#, RADiance, IF=I/F' $ , after='bunit' tmp = getBiasCols(lclhdr, complement=imgCols) tmp = getPOCRows(lclhdr, complement=imgRows) subWindow = imgOut[imgCols[0]+1:imgCols[1]-1, imgRows[0]+1:imgRows[1]-1 ,idata] fxaddpar, lclhdr, 'DATAMIN', min(subWindow), ' Minimum pixel in image array' fxaddpar, lclhdr, 'DATAMAX', max(subWindow), ' Maximum pixel in image array' fxaddpar, lclhdr, 'MEDPVAL', median(subWindow), ' Median pixel in image array' fxaddpar, lclhdr, 'STDPVAL', stddev(subWindow), ' Standard Deviation of pixel values' fxaddpar, lclhdr, 'PSATVAL', -999 fxaddpar, lclhdr, 'PSATNUM', -999 keyNames = ['BADPXCT','MISSPXCT','DESPIKCT','INTERPCT','PSATPXCT','SATPXCT','ASATPXCT','ULTCMPCT'] for j=0,7 do begin w = where((flags[*,*,idata] AND ishft(1,j)) NE 0,flgCnt) fxaddpar, lclhdr, keyNames[j],flgCnt endfor ;;; Extract FITS output file basename (i.e. excluding directory) ;;; and put into FILENAME and FILESDC keywords ;;;ps = path_sep_for_most_idl_versions() ps = path_sep() if size(fn,/type) eq size('',/type) then begin fitnametokens = strsplit(fn,ps,/extract) fitfilename = fitnametokens[n_elements(fitnametokens)-1] endif else begin fitfilename = 'UNKNOWN' endelse pdsmastfn = dical_getpdsmastfn(lclhdr,caltypeShort+'.fit') fxaddpar, lclhdr, 'FILENAME', pdsmastfn $ , ' FITS filename for PDS/MAST', before='FILENAMR' fxaddpar, lclhdr, 'FILESDC', fitfilename $ , ' FITS filename for SDC', after='FILENAMR' imgOutTmp = float(imgOut[*,*,idata]) check_FITS, imgOutTmp, lclhdr, /UPDATE, /FITS, /SILENT writeFITS, fn, imgOutTmp, lclhdr di_writefits_flagmap, fn, flags[*,*,idata] IF dtctr EQ 1 THEN BEGIN writeFITS, fn, float(specMap[*,*,0]), specHdr0, /append writeFITS, fn, float(specMap[*,*,1]), specHdr1, /append ENDIF writeFITS,fn,float(snr),snrHdr, /append if n_elements(rtnPY) gt 2L then writeFITS,fn,float(rtnPY),pyHdr, /append IF doverbose THEN print, 'Calibrated file written: '+fn ENDFOR ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;BTCarcich2005-04-14 ;;; Use saved quiet state ;;;!QUIET = 0 IF not doverbose THEN BEGIN !QUIET = savequiet ENDIF ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; RETURN END