PRO mk_xrt_comp3_unc, l_idx, l_da, m_idx, m_da, s_idx, s_da, c_idx, c_da, $ despike=despike, despot=despot, uncert_map=uncert_map, $ verbose=verbose ; ========================================================================= ;+ ; PROJECT: ; Solar-B / XRT ; ; NAME: ; ; MK_XRT_COMP3 ; ; CATEGORY: ; ; Data graphics. ; ; PURPOSE: ; Uncertainty array version of MK_XRT_COMPOSITE. ; Uncertainty arrays of the long/short images obtained thourgh ; XRT_PREP are made into a composite, according to the saturation ; maps of the image data. ; ; Combine a SINGLE long/medium/short exposure triple of images to ; increase the dynamic range and to replace saturation with ; unsaturated data. ; See Note #1. ; ; CALLING SEQUENCE: ; ; MK_XRT_COMP3_UNC, l_idx, l_da, m_idx, m_da, s_idx, s_da, $ ; c_idx, c_da, uncert_map=uncert_map, $ ; [,despike=despike] [,despot=despot] [,/verbose] ; ; 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. ; M_IDX - [Mandatory] (structure scalar) ; This is the "index" for the medium exposure image. ; Can be Level-0 or Level-1, as long as it matches M_DA. ; M_DA - [Mandatory] (2D number array, [Nx,Ny]) ; This is the 2D data array for the medium exposure image. ; Can be Level-0 or Level-1, as long as it matches M_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 not 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 of the long exposure image. ; Use this if the data have already been xrt_prepped and ; the map was generated or if xrt_prep is not to be used. ; Array must be same size as images to be composited. ; 1 represents the pixel to be replaced. ; MSAT_MAP - [Optional] (Byte array) Set to a boolian map of the ; saturated pixels of the medium exposure image. ; /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. ; UNCERT_MAP - [as a keyword] an array of uncertainties. ; ; EXAMPLES: ; ; Basic usage: ; IDL> mk_xrt_comp3_unc, l_idx, l_da, m_idx, m_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_comp3, l_idx, l_da, m_idx, m_da, s_idx, s_da, $ ; IDL> c_idx, c_da, despike=2, uncert_map=unc_map ; IDL> tvscl, xrtdisp(unc_map, /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. ; ; 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 = 'v2015-Mar-31' ;--- (A.Takeda) prepared as a modified version of ; MK_XRT_COMP3.PRO(v2014-Sep-12). ;- ; ========================================================================= progname = 'MK_XRT_COMP3' despike=keyword_set(despike) despot=keyword_set(despot) q_lsat_map=keyword_set(lsat_map) q_msat_map=keyword_set(msat_map) 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 m_idx1 = m_idx & m_da1 = m_da s_idx1 = s_idx & s_da1 = s_da if (q_lsat_map eq 1) then l_map = lsat_map if (q_msat_map eq 1) then m_map = msat_map ; determine if prepping is necessary. l_prep=get_history(l_idx1,'XRT_PREP',found=l_found) m_prep=get_history(m_idx1,'XRT_PREP',found=m_found) s_prep=get_history(s_idx1,'XRT_PREP',found=s_found) ; If you prep one you must prep all. if ((l_found ne s_found) or (l_found ne m_found)) then begin print, $ 'Long/Medium/Short images must have same type of processing... Returning.' return endif ; If you prep the long you must provide the saturated pixel map. if (l_found eq 1 and q_lsat_map eq 0) then begin print,'Must use the lsat_map keyword for xrt_prepped data... Returning.' return endif if (m_found eq 1 and q_msat_map eq 0) then begin print,'Must use the msat_map keyword for xrt_prepped data... Returning.' return endif ; Make sure that l_da & lsat_map dimensions match. if (l_found eq 1 or q_lsat_map eq 1) then begin sz_map=size(l_map,/dimensions) sz_lda=size(l_da1,/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 if (m_found eq 1 or q_msat_map eq 1) then begin sz_map=size(m_map,/dimensions) sz_mda=size(m_da1,/dimensions) ndim1=n_elements(m_da1) ndim2=n_elements(m_map) if (ndim1 ne ndim2) then begin print,'Error medium image and sat_map array 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) m_norm=get_history(m_idx1,'XRT_RENORMALIZE',found=nm_found) s_norm=get_history(s_idx1,'XRT_RENORMALIZE',found=ns_found) ; Perhaps they normalized them themselves. if (l_found eq 1 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 exposures.' endif if (m_found eq 1 and nm_found eq 0) then begin m_da1=temporary(m_da1)/get_xrt_expdur(m_idx1, /sec) m_idx1.e_etim=1000000 ; [microsec] for consistancy with XRT_PREP ; & returned index if (loud eq 1) then print,'Renormalizing the exposures.' endif if (s_found eq 1 and nl_found eq 0) then begin s_da1=temporary(s_da1)/get_xrt_expdur(s_idx1, /sec) endif ; Use xrt_prep on images... If needed. if (s_found eq 0) then begin if (loud eq 1) then print,'Prepping short for composite' ; s_map=1 s_unc=1 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 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, $ uncert_map=s_unc, $ n_grade_pixels=ng_pix, verbose=verbose, quiet=quiet n_sat_pixels=ng_pix(0) ; extract # of saturated pixels endif if (m_found eq 0) then begin if (loud eq 1) then print,'Prepping medium for composite' gm_map=1 ngm_pix=1 m_unc=1 xrt_prep, temporary(m_idx1), temporary(m_da1), m_idx1, m_da1, $ despike_despot=despike_despot, only_despot=only_despot, $ only_despike=only_despike, /float, /normalize, $ grade_map=gm_map,n_grade_pixels=ngm_pix, verbose=verbose, $ uncert_map=m_unc, $ quiet=quiet n_sat_m = ngm_pix[0] ; extract # of saturated pixels m_map = (fix(gm_map)/2. ne fix(fix(gm_map)/2.)) ; extract saturated pixels ; (odd values in gm_map) ;***** treatment of the MEDIUM 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(m_idx,'PIXCNT') eq 1 then $ if m_idx.PIXCNT lt n_elements(m_da) then begin ss_mfrm = where(m_da eq 0,nmiss) if nmiss gt 0 then m_map[ss_mfrm] = 1 endif endif else begin ; input is processed (Level 1) data ss_mfrm = where(m_da eq -999,nmiss) if nmiss gt 0 then m_map[ss_mfrm] = 1 endelse if (l_found eq 0) then begin if (loud eq 1) then print,'Prepping long for composite' gl_map=1 ngl_pix=1 l_unc=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=gl_map, n_grade_pixels=ngl_pix, verbose=verbose, $ uncert_map=l_unc, $ quiet=quiet n_sat_l=ngl_pix(0) ; extract # of saturated pixels l_map=(fix(gl_map)/2. ne fix(fix(gl_map)/2.)) ; extract saturated pixels ; (odd values in gl_map) ;***** 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 begin ss_mfrl = where(l_da eq 0,nmiss) if nmiss gt 0 then l_map[ss_mfrl] = 1 endif endif else begin ss_mfrl = where(l_da eq -999,nmiss) if nmiss gt 0 then l_map[ss_mfrl] = 1 endelse ;***** treatment of the SHORT missing frames here ***** if tag_exist(s_idx,'PIXCNT') eq 1 then $ if s_idx.PIXCNT lt n_elements(s_da) then begin ss_mfrs = where(s_da1 eq -999,nmiss) if nmiss gt 0 then s_da1[ss_mfrs] = m_da1[ss_mfrs] ; in case same pixel is bad in both short and medium exposures: ss_mfrs = where(s_da1 eq -999,nmiss) if nmiss gt 0 then s_da1[ss_mfrs] = l_da1[ss_mfrs] endif if (l_found ne 0) 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) m_da_out = xrt_despike2(m_da1,m_idx1) s_da_out = xrt_despike2(s_da1,s_idx1) l_da1 = l_da_out m_da1 = m_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,m_idx1,m_da1,m_idx_out,m_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 m_da1 = m_da_out m_idx1 = m_idx_out s_da1 = s_da_out s_idx1 = s_idx_out endif endif ; Composite creation ************************************ if q_lsat_map eq 1 then begin ; check for case of input satmap=0 if total(l_map gt 0) then n_sat_l = 1 else n_sat_l = 0 endif if q_msat_map eq 1 then begin ; check for case of input satmap=0 if total(m_map gt 0) then n_sat_m = 1 else n_sat_m = 0 endif ; Use sat_map to replace saturated pixels of long exposure. if (n_sat_m gt 0) then begin ; JDS 8/15/2014: added criterion that the short exposure data are > 0 ; to prevent replacing saturated pixels with bad data ssm = where((m_map eq 1) and (s_da1 > 0.),ns) if ns gt 0 then m_da1[ssm] = (s_da1[ssm] > min(m_da1)) if ns gt 0 then m_unc[ssm] = s_unc[ssm] endif else begin message,'No saturated pixels to replace - in middle image',/continue endelse if (n_sat_l gt 0) then begin ; JDS 8/15/2014: added criterion that the medium exposure data are > 0 ; to prevent replacing saturated pixels with bad data ssl = where((l_map eq 1) and (m_da1 > 0.),ns) if ns gt 0 then l_da1[ssl] = (m_da1[ssl] > min(l_da1)) if ns gt 0 then l_unc[ssl] = m_unc[ssl] c_da = l_da1 uncert_map = l_unc endif else begin c_da = l_da1 uncert_map = l_unc message,'No saturated pixels to replace - in long image',/continue endelse ;; ----- Added 11-June-2007 -----;; ; handeling the edge pixels ; using Saar's method of filtering whole image and then replacing, ; not just filtering bloom pixels. ; P.R. Jibben -He also calcuated shift better, only in vertical direction. ; Only doing this for long image siz = size(c_da) ny = siz[2] donut=[[1,1,1],[1,0,1],[1,1,1]] / 8. Result=CONVOL(c_da,donut) sl_map1=shift(l_map,[0,1]) sl_map1[*,0]=0 sl_map2=shift(fix(l_map),[0,-1]) sl_map2[*,ny-1]=0 bloomers=temporary(sl_map1)+temporary(sl_map2) bloomers=(bloomers eq 1) replace_pix=where(bloomers eq 1,cnt) bloom_tag='' if (cnt gt 0) then begin nbloom=strcompress(cnt,/rem) c_da[replace_pix]=Result[replace_pix] bloom_tag='Replaced ' + nbloom + ' bloom pixels in Composite.' endif ;generating the c_idx c_idx = l_idx1 ; adding new header keywords that provide easily accessible information ; on the filenames for the long and short exosure image files lfname = xrt_index2filename(l_idx) mfname = xrt_index2filename(m_idx) sfname = xrt_index2filename(s_idx) c_idx = add_tag2(c_idx,lfname,'LNGFNAME',index='OBS_MODE') c_idx = add_tag2(c_idx,mfname,'MEDFNAME',index='LNGFNAME') c_idx = add_tag2(c_idx,sfname,'SRTFNAME',index='MEDFNAME') ;; ----- Added 11-June-2007 ----- ;; ;Updating the index data level ;level2 data c_idx.data_lev=2 tagval0='' change_tag_value,c_idx,tagval0,'HISTORY' tagval0=[tagval0,progname + ' ' + progver + ': '] 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'] long_hist=get_history(l_idx1) medium_hist=get_history(m_idx1) short_hist=get_history(s_idx1) tagval0=[tagval0,bloom_tag,'Long image history --',long_hist, $ 'Medium image history --',medium_hist, $ 'Short image history --',short_hist] update_history,c_idx,tagval0,/noroutine RETURN END ;======================================================================