;+ Helper routine; main program starts below function IRIS_DUST_AMOEBA, shifts ;+ ; ; Helper routine called by IRIS_DUSTBUSTER when trying to optimize the alignment of the ; dust mask. This routine returns the total pixel value of all pixels that are included ; in the dust mask; it is called iteratively using IDL's AMOEBA minimizer to find the ; X and Y offsets that best align the dust mask with the dust in the image cube. ; ;- common DUSTBUSTER_FIT_AMOEBA_CB, l2data, bpx2, bpy2, bpxcen, bpycen, nx, ny, zaddress, maxshift nbpx2 = bpx2 + shifts[0] nbpy2 = bpy2 + shifts[1] maxradius = 880 ; more than this many arcsec from disk center and we ignore it bpaddress = nbpx2 + LONG(nx) * nbpy2 + zaddress keepers = WHERE( (bpx2 ge 0) and (bpx2 lt nx) and $ (bpy2 ge 0) and (bpy2 lt ny) and $ (ABS(bpxcen) le maxradius) and (ABS(bpycen) le maxradius) and $ (l2data[bpaddress] ne -200) and FINITE(l2data[bpaddress]), numkeep ) if numkeep eq 0 then begin BOX_MESSAGE, ['IRIS_DUST_AMOEBA: No bad pixels!', 'Returning 1e8'] RETURN, 1e8 endif bpaddress = bpaddress[keepers] bpaddress = bpaddress[UNIQ(bpaddress, SORT(bpaddress))] result = MEAN(l2data[bpaddress]) if (MAX(ABS(shifts)) gt (maxshift+1)) then result += 1e5 RETURN, result end function IRIS_DUSTBUSTER, l2index, l2data, $ bpaddress, clean_values, $ smoothsize = smoothsize, meanfill = meanfill, noblurfill = noblurfill, $ test = test, list_only = list_only, pstop = pstop, run = run, loud = loud, $ slow = slow, slit = slit, shift_badpix=shift_badpix, noshift = noshift ;+ ; ; Cleans up the dust on the level 2 FUV SJI movies by picking good pixels from ; frames that are adjacent in time (so it works better on rasters, not sit-and- ; stares, though it's safe to run on anything). ; ; INPUTS: ; l2index - Index structure from a L2 SJI FITS file (read in using read_iris_l2) ; l2data - Data cube from a L2 SJI FITS file (read in using read_iris_l2) ; Note that it should be of FLOAT type in order to avoid some funny ; features like popping in data that should be out of the frame ; ; OUTPUTS: ; bpaddress - An N-element LONG array specifying the (1-D) address of the ; bad pixels in the L2 data cube ; clean_values- An N-element FLOAT array giving the pixel values to poke into the ; bad pixels ; ; RETURN: ; clean_data - A data cube with the same dimensions as l2data, with the dust ; busted (note that if the /list_only keyword is set, then only ; a scalar integer is returned (the user can do the replacement ; separately) ; ; KEYWORDS: ; smoothsize - Set to the numebr of pixels over which to smooth the bad pixel ; mask. Default is not to use smoothing (it used to default to 6, ; back when we had a rounding error in the dust mask remapping) ; /meanfill - Plug unbustable dust specks (e.g. those at the edges of the ; first and last frames of the raster or all dust in sit-n-stares) ; with the overall median pixel value ; /slow - If set, then unbustable pixels are filled with a median-filtered ; value instead of a straight smoothed value. Takes much longer ; (minutes instead of seconds), but avoids an issue with snow ; off the limb ; /noblurfill - Unless this keyword is set, it will plug unbustable dust ; specks (e.g. those at the edges of the first and last frames ; of the raster or all dust in sit-n-stares with a locally ; appropriate pixel value found by blurring the image) ; /test - Set to generate some diagnostic plots/images (can set to an ; integer greater than one to run plots on different frames in ; the cube) ; /pstop - Stop for debugging ; /run - Print the elapsed run time ; /list_only - See description in { RETURN: clean_data } above ; /loud - Print some diagnostic information (not much) ; /slit - If set, then the slit is removed from the SJI as well as the dust ; /noshift - If set, then the dust pattern is not shifted to match the dust in the ; image. By default, the dust alignment is optimized using IDL's ; AMOEBA minimizer ; ; HISTORY: ; PB 2014-05-01 - Initial version ; PB 2014-12-22 - Updates including blur filling for sit-n-stares ; SLF 2015-10-21 - Adjust default SMOOTHSIZE for >=4 spatial x 2spectra summing (divide Paul defaults by 4) ; add .HISTORY record w/smoothsize applied ; PB 2020-06-12 - Added shifting of dust mask using AMOEBA; ; Fixed rounding error that caused gap in dust mask; ; Remove smoothing of dust mask by default ; ; EXAMPLE: ; (Note: it works best on many-step or coarse rasters; for sit-n-stares or ; narrow rasters, it fills in with spatial blurring instead of temporal ; blurring. It shouldn't ever really hurt anything, though...) ; Four-step sparse raster: ; IDL> cd, '/irisa/data/level2/2014/04/26/20140426_010036_3820257468/' ; IDL> read_iris_l2, 'iris_l2_20140426_010036_3820257468_SJI_1400_t000.fits', l2index, l2data ; % READFITS: Now reading 742 by 773 by 940 array ; IDL> cdata = IRIS_DUSTBUSTER(l2index, l2data, bpaddress, clean_values, /run) ; IRIS_DUSTBUSTER run time = 3.0 ; ; Four-step sparse raster: ; IDL> cd, '/irisa/data/level2/2014/04/22/20140422_184749_3820113694/' ; IDL> read_iris_l2, 'iris_l2_20140422_184749_3820113694_SJI_1400_t000.fits', l2index, l2data ; % READFITS: Now reading 738 by 387 by 512 array ; IDL> cdata = IRIS_DUSTBUSTER(l2index, l2data, bpaddress, clean_values, /run) ; IRIS_DUSTBUSTER run time = 1.2 ; IDL> ploop, alog10(l2data>1), alog10(cdata>1), min = 0, max = 4 ; ;- t0 = SYSTIME(/sec) if SIZE(l2data, /type) ne 4 then PRINT, 'IRIS_DUSTBUSTER works best on FLOAT data...' if KEYWORD_SET(meanfill) then noblurfill = 1 if n_elements(shift_badpix) eq 2 then begin xtweak=shift_badpix[0] ytweak=shift_badpix[1] endif else begin ; PB defaults xtweak = 0. ; Hard-coded manual adjustment of bad pixel mask alignment ytweak = -0.5 endelse ; Set up for the optimization of the bad pixel mask if not KEYWORD_SET(noshift) then begin common DUSTBUSTER_FIT_AMOEBA_CB, l2data_cb, bpx2, bpy2, bpxcen, bpycen, nx, ny, zaddress, maxshift maxshift = 7 ; Maximum allowable shift, in pixels endif ;+++++++++++++++++++++++++ ; Pull some information ; out of the header and ; set up defaults ;------------------------- img_path = l2index[0].tdesc1 imgpath = STRJOIN(STRSPLIT(/extract, /regex, img_path, '_')) sumspat = l2index[0].sumspat sumsptrl = l2index[0].sumsptrl imsize = SIZE(l2data) nx = imsize[1] ny = imsize[2] if imsize[0] eq 0 then nz = 1 else nz = imsize[3] ndat = nx * ny * nz offspat = (sumspat - 1.0) / (sumspat * 2.0) offsptrl = (sumsptrl - 1.0) / (sumsptrl * 2.0) exptime = l2index.exptime ; Set up blurring of bad pixel mask, depending on summing and img_path if not KEYWORD_SET(smoothsize) then begin case STRING(sumspat, ',', img_path, form = '(i1,a1,a8)') of '1,SJI_1330' : smoothsize = 0 '1,SJI_1400' : smoothsize = 0 '1,SJI_2796' : smoothsize = 0 '1,SJI_2832' : smoothsize = 4 '2,SJI_1330' : smoothsize = 0 '2,SJI_1400' : smoothsize = 0 '2,SJI_2796' : smoothsize = 0 '2,SJI_2832' : smoothsize = 6 '4,SJI_1330' : smoothsize = 0 '4,SJI_1400' : smoothsize = 0 '4,SJI_2796' : smoothsize = 0 '4,SJI_2832' : smoothsize = 8 '8,SJI_1330' : smoothsize = 0 '8,SJI_1400' : smoothsize = 0 '8,SJI_2796' : smoothsize = 0 '8,SJI_2832' : smoothsize = 10 else: smoothsize = 2 endcase smoothsize=smoothsize/([1,4])(sumspat ge 4 and sumsptrl ge 2) > 4 ; slf, 21-oct-2015 endif update_history,l2index,'SMOOTHSIZE='+strtrim(smoothsize,2),/caller ; Slit center coordinates in L2 frame slitx = l2index.sltpx1ix - 1 slity = l2index.sltpx2ix - 1 ;+++++++++++++++++++++++++ ; Read in the bad pixel ; list for the desired ; channel and time ;------------------------- dummyhead = CREATE_STRUCT('IMG_PATH', l2index[0].tdesc1, 'T_OBS', l2index[0].date_obs) badpix_list = IRIS_PREP_GET_BADPIX(dummyhead) if KEYWORD_SET(slit) then begin slitpix_list = IRIS_PREP_GET_SLITPIX(dummyhead) badpix_list = SETUNION(badpix_list, slitpix_list) if not KEYWORD_SET(noshift) then begin BOX_MESSAGE, ['IRIS_DUSTBUSTER: No mask shifting ', $ ' when slit removal is requested...', $ ' Setting /noshift'] noshift = 1 endif endif ; Expand the bad pixel mask by smoothing if smoothsize gt 0 then begin if smoothsize eq 1 then smoothsize = 4 mask = FLTARR(2072, 1096) mask[badpix_list] = 1 ; smask = GAUSS_SMOOTH(mask, smoothsize, /edge_mirr) gt 0 smask = SMOOTH(mask, smoothsize) gt 0 badpix_list = WHERE(smask gt 0) endif else smoothsize = 1 nbad = N_ELEMENTS(badpix_list) ;+++++++++++++++++++++++++ ; Map the bad pixel list ; to give X and Y distance ; from slit center ;------------------------- ; Bad pixel coordinates in SJI frame bpx = (badpix_list mod 2072) / sumsptrl bpy = (badpix_list / 2072) / sumspat bpx = bpx + xtweak ; Manual adjustment of bad pixel mask bpy = bpy + ytweak ; ; Slit center coordinates in SJI frame pdb = IRIS_MK_POINTDB(ver = pdbver) pdb_tags = TAG_NAMES(pdb) sji_wv_nm = STRMID(imgpath, 3, 3) if sji_wv_nm eq '160' then sji_wv_nm = 'MIR' if sji_wv_nm eq '500' then sji_wv_nm = 'FSI' slitx_0 = (pdb.(WHERE(pdb_tags eq 'CPX1_' + sji_wv_nm))-1)/sumsptrl + offsptrl slity_0 = (pdb.(WHERE(pdb_tags eq 'CPX2_' + sji_wv_nm))-1)/sumspat + offspat ; Bad pixel to slit center distance in SJI frame bpdx = bpx - slitx_0 bpdy = bpy - slity_0 ; Bad pixel to slit center distance in L2 frame slitd = SQRT(bpdx^2d + bpdy^2d) slitq = ATAN(bpdx, bpdy) scale = pdb.(WHERE(pdb_tags eq 'CDLT_' + sji_wv_nm)) roll = pdb.(WHERE(pdb_tags eq 'BE_' + sji_wv_nm)) sref = pdb.cdlt_1p5 rref = 0. mag = scale / sref roll = roll - rref slitx_2 = (slitd * mag) * SIN(slitq - roll * !pi / 180.) slity_2 = (slitd * mag) * COS(slitq - roll * !pi / 180.) ; Bad pixel coordinates in L2 frame ;bpx2 = ROUND( REBIN(slitx_2, nbad, nz) + REBIN(TRANSPOSE(slitx), nbad, nz) ) ;bpy2 = ROUND( REBIN(slity_2, nbad, nz) + REBIN(TRANSPOSE(slity), nbad, nz) ) ;zaddress = nx * LONG(ny) * REBIN(REFORM(LINDGEN(nz), 1, nz), nbad, nz) ; ; Note: The rounding that we used to use resulted in gaps in the bad pixel map in ; some cases. Instead we're going to combine floor and ceil versions, and later on ; drop overlapping bad pixels bpx2 = [FLOOR( REBIN(slitx_2, nbad, nz) + REBIN(TRANSPOSE(slitx), nbad, nz) ) , $ CEIL( REBIN(slitx_2, nbad, nz) + REBIN(TRANSPOSE(slitx), nbad, nz) ) ] bpy2 = [FLOOR( REBIN(slity_2, nbad, nz) + REBIN(TRANSPOSE(slity), nbad, nz) ) , $ CEIL( REBIN(slity_2, nbad, nz) + REBIN(TRANSPOSE(slity), nbad, nz) ) ] zaddress = [nx * LONG(ny) * REBIN(REFORM(LINDGEN(nz), 1, nz), nbad, nz), $ nx * LONG(ny) * REBIN(REFORM(LINDGEN(nz), 1, nz), nbad, nz) ] ; Bad pixel coordinates in helioprojective frame (arcsec x/y) bpxcen = bpx2*0 bpycen = bpy2*0 for frame = 0, N_ELEMENTS(l2index)-1 do begin bpxcen[*,frame] = (bpx2[*,frame] - l2index[frame].crpix1) * $ l2index[frame].cdelt1 + l2index[frame].crval1 bpycen[*,frame] = (bpy2[*,frame] - l2index[frame].crpix2) * $ l2index[frame].cdelt2 + l2index[frame].crval2 endfor ; Check for shifts in the bad pixel pattern if not KEYWORD_SET(noshift) then begin l2data_cb = l2data optimal_shift = ROUND(AMOEBA(1e-5, function_name="IRIS_DUST_AMOEBA", scale=1, p0=[0,0])) msg = 'IRIS_DUSTBUSTER: Optimal shift = ' + STRTRIM(optimal_shift, 2) if MAX(ABS(optimal_shift)) le maxshift then begin bpx2 = bpx2 + optimal_shift[0] bpy2 = bpy2 + optimal_shift[1] endif else begin msg = [msg, "Exceeds max shift; not applied..."] endelse BOX_MESSAGE, msg endif ;+++++++++++++++++++++++++ ; Look up the replacement ; pixel values for the ; bad pixels ;------------------------- ; Convert to 1-D address, and throw out anything out-of-bounds bpaddress = bpx2 + LONG(nx) * bpy2 + zaddress if KEYWORD_SET(pstop) then STOP keepers = WHERE( (bpx2 ge 0) and (bpx2 lt nx) and $ (bpy2 ge 0) and (bpy2 lt ny) and $ (l2data[bpaddress] ne -200) and FINITE(l2data[bpaddress]), numkeep ) if numkeep eq 0 then begin BOX_MESSAGE, ['IRIS_DUSTBUSTER: No bad pixels!', 'Returning input data cube...'] if KEYWORD_SET(pstop) then STOP RETURN, l2data endif bpaddress = bpaddress[keepers] bpsort = UNIQ(bpaddress, SORT(bpaddress)) bpaddress = bpaddress[bpsort] ; Note non-unique bpaddress due to smoothing and rounding numkeep = N_ELEMENTS(bpaddress) bphist = HISTOGRAM(bpaddress, min = 0) bpz = bpaddress / (nx * ny) ; Frame number bpexp = exptime[bpz] ; Find good values in adjacent (in time) pixels adjacents = LONG([-2,-1,1,2]) * nx * ny ; 1-d address offset for adjacent frames numadj = N_ELEMENTS(adjacents) candidates = REBIN(TRANSPOSE(adjacents), numkeep, numadj) + REBIN(bpaddress, numkeep, numadj) bad_candidates = WHERE( (candidates lt 0) or (candidates gt ndat) or $ (bphist[candidates] ne 0) or (l2data[candidates] lt -199), numbad) if numbad gt 0 then candidates[bad_candidates] = 0 replacement_values = FLOAT(l2data[candidates]) ; Need to float in case a bytarr is passed in, e.g. intscaled if numbad gt 0 then replacement_values[bad_candidates] = !values.f_nan ; Handle AEC by normalizing by exposure time candidate_frame = candidates / (nx * ny) ; frame number candidate_exp = exptime[candidate_frame] replacement_values = replacement_values / candidate_exp clean_values = MEDIAN(replacement_values, dim = 2) * bpexp ; Check for pixels where there are no good adjacent replacements; if there are ; any, then replace them with the global (raster) median goodmask2d = replacement_values eq replacement_values goodmask = TOTAL(goodmask2d, 2) no_adjacent = WHERE(goodmask eq 0, numno) if KEYWORD_SET(loud) then PRINT, numno, ' pixels with no good replacements' if (numno gt 0) and (KEYWORD_SET(meanfill)) then begin allmedian = MEDIAN(l2data>0) clean_values[no_adjacent] = allmedian endif if (numno gt 0) and (not KEYWORD_SET(noblurfill)) then begin blurdat = l2data ; Wasteful, but probably ok outframe = WHERE(l2data lt -199, numout) if numout gt 0 then blurdat[outframe] = !values.f_nan blurdat[bpaddress] = !values.f_nan blursmooth = smoothsize * 4. ; Decide how to do the blurring for replacement pixel values. If the image is ; too small, it can't be done at all. If it is small in Z, or the user requests ; /slow, then it is done with a frame-by-frame median value. Otherwise, use a ; global smooth value doblur = 1 if nz le blursmooth then doblur = 2 if KEYWORD_SET(slow) then doblur = 2 if MIN([nx,ny]) le blursmooth then doblur = 0 if KEYWORD_SET(pstop) then STOP case doblur of 2 : for i = 0, nz - 1 do blurdat[*,*,i] = MEDIAN(blurdat[*,*,i], blursmooth) > 0 1 : blurdat = SMOOTH(blurdat, blursmooth, /nan) > 0 else: PRINT, 'IRIS_DUSTBUSTER: cannot use blurfill (image too small?)' endcase clean_values[no_adjacent] = blurdat[bpaddress[no_adjacent]] endif ; If requested, replace the bad values and return the clean array if KEYWORD_SET(list_only) then begin cleandata = 1 endif else begin cleandata = l2data cleandata[bpaddress] = clean_values endelse ; Do some diagnostics to ensure that we're busting the right pixels if KEYWORD_SET(test) then begin TVLCT, rr, gg, bb, /get ; Blink an image and a mask IRIS_LCT, img_path inframe = l2data[*,*,test-1] outframe = cleandata[*,*,test-1] diffframe = 255 * (inframe ne outframe) if SIZE(l2data, /type) ne 1 then begin inframe = IRIS_INTSCALE(inframe, l2index) outframe = IRIS_INTSCALE(outframe, l2index) endif PLOOP, [[[inframe]], [[outframe]], [[diffframe]]], wait = 0.1, $ title = ['Input', 'Output', 'Input ne Output'] TVLCT, rr, gg, bb endif if KEYWORD_SET(run) then PRINT, STRING('IRIS_DUSTBUSTER run time =', $ SYSTIME(/sec) - t0, format = '(a40, f10.1)') if KEYWORD_SET(pstop) then STOP RETURN, cleandata end