;----------------------------------------------------------------------------- PRO fg_flatfield, index, data, index_out, data_out, mask, $ flatdir=flatdir, $ user_mask=user_mask, user_file=user_file, $ user_flat=user_flat, user_index=user_index, $ flat_image=flat_image, flat_index=flat_index, $ no_float=no_float, badthresh=badthresh, $ version=progver, name=prognam, $ loud=loud, verbose=verbose, $ retrn=retrn, qstop=qstop, $ update=update, run_time=run_time ;+ ; NAME: ; FG_FLATFIELD ; ; PURPOSE: ; Divide SOT/FG images by an appropriate flat field image. ; ; CATEGORY: ; FITS processing ; ; SAMPLE CALLING SEQUENCE: ; fg_flatfield, index, image, index_out, image_out, mask $ ; [, user_mask=user_mask, flatdir=flatdir $ ; , user_flat=user_flat, user_index=user_index $ ; , flat_image=flat_img, flat_index=flat_index $ ; , version=version, name=name, /update ] ; ; INPUTS: ; INDEX - The index structure for each image, including /image ; IMAGE - The data array (NX x NY) ; ; OUTPUTS: ; INDEX_OUT - The index structure for the flatfielded image. ; IMAGE_OUT - The original image normalized by the flatfield ; image, converted to floating point unless ; /NO_FLOAT is set. ; ; OUTPUTS (Optional): ; MASK - a 2-D image showing valid and invalid pixels. This array ; is useful when filling the bad pixels in fg_bad_pix. ; ; OPTIONAL INPUT KEYWORD PARAMETERS: ; USER_MASK - A user-supplied mask of dead and/or hot pixels. ; The MASK array returned by FG_SHIFT_PIX can be used. ; FLATDIR - Alternate directory for processed flat current fits images, ; default = 'concat_dir('$SOT_SSWDB','flats')' for flat images ; USER_FILE - If set, read in a user supplied flat current fits ; image as the first image with filename given by 'user_file'. ; This image is assumed to be sized, binned, and pixel aligned ; already. ; USER_FLAT - A user supplied flat image can be input in this ; keyword. ; USER_INDEX - The index structure associated with user_flat. ; Mandatory if user_flat is supplied. ; BADTHRESH - Minimum and maximumu threshold in a gain image. ; Default: [0.2,15] ; NO_FLOAT - If set, the flatfielded image is not converted to floating point. ; VERBOSE - Set for a few more messages; overrides /loud ; LOUD - Not set to minimize the messages ; QSTOP - Set to stop in this procedure for debuging purposes ; UPDATE - If set, update history record in INDEX_OUT ; ; OUTPUT KEYWORD PARAMETERS: ; FLAT_IMAGE - The flat image used in the process ; FLAT_INDEX - The index structure associated with flat_image ; VERSION - The current version number of the program ; NAME - The name of this routine ; RETRN - Flag for branching on unusual return ; ; PROCEDURE: ; ; 1. Find an appropriate flat in the database. ; 2. Create the gain image and multiply the image by it. ; ; Perform steps in floating point arithmetic. Return image in ; floating point for subsequent processing unless NO_FLOAT is ; set. ; ; MODIFICATION HISTORY: ;V0.1 09-Oct-2006. Created by T. Berger, LMSAL ;v0.9 16-Apr-2007. Major update by YK, NAOJ ;V1.0 24-Apr-2007. Minor typing corrections. TEB. ;v1.1 Update to handle multi-image NFI data arrays. 24-Oct-2007 TEB. ;v1.2a 2x2 fix by YK added. ;v1.2b Fix for unity flatfields returned by FG_GET_FLAT. TEB ; 4-Jan-2008. ;v1.23 Minor bug fix on line 173. TEB 7-Jan-2008. ;v1.3 Modify flatfields for binned images. TEB 28-Jul-2008. ;v1.4 Bug fix to handle binned images. 28-Oct-2008, YK ;v1.4.1 Bug fix in extraction code. TEB 4-Feb-2009. ;v1.4.2 inhibit extraction if no flat field found (avoid unity array subscript error) ;v1.4.3 2020-07-28 - GLS - Check that flat_mask and imask are same size as image. ; TO DO: ; - Implement common block (or other) mechanism to avoid ; re-reading flat images for every image in a series. ;- ;----------------------------------------------------------------------------- ; ; Program name and version: prognam = 'FG_FLATFIELD.PRO' progver = 'V1.4.3' prog_nam_ver = 'FG_FLATFIELD.PRO V1.4.3' ;Constants xmax = 4095 ymax = 2047 x0 = 0 y0 = 0 time0 = SYSTIME(1) retrn = 0 loud = KEYWORD_SET(loud) verbose = KEYWORD_SET(verbose) if (verbose eq 1) then loud = 1 if (loud) then MESSAGE,/info, 'Start '+prognam+' '+progver if KEYWORD_SET(no_float) then nofloat = 1 else nofloat = 0 ;if image not floating point then image = FLOAT(image) ; keep original image ;; data0 = data ;; case data_chk(image,/type) of ;; 2: if (KEYWORD_SET(float)) then image = FLOAT(TEMPORARY(image)) ;; 4: if (not KEYWORD_SET(float)) then image = FIX(ROUND(TEMPORARY(image))) ;; 5: if (not KEYWORD_SET(float)) then image = FIX(ROUND(TEMPORARY(image))) ;; else: MESSAGE,/INFO,'Warning: Image is a strange type...' ;; endcase ;image characteristics img_dateobs = gt_tagval(index[0],/date_obs) ; image date/time img_genid = fg_genid(gt_tagval(index[0],/gen_id)) ; GEN_ID has a 1:1 mapping to obs_type img_wave = fg_waveid(gt_tagval(index[0], /waveid)) ; wavelength of observable img_exptime = gt_tagval(index[0],/exptime) ; exposure time img_nx0 = gt_tagval(index[0],/naxis1) ; image size nx0 img_ny0 = gt_tagval(index[0],/naxis2) ; image size ny0 img_ix0 = gt_tagval(index[0],/fgccdix0) ; image x0 in CCD coordinates img_ix1 = gt_tagval(index[0],/fgccdix1) ; image x1 in CCD coordinates img_iy0 = gt_tagval(index[0],/fgccdiy0) ; image y0 in CCD coordinates img_iy1 = gt_tagval(index[0],/fgccdiy1) ; image y1 in CCD coordinates img_gain = gt_tagval(index[0],/camgain) ; camera gain img_amp = gt_tagval(index[0],/camamp) ; ccd amplifier used img_daca = gt_tagval(index[0],/camdaca) ; ccd DAC A img_dacb = gt_tagval(index[0],/camdacb) ; ccd DAC B img_psum = gt_tagval(index[0],/campsum) ; Parallel summing (summing in x direction) img_ssum = gt_tagval(index[0],/camssum) ; Serial summing (in y-direction) img_t_fgccd = gt_tagval(index[0],/t_fgccd) ; Temperature of the FG CCD img_t_fgceb = gt_tagval(index[0],/t_fgceb) ; Temperature of the FG camera electronics box img_binx = gt_tagval(index[0],/fgbinx) ; Binning in x-direction img_biny = gt_tagval(index[0],/fgbiny) ; Binning in y-direction ;BFI or NFI image? if gt_tagval(index[0],/waveid) le 6 then nfi = 0 else nfi = 1 ;full or partial FOV readout? fullfov=0 if gt_tagval(index[0],/roistart) eq 0 AND gt_tagval(index[0],/roistop) eq 1025 then fullfov = 1 if (verbose) then begin MESSAGE,/info,'Input image information.' MESSAGE,/info,' Date/Time: '+img_dateobs MESSAGE,/info,' Wavelength: '+img_wave+' Angstrom' MESSAGE,/info,' Gen_id: '+img_genid MESSAGE,/info,' Exposure time: '+STRTRIM(img_exptime,2)+' seconds' MESSAGE,/info,' CCD gain: '+STRTRIM(img_gain,2) MESSAGE,/info,' CCD amp: '+STRTRIM(img_amp,2) MESSAGE,/info,' CCD summing: '+STRTRIM(img_psum,2)+' x '+STRTRIM(img_ssum,2)+' (parallel x serial)' MESSAGE,/info,' CCD binning: '+STRTRIM(img_binx,2)+' x '+STRTRIM(img_biny,2)+' (X x Y)' MESSAGE,/info,' CEB temp: '+STRTRIM(img_t_fgceb,2)+' deg C' MESSAGE,/info,' CCD temp: '+STRTRIM(img_t_fgccd,2)+' deg C' end ;get appropriate flatfield image. Depends only on wavelength. new_flat = 1 if KEYWORD_SET(user_flat) then begin ;Flat image and index supplied by keywords... if (not KEYWORD_SET(user_index)) then begin MESSAGE,'Need to supply SSW index structure with flat image...' RETURN end flat_image = user_flat flat_index = user_index flat_exptime = gt_tagval(flat_index[0],/exptime) flat_nx0 = gt_tagval(flat_index[0],/naxis1) flat_ny0 = gt_tagval(flat_index[0],/naxis2) if (flat_nx0 ge img_nx0) and (flat_ny0 ge img_ny0) then new_flat=0 $ else begin MESSAGE,/INFO,'Flatfield too small for image. Returning...' index_out = -1 image_out = -1 retrn=1 RETURN end end if KEYWORD_SET(user_file) then begin ; read in user supplied flat image filename read_sot, user_file, flat_index, flat_image flat_exptime = gt_tagval(flat_index[0],/exptime) flat_nx0 = gt_tagval(flat_index[0],/naxis1) flat_ny0 = gt_tagval(flat_index[0],/naxis2) if (flat_nx0 ge img_nx0) and (flat_ny0 gt img_ny0) then new_flat=0 $ else begin MESSAGE,/INFO,'Flatfield too small for image. Returning...' index_out = -1 image_out = -1 retrn=1 RETURN end end if (new_flat) then begin ;find flat file first time thru if (loud) then begin MESSAGE, /info, 'Reading nearest flat file in time from database...' endif ; get flat field if N_ELEMENTS(flatdir) eq 0 then flatdir = concat_dir(getenv('SOT_SSWDB_FG_CAL'),'flats') fg_get_flat, index, flat_index, flat_image, $ flatdir=flatdir, deltaday=deltaday, retrn=retrn, flatfile=flatfile if deltaday eq -1 then unitarr=1 else unitarr=0 if (retrn eq 1) then begin MESSAGE, /CONT, 'No flat file found. Returning uncalibrated image.' image_out = flat_image index_out = index retrn = 1 RETURN end if (loud) and (N_ELEMENTS(flat_image) gt 0) and (deltaday ne -1) then begin MESSAGE, /info, 'Flat file found in the database. ' MESSAGE, /info, ' File name: '+flatfile MESSAGE, /info, 'Time difference: '+STRTRIM(deltaday,2)+ ' days' endif flat_exptime = gt_tagval(flat_index[0],/exptime) flat_nx0 = gt_tagval(flat_index[0],/naxis1) flat_ny0 = gt_tagval(flat_index[0],/naxis2) if (flat_nx0 lt img_nx0) or (flat_ny0 lt img_ny0) then begin MESSAGE,/INFO,'Flat file too small for image. Returning...' index_out = index data_out = data retrn=1 RETURN endif endif ;Have a flat image now. Load the tag values: flat_img_dateobs = gt_tagval(flat_index[0],/date_obs) ; flat date and time in YYYY-MM-DDThh:mm:ss.mmm format flat_img_date = STRMID(flat_img_dateobs,0,10) ; flat image date flat_img_time = STRMID(flat_img_dateobs,12,12) ; flat image time flat_amp = gt_tagval(flat_index[0],/camamp) ; flat ccd amplifier used flat_gain = gt_tagval(flat_index[0],/camgain) ; flat ccd gain flat_daca = gt_tagval(flat_index[0],/camdaca) ; flat DAC A flat_dacb = gt_tagval(flat_index[0],/camdacb) ; flat DAC B flat_psum = gt_tagval(flat_index[0],/campsum) ; flat parallel summing factor flat_ssum = gt_tagval(flat_index[0],/camssum) ; flat serial summing factor flat_t_fgccd = gt_tagval(flat_index[0],/t_fgccd) ; flat ccd temperature flat_x0 = gt_tagval(flat_index[0],/fgccdix0) ; flat x0 flat_y0 = gt_tagval(flat_index[0],/fgccdiy0) ; flat y0 flat_binx = gt_tagval(flat_index[0],/fgbinx) ; flat binning in x flat_biny = gt_tagval(flat_index[0],/fgbiny) ; flat binning in y flat_wave = fg_waveid(gt_tagval(flat_index[0],/waveid)) ; wavelength of observable if (img_wave ne flat_wave) then $ MESSAGE,/INFO,'WARNING: flatfield wavelength ('+STRTRIM(flat_wave,2)+') does not match image wavelength ('+STRTRIM(img_wave,2)+').' ; check FG_SHIFT_PIX has been applied to the data previous = get_history(index[0], caller='FG_SHIFT_PIX',found=found) if (found) and (not nfi) then begin ;NFI flats already passed through shift_pix in creation phase. fg_shift_pix, flat_index, flat_image, flat_index2, flat_image2, flat_mask flat_index = TEMPORARY(flat_index2) flat_image = TEMPORARY(flat_image2) endif else flat_mask = BYTE(flat_image)*0b+1b ;HELP, flat_index, flat_image, flat_index2, flat_image2, flat_mask & $ ; STOP, ' Stopping in fg_flatfield.pro after getting a flat image.' ;Extraction and summing and/or binning. ; Note that we assume a user-supplied flat is already extracted and summed properly. ; Summing should be done after extraction if not KEYWORD_SET(user_flat) then begin ;Corner coordinates in IMAGE ;orientation (not CCD orientation) x0i = ( (4095-img_ix1)>0<4095 ) x1i = ( (4095-img_ix0)>0<4095 ) y0i = ( (2047-img_iy1)>0<2047 ) y1i = ( (2047-img_iy0)>0<2047 ) x0=x0i x1=x1i y0=y0i y1=y1i ;handle camera summing on the flat images if flat_psum ge 2 then begin if x0 mod flat_psum eq 0 then x0 = x0/flat_psum $ else x0 = float(x0)/flat_psum if x1-1 mod flat_psum eq 0 then x1 = (x1-1)/flat_psum $ else x1 = float(x1-1)/flat_psum endif if flat_ssum ge 2 then begin if y0 mod flat_ssum eq 0 then y0 = y0/flat_ssum $ else y0 = float(y0)/flat_ssum if y1-1 mod flat_ssum eq 0 then y1 = (y1-1)/flat_ssum $ else y1 = float(y1-1)/flat_ssum endif ;handle binning in the same way, possibly cumulative to summing... if flat_binx ge 2 then begin if x0 mod flat_binx eq 0 then x0 = x0/flat_psum $ else x0 = float(x0)/flat_binx if x1-1 mod flat_binx eq 0 then x1 = (x1-1)/flat_binx $ else x1 = float(x1-1)/flat_binx endif if flat_biny ge 2 then begin if y0 mod flat_biny eq 0 then y0 = y0/flat_biny $ else y0 = float(y0)/flat_biny if y1-1 mod flat_biny eq 0 then y1 = (y1-1)/flat_biny $ else y1 = float(y1-1)/flat_biny endif ;Extraction of sub-area. If floating point values, use interpolation. if (size(x0,/type) le 3) and (size(x1,/type) le 3) $ and (size(y0,/type) le 3) and (size(y1,/type) le 3) then begin if min(flat_image) ne max(flat_image) then begin flat_image = flat_image[x0:x1, y0:y1] ;v1.4.1: removed if statement that was preventing extraction. flat_mask = flat_mask[x0:x1, y0:y1] endif else box_message,'unity flat? proceeding..' endif else begin nnx = ROUND(x1-x0)+1 nny = ROUND(y1-y0)+1 xarr = (FINDGEN(nnx)+x0)#REPLICATE(1,nny) yarr = REPLICATE(1,nnx)#(FINDGEN(nny)+y0) flat_image = INTERPOLATE(flat_image,xarr,yarr) flat_mask = INTERPOLATE(flat_mask,xarr,yarr) endelse szxy = SIZE(flat_image,/dimen) flsx = szxy[0] flsy = szxy[1] if (flat_psum*flat_binx lt img_psum*img_binx) then begin flsx = flsx*flat_psum*flat_binx/(img_psum*img_binx) flsy = flsy*flat_ssum*flat_biny/(img_ssum*img_biny) flat_image = REBIN(flat_image, flsx, flsy) flat_mask = REBIN(flat_mask, flsx, flsy) if (loud) then MESSAGE,/info,'WARNING: flatfield is being rebinned to match image summing condition.' end flat_image = FLOAT(flat_image) flat_image = flat_image > 0 ; Now have a valid flat image end ;Make gain image ss = WHERE(flat_image le 0, ns) if ns gt 0 then begin flat_image[ss] = 1.0 flat_mask[ss] = 0 endif ; gain = flat_mean/flat_image gain = 1./flat_image ;Kuhn-Lin flatfields are pre-normalized to average 1. ; Check for bad pixels. ; Default threshold is [0.2, 15] if N_ELEMENTS(badthresh) eq 0 then badthresh = [0.2, 15] ss = WHERE(gain lt badthresh[0] or gain gt badthresh[1], ns) if ns gt 0 then begin gain[ss] = 1.0 ;just in case... flat_mask[ss] = 0 endif if (loud) then begin IMAGE_STATISTICS, gain, max=gain_max, min=gain_min, mean=gain_mean, stddev=gain_stdev $ ,count=zcount MESSAGE, /info, 'Gain image statistics: ' MESSAGE, /info, ' Min = '+STRTRIM(gain_min,2) MESSAGE, /info, ' Max = '+STRTRIM(gain_max,2) MESSAGE, /info, ' Mean = '+STRTRIM(gain_mean,2) MESSAGE, /info, ' Std Deviation = '+STRTRIM(gain_stdev,2) end if N_ELEMENTS(user_mask) le 0 then user_mask = BYTE(data)*0b+1b mask = user_mask nexp = fg_num_eff_exp(index, quiet=(1-loud)) nimage = N_ELEMENTS(nexp) ; GLS - Added 2022-06-01 to handle 3D data arrays: ndim_data = (size(data))[0] if ndim_data eq 2 then nimage = 1 if ndim_data eq 3 then nimage = n_elements(data[0,0,*]) ;if n_elements(data[0,0,*]) gt 1 then nimage = n_elements(data[0,0,*]) ;some gen_id's [33--39] have I and V bit-shifted to avoid ;saturation. Must therefore re-scale these images after flat-fielding: scale_factor = FLOAT(MAX(nexp))/nexp ; GLS - Added 2022-06-06 to handle 3D data arrays: if n_elements(scale_factor) eq 1 and n_elements(data[0,0,*]) gt 1 then $ scale_factor = replicate(scale_factor[0],3) if (KEYWORD_SET(qstop)) then STOP if (nofloat) then data_out = data else data_out = data*0.0 index_out = index ;Loop through images in data array for i=0,nimage-1 do begin image = data[*,*,i] if ( (n_elements(gain) ne n_elements(image)) or (n_elements(flat_mask) ne n_elements(image)) ) then begin MESSAGE,/INFO,' FG_FLATFIELD: Problem with size of gain image or flat_mask image. Returning...' index_out = index data_out = data retrn=1 RETURN endif image_type = fg_image_type(index,i) imask = flat_mask*mask[*,*,i] ; For 2x2 flatfield in partial FOV, need to mask out the other bright rows: if NOT (fullfov) AND (img_ssum eq 2) then $ if (img_binx lt 2) then begin imask[img_nx0/2-1,*]=0 imask[img_nx0/2+1,*]=0 end image = image * gain * scale_factor[i] image *= imask if (image_type eq 'FG') or (image_type eq 'SI') then image >= 0 if (nofloat) then image = FIX(image) ;Check image statistics if (loud) then begin IMAGE_STATISTICS,image,min=outmin,max=outmax,mean=outmean,stddev=outstdev,count=outcount MESSAGE, /info, 'Flatfielded image statistics:' MESSAGE, /info, ' Min = '+STRTRIM(outmin,2) MESSAGE, /info, ' Max = '+STRTRIM(outmax,2) MESSAGE, /info, ' Mean = '+STRTRIM(outmean,2) MESSAGE, /info, ' Std Deviation = '+STRTRIM(outstdev,2) end data_out[*,*,i] = image mask[*,*,i] = imask endfor run_time = SYSTIME(1) - time0 if (loud) then MESSAGE,/info, 'FG_FLATFIELD took '+strtrim(run_time,2)+' seconds for 1 image' ; Update history record in the index if requested if KEYWORD_set(update) then begin tagval = 'Applied flatfield correction: '+flatfile update_history,index_out, tagval, version=progver if (loud) then MESSAGE, /INFO, 'HISTORY record updated: '+tagval endif if (KEYWORD_SET(qstop)) then STOP RETURN END