PRO mk_xrt_composite, l_idx, l_da, s_idx, s_da, c_idx, c_da, $ despike=despike, despot=despot, lsat_map=lsat_map, verbose=verbose, $ gmap_l=gmap_l, gmap_s=gmap_s, gmap_c=gmap_c, gc_idx=gc_idx ; ========================================================================= ;+ ; PROJECT: ; Solar-B / XRT ; ; NAME: ; ; MK_XRT_COMPOSITE ; ; CATEGORY: ; ; Data graphics. ; ; PURPOSE: ; ; Combine a SINGLE long/short exposure pair of images to increase the ; dynamic range and to replace saturation with unsaturated data. ; See Note #1. ; ; CALLING SEQUENCE: ; ; MK_XRT_COMPOSITE, l_idx, l_da, s_idx, s_da, c_idx, c_da ; [,despike=despike] [,despot=despot] ; [,lsat_map=lsat_map] [,/verbose] ; [,gmap_l=gmap_l, gmap_s=gmap_s] [, gmap_c=gmap_c] ; [,gc_idx=gc_idx] ; ; INPUTS: ; ; L_IDX - [Mandatory] (structure scalar) ; This is the "index" for the long exposure image. ; Can be Level-0 or Level-1, as long as it matches L_DA. ; L_DA - [Mandatory] (2D number array, [Nx,Ny]) ; This is the 2D data array for the long exposure image. ; Can be Level-0 or Level-1, as long as it matches L_IDX. ; S_IDX - [Mandatory] (structure scalar) ; This is the "index" for the short exposure image. ; Can be Level-0 or Level-1, as long as it matches S_DA. ; S_DA - [Mandatory] (2D number array, [Nx,Ny]) ; This is the 2D data array for the short exposure image. ; Can be Level-0 or Level-1, as long as it matches S_IDX. ; ; KEYWORDS: ; ; DESPIKE - [Optional] (Boolean or integer) Set if 1 pass of ; despike wanted to remove radiation-belt/cosmic-ray ; spikes, or set to 1-3 passes of despike. (Two passes ; works well even for heavy spiking.) Default = no ; despiking if the keyword is not set. ; (See for more information.) ; DESPOT - [Optional] (Boolean) If set, composite image will have ; cosmetic correction to contamination spots using ; xrt_spotcor [Note: currently we do only include the option ; to use a non-default threshold value.] ; LSAT_MAP - [Optional] (Byte array) Set to a boolian map of the ; location of saturated pixels to be replaced. This can ; be used if the data have already been xrt_prepped, though ; supplying grade maps is the preferred method. Array ; must be same size as images to be composited. 1 represents ; a pixel to be replaced, 0 indicates don't replace. ; GMAP_L - [Optional] (Byte array) Grade map for the long exposure ; image, generally generated by xrt_pixel_grade, which holds ; values that indicate the grade of the pixel, saturated, ; bleed, spot, etc.. If provided will be used to determine ; which pixels will be replaced by short exposure pixels ; (saturated and bleed pixels). Must match size of images ; GMAP_S - [Optional] (Byte array) Grade map for the short exposure ; image, generally generated by xrt_pixel_grade, which holds ; values that indicate the grade of the pixel, e.g. ; saturated, bleed, spot. Useful if an output grade map is ; going to be used which will hold the grades of the pixels ; depending on which image they come from. ; /VERBOSE - [Optional] (Boolean) If set, print out extra messages. ; ; OUTPUTS: ; ; C_IDX - [Mandatory] (structure scalar) ; This is the "index" for the composite image. It will ; have the same content as L_IDX (after it has been ; through and normalized), except that ; there will be an update to the HISTORY field. ; C_DA - [Mandatory] (2D number array, [Nx,Ny]) ; This is the 2D data array for the composite image. ; GMAP_C - [Optional] (Byte array) Grade map for composite image. ; If grade maps have been provided then this map uses pixel ; values in the input grade maps. If pixels in the long ; exposure image are replaced then the grades for those ; pixels correspond to those in the short exposure map. The ; non-replaced pixels have the corresponding grades of the ; long exposure pixels. Note: if lsat_map is input and the ; grade maps are not, then a pseudo grade map is created that ; contains only the saturated and bloom grades for pixels ; in the long exposure map. Since information is not ; available in this case for the short exposure image, this ; would allow a user to at least know which pixels were ; replaced in the composite image. ; GC_IDX - [Optional] (scalar structure) Output header structure ; created for the grade map. ; ; EXAMPLES: ; ; Basic usage: ; IDL> mk_xrt_composite, l_idx, l_da, s_idx, s_da, c_idx, c_da ; IDL> tvscl, xrtdisp(c_da, /log) ; ; Recommended to do despiking on Level-0 data: ; IDL> mk_xrt_composite, l_idx, l_da, s_idx, s_da, c_idx, c_da, $ ; IDL> despike=2 ; IDL> tvscl, xrtdisp(c_da, /log) ; ; COMMON BLOCKS: ; ; none ; ; NOTES: ; ; 1) A simple proceedure for combining long and short XRT images. ; This procedure simply generates a composite image for XRT ; data. Only makes 1 composite at a time. It checks to see ; if the images have been xrt_prepped, and preps them if ; they have not been. The index structure for the long image ; will be used for the combined index structure with the ; history updated with the short image information. ; 2) We have retained the old method of providing a saturation map ; though the preferred method is to provide grade maps as ; generated by xrt_prep. If a saturation map is provided then ; the returned grade map only includes saturation pixels and doesn't ; include dust/spots etc. and so is not a true grade map. If the ; input images are level 0 (i.e. not prepped) then the prepping ; creates the grade maps. ; ; CONTACT: ; ; Comments, feedback, and bug reports regarding this routine may be ; directed to this email address: ; xrt_manager ~at~ head.cfa.harvard.edu ; ; MODIFICATION HISTORY: ; progver = 'v2006-May-09' ;--- (P.Jibben) Written. progver = 'v2007-Jun-13' ;--- (P.Jibben, S.Saar) This will handle the ; bloom pixels in a special way without ; altering the other data. Also changed ; the way the history is saved, as well ; as changed the data_lev keyword to 2. progver = 'v2007-Jun-28' ;--- (P.Jibben) This will allow users to pass ; their own saturated pixel map so they are ; not forced to use xrt_prep. progver = 'v2007-Jul-13' ;--- (M.Weber) Do not replace saturated ; pixels with values below the ; minimum value of the long exposure. ; I.e., preserve minimum level of long ; exposure. progver = 'v2009-Oct-09' ;--- (S.Saar) Updated to reflect changes in ; the xrt_prep code regarding the saturated ; pixel array. Added check for case where ; input lsat_map=0. progver = 'v2011-Oct-20' ;--- (S.Saar) Removed note 2 (unnecessary, as program ; checks for normalization) progver = 'v2014-Mar-26' ;--- (J.Slavin) Added despot keyword. Also made it ; so that if despike keyword is given, despiking ; is done even if the data has been processed ; (i.e. it's Level-1 data) progver = 'v2014-Mar-31' ;--- (J.Slavin) Added two extra tags, LNGFNAME and ; SRTFNAME, the L1 file names for the long and ; short images respectively, to the composite ; image header index to allow easier access to ; those filenames progver = 'v2014-Aug-14' ;--- (J.Slavin) Added check for negative short ; exposure pixels when replacing long exposure ; pixels -- deals with issue of replacing ; saturated pixels with bad values progver = 'v2014-Sep-12' ;--- (J.Slavin) included code from mk_xrt_comp3.pro ; that includes bad pixels (i.e. from data ; dropouts -- missing frames) in saturation map ; so missing data replaced by short image data progver = 'v2016-Nov-21' ;--- (J.Slavin) fixed bug in constructing bloom/bleed ; pixel image that allowed pixels to be both bloom ; and saturated. Also added optional output ; parameter BLOOMAP that holds the map of ; bloom/bleed pixels progver = 'v2017-Feb-02' ;--- (J. Slavin) Removed BLOOMAP parameter in favor ; of GMAP_L, GMAP_S parameters that allow input ; of grade maps which include both saturation and ; bleed/bloom information. Also added gmap_c ; optional output parameter to return the grade ; map for the composite image. Changed treatment ; of bloom/bleed pixels - now just replace them ; with pixels from short exposure. progver = 'v2017-Feb-06' ;--- (J. Slavin) Fixed bug for case where saturation ; map is provided. Composite grade map should ; hold grades of pixels as replaced. For case ; where user supplies saturation map, the ; composite grade map holds saturation and bloom ; grades of the long exposure. ;- ; ========================================================================= progname = 'MK_XRT_COMPOSITE' despike = keyword_set(despike) despot = keyword_set(despot) q_lsat_map = keyword_set(lsat_map) q_gmap_l = keyword_set(gmap_l) q_gmap_s = keyword_set(gmap_s) ; This is not strictly needed, but the whole idea of providing grade maps is ; to create a composite grade map, which requires both maps if (q_gmap_l ne q_gmap_s) then begin print,'Both grade maps (for long and short exposure) must be ' $ + 'provided for grade maps to be used... Returning.' return endif loud=keyword_set(verbose) if (not keyword_set(verbose)) then quiet=1 ; Added this bit so that inputs are not overwritten. l_idx1 = l_idx l_da1 = l_da s_idx1 = s_idx s_da1 = s_da if q_lsat_map then l_map=lsat_map ; determine if prepping is necessary. l_prep = get_history(l_idx1,'XRT_PREP',found=l_found) s_prep = get_history(s_idx1,'XRT_PREP',found=s_found) if loud then begin if l_found then print,'Long exposure is prep''d' $ else print,'Long exposure is not prep''d' if s_found then print,'Short exposure is prep''d' $ else print,'Short exposure is not prep''d' endif ; If you prep one you must prep both. if (l_found ne s_found) then begin print,'Long/Short images must have same type of processing... Returning.' return endif ; If the images are prepped you must provide the saturated pixel map or ; the grade maps. if (l_found and (not q_lsat_map) and (not q_gmap_l)) then begin print,'Must provide either saturation map or grade maps for ' $ + 'xrt_prep''d data... Returning.' return endif ; Make sure that l_da & lsat_map dimensions match. if q_lsat_map then begin sz_map = size(l_map,/dimensions) sz_lda = size(lda1,/dimensions) ndim1 = n_elements(l_da1) ndim2 = n_elements(l_map) if (ndim1 ne ndim2) then begin print,'Error long image and sat_map array must have same ' $ + 'dimensions... Returning.' return endif endif ; Make sure that l_da & gmap_l & gmap_s dimensions match. if q_gmap_l then begin sz_mapl = size(gmap_l,/dimensions) sz_maps = size(gmap_s,/dimensions) sz_lda = size(lda1,/dimensions) ndim1 = n_elements(l_da1) ndim2 = n_elements(gmap_l) ndim3 = n_elements(gmap_s) if (ndim1 ne ndim2) or (ndim2 ne ndim3) then begin print,'Error: images, gmap_l and gmap_s must have same ' $ + 'dimensions... Returning.' return endif endif ; Need to make sure that images are normalized for exposure if xrt_prep is ; used. If xrt_prep not used, assume user normalized images themselves. l_norm = get_history(l_idx1,'XRT_RENORMALIZE',found=nl_found) s_norm = get_history(s_idx1,'XRT_RENORMALIZE',found=ns_found) ; Renormalize xrt_prep'd exposures if they haven't already been renormalized if (l_found and (nl_found eq 0)) then begin l_da1 = temporary(l_da1)/get_xrt_expdur(l_idx1, /sec) l_idx1.e_etim = 1000000 ; [microsec] for consistancy with XRT_PREP ; & returned index if (loud eq 1) then print,'Renormalizing the long exposure.' endif if (s_found and (ns_found eq 0)) then begin s_da1 = temporary(s_da1)/get_xrt_expdur(s_idx1, /sec) if (loud eq 1) then print,'Renormalizing the short exposure.' endif ; Use xrt_prep on images... If needed. smiss = 0 if (not s_found) then begin if loud then print,'Prepping short for composite' if despike and despot then begin despike_despot = 1 only_despike = 0 only_despot = 0 endif else if despike and (not despot) then begin despike_despot = 0 only_despike = 1 only_despot = 0 endif else if despot and (not despike) then begin despike_despot = 0 only_despike = 0 only_despot = 1 endif else begin despike_despot = 0 only_despike = 0 only_despot = 0 endelse gmap_s = 1 xrt_prep, temporary(s_idx1), temporary(s_da1), s_idx1, s_da1, $ despike_despot=despike_despot, only_despot=only_despot, $ only_despike=only_despike, /float, /normalize, grade_map=gmap_s, $ verbose=verbose, quiet=quiet q_gmap_s = 1B ;***** treatment of the SHORT missing frames ***** ; Note: missing frames show up as 0 values in level0 files and -999 in ; level1 files. It may be that despiking and/or despotting could change ; the value in the processed images to something other than -999 if tag_exist(s_idx,'PIXCNT') eq 1 then $ if s_idx.PIXCNT lt n_elements(s_da) then $ s_mfr = where(s_da eq 0,smiss) endif else begin if tag_exist(s_idx1,'PIXCNT') eq 1 then $ if s_idx.PIXCNT lt n_elements(s_da) then $ s_mfr = where(s_da eq -999,smiss) endelse lmiss = 0 if (not l_found) then begin if loud then print,'Prepping long for composite' gmap_l = 1 ; Note: we use keyword values for despike_despot, etc. as determined ; above for the short exposure (since if s_found eq 0 then l_found eq 0) xrt_prep, temporary(l_idx1), temporary(l_da1), l_idx1, l_da1, $ despike_despot=despike_despot, only_despot=only_despot, $ only_despike=only_despike, /float, /normalize, grade_map=gmap_l, $ verbose=verbose, quiet=quiet q_gmap_l = 1B ;***** treatment of the LONG missing frames ***** if tag_exist(l_idx,'PIXCNT') eq 1 then $ if l_idx.PIXCNT lt n_elements(l_da) then $ l_mfr = where(l_da eq 0,lmiss) endif else begin l_mfr = where(l_da eq -999,lmiss) endelse if l_found and (despike or despot) then begin ; assume here that if the long exposure was processed then ; all of them were -- also all assumed to be the same with ; regard to despiking and despotting if despike then begin l_da_out = xrt_despike2(l_da1,l_idx1,hist=lhist) s_da_out = xrt_despike2(s_da1,s_idx1,hist=shist) update_history,l_idx1,lhist,/noroutine update_history,s_idx1,shist,/noroutine l_da1 = l_da_out s_da1 = s_da_out endif if despot then begin xrt_spotcor,l_idx1,l_da1,l_idx_out,l_da_out,quiet=quiet xrt_spotcor,s_idx1,s_da1,s_idx_out,s_da_out,quiet=quiet l_da1 = l_da_out l_idx1 = l_idx_out s_da1 = s_da_out s_idx1 = s_idx_out endif endif ; Note: this overrides a saturation map provided by lsat_map if both ; grade maps and a saturation map are provided (which makes no sense) if q_gmap_l then begin ; use bit arithmetic to extract saturated pixels l_map = (gmap_l and 1B) ; use bit arithmetic to extract bleed/bloom pixels ; (Note: the "ne 0B" part causes map to contains 1s rather than 2s) bl_map = ((gmap_l and 2B) ne 0B) endif ; Flag missing frames for replacement in long exposure images ; For short exposure image, replace image values with those from long if lmiss gt 0 then l_map[l_mfr] = 1B if smiss gt 0 then s_da1[s_mfr] = l_da1[s_mfr] ; check if there are any saturated pixels n_sat_l = long(total(l_map)) ; In case the grade maps are not provided gmap_c just contains the ; saturation map - but will end up being mostly 0's once saturated pixels ; are replaced if q_gmap_l then gmap_c = gmap_l else gmap_c = l_map ; Use sat_map to replace saturated pixels of long exposure. if (n_sat_l gt 0) then begin ; make sure short exposure pixels are good (i.e. > 0) locations ; where we replace saturated pixels in long exposure image ss = where((l_map eq 1) and s_da1 gt 0.,nss) if nss gt 0 then begin l_da1[ss] = (s_da1[ss] > min(l_da1)) ; If grade maps not provided then leave grade map as saturation ; map if q_gmap_s then gmap_c[ss] = gmap_s[ss] endif sat_tag = progname + ': Replaced ' + strcompress(nss,/rem) $ + ' saturated pixels' c_da = l_da1 ; Handling of bloom/bleed pixels - changed from previous behavior - now ; simply replace bloom/bleed pixels in long exposure with pixels from ; short exposure the same as for saturated pixels if q_gmap_l then begin ss = where(bl_map and (s_da1 gt 0),nbs) if nbs gt 0 then begin ; use long exposure pixel value as floor - saturation will ; cause the pixel values to be too low in general c_da[ss] = (s_da1[ss] > min(l_da1)) gmap_c[ss] = gmap_s[ss] endif bloom_tag = progname + ': Replaced ' + strcompress(nbs,/rem) $ + ' bloom pixels' endif else begin ; Retaining this method for case where we don't have a grade map: ; Assume pixels on either side in y of saturated pixels are ; bloom/bleed pixels ; In this case grade map only includes saturation and bloom ; Taken almost verbatim from xrt_pixel_grade n1 = l_idx1.naxis1 n2 = l_idx1.naxis2 mask_edge = intarr(n1,n2) + 1 mask_edge[*,[0,n2-1]] = 0 ; mask to kill wrap-around bloom = mask_edge*(shift(l_map,0,1) + shift(l_map,0,-1)) bloom = (bloom*(l_map eq 0) < 1) ; ensure no double counts ; rpix = where(bloom eq 1 and (s_da1 gt 0),cnt) nbloom = strcompress(cnt,/rem) bloom_tag = progname + ': Replaced ' + nbloom + $ ' bloom pixels' if (cnt gt 0) then begin ; Note: we now simply replace bloom pixels with short exposure ; pixels (previously used convolved long exposure pixel values) c_da[rpix] = (s_da1[rpix] > min(l_da1)) ; Note: if only a saturation map is povided then output ; grade map holds the input saturation map plus bloom pixel ; grades gmap_c[rpix] = 2B endif endelse endif else begin c_da = l_da1 sat_tag = progname + ': No saturated pixels, long exposure used' bloom_tag = progname + ': No bloom pixels because no saturated pixels' message,'No saturated pixels to replace - returning long image', $ /continue endelse ;generating the c_idx c_idx = l_idx1 ; header keywords that provide easily accessible information ; on the filenames for the long and short exosure image files lfname = xrt_index2filename(l_idx) sfname = xrt_index2filename(s_idx) c_idx = add_tag2(c_idx,lfname,'LNGFNAME',index='OBS_MODE') c_idx = add_tag2(c_idx,sfname,'SRTFNAME',index='LNGFNAME') c_idx.data_lev = 2 ;tagval0 = '' ;change_tag_value,c_idx,tagval0,'HISTORY' ;tagval0 = [tagval0, progname + ' ' + progver + ': '] tagval0 = [progname + ': ' + progver] if q_gmap_l then tagval0 = [tagval0, $ progname + ': Saturated pixels taken from grade map' ] $ else tagval0 = [tagval0, progname + ': Saturated pixels taken ' $ + 'from saturation map'] tagvalg = tagval0 if despike then tagval0 = [tagval0, $ progname + ': despiking done on composite image using xrt_despike2'] if despot then tagval0 = [tagval0, $ progname + ': spot correction done on composite image' + $ ' using xrt_spotcor'] ; generate the output grade map header structure if q_gmap_l then begin tagvalg = [tagvalg, progname + ': Grade map created from ' + $ 'long and short exposure grade maps associated with ' + $ lfname + ' and ' + sfname] endif else begin tagvalg = [tagvalg, progname + ': Partial grade map including ' + $ 'only saturation and bloom grades created from user ' + $ 'provided long exposure saturation map.'] endelse long_hist = l_idx1.history short_hist = s_idx1.history tagval0 = [tagval0, sat_tag, bloom_tag, 'Long image history --', $ long_hist, 'Short image history --', short_hist] tagvalg = [tagvalg, sat_tag, bloom_tag, 'Long image history --', $ long_hist, 'Short image history --', short_hist] gc_idx = c_idx ;update_history,c_idx,tagval0,/noroutine ;update_history,gc_idx,tagvalg,/noroutine c_idx = rep_tag_value(c_idx,tagval0,'HISTORY') gc_idx = rep_tag_value(gc_idx,tagvalg,'HISTORY') return end